Commit ee4ec64f authored by Ignacio Pisso's avatar Ignacio Pisso
Browse files

initial commit from /home/ignacio/flexpart/NILU_footprint/PREPARE_FOOTPRINT_laptop/

parents
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
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
SHELL = /bin/bash
MAIN = prepare_footprint.exe
#
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 = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
# LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
# linker flags not necessary: no need for grib files
OBJECTS = prepare_footprint.o\
readheader.o \
readgrid.o \
juldate.o \
caldate.o \
$(MAIN): $(OBJECTS)
$(FC) *.o -o $(MAIN) $(FFLAGS)
$(OBJECTS): $(INCF)
clean:
rm *.o
program flexpartplot_latlon
parameter(nxmax=360,nymax=200,nzmax=3,maxtime=40)
parameter(maxpoint=5000,maxspec=88,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
*******************
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
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
else
filegrid=dirname(1:len)//'grid_time_nest_'//aday//atime//'_001'
fileout=dirname(1:len)//'footprint_nest_'//aday//atime
endif
write(*,*) filegrid
call readgrid(filegrid,itime,nxmax,numxgrid,nymax,numygrid,
+ nzmax,numzgrid,maxspec,nspec,maxageclass,nageclass,grid,
+ lage)
*****************************************
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'
else
fileout=dirname(1:len)//'footprint_total_nest'
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
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
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
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)
open(10,file=filename,form='unformatted',status='old')
read(10) ibdate,ibtime
write(*,*) ibdate,ibtime
read(10) loutstep,loutaver,loutsample
read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
read(10) numzgrid,(outheight(i),i=1,numzgrid)
read(10) jjjjmmdd,ihmmss
write(*,*) jjjjmmdd,ihmss
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
write(*,*) (lage(i),i=0,nageclass)
do 130 ix=1,numxgrid
130 read(10) (oro(ix,jy),jy=1,numygrid)
close(10)
write(*,*) (outheight(i),i=1,numzgrid)
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
****************************************