...
 
Commits (3)
......@@ -51,7 +51,8 @@ program flexpart
use par_mod
use com_mod
use conv_mod
use oh_mod
use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1
use class_gribfile
......@@ -136,10 +137,6 @@ program flexpart
endif
if (verbosity.gt.0) then
print*, 'nxmax=',nxmax
print*, 'nymax=',nymax
print*, 'nzmax=',nzmax
print*,'nxshift=',nxshift
write(*,*) 'call readpaths'
endif
call readpaths(pathfile)
......@@ -454,21 +451,16 @@ program flexpart
call timemanager(metdata_format)
if (verbosity.gt.0) then
! NIK 16.02.2005
do i=1,nspec
if (tot_inc_count(i).gt.0) then
write(*,*) '**********************************************'
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
write(*,*) '**********************************************'
endif
end do
write (*,*) 'timemanager> call wetdepo'
endif
do i=1,nspec
write(*,*) '**********************************************'
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
write(*,*) '**********************************************'
end do
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
&XPART MODEL RUN!'
......
......@@ -376,6 +376,7 @@ program flexpart
call readOHfield
endif
! Write basic information on the simulation to a file "header"
! and open files that are to be kept open throughout the simulation
!******************************************************************
......
......@@ -221,6 +221,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
! Determine the lower left corner and its distance to the current position
......@@ -244,7 +249,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
! eso: Temporary fix for particle exactly at north pole
if (jyp >= nymax) then
! write(*,*) 'WARNING: advance.f90 jyp >= nymax. xt,yt:',xt,yt
write(*,*) 'WARNING: advance.f90 jyp >= nymax. xt,yt:',xt,yt
jyp=jyp-1
end if
......@@ -538,9 +543,9 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
do nsp=1,nspec
if (xmass(nrelpoint,nsp).gt.eps2) exit
end do
if (nsp.gt.nspec) then
if (nsp.gt.nspec+1) then
! This should never happen
write(*,*) 'advance.f90: ERROR: could not find releasepoint'
write(*,*) 'advance.f90#1: ERROR: could not find releasepoint'
stop
end if
if (density(nsp).gt.0.) then
......@@ -708,9 +713,9 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
do nsp=1,nspec
if (xmass(nrelpoint,nsp).gt.eps2) exit
end do
if (nsp.gt.nspec) then
if (nsp.gt.nspec+1) then
! This should never happen
write(*,*) 'advance.f90: ERROR: could not find releasepoint'
write(*,*) 'advance.f90#2: ERROR: could not find releasepoint'
stop
end if
if (density(nsp).gt.0.) then
......@@ -918,9 +923,9 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
do nsp=1,nspec
if (xmass(nrelpoint,nsp).gt.eps2) exit
end do
if (nsp.gt.nspec) then
if (nsp.gt.nspec+1) then
! This should never happen
write(*,*) 'advance.f90: ERROR: could not find releasepoint'
write(*,*) 'advance.f90#3: ERROR: could not find releasepoint'
stop
end if
if (density(nsp).gt.0.) then
......
......@@ -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
......@@ -173,8 +176,11 @@ module com_mod
real :: vset(maxspec,ni),schmi(maxspec,ni),fract(maxspec,ni)
real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
real :: rm(maxspec),dryvel(maxspec)
real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
real :: use_NH3_loss(maxspec)
! se it is possible to associate a species with a second one to make transfer from gas to aerosol
integer :: spec_ass(maxspec)
real :: area_hour(maxspec,24),point_hour(maxspec,24)
real :: area_dow(maxspec,7),point_dow(maxspec,7)
......@@ -358,7 +364,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)
......@@ -380,6 +388,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
......@@ -577,7 +586,7 @@ module com_mod
integer :: numxgridn,numygridn
real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
!real outheight(maxzgrid),outheighthalf(maxzgrid)
logical :: NH3LOSS
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
& OHREA,ASSSPEC
logical :: DRYBKDEP,WETBKDEP
......@@ -602,7 +611,7 @@ module com_mod
! ASSSPEC .true., if there are two species asscoiated
! DRYBKDEP,WETBKDEP .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition
! (i.e. transfer of mass between these two occurs
! NH3LOSS .true., if the NH3 lifetime shall be reduced according to a given field
! if output for each releasepoint shall be created maxpointspec=number of releasepoints
......
......@@ -101,6 +101,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
......@@ -202,6 +206,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
......@@ -213,8 +220,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
......@@ -362,7 +375,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
......@@ -613,6 +626,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
......@@ -639,7 +695,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
......@@ -550,6 +563,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
......
......@@ -150,9 +150,12 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
do il=2,nz
if (height(il).gt.ztra1(jpart)) then
hz=il-1
! goto 26
exit
endif
end do
!26 continue
if (ngrid.eq.0) then
clouds_v=clouds(ix,jy,hz,n)
......@@ -199,19 +202,26 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
endif
!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp - 2 and 3 not used removed by SE
!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
! for now they are treated the same
grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
grfraction(2)=max(0.05,cc*(lfr(i)))
grfraction(3)=max(0.05,cc*(cfr(j)))
! 2) Computation of precipitation rate in sub-grid cell
!******************************************************
prec(1)=(lsp+convp)/grfraction(1)
prec(2)=(lsp)/grfraction(2)
prec(3)=(convp)/grfraction(3)
! 3) Computation of scavenging coefficients for all species
! Computation of wet deposition
!**********************************************************
if (ngrid.gt.0) then
act_temp=ttn(ix,jy,hz,n,ngrid)
else
......@@ -226,6 +236,7 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
! For gas: if positive below-cloud parameters (A or B), and dquer<=0
!******************************************************************
if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
blc_count(ks)=blc_count(ks)+1
wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
......@@ -259,6 +270,8 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
endif
! write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
endif ! gas or particle
! endif ! positive below-cloud scavenging parameters given in Species file
endif !end BELOW
......@@ -271,6 +284,7 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
! given in species file, or if gas and positive Henry's constant
if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
inc_count(ks)=inc_count(ks)+1
! write(*,*) 'Incloud: ',inc_count
! if negative coefficients (turned off) set to zero for use in equation
if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
if (in_aero(ks).lt.0.) in_aero(ks)=0.
......@@ -283,11 +297,7 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
else !parameterize cloudwater m2/m3
!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
! sec test
! cl=1E6*1E-7*prec(1)**0.3 !Sec GFS new
cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new, is also suitable for GFS
! cl=2E-7*prec(1)**0.36 !Andreas
! cl=1.6E-6*prec(1)**0.36 !Henrik
cl=1.6E-6*prec(1)**0.36
endif
!ZHG: Calculate the partition between liquid and water phase water.
......@@ -307,16 +317,24 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
!********
if (dquer(ks).gt.0.) then
S_i= frac_act/cl
! GAS
!****
else
cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
!REPLACE to switch old/ new scheme
! S_i=frac_act/cle
S_i=1/cle
endif ! gas or particle
! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
!SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
!OLD
if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
else
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
endif
endif ! positive in-cloud scavenging parameters given in Species file
endif !incloud
......
......@@ -86,6 +86,10 @@ subroutine getfields(itime,nstop,metdata_format)
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
......@@ -209,6 +213,29 @@ subroutine getfields(itime,nstop,metdata_format)
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
......@@ -60,10 +60,10 @@ ifeq ($(gcc), 4.9)
INCPATH1 = ${ROOT_DIR}/gcc-4.9.1/include
INCPATH2 = ${ROOT_DIR}/include
LIBPATH1 = ${ROOT_DIR}/lib
else
# Compiled libraries under user ~flexpart, gfortran v5.4
ROOT_DIR = /homevip/flexpart/
else
# Gfortran v5.4
ROOT_DIR = /homevip/flexpart/
F90 = /usr/bin/gfortran
MPIF90 = /usr/bin/mpifort
......@@ -87,7 +87,7 @@ VPATH = gributils/
## OPTIMIZATION LEVEL
O_LEV = 0 # [0,1,2,3,g,s,fast]
O_LEV = 2 # [0,1,2,3,g,s,fast]
O_LEV_DBG = g # [0,g]
## LIBRARIES
......@@ -96,7 +96,7 @@ LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(NCOPT) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) $(NCOPT) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
LDDEBUG = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
......@@ -118,7 +118,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 \
......@@ -126,6 +126,8 @@ OBJECTS_SERIAL = \
boundcond_domainfill.o \
redist.o \
concoutput_surf.o concoutput_surf_nest.o \
concoutput_inversion_nest.o \
concoutput_inversion.o \
getfields.o \
readwind_ecmwf.o
......@@ -167,6 +169,7 @@ interpol_all.o readpaths.o \
getrb.o obukhov.o \
getrc.o convmix.o \
getvdep.o readspecies.o \
read_hourly_NH3field.o \
interpol_misslev.o richardson.o \
scalev.o verttransform_ecmwf.o \
pbl_profile.o readOHfield.o \
......@@ -195,10 +198,10 @@ shift_field.o \
openreceptors.o \
readoutgrid_nest.o \
writeheader_nest.o writeheader_nest_surf.o \
wetdepokernel_nest.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
......@@ -270,9 +273,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
......@@ -319,6 +324,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
......@@ -364,6 +370,7 @@ readdepo.o: com_mod.o par_mod.o
readlanduse.o: com_mod.o par_mod.o
#readlanduse_int1.o: com_mod.o par_mod.o
readOHfield.o: com_mod.o oh_mod.o par_mod.o
read_hourly_NH3field.o: com_mod.o oh_mod.o par_mod.o
readoutgrid.o: com_mod.o outg_mod.o par_mod.o
readoutgrid_nest.o: com_mod.o outg_mod.o par_mod.o
readpartpositions.o: com_mod.o par_mod.o random_mod.o
......
......@@ -66,7 +66,9 @@ module netcdf_output_mod
drydep,wetdep,decay,weta_gas,wetb_gas, numbnests, &
ccn_aero,in_aero, & ! wetc_in,wetd_in, &
reldiff,henry,f0,density,dquer,dsigma,dryvel,&
weightmolar,ohcconst,ohdconst,vsetaver,&
! weightmolar,ohreact,spec_ass,kao,vsetaver,&
weightmolar,ohcconst,ohdconst,spec_ass,kao,vsetaver,&
use_NH3_loss, &
! for concoutput_netcdf and concoutput_nest_netcdf
nxmin1,nymin1,nz,oro,oron,rho,rhon,&
memind,xresoln,yresoln,xrn, xln, yrn,yln,nxn,nyn,&
......@@ -112,8 +114,6 @@ module netcdf_output_mod
logical, parameter :: write_vol = .false.
logical, parameter :: write_area = .false.
! coordinate transformation from internal to world coord
real :: xp1,yp1,xp2,yp2
contains
!****************************************************************
......@@ -512,7 +512,10 @@ subroutine writeheader_netcdf(lnest)
! call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
call nf90_err(nf90_put_att(ncid, sID, 'NH3_loss', use_NH3_loss(i)))
if (lnest) then
specIDn(i) = sID
......@@ -533,7 +536,10 @@ subroutine writeheader_netcdf(lnest)
! call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
call nf90_err(nf90_put_att(ncid, sID, 'NH3_loss', use_NH3_loss(i)))
if (lnest) then
specIDnppt(i) = sID
......@@ -669,14 +675,10 @@ subroutine writeheader_netcdf(lnest)
call nf90_err(nf90_put_var(ncid, relstartID, ireleasestart(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
xp1=xpoint1(i)*dx+xlon0
yp1=ypoint1(i)*dy+ylat0
xp2=xpoint2(i)*dx+xlon0
yp2=ypoint2(i)*dy+ylat0
call nf90_err(nf90_put_var(ncid, rellng1ID, xp1, (/i/)))
call nf90_err(nf90_put_var(ncid, rellng2ID, xp2, (/i/)))
call nf90_err(nf90_put_var(ncid, rellat1ID, yp1, (/i/)))
call nf90_err(nf90_put_var(ncid, rellat2ID, yp2, (/i/)))
call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))
call nf90_err(nf90_put_var(ncid, relpartID, npart(i), (/i/)))
......
......@@ -35,4 +35,9 @@ module oh_mod
real, dimension(360) :: lonjr
real, dimension(180) :: latjr
integer :: nxNH3,nyNH3,nzNH3,ntNH3
real, allocatable, dimension(:) :: lonNH3,latNH3,altNH3
double precision, allocatable, dimension(:) :: timeNH3
real, allocatable, dimension (:,:,:) :: NH3LOSS_field
end module oh_mod
......@@ -50,15 +50,32 @@ subroutine ohreaction(itime,ltsample,loutnext)
integer :: jpart,itime,ltsample,loutnext,ldeltat,j,k,ix,jy!,ijx,jjy
integer :: ngrid,interp_time,n,m,h,indz,i!,ia,il
integer :: jjjjmmdd,hhmmss,OHx,OHy,OHz
real, dimension(nzOH) :: altOHtop
integer :: jjjjmmdd,hhmmss,OHx,OHy,OHz,NH3x,Nh3y,NH3z
logical :: init=.true.
real, dimension(nzOH) :: alttopOH
real, dimension(nzNH3) :: alttopNH3
real :: xlon,ylat
real :: xtn,ytn
real :: restmass,ohreacted,oh_average
real :: ohrate,temp
real :: restmass,ohreacted,oh_average,nh3reachted,nh3_average
real :: ohrate,nh3lifetime,temp,nh3reacted
real, parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
real(kind=dp) :: jul
! quick sanity check
if (init) then
if (max_nh3_loss > 1.0 .or. max_nh3_loss < 0.0) then
write (*,*) 'Invalid max_nh3_loss'
stop
end if
if (max_nh3_gain <= 0.0) then
<