diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 52b55ad252..75f6d7e32c 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -49,6 +49,9 @@ module MOM_thickness_diffuse real :: kappa_smooth !< Vertical diffusivity used to interpolate more sensible values !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] logical :: thickness_diffuse !< If true, interfaces heights are diffused. + logical :: full_depth_khth_min !< If true, KHTH_MIN is enforced throughout the whole water column. + !! Otherwise, KHTH_MIN is only enforced at the surface. This parameter + !! is only available when KHTH_USE_EBT_STRUCT=True and KHTH_MIN>0. logical :: use_FGNV_streamfn !< If true, use the streamfunction formulation of !! Ferrari et al., 2010, which effectively emphasizes !! graver vertical modes by smoothing in the vertical. @@ -294,10 +297,18 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo ; enddo if (khth_use_ebt_struct) then - !$OMP do - do K=2,nz+1 ; do j=js,je ; do I=is-1,ie - KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) - enddo ; enddo ; enddo + if (CS%full_depth_khth_min) then + !$OMP do + do K=2,nz+1 ; do j=js,je ; do I=is-1,ie + KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + KH_u(I,j,K) = max(KH_u(I,j,K), CS%Khth_Min) + enddo ; enddo ; enddo + else + !$OMP do + do K=2,nz+1 ; do j=js,je ; do I=is-1,ie + KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + enddo ; enddo ; enddo + endif else !$OMP do do K=2,nz+1 ; do j=js,je ; do I=is-1,ie @@ -390,10 +401,18 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif if (khth_use_ebt_struct) then - !$OMP do - do K=2,nz+1 ; do J=js-1,je ; do i=is,ie - KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) - enddo ; enddo ; enddo + if (CS%full_depth_khth_min) then + !$OMP do + do K=2,nz+1 ; do J=js-1,je ; do i=is,ie + KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + KH_v(i,J,K) = max(KH_v(i,J,K), CS%Khth_Min) + enddo ; enddo ; enddo + else + !$OMP do + do K=2,nz+1 ; do J=js-1,je ; do i=is,ie + KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + enddo ; enddo ; enddo + endif else !$OMP do do K=2,nz+1 ; do J=js-1,je ; do i=is,ie @@ -2132,7 +2151,11 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) ! rotation [nondim]. 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. + logical :: khth_use_ebt_struct ! If true, uses the equivalent barotropic structure + ! as the vertical structure of thickness diffusivity. + ! Used to determine if FULL_DEPTH_KHTH_MIN should be + ! available. + integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. integer :: i, j CS%initialized = .true. @@ -2178,6 +2201,17 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) call get_param(param_file, mdl, "KHTH_MIN", CS%KHTH_Min, & "The minimum horizontal thickness diffusivity.", & default=0.0, units="m2 s-1", scale=US%m_to_L**2*US%T_to_s) + call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", khth_use_ebt_struct, & + "If true, uses the equivalent barotropic structure "//& + "as the vertical structure of thickness diffusivity.",& + default=.false., do_not_log=.true.) + if (khth_use_ebt_struct .and. CS%KHTH_Min>0.0) then + call get_param(param_file, mdl, "FULL_DEPTH_KHTH_MIN", CS%full_depth_khth_min, & + "If true, KHTH_MIN is enforced throughout the whole water column. "//& + "Otherwise, KHTH_MIN is only enforced at the surface. This parameter "//& + "is only available when KHTH_USE_EBT_STRUCT=True and KHTH_MIN>0.", & + default=.false.) + endif call get_param(param_file, mdl, "KHTH_MAX", CS%KHTH_Max, & "The maximum horizontal thickness diffusivity.", & default=0.0, units="m2 s-1", scale=US%m_to_L**2*US%T_to_s) diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 6f4e5d0f90..2e45803a60 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -52,8 +52,11 @@ module MOM_tracer_hor_diff real :: max_diff_CFL !< If positive, locally limit the along-isopycnal !! tracer diffusivity to keep the diffusive CFL !! locally at or below this value [nondim]. - logical :: KhTh_use_ebt_struct !< If true, uses the equivalent barotropic structure + logical :: KhTr_use_ebt_struct !< If true, uses the equivalent barotropic structure !! as the vertical structure of tracer diffusivity. + logical :: full_depth_khtr_min !< If true, KHTR_MIN is enforced throughout the whole water column. + !! Otherwise, KHTR_MIN is only enforced at the surface. This parameter + !! is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0. logical :: Diffuse_ML_interior !< If true, diffuse along isopycnals between !! the mixed layer and the interior. logical :: check_diffusive_CFL !< If true, automatically iterate the diffusion @@ -412,21 +415,40 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online enddo enddo enddo - if (CS%KhTh_use_ebt_struct) then - do K=2,nz+1 - do J=js-1,je - do i=is,ie - Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + if (CS%KhTr_use_ebt_struct) then + if (CS%full_depth_khtr_min) then + do K=2,nz+1 + do J=js-1,je + do i=is,ie + Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + Coef_y(i,J,K) = max(Coef_y(i,J,K), CS%KhTr_min) + enddo enddo enddo - enddo - do k=2,nz+1 - do j=js,je - do I=is-1,ie - Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + do k=2,nz+1 + do j=js,je + do I=is-1,ie + Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + Coef_x(I,j,K) = max(Coef_x(I,j,K), CS%KhTr_min) + enddo enddo enddo - enddo + else + do K=2,nz+1 + do J=js-1,je + do i=is,ie + Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + enddo + enddo + enddo + do k=2,nz+1 + do j=js,je + do I=is-1,ie + Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + enddo + enddo + enddo + endif endif do itt=1,num_itts @@ -468,7 +490,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online enddo enddo enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_ebt_struct) then do K=2,nz+1 do J=js-1,je do i=is,ie @@ -594,7 +616,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online do j=js,je ; do I=is-1,ie Kh_u(I,j,:) = G%mask2dCu(I,j)*Kh_u(I,j,1) enddo ; enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_ebt_struct) then do K=2,nz+1 do j=js,je do I=is-1,ie @@ -610,7 +632,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online do J=js-1,je ; do i=is,ie Kh_v(i,J,:) = G%mask2dCv(i,J)*Kh_v(i,J,1) enddo ; enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_ebt_struct) then do K=2,nz+1 do J=js-1,je do i=is,ie @@ -636,7 +658,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online (G%mask2dCv(i,J-1)+G%mask2dCv(i,J)) + 1.0e-37) Kh_h(i,j,:) = normalize*G%mask2dT(i,j)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & (Kh_v(i,J-1,1)+Kh_v(i,J,1))) - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_ebt_struct) then do K=2,nz+1 Kh_h(i,j,K) = normalize*G%mask2dT(i,j)*VarMix%ebt_struct(i,j,k-1)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & (Kh_v(i,J-1,1)+Kh_v(i,J,1))) @@ -1561,7 +1583,7 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic call get_param(param_file, mdl, "KHTR", CS%KhTr, & "The background along-isopycnal tracer diffusivity.", & units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) - call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, & + call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTr_use_ebt_struct, & "If true, uses the equivalent barotropic structure "//& "as the vertical structure of the tracer diffusivity.",& default=.false.) @@ -1573,6 +1595,13 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic call get_param(param_file, mdl, "KHTR_MIN", CS%KhTr_Min, & "The minimum along-isopycnal tracer diffusivity.", & units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) + if (CS%KhTr_use_ebt_struct .and. CS%KhTr_Min > 0.0) then + call get_param(param_file, mdl, "FULL_DEPTH_KHTR_MIN", CS%full_depth_khtr_min, & + "If true, KHTR_MIN is enforced throughout the whole water column. "//& + "Otherwise, KHTR_MIN is only enforced at the surface. This parameter "//& + "is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0.", & + default=.false.) + endif call get_param(param_file, mdl, "KHTR_MAX", CS%KhTr_Max, & "The maximum along-isopycnal tracer diffusivity.", & units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s)