From 1006c1e8df33e37dfb9dc404b258ed38e166ad93 Mon Sep 17 00:00:00 2001 From: Steve Goldhaber Date: Fri, 1 Sep 2023 13:29:46 +0200 Subject: [PATCH] Restore pp_trop_mam_oslo from cam_cesm2_1_rel_05-Nor_v1.0.5 Copy pp_trop_mam_oslo mechanism and pp source from cam_cesm2_1_rel_05-Nor_v1.0.5 Specify trop_mam_oslo for NF compsets Some cleanup in configure Fortran is identical to that created with the new chemistry preprocessor --- bld/config_files/definition.xml | 2 +- bld/configure | 10 +- cime_config/config_component.xml | 3 +- src/chemistry/pp_trop_mam_oslo/chem_mech.doc | 162 ++++++++ src/chemistry/pp_trop_mam_oslo/chem_mech.in | 109 +++++ src/chemistry/pp_trop_mam_oslo/chem_mods.F90 | 50 +++ src/chemistry/pp_trop_mam_oslo/m_rxt_id.F90 | 16 + src/chemistry/pp_trop_mam_oslo/m_spc_id.F90 | 33 ++ src/chemistry/pp_trop_mam_oslo/mo_adjrxt.F90 | 34 ++ src/chemistry/pp_trop_mam_oslo/mo_exp_sol.F90 | 79 ++++ src/chemistry/pp_trop_mam_oslo/mo_imp_sol.F90 | 392 ++++++++++++++++++ src/chemistry/pp_trop_mam_oslo/mo_indprd.F90 | 56 +++ .../pp_trop_mam_oslo/mo_lin_matrix.F90 | 74 ++++ .../pp_trop_mam_oslo/mo_lu_factor.F90 | 57 +++ .../pp_trop_mam_oslo/mo_lu_solve.F90 | 90 ++++ .../pp_trop_mam_oslo/mo_nln_matrix.F90 | 103 +++++ src/chemistry/pp_trop_mam_oslo/mo_phtadj.F90 | 24 ++ .../pp_trop_mam_oslo/mo_prod_loss.F90 | 97 +++++ .../pp_trop_mam_oslo/mo_rxt_rates_conv.F90 | 25 ++ src/chemistry/pp_trop_mam_oslo/mo_setrxt.F90 | 90 ++++ src/chemistry/pp_trop_mam_oslo/mo_sim_dat.F90 | 132 ++++++ 21 files changed, 1630 insertions(+), 8 deletions(-) create mode 100644 src/chemistry/pp_trop_mam_oslo/chem_mech.doc create mode 100644 src/chemistry/pp_trop_mam_oslo/chem_mech.in create mode 100644 src/chemistry/pp_trop_mam_oslo/chem_mods.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/m_rxt_id.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/m_spc_id.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_adjrxt.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_exp_sol.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_imp_sol.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_indprd.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_lin_matrix.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_lu_factor.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_lu_solve.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_nln_matrix.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_phtadj.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_prod_loss.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_rxt_rates_conv.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_setrxt.F90 create mode 100644 src/chemistry/pp_trop_mam_oslo/mo_sim_dat.F90 diff --git a/bld/config_files/definition.xml b/bld/config_files/definition.xml index 251669792a..d00eda8f14 100644 --- a/bld/config_files/definition.xml +++ b/bld/config_files/definition.xml @@ -224,7 +224,7 @@ User specified C compiler options to append to Makefile defaults. User specified Fortran compiler overrides Makefile default. - + Type of Fortran compiler. Used when -fc specifies a generic wrapper script such as mpif90 or ftn. diff --git a/bld/configure b/bld/configure index dce056db30..5df0dd23f2 100755 --- a/bld/configure +++ b/bld/configure @@ -154,13 +154,13 @@ OPTIONS -cosp_libdir Directory containing COSP library. -esmf_libdir Directory containing ESMF library and esmf.mk file. -fc User specified Fortran compiler. Overrides Makefile default. - -fc_type Type of Fortran compiler [pgi | intel | gnu | pathscale - | ibm | nag]. This argument is used in conjunction + -fc_type Type of Fortran compiler [intel | gnu]. + This argument is used in conjunction with the -fc argument when the name of the fortran - compiler refers to a wrapper script (e.g., mpif90 - or ftn). In this case the user needs to specify + compiler refers to a wrapper script (e.g., mpif90). + In this case the user needs to specify the type of Fortran compiler that is being invoked - by the wrapper script. Default: pgi + by the wrapper script. Default: gnu -fflags A string of user specified Fortran compiler flags. Appended to Makefile defaults. See -fopt to override optimization flags. -fopt A string of user specified Fortran compiler optimization flags. diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index 6326bd3c8a..02ca4efc54 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -145,8 +145,7 @@ -phys cam_dev -chem ghg_mam4 -chem trop_strat_mam5_vbs - -camnor - + -camnor -chem trop_mam_oslo -chem trop_mam7 -chem trop_strat_mam5_vbsext -chem trop_strat_mam5_ts2 diff --git a/src/chemistry/pp_trop_mam_oslo/chem_mech.doc b/src/chemistry/pp_trop_mam_oslo/chem_mech.doc new file mode 100644 index 0000000000..328b93a925 --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/chem_mech.doc @@ -0,0 +1,162 @@ + + + |--------------------------------------------------------------------------------------------------| + | | + | | + | This is a mam-oslo simulation with : | + | (1) comments to add | + | (2) comments to add | + | | + | | + |--------------------------------------------------------------------------------------------------| + + + Solution species + ( 1) SO2 + ( 2) H2SO4 + ( 3) DMS (CH3SCH3) + ( 4) H2O2 + ( 5) SO4_NA (H2SO4) + ( 6) SO4_A1 (H2SO4) + ( 7) SO4_A2 (NH4HSO4) + ( 8) SO4_AC (H2SO4) + ( 9) SO4_PR (H2SO4) + ( 10) BC_N (C) + ( 11) BC_AX (C) + ( 12) BC_NI (C) + ( 13) BC_A (C) + ( 14) BC_AI (C) + ( 15) BC_AC (C) + ( 16) OM_NI (C) + ( 17) OM_AI (C) + ( 18) OM_AC (C) + ( 19) DST_A2 (AlSiO5) + ( 20) DST_A3 (AlSiO5) + ( 21) SS_A1 (NaCl) + ( 22) SS_A2 (NaCl) + ( 23) SS_A3 (NaCl) + ( 24) SOA_NA (C10H16O2) + ( 25) SOA_A1 (C10H16O2) + ( 26) SOA_LV (C10H16O2) + ( 27) SOA_SV (C10H16O2) + ( 28) monoterp (C10H16) + ( 29) isoprene (C5H8) + ( 30) H2O + + + Invariant species + ( 1) M + ( 2) N2 + ( 3) O2 + ( 4) O3 + ( 5) OH + ( 6) NO3 + ( 7) HO2 + + + Column integrals + ( 1) O3 - 0.000E+00 + ( 2) O2 - 0.000E+00 + +Class List +========== + + Implicit + -------- + ( 1) DMS + ( 2) SO2 + ( 3) H2O2 + ( 4) SO4_NA + ( 5) SO4_A1 + ( 6) SO4_A2 + ( 7) SO4_AC + ( 8) SO4_PR + ( 9) BC_N + ( 10) BC_AX + ( 11) BC_NI + ( 12) BC_A + ( 13) BC_AI + ( 14) BC_AC + ( 15) OM_NI + ( 16) OM_AI + ( 17) OM_AC + ( 18) DST_A2 + ( 19) DST_A3 + ( 20) SS_A1 + ( 21) SS_A2 + ( 22) SS_A3 + ( 23) H2SO4 + ( 24) SOA_NA + ( 25) SOA_A1 + ( 26) SOA_LV + ( 27) SOA_SV + ( 28) monoterp + ( 29) isoprene + ( 30) H2O + + Photolysis + jh2o2 ( 1) H2O2 + hv -> (No products) rate = ** User defined ** ( 1) + + Reactions + usr_HO2_HO2 ( 1) HO2 + HO2 -> H2O2 rate = ** User defined ** ( 2) + ( 2) H2O2 + OH -> H2O + HO2 rate = 2.90E-12*exp( -160./t) ( 3) + ( 3) DMS + OH -> SO2 rate = 9.60E-12*exp( -234./t) ( 4) + ( 4) DMS + NO3 -> SO2 + {HNO3} rate = 1.90E-13*exp( 520./t) ( 5) + ( 5) SO2 + OH + M -> H2SO4 + M troe : ko=3.00E-31*(300/t)**3.30 ( 6) + ki=1.50E-12 + f=0.60 + usr_DMS_OH ( 6) DMS + OH -> .75*SO2 + .5*HO2 + 0.029*SOA_LV + 0.114*SOA_SV rate = ** User defined ** ( 7) + ( 7) monoterp + O3 -> .15*SOA_LV rate = 8.05E-16*exp( -640./t) ( 8) + ( 8) monoterp + OH -> .15*SOA_SV rate = 1.20E-11*exp( 440./t) ( 9) + ( 9) monoterp + NO3 -> .15*SOA_SV rate = 1.20E-12*exp( 490./t) ( 10) + ( 10) isoprene + O3 -> .05*SOA_SV rate = 1.03E-14*exp( -1995./t) ( 11) + ( 11) isoprene + OH -> .05*SOA_SV rate = 2.70E-11*exp( 390./t) ( 12) + ( 12) isoprene + NO3 -> .05*SOA_SV rate = 3.15E-12*exp( -450./t) ( 13) + +Heterogeneous loss species + +Extraneous prod/loss species + ( 1) SO2 (dataset) + ( 2) BC_NI (dataset) + ( 3) BC_AX (dataset) + ( 4) BC_N (dataset) + ( 5) OM_NI (dataset) + ( 6) SO4_PR (dataset) + ( 7) H2O (dataset) + + + Equation Report + + d(SO2)/dt = r3*OH*DMS + r4*NO3*DMS + .75*r6*OH*DMS + - r5*OH*M*SO2 + d(H2SO4)/dt = r5*OH*M*SO2 + d(DMS)/dt = - r3*OH*DMS - r4*NO3*DMS - r6*OH*DMS + d(H2O2)/dt = r1 + - j1*H2O2 - r2*OH*H2O2 + d(SO4_NA)/dt = 0 + d(SO4_A1)/dt = 0 + d(SO4_A2)/dt = 0 + d(SO4_AC)/dt = 0 + d(SO4_PR)/dt = 0 + d(BC_N)/dt = 0 + d(BC_AX)/dt = 0 + d(BC_NI)/dt = 0 + d(BC_A)/dt = 0 + d(BC_AI)/dt = 0 + d(BC_AC)/dt = 0 + d(OM_NI)/dt = 0 + d(OM_AI)/dt = 0 + d(OM_AC)/dt = 0 + d(DST_A2)/dt = 0 + d(DST_A3)/dt = 0 + d(SS_A1)/dt = 0 + d(SS_A2)/dt = 0 + d(SS_A3)/dt = 0 + d(SOA_NA)/dt = 0 + d(SOA_A1)/dt = 0 + d(SOA_LV)/dt = .029*r6*OH*DMS + .15*r7*O3*monoterp + d(SOA_SV)/dt = .114*r6*OH*DMS + .15*r8*OH*monoterp + .15*r9*NO3*monoterp + .05*r10*O3*isoprene + + .05*r11*OH*isoprene + .05*r12*NO3*isoprene + d(monoterp)/dt = - r7*O3*monoterp - r8*OH*monoterp - r9*NO3*monoterp + d(isoprene)/dt = - r10*O3*isoprene - r11*OH*isoprene - r12*NO3*isoprene + d(H2O)/dt = r2*OH*H2O2 diff --git a/src/chemistry/pp_trop_mam_oslo/chem_mech.in b/src/chemistry/pp_trop_mam_oslo/chem_mech.in new file mode 100644 index 0000000000..5897d33683 --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/chem_mech.in @@ -0,0 +1,109 @@ +BEGSIM + SPECIES + + Solution + SO2, H2SO4 + DMS -> CH3SCH3, H2O2 + SO4_NA->H2SO4, SO4_A1->H2SO4, SO4_A2->NH4HSO4 + SO4_AC->H2SO4, SO4_PR->H2SO4, BC_N->C + BC_AX->C, BC_NI->C, BC_A->C, BC_AI->C + BC_AC->C, OM_NI->C, OM_AI->C, OM_AC->C + DST_A2->AlSiO5, DST_A3->AlSiO5 + SS_A1->NaCl, SS_A2->NaCl, SS_A3->NaCl +* Approximate soa species with those of monoterpene oxidation products +* based on Paasonen et al. (2010); Taipale et al. (2008). + SOA_NA->C10H16O2, SOA_A1->C10H16O2 + SOA_LV ->C10H16O2, SOA_SV->C10H16O2 + monoterp -> C10H16, isoprene -> C5H8 + H2O + End Solution + + Fixed + M, N2, O2, O3, OH, NO3, HO2 + End Fixed + + Col-int + O3 = 0. + O2 = 0. + End Col-int + + End SPECIES + + Solution Classes + Explicit + End Explicit + Implicit + DMS, SO2, H2O2 + SO4_NA, SO4_A1, SO4_A2 + SO4_AC, SO4_PR, BC_N + BC_AX, BC_NI, BC_A, BC_AI + BC_AC, OM_NI, OM_AI, OM_AC + DST_A2, DST_A3 + SS_A1, SS_A2, SS_A3 , H2SO4 + SOA_NA, SOA_A1 + SOA_LV,SOA_SV, monoterp, isoprene + H2O + End Implicit + End Solution Classes + + CHEMISTRY + Photolysis + [jh2o2] H2O2 + hv -> + End Photolysis + + Reactions + [usr_HO2_HO2] HO2 + HO2 -> H2O2 + H2O2 + OH -> H2O + HO2 ; 2.9e-12, -160 + DMS + OH -> SO2 ; 9.6e-12, -234. + DMS + NO3 -> SO2 + HNO3 ; 1.9e-13, 520. + SO2 + OH + M -> H2SO4 + M ; 3.0e-31, 3.3, 1.5e-12, 0.0, 0.6 +* SOA has MW=168, and MSA=96, so to get correct MSA mass ==> factor of 96/168 = 0.57 +* Then account for 0.25 which is 0.25 MSA molec per DMS molec (the other 0.75 goes to SO2) +* Then 0.2 assumed yield for SOA_LV and 0.8 assumed yield for SOA_SV gives the coefficients below +* reaction rate from Chin et al 1996, JGR, vol 101, no D13 +* + [usr_DMS_OH] DMS + OH -> .75 * SO2 + .5 * HO2 + 0.029*SOA_LV + 0.114*SOA_SV +* +*cka: added organic vapor oxidation with constants from IUPAC below +* Assume a yield of 15% for SOA LV production from these reactions +* Assume a yield of 15 % for monoterpene and 5% for isoprene SOA SV production reactions +* SOA_LV: very low volatility, can nucleate or grow small particles (oxidation products from O3+monoterp) +* SOA_SV: rest of SOA formed + monoterp + O3 -> .15*SOA_LV ; 8.05e-16, -640. + monoterp + OH -> .15*SOA_SV ; 1.2e-11, 440. + monoterp + NO3 -> .15*SOA_SV ; 1.2e-12, 490. + isoprene + O3 -> .05*SOA_SV ; 1.03e-14, -1995. + isoprene + OH -> .05*SOA_SV ; 2.7e-11, 390. + isoprene + NO3 -> .05*SOA_SV ; 3.15e-12, -450. + End Reactions + + Heterogeneous + H2O2, SO2 + End Heterogeneous + + Ext Forcing + SO2 <- dataset + BC_NI <-dataset + BC_AX <-dataset + BC_N <-dataset + OM_NI <-dataset + SO4_PR <-dataset + H2O <- dataset + End Ext Forcing + + END CHEMISTRY + + SIMULATION PARAMETERS + + Version Options + model = cam + machine = intel + architecture = hybrid + vec_ftns = on + multitask = on + namemod = on + modules = on + End Version Options + + END SIMULATION PARAMETERS +ENDSIM diff --git a/src/chemistry/pp_trop_mam_oslo/chem_mods.F90 b/src/chemistry/pp_trop_mam_oslo/chem_mods.F90 new file mode 100644 index 0000000000..bcd985192b --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/chem_mods.F90 @@ -0,0 +1,50 @@ + module chem_mods +!-------------------------------------------------------------- +! ... Basic chemistry parameters and arrays +!-------------------------------------------------------------- + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none + save + integer, parameter :: phtcnt = 1, & ! number of photolysis reactions + rxntot = 13, & ! number of total reactions + gascnt = 12, & ! number of gas phase reactions + nabscol = 2, & ! number of absorbing column densities + gas_pcnst = 30, & ! number of "gas phase" species + nfs = 7, & ! number of "fixed" species + relcnt = 0, & ! number of relationship species + grpcnt = 0, & ! number of group members + nzcnt = 38, & ! number of non-zero matrix entries + extcnt = 7, & ! number of species with external forcing + clscnt1 = 0, & ! number of species in explicit class + clscnt2 = 0, & ! number of species in hov class + clscnt3 = 0, & ! number of species in ebi class + clscnt4 = 30, & ! number of species in implicit class + clscnt5 = 0, & ! number of species in rodas class + indexm = 1, & ! index of total atm density in invariant array + indexh2o = 0, & ! index of water vapor density + clsze = 1, & ! loop length for implicit chemistry + rxt_tag_cnt = 3, & + enthalpy_cnt = 0, & + nslvd = 0 + integer :: clscnt(5) = 0 + integer :: cls_rxt_cnt(4,5) = 0 + integer :: clsmap(gas_pcnst,5) = 0 + integer :: permute(gas_pcnst,5) = 0 + integer :: diag_map(clscnt4) = 0 + real(r8) :: adv_mass(gas_pcnst) = 0._r8 + real(r8) :: crb_mass(gas_pcnst) = 0._r8 + real(r8) :: fix_mass(max(1,nfs)) + real(r8), allocatable :: cph_enthalpy(:) + integer, allocatable :: cph_rid(:) + integer, allocatable :: num_rnts(:) + integer, allocatable :: rxt_tag_map(:) + real(r8), allocatable :: pht_alias_mult(:,:) + character(len=32), allocatable :: rxt_tag_lst(:) + character(len=16), allocatable :: pht_alias_lst(:,:) + character(len=16) :: inv_lst(max(1,nfs)) + character(len=16) :: extfrc_lst(max(1,extcnt)) + logical :: frc_from_dataset(max(1,extcnt)) + logical :: is_vector + logical :: is_scalar + character(len=16) :: slvd_lst(max(1,nslvd)) + end module chem_mods diff --git a/src/chemistry/pp_trop_mam_oslo/m_rxt_id.F90 b/src/chemistry/pp_trop_mam_oslo/m_rxt_id.F90 new file mode 100644 index 0000000000..278910603d --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/m_rxt_id.F90 @@ -0,0 +1,16 @@ + module m_rxt_id + implicit none + integer, parameter :: rid_jh2o2 = 1 + integer, parameter :: rid_usr_HO2_HO2 = 2 + integer, parameter :: rid_usr_DMS_OH = 7 + integer, parameter :: rid_r0003 = 3 + integer, parameter :: rid_r0004 = 4 + integer, parameter :: rid_r0005 = 5 + integer, parameter :: rid_r0006 = 6 + integer, parameter :: rid_r0008 = 8 + integer, parameter :: rid_r0009 = 9 + integer, parameter :: rid_r0010 = 10 + integer, parameter :: rid_r0011 = 11 + integer, parameter :: rid_r0012 = 12 + integer, parameter :: rid_r0013 = 13 + end module m_rxt_id diff --git a/src/chemistry/pp_trop_mam_oslo/m_spc_id.F90 b/src/chemistry/pp_trop_mam_oslo/m_spc_id.F90 new file mode 100644 index 0000000000..f288d118fa --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/m_spc_id.F90 @@ -0,0 +1,33 @@ + module m_spc_id + implicit none + integer, parameter :: id_SO2 = 1 + integer, parameter :: id_H2SO4 = 2 + integer, parameter :: id_DMS = 3 + integer, parameter :: id_H2O2 = 4 + integer, parameter :: id_SO4_NA = 5 + integer, parameter :: id_SO4_A1 = 6 + integer, parameter :: id_SO4_A2 = 7 + integer, parameter :: id_SO4_AC = 8 + integer, parameter :: id_SO4_PR = 9 + integer, parameter :: id_BC_N = 10 + integer, parameter :: id_BC_AX = 11 + integer, parameter :: id_BC_NI = 12 + integer, parameter :: id_BC_A = 13 + integer, parameter :: id_BC_AI = 14 + integer, parameter :: id_BC_AC = 15 + integer, parameter :: id_OM_NI = 16 + integer, parameter :: id_OM_AI = 17 + integer, parameter :: id_OM_AC = 18 + integer, parameter :: id_DST_A2 = 19 + integer, parameter :: id_DST_A3 = 20 + integer, parameter :: id_SS_A1 = 21 + integer, parameter :: id_SS_A2 = 22 + integer, parameter :: id_SS_A3 = 23 + integer, parameter :: id_SOA_NA = 24 + integer, parameter :: id_SOA_A1 = 25 + integer, parameter :: id_SOA_LV = 26 + integer, parameter :: id_SOA_SV = 27 + integer, parameter :: id_monoterp = 28 + integer, parameter :: id_isoprene = 29 + integer, parameter :: id_H2O = 30 + end module m_spc_id diff --git a/src/chemistry/pp_trop_mam_oslo/mo_adjrxt.F90 b/src/chemistry/pp_trop_mam_oslo/mo_adjrxt.F90 new file mode 100644 index 0000000000..fe6931f11d --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_adjrxt.F90 @@ -0,0 +1,34 @@ + module mo_adjrxt + private + public :: adjrxt + contains + subroutine adjrxt( rate, inv, m, ncol, nlev ) + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : nfs, rxntot + implicit none +!-------------------------------------------------------------------- +! ... dummy arguments +!-------------------------------------------------------------------- + integer, intent(in) :: ncol, nlev + real(r8), intent(in) :: inv(ncol,nlev,nfs) + real(r8), intent(in) :: m(ncol,nlev) + real(r8), intent(inout) :: rate(ncol,nlev,rxntot) +!-------------------------------------------------------------------- +! ... local variables +!-------------------------------------------------------------------- + real(r8) :: im(ncol,nlev) + im(:,:) = 1._r8 / m(:,:) + rate(:,:, 3) = rate(:,:, 3) * inv(:,:, 5) + rate(:,:, 4) = rate(:,:, 4) * inv(:,:, 5) + rate(:,:, 5) = rate(:,:, 5) * inv(:,:, 6) + rate(:,:, 7) = rate(:,:, 7) * inv(:,:, 5) + rate(:,:, 8) = rate(:,:, 8) * inv(:,:, 4) + rate(:,:, 9) = rate(:,:, 9) * inv(:,:, 5) + rate(:,:, 10) = rate(:,:, 10) * inv(:,:, 6) + rate(:,:, 11) = rate(:,:, 11) * inv(:,:, 4) + rate(:,:, 12) = rate(:,:, 12) * inv(:,:, 5) + rate(:,:, 13) = rate(:,:, 13) * inv(:,:, 6) + rate(:,:, 2) = rate(:,:, 2) * inv(:,:, 7) * inv(:,:, 7) * im(:,:) + rate(:,:, 6) = rate(:,:, 6) * inv(:,:, 5) * inv(:,:, 1) + end subroutine adjrxt + end module mo_adjrxt diff --git a/src/chemistry/pp_trop_mam_oslo/mo_exp_sol.F90 b/src/chemistry/pp_trop_mam_oslo/mo_exp_sol.F90 new file mode 100644 index 0000000000..cfde22391a --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_exp_sol.F90 @@ -0,0 +1,79 @@ +module mo_exp_sol + private + public :: exp_sol + public :: exp_sol_inti +contains + subroutine exp_sol_inti + use mo_tracname, only : solsym + use chem_mods, only : clscnt1, clsmap + use ppgrid, only : pver + use cam_history, only : addfld + implicit none + integer :: i,j + do i = 1,clscnt1 + j = clsmap(i,1) + call addfld( trim(solsym(j))//'_CHMP', (/ 'lev' /), 'I', '/cm3/s', 'chemical production rate' ) + call addfld( trim(solsym(j))//'_CHML', (/ 'lev' /), 'I', '/cm3/s', 'chemical loss rate' ) + enddo + end subroutine exp_sol_inti + subroutine exp_sol( base_sol, reaction_rates, het_rates, extfrc, delt, xhnm, ncol, lchnk, ltrop ) + !----------------------------------------------------------------------- + ! ... Exp_sol advances the volumetric mixing ratio + ! forward one time step via the fully explicit + ! Euler scheme + !----------------------------------------------------------------------- + use chem_mods, only : clscnt1, extcnt, gas_pcnst, clsmap, rxntot + use ppgrid, only : pcols, pver + use mo_prod_loss, only : exp_prod_loss + use mo_indprd, only : indprd + use shr_kind_mod, only : r8 => shr_kind_r8 + use cam_history, only : outfld + use mo_tracname, only : solsym + implicit none + !----------------------------------------------------------------------- + ! ... Dummy arguments + !----------------------------------------------------------------------- + integer, intent(in) :: ncol ! columns in chunck + integer, intent(in) :: lchnk ! chunk id + real(r8), intent(in) :: delt ! time step (s) + real(r8), intent(in) :: het_rates(ncol,pver,max(1,gas_pcnst)) ! het rates (1/cm^3/s) + real(r8), intent(in) :: reaction_rates(ncol,pver,rxntot) ! rxt rates (1/cm^3/s) + real(r8), intent(in) :: extfrc(ncol,pver,extcnt) ! "external insitu forcing" (1/cm^3/s) + real(r8), intent(in) :: xhnm(ncol,pver) + integer, intent(in) :: ltrop(pcols) ! chemistry troposphere boundary (index) + real(r8), intent(inout) :: base_sol(ncol,pver,gas_pcnst) ! working mixing ratios (vmr) + !----------------------------------------------------------------------- + ! ... Local variables + !----------------------------------------------------------------------- + integer :: i, k, l, m + real(r8), dimension(ncol,pver,clscnt1) :: & + prod, & + loss, & + ind_prd + real(r8), dimension(ncol,pver) :: wrk + !----------------------------------------------------------------------- + ! ... Put "independent" production in the forcing + !----------------------------------------------------------------------- + call indprd( 1, ind_prd, clscnt1, base_sol, extfrc, & + reaction_rates, ncol ) + !----------------------------------------------------------------------- + ! ... Form F(y) + !----------------------------------------------------------------------- + call exp_prod_loss( prod, loss, base_sol, reaction_rates, het_rates ) + !----------------------------------------------------------------------- + ! ... Solve for the mixing ratio at t(n+1) + !----------------------------------------------------------------------- + do m = 1,clscnt1 + l = clsmap(m,1) + do i = 1,ncol + do k = ltrop(i)+1,pver + base_sol(i,k,l) = base_sol(i,k,l) + delt * (prod(i,k,m) + ind_prd(i,k,m) - loss(i,k,m)) + end do + end do + wrk(:,:) = (prod(:,:,m) + ind_prd(:,:,m))*xhnm + call outfld( trim(solsym(l))//'_CHMP', wrk(:,:), ncol, lchnk ) + wrk(:,:) = (loss(:,:,m))*xhnm + call outfld( trim(solsym(l))//'_CHML', wrk(:,:), ncol, lchnk ) + end do + end subroutine exp_sol +end module mo_exp_sol diff --git a/src/chemistry/pp_trop_mam_oslo/mo_imp_sol.F90 b/src/chemistry/pp_trop_mam_oslo/mo_imp_sol.F90 new file mode 100644 index 0000000000..d885728ba4 --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_imp_sol.F90 @@ -0,0 +1,392 @@ +module mo_imp_sol + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : clscnt4, gas_pcnst, clsmap + use cam_logfile, only : iulog + implicit none + private + public :: imp_slv_inti, imp_sol + save + real(r8), parameter :: rel_err = 1.e-3_r8 + real(r8), parameter :: high_rel_err = 1.e-4_r8 + !----------------------------------------------------------------------- + ! Newton-Raphson iteration limits + !----------------------------------------------------------------------- + integer, parameter :: itermax = 11 + integer, parameter :: cut_limit = 5 + real(r8), parameter :: small = 1.e-40_r8 + real(r8) :: epsilon(clscnt4) + logical :: factor(itermax) +contains + subroutine imp_slv_inti + !----------------------------------------------------------------------- + ! ... Initialize the implict solver + !----------------------------------------------------------------------- + use mo_chem_utls, only : get_spc_ndx + implicit none + !----------------------------------------------------------------------- + ! ... Local variables + !----------------------------------------------------------------------- + integer :: m, ox_ndx, o3a_ndx + real(r8) :: eps(gas_pcnst) + factor(:) = .true. + eps(:) = rel_err + ox_ndx = get_spc_ndx( 'OX' ) + if( ox_ndx < 1 ) then + ox_ndx = get_spc_ndx( 'O3' ) + end if + if( ox_ndx > 0 ) then + eps(ox_ndx) = high_rel_err + end if + m = get_spc_ndx( 'NO' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'NO2' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'NO3' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'HNO3' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'HO2NO2' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'N2O5' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'OH' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'HO2' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + o3a_ndx = get_spc_ndx( 'O3A' ) + if( o3a_ndx > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'XNO' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'XNO2' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'XNO3' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'XHNO3' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'XHO2NO2' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'XNO2NO3' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + m = get_spc_ndx( 'NO2XNO3' ) + if( m > 0 ) then + eps(m) = high_rel_err + end if + do m = 1,clscnt4 + epsilon(m) = eps(clsmap(m,4)) + end do + end subroutine imp_slv_inti + subroutine imp_sol( base_sol, reaction_rates, het_rates, extfrc, delt, & + ncol,nlev, lchnk, prod_out, loss_out ) + !----------------------------------------------------------------------- + ! ... imp_sol advances the volumetric mixing ratio + ! forward one time step via the fully implicit euler scheme. + ! this source is meant for small l1 cache machines such as + ! the intel pentium and itanium cpus + !----------------------------------------------------------------------- + use chem_mods, only : rxntot, extcnt, nzcnt, permute, cls_rxt_cnt + use mo_tracname, only : solsym + use mo_lin_matrix, only : linmat + use mo_nln_matrix, only : nlnmat + use mo_lu_factor, only : lu_fac + use mo_lu_solve, only : lu_slv + use mo_prod_loss, only : imp_prod_loss + use mo_indprd, only : indprd + use time_manager, only : get_nstep + use perf_mod, only : t_startf, t_stopf + implicit none + !----------------------------------------------------------------------- + ! ... dummy args + !----------------------------------------------------------------------- + integer, intent(in) :: ncol ! columns in chunck + integer, intent(in) :: nlev + integer, intent(in) :: lchnk ! chunk id + real(r8), intent(in) :: delt ! time step (s) + real(r8), intent(in) :: reaction_rates(ncol,nlev,max(1,rxntot)) ! rxt rates (1/cm^3/s) + real(r8), intent(in) :: extfrc(ncol,nlev,max(1,extcnt)) ! external in-situ forcing (1/cm^3/s) + real(r8), intent(in) :: het_rates(ncol,nlev,max(1,gas_pcnst)) ! washout rates (1/s) + real(r8), intent(inout) :: base_sol(ncol,nlev,gas_pcnst) ! species mixing ratios (vmr) + real(r8), intent(out) :: prod_out(ncol,nlev,max(1,clscnt4)) + real(r8), intent(out) :: loss_out(ncol,nlev,max(1,clscnt4)) + !----------------------------------------------------------------------- + ! ... local variables + !----------------------------------------------------------------------- + integer :: nr_iter, & + lev, & + i, & + j, & + k, l, & + m + integer :: fail_cnt, cut_cnt, stp_con_cnt + integer :: nstep + real(r8) :: interval_done, dt, dti + real(r8) :: max_delta(max(1,clscnt4)) + real(r8) :: sys_jac(max(1,nzcnt)) + real(r8) :: lin_jac(max(1,nzcnt)) + real(r8), dimension(max(1,clscnt4)) :: & + solution, & + forcing, & + iter_invariant, & + prod, & + loss + real(r8) :: lrxt(max(1,rxntot)) + real(r8) :: lsol(max(1,gas_pcnst)) + real(r8) :: lhet(max(1,gas_pcnst)) + real(r8), dimension(ncol,nlev,max(1,clscnt4)) :: & + ind_prd + logical :: convergence + logical :: frc_mask, iter_conv + logical :: converged(max(1,clscnt4)) + solution(:) = 0._r8 + !----------------------------------------------------------------------- + ! ... class independent forcing + !----------------------------------------------------------------------- + if( cls_rxt_cnt(1,4) > 0 .or. extcnt > 0 ) then + call indprd( 4, ind_prd, clscnt4, base_sol, extfrc, & + reaction_rates, ncol ) + else + do m = 1,max(1,clscnt4) + ind_prd(:,:,m) = 0._r8 + end do + end if + level_loop : do lev = 1,nlev + column_loop : do i = 1,ncol + !----------------------------------------------------------------------- + ! ... transfer from base to local work arrays + !----------------------------------------------------------------------- + do m = 1,rxntot + lrxt(m) = reaction_rates(i,lev,m) + end do + if( gas_pcnst > 0 ) then + do m = 1,gas_pcnst + lhet(m) = het_rates(i,lev,m) + end do + end if + !----------------------------------------------------------------------- + ! ... time step loop + !----------------------------------------------------------------------- + dt = delt + cut_cnt = 0 + fail_cnt = 0 + stp_con_cnt = 0 + interval_done = 0._r8 + time_step_loop : do + dti = 1._r8 / dt + !----------------------------------------------------------------------- + ! ... transfer from base to local work arrays + !----------------------------------------------------------------------- + do m = 1,gas_pcnst + lsol(m) = base_sol(i,lev,m) + end do + !----------------------------------------------------------------------- + ! ... transfer from base to class array + !----------------------------------------------------------------------- + do k = 1,clscnt4 + j = clsmap(k,4) + m = permute(k,4) + solution(m) = lsol(j) + end do + !----------------------------------------------------------------------- + ! ... set the iteration invariant part of the function f(y) + !----------------------------------------------------------------------- + if( cls_rxt_cnt(1,4) > 0 .or. extcnt > 0 ) then + do m = 1,clscnt4 + iter_invariant(m) = dti * solution(m) + ind_prd(i,lev,m) + end do + else + do m = 1,clscnt4 + iter_invariant(m) = dti * solution(m) + end do + end if + !----------------------------------------------------------------------- + ! ... the linear component + !----------------------------------------------------------------------- + if( cls_rxt_cnt(2,4) > 0 ) then + call t_startf( 'lin_mat' ) + call linmat( lin_jac, lsol, lrxt, lhet ) + call t_stopf( 'lin_mat' ) + end if + !======================================================================= + ! the newton-raphson iteration for f(y) = 0 + !======================================================================= + iter_loop : do nr_iter = 1,itermax + !----------------------------------------------------------------------- + ! ... the non-linear component + !----------------------------------------------------------------------- + if( factor(nr_iter) ) then + call t_startf( 'nln_mat' ) + call nlnmat( sys_jac, lsol, lrxt, lin_jac, dti ) + call t_stopf( 'nln_mat' ) + !----------------------------------------------------------------------- + ! ... factor the "system" matrix + !----------------------------------------------------------------------- + call t_startf( 'lu_fac' ) + call lu_fac( sys_jac ) + call t_stopf( 'lu_fac' ) + end if + !----------------------------------------------------------------------- + ! ... form f(y) + !----------------------------------------------------------------------- + call t_startf( 'prod_loss' ) + call imp_prod_loss( prod, loss, lsol, lrxt, lhet ) + call t_stopf( 'prod_loss' ) + do m = 1,clscnt4 + forcing(m) = solution(m)*dti - (iter_invariant(m) + prod(m) - loss(m)) + end do + !----------------------------------------------------------------------- + ! ... solve for the mixing ratio at t(n+1) + !----------------------------------------------------------------------- + call t_startf( 'lu_slv' ) + call lu_slv( sys_jac, forcing ) + call t_stopf( 'lu_slv' ) + do m = 1,clscnt4 + solution(m) = solution(m) + forcing(m) + end do + !----------------------------------------------------------------------- + ! ... convergence measures + !----------------------------------------------------------------------- + if( nr_iter > 1 ) then + do k = 1,clscnt4 + m = permute(k,4) + if( abs(solution(m)) > 1.e-20_r8 ) then + max_delta(k) = abs( forcing(m)/solution(m) ) + else + max_delta(k) = 0._r8 + end if + end do + end if + !----------------------------------------------------------------------- + ! ... limit iterate + !----------------------------------------------------------------------- + where( solution(:) < 0._r8 ) + solution(:) = 0._r8 + endwhere + !----------------------------------------------------------------------- + ! ... transfer latest solution back to work array + !----------------------------------------------------------------------- + do k = 1,clscnt4 + j = clsmap(k,4) + m = permute(k,4) + lsol(j) = solution(m) + end do + !----------------------------------------------------------------------- + ! ... check for convergence + !----------------------------------------------------------------------- + converged(:) = .true. + if( nr_iter > 1 ) then + do k = 1,clscnt4 + m = permute(k,4) + frc_mask = abs( forcing(m) ) > small + if( frc_mask ) then + converged(k) = abs(forcing(m)) <= epsilon(k)*abs(solution(m)) + else + converged(k) = .true. + end if + end do + convergence = all( converged(:) ) + if( convergence ) then + exit + end if + end if + end do iter_loop + !----------------------------------------------------------------------- + ! ... check for newton-raphson convergence + !----------------------------------------------------------------------- + if( .not. convergence ) then + !----------------------------------------------------------------------- + ! ... non-convergence + !----------------------------------------------------------------------- + fail_cnt = fail_cnt + 1 + nstep = get_nstep() + write(iulog,'('' imp_sol: Time step '',1p,e21.13,'' failed to converge @ (lchnk,lev,col,nstep) = '',4i6)') & + dt,lchnk,lev,i,nstep + stp_con_cnt = 0 + if( cut_cnt < cut_limit ) then + cut_cnt = cut_cnt + 1 + if( cut_cnt < cut_limit ) then + dt = .5_r8 * dt + else + dt = .1_r8 * dt + end if + cycle time_step_loop + else + write(iulog,'('' imp_sol: Failed to converge @ (lchnk,lev,col,nstep,dt,time) = '',4i6,1p,2e21.13)') & + lchnk,lev,i,nstep,dt,interval_done+dt + do m = 1,clscnt4 + if( .not. converged(m) ) then + write(iulog,'(1x,a8,1x,1pe10.3)') solsym(clsmap(m,4)), max_delta(m) + end if + end do + end if + end if + !----------------------------------------------------------------------- + ! ... check for interval done + !----------------------------------------------------------------------- + interval_done = interval_done + dt + if( abs( delt - interval_done ) <= .0001_r8 ) then + if( fail_cnt > 0 ) then + write(iulog,*) 'imp_sol : @ (lchnk,lev,col) = ',lchnk,lev,i,' failed ',fail_cnt,' times' + end if + exit time_step_loop + else + !----------------------------------------------------------------------- + ! ... transfer latest solution back to base array + !----------------------------------------------------------------------- + if( convergence ) then + stp_con_cnt = stp_con_cnt + 1 + end if + do m = 1,gas_pcnst + base_sol(i,lev,m) = lsol(m) + end do + if( stp_con_cnt >= 2 ) then + dt = 2._r8*dt + stp_con_cnt = 0 + end if + dt = min( dt,delt-interval_done ) + ! write(iulog,'('' imp_sol: New time step '',1p,e21.13)') dt + end if + end do time_step_loop + !----------------------------------------------------------------------- + ! ... Transfer latest solution back to base array + !----------------------------------------------------------------------- + cls_loop: do k = 1,clscnt4 + j = clsmap(k,4) + m = permute(k,4) + base_sol(i,lev,j) = solution(m) + ! output diagnostics + prod_out(i,lev,k) = prod(k) + ind_prd(i,lev,k) + loss_out(i,lev,k) = loss(k) + end do cls_loop + end do column_loop + end do level_loop + end subroutine imp_sol +end module mo_imp_sol diff --git a/src/chemistry/pp_trop_mam_oslo/mo_indprd.F90 b/src/chemistry/pp_trop_mam_oslo/mo_indprd.F90 new file mode 100644 index 0000000000..f9ec6830fb --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_indprd.F90 @@ -0,0 +1,56 @@ + module mo_indprd + use shr_kind_mod, only : r8 => shr_kind_r8 + private + public :: indprd + contains + subroutine indprd( class, prod, nprod, y, extfrc, rxt, ncol ) + use chem_mods, only : gas_pcnst, extcnt, rxntot + use ppgrid, only : pver + implicit none +!-------------------------------------------------------------------- +! ... dummy arguments +!-------------------------------------------------------------------- + integer, intent(in) :: class + integer, intent(in) :: ncol + integer, intent(in) :: nprod + real(r8), intent(in) :: y(ncol,pver,gas_pcnst) + real(r8), intent(in) :: rxt(ncol,pver,rxntot) + real(r8), intent(in) :: extfrc(ncol,pver,extcnt) + real(r8), intent(inout) :: prod(ncol,pver,nprod) +!-------------------------------------------------------------------- +! ... "independent" production for Implicit species +!-------------------------------------------------------------------- + if( class == 4 ) then + prod(:,:,1) = 0._r8 + prod(:,:,2) = + extfrc(:,:,1) + prod(:,:,3) =rxt(:,:,2) + prod(:,:,4) = 0._r8 + prod(:,:,5) = 0._r8 + prod(:,:,6) = 0._r8 + prod(:,:,7) = 0._r8 + prod(:,:,8) = + extfrc(:,:,6) + prod(:,:,9) = + extfrc(:,:,4) + prod(:,:,10) = + extfrc(:,:,3) + prod(:,:,11) = + extfrc(:,:,2) + prod(:,:,12) = 0._r8 + prod(:,:,13) = 0._r8 + prod(:,:,14) = 0._r8 + prod(:,:,15) = + extfrc(:,:,5) + prod(:,:,16) = 0._r8 + prod(:,:,17) = 0._r8 + prod(:,:,18) = 0._r8 + prod(:,:,19) = 0._r8 + prod(:,:,20) = 0._r8 + prod(:,:,21) = 0._r8 + prod(:,:,22) = 0._r8 + prod(:,:,23) = 0._r8 + prod(:,:,24) = 0._r8 + prod(:,:,25) = 0._r8 + prod(:,:,26) = 0._r8 + prod(:,:,27) = 0._r8 + prod(:,:,28) = 0._r8 + prod(:,:,29) = 0._r8 + prod(:,:,30) = + extfrc(:,:,7) + end if + end subroutine indprd + end module mo_indprd diff --git a/src/chemistry/pp_trop_mam_oslo/mo_lin_matrix.F90 b/src/chemistry/pp_trop_mam_oslo/mo_lin_matrix.F90 new file mode 100644 index 0000000000..e4c7687ebb --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_lin_matrix.F90 @@ -0,0 +1,74 @@ + module mo_lin_matrix + private + public :: linmat + contains + subroutine linmat01( mat, y, rxt, het_rates ) +!---------------------------------------------- +! ... linear matrix entries for implicit species +!---------------------------------------------- + use chem_mods, only : gas_pcnst, rxntot, nzcnt + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none +!---------------------------------------------- +! ... dummy arguments +!---------------------------------------------- + real(r8), intent(in) :: y(gas_pcnst) + real(r8), intent(in) :: rxt(rxntot) + real(r8), intent(in) :: het_rates(max(1,gas_pcnst)) + real(r8), intent(inout) :: mat(nzcnt) + mat(1) = -( rxt(4) + rxt(5) + rxt(7) + het_rates(3) ) + mat(5) = -( rxt(6) + het_rates(1) ) + mat(2) = rxt(4) + rxt(5) + .750_r8*rxt(7) + mat(7) = -( rxt(1) + rxt(3) + het_rates(4) ) + mat(9) = -( het_rates(5) ) + mat(10) = -( het_rates(6) ) + mat(11) = -( het_rates(7) ) + mat(12) = -( het_rates(8) ) + mat(13) = -( het_rates(9) ) + mat(14) = -( het_rates(10) ) + mat(15) = -( het_rates(11) ) + mat(16) = -( het_rates(12) ) + mat(17) = -( het_rates(13) ) + mat(18) = -( het_rates(14) ) + mat(19) = -( het_rates(15) ) + mat(20) = -( het_rates(16) ) + mat(21) = -( het_rates(17) ) + mat(22) = -( het_rates(18) ) + mat(23) = -( het_rates(19) ) + mat(24) = -( het_rates(20) ) + mat(25) = -( het_rates(21) ) + mat(26) = -( het_rates(22) ) + mat(27) = -( het_rates(23) ) + mat(28) = -( het_rates(2) ) + mat(6) = rxt(6) + mat(29) = -( het_rates(24) ) + mat(30) = -( het_rates(25) ) + mat(31) = -( het_rates(26) ) + mat(3) = .029_r8*rxt(7) + mat(33) = .150_r8*rxt(8) + mat(32) = -( het_rates(27) ) + mat(4) = .114_r8*rxt(7) + mat(34) = .150_r8*rxt(9) + .150_r8*rxt(10) + mat(36) = .050_r8*rxt(11) + .050_r8*rxt(12) + .050_r8*rxt(13) + mat(35) = -( rxt(8) + rxt(9) + rxt(10) + het_rates(28) ) + mat(37) = -( rxt(11) + rxt(12) + rxt(13) + het_rates(29) ) + mat(38) = -( het_rates(30) ) + mat(8) = rxt(3) + end subroutine linmat01 + subroutine linmat( mat, y, rxt, het_rates ) +!---------------------------------------------- +! ... linear matrix entries for implicit species +!---------------------------------------------- + use chem_mods, only : gas_pcnst, rxntot, nzcnt + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none +!---------------------------------------------- +! ... dummy arguments +!---------------------------------------------- + real(r8), intent(in) :: y(gas_pcnst) + real(r8), intent(in) :: rxt(rxntot) + real(r8), intent(in) :: het_rates(max(1,gas_pcnst)) + real(r8), intent(inout) :: mat(nzcnt) + call linmat01( mat, y, rxt, het_rates ) + end subroutine linmat + end module mo_lin_matrix diff --git a/src/chemistry/pp_trop_mam_oslo/mo_lu_factor.F90 b/src/chemistry/pp_trop_mam_oslo/mo_lu_factor.F90 new file mode 100644 index 0000000000..703de44018 --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_lu_factor.F90 @@ -0,0 +1,57 @@ + module mo_lu_factor + private + public :: lu_fac + contains + subroutine lu_fac01( lu ) + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none +!----------------------------------------------------------------------- +! ... dummy args +!----------------------------------------------------------------------- + real(r8), intent(inout) :: lu(:) + lu(1) = 1._r8 / lu(1) + lu(2) = lu(2) * lu(1) + lu(3) = lu(3) * lu(1) + lu(4) = lu(4) * lu(1) + lu(5) = 1._r8 / lu(5) + lu(6) = lu(6) * lu(5) + lu(7) = 1._r8 / lu(7) + lu(8) = lu(8) * lu(7) + lu(9) = 1._r8 / lu(9) + lu(10) = 1._r8 / lu(10) + lu(11) = 1._r8 / lu(11) + lu(12) = 1._r8 / lu(12) + lu(13) = 1._r8 / lu(13) + lu(14) = 1._r8 / lu(14) + lu(15) = 1._r8 / lu(15) + lu(16) = 1._r8 / lu(16) + lu(17) = 1._r8 / lu(17) + lu(18) = 1._r8 / lu(18) + lu(19) = 1._r8 / lu(19) + lu(20) = 1._r8 / lu(20) + lu(21) = 1._r8 / lu(21) + lu(22) = 1._r8 / lu(22) + lu(23) = 1._r8 / lu(23) + lu(24) = 1._r8 / lu(24) + lu(25) = 1._r8 / lu(25) + lu(26) = 1._r8 / lu(26) + lu(27) = 1._r8 / lu(27) + lu(28) = 1._r8 / lu(28) + lu(29) = 1._r8 / lu(29) + lu(30) = 1._r8 / lu(30) + lu(31) = 1._r8 / lu(31) + lu(32) = 1._r8 / lu(32) + lu(35) = 1._r8 / lu(35) + lu(37) = 1._r8 / lu(37) + lu(38) = 1._r8 / lu(38) + end subroutine lu_fac01 + subroutine lu_fac( lu ) + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none +!----------------------------------------------------------------------- +! ... dummy args +!----------------------------------------------------------------------- + real(r8), intent(inout) :: lu(:) + call lu_fac01( lu ) + end subroutine lu_fac + end module mo_lu_factor diff --git a/src/chemistry/pp_trop_mam_oslo/mo_lu_solve.F90 b/src/chemistry/pp_trop_mam_oslo/mo_lu_solve.F90 new file mode 100644 index 0000000000..862191c56b --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_lu_solve.F90 @@ -0,0 +1,90 @@ + module mo_lu_solve + private + public :: lu_slv + contains + subroutine lu_slv01( lu, b ) + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : clscnt4, nzcnt + implicit none +!----------------------------------------------------------------------- +! ... Dummy args +!----------------------------------------------------------------------- + real(r8), intent(in) :: lu(:) + real(r8), intent(inout) :: b(:) +!----------------------------------------------------------------------- +! ... Local variables +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +! ... solve L * y = b +!----------------------------------------------------------------------- + b(2) = b(2) - lu(2) * b(1) + b(26) = b(26) - lu(3) * b(1) + b(27) = b(27) - lu(4) * b(1) + b(23) = b(23) - lu(6) * b(2) + b(30) = b(30) - lu(8) * b(3) + end subroutine lu_slv01 + subroutine lu_slv02( lu, b ) + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : clscnt4, nzcnt + implicit none +!----------------------------------------------------------------------- +! ... Dummy args +!----------------------------------------------------------------------- + real(r8), intent(in) :: lu(:) + real(r8), intent(inout) :: b(:) +!----------------------------------------------------------------------- +! ... Local variables +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +! ... solve L * y = b +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +! ... Solve U * x = y +!----------------------------------------------------------------------- + b(30) = b(30) * lu(38) + b(29) = b(29) * lu(37) + b(27) = b(27) - lu(36) * b(29) + b(28) = b(28) * lu(35) + b(27) = b(27) - lu(34) * b(28) + b(26) = b(26) - lu(33) * b(28) + b(27) = b(27) * lu(32) + b(26) = b(26) * lu(31) + b(25) = b(25) * lu(30) + b(24) = b(24) * lu(29) + b(23) = b(23) * lu(28) + b(22) = b(22) * lu(27) + b(21) = b(21) * lu(26) + b(20) = b(20) * lu(25) + b(19) = b(19) * lu(24) + b(18) = b(18) * lu(23) + b(17) = b(17) * lu(22) + b(16) = b(16) * lu(21) + b(15) = b(15) * lu(20) + b(14) = b(14) * lu(19) + b(13) = b(13) * lu(18) + b(12) = b(12) * lu(17) + b(11) = b(11) * lu(16) + b(10) = b(10) * lu(15) + b(9) = b(9) * lu(14) + b(8) = b(8) * lu(13) + b(7) = b(7) * lu(12) + b(6) = b(6) * lu(11) + b(5) = b(5) * lu(10) + b(4) = b(4) * lu(9) + b(3) = b(3) * lu(7) + b(2) = b(2) * lu(5) + b(1) = b(1) * lu(1) + end subroutine lu_slv02 + subroutine lu_slv( lu, b ) + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : clscnt4, nzcnt + implicit none +!----------------------------------------------------------------------- +! ... Dummy args +!----------------------------------------------------------------------- + real(r8), intent(in) :: lu(:) + real(r8), intent(inout) :: b(:) + call lu_slv01( lu, b ) + call lu_slv02( lu, b ) + end subroutine lu_slv + end module mo_lu_solve diff --git a/src/chemistry/pp_trop_mam_oslo/mo_nln_matrix.F90 b/src/chemistry/pp_trop_mam_oslo/mo_nln_matrix.F90 new file mode 100644 index 0000000000..3aa3cd3972 --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_nln_matrix.F90 @@ -0,0 +1,103 @@ + module mo_nln_matrix + use shr_kind_mod, only : r8 => shr_kind_r8 + private + public :: nlnmat + contains + subroutine nlnmat( mat, y, rxt, lmat, dti ) + use chem_mods, only : gas_pcnst, rxntot, nzcnt + implicit none +!---------------------------------------------- +! ... dummy arguments +!---------------------------------------------- + real(r8), intent(in) :: dti + real(r8), intent(in) :: lmat(nzcnt) + real(r8), intent(in) :: y(gas_pcnst) + real(r8), intent(in) :: rxt(rxntot) + real(r8), intent(inout) :: mat(nzcnt) + call nlnmat_finit( mat, lmat, dti ) + end subroutine nlnmat + subroutine nlnmat_finit( mat, lmat, dti ) + use chem_mods, only : gas_pcnst, rxntot, nzcnt + implicit none +!---------------------------------------------- +! ... dummy arguments +!---------------------------------------------- + real(r8), intent(in) :: dti + real(r8), intent(in) :: lmat(nzcnt) + real(r8), intent(inout) :: mat(nzcnt) +!---------------------------------------------- +! ... local variables +!---------------------------------------------- +!---------------------------------------------- +! ... complete matrix entries implicit species +!---------------------------------------------- + mat( 1) = lmat( 1) + mat( 2) = lmat( 2) + mat( 3) = lmat( 3) + mat( 4) = lmat( 4) + mat( 5) = lmat( 5) + mat( 6) = lmat( 6) + mat( 7) = lmat( 7) + mat( 8) = lmat( 8) + mat( 9) = lmat( 9) + mat( 10) = lmat( 10) + mat( 11) = lmat( 11) + mat( 12) = lmat( 12) + mat( 13) = lmat( 13) + mat( 14) = lmat( 14) + mat( 15) = lmat( 15) + mat( 16) = lmat( 16) + mat( 17) = lmat( 17) + mat( 18) = lmat( 18) + mat( 19) = lmat( 19) + mat( 20) = lmat( 20) + mat( 21) = lmat( 21) + mat( 22) = lmat( 22) + mat( 23) = lmat( 23) + mat( 24) = lmat( 24) + mat( 25) = lmat( 25) + mat( 26) = lmat( 26) + mat( 27) = lmat( 27) + mat( 28) = lmat( 28) + mat( 29) = lmat( 29) + mat( 30) = lmat( 30) + mat( 31) = lmat( 31) + mat( 32) = lmat( 32) + mat( 33) = lmat( 33) + mat( 34) = lmat( 34) + mat( 35) = lmat( 35) + mat( 36) = lmat( 36) + mat( 37) = lmat( 37) + mat( 38) = lmat( 38) + mat( 1) = mat( 1) - dti + mat( 5) = mat( 5) - dti + mat( 7) = mat( 7) - dti + mat( 9) = mat( 9) - dti + mat( 10) = mat( 10) - dti + mat( 11) = mat( 11) - dti + mat( 12) = mat( 12) - dti + mat( 13) = mat( 13) - dti + mat( 14) = mat( 14) - dti + mat( 15) = mat( 15) - dti + mat( 16) = mat( 16) - dti + mat( 17) = mat( 17) - dti + mat( 18) = mat( 18) - dti + mat( 19) = mat( 19) - dti + mat( 20) = mat( 20) - dti + mat( 21) = mat( 21) - dti + mat( 22) = mat( 22) - dti + mat( 23) = mat( 23) - dti + mat( 24) = mat( 24) - dti + mat( 25) = mat( 25) - dti + mat( 26) = mat( 26) - dti + mat( 27) = mat( 27) - dti + mat( 28) = mat( 28) - dti + mat( 29) = mat( 29) - dti + mat( 30) = mat( 30) - dti + mat( 31) = mat( 31) - dti + mat( 32) = mat( 32) - dti + mat( 35) = mat( 35) - dti + mat( 37) = mat( 37) - dti + mat( 38) = mat( 38) - dti + end subroutine nlnmat_finit + end module mo_nln_matrix diff --git a/src/chemistry/pp_trop_mam_oslo/mo_phtadj.F90 b/src/chemistry/pp_trop_mam_oslo/mo_phtadj.F90 new file mode 100644 index 0000000000..aaa43829fe --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_phtadj.F90 @@ -0,0 +1,24 @@ + module mo_phtadj + private + public :: phtadj + contains + subroutine phtadj( p_rate, inv, m, ncol, nlev ) + use chem_mods, only : nfs, phtcnt + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none +!-------------------------------------------------------------------- +! ... dummy arguments +!-------------------------------------------------------------------- + integer, intent(in) :: ncol, nlev + real(r8), intent(in) :: inv(ncol,nlev,max(1,nfs)) + real(r8), intent(in) :: m(ncol,nlev) + real(r8), intent(inout) :: p_rate(ncol,nlev,max(1,phtcnt)) +!-------------------------------------------------------------------- +! ... local variables +!-------------------------------------------------------------------- + integer :: k + real(r8) :: im(ncol,nlev) + do k = 1,nlev + end do + end subroutine phtadj + end module mo_phtadj diff --git a/src/chemistry/pp_trop_mam_oslo/mo_prod_loss.F90 b/src/chemistry/pp_trop_mam_oslo/mo_prod_loss.F90 new file mode 100644 index 0000000000..0cbb77be48 --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_prod_loss.F90 @@ -0,0 +1,97 @@ + module mo_prod_loss + use shr_kind_mod, only : r8 => shr_kind_r8 + private + public :: exp_prod_loss + public :: imp_prod_loss + contains + subroutine exp_prod_loss( prod, loss, y, rxt, het_rates ) + use ppgrid, only : pver + implicit none +!-------------------------------------------------------------------- +! ... dummy args +!-------------------------------------------------------------------- + real(r8), dimension(:,:,:), intent(out) :: & + prod, & + loss + real(r8), intent(in) :: y(:,:,:) + real(r8), intent(in) :: rxt(:,:,:) + real(r8), intent(in) :: het_rates(:,:,:) + end subroutine exp_prod_loss + subroutine imp_prod_loss( prod, loss, y, rxt, het_rates ) + use ppgrid, only : pver + implicit none +!-------------------------------------------------------------------- +! ... dummy args +!-------------------------------------------------------------------- + real(r8), dimension(:), intent(out) :: & + prod, & + loss + real(r8), intent(in) :: y(:) + real(r8), intent(in) :: rxt(:) + real(r8), intent(in) :: het_rates(:) +!-------------------------------------------------------------------- +! ... loss and production for Implicit method +!-------------------------------------------------------------------- + loss(1) = ( + rxt(4) + rxt(5) + rxt(7) + het_rates(3))* y(3) + prod(1) = 0._r8 + loss(2) = ( + rxt(6) + het_rates(1))* y(1) + prod(2) = (rxt(4) +rxt(5) +.750_r8*rxt(7))*y(3) + loss(3) = ( + rxt(1) + rxt(3) + het_rates(4))* y(4) + prod(3) = 0._r8 + loss(4) = ( + het_rates(5))* y(5) + prod(4) = 0._r8 + loss(5) = ( + het_rates(6))* y(6) + prod(5) = 0._r8 + loss(6) = ( + het_rates(7))* y(7) + prod(6) = 0._r8 + loss(7) = ( + het_rates(8))* y(8) + prod(7) = 0._r8 + loss(8) = ( + het_rates(9))* y(9) + prod(8) = 0._r8 + loss(9) = ( + het_rates(10))* y(10) + prod(9) = 0._r8 + loss(10) = ( + het_rates(11))* y(11) + prod(10) = 0._r8 + loss(11) = ( + het_rates(12))* y(12) + prod(11) = 0._r8 + loss(12) = ( + het_rates(13))* y(13) + prod(12) = 0._r8 + loss(13) = ( + het_rates(14))* y(14) + prod(13) = 0._r8 + loss(14) = ( + het_rates(15))* y(15) + prod(14) = 0._r8 + loss(15) = ( + het_rates(16))* y(16) + prod(15) = 0._r8 + loss(16) = ( + het_rates(17))* y(17) + prod(16) = 0._r8 + loss(17) = ( + het_rates(18))* y(18) + prod(17) = 0._r8 + loss(18) = ( + het_rates(19))* y(19) + prod(18) = 0._r8 + loss(19) = ( + het_rates(20))* y(20) + prod(19) = 0._r8 + loss(20) = ( + het_rates(21))* y(21) + prod(20) = 0._r8 + loss(21) = ( + het_rates(22))* y(22) + prod(21) = 0._r8 + loss(22) = ( + het_rates(23))* y(23) + prod(22) = 0._r8 + loss(23) = ( + het_rates(2))* y(2) + prod(23) =rxt(6)*y(1) + loss(24) = ( + het_rates(24))* y(24) + prod(24) = 0._r8 + loss(25) = ( + het_rates(25))* y(25) + prod(25) = 0._r8 + loss(26) = ( + het_rates(26))* y(26) + prod(26) =.029_r8*rxt(7)*y(3) +.150_r8*rxt(8)*y(28) + loss(27) = ( + het_rates(27))* y(27) + prod(27) = (.050_r8*rxt(11) +.050_r8*rxt(12) +.050_r8*rxt(13))*y(29) & + + (.150_r8*rxt(9) +.150_r8*rxt(10))*y(28) +.114_r8*rxt(7)*y(3) + loss(28) = ( + rxt(8) + rxt(9) + rxt(10) + het_rates(28))* y(28) + prod(28) = 0._r8 + loss(29) = ( + rxt(11) + rxt(12) + rxt(13) + het_rates(29))* y(29) + prod(29) = 0._r8 + loss(30) = ( + het_rates(30))* y(30) + prod(30) =rxt(3)*y(4) + end subroutine imp_prod_loss + end module mo_prod_loss diff --git a/src/chemistry/pp_trop_mam_oslo/mo_rxt_rates_conv.F90 b/src/chemistry/pp_trop_mam_oslo/mo_rxt_rates_conv.F90 new file mode 100644 index 0000000000..7f9000e78e --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_rxt_rates_conv.F90 @@ -0,0 +1,25 @@ +module mo_rxt_rates_conv + use shr_kind_mod, only : r8 => shr_kind_r8 + implicit none + private + public :: set_rates +contains + subroutine set_rates( rxt_rates, sol, ncol ) + real(r8), intent(inout) :: rxt_rates(:,:,:) + real(r8), intent(in) :: sol(:,:,:) + integer, intent(in) :: ncol + rxt_rates(:ncol,:, 1) = rxt_rates(:ncol,:, 1)*sol(:ncol,:, 4) ! rate_const*H2O2 + ! rate_const + rxt_rates(:ncol,:, 3) = rxt_rates(:ncol,:, 3)*sol(:ncol,:, 4) ! rate_const*OH*H2O2 + rxt_rates(:ncol,:, 4) = rxt_rates(:ncol,:, 4)*sol(:ncol,:, 3) ! rate_const*OH*DMS + rxt_rates(:ncol,:, 5) = rxt_rates(:ncol,:, 5)*sol(:ncol,:, 3) ! rate_const*NO3*DMS + rxt_rates(:ncol,:, 6) = rxt_rates(:ncol,:, 6)*sol(:ncol,:, 1) ! rate_const*OH*M*SO2 + rxt_rates(:ncol,:, 7) = rxt_rates(:ncol,:, 7)*sol(:ncol,:, 3) ! rate_const*OH*DMS + rxt_rates(:ncol,:, 8) = rxt_rates(:ncol,:, 8)*sol(:ncol,:, 28) ! rate_const*O3*monoterp + rxt_rates(:ncol,:, 9) = rxt_rates(:ncol,:, 9)*sol(:ncol,:, 28) ! rate_const*OH*monoterp + rxt_rates(:ncol,:, 10) = rxt_rates(:ncol,:, 10)*sol(:ncol,:, 28) ! rate_const*NO3*monoterp + rxt_rates(:ncol,:, 11) = rxt_rates(:ncol,:, 11)*sol(:ncol,:, 29) ! rate_const*O3*isoprene + rxt_rates(:ncol,:, 12) = rxt_rates(:ncol,:, 12)*sol(:ncol,:, 29) ! rate_const*OH*isoprene + rxt_rates(:ncol,:, 13) = rxt_rates(:ncol,:, 13)*sol(:ncol,:, 29) ! rate_const*NO3*isoprene + end subroutine set_rates +end module mo_rxt_rates_conv diff --git a/src/chemistry/pp_trop_mam_oslo/mo_setrxt.F90 b/src/chemistry/pp_trop_mam_oslo/mo_setrxt.F90 new file mode 100644 index 0000000000..b9909e2e6b --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_setrxt.F90 @@ -0,0 +1,90 @@ + + module mo_setrxt + + use shr_kind_mod, only : r8 => shr_kind_r8 + + private + public :: setrxt + public :: setrxt_hrates + + contains + + subroutine setrxt( rate, temp, m, ncol ) + + use ppgrid, only : pver, pcols + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : rxntot + use mo_jpl, only : jpl + + implicit none + +!------------------------------------------------------- +! ... dummy arguments +!------------------------------------------------------- + integer, intent(in) :: ncol + real(r8), intent(in) :: temp(pcols,pver) + real(r8), intent(in) :: m(ncol,pver) + real(r8), intent(inout) :: rate(ncol,pver,rxntot) + +!------------------------------------------------------- +! ... local variables +!------------------------------------------------------- + integer :: n + real(r8) :: itemp(ncol,pver) + real(r8) :: exp_fac(ncol,pver) + real(r8) :: ko(ncol,pver) + real(r8) :: kinf(ncol,pver) + + itemp(:ncol,:) = 1._r8 / temp(:ncol,:) + n = ncol*pver + rate(:,:,3) = 2.9e-12_r8 * exp( -160._r8 * itemp(:,:) ) + rate(:,:,4) = 9.6e-12_r8 * exp( -234._r8 * itemp(:,:) ) + rate(:,:,5) = 1.9e-13_r8 * exp( 520._r8 * itemp(:,:) ) + rate(:,:,8) = 8.05e-16_r8 * exp( -640._r8 * itemp(:,:) ) + rate(:,:,9) = 1.2e-11_r8 * exp( 440._r8 * itemp(:,:) ) + rate(:,:,10) = 1.2e-12_r8 * exp( 490._r8 * itemp(:,:) ) + rate(:,:,11) = 1.03e-14_r8 * exp( -1995._r8 * itemp(:,:) ) + rate(:,:,12) = 2.7e-11_r8 * exp( 390._r8 * itemp(:,:) ) + rate(:,:,13) = 3.15e-12_r8 * exp( -450._r8 * itemp(:,:) ) + + itemp(:,:) = 300._r8 * itemp(:,:) + + ko(:,:) = 3.0e-31_r8 * itemp(:,:)**3.3_r8 + kinf(:,:) = 1.5e-12_r8 + call jpl( rate(1,1,6), m, 0.6_r8, ko, kinf, n ) + + end subroutine setrxt + + + subroutine setrxt_hrates( rate, temp, m, ncol, kbot ) + + use ppgrid, only : pver, pcols + use shr_kind_mod, only : r8 => shr_kind_r8 + use chem_mods, only : rxntot + use mo_jpl, only : jpl + + implicit none + +!------------------------------------------------------- +! ... dummy arguments +!------------------------------------------------------- + integer, intent(in) :: ncol + integer, intent(in) :: kbot + real(r8), intent(in) :: temp(pcols,pver) + real(r8), intent(in) :: m(ncol,pver) + real(r8), intent(inout) :: rate(ncol,pver,rxntot) + +!------------------------------------------------------- +! ... local variables +!------------------------------------------------------- + integer :: n + real(r8) :: itemp(ncol,kbot) + real(r8) :: exp_fac(ncol,kbot) + real(r8) :: ko(ncol,kbot) + real(r8) :: kinf(ncol,kbot) + real(r8) :: wrk(ncol,kbot) + + + end subroutine setrxt_hrates + + end module mo_setrxt diff --git a/src/chemistry/pp_trop_mam_oslo/mo_sim_dat.F90 b/src/chemistry/pp_trop_mam_oslo/mo_sim_dat.F90 new file mode 100644 index 0000000000..ff81ecfc5d --- /dev/null +++ b/src/chemistry/pp_trop_mam_oslo/mo_sim_dat.F90 @@ -0,0 +1,132 @@ + + module mo_sim_dat + + private + public :: set_sim_dat + + contains + + subroutine set_sim_dat + + use chem_mods, only : clscnt, cls_rxt_cnt, clsmap, permute, adv_mass, fix_mass, crb_mass + use chem_mods, only : diag_map + use chem_mods, only : phtcnt, rxt_tag_cnt, rxt_tag_lst, rxt_tag_map + use chem_mods, only : pht_alias_lst, pht_alias_mult + use chem_mods, only : extfrc_lst, inv_lst, slvd_lst + use chem_mods, only : enthalpy_cnt, cph_enthalpy, cph_rid, num_rnts, rxntot + use cam_abortutils,only : endrun + use mo_tracname, only : solsym + use chem_mods, only : frc_from_dataset + use chem_mods, only : is_scalar, is_vector + use shr_kind_mod, only : r8 => shr_kind_r8 + use cam_logfile, only : iulog + + implicit none + +!-------------------------------------------------------------- +! ... local variables +!-------------------------------------------------------------- + integer :: ios + + is_scalar = .true. + is_vector = .false. + + clscnt(:) = (/ 0, 0, 0, 30, 0 /) + + cls_rxt_cnt(:,4) = (/ 1, 12, 0, 30 /) + + solsym(: 30) = (/ 'SO2 ','H2SO4 ','DMS ','H2O2 ','SO4_NA ', & + 'SO4_A1 ','SO4_A2 ','SO4_AC ','SO4_PR ','BC_N ', & + 'BC_AX ','BC_NI ','BC_A ','BC_AI ','BC_AC ', & + 'OM_NI ','OM_AI ','OM_AC ','DST_A2 ','DST_A3 ', & + 'SS_A1 ','SS_A2 ','SS_A3 ','SOA_NA ','SOA_A1 ', & + 'SOA_LV ','SOA_SV ','monoterp ','isoprene ','H2O ' /) + + adv_mass(: 30) = (/ 64.064800_r8, 98.078400_r8, 62.132400_r8, 34.013600_r8, 98.078400_r8, & + 98.078400_r8, 115.107340_r8, 98.078400_r8, 98.078400_r8, 12.011000_r8, & + 12.011000_r8, 12.011000_r8, 12.011000_r8, 12.011000_r8, 12.011000_r8, & + 12.011000_r8, 12.011000_r8, 12.011000_r8, 135.064039_r8, 135.064039_r8, & + 58.442468_r8, 58.442468_r8, 58.442468_r8, 168.227200_r8, 168.227200_r8, & + 168.227200_r8, 168.227200_r8, 136.228400_r8, 68.114200_r8, 18.014200_r8 /) + + crb_mass(: 30) = (/ 0.000000_r8, 0.000000_r8, 24.022000_r8, 0.000000_r8, 0.000000_r8, & + 0.000000_r8, 0.000000_r8, 0.000000_r8, 0.000000_r8, 12.011000_r8, & + 12.011000_r8, 12.011000_r8, 12.011000_r8, 12.011000_r8, 12.011000_r8, & + 12.011000_r8, 12.011000_r8, 12.011000_r8, 0.000000_r8, 0.000000_r8, & + 0.000000_r8, 0.000000_r8, 0.000000_r8, 120.110000_r8, 120.110000_r8, & + 120.110000_r8, 120.110000_r8, 120.110000_r8, 60.055000_r8, 0.000000_r8 /) + + fix_mass(: 7) = (/ 0.00000000_r8, 28.0134800_r8, 31.9988000_r8, 47.9982000_r8, 17.0068000_r8, & + 62.0049400_r8, 33.0062000_r8 /) + + clsmap(: 30,4) = (/ 3, 1, 4, 5, 6, 7, 8, 9, 10, 11, & + 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, & + 22, 23, 2, 24, 25, 26, 27, 28, 29, 30 /) + + permute(: 30,4) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, & + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, & + 21, 22, 23, 24, 25, 26, 27, 28, 29, 30 /) + + diag_map(: 30) = (/ 1, 5, 7, 9, 10, 11, 12, 13, 14, 15, & + 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, & + 26, 27, 28, 29, 30, 31, 32, 35, 37, 38 /) + + extfrc_lst(: 7) = (/ 'SO2 ','BC_NI ','BC_AX ','BC_N ','OM_NI ', & + 'SO4_PR ','H2O ' /) + + frc_from_dataset(: 7) = (/ .true., .true., .true., .true., .true., & + .true., .true. /) + + inv_lst(: 7) = (/ 'M ', 'N2 ', 'O2 ', 'O3 ', 'OH ', & + 'NO3 ', 'HO2 ' /) + + if( allocated( rxt_tag_lst ) ) then + deallocate( rxt_tag_lst ) + end if + allocate( rxt_tag_lst(rxt_tag_cnt),stat=ios ) + if( ios /= 0 ) then + write(iulog,*) 'set_sim_dat: failed to allocate rxt_tag_lst; error = ',ios + call endrun + end if + if( allocated( rxt_tag_map ) ) then + deallocate( rxt_tag_map ) + end if + allocate( rxt_tag_map(rxt_tag_cnt),stat=ios ) + if( ios /= 0 ) then + write(iulog,*) 'set_sim_dat: failed to allocate rxt_tag_map; error = ',ios + call endrun + end if + rxt_tag_lst( 1: 3) = (/ 'jh2o2 ', 'usr_HO2_HO2 ', & + 'usr_DMS_OH ' /) + rxt_tag_map(:rxt_tag_cnt) = (/ 1, 2, 7 /) + if( allocated( pht_alias_lst ) ) then + deallocate( pht_alias_lst ) + end if + allocate( pht_alias_lst(phtcnt,2),stat=ios ) + if( ios /= 0 ) then + write(iulog,*) 'set_sim_dat: failed to allocate pht_alias_lst; error = ',ios + call endrun + end if + if( allocated( pht_alias_mult ) ) then + deallocate( pht_alias_mult ) + end if + allocate( pht_alias_mult(phtcnt,2),stat=ios ) + if( ios /= 0 ) then + write(iulog,*) 'set_sim_dat: failed to allocate pht_alias_mult; error = ',ios + call endrun + end if + pht_alias_lst(:,1) = (/ ' ' /) + pht_alias_lst(:,2) = (/ ' ' /) + pht_alias_mult(:,1) = (/ 1._r8 /) + pht_alias_mult(:,2) = (/ 1._r8 /) + allocate( num_rnts(rxntot-phtcnt),stat=ios ) + if( ios /= 0 ) then + write(iulog,*) 'set_sim_dat: failed to allocate num_rnts; error = ',ios + call endrun + end if + num_rnts(:) = (/ 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, & + 2, 2 /) + + end subroutine set_sim_dat + + end module mo_sim_dat