RRTMGP mods for MT configurations (ESCOMP#266)

Originator(s): JulioTBacmeister peverwhee cacraigucar

Description (include issue title and the keyword ['closes', 'fixes',
'resolves'] and issue number):
No issue here, but in ESCOMP/CAM
ESCOMP/CAM#1279: RRTMGP error with mpas and
L93 mid-top configuration

List all namelist files that were added or changed:
- Added p_top_for_rrtmgp which is supplied by the host model's radiation
namelist

List all files eliminated and why: N/A

List all files added and what they do: 
A       to_be_ccppized/nlte_fomichev.F90
              - Copied from ESCOMP/CAM
              - 
List all existing files that have been modified, and describe the
changes:
(Helpful git command: `git diff --name-status
development...<your_branch_name>`)
M       schemes/rrtmgp/utils/calculate_net_heating.F90
              - modifications to fix the MT configuration

List all automated tests that failed, as well as an explanation for why
they weren't fixed: N/A

Is this an answer-changing PR? If so, is it a new physics package,
algorithm change, tuning change, etc?
     - Answer changing for MT only in ESCOMP/CAM

If yes to the above question, describe how this code was validated with
the new/modified features: @JulioTBacmeister made scientific runs

---------

Co-authored-by: peverwhee <courtneyp@ucar.edu>
Co-authored-by: Julio Bacmeister <juliob@ucar.edu>

Author: 3 people

schemes/rrtmgp/rrtmgp_inputs.F90

@@ -13,7 +13,8 @@ module rrtmgp_inputs
   subroutine rrtmgp_inputs_init(ktopcam, ktoprad, nlaycam, sw_low_bounds, sw_high_bounds, nswbands,         &
                    pref_edge, nlay, pver, pverp, kdist_sw, kdist_lw, qrl, is_first_step, use_rad_dt_cosz,   &
                    timestep_size, nstep, iradsw, dt_avg, irad_always, is_first_restart_step, is_root,       &
-                   nlwbands, nradgas, gasnamelength, iulog, idx_sw_diag, idx_nir_diag, idx_uv_diag,         &
+                   p_top_for_rrtmgp,                                                                        &
+                   nlwbands, nradgas, gasnamelength, idx_sw_diag, idx_nir_diag, idx_uv_diag,                &
                    idx_sw_cloudsim, idx_lw_diag, idx_lw_cloudsim, nswgpts, nlwgpts, nlayp,                  &
                    nextsw_cday, current_cal_day, band2gpt_sw, errmsg, errflg)
      use ccpp_kinds,             only: kind_phys
@@ -26,10 +27,10 @@ subroutine rrtmgp_inputs_init(ktopcam, ktoprad, nlaycam, sw_low_bounds, sw_high_
      integer,                         intent(in) :: nradgas                ! Number of radiatively active gases
      integer,                         intent(in) :: pverp                  ! Number of vertical interfaces
      integer,                         intent(in) :: pver                   ! Number of vertical layers
-     integer,                         intent(in) :: iradsw                 ! Freq. of shortwave radiation calc in time steps (positive) or hours (negative).
+     integer,                         intent(in) :: iradsw                 ! Freq. of shortwave radiation calc in time steps
+                                                                           ! (positive) or hours (negative).
      integer,                         intent(in) :: timestep_size          ! Timestep size (s)
      integer,                         intent(in) :: nstep                  ! Current timestep number
-     integer,                         intent(in) :: iulog                  ! Logging unit
      integer,                         intent(in) :: gasnamelength          ! Length of all of the gas_list entries
      real(kind_phys),                 intent(in) :: current_cal_day        ! Current calendar day
      real(kind_phys), dimension(:),   intent(in) :: pref_edge              ! Reference pressures (interfaces) (Pa)
@@ -39,6 +40,7 @@ subroutine rrtmgp_inputs_init(ktopcam, ktoprad, nlaycam, sw_low_bounds, sw_high_
      logical,                         intent(in) :: is_first_restart_step  ! Flag for whether this is the first restart step (.true. = yes)
      logical,                         intent(in) :: use_rad_dt_cosz        ! Use adjusted radiation timestep for cosz calculation
      logical,                         intent(in) :: is_root                ! Flag for whether this is the root task
+     real(kind_phys),                 intent(in) :: p_top_for_rrtmgp       ! Top pressure to use for RRTMGP
 
      ! Outputs
      integer,                         intent(out) :: ktopcam               ! Index in CAM arrays of top level (layer or interface) at which RRTMGP is active
@@ -81,7 +83,7 @@ subroutine rrtmgp_inputs_init(ktopcam, ktoprad, nlaycam, sw_low_bounds, sw_high_
      ! pressure interfaces below 1 Pa.  When the entire model atmosphere is
      ! below 1 Pa then an extra layer is added to the top of the model for
      ! the purpose of the radiation calculation.
-     nlay = count( pref_edge(:) > 1._kind_phys ) ! pascals (0.01 mbar)
+     nlay = count( pref_edge(:) > p_top_for_rrtmgp )
      nlayp = nlay + 1
 
      if (nlay == pverp) then
@@ -90,16 +92,6 @@ subroutine rrtmgp_inputs_init(ktopcam, ktoprad, nlaycam, sw_low_bounds, sw_high_
         ktopcam = 1
         ktoprad = 2
         nlaycam = pver
-     else if (nlay == (pverp-1)) then
-        ! Special case nlay == (pverp-1) -- topmost interface outside bounds (CAM MT config), treat as if it is ok.
-        ktopcam = 1
-        ktoprad = 2
-        nlaycam = pver
-        nlay = nlay+1 ! reassign the value so later code understands to treat this case like nlay==pverp
-        if (is_root) then
-           write(iulog,*) 'RADIATION_INIT: Special case of 1 model interface at p < 1Pa. Top layer will be INCLUDED in radiation calculation.'
-           write(iulog,*) 'RADIATION_INIT: nlay = ',nlay, ' same as pverp: ',nlay==pverp
-        end if
      else
         ! nlay < pverp.  nlay layers are used in radiation calcs, and they are
         ! all CAM layers.
@@ -254,33 +246,49 @@ subroutine rrtmgp_inputs_run(dosw, dolw, snow_associated, graupel_associated, &
      ! to be consistent with t_sfc.
      emis_sfc(:,:) = 1._kind_phys
 
+     !-------------------------------------------------------------------------
+     ! RRTMGP enforces  P > 1 Pa for validity.
+     ! In radiation.F90 we count layers based on P_ref > 10 Pa to safely account
+     ! for possible situations in MPAS (z-based vert. coordinate) in which
+     ! full 3D pressure could be significanlty below min(P_ref).
+     !
+     ! If 
+     ! 1) entire vertical domain has P_ref> 10Pa (e.g. CAM7 LT) then
+     !    nlay = pverp
+     !    ktoprad = 2
+     !    ktopcam = 1
+     ! 2) min(P_ref) < 10Pa (e.g. CAM7 MT) then
+     !    nlay < pverp
+     !    ktoprad = 1
+     !    ktopcam = pver - nlay + 1
+     !
      ! Level ordering is the same for both CAM and RRTMGP (top to bottom)
+     ! Note in Case 2, tops of {t,pint,pmid}_rad start at a 'valid' level,
+     ! i.e. pref > p_top_for_rrtmgp
+     !-------------------------------------------------------
      t_rad(:,ktoprad:) = t(:,ktopcam:)
      pmid_rad(:,ktoprad:) = pmid(:,ktopcam:)
      pint_rad(:,ktoprad:) = pint(:,ktopcam:)
 
-     ! Add extra layer values if needed.
+     ! Deal with vertical grid for RRTMGP 
      if (nlay == pverp) then
-        t_rad(:,1) = t(:,1)
+        ! This case is the CAM7 LT situation, i.e., all model layers are
+        ! within RRTMGP's range of valid pressures - (Case 1 above)
+        t_rad(:,1)      = t(:,1)
         ! The top reference pressure from the RRTMGP coefficients datasets is 1.005183574463 Pa
         ! Set the top of the extra layer just below that.
         pint_rad(:,1) = 1.01_kind_phys
-
-        ! next interface down in LT will always be > 1Pa
-        ! but in MT we apply adjustment to have it be 1.02 Pa if it was too high
-        where (pint_rad(:,2) <= pint_rad(:,1)) pint_rad(:,2) = pint_rad(:,1)+0.01_kind_phys
-
-        ! set the highest pmid (in the "extra layer") to the midpoint (guarantees > 1Pa)
+        ! set the highest pmid (in the "extra layer") to the midpoint (guarantees > 1Pa) 
         pmid_rad(:,1)   = pint_rad(:,1) + 0.5_kind_phys * (pint_rad(:,2) - pint_rad(:,1))
-
-        ! For case of CAM MT, also ensure pint_rad(:,2) > pint_rad(:,1) & pmid_rad(:,2) > max(pmid_rad(:,1), min_pressure)
-        where (pmid_rad(:,2) <= kdist_sw%gas_props%get_press_min()) pmid_rad(:,2) = pint_rad(:,2) + 0.01_kind_phys
      else
-        ! nlay < pverp, thus the 1 Pa level is within a CAM layer.  Assuming the top interface of
-        ! this layer is at a pressure < 1 Pa, we need to adjust the top of this layer so that it
-        ! is within the valid pressure range of RRTMGP (otherwise RRTMGP issues an error).  Then
-        ! set the midpoint pressure halfway between the interfaces.
+        ! nlay < pverp : model min(pref) < p_top_for_rrtmgp  (Case 2 above)
+        ! min(pref) could be 9.999 or 0.0999  
+        ! Assuming the top interface of this layer is at a pressure < 1 Pa, we need to adjust
+        ! so that it is within the valid pressure range of RRTMGP (otherwise RRTMGP issues
+        ! an error).  Then we set the midpoint pressure halfway between the interfaces.
         pint_rad(:,1) = 1.01_kind_phys
+        ! The following *should* work since pint_rad is all in valid range.
+        ! Need to think about possible edge cases ... (jtb 07/31/25)
         pmid_rad(:,1) = 0.5_kind_phys * (pint_rad(:,1) + pint_rad(:,2))
      end if
 

schemes/rrtmgp/utils/calculate_net_heating.F90

@@ -23,8 +23,9 @@ module calculate_net_heating
 !> \section arg_table_calculate_net_heating_run Argument Table
 !! \htmlinclude calculate_net_heating_run.html
 !!
-subroutine calculate_net_heating_run(ncol, rad_heat, qrl_prime, qrs_prime, fsns, fsnt, &
-                flns, flnt, is_offline_dyn, net_flx, errmsg, errflg)
+subroutine calculate_net_heating_run(ncol, rad_heat, qrl_prime, qrs_prime,     &
+                is_offline_dyn, fsns, fsnt, flns, flnt, net_flx, &
+                errmsg, errflg)
 !-----------------------------------------------------------------------
 ! Compute net radiative heating from qrs and qrl, and the associated net
 ! boundary flux.
@@ -52,6 +53,7 @@ subroutine calculate_net_heating_run(ncol, rad_heat, qrl_prime, qrs_prime, fsns,
    ! Set error variables
    errmsg = ''
    errflg = 0
+
    if (.not. is_offline_dyn) then
       rad_heat(:,:) = (qrs_prime(:,:) + qrl_prime(:,:))
    end if
@@ -62,5 +64,4 @@ subroutine calculate_net_heating_run(ncol, rad_heat, qrl_prime, qrs_prime, fsns,
 
 end subroutine calculate_net_heating_run
 
-!================================================================================================
 end module calculate_net_heating