diff --git a/data/cross_sections/BRO_JPL06.nc b/data/cross_sections/BRO_JPL06.nc new file mode 100644 index 00000000..5ad156c8 Binary files /dev/null and b/data/cross_sections/BRO_JPL06.nc differ diff --git a/data/cross_sections/CL2O2_JPL10.nc b/data/cross_sections/CL2O2_JPL10.nc new file mode 100644 index 00000000..1656cee4 Binary files /dev/null and b/data/cross_sections/CL2O2_JPL10.nc differ diff --git a/data/cross_sections/CLO_JPL06.nc b/data/cross_sections/CLO_JPL06.nc new file mode 100644 index 00000000..9b0c4281 Binary files /dev/null and b/data/cross_sections/CLO_JPL06.nc differ diff --git a/data/cross_sections/HNO3_JPL06.nc b/data/cross_sections/HNO3_JPL06.nc new file mode 100644 index 00000000..f56b8dd4 Binary files /dev/null and b/data/cross_sections/HNO3_JPL06.nc differ diff --git a/examples/ts1_tsmlt.json b/examples/ts1_tsmlt.json index b00663f3..fb668cc0 100644 --- a/examples/ts1_tsmlt.json +++ b/examples/ts1_tsmlt.json @@ -321,7 +321,7 @@ "__reaction": "HNO3 + hv -> OH + NO2", "cross section": { "netcdf files": [ - { "file path": "data/cross_sections/HNO3_1.nc" } + { "file path": "data/cross_sections/HNO3_JPL06.nc" } ], "type": "HNO3+hv->OH+NO2" }, @@ -555,11 +555,7 @@ "cross section": { "netcdf files": [ { - "file path": "data/cross_sections/BrO_1.nc", - "interpolator": { - "type": "fractional target", - "fold in": true - } + "file path": "data/cross_sections/BRO_JPL06.nc" } ], "type": "base" @@ -818,7 +814,7 @@ "__comments": "TODO - this doesn't exactly match the products in TS1", "cross section": { "netcdf files": [ - { "file path": "data/cross_sections/ClOOCl_1.nc" } + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } ], "type": "base" }, @@ -840,6 +836,20 @@ "type": "ClO+hv->Cl+O(1D)" } }, + { + "name": "jclo_o3p", + "__reaction": "ClO + hv -> Cl + O", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CLO_JPL06.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, { "name": "jclono2_a", "__reaction": "ClONO2 + hv -> Cl + NO3", diff --git a/test/data/xsqy.doug.config.json b/test/data/xsqy.doug.config.json index 901c9c2d..971775d8 100644 --- a/test/data/xsqy.doug.config.json +++ b/test/data/xsqy.doug.config.json @@ -109,5 +109,126 @@ }, "label": "CH2Br2 + hv -> 2Br", "tolerance": 1.0e-3 + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/BRO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "BRO + hv -> Br + O", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask" : [ { "index": 62 }, { "index": 86 }] + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/BRO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "BRO + hv -> Br + O", + "__note": "second test: including edges of interpolation with relaxed tolerance", + "tolerance": 5.0e-3 + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "Cl2O2 + hv -> Cl + ClOO", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask" : [ { "index": 34 }, { "index": 97 } ] + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "Cl2O2 + hv -> Cl + ClOO", + "__note": "second test: including edges of interpolation with relaxed tolerance", + "tolerance": 1.0e-3 + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CLO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "ClO + hv -> Cl + O", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask": [ { "index": 51 }, { "index": 71 }] + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CLO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "ClO + hv -> Cl + O", + "__note": "second test: including edges of interpolation with relaxed tolerance", + "tolerance": 1.0e-3 + }, + { + "cross section": { + "type": "HNO3+hv->OH+NO2", + "netcdf files": [ + { "file path": "data/cross_sections/HNO3_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "HNO3 + hv -> OH + NO2", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask": [ { "index": 30 }, { "index": 79 } ] + }, + { + "cross section": { + "type": "HNO3+hv->OH+NO2", + "netcdf files": [ + { "file path": "data/cross_sections/HNO3_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "HNO3 + hv -> OH + NO2", + "__note": "second test: including lower edge of interpolation with relaxed tolerance (upper edge is a very small value with large relative difference)", + "tolerance": 1.0e-3, + "mask": [ { "index": 79 } ] } ] diff --git a/test/unit/tuv_doug/JCALC/CMakeLists.txt b/test/unit/tuv_doug/JCALC/CMakeLists.txt index 93fd2604..6e70ac66 100644 --- a/test/unit/tuv_doug/JCALC/CMakeLists.txt +++ b/test/unit/tuv_doug/JCALC/CMakeLists.txt @@ -3,8 +3,12 @@ target_sources(tuv_doug PRIVATE + XSQY_BRO.f XSQY_CH2BR2.f + XSQY_CL2O2.f + XSQY_CLO.f XSQY_H2O.f + XSQY_HNO3.f XSQY_N2O5.f ) diff --git a/test/unit/tuv_doug/JCALC/XSQY_BRO.f b/test/unit/tuv_doug/JCALC/XSQY_BRO.f new file mode 100644 index 00000000..df054c6e --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_BRO.f @@ -0,0 +1,116 @@ + subroutine XSQY_BRO(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield): ! +! BrO + hv -> Br + O ! +! cross section: JPL06 ! +! quantum yield: is 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/27/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, ierr, iz + real x1(kdata) + real y1(kdata) + real yg(kw) + real qy + +!---------------------------------------------- +! ... jlabel(j) = 'BRO + hv -> Br + O' +!---------------------------------------------- + j = j+1 + jlabel(j) = 'BRO + hv -> Br + O' + +!---------------------------------------------------- +! ... cross sections from JPL06 recommendation +!---------------------------------------------------- +! ... 0.5nm resolution JPL06. + open(kin,file=TRIM(pn)//'XS_BRO_JPL06.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + close(kin) + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + call inter2(nw,wl,yg, n,x1, y1,ierr) + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + +!------------------------------------------------------- +! ... quantum yield (assumed) to be unity (JPL06) +!------------------------------------------------------- + qy = 1.0 + + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy * yg(iw) + + enddo + enddo + +!------------------------------------------------------- +! ... Check Routine +! do iw = 61, 87 +! print*, iw, wc(iw), (qy * yg(iw)) +! enddo +! stop +!------------------------------------------------------- + + end subroutine XSQY_BRO diff --git a/test/unit/tuv_doug/JCALC/XSQY_CL2O2.f b/test/unit/tuv_doug/JCALC/XSQY_CL2O2.f new file mode 100644 index 00000000..87abe738 --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_CL2O2.f @@ -0,0 +1,117 @@ + subroutine XSQY_CL2O2(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield) for hcl photolysis: ! +! cl2o2 + hv -> cl + cloo ! +! cross section: from JPL10 ! +! ! +!-----------------------------------------------------------------------------! +! 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: ! +! 01/13/2012 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=600) + integer i, iw, n, idum, ierr, iz + real x1(kdata) + real y1(kdata) + real yg(kw) + real qy + +!---------------------------------------------- +! ... jlabel(j) = 'cl2o2 -> cl + cloo' +!---------------------------------------------- + j = j+1 + jlabel(j) = 'Cl2O2 + hv -> Cl + ClOO' +! print*,jlabel(j) +!---------------------------------------------------- +! ... cross sections +!---------------------------------------------------- + open(kin,file= + $ TRIM(pn)//'XS_CL2O2_JPL10_500nm.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + close(kin) + +!---------------------------------------------------------- +! do i = 1, n +! print*, i, x1(i), y1(i) +! enddo +! stop +!---------------------------------------------------------- + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + call inter2(nw,wl,yg,n,x1,y1,ierr) + + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + +!-------------------------------------------------------------- +! ... quantum yield assumed to be 1.0 +!-------------------------------------------------------------- + qy = 1.0 + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy * yg(iw) + enddo + enddo +!---------------------------------------------------------- +! do iw = 28, 99 +! print*, iw, wc(iw), yg(iw) +! enddo +! stop +!---------------------------------------------------------- + end subroutine XSQY_CL2O2 diff --git a/test/unit/tuv_doug/JCALC/XSQY_CLO.f b/test/unit/tuv_doug/JCALC/XSQY_CLO.f new file mode 100644 index 00000000..2dab1cae --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_CLO.f @@ -0,0 +1,115 @@ + subroutine XSQY_CLO(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield): ! +! ClO + hv -> Cl + O ! +! cross section: JPL06 ! +! quantum yield: is 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/27/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, ierr, iz + real x1(kdata) + real y1(kdata) + real yg(kw) + real qy + +!---------------------------------------------- +! ... jlabel(j) = 'ClO + hv -> Cl + O' +!---------------------------------------------- + j = j+1 + jlabel(j) = 'ClO + hv -> Cl + O' + +!---------------------------------------------------- +! ... cross sections from JPL06 recommendation +!---------------------------------------------------- + open(kin,file=TRIM(pn)//'XS_CLO_JPL06.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + close(kin) + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + + call inter2(nw,wl,yg,n,x1,y1,ierr) + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif +!------------------------------------------------------- +! ... quantum yield (assumed) to be unity (JPL06) +!------------------------------------------------------- + qy = 1.0 + + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy * yg(iw) + + enddo + enddo +!------------------------------------------------------- +! ... Check routine (no temperature dependence +! print*,'jclo' +! do iw = 30, 72 +! print*, iw, wc(iw), (qy * yg(iw)) +! enddo +! stop +!------------------------------------------------------- + + end subroutine XSQY_CLO diff --git a/test/unit/tuv_doug/JCALC/XSQY_HNO3.f b/test/unit/tuv_doug/JCALC/XSQY_HNO3.f new file mode 100644 index 00000000..67792b40 --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_HNO3.f @@ -0,0 +1,140 @@ + subroutine XSQY_HNO3(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product of (cross section) x (quantum yield) for photolysis ! +! hno3 + hv -> oh + no2 ! +! cross section: burkholder et al., 1993 (and JPL06) ! +! 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: ! +! 05/98 original, adapted from former jspec1 subroutine ! +! 01/15/08 minor update,dek ! +!-----------------------------------------------------------------------------! + 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=100) + integer n1, n2 + integer i, iw, n, idum, iz + integer ierr + real x1 (kdata), x2 (kdata) + real y1 (kdata), y2 (kdata) + real yg1(kw), yg2(kw) + real yg( kw) + real tin(nz) + +!---------------------------------------------- +! ... tin set to tlev +!---------------------------------------------- + tin(:) = tlev(:) + +!---------------------------------------------- +! ... jlabel(j) = 'HNO3 -> OH + NO2 +!---------------------------------------------- + j = j + 1 + jlabel(j) = 'HNO3 + hv -> OH + NO2' + +!----------------------------------------------------------------------- +! ... hno3 cross section parameters from burkholder et al. 1993 +!----------------------------------------------------------------------- + open(kin,file=TRIM(pn)//'XS_HNO3_JPL06.txt',status='old') + +!... read lambda and cross sections + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + +!... read lambda and T-dep coeff. + read(kin,*) + do i = 1, n + read(kin,*) x2(i), y2(i) + enddo + close(kin) + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + call inter2(nw,wl,yg1,n,x1,y1,ierr) + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + + n= 80 + call addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) + call addpnt(x2,y2,kdata,n, 0.,0.) + call addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),0.) + call addpnt(x2,y2,kdata,n, 1.e+38,0.) + call inter2(nw,wl,yg2,n,x2,y2,ierr) + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + +!-------------------------------------------------- +! ... quantum yield = 1 +! correct for temperature dependence +!-------------------------------------------------- + do iw = 1, nw - 1 + do iz = 1, nz + sq(j,iz,iw) = yg1(iw) + $ * exp( yg2(iw)/1.e3*(tin(iz)-298.) ) + enddo + enddo + +!------------------------------------------------------- +! ... Check routine (no temperature dependence +! iz = 1 +! do iw = 29, 79 +! print*, iw, wc(iw), sq(j,iz,iw) +! enddo +! stop +!------------------------------------------------------- + + end subroutine XSQY_HNO3 diff --git a/test/unit/tuv_doug/data_sets.F90 b/test/unit/tuv_doug/data_sets.F90 index 434014ee..f62da19c 100644 --- a/test/unit/tuv_doug/data_sets.F90 +++ b/test/unit/tuv_doug/data_sets.F90 @@ -45,9 +45,10 @@ subroutine test_data_set( ) class(cross_section_t), pointer :: cross_section class(quantum_yield_t), pointer :: quantum_yield - character(len=*), parameter :: Iam = "H2O cross section test" + character(len=*), parameter :: Iam = "Doug's cross section tests" type(config_t) :: config, config_pair, cs_config, qy_config - class(iterator_t), pointer :: iter + type(config_t) :: mask_points_config, mask_point_config + class(iterator_t), pointer :: iter, mask_points_iter type(string_t) :: cs_type_name, qy_type_name, label character, allocatable :: buffer(:) integer :: pos, pack_size @@ -59,7 +60,9 @@ subroutine test_data_set( ) class(profile_t), pointer :: air, temperature class(grid_t), pointer :: wavelength real(kind=dk) :: tolerance + integer, allocatable :: mask_points(:) integer :: i + logical :: found ! Load grids based on Doug's TUV grids => get_grids( ) @@ -84,6 +87,22 @@ subroutine test_data_set( ) call config_pair%get( "label", label, Iam ) call config_pair%get( "tolerance", tolerance, Iam, & default = 1.0e-6_dk ) + call config_pair%get( "mask", mask_points_config, Iam, & + found = found ) + if( found ) then + mask_points_iter => mask_points_config%get_iterator( ) + allocate( mask_points( mask_points_config%number_of_children( ) ) ) + do i = 1, size( mask_points ) + call assert( 564855121, mask_points_iter%next( ) ) + call mask_points_config%get( mask_points_iter, mask_point_config, & + Iam ) + call mask_point_config%get( "index", mask_points( i ), Iam ) + end do + call assert( 888375064, .not. mask_points_iter%next( ) ) + deallocate( mask_points_iter ) + else + allocate( mask_points(0) ) + end if ! Load and test cross section if( musica_mpi_rank( comm ) == 0 ) then @@ -144,6 +163,8 @@ subroutine test_data_set( ) ! Skip first two bins because Lyman-Alpha bins are different in ! Doug's version of TUV-x. Data sets were adapted to have Lyman-Alpha ! specific data go into the TUV-x Lyman-Alpha bin 121.4-121.9 nm + ! Also skip any points explicitly masked in the configuration + tuvx_xsqy(:,mask_points(:)) = doug_xsqy(:,mask_points(:)) call check_values( 377150482, tuvx_xsqy(:,3:), & real( doug_xsqy(:,3:), kind=dk ), tolerance ) @@ -152,6 +173,7 @@ subroutine test_data_set( ) deallocate( cross_section_data ) deallocate( quantum_yield_data ) deallocate( tuvx_xsqy ) + deallocate( mask_points ) end do diff --git a/test/unit/tuv_doug/driver.F90 b/test/unit/tuv_doug/driver.F90 index cbdd8a7c..1dea6336 100644 --- a/test/unit/tuv_doug/driver.F90 +++ b/test/unit/tuv_doug/driver.F90 @@ -141,6 +141,18 @@ subroutine calculate( label, temperature, air_density, xsqy ) case( "CH2Br2 + hv -> 2Br" ) call XSQY_CH2BR2(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "BRO + hv -> Br + O" ) + call XSQY_BRO(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "Cl2O2 + hv -> Cl + ClOO" ) + call XSQY_CL2O2(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "ClO + hv -> Cl + O" ) + call XSQY_CLO(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case ("HNO3 + hv -> OH + NO2" ) + call XSQY_HNO3(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) case default call die( 946669022 ) end select diff --git a/tool/data_conversion/photo.config.json b/tool/data_conversion/photo.config.json new file mode 100644 index 00000000..5ee4685d --- /dev/null +++ b/tool/data_conversion/photo.config.json @@ -0,0 +1,44 @@ +{ + "photoreactions": [ + { + "molecule": "BRO", + "cross-sections": [ + { + "filespec": "XS_BRO_JPL06.txt", + "nPreSkip": 17, + "nRead": 198 + } + ] + }, + { + "molecule": "CL2O2", + "cross-sections": [ + { + "filespec": "XS_CL2O2_JPL10_500nm.txt", + "nPreSkip": 32, + "nRead": 521 + } + ] + }, + { + "molecule": "CLO", + "cross-sections": [ + { + "filespec": "XS_CLO_JPL06.txt", + "nPreSkip": 19, + "nRead": 36 + } + ] + }, + { + "molecule": "HNO3", + "cross-sections": [ + { + "filespec": "XS_HNO3_JPL06.txt", + "nPreSkip": 26, + "nRead": 80 + } + ] + } + ] +} diff --git a/tool/data_conversion/text_to_netcdf.py b/tool/data_conversion/text_to_netcdf.py new file mode 100644 index 00000000..bdf99226 --- /dev/null +++ b/tool/data_conversion/text_to_netcdf.py @@ -0,0 +1,61 @@ +#!/Library/Frameworks/Python.framework/Versions/3.8/bin/python3 + +import numpy as np +import sys +import json +from xsqy_subs import xform_to_netCDF + +#----------------------------------------------------- +# json config file is only possible argument +#----------------------------------------------------- +if( len(sys.argv) > 2 ): + print(f'\n{sys.argv[0]}: requires one or no arguments') + sys.exit( -1) +elif( len(sys.argv) == 2 ): + filespec = sys.argv[1] +else: + filespec = 'photo.config.json' + + +#----------------------------------------------------- +# open json photo config file +#----------------------------------------------------- +#ilespec = 'photo.config.tst.json' +try: + fp = open(filespec,'r') +except: + print(f"Failed to open {filespec}") + sys.exit(-1) + +#----------------------------------------------------- +# transfer config file into dictionary +#----------------------------------------------------- +try: + photDict = json.load(fp) +except: + print(f"Failed to load json file {filespec}") + sys.exit(-1) + +#----------------------------------------------------- +# done with json input file +#----------------------------------------------------- +fp.close() + +#----------------------------------------------------- +# loop through photo reactions in dictionary +#----------------------------------------------------- +list = photDict['photoreactions'] +molecule = "" +for rxt in list: +#----------------------------------------------------- +# call xform_to_netCDF +#----------------------------------------------------- + if( molecule != rxt['molecule']): + nFile = 1 + molecule = rxt['molecule'] + else: + nFile += 1 + xform_to_netCDF(nFile,rxt,'./') + +print('\n') +print(f'\nThere are {len(list)} photoreactions in {filespec}') diff --git a/tool/data_conversion/xsqy_subs.py b/tool/data_conversion/xsqy_subs.py new file mode 100644 index 00000000..d7f53404 --- /dev/null +++ b/tool/data_conversion/xsqy_subs.py @@ -0,0 +1,218 @@ +import numpy as np +import sys +from netCDF4 import Dataset +import netCDF4 as ncd +from datetime import datetime as dt + +""" +Function to read the data file(s) +""" +def read_data_file(data_dictionary,dataTray): + + InpFileSpec = data_dictionary['filespec'] + try: + InpFile = open(InpFileSpec,'r') + except: + print(f'Failed to open data file {InpFileSpec}') + sys.exit(-3) + + print(f'Opened data file {InpFileSpec}') + nLines = len(InpFile.readlines()) + InpFile.seek(0) + nskipHdr = data_dictionary['nPreSkip'] if 'nPreSkip' in data_dictionary else 0 + nRead = data_dictionary['nRead'] if 'nRead' in data_dictionary else 0 + + header = '' +# if header lines exist then read them + if( nskipHdr > 0 ): + for ndx in range(nskipHdr): + header += InpFile.readline() + InpFile.seek(0) + + nskipHdr = abs(nskipHdr) + nskipEnd = nLines - (nskipHdr + nRead) + print(f'nLines,nskipHdr,nRead,nskipEnd = {nLines},{nskipHdr},{nRead},{nskipEnd}') + try: + data = np.genfromtxt(InpFile,dtype='float64',skip_header=nskipHdr,skip_footer=nskipEnd,comments=None) + print(f'Read cross section file {InpFileSpec}') + except: + print(f'Failed to read data file {InpFileSpec}') + sys.exit(-2) + + try: + dataTray.append(data) + except: + print('Failed to append data array to dataTray') + sys.exit(-2) + + InpFile.close() + print(f'Closed data file {InpFileSpec}') + return(header) + +""" +Function to write the netCDF file +""" +def stuff_netCDF_file(ncFile,interpolationTemps,dataTray,hasLambdaGrid,InpFileSpecs,var_typ,Headers): + + ndataVars = len(dataTray) + print(f'ndataVars = {ndataVars}') + + for dataVarNdx in range(ndataVars): + nparameterRow,nparameterCol = np.shape(dataTray[dataVarNdx]) + if( hasLambdaGrid ): + nparameterCol -= 1 + ntemperature = min( len(interpolationTemps),nparameterCol ) + print(f'data array is ({nparameterRow},{nparameterCol})') + DataTag = var_typ + "_parameters" +# define dimensions + RowDimName = 'bins' + ColDimName = 'parameters' + TempDimName = 'temperatures' + print(f'Variable type = {var_typ}') + print(f'RowDimName,ColDimName = {RowDimName},{ColDimName}') + ncFile.createDimension(RowDimName,nparameterRow) + ncFile.createDimension(ColDimName,nparameterCol) + ncFile.createDimension(TempDimName,ntemperature) +# create wavelength grid + if( hasLambdaGrid ): + Var = ncFile.createVariable('wavelength',np.float64,(RowDimName)) + Var.units = 'nm' +# write wavelength grid + Var[:] = dataTray[dataVarNdx][:,0] +# create interpolation temperature array + if( len(interpolationTemps) > 0 ): + Var = ncFile.createVariable('temperature',np.float64,(TempDimName)) + Var.units = 'K' +# write interpolation temperature array + Var[:] = interpolationTemps[:ntemperature] +# create cross section or quantum yield data array + Var = ncFile.createVariable(DataTag,np.float64,(ColDimName,RowDimName)) + Var.hdr = Headers[dataVarNdx] +# write data array + if( hasLambdaGrid ): + Var[:,:] = np.transpose(dataTray[dataVarNdx][:,1:]) + else: + Var[:,:] = np.transpose(dataTray[dataVarNdx]) + if( var_typ == 'cross_section' ): + if( hasLambdaGrid ): + Var.units = 'cm^2' + else: + Var.units = 'see source code' + else: + Var.units = 'fraction' + +# global attributes + version = '1.0' + ncFile.Author = 'TUV Data Xformer ' + version + now = dt.now() + ncFile.created = now.strftime("%Y-%m-%d %H:%M:%S") + if( var_typ == 'cross_section' ): + ncFile.title = 'Cross section parameters' + else: + ncFile.title = 'Quantum yield parameters' + ncFile.file = InpFileSpecs + +""" +Transform ascii data file(s) to netCDF counterpart +""" +def xform_to_netCDF(nFile,phtDictionary,ncd_path): + +# cross section + if( 'cross-sections' in phtDictionary ): +# form netCDF file for cross sections + molecule = phtDictionary['molecule'] + + print(f'\nProcessing {molecule} cross section') + xsects = phtDictionary['cross-sections'] + nxsects = len(xsects) + print(f' There are {nxsects} cross section files') + dataTray = [] + interpolationTemps = [] + Headers = [] + +# loop over ascii input data files + for xsect in xsects: + ncdFilespec = ncd_path + '/' + molecule + '_cross_section_' + str(nFile) + '.nc' +# create the netcdf dataset + print(f'\nCreating netCDF file {ncdFilespec}') + + try: + ncFile = Dataset(ncdFilespec,mode='w',format='NETCDF4_CLASSIC') + except: + print(f'Failed to create netCDF4 dataset {ncdFilespec}') + sys.exit(-1) + +# 1st data column wavelength grid? + if( 'has lambda grid' in xsect ): + hasLambdaGrid = xsect['has lambda grid'] + else: + hasLambdaGrid = True + +# interpolation temperatures? + if( 'interpolation temperature' in xsect ): + interpolationTemps = xsect['interpolation temperature'] + print("\nInterpolation temperatures:") + print(interpolationTemps) + + header = '' + header = read_data_file(xsect,dataTray) + Headers.append(header) + + InpFileSpecs = xsect['filespec'] + + print(f'\nThere are {len(dataTray)} arrays in dataTray') + print(f'\nShape data array is {dataTray[0].shape}') + stuff_netCDF_file(ncFile,interpolationTemps,dataTray,hasLambdaGrid,InpFileSpecs,'cross_section',Headers) +# close netcdf file + ncFile.close() + + print('\n') + +# quantum yield + if( 'quantum-yields' in phtDictionary ): +# form netCDF file for cross sections + molecule = phtDictionary['molecule'] + + print(f'\nProcessing {molecule} quantum yield') + qylds = phtDictionary['quantum-yields'] + nqylds = len(qylds) + print(f' There are {nqylds} quantum yield files') + dataTray = [] + interpolationTemps = [] + Headers = [] + +# loop over ascii input data files + for qyld in qylds: + ncdFilespec = ncd_path + '/' + molecule + '_quantum_yield_' + str(nFile) + '.nc' +# create the netcdf dataset + print(f'\nCreating netCDF file {ncdFilespec}') + + try: + ncFile = Dataset(ncdFilespec,mode='w',format='NETCDF4_CLASSIC') + except: + print(f'Failed to create netCDF4 dataset {ncdFilespec}') + sys.exit(-1) + +# 1st data column wavelength grid? + if( 'has lambda grid' in qyld ): + hasLambdaGrid = qyld['has lambda grid'] + else: + hasLambdaGrid = True + +# interpolation temperatures? + if( 'interpolation temperature' in qyld ): + interpolationTemps = qyld['interpolation temperature'] + print("\nInterpolation temperatures:") + print(interpolationTemps) + + header = '' + header = read_data_file(qyld,dataTray) + Headers.append(header) + + InpFileSpecs = qyld['filespec'] + + print(f'\nThere are {len(dataTray)} arrays in dataTray') + print(f'\nShape data array is {dataTray[0].shape}') + stuff_netCDF_file(ncFile,interpolationTemps,dataTray,hasLambdaGrid,InpFileSpecs,'quantum_yield',Headers) +# close netcdf file + ncFile.close()