Skip to content

Reversion of MOM_mixed_layer_restrat growth_time #358

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

Merged
Hallberg-NOAA merged 1 commit into from
Apr 29, 2023
Merged

Reversion of MOM_mixed_layer_restrat growth_time #358

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
89 changes: 81 additions & 8 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 @@ -159,6 +159,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
real :: h_vel ! htot interpolated onto velocity points [Z ~> m] (not H).
real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
real :: u_star ! surface friction velocity, interpolated to velocity points [Z T-1 ~> m s-1].
real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
real :: timescale ! mixing growth timescale [T ~> s]
real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
real :: h_neglect ! tiny thickness usually lost in roundoff so can be neglected [H ~> m or kg m-2]
Expand Down Expand Up @@ -194,6 +195,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
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 :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1]
real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times
! pi squared [nondim]
logical :: line_is_empty, keep_going, res_upscale
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
Expand All @@ -205,6 +208,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
covTS(:) = 0.0 !!Functionality not implemented yet; in future, should be passed in tv
varS(:) = 0.0

vonKar_x_pi2 = CS%vonKar * 9.8696

if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
"An equation of state must be used with this module.")
if (.not. allocated(VarMix%Rd_dx_h) .and. CS%front_length > 0.) &
Expand Down Expand Up @@ -316,7 +321,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var

p0(:) = 0.0
EOSdom(:) = EOS_domain(G%HI, halo=1)
!$OMP parallel default(shared) private(rho_ml,h_vel,u_star,absf,timescale, &
!$OMP parallel default(shared) private(rho_ml,h_vel,u_star,absf,mom_mixrate,timescale, &
!$OMP line_is_empty, keep_going,res_scaling_fac, &
!$OMP a,IhTot,b,Ihtot_slow,zpb,hAtVel,zpa,dh) &
!$OMP firstprivate(uDml,vDml,uDml_slow,vDml_slow)
Expand Down Expand Up @@ -379,21 +384,40 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
!$OMP do
do j=js,je ; do I=is-1,ie
u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j)))

absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
! If needed, res_scaling_fac = min( ds, L_d ) / l_f
if (res_upscale) res_scaling_fac = &
( sqrt( 0.5 * ( G%dxCu(I,j)**2 + G%dyCu(I,j)**2 ) ) * I_LFront ) &
* min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i+1,j) ) )

! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
! momentum mixing rate: pi^2*visc/h_ml^2
h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect) * GV%H_to_Z
timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

! NOTE: growth_time changes answers on some systems, see below.
! timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

mom_mixrate = vonKar_x_pi2*u_star**2 / &
(absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
timescale = timescale * CS%ml_restrat_coef

if (res_upscale) timescale = timescale * res_scaling_fac
uDml(I) = timescale * G%OBCmaskCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * &
(Rml_av_fast(i+1,j)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H)

! As above but using the slow filtered MLD
h_vel = 0.5*((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect) * GV%H_to_Z
timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef2)

! NOTE: growth_time changes answers on some systems, see below.
! timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef2)

mom_mixrate = vonKar_x_pi2*u_star**2 / &
(absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
timescale = timescale * CS%ml_restrat_coef2

if (res_upscale) timescale = timescale * res_scaling_fac
uDml_slow(I) = timescale * G%OBCmaskCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * &
(Rml_av_slow(i+1,j)-Rml_av_slow(i,j)) * (h_vel**2 * GV%Z_to_H)
Expand Down Expand Up @@ -447,21 +471,40 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
!$OMP do
do J=js-1,je ; do i=is,ie
u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)))

absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
! If needed, res_scaling_fac = min( ds, L_d ) / l_f
if (res_upscale) res_scaling_fac = &
( sqrt( 0.5 * ( (G%dxCv(i,J))**2 + (G%dyCv(i,J))**2 ) ) * I_LFront ) &
* min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i,j+1) ) )

! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
! momentum mixing rate: pi^2*visc/h_ml^2
h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect) * GV%H_to_Z
timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

! NOTE: growth_time changes answers on some systems, see below.
! timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

mom_mixrate = vonKar_x_pi2*u_star**2 / &
(absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
timescale = timescale * CS%ml_restrat_coef

if (res_upscale) timescale = timescale * res_scaling_fac
vDml(i) = timescale * G%OBCmaskCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * &
(Rml_av_fast(i,j+1)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H)

! As above but using the slow filtered MLD
h_vel = 0.5*((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect) * GV%H_to_Z
timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef2)

! NOTE: growth_time changes answers on some systems, see below.
! timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef2)

mom_mixrate = vonKar_x_pi2*u_star**2 / &
(absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
timescale = timescale * CS%ml_restrat_coef2

if (res_upscale) timescale = timescale * res_scaling_fac
vDml_slow(i) = timescale * G%OBCmaskCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * &
(Rml_av_slow(i,j+1)-Rml_av_slow(i,j)) * (h_vel**2 * GV%Z_to_H)
Expand Down Expand Up @@ -608,6 +651,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
real :: h_vel ! htot interpolated onto velocity points [Z ~> m]. (The units are not H.)
real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
real :: u_star ! surface friction velocity, interpolated to velocity points [Z T-1 ~> m s-1].
real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times
! pi squared [nondim]
real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
real :: timescale ! mixing growth timescale [T ~> s]
real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
real :: h_neglect ! tiny thickness usually lost in roundoff and can be neglected [H ~> m or kg m-2]
Expand Down Expand Up @@ -642,6 +688,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
uDml(:) = 0.0 ; vDml(:) = 0.0
I4dt = 0.25 / dt
g_Rho0 = GV%g_Earth / GV%Rho0
vonKar_x_pi2 = CS%vonKar * 9.8696
use_EOS = associated(tv%eqn_of_state)
h_neglect = GV%H_subroundoff
dz_neglect = GV%H_subroundoff*GV%H_to_Z
Expand All @@ -657,7 +704,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)

p0(:) = 0.0
EOSdom(:) = EOS_domain(G%HI, halo=1)
!$OMP parallel default(shared) private(Rho0,h_vel,u_star,absf,timescale, &
!$OMP parallel default(shared) private(Rho0,h_vel,u_star,absf,mom_mixrate,timescale, &
!$OMP I2htot,z_topx2,hx2,a) &
!$OMP firstprivate(uDml,vDml)
!$OMP do
Expand Down Expand Up @@ -689,8 +736,19 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * GV%H_to_Z

u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j)))

absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

! NOTE: growth_time changes answers on some systems, see below.
! timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
! momentum mixing rate: pi^2*visc/h_ml^2
mom_mixrate = vonKar_x_pi2*u_star**2 / &
(absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)

timescale = timescale * CS%ml_restrat_coef
! timescale = timescale*(2?)*(L_def/L_MLI) * min(EKE/MKE,1.0 + (G%dyCv(i,j)/L_def)**2)

uDml(I) = timescale * G%OBCmaskCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * &
Expand Down Expand Up @@ -729,8 +787,19 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * GV%H_to_Z

u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)))

absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

! NOTE: growth_time changes answers on some systems, see below.
! timescale = growth_time(u_star, h_vel, absf, dz_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)

! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
! momentum mixing rate: pi^2*visc/h_ml^2
mom_mixrate = vonKar_x_pi2*u_star**2 / &
(absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)

timescale = timescale * CS%ml_restrat_coef
! timescale = timescale*(2?)*(L_def/L_MLI) * min(EKE/MKE,1.0 + (G%dyCv(i,j)/L_def)**2)

vDml(i) = timescale * G%OBCmaskCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * &
Expand Down Expand Up @@ -799,6 +868,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)

end subroutine mixedlayer_restrat_BML

! 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.

!> 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 [Z T-1 ~> m s-1]
Expand Down Expand Up @@ -945,6 +1017,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
"used in the MLE scheme. This simply multiplies MLD wherever used.",&
units="nondim", default=1.0)
endif

call get_param(param_file, mdl, "KV_RESTRAT", CS%Kv_restrat, &
"A small viscosity that sets a floor on the momentum mixing rate during "//&
"restratification. If this is positive, it will prevent some possible "//&
Expand Down