...
 
Commits (8)
&COMMAND
LDIRECT= 1,
IBDATE= 20140907,
IBTIME= 060000,
IEDATE= 20140907,
IETIME= 120000,
LOUTSTEP= 3600,
LOUTAVER= 3600,
LOUTSAMPLE= 900,
ITSPLIT= 99999999,
LSYNCTIME= 900,
CTL= -5.0000000 ,
IFINE= 4,
IOUT= 9,
IPOUT= 0,
LSUBGRID= 1,
LCONVECTION= 1,
LAGESPECTRA= 0,
IPIN= 0,
IOUTPUTFOREACHRELEASE= 0,
IFLUX= 0,
MDOMAINFILL= 0,
IND_SOURCE= 1,
IND_RECEPTOR= 1,
MQUASILAG= 0,
NESTED_OUTPUT= 1,
LINIT_COND= 0,
SURF_ONLY= 0,
CBLFLAG= 0,
OHFIELDS_PATH= "../../flexin/",
LDIRECT= -1,
IBDATE= 20101222,
IBTIME= 000000,
IEDATE= 20110103,
IETIME= 000000,
LOUTSTEP= 86400,
LOUTAVER= 86400,
LOUTSAMPLE= 1800,
ITSPLIT= 999999999,
LSYNCTIME= 600,
CTL= -5.00,
IFINE= 4,
IOUT= 1,
IPOUT= 2,
LSUBGRID= 1,
LCONVECTION= 1,
LAGESPECTRA= 1,
IPIN= 0,
IOUTPUTFOREACHRELEASE= 1,
IFLUX= 0,
MDOMAINFILL= 0,
IND_SOURCE= 1,
IND_RECEPTOR= 2,
MQUASILAG= 0,
NESTED_OUTPUT= 0,
LINIT_COND= 2,
SURF_ONLY= 1,
LINVERSIONOUT= 1,
OHFIELDS_PATH= "/xnilu_wrk/rlt/NILU/FLEXDATA/",
/
&SPECIES_PARAMS
PSPECIES="C-CO2 ",
PWEIGHTMOLAR= 12.0000000 ,
/
grep Tra SPEC*
SPECIES_001: PSPECIES=TRACER
SPECIES_002: PSPECIES=O3
SPECIES_003: PSPECIES=NO
SPECIES_004: PSPECIES=NO2
SPECIES_005: PSPECIES=HNO3
SPECIES_006: PSPECIES=HNO2
SPECIES_007: PSPECIES=H2O2
SPECIES_009: PSPECIES=HCHO
SPECIES_010: PSPECIES=PAN
SPECIES_011: PSPECIES=NH3
SPECIES_012: PSPECIES=SO4-aero
SPECIES_013: PSPECIES=NO3-aero
SPECIES_014: PSPECIES=I2-131
SPECIES_015: PSPECIES=I-131
SPECIES_016: PSPECIES=Cs-137
SPECIES_017: PSPECIES=Y-91
SPECIES_018: PSPECIES=Ru-106
SPECIES_019: PSPECIES=Kr-85
SPECIES_020: PSPECIES=Sr-90
SPECIES_021: PSPECIES=Xe-133
SPECIES_022: PSPECIES=CO
SPECIES_023: PSPECIES=SO2dryOH
SPECIES_024: PSPECIES=AIRTRACER
SPECIES_025: PSPECIES=AERO-TRACE
SPECIES_026: PSPECIES=CH4
SPECIES_027: PSPECIES=C2H6
SPECIES_028: PSPECIES=C3H8
SPECIES_029: PSPECIES=C-CO2
SPECIES_031: PSPECIES=PCB28
SPECIES_034: PSPECIES=G-HCH
......@@ -220,6 +220,11 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
endif
ixp=ix+1
jyp=jy+1
! fix (Espen 02.05.2017)
if(jyp>=nymax) then
write(*,*) 'WARNING: advance.f90 jyp>nymax'
jyp=jyp-1
endif
! Compute maximum mixing height around particle position
......
......@@ -120,6 +120,9 @@ module com_mod
integer :: lnetcdfout
! lnetcdfout 1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
integer :: linversionout
! linversionout 1 for one grid_time output file for each release containing all timesteps
integer :: nageclass,lage(maxageclass)
! nageclass number of ageclasses for the age spectra calculation
......@@ -356,7 +359,9 @@ module com_mod
real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !liquid [kg/kg]
real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice [kg/kg]
real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !combined [m3/m3]
! RLT add pressure and dry air density
real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
......@@ -378,6 +383,7 @@ module com_mod
! uu,vv,ww [m/2] wind components in x,y and z direction
! uupol,vvpol [m/s] wind components in polar stereographic projection
! tt [K] temperature data
! prs air pressure
! qv specific humidity data
! pv (pvu) potential vorticity
! rho [kg/m3] air density
......
......@@ -100,6 +100,11 @@ subroutine conccalc(itime,weight)
jy=int(ytra1(i))
ixp=ix+1
jyp=jy+1
! fix (Espen 02.05.2017)
if(jyp>=nymax) then
write(*,*) 'WARNING: conccalc.f90 jyp>nymax'
jyp=jyp-1
endif
ddx=xtra1(i)-real(ix)
ddy=ytra1(i)-real(jy)
rddx=1.-ddx
......
......@@ -71,6 +71,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -105,6 +108,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! mind eso:added to ensure identical results between 2&3-fields versions
character(LEN=8),save :: file_stat='REPLACE'
logical :: ldates_file
logical :: lexist
integer :: ierr
character(LEN=100) :: dates_char
......@@ -200,6 +204,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
end do
end do
end do
......@@ -211,8 +218,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -351,7 +364,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
if ((ldirect.eq.1).and.(WETDEP)) then
do jy=0,numygrid-1
do ix=0,numxgrid-1
!oncentraion greater zero
!concentraion greater zero
if (wetgrid(ix,jy).gt.smallnum) then
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
......@@ -594,6 +607,49 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgrid+kz*numxgrid*numygrid
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal
......@@ -620,7 +676,23 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
endif
! RLT Aug 2017
! Write out conversion factor for dry air
if (numreceptor.gt.0) then
inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='new',action='write')
endif
write(unitoutfactor) itime
write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
close(unitoutfactor)
endif
! Reinitialization of grid
!*************************
......
This diff is collapsed.
This diff is collapsed.
......@@ -69,6 +69,9 @@ subroutine concoutput_nest(itime,outnum)
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -95,6 +98,7 @@ subroutine concoutput_nest(itime,outnum)
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
integer :: mind
! mind eso:added to ensure identical results between 2&3-fields versions
......@@ -163,6 +167,9 @@ subroutine concoutput_nest(itime,outnum)
! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
end do
end do
end do
......@@ -174,8 +181,14 @@ subroutine concoutput_nest(itime,outnum)
jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -539,6 +552,47 @@ subroutine concoutput_nest(itime,outnum)
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,numzgrid
do jy=0,numygridn-1
do ix=0,numxgridn-1
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgridn+kz*numxgridn*numygridn
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
! Reinitialization of grid
......
......@@ -71,6 +71,9 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -100,6 +103,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
if (verbosity.eq.1) then
......@@ -179,6 +183,9 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
rho(iix,jjy,kzz-1,2)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
end do
end do
end do
......@@ -189,8 +196,14 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,2)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,2)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -459,6 +472,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
endif ! concentration output
......@@ -503,6 +517,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
......@@ -595,6 +610,49 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,1
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgrid+kz*numxgrid*numygrid
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal
......@@ -621,7 +679,23 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
endif
! RLT Aug 2017
! Write out conversion factor for dry air
if (numreceptor.gt.0) then
inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='new',action='write')
endif
write(unitoutfactor) itime
write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
close(unitoutfactor)
endif
! Reinitialization of grid
!*************************
......
......@@ -69,6 +69,9 @@ subroutine concoutput_surf_nest(itime,outnum)
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -95,7 +98,7 @@ subroutine concoutput_surf_nest(itime,outnum)
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
! Determine current calendar date, needed for the file name
!**********************************************************
......@@ -158,18 +161,27 @@ subroutine concoutput_surf_nest(itime,outnum)
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
rho(iix,jjy,kzz-1,2)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
end do
end do
end do
do i=1,numreceptor
xl=xreceptor(i)
yl=yreceptor(i)
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,2)
end do
do i=1,numreceptor
xl=xreceptor(i)
yl=yreceptor(i)
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,2)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,2)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -316,8 +328,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! Dry deposition
sp_count_i=0
......@@ -353,8 +365,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
......@@ -398,8 +410,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
else
! write full vertical resolution
......@@ -439,8 +451,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
endif ! surf_only
......@@ -486,8 +498,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! Dry deposition
......@@ -525,8 +537,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! Mixing ratios
......@@ -569,8 +581,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
else
! write full vertical resolution
......@@ -610,8 +622,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
endif ! surf_only
endif ! output for ppt
......@@ -624,6 +636,48 @@ subroutine concoutput_surf_nest(itime,outnum)
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,1
do jy=0,numygridn-1
do ix=0,numxgridn-1
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgridn+kz*numxgridn*numygridn
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
! Reinitialization of grid
......
......@@ -79,6 +79,10 @@ subroutine getfields(itime,nstop)
real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
! RLT added partial pressure water vapor
real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
integer :: kz, ix
character(len=100) :: rowfmt
integer :: indmin = 1
......@@ -170,6 +174,29 @@ subroutine getfields(itime,nstop)
endif
! RLT calculate dry air density
pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
rho_dry=(prs-pwater)/(r_air*tt)
! test density
! write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
! if(itime.eq.0) then
! open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
! do kz=1,nzmax
! do ix=1,nxmax
! write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
! end do
! end do
! close(500)
! open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
! do kz=1,nzmax
! do ix=1,nxmax
! write(500,fmt=rowfmt) rho(ix,:,kz,1)
! end do
! end do
! close(500)
! endif
lwindinterv=abs(memtime(2)-memtime(1))
if (lwindinterv.gt.idiffmax) nstop=3
......
......@@ -41,16 +41,25 @@ subroutine gethourlyOH(itime)
implicit none
integer :: itime
integer :: ix,jy,kz,m1,m2
integer :: ix,jy,kz,m1,m2,d1,d2,m,dd,ndays
integer :: ijx,jjy
integer :: jjjjmmdd,hhmmss
real :: sza,jrate,photo_O1D,zenithangle
real :: sza,jrate,photo_O1D,zenithangle,OH_tmp,jrate_tmp
real(kind=dp) :: jul1,jul2
integer, dimension(12) :: monthdays
monthdays=(/31,28,31,30,31,30,31,31,30,31,30,31/)
! print*, 'itime: ',itime
! print*, 'memOHtime(1):',memOHtime(1)
! print*, 'memOHtime(2):',memOHtime(2)
!! testing OH
! open(1000,file=path(2)(1:length(2))//'OH_50N_0E.txt',action='write',status='old',access='append')
! open(999,file=path(2)(1:length(2))//'jrate_50N_0E.txt',action='write',status='old',access='append')
! open(998,file=path(2)(1:length(2))//'OH_50S_0E.txt',action='write',status='old',access='append')
! open(997,file=path(2)(1:length(2))//'jrate_50S_0E.txt',action='write',status='old',access='append')
! Check hourly OH field is available for the current time step
!**************************************************************
......@@ -79,9 +88,13 @@ subroutine gethourlyOH(itime)
jul2=bdate+memOHtime(2)/86400._dp ! date for next hour
call caldate(jul2,jjjjmmdd,hhmmss)
m2=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
dd=(jjjjmmdd-(jjjjmmdd/100)*100)
ndays=monthdays(m2)
! print*, 'jul2:',jul2
! print*, 'm2:',m2
! print*, 'dd:',dd
! print*, 'ndays:',ndays
do kz=1,nzOH
do jy=1,nyOH
......@@ -92,19 +105,33 @@ subroutine gethourlyOH(itime)
sza=zenithangle(latOH(jy),lonOH(ix),jul2)
! calculate J(O1D) (jrate)
jrate=photo_O1D(sza)
! interpolate monthly fields to current day
if(dd.ge.(ndays/2)) then
m=min((m2+1),12)
OH_tmp=OH_field(ix,jy,kz,m2)+(OH_field(ix,jy,kz,m)-OH_field(ix,jy,kz,m2))*real(dd-ndays/2)/real(ndays)
jrate_tmp=jrate_average(ijx,jjy,m2)+(jrate_average(ijx,jjy,m)-jrate_average(ijx,jjy,m2))*real(dd-ndays/2)/real(ndays)
else
m=max((m2-1),1)
OH_tmp=OH_field(ix,jy,kz,m)+(OH_field(ix,jy,kz,m2)-OH_field(ix,jy,kz,m))*real(dd+ndays/2)/real(ndays)
jrate_tmp=jrate_average(ijx,jjy,m)+(jrate_average(ijx,jjy,m2)-jrate_average(ijx,jjy,m))*real(dd+ndays/2)/real(ndays)
endif
! apply hourly correction to OH
if(jrate_average(ijx,jjy,m2).gt.0.) then
OH_hourly(ix,jy,kz,2)=OH_field(ix,jy,kz,m2)*jrate/jrate_average(ijx,jjy,m2)
if(jrate_tmp.gt.0.) then
OH_hourly(ix,jy,kz,2)=OH_tmp*jrate/jrate_tmp
else
OH_hourly(ix,jy,kz,2)=0.
endif
!! for testing !!
! if(jy.eq.36.and.ix.eq.36.and.kz.eq.1) then
! write(999,fmt='(F6.3)') jrate/jrate_average(ijx,jjy,m2)
! endif
! if(jy.eq.11.and.ix.eq.36.and.kz.eq.1) then
! write(998,fmt='(F6.3)') jrate/jrate_average(ijx,jjy,m2)
! endif
!! testing OH
! if(jy.eq.36.and.ix.eq.37.and.kz.eq.1) then
! write(999,fmt='(E12.6)') OH_hourly(ix,jy,kz,2)
! write(1000,fmt='(E12.6)') OH_tmp
! write(999,fmt='(E12.6)') jrate/jrate_tmp
! endif
! if(jy.eq.11.and.ix.eq.37.and.kz.eq.1) then
! write(998,fmt='(E12.6)') OH_hourly(ix,jy,kz,2)
! write(998,fmt='(E12.6)') OH_tmp
! write(997,fmt='(E12.6)') jrate/jrate_tmp
! endif
end do
end do
end do
......@@ -118,15 +145,18 @@ subroutine gethourlyOH(itime)
call caldate(jul1,jjjjmmdd,hhmmss)
m1=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
memOHtime(1)=0.
d1=(jjjjmmdd-(jjjjmmdd/100)*100)
jul2=bdate+ldirect*real(1./24.,kind=dp) ! date for next hour
call caldate(jul2,jjjjmmdd,hhmmss)
m2=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
memOHtime(2)=ldirect*3600.
d2=(jjjjmmdd-(jjjjmmdd/100)*100)
! print*, 'jul1:',jul1
! print*, 'jul2:',jul2
! print*, 'm1,m2:',m1,m2
! print*, 'd1,d2:',d1,d2
do kz=1,nzOH
do jy=1,nyOH
......@@ -137,9 +167,18 @@ subroutine gethourlyOH(itime)
sza=zenithangle(latOH(jy),lonOH(ix),jul1)
! calculate J(O1D) (jrate), beginning
jrate=photo_O1D(sza)
! interpolate monthly fields to current day
ndays=monthdays(m1)
if(d1.ge.(ndays/2)) then
m=min((m1+1),12)
OH_tmp=OH_field(ix,jy,kz,m1)+(OH_field(ix,jy,kz,m)-OH_field(ix,jy,kz,m1))*real(d1-ndays/2)/real(ndays)
else
m=max((m1-1),1)
OH_tmp=OH_field(ix,jy,kz,m)+(OH_field(ix,jy,kz,m1)-OH_field(ix,jy,kz,m))*real(d1+ndays/2)/real(ndays)
endif
! apply hourly correction to OH
if(jrate_average(ijx,jjy,m1).gt.0.) then
OH_hourly(ix,jy,kz,1)=OH_field(ix,jy,kz,m1)*jrate/jrate_average(ijx,jjy,m1)
OH_hourly(ix,jy,kz,1)=OH_tmp*jrate/jrate_average(ijx,jjy,m1)
else
OH_hourly(ix,jy,kz,1)=0.
endif
......@@ -147,9 +186,18 @@ subroutine gethourlyOH(itime)
sza=zenithangle(latOH(jy),lonOH(ix),jul2)
! calculate J(O1D) (jrate), after 1-hour
jrate=photo_O1D(sza)
! interpolate monthly fields to current day
ndays=monthdays(m2)
if(d2.ge.(ndays/2)) then
m=min((m2+1),12)
OH_tmp=OH_field(ix,jy,kz,m2)+(OH_field(ix,jy,kz,m)-OH_field(ix,jy,kz,m2))*real(d2-ndays/2)/real(ndays)
else
m=max((m2-1),1)
OH_tmp=OH_field(ix,jy,kz,m)+(OH_field(ix,jy,kz,m2)-OH_field(ix,jy,kz,m))*real(d2+ndays/2)/real(ndays)
endif
! apply hourly correction to OH
if(jrate_average(ijx,jjy,m2).gt.0.) then
OH_hourly(ix,jy,kz,2)=OH_field(ix,jy,kz,m2)*jrate/jrate_average(ijx,jjy,m2)
OH_hourly(ix,jy,kz,2)=OH_tmp*jrate/jrate_average(ijx,jjy,m2)
else
OH_hourly(ix,jy,kz,2)=0.
endif
......@@ -159,5 +207,10 @@ subroutine gethourlyOH(itime)
endif
!! testing OH
! close(999)
! close(998)
end subroutine gethourlyOH
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine initial_cond_output_inversion(itime)
! i
!*****************************************************************************
! *
! Output of the initial condition sensitivity field. *
! *
! Author: A. Stohl *
! *
! 24 May 1995 *
! *
! 13 April 1999, Major update: if output size is smaller, dump output *
! in sparse matrix format; additional output of *
! uncertainty *
! *
! 05 April 2000, Major update: output of age classes; output for backward*
! runs is time spent in grid cell times total mass of *
! species. *
! *
! 17 February 2002, Appropriate dimensions for backward and forward runs *
! are now specified in file par_mod *
! *
! June 2006, write grid in sparse matrix with a single write command *
! in order to save disk space *
! *
! 2008 new sparse matrix format *
! *
!*****************************************************************************
! *
! Variables: *
! ncells number of cells with non-zero concentrations *
! sparse .true. if in sparse matrix format, else .false. *
! *
!*****************************************************************************
use unc_mod
use point_mod
use outg_mod
use par_mod
use com_mod
implicit none
integer :: itime,i,ix,jy,kz,ks,kp,sp_count_i,sp_count_r
integer :: jjjjmmdd, ihmmss
real(kind=dp) :: jul
real :: sp_fact,fact_recept
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
logical :: sp_zer,lexist
logical,save :: listart=.true.
logical,save,allocatable,dimension(:) :: listartrel
character :: adate*8,atime*6
character :: areldate*8,areltime*6
character(len=3) :: anspec
if(listart) then
allocate(listartrel(maxpointspec_act))
listartrel(:)=.true.
endif
print*, 'listartrel = ',listartrel
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
do ks=1,nspec
write(anspec,'(i3.3)') ks
do kp=1,maxpointspec_act
! calculate date of release
jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp ! this is the current day
call caldate(jul,jjjjmmdd,ihmmss)
write(areldate,'(i8.8)') jjjjmmdd
write(areltime,'(i6.6)') ihmmss
print*, areldate//areltime
! calculate date of field
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,ihmmss)
write(adate,'(i8.8)') jjjjmmdd
write(atime,'(i6.6)') ihmmss
print*, adate//atime
inquire(file=path(2)(1:length(2))//'grid_initial_'//areldate// &
areltime//'_'//anspec,exist=lexist)
if(lexist.and..not.listartrel(kp)) then
! open and append to existing file
open(97,file=path(2)(1:length(2))//'grid_initial_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
else
! open new file
open(97,file=path(2)(1:length(2))//'grid_initial_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='replace',action='write')
endif
write(97) jjjjmmdd
write(97) ihmmss
listartrel(kp)=.false.
if (ind_rel.eq.1) then
fact_recept=rho_rel(kp)
else
fact_recept=1.
endif
!*******************************************************************
! Generate output: may be in concentration (ng/m3) or in mixing
! ratio (ppt) or both
! Output the position and the values alternated multiplied by
! 1 or -1, first line is number of values, number of positions
! For backward simulations, the unit is seconds, stored in grid_time
!*******************************************************************
! Write out dummy "wet and dry deposition" fields, to keep same format
! as for concentration output
! sp_count_i=0
! sp_count_r=0
! write(97) sp_count_i
! write(97) (sparse_dump_i(i),i=1,sp_count_i)
! write(97) sp_count_r
! write(97) (sparse_dump_r(i),i=1,sp_count_r)
! write(97) sp_count_i
! write(97) (sparse_dump_i(i),i=1,sp_count_i)
! write(97) sp_count_r
! write(97) (sparse_dump_r(i),i=1,sp_count_r)
! Write out sensitivity to initial conditions
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (init_cond(ix,jy,kz,ks,kp).gt.smallnum) then
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgrid+kz*numxgrid*numygrid
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)=sp_fact* &
init_cond(ix,jy,kz,ks,kp)/xmass(kp,ks)*fact_recept
else ! concentration is zero
sp_zer=.true.
endif
end do
end do
end do
write(97) sp_count_i
write(97) (sparse_dump_i(i),i=1,sp_count_i)
write(97) sp_count_r
write(97) (sparse_dump_r(i),i=1,sp_count_r)
close(97)
end do
end do
! reset listart
if (listart) then
listart=.false.
endif
end subroutine initial_cond_output_inversion
......@@ -50,13 +50,12 @@ ifeq ($(gcc), 4.9)
LIBPATH1 = ${ROOT_DIR}/lib
else
# Default: System libraries at NILU, gfortran v4.6
# Gfortran v5.4
F90 = /usr/bin/gfortran
MPIF90 = /usr/bin/mpif90.openmpi
INCPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/include
MPIF90 = /usr/bin/mpifort
INCPATH2 = /usr/include
LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
INCPATH1 = /homevip/flexpart/gcc-5.4.0/include
LIBPATH1 = /homevip/flexpart/gcc-5.4.0/lib
endif
......@@ -90,7 +89,7 @@ mpi_mod.o
OBJECTS_SERIAL = \
releaseparticles.o partoutput.o \
conccalc.o \
init_domainfill.o concoutput.o \
init_domainfill.o concoutput.o \
timemanager.o FLEXPART.o \
readpartpositions.o \
partoutput_short.o \
......@@ -98,6 +97,8 @@ OBJECTS_SERIAL = \
boundcond_domainfill.o \
redist.o \
concoutput_surf.o concoutput_surf_nest.o \
concoutput_inversion_nest.o \
concoutput_inversion.o \
getfields.o
## For MPI version
......@@ -191,7 +192,7 @@ writeheader_nest.o writeheader_nest_surf.o \
wetdepokernel_nest.o \
drydepokernel_nest.o zenithangle.o \
ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o \
initial_cond_calc.o initial_cond_output.o initial_cond_output_inversion.o \
dynamic_viscosity.o get_settling.o \
initialize_cbl_vel.o re_initialize_particle.o \
cbl.o netcdf_output_mod.o
......@@ -293,9 +294,11 @@ com_mod.o: par_mod.o
conccalc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
conccalc_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o unc_mod.o
concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o mean_mod.o
concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_inversion_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
concoutput_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o mean_mod.o
concoutput_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
......@@ -342,6 +345,7 @@ init_domainfill_mpi.o: com_mod.o mpi_mod.o par_mod.o point_mod.o \
random_mod.o
initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
initial_cond_output.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
initialize.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o random_mod.o
initialize_cbl_vel.o: com_mod.o par_mod.o random_mod.o
interpol_all.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o
......
......@@ -36,6 +36,8 @@ module outg_mod
real,allocatable, dimension (:,:,:) :: areaeast
real,allocatable, dimension (:,:,:) :: areanorth
real,allocatable, dimension (:,:,:) :: densityoutgrid
real,allocatable, dimension (:,:,:) :: densitydrygrid ! added RLT
real,allocatable, dimension (:,:,:) :: factor_drygrid ! added RLT
real,allocatable, dimension (:,:,:) :: factor3d
real,allocatable, dimension (:,:,:) :: grid
real(sp),allocatable, dimension (:,:) :: wetgrid
......
......@@ -251,6 +251,13 @@ subroutine outgrid_init
allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
! RLT
allocate(densitydrygrid(0:max(numxgrid,numxgridn)-1, &
0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
allocate(factor_drygrid(0:max(numxgrid,numxgridn)-1, &
0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
allocate(factor3d(0:max(numxgrid,numxgridn)-1, &
0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
......
......@@ -79,6 +79,9 @@ module par_mod
real,parameter :: pi=3.14159265, r_earth=6.371e6, r_air=287.05, ga=9.81
real,parameter :: cpa=1004.6, kappa=0.286, pi180=pi/180., vonkarman=0.4
! additional constants RLT Aug-2017
real,parameter :: rgas=8.31447
real,parameter :: r_water=461.495
! pi number "pi"
! pi180 pi/180.
......@@ -88,6 +91,8 @@ module par_mod
! cpa specific heat for dry air
! kappa exponent of formula for potential temperature
! vonkarman von Karman constant
! rgas universal gas constant [J/mol/K]
! r_water specific gas constant for water vapor [J/kg/K]
real,parameter :: karman=0.40, href=15., convke=2.0
real,parameter :: hmixmin=100., hmixmax=4500., turbmesoscale=0.16
......@@ -260,6 +265,8 @@ module par_mod
integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
integer,parameter :: unitboundcond=89
integer,parameter :: unittmp=101
! RLT
integer,parameter :: unitoutfactor=102
!******************************************************
! integer code for missing values, used in wet scavenging (PS, 2012)
......
......@@ -111,6 +111,7 @@ subroutine readcommand
lnetcdfout, &
surf_only, &
cblflag, &
linversionout, &
ohfields_path
! Presetting namelist command
......@@ -143,6 +144,7 @@ subroutine readcommand
lnetcdfout=0
surf_only=0
cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine
linversionout=0
ohfields_path="../../flexin/"
! Open the command file and read user options
......@@ -232,6 +234,8 @@ subroutine readcommand
if (old) call skplin(3,unitcommand)
read(unitcommand,*) linit_cond
if (old) call skplin(3,unitcommand)
read(unitcommand,*) linversionout !added by RT
if (old) call skplin(3,unitcommand)
read(unitcommand,*) surf_only
! Removed for backwards compatibility.
! if (old) call skplin(3,unitcommand) !added by mc
......@@ -448,6 +452,16 @@ subroutine readcommand
stop
endif
! Inversion output format only for backward runs
!*****************************************************************************
if ((linversionout.eq.1).and.(ldirect.eq.1)) then
write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
write(*,*) '#### INVERSION OUTPUT FORMAT ONLY FOR ####'
write(*,*) '#### BACKWARD RUNS ####'
stop
endif
! For domain-filling trajectories, a plume centroid trajectory makes no sense,
! For backward runs, only residence time output (iout=1) or plume trajectories (iout=4),
......
......@@ -114,6 +114,7 @@ subroutine timemanager
real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
real :: xold,yold,zold,xmassfract
real, parameter :: e_inv = 1.0/exp(1.0)
logical :: lexist
!double precision xm(maxspec,maxpointspec_act),
! + xm_depw(maxspec,maxpointspec_act),
! + xm_depd(maxspec,maxpointspec_act)
......@@ -143,7 +144,16 @@ subroutine timemanager
! species_lifetime(:,:)=0
! print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
!! RLT 2018
!! testing OH
! if(OHREA) then
! open(1000,file=path(2)(1:length(2))//'OH_50N_0E.txt',action='write',status='new')
! open(999,file=path(2)(1:length(2))//'jrate_50N_0E.txt',action='write',status='new')
! open(998,file=path(2)(1:length(2))//'OH_50S_0E.txt',action='write',status='new')
! open(997,file=path(2)(1:length(2))//'jrate_50S_0E.txt',action='write',status='new')
! endif
write(*,46) float(itime)/3600,itime,numpart
if (verbosity.gt.0) then
......@@ -348,6 +358,11 @@ subroutine timemanager
! Check whether concentrations are to be calculated
!**************************************************
print*, 'itime:',itime
print*, 'loutstart:',loutstart
print*, 'loutend:',loutend
print*, 'loutnext:',loutnext
if ((ldirect*itime.ge.ldirect*loutstart).and. &
(ldirect*itime.le.ldirect*loutend)) then ! add to grid
if (mod(itime-loutstart,loutsample).eq.0) then
......@@ -391,7 +406,16 @@ subroutine timemanager
if (lnetcdfout.eq.1) then
call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
else
call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
if (linversionout.eq.1) then
call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
if (verbosity.eq.1) then
print*,'called concoutput_inversion'
call system_clock(count_clock)
write(*,*) 'system clock',count_clock - count_clock0
endif
else
call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
endif
if (verbosity.eq.1) then
print*,'called concoutput_surf '
call system_clock(count_clock)
......@@ -404,8 +428,12 @@ subroutine timemanager
if (lnetcdfout.eq.0) then
if (surf_only.ne.1) then
call concoutput_nest(itime,outnum)
else
call concoutput_surf_nest(itime,outnum)
else
if(linversionout.eq.1) then
call concoutput_inversion_nest(itime,outnum)
else
call concoutput_surf_nest(itime,outnum)
endif
endif
else
if (surf_only.ne.1) then
......@@ -668,7 +696,13 @@ subroutine timemanager
if (ipout.eq.2) call partoutput(itime) ! dump particle positions
if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field
if (linit_cond.ge.1) then
if(linversionout.eq.1) then
call initial_cond_output_inversion(itime) ! dump initial cond. field
else
call initial_cond_output(itime) ! dump initial cond. fielf
endif
endif
!close(104)
......
......@@ -75,6 +75,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
! RLT added pressure
real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: prsh
real,dimension(0:nxmax-1,0:nymax-1) :: tvold,pold,pint,tv
real,dimension(0:nymax-1) :: cosf
......@@ -185,17 +187,19 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! Loop over the whole grid
!*************************
do jy=0,nymin1
do ix=0,nxmin1
tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ &
ps(ix,jy,1,n))
enddo
enddo
pold=ps(:,:,1,n)
uvzlev(:,:,1)=0.
wzlev(:,:,1)=0.
rhoh(:,:,1)=pold/(r_air*tvold)
! RLT add pressure
prsh(:,:,1)=ps(:,:,1,n)
! Compute heights of eta levels
......@@ -203,6 +207,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
do kz=2,nuvz
pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
! RLT add pressure
prsh(:,:,kz)=pint
tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
......@@ -254,6 +260,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!hg
pv(:,:,1,n)=pvh(:,:,1)
rho(:,:,1,n)=rhoh(:,:,1)
! RLT add pressure
prs(:,:,1,n)=prsh(:,:,1)
uu(:,:,nz,n)=uuh(:,:,nuvz)
vv(:,:,nz,n)=vvh(:,:,nuvz)
......@@ -267,7 +275,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!hg
pv(:,:,nz,n)=pvh(:,:,nuvz)
rho(:,:,nz,n)=rhoh(:,:,nuvz)
! RLT
prs(:,:,nz,n)=prsh(:,:,nuvz)
kmin=2
idx=kmin
......@@ -287,6 +296,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!hg
pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
! RLT
prs(ix,jy,iz,n)=prs(ix,jy,nz,n)
else
innuvz: do kz=idx(ix,jy),nuvz
if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
......@@ -320,6 +331,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!hg
pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
! RLT add pressure
prs(ix,jy,iz,n)=(prsh(ix,jy,kz-1)*dz2+prsh(ix,jy,kz)*dz1)/dz
endif
enddo
enddo
......
......@@ -87,7 +87,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
!ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
real :: frac_act, liq_frac, dquer_m
real :: frac_act, liq_frac, ice_frac, dquer_m
integer :: blc_count, inc_count
real :: Si_dummy, wetscav_dummy
......@@ -335,16 +335,25 @@ subroutine wetdepo(itime,ltsample,loutnext)
cl=1.6E-6*prec(1)**0.36
endif
!ZHG: Calculate the partition between liquid and water phase water.
!ZHG: Calculate the partition between liquid and water phase water.
if (act_temp .le. 253.) then
liq_frac=0
ice_frac=1
else if (act_temp .ge. 273.) then
liq_frac=1
ice_frac=0
else
liq_frac =((act_temp-273.)/(273.-253.))**2.
! sec bugfix after FLEXPART paper review, liq_frac was 1-liq_frac
! IP bugfix v10.4, calculate ice_frac and liq_frac
ice_frac= ((act_temp-273.)/(273.-253.))**2.
!liq_frac = 1-ice_frac !((act_temp-253.)/(273.-253.))**2.
liq_frac=max(0.,1.-ice_frac)
end if
! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
frac_act =