add new reading routines and wrap programs for header v10 in f77 fixed and f90 free comments

files originally from 16/8/2018

README

+aim: read output from flexpart default and parameter sweep with fortran
+
+history: 
+2018-08-16
+cloned from /Users/ignacio/todo/doing/flexpart/2018-06-03_Sun22/PREPARE_FOOTPRINT
+
+modified read header and read grid from f77 to f90
+
+added
+
+readheader10.f
+ from readheader10.f modified to be more verbose
+readheader_stdout.f
+ wrapper to call the above function and display the read variables
+
+readheader1090.f90
+readheader_stdout90.f90
+ translated into free f90 format from above
+
+
+
+examples from before
+sum_avg
+ compilation:
+	make sum_avg	
+ run:
+	# sum_avg needs dir to be in file
+	./sum_avg conc
+
+flexout_sum_max
+ compilation:
+	make flexout_sum_max	
+ run:
+	./flexout_sum_max /Users/ignacio/repos/flex_cases/output/ conc
+
+
+read and display header
+
+f77
+make readheader_stdout
+./readheader_stdout 
+f90
+make readheader_stdout90
+./readheader_stdout90
+
+read and display grid
+
+
+
+read and display particles
+
+read and verify a fingerprint
+
+particles: some central measure or other moments
+
+
+2018-08-22
+can read (header and grids) group 1 of standard use cases
+
+08.22 22.58 ignacio@Kjeller:~/repos/flex_read_fortran$./flexout_sum_max /Users/ignacio/repos/flex_gen_input/output_2/ pptv
+08.22 22.58 ignacio@Kjeller:~/repos/flex_read_fortran$./flexout_sum_max /Users/ignacio/repos/flex_gen_input/output_3/ conc
+08.22 22.58 ignacio@Kjeller:~/repos/flex_read_fortran$./flexout_sum_max /Users/ignacio/repos/flex_gen_input/output_3/ pptv
+08.22 22.58 ignacio@Kjeller:~/repos/flex_read_fortran$./flexout_sum_max /Users/ignacio/repos/flex_gen_input/output_5/ conc
+ 
+
+2018-08-23
+
+compare with matlab
+
+case 1A _2
+
+./flexout_sum_max /Users/ignacio/repos/flex_cases_reference/output_2/ pptv
+ sum_sum_grid_out=   90265168.0    
+ max_max_grid_out=   28875172.0   
+
+with matlab:
+% diagnostics
+% sumgrid=90265170.9066
+% maxgrid=28875172
+
+CAse 1A _3 
+
+./flexout_sum_max /Users/ignacio/repos/flex_cases_reference/output_3/ pptv
+ sum_sum_grid_out=   90265168.0    
+ max_max_grid_out=   28875172.0   
+
+./flexout_sum_max /Users/ignacio/repos/flex_cases_reference/output_3/ conc
+ sum_sum_grid_out=   101040392.    
+ max_max_grid_out=   32318520.0    
+
+
+./flexout_sum_max /Users/ignacio/repos/flex_cases_reference/output_5/ conc
+ sum_sum_grid_out=   101040392.    
+ max_max_grid_out=   32318520.0 

dirlist

+/Users/ignacio/repos/flex_cases/output/

makefile

@@ -34,5 +34,12 @@ sum_avg:readheader.o readgrid.o juldate.o caldate.o sum_avg.o
 	$(FC) readheader.o readgrid.o juldate.o caldate.o sum_avg.o -o sum_avg
 flexout_sum_max:readheader.o readgrid.o juldate.o caldate.o flexout_sum_max.o 
 	$(FC) readheader.o readgrid.o juldate.o caldate.o flexout_sum_max.o -o flexout_sum_max
+readheader_stdout:readheader10.o juldate.o caldate.o readheader_stdout.o 
+	$(FC) readheader10.o juldate.o caldate.o readheader_stdout.o -o readheader_stdout
+%.o: %.f90
+	+$(FC) -c $<	
+readheader_stdout90:readheader1090.o juldate.o caldate.o readheader_stdout90.o 
+	$(FC) readheader1090.o juldate.o caldate.o readheader_stdout90.o -o readheader_stdout90
+
 clean:
 	        rm *.o

readheader10.f

+      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)
+      !real 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 
+      print*,'subroutine readheader10f'
+      print*,'! Open header output file: ', filename
+      open(10,file=filename,form='unformatted',status='old')
+
+      
+      !l76
+      print*,'ibdate,ibtime fwd. iedate,ietime bwd'
+      read(10) ibdate,ibtime, flexversion
+      write(*,*) 'date,time=',ibdate,ibtime
+      print*, 'flexversion=',flexversion
+
+      print*,'finfo on output interval, averaging time, sampling time'
+      !l85
+      read(10) loutstep,loutaver,loutsample
+      print*,'loutstep, aver, sample=',loutstep,loutaver,loutsample
+
+      print*, 'information on output grid setup'
+      read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+      write(*,*) 'outlon0,outlat0,numxgrid,numygrid,dxout,dyout'
+      write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+      read(10) numzgrid,(outheight(i),i=1,numzgrid)
+      write(*,*) 'numzgrid=', numzgrid
+      write(*,*) 'outheight=', outheight
+
+      read(10) jjjjmmdd,ihmmss
+      write(*,*) 'jjjjmmdd,ihmmss=',jjjjmmdd,ihmss
+
+      print*, 'information on species'
+      !l102  3*nspec,maxpointspec_act implicit? 
+      read(10) nspec, maxpointspec_act
+      print*,'3*nspec,maxpointspec_act=',nspec,maxpointspec_act
+
+      nspec=nspec/3
+      do 8 n=1,nspec
+        read(10) numzgrid,species(n)
+        print*,'numzgrid,species(n)',numzgrid,species(n)  
+        read(10) numzgrid,species(n)
+        print*,'numzgrid,species(n)',numzgrid,species(n)  
+        read(10) numzgrid,species(n)
+8       print*,'numzgrid,species(n)',numzgrid,species(n)  
+
+
+      print*, 'information on release points:'
+      read(10) numpoint !l113
+      print*,'numpoint=',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)
+      print*,'! Write information on some model switches'
+      read(10) method,lsubgrid,lconvection,ind_source,ind_receptor
+      print*,'method,lsubgrid,lconvection,ind_source,ind_receptor'
+      print*,method,lsubgrid,lconvection,ind_source,ind_receptor
+
+      print*,'!Write age class information'
+      read(10) nageclass,(lage(i),i=1,nageclass)
+      lage(0)=0
+      print*,'nageclass=',nageclass
+      write(*,*) 'lage(0:nageclass)=', (lage(i),i=0,nageclass)
+
+      print*,'!Write topography to output file' 
+      do 130 ix=1,numxgrid
+130     read(10) (oro(ix,jy),jy=1,numygrid)
+      print*,'(oro(ix,jy)'
+      print*,'numxgrid,numygrid=',numxgrid,numygrid
+      print*,'nxmax,nymax=',nxmax,nymax
+      
+      close(10)
+
+
+
+      !write(*,*) (outheight(i),i=1,numzgrid)
+
+      !if (loutstep.lt.0) nspec=numpoint
+
+
+      print*,' !if (loutstep.lt.0) nspec=numpoint'
+      print*,'outheight(0)=0. ???'
+ 
+      print*,'c Calculate height, which is outheight plus topography'
+
+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
+
+      print*,'C Determine area of each output grid box'
+
+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
+
+      print*,'end subroutine readheader10f'      
+      end

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
+    print*,'subroutine readheader10f'
+    print*,'! Open header output file: ', filename
+    open(10,file=filename,form='unformatted',status='old')
+
+          
+! 76
+    print*,'ibdate,ibtime fwd. iedate,ietime bwd'
+    read(10) ibdate,ibtime, flexversion
+    write(*,*) 'date,time=',ibdate,ibtime
+    print*, 'flexversion=',flexversion
+
+    print*,'finfo on output interval, averaging time, sampling time'
+! 85
+    read(10) loutstep,loutaver,loutsample
+    print*,'loutstep, aver, sample=',loutstep,loutaver,loutsample
+
+    print*, 'information on output grid setup'
+    read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+    write(*,*) 'outlon0,outlat0,numxgrid,numygrid,dxout,dyout'
+    write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
+    read(10) numzgrid,(outheight(i),i=1,numzgrid)
+    write(*,*) 'numzgrid=', numzgrid
+    write(*,*) 'outheight=', outheight
+
+    read(10) jjjjmmdd,ihmmss
+    write(*,*) 'jjjjmmdd,ihmmss=',jjjjmmdd,ihmss
+
+    print*, 'information on species'
+! 102  3*nspec,maxpointspec_act implicit?
+    read(10) nspec, maxpointspec_act
+    print*,'3*nspec,maxpointspec_act=',nspec,maxpointspec_act
+
+    nspec=nspec/3
+    do 8 n=1,nspec
+        read(10) numzgrid,species(n)
+        print*,'numzgrid,species(n)',numzgrid,species(n)
+        read(10) numzgrid,species(n)
+        print*,'numzgrid,species(n)',numzgrid,species(n)
+        read(10) numzgrid,species(n)
+        print*,'numzgrid,species(n)',numzgrid,species(n)
+    8 END DO
+
+
+    print*, 'information on release points:'
+    read(10) numpoint !l113
+    print*,'numpoint=',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)
+            read(10) xmass(i,j)
+    13 END DO
+    print*,'! Write information on some model switches'
+    read(10) method,lsubgrid,lconvection,ind_source,ind_receptor
+    print*,'method,lsubgrid,lconvection,ind_source,ind_receptor'
+    print*,method,lsubgrid,lconvection,ind_source,ind_receptor
+
+    print*,'!Write age class information'
+    read(10) nageclass,(lage(i),i=1,nageclass)
+    lage(0)=0
+    print*,'nageclass=',nageclass
+    write(*,*) 'lage(0:nageclass)=', (lage(i),i=0,nageclass)
+
+    print*,'!Write topography to output file'
+    do 130 ix=1,numxgrid
+        read(10) (oro(ix,jy),jy=1,numygrid)
+    130 END DO
+    print*,'(oro(ix,jy)'
+    print*,'numxgrid,numygrid=',numxgrid,numygrid
+    print*,'nxmax,nymax=',nxmax,nymax
+          
+    close(10)
+
+
+
+! rite(*,*) (outheight(i),i=1,numzgrid)
+
+! f (loutstep.lt.0) nspec=numpoint
+
+
+    print*,' !if (loutstep < 0) nspec=numpoint'
+    print*,'outheight(0)=0. ???'
+     
+    print*,'c Calculate height, which is outheight plus topography'
+
+! Calculate height, which is outheight plus topography
+!*****************************************************
+
+    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
+
+    print*,'C Determine area of each output grid box'
+
+! 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
+
+    print*,'end subroutine readheader10f'
+    end subroutine readheader10f

readheader_stdout.f

+      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,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
+
+c      do 50 mother_or_nest=1,2
+      do 50 mother_or_nest=1,1
+
+**********************************
+C Read FLEXPART header information
+**********************************
+
+        if (mother_or_nest.eq.1) then
+          filename=dirname(1:len)//'header'
+        else if (mother_or_nest.eq.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,nxmax,numxgrid,nymax,numygrid,nzmax='  
+        print*,filename,nxmax,numxgrid,nymax,numygrid,nzmax  
+        !print*,'numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime'
+        !print*,numzgrid,outlon0,outlat0,dxout,dyout,outheight,ibdate,ibtime
+        print*,''
+        !print*,
+        print*,''
+        !print*,
+        print*,''
+        !print*,
+
+        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        continue
+
+
+!      print*, 'about to exit' 
+!      exit
+
+50      continue
+
+      if (.not.dirnameisarg) then
+      goto 200 !read(15, dirlist
+199   close(15) !dirlist
+      endif 
+
+         !max_max_grid=max(max_max_grid,max_grid)
+         !sum_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

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.