Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

make MLD_003 computation available to mixedlayer_restrat_Bodner #327

Merged
merged 2 commits into from
Jan 7, 2025
Merged

make MLD_003 computation available to mixedlayer_restrat_Bodner #327

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1709248
make MLD_003 computation available to mixedlayer_restrat_Bodner
alperaltuntas Dec 20, 2024
ceb51f6
add missing params to Bodner MLD detection and rename new subroutine …
alperaltuntas Dec 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
175 changes: 116 additions & 59 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Open in desktop
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
Open in desktop
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
Loading