diff --git a/tests/flex_read_fortran/README b/tests/flex_read_fortran/README
new file mode 100644
index 0000000000000000000000000000000000000000..189e68701d3354e0cfef58da361ef7c2ccfd64f2
--- /dev/null
+++ b/tests/flex_read_fortran/README
@@ -0,0 +1,12 @@
+
+configure the path to netcdf in makefile
+and build
+
+make
+
+equivalent to
+
+make printgrid
+make printgrid_netcdf
+make printpart
+
diff --git a/tests/flex_read_fortran/caldate.f b/tests/flex_read_fortran/caldate.f
new file mode 100644
index 0000000000000000000000000000000000000000..07ce4f2fcebd0ec38d4124e5b2349dfa4e3d725b
--- /dev/null
+++ b/tests/flex_read_fortran/caldate.f
@@ -0,0 +1,68 @@
+ 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
diff --git a/tests/flex_read_fortran/flex_read_compare2.f90 b/tests/flex_read_fortran/flex_read_compare2.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6d7961c4ee5a1aa31fb47dcf8c07a72f45a11847
--- /dev/null
+++ b/tests/flex_read_fortran/flex_read_compare2.f90
@@ -0,0 +1,248 @@
+ 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
diff --git a/tests/flex_read_fortran/juldate.f b/tests/flex_read_fortran/juldate.f
new file mode 100644
index 0000000000000000000000000000000000000000..90699acc393597d96a314ef8bcf08c3211ced179
--- /dev/null
+++ b/tests/flex_read_fortran/juldate.f
@@ -0,0 +1,58 @@
+ 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
diff --git a/tests/flex_read_fortran/makefile b/tests/flex_read_fortran/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..b030f249917273900d5aecbe59e338b3cdf6cd06
--- /dev/null
+++ b/tests/flex_read_fortran/makefile
@@ -0,0 +1,46 @@
+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
diff --git a/tests/flex_read_fortran/printgrid.f90 b/tests/flex_read_fortran/printgrid.f90
new file mode 100644
index 0000000000000000000000000000000000000000..60e2f0e2852a6d74cc6f5b1128a9677191dc5ab3
--- /dev/null
+++ b/tests/flex_read_fortran/printgrid.f90
@@ -0,0 +1,355 @@
+ 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
diff --git a/tests/flex_read_fortran/printgrid_ncf.f90 b/tests/flex_read_fortran/printgrid_ncf.f90
new file mode 100644
index 0000000000000000000000000000000000000000..7b0d447530d8fcbedfe0e40b37565d55f27caada
--- /dev/null
+++ b/tests/flex_read_fortran/printgrid_ncf.f90
@@ -0,0 +1,466 @@
+ 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
+
diff --git a/tests/flex_read_fortran/printpart.f90 b/tests/flex_read_fortran/printpart.f90
new file mode 100644
index 0000000000000000000000000000000000000000..342637001fb83fb79205de04f69c831be4d93e5f
--- /dev/null
+++ b/tests/flex_read_fortran/printpart.f90
@@ -0,0 +1,83 @@
+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
diff --git a/tests/flex_read_fortran/readgrid.f b/tests/flex_read_fortran/readgrid.f
new file mode 100644
index 0000000000000000000000000000000000000000..90812aab6a4cb7db8e2984c6ea2f8ddb1544979f
--- /dev/null
+++ b/tests/flex_read_fortran/readgrid.f
@@ -0,0 +1,133 @@
+ 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
diff --git a/tests/flex_read_fortran/readgrid_ncf.f b/tests/flex_read_fortran/readgrid_ncf.f
new file mode 100644
index 0000000000000000000000000000000000000000..90812aab6a4cb7db8e2984c6ea2f8ddb1544979f
--- /dev/null
+++ b/tests/flex_read_fortran/readgrid_ncf.f
@@ -0,0 +1,133 @@
+ 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
diff --git a/tests/flex_read_fortran/readheader.f b/tests/flex_read_fortran/readheader.f
new file mode 100644
index 0000000000000000000000000000000000000000..0c04759a836ef7e1dc0e9f65a95295ace23fbc7e
--- /dev/null
+++ b/tests/flex_read_fortran/readheader.f
@@ -0,0 +1,124 @@
+ 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
diff --git a/tests/flex_read_fortran/readheader1090.f90 b/tests/flex_read_fortran/readheader1090.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c75646f6a579e1e4144b3f1bbfaba03bfc9d2f35
--- /dev/null
+++ b/tests/flex_read_fortran/readheader1090.f90
@@ -0,0 +1,253 @@
+ 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
diff --git a/tests/flex_read_fortran/readheader_ncf.f b/tests/flex_read_fortran/readheader_ncf.f
new file mode 100644
index 0000000000000000000000000000000000000000..5e0c32a339b34d935de70c684c2e8cee584f4ad8
--- /dev/null
+++ b/tests/flex_read_fortran/readheader_ncf.f
@@ -0,0 +1,131 @@
+ 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
diff --git a/tests/flex_read_fortran/readheader_stdout90.f90 b/tests/flex_read_fortran/readheader_stdout90.f90
new file mode 100644
index 0000000000000000000000000000000000000000..e3b8316ef853e373272e767a0c3770f16e366b16
--- /dev/null
+++ b/tests/flex_read_fortran/readheader_stdout90.f90
@@ -0,0 +1,183 @@
+ 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
diff --git a/tests/flex_read_fortran/readpart.f b/tests/flex_read_fortran/readpart.f
new file mode 100644
index 0000000000000000000000000000000000000000..4704da8ab6842e244958b42f52ae3d7955f06a32
--- /dev/null
+++ b/tests/flex_read_fortran/readpart.f
@@ -0,0 +1,97 @@
+ 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
diff --git a/tests/flex_read_fortran/work_arrays.f90 b/tests/flex_read_fortran/work_arrays.f90
new file mode 100644
index 0000000000000000000000000000000000000000..83e70f9a36fd18988c57798d7e281f27084630a1
--- /dev/null
+++ b/tests/flex_read_fortran/work_arrays.f90
@@ -0,0 +1,31 @@
+!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
+
diff --git a/tests/read_examples/config.10.4.NILU b/tests/read_examples/config.10.4.NILU
new file mode 100644
index 0000000000000000000000000000000000000000..b864cd59558f9d145b4b10493c42b6a25f69b22e
--- /dev/null
+++ b/tests/read_examples/config.10.4.NILU
@@ -0,0 +1,7 @@
+
+module purge
+module load eccodes/2.28.0-gcc9.4
+export FLEXPART=/home/ignacio/repos/flexpart/src/FLEXPART
+#export printgrid=../flex_read_fortran/printgrid
+
+
diff --git a/tests/read_examples/declare_examples b/tests/read_examples/declare_examples
new file mode 100644
index 0000000000000000000000000000000000000000..39433eddbd1ec23936ed974b85b94de5e7a13549
--- /dev/null
+++ b/tests/read_examples/declare_examples
@@ -0,0 +1,109 @@
+#!/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[@]}" )
+
+
diff --git a/tests/read_examples/read_examples.sh b/tests/read_examples/read_examples.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b26d6f4656a16c335f71c29ac1f158b9c9573d49
--- /dev/null
+++ b/tests/read_examples/read_examples.sh
@@ -0,0 +1,103 @@
+#!/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
diff --git a/tests/read_examples/set_examples_all b/tests/read_examples/set_examples_all
new file mode 100755
index 0000000000000000000000000000000000000000..9d4679c05d204087372ef9703545b1ae45c765ec
--- /dev/null
+++ b/tests/read_examples/set_examples_all
@@ -0,0 +1,36 @@
+#!/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
diff --git a/tests/read_examples/set_examples_i.sh b/tests/read_examples/set_examples_i.sh
new file mode 100755
index 0000000000000000000000000000000000000000..da0d07a44c93d6ae92a059696b03d3214d34f603
--- /dev/null
+++ b/tests/read_examples/set_examples_i.sh
@@ -0,0 +1,42 @@
+#!/bin/bash
+# display in 83 columns
+source ./declare_examples
+
+i=$1
+
+
+examples_value=$(echo "examples"$i)
+#echo $examples_value
+ #printf -v "examplesi" "%s" "$examples_value"
+ #echo $examplesi
+examples="${examples_value}[@]"
+#echo examples 1: $examples
+ #export examples=("${examplesi[@]}")
+export examples=("${!examples}")
+ #echo ${examples[@]}
+ #echo examples 2: $examples[@]
+#echo examples 2: ${examples[@] }
+echo examples: ${examples[@] }
+
+array_name=$(echo "units"$i)
+units="${array_name}[@]"
+export units=("${!units}")
+echo units: ${units[@] }
+
+#export formats=("${formati[@]}")
+
+array_name=$(echo "format"$i)
+target_array="${array_name}[@]"
+export formats=("${!target_array}")
+echo formats: ${formats[@] }
+
+array_name=$(echo "filens"$i)
+target_array="${array_name}[@]"
+export filens=("${!target_array}")
+echo filens : ${filens[@] }
+
+#source ./display_examples.sh
+# output eg:
+#3-1_part 20120101100000 part
+#3-2_part_end end part
+#3-3_part_bwd 20120101080000 part