Add namelist control around the qsat calculation for atm-ocn fluxes

.github/workflows/srt.yml

@@ -19,7 +19,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [ 3.8, 3.11, 3.x ]
+        python-version: [ "3.10", 3.12, 3.x ]
     env:
       CC: mpicc
       FC: mpifort

cesm/flux_atmocn/flux_atmocn_COARE_mod.F90

@@ -45,6 +45,7 @@ subroutine flux_atmOcn_COARE(                       &
        ts, mask, seq_flux_atmocn_minwind,             &
        sen, lat, lwup, evap,                          &
        taux ,tauy, tref, qref,                        &
+       aofluxes_use_shr_wv_sat,                       &
        duu10n, ugust_out, u10res,                     &
        ustar_sv, re_sv, ssq_sv)
 
@@ -53,6 +54,7 @@ subroutine flux_atmOcn_COARE(                       &
     real(R8) , intent(in) :: spval
     integer  , intent(in) :: nMax        ! data vector length
     integer  , intent(in) :: mask (nMax) ! ocn domain mask       0 <=> out of domain
+    logical  , intent(in) :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
     real(R8) , intent(in) :: zbot (nMax) ! atm level height      (m)
     real(R8) , intent(in) :: ubot (nMax) ! atm u wind            (m/s)
     real(R8) , intent(in) :: vbot (nMax) ! atm v wind            (m/s)
@@ -107,6 +109,8 @@ subroutine flux_atmOcn_COARE(                       &
     real(R8)    :: hsb,hlb         ! sens & lat heat flxs at zbot
     real(R8)    :: tau             ! stress at zbot
     real(R8)    :: trf,qrf,urf,vrf ! reference-height quantities
+    real(r8)    :: esat_val        ! value of esat (saturation vapor pressure) at this point
+    real(r8)    :: qsat_val        ! value of qsat (saturation specific humidity) at this point
 
     !--- local functions --------------------------------
     real(R8)    :: qsat   ! function: the saturation humididty of air (kg/m^3)
@@ -116,6 +120,9 @@ subroutine flux_atmOcn_COARE(                       &
     real(R8)    :: tdiff(nMax) ! tbot - ts
     real(R8)    :: vscl
 
+    !--- functions ---
+    qsat(Tk)   = 640380.0_R8 / exp(5107.4_R8/Tk)
+
     !--- formats ----------------------------------------
     character(*),parameter :: subName = '(flux_atmOcn_COARE) '
     character(*),parameter ::   F00 = "('(flux_atmOcn_COARE) ',4a)"
@@ -145,8 +152,15 @@ subroutine flux_atmOcn_COARE(                       &
              endif
           endif
 
-          call shr_wv_sat_qsat_liquid(ts(n), pslv(n), qsat, ssq)
-          ssq = 0.98_R8 * ssq   ! sea surf hum (kg/kg)
+          if (aofluxes_use_shr_wv_sat) then
+             ! This version uses a qsat calculation method consistent with what's used in CAM
+             call shr_wv_sat_qsat_liquid(ts(n), pslv(n), esat_val, qsat_val)
+             ssq = 0.98_R8 * qsat_val   ! sea surf hum (kg/kg)
+          else
+             ! This version uses the qsat calculation method that was used for many years,
+             ! prior to Aug 2025, and which is still being used by default in NorESM
+             ssq = 0.98_R8 * qsat(ts(n)) / rbot(n)   ! sea surf hum (kg/kg)
+          end if
 
           call cor30a(ubot(n),vbot(n),tbot(n),qbot(n),rbot(n), & ! in atm params
                us(n),vs(n),ts(n),ssq,                          & ! in surf params

cesm/flux_atmocn/flux_atmocn_Diurnal_mod.F90

@@ -26,7 +26,6 @@ module flux_atmocn_diurnal_mod
    use shr_const_mod,         only : shr_const_ocn_ref_sal, shr_const_zsrflyr, shr_const_rgas
    use shr_sys_mod,           only : shr_sys_abort
    use flux_atmocn_COARE_mod, only : cor30a
-   use shr_wv_sat_mod,        only : shr_wv_sat_qsat_liquid ! use saturation calculation consistent with CAM
 
    implicit none
    private
@@ -236,7 +235,8 @@ subroutine flux_atmOcn_diurnal(                      &
     real(R8)    :: tdiff(nMax)               ! tbot - ts
     real(R8)    :: vscl
 
-    ! NOTE: this should use the shr_wv_sat_qsat_liquid if this routine is ever used in production
+    ! NOTE: this should use the shr_wv_sat_qsat_liquid if this routine is ever used in
+    ! production (see https://github.com/ESCOMP/CMEPS/issues/624)
     qsat(Tk)   = 640380.0_R8 / exp(5107.4_R8/Tk)
     cdn(Umps)  =   0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps
     psimhu(xd) = log((1.0_R8+xd*(2.0_R8+xd))*(1.0_R8+xd*xd)/8.0_R8) - 2.0_R8*atan(xd) + 1.571_R8
@@ -354,10 +354,9 @@ subroutine flux_atmOcn_diurnal(                      &
              speed(n)   = 0.0_R8
           endif
 
-          ! This should be changed to use the subroutine below
+          ! This should be changed to use shr_wv_sat_qsat_liquid (see
+          ! https://github.com/ESCOMP/CMEPS/issues/624)
           ssq = 0.98_R8 * qsat(tBulk(n)) / rbot(n)   ! sea surf hum (kg/kg)
-          ! call shr_wv_sat_qsat_liquid(tBulk(n), pslv(n), qsat, ssq)
-          ! ssq    = 0.98_R8 * ssq                   ! sea surf hum (kg/kg)
 
           delt   = thbot(n) - tBulk(n)               ! pot temp diff (K)
           delq   = qbot(n) - ssq                     ! spec hum dif (kg/kg)
@@ -503,10 +502,9 @@ subroutine flux_atmOcn_diurnal(                      &
 
              !--need to update ssq,delt,delq as function of tBulk ----
 
-             ! This should be changed to use the subroutine below
+             ! This should be changed to use shr_wv_sat_qsat_liquid (see
+             ! https://github.com/ESCOMP/CMEPS/issues/624)
              ssq = 0.98_R8 * qsat(tBulk(n)) / rbot(n)      ! sea surf hum (kg/kg)
-             ! call shr_wv_sat_qsat_liquid(tBulk(n), pslv(n), qsat, ssq)
-             ! ssq    = 0.98_R8 * ssq                      ! sea surf hum (kg/kg)
 
              delt   = thbot(n) - tBulk(n)                  ! pot temp diff (K)
              delq   = qbot(n) - ssq                        ! spec hum dif (kg/kg)

cesm/flux_atmocn/flux_atmocn_Large.F90

@@ -42,7 +42,8 @@ subroutine flux_atmOcn_large(              &
        ts, mask,  seq_flux_atmocn_minwind,   &
        sen, lat, lwup, evap,                 &
        taux, tauy, tref, qref,               &
-       add_gusts, duu10n, ugust_out, u10res, &
+       add_gusts, aofluxes_use_shr_wv_sat,   &
+       duu10n, ugust_out, u10res, &
        ustar_sv, re_sv, ssq_sv)
 
     !--- input arguments --------------------------------
@@ -51,6 +52,7 @@ subroutine flux_atmOcn_large(              &
     integer  ,intent(in) :: nMax        ! data vector length
     integer  ,intent(in) :: mask (nMax) ! ocn domain mask       0 <=> out of domain
     logical  ,intent(in) :: add_gusts
+    logical  ,intent(in) :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
     real(R8) ,intent(in) :: zbot (nMax) ! atm level height (m)
     real(R8) ,intent(in) :: ubot (nMax) ! atm u wind (m/s)
     real(R8) ,intent(in) :: vbot (nMax) ! atm v wind (m/s)
@@ -121,6 +123,8 @@ subroutine flux_atmOcn_large(              &
     real(R8) :: cp     ! specific heat of moist air
     real(R8) :: fac    ! vertical interpolation factor
     real(R8) :: wind0  ! resolved large-scale 10m wind (no gust added)
+    real(r8) :: esat_val ! value of esat (saturation vapor pressure) at this point
+    real(r8) :: qsat_val ! value of qsat (saturation specific humidity) at this point
 
     !--- local functions --------------------------------
     real(R8) :: qsat   ! function: the saturation humididty of air (kg/m^3)
@@ -143,6 +147,8 @@ subroutine flux_atmOcn_large(              &
     real(R8) :: gprec ! convective rainfall argument for ugust
     ! -------------------------------------------------------------------------
 
+    qsat(Tk)   = 640380.0_R8 / exp(5107.4_R8/Tk)
+
     ! Large and Yeager 2009
     cdn(Umps)  =  0.0027_R8 / min(33.0000_R8,Umps) + 0.000142_R8 + &
          0.0000764_R8 * min(33.0000_R8,Umps) - 3.14807e-13_r8 * min(33.0000_R8,Umps)**6
@@ -208,8 +214,15 @@ subroutine flux_atmOcn_large(              &
              endif
           endif
 
-          call shr_wv_sat_qsat_liquid(ts(n), pslv(n), qsat, ssq)
-          ssq    = 0.98_R8 * ssq                     ! sea surf hum (kg/kg)
+          if (aofluxes_use_shr_wv_sat) then
+             ! This version uses a qsat calculation method consistent with what's used in CAM
+             call shr_wv_sat_qsat_liquid(ts(n), pslv(n), esat_val, qsat_val)
+             ssq    = 0.98_R8 * qsat_val                ! sea surf hum (kg/kg)
+          else
+             ! This version uses the qsat calculation method that was used for many years,
+             ! prior to Aug 2025, and which is still being used by default in NorESM
+             ssq    = 0.98_R8 * qsat(ts(n)) / rbot(n)   ! sea surf hum (kg/kg)
+          end if
           delt   = thbot(n) - ts(n)                  ! pot temp diff (K)
           delq   = qbot(n) - ssq                     ! spec hum dif (kg/kg)
           alz    = log(zbot(n)/zref)

cesm/flux_atmocn/flux_atmocn_driver_mod.F90

@@ -25,14 +25,16 @@ subroutine flux_atmOcn_driver(logunit, nMax, &
        sen, lat, lwup, evap,                   &
        taux, tauy, tref, qref,                 &
        ocn_surface_flux_scheme,                &
-       add_gusts, duu10n, ugust_out, u10res,   &
+       add_gusts, aofluxes_use_shr_wv_sat,     &
+       duu10n, ugust_out, u10res,              &
        ustar_sv, re_sv, ssq_sv, missval)
 
     !--- input arguments --------------------------------
     integer  , intent(in) :: logunit
     integer  , intent(in) :: nMax        ! data vector length
     integer  , intent(in) :: mask (nMax) ! ocn domain mask       0 <=> out of domain
     logical  , intent(in) :: add_gusts
+    logical  , intent(in) :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
     real(R8) , intent(in) :: zbot (nMax) ! atm level height (m)
     real(R8) , intent(in) :: ubot (nMax) ! atm u wind (m/s)
     real(R8) , intent(in) :: vbot (nMax) ! atm v wind (m/s)
@@ -94,7 +96,8 @@ subroutine flux_atmOcn_driver(logunit, nMax, &
             ts, mask,  seq_flux_atmocn_minwind,   &
             sen, lat, lwup,  evap,                &
             taux, tauy, tref, qref,               &
-            add_gusts, duu10n, ugust_out, u10res, &
+            add_gusts, aofluxes_use_shr_wv_sat,   &
+            duu10n, ugust_out, u10res, &
             ustar_sv=ustar_sv, re_sv=re_sv, ssq_sv=ssq_sv)
 
     else if (ocn_surface_flux_scheme == ocn_flux_scheme_coare) then
@@ -107,6 +110,7 @@ subroutine flux_atmOcn_driver(logunit, nMax, &
             ts, mask,  seq_flux_atmocn_minwind, &
             sen, lat, lwup, evap,               &
             taux, tauy, tref, qref,             &
+            aofluxes_use_shr_wv_sat,            &
             duu10n, ugust_out, u10res,          &
             ustar_sv=ustar_sv, re_sv=re_sv, ssq_sv=ssq_sv)
 

cime_config/namelist_definition_drv.xml

@@ -938,11 +938,33 @@
     <group>MED_attributes</group>
     <desc>
       atm/ocn flux calculation scheme
+
+      0: Large and Pond
+      1: COARE algorithm
+      2: UA algorithm
     </desc>
     <values>
       <value>0</value>
     </values>
   </entry>
+  <entry id="aofluxes_use_shr_wv_sat">
+    <type>logical</type>
+    <category>control</category>
+    <group>MED_attributes</group>
+    <desc>
+      If true, use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations.
+
+      If false, use the older inline calculation of qsat, which uses a different
+      formulation.
+
+      (Currently only relevant for ocn_surface_flux_scheme = 0 or 1.)
+    </desc>
+    <values>
+      <!-- Once CIME_MODEL distinguishes between NorESM and CESM, this value will differ
+      between these two. -->
+      <value>.true.</value>
+    </values>
+  </entry>
   <entry id="add_gusts">
     <type>logical</type>
     <category>control</category>

mediator/med_phases_aofluxes_mod.F90

@@ -80,6 +80,7 @@ module med_phases_aofluxes_mod
   logical :: compute_atm_thbot
   integer :: ocn_surface_flux_scheme ! use case
   logical :: add_gusts
+  logical :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
 
   character(len=CS), pointer :: fldnames_ocn_in(:)
   character(len=CS), pointer :: fldnames_atm_in(:)
@@ -419,6 +420,20 @@ subroutine med_aofluxes_init(gcomp, aoflux_in, aoflux_out, rc)
        add_gusts = .false.
     end if
 
+    call NUOPC_CompAttributeGet(gcomp, name='aofluxes_use_shr_wv_sat', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
+    if (chkerr(rc,__LINE__,u_FILE_u)) return
+    if (isPresent .and. isSet) then
+       read(cvalue,*) aofluxes_use_shr_wv_sat
+    else
+       aofluxes_use_shr_wv_sat = .false.
+    end if
+#ifdef CESMCOUPLED
+    if (maintask) then
+       write(logunit,*)
+       write(logunit,'(a,l7)') trim(subname)//' aofluxes_use_shr_wv_sat = ', aofluxes_use_shr_wv_sat
+    end if
+#endif
+
     ! bottom level potential temperature and/or botom level density
     ! will need to be computed if not received from the atm
     if (FB_fldchk(is_local%Wrap%FBImp(Compatm,Compatm), 'Sa_ptem', rc=rc)) then
@@ -1082,7 +1097,8 @@ subroutine med_aofluxes_update(gcomp, aoflux_in, aoflux_out, rc)
          sen=aoflux_out%sen, lat=aoflux_out%lat, lwup=aoflux_out%lwup, evap=aoflux_out%evap, &
          taux=aoflux_out%taux, tauy=aoflux_out%tauy, tref=aoflux_out%tref, qref=aoflux_out%qref, &
          ocn_surface_flux_scheme=ocn_surface_flux_scheme, &
-         add_gusts=add_gusts, duu10n=aoflux_out%duu10n, ugust_out = aoflux_out%ugust_out, u10res = aoflux_out%u10res, &
+         add_gusts=add_gusts, aofluxes_use_shr_wv_sat=aofluxes_use_shr_wv_sat, &
+         duu10n=aoflux_out%duu10n, ugust_out = aoflux_out%ugust_out, u10res = aoflux_out%u10res, &
          ustar_sv=aoflux_out%ustar, re_sv=aoflux_out%re, ssq_sv=aoflux_out%ssq, missval=0.0_r8)
 
 #else