Commit 8c7a4e2d authored by Ignacio Pisso's avatar Ignacio Pisso
Browse files

branch flex_read for paper 10.3

parent d79e1653
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()
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
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
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)
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
max_max_grid=0
max_grid=0
sum_grid=0
sum_sum_grid=0
sizegrid_m=0
!print*, 'init sizegrid_m=', sizegrid_m !=0
do 70 k=1,nspec
do 70 ix=0,numxgrid-1
do 70 jy=0,numygrid-1
70 footprint_total(ix,jy,k)=0.
*******************************************
C Loop about all times, given in file dates
*******************************************
print*,'> open dates'
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 (mother_or_nest.eq.1) then
!filegrid=dirname(1:len)//'/grid_time_'//aday//atime//'_001'
filegrid=dirname(1:len)//'/grid_conc_'//aday//atime//'_001'
!filegrid=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
!write(*,*) filegrid
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)
****************************************
C Display read values
****************************************
!sizegrid_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
!print*, 'sum sum(grid)=', sum_sum_grid
!print*, 'max maxval(grid)=', max_max_grid
print*, 'sum_sum_grid=', sum_sum_grid
print*, 'max_max_grid=', max_max_grid
!print*, 'sizegrid_m=', sizegrid_m
if (0.eq.1) then
*****************************************
C Write out footprint files for each date
*****************************************
open(97,file=fileout,form='unformatted')
write(97) itime
C 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 continue
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).gt.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 continue
C 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 continue
close(97)
endif ! Write out footprint files for each date
goto 100 !read dates
99 close(20)
C Dump time-integrated footprint
********************************
if (0.eq.1) then
if (mother_or_nest.eq.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).gt.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 continue
C 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 continue
close(97)
endif ! write footprint total
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
SHELL = /bin/bash
MAIN = prepare_footprint.exe
#MAIN = sum_avg
MAIN = prepare_footprint
#
FC = gfortran
INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
FFLAGS = -g -fcheck=all -Wall -fbacktrace -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
#LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
#FFLAGS = -g -fcheck=all -Wall -fbacktrace -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
FFLAGS = -g -fcheck=all -Wall -fbacktrace -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4
# LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
# linker flags not necessary: no need for grib files
......@@ -16,10 +18,21 @@ readgrid.o \
juldate.o \
caldate.o \
OBJECTS_AUX = readheader.o \
readgrid.o \
juldate.o \
caldate.o \
$(MAIN): $(OBJECTS)
$(FC) *.o -o $(MAIN) $(FFLAGS)
$(OBJECTS): $(INCF)
#$(OBJECTS):
prepare_footprint_modif:readheader.o readgrid.o juldate.o caldate.o prepare_footprint_modif.o
$(FC) *.o -o prepare_footprint_modif
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
clean:
rm *.o
program flexpartplot_latlon
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
character filename*150,fileout*150,aday*8,atime*6
character*250 filegrid
logical sp_zer
real sizegrid_m=0
*******************
print*, " Read CONTROL file "
*******************
open(15,file='dirlist')
200 read(15,'(a)',end=199) dirname
len=index(dirname,' ')-1
do 50 mother_or_nest=1,2
**********************************
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
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)
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))
30 continue
sizegrid_m=0
print*, 'init sizegrid_m=', sizegrid_m
do 70 k=1,nspec
do 70 ix=0,numxgrid-1
do 70 jy=0,numygrid-1
70 footprint_total(ix,jy,k)=0.
*******************************************
C Loop about all times, given in file dates
*******************************************
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 (mother_or_nest.eq.1) then
filegrid=dirname(1:len)//'/grid_time_'//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
write(*,*) filegrid
call readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid,
+ nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid,
+ lage)
****************************************
C Display read values
****************************************
sizegrid_m=sizegrid_m+sum(grid)
write(*,*) 'sum(grid)', sum(grid)
print*, 'sizegrid_m=', sizegrid_m
*****************************************
C Write out footprint files for each date
*****************************************
open(97,file=fileout,form='unformatted')
write(97) itime
C 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 continue
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).gt.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 continue
C 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 continue
close(97)
goto 100
99 close(20)
C Dump time-integrated footprint
********************************
if (mother_or_nest.eq.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).gt.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 continue
C 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 continue
close(97)
50 continue
goto 200
199 close(15)
end
......@@ -17,7 +17,8 @@
open(10,file=filegrid,form='unformatted',status='old')
read(10,end=99) itime
write(*,*) itime
!write(*,*) itime
!print*, 'readgrid> itime=', itime
C Loop about all species
************************
......
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)