diff --git a/.github/workflows/srt.yml b/.github/workflows/srt.yml index 3afe2b7d1..361c561fd 100644 --- a/.github/workflows/srt.yml +++ b/.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 diff --git a/cesm/flux_atmocn/flux_atmocn_COARE_mod.F90 b/cesm/flux_atmocn/flux_atmocn_COARE_mod.F90 index 824d4097a..15e38d02a 100644 --- a/cesm/flux_atmocn/flux_atmocn_COARE_mod.F90 +++ b/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 diff --git a/cesm/flux_atmocn/flux_atmocn_Diurnal_mod.F90 b/cesm/flux_atmocn/flux_atmocn_Diurnal_mod.F90 index ed0dd9a4a..8e4106409 100644 --- a/cesm/flux_atmocn/flux_atmocn_Diurnal_mod.F90 +++ b/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) diff --git a/cesm/flux_atmocn/flux_atmocn_Large.F90 b/cesm/flux_atmocn/flux_atmocn_Large.F90 index d58d512ba..b0df413db 100644 --- a/cesm/flux_atmocn/flux_atmocn_Large.F90 +++ b/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) diff --git a/cesm/flux_atmocn/flux_atmocn_driver_mod.F90 b/cesm/flux_atmocn/flux_atmocn_driver_mod.F90 index 82a2b97d8..6bca47e05 100644 --- a/cesm/flux_atmocn/flux_atmocn_driver_mod.F90 +++ b/cesm/flux_atmocn/flux_atmocn_driver_mod.F90 @@ -25,7 +25,8 @@ 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 -------------------------------- @@ -33,6 +34,7 @@ subroutine flux_atmOcn_driver(logunit, nMax, & 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) diff --git a/cime_config/namelist_definition_drv.xml b/cime_config/namelist_definition_drv.xml index e33d9455f..b902ba356 100644 --- a/cime_config/namelist_definition_drv.xml +++ b/cime_config/namelist_definition_drv.xml @@ -938,11 +938,33 @@ MED_attributes atm/ocn flux calculation scheme + + 0: Large and Pond + 1: COARE algorithm + 2: UA algorithm 0 + + logical + control + MED_attributes + + 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.) + + + + .true. + + logical control diff --git a/mediator/med_phases_aofluxes_mod.F90 b/mediator/med_phases_aofluxes_mod.F90 index c0d0192b5..64de18785 100644 --- a/mediator/med_phases_aofluxes_mod.F90 +++ b/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