diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index fd5cfbeee3..4bac09b7cc 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -766,7 +766,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) & - .or. allocated(MEKE%Le)) then + .or. allocated(MEKE%Le)) then call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_Kh, G%Domain) call cpu_clock_end(CS%id_clock_pass) @@ -1570,7 +1570,7 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain) if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) & - .or. allocated(MEKE%Le)) & + .or. allocated(MEKE%Le)) & call do_group_pass(CS%pass_Kh, G%Domain) end function MEKE_init diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index fd8ba24afb..27fe0d335c 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -129,7 +129,7 @@ module MOM_lateral_mixing_coeffs real, allocatable :: khtr_struct(:,:,:) !< Vertical structure function used in tracer diffusivity [nondim] real, allocatable :: kdgl90_struct(:,:,:) !< Vertical structure function used in GL90 diffusivity [nondim] real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0. - real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 0.0 + real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 1.0 logical :: BS_use_sqg_struct !< If true, use sqg_stuct for backscatter vertical structure. @@ -556,12 +556,14 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] real, dimension(SZI_(G), SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1] - real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2] - real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m] + real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2] + real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m] + real :: f_subround ! The minimal resolved value of Coriolis parameter to prevent division by zero [T-1 ~> s-1] real, dimension(SZI_(G), SZJ_(G)) :: Le ! Eddy length scale [m] integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + f_subround = 1.0e-40 * US%s_to_T if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//& "Module must be initialized before it is used.") @@ -570,29 +572,33 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v,dzu=dzu, dzv=dzv, & dzSxN=dzSxN, dzSyN=dzSyN, halo=1) - - do j=js,je ; do i=is,ie - CS%sqg_struct(i,j,1) = 1.0 - enddo ; enddo - if (allocated(MEKE%Le)) then - do j=js,je ; do i=is,ie - Le(i,j) = MEKE%Le(i,j) - f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & - (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), 1.0e-8 * US%s_to_T) - enddo ; enddo + + if (CS%sqg_expo<=0.) then + CS%sqg_struct(:,:,:) = 1. else do j=js,je ; do i=is,ie - Le(i,j) = sqrt(G%areaT(i,j)) - f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & - (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), 1.0e-8 * US%s_to_T) + CS%sqg_struct(i,j,1) = 1.0 enddo ; enddo + if (allocated(MEKE%Le)) then + do j=js,je ; do i=is,ie + Le(i,j) = MEKE%Le(i,j) + f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), f_subround) + enddo ; enddo + else + do j=js,je ; do i=is,ie + Le(i,j) = sqrt(G%areaT(i,j)) + f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), f_subround) + enddo ; enddo + endif + do k=2,nz ; do j=js,je ; do i=is,ie + N2 = max(0.25 * ((N2_u(I-1,j,k) + N2_u(I,j,k)) + (N2_v(i,J-1,k) + N2_v(i,J,k))), 0.0) + dzc = 0.25 * ((dzu(I-1,j,k) + dzu(I,j,k)) + (dzv(i,J-1,k) + dzv(i,J,k))) + CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * & + exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j)))) + enddo ; enddo ; enddo endif - do k=2,nz ; do j=js,je ; do i=is,ie - N2 = max(0.25 * ((N2_u(I-1,j,k) + N2_u(I,j,k)) + (N2_v(i,J-1,k) + N2_v(i,J,k))), 0.0) - dzc = 0.25 * ((dzu(I-1,j,k) + dzu(I,j,k)) + (dzv(i,J-1,k) + dzv(i,J,k))) - CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * & - exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j)))) - enddo ; enddo ; enddo if (query_averaging_enabled(CS%diag)) then @@ -1393,7 +1399,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, & "Power to raise EBT vertical structure to when backscatter "// & "has vertical structure.", units="nondim", default=0.0) - call get_param(param_file, mdl, "BS_USE_SQG", CS%BS_use_sqg_struct, & + call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", CS%BS_use_sqg_struct, & "If true, the SQG vertical structure is used for backscatter "//& "on the condition that BS_EBT_power=0", & default=.false.) @@ -1480,26 +1486,26 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (CS%BS_EBT_power>0. .and. CS%BS_use_sqg_struct) then call MOM_error(FATAL, & - "calc_resoln_function: BS_EBT_POWER>0. & - and BS_USE_SQG=True cannot be set together") + "calc_resoln_function: BS_EBT_POWER>0. & + and BS_USE_SQG=True cannot be set together") endif if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) then call MOM_error(FATAL, & - "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT & - and KHTH_USE_SQG_STRUCT can be true") + "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT & + and KHTH_USE_SQG_STRUCT can be true") endif if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) then call MOM_error(FATAL, & - "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT & - and KHTR_USE_SQG_STRUCT can be true") + "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT & + and KHTR_USE_SQG_STRUCT can be true") endif if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) then call MOM_error(FATAL, & - "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT & - and KD_GL90_USE_SQG_STRUCT can be true") + "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT & + and KD_GL90_USE_SQG_STRUCT can be true") endif if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then @@ -1616,7 +1622,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, & 'Vertical structure of SQG mode', 'nondim') if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & - .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then + .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif @@ -1890,7 +1896,7 @@ subroutine VarMix_end(CS) type(VarMix_CS), intent(inout) :: CS if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & - .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct) + .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct) if (allocated(CS%sqg_struct)) deallocate(CS%sqg_struct) if (allocated(CS%BS_struct)) deallocate(CS%BS_struct) if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) deallocate(CS%khth_struct)