Extended non-LTE CO2 cooling parameterization

.gitmodules

@@ -107,6 +107,14 @@
         fxtag = clubb_4ncar_20240605_73d60f6_gpufixes_posinf
         fxDONOTUSEurl = https://github.com/larson-group/clubb_release
 
+[submodule "src/physics/ext_co2_cooling"]
+	path = src/physics/ext_co2_cooling
+	url = https://github.com/fedef17/CO2_cool_fort.git
+	fxrequired = AlwaysRequired
+	fxsparse = ../.ext_co2_cooling_sparse_checkout
+	fxtag = v1.1
+	fxDONOTUSEurl = https://github.com/fedef17/CO2_cool_fort.git
+
 [submodule "rtm"]
 path = components/rtm
 url = https://github.com/ESCOMP/RTM

bld/config_files/definition.xml

@@ -56,6 +56,9 @@ Model top specifier - set by compset definition
 <entry id="waccmx" valid_values="0,1" value="0">
 Option to turn on waccmx thermosphere/ionosphere extension: 0 => no, 1 => yes
 </entry>
+<entry id="ext_co2_cool" valid_values="0,1" value="0">
+Option to turn on extended non-LTE CO2 cooling: 0 => no, 1 => yes
+</entry>
 <entry id="ionosphere" valid_values="none,wxie" value="none">
 Ionosphere model used in WACCMX.
 </entry>

bld/configure

@@ -80,6 +80,7 @@ OPTIONS
      -cpl               Coupling framework [mct | nuopc].  Default: mct.
      -dyn <name>        Dynamical core option: [fv | se | fv3 | mpas].  Default: fv.
      -edit_chem_mech    Invokes CAMCHEM_EDITOR to allow the user to edit the chemistry mechanism file
+     -ext_co2_cool      Switch on extended non-LTE CO2 cooling
      -hgrid <name>      Specify horizontal grid.  Use nlatxnlon for spectral grids;
                         dlatxdlon for fv grids (dlat and dlon are the grid cell size
                         in degrees for latitude and longitude respectively); nexnp for
@@ -300,6 +301,7 @@ GetOptions(
     "version"                   => \$opts{'version'},
     "waccm_phys"                => \$opts{'waccm_phys'},
     "waccmx"                    => \$opts{'waccmx'},
+    "ext_co2_cool"              => \$opts{'ext_co2_cool'},
 )  or usage();
 
 # Give usage message.
@@ -706,6 +708,12 @@ if (defined $opts{'waccmx'}) {
 my $waccmx = $cfg_ref->get('waccmx');
 my $ionos  = $cfg_ref->get('ionosphere');
 
+# Extended CO2 non-LTE option
+if (defined $opts{'ext_co2_cool'}) {
+    $cfg_ref->set('ext_co2_cool', $opts{'ext_co2_cool'});
+}
+my $ext_co2_cool = $cfg_ref->get('ext_co2_cool');
+
 #-----------------------------------------------------------------------------------------------
 
 # Prognostic species package(s)
@@ -1768,6 +1776,11 @@ if ($waccmx) {
     $cfg_cppdefs .= ' -DWACCMX_PHYS';
 }
 
+# Extended Non-LTE CO2 cooling
+if ($ext_co2_cool) {
+    $cfg_cppdefs .= ' -DEXT_CO2_COOL';
+}
+
 # PIO
 my $pio2_build = (defined $opts{'pio2'}) ? 1 : 0;
 
@@ -1967,6 +1980,7 @@ sub write_filepath
     my $spmd          = $cfg_ref->get('spmd');
     my $offline_drv   = $cfg_ref->get('offline_drv');
     my $inic_val      = $cfg_ref->get('analytic_ic');
+    my $ext_co2_cool  = $cfg_ref->get('ext_co2_cool');
 
     # Root directory
     my $camsrcdir = $cfg_ref->get('cam_dir');
@@ -2058,6 +2072,10 @@ sub write_filepath
 	print $fh "$camsrcdir/src/physics/ali_arms/subs\n";
 	print $fh "$camsrcdir/src/physics/ali_arms/include\n";
     }
+    if ($ext_co2_cool) {
+        print $fh "$camsrcdir/src/physics/ext_co2_cooling/source/modules\n";
+    }	
+
     print $fh "$camsrcdir/src/ionosphere\n";
 
     print $fh "$camsrcdir/src/chemistry/mozart\n";

bld/namelist_files/namelist_defaults_cam.xml

@@ -182,7 +182,7 @@
 <ncdata hgrid="ne16np4"  nlev="126" waccmx="1"                                 >atm/waccm/ic/waccmx_ne16np4_126L_c200108.nc</ncdata>
 <ncdata hgrid="ne16np4"  nlev="126" waccmx="1" aquaplanet="1"                  >atm/waccm/ic/waccmx_aqua_ne16np4_126L_c191108.nc</ncdata>
 <ncdata hgrid="ne16np4"  nlev="126" waccmx="1" aquaplanet="1" ionosphere="none">atm/waccm/ic/waccmx4_neutral_aquap_ne16np4_126lev_c200827.nc</ncdata>
-<ncdata hgrid="ne16np4"  nlev="130" waccmx="1"     >atm/waccm/ic/waccmx_ne16pg3_L130_c250318.nc</ncdata>
+<ncdata hgrid="ne16np4"  nlev="130" waccmx="1"     >atm/waccm/ic/FX2000_ne16pg3_spinup01.cam.i.0002-01-01_c251215.nc</ncdata>
 
 <ncdata hgrid="ne30np4"  nlev="26"                 >atm/cam/inic/se/ape_topo_cam4_ne30np4_L26_c171020.nc</ncdata>
 <ncdata hgrid="ne30np4"  nlev="26"  aquaplanet="1" >atm/cam/inic/se/ape_cam4_ne30np4_L26_c170417.nc</ncdata>

bld/namelist_files/namelist_definition.xml

@@ -7994,6 +7994,13 @@ Call ALI-ARMS every X timesteps
 Default: 1
 </entry>
 
+<entry id="nlte_use_extco2" type="logical" category="waccm_phys"
+       group="radheat_nl" valid_values="" >
+If TRUE, then use the extended CO2 scheme of Lopez-Puertas et al. 2024 as an alernative method of computing non-LTE
+CO2 cooling rates in the upper atmosphere. (DOI:10.5194/gmd-17-4401-2024)
+Default: FALSE
+</entry>
+
 <entry id="qbo_cyclic" type="logical" category="waccm_phys"
        group="qbo_nl" valid_values="" >
 TRUE implies assume cyclic qbo data.

cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/shell_commands

@@ -1,2 +1,3 @@
 ./xmlchange ROF_NCPL=\$ATM_NCPL
 ./xmlchange GLC_NCPL=\$ATM_NCPL
+./xmlchange --append CAM_CONFIG_OPTS="-ext_co2_cool"

cime_config/testdefs/testmods_dirs/cam/outfrq9s_ctem/user_nl_cam

@@ -9,6 +9,9 @@ ctem_diags_numlats = 90
 
 fincl2 = 'THtem','Utem','Vtem','Wtem','VTHtem','WTHtem','UVtem', 'UWtem', 'MSKtem', 'PStem'
 fincl3 = 'Uzm', 'Vzm','Wzm', 'THzm','VTHzm','WTHzm','UVzm', 'UWzm', 'PSzm'
+fincl4 = 'QRL_TOT','QRLNLTE','QCO2','QO3','QHC2S','QNO','QO3P',
+         'QCO2ext','TCO2ext','PCO2ext','CO2_ext','N2_ext','O_ext','O2_ext'
 
 do_circulation_diags = .false.
 fincl7 = ' '
+nlte_use_extco2 = .true.

src/physics/.ext_co2_cooling_sparse_checkout

@@ -0,0 +1 @@
+source/modules

src/physics/waccm/nlte_extco2.F90

@@ -0,0 +1,141 @@
+!--------------------------------------------------------------------------------
+! Ref:
+!  Lopez-Puertas, M., Fabiano, F., Fomichev, V., Funke, B., and Marsh, D. R.:
+!  An improved and extended parameterization of the CO2 15 um cooling in the middle
+!  and upper atmosphere (CO2_cool_fort-1.0), Geosci. Model Dev., 17, 4401-4432,
+!  https://doi.org/10.5194/gmd-17-4401-2024, 2024.
+!--------------------------------------------------------------------------------
+module nlte_extco2
+
+  use ppgrid,       only: pcols, pver
+  use shr_kind_mod, only: r8 => shr_kind_r8
+  use cam_logfile,  only: iulog
+  use spmd_utils,   only: masterproc
+  use cam_history,  only: add_default, addfld, outfld
+
+  implicit none
+
+  private
+  public :: nlte_extco2_init
+  public :: nlte_extco2_hrate
+
+  real(r8) :: o1_mw_inv  = -huge(1.0_r8)       ! O molecular weight (inverse)
+  real(r8) :: o2_mw_inv  = -huge(1.0_r8)       ! O2 molecular weight (inverse)
+  real(r8) :: co2_mw_inv = -huge(1.0_r8)       ! CO2 molecular weight (inverse)
+  real(r8) :: n2_mw_inv  = -huge(1.0_r8)       ! N2 molecular weight (inverse)
+
+contains
+
+  !------------------------------------------------------------------------------
+  !------------------------------------------------------------------------------
+  subroutine nlte_extco2_init(co2_mw,n2_mw,o1_mw,o2_mw)
+    use ref_pres, only: pref_mid
+
+    ! Input variables
+    real(r8), intent(in) :: o1_mw     ! O molecular weight
+    real(r8), intent(in) :: o2_mw     ! O2 molecular weight
+    real(r8), intent(in) :: co2_mw    ! CO2 molecular weight
+    real(r8), intent(in) :: n2_mw     ! N2 molecular weight
+
+#ifdef EXT_CO2_COOL
+
+    co2_mw_inv = 1._r8/co2_mw
+    o1_mw_inv  = 1._r8/o1_mw
+    o2_mw_inv  = 1._r8/o2_mw
+    n2_mw_inv  = 1._r8/n2_mw
+
+    call addfld ('QCO2ext',   (/ 'lev' /), 'A','K/s','Extended CO2 cooling')
+    call addfld ('TCO2ext',   (/ 'lev' /), 'A','K','Temp used in Extended CO2 cooling')
+    call addfld ('PCO2ext',   (/ 'lev' /), 'A','hPa','Press used in Extended CO2 cooling')
+    call addfld ('CO2_ext',   (/ 'lev' /), 'A','mol/mol','CO2 vmr used in Extended CO2 cooling')
+    call addfld ('N2_ext',    (/ 'lev' /), 'A','mol/mol','N2 vmr used in Extended CO2 cooling')
+    call addfld ('O_ext',     (/ 'lev' /), 'A','mol/mol','O vmr used in Extended CO2 cooling')
+    call addfld ('O2_ext',    (/ 'lev' /), 'A','mol/mol','O2 vmr used in Extended CO2 cooling')
+
+#endif
+  end subroutine nlte_extco2_init
+
+
+  !------------------------------------------------------------------------------
+  !------------------------------------------------------------------------------
+  subroutine nlte_extco2_hrate(lchnk, ncol, temp, pres, co2mmr, n2mmr, ommr, o2mmr, &
+                               co2cooling)
+    use air_composition, only: mbarv
+#ifdef EXT_CO2_COOL
+    use co2cool, only: co2_nlte_cool
+#endif
+    ! Input variables
+    integer, intent(in) :: ncol           ! number of atmospheric columns
+    integer, intent(in) :: lchnk          ! chunk identifier
+
+    real(r8), intent(in) :: temp(:,:)     ! temperature (K)
+    real(r8), intent(in) :: pres(:,:)     ! pressure (Pa)
+    real(r8), intent(in) :: co2mmr(:,:)   ! CO2 mass mixing ratio
+    real(r8), intent(in) :: n2mmr(:,:)    ! N2 mass mixing ratio
+    real(r8), intent(in) :: ommr(:,:)     ! O mass mixing ratio
+    real(r8), intent(in) :: o2mmr(:,:)    ! O2 mass mixing ratio
+
+    ! Output
+    real(r8), intent(out) :: co2cooling(:,:)   ! K sec-1
+
+    ! Local vars
+    real(r8) :: co2vmr(pver)
+    real(r8) :: ovmr(pver)
+    real(r8) :: n2vmr(pver)
+    real(r8) :: o2vmr(pver)
+    real(r8) :: heatrate(pver) ! K/day
+    real(r8) :: hPa(pver)
+
+    real(r8) :: surf_temp ! K
+    integer :: icol
+
+    real(r8) :: tempout(pcols,pver)
+    real(r8) :: presout(pcols,pver)
+    real(r8) :: co2_out(pcols,pver)
+    real(r8) :: n2_out(pcols,pver)
+    real(r8) :: o_out(pcols,pver)
+    real(r8) :: o2_out(pcols,pver)
+
+    real(r8), parameter :: day_per_sec = 1._r8/86400._r8
+
+    co2cooling = 0._r8
+
+#ifdef EXT_CO2_COOL
+    do icol=1,ncol
+
+       ! Convert to VMR from mmr
+       co2vmr(:) = mbarv(icol,:,lchnk) * co2mmr(icol,:) * co2_mw_inv
+       ovmr(:)   = mbarv(icol,:,lchnk) * ommr(icol,:) * o1_mw_inv
+       n2vmr(:)  = mbarv(icol,:,lchnk) * n2mmr(icol,:) * n2_mw_inv
+       o2vmr(:)  = mbarv(icol,:,lchnk) * o2mmr(icol,:) * o2_mw_inv
+       hPa(:)    = pres(icol,:) * 1.e-2_r8 ! Pa --> hPa
+       surf_temp = temp(icol,pver)
+
+       tempout(icol,:) = temp(icol,:)
+       presout(icol,:) = hPa(:)
+       co2_out(icol,:) = co2vmr(:)
+       o_out(icol,:)   = ovmr(:)
+       o2_out(icol,:)  = o2vmr(:)
+       n2_out(icol,:)  = n2vmr(:)
+
+       heatrate(:) = 0._r8
+
+       ! Units: Temperature in K, pressure in hPa, vmrs in mol/mol, heating rate in K/day
+       call co2_nlte_cool(temp(icol,:), hPa, co2vmr, ovmr, o2vmr, n2vmr, pver, &
+                          surf_temp, heatrate) !   (K day-1)
+
+       co2cooling(icol,:) = heatrate(:) * day_per_sec ! K day-1 --> K sec-1
+
+    end do
+
+    call outfld ('QCO2ext', co2cooling(:ncol,:), ncol, lchnk)
+    call outfld ('TCO2ext', tempout(:ncol,:), ncol, lchnk)
+    call outfld ('PCO2ext', presout(:ncol,:), ncol, lchnk)
+    call outfld ('CO2_ext', co2_out(:ncol,:), ncol, lchnk)
+    call outfld ('O_ext',   o_out(:ncol,:),   ncol, lchnk)
+    call outfld ('O2_ext',  o2_out(:ncol,:),  ncol, lchnk)
+    call outfld ('N2_ext',  n2_out(:ncol,:),  ncol, lchnk)
+#endif
+  end subroutine nlte_extco2_hrate
+
+end module nlte_extco2

src/physics/waccm/nlte_lw.F90

@@ -11,9 +11,11 @@ module nlte_lw
 
   use nlte_fomichev,      only: nlte_fomichev_init, nlte_fomichev_calc, nocooling, o3pcooling
   use nlte_aliarms,       only: nlte_aliarms_init, nlte_aliarms_calc
+  use nlte_extco2,        only: nlte_extco2_init, nlte_extco2_hrate
 
   use waccm_forcing,      only: waccm_forcing_init, waccm_forcing_adv,  get_cnst
   use cam_logfile,        only: iulog
+  use cam_abortutils,     only: endrun
 
   implicit none
   private
@@ -36,6 +38,8 @@ module nlte_lw
   logical :: nlte_use_aliarms = .false.
   integer :: nlte_aliarms_every_X = 0
 
+  logical :: nlte_use_extco2 = .false.
+
   logical :: use_data_o3
   logical :: use_waccm_forcing = .false.
 
@@ -68,7 +72,8 @@ end subroutine nlte_register
 
 !================================================================================================
 
-  subroutine nlte_init (pref_mid, max_pressure_lw, nlte_use_mo_in, nlte_limit_co2, nlte_use_aliarms_in, nlte_aliarms_every_X_in)
+  subroutine nlte_init (pref_mid, max_pressure_lw, nlte_use_mo_in, nlte_limit_co2, nlte_use_aliarms_in, &
+       nlte_aliarms_every_X_in , nlte_use_extco2_in)
 !
 ! Initialize the nlte parameterizations and tgcm forcing data, if required
 !------------------------------------------------------------------------
@@ -84,6 +89,7 @@ subroutine nlte_init (pref_mid, max_pressure_lw, nlte_use_mo_in, nlte_limit_co2,
     logical,          intent(in) :: nlte_limit_co2
     logical,          intent(in) :: nlte_use_aliarms_in
     integer,          intent(in) :: nlte_aliarms_every_X_in
+    logical,          intent(in) :: nlte_use_extco2_in
 
 
     real(r8) :: o1_mw = -huge(1.0_r8)      ! O molecular weight
@@ -106,6 +112,11 @@ subroutine nlte_init (pref_mid, max_pressure_lw, nlte_use_mo_in, nlte_limit_co2,
     nlte_use_mo          = nlte_use_mo_in
     nlte_use_aliarms     = nlte_use_aliarms_in
     nlte_aliarms_every_X = nlte_aliarms_every_X_in
+    nlte_use_extco2      = nlte_use_extco2_in
+
+    if (nlte_use_aliarms .and. nlte_use_extco2) then
+       call endrun('nlte_init: cannot set both nlte_use_aliarms and nlte_use_extco2 set to TRUE')
+    end if
 
     ! ask rad_constituents module whether the O3 used in the climate
     ! calculation is from data
@@ -178,6 +189,8 @@ subroutine nlte_init (pref_mid, max_pressure_lw, nlte_use_mo_in, nlte_limit_co2,
 ! Initialize ALI-ARMS parameterization
     if (nlte_use_aliarms) then
        call nlte_aliarms_init (max_pressure_lw,co2_mw,n2_mw,o1_mw,o2_mw)
+    else if (nlte_use_extco2) then
+       call nlte_extco2_init(co2_mw,n2_mw,o1_mw,o2_mw)
     end if
 
 ! Initialize waccm forcing data
@@ -288,6 +301,7 @@ subroutine nlte_tend(state, pbuf, qrlf)
 
     integer :: k
     integer :: nstep
+    real(r8) :: co2cooling(pcols,pver)
 
 !------------------------------------------------------------------------
 
@@ -351,6 +365,16 @@ subroutine nlte_tend(state, pbuf, qrlf)
 
        call t_stopf('nlte_aliarms_calc')
 
+    else if (nlte_use_extco2) then
+
+       call t_startf('nlte_extco2_hrate')
+
+       call nlte_extco2_hrate(lchnk, ncol, state%t, state%pmid, xco2mmr, xn2mmr, xommr, xo2mmr, co2cooling)
+
+       qrlf(:ncol,:) = o3cool(:ncol,:) + co2cooling(:ncol,:) * cpairv(:ncol,:,lchnk)
+
+       call t_stopf('nlte_extco2_hrate')
+
     else
        qrlf(:ncol,:) = qrlfomichev(:ncol,:)
     end if

src/physics/waccm/radheat.F90

@@ -46,6 +46,7 @@ module radheat
   logical :: nlte_limit_co2 = .false. ! if true apply upper limit to co2 in the Fomichev scheme
   logical :: nlte_use_aliarms = .false. ! If true, use ALI-ARMS for the cooling rate calculation
   integer :: nlte_aliarms_every_X = 1 ! Call aliarms every X times radiation is called
+  logical :: nlte_use_extco2 = .false. ! If true, use the extended CO2 scheme of Lopez-Puertas et al. 2024
   real(r8), public :: p_top_for_equil_rad = 0._r8 ! Pressure top for blending layer
 
 ! Private variables for merging heating rates
@@ -88,7 +89,8 @@ subroutine radheat_readnl(nlfile)
     integer :: unitn, ierr
     character(len=*), parameter :: subname = 'radheat_readnl'
 
-    namelist /radheat_nl/ nlte_use_mo, nlte_limit_co2, nlte_use_aliarms,nlte_aliarms_every_X, p_top_for_equil_rad
+    namelist /radheat_nl/ nlte_use_mo, nlte_limit_co2, nlte_use_aliarms,nlte_aliarms_every_X, &
+                          nlte_use_extco2, p_top_for_equil_rad
 
     if (masterproc) then
        unitn = getunit()
@@ -113,9 +115,17 @@ subroutine radheat_readnl(nlfile)
     if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_use_aliarms")
     call mpi_bcast (nlte_aliarms_every_X, 1, mpi_integer, masterprocid, mpicom, ierr)
     if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_aliarms_every_X")
+    call mpi_bcast (nlte_use_extco2, 1, mpi_logical, masterprocid, mpicom, ierr)
+    if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: nlte_use_extco2")
     call mpi_bcast(p_top_for_equil_rad, 1, mpi_real8, masterprocid, mpicom, ierr)
     if (ierr /= 0) call endrun("radheat_readnl: FATAL: mpi_bcast: p_top_for_equil_rad")
 
+#ifndef EXT_CO2_COOL
+    if (nlte_use_extco2) then
+       call endrun(subname // ':: model not configured for extended non-LTE CO2 coooling')
+    end if
+#endif
+
     ! Have waccm_forcing read its namelist as well.
     call waccm_forcing_readnl(nlfile)
 
@@ -246,7 +256,8 @@ subroutine radheat_init(pref_mid)
     end if
 
     if (waccm_heating) then
-       call nlte_init(pref_mid, max_pressure_lw, nlte_use_mo, nlte_limit_co2, nlte_use_aliarms,nlte_aliarms_every_X)
+       call nlte_init(pref_mid, max_pressure_lw, nlte_use_mo, nlte_limit_co2, nlte_use_aliarms, &
+                      nlte_aliarms_every_X, nlte_use_extco2)
     endif
 
 ! Add history variables to master field list