Skip to content

Commit

Permalink
Merge pull request #37 from NCAR/develop-35-ch4-quantum-yield
Browse files Browse the repository at this point in the history
Add TS1/TSMLT datasets, tests, and example configurations
  • Loading branch information
mattldawson authored Jan 19, 2024
2 parents 538f397 + a083881 commit dc042ed
Show file tree
Hide file tree
Showing 8 changed files with 436 additions and 11 deletions.
Binary file added data/cross_sections/CH3CL_JPL06.nc
Binary file not shown.
Binary file added data/cross_sections/SO2_Mills.nc
Binary file not shown.
45 changes: 34 additions & 11 deletions examples/ts1_tsmlt.json
Original file line number Diff line number Diff line change
Expand Up @@ -479,14 +479,37 @@
"name": "jacet",
"__reaction": "CH3COCH3 + hv -> CH3CO + CH3",
"cross section": {
"netcdf files": [
{ "file path": "data/cross_sections/CH3COCH3_1.nc" }
],
"type": "CH3COCH3+hv->CH3CO+CH3"
},
"quantum yield": {
"type": "CH3COCH3+hv->CH3CO+CH3"
}
"type": "temperature based",
"netcdf file": "data/cross_sections/CH3CL_JPL06.nc",
"parameterization": {
"AA": [ -299.80, 5.1047, -3.3630e-2, 9.5805e-5, -1.0135e-7 ],
"BB": [ -7.1727, 1.4837e-1, -1.1463e-3, 3.9188e-6, -4.9994e-9 ],
"lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ],
"minimum wavelength": 174.1,
"maximum wavelength": 216.0,
"base temperature": 273.0,
"base wavelength": 0.0,
"logarithm": "base 10",
"temperature ranges": [
{
"maximum": 209.999999999999,
"fixed value": 210.0
},
{
"minimum": 210,
"maximum": 300
},
{
"minimum": 300.00000000001,
"fixed value": 300.0
}
]
}
},
"quantum yield": {
"type": "base",
"constant value": 1.0
}
},
{
"name": "jmgly",
Expand Down Expand Up @@ -1506,7 +1529,7 @@
},
"quantum yield": {
"type": "base",
"constant value": 0.55
"constant value": 0.45
}
},
{
Expand All @@ -1520,7 +1543,7 @@
},
"quantum yield": {
"type": "base",
"constant value": 0.45
"constant value": 0.55
}
},
{
Expand Down Expand Up @@ -1627,7 +1650,7 @@
"cross section": {
"type": "base",
"netcdf files": [
{ "file path": "data/cross_sections/SO2_1.nc" }
{ "file path": "data/cross_sections/SO2_Mills.nc" }
]
},
"quantum yield": {
Expand Down
50 changes: 50 additions & 0 deletions test/data/xsqy.doug.config.json
Original file line number Diff line number Diff line change
Expand Up @@ -971,5 +971,55 @@
"__note": "second test: including edges of interpolation with relaxed tolerance",
"tolerance": 1.0e-3,
"mask": [ { "index": 79 } ]
},
{
"cross section": {
"type": "temperature based",
"netcdf file": "data/cross_sections/CH3CL_JPL06.nc",
"parameterization": {
"AA": [ -299.80, 5.1047, -3.3630e-2, 9.5805e-5, -1.0135e-7 ],
"BB": [ -7.1727, 1.4837e-1, -1.1463e-3, 3.9188e-6, -4.9994e-9 ],
"lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ],
"minimum wavelength": 174.1,
"maximum wavelength": 216.0,
"base temperature": 273.0,
"base wavelength": 0.0,
"logarithm": "base 10",
"temperature ranges": [
{
"maximum": 209.999999999999,
"fixed value": 210.0
},
{
"minimum": 210,
"maximum": 300
},
{
"minimum": 300.00000000001,
"fixed value": 300.0
}
]
}
},
"quantum yield": {
"type": "base",
"constant value": 1.0
},
"label": "CH3Cl + hv -> Cl",
"tolerance": 5.0e-3
},
{
"cross section": {
"type": "base",
"netcdf files": [
{ "file path": "data/cross_sections/SO2_Mills.nc" }
]
},
"quantum yield": {
"type": "base",
"constant value": 1.0
},
"label": "SO2 + hv -> SO + O",
"tolerance": 1.0e-4
}
]
2 changes: 2 additions & 0 deletions test/unit/tuv_doug/JCALC/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ target_sources(tuv_doug
XSQY_CFCL3.f
XSQY_CH2BR2.f
XSQY_CH3BR.f
XSQY_CH3CL.f
XSQY_CHBR3.f
XSQY_CL2O2.f
XSQY_CLO.f
Expand All @@ -24,6 +25,7 @@ target_sources(tuv_doug
XSQY_HNO3.f
XSQY_HO2NO2.f
XSQY_N2O5.f
XSQY_SO2.f
)

################################################################################
236 changes: 236 additions & 0 deletions test/unit/tuv_doug/JCALC/XSQY_CH3CL.f
Original file line number Diff line number Diff line change
@@ -0,0 +1,236 @@
subroutine XSQY_CH3CL(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn)
!-----------------------------------------------------------------------------!
! purpose: !
! provide product (cross section) x (quantum yield) for ch3cl photolysis: !
! ch3cl + hv -> products !
! cross section: from JPL06 recommendation !
! quantum yield: assumed to be unity !
!-----------------------------------------------------------------------------!
! parameters: !
! nw - integer, number of specified intervals + 1 in working (i) !
! wavelength grid !
! wl - real, vector of lower limits of wavelength intervals in (i) !
! working wavelength grid !
! wc - real, vector of center points of wavelength intervals in (i) !
! working wavelength grid !
! nz - integer, number of altitude levels in working altitude grid (i) !
! tlev - real, temperature (k) at each specified altitude level (i) !
! airlev - real, air density (molec/cc) at each altitude level (i) !
! j - integer, counter for number of weighting functions defined (io) !
! sq - real, cross section x quantum yield (cm^2) for each (o) !
! photolysis reaction defined, at each defined wavelength and !
! at each defined altitude level !
! jlabel - character*60, string identifier for each photolysis reaction (o) !
! defined !
!-----------------------------------------------------------------------------!
! edit history: !
! 07/30/07 Doug Kinnison !
!-----------------------------------------------------------------------------!
implicit none
include 'params'

!-----------------------------------------------------------------------------!
! ... input !
!-----------------------------------------------------------------------------!
real, intent(in) :: wl(kw)
real, intent(in) :: wc(kw)
real, intent(in) :: tlev(kz)
real, intent(in) :: airlev(kz)

integer, intent(in) :: nz
integer, intent(in) :: nw

character*80, intent(in) :: pn
character*60, intent(out) :: jlabel(kj)
real, intent(out) :: sq(kj,kz,kw)

!-----------------------------------------------------------------------------!
! ... input/output !
!-----------------------------------------------------------------------------!
integer, intent(inout) :: j

!-----------------------------------------------------------------------------!
! ... local !
!-----------------------------------------------------------------------------!
integer kdata
parameter(kdata=300)
integer i, iw, n, idum, nloop, n1
integer ierr, iz, iwc, icnt
real x1 (kdata), y1 (kdata)
real xin (kdata), yin (kdata)
real wctmp(kdata), wcb (kdata)
real ytmp (nz,kdata), ycomb(nz,kdata)
real ytd (nz,kw), tin (nz)
real AA(5), BB(5), lp(5)
real yg1 (kw)
real qy, ysave

AA(1) = -299.80
AA(2) = 5.1047
AA(3) = -3.3630e-2
AA(4) = 9.5805e-5
AA(5) = -1.0135e-7

BB(1) = -7.1727
BB(2) = 1.4837e-1
BB(3) = -1.1463e-3
BB(4) = 3.9188e-6
BB(5) = -4.9994e-9

lp(1) = 0.0
lp(2) = 1.0
lp(3) = 2.0
lp(4) = 3.0
lp(5) = 4.0

!----------------------------------------------
! ... tin set to tlev
!----------------------------------------------
tin(:) = tlev(:)

!----------------------------------------------
! ... jlabel(j) = 'CH3Cl + hv -> Cl'
!----------------------------------------------
j = j+1
jlabel(j) = 'CH3Cl + hv -> Cl'

!----------------------------------------------
! Derive temperature dependence
!----------------------------------------------
! Temperature dependence good between 210-300K
! and 174 nm-216 nm.
!----------------------------------------------
iwc = 1
ytmp(:,:)= 0.0

do iw = 1, nw-1

IF ((wc(iw) .GE. 174.) .AND. (wc(iw) .LE.216.)) THEN

do iz = 1, nz

IF (tin(iz) .LT. 210.) THEN
do nloop = 1, 5
ytmp(iz,iwc) = ytmp(iz,iwc)
& + AA(nloop)* (wc(iw)**lp(nloop))
& + (210.0-273.0)*BB(nloop)*wc(iw)**lp(nloop)
enddo
wctmp(iwc) = wc(iw)
ENDIF

IF ((tin(iz) .GE. 210.).AND.(tin(iz) .LE. 300.)) THEN
do nloop = 1,5

ytmp(iz,iwc) = ytmp(iz,iwc)
& + AA(nloop)* (wc(iw)**lp(nloop))
& + (tin(iz)-273.0)*BB(nloop)*wc(iw)**lp(nloop)
enddo
wctmp(iwc) = wc(iw)
ENDIF

IF (tin(iz) .GT. 300.) THEN
do nloop = 1, 5
ytmp(iz,iwc) = ytmp(iz,iwc)
& + AA(nloop)* (wc(iw)**lp(nloop))
& + (300.0-273.0)*BB(nloop)*wc(iw)**lp(nloop)
enddo
wctmp(iwc) = wc(iw)
ENDIF
enddo
iwc = iwc+ 1

ENDIF

enddo

!----------------------------------------------
! ... For wavelengths >216 nm and <174 nm
!----------------------------------------------
open(kin,file=TRIM(pn)//'XS_CH3CL_JPL06.txt',status='old')

read(kin,*) idum, n
do i = 1, idum-2
read(kin,*)
enddo

do i = 1, n
read(kin,*) xin(i), yin(i)
enddo
close(kin)

!----------------------------------------------
! ... Combine cross sections
!----------------------------------------------
do iz = 1, nz
icnt = 1

! ... < 174nm
do i = 1, n
IF (xin(i) .LT. 174.1) THEN
ycomb(iz,icnt) = yin(i)
wcb (icnt) = xin(i)
icnt = icnt + 1
ENDIF
enddo
! ... 174-216 nm
do i = 1, iwc-1
ycomb(iz,icnt) = 10**(ytmp(iz,i))
wcb (icnt) = wctmp(i)
icnt = icnt+1
enddo
! ... >216nm
do i = 1, n
IF (xin(i) .GT. 216.) THEN
ycomb(iz,icnt) = yin(i)
wcb (icnt) = xin(i)
icnt = icnt+1
ENDIF
enddo
enddo

!----------------------------------------------
! ... interpolate
!----------------------------------------------
do iz = 1, nz
n1 = icnt-1
y1 = ycomb(iz,:)
x1 = wcb
!----------------------------------------------
! Check routine
! do iw = 1, icnt-1
! print*, iw, wcb(iw), ycomb(iz,iw), tin(iz)
! enddo
! stop
!----------------------------------------------
call addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
call addpnt(x1,y1,kdata,n1, 0.,0.)
call addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
call addpnt(x1,y1,kdata,n1, 1e38,0.)
call inter2(nw,wl,yg1,n1,x1,y1,ierr)
ytd(iz,:) = yg1(:)

if (ierr .ne. 0) then
write(*,*) ierr, jlabel(j)
stop
endif
enddo

!----------------------------------------------
! Check routine
! iz = 1
! do iw = 19, 64
! print*, iw, wc(iw), ytd(iz,iw), tin(iz)
! enddo
! stop
!----------------------------------------------
!----------------------------------------------
! ...quantum yield assumed to be unity
!----------------------------------------------
qy = 1.
do iw = 1, nw-1
do iz = 1, nz
sq(j,iz,iw) = qy * ytd(iz,iw)
enddo
enddo

end subroutine XSQY_CH3CL
Loading

0 comments on commit dc042ed

Please sign in to comment.