text:
c     Program to amalgamate USAF files in second update
c     AJ_Kettle, 21Apr2021

c     size of source directories
c     2019: 54 750 032
c     2020: 45 882 528
c     2021:  6 624 176

c     22Apr2021: full stnlist for 2019-2021 15577 stn
c     27Apr2021: full run amalgamation files; c/c=16.4/117.0min
c     27Apr2021: full run amalgamation files; c/c=16.14/86.30min

c************************************************************************
c2-------10--------20--------30--------40--------50--------60--------70--------80
      IMPLICIT NONE
c************************************************************************
c     Declare variables

      INTEGER             :: i,j,k,ii,jj,kk

      REAL                :: f_time_st,f_time_en

      REAL                :: f_deltime_cpu_s
      REAL                :: f_deltime_clock_s

      CHARACTER(LEN=8)    :: s_date_st
      CHARACTER(LEN=10)   :: s_time_st
      CHARACTER(LEN=5)    :: s_zone_st
      INTEGER             :: i_values_st(8)
      CHARACTER(LEN=8)    :: s_date_en
      CHARACTER(LEN=10)   :: s_time_en
      CHARACTER(LEN=5)    :: s_zone_en
      INTEGER             :: i_values_en(8)

      REAL                :: f_clock_st_s,f_clock_en_s
c*****
      CHARACTER(LEN=300)  :: s_directory_filelist
      CHARACTER(LEN=300)  :: s_directory_dataroot1
      CHARACTER(LEN=300)  :: s_directory_dataroot2
      CHARACTER(LEN=300)  :: s_directory_dataroot3

      CHARACTER(LEN=300)  :: s_directory_assemstns
      CHARACTER(LEN=300)  :: s_directory_diagnostics

      CHARACTER(LEN=300)  :: s_subdirectory_sorted

      CHARACTER(LEN=300)  :: s_filelist
      CHARACTER(LEN=300)  :: s_file_liststn
      CHARACTER(LEN=300)  :: s_file_linecount

c*****
      INTEGER,PARAMETER   :: l_rgh_stn=20000
      INTEGER             :: l_stn2019
      INTEGER             :: l_stn2020
      INTEGER             :: l_stn2021
      CHARACTER(LEN=32)   :: s_vec_stnlist2019(l_rgh_stn)
      CHARACTER(LEN=32)   :: s_vec_stnlist2020(l_rgh_stn)
      CHARACTER(LEN=32)   :: s_vec_stnlist2021(l_rgh_stn)
c*****
      INTEGER             :: l_amal
      CHARACTER(LEN=32)   :: s_vec_stnlist_amal(l_rgh_stn)
      INTEGER             :: i_mat_stnlist_flag(l_rgh_stn,3)
c*****
      CHARACTER(LEN=300)  :: s_command
      INTEGER             :: io
c************************************************************************
      print*,'start am1'
c************************************************************************
c     Find start time
      CALL CPU_TIME(f_time_st)

c     Find date & time
      CALL DATE_AND_TIME(s_date_st,s_time_st,s_zone_st,i_values_st)
c************************************************************************
c************************************************************************
c     Directory of files

c     Note root is in GWS for update1 & scratch-pw for update2
      s_directory_dataroot1=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/'//
     + 'Y2019/'
      s_directory_dataroot2=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/'//
     + 'Y2020/'
      s_directory_dataroot3=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/'//
     + 'Y2021/'

      s_directory_assemstns=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/Y_amal/'//
     + 'Data_5sorted/'
      s_directory_diagnostics=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/Y_amal/'//
     + 'Data_5diag/'

      s_subdirectory_sorted='Data_4sorted/'

      s_filelist='z_filelist.dat'
      s_file_liststn='list_stn_2019_2020_2021.dat'
      s_file_linecount='linecount_2019_2020_2021.dat'

c      s_file_diag_fullstn
c************************************************************************
c     Erase files output directory

      s_command=
     +  'rm '//TRIM(s_directory_assemstns)//'*'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

      s_command=
     +  'rm '//TRIM(s_directory_diagnostics)//'*'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

c      STOP 'am1'
c************************************************************************
c     2019 - Create datafile list
      s_command=
     +  '(cd '//TRIM(s_directory_dataroot1)//
     +  TRIM(s_subdirectory_sorted)//
     +  ' && ls) > '//
     +  TRIM(s_filelist)
c      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

c     Read in data file lists
      CALL readin_stnlist(s_filelist,
     +  l_rgh_stn,l_stn2019,s_vec_stnlist2019)

c     Remove data file lists
      s_command='rm '//TRIM(s_filelist)
c      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)
c*****
c     2020 - Create datafile list
      s_command=
     +  '(cd '//TRIM(s_directory_dataroot2)//
     +  TRIM(s_subdirectory_sorted)//
     +  ' && ls) > '//
     +  TRIM(s_filelist)
c      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

c     Read in data file lists
      CALL readin_stnlist(s_filelist,
     +  l_rgh_stn,l_stn2020,s_vec_stnlist2020)

c     Remove data file lists
      s_command='rm '//TRIM(s_filelist)
c      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)
c*****
c     2021 - Create datafile list
      s_command=
     +  '(cd '//TRIM(s_directory_dataroot3)//
     +  TRIM(s_subdirectory_sorted)//
     +  ' && ls) > '//
     +  TRIM(s_filelist)
c      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

c     Read in data file lists
      CALL readin_stnlist(s_filelist,
     +  l_rgh_stn,l_stn2021,s_vec_stnlist2021)

c     Remove data file lists
      s_command='rm '//TRIM(s_filelist)
c      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)
c*****
c************************************************************************
c     Create diagnostics files
      CALL find_single_bigvector2(l_rgh_stn,
     +  l_stn2019,s_vec_stnlist2019,l_stn2020,s_vec_stnlist2020,
     +  l_stn2021,s_vec_stnlist2021,
     +  l_amal,
     +  s_vec_stnlist_amal,i_mat_stnlist_flag)
c************************************************************************
c     Export station diagnostics list
      CALL export_stnlist(s_date_st,s_time_st,
     +  s_directory_diagnostics,s_file_liststn,
     +  l_rgh_stn,l_amal,
     +  s_vec_stnlist_amal,i_mat_stnlist_flag)
c************************************************************************
c     Create amalgamated set of data files

      CALL combine_3datasets(s_date_st,s_time_st,
     +  s_directory_dataroot1,
     +  s_directory_dataroot2,s_directory_dataroot3,
     +  s_subdirectory_sorted,
     +  s_directory_assemstns,s_directory_diagnostics,
     +  s_file_linecount,
     +  l_rgh_stn,l_amal,
     +  s_vec_stnlist_amal,i_mat_stnlist_flag)

c************************************************************************
c************************************************************************
c     Find end time
      CALL CPU_TIME(f_time_en)

c     Find date & time
      CALL DATE_AND_TIME(s_date_en,s_time_en,s_zone_en,i_values_en)
c************************************************************************
      print*,'s_date_en=',TRIM(s_date_en)
      print*,'s_time_en=',TRIM(s_time_en)
      print*,'i_values_st',(i_values_st(i),i=1,8)
      print*,'i_values_en',(i_values_en(i),i=1,8)

      f_clock_st_s=FLOAT(i_values_st(3))*24.0*60.0*60.0+
     +  FLOAT(i_values_st(5))*60.0*60.0+FLOAT(i_values_st(6))*60.0+
     +  FLOAT(i_values_st(7))
      f_clock_en_s=FLOAT(i_values_en(3))*24.0*60.0*60.0+
     +  FLOAT(i_values_en(5))*60.0*60.0+FLOAT(i_values_en(6))*60.0+
     +  FLOAT(i_values_en(7))

      print*,'f_clock_st_s,f_clock_en_s=',f_clock_st_s,f_clock_en_s

      f_deltime_clock_s=f_clock_en_s-f_clock_st_s
c************************************************************************
      f_deltime_cpu_s=f_time_en-f_time_st
      print*,'f_deltime_cpu_s,f_deltime_cpu_min=',
     +  f_deltime_cpu_s,f_deltime_cpu_s/60.0
      print*,'f_deltime_clock_s,f_deltime_clock_min=',
     +  f_deltime_clock_s,f_deltime_clock_s/60.0
c************************************************************************
      print*,'successfully finished am1'

      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


SUBROUTINE DiffXEven(var,varX,order,Nx,Ny,Nz,DiffStencils,alpha)
!
! Explicit order differentiation in X-direction.
!
! No conditional loops.
!
! By Allan P. Engsig-Karup.
USE Precision
USE DataTypes
IMPLICIT NONE
TYPE (Diff_def) :: DiffStencils
INTEGER :: Nx, Ny, Nz, rank, order, alpha, i,j,k
REAL(KIND=long), DIMENSION(Nz,Nx,Ny) :: var, varX
REAL(KIND=long), DIMENSION(2*alpha+1) :: Stencil
rank = 2*alpha+1
DO i = 1, alpha
	Stencil = DiffStencils%StencilX(i,1:rank,order)
	DO j = 1, Ny
		DO k = 1, Nz
			varX(k,i,j) = DOT_PRODUCT(Stencil,var(k,1:rank,j))
		END DO
	END DO
END DO
DO i = alpha+1, Nx-alpha
	Stencil = DiffStencils%StencilX(i,1:rank,order)
	DO j = 1, Ny
		DO k = 1, Nz
			varX(k,i,j) = DOT_PRODUCT(Stencil,var(k,i-alpha:i+alpha,j))
		END DO
	END DO
END DO
DO i = Nx-alpha+1, Nx
	Stencil = DiffStencils%StencilX(i,1:rank,order)
	DO j = 1, Ny
		DO k = 1, Nz
			varX(k,i,j) = DOT_PRODUCT(Stencil,var(k,Nx+1-rank:Nx,j))
		END DO
	END DO
END DO
END SUBROUTINE DiffXEven


text:
      subroutine readgmv_binary(ifile,ierror)
c
c #####################################################################
c
c     purpose -
c
c        create a gmv formatted file from the current mesh object
c
c     input arguments -
c
c        none
c
c     output arguments -
c
c        none
c
c     change history -
c
C        $Log: readgmv_binary.f,v $
C        Revision 2.00  2007/11/09 20:03:59  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.9   12 Jun 2003 11:53:06   tam
CPVCS    initialized ctype and cmoattnam to 32 blanks for opt SGI
CPVCS
CPVCS       Rev 1.8   24 Mar 2003 14:38:00   gable
CPVCS    Recognize AMR mesh variables itetpar, itetkid, itetlev as VINT.
CPVCS
CPVCS       Rev 1.7   18 Jul 2001 15:44:50   dcg
CPVCS    get rid of all lines having to do with random io calls
CPVCS    i.e call cposition, iadr calculation
CPVCS
CPVCS       Rev 1.5   29 Nov 2000 10:08:56   dcg
CPVCS    check if geometry exists
CPVCS
CPVCS       Rev 1.4   Wed Apr 12 15:06:16 2000   dcg
CPVCS    change nulls to blanks after cread ( nulls come from c io)
CPVCS
CPVCS       Rev 1.3   Wed Apr 05 13:34:56 2000   nnc
CPVCS    Minor source modifications required by the Absoft compiler.
CPVCS
CPVCS       Rev 1.2   Wed Feb 02 13:43:08 2000   dcg
CPVCS
CPVCS       Rev 1.1   Mon Jan 31 13:31:08 2000   dcg
CPVCS
CPVCS       Rev 1.0   Mon Jan 31 13:25:36 2000   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.16   Tue Nov 30 17:08:12 1999   jtg
CPVCS    slightly cleaned up, nothing major
CPVCS
CPVCS       Rev 1.15   Thu Aug 05 14:03:34 1999   dcg
CPVCS    allow line type mesh object elements
CPVCS
CPVCS       Rev 1.14   Tue Jul 06 19:29:40 1999   jtg
CPVCS    modifed call to geniee (put at end)
CPVCS
CPVCS       Rev 1.13   Fri Jan 22 15:35:16 1999   dcg
CPVCS    change itetclr(ielements)=imat to itetclr(ielements)=1 (imat was not defined)
CPVCS
CPVCS       Rev 1.12   Fri Oct 31 10:49:46 1997   dcg
CPVCS    declare ipcmoprm as a pointer
CPVCS
CPVCS       Rev 1.11   Thu May 01 14:14:20 1997   dcg
CPVCS    read polygon data correctly (4 bytes at a time)
CPVCS
CPVCS       Rev 1.10   Mon Apr 14 16:57:48 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.9   Sun Feb 23 10:33:24 1997   het
CPVCS    Change to read binary AMR grids.
CPVCS
CPVCS       Rev 1.8   Mon Dec 02 08:54:18 1996   het
CPVCS    Read element size variable fields.
CPVCS
CPVCS       Rev 1.7   Wed Nov 13 13:48:48 1996   het
CPVCS    read connectivity for quad and tri GMV cells.
c
c ######################################################################
c
      implicit none
C
C args
      character*(*) ifile
      integer ierror
C
      integer charlength
      parameter (charlength=32)
C
      include "local_element.h"
      include 'geom_lg.h'
C
      pointer (ipimt1, imt1)
      integer imt1(*)
C
      pointer (ipitetclr, itetclr(*))
      pointer (ipitettyp, itettyp(*))
      pointer (ipitetoff, itetoff(*))
      pointer (ipjtetoff, jtetoff(*))
      pointer (ipitet, itet1(*))
      pointer (ipjtet, jtet1(*))
      integer itetclr,itettyp,itetoff,jtetoff,itet1,jtet1

C

      pointer (ipxic, xic)
      pointer (ipyic, yic)
      pointer (ipzic, zic)
      real*8 xic(*), yic(*), zic(*)

      pointer (ipout,out)
      real*8 out(*)
C
      real*4 xtime4
      pointer (ipxtemp4, xtemp4)
      real*4 xtemp4(*)
      pointer (ipitemp4, itemp4)
      integer*4 itemp4(*)

      integer*4 liste(100)

      integer*4 iunit4
      integer*4 nelements4,nnodes4,iflag4,maxclrpoint4,ihcycle4
C
      pointer (ipiarray, iarray)
      integer iarray(*)
      pointer (ipxarray, xarray)
      real*8 xarray(*)

      real*8 a,b,c,d,e,f,crosx,crosy,crosz,time,rout
      pointer (ipvels,vels)
      real*8 vels(3,*)

C     for defining sizes in gmv binary files
      integer  IEEE, ASCII, IEEEI4R4, IEEEI4R8, IEEEI8R4,
     * IEEEI8R8, IECXI4R4, IECXI4R8, IECXI8R4, IECXI8R8

      integer ihcycle, ivoronoi2d,ivoronoi3d,ipolydata,
     *  iflag_all,it,ird,lenc,index,lenout,icharlnb,
     *  ics,iout,iflag,itype,ifelmnew,nee,nsdtopo,nsdgeom,
     *  i,i1,i2,i3,i4,i5,i6,i7,i8,icount,ifound,
     *  itoff,jtoff,nelements,ftype,
     *  length,nnodes,iadr,ibytes,icharlnf,len1,maxclrpoint,
     *  iunit,ihybrid,ilen,ityp,ierr,ielements,
     *  j,leni,itp,lout,ierror_return,icscode,ierrwrt,
     *  ier
      integer nsd, nen, nef

      integer nbyte_total, idebug

      logical isgeom
c
      character*32 isubname,cmo,cmoattnam,ctype,cvelnm,
     *         clen,cinter,cpers,cio,crank
      character*32  cmotype,cout,geom_name,c32blank
      character*32  iword1, iword2
      character*132 iline
      character*132 file, Fifile
      character*132 logmess
      character*8092 cbuff
c
      data ivoronoi2d / 1 /
      data ivoronoi3d / 0 /
      data ipolydata / 1 /
      data iflag_all / 0 /

      data IEEE / 0 /
      data ASCII / 1 /
      data IEEEI4R4 / 0 /
      data IEEEI4R8 / 2 /
      data IEEEI8R4 / 3 /
      data IEEEI8R8 / 4 /
      data IECXI4R4 / 5 /
      data IECXI4R8 / 6 /
      data IECXI8R4 / 7 /
      data IECXI8R8 / 8 /

c
      crosx(a,b,c,d,e,f)=b*f-c*e
      crosy(a,b,c,d,e,f)=c*d-a*f
      crosz(a,b,c,d,e,f)=a*e-b*d
c
C###############################################################
C begin

      isubname="readgmv_binary"
c     initialize strings going into cread
      c32blank="                                "
      cmoattnam=c32blank
      ctype=c32blank
      nbyte_total = 0
      idebug = 0
C
      len1=icharlnf(ifile)
      Fifile=ifile(1:len1) // 'F'
      file=ifile(1:len1) // char(0)
      iunit=-1
      call hassign(iunit,Fifile,ierror)
      if (iunit.lt.0 .or. ierror.lt.0) then
        call x3d_error(isubname,'hassign bad file unit')
        write(logmess,*) 'WARNING: can not get file no: '
     *   // Fifile
        call writloga('default',0,logmess,0,ierror)
        goto 9998
      endif
      iunit4=iunit
      call cassignr(iunit4,file,ierror)
C
C tam
C     Removed calls to termcode since we really
C     do not want to stop the program
C     we just want to exit the routine trying to read gmv
C     also moved error reporting from cread to end with rest
C     instead of using a goto statement back up to 98 here.

C     Look for first header word 'gmvinput'
      iadr=0
      ibytes=8

      nbyte_total = nbyte_total + ibytes
      call cread(iunit4,iword1,ibytes,ierror)
      if (ierror .ne.0) go to 9998

      if(iword1(1:8).eq.'gmvinput') then

C       Check second header word for gmv format tag
C       These tags are found in various gmv formats
C       and can be used to recognize and parse new gmv formats
C       as described in sample I/O files found on GMV web site
C       Original lagrit type is IEEE / 0 / same as IEEEI4R4 / 0 /

        ibytes=8
        nbyte_total = nbyte_total + ibytes
        call cread(iunit,iword2,ibytes,ierror)
        if(ierror.ne.0) go to 9998

        ftype = -1
        if (iword2(1:5).eq. 'ascii') ftype=ASCII
        if (iword2(1:6).eq. ' ascii') ftype=ASCII
        if (iword2(1:7).eq. '  ascii') ftype=ASCII
        if (iword2(1:8).eq. '   ascii') ftype=ASCII
        if (iword2(1:4).eq. 'ieee') ftype=IEEEI4R4
        if (iword2(1:5).eq.' ieee') ftype=IEEEI4R4
        if (iword2(1:8).eq.'ieeei4r4') ftype=IEEEI4R4
        if (iword2(1:9).eq.' ieeei4r4') ftype=IEEEI4R4
        if (iword2(1:8).eq.'ieeei4r8') ftype=IEEEI4R8
        if (iword2(1:9).eq.' ieeei4r8') ftype=IEEEI4R8
        if (iword2(1:8).eq.'ieeei8r4') ftype=IEEEI8R4
        if (iword2(1:9).eq.' ieeei8r4') ftype=IEEEI8R4
        if (iword2(1:8).eq.'ieeei8r8') ftype=IEEEI8R8
        if (iword2(1:9).eq.' ieeei8r8') ftype=IEEEI8R8
        if (iword2(1:8).eq.'iecxi4r4') ftype=IECXI4R4
        if (iword2(1:9).eq.' iecxi4r4') ftype=IECXI4R4
        if (iword2(1:8).eq.'iecxi4r8') ftype=IECXI4R8
        if (iword2(1:9).eq.' iecxi4r8') ftype=IECXI4R8
        if (iword2(1:8).eq.'iecxi8r4') ftype=IECXI8R4
        if (iword2(1:9).eq.' iecxi8r4') ftype=IECXI8R4
        if (iword2(1:8).eq.'iecxi8r8') ftype=IECXI8R8
        if (iword2(1:9).eq.' iecxi8r8') ftype=IECXI8R8

C       FORMAT TYPE NOT FOUND
        if (ftype .eq. -1) then
          write(logmess,'(a,a)')
     *    "WARNING: Invalid GMV file type: ", iword2(1:8)
          call writloga('default',0,logmess,0,ierror)
          write(logmess,'(a,a)')"Must be of file type: ",
     *    " ascii, ieee, or iecx "
          call writloga('default',0,logmess,0,ierror)

          call cclose(iunit4)
          goto 9999

C       FORMAT TYPE NOT SUPPORTED
        elseif (ftype .gt. IEEEI4R4 ) then
          write(logmess,'(a,a)')"Unsupported GMV file type: ",
     *    iword2(1:8)
          call writloga('default',0,logmess,0,ierror)
          call cclose(iunit4)
          goto 9999

C       FORMAT TYPE RECOGNIZED
        else
          write(logmess,'(a,a,a)') "Reading GMV binary file "
     >    ,iword1(1:8),iword2(1:8)
          call writloga('default',0,logmess,0,ierror)
        endif

C     IF word 1 is not gmvinput
      else
         write(logmess,'(a)')
     *   "Not a GMV file, gmvinput token not found."
         call writloga('default',0,logmess,0,ierror)
         call cclose(iunit4)
         goto 9999
      endif
C
C     done checking header for tag gmvinput and format ieee
C     continue reading file
      ihybrid=0

C     foreach word1 keyword, read keyword, then fill cmo
C     return here to read each word flag

 100  continue

      ibytes=8
      nbyte_total = nbyte_total + ibytes
      call cread(iunit4,iword1,ibytes,ierror)
      if(ierror.ne.0) go to 9998

      if(iword1(1:6).eq.'endgmv') goto 9999
C
      if(iword1(1:5).eq.'nodes') then
C
         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,nnodes4,ibytes,ierror)
         if(ierror.ne.0) go to 9998
         nnodes = nnodes4
C
         call cmo_get_name(cmo,icscode)
         call cmo_set_info('nnodes',cmo,nnodes,1,1,icscode)
C
         if(nnodes.gt.0) then
C
            call cmo_newlen(cmo,icscode)
C
            call cmo_get_info('xic',cmo,ipxic,ilen,ityp,ierr)
            call cmo_get_info('yic',cmo,ipyic,ilen,ityp,ierr)
            call cmo_get_info('zic',cmo,ipzic,ilen,ityp,ierr)
C
            length=nnodes
            call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,icscode)
C
            ibytes=4*nnodes
            nbyte_total = nbyte_total + ibytes
            call cread(iunit4,xtemp4,ibytes,ierror)
            if(ierror.ne.0) go to 9998
            do i=1,nnodes
               xic(i)=xtemp4(i)
            enddo

            ibytes=4*nnodes
            nbyte_total = nbyte_total + ibytes
            call cread(iunit4,xtemp4,ibytes,ierror)
            if (ierror.ne.0 ) go to 9998
            do i=1,nnodes
               yic(i)=xtemp4(i)
            enddo

            ibytes=4*nnodes
            nbyte_total = nbyte_total + ibytes
            call cread(iunit4,xtemp4,ibytes,ierror)
            if (ierror.ne.0 ) go to 9998
            do i=1,nnodes
               zic(i)=xtemp4(i)
            enddo
C
            call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)
C
         endif
         goto 100

      elseif(iword1(1:5).eq.'cells') then
C
         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,nelements4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         nelements=nelements4
C
         call cmo_get_name(cmo,icscode)
         call cmo_set_info('nnodes',cmo,nnodes,1,1,icscode)
         call cmo_set_info('nelements',cmo,nelements,1,1,icscode)
C
         if(nelements.gt.0) then
            itoff=0
            jtoff=0
            do ielements=1,nelements

               ibytes=8
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,ctype,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               call nulltoblank_lg(ctype,charlength)
               ifound=0

               if(ctype(1:3).eq.'lin') then

                  ibytes=4*(1+nelmnen(ifelmlin))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i1=liste(2)
                  i2=liste(3)
                  do i=2,3
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,4i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,3)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
               elseif(ctype(1:3).eq.'tri') then

                  ibytes=4*(1+nelmnen(ifelmtri))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i1=liste(2)
                  i2=liste(3)
                  i3=liste(4)
                  do i=2,4
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,4i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,4)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
               elseif(ctype(1:4).eq.'quad') then

                  ibytes=4*(1+nelmnen(ifelmqud))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i1=liste(2)
                  i2=liste(3)
                  i3=liste(4)
                  i4=liste(5)
                  do i=2,5
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,4i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,5)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
               elseif(ctype(1:3).eq.'tet') then

                  ibytes=4*(1+nelmnen(ifelmtet))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i1=liste(2)
                  i2=liste(3)
                  i3=liste(4)
                  i4=liste(5)
                  do i=2,5
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,4i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,5)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
                  if(i1.eq.i4) ctype='tri'
               elseif(ctype(1:3).eq.'pyr') then

                  ibytes=4*(1+nelmnen(ifelmpyr))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i5=liste(2)
                  i1=liste(3)
                  i4=liste(4)
                  i3=liste(5)
                  i2=liste(6)
                  do i=2,6
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,9i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,6)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
               elseif(ctype(1:3).eq.'pri') then

                  ibytes=4*(1+nelmnen(ifelmpri))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i1=liste(2)
                  i2=liste(3)
                  i3=liste(4)
                  i4=liste(5)
                  i5=liste(6)
                  i6=liste(7)
                  do i=2,7
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,9i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,7)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
               elseif(ctype(1:3).eq.'hex') then

                  ibytes=4*(1+nelmnen(ifelmhex))
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,liste,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998

                  icount=liste(1)
                  i1=liste(2)
                  i2=liste(3)
                  i3=liste(4)
                  i4=liste(5)
                  i5=liste(6)
                  i6=liste(7)
                  i7=liste(8)
                  i8=liste(9)
                  do i=2,9
                     if(liste(i).lt.1.or.liste(i).gt.nnodes) then
                        write(logmess,'(a,9i10)') 'Illegal NN: ',
     *                                            ielements,
     *                                            (liste(j),j=2,9)
                        call writloga('default',0,logmess,0,ierror)
                     endif
                  enddo
                  if(i1.eq.i5.and.i2.eq.i6 .and.
     *                  i3.eq.i7.and.i4.eq.i8) ctype='quad'
               else
                  write(logmess,'(a,a)')
     *               'Unsupported GMV element type: ',ctype
                  call writloga('default',0,logmess,0,ierrwrt)
               endif
               if(ifound.eq.0.and.ielements.eq.1) then
                  if(ctype(1:3).eq.'lin') then
                        cmotype='lin'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=2
                        nen=nelmnen(ifelmlin)
                        nef=nelmnef(ifelmlin)
                        nee=nelmnee(ifelmlin)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  elseif(ctype(1:3).eq.'tri') then
                        cmotype='tri'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=2
                        nen=nelmnen(ifelmtri)
                        nef=nelmnef(ifelmtri)
                        nee=nelmnee(ifelmtri)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  elseif(ctype(1:4).eq.'quad') then
                        cmotype='quad'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=2
                        nen=nelmnen(ifelmqud)
                        nef=nelmnef(ifelmqud)
                        nee=nelmnee(ifelmqud)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  elseif(ctype(1:3).eq.'tet') then
                        cmotype='tet'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=3
                        nen=nelmnen(ifelmtet)
                        nef=nelmnef(ifelmtet)
                        nee=nelmnee(ifelmtet)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  elseif(ctype(1:3).eq.'pyr') then
                        cmotype='pyr'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=3
                        nen=nelmnen(ifelmpyr)
                        nef=nelmnef(ifelmpyr)
                        nee=nelmnee(ifelmpyr)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  elseif(ctype(1:3).eq.'pri') then
                        cmotype='pri'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=3
                        nen=nelmnen(ifelmpri)
                        nef=nelmnef(ifelmpri)
                        nee=nelmnee(ifelmpri)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  elseif(ctype(1:3).eq.'hex') then
                        cmotype='hex'
                        nsd=3
                        nsdgeom=3
                        nsdtopo=3
                        nen=nelmnen(ifelmhex)
                        nef=nelmnef(ifelmhex)
                        nee=nelmnee(ifelmhex)
                        call cmo_set_info('nnodes',cmo,
     *                                    nnodes,1,1,ierror)
                        call cmo_set_info('nelements',cmo,
     *                                    nelements,1,1,ierror)
                        call cmo_set_info('ndimensions_geom',cmo,
     *                                    nsdgeom,1,1,ierror)
                        call cmo_set_info('ndimensions_topo',cmo,
     *                                    nsdtopo,1,1,ierror)
                        call cmo_set_info('nodes_per_element',cmo,
     *                                    nen,1,1,ierror)
                        call cmo_set_info('faces_per_element',cmo,
     *                                    nef,1,1,ierror)
                        call cmo_set_info('edges_per_element',cmo,
     *                                    nee,1,1,ierror)
                        call cmo_newlen(cmo,ierror)
                        call cmo_get_info('itetclr',cmo,
     *                                    ipitetclr,leni,itp,ier)
                        call cmo_get_info('itettyp',cmo,
     *                                    ipitettyp,leni,itp,ier)
                        call cmo_get_info('itetoff',cmo,
     *                                    ipitetoff,leni,itp,ier)
                        call cmo_get_info('jtetoff',cmo,
     *                                    ipjtetoff,leni,itp,ier)
                        call cmo_get_info('itet',cmo,
     *                                    ipitet,leni,itp,ier)
                        call cmo_get_info('jtet',cmo,
     *                                    ipjtet,leni,itp,ier)
                  endif
               endif
               if(ctype(1:3).eq.'tri') then
                  ifelmnew=ifelmtri
               elseif(ctype(1:4).eq.'line') then
                  ifelmnew=ifelmlin
               elseif(ctype(1:4).eq.'quad') then
                  ifelmnew=ifelmqud
               elseif(ctype(1:3).eq.'tet') then
                  ifelmnew=ifelmtet
               elseif(ctype(1:3).eq.'pyr') then
                  ifelmnew=ifelmpyr
               elseif(ctype(1:3).eq.'pri') then
                  ifelmnew=ifelmpri
               elseif(ctype(1:3).eq.'hex') then
                  ifelmnew=ifelmhex
               endif
               if(ihybrid.eq.0 .and.
     *            (nelmnen(ifelmnew).ne.nen .or.
     *             nelmnef(ifelmnew).ne.nef)) then
                  ihybrid=1
                  nen=nelmnen(ifelmhyb)
                  nef=nelmnef(ifelmhyb)
                  nee=nelmnee(ifelmhyb)
                  call cmo_set_info('nodes_per_element',cmo,
     *                              nen,1,1,ierror)
                  call cmo_set_info('faces_per_element',cmo,
     *                              nef,1,1,ierror)
                  call cmo_set_info('edges_per_element',cmo,
     *                              nee,1,1,ierror)
                  call cmo_newlen(cmo,ierror)
                  call cmo_get_info('itetclr',cmo,
     *                              ipitetclr,leni,itp,ier)
                  call cmo_get_info('itettyp',cmo,
     *                              ipitettyp,leni,itp,ier)
                  call cmo_get_info('itetoff',cmo,
     *                              ipitetoff,leni,itp,ier)
                  call cmo_get_info('jtetoff',cmo,
     *                              ipjtetoff,leni,itp,ier)
                  call cmo_get_info('itet',cmo,
     *                              ipitet,leni,itp,ier)
                  call cmo_get_info('jtet',cmo,
     *                              ipjtet,leni,itp,ier)
               endif
               if(ctype(1:3).eq.'lin') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmlin
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  itoff=itoff+nelmnen(ifelmlin)
                  jtoff=jtoff+nelmnef(ifelmlin)
               elseif(ctype(1:3).eq.'tri') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmtri
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  itet1(itoff+3)=i3
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  jtet1(jtoff+3)=0
                  itoff=itoff+nelmnen(ifelmtri)
                  jtoff=jtoff+nelmnef(ifelmtri)
               elseif(ctype(1:4).eq.'quad') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmqud
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  itet1(itoff+3)=i3
                  itet1(itoff+4)=i4
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  jtet1(jtoff+3)=0
                  jtet1(jtoff+4)=0
                  itoff=itoff+nelmnen(ifelmqud)
                  jtoff=jtoff+nelmnef(ifelmqud)
               elseif(ctype(1:3).eq.'tet') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmtet
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  itet1(itoff+3)=i3
                  itet1(itoff+4)=i4
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  jtet1(jtoff+3)=0
                  jtet1(jtoff+4)=0
                  itoff=itoff+nelmnen(ifelmtet)
                  jtoff=jtoff+nelmnef(ifelmtet)
               elseif(ctype(1:3).eq.'pyr') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmpyr
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  itet1(itoff+3)=i3
                  itet1(itoff+4)=i4
                  itet1(itoff+5)=i5
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  jtet1(jtoff+3)=0
                  jtet1(jtoff+4)=0
                  jtet1(jtoff+5)=0
                  itoff=itoff+nelmnen(ifelmpyr)
                  jtoff=jtoff+nelmnef(ifelmpyr)
               elseif(ctype(1:3).eq.'pri') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmpri
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  itet1(itoff+3)=i3
                  itet1(itoff+4)=i4
                  itet1(itoff+5)=i5
                  itet1(itoff+6)=i6
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  jtet1(jtoff+3)=0
                  jtet1(jtoff+4)=0
                  jtet1(jtoff+5)=0
                  itoff=itoff+nelmnen(ifelmpri)
                  jtoff=jtoff+nelmnef(ifelmpri)
               elseif(ctype(1:3).eq.'hex') then
                  itetclr(ielements)=1
                  itettyp(ielements)=ifelmhex
                  itetoff(ielements)=itoff
                  jtetoff(ielements)=jtoff
                  itet1(itoff+1)=i1
                  itet1(itoff+2)=i2
                  itet1(itoff+3)=i3
                  itet1(itoff+4)=i4
                  itet1(itoff+5)=i5
                  itet1(itoff+6)=i6
                  itet1(itoff+7)=i7
                  itet1(itoff+8)=i8
                  jtet1(jtoff+1)=0
                  jtet1(jtoff+2)=0
                  jtet1(jtoff+3)=0
                  jtet1(jtoff+4)=0
                  jtet1(jtoff+5)=0
                  jtet1(jtoff+6)=0
                  itoff=itoff+nelmnen(ifelmhex)
                  jtoff=jtoff+nelmnef(ifelmhex)
               else

               write(logmess,'(a,i10,i5)') "Bad element type: ",
     *         ielements,itype
               call writloga('default',0,logmess,0,ierror)

               endif
            enddo
CCC         call geniee_cmo(cmo)
CCCC        ! moved to end
         endif
         goto 100
C
      elseif(iword1(1:8).eq.'velocity') then

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,iflag4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         iflag=iflag4

C        if iflag is 1
C        create and assign vector attribute vels
         if(iflag.eq.1) then
            call cmo_get_name(cmo,icscode)
            call cmo_get_info('nnodes',cmo,nnodes,ilen,itype,icscode)
            if(nnodes.gt.0) then
               call cmo_get_attinfo('velname',cmo,iout,rout,cvelnm,
     *                        ipout,lout,itype,ierror_return)

               if(ier.ne.0) cvelnm='vels'
               call cmo_get_info(cvelnm,cmo,ipvels,ilen,ityp,ierr)
               if(ierr.ne.0) then
                 cbuff ='cmo/addatt/-def-/vels/VDOUBLE/vector/' //
     *                  'nnodes/linear/permanent/gxa/0.0 ;  ' //
     *                  'finish'
                 call dotask(cbuff,ier)
                 call cmo_get_info('vels',cmo,ipvels,ilen,ityp,ierr)
               endif
               length=nnodes
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nnodes
                  vels(1,i)=xtemp4(i)
               enddo

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nnodes
                  vels(2,i)=xtemp4(i)
               enddo

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nnodes
                  vels(3,i)=xtemp4(i)
               enddo

               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

            endif

C        else iflag is not 1 (should be 0)
C        data were written fgmvwritevelocitydata(igmvtype,u,v,w)
C        note that values are read, but are not copied into cmo
         else
            call cmo_get_name(cmo,icscode)
            call cmo_get_info('nelements',cmo,
     *                        nelements,ilen,itype,icscode)

            if(nelements.gt.0) then
               length=nelements
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998

               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

            endif
         endif
         goto 100
C
C     read and fill each of the attribute arrays
C     iflag 1 is node and iflag 0 is element
C     keyword endvars will end the list
      elseif(iword1(1:8).eq.'variable') then
C
 110     continue

         ibytes=8
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,cmoattnam,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         call nulltoblank_lg(cmoattnam,charlength)
         if(cmoattnam(1:7).eq.'endvars') goto 100

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,iflag4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         iflag=iflag4

C        node attributes
         if(iflag.eq.1) then
            call cmo_get_info('nnodes',cmo,
     *                        nnodes,ilen,itype,icscode)
            len1=icharlnb(cmoattnam)
            do j=1,len1
               if(cmoattnam(j:j).eq.' ') then
                  cmoattnam(j:j)='_'
               endif
            enddo
            call mmfindbk(cmoattnam,cmo,ipout,lenout,icscode)
            if(icscode.eq.0) then
               call cmo_get_attparam(cmoattnam,cmo,index,ctype,crank,
     *          clen,cinter,cpers,cio,ierror_return)

               lenc=icharlnf(ctype)
               length=nnodes
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998

               if(ctype(1:lenc).eq.'VINT') then
                  call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
                  do i=1,nnodes
                     iarray(i)=xtemp4(i)
                  enddo
               elseif(ctype(1:lenc).eq.'VDOUBLE') then
                  call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
                  do i=1,nnodes
                     xarray(i)=xtemp4(i)
                  enddo
               endif
               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

            else
               cbuff='cmo/addatt/' //
     *               cmo(1:icharlnf(cmo)) //
     *               '/' //
     *               cmoattnam(1:icharlnf(cmoattnam)) //
     *               '/VDOUBLE' //
     *               '/scalar/nnodes/linear/permanent/gxa/0.0' //
     *               ' ; finish '
               call dotask(cbuff,ierror)
               call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
               length=nnodes
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nnodes
                  xarray(i)=xtemp4(i)
               enddo

               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)
            endif

C        element attributes
         elseif(iflag.eq.0) then
            call cmo_get_info('nelements',cmo,
     *                        nelements,ilen,itype,icscode)
            len1=icharlnb(cmoattnam)
            do j=1,len1
               if(cmoattnam(j:j).eq.' ') then
                  cmoattnam(j:j)='_'
               endif
            enddo
            call mmfindbk(cmoattnam,cmo,ipout,lenout,icscode)
            if(icscode.eq.0) then
               call cmo_get_attparam(cmoattnam,cmo,index,ctype,crank,
     *         clen,cinter,cpers,cio,ierror_return)

               lenc=icharlnf(ctype)
               length=nelements
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998

               if(ctype(1:lenc).eq.'VINT') then
                  call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
                  do i=1,nelements
                     iarray(i)=xtemp4(i)
                  enddo
               elseif(ctype(1:lenc).eq.'VDOUBLE') then
                  call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
                  do i=1,nelements
                     xarray(i)=xtemp4(i)
                  enddo
               endif
               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

            else
C
C    Special case for AMR variables
C
              if((cmoattnam(1:7).eq.'itetpar') .or.
     *           (cmoattnam(1:7).eq.'itetkid') .or.
     *           (cmoattnam(1:7).eq.'itetlev'))then
               cbuff='cmo/addatt/' //
     *               cmo(1:icharlnf(cmo)) //
     *               '/' //
     *               cmoattnam(1:icharlnf(cmoattnam)) //
     *               '/VINT' //
     *               '/scalar/nelements/linear/permanent/gxa/0.0' //
     *               ' ; finish '
               call dotask(cbuff,ierror)
               call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
               length=nelements
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nelements
                  iarray(i)=nint(xtemp4(i))
               enddo

               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)
C
C    END Special case for AMR variables
C
              else
               cbuff='cmo/addatt/' //
     *               cmo(1:icharlnf(cmo)) //
     *               '/' //
     *               cmoattnam(1:icharlnf(cmoattnam)) //
     *               '/VDOUBLE' //
     *               '/scalar/nelements/linear/permanent/gxa/0.0' //
     *               ' ; finish '
               call dotask(cbuff,ierror)
               call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
               length=nelements
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nelements
                  xarray(i)=xtemp4(i)
               enddo

               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)
              endif
            endif
         else
            call cmo_get_info('nelements',cmo,
     *                        nelements,ilen,itype,icscode)
            len1=icharlnb(cmoattnam)
            do j=1,len1
               if(cmoattnam(j:j).eq.' ') then
                  cmoattnam(j:j)='_'
               endif
            enddo
            call mmfindbk(cmoattnam,cmo,ipout,lenout,icscode)
            if(icscode.eq.0) then

               call cmo_get_attparam(cmoattnam,cmo,index,ctype,crank,
     *          clen,cinter,cpers,cio,ierror_return)

               lenc=icharlnf(ctype)
               length=nelements
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998

               if(ctype(1:lenc).eq.'VINT') then
                  call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
                  do i=1,nelements
                     iarray(i)=xtemp4(i)
                  enddo
               elseif(ctype(1:lenc).eq.'VDOUBLE') then
                  call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
                  do i=1,nelements
                     xarray(i)=xtemp4(i)
                  enddo
               endif
               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)
            else
               cbuff='cmo/addatt/' //
     *               cmo(1:icharlnf(cmo)) //
     *               '/' //
     *               cmoattnam(1:icharlnf(cmoattnam)) //
     *               '/VDOUBLE' //
     *               '/scalar/nelements/linear/permanent/gxa/0.0' //
     *               ' ; finish '
               call dotask(cbuff,ierror)
               call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
               length=nelements
               call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,xtemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do i=1,nelements
                  xarray(i)=xtemp4(i)
               enddo

               call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

            endif
         endif
         goto 110

C
      elseif(iword1(1:5).eq.'flags') then
 120     continue

         ibytes=8
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,cmoattnam,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         call nulltoblank_lg(cmoattnam,charlength)
         if(cmoattnam(1:7).eq.'endflag') goto 100

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,maxclrpoint4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         maxclrpoint=maxclrpoint4

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,iflag4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         iflag=iflag4

         do i=1,maxclrpoint
            iword2=' '
            ibytes=8
            nbyte_total = nbyte_total + ibytes
            call cread(iunit4,iword2,ibytes,ierror)
            if (ierror.ne.0 ) go to 9998
         enddo

         if(iflag.eq.1) then
            call mmfindbk(cmoattnam,cmo,ipout,lenout,icscode)
            if(icscode.eq.0) then
               call cmo_get_attparam(cmoattnam,cmo,index,ctype,crank,
     *         clen,cinter,cpers,cio,ierror_return)

               lenc=icharlnf(ctype)
               if(ctype(1:lenc).eq.'VINT') then
                  call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
                  length=nnodes
                  call mmgetblk('itemp4',isubname,ipitemp4,length,1
     &               ,icscode)

                  ibytes=4*nnodes
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,itemp4,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998
                  do it=1,nnodes
                      iarray(it)=itemp4(it)
                  enddo
                  call mmrelblk('itemp4',isubname,ipitemp4,icscode)

               elseif(ctype(1:lenc).eq.'VDOUBLE') then
                  call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
                  length=nnodes
                  call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

                  ibytes=4*nnodes
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,xtemp4,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998
                  do i=1,nnodes
                     xarray(i)=xtemp4(i)
                  enddo
                  call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

               endif
            else
               cbuff='cmo/addatt/' //
     *               cmo(1:icharlnf(cmo)) //
     *               '/' //
     *               cmoattnam(1:icharlnf(cmoattnam)) //
     *               '/VINT' //
     *               '/scalar/nnodes/linear/permanent/gxa/0.0' //
     *               ' ; finish '
               call dotask(cbuff,ierror)
               call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
               length=nnodes
               call mmgetblk('itemp4',isubname,ipitemp4,length,1
     &            ,icscode)

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,itemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do it=1,nnodes
                  iarray(it)=itemp4(it)
               enddo
               call mmrelblk('itemp4',isubname,ipitemp4,icscode)

            endif
         elseif(iflag.eq.0) then
            call mmfindbk(cmoattnam,cmo,ipout,lenout,icscode)
            if(icscode.eq.0) then
               call cmo_get_attparam(cmoattnam,cmo,index,ctype,crank,
     *         clen,cinter,cpers,cio,ierror_return)
               lenc=icharlnf(ctype)
               if(ctype(1:lenc).eq.'VINT') then
                  call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
                  length=nelements
                  call mmgetblk('itemp4',isubname,ipitemp4,length,1
     &               ,icscode)

                  ibytes=4*nelements
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,itemp4,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998
                  do it=1,nelements
                     iarray(it)=itemp4(it)
                  enddo
                  call mmrelblk('itemp4',isubname,ipitemp4,icscode)

               elseif(ctype(1:lenc).eq.'VDOUBLE') then
                  call mmfindbk(cmoattnam,cmo,ipxarray,lenout,icscode)
                  length=nelements
                  call mmgetblk('xtemp4',isubname,ipxtemp4,length,2,ics)

                  ibytes=4*nelements
                  nbyte_total = nbyte_total + ibytes
                  call cread(iunit4,xtemp4,ibytes,ierror)
                  if (ierror.ne.0 ) go to 9998
                  do i=1,nelements
                     xarray(i)=xtemp4(i)
                  enddo
                  call mmrelblk('xtemp4',isubname,ipxtemp4,icscode)

               endif
            else
               cbuff='cmo/addatt/' //
     *               cmo(1:icharlnf(cmo)) //
     *               '/' //
     *               cmoattnam(1:icharlnf(cmoattnam)) //
     *               '/VINT' //
     *               '/scalar/nelements/linear/permanent/gxa/0.0' //
     *               ' ; finish '
               call dotask(cbuff,ierror)
               call mmfindbk(cmoattnam,cmo,ipiarray,lenout,icscode)
               length=nelements
               call mmgetblk('itemp4',isubname,ipitemp4,length,1
     &            ,icscode)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,itemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do it=1,nelements
                   iarray(it)=itemp4(it)
               enddo
               call mmrelblk('itemp4',isubname,ipitemp4,icscode)

            endif
         endif
C
         goto 120
      elseif(iword1(1:8).eq.'material') then
      call cmo_get_attinfo('geom_name',cmo,iout,rout,geom_name,
     *                        ipout,lout,itype,ierror)
      isgeom=.true.
      call mmfindbk('cmregs',geom_name,ipcmregs,length,ierror)
      if(ierror.ne.0) then
             isgeom=.false.
             write(logmess,'(a)') ' no geometry for mesh object'
             call writloga('default',0,logmess,0,ierrwrt)
      endif
      call mmfindbk('matregs',geom_name,ipmatregs,length,ierror)
      if(ierror.ne.0) isgeom=.false.
C
         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,maxclrpoint4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         maxclrpoint=maxclrpoint4

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,iflag4,ibytes,ierror)
         iflag=iflag4

         if (ierror.ne.0 ) go to 9998
         do i=1,maxclrpoint
            iword2=' '
            ibytes=8
            nbyte_total = nbyte_total + ibytes
            call cread(iunit4,iword2,ibytes,ierror)
            if (ierror.ne.0 ) go to 9998
            call nulltoblank_lg(iword2,charlength)
            len1=icharlnb(iword2)
            do j=1,len1
               if(iword2(j:j).eq.' ') then
                  iword2(j:j)='_'
               endif
            enddo

            if(isgeom) then
              do j=1,nmregs
                if(cmregs(j).eq.iword2(1:len1)) matregs(j)=i
              enddo
            endif
         enddo

         if(iflag.eq.1) then
            call cmo_get_info('nnodes',cmo,nnodes,ilen,itype,icscode)
            if(nnodes.gt.0) then
               call cmo_get_info('nelements',cmo,
     *                           nelements,ilen,itype,icscode)
               call cmo_get_info('imt1',cmo,ipimt1,ilen,ityp,ierr)
               call cmo_get_info('itetclr',cmo,ipitetclr,ilen,ityp,ierr)
               call cmo_get_info('itetoff',cmo,
     *                           ipitetoff,leni,itp,ier)
               call cmo_get_info('itet',cmo,
     *                           ipitet,leni,itp,ier)
               length=nnodes
               call mmgetblk('itemp4',isubname,ipitemp4,length,1
     &            ,icscode)

               ibytes=4*nnodes
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,itemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do it=1,nnodes
                  imt1(it)=itemp4(it)
               enddo
               do it=1,nelements
                  itetclr(it)=imt1(itet1(itetoff(it)+1))
               enddo
               call mmrelblk('itemp4',isubname,ipitemp4,icscode)

            endif
         elseif(iflag.eq.0) then
            call cmo_get_info('nelements',cmo,
     *                        nelements,ilen,itype,icscode)
            call cmo_get_info('itetclr',cmo,ipitetclr,ilen,ityp,ierr)
            if(nelements.gt.0) then
               length=nelements
               call mmgetblk('itemp4',isubname,ipitemp4,length,1
     &            ,icscode)

               ibytes=4*nelements
               nbyte_total = nbyte_total + ibytes
               call cread(iunit4,itemp4,ibytes,ierror)
               if (ierror.ne.0 ) go to 9998
               do it=1,nelements
                  itetclr(it)=itemp4(it)
               enddo
               call mmrelblk('itemp4',isubname,ipitemp4,icscode)

            endif
         endif
         goto 100
C
      elseif(iword1(1:8).eq.'polygons') then
         len1=icharlnf(iword1)
         ird=1
         ibytes=4
         dowhile(ird.eq.1)

             nbyte_total = nbyte_total + ibytes
             call cread(iunit4,iword1,ibytes,ierror)
             if (ierror.ne.0 ) go to 9998

             if(iword1(1:4).eq.'endp') then
                nbyte_total = nbyte_total + ibytes
                call cread(iunit4,iword1,ibytes,ierror)
                if (ierror.ne.0 ) go to 9998
                if(iword1(1:3).eq.'oly') ird=0
             endif

         enddo
         goto 100
C
      elseif(iword1(1:7).eq.'cycleno') then

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,ihcycle4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         ihcycle = ihcycle4

         call set_global('ihcycle',
     *                   ihcycle,rout,cout,1,icscode)
         if (icscode .ne. 0) call x3d_error(isubname,'set_global')
         goto 100

      elseif(iword1(1:8).eq.'probtime') then

         ibytes=4
         nbyte_total = nbyte_total + ibytes
         call cread(iunit4,xtime4,ibytes,ierror)
         if (ierror.ne.0 ) go to 9998
         time=xtime4

         call set_global('time',iout,
     *                   time,cout,2,icscode)
         if (icscode .ne. 0) call x3d_error(isubname,'set_global')
         goto 100

C     this is not always an error - skip and continue
C     often encounter spaces or ?? characters
C     need to check that all tags are parsed correctly
C     and that write gmv is writing binary correctly
C     check variable polygons cycleno probtime
      else

C        write(logmess,9000) iword1(1:icharlnf(iline))
C        call writloga('default',0,logmess,0,ierr)
C9000    format('GMV binary input not recognized:  ',a)
         goto 100
      endif
C

 9998 continue

      if (ierror .ne. 0) then
         write(logmess,98) file(1:icharlnf(file)),nbyte_total
   98    format('ERROR reading GMV binary file: ',a32,
     >   " near byte ",i10)
         call writloga('default',0,logmess,0,ierror)
      endif

 9999 continue

      call cmo_get_info('nnodes',cmo,
     *                   nnodes,ilen,itype,ierror)
      call cmo_get_info('nelements',cmo,
     *                   nelements,ilen,itype,icscode)

      if (ierror.ne. 0) then
         write(logmess,'(a,a)')
     *   "Error: Problem setting up cmo ",cmo
         call writloga('default',1,logmess,0,ierr)

      else if (nelements.gt.0) then
          call dotask('geniee; finish',ierror)
      endif

      write(logmess,90) nnodes,nelements
   90 format('Nodes: ',i10 ' Elements: ',i10)
      call writloga('default',0,logmess,0,ierror)

      write(logmess,91) nbyte_total,file(1:icharlnf(file))
   91 format('Done reading GMV binary total bytes: ',i10,
     >" file:  ",a)
      call writloga('default',0,logmess,1,ierror)

      return
      end

ext: f


subroutine qinit(maxmx, maxmy, meqn, mbc, mx, my, x_lower, y_lower, dx, dy, q)

! Sets initial condition.

    !$ use omp_lib
    implicit none

    integer, intent(in) :: maxmx, maxmy, meqn, mbc, mx, my
    double precision, intent(in) :: x_lower, y_lower, dx, dy
    double precision, intent(out) :: q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)

    integer :: ix, iy
    double precision :: x, y
    double precision, external :: true_solution

    !$omp parallel do private(x, y, ix)
    do iy = 1, my
        y = y_lower + (iy-0.5d0)*dy
        do ix = 1, mx
            x = x_lower + (ix-0.5d0)*dx
            call set_true_solution(x, y, 0.d0, q(ix, iy, :))
        end do
    end do

end subroutine qinit


! { dg-do run }
      SUBROUTINE WORK1()
      END SUBROUTINE WORK1
      SUBROUTINE WORK2()
      END SUBROUTINE WORK2
      PROGRAM A10
!$OMP PARALLEL
!$OMP SINGLE
        print *, "Beginning work1."
!$OMP END SINGLE
        CALL WORK1()
!$OMP SINGLE
        print *, "Finishing work1."
!$OMP END SINGLE
!$OMP SINGLE
        print *, "Finished work1 and beginning work2."
!$OMP END SINGLE NOWAIT
        CALL WORK2()
!$OMP END PARALLEL
      END PROGRAM A10


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     #############################################################
c     ##                                                         ##
c     ##  subroutine makexyz  --  convert internal to Cartesian  ##
c     ##                                                         ##
c     #############################################################
c
c
c     "makexyz" generates a complete set of Cartesian coordinates
c     for a full structure from the internal coordinate values
c
c
      subroutine makexyz
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'zcoord.i'
      integer i,chiral
      integer ia,ib,ic
      real*8 bond
      real*8 angle1
      real*8 angle2
c
c
c     loop over each atom in turn, finding its coordinates
c
      do i = 1, n
         ia = iz(1,i)
         ib = iz(2,i)
         ic = iz(3,i)
         chiral = iz(4,i)
         bond = zbond(i)
         angle1 = zang(i)
         angle2 = ztors(i)
         call xyzatm (i,ia,bond,ib,angle1,ic,angle2,chiral)
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
! posterior samples code

module posterior
  use utils1
  use xmeans_clstr

  implicit none

  double precision, dimension(:,:,:), allocatable :: evdatp
  integer, dimension(:), allocatable :: nbranchp,nPtPerNode,ncon,nSamp,clstrdNode
  double precision, dimension(:,:), allocatable :: branchp
  double precision, dimension(:,:), allocatable :: pts,pts2,consP,unconsP,pwt,pNwt
  logical, dimension(:), allocatable :: check
  integer nClst,nUncon
  double precision, dimension(:,:), allocatable :: stMu,stSigma

contains

!------------------------------------------------------------------------

 subroutine pos_samp(Ztol,nIter,root,nLpt,ndim,nCdim,nPar,multimodal,outfile,globZ,globinfo,ic_n,ic_Z,ic_info,ic_reme, &
 ic_vnow,ic_npt,ic_nBrnch,ic_brnch,phyP,l,evDataAll,IS,IS_Z,dumper,context)
 !subroutine pos_samp
  	implicit none

	double precision Ztol !null evidence
	integer nIter !globff (total no. replacements)
	character(len=root_max_len)root !base root
	integer nLpt !no. of live points
	integer ndim !dimensionality
	integer nCdim !no. of parameters to cluster on
	integer nPar !total no. of parameters to save
	logical multimodal
	logical outfile !write output files?
	integer i,j,k,i1,ios,ic_n,m,indx
	integer ic_npt(ic_n),ic_nptloc(ic_n),ic_nBrnch(ic_n)
	double precision ic_Z(ic_n),ic_info(ic_n),ic_vnow(ic_n),ic_brnch(ic_n,ic_n),phyP(nPar,nLpt),l(nLpt),evDataAll(:)
	logical IS !importance sampling?
	double precision IS_Z(2) !importance sampling log(Z) estimate & its standard deviation
	logical ic_reme(ic_n)
  	character(len=root_max_len) evfile,livefile,postfile,resumefile
      	character(len=root_max_len) sepFile,statsFile,postfile4,strictSepFile,summaryFile
  	character(len=32) fmt,fmt2
      	logical l1
      	double precision d1,d2,urv
      	double precision ltmp(nPar+2),ic_zloc(ic_n),ic_infoloc(ic_n),gzloc,ginfoloc
	integer phyID(nLpt)
	integer context

      	!posterior info
      	double precision lognpost,globZ,globInfo,gZold,maxWt
      	integer npost !no. of equally weighted posterior samples
      	double precision, dimension(:,:), allocatable :: wt,temp !probability weight of each posterior sample
      	double precision, dimension(:), allocatable :: ic_zold,llike !local evidence

	! parameters for dumper
	double precision, pointer :: physLive(:,:), posterior(:,:), paramConstr(:)
	double precision maxLogLike, logZ, INSlogZ, logZerr
	integer nSamples

	INTERFACE
		!the user dumper function
    		subroutine dumper(nSamples, nlive, nPar, physLive, posterior, paramConstr, maxLogLike, logZ, INSlogZ, logZerr, context_pass)
			integer nSamples, nlive, nPar, context_pass
			double precision, pointer :: physLive(:,:), posterior(:,:), paramConstr(:)
			double precision maxLogLike, logZ, INSlogZ, logZerr
		end subroutine dumper
	end INTERFACE

	ic_zloc=ic_z
	ic_infoloc=ic_info
	gzloc=globz
	ginfoloc=globinfo
	ic_nptloc=ic_npt

	!Ztol=-1.d99
      	!file names
      	postfile=TRIM(root)//'.txt'
      	statsFile=TRIM(root)//'stats.dat'
      	postfile4=TRIM(root)//'post_equal_weights.dat'
      	strictSepFile=TRIM(root)//'post_separate_strict.dat'
      	sepFile=TRIM(root)//'post_separate.dat'
      	evfile=TRIM(root)//'ev.dat'
      	livefile = TRIM(root)//'phys_live.points'
      	resumefile = TRIM(root)//'resume.dat'
      	summaryFile = TRIM(root)//'summary.txt'

      	allocate(branchp(0:ic_n,ic_n),evdatp(ic_n,nIter,nPar+2),wt(ic_n,nIter))
      	allocate(nbranchp(0:ic_n),nPtPerNode(ic_n))
      	allocate(pts2(ndim+1,nIter),pts(nPar+2,nIter),consP(nIter,nPar+2),unconsP(nIter,nPar+2),pwt(nIter,ic_n), &
	pNwt(nIter,ic_n))
      	allocate(ncon(ic_n),nSamp(ic_n),clstrdNode(ic_n))
	allocate(check(ic_n),ic_zold(ic_n),temp(ic_n,3))
	allocate(stMu(ic_n,nPar),stSigma(ic_n,nPar),llike(ic_n))

	nPtPerNode=0
	clstrdNode=0
      	nbranchp=0
	check=.false.
	nSamp=0

	nbranchp(1:ic_n)=ic_nbrnch(1:ic_n)
	branchp(1:ic_n,:)=ic_brnch(1:ic_n,:)

	!set the membership of live points to their coresponding nodes
	k=0
	do i=1,ic_n
		if(ic_nptloc(i)==0) cycle
		phyID(k+1:k+ic_nptloc(i))=i
		k=k+ic_nptloc(i)
	enddo

	!add in the contribution of the remaining live points to the evidence
	i1=1
	do i=1,nLpt
		ltmp(1:nPar)=phyP(1:nPar,i)
		ltmp(nPar+1)=l(i)
		i1=phyID(i)

		d1=ltmp(nPar+1)+log(ic_vnow(i1)/dble(ic_nptloc(i1)))

		!global evidence & info
		gZold=gzloc
		gzloc=LogSumExp(gzloc,d1)
!		ginfoloc=exp(d1-gzloc)*ltmp(nPar+1)+exp(gZold-gzloc)*(ginfoloc+gZold)-gzloc
		ginfoloc=ginfoloc*exp(gzold-gzloc)+exp(d1-gzloc)*ltmp(nPar+1)

		!local evidence & info
!		gZold=ic_Z(i1)
		ic_zold(i1)=ic_zloc(i1)
		ic_zloc(i1)=LogSumExp(ic_zloc(i1),d1)
!		ic_infoloc(i1)=exp(d1-ic_zloc(i1))*ltmp(nPar+1)+exp(gZold-ic_zloc(i1))*(ic_infoloc(i1)+gZold)-ic_zloc(i1)
		ic_infoloc(i1)=ic_infoloc(i1)*exp(ic_zold(i1)-ic_zloc(i1))+exp(d1-ic_zloc(i1))*ltmp(nPar+1)
	enddo

	ginfoloc=ginfoloc-gzloc
	ic_infoloc(1:ic_n)=ic_infoloc(1:ic_n)-ic_zloc(1:ic_n)

      	!make the top level branch
      	nbranchp(0)=1
      	branchp(0,1)=1.d0

      	lognpost=0.d0
	if(outfile) then
      		!read the ev.dat file & calculate the probability weights
      		open(unit=55,file=evfile,status='old')
	endif
	i=0
    	do
		if(outfile) then
    			read(55,*,IOSTAT=ios) ltmp(1:nPar+2),i1

			!end of file?
        		if(ios<0) then
				close(55)
				exit
			endif
		else
			i=i+1
			if(i*(nPar+3)>size(evDataAll)) exit
			ltmp(1:nPar+2)=evDataAll((i-1)*(nPar+3)+1:i*(nPar+3)-1)
			i1=int(evDataAll(i*(nPar+3)))
		endif

		nPtPerNode(i1)=nPtPerNode(i1)+1
		evdatp(i1,nPtPerNode(i1),1:nPar+2)=ltmp(1:nPar+2)

		!probability weight
		wt(i1,nPtPerNode(i1))=exp(evdatp(i1,nPtPerNode(i1),nPar+1)+evdatp(i1,nPtPerNode(i1),nPar+2)-gzloc)
		if(wt(i1,nPtPerNode(i1))>0.d0) lognpost=lognpost-wt(i1,nPtPerNode(i1))*log(wt(i1,nPtPerNode(i1)))
    	enddo

	allocate(posterior(nIter, nPar+2), physLive(nLpt, nPar+1), paramConstr(4*nPar))
	paramConstr(1:2*nPar) = 0d0
	maxLogLike = -huge(1d0)
	logZ = gzloc
	nSamples = nIter
	maxWt = 0d0
	m = 0

 	!read in final remaining points & calculate their probability weights
     	do i=1,nLpt
		ltmp(1:nPar)=phyP(1:nPar,i)
		ltmp(nPar+1)=l(i)
		i1=phyID(i)

		physLive(i,1:nPar+1) = ltmp(1:nPar+1)
		if( ltmp(nPar+1) > maxLogLike ) then
			maxLogLike = ltmp(nPar+1)
			indx = i
		endif
		nPtPerNode(i1)=nPtPerNode(i1)+1
            	evdatp(i1,nPtPerNode(i1),1:nPar+1)=ltmp(1:nPar+1)
            	evdatp(i1,nPtPerNode(i1),nPar+2)=log(ic_vnow(i1)/dble(ic_nptloc(i1)))
		wt(i1,nPtPerNode(i1))=exp(evdatp(i1,nPtPerNode(i1),nPar+1)+evdatp(i1,nPtPerNode(i1),nPar+2)-gzloc)
		if(wt(i1,nPtPerNode(i1))>0.d0) lognpost=lognpost-wt(i1,nPtPerNode(i1))*log(wt(i1,nPtPerNode(i1)))
	enddo
	! global maxlike parameters
	paramConstr(nPar*2+1:nPar*3) = physLive(indx,1:nPar)

      	!no. of equally weighted posterior samples
      	npost=nint(exp(lognpost))

      	!write the global posterior files
	if(outfile) then
      		open(55,file=postfile,form='formatted',status='replace')
      		open(56,file=postfile4,form='formatted',status='replace')
      		write(fmt,'(a,i4,a)')  '(',nPar+2,'E28.18)'
      		write(fmt2,'(a,i4,a)')  '(',nPar+1,'E28.18)'
	endif
      	do i=1,ic_n
      		do j=1,nPtPerNode(i)
			m = m + 1
			posterior(m, 1:nPar+1) = evdatp(i, j, 1:nPar+1)
			posterior(m, nPar+2) = wt(i,j)
			if( wt(i,j) > maxWt ) then
				indx = m
				maxWt = wt(i,j)
			endif
			! global paramater means
			paramConstr(1:nPar) = paramConstr(1:nPar) + evdatp(i, j, 1:nPar) * wt(i,j)
			! global paramater standard deviations
			paramConstr(nPar+1:2*nPar) = paramConstr(nPar+1:2*nPar) + (evdatp(i, j, 1:nPar)**2.0) * wt(i,j)

            		if(wt(i,j)>1.d-99) then
            			if(outfile) write(55,fmt) wt(i,j),-2.d0*evdatp(i,j,nPar+1),evdatp(i,j,1:nPar)

      				!find the multiplicity
      				d1=wt(i,j)*npost
            			!calculate the integer part of multiplicity
        	    		k=int(d1)
            			!calculate the remaining part of multiplicity
            			d2=d1-dble(k)
            			!increase the multiplicity by one with probability d2
				urv=ranmarns(0)
        	    		if(urv<=d2) k=k+1
				if(outfile) then
      					do i1=1,k
            					write(56,fmt) evdatp(i,j,1:nPar+1)
            				enddo
				endif
			endif
		enddo
	enddo

	! global paramater standard deviations
	paramConstr(nPar+1:2*nPar) = sqrt( paramConstr(nPar+1:2*nPar) - paramConstr(1:nPar)**2.0 )
	! global MAP parameters
	paramConstr(nPar*3+1:nPar*4) = posterior(indx,1:nPar)

	if(outfile) then
      		close(55)
      		close(56)

		open(unit=57,file=statsFile,form='formatted',status='replace')
      		!stats file
		write(57,'(a,E28.18,a,E28.18)')		"Nested Sampling Global Log-Evidence           :",gzloc,"  +/-",sqrt(ginfoloc/dble(nLpt))
		if( IS ) write(57,'(a,E28.18,a,E28.18)')"Nested Importance Sampling Global Log-Evidence:",IS_Z(1),"  +/-",IS_Z(2)

		!now the separated posterior samples

	      	!generate the point set to be used by the constrained clustering algorithm
	      	nUncon=0
	      	nClst=0
		i=0
		j=1
	      	call genPoints(i,j,nPar,ic_reme)

		do i=1,nClst
			temp(i,1)=ic_zloc(clstrdNode(i))
			temp(i,2)=ic_infoloc(clstrdNode(i))
			temp(i,3)=ic_nptloc(clstrdNode(i))
		enddo
		ic_zloc(1:nClst)=temp(1:nClst,1)
		ic_infoloc(1:nClst)=temp(1:nClst,2)
		ic_nptloc(1:nClst)=temp(1:nClst,3)

		!now arrange the point set, constrained points first, unconstrained later
		pNwt=0.d0
		pwt=0.d0
	      	k=0
	      	do i=1,nClst
			llike(i)=minval(consP(k+1:k+nCon(i),nPar+1))
	      		do j=1,nCon(i)
	      			k=k+1
	          		pts(1:nPar+2,k)=consP(k,1:nPar+2)
				pwt(k,i)=1.d0
			enddo
		enddo
	      	do i=1,nUncon
	      		k=k+1
	      		pts(1:nPar+2,k)=unconsP(i,1:nPar+2)
		enddo

		i1=0
	      	do
			if(.false.) then
				i1=i1+1

				do i=1,k
					pts2(1:ndim,i)=pts(1:ndim,i)
					pts2(ndim+1,i)=pts(nPar+1,i)
				enddo

	      			!perform constrained clustering
				i=nCdim
				l1=.true.

	      			call GaussMixExpMaxLike(i,nClst,k,nCon(1:nClst),.true.,pts2(1:i,1:k), &
	      			pts(nPar+1:nPar+2,1:k),pwt(1:k,1:nClst),pNwt(1:k,1:nClst),ic_zloc(1:nClst), &
				llike(1:nClst),l1)
			endif

			!calculate cluster properties
			do i=1,nClst
				call rGaussProp(k,nCdim,pts(1:nCdim,1:k),pts(nPar+1:nPar+2,1:k), &
				pwt(1:k,i),pNwt(1:k,i),stMu(i,1:nCdim),stSigma(i,1:nCdim),ic_zloc(i))
			enddo

			i=sum(nCon(1:nClst))
			j=nCdim

			!if(.not.merge(nClst,j,nPar,i,k,nCon(1:nClst),pts(:,1:k),stMu(1:nClst,1:j), &
			!stSigma(1:nClst,1:j),ic_zloc(1:nClst),Ztol,pwt(1:k,1:nClst),pNwt(1:k,1:nClst), &
			!llike(1:nClst))) exit
			exit
		enddo

		!open the output file
		if(multimodal .and. outfile) then
			!open(unit=56,file=strictSepFile,form='formatted',status='replace')
	      		open(unit=55,file=sepFile,form='formatted',status='replace')
			write(57,'(a)')
	      		write(57,'(a)')"Local Mode Properties"
	      		write(57,'(a)')"-------------------------------------------"
	      		write(57,'(a)')
			write(57,'(a,i12)')"Total Modes Found:",nClst
		endif

		open(unit=58,file=summaryFile,status='unknown')
		if( IS ) then
			write(fmt,'(a,i4,a)')  '(',nPar*4+4,'E28.18)'
			write(58,fmt)paramConstr(1:nPar),paramConstr(nPar+1:2*nPar),paramConstr(nPar*2+1:nPar*3), &
			paramConstr(nPar*3+1:nPar*4),globZ,maxLogLike,IS_Z
		else
			write(fmt,'(a,i4,a)')  '(',nPar*4+2,'E28.18)'
			write(58,fmt)paramConstr(1:nPar),paramConstr(nPar+1:2*nPar),paramConstr(nPar*2+1:nPar*3), &
			paramConstr(nPar*3+1:nPar*4),globZ,maxLogLike
		endif
		call genSepFiles(k,nPar,nClst,Ztol,pts,pNwt(1:k,1:nClst),nCon(1:nClst),ic_zloc(1:nClst), &
		ic_infoloc(1:nClst), ic_nptloc(1:nClst),55,56,57,58,multimodal)
		close(58)

	      	if(multimodal) close(55)
	      	close(57)
	endif

	!error on global evidence
	logZerr=sqrt(ginfoloc/dble(nLpt))

	! call the dumper
	INSlogZ = logZ
	if( IS ) INSlogZ = IS_Z(1)
	call dumper(nSamples, nLpt, nPar, physLive, posterior, paramConstr, maxLogLike, logZ, INSlogZ, logZerr, context)

      	deallocate(branchp,evdatp,wt)
      	deallocate(nbranchp,nPtPerNode)
      	deallocate(pts2,pts,consP,unconsP,pwt,pNwt)
      	deallocate(ncon,nSamp,clstrdNode)
	deallocate(check,ic_zold,temp)
	deallocate(stMu,stSigma,llike)
	deallocate(posterior, physLive, paramConstr)

  end subroutine pos_samp

!------------------------------------------------------------------------

  recursive subroutine genPoints(br,brNum,nPar,ic_reme)

  	implicit none

      	integer br !branch to be analyzed
      	integer brNum !branching no. of the branch to be analyzed
	integer nPar !dimensionality
	logical ic_reme(:)
      	!work variables
      	integer i,j,k,i1,i2,node

      	node=int(branchp(br,brNum))

      	!find starting node
      	i1=1

      	!find out the ending position
      	i2=nPtPerNode(node)

      	if(nbranchp(node)==0 .and. .not.ic_reme(node) .and. nPtPerNode(node)>0) then
      		!add the points to the constrained point set if encountered the leaf
		!& calculate the means & sigmas
           	nClst=nClst+1
           	nCon(nClst)=i2-i1+1
            	j=sum(nCon(1:nClst-1))
		stMu(nClst,:)=0.d0
		stSigma(nClst,:)=0.d0
		ClstrdNode(nClst)=node
            	do i=i1,i2
            		j=j+1
            		consP(j,1:nPar+2)=evdatp(node,i,1:nPar+2)
			stMu(nClst,1:nPar)=stMu(nClst,1:nPar)+evdatp(node,i,1:nPar)
			stSigma(nClst,1:nPar)=stSigma(nClst,1:nPar)+evdatp(node,i,1:nPar)**2
            	enddo
		stMu(nClst,1:nPar)=stMu(nClst,1:nPar)/dble(nCon(nClst))
		stSigma(nClst,1:nPar)=stSigma(nClst,1:nPar)/dble(nCon(nClst))
		stSigma(nClst,1:nPar)=sqrt(stSigma(nClst,1:nPar)-stMu(nClst,1:nPar)**2)
	elseif(.not.(nbranchp(node)==0)) then
		if(.not.check(node)) then
            		!add the points to the un-constrained point set
            		do i=i1,i2
            			nUncon=nUncon+1
            			unconsP(nUncon,1:nPar+2)=evdatp(node,i,1:nPar+2)
            		enddo
			check(node)=.true.
		endif
            	!now parse the daughter branches
      		do i=1,nbranchp(node)
            		k=node
            		i1=i
			call genPoints(k,i1,nPar,ic_reme)
		enddo
	endif

  end subroutine genPoints

!------------------------------------------------------------------------

  subroutine genSepFiles(npt,nPar,nCls,Ztol,pt,pwt,nCon,locZ,locInfo,locNpt,funit1,funit2,funit3,funit4,multimodal)

  	implicit none

	!input variables
	integer npt !total no. of points
	integer nPar !dimensionality
	integer nCls !no. of modes
	double precision pt(nPar+2,npt) !points
	integer nCon(nCls) !no. of constrained points
	double precision pwt(npt,nCls) !ptrobability weights
	double precision locZ(nCls), locInfo(nCls) !local evidence
	integer locNpt(nCls)
	double precision Ztol
	integer funit1 !file having strictly separated samples
	integer funit2 !file having separated samples
	integer funit3 !stats file
	integer funit4 !summary file
	logical multimodal
	!work variables
	integer i,j,k,indx(1),nliveP
	double precision d1,d2,d3,mean(nCls,nPar),sigma(nCls,nPar),maxLike(nPar),MAP(nPar)
	double precision old_slocZ,sinfo,slocZ
	character*30 fmt,stfmt

	nliveP=sum(locNpt(1:nCls))
	write(stfmt,'(a,i4,a)')  '(',nPar*4+2,'E28.18)'

	!calculate the weights including the posterior component
	do i=1,nCls
		k=sum(nCon(1:i-1))
		!first calculate the evidence & info for points strictly lying the cluster
		slocZ=-huge(1.d0)*epsilon(1.d0) !logZero
		sinfo=0.d0
		do j=k+1,k+nCon(i)
			!local evidence
			old_slocZ=slocZ
			slocZ=logSumExp(slocZ,pt(nPar+1,j)+pt(nPar+2,j))
			!local info
			sinfo=exp(pt(nPar+1,j)+pt(nPar+2,j)-slocZ)*pt(nPar+1,j) &
				+exp(old_slocZ-slocZ)*(sinfo+old_slocZ)-slocZ
		enddo

		if(locZ(i)<Ztol) cycle

		mean(i,:)=0.d0
		sigma(i,:)=0.d0
		nSamp(i)=0

		!normalize
		d1=sum(pwt(1:npt,i))
		pwt(1:npt,i)=pwt(1:npt,i)/d1

		!insert two blank line
!		write(funit2,*)
!            	write(funit2,*)
		if(multimodal) then
			write(funit1,*)
            		write(funit1,*)
		endif
		do j=1,npt
			mean(i,1:nPar)=mean(i,1:nPar)+pt(1:nPar,j)*pwt(j,i)
			sigma(i,1:nPar)=sigma(i,1:nPar)+(pt(1:nPar,j)**2)*pwt(j,i)

			if(multimodal) then
				!write the strictly separate file
      				write(fmt,'(a,i4,a)')  '(',nPar+2,'E28.18)'
				!strictly separate points
!				if(j>k .and. j<k+nCon(i)+1) then
!					!probability weight
!					swt=exp(pt(nPar+1,j)+pt(nPar+2,j)-slocZ)
!					if(swt>1.d-99) then
!						write(funit2,fmt)swt,-2.d0*pt(nPar+1,j),pt(1:nPar,j)
!					endif
!				endif

				!write the separate file
				if(pwt(j,i)>1.d-99) then
					nSamp(i)=nSamp(i)+1
					write(funit1,fmt)pwt(j,i),-2.d0*pt(nPar+1,j),pt(1:nPar,j)
				endif
			endif
		enddo
		sigma(i,1:nPar)=sqrt(sigma(i,1:nPar)-mean(i,1:nPar)**2.)

		!stMu(i,:)=mean(i,:)
		!stSigma(i,:)=sigma(i,:)

		!find maxLike parameters
		k=sum(nCon(1:i-1))
		indx=maxloc(pt(nPar+1,k+1:k+nCon(i)))
		maxLike(1:nPar)=pt(1:nPar,indx(1)+k)
		d2=pt(nPar+1,indx(1)+k)

		!find MAP parameters
		indx=maxloc(pwt(1:npt,i))
		MAP(1:nPar)=pt(1:nPar,indx(1))

		!write the stats file
		if(multimodal) then
			write(funit3,*)
			write(funit3,*)
			write(funit3,'(a,i4)')'Mode',i
			d3=(nliveP-locNpt(i))*sinfo/locInfo(i)+locNpt(i)
			write(funit3,'(a,E28.18,a,E28.18)')"Strictly Local Log-Evidence",slocZ," +/-",sqrt(sinfo/locNpt(i))
			write(funit3,'(a,E28.18,a,E28.18)')"Local Log-Evidence",locZ(i)," +/-",sqrt(locInfo(i)/d3)
     		endif
		write(funit3,'(a)')""
		write(funit3,'(a)')"Dim No.       Mean        Sigma"
           	do j=1,nPar
           		!write(funit3,'(i4,2E28.18)')j,stMu(i,j),stSigma(i,j)
           		write(funit3,'(i4,2E28.18)')j,mean(i,j),sigma(i,j)
           	enddo
            	write(funit3,'(a)')""
            	write(funit3,'(a)')"Maximum Likelihood Parameters"
            	write(funit3,'(a)')"Dim No.        Parameter"
            	do j=1,nPar
           		write(funit3,'(i4,1E28.18)')j,maxLike(j)
           	enddo
      		write(funit3,'(a)')""
           	write(funit3,'(a)')"MAP Parameters"
           	write(funit3,'(a)')"Dim No.        Parameter"
           	do j=1,nPar
           		write(funit3,'(i4,1E28.18)')j,MAP(j)
           	enddo
		write(funit4,stfmt)mean(i,1:nPar),sigma(i,1:nPar),maxLike(1:nPar),MAP(1:nPar),locZ(i),d2
	enddo

  end subroutine genSepFiles

!------------------------------------------------------------------------

 logical function merge(n,nCdim,nPar,npt,gnpt,nptx,pt,mean,sigma,locEv,nullEv,wt,nWt,llike)

 	implicit none

	!input variables
	integer n !no. of clusters
	integer nCdim,nPar !dimensionality
	integer npt !total no. of constrained points
	integer gnpt !total no. of  points
	integer nptx(n) !no. of points in each cluster
	double precision pt(nPar+2,gnpt) !points
	double precision mean(n,nCdim),sigma(n,nCdim)
	double precision locEv(n) !local evidence
	double precision nullEv !null evidence
	double precision wt(gnpt,n)
	double precision nWt(gnpt,n),llike(n)

	!work variables
	integer i,j,k,l,m
	double precision tP(nPar+2,npt),tWt(npt,n)
	logical check(n),flag

	flag = .false.
	merge=.false.
	check=.false.
	do i=1,n
		if(nptx(i)==0 .or. locEv(i)<nullEv) cycle
		do
			do j=1,n
				if(i==j .or. nptx(j)==0 .or. locEv(j)<nullEv) cycle
				flag=.false.

				!merge required?
				l=0
				do k=1,nCdim
					if(abs(mean(i,k)-mean(j,k))<=1.d0*sigma(i,k)) then
						l=l+1
					else
						exit
					endif
				enddo

				!yes, then merge the modes
				if(l==nCdim) then
					flag=.true.
					merge=.true.

					!re-arrange the constrained points & weights
					wt(:,i)=wt(:,i)+wt(:,j)

					!lowest likelihood
					llike(i)=max(llike(i),llike(j))

					l=sum(nptx(1:j-1))
					m=sum(nptx(1:i-1))
					tP(:,1:nptx(j))=pt(:,l+1:l+nptx(j))
					tWt(1:nptx(j),:)=wt(l+1:l+nptx(j),:)
					if(i<j) then
						pt(:,m+nptx(i)+nptx(j)+1:l+nptx(j))=pt(:,m+nptx(i)+1:l)
						pt(:,m+nptx(i)+1:m+nptx(i)+nptx(j))=tP(:,1:nptx(j))
						wt(m+nptx(i)+nptx(j)+1:l+nptx(j),:)=wt(m+nptx(i)+1:l,:)
						wt(m+nptx(i)+1:m+nptx(i)+nptx(j),:)=tWt(1:nptx(j),:)
					else
						pt(:,l+1:m-nptx(j))=pt(:,l+nptx(j)+1:m)
						pt(:,m-nptx(j)+1:m)=tP(:,1:nptx(j))
						wt(l+1:m-nptx(j),:)=wt(l+nptx(j)+1:m,:)
						wt(m-nptx(j)+1:m,:)=tWt(1:nptx(j),:)
					endif
					nptx(i)=nptx(i)+nptx(j)
					nptx(j)=0

					!recalculate the means & sigmas
					call rGaussProp(gnpt,nCdim,pt(1:nCdim,:),pt(nPar+1:nPar+2,:),wt(:,i), &
					nWt(:,i),mean(i,:),sigma(i,:),locEv(i))

					exit
				endif
			enddo

			if(.not.flag) exit
		enddo
	enddo
	j=0
	do i=1,n
		if(nptx(i)>0) then
			j=j+1
			nptx(j)=nptx(i)
			locEv(j)=locEv(i)
			mean(j,:)=mean(i,:)
			sigma(j,:)=sigma(i,:)
			wt(:,j)=wt(:,i)
			nWt(:,j)=nWt(:,i)
			llike(j)=llike(i)
		endif
	enddo
	n=j

 end function merge

!------------------------------------------------------------------------

  subroutine rGaussProp(n,d,p,LX,wt,nWt,mean,sigma,Z)

	implicit none

      	!input variables
      	integer n !no. of points
	integer d !dimensionality
	double precision p(d,n) !points
	double precision LX(2,n) !log-like & log-dX of points
	double precision wt(n)
	double precision nWt(n)

      	!output variables
      	double precision mean(d) !mean
      	double precision sigma(d) !standard deviations
      	double precision Z !local evidence

      	!work variables
      	integer i

	nWt=wt
	!calculate the evidence
!      	Z=-huge(1.d0)*epsilon(1.d0) !logZero
!	do i=1,n
!		if(nWt(i)>0.d0) Z=logSumExp(Z,LX(1,i)+LX(2,i)+log(nWt(i)))
!	enddo

	!now calculate the normalized posterior probabilty weights
	do i=1,n
		if(nWt(i)>0.d0) nWt(i)=nWt(i)*exp(LX(1,i)+LX(2,i)-Z)
	enddo

	!mean & sigma
      	mean=0.d0
	sigma=0.d0
      	do i=1,n
            	mean(1:d)=mean(1:d)+p(1:d,i)*nWt(i)
            	sigma(1:d)=sigma(1:d)+(p(1:d,i)**2)*nWt(i)
	enddo
	sigma(1:d)=sqrt(sigma(1:d)-mean(1:d)**2)

  end subroutine rGaussProp

!----------------------------------------------------------------------

end module posterior

max_stars_repo_name: annis/cosmosis

max_issues_repo_name: annis/cosmosis

max_forks_repo_name: annis/cosmosis


! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

! Test 8.5.18 constraints on the VALUE attribute

module m
  type :: hasCoarray
    real :: coarray[*]
  end type
 contains
  !ERROR: VALUE attribute may apply only to a dummy data object
  subroutine C863(notData,assumedSize,coarray,coarrayComponent)
    external :: notData
    !ERROR: VALUE attribute may apply only to a dummy argument
    real, value :: notADummy
    value :: notData
    !ERROR: VALUE attribute may not apply to an assumed-size array
    real, value :: assumedSize(10,*)
    !ERROR: VALUE attribute may not apply to a coarray
    real, value :: coarray[*]
    !ERROR: VALUE attribute may not apply to a type with a coarray ultimate component
    type(hasCoarray), value :: coarrayComponent
  end subroutine
  subroutine C864(allocatable, inout, out, pointer, volatile)
    !ERROR: VALUE attribute may not apply to an ALLOCATABLE
    real, value, allocatable :: allocatable
    !ERROR: VALUE attribute may not apply to an INTENT(IN OUT) argument
    real, value, intent(in out) :: inout
    !ERROR: VALUE attribute may not apply to an INTENT(OUT) argument
    real, value, intent(out) :: out
    !ERROR: VALUE attribute may not apply to a POINTER
    real, value, pointer :: pointer
    !ERROR: VALUE attribute may not apply to a VOLATILE
    real, value, volatile :: volatile
  end subroutine
  subroutine C865(optional) bind(c)
    !ERROR: VALUE attribute may not apply to an OPTIONAL in a BIND(C) procedure
    real, value, optional :: optional
  end subroutine
end module