Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into test_mle_diag_v3
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Feb 9, 2024
2 parents cd69215 + 4757bcb commit dd160e0
Show file tree
Hide file tree
Showing 3 changed files with 78 additions and 32 deletions.
9 changes: 4 additions & 5 deletions src/diagnostics/MOM_wave_speed.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,10 +197,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N
enddo ; enddo ; enddo
endif

g_Rho0 = GV%g_Earth*GV%H_to_Z / GV%Rho0
nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
H_to_pres = GV%H_to_RZ * GV%g_Earth
! Note that g_Rho0 = H_to_pres / GV%Rho0**2
nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
if (.not.nonBous) g_Rho0 = GV%g_Earth*GV%H_to_Z / GV%Rho0
use_EOS = associated(tv%eqn_of_state)

better_est = CS%better_cg1_est
Expand Down Expand Up @@ -900,9 +900,9 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
is = G%isc - halo ; ie = G%iec + halo ; js = G%jsc - halo ; je = G%jec + halo
endif

g_Rho0 = GV%g_Earth * GV%H_to_Z / GV%Rho0
H_to_pres = GV%H_to_RZ * GV%g_Earth
nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
H_to_pres = GV%H_to_RZ * GV%g_Earth
if (.not.nonBous) g_Rho0 = GV%g_Earth * GV%H_to_Z / GV%Rho0
use_EOS = associated(tv%eqn_of_state)

if (CS%c1_thresh < 0.0) &
Expand Down Expand Up @@ -1057,7 +1057,6 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
enddo
endif
endif
cg1_est = g_Rho0 * drxh_sum
else ! Not use_EOS
drxh_sum = 0.0 ; dSpVxh_sum = 0.0
if (better_est) then
Expand Down
98 changes: 73 additions & 25 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ module MOM_mixed_layer_restrat
real :: nstar !< The n* value used to estimate the turbulent vertical momentum flux [nondim]
real :: min_wstar2 !< The minimum lower bound to apply to the vertical momentum flux,
!! w'u', in the Bodner et al., restratification parameterization
!! [m2 s-2]. This avoids a division-by-zero in the limit when u*
!! [Z2 T-2 ~> m2 s-2]. This avoids a division-by-zero in the limit when u*
!! and the buoyancy flux are zero.
real :: BLD_growing_Tfilt !< The time-scale for a running-mean filter applied to the boundary layer
!! depth (BLD) when the BLD is deeper than the running mean [T ~> s].
Expand All @@ -82,6 +82,11 @@ module MOM_mixed_layer_restrat
real :: MLD_growing_Tfilt !< The time-scale for a running-mean filter applied to the time-filtered
!! MLD, when the latter is deeper than the running mean [T ~> s].
!! A value of 0 instantaneously sets the running mean to the current value of MLD.
integer :: answer_date !< The vintage of the order of arithmetic and expressions in the
!! mixed layer restrat calculations. Values below 20240201 recover
!! the answers from the end of 2023, while higher values use the new
!! cuberoot function in the Bodner code to avoid needing to undo
!! dimensional rescaling.

logical :: debug = .false. !< If true, calculate checksums of fields for debugging.

Expand Down Expand Up @@ -281,7 +286,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
!! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA
pRef_MLD(:) = 0.
EOSdom(:) = EOS_domain(G%HI, halo=1)
do j = js-1, je+1
do j=js-1,je+1
dK(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer
if (CS%use_Stanley_ML) then
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, tv%varT(:,j,1), covTS, varS, &
Expand All @@ -291,7 +296,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
endif
deltaRhoAtK(:) = 0.
MLD_fast(:,j) = 0.
do k = 2, nz
do k=2,nz
dKm1(:) = dK(:) ! Depth of center of layer K-1
dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Depth of center of layer K
! Mixed-layer depth, using sigma-0 (surface reference pressure)
Expand All @@ -302,10 +307,10 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
else
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, EOSdom)
endif
do i = is-1,ie+1
do i=is-1,ie+1
deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
enddo
do i = is-1, ie+1
do i=is-1,ie+1
ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i)
if ((MLD_fast(i,j)==0.) .and. (ddRho>0.) .and. &
(deltaRhoAtKm1(i)<CS%MLE_density_diff) .and. (deltaRhoAtK(i)>=CS%MLE_density_diff)) then
Expand All @@ -314,15 +319,15 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
endif
enddo ! i-loop
enddo ! k-loop
do i = is-1, ie+1
do i=is-1,ie+1
MLD_fast(i,j) = CS%MLE_MLD_stretch * MLD_fast(i,j)
if ((MLD_fast(i,j)==0.) .and. (deltaRhoAtK(i)<CS%MLE_density_diff)) &
MLD_fast(i,j) = dK(i) ! Assume mixing to the bottom
enddo
enddo ! j-loop
elseif (CS%MLE_use_PBL_MLD) then
if (GV%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
do j = js-1, je+1 ; do i = is-1, ie+1
do j=js-1,je+1 ; do i=is-1,ie+1
MLD_fast(i,j) = CS%MLE_MLD_stretch * GV%Z_to_H * MLD_in(i,j)
enddo ; enddo
else ! The fully non-Boussinesq conversion between height in MLD_in and thickness.
Expand Down Expand Up @@ -358,7 +363,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
endif
aFac = CS%MLE_MLD_decay_time / ( dt + CS%MLE_MLD_decay_time )
bFac = dt / ( dt + CS%MLE_MLD_decay_time )
do j = js-1, je+1 ; do i = is-1, ie+1
do j=js-1,je+1 ; do i=is-1,ie+1
! Expression bFac*MLD_fast(i,j) + aFac*CS%MLD_filtered(i,j) is the time-filtered
! (running mean) of MLD. The max() allows the "running mean" to be reset
! instantly to a deeper MLD.
Expand All @@ -375,15 +380,15 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
endif
aFac = CS%MLE_MLD_decay_time2 / ( dt + CS%MLE_MLD_decay_time2 )
bFac = dt / ( dt + CS%MLE_MLD_decay_time2 )
do j = js-1, je+1 ; do i = is-1, ie+1
do j=js-1,je+1 ; do i=is-1,ie+1
! Expression bFac*MLD_fast(i,j) + aFac*CS%MLD_filtered(i,j) is the time-filtered
! (running mean) of MLD. The max() allows the "running mean" to be reset
! instantly to a deeper MLD.
CS%MLD_filtered_slow(i,j) = max( MLD_fast(i,j), bFac*MLD_fast(i,j) + aFac*CS%MLD_filtered_slow(i,j) )
MLD_slow(i,j) = CS%MLD_filtered_slow(i,j)
enddo ; enddo
else
do j = js-1, je+1 ; do i = is-1, ie+1
do j=js-1,je+1 ; do i=is-1,ie+1
MLD_slow(i,j) = MLD_fast(i,j)
enddo ; enddo
endif
Expand Down Expand Up @@ -822,8 +827,8 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
real :: g_Rho0 ! G_Earth/Rho0 times a thickness conversion factor
! [L2 H-1 T-2 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2]
real :: h_vel ! htot interpolated onto velocity points [H ~> m or kg m-2]
real :: w_star3 ! Cube of turbulent convective velocity [m3 s-3]
real :: u_star3 ! Cube of surface fruction velocity [m3 s-3]
real :: w_star3 ! Cube of turbulent convective velocity [Z3 T-3 ~> m3 s-3]
real :: u_star3 ! Cube of surface friction velocity [Z3 T-3 ~> m3 s-3]
real :: r_wpup ! reciprocal of vertical momentum flux [T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1]
real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
real :: f2_h ! Coriolis parameter at h-points squared [T-2 ~> s-2]
Expand All @@ -842,6 +847,10 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
real :: muza ! mu(z) at top of the layer [nondim]
real :: dh ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2]
real :: res_scaling_fac ! The resolution-dependent scaling factor [nondim]
real :: Z3_T3_to_m3_s3 ! Conversion factors to undo scaling and permit terms to be raised to a
! fractional power [T3 m3 Z-3 s-3 ~> 1]
real :: m2_s2_to_Z2_T2 ! Conversion factors to restore scaling after a term is raised to a
! fractional power [Z2 s2 T-2 m-2 ~> 1]
real, parameter :: two_thirds = 2./3. ! [nondim]
logical :: line_is_empty, keep_going
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
Expand Down Expand Up @@ -889,7 +898,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
! Apply time filter to BLD (to remove diurnal cycle) to obtain "little h".
! "little h" is representative of the active mixing layer depth, used in B22 formula (eq 27).
if (GV%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
do j = js-1, je+1 ; do i = is-1, ie+1
do j=js-1,je+1 ; do i=is-1,ie+1
little_h(i,j) = rmean2ts(GV%Z_to_H*BLD(i,j), CS%MLD_filtered(i,j), &
CS%BLD_growing_Tfilt, CS%BLD_decaying_Tfilt, dt)
CS%MLD_filtered(i,j) = little_h(i,j)
Expand Down Expand Up @@ -920,21 +929,49 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
endif

! Calculate "big H", representative of the mixed layer depth, used in B22 formula (eq 27).
do j = js-1, je+1 ; do i = is-1, ie+1
do j=js-1,je+1 ; do i=is-1,ie+1
big_H(i,j) = rmean2ts(little_h(i,j), CS%MLD_filtered_slow(i,j), &
CS%MLD_growing_Tfilt, CS%MLD_decaying_Tfilt, dt)
CS%MLD_filtered_slow(i,j) = big_H(i,j)
enddo ; enddo

! Estimate w'u' at h-points
do j = js-1, je+1 ; do i = is-1, ie+1
w_star3 = max(0., -bflux(i,j)) * BLD(i,j) & ! (this line in Z3 T-3 ~> m3 s-3)
* ( ( US%Z_to_m * US%s_to_T )**3 ) ! [m3 T3 Z-3 s-3 ~> 1]
u_star3 = ( US%Z_to_m * US%s_to_T * U_star_2d(i,j) )**3 ! m3 s-3
wpup(i,j) = max( CS%min_wstar2, & ! The max() avoids division by zero later
( CS%mstar * u_star3 + CS%nstar * w_star3 )**two_thirds ) & ! (this line m2 s-2)
* ( US%m_to_L * GV%m_to_H * US%T_to_s**2 ) ! [L H s2 m-2 T-2 ~> 1 or kg m-3]
! We filter w'u' with the same time scales used for "little h"
! Estimate w'u' at h-points, with a floor to avoid division by zero later.
if (allocated(tv%SpV_avg) .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
do j=js-1,je+1 ; do i=is-1,ie+1
! This expression differs by a factor of 1. / (Rho_0 * SpV_avg) compared with the other
! expressions below, and it is invariant to the value of Rho_0 in non-Boussinesq mode.
wpup(i,j) = max((cuberoot( CS%mstar * U_star_2d(i,j)**3 + &
CS%nstar * max(0., -bflux(i,j)) * BLD(i,j) ))**2, CS%min_wstar2) * &
( US%Z_to_L * GV%RZ_to_H / tv%SpV_avg(i,j,1))
! The final line above converts from [Z2 T-2 ~> m2 s-2] to [L H T-2 ~> m2 s-2 or Pa].
! Some rescaling factors and the division by specific volume compensating for other
! factors that are in find_ustar_mech, and others effectively converting the wind
! stresses from [R L Z T-2 ~> Pa] to [L H T-2 ~> m2 s-2 or Pa]. The rescaling factors
! and density being applied to the buoyancy flux are not so neatly explained because
! fractional powers cancel out or combine with terms in the definitions of BLD and
! bflux (such as SpV_avg**-2/3 combining with other terms in bflux to give the thermal
! expansion coefficient) and because the specific volume does vary within the mixed layer.
enddo ; enddo
elseif (CS%answer_date < 20240201) then
Z3_T3_to_m3_s3 = (US%Z_to_m * US%s_to_T)**3
m2_s2_to_Z2_T2 = (US%m_to_Z * US%T_to_s)**2
do j=js-1,je+1 ; do i=is-1,ie+1
w_star3 = max(0., -bflux(i,j)) * BLD(i,j) ! In [Z3 T-3 ~> m3 s-3]
u_star3 = U_star_2d(i,j)**3 ! In [Z3 T-3 ~> m3 s-3]
wpup(i,j) = max(m2_s2_to_Z2_T2 * (Z3_T3_to_m3_s3 * ( CS%mstar * u_star3 + CS%nstar * w_star3 ) )**two_thirds, &
CS%min_wstar2) * &
( US%Z_to_L * US%Z_to_m * GV%m_to_H ) ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2]
enddo ; enddo
else
do j=js-1,je+1 ; do i=is-1,ie+1
w_star3 = max(0., -bflux(i,j)) * BLD(i,j) ! In [Z3 T-3 ~> m3 s-3]
wpup(i,j) = max( (cuberoot(CS%mstar * U_star_2d(i,j)**3 + CS%nstar * w_star3))**2, CS%min_wstar2 ) * &
( US%Z_to_L * US%Z_to_m * GV%m_to_H ) ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2]
enddo ; enddo
endif

! We filter w'u' with the same time scales used for "little h"
do j=js-1,je+1 ; do i=is-1,ie+1
wpup(i,j) = rmean2ts(wpup(i,j), CS%wpup_filtered(i,j), &
CS%BLD_growing_Tfilt, CS%BLD_decaying_Tfilt, dt)
CS%wpup_filtered(i,j) = wpup(i,j)
Expand Down Expand Up @@ -1483,7 +1520,7 @@ end subroutine mixedlayer_restrat_BML
!> Return the growth timescale for the submesoscale mixed layer eddies in [T ~> s]
real function growth_time(u_star, hBL, absf, h_neg, vonKar, Kv_rest, restrat_coef)
real, intent(in) :: u_star !< Surface friction velocity in thickness-based units [H T-1 ~> m s-1 or kg m-2 s-1]
real, intent(in) :: hBL !< Boundary layer thickness including at least a neglible
real, intent(in) :: hBL !< Boundary layer thickness including at least a negligible
!! value to keep it positive definite [H ~> m or kg m-2]
real, intent(in) :: absf !< Absolute value of the Coriolis parameter [T-1 ~> s-1]
real, intent(in) :: h_neg !< A tiny thickness that is usually lost in roundoff so can be
Expand Down Expand Up @@ -1537,6 +1574,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
real :: ustar_min_dflt ! The default value for RESTRAT_USTAR_MIN [Z T-1 ~> m s-1]
real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
! temperature variance [nondim]
integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags
! This include declares and sets the variable "version".
# include "version_variable.h"
integer :: i, j
Expand Down Expand Up @@ -1605,13 +1643,23 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
"BLD, when the latter is shallower than the running mean. A value of 0 "//&
"instantaneously sets the running mean to the current value filtered BLD.", &
units="s", default=0., scale=US%s_to_T)
call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
"This sets the default value for the various _ANSWER_DATE parameters.", &
default=99991231)
call get_param(param_file, mdl, "ML_RESTRAT_ANSWER_DATE", CS%answer_date, &
"The vintage of the order of arithmetic and expressions in the mixed layer "//&
"restrat calculations. Values below 20240201 recover the answers from the end "//&
"of 2023, while higher values use the new cuberoot function in the Bodner code "//&
"to avoid needing to undo dimensional rescaling.", &
default=default_answer_date, &
do_not_log=.not.(CS%use_Bodner.and.(GV%Boussinesq.or.GV%semi_Boussinesq)))
call get_param(param_file, mdl, "MIN_WSTAR2", CS%min_wstar2, &
"The minimum lower bound to apply to the vertical momentum flux, w'u', "//&
"in the Bodner et al., restratification parameterization. This avoids "//&
"a division-by-zero in the limit when u* and the buoyancy flux are zero. "//&
"The default is less than the molecular viscosity of water times the Coriolis "//&
"parameter a micron away from the equator.", &
units="m2 s-2", default=1.0e-24) ! This parameter stays in MKS units.
units="m2 s-2", default=1.0e-24, scale=US%m_to_Z**2*US%T_to_s**2)
call get_param(param_file, mdl, "TAIL_DH", CS%MLE_tail_dh, &
"Fraction by which to extend the mixed-layer restratification "//&
"depth used for a smoother stream function at the base of "//&
Expand Down
3 changes: 1 addition & 2 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -773,11 +773,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

I4dt = 0.25 / dt
I_slope_max2 = 1.0 / (CS%slope_max**2)
G_scale = GV%g_Earth * GV%H_to_Z

h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 ; hn_2 = 0.5*h_neglect
dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect**2
G_rho0 = GV%g_Earth / GV%Rho0
if (GV%Boussinesq) G_rho0 = GV%g_Earth / GV%Rho0
N2_floor = CS%N2_floor * US%Z_to_L**2

use_EOS = associated(tv%eqn_of_state)
Expand Down

0 comments on commit dd160e0

Please sign in to comment.