Skip to content

Commit 35f4fc1

Browse files
billsacksbillsacks
authored
Merge pull request #625 from billsacks/shr_flux_namelist
Add namelist control around the qsat calculation for atm-ocn fluxes ### Description of changes The calculation of qsat in the atm-ocn fluxes was [recently changed](#589) to use the shr_wv_sat module. However, NorESM would like to continue to use the old formulation for now. This PR reintroduces the old code and introduces a new config variable that controls whether to use the old or new formulation for qsat. The goal here is to allow NorESM and CESM to use the same code, so that Mariana doesn't need to maintain a separate fork of CMEPS for NorESM. For now, NorESM will need a one-line difference in the namelist definition file, setting the default value of aofluxes_use_shr_wv_sat to false instead of true. Note that this only applies to ocn_surface_flux_scheme = 0 (Large) or 1 (COARE) (but see #624 regarding extending this at some point). This PR also renames some arguments in the calls to shr_wv_sat_qsat_liquid to better match the argument names in that subroutine. ### Specific notes Contributors other than yourself, if any: CMEPS Issues Fixed (include github issue #): Are changes expected to change answers? (specify if bfb, different at roundoff, more substantial) No - bfb Any User Interface Changes (namelist or namelist defaults changes)? Adds new namelist flag, aofluxes_use_shr_wv_sat. ### Testing performed Ran these tests in the context of cesm3_0_alpha08a (but with share code updated to the latest, since that's a dependency of the latest CMEPS): ``` ERS_D_Ld3.ne30pg3_t232.B1850C_LTso.derecho_intel.allactive-defaultio SMS_Ld1.ne30_ne30_mg17.FC2010climo.derecho_intel.cam-outfrq1d ERP_D_Ln9.ne30pg3_ne30pg3_mt232.F1850C_MTso.derecho_intel.cam-outfrq9s ``` I ran 4 versions of these tests: - ocn_surface_flux_scheme = 0, aofluxes_use_shr_wv_sat = true (out-of-the-box defaults): bit-for-bit - ocn_surface_flux_scheme = 1, aofluxes_use_shr_wv_sat = true, comparing against baselines with ocn_surface_flux_scheme = 1: bit-for-bit - ocn_surface_flux_scheme = 0, aofluxes_use_shr_wv_sat = false: changes answers, as expected - ocn_surface_flux_scheme = 1, aofluxes_use_shr_wv_sat = false, comparing against baselines with ocn_surface_flux_scheme = 1: changes answers as expected
2 parents 8d3ff3a + 4a9842a commit 35f4fc1

7 files changed

Lines changed: 84 additions & 17 deletions

File tree

‎.github/workflows/srt.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ jobs:
1919
strategy:
2020
fail-fast: false
2121
matrix:
22-
python-version: [ 3.8, 3.11, 3.x ]
22+
python-version: [ "3.10", 3.12, 3.x ]
2323
env:
2424
CC: mpicc
2525
FC: mpifort

‎cesm/flux_atmocn/flux_atmocn_COARE_mod.F90‎

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,7 @@ subroutine flux_atmOcn_COARE( &
4545
ts, mask, seq_flux_atmocn_minwind, &
4646
sen, lat, lwup, evap, &
4747
taux ,tauy, tref, qref, &
48+
aofluxes_use_shr_wv_sat, &
4849
duu10n, ugust_out, u10res, &
4950
ustar_sv, re_sv, ssq_sv)
5051

@@ -53,6 +54,7 @@ subroutine flux_atmOcn_COARE( &
5354
real(R8) , intent(in) :: spval
5455
integer , intent(in) :: nMax ! data vector length
5556
integer , intent(in) :: mask (nMax) ! ocn domain mask 0 <=> out of domain
57+
logical , intent(in) :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
5658
real(R8) , intent(in) :: zbot (nMax) ! atm level height (m)
5759
real(R8) , intent(in) :: ubot (nMax) ! atm u wind (m/s)
5860
real(R8) , intent(in) :: vbot (nMax) ! atm v wind (m/s)
@@ -107,6 +109,8 @@ subroutine flux_atmOcn_COARE( &
107109
real(R8) :: hsb,hlb ! sens & lat heat flxs at zbot
108110
real(R8) :: tau ! stress at zbot
109111
real(R8) :: trf,qrf,urf,vrf ! reference-height quantities
112+
real(r8) :: esat_val ! value of esat (saturation vapor pressure) at this point
113+
real(r8) :: qsat_val ! value of qsat (saturation specific humidity) at this point
110114

111115
!--- local functions --------------------------------
112116
real(R8) :: qsat ! function: the saturation humididty of air (kg/m^3)
@@ -116,6 +120,9 @@ subroutine flux_atmOcn_COARE( &
116120
real(R8) :: tdiff(nMax) ! tbot - ts
117121
real(R8) :: vscl
118122

123+
!--- functions ---
124+
qsat(Tk) = 640380.0_R8 / exp(5107.4_R8/Tk)
125+
119126
!--- formats ----------------------------------------
120127
character(*),parameter :: subName = '(flux_atmOcn_COARE) '
121128
character(*),parameter :: F00 = "('(flux_atmOcn_COARE) ',4a)"
@@ -145,8 +152,15 @@ subroutine flux_atmOcn_COARE( &
145152
endif
146153
endif
147154

148-
call shr_wv_sat_qsat_liquid(ts(n), pslv(n), qsat, ssq)
149-
ssq = 0.98_R8 * ssq ! sea surf hum (kg/kg)
155+
if (aofluxes_use_shr_wv_sat) then
156+
! This version uses a qsat calculation method consistent with what's used in CAM
157+
call shr_wv_sat_qsat_liquid(ts(n), pslv(n), esat_val, qsat_val)
158+
ssq = 0.98_R8 * qsat_val ! sea surf hum (kg/kg)
159+
else
160+
! This version uses the qsat calculation method that was used for many years,
161+
! prior to Aug 2025, and which is still being used by default in NorESM
162+
ssq = 0.98_R8 * qsat(ts(n)) / rbot(n) ! sea surf hum (kg/kg)
163+
end if
150164

151165
call cor30a(ubot(n),vbot(n),tbot(n),qbot(n),rbot(n), & ! in atm params
152166
us(n),vs(n),ts(n),ssq, & ! in surf params

‎cesm/flux_atmocn/flux_atmocn_Diurnal_mod.F90‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,6 @@ module flux_atmocn_diurnal_mod
2626
use shr_const_mod, only : shr_const_ocn_ref_sal, shr_const_zsrflyr, shr_const_rgas
2727
use shr_sys_mod, only : shr_sys_abort
2828
use flux_atmocn_COARE_mod, only : cor30a
29-
use shr_wv_sat_mod, only : shr_wv_sat_qsat_liquid ! use saturation calculation consistent with CAM
3029

3130
implicit none
3231
private
@@ -236,7 +235,8 @@ subroutine flux_atmOcn_diurnal( &
236235
real(R8) :: tdiff(nMax) ! tbot - ts
237236
real(R8) :: vscl
238237

239-
! NOTE: this should use the shr_wv_sat_qsat_liquid if this routine is ever used in production
238+
! NOTE: this should use the shr_wv_sat_qsat_liquid if this routine is ever used in
239+
! production (see https://github.com/ESCOMP/CMEPS/issues/624)
240240
qsat(Tk) = 640380.0_R8 / exp(5107.4_R8/Tk)
241241
cdn(Umps) = 0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps
242242
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( &
354354
speed(n) = 0.0_R8
355355
endif
356356

357-
! This should be changed to use the subroutine below
357+
! This should be changed to use shr_wv_sat_qsat_liquid (see
358+
! https://github.com/ESCOMP/CMEPS/issues/624)
358359
ssq = 0.98_R8 * qsat(tBulk(n)) / rbot(n) ! sea surf hum (kg/kg)
359-
! call shr_wv_sat_qsat_liquid(tBulk(n), pslv(n), qsat, ssq)
360-
! ssq = 0.98_R8 * ssq ! sea surf hum (kg/kg)
361360

362361
delt = thbot(n) - tBulk(n) ! pot temp diff (K)
363362
delq = qbot(n) - ssq ! spec hum dif (kg/kg)
@@ -503,10 +502,9 @@ subroutine flux_atmOcn_diurnal( &
503502

504503
!--need to update ssq,delt,delq as function of tBulk ----
505504

506-
! This should be changed to use the subroutine below
505+
! This should be changed to use shr_wv_sat_qsat_liquid (see
506+
! https://github.com/ESCOMP/CMEPS/issues/624)
507507
ssq = 0.98_R8 * qsat(tBulk(n)) / rbot(n) ! sea surf hum (kg/kg)
508-
! call shr_wv_sat_qsat_liquid(tBulk(n), pslv(n), qsat, ssq)
509-
! ssq = 0.98_R8 * ssq ! sea surf hum (kg/kg)
510508

511509
delt = thbot(n) - tBulk(n) ! pot temp diff (K)
512510
delq = qbot(n) - ssq ! spec hum dif (kg/kg)

‎cesm/flux_atmocn/flux_atmocn_Large.F90‎

Lines changed: 16 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,8 @@ subroutine flux_atmOcn_large( &
4242
ts, mask, seq_flux_atmocn_minwind, &
4343
sen, lat, lwup, evap, &
4444
taux, tauy, tref, qref, &
45-
add_gusts, duu10n, ugust_out, u10res, &
45+
add_gusts, aofluxes_use_shr_wv_sat, &
46+
duu10n, ugust_out, u10res, &
4647
ustar_sv, re_sv, ssq_sv)
4748

4849
!--- input arguments --------------------------------
@@ -51,6 +52,7 @@ subroutine flux_atmOcn_large( &
5152
integer ,intent(in) :: nMax ! data vector length
5253
integer ,intent(in) :: mask (nMax) ! ocn domain mask 0 <=> out of domain
5354
logical ,intent(in) :: add_gusts
55+
logical ,intent(in) :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
5456
real(R8) ,intent(in) :: zbot (nMax) ! atm level height (m)
5557
real(R8) ,intent(in) :: ubot (nMax) ! atm u wind (m/s)
5658
real(R8) ,intent(in) :: vbot (nMax) ! atm v wind (m/s)
@@ -121,6 +123,8 @@ subroutine flux_atmOcn_large( &
121123
real(R8) :: cp ! specific heat of moist air
122124
real(R8) :: fac ! vertical interpolation factor
123125
real(R8) :: wind0 ! resolved large-scale 10m wind (no gust added)
126+
real(r8) :: esat_val ! value of esat (saturation vapor pressure) at this point
127+
real(r8) :: qsat_val ! value of qsat (saturation specific humidity) at this point
124128

125129
!--- local functions --------------------------------
126130
real(R8) :: qsat ! function: the saturation humididty of air (kg/m^3)
@@ -143,6 +147,8 @@ subroutine flux_atmOcn_large( &
143147
real(R8) :: gprec ! convective rainfall argument for ugust
144148
! -------------------------------------------------------------------------
145149

150+
qsat(Tk) = 640380.0_R8 / exp(5107.4_R8/Tk)
151+
146152
! Large and Yeager 2009
147153
cdn(Umps) = 0.0027_R8 / min(33.0000_R8,Umps) + 0.000142_R8 + &
148154
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( &
208214
endif
209215
endif
210216

211-
call shr_wv_sat_qsat_liquid(ts(n), pslv(n), qsat, ssq)
212-
ssq = 0.98_R8 * ssq ! sea surf hum (kg/kg)
217+
if (aofluxes_use_shr_wv_sat) then
218+
! This version uses a qsat calculation method consistent with what's used in CAM
219+
call shr_wv_sat_qsat_liquid(ts(n), pslv(n), esat_val, qsat_val)
220+
ssq = 0.98_R8 * qsat_val ! sea surf hum (kg/kg)
221+
else
222+
! This version uses the qsat calculation method that was used for many years,
223+
! prior to Aug 2025, and which is still being used by default in NorESM
224+
ssq = 0.98_R8 * qsat(ts(n)) / rbot(n) ! sea surf hum (kg/kg)
225+
end if
213226
delt = thbot(n) - ts(n) ! pot temp diff (K)
214227
delq = qbot(n) - ssq ! spec hum dif (kg/kg)
215228
alz = log(zbot(n)/zref)

‎cesm/flux_atmocn/flux_atmocn_driver_mod.F90‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -25,14 +25,16 @@ subroutine flux_atmOcn_driver(logunit, nMax, &
2525
sen, lat, lwup, evap, &
2626
taux, tauy, tref, qref, &
2727
ocn_surface_flux_scheme, &
28-
add_gusts, duu10n, ugust_out, u10res, &
28+
add_gusts, aofluxes_use_shr_wv_sat, &
29+
duu10n, ugust_out, u10res, &
2930
ustar_sv, re_sv, ssq_sv, missval)
3031

3132
!--- input arguments --------------------------------
3233
integer , intent(in) :: logunit
3334
integer , intent(in) :: nMax ! data vector length
3435
integer , intent(in) :: mask (nMax) ! ocn domain mask 0 <=> out of domain
3536
logical , intent(in) :: add_gusts
37+
logical , intent(in) :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
3638
real(R8) , intent(in) :: zbot (nMax) ! atm level height (m)
3739
real(R8) , intent(in) :: ubot (nMax) ! atm u wind (m/s)
3840
real(R8) , intent(in) :: vbot (nMax) ! atm v wind (m/s)
@@ -94,7 +96,8 @@ subroutine flux_atmOcn_driver(logunit, nMax, &
9496
ts, mask, seq_flux_atmocn_minwind, &
9597
sen, lat, lwup, evap, &
9698
taux, tauy, tref, qref, &
97-
add_gusts, duu10n, ugust_out, u10res, &
99+
add_gusts, aofluxes_use_shr_wv_sat, &
100+
duu10n, ugust_out, u10res, &
98101
ustar_sv=ustar_sv, re_sv=re_sv, ssq_sv=ssq_sv)
99102

100103
else if (ocn_surface_flux_scheme == ocn_flux_scheme_coare) then
@@ -107,6 +110,7 @@ subroutine flux_atmOcn_driver(logunit, nMax, &
107110
ts, mask, seq_flux_atmocn_minwind, &
108111
sen, lat, lwup, evap, &
109112
taux, tauy, tref, qref, &
113+
aofluxes_use_shr_wv_sat, &
110114
duu10n, ugust_out, u10res, &
111115
ustar_sv=ustar_sv, re_sv=re_sv, ssq_sv=ssq_sv)
112116

‎cime_config/namelist_definition_drv.xml‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -938,11 +938,33 @@
938938
<group>MED_attributes</group>
939939
<desc>
940940
atm/ocn flux calculation scheme
941+
942+
0: Large and Pond
943+
1: COARE algorithm
944+
2: UA algorithm
941945
</desc>
942946
<values>
943947
<value>0</value>
944948
</values>
945949
</entry>
950+
<entry id="aofluxes_use_shr_wv_sat">
951+
<type>logical</type>
952+
<category>control</category>
953+
<group>MED_attributes</group>
954+
<desc>
955+
If true, use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations.
956+
957+
If false, use the older inline calculation of qsat, which uses a different
958+
formulation.
959+
960+
(Currently only relevant for ocn_surface_flux_scheme = 0 or 1.)
961+
</desc>
962+
<values>
963+
<!-- Once CIME_MODEL distinguishes between NorESM and CESM, this value will differ
964+
between these two. -->
965+
<value>.true.</value>
966+
</values>
967+
</entry>
946968
<entry id="add_gusts">
947969
<type>logical</type>
948970
<category>control</category>

‎mediator/med_phases_aofluxes_mod.F90‎

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,7 @@ module med_phases_aofluxes_mod
8080
logical :: compute_atm_thbot
8181
integer :: ocn_surface_flux_scheme ! use case
8282
logical :: add_gusts
83+
logical :: aofluxes_use_shr_wv_sat ! use shr_wv_sat_mod to calculate qsat for atm-ocn flux calculations
8384

8485
character(len=CS), pointer :: fldnames_ocn_in(:)
8586
character(len=CS), pointer :: fldnames_atm_in(:)
@@ -419,6 +420,20 @@ subroutine med_aofluxes_init(gcomp, aoflux_in, aoflux_out, rc)
419420
add_gusts = .false.
420421
end if
421422

423+
call NUOPC_CompAttributeGet(gcomp, name='aofluxes_use_shr_wv_sat', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
424+
if (chkerr(rc,__LINE__,u_FILE_u)) return
425+
if (isPresent .and. isSet) then
426+
read(cvalue,*) aofluxes_use_shr_wv_sat
427+
else
428+
aofluxes_use_shr_wv_sat = .false.
429+
end if
430+
#ifdef CESMCOUPLED
431+
if (maintask) then
432+
write(logunit,*)
433+
write(logunit,'(a,l7)') trim(subname)//' aofluxes_use_shr_wv_sat = ', aofluxes_use_shr_wv_sat
434+
end if
435+
#endif
436+
422437
! bottom level potential temperature and/or botom level density
423438
! will need to be computed if not received from the atm
424439
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)
10821097
sen=aoflux_out%sen, lat=aoflux_out%lat, lwup=aoflux_out%lwup, evap=aoflux_out%evap, &
10831098
taux=aoflux_out%taux, tauy=aoflux_out%tauy, tref=aoflux_out%tref, qref=aoflux_out%qref, &
10841099
ocn_surface_flux_scheme=ocn_surface_flux_scheme, &
1085-
add_gusts=add_gusts, duu10n=aoflux_out%duu10n, ugust_out = aoflux_out%ugust_out, u10res = aoflux_out%u10res, &
1100+
add_gusts=add_gusts, aofluxes_use_shr_wv_sat=aofluxes_use_shr_wv_sat, &
1101+
duu10n=aoflux_out%duu10n, ugust_out = aoflux_out%ugust_out, u10res = aoflux_out%u10res, &
10861102
ustar_sv=aoflux_out%ustar, re_sv=aoflux_out%re, ssq_sv=aoflux_out%ssq, missval=0.0_r8)
10871103

10881104
#else

0 commit comments

Comments
 (0)