diff --git a/.gitignore b/.gitignore index ae916806c6..130f252e2c 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ src/dynamics/fv3/atmos_cubed_sphere libraries/FMS libraries/mct libraries/parallelio +libraries/tuv-x src/atmos_phys src/dynamics/mpas/dycore share @@ -23,6 +24,9 @@ src/hemco buildnmlc buildcppc +# ignore ide configs +.vscode + # Ignore editor temporaries and backups *~ .#* diff --git a/Externals_CAM.cfg b/Externals_CAM.cfg index 03698dc3f6..8de90f77a0 100644 --- a/Externals_CAM.cfg +++ b/Externals_CAM.cfg @@ -71,6 +71,13 @@ sparse = ../.mpas_sparse_checkout hash = b8c33daa required = True +[TUV-x] +local_path = libraries/tuv-x +protocol = git +repo_url = https://github.com/NCAR/tuv-x.git +hash = ccda76c0553064a9b7b0eba73162ddeee6d8eaff +required = True + [hemco] local_path = src/hemco tag = hemco-cesm1_2_1_hemco3_6_3_cesm diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 286bf4f953..912c2848b2 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -7352,6 +7352,17 @@ Maximum zenith angle (degrees) used for photolysis. Default: set by build-namelist. + +Filepath of TUV-X configuration specification. +Default: NONE + + + +Switch to turn on TUV-X photolysis. +Default: FALSE + diff --git a/cime_config/buildlib b/cime_config/buildlib index 73db5db3dd..0b125e032b 100755 --- a/cime_config/buildlib +++ b/cime_config/buildlib @@ -7,7 +7,7 @@ create the cam library # pylint: disable=unused-wildcard-import, bad-whitespace, too-many-locals # pylint: disable=invalid-name import sys, os, filecmp, shutil, imp - +from glob import glob _CIMEROOT = os.environ.get("CIMEROOT") if _CIMEROOT is None: @@ -21,6 +21,7 @@ from CIME.case import Case from CIME.utils import run_sub_or_cmd, expect, run_cmd from CIME.buildlib import parse_input from CIME.build import get_standard_makefile_args +from CIME.XML.env_build import EnvBuild logger = logging.getLogger(__name__) @@ -113,12 +114,151 @@ def _build_cam(caseroot, libroot, bldroot): logger.info("%s: \n\n output:\n %s \n\n err:\n\n%s\n", cmd, out, err) expect(rc == 0, "Command %s failed with rc=%s" % (cmd, rc)) +############################################################################### +def _run_cmd(command, working_dir): +############################################################################### + + rc, out, err = run_cmd(command, from_dir=working_dir, verbose=True) + expect(rc == 0, "Command {} failed with rc={}".format(command, rc)) + +############################################################################### +def _cmake_default_args(caseroot): +############################################################################### +# Returns a dictionary of CMake variables based on the Macros.cmake file for +# the build. + + build = EnvBuild(case_root=caseroot) + with Case(caseroot) as case: + macro_path = os.path.abspath(os.path.join(caseroot, "cmake_macros", "")) + args = "-DCONVERT_TO_MAKE=ON " + args += "-DCASEROOT={} ".format(caseroot) + args += "-DCOMPILER={} ".format(build.get_value("COMPILER")) + args += "-DOS={} ".format(build.get_value("OS")) + args += "-DMACH={} ".format(case.get_value("MACH")) + args += "-DCMAKE_C_COMPILER_WORKS=1 " + args += "-DCMAKE_Fortran_COMPILER_WORKS=1 " + args += "-DCMAKE_CXX_COMPILER_WORKS=1 " + args += "-DDEBUG={} ".format(build.get_value("DEBUG")) + cmd = "cmake {} .".format(args) + rc, out, err = run_cmd(cmd, combine_output=True, from_dir=macro_path) + expect(rc == 0, "Command {} failed with rc={} out={} err={}".format(cmd, rc, out, err)) + + arg_dict = {} + for line in out.splitlines(): + if ":=" in line: + key, val = line.split(":=") + arg_dict[key.replace('CIME_SET_MAKEFILE_VAR','').strip()] = val.strip() + + return arg_dict + +############################################################################### +def _build_tuvx(caseroot, libroot, bldroot): +############################################################################### +# Builds the TUV-x library and updates the case variables used to set the +# include paths and linked libraries + + build = EnvBuild(case_root=caseroot) + with Case(caseroot) as case: + bldpath = os.path.join(bldroot, "tuv-x") + if not os.path.exists(bldpath): + os.makedirs(bldpath) + srcpath = os.path.abspath(os.path.join(case.get_value("COMP_ROOT_DIR_ATM"), \ + "libraries", "tuv-x", "")) + logger.info("Building TUV-x in {} from source in {}\n".format(bldpath, srcpath)) + + arg_dict = _cmake_default_args(caseroot) + if build.get_value("MPILIB") == "mpi-serial": + cmake_args = "-DCMAKE_Fortran_COMPILER={} ".format(arg_dict["SFC"]) + else: + cmake_args = "-DCMAKE_Fortran_COMPILER={} ".format(arg_dict["MPIFC"]) + cmake_args += "-DTUVX_ENABLE_MPI:BOOL=TRUE " + if case.get_value("DEBUG"): + cmake_args += "-DCMAKE_BUILD_TYPE=Debug " + else: + cmake_args += "-DCMAKE_BUILD_TYPE=Release " + cmake_args += "-DCMAKE_C_COMPILER_WORKS=1 " + cmake_args += "-DCMAKE_CXX_COMPILER_WORKS=1 " + cmake_args += "-DTUVX_ENABLE_TESTS=OFF " + cmake_args += "-DTUVX_ENABLE_COVERAGE=OFF " + cmake_args += "-DCMAKE_Fortran_FLAGS='{}' ".format(arg_dict["FFLAGS"]) + cmake_args += "-DCMAKE_INSTALL_PREFIX='{}' ".format(libroot) + cmake_args += "-DTUVX_INSTALL_INCLUDE_DIR='{}' ".format(_tuvx_include_dir(libroot)) + cmake_args += "-DTUVX_INSTALL_MOD_DIR='{}' ".format(_tuvx_include_dir(libroot)) + cmake_args += srcpath + + _run_cmd("cmake {}".format(cmake_args), bldpath) + _run_cmd(case.get_value("GMAKE"), bldpath) + _run_cmd("{} install".format(case.get_value("GMAKE")), bldpath) + + # add TUV-x to include paths + incldir = os.environ.get('USER_INCLDIR') + if incldir is None: + incldir = '' + os.environ['USER_INCLDIR'] = incldir + \ + " -I{} ".format(_tuvx_include_dir(libroot)) + + # create simlink to library in folder CIME expects libraries to be in + dst = os.path.join(libroot, "libtuvx.a") + if os.path.isfile(dst): + os.remove(dst) + os.symlink(_tuvx_lib_path(libroot), dst) + +############################################################################### +def _tuvx_include_dir(libroot): +############################################################################### +# Returns the path to the TUV-x include directory + + coreinc = os.path.join(libroot, "include", "") + expect(os.path.exists(coreinc), \ + "TUV-x include directory not found at {}".format(coreinc)) + return coreinc + +############################################################################### +def _tuvx_lib_path(libroot): +############################################################################### +# Returns the path to the TUV-x library + + corelib = os.path.join(_tuvx_install_dir(libroot), "lib64", "libtuvx.a") + if not os.path.exists(corelib): + corelib = os.path.join(_tuvx_install_dir(libroot), "lib", "libtuvx.a") + expect(os.path.exists(corelib), \ + "TUV-x library not found at {}".format(corelib)) + return corelib + +############################################################################### +def _tuvx_install_dir(libroot): +############################################################################### +# Returns the path to the TUV-x install directory + + corepaths = glob(os.path.join(libroot, "tuvx*")) + expect(len(corepaths)>0, \ + "TUV-x not found at {}".format(libroot)) + expect(len(corepaths)<2, \ + "Multiple TUV-x versions found at {}".format(libroot)) + expect(os.path.exists(corepaths[0]), \ + "TUV-x install directory not found at {}".format(corepaths[0])) + return corepaths[0] + +############################################################################### +def _tuvx_package_dir(libroot): +############################################################################### +# Returns the path to the TUV-x CMake package + + paths = glob(os.path.join(libroot, "tuvx*", "cmake", "tuvx*" )) + expect(len(paths)>0, \ + "TUV-x package not found at {}".format(libroot)) + expect(len(paths)<2, \ + "Multiple TUV-x versions found at {}".format(libroot)) + expect(os.path.exists(paths[0]), \ + "TUV-x package directory not found at {}".format(paths[0])) + return paths[0] ############################################################################### def _main_func(): caseroot, libroot, bldroot = parse_input(sys.argv) + _build_tuvx(caseroot, libroot, bldroot) _build_cam(caseroot, libroot, bldroot) diff --git a/cime_config/buildnml b/cime_config/buildnml index 28e3e8198c..bebcd94426 100755 --- a/cime_config/buildnml +++ b/cime_config/buildnml @@ -211,6 +211,17 @@ def buildnml(case, caseroot, compname): if (os.path.isfile(file1)) and (not os.path.isfile(file2)): shutil.copy(file1,file2) + # Temporary copy of TUV-x data for development + dest_data = os.path.join(rundir, "data") + if os.path.exists(dest_data): + shutil.rmtree(dest_data) + shutil.copytree(os.path.join(srcroot, "libraries", "tuv-x", "data"), \ + dest_data) + shutil.copy2(os.path.join(srcroot, "cime_config", "tuvx_MOZART.json"), \ + os.path.join(rundir, "tuvx_MOZART.json")) + shutil.copy2(os.path.join(srcroot, "cime_config", "tuvx_MOZART_TS1.json"), \ + os.path.join(rundir, "tuvx_MOZART_TS1.json")) + ############################################################################### def _main_func(): diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index 0dc9a4d10b..fc1ca74b8b 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -203,6 +203,18 @@ + + char + + -ltuvx -lstdc++ + build_component_cam + env_build.xml + + CAM linked libraries. The libraries are built by CAM's buildlib script and should be included + during linking of the model executable. + + + char diff --git a/cime_config/tuvx_MOZART.json b/cime_config/tuvx_MOZART.json new file mode 100644 index 0000000000..8129266467 --- /dev/null +++ b/cime_config/tuvx_MOZART.json @@ -0,0 +1,540 @@ +{ + "__description": "TUV-x configuration for the MOZART chemical mechanism", + "O2 absorption" : { + "cross section parameters file": "data/cross_sections/O2_parameters.txt" + }, + "grids": [ + ], + "profiles": [ + ], + "radiative transfer": { + "solver": { + "type": "discrete ordinate", + "number of streams": 4 + }, + "cross sections": [ + { + "name": "air", + "type": "air" + }, + { + "name": "O3", + "netcdf files": [ + { "file path": "data/cross_sections/O3_1.nc" }, + { "file path": "data/cross_sections/O3_2.nc" }, + { "file path": "data/cross_sections/O3_3.nc" }, + { "file path": "data/cross_sections/O3_4.nc" } + ], + "type": "O3" + }, + { + "name": "O2", + "netcdf files": [ + { + "file path": "data/cross_sections/O2_1.nc", + "lower extrapolation": { "type": "boundary" }, + "interpolator": { "type": "fractional target" } + } + ], + "type": "base" + } + ], + "radiators": [ + { + "name": "air", + "type": "base", + "treat as air": true, + "cross section": "air", + "vertical profile": "air", + "vertical profile units": "molecule cm-3" + }, + { + "name": "O2", + "type": "base", + "cross section": "O2", + "vertical profile": "O2", + "vertical profile units": "molecule cm-3" + }, + { + "name": "O3", + "type": "base", + "cross section": "O3", + "vertical profile": "O3", + "vertical profile units": "molecule cm-3" + } + ] + }, + "photolysis": { + "reactions": [ + { + "name": "jo2", + "__reaction": "O2 + hv -> O + O", + "cross section": { + "netcdf files": [ + { + "file path": "data/cross_sections/O2_1.nc", + "lower extrapolation": { "type": "boundary" }, + "interpolator": { "type": "fractional target" } + } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jo1d", + "__reaction": "O3 + hv -> O2 + O(1D)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/O3_1.nc" }, + { "file path": "data/cross_sections/O3_2.nc" }, + { "file path": "data/cross_sections/O3_3.nc" }, + { "file path": "data/cross_sections/O3_4.nc" } + ], + "type": "O3" + }, + "quantum yield": { + "type": "O3+hv->O2+O(1D)" + } + }, + { + "name": "jo3p", + "__reaction": "O3 + hv -> O2 + O(3P)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/O3_1.nc" }, + { "file path": "data/cross_sections/O3_2.nc" }, + { "file path": "data/cross_sections/O3_3.nc" }, + { "file path": "data/cross_sections/O3_4.nc" } + ], + "type": "O3" + }, + "quantum yield": { + "type": "O3+hv->O2+O(3P)" + } + }, + { + "name": "jn2o", + "__reaction": "N2O + hv -> N2 + O(1D)", + "cross section": { + "type": "N2O+hv->N2+O(1D)" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jno2", + "__reaction": "NO2 + hv -> NO + O(3P)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/NO2_1.nc" } + ], + "type": "NO2 tint" + }, + "quantum yield": { + "netcdf files": ["data/quantum_yields/NO2_1.nc"], + "type": "NO2 tint", + "lower extrapolation": { "type": "boundary" } + } + }, + { + "name": "jn2o5", + "__reaction": "N2O5 + hv -> NO2 + NO3", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/N2O5_1.nc" }, + { "file path": "data/cross_sections/N2O5_2.nc" } + ], + "type": "N2O5+hv->NO2+NO3" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhno3", + "__reaction": "HNO3 + hv -> OH + NO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HNO3_1.nc" } + ], + "type": "HNO3+hv->OH+NO2" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jno3_a", + "__reaction": "NO3 + hv -> NO2 + O(3P)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/NO3_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/NO3-NO2+O(3P)_1.nc" + ], + "type": "tint", + "lower extrapolation": { + "type": "constant", + "value": 1.0 + } + } + }, + { + "name": "jno3_b", + "__reaction": "NO3 + hv -> NO + O2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/NO3_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/NO3-NO+O2_1.nc" + ], + "type": "tint" + } + }, + { + "name": "jch3ooh", + "__reaction": "CH3OOH + hv -> CH3O + OH", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3OOH_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jch2o_a", + "__reaction": "CH2O + hv -> H + HCO", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH2O_1.nc" } + ], + "type": "CH2O" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CH2O_1.nc" + ], + "type": "base", + "lower extrapolation": { + "type": "boundary" + } + } + }, + { + "name": "jch2o_b", + "__reaction": "CH2O + hv -> H2 + CO", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH2O_1.nc" } + ], + "type": "CH2O" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CH2O_1.nc" + ], + "type": "CH2O", + "lower extrapolation": { + "type": "boundary" + } + } + }, + { + "name": "jh2o2", + "__reaction": "H2O2 + hv -> OH + OH", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/H2O2_1.nc" } + ], + "type": "H2O2+hv->OH+OH" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jch3cho", + "__reaction": "CH3CHO + hv -> CH3 + HCO", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3CHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CH3CHO_1.nc" + ], + "type": "CH3CHO+hv->CH3+HCO" + } + }, + { + "name": "jpan", + "__reaction": "PAN + hv -> 0.6*CH3CO3 + 0.6*NO2 + 0.4*CH3O2 + 0.4*NO3 + 0.4*CO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/PAN_1.nc" } + ], + "type": "CH3ONO2+hv->CH3O+NO2" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jmvk", + "__reaction": "MVK + hv -> 0.7*C3H6 + 0.7*CO + 0.3*CH3O2 + 0.3*CH3CO3", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/MVK_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "MVK+hv->Products" + } + }, + { + "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" + } + }, + { + "name": "jmgly", + "__reaction": "CH3COCHO + hv -> CH3CO3 + CO + HO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3COCHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "CH3COCHO+hv->CH3CO+HCO" + } + }, + { + "name": "jglyald", + "__reaction": "GLYALD + hv -> 2*HO2 + CO + CH2O", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HOCH2CHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jglyoxal", + "__reaction": "GLYOXAL + hv -> 2*CO + 2*HO2", + "__comments": "TODO the products of this reaction don't exactly match", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CHOCHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CHOCHO-H2_CO_CO_1.nc" + ], + "type": "base", + "lower extrapolation": { + "type": "boundary" + } + } + }, + { + "name": "jho2no2_a", + "__reaction": "HNO4 + hv -> OH + NO3", + "__comments": "TODO Doug's data sets have special temperature dependence - need new type?", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HNO4_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.2 + } + }, + { + "name": "jho2no2_b", + "__reaction": "HNO4 + hv -> HO2 + NO2", + "__comments": "TODO Doug's data sets have special temperature dependence - need new type?", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HNO4_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.8 + } + }, + { + "name": "jmacr_a", + "__reaction": "CH2=C(CH3)CHO->1.34HO2+0.66MCO3+1.34CH2O+CH3CO3", + "__comments": "Methacrolein photolysis channel 1", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/Methacrolein_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.005 + } + }, + { + "name": "jmacr_b", + "__reaction": "CH2=C(CH3)CHO->0.66OH+1.34CO", + "__comments": "Methacrolein photolysis channel 2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/Methacrolein_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.005 + } + }, + { + "name": "jhyac", + "__reaction": "CH2(OH)COCH3->CH3CO3+HO2+CH2O", + "__comments": "hydroxy acetone", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/Hydroxyacetone_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.65 + } + } + ] + }, + "__CAM options": { + "aliasing": { + "default matching": "backup", + "pairs": [ + { + "to": "jpooh", + "__reaction": "POOH (C3H6OHOOH) + hv -> CH3CHO + CH2O + HO2 + OH", + "from": "jch3ooh" + }, + { + "to": "jch3co3h", + "__reaction": "CH3COOOH + hv -> CH3O2 + OH + CO2", + "from": "jh2o2", + "scale by": 0.28 + }, + { + "to": "jmpan", + "__reaction": "MPAN + hv -> MCO3 + NO2", + "from": "jpan" + }, + { + "to": "jc2h5ooh", + "__reaction": "C2H5OOH + hv -> CH3CHO + HO2 + OH", + "from": "jch3ooh" + }, + { + "to": "jc3h7ooh", + "__reaction": "C3H7OOH + hv -> 0.82*CH3COCH3 + OH + HO2", + "from": "jch3ooh" + }, + { + "to": "jrooh", + "__reaction": "ROOH + hv -> CH3CO3 + CH2O + OH", + "from": "jch3ooh" + }, + { + "to": "jxooh", + "__reaction": "XOOH + hv -> OH", + "from": "jch3ooh" + }, + { + "to": "jonitr", + "__reaction": "ONITR + hv -> NO2", + "from": "jch3cho" + }, + { + "to": "jisopooh", + "__reaction": "ISOPOOH + hv -> 0.402*MVK + 0.288*MACR + 0.69*CH2O + HO2", + "from": "jch3ooh" + }, + { + "to": "jmek", + "__reaction": "MEK + hv -> CH3CO3 + C2H5O2", + "from": "jacet" + }, + { + "to": "jbigald", + "__reaction": "BIGALD + hv -> .45*CO + .13*GLYOXAL + .56*HO2 + .13*CH3CO3 + .18*CH3COCHO", + "from": "jno2", + "scale by": 0.2 + }, + { + "to": "jalkooh", + "__reaction": "ALKOOH + hv -> .4*CH3CHO + .1*CH2O + .25*CH3COCH3 + .9*HO2 + .8*MEK + OH", + "from": "jch3ooh" + }, + { + "to": "jmekooh", + "__reaction": "MEKOOH + hv -> OH + CH3CO3 + CH3CHO", + "from": "jch3ooh" + }, + { + "to": "jtolooh", + "__reaction": "TOLOOH + hv -> OH + .45*GLYOXAL + .45*CH3COCHO + .9*BIGALD", + "from": "jch3ooh" + }, + { + "to": "jterpooh", + "__reaction": "TERPOOH + hv -> OH + .1*CH3COCH3 + HO2 + MVK + MACR", + "from": "jch3ooh" + } + ] + } + } +} diff --git a/cime_config/tuvx_MOZART_TS1.json b/cime_config/tuvx_MOZART_TS1.json new file mode 100644 index 0000000000..82b5d0175c --- /dev/null +++ b/cime_config/tuvx_MOZART_TS1.json @@ -0,0 +1,2006 @@ +{ + "__description": "TUV-x configuration for the MOZART-TS1 and MOZART-TSMLT chemical mechanisms", + "O2 absorption" : { + "cross section parameters file": "data/cross_sections/O2_parameters.txt" + }, + "grids": [ + ], + "profiles": [ + ], + "radiative transfer": { + "solver": { + "type": "discrete ordinate", + "number of streams": 4 + }, + "cross sections": [ + { + "name": "air", + "type": "air" + }, + { + "name": "O3", + "netcdf files": [ + { "file path": "data/cross_sections/O3_1.nc" }, + { "file path": "data/cross_sections/O3_2.nc" }, + { "file path": "data/cross_sections/O3_3.nc" }, + { "file path": "data/cross_sections/O3_4.nc" } + ], + "type": "O3" + }, + { + "name": "O2", + "netcdf files": [ + { + "file path": "data/cross_sections/O2_1.nc", + "lower extrapolation": { "type": "boundary" }, + "interpolator": { "type": "fractional target" } + } + ], + "type": "base" + } + ], + "radiators": [ + { + "name": "air", + "type": "base", + "treat as air": true, + "cross section": "air", + "vertical profile": "air", + "vertical profile units": "molecule cm-3" + }, + { + "name": "O2", + "type": "base", + "cross section": "O2", + "vertical profile": "O2", + "vertical profile units": "molecule cm-3" + }, + { + "name": "O3", + "type": "base", + "cross section": "O3", + "vertical profile": "O3", + "vertical profile units": "molecule cm-3" + } + ] + }, + "photolysis": { + "reactions": [ + { + "name": "jo2_a", + "__reaction": "O2 + hv -> O + O1D", + "cross section": { + "apply O2 bands": true, + "netcdf files": [ + { + "file path": "data/cross_sections/O2_1.nc", + "lower extrapolation": { "type": "boundary" }, + "interpolator": { "type": "fractional target" } + } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0, + "override bands": [ + { + "band": "lyman-alpha", + "value": 0.53 + }, + { + "band": "schumann-runge continuum", + "value": 1.0 + } + ] + }, + "heating" : { + "energy term": 175.05 + } + }, + { + "name": "jo2_b", + "__reaction": "O2 + hv -> O + O", + "cross section": { + "apply O2 bands": true, + "netcdf files": [ + { + "file path": "data/cross_sections/O2_1.nc", + "lower extrapolation": { "type": "boundary" }, + "interpolator": { "type": "fractional target" } + } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0, + "override bands": [ + { + "band": "lyman-alpha", + "value": 0.47 + }, + { + "band": "schumann-runge continuum", + "value": 0.0 + } + ] + }, + "heating" : { + "energy term": 242.37 + } + }, + { + "name": "jo3_a", + "__reaction": "O3 + hv -> O2 + O(1D)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/O3_1.nc" }, + { "file path": "data/cross_sections/O3_2.nc" }, + { "file path": "data/cross_sections/O3_3.nc" }, + { "file path": "data/cross_sections/O3_4.nc" } + ], + "type": "O3" + }, + "quantum yield": { + "type": "O3+hv->O2+O(1D)" + }, + "heating" : { + "energy term": 310.32 + } + }, + { + "name": "jo3_b", + "__reaction": "O3 + hv -> O2 + O(3P)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/O3_1.nc" }, + { "file path": "data/cross_sections/O3_2.nc" }, + { "file path": "data/cross_sections/O3_3.nc" }, + { "file path": "data/cross_sections/O3_4.nc" } + ], + "type": "O3" + }, + "quantum yield": { + "type": "O3+hv->O2+O(3P)" + }, + "heating" : { + "energy term": 1179.87 + } + }, + { + "name": "jn2o", + "__reaction": "N2O + hv -> N2 + O(1D)", + "cross section": { + "type": "N2O+hv->N2+O(1D)" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jno2", + "__reaction": "NO2 + hv -> NO + O(3P)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/NO2_1.nc" } + ], + "type": "NO2 tint" + }, + "quantum yield": { + "netcdf files": ["data/quantum_yields/NO2_1.nc"], + "type": "NO2 tint", + "lower extrapolation": { "type": "boundary" } + } + }, + { + "name": "jn2o5_a", + "__reaction": "N2O5 + hv -> NO2 + NO3", + "cross section": { + "type":"temperature based", + "netcdf file": "data/cross_sections/N2O5_JPL06.nc", + "parameterization": { + "type": "HARWOOD", + "aa": [ -18.27, -18.42, -18.59, -18.72, -18.84, + -18.90, -18.93, -18.87, -18.77, -18.71, + -18.31, -18.14, -18.01, -18.42, -18.59, + -18.13 ], + "bb": [ -91.0, -104.0, -112.0, -135.0, -170.0, + -226.0, -294.0, -388.0, -492.0, -583.0, + -770.0, -885.0, -992.0, -949.0, -966.0, + -1160.0 ], + "base temperature": 0.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "minimum wavelength": 260.0, + "maximum wavelength": 410.0, + "temperature ranges": [ + { + "maximum": 199.999999999999, + "fixed value": 200 + }, + { + "minimum": 200, + "maximum": 295 + }, + { + "minimum": 295.00000000001, + "fixed value": 295.0 + } + ] + }, + "parameterization wavelength grid": { + "name": "custom wavelengths", + "type": "from config file", + "units": "nm", + "values": [ + 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, + 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, + 375.0, 385.0, 395.0, 405.0, 415.0 + ] + } + }, + "quantum yield": { + "type": "Taylor series", + "constant value": 0.0, + "coefficients": [ -2.832441, 0.012809638 ], + "override bands": [ + { + "band": "range", + "minimum wavelength": 300.0, + "value": 1.0 + } + ] + } + }, + { + "name": "jn2o5_b", + "__reaction": "N2O5 + hv -> NO + O + NO3", + "cross section": { + "type":"temperature based", + "netcdf file": "data/cross_sections/N2O5_JPL06.nc", + "parameterization": { + "type": "HARWOOD", + "aa": [ -18.27, -18.42, -18.59, -18.72, -18.84, + -18.90, -18.93, -18.87, -18.77, -18.71, + -18.31, -18.14, -18.01, -18.42, -18.59, + -18.13 ], + "bb": [ -91.0, -104.0, -112.0, -135.0, -170.0, + -226.0, -294.0, -388.0, -492.0, -583.0, + -770.0, -885.0, -992.0, -949.0, -966.0, + -1160.0 ], + "base temperature": 0.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "minimum wavelength": 260.0, + "maximum wavelength": 410.0, + "temperature ranges": [ + { + "maximum": 199.999999999999, + "fixed value": 200 + }, + { + "minimum": 200, + "maximum": 295 + }, + { + "minimum": 295.00000000001, + "fixed value": 295.0 + } + ] + }, + "parameterization wavelength grid": { + "name": "custom wavelengths", + "type": "from config file", + "units": "nm", + "values": [ + 255.0, 265.0, 275.0, 285.0, 295.0, 305.0, + 315.0, 325.0, 335.0, 345.0, 355.0, 365.0, + 375.0, 385.0, 395.0, 405.0, 415.0 + ] + } + }, + "quantum yield": { + "type": "Taylor series", + "constant value": 0.0, + "coefficients": [ 3.832441, -0.012809638 ], + "override bands": [ + { + "band": "range", + "minimum wavelength": 300.0, + "value": 0.0 + } + ] + } + }, + { + "name": "jhno3", + "__reaction": "HNO3 + hv -> OH + NO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HNO3_JPL06.nc" } + ], + "type": "HNO3+hv->OH+NO2" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jno3_a", + "__reaction": "NO3 + hv -> NO2 + O(3P)", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/NO3_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/NO3-NO2+O(3P)_1.nc" + ], + "type": "tint", + "lower extrapolation": { + "type": "constant", + "value": 1.0 + } + } + }, + { + "name": "jno3_b", + "__reaction": "NO3 + hv -> NO + O2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/NO3_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/NO3-NO+O2_1.nc" + ], + "type": "tint" + } + }, + { + "name": "jch3ooh", + "__reaction": "CH3OOH + hv -> CH3O + OH", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3OOH_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jch2o_a", + "__reaction": "CH2O + hv -> H + HCO", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH2O_1.nc" } + ], + "type": "CH2O" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CH2O_1.nc" + ], + "type": "base", + "lower extrapolation": { + "type": "boundary" + } + } + }, + { + "name": "jch2o_b", + "__reaction": "CH2O + hv -> H2 + CO", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH2O_1.nc" } + ], + "type": "CH2O" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CH2O_1.nc" + ], + "type": "CH2O", + "lower extrapolation": { + "type": "boundary" + } + } + }, + { + "name": "jh2o2", + "__reaction": "H2O2 + hv -> OH + OH", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/H2O2_1.nc" } + ], + "type": "H2O2+hv->OH+OH" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jch3cho", + "__reaction": "CH3CHO + hv -> CH3 + HCO", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3CHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CH3CHO_1.nc" + ], + "type": "CH3CHO+hv->CH3+HCO" + } + }, + { + "name": "jpan", + "__reaction": "PAN + hv -> 0.6*CH3CO3 + 0.6*NO2 + 0.4*CH3O2 + 0.4*NO3 + 0.4*CO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/PAN_1.nc" } + ], + "type": "CH3ONO2+hv->CH3O+NO2" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jmvk", + "__reaction": "MVK + hv -> 0.7*C3H6 + 0.7*CO + 0.3*CH3O2 + 0.3*CH3CO3", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/MVK_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "MVK+hv->Products" + } + }, + { + "name": "jacet", + "__reaction": "CH3COCH3 + hv -> CH3CO + CH3", + "cross section": { + "type": "temperature based", + "parameterization": { + "type": "TAYLOR_SERIES", + "netcdf file": { + "file path": "data/cross_sections/ACETONE_JPL06.nc" + }, + "base temperature": 0.0, + "temperature ranges": [ + { + "maximum": 234.999999999999, + "fixed value": 235.0 + }, + { + "minimum": 235.0, + "maximum": 298.0 + }, + { + "minimum": 298.00000000001, + "fixed value": 298.0 + } + ] + } + }, + "quantum yield": { + "type": "CH3COCH3+hv->CH3CO+CH3", + "branch": "CO+CH3CO", + "low wavelength value": 1, + "minimum temperature": 218, + "maximum temperature": 295 + } + }, + { + "name": "jmgly", + "__reaction": "CH3COCHO + hv -> CH3CO3 + CO + HO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3COCHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "CH3COCHO+hv->CH3CO+HCO" + } + }, + { + "name": "jglyald", + "__reaction": "GLYALD + hv -> 2*HO2 + CO + CH2O", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HOCH2CHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.5 + } + }, +{ + "name": "jglyoxal", + "__reaction": "GLYOXAL + hv -> 2*CO + 2*HO2", + "__comments": "TODO the products of this reaction don't exactly match", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CHOCHO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "netcdf files": [ + "data/quantum_yields/CHOCHO-H2_CO_CO_1.nc" + ], + "type": "base", + "lower extrapolation": { + "type": "boundary" + } + } + }, + { + "name": "jbrcl", + "__reaction": "BrCl + hv -> Br + Cl", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/BrCl_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jbro", + "__reaction": "BrO + hv -> Br + O", + "cross section": { + "netcdf files": [ + { + "file path": "data/cross_sections/BRO_JPL06.nc" + } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jbrono2_a", + "__reaction": "BrONO2 + hv -> Br + NO3", + "cross section": { + "type": "temperature based", + "parameterization": { + "type": "TAYLOR_SERIES", + "netcdf file": { + "file path": "data/cross_sections/BRONO2_JPL06.nc" + }, + "base temperature": 296.0, + "temperature ranges": [ + { + "maximum": 199.999999999999, + "fixed value": 200.0 + }, + { + "minimum": 200.0, + "maximum": 296.0 + }, + { + "minimum": 296.00000000001, + "fixed value": 296.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.85 + } + }, + { + "name": "jbrono2_b", + "__reaction": "BrONO2 + hv -> BrO + NO2", + "cross section": { + "type": "temperature based", + "parameterization": { + "type": "TAYLOR_SERIES", + "netcdf file": { + "file path": "data/cross_sections/BRONO2_JPL06.nc" + }, + "base temperature": 296.0, + "temperature ranges": [ + { + "maximum": 199.999999999999, + "fixed value": 200.0 + }, + { + "minimum": 200.0, + "maximum": 296.0 + }, + { + "minimum": 296.00000000001, + "fixed value": 296.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.15 + } + }, + { + "name": "jccl4", + "__reaction": "CCl4 + hv -> Products", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CCl4_1.nc" } + ], + "type": "CCl4+hv->Products" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcf2clbr", + "__reaction": "CF2BrCl + hv -> Products", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CF2BrCl_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcf3br", + "__reaction": "CF3Br + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/H1301_JPL06.nc", + "parameterization": { + "AA": [ 62.563, -2.0068, 1.6592e-2, -5.6465e-5, 6.7459e-8 ], + "BB": [ -9.1755e-1, 1.8575e-2, -1.3857e-4, 4.5066e-7, -5.3803e-10 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 178.0, + "maximum wavelength": 280.0, + "base temperature": 273.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "temperature ranges": [ + { + "maximum": 209.999999999999, + "fixed value": 210.0 + }, + { + "minimum": 210.0, + "maximum": 300.0 + }, + { + "minimum": 300.00000000001, + "fixed value": 300.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcfcl3", + "__reaction": "CCl3F + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CFCL3_JPL06.nc", + "parameterization": { + "AA": [ -84.611, 7.9551e-1, -2.0550e-3, -4.4812e-6, 1.5838e-8 ], + "BB": [ -5.7912, 1.1689e-1, -8.8069e-4, 2.9335e-6, -3.6421e-9 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 174.1, + "maximum wavelength": 230.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": "jcfc113", + "__reaction": "CFC-113 + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CFC113_JPL06.nc", + "parameterization": { + "AA": [ -1087.9, 20.004, -1.3920e-1, 4.2828e-4, -4.9384e-7 ], + "BB": [ 12.493, -2.3937e-1, 1.7142e-3, -5.4393e-6, 6.4548e-9 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 182.0, + "maximum wavelength": 230.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": "jcfc114", + "__reaction": "CFC-114 + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CFC114_JPL10.nc", + "parameterization": { + "AA": [ -160.50, 2.4807, -1.5202e-2, 3.8412e-5, -3.4373e-8 ], + "BB": [ -1.5296, 3.5248e-2, -2.9951e-4, 1.1129e-6, -1.5259e-9 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 172.0, + "maximum wavelength": 220.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": "jcfc115", + "__reaction": "CFC-115 + hv -> Products", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CFC115_JPL10.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcf2cl2", + "__reaction": "CCl2F2 + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CF2CL2_JPL06.nc", + "parameterization": { + "AA": [ -43.8954569, -2.403597e-1, -4.2619e-4, 9.8743e-6, 0.0 ], + "BB": [ 4.8438e-3, 4.96145e-4, -5.6953e-6, 0.0, 0.0 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 200.0, + "maximum wavelength": 231.0, + "base temperature": 296.0, + "base wavelength": 200.0, + "logarithm": "natural", + "temperature ranges": [ + { + "maximum": 219.999999999999, + "fixed value": 220.0 + }, + { + "minimum": 220, + "maximum": 296 + }, + { + "minimum": 296.00000000001, + "fixed value": 296.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jch2br2", + "__reaction": "CH2BR2 + hv -> 2*BR", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CH2BR2_1.nc", + "parameterization": { + "AA": [ -70.211776, 1.940326e-1, 2.726152e-3, -1.695472e-5, 2.500066e-8 ], + "BB": [ 2.899280, -4.327724e-2, 2.391599e-4, -5.807506e-7, 5.244883e-10 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 210.0, + "maximum wavelength": 290.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": "jch3br", + "__reaction": "CH3Br + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CH3BR_JPL06.nc", + "parameterization": { + "AA": [ 46.520, -1.4580, 1.1469e-2, -3.7627e-5, 4.3264e-8 ], + "BB": [ 9.3408e-1, -1.6887e-2, 1.1487e-4, -3.4881e-7, 3.9945e-10 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 200.0, + "maximum wavelength": 280.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": "jch3ccl3", + "__reaction": "CH3CCl3+hv->Products", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CH3CCl3_1.nc" } + ], + "type": "tint" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jch3cl", + "__reaction": "CH3Cl + hv -> Products", + "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 + } + }, + { + "name": "jchbr3", + "__reaction": "CHBr3 + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/CHBR3_JPL10.nc", + "parameterization": { + "AA": [ -32.6067, 0.10308, 6.39e-5, -7.7392e-7, -2.2513e-9, 6.1376e-12 ], + "BB": [ 0.1582, -0.0014758, 3.8058e-6, 9.187e-10, -1.0772e-11, 0.0 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 ], + "minimum wavelength": 260.0, + "maximum wavelength": 362.0, + "base temperature": 296.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "invert temperature offset": true, + "temperature ranges": [ + { + "maximum": 259.999999999999, + "fixed value": 260.0 + }, + { + "minimum": 260.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcl2", + "__reaction": "Cl2 + hv -> Cl + Cl", + "cross section": { + "type": "Cl2+hv->Cl+Cl" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcl2o2", + "__reaction": "ClOOCl + hv -> Cl + ClOO", + "__comments": "TODO - this doesn't exactly match the products in TS1", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jclo", + "__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", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/ClONO2_1.nc" } + ], + "type": "ClONO2" + }, + "quantum yield": { + "type": "ClONO2+hv->Cl+NO3" + } + }, + { + "name": "jclono2_b", + "__reaction": "ClONO2 + hv -> ClO + NO2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/ClONO2_1.nc" } + ], + "type": "ClONO2" + }, + "quantum yield": { + "type": "ClONO2+hv->ClO+NO2" + } + }, + { + "name": "jcof2", + "__reaction": "CF2O + hv -> Products", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CF2O_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jcofcl", + "__reaction": "CClFO + hv -> Products", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CClFO_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jh2402", + "__reaction": "H2402 + hv -> 2*BR + 2*COF2", + "__comments": "TUV data set name CF2BrCF2Br", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/H2402_JPL06.nc", + "parameterization": { + "AA": [ 34.026, -1.152616, 8.959798e-3, -2.9089e-5, 3.307212e-8 ], + "BB": [ 4.010664e-1, -8.358968e-3, 6.415741e-5, -2.157554e-7, 2.691871e-10 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 190.0, + "maximum wavelength": 290.0, + "base temperature": 273.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "temperature ranges": [ + { + "maximum": 209.999999999999, + "fixed value": 210.0 + }, + { + "minimum": 210.0, + "maximum": 300.0 + }, + { + "minimum": 300.00000000001, + "fixed value": 300.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhcfc141b", + "__reaction": "HCFC-141b + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/HCFC141b_JPL10.nc", + "parameterization": { + "AA": [ -682.913042, 12.122290, -8.187699e-2, 2.437244e-4, -2.719103e-7 ], + "BB": [ 4.074747, -8.053899e-2, 5.946552e-4, -1.945048e-6, 2.380143e-9 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 172.0, + "maximum wavelength": 240.0, + "base temperature": 273.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "temperature ranges": [ + { + "maximum": 209.999999999999, + "fixed value": 210.0 + }, + { + "minimum": 210.0, + "maximum": 300.0 + }, + { + "minimum": 300.00000000001, + "fixed value": 300.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhcfc142b", + "__reaction": "HCFC-142b + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/HCFC142b_JPL10.nc", + "parameterization": { + "AA": [ -328.092008, 6.342799, -4.810362e-2, 1.611991e-4, -2.042613e-7 ], + "BB": [ 4.289533e-1, -9.042817e-3, 7.018009e-5, -2.389064e-7, 3.039799e-10 ], + "lp": [ 0.0, 1.0, 2.0, 3.0, 4.0 ], + "minimum wavelength": 172.0, + "maximum wavelength": 230.0, + "base temperature": 273.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "temperature ranges": [ + { + "maximum": 209.999999999999, + "fixed value": 210.0 + }, + { + "minimum": 210.0, + "maximum": 300.0 + }, + { + "minimum": 300.00000000001, + "fixed value": 300.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhcfc22", + "__reaction": "HCFC-22 + hv -> Products", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/HCFC22_JPL06.nc", + "parameterization wavelength grid": { + "name": "custom wavelengths", + "type": "from config file", + "units": "nm", + "values": [ + 169.0, 171.0, 173.0, 175.0, 177.0, 179.0, 181.0, 183.0, 185.0, + 187.0, 189.0, 191.0, 193.0, 195.0, 197.0, 199.0, 201.0, 203.0, + 205.0, 207.0, 209.0, 211.0, 213.0, 215.0, 217.0, 219.0, 221.0 + ] + }, + "parameterization": { + "AA": [ -106.029, 1.5038, -8.2476e-3, 1.4206e-5 ], + "BB": [ -1.3399e-1, 2.7405e-3, -1.8028e-5, 3.8504e-8 ], + "lp": [ 0.0, 1.0, 2.0, 3.0 ], + "minimum wavelength": 174.0, + "maximum wavelength": 204.0, + "base temperature": 273.0, + "base wavelength": 0.0, + "logarithm": "base 10", + "temperature ranges": [ + { + "maximum": 209.999999999999, + "fixed value": 210.0 + }, + { + "minimum": 210.0, + "maximum": 300.0 + }, + { + "minimum": 300.00000000001, + "fixed value": 300.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhcl", + "__reaction": "HCl + hv -> H + Cl", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HCl_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhobr", + "__reaction": "HOBr + hv -> OH + Br", + "cross section": { + "type": "HOBr+hv->OH+Br" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhocl", + "__reaction": "HOCl + hv -> HO + Cl", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/HOCl_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "joclo", + "__reaction": "OClO + hv -> Products", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/OClO_1.nc" }, + { "file path": "data/cross_sections/OClO_2.nc" }, + { "file path": "data/cross_sections/OClO_3.nc" } + ], + "type": "OClO+hv->Products" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jho2no2_a", + "__reaction": "HNO4 + hv -> OH + NO3", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/HO2NO2_JPL06.nc", + "parameterization": { + "type": "BURKHOLDER", + "netcdf file": { + "file path": "data/cross_sections/HO2NO2_temp_JPL06.nc" + }, + "A": -988.0, + "B": 0.69, + "temperature ranges": [ + { + "maximum": 279.999999999999, + "fixed value": 280.0 + }, + { + "minimum": 280.0, + "maximum": 350.0 + }, + { + "minimum": 350.00000000001, + "fixed value": 350.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.30, + "override bands": [ + { + "band": "range", + "minimum wavelength": 200.0, + "value": 0.20 + } + ] + } + }, + { + "name": "jho2no2_b", + "__reaction": "HNO4 + hv -> HO2 + NO2", + "cross section": { + "type": "temperature based", + "netcdf file": "data/cross_sections/HO2NO2_JPL06.nc", + "parameterization": { + "type": "BURKHOLDER", + "netcdf file": { + "file path": "data/cross_sections/HO2NO2_temp_JPL06.nc" + }, + "A": -988.0, + "B": 0.69, + "temperature ranges": [ + { + "maximum": 279.999999999999, + "fixed value": 280.0 + }, + { + "minimum": 280.0, + "maximum": 350.0 + }, + { + "minimum": 350.00000000001, + "fixed value": 350.0 + } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.70, + "override bands": [ + { + "band": "range", + "minimum wavelength": 200.0, + "value": 0.80 + } + ] + } + }, + { + "name": "jmacr_a", + "__reaction": "CH2=C(CH3)CHO->1.34HO2+0.66MCO3+1.34CH2O+CH3CO3", + "__comments": "Methacrolein photolysis channel 1", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/Methacrolein_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.005 + } + }, + { + "name": "jmacr_b", + "__reaction": "CH2=C(CH3)CHO->0.66OH+1.34CO", + "__comments": "Methacrolein photolysis channel 2", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/Methacrolein_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.005 + } + }, + { + "name": "jhyac", + "__reaction": "CH2(OH)COCH3->CH3CO3+HO2+CH2O", + "__comments": "hydroxy acetone TODO: the products of this reaction differ from standalone TUV-x", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/Hydroxyacetone_1.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 0.65 + } + }, + { + "name": "jh2o_a", + "__reaction": "H2O + hv -> OH + H", + "cross section": { + "type": "base", + "merge data": true, + "netcdf files": [ + { + "file path": "data/cross_sections/H2O_1.nc", + "zero above": 183.0 + }, + { + "file path": "data/cross_sections/H2O_2.nc", + "zero below": 183.00000000001, + "zero above": 190.0 + }, + { + "file path": "data/cross_sections/H2O_3.nc", + "zero below": 190.00000000001 + } + ] + }, + "quantum yield" : { + "type": "base", + "netcdf files": [ "data/quantum_yields/H2O_H_OH.nc" ] + } + }, + { + "name": "jh2o_b", + "__reaction": "H2O + hv -> H2 + O1D", + "cross section": { + "type": "base", + "merge data": true, + "netcdf files": [ + { + "file path": "data/cross_sections/H2O_1.nc", + "zero above": 183.0 + }, + { + "file path": "data/cross_sections/H2O_2.nc", + "zero below": 183.00000000001, + "zero above": 190.0 + }, + { + "file path": "data/cross_sections/H2O_3.nc", + "zero below": 190.00000000001 + } + ] + }, + "quantum yield" : { + "type": "base", + "netcdf files": [ "data/quantum_yields/H2O_H2_O1D.nc" ] + } + }, + { + "name": "jh2o_c", + "__reaction": "H2O + hv -> 2*H + O", + "cross section": { + "type": "base", + "merge data": true, + "netcdf files": [ + { + "file path": "data/cross_sections/H2O_1.nc", + "zero above": 183.0 + }, + { + "file path": "data/cross_sections/H2O_2.nc", + "zero below": 183.00000000001, + "zero above": 190.0 + }, + { + "file path": "data/cross_sections/H2O_3.nc", + "zero below": 190.00000000001 + } + ] + }, + "quantum yield" : { + "type": "base", + "netcdf files": [ "data/quantum_yields/H2O_2H_O3P.nc" ] + } + }, + { + "name": "jch4_a", + "__reaction": "CH4 + hv -> H + CH3O2", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CH4_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 0.45 + } + }, + { + "name": "jch4_b", + "__reaction": "CH4 + hv -> 1.44*H2 + 0.18*CH2O + 0.18*O + 0.33*OH + 0.33*H + 0.44*CO2 + 0.38*CO + 0.05*H2O", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CH4_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 0.55 + } + }, + { + "name": "jco2", + "__reaction": "CO2 + hv -> CO + O", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CO2_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhbr", + "__reaction": "HBR + hv -> BR + H", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/HBr_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jhf", + "__reaction": "HF + hv -> H + F", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/HF_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jsf6", + "__reaction": "SF6 + hv -> sink", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/SF6_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jh2so4", + "__reaction": "H2SO4 + hv -> SO3 + H2O", + "cross section": { + "type": "base", + "data": { + "default value": 0.0, + "point values": [ + { "wavelength": 121.65, "value": 6.3e-17 }, + { "wavelength": 525.0, "value": 1.43e-26 }, + { "wavelength": 625.0, "value": 1.8564e-25 }, + { "wavelength": 725.0, "value": 3.086999e-24 } + ] + } + }, + "quantum yield": { + "type": "H2SO4 Mills", + "netcdf files": [ + "data/quantum_yields/H2SO4_mills.nc" + ], + "parameterized wavelengths": [ + 525, + 625, + 725 + ], + "collision interval s": [ + 1.1e-9, + 8.9e-9, + 1.7e-7 + ], + "molecular diameter m": 4.18e-10, + "molecular weight kg mol-1": 98.078479e-3 + } + }, + { + "name": "jocs", + "__reaction": "OCS + hv -> S + CO", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/OCS_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jso", + "__reaction": "SO + hv -> S + O", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/SO_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jso2", + "__reaction": "SO2 + hv -> SO + O", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/SO2_Mills.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jso3", + "__reaction": "SO3 + hv -> SO2 + O", + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/SO3_1.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, + { + "name": "jno_i", + "__reaction": "NO + hv -> NOp + e", + "cross section": { + "type": "base", + "data": { + "default value": 0.0, + "point values": [ + { "wavelength": 121.65, "value": 2.0e-18 } + ] + } + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + } + ] + }, + "__CAM options": { + "disable clouds": true, + "disable aerosols": true, + "aliasing": { + "default matching": "backup", + "pairs": [ + { + "to": "jalknit", + "__reaction": "ALKNIT + hv -> NO2 + 0.4*CH3CHO + 0.1*CH2O + 0.25*CH3COCH3 + HO2 + 0.8*MEK", + "from": "jch3ooh" + }, + { + "to": "jpooh", + "__reaction": "POOH (C3H6OHOOH) + hv -> CH3CHO + CH2O + HO2 + OH", + "from": "jch3ooh" + }, + { + "to": "jch3co3h", + "__reaction": "CH3COOOH + hv -> CH3O2 + OH + CO2", + "from": "jh2o2", + "scale by": 0.28 + }, + { + "to": "jmpan", + "__reaction": "MPAN + hv -> MCO3 + NO2", + "from": "jpan" + }, + { + "to": "jc2h5ooh", + "__reaction": "C2H5OOH + hv -> CH3CHO + HO2 + OH", + "from": "jch3ooh" + }, + { + "to": "jc3h7ooh", + "__reaction": "C3H7OOH + hv -> 0.82*CH3COCH3 + OH + HO2", + "from": "jch3ooh" + }, + { + "to": "jc6h5ooh", + "__reaction": "C6H5OOH + hv -> PHENO + OH", + "from": "jch3ooh" + }, + { + "to": "jeooh", + "__reaction": "EOOH + hv -> EO + OH", + "from": "jch3ooh" + }, + { + "to": "jrooh", + "__reaction": "ROOH + hv -> CH3CO3 + CH2O + OH", + "from": "jch3ooh" + }, + { + "to": "jxooh", + "__reaction": "XOOH + hv -> OH", + "from": "jch3ooh" + }, + { + "to": "jonitr", + "__reaction": "ONITR + hv -> NO2", + "from": "jch3cho" + }, + { + "to": "jisopooh", + "__reaction": "ISOPOOH + hv -> 0.402*MVK + 0.288*MACR + 0.69*CH2O + HO2", + "from": "jch3ooh" + }, + { + "to": "jmek", + "__reaction": "MEK + hv -> CH3CO3 + C2H5O2", + "from": "jacet" + }, + { + "to": "jalkooh", + "__reaction": "ALKOOH + hv -> .4*CH3CHO + .1*CH2O + .25*CH3COCH3 + .9*HO2 + .8*MEK + OH", + "from": "jch3ooh" + }, + { + "to": "jbenzooh", + "__reaction": "BENZOOH + hv -> OH + GLYOXAL + 0.5*BIGALD1 + HO2", + "from": "jch3ooh" + }, + { + "to": "jbepomuc", + "__reaction": "BEPOMUC + hv -> BIGALD1 + 1.5*HO2 + 1.5*CO", + "from": "jno2", + "scale by": 0.1 + }, + { + "to": "jbigald", + "__reaction": "BIGALD + hv -> 0.45*CO + 0.13*GLYOXAL + 0.56*HO2 + 0.13*CH3CO3 + 0.18*CH3COCHO", + "from": "jno2", + "scale by": 0.2 + }, + { + "to": "jbigald1", + "__reaction": "BIGALD1 + hv -> 0.6*MALO2 + HO2", + "from": "jno2", + "scale by": 0.14 + }, + { + "to": "jbigald2", + "__reaction": "BIGALD2 + hv -> 0.6*HO2 + 0.6*DICARBO2", + "from": "jno2", + "scale by": 0.2 + }, + { + "to": "jbigald3", + "__reaction": "BIGALD3 + hv -> 0.6*HO2 + 0.6*CO + 0.6*MDIALO2", + "from": "jno2", + "scale by": 0.2 + }, + { + "to": "jbigald4", + "__reaction": "BIGALD4 + hv -> HO2 + CO + CH3COCHO + CH3CO3", + "from": "jno2", + "scale by": 0.006 + }, + { + "to": "jbzooh", + "__reaction": "BZOOH + hv -> BZALD + OH + HO2", + "from": "jch3ooh" + }, + { + "to": "jmekooh", + "__reaction": "MEKOOH + hv -> OH + CH3CO3 + CH3CHO", + "from": "jch3ooh" + }, + { + "to": "jtolooh", + "__reaction": "TOLOOH + hv -> OH + .45*GLYOXAL + .45*CH3COCHO + .9*BIGALD", + "from": "jch3ooh" + }, + { + "to": "jterpooh", + "__reaction": "TERPOOH + hv -> OH + .1*CH3COCH3 + HO2 + MVK + MACR", + "from": "jch3ooh" + }, + { + "to": "jhonitr", + "__reaction": "HONITR + hv -> NO2 + 0.67*HO2 + 0.33*CH3CHO + 0.33*CH2O + 0.33*CO + 0.33*GLYALD + 0.33*CH3CO3 + 0.17*HYAC + 0.17*CH3COCH3", + "from": "jch2o_a" + }, + { + "to": "jhpald", + "__reaction": "HPALD + hv -> BIGALD3 + OH + HO2", + "from": "jno2", + "scale by": 0.006 + }, + { + "to": "jisopnooh", + "__reaction": "ISOPNOOH + hv -> NO2 + HO2 + ISOPOOH", + "from": "jch3ooh" + }, + { + "to": "jnc4cho", + "__reaction": "NC4CHO + hv -> BIGALD3 + NO2 + HO2", + "from": "jch2o_a" + }, + { + "to": "jnoa", + "__reaction": "NOA + hv -> NO2 + CH2O + CH3CO3", + "from": "jch2o_a" + }, + { + "to": "jnterpooh", + "__reaction": "NTERPOOH + hv -> TERPROD1 + NO2 + OH", + "from": "jch3ooh" + }, + { + "to": "jphenooh", + "__reaction": "PHENOOH + hv -> OH + HO2 + 0.7*GLYOXAL", + "from": "jch3ooh" + }, + { + "to": "jtepomuc", + "__reaction": "TEPOMUC + hv -> 0.5*CH3CO3 + HO2 + 1.5*CO", + "from": "jno2", + "scale by": 0.1 + }, + { + "to": "jterp2ooh", + "__reaction": "TERP2OOH + hv -> OH + 0.375*CH2O + 0.3*CH3COCH3 + 0.25*CO + CO2 + TERPROD2 + HO2 + 0.25*GLYALD", + "from": "jch3ooh" + }, + { + "to": "jterpnit", + "__reaction": "TERPNIT + hv -> TERPROD1 + NO2 + HO2", + "from": "jch3ooh" + }, + { + "to": "jterprd1", + "__reaction": "TERPROD1 + hv -> HO2 + CO + TERPROD2", + "from": "jch3cho" + }, + { + "to": "jterprd2", + "__reaction": "TERPROD2 + hv -> 0.15*RO2 + 0.68*CH2O + 0.8*CO2 + 0.5*CH3COCH3 + 0.65*CH3CO3 + 1.2*HO2 + 1.7*CO", + "from": "jch3cho" + }, + { + "to": "jxylenooh", + "__reaction": "XYLENOOH + hv -> OH + HO2 + 0.34*GLYOXAL + 0.54*CH3COCHO + 0.06*BIGALD1 + 0.2*BIGALD2 + 0.15*BIGALD3 + 0.21*BIGALD4", + "from": "jch3ooh" + }, + { + "to": "jxylolooh", + "__reaction": "XYLOLOOH + hv -> OH + 0.17*GLYOXAL + 0.51*CH3COCHO + HO2", + "from": "jch3ooh" + }, + { + "to": "jsoa1_a1", + "__reaction": "soa1_a1 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa1_a2", + "__reaction": "soa1_a2 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa2_a1", + "__reaction": "soa2_a1 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa2_a2", + "__reaction": "soa2_a2 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa3_a1", + "__reaction": "soa3_a1 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa3_a2", + "__reaction": "soa3_a2 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa4_a1", + "__reaction": "soa4_a1 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa4_a2", + "__reaction": "soa4_a2 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa5_a1", + "__reaction": "soa5_a1 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + }, + { + "to": "jsoa5_a2", + "__reaction": "soa5_a2 + hv -> Products", + "from": "jno2", + "scale by": 0.0004 + } + ] + } + } +} diff --git a/src/chemistry/mozart/chemistry.F90 b/src/chemistry/mozart/chemistry.F90 index ff42e870d9..9557907f1f 100644 --- a/src/chemistry/mozart/chemistry.F90 +++ b/src/chemistry/mozart/chemistry.F90 @@ -162,6 +162,7 @@ subroutine chem_register use short_lived_species, only : slvd_index, short_lived_map=>map, register_short_lived_species use cfc11star, only : register_cfc11star use mo_photo, only : photo_register + use mo_tuvx, only : tuvx_register, tuvx_active use mo_aurora, only : aurora_register use aero_model, only : aero_model_register use physics_buffer, only : pbuf_add_field, dtype_r8 @@ -311,9 +312,12 @@ subroutine chem_register call register_cfc11star() if ( waccmx_is('ionosphere') ) then - call photo_register() + if( .not. tuvx_active ) then + call photo_register( ) + end if call aurora_register() endif + call tuvx_register( ) ! add fields to pbuf needed by aerosol models call aero_model_register() @@ -340,6 +344,7 @@ subroutine chem_readnl(nlfile) use mo_sulf, only: sulf_readnl use species_sums_diags,only: species_sums_readnl use ocean_emis, only: ocean_emis_readnl + use mo_tuvx, only: tuvx_readnl ! args @@ -551,6 +556,7 @@ subroutine chem_readnl(nlfile) call sulf_readnl(nlfile) call species_sums_readnl(nlfile) call ocean_emis_readnl(nlfile) + call tuvx_readnl(nlfile) end subroutine chem_readnl @@ -1035,6 +1041,7 @@ subroutine chem_timestep_init(phys_state,pbuf2d) use mo_aurora, only : aurora_timestep_init use mo_photo, only : photo_timestep_init + use mo_tuvx, only : tuvx_active, tuvx_timestep_init use cfc11star, only : update_cfc11star use physics_buffer, only : physics_buffer_desc @@ -1106,6 +1113,11 @@ subroutine chem_timestep_init(phys_state,pbuf2d) !----------------------------------------------------------------------------- call photo_timestep_init( calday ) + !----------------------------------------------------------------------------- + ! ... setup the TUV-x profiles for this timestep + !----------------------------------------------------------------------------- + if( tuvx_active ) call tuvx_timestep_init( ) + call update_cfc11star( pbuf2d, phys_state ) ! Galatic Cosmic Rays ... @@ -1258,7 +1270,7 @@ subroutine chem_timestep_tend( state, ptend, cam_in, cam_out, dt, pbuf, fh2o) fsds, cam_in%ts, cam_in%asdir, cam_in%ocnfrac, cam_in%icefrac, & cam_out%precc, cam_out%precl, cam_in%snowhland, ghg_chem, state%latmapback, & drydepflx, wetdepflx, cam_in%cflx, cam_in%fireflx, cam_in%fireztop, & - nhx_nitrogen_flx, noy_nitrogen_flx, use_hemco, ptend%q, pbuf ) + nhx_nitrogen_flx, noy_nitrogen_flx, use_hemco, ptend%q, pbuf, state ) if (associated(cam_out%nhx_nitrogen_flx)) then cam_out%nhx_nitrogen_flx(:ncol) = nhx_nitrogen_flx(:ncol) endif @@ -1333,10 +1345,12 @@ subroutine chem_final() use mee_ionization, only: mee_ion_final use rate_diags, only: rate_diags_final use species_sums_diags, only: species_sums_final + use mo_tuvx, only: tuvx_finalize call mee_ion_final() call rate_diags_final() call species_sums_final() + call tuvx_finalize() end subroutine chem_final diff --git a/src/chemistry/mozart/mo_chemini.F90 b/src/chemistry/mozart/mo_chemini.F90 index 0f4005be96..a13de223b9 100644 --- a/src/chemistry/mozart/mo_chemini.F90 +++ b/src/chemistry/mozart/mo_chemini.F90 @@ -48,6 +48,7 @@ subroutine chemini & use mo_srf_emissions, only : srf_emissions_inti use mo_sulf, only : sulf_inti use mo_photo, only : photo_inti + use mo_tuvx, only : tuvx_init, tuvx_active use mo_drydep, only : drydep_inti use mo_imp_sol, only : imp_slv_inti use mo_exp_sol, only : exp_sol_inti @@ -196,11 +197,17 @@ subroutine chemini & call euvac_init (euvac_file) call photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & - photon_file, electron_file, & - exo_coldens_file, photo_max_zen ) - + photon_file, electron_file, exo_coldens_file, photo_max_zen ) if (masterproc) write(iulog,*) 'chemini: after photo_inti on node ',iam + !----------------------------------------------------------------------- + ! ... initialize the TUV-x photolysis rate constant calculator + !----------------------------------------------------------------------- + if( tuvx_active ) then + call tuvx_init( photon_file, electron_file, photo_max_zen, pbuf2d ) + if (masterproc) write(iulog,*) 'chemini: after tuvx_init on node ',iam + end if + !----------------------------------------------------------------------- ! ... initialize ion production !----------------------------------------------------------------------- diff --git a/src/chemistry/mozart/mo_gas_phase_chemdr.F90 b/src/chemistry/mozart/mo_gas_phase_chemdr.F90 index 68657d0739..c894e99c74 100644 --- a/src/chemistry/mozart/mo_gas_phase_chemdr.F90 +++ b/src/chemistry/mozart/mo_gas_phase_chemdr.F90 @@ -259,8 +259,9 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, & delt, ps, & fsds, ts, asdir, ocnfrac, icefrac, & precc, precl, snowhland, ghg_chem, latmapback, & - drydepflx, wetdepflx, cflx, fire_sflx, fire_ztop, nhx_nitrogen_flx, noy_nitrogen_flx, & - use_hemco, qtend, pbuf) + drydepflx, wetdepflx, cflx, fire_sflx, fire_ztop, & + nhx_nitrogen_flx, noy_nitrogen_flx, & + use_hemco, qtend, pbuf, state) !----------------------------------------------------------------------- ! ... Chem_solver advances the volumetric mixing ratio @@ -272,6 +273,7 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, & use chem_mods, only : nabscol, nfs, indexm, clscnt4 use physconst, only : rga, gravit use mo_photo, only : set_ub_col, setcol, table_photo + use mo_tuvx, only : tuvx_get_photo_rates, tuvx_active use mo_exp_sol, only : exp_sol use mo_imp_sol, only : imp_sol use mo_setrxt, only : setrxt @@ -304,6 +306,7 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, & use mo_chm_diags, only : chm_diags, het_diags use perf_mod, only : t_startf, t_stopf use gas_wetdep_opts, only : gas_wetdep_method + use physics_types, only : physics_state use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx use infnan, only : nan, assignment(=) use rate_diags, only : rate_diags_calc, rate_diags_o3s_loss @@ -363,6 +366,7 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, & logical, intent(in) :: use_hemco ! use Harmonized Emissions Component (HEMCO) type(physics_buffer_desc), pointer :: pbuf(:) + type(physics_state), target, intent(in) :: state !----------------------------------------------------------------------- ! ... Local variables @@ -811,12 +815,22 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, & call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr , & delta, esfact ) - !----------------------------------------------------------------- - ! ... lookup the photolysis rates from table - !----------------------------------------------------------------- - call table_photo( reaction_rates, pmid, pdel, tfld, zmid, zint, & - col_dens, zen_angle, asdir, cwat, cldfr, & - esfact, vmr, invariants, ncol, lchnk, pbuf ) + if (tuvx_active) then + !----------------------------------------------------------------- + ! ... get calculated photolysis rates from TUV-x + !----------------------------------------------------------------- + call tuvx_get_photo_rates( state, pbuf, ncol, lchnk, zmid, zint, & + tfld, ts, invariants, vmr, col_delta, & + asdir, zen_angle, esfact, pdel, cldfr,& + cwat, reaction_rates(:,:,1:phtcnt) ) + else + !----------------------------------------------------------------- + ! ... lookup the photolysis rates from table + !----------------------------------------------------------------- + call table_photo( reaction_rates, pmid, pdel, tfld, zmid, zint, & + col_dens, zen_angle, asdir, cwat, cldfr, & + esfact, vmr, invariants, ncol, lchnk, pbuf ) + endif do i = 1,phtcnt call outfld( tag_names(i), reaction_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk ) diff --git a/src/chemistry/mozart/mo_jeuv.F90 b/src/chemistry/mozart/mo_jeuv.F90 index c13255250e..b893ad9f20 100644 --- a/src/chemistry/mozart/mo_jeuv.F90 +++ b/src/chemistry/mozart/mo_jeuv.F90 @@ -80,6 +80,7 @@ subroutine jeuv_init (photon_file, electron_file, indexer) character(len=2) :: mstring character(len=7) :: jstring logical :: do_jeuv + logical, save :: is_initialized = .false. do_jeuv=.false. do_heating=.false. @@ -97,6 +98,11 @@ subroutine jeuv_init (photon_file, electron_file, indexer) endif enddo + ! If the module has already been initialized, return after + ! computing index map + if( is_initialized ) return + is_initialized = .true. + if (.not.do_jeuv) return if (solar_euv_data_active) then diff --git a/src/chemistry/mozart/mo_jshort.F90 b/src/chemistry/mozart/mo_jshort.F90 index aa47dffb31..65efa66de4 100644 --- a/src/chemistry/mozart/mo_jshort.F90 +++ b/src/chemistry/mozart/mo_jshort.F90 @@ -27,6 +27,7 @@ module mo_jshort public :: jshort public :: sphers public :: slant_col + public :: calc_jno public :: nj save diff --git a/src/chemistry/mozart/mo_tuvx.F90 b/src/chemistry/mozart/mo_tuvx.F90 new file mode 100644 index 0000000000..8791053431 --- /dev/null +++ b/src/chemistry/mozart/mo_tuvx.F90 @@ -0,0 +1,1947 @@ +!---------------------------------------------------------------------- +! Wrapper for TUV-x photolysis rate constant calculator +!---------------------------------------------------------------------- +module mo_tuvx + + use musica_map, only : map_t + use musica_string, only : string_t + use ppgrid, only : pver, & ! number of vertical layers + pverp ! number of vertical interfaces (pver + 1) + use shr_kind_mod, only : r8 => shr_kind_r8, cl=>shr_kind_cl + use tuvx_core, only : core_t + use tuvx_grid_from_host, only : grid_updater_t + use tuvx_profile_from_host, only : profile_updater_t + use tuvx_radiator_from_host, only : radiator_updater_t + + implicit none + + private + + public :: tuvx_readnl + public :: tuvx_register + public :: tuvx_init + public :: tuvx_timestep_init + public :: tuvx_get_photo_rates + public :: tuvx_finalize + public :: tuvx_active + + ! Inidices for grid updaters + integer, parameter :: NUM_GRIDS = 2 ! number of grids that CAM will update at runtime + integer, parameter :: GRID_INDEX_HEIGHT = 1 ! Height grid index + integer, parameter :: GRID_INDEX_WAVELENGTH = 2 ! Wavelength grid index + + ! Indices for profile updaters + integer, parameter :: NUM_PROFILES = 8 ! number of profiles that CAM will update at runtime + integer, parameter :: PROFILE_INDEX_TEMPERATURE = 1 ! Temperature profile index + integer, parameter :: PROFILE_INDEX_ALBEDO = 2 ! Surface albedo profile index + integer, parameter :: PROFILE_INDEX_ET_FLUX = 3 ! Extraterrestrial flux profile index + integer, parameter :: PROFILE_INDEX_AIR = 4 ! Air density profile index + integer, parameter :: PROFILE_INDEX_O3 = 5 ! Ozone profile index + integer, parameter :: PROFILE_INDEX_O2 = 6 ! Molecular oxygen profile index + integer, parameter :: PROFILE_INDEX_SO2 = 7 ! Sulfur dioxide profile index + integer, parameter :: PROFILE_INDEX_NO2 = 8 ! Nitrogen dioxide profile index + + ! Indices for radiator updaters + integer, parameter :: NUM_RADIATORS = 2 ! number of radiators that CAM will update at runtime + integer, parameter :: RADIATOR_INDEX_AEROSOL = 1 ! Aerosol radiator index + integer, parameter :: RADIATOR_INDEX_CLOUDS = 2 ! Cloud radiator index + + ! Definition of the MS93 wavelength grid TODO add description of this + integer, parameter :: NUM_BINS_MS93 = 4 + real(kind=r8), parameter :: WAVELENGTH_EDGES_MS93(NUM_BINS_MS93+1) = & + (/ 181.6_r8, 183.1_r8, 184.6_r8, 190.2_r8, 192.5_r8 /) + + ! Heating rate indices + integer :: number_of_heating_rates = 0 ! number of heating rates in TUV-x + integer :: index_cpe_jo2_a = -1 ! index for jo2_a in heating rate array + integer :: index_cpe_jo2_b = -1 ! index for jo2_b in heating rate array + integer :: index_cpe_jo3_a = -1 ! index for jo3_a in heating rate array + integer :: index_cpe_jo3_b = -1 ! index for jo3_b in heating rate array + integer :: cpe_jo2_a_pbuf_index = -1 ! index in physics buffer for jo2_a heating rate + integer :: cpe_jo2_b_pbuf_index = -1 ! index in physics buffer for jo2_b heating rate + integer :: cpe_jo3_a_pbuf_index = -1 ! index in physics buffer for jo3_a heating rate + integer :: cpe_jo3_b_pbuf_index = -1 ! index in physics buffer for jo3_b heating rate + + ! Information needed to access CAM species state data + logical :: is_fixed_N2 = .false. ! indicates whether N2 concentrations are fixed + logical :: is_fixed_O = .false. ! indicates whether O concentrations are fixed + logical :: is_fixed_O2 = .false. ! indicates whether O2 concentrations are fixed + logical :: is_fixed_O3 = .false. ! indicates whether O3 concentrations are fixed + logical :: is_fixed_NO = .false. ! indicates whether NO concentrations are fixed + integer :: index_N2 = 0 ! index for N2 in concentration array + integer :: index_O = 0 ! index for O in concentration array + integer :: index_O2 = 0 ! index for O2 in concentration array + integer :: index_O3 = 0 ! index for O3 in concentration array + integer :: index_NO = 0 ! index for NO in concentration array + + ! Information needed to access aerosol and cloud optical properties + logical :: do_aerosol = .false. ! indicates whether aerosol optical properties + ! are available and should be used in radiative + ! transfer calculations + logical :: do_clouds = .false. ! indicates whether cloud optical properties + ! should be calculated and used in radiative + ! transfer calculations + + ! Information needed to set extended-UV photo rates + logical :: do_euv = .false. ! Indicates whether to calculate + ! extended-UV photo rates + integer :: ion_rates_pbuf_index = 0 ! Index in physics buffer for + ! ionization rates + + ! Information needed to do special NO photolysis rate calculation + logical :: do_jno = .false. ! Indicates whether to calculate jno + integer :: jno_index = 0 ! Index in tuvx_ptr::photo_rates_ array for jno + + ! Cutoff solar zenith angle for doing photolysis rate calculations [degrees] + real(r8) :: max_sza = 0.0_r8 + + ! TODO how should these paths be set and communicated to this wrapper? + character(len=*), parameter :: wavelength_config_path = & + "data/grids/wavelength/cam.csv" + logical, parameter :: enable_diagnostics = .true. + + ! TUV-x calculator for each OMP thread + type :: tuvx_ptr + type(core_t), pointer :: core_ => null( ) ! TUV-x calculator + integer :: n_photo_rates_ = 0 ! number of photo reactions in TUV-x + integer :: n_euv_rates_ = 0 ! number of extreme-UV rates + integer :: n_special_rates_ = 0 ! number of special photo rates + integer :: n_wavelength_bins_ = 0 ! number of wavelength bins in TUV-x + type(map_t) :: photo_rate_map_ ! map between TUV-x and CAM + ! photo rate constant arrays + type(grid_updater_t) :: grids_(NUM_GRIDS) ! grid updaters + type(profile_updater_t) :: profiles_(NUM_PROFILES) ! profile updaters + type(radiator_updater_t) :: radiators_(NUM_RADIATORS) ! radiator updaters + real(r8) :: height_delta_(pver+1) ! change in height in each + ! vertical layer (km) + real(r8), allocatable :: wavelength_edges_(:) ! TUV-x wavelength bin edges (nm) + ! on the TUV-x wavelength grid + real(r8) :: et_flux_ms93_(NUM_BINS_MS93) ! extraterrestrial flux on the MS93 grid + ! [photon cm-2 nm-1 s-1] + end type tuvx_ptr + type(tuvx_ptr), allocatable :: tuvx_ptrs(:) + + ! Diagnostic photolysis rate constant output + type :: diagnostic_t + character(len=:), allocatable :: name_ ! Name of the output field + integer :: index_ ! index of the photolysis rate constant from TUV-x + end type diagnostic_t + type(diagnostic_t), allocatable :: diagnostics(:) + + ! namelist options + character(len=cl) :: tuvx_config_path = 'NONE' ! absolute path to TUVX configuration file + logical, protected :: tuvx_active = .false. + +!================================================================================================ +contains +!================================================================================================ + + !----------------------------------------------------------------------- + ! registers fields in the physics buffer + !----------------------------------------------------------------------- + subroutine tuvx_register( ) + + use mo_jeuv, only : nIonRates + use physics_buffer, only : pbuf_add_field, dtype_r8 + use ppgrid, only : pcols ! maximum number of columns + + if( .not. tuvx_active ) return + + ! add photo-ionization rates to physics buffer for WACCMX Ionosphere module + call pbuf_add_field( 'IonRates', 'physpkg', dtype_r8, (/ pcols, pver, nIonRates /), & + ion_rates_pbuf_index ) ! Ionization rates for O+, O2+, N+, N2+, NO+ + + call pbuf_add_field( 'CPE_jO2a', 'global', dtype_r8, (/ pcols, pver /), cpe_jo2_a_pbuf_index ) + call pbuf_add_field( 'CPE_jO2b', 'global', dtype_r8, (/ pcols, pver /), cpe_jo2_b_pbuf_index ) + call pbuf_add_field( 'CPE_jO3a', 'global', dtype_r8, (/ pcols, pver /), cpe_jo3_a_pbuf_index ) + call pbuf_add_field( 'CPE_jO3b', 'global', dtype_r8, (/ pcols, pver /), cpe_jo3_b_pbuf_index ) + + end subroutine tuvx_register + +!================================================================================================ + + !----------------------------------------------------------------------- + ! read namelist options + !----------------------------------------------------------------------- + subroutine tuvx_readnl(nlfile) + +#ifdef HAVE_MPI + use mpi +#endif + use cam_abortutils, only : endrun + use cam_logfile, only : iulog ! log file output unit + use namelist_utils, only : find_group_name + use spmd_utils, only : mpicom, is_main_task => masterproc, & + main_task => masterprocid + + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + integer :: unitn, ierr + character(len=*), parameter :: subname = 'tuvx_readnl' + + ! =================== + ! Namelist definition + ! =================== + namelist /tuvx_opts/ tuvx_active, tuvx_config_path + + ! ============= + ! Read namelist + ! ============= + if (is_main_task) then + open( newunit=unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'tuvx_opts', status=ierr) + if (ierr == 0) then + read(unitn, tuvx_opts, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading namelist') + end if + end if + close(unitn) + end if + + ! ============================ + ! Broadcast namelist variables + ! ============================ +#ifdef HAVE_MPI + call mpi_bcast(tuvx_config_path, len(tuvx_config_path), mpi_character, main_task, mpicom, ierr) + call mpi_bcast(tuvx_active, 1, mpi_logical, main_task, mpicom, ierr) +#endif + + if (tuvx_active .and. tuvx_config_path == 'NONE') then + call endrun(subname // ' : must set tuvx_config_path when TUV-X is active') + end if + + if (is_main_task) then + write(iulog,*) 'tuvx_readnl: tuvx_config_path = ', trim(tuvx_config_path) + write(iulog,*) 'tuvx_readnl: tuvx_active = ', tuvx_active + end if + + end subroutine tuvx_readnl + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Initializes TUV-x for photolysis calculations + !----------------------------------------------------------------------- + subroutine tuvx_init( photon_file, electron_file, max_solar_zenith_angle, pbuf2d ) + +#ifdef HAVE_MPI + use mpi +#endif + use cam_history, only : addfld + use cam_logfile, only : iulog ! log file output unit + use infnan, only : nan, assignment(=) + use mo_chem_utls, only : get_spc_ndx, get_inv_ndx + use mo_jeuv, only : neuv ! number of extreme-UV rates + use musica_assert, only : assert_msg, die_msg + use musica_config, only : config_t + use musica_mpi, only : musica_mpi_rank, & + musica_mpi_pack_size, & + musica_mpi_pack, & + musica_mpi_unpack + use musica_string, only : string_t, to_char + use physics_buffer, only : physics_buffer_desc + use physics_buffer, only : pbuf_set_field + use ppgrid, only : pcols ! maximum number of columns + use shr_const_mod, only : pi => shr_const_pi + use solar_irrad_data, only : has_spectrum + use spmd_utils, only : main_task => masterprocid, & + is_main_task => masterproc, & + mpicom + use tuvx_grid, only : grid_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_profile_warehouse, only : profile_warehouse_t + use tuvx_radiator_warehouse, only : radiator_warehouse_t + use time_manager, only : is_first_step + + character(len=*), intent(in) :: photon_file ! photon file used in extended-UV module setup + character(len=*), intent(in) :: electron_file ! electron file used in extended-UV module setup + real(r8), intent(in) :: max_solar_zenith_angle ! cutoff solar zenith angle for + ! photo rate calculations [degrees] + type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer + + character(len=*), parameter :: my_name = "TUV-x wrapper initialization" + class(core_t), pointer :: core + character, allocatable :: buffer(:) + type(string_t) :: config_path + type(config_t) :: tuvx_config, cam_config, map_config + type(map_t) :: map + class(grid_t), pointer :: height + class(grid_t), pointer :: wavelength + class(grid_warehouse_t), pointer :: cam_grids + class(profile_warehouse_t), pointer :: cam_profiles + class(radiator_warehouse_t), pointer :: cam_radiators + integer :: pack_size, pos, i_core, i_err + logical :: disable_aerosols, disable_clouds + type(string_t) :: required_keys(1), optional_keys(2) + logical, save :: is_initialized = .false. + + type(string_t), allocatable :: labels(:) + character(len=16) :: label + integer :: i + real(r8) :: nanval + + if( .not. tuvx_active ) return + if( is_initialized ) return + is_initialized = .true. + + nanval=nan + + call pbuf_set_field( pbuf2d, ion_rates_pbuf_index, nanval ) + + if( is_first_step( ) ) then + call pbuf_set_field( pbuf2d, cpe_jo2_a_pbuf_index, 0.0_r8 ) + call pbuf_set_field( pbuf2d, cpe_jo2_b_pbuf_index, 0.0_r8 ) + call pbuf_set_field( pbuf2d, cpe_jo3_a_pbuf_index, 0.0_r8 ) + call pbuf_set_field( pbuf2d, cpe_jo3_b_pbuf_index, 0.0_r8 ) + end if + + if( is_main_task ) write(iulog,*) "Beginning TUV-x Initialization" + + config_path = trim(tuvx_config_path) + + ! =============================== + ! CAM TUV-x configuration options + ! =============================== + required_keys(1) = "aliasing" + optional_keys(1) = "disable aerosols" + optional_keys(2) = "disable clouds" + +#ifndef HAVE_MPI + call assert_msg( 113937299, is_main_task, "Multiple tasks present without " & + //"MPI support enabled for TUV-x" ) +#endif + + ! =============================================================== + ! set the maximum solar zenith angle to calculate photo rates for + ! =============================================================== + max_sza = max_solar_zenith_angle + if( max_sza <= 0.0_r8 .or. max_sza > 180.0_r8 ) then + call die_msg( 723815691, "TUV-x max solar zenith angle must be between 0 and 180 degress" ) + end if + + ! ================================= + ! initialize the extended-UV module + ! ================================= + call initialize_euv( photon_file, electron_file, do_euv ) + + ! ========================================================================== + ! create the set of TUV-x grids and profiles that CAM will update at runtime + ! ========================================================================== + cam_grids => get_cam_grids( wavelength_config_path ) + cam_profiles => get_cam_profiles( cam_grids ) + cam_radiators => get_cam_radiators( cam_grids ) + + ! ================================================================== + ! construct a core and a map between TUV-x and CAM photolysis arrays + ! on the primary process and pack them onto an MPI buffer + ! ================================================================== + if( is_main_task ) then + call tuvx_config%from_file( config_path%to_char( ) ) + call tuvx_config%get( "__CAM options", cam_config, my_name ) + call assert_msg( 973680295, & + cam_config%validate( required_keys, optional_keys ), & + "Bad configuration for CAM TUV-x options." ) + call cam_config%get( "disable aerosols", disable_aerosols, my_name, & + default = .false. ) + call cam_config%get( "disable clouds", disable_clouds, my_name, & + default = .false. ) + call cam_config%get( "aliasing", map_config, my_name ) + core => core_t( config_path, cam_grids, cam_profiles, cam_radiators ) + call set_photo_rate_map( core, map_config, do_euv, do_jno, jno_index, map ) + pack_size = core%pack_size( mpicom ) + & + map%pack_size( mpicom ) + & + musica_mpi_pack_size( do_jno, mpicom ) + & + musica_mpi_pack_size( jno_index, mpicom ) + & + musica_mpi_pack_size( disable_aerosols, mpicom ) + & + musica_mpi_pack_size( disable_clouds, mpicom ) + allocate( buffer( pack_size ) ) + pos = 0 + call core%mpi_pack( buffer, pos, mpicom ) + call map%mpi_pack( buffer, pos, mpicom ) + call musica_mpi_pack( buffer, pos, do_jno, mpicom ) + call musica_mpi_pack( buffer, pos, jno_index, mpicom ) + call musica_mpi_pack( buffer, pos, disable_aerosols, mpicom ) + call musica_mpi_pack( buffer, pos, disable_clouds, mpicom ) + deallocate( core ) + end if + +#ifdef HAVE_MPI + ! ==================================================== + ! broadcast the core and map data to all MPI processes + ! ==================================================== + call mpi_bcast( pack_size, 1, MPI_INTEGER, main_task, mpicom, i_err ) + if( i_err /= MPI_SUCCESS ) then + write(iulog,*) "TUV-x MPI int bcast error" + call mpi_abort( mpicom, 1, i_err ) + end if + if( .not. is_main_task ) allocate( buffer( pack_size ) ) + call mpi_bcast( buffer, pack_size, MPI_CHARACTER, main_task, mpicom, i_err ) + if( i_err /= MPI_SUCCESS ) then + write(iulog,*) "TUV-x MPI char array bcast error" + call mpi_abort( mpicom, 1, i_err ) + end if +#endif + + ! ================================================================ + ! unpack the core and map for each OMP thread on every MPI process + ! ================================================================ + allocate( tuvx_ptrs( max_threads( ) ) ) + do i_core = 1, size( tuvx_ptrs ) + associate( tuvx => tuvx_ptrs( i_core ) ) + allocate( tuvx%core_ ) + pos = 0 + call tuvx%core_%mpi_unpack( buffer, pos, mpicom ) + call tuvx%photo_rate_map_%mpi_unpack( buffer, pos, mpicom ) + call musica_mpi_unpack( buffer, pos, do_jno, mpicom ) + call musica_mpi_unpack( buffer, pos, jno_index, mpicom ) + call musica_mpi_unpack( buffer, pos, disable_aerosols, mpicom ) + call musica_mpi_unpack( buffer, pos, disable_clouds, mpicom ) + + ! =================================================================== + ! Set up connections between CAM and TUV-x input data for each thread + ! =================================================================== + call create_updaters( tuvx, cam_grids, cam_profiles, cam_radiators, & + disable_aerosols, disable_clouds ) + + ! ============================================ + ! Save the number of photolysis rate constants + ! ============================================ + tuvx%n_photo_rates_ = tuvx%core_%number_of_photolysis_reactions( ) + if( do_euv ) tuvx%n_euv_rates_ = neuv + if( do_jno ) tuvx%n_special_rates_ = tuvx%n_special_rates_ + 1 + + end associate + end do + deallocate( cam_grids ) + deallocate( cam_profiles ) + deallocate( cam_radiators ) + + ! ============================================= + ! Get index info for CAM species concentrations + ! ============================================= + index_N2 = get_inv_ndx( 'N2' ) + is_fixed_N2 = index_N2 > 0 + if( .not. is_fixed_N2 ) index_N2 = get_spc_ndx( 'N2' ) + index_O = get_inv_ndx( 'O' ) + is_fixed_O = index_O > 0 + if( .not. is_fixed_O ) index_O = get_spc_ndx( 'O' ) + index_O2 = get_inv_ndx( 'O2' ) + is_fixed_O2 = index_O2 > 0 + if( .not. is_fixed_O2 ) index_O2 = get_spc_ndx( 'O2' ) + index_O3 = get_inv_ndx( 'O3' ) + is_fixed_O3 = index_O3 > 0 + if( .not. is_fixed_O3 ) index_O3 = get_spc_ndx( 'O3' ) + index_NO = get_inv_ndx( 'NO' ) + is_fixed_NO = index_NO > 0 + if( .not. is_fixed_NO ) index_NO = get_spc_ndx( 'NO' ) + + ! ==================================================== + ! make sure extraterrestrial flux values are available + ! ==================================================== + call assert_msg( 170693514, has_spectrum, & + "Solar irradiance spectrum needed for TUV-x" ) + + ! ============================================ + ! set up diagnostic output of photolysis rates + ! ============================================ + call initialize_diagnostics( tuvx_ptrs( 1 ) ) + + ! =============================== + ! set up map to CAM heating rates + ! =============================== + labels = tuvx_ptrs(1)%core_%heating_rate_labels( ) + do i = 1, size( labels ) + label = trim( labels( i )%to_char( ) ) + select case( label ) + case( 'jo2_a' ) + index_cpe_jo2_a = i + number_of_heating_rates = number_of_heating_rates + 1 + call addfld('CPE_jO2a',(/ 'lev' /), 'A', 'joules sec-1', & + trim(label)//' chemical potential energy') + case( 'jo2_b' ) + index_cpe_jo2_b = i + number_of_heating_rates = number_of_heating_rates + 1 + call addfld('CPE_jO2b',(/ 'lev' /), 'A', 'joules sec-1', & + trim(label)//' chemical potential energy') + case( 'jo3_a' ) + index_cpe_jo3_a = i + number_of_heating_rates = number_of_heating_rates + 1 + call addfld('CPE_jO3a',(/ 'lev' /), 'A', 'joules sec-1', & + trim(label)//' chemical potential energy') + case( 'jo3_b' ) + index_cpe_jo3_b = i + number_of_heating_rates = number_of_heating_rates + 1 + call addfld('CPE_jO3b',(/ 'lev' /), 'A', 'joules sec-1', & + trim(label)//' chemical potential energy') + end select + end do + call assert_msg( 398372957, & + number_of_heating_rates == size( labels ), & + "TUV-x heating rate mismatch. Expected "// & + trim( to_char( size( labels ) )// & + " rates, but only matched "// & + trim( to_char( number_of_heating_rates ) )//"." ) ) + + if( is_main_task ) call log_initialization( labels ) + + end subroutine tuvx_init + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Updates TUV-x profiles that depend on time but not space + !----------------------------------------------------------------------- + subroutine tuvx_timestep_init( ) + + integer :: i_thread + + if( .not. tuvx_active ) return + + do i_thread = 1, size( tuvx_ptrs ) + associate( tuvx => tuvx_ptrs( i_thread ) ) + call set_et_flux( tuvx ) + end associate + end do + + end subroutine tuvx_timestep_init + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Calculates and returns photolysis rate constants + !----------------------------------------------------------------------- + subroutine tuvx_get_photo_rates( state, pbuf, ncol, lchnk, height_mid, & + height_int, temperature_mid, surface_temperature, fixed_species_conc, & + species_vmr, exo_column_conc, surface_albedo, solar_zenith_angle, & + earth_sun_distance, pressure_delta, cloud_fraction, liquid_water_content, & + photolysis_rates ) + + use cam_history, only : outfld + use cam_logfile, only : iulog ! log info output unit + use chem_mods, only : phtcnt, & ! number of photolysis reactions + gas_pcnst, & ! number of non-fixed species + nfs, & ! number of fixed species + nabscol ! number of absorbing species (radiators) + use physics_types, only : physics_state + use physics_buffer, only : physics_buffer_desc + use physics_buffer, only : pbuf_get_field + use ppgrid, only : pcols ! maximum number of columns + use shr_const_mod, only : pi => shr_const_pi + use spmd_utils, only : main_task => masterprocid, & + is_main_task => masterproc, & + mpicom + + type(physics_state), target, intent(in) :: state + type(physics_buffer_desc), pointer, intent(inout) :: pbuf(:) + integer, intent(in) :: ncol ! number of active columns on this thread + integer, intent(in) :: lchnk ! identifier for this thread + real(r8), intent(in) :: height_mid(ncol,pver) ! height at mid-points (km) + real(r8), intent(in) :: height_int(ncol,pver+1) ! height at interfaces (km) + real(r8), intent(in) :: temperature_mid(pcols,pver) ! midpoint temperature (K) + real(r8), intent(in) :: surface_temperature(pcols) ! surface temperature (K) + real(r8), intent(in) :: fixed_species_conc(ncol,pver,max(1,nfs)) ! fixed species densities + ! (molecule cm-3) + real(r8), intent(in) :: species_vmr(ncol,pver,max(1,gas_pcnst)) ! species volume mixing + ! ratios (mol mol-1) + real(r8), intent(in) :: exo_column_conc(ncol,0:pver,max(1,nabscol)) ! layer column densities + ! (molecule cm-2) + real(r8), intent(in) :: surface_albedo(pcols) ! surface albedo (unitless) + real(r8), intent(in) :: solar_zenith_angle(ncol) ! solar zenith angle (radians) + real(r8), intent(in) :: earth_sun_distance ! Earth-Sun distance (AU) + real(r8), intent(in) :: pressure_delta(pcols,pver) ! pressure delta about midpoints (Pa) + real(r8), intent(in) :: cloud_fraction(ncol,pver) ! cloud fraction (unitless) + real(r8), intent(in) :: liquid_water_content(ncol,pver) ! liquid water content (kg/kg) + real(r8), intent(inout) :: photolysis_rates(ncol,pver,phtcnt) ! photolysis rate + ! constants (1/s) + + integer :: i_col ! column index + integer :: i_level ! vertical level index + real(r8) :: sza ! solar zenith angle [degrees] + real(r8) :: cpe_rates(ncol,pverp+1,number_of_heating_rates) ! heating rates from TUV-x + real(r8), pointer :: cpe_jo2_a(:,:) ! heating rate for jo2_a in physics buffer + real(r8), pointer :: cpe_jo2_b(:,:) ! heating rate for jo2_b in physics buffer + real(r8), pointer :: cpe_jo3_a(:,:) ! heating rate for jo3_a in physics buffer + real(r8), pointer :: cpe_jo3_b(:,:) ! heating rate for jo3_b in physics buffer + + ! working arrays + real(r8), allocatable :: photo_rates(:,:,:) ! calculated photo rate constants (column, level, reaction) [s-1] + real(r8), allocatable :: optical_depth(:,:,:) ! aerosol optical depth (column, level, wavelength) [unitless] + real(r8), allocatable :: single_scattering_albedo(:,:,:) ! aerosol single scattering albedo (column, level, wavelength) [unitless] + real(r8), allocatable :: asymmetry_factor(:,:,:) ! aerosol asymmetry factor (column, level, wavelength) [unitless] + + if( .not. tuvx_active ) return + + call pbuf_get_field(pbuf, cpe_jo2_a_pbuf_index, cpe_jo2_a) + call pbuf_get_field(pbuf, cpe_jo2_b_pbuf_index, cpe_jo2_b) + call pbuf_get_field(pbuf, cpe_jo3_a_pbuf_index, cpe_jo3_a) + call pbuf_get_field(pbuf, cpe_jo3_b_pbuf_index, cpe_jo3_b) + + cpe_rates(:,:,:) = 0.0_r8 + cpe_jo2_a(:,:) = 0.0_r8 + cpe_jo2_b(:,:) = 0.0_r8 + cpe_jo3_a(:,:) = 0.0_r8 + cpe_jo3_b(:,:) = 0.0_r8 + + associate( tuvx => tuvx_ptrs( thread_id( ) ) ) + + allocate( photo_rates( pcols, pver+2, tuvx%n_photo_rates_ + tuvx%n_euv_rates_ & + + tuvx%n_special_rates_ ) ) + allocate( optical_depth( pcols, pver+1, tuvx%n_wavelength_bins_ ) ) + allocate( single_scattering_albedo( pcols, pver+1, tuvx%n_wavelength_bins_ ) ) + allocate( asymmetry_factor( pcols, pver+1, tuvx%n_wavelength_bins_ ) ) + photo_rates(:,:,:) = 0.0_r8 + + ! ============================================== + ! set aerosol optical properties for all columns + ! ============================================== + call get_aerosol_optical_properties( tuvx, state, pbuf, optical_depth, & + single_scattering_albedo, asymmetry_factor ) + + do i_col = 1, ncol + + ! =================================== + ! skip columns in near total darkness + ! =================================== + sza = solar_zenith_angle(i_col) * 180.0_r8 / pi + if( sza < 0.0_r8 .or. sza > max_sza ) cycle + + ! =================== + ! update grid heights + ! =================== + call set_heights( tuvx, i_col, ncol, height_mid, height_int ) + + ! ======================================= + ! set conditions for this column in TUV-x + ! ======================================= + call set_temperatures( tuvx, i_col, temperature_mid, surface_temperature ) + call set_surface_albedo( tuvx, i_col, surface_albedo ) + call set_radiator_profiles( tuvx, i_col, ncol, fixed_species_conc, & + species_vmr, exo_column_conc, & + pressure_delta(1:ncol,:), cloud_fraction, & + liquid_water_content, optical_depth, & + single_scattering_albedo, asymmetry_factor) + + ! =================================================== + ! Calculate photolysis rate constants for this column + ! =================================================== + call tuvx%core_%run( solar_zenith_angle = sza, & + earth_sun_distance = earth_sun_distance, & + photolysis_rate_constants = & + photo_rates(i_col,:,1:tuvx%n_photo_rates_), & + heating_rates = cpe_rates(i_col,:,:) ) + + ! ============================== + ! Calculate the extreme-UV rates + ! ============================== + if( do_euv ) then + associate( euv_begin => tuvx%n_photo_rates_ + 1, & + euv_end => tuvx%n_photo_rates_ + tuvx%n_euv_rates_ ) + call calculate_euv_rates( sza, & + fixed_species_conc(i_col,:,:), & + species_vmr(i_col,:,:), & + height_mid(i_col,:), & + height_int(i_col,:), & + photo_rates(i_col,2:pver+1,euv_begin:euv_end) ) + end associate + end if + + ! ============================= + ! Calculate special photo rates + ! ============================= + if( do_jno ) then + call calculate_jno( sza, & + tuvx%et_flux_ms93_, & + fixed_species_conc(i_col,:,:), & + species_vmr(i_col,:,:), & + height_int(i_col,:), & + photo_rates(i_col,2:pver+1,jno_index) ) + end if + end do + + ! ===================== + ! Filter negative rates + ! ===================== + photo_rates(:,:,:) = max( 0.0_r8, photo_rates(:,:,:) ) + + ! ============================================ + ! Return the photolysis rates on the CAM grids + ! ============================================ + do i_col = 1, ncol + do i_level = 1, pver + call tuvx%photo_rate_map_%apply( photo_rates(i_col,pver-i_level+2,:), & + photolysis_rates(i_col,i_level,:) ) + end do + end do + + call output_diagnostics( tuvx, ncol, lchnk, photo_rates ) + + end associate + + if (index_cpe_jo2_a>0) then + do i_level = 1, pver + cpe_jo2_a(:ncol,i_level) = & + 0.5_r8 * ( cpe_rates(:ncol,pver-i_level+2,index_cpe_jo2_a) & + + cpe_rates(:ncol,pver-i_level+1,index_cpe_jo2_a) ) + end do + call outfld('CPE_jO2a', cpe_jo2_a(:ncol,:), ncol, lchnk ) + end if + if (index_cpe_jo2_b>0) then + do i_level = 1, pver + cpe_jo2_b(:ncol,i_level) = & + 0.5_r8 * ( cpe_rates(:ncol,pver-i_level+2,index_cpe_jo2_b) & + + cpe_rates(:ncol,pver-i_level+1,index_cpe_jo2_b) ) + end do + call outfld('CPE_jO2b', cpe_jo2_b(:ncol,:), ncol, lchnk ) + end if + + if (index_cpe_jo3_a>0) then + do i_level = 1, pver + cpe_jo3_a(:ncol,i_level) = & + 0.5_r8 * ( cpe_rates(:ncol,pver-i_level+2,index_cpe_jo3_a) & + + cpe_rates(:ncol,pver-i_level+1,index_cpe_jo3_a) ) + end do + call outfld('CPE_jO3a', cpe_jo3_a(:ncol,:), ncol, lchnk ) + end if + if (index_cpe_jo3_b>0) then + do i_level = 1, pver + cpe_jo3_b(:ncol,i_level) = & + 0.5_r8 * ( cpe_rates(:ncol,pver-i_level+2,index_cpe_jo3_b) & + + cpe_rates(:ncol,pver-i_level+1,index_cpe_jo3_b) ) + end do + call outfld('CPE_jO3b', cpe_jo3_b(:ncol,:), ncol, lchnk ) + end if + + end subroutine tuvx_get_photo_rates + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Cleans up memory associated with TUV-x calculators + !----------------------------------------------------------------------- + subroutine tuvx_finalize( ) + + integer :: i_core + + if( allocated( tuvx_ptrs ) ) then + do i_core = 1, size( tuvx_ptrs ) + associate( tuvx => tuvx_ptrs( i_core ) ) + if( associated( tuvx%core_ ) ) deallocate( tuvx%core_ ) + end associate + end do + end if + + end subroutine tuvx_finalize + +!================================================================================================ +!================================================================================================ +! +! Support functions +! +!================================================================================================ +!================================================================================================ + + !----------------------------------------------------------------------- + ! Returns the id of the current OpenMP thread, or 1 if not + ! using OpenMP (1 <= id <= max_threads()) + !----------------------------------------------------------------------- + integer function thread_id( ) +#ifdef _OPENMP + use omp_lib, only : omp_get_thread_num +#endif + +#ifdef _OPENMP + thread_id = 1 + omp_get_thread_num( ) +#else + thread_id = 1 +#endif + + end function thread_id + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Returns the number of threads available for calculations at + ! runtime + !----------------------------------------------------------------------- + integer function max_threads( ) +#ifdef _OPENMP + use omp_lib, only : omp_get_max_threads +#endif + +#ifdef _OPENMP + max_threads = omp_get_max_threads( ) +#else + max_threads = 1 +#endif + + end function max_threads + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Prints initialization conditions to the log file + !----------------------------------------------------------------------- + subroutine log_initialization( heating_rate_labels ) + + use cam_logfile, only : iulog ! log info output unit + use musica_string, only : to_char + use spmd_utils, only : main_task => masterprocid, & + is_main_task => masterproc + + type(string_t), intent(in) :: heating_rate_labels(:) ! heating rate labels + + integer :: i + + if( is_main_task ) then + write(iulog,*) "Initialized TUV-x" +#ifdef HAVE_MPI + write(iulog,*) " - with MPI support on task "//trim( to_char( main_task ) ) +#else + write(iulog,*) " - without MPI support" +#endif +#ifdef _OPENMP + write(iulog,*) " - with OpenMP support for "// & + trim( to_char( max_threads( ) ) )//" threads, on thread " & + //trim( to_char( thread_id( ) ) ) +#else + write(iulog,*) " - without OpenMP support" +#endif + write(iulog,*) " - with configuration file: '"//trim( tuvx_config_path )//"'" + if( do_aerosol ) then + write(iulog,*) " - with on-line aerosols" + else + write(iulog,*) " - without on-line aerosols" + end if + if( do_clouds ) then + write(iulog,*) " - with on-line clouds" + else + write(iulog,*) " - without on-line clouds" + end if + if( index_N2 > 0 ) write(iulog,*) " - including N2" + if( index_O > 0 ) write(iulog,*) " - including O" + if( index_O2 > 0 ) write(iulog,*) " - including O2" + if( index_O3 > 0 ) write(iulog,*) " - including O3" + if( index_NO > 0 ) write(iulog,*) " - including NO" + if( do_euv ) write(iulog,*) " - doing Extreme-UV calculations" + if( do_jno ) write(iulog,*) " - including special jno rate calculation" + write(iulog,*) " - max solar zenith angle [degrees]:", max_sza + if( size( heating_rate_labels ) > 0 ) then + write(iulog,*) " - with heating rates:" + do i = 1, size( heating_rate_labels ) + write(iulog,*) " - "//trim( heating_rate_labels( i )%to_char( ) ) + end do + end if + end if + + end subroutine log_initialization + +!================================================================================================ + + !----------------------------------------------------------------------- + ! initializes the external extreme-UV module + !----------------------------------------------------------------------- + subroutine initialize_euv( photon_file, electron_file, do_euv ) + + use chem_mods, only : phtcnt ! number of CAM-Chem photolysis reactions + use mo_jeuv, only : jeuv_init ! extreme-UV initialization + + character(len=*), intent(in) :: photon_file ! photon file used in extended-UV module + ! setup + character(len=*), intent(in) :: electron_file ! electron file used in extended-UV + ! module setup + logical, intent(out) :: do_euv ! indicates whether extreme-UV + ! calculations are needed + + integer, allocatable :: euv_index_map(:) + + allocate( euv_index_map( phtcnt ) ) + euv_index_map(:) = 0 + call jeuv_init( photon_file, electron_file, euv_index_map ) + do_euv = any( euv_index_map(:) > 0 ) + + end subroutine initialize_euv + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Registers fields for diagnostic output + !----------------------------------------------------------------------- + subroutine initialize_diagnostics( this ) + + use cam_history, only : addfld + use musica_assert, only : assert + use musica_string, only : string_t + + type(tuvx_ptr), intent(in) :: this + + type(string_t), allocatable :: labels(:), all_labels(:) + integer :: i_label + + if( .not. enable_diagnostics ) then + allocate( diagnostics( 0 ) ) + return + end if + + ! ========================================================== + ! add output for specific photolysis reaction rate constants + ! ========================================================== + labels = this%core_%photolysis_reaction_labels( ) + allocate( all_labels( size( labels ) + this%n_special_rates_ ) ) + all_labels( 1 : size( labels ) ) = labels(:) + i_label = size( labels ) + 1 + if( do_jno ) then + all_labels( i_label ) = "jno" + i_label = i_label + 1 + end if + call assert( 522515214, i_label == size( all_labels ) + 1 ) + allocate( diagnostics( size( all_labels ) ) ) + do i_label = 1, size( all_labels ) + diagnostics( i_label )%name_ = trim( all_labels( i_label )%to_char( ) ) + diagnostics( i_label )%index_ = i_label + call addfld( "tuvx_"//diagnostics( i_label )%name_, (/ 'lev' /), 'A', 'sec-1', & + 'photolysis rate constant' ) + end do + + end subroutine initialize_diagnostics + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Sets up a map between the TUV-x and CAM photolysis rate arrays + !----------------------------------------------------------------------- + subroutine set_photo_rate_map( core, config, do_euv, do_jno, jno_index, map ) + + use cam_logfile, only : iulog ! log info output unit + use chem_mods, only : phtcnt, & ! number of photolysis reactions + rxt_tag_lst ! labels for all chemical reactions + ! NOTE photolysis reactions are + ! expected to appear first + use mo_jeuv, only : neuv ! number of extreme-UV rates + use musica_config, only : config_t + use musica_string, only : string_t + + type(core_t), intent(in) :: core ! TUV-x core + type(config_t), intent(inout) :: config ! CAM<->TUV-x map configuration + logical, intent(in) :: do_euv ! indicates whether to include + ! extreme-UV rates in the mapping + logical, intent(out) :: do_jno ! indicates whether jno should be + ! calculated + integer, intent(out) :: jno_index ! index for jno in source photo + ! rate array + type(map_t), intent(out) :: map + + integer :: i_label, i_start, i_end + type(string_t) :: str_label + type(string_t), allocatable :: tuvx_labels(:), euv_labels(:), special_labels(:), & + all_labels(:), cam_labels(:) + + ! ================== + ! MOZART photo rates + ! ================== + allocate( cam_labels( phtcnt ) ) + do i_label = 1, phtcnt + cam_labels( i_label ) = trim( rxt_tag_lst( i_label ) ) + end do + + ! ================= + ! TUV-x photo rates + ! ================= + tuvx_labels = core%photolysis_reaction_labels( ) + + ! ====================== + ! Extreme-UV photo rates + ! ====================== + if( do_euv ) then + allocate( euv_labels( neuv ) ) + do i_label = 1, neuv + str_label = i_label + euv_labels( i_label ) = "jeuv_"//str_label + end do + else + allocate( euv_labels(0) ) + end if + + ! =============================== + ! Special photo rate calculations + ! =============================== + do_jno = .false. + jno_index = 0 + do i_label = 1, size( cam_labels ) + if( cam_labels( i_label ) == "jno" ) then + do_jno = .true. + exit + end if + end do + if( do_jno ) then + allocate( special_labels(1) ) + special_labels(1) = "jno" + else + allocate( special_labels(0) ) + end if + + ! ========================== + ! Combine photo rate sources + ! ========================== + allocate( all_labels( size( tuvx_labels ) + size( euv_labels ) + & + size( special_labels ) ) ) + i_end = 0 + if( size( tuvx_labels ) > 0 ) then + i_start = i_end + 1 + i_end = i_start + size( tuvx_labels ) - 1 + all_labels( i_start : i_end ) = tuvx_labels(:) + end if + if( size( euv_labels ) > 0 ) then + i_start = i_end + 1 + i_end = i_start + size( euv_labels ) - 1 + all_labels( i_start : i_end ) = euv_labels(:) + end if + if( size( special_labels ) > 0 ) then + i_start = i_end + 1 + i_end = i_start + size( special_labels ) - 1 + all_labels( i_start : i_end ) = special_labels(:) + jno_index = i_start + end if + + ! ========== + ! Create map + ! ========== + map = map_t( config, all_labels, cam_labels ) + write(iulog,*) + write(iulog,*) "TUV-x --> CAM-Chem photolysis rate constant map" + call map%print( all_labels, cam_labels, iulog ) + + end subroutine set_photo_rate_map + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Outputs diagnostic information for the current time step + !----------------------------------------------------------------------- + subroutine output_diagnostics( this, ncol, lchnk, photo_rates ) + + use cam_history, only : outfld + + type(tuvx_ptr), intent(in) :: this + integer, intent(in) :: ncol ! number of active columns on this thread + integer, intent(in) :: lchnk ! identifier for this thread + real(r8), intent(in) :: photo_rates(:,:,:) ! photo rate constants + ! (column, level, reaction) [s-1] + + integer :: i_diag + + if( .not. enable_diagnostics ) return + + do i_diag = 1, size( diagnostics ) + associate( diag => diagnostics( i_diag ) ) + call outfld( "tuvx_"//diag%name_, photo_rates(:ncol,pver+1:2:-1,diag%index_), & + ncol, lchnk ) + end associate + end do + + end subroutine output_diagnostics + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Creates and loads a grid warehouse with grids that CAM will + ! update at runtime + !----------------------------------------------------------------------- + function get_cam_grids( wavelength_path ) result( grids ) + + use musica_assert, only : assert_msg + use musica_config, only : config_t + use tuvx_grid, only : grid_t + use tuvx_grid_factory, only : grid_builder + use tuvx_grid_from_host, only : grid_from_host_t + use tuvx_grid_warehouse, only : grid_warehouse_t + + character(len=*), intent(in) :: wavelength_path ! path to the wavelength data file + class(grid_warehouse_t), pointer :: grids ! collection of grids to be updated by CAM + + character(len=*), parameter :: my_name = "CAM grid creator" + class(grid_t), pointer :: host_grid + type(config_t) :: config + + grids => grid_warehouse_t( ) + + ! ========================= + ! heights above the surface + ! ========================= + host_grid => grid_from_host_t( "height", "km", pver+1 ) + call grids%add( host_grid ) + deallocate( host_grid ) + + ! ============================================ + ! wavelengths (will not be updated at runtime) + ! ============================================ + call config%empty( ) + call config%add( "type", "from csv file", my_name ) + call config%add( "name", "wavelength", my_name ) + call config%add( "units", "nm", my_name ) + call config%add( "file path", wavelength_path, my_name ) + host_grid => grid_builder( config ) + call grids%add( host_grid ) + deallocate( host_grid ) + + end function get_cam_grids + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Creates and loads a profile warehouse with profiles that CAM + ! will update at runtime + !----------------------------------------------------------------------- + function get_cam_profiles( grids ) result( profiles ) + + use tuvx_grid, only : grid_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_profile_from_host, only : profile_from_host_t + use tuvx_profile_warehouse, only : profile_warehouse_t + + class(grid_warehouse_t), intent(in) :: grids ! CAM grids used in TUV-x + class(profile_warehouse_t), pointer :: profiles ! collection of profiles to be updated by CAM + + class(profile_from_host_t), pointer :: host_profile + class(grid_t), pointer :: height, wavelength + + profiles => profile_warehouse_t( ) + height => grids%get_grid( "height", "km" ) + wavelength => grids%get_grid( "wavelength", "nm" ) + + ! ================================== + ! Temperature profile on height grid + ! ================================== + host_profile => profile_from_host_t( "temperature", "K", height%size( ) ) + call profiles%add( host_profile ) + deallocate( host_profile ) + + ! ================================= + ! Surface albedo on wavelength grid + ! ================================= + host_profile => profile_from_host_t( "surface albedo", "none", wavelength%size( ) ) + call profiles%add( host_profile ) + deallocate( host_profile ) + + ! ======================================== + ! Extraterrestrial flux on wavelength grid + ! ======================================== + host_profile => profile_from_host_t( "extraterrestrial flux", "photon cm-2 s-1", & + wavelength%size( ) ) + call profiles%add( host_profile ) + deallocate( host_profile ) + + ! =========== + ! Air profile + ! =========== + host_profile => profile_from_host_t( "air", "molecule cm-3", height%size( ) ) + call profiles%add( host_profile ) + deallocate( host_profile ) + + ! ========== + ! O3 profile + ! ========== + host_profile => profile_from_host_t( "O3", "molecule cm-3", height%size( ) ) + call profiles%add( host_profile ) + deallocate( host_profile ) + + ! ========== + ! O2 profile + ! ========== + host_profile => profile_from_host_t( "O2", "molecule cm-3", height%size( ) ) + call profiles%add( host_profile ) + deallocate( host_profile ) + + deallocate( height ) + deallocate( wavelength ) + + end function get_cam_profiles + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Creates and loads a radiator warehouse with radiators that CAM will + ! update at runtime + !----------------------------------------------------------------------- + function get_cam_radiators( grids ) result( radiators ) + + use tuvx_grid, only : grid_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_radiator_from_host, only : radiator_from_host_t + use tuvx_radiator_warehouse, only : radiator_warehouse_t + + type(grid_warehouse_t), intent(in) :: grids ! CAM grids used in TUV-x + class(radiator_warehouse_t), pointer :: radiators ! collection of radiators to be updated by CAM + + class(radiator_from_host_t), pointer :: host_radiator + class(grid_t), pointer :: height, wavelength + + radiators => radiator_warehouse_t( ) + height => grids%get_grid( "height", "km" ) + wavelength => grids%get_grid( "wavelength", "nm" ) + + ! ================ + ! Aerosol radiator + ! ================ + host_radiator => radiator_from_host_t( "aerosol", height, wavelength ) + call radiators%add( host_radiator ) + deallocate( host_radiator ) + + ! ============== + ! Cloud radiator + ! ============== + host_radiator => radiator_from_host_t( "clouds", height, wavelength ) + call radiators%add( host_radiator ) + deallocate( host_radiator ) + + deallocate( height ) + deallocate( wavelength ) + + end function get_cam_radiators + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Forms connections between CAM and TUV-x data structures + ! + ! - creates updaters for each grid and profile that CAM will use + ! to update TUV-x at each timestep. + ! - allocates working arrays when needed for interpolation between + ! grids + ! - initializes CAM modules if necessary and sets parameters needed + ! for runtime access of CAM data + ! + !----------------------------------------------------------------------- + subroutine create_updaters( this, grids, profiles, radiators, disable_aerosols, & + disable_clouds ) + + use ppgrid, only : pcols ! maximum number of columns + use rad_constituents, only : rad_cnst_get_info + use musica_assert, only : assert + use tuvx_grid, only : grid_t + use tuvx_grid_from_host, only : grid_updater_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_profile, only : profile_t + use tuvx_profile_from_host, only : profile_updater_t + use tuvx_profile_warehouse, only : profile_warehouse_t + use tuvx_radiator, only : radiator_t + use tuvx_radiator_from_host, only : radiator_updater_t + use tuvx_radiator_warehouse, only : radiator_warehouse_t + + class(tuvx_ptr), intent(inout) :: this + class(grid_warehouse_t), intent(in) :: grids + class(profile_warehouse_t), intent(in) :: profiles + class(radiator_warehouse_t), intent(in) :: radiators + logical, intent(in) :: disable_aerosols + logical, intent(in) :: disable_clouds + + class(grid_t), pointer :: height, wavelength + class(profile_t), pointer :: host_profile + class(radiator_t), pointer :: host_radiator + integer :: n_modes + logical :: found + + ! ============= + ! Grid updaters + ! ============= + + height => grids%get_grid( "height", "km" ) + this%grids_( GRID_INDEX_HEIGHT ) = this%core_%get_updater( height, found ) + call assert( 213798815, found ) + + ! ============================================ + ! wavelength grid cannot be updated at runtime + ! ============================================ + wavelength => grids%get_grid( "wavelength", "nm" ) + this%n_wavelength_bins_ = wavelength%size( ) + allocate( this%wavelength_edges_( this%n_wavelength_bins_ + 1 ) ) + this%wavelength_edges_(:) = wavelength%edge_(:) + + deallocate( height ) + deallocate( wavelength ) + + ! ================ + ! Profile updaters + ! ================ + + host_profile => profiles%get_profile( "temperature", "K" ) + this%profiles_( PROFILE_INDEX_TEMPERATURE ) = this%core_%get_updater( host_profile, found ) + call assert( 418735162, found ) + deallocate( host_profile ) + + host_profile => profiles%get_profile( "surface albedo", "none" ) + this%profiles_( PROFILE_INDEX_ALBEDO ) = this%core_%get_updater( host_profile, found ) + call assert( 720785186, found ) + deallocate( host_profile ) + + host_profile => profiles%get_profile( "extraterrestrial flux", "photon cm-2 s-1" ) + this%profiles_( PROFILE_INDEX_ET_FLUX ) = this%core_%get_updater( host_profile, found ) + call assert( 550628282, found ) + deallocate( host_profile ) + + host_profile => profiles%get_profile( "air", "molecule cm-3" ) + this%profiles_( PROFILE_INDEX_AIR ) = this%core_%get_updater( host_profile, found ) + call assert( 380471378, found ) + deallocate( host_profile ) + + host_profile => profiles%get_profile( "O3", "molecule cm-3" ) + this%profiles_( PROFILE_INDEX_O3 ) = this%core_%get_updater( host_profile, found ) + call assert( 210314474, found ) + deallocate( host_profile ) + + host_profile => profiles%get_profile( "O2", "molecule cm-3" ) + this%profiles_( PROFILE_INDEX_O2 ) = this%core_%get_updater( host_profile, found ) + call assert( 105165970, found ) + deallocate( host_profile ) + + ! ================= + ! radiator updaters + ! ================= + + ! ==================================================================== + ! determine if aerosol optical properties will be available, and if so + ! intialize the aerosol optics module + ! ==================================================================== + call rad_cnst_get_info( 0, nmodes = n_modes ) + if( n_modes > 0 .and. .not. do_aerosol .and. .not. disable_aerosols ) then + do_aerosol = .true. + ! TODO update to use new aerosol_optics class + ! call modal_aer_opt_init( ) + end if + host_radiator => radiators%get_radiator( "aerosol" ) + this%radiators_( RADIATOR_INDEX_AEROSOL ) = & + this%core_%get_updater( host_radiator, found ) + call assert( 675200430, found ) + nullify( host_radiator ) + + ! ===================================== + ! get an updater for the cloud radiator + ! ===================================== + do_clouds = .not. disable_clouds + host_radiator => radiators%get_radiator( "clouds" ) + this%radiators_( RADIATOR_INDEX_CLOUDS ) = & + this%core_%get_updater( host_radiator, found ) + call assert( 993715720, found ) + nullify( host_radiator ) + + end subroutine create_updaters + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Sets the height values in TUV-x for the given column + ! + ! NOTE: This function must be called before updating any profile data. + ! + ! CAM to TUV-x height grid mapping + ! + ! TUV-x heights are "bottom-up" and require atmospheric constituent + ! concentrations at interfaces. Therefore, CAM mid-points are used as + ! TUV-x grid interfaces, with an additional layer introduced between + ! the surface and the lowest CAM mid-point, and a layer at the + ! top of the TUV-x grid to hold species densities above the top CAM + ! mid-point. + ! + ! ---- (interface) ===== (mid-point) + ! + ! CAM TUV-x + ! ------(top)------ i_int = 1 -------(top)------ i_int = pver + 2 + ! ************************ (exo values) ***************************** + ! ================== i_mid = pver + 1 + ! ================= i_mid = 1 ------------------ i_int = pver + 1 + ! ----------------- i_int = 2 ================== i_mid = pver + ! ------------------ i_int = pver + ! || + ! || || + ! || + ! ----------------- i_int = pver + ! ================= i_imd = pver ------------------ i_int = 2 + ! ================== i_mid = 1 + ! -----(ground)---- i_int = pver+1 -----(ground)----- i_int = 1 + ! + !----------------------------------------------------------------------- + subroutine set_heights( this, i_col, ncol, height_mid, height_int ) + + use ppgrid, only : pcols ! maximum number of columns + + class(tuvx_ptr), intent(inout) :: this ! TUV-x calculator + integer, intent(in) :: i_col ! column to set conditions for + integer, intent(in) :: ncol ! number of colums to calculated photolysis for + real(r8), intent(in) :: height_mid(ncol,pver) ! height above the surface at mid-points (km) + real(r8), intent(in) :: height_int(ncol,pver+1) ! height above the surface at interfaces (km) + + integer :: i_level + real(r8) :: edges(pver+2) + real(r8) :: mid_points(pver+1) + + edges(1) = height_int(i_col,pver+1) + edges(2:pver+1) = height_mid(i_col,pver:1:-1) + edges(pver+2) = height_int(i_col,1) + mid_points(1) = ( height_mid(i_col,pver) - height_int(i_col,pver+1) ) * 0.5_r8 & + + height_int(i_col,pver+1) + mid_points(2:pver) = height_int(i_col,pver:2:-1) + mid_points(pver+1) = 0.5_r8 * ( edges(pver+1) + edges(pver+2) ) + call this%grids_( GRID_INDEX_HEIGHT )%update( edges = edges, mid_points = mid_points ) + this%height_delta_(1:pver+1) = edges(2:pver+2) - edges(1:pver+1) + + end subroutine set_heights + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Sets the temperatures in TUV-x for the given column + ! + ! See description of `set_heights` for CAM <-> TUV-x vertical grid + ! mapping. + !----------------------------------------------------------------------- + subroutine set_temperatures( this, i_col, temperature_mid, surface_temperature ) + + use ppgrid, only : pcols ! maximum number of columns + + class(tuvx_ptr), intent(inout) :: this ! TUV-x calculator + integer, intent(in) :: i_col ! column to set conditions for + real(r8), intent(in) :: temperature_mid(pcols,pver) ! midpoint temperature (K) + real(r8), intent(in) :: surface_temperature(pcols) ! surface temperature (K) + + real(r8) :: edges(pver+2) + + edges(1) = surface_temperature(i_col) + edges(2:pver+1) = temperature_mid(i_col,pver:1:-1) + edges(pver+2) = temperature_mid(i_col,1) ! Use upper mid-point temperature for top edge + call this%profiles_( PROFILE_INDEX_TEMPERATURE )%update( edge_values = edges ) + + end subroutine set_temperatures + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Sets the surface albedo in TUV-x for the given column + ! + ! CAM uses a single value for surface albedo at all wavelengths + !----------------------------------------------------------------------- + subroutine set_surface_albedo( this, i_col, surface_albedo ) + + use ppgrid, only : pcols ! maximum number of columns + + class(tuvx_ptr), intent(inout) :: this ! TUV-x calculator + integer, intent(in) :: i_col ! column to set conditions for + real(r8), intent(in) :: surface_albedo(pcols) ! surface albedo (unitless) + + real(r8) :: albedos(this%n_wavelength_bins_ + 1) + + albedos(:) = surface_albedo( i_col ) + call this%profiles_( PROFILE_INDEX_ALBEDO )%update( & + edge_values = albedos(:) ) + + end subroutine set_surface_albedo + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Sets the extraterrestrial flux in TUV-x for the given column + ! + ! Extraterrestrial flux is read from data files and interpolated to the + ! TUV-x wavelength grid. CAM ET Flux values are multiplied by the + ! width of the wavelength bins to get the TUV-x units of photon cm-2 s-1 + ! + ! NOTE: TUV-x only uses mid-point values for ET Flux + !----------------------------------------------------------------------- + subroutine set_et_flux( this ) + + use mo_util, only : rebin + use solar_irrad_data, only : nbins, & ! number of wavelength bins + we, & ! wavelength bin edges + sol_etf ! extraterrestrial flux + ! (photon cm-2 nm-1 s-1) + + class(tuvx_ptr), intent(inout) :: this ! TUV-x calculator + + real(r8) :: et_flux_orig(nbins), tuv_mid_values(this%n_wavelength_bins_), & + tuv_edge_values(this%n_wavelength_bins_ + 1) + integer :: i_bin + + ! =============================================== + ! regrid normalized flux to TUV-x wavelength grid + !================================================ + et_flux_orig(:) = sol_etf(:) + call rebin( nbins, this%n_wavelength_bins_, we, this%wavelength_edges_, & + et_flux_orig, tuv_mid_values ) + + ! ======================================================== + ! convert normalized flux to flux on TUV-x wavelength grid + ! ======================================================== + tuv_mid_values(:) = tuv_mid_values(:) * & + ( this%wavelength_edges_(2:this%n_wavelength_bins_+1) - & + this%wavelength_edges_(1:this%n_wavelength_bins_) ) + + ! ==================================== + ! estimate unused edge values for flux + ! ==================================== + tuv_edge_values(1) = tuv_mid_values(1) - & + ( tuv_mid_values(2) - tuv_mid_values(1) ) * 0.5_r8 + do i_bin = 2, this%n_wavelength_bins_ + tuv_edge_values(i_bin) = tuv_mid_values(i_bin-1) + & + ( tuv_mid_values(i_bin) - tuv_mid_values(i_bin-1) ) * 0.5_r8 + end do + tuv_edge_values(this%n_wavelength_bins_+1) = & + tuv_mid_values(this%n_wavelength_bins_) + & + ( tuv_mid_values(this%n_wavelength_bins_) - & + tuv_mid_values(this%n_wavelength_bins_-1) ) * 0.5_r8 + + ! ============================ + ! update TUV-x ET flux profile + ! ============================ + call this%profiles_( PROFILE_INDEX_ET_FLUX )%update( & + mid_point_values = tuv_mid_values, & + edge_values = tuv_edge_values) + + ! ====================================================================== + ! rebin extraterrestrial flux to MS93 grid for use with jno calculations + ! ====================================================================== + call rebin( nbins, NUM_BINS_MS93, we, WAVELENGTH_EDGES_MS93, et_flux_orig, & + this%et_flux_ms93_ ) + + end subroutine set_et_flux + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Sets the profiles of optically active atmospheric constituents + ! in TUV-x for the given column + ! + ! See `set_height` for a description of the CAM <-> TUV-x vertical grid + ! mapping. + ! + ! Above layer densities are calculated using a scale height for air + ! and pre-calculated values for O2 and O3 + !----------------------------------------------------------------------- + subroutine set_radiator_profiles( this, i_col, ncol, fixed_species_conc, species_vmr, & + exo_column_conc, delta_pressure, cloud_fraction, liquid_water_content, & + optical_depth, single_scattering_albedo, asymmetry_factor ) + + use chem_mods, only : gas_pcnst, & ! number of non-fixed species + nfs, & ! number of fixed species + nabscol, & ! number of absorbing species (radiators) + indexm ! index for air density in fixed species array + use ppgrid, only : pcols ! maximum number of columns + + class(tuvx_ptr), intent(inout) :: this ! TUV-x calculator + integer, intent(in) :: i_col ! column to set conditions for + integer, intent(in) :: ncol ! number of columns + real(r8), intent(in) :: fixed_species_conc(ncol,pver,max(1,nfs)) ! fixed species densities + ! (molecule cm-3) + real(r8), intent(in) :: species_vmr(ncol,pver,max(1,gas_pcnst)) ! species volume mixing + ! ratios (mol mol-1) + real(r8), intent(in) :: exo_column_conc(ncol,0:pver,max(1,nabscol)) ! above column densities + ! (molecule cm-2) + real(r8), intent(in) :: delta_pressure(ncol,pver) ! pressure delta about midpoints (Pa) + real(r8), intent(in) :: cloud_fraction(ncol,pver) ! cloud fraction (unitless) + real(r8), intent(in) :: liquid_water_content(ncol,pver) ! liquid water content (kg/kg) + real(r8), intent(in) :: optical_depth(pcols, pver+1, this%n_wavelength_bins_) ! aerosol optical depth [unitless] + real(r8), intent(in) :: single_scattering_albedo(pcols, pver+1, this%n_wavelength_bins_) ! single scattering albedo [unitless] + real(r8), intent(in) :: asymmetry_factor(pcols, pver+1, this%n_wavelength_bins_) ! asymmetry factor [unitless] + + integer :: i_level + real(r8) :: tmp(pver) + real(r8) :: tau(pver+1, this%n_wavelength_bins_) + real(r8) :: edges(pver+2), densities(pver+1) + real(r8) :: exo_val + real(r8), parameter :: rgrav = 1.0_r8 / 9.80616_r8 ! reciprocal of acceleration by gravity (s/m) + real(r8), parameter :: km2cm = 1.0e5_r8 ! conversion from km to cm + + ! =========== + ! air profile + ! =========== + edges(1) = fixed_species_conc(i_col,pver,indexm) + edges(2:pver+1) = fixed_species_conc(i_col,pver:1:-1,indexm) + edges(pver+2) = fixed_species_conc(i_col,1,indexm) ! use upper mid-point value for top edge + densities(1:pver+1) = this%height_delta_(1:pver+1) * km2cm * & + sqrt(edges(1:pver+1)) * sqrt(edges(2:pver+2)) + call this%profiles_( PROFILE_INDEX_AIR )%update( & + edge_values = edges, layer_densities = densities, & + scale_height = 8.01_r8 ) ! scale height in [km] + + ! ========== + ! O2 profile + ! ========== + if( is_fixed_O2 ) then + edges(1) = fixed_species_conc(i_col,pver,index_O2) + edges(2:pver+1) = fixed_species_conc(i_col,pver:1:-1,index_O2) + edges(pver+2) = fixed_species_conc(i_col,1,index_O2) + else if( index_O2 > 0 ) then + edges(1) = species_vmr(i_col,pver,index_O2) * & + fixed_species_conc(i_col,pver,indexm) + edges(2:pver+1) = species_vmr(i_col,pver:1:-1,index_O2) * & + fixed_species_conc(i_col,pver:1:-1,indexm) + edges(pver+2) = species_vmr(i_col,1,index_O2) * & + fixed_species_conc(i_col,1,indexm) + else + edges(:) = 0.0_r8 + end if + densities(1:pver+1) = this%height_delta_(1:pver+1) * km2cm * & + sqrt(edges(1:pver+1)) * sqrt(edges(2:pver+2)) + call this%profiles_( PROFILE_INDEX_O2 )%update( & + edge_values = edges, layer_densities = densities, & + scale_height = 7.0_r8 ) + + ! ========== + ! O3 profile + ! ========== + if( is_fixed_O3 ) then + edges(1) = fixed_species_conc(i_col,pver,index_O3) + edges(2:pver+1) = fixed_species_conc(i_col,pver:1:-1,index_O3) + edges(pver+2) = fixed_species_conc(i_col,1,index_O3) + else if( index_O3 > 0 ) then + edges(1) = species_vmr(i_col,pver,index_O3) * & + fixed_species_conc(i_col,pver,indexm) + edges(2:pver+1) = species_vmr(i_col,pver:1:-1,index_O3) * & + fixed_species_conc(i_col,pver:1:-1,indexm) + edges(pver+2) = species_vmr(i_col,1,index_O3) * & + fixed_species_conc(i_col,1,indexm) + else + edges(:) = 0.0_r8 + end if + if( nabscol >= 1 ) then + densities(1) = 0.5_r8 * exo_column_conc(i_col,pver,1) + densities(2:pver) = 0.5_r8 * ( exo_column_conc(i_col,pver-1:1:-1,1) & + + exo_column_conc(i_col,pver:2:-1,1) ) + densities(pver+1) = exo_column_conc(i_col,0,1) & + + 0.5_r8 * exo_column_conc(i_col,1,1) + call this%profiles_( PROFILE_INDEX_O3 )%update( & + edge_values = edges, layer_densities = densities, & + exo_density = exo_column_conc(i_col,0,1) ) + else + densities(1:pver+1) = this%height_delta_(1:pver+1) * km2cm * & + ( edges(1:pver+1) + edges(2:pver+2) ) * 0.5_r8 + call this%profiles_( PROFILE_INDEX_O3 )%update( & + edge_values = edges, layer_densities = densities, & + scale_height = 7.0_r8 ) + end if + + ! =============== + ! aerosol profile + ! =============== + if( do_aerosol ) then + call this%radiators_( RADIATOR_INDEX_AEROSOL )%update( & + optical_depths = optical_depth(i_col,:,:), & + single_scattering_albedos = single_scattering_albedo(i_col,:,:), & + asymmetry_factors = asymmetry_factor(i_col,:,:) ) + end if + + ! ============= + ! cloud profile + ! ============= + if( do_clouds ) then + ! =================================================== + ! estimate cloud optical depth as: + ! liquid_water_path * 0.155 * cloud_fraction^(1.5) + ! =================================================== + associate( clouds => cloud_fraction(i_col,:) ) + where( clouds(:) /= 0.0_r8 ) + tmp(:) = ( rgrav * liquid_water_content(i_col,:) * delta_pressure(i_col,:) & + * 1.0e3_r8 / clouds(:) ) * 0.155_r8 * clouds(:)**1.5_r8 + elsewhere + tmp(:) = 0.0_r8 + end where + end associate + do i_level = 1, pver + tau(i_level,:) = tmp(pver-i_level+1) + end do + tau(pver+1,:) = 0.0_r8 + call this%radiators_( RADIATOR_INDEX_CLOUDS )%update( optical_depths = tau ) + end if + + end subroutine set_radiator_profiles + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Updates working arrays of aerosol optical properties for all + ! columns from the aerosol package + !----------------------------------------------------------------------- + subroutine get_aerosol_optical_properties( this, state, pbuf, optical_depth, & + single_scattering_albedo, asymmetry_factor ) + + use aer_rad_props, only : aer_rad_props_sw + use mo_util, only : rebin + use physics_types, only : physics_state + use physics_buffer, only : physics_buffer_desc + use ppgrid, only : pcols ! maximum number of columns + use radconstants, only : nswbands, & ! Number of CAM shortwave radiation bands + get_sw_spectral_boundaries + + class(tuvx_ptr), intent(inout) :: this ! TUV-x calculator + type(physics_state), target, intent(in) :: state + type(physics_buffer_desc), pointer, intent(inout) :: pbuf(:) + real(r8), intent(out) :: optical_depth(pcols,pver+1,this%n_wavelength_bins_) ! aerosol optical depth [unitless] + real(r8), intent(out) :: single_scattering_albedo(pcols,pver+1,this%n_wavelength_bins_) ! aerosol single scattering albedo [unitless] + real(r8), intent(out) :: asymmetry_factor(pcols,pver+1,this%n_wavelength_bins_) ! aerosol asymmetry factor [unitless] + + real(r8) :: wavelength_edges(nswbands+1) ! CAM radiation wavelength grid edges [nm] + real(r8) :: aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth + real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau + real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau + real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau + real(r8) :: low_bound(nswbands) ! lower bound of CAM wavenumber bins + real(r8) :: high_bound(nswbands) ! upper bound of CAM wavenumber bins + integer :: n_night ! number of night columns + integer :: idx_night(pcols) ! indices of night columns + integer :: n_tuvx_bins ! number of TUV-x wavelength bins + integer :: i_col, i_level ! column and level indices + + ! =================================================================== + ! return default optical properties if no aerosol module is available + ! =================================================================== + if( .not. do_aerosol ) then + optical_depth(:,:,:) = 0.0_r8 + single_scattering_albedo(:,:,:) = 0.0_r8 + asymmetry_factor(:,:,:) = 0.0_r8 + return + end if + + ! TODO just assume all daylight columns for now + ! can adjust later if necessary + n_night = 0 + idx_night(:) = 0 + + ! =========================================================== + ! get aerosol optical properties on native CAM radiation grid + ! =========================================================== + call aer_rad_props_sw( 0, state, pbuf, n_night, idx_night, & + aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f ) + + ! ========================================================================= + ! Convert CAM wavenumber grid to wavelength grid and re-order optics arrays + ! NOTE: CAM wavenumber grid is continuous and increasing, except that the + ! last bin should be moved to the just before the first bin (!?!) + ! ========================================================================= + call get_sw_spectral_boundaries( low_bound, high_bound, 'nm' ) + wavelength_edges(1:nswbands-1) = low_bound( nswbands-1:1:-1) + wavelength_edges(nswbands ) = low_bound( nswbands ) + wavelength_edges(nswbands+1 ) = high_bound(nswbands ) + call reorder_optics_array( aer_tau ) + call reorder_optics_array( aer_tau_w ) + call reorder_optics_array( aer_tau_w_g ) + + ! ============================================================= + ! regrid optical properties to TUV-x wavelength and height grid + ! ============================================================= + ! TODO is this the correct regridding scheme to use? + n_tuvx_bins = this%n_wavelength_bins_ + do i_col = 1, pcols + do i_level = 1, pver + 1 + call rebin(nswbands, n_tuvx_bins, wavelength_edges, this%wavelength_edges_, & + aer_tau(i_col,pver+1-i_level,:), optical_depth(i_col,i_level,:)) + call rebin(nswbands, n_tuvx_bins, wavelength_edges, this%wavelength_edges_, & + aer_tau_w(i_col,pver+1-i_level,:), single_scattering_albedo(i_col,i_level,:)) + call rebin(nswbands, n_tuvx_bins, wavelength_edges, this%wavelength_edges_, & + aer_tau_w_g(i_col,pver+1-i_level,:), asymmetry_factor(i_col,i_level,:)) + end do + end do + + ! ================================================================ + ! back-calculate the single scattering albedo and asymmetry factor + ! ================================================================ + associate( tau => optical_depth, & + omega => single_scattering_albedo, & + g => asymmetry_factor ) + where(omega > 0.0_r8 .and. g > 0.0_r8) + g = g / omega + elsewhere + g = 0.0_r8 + end where + where(tau > 0.0_r8 .and. omega > 0.0_r8) + omega = omega / tau + elsewhere + omega = 0.0_r8 + end where + end associate + + end subroutine get_aerosol_optical_properties + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Reorders elements of an optical property array for conversion + ! from wavenumber to wavelength grid and out-of-order final + ! element in CAM wavenumber grid + !----------------------------------------------------------------------- + subroutine reorder_optics_array( optics_array ) + + use ppgrid, only : pcols ! maximum number of columns + use radconstants, only : nswbands ! Number of CAM shortwave radiation bands + + real(r8), intent(inout) :: optics_array(pcols,0:pver,nswbands) ! optics array to reorder + + real(r8) :: working(pcols,0:pver,nswbands) ! working array + + working(:,:,:) = optics_array(:,:,:) + optics_array(:,:,1:nswbands-1) = working(:,:,nswbands-1:1:-1) + + end subroutine reorder_optics_array + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Calculates extreme-UV ionization rates + ! + ! NOTE This never includes an above-column layer + !----------------------------------------------------------------------- + subroutine calculate_euv_rates( solar_zenith_angle, fixed_species_conc, & + species_vmr, height_mid, height_int, euv_rates ) + + use chem_mods, only : gas_pcnst, & ! number of non-fixed species + nfs, & ! number of fixed species + indexm ! index for air density in fixed species array + use mo_jeuv, only : jeuv, neuv ! number of extreme-UV rates + use ref_pres, only : ptop_ref ! pressure at the top of the column (Pa) + + real(r8), intent(in) :: solar_zenith_angle ! degrees + real(r8), intent(in) :: fixed_species_conc(pver,max(1,nfs)) ! fixed species densities + ! (molecule cm-3) + real(r8), intent(in) :: species_vmr(pver,max(1,gas_pcnst)) ! species volume mixing + ! ratios (mol mol-1) + real(r8), intent(in) :: height_mid(pver) ! height at mid-points (km) + real(r8), intent(in) :: height_int(pver+1) ! height at interfaces (km) + real(r8), intent(out) :: euv_rates(pver,neuv) ! calculated extreme-UV rates + + real(r8) :: o_dens(pver), o2_dens(pver), n2_dens(pver), height_arg(pver) + + ! ========== + ! N2 density + ! ========== + if( is_fixed_N2 ) then + n2_dens(:) = fixed_species_conc(:pver,index_N2) + else + n2_dens(:) = species_vmr(:pver,index_N2) * fixed_species_conc(:pver,indexm) + end if + + ! ========= + ! O density + ! ========= + if( is_fixed_O ) then + o_dens(:) = fixed_species_conc(:pver,index_O) + else + o_dens(:) = species_vmr(:pver,index_O) * fixed_species_conc(:pver,indexm) + end if + + ! ========== + ! O2 density + ! ========== + if( is_fixed_O2 ) then + o2_dens(:) = fixed_species_conc(:pver,index_O2) + else + o2_dens(:) = species_vmr(:pver,index_O2) * fixed_species_conc(:pver,indexm) + end if + + ! ======================= + ! special height argument + ! ======================= + height_arg(:) = height_mid(:) + + call jeuv( pver, solar_zenith_angle, o_dens, o2_dens, n2_dens, height_arg, & + euv_rates(pver:1:-1,:) ) + + end subroutine calculate_euv_rates + +!================================================================================================ + + !----------------------------------------------------------------------- + ! Calculates NO photolysis rates + ! + ! NOTE: Always includes an above-column layer + !----------------------------------------------------------------------- + subroutine calculate_jno( solar_zenith_angle, et_flux, fixed_species_conc, species_vmr, & + height_int, jno ) + + use chem_mods, only : gas_pcnst, & ! number of non-fixed species + nfs, & ! number of fixed species + indexm ! index for air density in fixed species array + use mo_jshort, only : sphers, slant_col, calc_jno + use ref_pres, only : ptop_ref ! pressure at the top of the column (Pa) + + real(r8), intent(in) :: solar_zenith_angle ! degrees + real(r8), intent(in) :: et_flux(NUM_BINS_MS93) ! extraterrestrial flux MS93 grid + ! (photon cm-2 nm-1 s-1) + real(r8), intent(in) :: fixed_species_conc(pver,max(1,nfs)) ! fixed species densities + ! (molecule cm-3) + real(r8), intent(in) :: species_vmr(pver,max(1,gas_pcnst)) ! species volume mixing + ! ratios (mol mol-1) + real(r8), intent(in) :: height_int(pver+1) ! height at interfaces (km) + real(r8), intent(out) :: jno(pver) ! calculated NO rate + + ! species column densities (molecule cm-3) + real(kind=r8) :: n2_dens(pver+1), o2_dens(pver+1), o3_dens(pver+1), no_dens(pver+1) + ! species slant column densities (molecule cm-2) + real(kind=r8) :: o2_slant(pver+1), o3_slant(pver+1), no_slant(pver+1) + ! working photo rate array + real(kind=r8) :: work_jno(pver+1) + ! parameters needed to calculate slant column densities + ! (see sphers routine description for details) + integer :: nid(pver+1) + real(kind=r8) :: dsdh(0:pver+1,pver+1) + ! layer thickness (cm) + real(kind=r8) :: delz(pver+1) + ! conversion from km to cm + real(kind=r8), parameter :: km2cm = 1.0e5_r8 + + ! ========== + ! N2 density + ! ========== + if( is_fixed_N2 ) then + n2_dens(2:) = fixed_species_conc(:pver,index_N2) + else + n2_dens(2:) = species_vmr(:pver,index_N2) * fixed_species_conc(:pver,indexm) + end if + n2_dens(1) = n2_dens(2) * 0.9_r8 + + ! ========== + ! O2 density + ! ========== + if( is_fixed_O2 ) then + o2_dens(2:) = fixed_species_conc(:pver,index_O2) + else + o2_dens(2:) = species_vmr(:pver,index_O2) * fixed_species_conc(:pver,indexm) + end if + o2_dens(1) = o2_dens(2) * 7.0_r8 / ( height_int(1) - height_int(2) ) + + ! ========== + ! O3 density + ! ========== + if( is_fixed_O3 ) then + o3_dens(2:) = fixed_species_conc(:pver,index_O3) + else + o3_dens(2:) = species_vmr(:pver,index_O3) * fixed_species_conc(:pver,indexm) + end if + o3_dens(1) = o3_dens(2) * 7.0_r8 / ( height_int(1) - height_int(2) ) + + ! ========== + ! NO density + ! ========== + if( is_fixed_NO ) then + no_dens(2:) = fixed_species_conc(:pver,index_NO) + else + no_dens(2:) = species_vmr(:pver,index_NO) * fixed_species_conc(:pver,indexm) + end if + no_dens(1) = no_dens(2) * 0.9_r8 + + ! ================================ + ! calculate slant column densities + ! ================================ + call sphers( pver+1, height_int, solar_zenith_angle, dsdh, nid ) + delz(1:pver) = km2cm * ( height_int(1:pver) - height_int(2:pver+1) ) + call slant_col( pver+1, delz, dsdh, nid, o2_dens, o2_slant ) + call slant_col( pver+1, delz, dsdh, nid, o3_dens, o3_slant ) + call slant_col( pver+1, delz, dsdh, nid, no_dens, no_slant ) + + ! ========================================= + ! calculate the NO photolysis rate constant + ! ========================================= + call calc_jno( pver+1, et_flux, n2_dens, o2_slant, o3_slant, no_slant, work_jno ) + jno(:) = work_jno(pver:1:-1) + + end subroutine calculate_jno + +!================================================================================================ + +end module mo_tuvx diff --git a/src/chemistry/mozart/mo_waccm_hrates.F90 b/src/chemistry/mozart/mo_waccm_hrates.F90 index 4f368c749f..322360d94c 100644 --- a/src/chemistry/mozart/mo_waccm_hrates.F90 +++ b/src/chemistry/mozart/mo_waccm_hrates.F90 @@ -25,6 +25,11 @@ module mo_waccm_hrates logical :: has_hrates integer :: ele_temp_ndx, ion_temp_ndx + integer :: cpe_jo2a_ndx = -1 + integer :: cpe_jo2b_ndx = -1 + integer :: cpe_jo3a_ndx = -1 + integer :: cpe_jo3b_ndx = -1 + contains subroutine init_hrates( ) @@ -81,9 +86,34 @@ subroutine init_hrates( ) attr = 'total jo2 euv photolysis rate' call addfld( 'JO2_EUV', (/ 'lev' /), 'I', '/s', trim(attr) ) + attr = 'O2 + hv -> O1D + O3P tuvx photo-chem heating rate' + call addfld( 'QRS_O2A_tuvx', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + attr = 'O2 + hv -> O3P + O3P tuvx photo-chem heating rate' + call addfld( 'QRS_O2B_tuvx', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + attr = 'O3 + hv -> O1D + O2_1S tuvx photo-chem heating rate' + call addfld( 'QRS_O3A_tuvx', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + attr = 'O3 + hv -> O3P + O2 tuvx photo-chem heating rate' + call addfld( 'QRS_O3B_tuvx', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + + + attr = 'O2 + hv -> O1D + O3P table photo-chem heating rate' + call addfld( 'QRS_O2A_tabl', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + attr = 'O2 + hv -> O3P + O3P table photo-chem heating rate' + call addfld( 'QRS_O2B_tabl', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + attr = 'O3 + hv -> O1D + O2_1S table photo-chem heating rate' + call addfld( 'QRS_O3A_tabl', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + attr = 'O3 + hv -> O3P + O2 table photo-chem heating rate' + call addfld( 'QRS_O3B_tabl', (/ 'lev' /), 'I', 'K/s', trim(attr) ) + + ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index + cpe_jo2a_ndx = pbuf_get_index('CPE_jO2a',errcode=err) + cpe_jo2b_ndx = pbuf_get_index('CPE_jO2b',errcode=err) + cpe_jo3a_ndx = pbuf_get_index('CPE_jO3a',errcode=err) + cpe_jo3b_ndx = pbuf_get_index('CPE_jO3b',errcode=err) + end subroutine init_hrates subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) @@ -124,6 +154,9 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) use phys_control, only : waccmx_is use orbit, only : zenith + use physconst, only : avogad + use mo_tuvx, only : tuvx_active + !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- @@ -161,6 +194,14 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) real(r8) :: qrl_col(pver,4) ! column thermal heating > 200nm real(r8) :: qrs(ncol,pver,4) ! chunk thermal heating < 200nm real(r8) :: qrl(ncol,pver,4) ! chunk thermal heating > 200nm + + real(r8) :: qr_jo2a_tuvx(ncol,pver) ! heating rates (K/s) + real(r8) :: qr_jo2b_tuvx(ncol,pver) ! heating rates (K/s) + real(r8) :: qr_jo3a_tuvx(ncol,pver) ! heating rates (K/s) + real(r8) :: qr_jo3b_tuvx(ncol,pver) ! heating rates (K/s) + + real(r8) :: hfactor(pver) + real(r8) :: euv_hrate_col(pver) ! column euv thermal heating rate real(r8) :: co2_hrate_col(pver) ! column co2 nir heating rate real(r8) :: euv_hrate(ncol,pver) ! chunk euv thermal heating rate @@ -196,6 +237,11 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) real(r8), pointer :: ele_temp_fld(:,:) ! electron temperature pointer real(r8), pointer :: ion_temp_fld(:,:) ! ion temperature pointer + real(r8), pointer :: cpe_jo2a(:,:) ! chemical potential energy + real(r8), pointer :: cpe_jo2b(:,:) ! chemical potential energy + real(r8), pointer :: cpe_jo3a(:,:) ! chemical potential energy + real(r8), pointer :: cpe_jo3b(:,:) ! chemical potential energy + if ( ele_temp_ndx>0 .and. ion_temp_ndx>0 ) then call pbuf_get_field(pbuf, ele_temp_ndx, ele_temp_fld) call pbuf_get_field(pbuf, ion_temp_ndx, ion_temp_fld) @@ -204,6 +250,11 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) ion_temp_fld => state%t endif + if (cpe_jo2a_ndx > 0) call pbuf_get_field(pbuf, cpe_jo2a_ndx, cpe_jo2a) + if (cpe_jo2b_ndx > 0) call pbuf_get_field(pbuf, cpe_jo2b_ndx, cpe_jo2b) + if (cpe_jo3a_ndx > 0) call pbuf_get_field(pbuf, cpe_jo3a_ndx, cpe_jo3a) + if (cpe_jo3b_ndx > 0) call pbuf_get_field(pbuf, cpe_jo3b_ndx, cpe_jo3b) + qrs_tot(:ncol,:) = 0._r8 if (.not. has_hrates) return @@ -343,6 +394,12 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) qrl(:,k,m) = 0._r8 end do end do + + qr_jo2a_tuvx = 0._r8 + qr_jo2b_tuvx = 0._r8 + qr_jo3a_tuvx = 0._r8 + qr_jo3b_tuvx = 0._r8 + do k = 1,pver euv_hrate(:,k) = 0._r8 co2_hrate(:,k) = 0._r8 @@ -372,18 +429,29 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) mw(:) = mbar(i,:) cparg(:) = cpair(i,:) do_diag = .false. + call jshort( pver, sza, o2_line, o3_line, o2cc, & o3cc, tline, zarg, mw, qrs_col, & cparg, lchnk, i, co2cc, scco2, do_diag ) - call jlong( pver, sza, eff_alb, parg, tline, & - mw, o2_line, o3_line, colo3, qrl_col, & - cparg, kbot_hrates ) - do m = 1,4 - qrs(i,pver:1:-1,m) = qrs_col(:,m) * esfact - end do - do m = 2,4 - qrl(i,:,m) = qrl_col(:,m) * esfact - end do + + if (tuvx_active) then + hfactor(:) = avogad/(cparg(:)*mw(:)) + qr_jo2a_tuvx(i,:) = cpe_jo2a(i,:) * hfactor(:) * o2_line(:) * esfact + qr_jo2b_tuvx(i,:) = cpe_jo2b(i,:) * hfactor(:) * o2_line(:) * esfact + qr_jo3a_tuvx(i,:kbot_hrates) = cpe_jo3a(i,:kbot_hrates) * hfactor(:kbot_hrates) * o3_line(:kbot_hrates) * esfact + qr_jo3b_tuvx(i,:kbot_hrates) = cpe_jo3b(i,:kbot_hrates) * hfactor(:kbot_hrates) * o3_line(:kbot_hrates) * esfact + else + call jlong( pver, sza, eff_alb, parg, tline, & + mw, o2_line, o3_line, colo3, qrl_col, & + cparg, kbot_hrates ) + do m = 1,4 + qrs(i,pver:1:-1,m) = qrs_col(:,m) * esfact + end do + do m = 2,4 + qrl(i,:,m) = qrl_col(:,m) * esfact + end do + end if + call heuv( pver, sza, occ, o2cc, n2cc, & o_line, o2_line, n2_line, cparg, mw, & zarg, euv_hrate_col, kbot_hrates ) @@ -403,6 +471,8 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) write(iulog,*) '===================================' #endif co2_hrate(i,:kbot_hrates) = co2_hrate_col(:kbot_hrates) * esfact * daypsec + + end if end do column_loop @@ -414,10 +484,23 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) call outfld( 'QRS_LO2B', qrl(:,:,2), ncol, lchnk ) call outfld( 'QRS_LO3A', qrl(:,:,3), ncol, lchnk ) call outfld( 'QRS_LO3B', qrl(:,:,4), ncol, lchnk ) + call outfld( 'QRS_LO3', qrl(:,:,3)+qrl(:,:,4), ncol, lchnk ) call outfld( 'QRS_EUV', euv_hrate(:,:), ncol, lchnk ) call outfld( 'QRS_CO2NIR', co2_hrate(:,:), ncol, lchnk ) + if (tuvx_active) then + call outfld( 'QRS_O2A_tuvx', qr_jo2a_tuvx(:ncol,:), ncol, lchnk ) + call outfld( 'QRS_O2B_tuvx', qr_jo2b_tuvx(:ncol,:), ncol, lchnk ) + call outfld( 'QRS_O3A_tuvx', qr_jo3a_tuvx(:ncol,:), ncol, lchnk ) + call outfld( 'QRS_O3B_tuvx', qr_jo3b_tuvx(:ncol,:), ncol, lchnk ) + else + call outfld( 'QRS_O2A_tabl', qrs(:,:,1)+qrl(:,:,1), ncol, lchnk ) + call outfld( 'QRS_O2B_tabl', qrs(:,:,2)+qrl(:,:,2), ncol, lchnk ) + call outfld( 'QRS_O3A_tabl', qrs(:,:,3)+qrl(:,:,3), ncol, lchnk ) + call outfld( 'QRS_O3B_tabl', qrs(:,:,4)+qrl(:,:,4), ncol, lchnk ) + endif + !----------------------------------------------------------------------- ! ... chemical pot heating rate !----------------------------------------------------------------------- @@ -436,7 +519,7 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) call outfld( 'QRS_AUR', aur_hrate(:,:), ncol, lchnk ) !----------------------------------------------------------------------- -! ... airglow heating rate +! ... airglow heating rate -- for diagnostic output only within airglow !----------------------------------------------------------------------- call airglow( aghrate, vmr(1,1,id_o2_1s), vmr(1,1,id_o2_1d), vmr(1,1,id_o1d), reaction_rates, cpair, & ncol, lchnk ) @@ -444,11 +527,17 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) !----------------------------------------------------------------------- ! ... form total heating rate !----------------------------------------------------------------------- - do k = 1,kbot_hrates - qrs_tot(:ncol,k) = qrs(:,k,1) + qrs(:,k,2) + qrs(:,k,3) + qrs(:,k,4) & - + qrl(:,k,1) + qrl(:,k,2) + qrl(:,k,3) + qrl(:,k,4) - end do + if (tuvx_active) then + qrs_tot(:ncol,:kbot_hrates) = qr_jo2a_tuvx(:ncol,:kbot_hrates) + qr_jo2b_tuvx(:ncol,:kbot_hrates) & + + qr_jo3a_tuvx(:ncol,:kbot_hrates) + qr_jo3b_tuvx(:ncol,:kbot_hrates) + else + do k = 1,kbot_hrates + qrs_tot(:ncol,k) = qrs(:,k,1) + qrs(:,k,2) + qrs(:,k,3) + qrs(:,k,4) & + + qrl(:,k,1) + qrl(:,k,2) + qrl(:,k,3) + qrl(:,k,4) + end do + end if call outfld( 'QTHERMAL', qrs_tot, pcols, lchnk ) + do k = 1,kbot_hrates qrs_tot(:ncol,k) = qrs_tot(:ncol,k) & + cphrate(:,k) + euv_hrate(:,k) + aur_hrate(:,k) + co2_hrate(:,k)