Skip to content

Commit

Permalink
make MLD_003 computation available to mixedlayer_restrat_Bodner (#327)
Browse files Browse the repository at this point in the history 
This PR introduces the calculate_mld_003 subroutine, previously a code black in mixedlayer_restrat_OM4, so as to make it available in mixedlayer_restrat_bodner as well. The newly added flag BODNER_USE_MLD_003 determines whether mld_003 is indeed used in mixedlayer_restrat_bodner.

* make MLD_003 computation available to mixedlayer_restrat_Bodner
* add missing params to Bodner MLD detection and rename new subroutine and runtime parameter
  • Loading branch information
@alperaltuntas
alperaltuntas authored Jan 7, 2025
1 parent acadc6e commit bd33548
Show file tree
Hide file tree
Showing 2 changed files with 122 additions and 62 deletions.
175 changes: 116 additions & 59 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ module MOM_mixed_layer_restrat
logical :: MLE_use_PBL_MLD !< If true, use the MLD provided by the PBL parameterization.
!! if false, MLE will calculate a MLD based on a density difference
!! based on the parameter MLE_DENSITY_DIFF.
logical :: Bodner_detect_MLD !< If true, detect the MLD based on given density difference criterion
!! (MLE_DENSITY_DIFF) in the Bodner et al. parameterization.
real :: vonKar !< The von Karman constant as used for mixed layer viscosity [nondim]
real :: MLE_MLD_decay_time !< Time-scale to use in a running-mean when MLD is retreating [T ~> s].
real :: MLE_MLD_decay_time2 !< Time-scale to use in a running-mean when filtered MLD is retreating [T ~> s].
Expand Down Expand Up @@ -244,14 +246,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix,
real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! Meridional restratification timescale [T ~> s], stored for diagnostics.
real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK ! Densities and density differences [R ~> kg m-3]
real, dimension(SZI_(G)) :: dK, dKm1 ! Depths of layer centers [H ~> m or kg m-2].
real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
! densities [R L2 T-2 ~> Pa].
real, dimension(SZI_(G)) :: covTS, & ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
varS ! SGS S variance in Stanley param; currently 0 [S2 ~> ppt2]
real :: aFac, bFac ! Nondimensional ratios [nondim]
real :: ddRho ! A density difference [R ~> kg m-3]
real :: hAtVel ! Thickness at the velocity points [H ~> m or kg m-2]
real :: zpa ! Fractional position within the mixed layer of the interface above a layer [nondim]
real :: zpb ! Fractional position within the mixed layer of the interface below a layer [nondim]
Expand Down Expand Up @@ -286,48 +283,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix,
call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1, H_T_units=.true.)

if (CS%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD.
!! 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
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, &
rhoSurf, tv%eqn_of_state, EOSdom)
else
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom)
endif
deltaRhoAtK(:) = 0.
MLD_fast(:,j) = 0.
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)
deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
if (CS%use_Stanley_ML) then
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, tv%varT(:,j,k), covTS, varS, &
deltaRhoAtK, tv%eqn_of_state, EOSdom)
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
deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
enddo
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
aFac = ( CS%MLE_density_diff - deltaRhoAtKm1(i) ) / ddRho
MLD_fast(i,j) = dK(i) * aFac + dKm1(i) * (1. - aFac)
endif
enddo ! i-loop
enddo ! k-loop
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
call detect_mld(h, tv, MLD_fast, G, GV, CS)
elseif (CS%MLE_use_PBL_MLD) then
do j=js-1,je+1 ; do i=is-1,ie+1
MLD_fast(i,j) = CS%MLE_MLD_stretch * h_MLD(i,j)
Expand Down Expand Up @@ -787,6 +743,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
real, dimension(SZI_(G),SZJ_(G)) :: &
little_h, & ! "Little h" representing active mixing layer depth [H ~> m or kg m-2]
big_H, & ! "Big H" representing the mixed layer depth [H ~> m or kg m-2]
mld, & ! The mixed layer depth returned by detect_mld [H ~> m or kg m-2]
htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
buoy_av, & ! g_Rho0 times the average mixed layer density or G_Earth
! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
Expand Down Expand Up @@ -853,14 +810,19 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d

if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"An equation of state must be used with this module.")
if (.not.CS%MLE_use_PBL_MLD) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"To use the Bodner et al., 2023, MLE parameterization, MLE_USE_PBL_MLD must be True.")
if (CS%MLE_density_diff > 0.) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"MLE_density_diff is +ve and should not be in mixedlayer_restrat_Bodner.")
if (.not.associated(bflux)) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"Surface buoyancy flux was not associated.")
if (CS%MLE_use_PBL_MLD) then
if (CS%MLE_density_diff > 0.) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"MLE_density_diff is +ve and should not be in mixedlayer_restrat_Bodner.")
if (.not.associated(bflux)) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"Surface buoyancy flux was not associated.")
else
if (.not.CS%Bodner_detect_MLD) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
"To use the Bodner et al., 2023, MLE parameterization, either MLE_USE_PBL_MLD or "// &
"Bodner_detect_MLD must be True.")
endif

call pass_var(bflux, G%domain, halo=1)
if (associated(bflux)) &
call pass_var(bflux, G%domain, halo=1)

! Extract the friction velocity from the forcing type.
call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)
Expand All @@ -887,9 +849,19 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
enddo ; enddo

! Calculate "big H", representative of the mixed layer depth, used in B22 formula (eq 27).
if (CS%Bodner_detect_MLD) then
call detect_mld(h, tv, MLD, G, GV, CS)
do j=js-1,je+1 ; do i=is-1,ie+1
big_H(i,j) = rmean2ts(MLD(i,j), CS%MLD_filtered_slow(i,j), &
CS%MLD_growing_Tfilt, CS%MLD_decaying_Tfilt, dt)
enddo ; enddo
else
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)
enddo ; enddo
endif
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

Expand Down Expand Up @@ -1485,6 +1457,76 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)

end subroutine mixedlayer_restrat_BML

!> Detects the mixed layer depth using a density difference criterion (MLE_DENSITY_DIFF)
subroutine detect_mld(h, tv, MLD_fast, G, GV, CS)
type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: MLD_fast !< detected mixed layer depth [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables structure

! Local variables
real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
! densities [R L2 T-2 ~> Pa].
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK ! Densities and density differences [R ~> kg m-3]
real, dimension(SZI_(G)) :: dK, dKm1 ! Depths of layer centers [H ~> m or kg m-2].
real :: ddRho ! A density difference [R ~> kg m-3]
real :: aFac ! A nondimensional ratio [nondim]
real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [S2 ~> ppt2]
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, j, k, is, ie, js, je, nz

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

covTS(:) = 0.0 ! Might be in tv% in the future. Not implemented for the time being.
varS(:) = 0.0 ! Ditto.

!! 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
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, &
rhoSurf, tv%eqn_of_state, EOSdom)
else
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom)
endif
deltaRhoAtK(:) = 0.
MLD_fast(:,j) = 0.
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)
deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
if (CS%use_Stanley_ML) then
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, tv%varT(:,j,k), covTS, varS, &
deltaRhoAtK, tv%eqn_of_state, EOSdom)
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
deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
enddo
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
aFac = ( CS%MLE_density_diff - deltaRhoAtKm1(i) ) / ddRho
MLD_fast(i,j) = dK(i) * aFac + dKm1(i) * (1. - aFac)
endif
enddo ! i-loop
enddo ! k-loop
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
end subroutine detect_mld

! NOTE: This function appears to change answers on some platforms, so it is
! currently unused in the model, but we intend to introduce it in the future.

Expand Down Expand Up @@ -1643,9 +1685,15 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
"If true, the MLE parameterization will use the mixed-layer "//&
"depth provided by the active PBL parameterization. If false, "//&
"MLE will estimate a MLD based on a density difference with the "//&
"surface using the parameter MLE_DENSITY_DIFF.", default=.false.)
if (.not.CS%MLE_use_PBL_MLD) call MOM_error(FATAL, "mixedlayer_restrat_init: "// &
"To use MLE%USE_BODNER23=True then MLE_USE_PBL_MLD must be True.")
"surface using the parameter MLE_DENSITY_DIFF, unless "//&
"BODNER_DETECT_MLD is true.", default=.false.)
call get_param(param_file, mdl, "BODNER_DETECT_MLD", CS%Bodner_detect_MLD, &
"If true, the Bodner parameterization will use the mixed-layer depth "//&
"detected via the density difference criterion MLE_DENSITY_DIFF.", default=.false.)
if (.not.(CS%MLE_use_PBL_MLD.or.CS%Bodner_detect_MLD)) call MOM_error(FATAL, "mixedlayer_restrat_init: "// &
"To use MLE%USE_BODNER23=True then MLE_USE_PBL_MLD or BODNER_DETECT_MLD must be true.")
if (CS%MLE_use_PBL_MLD.and.CS%Bodner_detect_MLD) call MOM_error(FATAL, "mixedlayer_restrat_init: "// &
"MLE_USE_PBL_MLD and BODNER_DETECT_MLD cannot both be true.")
else
call closeParameterBlock(param_file) ! The remaining parameters do not have MLE% prepended
endif
Expand Down Expand Up @@ -1731,6 +1779,15 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
"restratification module. This can be tiny, but if this is greater than 0, "//&
"it will prevent divisions by zero when f and KV_RESTRAT are zero.", &
units="m s-1", default=US%Z_to_m*US%s_to_T*ustar_min_dflt, scale=GV%m_to_H*US%T_to_s)
elseif (CS%Bodner_detect_MLD) then
call get_param(param_file, mdl, "MLE_DENSITY_DIFF", CS%MLE_density_diff, &
"Density difference used to detect the mixed-layer "//&
"depth used for the mixed-layer eddy parameterization "//&
"by Fox-Kemper et al. (2010)", units="kg/m3", default=0.03, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "MLE_MLD_STRETCH", CS%MLE_MLD_stretch, &
"A scaling coefficient for stretching/shrinking the MLD "//&
"used in the MLE scheme. This simply multiplies MLD wherever used.",&
units="nondim", default=1.0)
endif

CS%diag => diag
Expand Down
9 changes: 6 additions & 3 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ module MOM_set_visc
use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_specific_vol_derivs
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_file_parser, only : openParameterBlock, closeParameterBlock
use MOM_forcing_type, only : forcing, mech_forcing, find_ustar
use MOM_grid, only : ocean_grid_type
use MOM_hor_index, only : hor_index_type
Expand Down Expand Up @@ -2783,8 +2784,12 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C
default=.false., do_not_log=.true.)
call get_param(param_file, mdl, "USE_IDEAL_AGE_TRACER", use_ideal_age, &
default=.false., do_not_log=.true.)
call openParameterBlock(param_file, 'MLE', do_not_log=.true.)
call get_param(param_file, mdl, "USE_BODNER23", MLE_use_Bodner, &
default=.false., do_not_log=.true.)
call closeParameterBlock(param_file)

if (MLE_use_PBL_MLD) then
if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
endif
if ((hfreeze >= 0.0) .or. MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. &
Expand All @@ -2804,8 +2809,6 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C
endif

! visc%sfc_buoy_flx is used to communicate the state of the (e)PBL or KPP to the rest of the model
call get_param(param_file, mdl, "MLE%USE_BODNER23", MLE_use_Bodner, &
default=.false., do_not_log=.true.)
if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
call safe_alloc_ptr(visc%sfc_buoy_flx, isd, ied, jsd, jed)
call register_restart_field(visc%sfc_buoy_flx, "SFC_BFLX", .false., restart_CS, &
Expand Down

0 comments on commit bd33548

Please sign in to comment.