Skip to content
Snippets Groups Projects
Commit 2464aef4 authored by Rona Thompson's avatar Rona Thompson
Browse files

Update to chemistry_mod to use new hourly interpolation scheme for OH

parent 2861005b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 78 additions and 186 deletions
src/chemistry_mod.f90
+ 78
− 186
View file @ 2464aef4
......@@ -31,8 +31,9 @@ module chemistry_mod
real, allocatable, dimension(:,:,:,:) :: clfield_cur ! chemical loss fields current hour
integer, dimension(2) :: memCLtime ! time of fields in memory (sec)
integer :: curCLhour ! current hour since start of simulation
real(kind=dp) :: memCLday ! current day of simulation
real(kind=dp), dimension(2) :: CL_time ! julian date of fields in memory
real, allocatable, dimension(:,:,:) :: jrate_average ! monthly average actinic flux
real, allocatable, dimension(:,:) :: jrate_average ! daily average actinic flux
real, allocatable, dimension(:) :: lonjr, latjr
private :: zenithangle, photo_O1D
......@@ -219,25 +220,6 @@ module chemistry_mod
error stop
endif
if (reag_hourly(nr).gt.0) then
! Read average jrates and store in 2nd position
!**********************************************
write(jr_name,'(A)') trim(reag_path(nr))//'jrate_average.nc'
inquire(file=jr_name,exist=lexist)
if (lexist) then
write(*,*) 'Reading jrate field: ',jr_name
call readjrate(jr_name, memid, mm)
else
write(*,*) '#### FLEXPART ERROR ####'
write(*,*) '#### JRATE_AVERAGE NOT FOUND ####'
write(*,*) '#### '//trim(jr_name)//' ####'
error stop
endif
endif
end do ! nreagent
else
......@@ -278,8 +260,6 @@ module chemistry_mod
CL_time(memid)=jdmid
endif
endif
!! testing
print*, 'getchemfield: memid, jd, jdmid, CL_time = ',memid,jd,jdmid,CL_time(1)
else
call caldate(jd, jjjjmmdd, hhmmss)
eomday=calceomday(jjjjmmdd/100)
......@@ -314,24 +294,6 @@ module chemistry_mod
error stop
endif
! Read average jrates
!********************
if (reag_hourly(nr).gt.0) then
! Read average jrates into memory
write(jr_name,'(A)') trim(reag_path(nr))//'jrate_average.nc'
inquire(file=jr_name,exist=lexist)
if (lexist) then
write(*,*) 'Reading jrate field: ',jr_name
call readjrate(jr_name, memid, mm)
else
write(*,*) '#### FLEXPART ERROR ####'
write(*,*) '#### JRATE_AVERAGE NOT FOUND ####'
write(*,*) '#### '//trim(jr_name)//' ####'
error stop
endif
endif ! reag_hourly
end do ! nreagent
end do ! memid
......@@ -452,44 +414,85 @@ module chemistry_mod
integer :: itime, curhour, interp_time
integer :: nr, kz, ix, jy, ixm, jym, ixp, jyp, indz, indzm, ii
real :: dt1, dt2, dtt, sza, jrate, jrate_cur, r
real :: dt1, dt2, dtt, sza, jrate, r
real, dimension(2) :: r1
real :: dz1, dz2, dz, ddx, ddy, rddx, rddy, p1, p2, p3, p4
integer :: jjjjmmdd, hhmmss
integer :: jrx, jry
real(kind=dp) :: jul
integer :: jjjjmmdd, hhmmss, hh
real(kind=dp) :: jul, curday
real, parameter :: avog = 6.02214e23 ! Avogadro constant (1/mol)
!! testing
! character(len=4) :: atime
! character(len=20) :: frmt
! real, dimension(nxjr,nyjr) :: jscalar, sza_grid
character(len=4) :: atime
real, dimension(nxCL(1),nyCL(1)) :: jscalar
character(len=20) :: frmt
!! testing
! jscalar(:,:)=1.
! sza_grid(:,:)=0.
jscalar(:,:)=0.
! current hour of simulation
curhour=itime/3600
write(*,*) 'getchemhourly: curhour, curCLhour = ',curhour, curCLhour
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,hhmmss)
if ((ldirect*curCLhour.eq.ldirect*curhour).and.(ldirect*itime.gt.0)) then
! chemical field is for current hour -> don't do anything
return
else
! interpolate chemical field
!*****************************
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,hhmmss)
! interpolate to middle of hour
curCLhour=curhour
interp_time=curhour*3600+1800
dt1=float(interp_time-memCLtime(1))
dt2=float(memCLtime(2)-interp_time)
dtt=1./(dt1+dt2)
if (any(reag_hourly.gt.0)) then
! interpolate using rate of O3 + hv -> O1D
! check if daily mean rate exists for current day
! note: assume that all chemical fields using hourly interpolation
! are on the same grid
do nr=1,nreagent
if (reag_hourly(nr).gt.0) exit
end do
curday=bdate+real(itime,kind=dp)/86400._dp
curday=floor(curday)
write(*,*) 'getchemhourly: curday, memCLday = ',curday, memCLday
if (memCLday.ne.curday) then
if (.not.allocated(jrate_average)) then
allocate(jrate_average(nxCL(1),nyCL(1)))
endif
memCLday=curday
jrate_average(:,:)=0.
!$OMP PARALLEL &
!$OMP PRIVATE(jy,ix,hh,sza,jrate) &
!$OMP REDUCTION(+:jrate_average)
!$OMP DO
do jy=1,nyCL(nr)
do ix = 1,nxCL(nr)
do hh=0,23
sza=zenithangle(latCL(jy,nr),lonCL(ix,nr),jjjjmmdd,hh*10000)
jrate=photo_O1D(sza)
jrate_average(ix,jy)=jrate_average(ix,jy)+jrate
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
jrate_average=jrate_average/24.
end if
!! testing
print*, 'getchemhourly: dt1, dt2, dtt = ',dt1,dt2,dtt
! print*, 'getchemhourly: range(jrate_average) = ',minval(jrate_average),maxval(jrate_average)
endif
! loop over reagents
do nr=1,nreagent
write(*,*) 'Interpolating fields for reagent: ',reagents(nr)
clfield_cur(:,:,:,nr)=(dt2*CL_field(:,:,:,nr,1) + dt1*CL_field(:,:,:,nr,2))*dtt
if (reag_unit(nr).eq.'mol/mol') then
! convert to molecule/cm3
!$OMP PARALLEL &
......@@ -522,11 +525,6 @@ module chemistry_mod
p2=ddx*rddy
p3=rddx*ddy
p4=ddx*ddy
!! testing
! if ((ix.lt.5).and.(jy.gt.10.and.jy.lt.20).and.kz.eq.1) then
! print*, 'getchemhourly: lonCL, xlon, latCL, ylat = ',lonCL(ix,nr), xlon0+ixm*dx, latCL(jy,nr), ylat0+jyp*dy
! print*, 'getchemhourly: altCL, height = ',altCL(kz,nr), (dz2*height(indzm)+dz1*height(indzm+1))*dz
! endif
! take density from first field (accurate enough)
do ii=1,2
indz=indzm+ii-1
......@@ -545,38 +543,24 @@ module chemistry_mod
!$OMP END DO
!$OMP END PARALLEL
endif
if (reag_hourly(nr).gt.0) then
! use actinic flux (jrate) for interpolation
! interpolate using rate of O3 + hv -> O1D
!$OMP PARALLEL &
!$OMP PRIVATE(kz,jy,jry,ix,jrx,sza,jrate,jrate_cur)
!$OMP PRIVATE(kz,jy,ix,sza,jrate)
!$OMP DO
do kz=1,nzCL(nr)
do jy=1,nyCL(nr)
! jrate_average dimensions given as grid midpoints
jry=int((latCL(jy,nr)-(latjr(1)-0.5*dyjr))/dyjr)+1
jry=min(jry,nyCL(nr))
!! testing
! if (kz.eq.1.and.jy.lt.10) print*, 'latCL, latjr = ',latCL(jy), latjr(jry)
do ix=1,nxCL(nr)
! jrate_average dimensions given as grid midpoints
jrx=int((lonCL(ix,nr)-(lonjr(1)-0.5*dxjr))/dxjr)+1
jrx=min(jrx,nxCL(nr))
!! testing
! if (kz.eq.1.and.jy.lt.10.and.ix.lt.10) print*, 'lonCL, lonjr = ',lonCL(ix), lonjr(jrx)
! solar zenith angle
sza=zenithangle(latjr(jry),lonjr(jrx),jjjjmmdd,hhmmss)
sza=zenithangle(latCL(jy,nr),lonCL(ix,nr),jjjjmmdd,hhmmss)
! calculate J(O1D) (jrate)
jrate=photo_O1D(sza)
jrate_cur=(dt2*jrate_average(jrx,jry,1) + &
dt1*jrate_average(jrx,jry,2))*dtt
! apply hourly correction to chem field
if(jrate_cur.gt.0.) then
clfield_cur(ix,jy,kz,nr)=clfield_cur(ix,jy,kz,nr)*jrate/jrate_cur
! hourly correction
if (jrate_average(ix,jy).ne.0) then
clfield_cur(ix,jy,kz,nr)=clfield_cur(ix,jy,kz,nr)*jrate/jrate_average(ix,jy)
!! testing
! if (kz.eq.1) then
! jscalar(ix,jy)=jrate/jrate_cur
! sza_grid(ix,jy)=sza
! endif
! jscalar(ix,jy)=jrate/jrate_average(ix,jy)
endif
end do
end do
......@@ -588,30 +572,18 @@ module chemistry_mod
endif ! curhour
!! testing
! write(frmt,fmt='(A,I4,A)') '(',nxCL,'(E14.6))'
! write(atime,fmt='(I4.4)') curhour
! open(600,file=path(2)(1:length(2))//'clfield_'//atime//'.txt',action='write',status='replace')
! do kz=1,nzCL
! do jy=1,nyCL
! write(600,fmt=frmt) clfield_cur(:,jy,kz,1)
! end do
! end do
! close(600)
! write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
! write(frmt,fmt='(A,I4,A)') '(',nxCL(1),'(E14.6))'
! open(600,file=path(2)(1:length(2))//'jscalar_'//atime//'.txt',action='write',status='replace')
! do jy=1,nyjr
! do jy=1,nyCL(1)
! write(600,fmt=frmt) jscalar(:,jy)
! end do
! close(600)
! write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
! open(600,file=path(2)(1:length(2))//'sza_'//atime//'.txt',action='write',status='replace')
! do jy=1,nyjr
! write(600,fmt=frmt) sza_grid(:,jy)
! end do
! close(600)
!!
end subroutine getchemhourly
subroutine chemreaction(itime)
!*****************************************************************************
......@@ -817,73 +789,13 @@ module chemistry_mod
end subroutine chemreaction
subroutine readjrate(jr_name,memid,mm)
!*****************************************************************************
! *
! Reads the average actinic flux rates into memory *
! *
! Author: Rona Thompson, Sep 2023 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! jr_name name of the file *
! memid time index to chemical field variable *
! mm month to read *
! *
!*****************************************************************************
implicit none
character(len=256) :: jr_name
integer :: memid,mm
integer :: ncid,dimid,varid
! Read netcdf file
!******************
! open file
call nf90_err( nf90_open(trim(jr_name),nf90_NOWRITE,ncid) )
! longitude
call nf90_err( nf90_inq_dimid(ncid,'longitude',dimid) )
call nf90_err( nf90_inquire_dimension(ncid,dimid,len=nxjr) )
if (.not.allocated(lonjr)) allocate(lonjr(nxjr))
call nf90_err( nf90_inq_varid(ncid,'longitude',varid) )
call nf90_err( nf90_get_var(ncid,varid,lonjr) )
dxjr=abs(lonjr(2)-lonjr(1))
! latitude
call nf90_err( nf90_inq_dimid(ncid,'latitude',dimid) )
call nf90_err( nf90_inquire_dimension(ncid,dimid,len=nyjr) )
if (.not.allocated(latjr)) allocate(latjr(nyjr))
call nf90_err( nf90_inq_varid(ncid,'latitude',varid) )
call nf90_err( nf90_get_var(ncid,varid,latjr) )
dyjr=abs(latjr(2)-latjr(1))
! jrate field
call nf90_err( nf90_inq_varid(ncid,'jrate',varid) )
if (.not.allocated(jrate_average)) then
allocate(jrate_average(nxjr,nyjr,2))
endif
call nf90_err( nf90_get_var(ncid,varid,jrate_average(:,:,memid),start=(/1,1,mm/),count=(/nxjr,nyjr,1/)) )
! close file
call nf90_err( nf90_close(ncid) )
return
end subroutine readjrate
function photo_O1D(sza) result(jrate)
!*****************************************************************************
! *
! Calculates J(O1D) photolysis rate based on solar zenith angle *
! *
! Author: A. Stohl, Nov 2014 *
! Ref: Hamer, P. D., et al. Atmos. Chem. Phys. 2015
! *
!*****************************************************************************
! *
......@@ -898,38 +810,18 @@ module chemistry_mod
implicit none
real, intent(in) :: sza
real :: jrate
integer :: iz,ik
real :: z1,z2,zg,f1,f2,dummy
real :: photo_NO2
integer, parameter :: nzenith=11
real :: sun, jrate
real, parameter :: pi=3.1415927
real, dimension(nzenith) :: zangle,fact_photo
! zangle: zenith angles for which fact_photo is tabulated
! fact_photo: conversion of photolysis rate of NO2 to photolysis
! rate of O3 into O1D as a function of solar zenith angle
zangle=(/0.,10.,20.,30.,40.,50.,60.,70.,78.,86.,90.0001/)
fact_photo=(/0.4616E-02,0.4478E-02,0.4131E-02,0.3583E-02,0.2867E-02,&
&0.2081E-02,0.1235E-02,0.5392E-03,0.2200E-03,0.1302E-03,0.0902E-03/)
if (sza.lt.90.) then
do iz=1,nzenith-1
if(sza.ge.zangle(iz)) ik=iz
end do
z1=1./cos(zangle(ik)*pi/180.)
z2=1./cos(zangle(ik+1)*pi/180.)
zg=1./cos(sza*pi/180.)
dummy=(zg-z1)/(z2-z1)
f1=alog(fact_photo(ik))
f2=alog(fact_photo(ik+1))
photo_NO2=1.45e-2*exp(-0.4/cos(sza*pi/180.))
jrate=photo_NO2*exp(f1+(f2-f1)*dummy)
sun=sza*pi/180.
if (sun.le.(pi/2.)) then
sun=cos(sun)
else
jrate=0.
sun=0.
endif
jrate=(8.978e-5*(sun**1.436)*exp(-0.936/sun))
end function photo_O1D
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment