Add forder read_examples/ containing bash scripts set to display the output of...

Add forder read_examples/ containing bash scripts set to display the output of the published FLEXPART 10.4 examples/ runs. Add also the folder flex_read_fortran containing the source code for reading the contants of the FLEXPART output files (binary and netcdf)

tests/flex_read_fortran/README

+
+configure the path to netcdf in makefile
+and build
+
+make 
+
+equivalent to 
+
+make printgrid
+make printgrid_netcdf
+make printpart
+

tests/flex_read_fortran/caldate.f

+      SUBROUTINE CALDATE(JULDATE,YYYYMMDD,HHMISS)
+c                           i       o       o
+* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+*                                                                             *
+*     Calculates the Gregorian date from the Julian date                      *
+*                                                                             *
+*     AUTHOR: Andreas Stohl (21 January 1994), adapted from Numerical Recipes *
+*                                                                             *
+*     Variables:                                                              *
+*     DD             Day                                                      *
+*     HH             Hour                                                     *
+*     HHMISS         Hour, Minute, Second                                     *
+*     JA,JB,JC,JD,JE help variables                                           *
+*     JALPHA         help variable                                            *
+*     JULDATE        Julian Date                                              *
+*     JULDAY         help variable                                            *
+*     MI             Minute                                                   *
+*     MM             Month                                                    *
+*     SS             Seconds                                                  *
+*     YYYY           Year                                                     *
+*     YYYYMMDD       Year, Month, Day                                         *
+*                                                                             *
+*     Constants:                                                              *
+*     IGREG          help constant                                            *
+*                                                                             *
+* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+
+      IMPLICIT NONE
+
+      INTEGER YYYYMMDD,YYYY,MM,DD,HHMISS,HH,MI,SS
+      INTEGER JULDAY,JA,JB,JC,JD,JE,IGREG,JALPHA
+      DOUBLE PRECISION JULDATE
+      PARAMETER (IGREG=2299161)
+
+      JULDAY=INT(JULDATE)
+      IF(JULDAY.GE.IGREG)THEN
+        JALPHA=INT(((JULDAY-1867216)-0.25)/36524.25)
+        JA=JULDAY+1+JALPHA-INT(0.25*JALPHA)
+      ELSE
+        JA=JULDAY
+      ENDIF
+      JB=JA+1524
+      JC=INT(6680.+((JB-2439870)-122.1)/365.25)
+      JD=365*JC+INT(0.25*JC)
+      JE=INT((JB-JD)/30.6001)
+      DD=JB-JD-INT(30.6001*JE)
+      MM=JE-1
+      IF (MM.GT.12) MM=MM-12
+      YYYY=JC-4715
+      IF (MM.GT.2) YYYY=YYYY-1
+      IF (YYYY.LE.0) YYYY=YYYY-1
+
+      YYYYMMDD=10000*YYYY+100*MM+DD
+      HH=INT(24.*(JULDATE-FLOAT(JULDAY)))
+      MI=INT(1440.*(JULDATE-FLOAT(JULDAY))-60.*FLOAT(HH))
+      SS=NINT(86400.*(JULDATE-FLOAT(JULDAY))-3600.*FLOAT(HH))
+     +-60.*FLOAT(MI)
+      IF (SS.EQ.60) THEN  ! 60 seconds = 1 minute
+        SS=0
+        MI=MI+1
+      ENDIF
+      IF (MI.EQ.60) THEN
+        MI=0
+        HH=HH+1
+      ENDIF
+      HHMISS=10000*HH+100*MI+SS
+
+      END

tests/flex_read_fortran/flex_read_compare2.f90

+    program flex_read_compare
+
+    parameter(nxmax=720,nymax=360,nzmax=40,maxtime=10)
+    parameter(maxpoint=1,maxspec=1,maxageclass=1)
+
+    double precision :: jul,julstart,juldate
+    real :: outheight(nzmax)
+    integer :: ireleasestart(maxpoint),ireleaseend(maxpoint)
+    integer :: lage(0:maxageclass)
+    real :: xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+    real :: zpoint2(maxpoint)
+    integer :: ireldate(maxspec),ireltime(maxspec)
+
+    real :: grid(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+    real :: grid_reference(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+    real :: grid_diff(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+    real :: footprint(0:nxmax-1,0:nymax-1)
+    real :: footprint_total(0:nxmax-1,0:nymax-1,maxspec)
+    real :: area(0:nxmax-1,0:nymax-1)
+    real :: heightnn(0:nxmax-1,0:nymax-1,0:nzmax)
+    real :: smallnum,r_sparse_dump(5000000)
+    integer :: i_sparse_dump(5000000)
+    parameter(smallnum=1.e-38)
+
+    character dirname*120, dirname_reference*120 
+    logical ::  dirnameisarg
+    character filename_reference*150,filename*150,fileout*150,aday*8,atime*6
+    character(250) :: filegrid, filegrid_reference 
+          
+    logical :: sp_zer
+
+    real :: sizegrid_m
+    real :: sum_grid, sum_sum_grid
+    real :: max_grid, max_max_grid
+    real :: maxdiff, diff_tolerance 
+    character(4) :: grid_magnitude
+
+!    grid_magnitude='time'
+    grid_magnitude='conc'
+!    grid_magnitude='pptv'
+    maxdiff=0.0  
+    
+    dirnameisarg= .FALSE. 
+    ! print*, iargc()
+    select case (iargc())
+    case (3)
+      call getarg(1,dirname)
+      len=index(dirname,' ')-1
+      call getarg(2,dirname_reference)
+      len_reference=index(dirname_reference,' ')-1
+      call getarg(3,grid_magnitude)
+    case (2)
+      !dirnameisarg= .TRUE. 
+      call getarg(1,dirname)
+      len=index(dirname,' ')-1
+      call getarg(2,grid_magnitude)
+    case (1)
+      call getarg(1,grid_magnitude)
+    case (0)
+      grid_magnitude='time'
+    print*,'Error: argument absent'
+    stop 2
+    end select
+    print*, "grid_magnitude=", grid_magnitude
+
+!    if ( .NOT. dirnameisarg) then
+    !******************
+!!        print*, " Read dirlist "
+    !******************
+!!        open(15,file='dirlist')
+!!        200 read(15,'(a)',end=199) dirname_reference
+!!       len_reference=index(dirname_reference,' ')-1
+!    endif
+
+
+     print*, 'dirname input: ', dirname
+     print*, 'dirname reference: ', dirname_reference
+     
+
+!      do 50 mother_or_nest=1,2
+    do 50 mother_or_nest=1,1
+
+    !*********************************
+    ! Read FLEXPART header information
+    !*********************************
+
+        if (mother_or_nest == 1) then
+            filename=dirname_reference(1:len_reference)//'header'
+        else if (mother_or_nest == 2) then
+            filename=dirname_reference(1:len_reference)//'header_nest'
+        endif
+        call readheader(filename,nxmax,numxgrid,nymax,numygrid,nzmax, &
+        numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime, &
+        loutstep,maxspec,nspec,maxageclass,nageclass,lage,ireleasestart, &
+        ireleaseend,maxpoint,numpoint,xpoint,ypoint,zpoint1,zpoint2, &
+        heightnn,area)
+
+        print*, 'ireleasestart, ireleaseend ',ireleasestart, ireleaseend
+
+
+        if (mother_or_nest == 1) then
+            filename=dirname(1:len)//'header'
+        else if (mother_or_nest == 2) then
+            filename=dirname(1:len)//'header_nest'
+        endif
+        call readheader(filename,nxmax,numxgrid,nymax,numygrid,nzmax, &
+        numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime, &
+        loutstep,maxspec,nspec,maxageclass,nageclass,lage,ireleasestart, &
+        ireleaseend,maxpoint,numpoint,xpoint,ypoint,zpoint1,zpoint2, &
+        heightnn,area)
+
+        print*, 'ireleasestart, ireleaseend ',ireleasestart, ireleaseend
+        !STOP
+
+
+        write(*,*) 'nspec',nspec
+        julstart=juldate(ibdate,ibtime)
+        do 30 k=1,nspec
+            jul=julstart+dble(float(ireleaseend(k))/86400.)
+            call caldate(jul,ireldate(k),ireltime(k))
+            print*,'ireldate',ireldate(k)
+            print*,'ireltime',ireltime(k)
+        30 END DO
+
+        max_max_grid=0
+        max_grid=0
+        sum_grid=0
+        sum_sum_grid=0
+        sizegrid_m=0
+
+    ! rint*, 'init sizegrid_m=', sizegrid_m !=0
+
+        do 70 k=1,nspec
+            do 70 ix=0,numxgrid-1
+                do 70 jy=0,numygrid-1
+                    footprint_total(ix,jy,k)=0.
+        70 END DO
+
+    !******************************************
+    ! Loop about all times, given in file dates
+    !******************************************
+        !print*,'> open dates'
+        
+        open(20,file=dirname(1:len)//'dates',form='formatted', &
+        status='old')
+        it=0
+        print*,'> loop on dates'
+        100 read(20,'(i8,i6)',end=99) inday,intime
+        write(aday,'(i8.8)') inday
+        write(atime,'(i6.6)') intime
+        it=it+1
+                 
+        if (mother_or_nest == 1) then
+            !filegrid=dirname(1:len)//'/grid_conc_'//aday//atime//'_001'
+            filegrid=dirname(1:len)//'/grid_'//grid_magnitude//'_' &
+            //aday//atime//'_001'
+
+            filegrid_reference=dirname_reference(1:len_reference)//'/grid_'//grid_magnitude//'_' &
+            //aday//atime//'_001'
+
+            fileout=dirname(1:len)//'footprint_'//aday//atime//'_001'
+        else
+            filegrid=dirname(1:len)//'grid_time_nest_'//aday//atime//'_001'
+            fileout=dirname(1:len)//'footprint_nest_'//aday//atime//'_001'
+        endif
+
+    ! rite(*,*) filegrid
+        print*, 'read filegrid input=', filegrid
+        call readgrid(filegrid_reference,itime,nxmax,numxgrid,nymax,numygrid, &
+        nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid_reference, &
+        lage)
+
+        print*, 'read filegrid reference=', filegrid_reference
+        call readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid, &
+        nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid, &
+        lage)
+
+
+
+
+        grid_diff=grid_reference-grid
+
+        print*, 'sum(reference - input)=', sum(grid_diff) 
+        print*, 'sum(reference - input)=', sum(grid_reference-grid) 
+
+        maxdiff = max(maxdiff,  sum(grid_reference-grid) )
+        print*,'maxdiff =',maxdiff
+        
+
+        !diff_tolerance= max(sum(grid_reference)/sum(grid) , sum(grid)/sum(grid_reference)  )
+        !if  (sum(abs(grid_diff)) > 1 ) then
+        if  (sum(abs(grid_diff)) > sum(abs(grid))/100  ) then
+           print*, 'fail! sum(abs(grid_diff)) > 1 '
+           stop 1
+        endif
+
+
+        
+        !print*, 'stop!'
+        !stop
+
+
+        print*, 'read grid:'
+        print*, 'maxval(grid)=', maxval(grid)
+
+        print*, 'sum(grid)=', sum(grid)
+    !***************************************
+    ! Display read values
+    !***************************************
+    ! izegrid_m=sizegrid_m+sum(grid)
+
+
+        sum_grid=sum(grid)
+        max_grid=maxval(grid)
+
+        print*, 'max_max_grid,max_grid', max_max_grid, max_grid
+        max_max_grid=max(max_max_grid,max_grid)
+        sum_sum_grid=sum_sum_grid+sum_grid
+
+        print*, 'sum_grid=',sum_grid !sum(grid)
+        print*, 'max_grid=', max_grid !  maxval(grid)
+        print*, 'part total', it, aday, atime
+    ! rint*, 'sum sum(grid)=', sum_sum_grid
+    ! rint*, 'max maxval(grid)=', max_max_grid
+        print*, 'sum_sum_grid=', sum_sum_grid
+        print*, 'max_max_grid=', max_max_grid
+    ! rint*, 'sizegrid_m=', sizegrid_m
+
+
+        goto 100 !read dates
+        99 close(20)
+
+
+    50 END DO
+
+!    if ( .NOT. dirnameisarg) then
+!!        goto 200 !read(15, dirlist
+!!        199 close(15) !dirlist
+!    endif
+
+
+    print*, 'for all times'
+    print*, 'sum sum(grid)=', sum_sum_grid
+    print*, 'max maxval(grid)=', max_max_grid
+
+!    print*, 'sum_sum_grid_out=',sum_sum_grid
+!    print*, 'max_max_grid_out=',max_max_grid
+    END PROGRAM

tests/flex_read_fortran/juldate.f

+      FUNCTION JULDATE(YYYYMMDD,HHMISS)
+* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+*                                                                             *
+*     Calculates the Julian date                                              *
+*                                                                             *
+*     AUTHOR: Andreas Stohl (15 October 1993)                                 *
+*                                                                             *
+*     Variables:                                                              *
+*     DD             Day                                                      *
+*     HH             Hour                                                     *
+*     HHMISS         Hour, minute + second                                    *
+*     JA,JM,JY       help variables                                           *
+*     JULDATE        Julian Date                                              *
+*     JULDAY         help variable                                            *
+*     MI             Minute                                                   *
+*     MM             Month                                                    *
+*     SS             Second                                                   *
+*     YYYY           Year                                                     *
+*     YYYYMMDDHH     Date and Time                                            *
+*                                                                             *
+*     Constants:                                                              *
+*     IGREG          help constant                                            *
+*                                                                             *
+* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+
+      IMPLICIT NONE
+
+      INTEGER YYYYMMDD,YYYY,MM,DD,HH,MI,SS,HHMISS
+      INTEGER JULDAY,JY,JM,JA,IGREG
+      DOUBLE PRECISION JULDATE
+      PARAMETER (IGREG=15+31*(10+12*1582))
+
+      YYYY=YYYYMMDD/10000
+      MM=(YYYYMMDD-10000*YYYY)/100
+      DD=YYYYMMDD-10000*YYYY-100*MM
+      HH=HHMISS/10000
+      MI=(HHMISS-10000*HH)/100
+      SS=HHMISS-10000*HH-100*MI
+
+      IF (YYYY.EQ.0) STOP 'There is no Year Zero.'
+      IF (YYYY.LT.0) YYYY=YYYY+1
+      IF (MM.GT.2) THEN
+        JY=YYYY
+        JM=MM+1
+      ELSE
+        JY=YYYY-1
+        JM=MM+13
+      ENDIF
+      JULDAY=INT(365.25*JY)+INT(30.6001*JM)+DD+1720995
+      IF (DD+31*(MM+12*YYYY).GE.IGREG) THEN
+        JA=INT(0.01*JY)
+        JULDAY=JULDAY+2-JA+INT(0.25*JA)
+      ENDIF
+
+      JULDATE=DBLE(FLOAT(JULDAY))+DBLE(FLOAT(HH)/24.)+
+     +DBLE(FLOAT(MI)/1440.)+DBLE(FLOAT(SS)/86400.)
+
+      END

tests/flex_read_fortran/makefile

+SHELL = /bin/bash
+
+FC       = gfortran
+
+INCPATH2 = /opt/tools/netcdf-fortran/4.6.0-gcc9.4/include 
+LIBPATH2 = /opt/tools/netcdf-fortran/4.6.0-gcc9.4/lib 
+ 
+ 
+FFLAGS   = -I$(INCPATH2) 
+ 
+LIBS = -lnetcdff 
+ 
+LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -Wl,-rpath,$(LIBPATH2) $(LIBS)  
+
+
+print:  printgrid_ncf printgrid printpart
+
+all: flex_read_compare2 printheader printgrid printpart
+
+%.o: %.f90
+	+$(FC) -c $(FFLAGS) $<		
+#	+$(FC) -c $<		
+flex_read_compare2:readheader.o readgrid.o juldate.o caldate.o flex_read_compare2.o 
+	$(FC) readheader.o readgrid.o juldate.o caldate.o flex_read_compare2.o -o flex_read_compare2
+printheader:readheader1090.o juldate.o caldate.o readheader_stdout90.o 
+	$(FC) readheader1090.o juldate.o caldate.o readheader_stdout90.o -o printheader
+printgrid:readheader.o readgrid.o juldate.o caldate.o printgrid.o 
+	$(FC) readheader.o readgrid.o juldate.o caldate.o printgrid.o -o printgrid
+printpart:work_arrays.o readpart.o printpart.o
+#	$(FC) -c work_arrays.f90 readpart.f printpart.f90 
+	$(FC) work_arrays.o readpart.o printpart.o -o printpart
+
+#printgrid_ncf:readheader.o readheader_ncf.o readgrid.o juldate.o caldate.o printgrid_ncf.o 
+#	$(FC) readheader.o readheader_ncf.o readgrid.o juldate.o caldate.o printgrid_ncf.o -o printgrid_ncf $(LDFLAGS)
+printgrid_ncf: readheader_ncf.o readgrid_ncf.o juldate.o caldate.o printgrid_ncf.o 
+	$(FC)  readheader_ncf.o readgrid_ncf.o juldate.o caldate.o printgrid_ncf.o -o printgrid_ncf $(LDFLAGS)
+clean:
+	        rm *.o
+		rm *.mod
+reset: 
+	rm  *.o *.mod printpart printgrid printgrid_ncf printheader flex_read_compare2
+
+output=../../output/
+test: 
+	./printheader $(output)
+	./printgrid $(output) conc

tests/flex_read_fortran/printgrid.f90

+    program flexpartread
+
+    parameter(nxmax=720,nymax=360,nzmax=40,maxtime=10)
+    parameter(maxpoint=1,maxspec=1,maxageclass=1)
+
+    double precision :: jul,julstart,juldate
+    real :: outheight(nzmax)
+    integer :: ireleasestart(maxpoint),ireleaseend(maxpoint)
+    integer :: lage(0:maxageclass)
+    real :: xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+    real :: zpoint2(maxpoint)
+    integer :: ireldate(maxspec),ireltime(maxspec)
+
+    real :: grid(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+    real :: footprint(0:nxmax-1,0:nymax-1)
+    real :: footprint_total(0:nxmax-1,0:nymax-1,maxspec)
+    real :: area(0:nxmax-1,0:nymax-1)
+    real :: heightnn(0:nxmax-1,0:nymax-1,0:nzmax)
+    real :: smallnum,r_sparse_dump(5000000)
+    integer :: i_sparse_dump(5000000)
+    parameter(smallnum=1.e-38)
+
+    character dirname*120
+    logical ::  dirnameisarg
+    character filename*150,fileout*150,aday*8,atime*6
+    character(250) :: filegrid
+          
+    logical :: sp_zer
+    logical :: debug=.FALSE.
+    
+    real :: sizegrid_m
+    real :: sum_grid, sum_sum_grid
+    real :: max_grid, max_max_grid
+    !character(4) :: grid_magnitude
+    character(6) :: grid_magnitude
+
+        max_max_grid=0
+        max_grid=0
+        sum_grid=0
+        sum_sum_grid=0
+        sizegrid_m=0
+
+
+    grid_magnitude='time'
+    grid_magnitude='conc'
+    grid_magnitude='pptv'
+
+    dirnameisarg= .FALSE. 
+
+    if (debug) then
+      print*, 'iargc', iargc()
+    endif
+
+    select case (iargc())
+    case (3)
+    debug=.TRUE.
+    dirnameisarg= .TRUE. 
+    call getarg(1,dirname)
+    len=index(dirname,' ')-1
+    call getarg(2,grid_magnitude)
+
+    case (2)
+    dirnameisarg= .TRUE. 
+    call getarg(1,dirname)
+    len=index(dirname,' ')-1
+    call getarg(2,grid_magnitude)
+    case (1)
+    call getarg(1,grid_magnitude)
+    case (0)
+    grid_magnitude='time'
+    end select
+    if (debug) then
+    print*, grid_magnitude
+    endif
+    
+    if ( .NOT. dirnameisarg) then
+    !******************
+        print*, " Read CONTROL file "
+    !******************
+        open(15,file='dirlist')
+        200 read(15,'(a)',end=199) dirname
+        len=index(dirname,' ')-1
+    endif
+
+!      do 50 mother_or_nest=1,2
+    do 50 mother_or_nest=1,1
+
+    !*********************************
+    ! Read FLEXPART header information
+    !*********************************
+
+        if (mother_or_nest == 1) then
+            filename=dirname(1:len)//'header'
+        else if (mother_or_nest == 2) then
+            filename=dirname(1:len)//'header_nest'
+        endif
+        call readheader(filename,nxmax,numxgrid,nymax,numygrid,nzmax, &
+        numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime, &
+        loutstep,maxspec,nspec,maxageclass,nageclass,lage,ireleasestart, &
+        ireleaseend,maxpoint,numpoint,xpoint,ypoint,zpoint1,zpoint2, &
+        heightnn,area)
+
+        if (debug) then
+          !write(*,*) 'nspec',nspec
+          write(*,*) 'simulation starts ', ibdate,ibtime
+          write(*,*) 'releases:'
+        endif
+
+        julstart=juldate(ibdate,ibtime)
+        do 30 k=1,nspec
+            jul=julstart+dble(float(ireleaseend(k))/86400.)
+            call caldate(jul,ireldate(k),ireltime(k))
+!            print*,'ireldate',ireldate(k)
+!             print*,'ireltime',ireltime(k)
+             if (debug) then 
+             print*, k, 'ireldate, ireltime',ireldate(k), ireltime(k)
+             endif
+             call caldate(julstart+dble(float(ireleasestart(k))/86400.) &
+             ,irelds,irelts)
+             if (debug) then
+             print*, k,irelds,irelts 
+             print*, k,xpoint(k),ypoint(k),zpoint1(k),zpoint2(k) 
+             !print*, k,irelds,irelts 
+             endif 
+
+        30 END DO
+        
+
+    ! rint*, 'init sizegrid_m=', sizegrid_m !=0
+
+        do 70 k=1,nspec
+            do 70 ix=0,numxgrid-1
+                do 70 jy=0,numygrid-1
+                    footprint_total(ix,jy,k)=0.
+        70 END DO
+
+    !******************************************
+    ! Loop about all times, given in file dates
+    !******************************************
+        if (debug) then
+          print*,'> open dates'
+        endif
+
+        open(20,file=dirname(1:len)//'dates',form='formatted', &
+        status='old')
+        it=0
+        100 read(20,'(i8,i6)',end=99) inday,intime
+        write(aday,'(i8.8)') inday
+        write(atime,'(i6.6)') intime
+        it=it+1
+
+        if (debug) then
+!         print*, 'it, aday, atime', it, aday, atime
+         print*, 'i, day, time', it, inday, intime
+         endif
+                 
+        if (mother_or_nest == 1) then
+                        
+        ! ilegrid=dirname(1:len)//'/grid_time_'//aday//atime//'_001'
+            filegrid=dirname(1:len)//'/grid_conc_'//aday//atime//'_001'
+        ! ilegrid=dirname(1:len)//'/grid_pptv_'//aday//atime//'_001'
+                        
+    !        filegrid=dirname(1:len)//'/grid_'//grid_magnitude//'_' &
+            filegrid=dirname(1:len)//'/grid_'//trim(grid_magnitude)//'_' &
+            //aday//atime//'_001'
+            fileout=dirname(1:len)//'footprint_'//aday//atime//'_001'
+        else
+            filegrid=dirname(1:len)//'grid_time_nest_'//aday//atime//'_001'
+            fileout=dirname(1:len)//'footprint_nest_'//aday//atime//'_001'
+        endif
+
+    !    print*, 'filegrid=', filegrid
+
+        call readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid, &
+        nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid, &
+        lage)
+
+        !print*, 'read grid:'
+        !print*, 'maxval(grid)=', maxval(grid)
+        !print*, 'sum(grid)=', sum(grid)
+        !print*, 'max, sum                 ', maxval(grid), sum(grid)
+
+    !***************************************
+    ! Display read values
+    !***************************************
+    ! izegrid_m=sizegrid_m+sum(grid)
+
+
+        sum_grid=sum(grid)
+        max_grid=maxval(grid)
+
+        if (debug) then
+          print*, 'max, sum                 ', max_grid, sum_grid
+        endif
+
+        max_max_grid=max(max_max_grid,max_grid)
+        sum_sum_grid=sum_sum_grid+sum_grid
+!        print*, 'max_max_grid,max_grid', max_max_grid, max_grid
+        if (debug) then
+          print*, 'max_max_grid,sum_sum_grid', max_max_grid, sum_sum_grid
+          !print*,
+        endif 
+!        print*, 'sum_grid=',sum_grid !sum(grid)
+!        print*, 'max_grid=', max_grid !  maxval(grid)
+!        print*, 'part total', it, aday, atime
+!        print*, 'sum_sum_grid=', sum_sum_grid
+!        print*, 'max_max_grid=', max_max_grid
+
+        if (0 == 1) then
+        !****************************************
+        ! Write out footprint files for each date
+        !****************************************
+                   
+            open(97,file=fileout,form='unformatted')
+            write(97) itime
+
+        ! Add contributions from this time step to gridded field
+        !*******************************************************
+
+            do 80 k=1,nspec
+                do 202 ix=0,numxgrid-1
+                    do 202 jy=0,numygrid-1
+                        footprint(ix,jy)=0.
+                        do 202 nage=1,nageclass
+                            footprint(ix,jy)=footprint(ix,jy)+ &
+                            grid(ix,jy,1,k,nage)
+                            footprint_total(ix,jy,k)=footprint_total(ix,jy,k)+ &
+                            grid(ix,jy,1,k,nage)
+                202 END DO
+
+
+                n_sp_count_i=0
+                n_sp_count_r=0
+                sp_fact=-1.
+                sp_zer= .TRUE. 
+                do 135 jy=0,numygrid-1
+                    do 135 ix=0,numxgrid-1
+                        if (footprint(ix,jy) > smallnum) then
+                            if (sp_zer) then ! first non zero value
+                                n_sp_count_i=n_sp_count_i+1
+                                i_sparse_dump(n_sp_count_i)= &
+                                ix+jy*numxgrid+1*numxgrid*numygrid   ! vertical level is 1
+                                sp_zer= .FALSE. 
+                                sp_fact=sp_fact*(-1)
+                            endif
+                            n_sp_count_r=n_sp_count_r+1
+                            r_sparse_dump(n_sp_count_r)=sp_fact*footprint(ix,jy)
+                        else ! concentration is zero
+                            sp_zer= .TRUE. 
+                        endif
+                135 END DO
+
+            ! First, zero wet and dry depo
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+
+                write(97) n_sp_count_i
+                write(97) (i_sparse_dump(i),i=1,n_sp_count_i)
+                write(97) n_sp_count_r
+                write(97) (r_sparse_dump(i),i=1,n_sp_count_r)
+
+            80 END DO
+            close(97)
+        endif ! Write out footprint files for each date
+
+        goto 100 !read dates
+        99 close(20)
+
+
+        if  (0 == 1) then
+        ! relic
+    ! Dump time-integrated footprint
+    !*******************************
+            if (mother_or_nest == 1) then
+                fileout=dirname(1:len)//'footprint_total'//'_001'
+            else
+                fileout=dirname(1:len)//'footprint_total_nest'//'_001'
+            endif
+
+            open(97,file=fileout,form='unformatted')
+            write(97) itime
+
+            do 180 k=1,nspec
+                n_sp_count_i=0
+                n_sp_count_r=0
+                sp_fact=-1.
+                sp_zer= .TRUE. 
+                do 145 jy=0,numygrid-1
+                    do 145 ix=0,numxgrid-1
+                        if (footprint_total(ix,jy,k) > smallnum) then
+                            if (sp_zer) then ! first non zero value
+                                n_sp_count_i=n_sp_count_i+1
+                                i_sparse_dump(n_sp_count_i)= &
+                                ix+jy*numxgrid+1*numxgrid*numygrid   ! vertical level is 1
+                                sp_zer= .FALSE. 
+                                sp_fact=sp_fact*(-1)
+                            endif
+                            n_sp_count_r=n_sp_count_r+1
+                            r_sparse_dump(n_sp_count_r)= &
+                            sp_fact*footprint_total(ix,jy,k)
+                        else ! concentration is zero
+                            sp_zer= .TRUE. 
+                        endif
+                145 END DO
+
+            ! First, zero wet and dry depo
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+
+                write(97) n_sp_count_i
+                write(97) (i_sparse_dump(i),i=1,n_sp_count_i)
+                write(97) n_sp_count_r
+                write(97) (r_sparse_dump(i),i=1,n_sp_count_r)
+
+            180 END DO
+            close(97)
+        endif ! write footprint total
+
+    50 END DO
+
+    if ( .NOT. dirnameisarg) then
+        goto 200 !read(15, dirlist
+        199 close(15) !dirlist
+    endif
+
+! ax_max_grid=max(max_max_grid,max_grid)
+! um_sum_grid=sum_sum_grid+sum_grid
+    if (debug) then
+
+    print*, 'for all times'
+    print*, 'sum sum(grid)=', sum_sum_grid
+    print*, 'max maxval(grid)=', max_max_grid
+    endif 
+    !if (debug) then
+
+    !print*, 'max:', max_max_grid, 'sum:', sum_sum_grid 
+    print*, 'max:', max_max_grid, 'mean:', sum_sum_grid/numxgrid/numygrid/numzgrid
+    !print*, numxgrid*numygrid*numzgrid
+    !print*, nxmax,numxgrid,nymax,numygrid,nzmax, numzgrid 
+
+!    print*, 'sum_sum_grid_out=',sum_sum_grid
+!    print*, 'max_max_grid_out=',max_max_grid
+    END PROGRAM

tests/flex_read_fortran/printgrid_ncf.f90

+    program flexpartread
+
+    use netcdf
+
+    parameter(nxmax=720,nymax=360,nzmax=40,maxtime=10)
+    parameter(maxpoint=1,maxspec=1,maxageclass=1)
+
+    double precision :: jul,julstart,juldate
+    real :: outheight(nzmax)
+    integer :: ireleasestart(maxpoint),ireleaseend(maxpoint)
+    integer :: lage(0:maxageclass)
+    real :: xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+    real :: zpoint2(maxpoint)
+    integer :: ireldate(maxspec),ireltime(maxspec)
+
+    real :: grid(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+    real :: footprint(0:nxmax-1,0:nymax-1)
+    real :: footprint_total(0:nxmax-1,0:nymax-1,maxspec)
+    real :: area(0:nxmax-1,0:nymax-1)
+    real :: heightnn(0:nxmax-1,0:nymax-1,0:nzmax)
+    real :: smallnum,r_sparse_dump(5000000)
+    integer :: i_sparse_dump(5000000)
+    parameter(smallnum=1.e-38)
+
+    character dirname*120
+    logical ::  dirnameisarg
+    character filename*150,fileout*150,aday*8,atime*6
+    character(250) :: filegrid
+    character(250) :: grid_filename
+          
+    logical :: sp_zer
+    logical :: debug=.FALSE.
+    
+    real :: sizegrid_m
+    real :: sum_grid, sum_sum_grid
+    real :: max_grid, max_max_grid
+    character(4) :: grid_magnitude
+
+    integer :: status, ncid, concvarid 
+    integer :: agedimid, latdimid, londimid, pointdimid, timedimid, heightdimid
+    integer ::  nlats, nlons, nheigths, ntimes, nages, npoints
+
+    real, allocatable, dimension(:,:,:,:,:,:) :: concvalue
+
+  
+        max_max_grid=0
+        max_grid=0
+        sum_grid=0
+        sum_sum_grid=0
+        sizegrid_m=0
+
+
+    grid_magnitude='time'
+    grid_magnitude='conc'
+    grid_magnitude='pptv'
+
+    dirnameisarg= .FALSE. 
+
+    if (debug) then
+      print*, 'iargc', iargc()
+    endif
+
+    select case (iargc())
+    case (3)
+    !debug=.TRUE.
+    dirnameisarg= .TRUE. 
+    call getarg(1,dirname)
+    len=index(dirname,' ')-1
+    call getarg(2,grid_magnitude)
+    call getarg(3,grid_filename)
+    case (2)
+    dirnameisarg= .TRUE. 
+    call getarg(1,dirname)
+    len=index(dirname,' ')-1
+    call getarg(2,grid_magnitude)
+    case (1)
+    call getarg(1,grid_magnitude)
+    case (0)
+    grid_magnitude='time'
+    end select
+    if (debug) then
+    print*, grid_magnitude
+    endif
+    
+    if ( .NOT. dirnameisarg) then
+    !******************
+        print*, " Read dirname file "
+    !******************
+        open(15,file='dirlist')
+
+        200 read(15,'(a)',end=199) dirname
+        len=index(dirname,' ')-1
+    endif
+
+!      do 50 mother_or_nest=1,2
+    do 50 mother_or_nest=1,1
+
+    !*********************************
+    ! Read FLEXPART header information
+    !*********************************
+
+
+        !filename=dirname(1:len)//'grid_conc_20190225190000.nc'
+        filename=dirname(1:len)//trim(grid_filename) !/grid_magnitude
+  
+        !print*,'test read ncf' 
+        !print*, 'filename=', filename  
+
+        status = nf90_open(filename, NF90_NOWRITE, ncid) 
+        !print*,'status:',status
+        call nf90_err(status) 
+      
+
+       ! nf90_inq_dimid(ncid, "latitude", latdimid) londimid, pointdimid, timedimid, heightdimid
+
+      status =  nf90_inq_dimid(ncid, "latitude", latdimid)
+      status =  nf90_inq_dimid(ncid, "longitude", londimid)
+      status =   nf90_inq_dimid(ncid, "height", heightdimid)
+      status =   nf90_inq_dimid(ncid, "time", timedimid)
+      status =    nf90_inq_dimid(ncid, "nageclass", agedimid)
+      status =    nf90_inq_dimid(ncid, "pointspec", pointdimid)
+
+       status = nf90_inquire_dimension(ncid, latdimid, len = nlats)
+
+       status = nf90_inquire_dimension(ncid, londimid, len = nlons)
+
+
+       status = nf90_inquire_dimension(ncid, heightdimid, len = nheigths)
+       status = nf90_inquire_dimension(ncid, timedimid, len = ntimes)
+       status = nf90_inquire_dimension(ncid, agedimid, len = nages)
+       status = nf90_inquire_dimension(ncid, pointdimid, len = npoints)
+
+       !print*, nlats, nlons, nheigths, ntimes, nages, npoints
+
+       !allocate(concvalue(nages,npoints,ntimes,nheigths,nlons,nlats))
+       allocate(concvalue(nlons,nlats,nheigths,ntimes, npoints, nages))
+
+
+!dimensions:
+!	time = UNLIMITED ; // (2 currently)
+!	longitude = 10 ;
+!	latitude = 10 ;
+!	height = 4 ;
+!	numspec = 1 ;
+!	pointspec = 1 ;
+!	nageclass = 1 ;
+!	nchar = 45 ;
+!	ncharrec = 16 ;
+!	numpoint = 1 ;
+ 
+       !	float spec001_mr(nageclass, pointspec, time, height, latitude, longitude) ;
+
+    if (grid_magnitude(1:4).eq."conc") then
+     
+        status = nf90_inq_varid(ncid, "spec001_mr", concvarid)
+   elseif (grid_magnitude(1:4).eq."pptv") then  
+        status = nf90_inq_varid(ncid, "spec001_pptv", concvarid)
+   elseif (grid_magnitude(1:4).eq."time") then  
+        status = nf90_inq_varid(ncid, "spec001_mr", concvarid)
+    else
+         print*, grid_magnitude 
+         stop 'unit not supported'
+   endif
+       
+
+        !print*, concvarid
+
+         status = nf90_get_var(ncid, concvarid, concvalue) 
+         call nf90_err(status)
+
+
+        status = nf90_close(ncid) 
+        call nf90_err(status) 
+
+
+
+        sum_grid=sum(concvalue)
+        max_grid=maxval(concvalue)
+
+     
+          !print*, 'max, sum, ave :                 ', max_grid, sum_grid, sum_grid/nlons/nlats/nheigths
+          !binary format
+          print*, 'max:', max_grid, 'mean:', sum_grid/nlons/nlats/nheigths ! sum_grid/numxgrid/numygrid/numzgrid
+
+
+      stop
+
+
+        if (mother_or_nest == 1) then
+            filename=dirname(1:len)//'header'
+        else if (mother_or_nest == 2) then
+            filename=dirname(1:len)//'header_nest'
+        endif
+        call readheader(filename,nxmax,numxgrid,nymax,numygrid,nzmax, &
+        numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime, &
+        loutstep,maxspec,nspec,maxageclass,nageclass,lage,ireleasestart, &
+        ireleaseend,maxpoint,numpoint,xpoint,ypoint,zpoint1,zpoint2, &
+        heightnn,area)
+
+        if (debug) then
+          !write(*,*) 'nspec',nspec
+          write(*,*) 'simulation starts ', ibdate,ibtime
+          write(*,*) 'releases:'
+        endif
+
+        julstart=juldate(ibdate,ibtime)
+        do 30 k=1,nspec
+            jul=julstart+dble(float(ireleaseend(k))/86400.)
+            call caldate(jul,ireldate(k),ireltime(k))
+!            print*,'ireldate',ireldate(k)
+!             print*,'ireltime',ireltime(k)
+             if (debug) then 
+             print*, k, 'ireldate, ireltime',ireldate(k), ireltime(k)
+             endif
+             call caldate(julstart+dble(float(ireleasestart(k))/86400.) &
+             ,irelds,irelts)
+             if (debug) then
+             print*, k,irelds,irelts 
+             print*, k,xpoint(k),ypoint(k),zpoint1(k),zpoint2(k) 
+             !print*, k,irelds,irelts 
+             endif 
+
+        30 END DO
+        
+
+    ! rint*, 'init sizegrid_m=', sizegrid_m !=0
+
+        do 70 k=1,nspec
+            do 70 ix=0,numxgrid-1
+                do 70 jy=0,numygrid-1
+                    footprint_total(ix,jy,k)=0.
+        70 END DO
+
+    !******************************************
+    ! Loop about all times, given in file dates
+    !******************************************
+        if (debug) then
+          print*,'> open dates'
+        endif
+
+        open(20,file=dirname(1:len)//'dates',form='formatted', &
+        status='old')
+        it=0
+        100 read(20,'(i8,i6)',end=99) inday,intime
+        write(aday,'(i8.8)') inday
+        write(atime,'(i6.6)') intime
+        it=it+1
+
+        if (debug) then
+!         print*, 'it, aday, atime', it, aday, atime
+         print*, 'i, day, time', it, inday, intime
+         endif
+                 
+        if (mother_or_nest == 1) then
+                        
+        ! ilegrid=dirname(1:len)//'/grid_time_'//aday//atime//'_001'
+            filegrid=dirname(1:len)//'/grid_conc_'//aday//atime//'_001'
+        ! ilegrid=dirname(1:len)//'/grid_pptv_'//aday//atime//'_001'
+                        
+            filegrid=dirname(1:len)//'/grid_'//grid_magnitude//'_' &
+            //aday//atime//'_001'
+            fileout=dirname(1:len)//'footprint_'//aday//atime//'_001'
+        else
+            filegrid=dirname(1:len)//'grid_time_nest_'//aday//atime//'_001'
+            fileout=dirname(1:len)//'footprint_nest_'//aday//atime//'_001'
+        endif
+
+        print*, 'filegrid=', filegrid
+
+        call readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid, &
+        nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid, &
+        lage)
+
+        print*, 'read grid:'
+        print*, 'maxval(grid)=', maxval(grid)
+        print*, 'sum(grid)=', sum(grid)
+        print*, 'max, sum                 ', maxval(grid), sum(grid)
+
+    !***************************************
+    ! Display read values
+    !***************************************
+    ! izegrid_m=sizegrid_m+sum(grid)
+
+
+        sum_grid=sum(grid)
+        max_grid=maxval(grid)
+
+        if (debug) then
+          print*, 'max, sum                 ', max_grid, sum_grid
+        endif
+
+        max_max_grid=max(max_max_grid,max_grid)
+        sum_sum_grid=sum_sum_grid+sum_grid
+!        print*, 'max_max_grid,max_grid', max_max_grid, max_grid
+        if (debug) then
+          print*, 'max_max_grid,sum_sum_grid', max_max_grid, sum_sum_grid
+          print*, ' ' 
+        endif 
+!        print*, 'sum_grid=',sum_grid !sum(grid)
+!        print*, 'max_grid=', max_grid !  maxval(grid)
+!        print*, 'part total', it, aday, atime
+!        print*, 'sum_sum_grid=', sum_sum_grid
+!        print*, 'max_max_grid=', max_max_grid
+
+        if (0 == 1) then
+        !****************************************
+        ! Write out footprint files for each date
+        !****************************************
+                   
+            open(97,file=fileout,form='unformatted')
+            write(97) itime
+
+        ! Add contributions from this time step to gridded field
+        !*******************************************************
+
+            do 80 k=1,nspec
+                do 202 ix=0,numxgrid-1
+                    do 202 jy=0,numygrid-1
+                        footprint(ix,jy)=0.
+                        do 202 nage=1,nageclass
+                            footprint(ix,jy)=footprint(ix,jy)+ &
+                            grid(ix,jy,1,k,nage)
+                            footprint_total(ix,jy,k)=footprint_total(ix,jy,k)+ &
+                            grid(ix,jy,1,k,nage)
+                202 END DO
+
+
+                n_sp_count_i=0
+                n_sp_count_r=0
+                sp_fact=-1.
+                sp_zer= .TRUE. 
+                do 135 jy=0,numygrid-1
+                    do 135 ix=0,numxgrid-1
+                        if (footprint(ix,jy) > smallnum) then
+                            if (sp_zer) then ! first non zero value
+                                n_sp_count_i=n_sp_count_i+1
+                                i_sparse_dump(n_sp_count_i)= &
+                                ix+jy*numxgrid+1*numxgrid*numygrid   ! vertical level is 1
+                                sp_zer= .FALSE. 
+                                sp_fact=sp_fact*(-1)
+                            endif
+                            n_sp_count_r=n_sp_count_r+1
+                            r_sparse_dump(n_sp_count_r)=sp_fact*footprint(ix,jy)
+                        else ! concentration is zero
+                            sp_zer= .TRUE. 
+                        endif
+                135 END DO
+
+            ! First, zero wet and dry depo
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+
+                write(97) n_sp_count_i
+                write(97) (i_sparse_dump(i),i=1,n_sp_count_i)
+                write(97) n_sp_count_r
+                write(97) (r_sparse_dump(i),i=1,n_sp_count_r)
+
+            80 END DO
+            close(97)
+        endif ! Write out footprint files for each date
+
+        goto 100 !read dates
+        99 close(20)
+
+
+        if  (0 == 1) then
+        ! relic
+    ! Dump time-integrated footprint
+    !*******************************
+            if (mother_or_nest == 1) then
+                fileout=dirname(1:len)//'footprint_total'//'_001'
+            else
+                fileout=dirname(1:len)//'footprint_total_nest'//'_001'
+            endif
+
+            open(97,file=fileout,form='unformatted')
+            write(97) itime
+
+            do 180 k=1,nspec
+                n_sp_count_i=0
+                n_sp_count_r=0
+                sp_fact=-1.
+                sp_zer= .TRUE. 
+                do 145 jy=0,numygrid-1
+                    do 145 ix=0,numxgrid-1
+                        if (footprint_total(ix,jy,k) > smallnum) then
+                            if (sp_zer) then ! first non zero value
+                                n_sp_count_i=n_sp_count_i+1
+                                i_sparse_dump(n_sp_count_i)= &
+                                ix+jy*numxgrid+1*numxgrid*numygrid   ! vertical level is 1
+                                sp_zer= .FALSE. 
+                                sp_fact=sp_fact*(-1)
+                            endif
+                            n_sp_count_r=n_sp_count_r+1
+                            r_sparse_dump(n_sp_count_r)= &
+                            sp_fact*footprint_total(ix,jy,k)
+                        else ! concentration is zero
+                            sp_zer= .TRUE. 
+                        endif
+                145 END DO
+
+            ! First, zero wet and dry depo
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+                write(97) 0
+                write(97)
+
+                write(97) n_sp_count_i
+                write(97) (i_sparse_dump(i),i=1,n_sp_count_i)
+                write(97) n_sp_count_r
+                write(97) (r_sparse_dump(i),i=1,n_sp_count_r)
+
+            180 END DO
+            close(97)
+        endif ! write footprint total
+
+    50 END DO
+
+    if ( .NOT. dirnameisarg) then
+        goto 200 !read(15, dirlist
+        199 close(15) !dirlist
+    endif
+
+! ax_max_grid=max(max_max_grid,max_grid)
+! um_sum_grid=sum_sum_grid+sum_grid
+    if (debug) then
+
+    print*, 'for all times'
+    print*, 'sum sum(grid)=', sum_sum_grid
+    print*, 'max maxval(grid)=', max_max_grid
+    endif 
+    !if (debug) then
+
+    !print*, 'max:', max_max_grid, 'sum:', sum_sum_grid 
+    print*, 'max:', max_max_grid, 'mean:', sum_sum_grid/numxgrid/numygrid/numzgrid
+    !print*, numxgrid*numygrid*numzgrid
+    !print*, nxmax,numxgrid,nymax,numygrid,nzmax, numzgrid 
+
+!    print*, 'sum_sum_grid_out=',sum_sum_grid
+!    print*, 'max_max_grid_out=',max_max_grid
+    END PROGRAM
+
+
+subroutine nf90_err(status)
+  !****************************************************************
+  ! netcdf error message handling
+  !****************************************************************
+  use netcdf
+  implicit none
+  integer, intent (in) :: status
+   if(status /= nf90_noerr) then
+      print *, trim(nf90_strerror(status))
+      stop 'Stopped'
+    end if
+end subroutine nf90_err
+

tests/flex_read_fortran/printpart.f90

+program printpart
+
+use work_arrays
+
+CHARACTER(len=150)     :: flexout
+CHARACTER (len=14), dimension(:), ALLOCATABLE  :: dates
+CHARACTER (len=14) :: date
+integer numpart !,nspec
+real :: mx,my,mz
+
+
+
+  flexout='/Users/ignacio/repos/flexpart/tests/examples/output_3-part1/'
+  !date='20120101080000'
+  !date='20120101090000'
+  date='20120101100000'
+  numpart=-1
+  nspec=1
+
+select case (iargc())
+case (4)
+  call getarg(1,flexout)  
+  call getarg(2,date)  
+  !call getarg(3,nspec)  
+  !call getarg(4,numpart)  
+case (3)
+  call getarg(1,flexout)  
+  call getarg(2,date)  
+  !call getarg(3,nspec)  
+case (2)
+  call getarg(1,flexout)  
+  call getarg(2,date)  
+case (1)
+  call getarg(1,flexout) 
+  date='0000'
+case (0)
+  print*,'0 arguments, default'
+  print*, flexout
+  print*, date
+
+end select
+
+!CALL readpart2(flexout,dates(idate),numpart)
+!CALL readpart2(flexout,date,numpart)
+CALL readpart(flexout,date,numpart)
+
+if (.false.) then
+print*, 'numpart read:', numpart
+numpart=5
+endif
+
+ALLOCATE(vnpoint(numpart),vitramem(numpart) )
+ALLOCATE(vxlon(numpart),vylat(numpart) ,vztra1(numpart), &
+          vtri(numpart),vhmixi(numpart),vtopo(numpart) , &
+          vpvi(numpart),vtti(numpart),                   &
+          vqvi(numpart),vrhoi(numpart)) 
+ALLOCATE(xmass0(numpart,nspec))
+
+CALL readpart(flexout,date,numpart)
+
+if (.false.) then
+print*, vnpoint
+print*, vxlon
+print*, vylat
+print*, vztra1
+endif
+
+mx = SUM(vxlon)/SIZE(vxlon) 
+my = SUM(vylat)/SIZE(vylat) 
+mz = SUM(vztra1)/SIZE(vztra1) 
+
+if (.false.) then
+print*, 'mean x = ', mx
+print*, 'mean y = ', my
+print*, 'mean z = ', mz
+endif
+
+print*, mx,my,mz
+
+
+!print*, SUM((vxlon- )**2)/SIZE(Array)
+
+end program

tests/flex_read_fortran/readgrid.f

+      subroutine readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid,
+     +nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid,lage)
+
+      real grid(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+      integer lage(0:maxageclass)
+
+***************** new sparse output
+      real smallnum
+      integer fact,ii,ir
+      real sparse_dump_r(nxmax*nymax*nzmax)
+      integer sparse_dump_i(nxmax*nymax*nzmax)
+      integer sp_count_i,sp_count_r
+***************** new sparse output
+
+
+      character*250 filegrid
+
+      open(10,file=filegrid,form='unformatted',status='old')
+      read(10,end=99) itime 
+      !write(*,*) itime
+      !print*, 'readgrid> itime=', itime
+
+C Loop about all species
+************************
+
+      do 38 k=1,nspec
+
+C Loop about all age classes
+****************************
+
+        do 38 nage=1,nageclass
+
+C Initialize age class fields
+*****************************
+
+        do 63 ix=0,numxgrid-1
+          do 63 jy=0,numygrid-1
+            do 63 kz=1,numzgrid
+63            grid(ix,jy,kz,k,nage)=0.
+
+
+C Read wet deposition
+*********************
+
+        fact=1
+        read(10) sp_count_i
+        read(10) (sparse_dump_i(ix),ix=1,sp_count_i)
+        read(10) sp_count_r
+        read(10) (sparse_dump_r(ix),ix=1,sp_count_r)
+
+c       ii=0
+c       do 32 ir=1,sp_count_r
+c         if ((sparse_dump_r(ir)*fact).gt.smallnum) then
+c            ii=ii+1
+c            n=sparse_dump_i(ii)
+c            fact=fact*(-1.)
+c         else
+c            pos=pos+1
+c         endif
+
+c         jy=n/numxgrid
+c         ix=n-numxgrid*jy
+c         wetgrid(ix,jy,k,nage)=sparse_dump_r(ir)
+
+c2      continue
+
+C Read dry deposition
+*********************
+
+        fact=1
+        read(10) sp_count_i
+        read(10) (sparse_dump_i(ix),ix=1,sp_count_i)
+        read(10) sp_count_r
+        read(10) (sparse_dump_r(ix),ix=1,sp_count_r)
+
+c       ii=0
+c       do 36 ir=1,sp_count_r
+c         if ((sparse_dump_r(ir)*fact).gt.smallnum) then
+c            ii=ii+1
+c            n=sparse_dump_i(ii)
+c            fact=fact*(-1.)
+c         else
+c            pos=pos+1
+c         endif
+
+c         jy=n/numxgrid
+c         ix=n-numxgrid*jy
+c         drygrid(ix,jy,k,nage)=sparse_dump_r(ir)
+
+c6      continue
+
+
+C Read concentrations
+*********************
+
+        fact=1
+        read(10) sp_count_i
+        read(10) (sparse_dump_i(ix),ix=1,sp_count_i)
+        read(10) sp_count_r
+        read(10) (sparse_dump_r(ix),ix=1,sp_count_r)
+c       write(*,*) sp_count_i,sp_count_r
+c       write(*,*) (sparse_dump_i(ix),ix=1,sp_count_i)
+c       write(*,*) (sparse_dump_r(ix),ix=1,sp_count_r)
+
+        ii=0
+        do 47 ir=1,sp_count_r
+          if ((sparse_dump_r(ir)*fact).gt.smallnum) then
+             ii=ii+1
+             n=sparse_dump_i(ii)
+             fact=fact*(-1.)
+          else
+             n=n+1
+          endif
+
+          kz=n/(numxgrid*numygrid)
+          jy=(n-kz*numxgrid*numygrid)/numxgrid
+          ix=n-numxgrid*numygrid*kz-numxgrid*jy
+          grid(ix,jy,kz,k,nage)=abs(sparse_dump_r(ir))
+
+47      continue
+
+C End species loop, end age class loop
+**************************************
+
+138       continue
+38      continue
+
+      close(10)
+
+99    continue
+
+
+      end

tests/flex_read_fortran/readgrid_ncf.f

+      subroutine readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid,
+     +nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid,lage)
+
+      real grid(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+      integer lage(0:maxageclass)
+
+***************** new sparse output
+      real smallnum
+      integer fact,ii,ir
+      real sparse_dump_r(nxmax*nymax*nzmax)
+      integer sparse_dump_i(nxmax*nymax*nzmax)
+      integer sp_count_i,sp_count_r
+***************** new sparse output
+
+
+      character*250 filegrid
+
+      open(10,file=filegrid,form='unformatted',status='old')
+      read(10,end=99) itime 
+      !write(*,*) itime
+      !print*, 'readgrid> itime=', itime
+
+C Loop about all species
+************************
+
+      do 38 k=1,nspec
+
+C Loop about all age classes
+****************************
+
+        do 38 nage=1,nageclass
+
+C Initialize age class fields
+*****************************
+
+        do 63 ix=0,numxgrid-1
+          do 63 jy=0,numygrid-1
+            do 63 kz=1,numzgrid
+63            grid(ix,jy,kz,k,nage)=0.
+
+
+C Read wet deposition
+*********************
+
+        fact=1
+        read(10) sp_count_i
+        read(10) (sparse_dump_i(ix),ix=1,sp_count_i)
+        read(10) sp_count_r
+        read(10) (sparse_dump_r(ix),ix=1,sp_count_r)
+
+c       ii=0
+c       do 32 ir=1,sp_count_r
+c         if ((sparse_dump_r(ir)*fact).gt.smallnum) then
+c            ii=ii+1
+c            n=sparse_dump_i(ii)
+c            fact=fact*(-1.)
+c         else
+c            pos=pos+1
+c         endif
+
+c         jy=n/numxgrid
+c         ix=n-numxgrid*jy
+c         wetgrid(ix,jy,k,nage)=sparse_dump_r(ir)
+
+c2      continue
+
+C Read dry deposition
+*********************
+
+        fact=1
+        read(10) sp_count_i
+        read(10) (sparse_dump_i(ix),ix=1,sp_count_i)
+        read(10) sp_count_r
+        read(10) (sparse_dump_r(ix),ix=1,sp_count_r)
+
+c       ii=0
+c       do 36 ir=1,sp_count_r
+c         if ((sparse_dump_r(ir)*fact).gt.smallnum) then
+c            ii=ii+1
+c            n=sparse_dump_i(ii)
+c            fact=fact*(-1.)
+c         else
+c            pos=pos+1
+c         endif
+
+c         jy=n/numxgrid
+c         ix=n-numxgrid*jy
+c         drygrid(ix,jy,k,nage)=sparse_dump_r(ir)
+
+c6      continue
+
+
+C Read concentrations
+*********************
+
+        fact=1
+        read(10) sp_count_i
+        read(10) (sparse_dump_i(ix),ix=1,sp_count_i)
+        read(10) sp_count_r
+        read(10) (sparse_dump_r(ix),ix=1,sp_count_r)
+c       write(*,*) sp_count_i,sp_count_r
+c       write(*,*) (sparse_dump_i(ix),ix=1,sp_count_i)
+c       write(*,*) (sparse_dump_r(ix),ix=1,sp_count_r)
+
+        ii=0
+        do 47 ir=1,sp_count_r
+          if ((sparse_dump_r(ir)*fact).gt.smallnum) then
+             ii=ii+1
+             n=sparse_dump_i(ii)
+             fact=fact*(-1.)
+          else
+             n=n+1
+          endif
+
+          kz=n/(numxgrid*numygrid)
+          jy=(n-kz*numxgrid*numygrid)/numxgrid
+          ix=n-numxgrid*numygrid*kz-numxgrid*jy
+          grid(ix,jy,kz,k,nage)=abs(sparse_dump_r(ir))
+
+47      continue
+
+C End species loop, end age class loop
+**************************************
+
+138       continue
+38      continue
+
+      close(10)
+
+99    continue
+
+
+      end

tests/flex_read_fortran/readheader.f

+      subroutine readheader(filename,nxmax,numxgrid,nymax,numygrid,
+     +nzmax,numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,
+     +ibtime,loutstep,maxspec,nspec,maxageclass,nageclass,lage,
+     +ireleasestart,ireleaseend,maxpoint,numpoint,xpoint,ypoint,
+     +zpoint1,zpoint2,heightnn,area)
+
+      parameter(pi=3.14159265,r_earth=6.371e6,pih=pi/180.)
+
+      character filename*150,compoint(maxpoint)*45,species(maxspec)*7
+      real outheight(nzmax),heightnn(nxmax,nymax,0:nzmax)
+      integer ireleasestart(maxpoint),ireleaseend(maxpoint)
+      integer npart(maxpoint),kind(maxpoint),lage(0:maxageclass)
+      real xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+      real zpoint2(maxpoint),xmass(maxpoint,maxspec)
+      real oro(nxmax,nymax),area(nxmax,nymax)
+      logical debug
+
+      debug=.FALSE.
+
+      open(10,file=filename,form='unformatted',status='old')
+      read(10) ibdate,ibtime
+
+      if (debug) then 
+      write(*,*) ibdate,ibtime
+      endif
+
+      read(10) loutstep,loutaver,loutsample
+      read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+
+      if (debug) then
+      write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+      endif
+
+      read(10) numzgrid,(outheight(i),i=1,numzgrid)
+      read(10) jjjjmmdd,ihmmss
+
+      if (debug) then
+      write(*,*) jjjjmmdd,ihmss
+      endif
+      
+      read(10) nspec
+      nspec=nspec/3
+      do 8 n=1,nspec
+        read(10) numzgrid,species(n)
+        read(10) numzgrid,species(n)
+8       read(10) numzgrid,species(n)
+
+
+      read(10) numpoint
+      do 13 i=1,numpoint
+        read(10) ireleasestart(i),ireleaseend(i)
+
+        read(10) xpoint(i),ypoint(i),xp2,yp2,zpoint1(i),zpoint2(i)
+        read(10) npart(i),kind(i)
+        read(10) compoint(i)
+        do 13 j=1,nspec
+          read(10)
+          read(10)
+13        read(10) xmass(i,j)
+      read(10) method
+      read(10) nageclass,(lage(i),i=1,nageclass)
+      lage(0)=0
+
+      if (debug) then
+      write(*,*) (lage(i),i=0,nageclass)
+      endif
+
+      do 130 ix=1,numxgrid
+130     read(10) (oro(ix,jy),jy=1,numygrid)
+
+      close(10)
+      if (debug) then
+      write(*,*) (outheight(i),i=1,numzgrid)
+      endif
+      
+      if (loutstep.lt.0) nspec=numpoint
+
+
+c Calculate height, which is outheight plus topography
+******************************************************
+
+      do 150 ix=1,numxgrid
+        do 150 jy=1,numygrid
+          if (ltopo.eq.1) then
+            heightnn (ix,jy,0) = oro(ix,jy)
+          else
+            heightnn (ix,jy,0) = 0.
+          endif
+          do 150 i=1,numzgrid
+            if (ltopo.eq.1) then
+              heightnn (ix,jy,i) = outheight(i) + oro(ix,jy)
+            else
+              heightnn (ix,jy,i) = outheight(i)
+            endif
+150         continue
+
+
+C Determine area of each output grid box
+****************************************
+
+      do 140 jy=1,numygrid
+        ylata=outlat0+(float(jy-1)+0.5)*dyout
+        ylatp=ylata+0.5*dyout
+        ylatm=ylata-0.5*dyout
+        if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
+          hzone=dyout*r_earth*pih
+        else
+          cosfact=cos(ylata*pih)*r_earth
+          cosfactp=cos(ylatp*pih)*r_earth
+          cosfactm=cos(ylatm*pih)*r_earth
+          if (cosfactp.lt.cosfactm) then
+            hzone=sqrt(r_earth**2-cosfactp**2)-
+     +      sqrt(r_earth**2-cosfactm**2)
+          else
+            hzone=sqrt(r_earth**2-cosfactm**2)-
+     +      sqrt(r_earth**2-cosfactp**2)
+          endif
+        endif
+        gridarea=2.*pi*r_earth*hzone*dxout/360.
+        do 140 ix=1,numxgrid
+140       area(ix,jy)=gridarea
+
+
+      end

tests/flex_read_fortran/readheader1090.f90

+    subroutine readheader10f(filename, &
+    nxmax,numxgrid,nymax,numygrid,nzmax,numzgrid, &
+    outlon0,outlat0,dxout,dyout,outheight,ibdate, &
+    ibtime,loutstep,maxspec,nspec,maxageclass,nageclass,lage, &
+    ireleasestart,ireleaseend,maxpoint,numpoint,xpoint,ypoint, &
+    zpoint1,zpoint2,heightnn,area)
+
+    parameter(pi=3.14159265,r_earth=6.371e6,pih=pi/180.)
+
+    character filename*150,compoint(maxpoint)*45,species(maxspec)*7
+    real :: outheight(nzmax),heightnn(nxmax,nymax,0:nzmax)
+! eal outheight(0:nzmax),heightnn(nxmax,nymax,0:nzmax) !v8.2
+    integer :: ireleasestart(maxpoint),ireleaseend(maxpoint)
+    integer :: npart(maxpoint),kind(maxpoint),lage(0:maxageclass)
+    real :: xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+    real :: zpoint2(maxpoint),xmass(maxpoint,maxspec)
+    real :: oro(nxmax,nymax),area(nxmax,nymax)
+          
+    character(len=256) :: pathfile, flexversion
+    
+    integer :: verbosity
+
+    verbosity=1
+
+    if (verbosity.gt.0) then
+    print*,'subroutine readheader10f'
+    endif 
+    
+    if (verbosity.gt.0) then
+        print*,'intent in:'
+        print*,'nxmax=', nxmax
+        print*,'nymax=', nymax
+        print*,'nzmax=', nzmax    
+    endif
+ 
+
+
+
+    open(10,file=filename,form='unformatted',status='old')
+
+          
+! 76
+    if (verbosity.gt.0) then
+    print*,'! Open header output file: '
+    print*, '**********************************************'
+    print*, filename    
+    endif
+    
+    read(10) ibdate,ibtime, flexversion
+
+    if (verbosity.gt.0) then
+    write(*,*) '! Write the header information'
+    print*, '**********************************************'
+    print*,'(ibdate,ibtime if fwd. -- iedate,ietime if bwd.)'
+    write(*,*) 'date,time=',ibdate,ibtime
+    print*, 'flexversion=',flexversion
+    endif
+    
+! 85
+    read(10) loutstep,loutaver,loutsample
+    if (verbosity.gt.0) then
+    print*,' info on output interval, averaging time, sampling time' 
+    print*, '**********************************************'
+    print*,'loutstep, aver, sample=',loutstep,loutaver,loutsample
+    endif
+
+!87
+    read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+    if (verbosity.gt.0) then
+    print*, 'information on output grid setup'
+    print*, '**********************************************'
+    write(*,*) 'outlon0,outlat0,numxgrid,numygrid,dxout,dyout'
+    write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+    endif
+
+    read(10) numzgrid,(outheight(i),i=1,numzgrid)
+    if (verbosity.gt.0) then  
+    write(*,*) 'numzgrid=', numzgrid
+!    write(*,*) 'outheight=', outheight
+        print*,'outheight='
+        do i=1,numzgrid
+        print*, i, outheight(i)
+        end do   
+    endif 
+
+    read(10) jjjjmmdd,ihmmss
+    if (verbosity.gt.0) then
+    write(*,*) 'jjjjmmdd,ihmmss=',jjjjmmdd,ihmss
+    endif
+
+  ! Write number of species, and name for each species (+extra name for depositions)
+  ! Indicate the dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+  ! concentration fields
+  !*****************************************************************************
+    if (verbosity.gt.0) then
+    print*, 'information on species'
+    print*, '**********************************************'    
+    endif 
+! 102  3*nspec,maxpointspec_act implicit?
+    read(10) nspec, maxpointspec_act
+    if (verbosity.gt.0) then
+    print*,'3*nspec,maxpointspec_act=',nspec,maxpointspec_act
+    endif
+    nspec=nspec/3
+    do 8 n=1,nspec
+        read(10) numzgrid,species(n)
+        if (verbosity.gt.0) then
+        print*,'numzgrid,species(n)',numzgrid,species(n)
+        endif
+        read(10) numzgrid,species(n)
+        if (verbosity.gt.0) then
+        print*,'numzgrid,species(n)',numzgrid,species(n)
+        endif
+        read(10) numzgrid,species(n)
+        if (verbosity.gt.0) then
+        print*,'numzgrid,species(n)',numzgrid,species(n)
+        endif
+    8 END DO
+
+    if (verbosity.gt.0) then
+    print*, 'information on release points:'
+    print*, '**********************************************'    
+    endif
+    
+    read(10) numpoint !l113
+    if (verbosity.gt.0) then
+    print*,'numpoint=',numpoint
+    endif
+
+    do 13 i=1,numpoint
+        read(10) ireleasestart(i),ireleaseend(i)
+        read(10) xpoint(i),ypoint(i),xp2,yp2,zpoint1(i),zpoint2(i)
+        read(10) npart(i),kind(i)
+        read(10) compoint(i)
+        do 13 j=1,nspec
+            read(10)
+            read(10)
+            read(10) xmass(i,j)
+    13 END DO
+    
+    
+
+
+    if (verbosity.gt.0) then
+    print*,'! Write information on some model switches'
+    print*, '**********************************************'
+    
+    endif
+    read(10) method,lsubgrid,lconvection,ind_source,ind_receptor
+    if (verbosity.gt.0) then
+    print*,'method,lsubgrid,lconvection,ind_source,ind_receptor'
+    print*,method,lsubgrid,lconvection,ind_source,ind_receptor
+    endif
+
+
+  ! Write age class information
+  !****************************
+    read(10) nageclass,(lage(i),i=1,nageclass)
+    lage(0)=0
+    if (verbosity.gt.0) then
+    print*,'!Write age class information'
+    print*, '**********************************************'
+    endif
+    if (verbosity.gt.0) then
+    print*,'nageclass=',nageclass
+    write(*,*) 'lage(0:nageclass)=', (lage(i),i=0,nageclass)
+    endif
+    
+
+  ! Write topography to output file
+  !********************************
+
+    
+    do 130 ix=1,numxgrid
+        read(10) (oro(ix,jy),jy=1,numygrid)
+    130 END DO
+    
+    if (verbosity.gt.0) then
+    print*,'!Write topography to output file' 
+    print*, '**********************************************'
+    
+    print*,'(oro(ix,jy)'
+    print*,'numxgrid,numygrid=',numxgrid,numygrid
+    print*,'nxmax,nymax=',nxmax,nymax
+    endif      
+    close(10)
+
+
+
+! rite(*,*) (outheight(i),i=1,numzgrid)
+! f (loutstep.lt.0) nspec=numpoint
+
+    if (verbosity.gt.0) then
+    print*,' !if (loutstep < 0) nspec=numpoint'
+    print*,'outheight(0)=0. ???'
+    endif
+    
+! Calculate height, which is outheight plus topography
+!*****************************************************
+
+    if (verbosity.gt.0) then     
+    print*,'Calculate height, which is outheight plus topography'
+    print*, '**********************************************'
+    endif 
+
+    do 150 ix=1,numxgrid
+        do 150 jy=1,numygrid
+            if (ltopo == 1) then
+                heightnn (ix,jy,0) = oro(ix,jy)
+            else
+                heightnn (ix,jy,0) = 0.
+            endif
+            do 150 i=1,numzgrid
+                if (ltopo == 1) then
+                    heightnn (ix,jy,i) = outheight(i) + oro(ix,jy)
+                else
+                    heightnn (ix,jy,i) = outheight(i)
+                endif
+    150 END DO
+
+    if (verbosity.gt.0) then
+    print*,'C Determine area of each output grid box'
+    endif
+
+! Determine area of each output grid box
+!***************************************
+
+    do 140 jy=1,numygrid
+        ylata=outlat0+(float(jy-1)+0.5)*dyout
+        ylatp=ylata+0.5*dyout
+        ylatm=ylata-0.5*dyout
+        if ((ylatm < 0) .AND. (ylatp > 0.)) then
+            hzone=dyout*r_earth*pih
+        else
+            cosfact=cos(ylata*pih)*r_earth
+            cosfactp=cos(ylatp*pih)*r_earth
+            cosfactm=cos(ylatm*pih)*r_earth
+            if (cosfactp < cosfactm) then
+                hzone=sqrt(r_earth**2-cosfactp**2)- &
+                sqrt(r_earth**2-cosfactm**2)
+            else
+                hzone=sqrt(r_earth**2-cosfactm**2)- &
+                sqrt(r_earth**2-cosfactp**2)
+            endif
+        endif
+        gridarea=2.*pi*r_earth*hzone*dxout/360.
+        do 140 ix=1,numxgrid
+            area(ix,jy)=gridarea
+    140 END DO
+    if (verbosity.gt.0) then  
+    print*,'end subroutine readheader10f'
+    endif
+    end subroutine readheader10f

tests/flex_read_fortran/readheader_ncf.f

+      subroutine readheader(filename,nxmax,numxgrid,nymax,numygrid,
+     +nzmax,numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,
+     +ibtime,loutstep,maxspec,nspec,maxageclass,nageclass,lage,
+     +ireleasestart,ireleaseend,maxpoint,numpoint,xpoint,ypoint,
+     +zpoint1,zpoint2,heightnn,area)
+
+      use netcdf
+
+      parameter(pi=3.14159265,r_earth=6.371e6,pih=pi/180.)
+
+      character filename*150,compoint(maxpoint)*45,species(maxspec)*7
+      real outheight(nzmax),heightnn(nxmax,nymax,0:nzmax)
+      integer ireleasestart(maxpoint),ireleaseend(maxpoint)
+      integer npart(maxpoint),kind(maxpoint),lage(0:maxageclass)
+      real xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+      real zpoint2(maxpoint),xmass(maxpoint,maxspec)
+      real oro(nxmax,nymax),area(nxmax,nymax)
+      logical debug
+
+      debug=.FALSE.
+
+
+      print*, 'filename=', filename
+
+      !stop
+
+      open(10,file=filename,form='unformatted',status='old')
+      read(10) ibdate,ibtime
+
+      if (debug) then 
+      write(*,*) ibdate,ibtime
+      endif
+
+      read(10) loutstep,loutaver,loutsample
+      read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+
+      if (debug) then
+      write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+      endif
+
+      read(10) numzgrid,(outheight(i),i=1,numzgrid)
+      read(10) jjjjmmdd,ihmmss
+
+      if (debug) then
+      write(*,*) jjjjmmdd,ihmss
+      endif
+      
+      read(10) nspec
+      nspec=nspec/3
+      do 8 n=1,nspec
+        read(10) numzgrid,species(n)
+        read(10) numzgrid,species(n)
+8       read(10) numzgrid,species(n)
+
+
+      read(10) numpoint
+      do 13 i=1,numpoint
+        read(10) ireleasestart(i),ireleaseend(i)
+
+        read(10) xpoint(i),ypoint(i),xp2,yp2,zpoint1(i),zpoint2(i)
+        read(10) npart(i),kind(i)
+        read(10) compoint(i)
+        do 13 j=1,nspec
+          read(10)
+          read(10)
+13        read(10) xmass(i,j)
+      read(10) method
+      read(10) nageclass,(lage(i),i=1,nageclass)
+      lage(0)=0
+
+      if (debug) then
+      write(*,*) (lage(i),i=0,nageclass)
+      endif
+
+      do 130 ix=1,numxgrid
+130     read(10) (oro(ix,jy),jy=1,numygrid)
+
+      close(10)
+      if (debug) then
+      write(*,*) (outheight(i),i=1,numzgrid)
+      endif
+      
+      if (loutstep.lt.0) nspec=numpoint
+
+
+c Calculate height, which is outheight plus topography
+******************************************************
+
+      do 150 ix=1,numxgrid
+        do 150 jy=1,numygrid
+          if (ltopo.eq.1) then
+            heightnn (ix,jy,0) = oro(ix,jy)
+          else
+            heightnn (ix,jy,0) = 0.
+          endif
+          do 150 i=1,numzgrid
+            if (ltopo.eq.1) then
+              heightnn (ix,jy,i) = outheight(i) + oro(ix,jy)
+            else
+              heightnn (ix,jy,i) = outheight(i)
+            endif
+150         continue
+
+
+C Determine area of each output grid box
+****************************************
+
+      do 140 jy=1,numygrid
+        ylata=outlat0+(float(jy-1)+0.5)*dyout
+        ylatp=ylata+0.5*dyout
+        ylatm=ylata-0.5*dyout
+        if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
+          hzone=dyout*r_earth*pih
+        else
+          cosfact=cos(ylata*pih)*r_earth
+          cosfactp=cos(ylatp*pih)*r_earth
+          cosfactm=cos(ylatm*pih)*r_earth
+          if (cosfactp.lt.cosfactm) then
+            hzone=sqrt(r_earth**2-cosfactp**2)-
+     +      sqrt(r_earth**2-cosfactm**2)
+          else
+            hzone=sqrt(r_earth**2-cosfactm**2)-
+     +      sqrt(r_earth**2-cosfactp**2)
+          endif
+        endif
+        gridarea=2.*pi*r_earth*hzone*dxout/360.
+        do 140 ix=1,numxgrid
+140       area(ix,jy)=gridarea
+
+
+      end

tests/flex_read_fortran/readheader_stdout90.f90

+    program flexpartread
+
+    parameter(nxmax=720,nymax=360,nzmax=40,maxtime=10)
+    parameter(maxpoint=1,maxspec=1,maxageclass=1)
+
+    double precision :: jul,julstart,juldate
+    real :: outheight(nzmax)
+    integer :: ireleasestart(maxpoint),ireleaseend(maxpoint)
+    integer :: lage(0:maxageclass)
+    real :: xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
+    real :: zpoint2(maxpoint)
+    integer :: ireldate(maxspec),ireltime(maxspec)
+
+    real :: grid(0:nxmax-1,0:nymax-1,nzmax,maxspec,maxageclass)
+    real :: footprint(0:nxmax-1,0:nymax-1)
+    real :: footprint_total(0:nxmax-1,0:nymax-1,maxspec)
+    real :: area(0:nxmax-1,0:nymax-1)
+    real :: heightnn(0:nxmax-1,0:nymax-1,0:nzmax)
+    real :: smallnum,r_sparse_dump(5000000)
+    integer :: i_sparse_dump(5000000)
+    parameter(smallnum=1.e-38)
+
+    character dirname*120
+    logical ::  dirnameisarg
+    character filename*150,fileout*150,aday*8,atime*6
+    character(250) :: filegrid
+          
+    logical :: sp_zer
+
+    real :: sizegrid_m
+    real :: sum_grid, sum_sum_grid
+    real :: max_grid, max_max_grid
+    character(4) :: grid_magnitude
+
+    grid_magnitude='time'
+    grid_magnitude='conc'
+    grid_magnitude='pptv'
+
+    dirnameisarg= .FALSE. 
+    print*, 'iargc=', iargc()
+    select case (iargc())
+    case (2)
+    dirnameisarg= .TRUE. 
+    call getarg(1,dirname)
+    len=index(dirname,' ')-1
+    call getarg(2,grid_magnitude)
+    case (1)
+    call getarg(1,dirname)
+    len=index(dirname,' ')-1
+    !call getarg(1,grid_magnitude)
+    case (0)
+    grid_magnitude='time'
+    end select
+    print*, 'grid_magnitude=', grid_magnitude
+
+ !   if ( .NOT. dirnameisarg) then
+ !   !******************
+ !       print*, " Read CONTROL file "
+    !******************
+ !       open(15,file='dirlist')
+ !       200 read(15,'(a)',end=199) dirname
+ !       len=index(dirname,' ')-1
+ !   endif
+
+!      do 50 mother_or_nest=1,2
+    do 50 mother_or_nest=1,1
+
+    !*********************************
+    ! Read FLEXPART header information
+    !*********************************
+
+        if (mother_or_nest == 1) then
+            filename=dirname(1:len)//'header'
+        else if (mother_or_nest == 2) then
+            filename=dirname(1:len)//'header_nest'
+        endif
+        print*, 'call readheader10f for', filename
+    !        call readheader(filename,nxmax,numxgrid,nymax,numygrid,nzmax,
+        call readheader10f(filename,nxmax,numxgrid,nymax,numygrid,nzmax, &
+        numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime, &
+        loutstep,maxspec,nspec,maxageclass,nageclass,lage,ireleasestart, &
+        ireleaseend,maxpoint,numpoint,xpoint,ypoint,zpoint1,zpoint2, &
+        heightnn,area)
+        print*, 'print readheader10f output'
+        print*,'! Write the header information'
+                
+        print*,'! Write info on output interval, averaging time, samplin &
+        g time'
+        print*,'! Write information on output grid setup'
+        print*,'! Write number of species, and name for each species '
+        print*,'dimension of the fields (1 for depo, numzgrid for conc)'
+        print*,'! Write information on release points:'
+        print*,'! Write information on some model switches'
+        print*,'! Write age class information'
+        print*,'! Write topography to output file'
+
+                
+        print*,'filename'
+        !print*,
+
+        print*,'nxmax,numxgrid,nymax,numygrid='
+        print*,filename,nxmax,numxgrid,nymax,numygrid
+        
+        print*,'nzmax, numzgrid='
+        print*,nzmax,numzgrid
+
+        print*,'outlon0,outlat0,dxout,dyout,' 
+        print*,outlon0,outlat0,dxout,dyout 
+        
+        print*,' outheight='
+        do i=1,numzgrid
+        print*, i, outheight(i)
+        end do
+
+        print*, 'ibdate,ibtime, loutstep'
+        print*, ibdate,ibtime, loutstep
+
+!  ,,,ireleasestart, &
+!        ireleaseend,maxpoint,numpoint,xpoint,ypoint,zpoint1,zpoint2, &
+!       heightnn,area)
+
+
+        print*,'maxspec,nspec'
+        print*,maxspec,nspec
+    ! rint*,
+        print*,'maxageclass,nageclass'
+        print*, maxageclass,nageclass
+    ! rint*,
+        print*,'lage'
+        print*, lage 
+
+        print*,'ireleasestart,   ireleaseend,'
+        print*, ireleasestart,   ireleaseend 
+
+        print*,'maxpoint,numpoint'
+        print*, maxpoint,numpoint 
+
+        print*,'xpoint,ypoint'
+        print*, xpoint,ypoint 
+
+        print*,'zpoint1,zpoint2'
+        print*, zpoint1,zpoint2  
+
+        print*,'heightnn'
+        print*, size(heightnn) 
+
+!        print*,'xpoint,ypoint,zpoint1,zpoint2,'
+!        print*, xpoint,ypoint,zpoint1,zpoint2 
+
+    ! rint*,
+
+        write(*,*) 'nspec',nspec
+        julstart=juldate(ibdate,ibtime)
+        do 30 k=1,nspec
+            print*, 'k', k 
+            jul=julstart+dble(float(ireleaseend(k))/86400.)
+            call caldate(jul,ireldate(k),ireltime(k))
+            print*,'ireldate',ireldate(k)
+            print*,'ireltime',ireltime(k)
+        30 END DO
+
+
+         
+    !      print*, 'about to exit'
+    !      exit
+
+    50 END DO
+
+ !   if ( .NOT. dirnameisarg) then
+ !       goto 200 !read(15, dirlist
+ !       199 close(15) !dirlist
+ !   endif
+
+! ax_max_grid=max(max_max_grid,max_grid)
+! um_sum_grid=sum_sum_grid+sum_grid
+
+!         print*, 'for all times'
+!         print*, 'sum sum(grid)=', sum_sum_grid
+!         print*, 'max maxval(grid)=', max_max_grid
+
+!         print*, 'sum_sum_grid_out=',sum_sum_grid
+!         print*, 'max_max_grid_out=',max_max_grid
+    END PROGRAM

tests/flex_read_fortran/readpart.f

+       SUBROUTINE readpart(flexout,datetime,numpart)
+       use work_arrays
+       implicit none
+c      version 2 of 10 columns of reals
+c      include '/home/ignacio/flexpart/FLEXPART8.1/includepar'
+c      include '/home/ignacio/flexpart/FLEXPART8.1/includecom'
+      DOUBLE PRECISION jul
+!c      integer itime,i,j,jjjjmmdd,ihmmss itime in work_arrays
+      integer i,j,jjjjmmdd,ihmmss
+      integer ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+      real xlon,ylat
+      real dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
+      real topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
+      real tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
+      real tr(2),tri
+      character datetime*14 , flexout*150, filename*150     !atime*6,path*90,directory*90
+      integer unitpartout
+      integer iargc
+      integer npoint
+      real ztra1
+      integer itramem(1)
+c      real topo,pvi,qvi,rhoi,hmixi,tri,tti,
+      real xmass(1,1) 
+      !integer nspec,p
+      integer p
+      integer numpart
+      integer ios 
+      logical verbose
+
+      verbose=.false. ! .true.
+      unitpartout=1
+      nspec=1
+
+      filename=trim(flexout)//'partposit_'//datetime
+
+      !print*,filename
+
+      if ( datetime(1:1) == '0' ) then
+        filename=trim(flexout)
+        !print*,filename
+        !stop 42
+      endif
+
+!      i=0
+!      open(unitpartout,file=trim(flexout)//
+!c     + 'part2_'//datetime,
+!     + 'partposit_'//datetime,
+!     + form='unformatted',IOSTAT=ios)
+!
+      open(unitpartout,file=trim(filename),
+     + form='unformatted',IOSTAT=ios)
+
+     
+      read(unitpartout) itime
+      if ( verbose ) then
+        print *, 'readpart2> itime :' , itime
+      endif
+      if ( numpart==-1 ) then
+        do while ( ios == 0 )
+!          read(unitpartout, IOSTAT=ios)  xlon,ylat,ztra1,tri,hmixi
+          read(unitpartout, IOSTAT=ios)  npoint, xlon,ylat,ztra1, !,tri,hmixi
+     +    itramem,topo,pvi,qvi,rhoi,hmixi,tri,tti, xmass
+
+        i=i+1
+          if ( verbose ) then
+            print *, 'readpart2> i,ios',  i,ios,npoint, xlon,ylat,ztra1 
+          endif
+!c      print *, 'ios', ios
+!c          read(unitpartout) npoint,xlon,ylat,ztra1,
+!c     +    itramem(i),topo,pvi,qvi,rhoi,hmixi,tri,tti,
+!c     +    (xmass(i,j),j=1,nspec)
+!npoint,xlon,ylat,ztra1,
+!c     +    itramem(i),topo,pvi,qvi,rhoi,hmixi,tri,tti,
+!c     +    xmass
+        end do
+!          numpart=i-1
+          numpart=i-2
+      else !numpart~=-1
+        do i=1,numpart
+!          read(unitpartout, IOSTAT=ios)  xlon,ylat,ztra1,tri,hmixi
+!     +    ,topo,pvi,tti,qvi,rhoi
+        read(unitpartout, IOSTAT=ios) vnpoint(i), 
+     +    vxlon(i),vylat(i), vztra1(i),
+     +    vitramem(i),
+     +    vtopo(i),vpvi(i),vqvi(i),vrhoi(i),vhmixi(i), 
+     +    vtri(i), vtti(i),
+     +    (xmass0(i,j),j=1,nspec)
+!        print *, 'readpart2> vztra1(i)',vztra1(i) ,i        
+!c       vlat(i)=ylat
+!c       interpolate OH from the field (needs month, x, y ,z )
+!c        OH=
+!c          print *, i,   xlon,ylat,ztra1,tri,hmixi
+!c     + ,topo,pvi,tti,qvi,rhoi
+        end do
+      endif ! numpart==-1 
+      close(unitpartout)
+      end subroutine

tests/flex_read_fortran/work_arrays.f90

+!arrays changing with time, because of reading or transformation
+   MODULE work_arrays
+      INTEGER itime
+      INTEGER ::  nspec=1
+      REAL   , DIMENSION(:,:),ALLOCATABLE :: xmass0(:,:), xmass1(:,:) , xmass1_injected(:,:), xmass1_remaining(:,:)
+      REAL   , DIMENSION(:,:),ALLOCATABLE :: xmass1_SG,xmass1_PG
+      REAL   , DIMENSION(:,:),ALLOCATABLE :: xmass2 !mass of artificial tracer 
+      integer, DIMENSION(:),  ALLOCATABLE :: vnpoint(:),vitramem(:)
+      REAL   , DIMENSION(:),  ALLOCATABLE :: vxlon(:),vylat(:),vztra1(:),vtri(:),vhmixi(:),vtopo(:),vpvi(:),vtti(:),vqvi(:),vrhoi(:)
+      REAL   , DIMENSION(:),  ALLOCATABLE ::  vxlon2,vxlon3
+      INTEGER, DIMENSION(:),  ALLOCATABLE :: i_reach_strat(:), i_reach_400(:), i_reach_2PVU(:)
+      INTEGER, DIMENSION(:),  ALLOCATABLE :: reach_strat(:), reach_400(:), reach_2PVU(:), flag_reach_strat(:), flag_over_400(:)
+!      INTEGER, DIMENSION(:),  ALLOCATABLE :: ireach_strat(:), ireach_400(:), ireach_2PVU(:)
+      INTEGER, DIMENSION(:),  ALLOCATABLE :: secs2_400(:), secs2_380(:), secs2_360(:), secs2_340(:)
+      INTEGER, DIMENSION(:),  ALLOCATABLE ::  secs2_conventional(:), secs2_trop_FLEX(:),secs2_2PVU380(:),secs2_35PVU380(:)
+      INTEGER ,DIMENSION(:),  allocatable :: ETOW(:),TOW(:),LMS(:),TTL(:),UT(:),STRAT(:),OW(:),MW(:)     
+      LOGICAL, DIMENSION(:),  ALLOCATABLE :: l_reached_strat(:),l_reaching_strat(:),l_still_in_trop
+      REAL   , DIMENSION(:),  ALLOCATABLE :: vtheta,vpres
+      REAL   , DIMENSION(:),  ALLOCATABLE :: lat_in_strat,lat_went_strat
+      REAL   , DIMENSION(:,:),ALLOCATABLE :: out_matrix 
+      INTEGER, DIMENSION(:,:,:),  ALLOCATABLE :: this_source
+      REAL   , DIMENSION(:),  ALLOCATABLE :: lon_down, lat_down
+      INTEGER  :: nlon_down, nlat_down
+      REAL   , DIMENSION(:,:),ALLOCATABLE ::   chi_30d_strat
+      INTEGER, DIMENSION(:,:),  ALLOCATABLE :: N_1, N_0
+
+     ! nspec=1
+
+      
+   END MODULE
+

tests/read_examples/config.10.4.NILU

+
+module purge 
+module load eccodes/2.28.0-gcc9.4 
+export FLEXPART=/home/ignacio/repos/flexpart/src/FLEXPART 
+#export  printgrid=../flex_read_fortran/printgrid
+
+

tests/read_examples/declare_examples

+#!/bin/bash
+
+
+#group 1 (6+2*1=8) IOUT=3 double output
+#declare -a  examples1=('1-1_1' '1-2_2' '1-3_3' '1-4_3')
+#declare -a  units1=('conc' 'pptv' 'conc' 'pptv')
+#declare -a  examples1=('1-1' '1-2' '1-3' '1-3' '1-5')
+#declare -a  examples1=('1-1_1' '1-2_2' '1-3_3' '1-3_3' '1-5_5')
+#declare -a  examples1=('1-1_1' '1-2_2' '1-3_3' '1-3_3' '1-5_5' '1-5_9' '1-5_10' '1-5_11')
+#declare -a  units1=('conc' 'pptv' 'conc' 'pptv' 'conc')
+declare -a examples1=('1-1_1' '1-2_2' '1-3_3' '1-3_3' '1-5_5' '1-6_9' '1-7_10' '1-8_11')
+declare -a    units1=('conc'  'pptv'  'conc'  'pptv'  'conc'   'conc'     'pptv'      'pptv')
+declare -a   format1=('cm'    'cm'    'cm'    'cm'    'cm'     'nc'    'nc'     'nc')
+declare -a   filens1=('cm'    'cm'    'cm'    'cm'    'cm'     'grid_conc_20120101060000.nc'    'grid_conc_20120101060000.nc'     'grid_conc_20120101060000.nc')
+
+#group 2
+#declare -a  examples2=('2-1_bwd' '2-2_bwd5')
+#declare -a  units2=('time' 'time')
+#group 2 (3)
+declare -a  examples2=('2-1_bwd' '2-2_bwd5' '2-3_bwd_nc')
+declare -a  units2=('time' 'time' 'time')
+declare -a  format2=('cm' 'cm' 'nc')
+declare -a  filens2=('cm' 'cm' 'grid_time_20120101120000.nc')
+
+
+#group 3 (3)
+declare -a  examples3=('3-part1' '3-part2' '3-part_bwd1')
+#declare -a  examples3=('3-1_part' '3-2_part_end' '3-3_part_bwd')
+#declare -a  units3=('part' 'part' 'part')
+declare -a  units3=('20120101100000' 'end' '20120101080000')
+declare -a  format3=('part' 'part' 'part')
+declare -a  filens3=('part' 'part' 'part')
+
+
+#group 4 (6)
+declare -a examples4=('4-1_ind_1_2' '4-2_ind_2_1' '4-3_ind_2_2' '4-4_bwd_ind_1_2' '4-5_bwd_ind_2_1' '4-6_bwd_ind_2_2' )
+declare -a    units4=('conc'         'conc'        'conc'        'time'           'time'            'time')
+declare -a   format4=('cm' 'cm' 'cm' 'cm' 'cm' 'cm')
+declare -a   filens4=('cm' 'cm' 'cm' 'cm' 'cm' 'cm')
+
+
+#group 5 (4)
+#declare -a  examples5=('5-1_specNO' '5-2_specAERO-TRACE' '5-3_specCO' '5-4_specBC' '5-5_bwd_specNO' )
+#declare -a     units5=('conc'       'conc'               'conc'       'conc'       'time')
+declare -a  examples5=('5-1_specNO' '5-2_specCO' '5-3_specAERO-TRACE' '5-4_specBC' )
+declare -a  units5=('conc' 'conc' 'conc'  'conc')
+declare -a    format5=('cm' 'cm' 'cm' 'cm')
+declare -a    filens5=('cm' 'cm' 'cm' 'cm')
+
+#group 6 (6)
+declare -a  examples6=('6-1_nested' '6-2_nested_mr'  '6-3_nested_bwd' '6-4_nested_nc' '6-5_nested_mr_nc' '6-6_nested_bwd_nc')
+declare -a  units6=('conc' 'pptv'  'time'  'conc' 'pptv' 'time')
+#declare -a  examples6=('6-1_nested' '6-2_nested_bwd' '6-3_nested1' '6-4_nested1_bwd' )
+#declare -a     units6=('conc'        'time'           'conc'        'time')
+declare -a    format6=('cm' 'cm' 'cm' 'nc' 'nc' 'nc')
+declare -a    filens6=('cm' 'cm' 'cm' 'grid_conc_20120101060000_nest.nc' 'grid_conc_20120101060000.nc' 'grid_time_20120101120000_nest.nc')
+
+#group 7 (1)
+declare -a  examples7=('7-1_DOMAINFILL')
+declare -a  units7=('conc')
+declare -a  format7=('cm')
+declare -a  filens7=('cm')
+
+
+#group 8 (3)
+#options_8-1_init_cond/
+#options_8-2_init_cond_surf/
+#options_8-3_init_cond_ind_2_1/
+#options_8-4_init_cond_ind_1_2/
+#options_8-5_init_cond_ind_2_2/
+#declare -a  examples8=('8-1_init_cond' '8-2_init_cond_surf' '8-3_init_cond_ind_2_1' '8-4_init_cond_ind_1_2' '8-5_init_cond_ind_2_2' )
+#declare -a  units8=('time' 'time' 'time'  'time' 'time')
+declare -a  examples8=('8-1_init_cond' '8-2_init_cond_surf' '8-3_init_cond_ind_2_1' )
+#declare -a  examples8=('8-1_init_cond' '8-2_init_cond_ind_1_2'  '8-3_init_cond_surf' )
+declare -a  units8=('time' 'time' 'time' )
+declare -a  format8=('cm' 'cm' 'cm')
+declare -a  filens8=('cm' 'cm' 'cm')
+
+#suffix=_8-1_init_cond
+#suffix=_8-2_init_cond_surf
+#suffix=_8-3_init_cond_ind_2_1
+
+#group 9 (2)
+#output_9-1_CBLFLAG/
+#output_9-2_CBLFLAG_bwd/
+declare -a  examples9=('9-1_CBLFLAG' '9-2_CBLFLAG_bwd')
+declare -a  units9=('conc' 'time')
+declare -a  format9=('cm' 'cm')
+declare -a  filens9=('cm' 'cm')
+
+#group 10 (2)
+declare -a  examples10=('10-1_bwd_ind_1_3' '10-2_bwd_ind_1_4')
+declare -a  units10=('wetdep' 'drydep')
+declare -a  filens10=('cm' 'cm')
+declare -a  format10=('cm' 'cm')
+#declare -a  format10=('cmd' 'cmd')
+
+#group 11 (1)
+declare -a  examples11=('11-1_Volc')
+declare -a  units11=('conc')
+declare -a  format11=('cm')
+declare -a  filens11=('cm')
+
+
+
+#examples=("${examples1[@]}" "${examples2[@]}" "${examples4[@]}" "${examples5[@]}" "${examples6[@]}" "${examples7[@]}" "${examples8[@]}" "${examples9[@]}")
+#   units=("${units1[@]}"    "${units2[@]}"    "${units4[@]}"    "${units5[@]}"    "${units6[@]}"    "${units7[@]}"    "${units8[@]}"    "${units9[@]}" )
+
+

tests/read_examples/read_examples.sh

+#!/bin/bash
+
+# export printgrid=$HOME/repos/flexpart/postproc/flex_read_fortran/printgrid
+
+# default path to examples
+#examples_default=./
+#examples_default=../../examples
+examples_default=../examples
+# environment override on default
+#examples_int=${1:-$examples_default}
+examples_env=${examples:-$examples_default}
+#internal (argument override)
+examples_int=${1:-$examples_env}
+
+
+#printgrid_default=printgrid
+#printgrid_default=/Users/ignacio/repos/flexpart/postproc/flex_read_fortran/printgrid
+printgrid_default=../flex_read_fortran/printgrid
+# todo test exixtence of executable
+# alternative: printgrid_default=$HOME/repos/flexpart/postprocess/read_fortran/printgrid
+#overrun by input from ENV
+printgrid_env=${printgrid:-$printgrid_default}
+# overrun by command line arguments 
+printgrid_int=${2:-$printgrid_env}
+
+#printpart_default=../../postprocess/flex_read_fortran_part/printpart
+printpart_default=../flex_read_fortran/printpart
+#overrun by input from ENV
+printpart_env=${printpart:-$printpart_default}
+# overrun by command line arguments 
+printpart_int=${2:-$printpart_env}
+
+
+
+echo default example
+stdout=$($printgrid_int $examples_int/output/ conc)
+echo default $stdout
+
+set_examples_default=set_examples_1.sh
+echo default example variations $set_examples_default
+
+set_examples_env=${set_examples:-$set_examples_default} 
+
+#source ./set_examples_3.sh
+#source ./set_examples_1.sh
+#source ./set_examples_all
+source ./${set_examples_env}
+
+echo 'echo ${examples[@]}' 
+echo ${examples[@]} 
+
+iend=$(( ${#examples[@]} - 1 ))
+#echo $(( iend + 1 )) 
+echo loop on $(( iend + 1  )) examples  
+ii=0
+for i in `seq 0 $iend`;
+do 
+  ii=$(( ii + 1 ))
+  #echo $i
+  #echo $example
+  b=${example:0:1}
+  #echo $b
+#  if [ $i == 3 ]; then  ii=3 ; fi
+  if [ $i == 3 ] && [ $b == 1 ]; then  ii=3 ; fi
+
+  example=${examples[$i]}
+  unit=${units[$i]}
+  format=${formats[$i]}
+#    if [ $unit == "x" ]; then
+#      echo $i '&' "not gridded binary"
+#    else
+    if [ $format == "cm" ]; then
+      #echo $printgrid_int $examples_int/output_${example}/ ${unit}
+      #stdout=$($printgrid_int $examples_int/output_1-1/ conc)
+      stdout=$($printgrid_int $examples_int/output_${example}/ ${unit})
+      #echo $i '&' ${example}: $stdout
+      echo $ii '&' ${example} '&' $stdout '\\'
+    elif [ $format == "part" ]; then
+      #echo $printpart_int $examples_int/output_${example}/ ${unit}
+      stdout=$($printpart_int $examples_int/output_${example}/ ${unit})
+      echo $ii '&' ${example} '& xyz:' $stdout '\\'     
+    elif [ $format == "nc" ]; then
+      grid_conc_default="grid_conc_20190225190000.nc"
+      grid_conc_env=${grid_conc:-$grid_conc_default} 
+      #echo ${printgrid_int}_ncf $examples_int/output_${example}/ ${unit} $grid_conc_env
+      filen=${filens[$i]}
+      #echo $filen
+      stdout=$(${printgrid_int}_ncf $examples_int/output_${example}/ ${unit} $filen)
+      #stdout=$(${printgrid_int}_ncf $examples_int/output_${example}/ ${unit} $grid_conc_env)
+      #echo $stdout
+      echo $ii '&' ${example} '&' $stdout '\\'       
+    else
+      echo $ii '&' ${example} '&' $format "not yet supported"
+    fi
+done
+
+exit
+$printgrid_int $examples_int/output_2/ pptv
+
+$printgrid_int $examples_int/output_bwd/ time
+
+
+exit 

tests/read_examples/set_examples_all

+#!/bin/bash
+# display in 83 columns
+source ./declare_examples
+
+examples=("${examples1[@]}" "${examples2[@]}" "${examples3[@]}"  "${examples4[@]}"           "${examples5[@]}" "${examples6[@]}" "${examples7[@]}" "${examples8[@]}"            "${examples9[@]}" "${examples10[@]}" "${examples11[@]}")                 
+units=( "${units1[@]}" "${units2[@]}" "${units3[@]}" "${units4[@]}" "${units5[@]}"         "${units6[@]}" "${units7[@]}" "${units8[@]}" "${units9[@]}" "${units10[@]}"       "${units11[@]}")
+formats=("${format1[@]}" "${format2[@]}" "${format3[@]}" "${format4[@]}"                    "${format5[@]}" "${format6[@]}" "${format7[@]}" "${format8[@]}"                    "${format9[@]}" "${format10[@]}" "${format11[@]}")
+
+filens=("${filens1[@]}" "${filens2[@]}" "${filens3[@]}" "${filens4[@]}"                    "${filens5[@]}" "${filens6[@]}" "${filens7[@]}" "${filens8[@]}"                    "${filens9[@]}" "${filens10[@]}" "${filens11[@]}")
+
+# source ./display_examples.sh
+
+#moved to 
+# ./display_examples.sh
+#iend=$(( ${#examples[@]} - 1 ))
+#echo 'iend + 1 = ' $(( iend + 1 )) 
+#for i in `seq 0 $iend`;
+#do 
+#  #echo $i
+#  example=${examples[$i]}
+#  unit=${units[$i]}
+#  format=${formats[$i]}
+
+#  echo $example $unit $format
+#    if [ $unit == "x" ]; then
+#      echo $i '&' "not gridded binary"
+#    else
+#      #echo $printgrid_int $examples_int/output_${example}/ ${unit}
+      #stdout=$($printgrid_int $examples_int/output_1-1/ conc)
+#      stdout=$($printgrid_int $examples_int/output_${example}/ ${unit})
+      #echo $i '&' ${example}: $stdout
+#      echo $i '&' ${example} '&' $stdout '\\'
+#    fi
+#done
+
+#exit