diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 1119b5660f..880af32fa9 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -45,6 +45,15 @@ module MOM_lateral_mixing_coeffs !! for layer thicknesses. In addition, masking at coastlines was not !! used which introduced potential restart issues. This flag will be !! deprecated in a future release. + logical :: calculate_cg1 !< If true, calls wave_speed() to calculate the first + !! baroclinic wave speed and populate CS%cg1. + !! This parameter is set depending on other parameters. + logical :: calculate_Rd_dx !< If true, calculates Rd/dx and populate CS%Rd_dx_h. + !! This parameter is set depending on other parameters. + logical :: calculate_res_fns !< If true, calculate all the resolution factors. + !! This parameter is set depending on other parameters. + logical :: calculate_Eady_growth_rate !< If true, calculate all the Eady growth rate. + !! This parameter is set depending on other parameters. real, dimension(:,:), pointer :: & SN_u => NULL(), & !< S*N at u-points (s^-1) SN_v => NULL(), & !< S*N at v-points (s^-1) @@ -139,10 +148,48 @@ subroutine calc_resoln_function(h, tv, G, GV, CS) if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "calc_resoln_function:"// & "Module must be initialized before it is used.") - if (.not. (CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. & - CS%Resoln_scaled_KhTr)) return - if (.not. ASSOCIATED(CS%cg1)) call MOM_error(FATAL, & - "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.") + if (CS%calculate_cg1) then + if (.not. ASSOCIATED(CS%cg1)) call MOM_error(FATAL, & + "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.") + if (CS%khth_use_ebt_struct) then + if (.not. ASSOCIATED(CS%ebt_struct)) call MOM_error(FATAL, & + "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.") + if (CS%Resoln_use_ebt) then + ! Both resolution fn and vertical structure are using EBT + call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp, modal_structure=CS%ebt_struct) + else + ! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode + call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp, modal_structure=CS%ebt_struct, use_ebt_mode=.true.) + call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp) + endif + call pass_var(CS%ebt_struct, G%Domain) + else + call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp) + endif + + call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain) + call do_group_pass(CS%pass_cg1, G%Domain) + endif + + ! Calculate and store the ratio between deformation radius and grid-spacing + ! at h-points (non-dimensional). + if (CS%calculate_rd_dx) then + if (.not. ASSOCIATED(CS%Rd_dx_h)) call MOM_error(FATAL, & + "calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.") +!$OMP parallel default(none) shared(is,ie,js,je,CS) +!$OMP do + do j=js-1,je+1 ; do i=is-1,ie+1 + CS%Rd_dx_h(i,j) = CS%cg1(i,j) / & + (sqrt(CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j))) + enddo ; enddo +!$OMP end parallel + if (query_averaging_enabled(CS%diag)) then + if (CS%id_Rd_dx > 0) call post_data(CS%id_Rd_dx, CS%Rd_dx_h, CS%diag) + endif + endif + + if (.not. CS%calculate_res_fns) return + if (.not. ASSOCIATED(CS%Res_fn_h)) call MOM_error(FATAL, & "calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.") if (.not. ASSOCIATED(CS%Res_fn_q)) call MOM_error(FATAL, & @@ -168,23 +215,6 @@ subroutine calc_resoln_function(h, tv, G, GV, CS) if (.not. ASSOCIATED(CS%beta_dx2_v)) call MOM_error(FATAL, & "calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.") - if (CS%khth_use_ebt_struct) then - if (CS%Resoln_use_ebt) then - ! Both resolution fn and vertical structure are using EBT - call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp, modal_structure=CS%ebt_struct) - else - ! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode - call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp, modal_structure=CS%ebt_struct, use_ebt_mode=.true.) - call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp) - endif - call pass_var(CS%ebt_struct, G%Domain) - else - call wave_speed(h, tv, G, GV, CS%cg1, CS%wave_speed_CSp) - endif - - call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain) - call do_group_pass(CS%pass_cg1, G%Domain) - ! Do this calculation on the extent used in MOM_hor_visc.F90, and ! MOM_tracer.F90 so that no halo update is needed. @@ -336,19 +366,10 @@ subroutine calc_resoln_function(h, tv, G, GV, CS) enddo ; enddo endif endif - - ! Calculate and store the ratio between deformation radius and grid-spacing - ! at h-points (non-dimensional). -!$OMP do - do j=js-1,je+1 ; do i=is-1,ie+1 - CS%Rd_dx_h(i,j) = CS%cg1(i,j) / & - (sqrt(CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j))) - enddo ; enddo !$OMP end parallel if (query_averaging_enabled(CS%diag)) then if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag) - if (CS%id_Rd_dx > 0) call post_data(CS%id_Rd_dx, CS%Rd_dx_h, CS%diag) endif end subroutine calc_resoln_function @@ -371,8 +392,8 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, CS) if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//& "Module must be initialized before it is used.") - call find_eta(h, tv, GV%g_Earth, G, GV, e, halo_size=2) - if (CS%use_variable_mixing) then + if (CS%calculate_Eady_growth_rate) then + call find_eta(h, tv, GV%g_Earth, G, GV, e, halo_size=2) if (CS%use_stored_slopes) then call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, & CS%slope_x, CS%slope_y, N2_u, N2_v, 1) @@ -389,10 +410,8 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, CS) if (CS%id_SN_v > 0) call post_data(CS%id_SN_v, CS%SN_v, CS%diag) if (CS%id_L2u > 0) call post_data(CS%id_L2u, CS%L2u, CS%diag) if (CS%id_L2v > 0) call post_data(CS%id_L2v, CS%L2v, CS%diag) - if (CS%use_stored_slopes) then - if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag) - if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag) - endif + if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag) + if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag) endif end subroutine calc_slope_functions @@ -407,7 +426,7 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in) :: N2_u !< Brunt-Vaisala frequency at u-points (1/s2) real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: slope_y !< Meridional isoneutral slope real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: N2_v !< Brunt-Vaisala frequency at v-points (1/s2) - type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients + type(VarMix_CS), pointer :: CS !< Variable mixing coefficients ! Local variables real :: E_x(SZIB_(G), SZJ_(G)) ! X-slope of interface at u points (for diagnostics) real :: E_y(SZI_(G), SZJB_(G)) ! Y-slope of interface at u points (for diagnostics) @@ -425,9 +444,9 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) real :: S2_u(SZIB_(G), SZJ_(G)) real :: S2_v(SZI_(G), SZJB_(G)) - if (LOC(CS)==0) call MOM_error(FATAL, "calc_slope_function:"// & + if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "calc_slope_function:"// & "Module must be initialized before it is used.") - if (.not. CS%use_variable_mixing) return + if (.not. CS%calculate_Eady_growth_rate) return if (.not. ASSOCIATED(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_u is not associated with use_variable_mixing.") if (.not. ASSOCIATED(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:"// & @@ -452,17 +471,6 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) ! calculate the first-mode gravity wave speed and then blend the equatorial ! and midlatitude deformation radii, using calc_resoln_function as a template. - ! Set the length scale at u-points. -!$OMP do - do j=js,je ; do I=is-1,ie - CS%L2u(I,j) = CS%Visbeck_L_scale**2 - enddo ; enddo - ! Set length scale at v-points -!$OMP do - do J=js-1,je ; do i=is,ie - CS%L2v(i,J) = CS%Visbeck_L_scale**2 - enddo ; enddo - !$OMP do do j = js,je do I=is-1,ie @@ -564,9 +572,9 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) endif if (CS%debug) then - call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) - call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", N2_u, N2_v, G%HI) - call uvchksum("calc_Visbeck_coeffs SN_[uv]", CS%SN_u, CS%SN_v, G%HI) + call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) + call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", N2_u, N2_v, G%HI) + call uvchksum("calc_Visbeck_coeffs SN_[uv]", CS%SN_u, CS%SN_v, G%HI) endif end subroutine calc_Visbeck_coeffs @@ -600,7 +608,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "calc_slope_function:"// & "Module must be initialized before it is used.") - if (.not. CS%use_variable_mixing) return + if (.not. CS%calculate_Eady_growth_rate) return if (.not. ASSOCIATED(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_u is not associated with use_variable_mixing.") if (.not. ASSOCIATED(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:"// & @@ -625,16 +633,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) ! calculate the first-mode gravity wave speed and then blend the equatorial ! and midlatitude deformation radii, using calc_resoln_function as a template. - ! Set the length scale at u-points. -!$OMP do - do j=js,je ; do I=is-1,ie - CS%L2u(I,j) = CS%Visbeck_L_scale**2 - enddo ; enddo - ! Set length scale at v-points -!$OMP do - do J=js-1,je ; do i=is,ie - CS%L2v(i,J) = CS%Visbeck_L_scale**2 - enddo ; enddo !$OMP do do k=nz,CS%VarMix_Ktop,-1 @@ -730,12 +728,11 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) type(VarMix_CS), pointer :: CS !< Variable mixing coefficients ! Local variables real :: KhTr_Slope_Cff, KhTh_Slope_Cff, oneOrTwo, N2_filter_depth + real :: KhTr_passivity_coeff real, parameter :: absurdly_small_freq2 = 1e-34 ! A miniscule frequency ! squared that is used to avoid division by 0, in s-2. This ! value is roughly (pi / (the age of the universe) )^2. - logical :: use_variable_mixing, Gill_equatorial_Ld, use_stored_slopes - logical :: Resoln_scaled_Kh, Resoln_scaled_KhTh, Resoln_scaled_KhTr - logical :: Resoln_use_ebt, khth_use_ebt_struct + logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name. @@ -752,33 +749,39 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) return endif + allocate(CS) + in_use = .false. ! Set to true to avoid deallocating + CS%diag => diag ! Diagnostics pointer + CS%calculate_cg1 = .false. + CS%calculate_Rd_dx = .false. + CS%calculate_res_fns = .false. + CS%calculate_Eady_growth_rate = .false. + ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") - ! This first set of parameters are read into local variables first, in case - ! the control structure should not be allocated. - call get_param(param_file, mdl, "USE_VARIABLE_MIXING", use_variable_mixing,& + call get_param(param_file, mdl, "USE_VARIABLE_MIXING", CS%use_variable_mixing,& "If true, the variable mixing code will be called. This \n"//& "allows diagnostics to be created even if the scheme is \n"//& "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, \n"//& "this is set to true regardless of what is in the \n"//& "parameter file.", default=.false.) - call get_param(param_file, mdl, "RESOLN_SCALED_KH", Resoln_scaled_Kh, & + call get_param(param_file, mdl, "RESOLN_SCALED_KH", CS%Resoln_scaled_Kh, & "If true, the Laplacian lateral viscosity is scaled away \n"//& "when the first baroclinic deformation radius is well \n"//& "resolved.", default=.false.) - call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", Resoln_scaled_KhTh, & + call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", CS%Resoln_scaled_KhTh, & "If true, the interface depth diffusivity is scaled away \n"//& "when the first baroclinic deformation radius is well \n"//& "resolved.", default=.false.) - call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", Resoln_scaled_KhTr, & + call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", CS%Resoln_scaled_KhTr, & "If true, the epipycnal tracer diffusivity is scaled \n"//& "away when the first baroclinic deformation radius is \n"//& "well resolved.", default=.false.) - call get_param(param_file, mdl, "RESOLN_USE_EBT", Resoln_use_ebt, & + call get_param(param_file, mdl, "RESOLN_USE_EBT", CS%Resoln_use_ebt, & "If true, uses the equivalent barotropic wave speed instead\n"//& "of first baroclinic wave for calculating the resolution fn.",& default=.false.) - call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", khth_use_ebt_struct, & + call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, & "If true, uses the equivalent barotropic structure\n"//& "as the vertical structure of thickness diffusivity.",& default=.false.) @@ -790,38 +793,34 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "The nondimensional coefficient in the Visbeck formula \n"//& "for the epipycnal tracer diffusivity", units="nondim", & default=0.0) - call get_param(param_file, mdl, "USE_STORED_SLOPES", use_stored_slopes,& + call get_param(param_file, mdl, "USE_STORED_SLOPES", CS%use_stored_slopes,& "If true, the isopycnal slopes are calculated once and\n"//& "stored for re-use. This uses more memory but avoids calling\n"//& "the equation of state more times than should be necessary.", & default=.false.) - if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) use_variable_mixing = .true. - - if (use_variable_mixing .or. Resoln_scaled_Kh .or. Resoln_scaled_KhTh .or. & - Resoln_scaled_KhTr .or. use_stored_slopes .or. khth_use_ebt_struct) then - allocate(CS) - CS%diag => diag ! Diagnostics pointer - CS%Resoln_scaled_Kh = Resoln_scaled_Kh - CS%Resoln_scaled_KhTh = Resoln_scaled_KhTh - CS%Resoln_scaled_KhTr = Resoln_scaled_KhTr - CS%Resoln_use_ebt = Resoln_use_ebt - CS%khth_use_ebt_struct = khth_use_ebt_struct - CS%use_variable_mixing = use_variable_mixing - CS%use_stored_slopes = use_stored_slopes - call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.) - else - return - endif - if (Resoln_use_ebt .or. khth_use_ebt_struct) then + call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, & + default=.false., do_not_log=.true.) + CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn + call get_param(param_file, mdl, "USE_MEKE", use_MEKE, & + default=.false., do_not_log=.true.) + CS%calculate_Eady_growth_rate = CS%calculate_Eady_growth_rate .or. use_MEKE + call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", KhTr_passivity_coeff, & + default=0., do_not_log=.true.) + CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (KhTr_passivity_coeff>0.) + + call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.) + + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct) then + in_use = .true. call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, & "The depth below which N2 is monotonized to avoid stratification\n"//& "artifacts from altering the equivalent barotropic mode structure.",& units='m', default=2000.) - endif - if (khth_use_ebt_struct) then allocate(CS%ebt_struct(isd:ied,jsd:jed,G%ke)) ; CS%ebt_struct(:,:,:) = 0.0 endif - if (use_variable_mixing) then + + if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then + CS%calculate_Eady_growth_rate = .true. call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", CS%Visbeck_S_max, & "If non-zero, is an upper bound on slopes used in the\n"// & "Visbeck formula for diffusivity. This does not affect the\n"// & @@ -829,8 +828,8 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) units="nondim", default=0.0) endif -! Allocate CS and memory if (CS%use_stored_slopes) then + in_use = .true. allocate(CS%slope_x(IsdB:IedB,jsd:jed,G%ke+1)) ; CS%slope_x(:,:,:) = 0.0 allocate(CS%slope_y(isd:ied,JsdB:JedB,G%ke+1)) ; CS%slope_y(:,:,:) = 0.0 call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & @@ -839,67 +838,61 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) default=1.0e-6) endif - if (CS%use_variable_mixing) then + if (CS%calculate_Eady_growth_rate) then + in_use = .true. allocate(CS%SN_u(IsdB:IedB,jsd:jed)) ; CS%SN_u(:,:) = 0.0 allocate(CS%SN_v(isd:ied,JsdB:JedB)) ; CS%SN_v(:,:) = 0.0 - allocate(CS%L2u(IsdB:IedB,jsd:jed)) ; CS%L2u(:,:) = 0.0 - allocate(CS%L2v(isd:ied,JsdB:JedB)) ; CS%L2v(:,:) = 0.0 - call MOM_mesg("VarMix_init: memory allocated for use_variable_mixing", 5) - - ! More run-time parameters + CS%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, Time, & + 'Inverse eddy time-scale, S*N, at u-points', 's^-1') + CS%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, Time, & + 'Inverse eddy time-scale, S*N, at v-points', 's^-1') call get_param(param_file, mdl, "VARMIX_KTOP", CS%VarMix_Ktop, & "The layer number at which to start vertical integration \n"//& "of S*N for purposes of finding the Eady growth rate.", & units="nondim", default=2) + endif + + if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then + in_use = .true. call get_param(param_file, mdl, "VISBECK_L_SCALE", CS%Visbeck_L_scale, & "The fixed length scale in the Visbeck formula.", units="m", & default=0.0) + allocate(CS%L2u(IsdB:IedB,jsd:jed)) ; CS%L2u(:,:) = CS%Visbeck_L_scale**2 + allocate(CS%L2v(isd:ied,JsdB:JedB)) ; CS%L2v(:,:) = CS%Visbeck_L_scale**2 - ! Register fields for output from this module. - CS%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, Time, & - 'Inverse eddy time-scale, S*N, at u-points', 's^-1') - CS%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, Time, & - 'Inverse eddy time-scale, S*N, at v-points', 's^-1') CS%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, Time, & 'Length scale squared for mixing coefficient, at u-points', 'm^2') CS%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, Time, & 'Length scale squared for mixing coefficient, at v-points', 'm^2') + endif - if (CS%use_stored_slopes) then - CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, & + if (CS%use_stored_slopes) then + CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, & 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', 's^-2') - CS%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, Time, & + CS%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, Time, & 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', 's^-2') - CS%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, Time, & + CS%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, Time, & 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', 's^-2') - CS%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, Time, & + CS%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, Time, & 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 's^-2') - endif endif - if (CS%Resoln_scaled_Kh .or. Resoln_scaled_KhTh .or. Resoln_scaled_KhTr) then - call wave_speed_init(CS%wave_speed_CSp, use_ebt_mode=Resoln_use_ebt, mono_N2_depth=N2_filter_depth) - - ! Allocate and initialize various arrays. + if (CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. CS%Resoln_scaled_KhTr) then + CS%calculate_Rd_dx = .true. + CS%calculate_res_fns = .true. allocate(CS%Res_fn_h(isd:ied,jsd:jed)) ; CS%Res_fn_h(:,:) = 0.0 allocate(CS%Res_fn_q(IsdB:IedB,JsdB:JedB)) ; CS%Res_fn_q(:,:) = 0.0 allocate(CS%Res_fn_u(IsdB:IedB,jsd:jed)) ; CS%Res_fn_u(:,:) = 0.0 allocate(CS%Res_fn_v(isd:ied,JsdB:JedB)) ; CS%Res_fn_v(:,:) = 0.0 - allocate(CS%cg1(isd:ied,jsd:jed)) ; CS%cg1(:,:) = 0.0 - allocate(CS%beta_dx2_h(isd:ied,jsd:jed)) ; CS%beta_dx2_h(:,:) = 0.0 allocate(CS%beta_dx2_q(IsdB:IedB,JsdB:JedB)) ; CS%beta_dx2_q(:,:) = 0.0 allocate(CS%beta_dx2_u(IsdB:IedB,jsd:jed)) ; CS%beta_dx2_u(:,:) = 0.0 allocate(CS%beta_dx2_v(isd:ied,JsdB:JedB)) ; CS%beta_dx2_v(:,:) = 0.0 - allocate(CS%f2_dx2_h(isd:ied,jsd:jed)) ; CS%f2_dx2_h(:,:) = 0.0 allocate(CS%f2_dx2_q(IsdB:IedB,JsdB:JedB)) ; CS%f2_dx2_q(:,:) = 0.0 allocate(CS%f2_dx2_u(IsdB:IedB,jsd:jed)) ; CS%f2_dx2_u(:,:) = 0.0 allocate(CS%f2_dx2_v(isd:ied,JsdB:JedB)) ; CS%f2_dx2_v(:,:) = 0.0 - allocate(CS%Rd_dx_h(isd:ied,jsd:jed)) ; CS%Rd_dx_h(:,:) = 0.0 CS%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, Time, & 'Resolution function for scaling diffusivities', 'Nondim') - CS%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, Time, & - 'Ratio between deformation radius and grid spacing', 'Nondim') call get_param(param_file, mdl, "KH_RES_SCALE_COEF", CS%Res_coef_khth, & "A coefficient that determines how KhTh is scaled away if \n"//& @@ -950,22 +943,11 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "Pedlosky's definition. These definitions differ by a factor\n"//& "of 2 infront of the beta term in the denominator. Gill's"//& "is the more appropriate definition.\n", default=.false.) - - ! Pre-calculate several static expressions for later use. - if (Gill_equatorial_Ld) then; oneOrTwo = 2.0 - else; oneOrTwo = 1.0; endif - - do j=js-1,je+1 ; do i=is-1,ie+1 - CS%f2_dx2_h(i,j) = (G%dxT(i,j)**2 + G%dyT(i,j)**2) * & - max(0.25 * ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & - (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)), & - absurdly_small_freq2) - CS%beta_dx2_h(i,j) = oneOrTwo * (G%dxT(i,j)**2 + G%dyT(i,j)**2) * (sqrt(0.5 * & - ( (((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + & - ((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2) + & - (((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2 + & - ((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2) ) )) - enddo ; enddo + if (Gill_equatorial_Ld) then + oneOrTwo = 2.0 + else + oneOrTwo = 1.0 + endif do J=js-1,Jeq ; do I=is-1,Ieq CS%f2_dx2_q(I,J) = (G%dxBu(I,J)**2 + G%dyBu(I,J)**2) * & @@ -1001,6 +983,43 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) endif + ! Resolution %Rd_dx_h + CS%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, Time, & + 'Ratio between deformation radius and grid spacing', 'Nondim') + CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (CS%id_Rd_dx>0) + + if (CS%calculate_Rd_dx) then + CS%calculate_cg1 = .true. ! We will need %cg1 + allocate(CS%Rd_dx_h(isd:ied,jsd:jed)) ; CS%Rd_dx_h(:,:) = 0.0 + allocate(CS%beta_dx2_h(isd:ied,jsd:jed)); CS%beta_dx2_h(:,:) = 0.0 + allocate(CS%f2_dx2_h(isd:ied,jsd:jed)) ; CS%f2_dx2_h(:,:) = 0.0 + do j=js-1,je+1 ; do i=is-1,ie+1 + CS%f2_dx2_h(i,j) = (G%dxT(i,j)**2 + G%dyT(i,j)**2) * & + max(0.25 * ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & + (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)), & + absurdly_small_freq2) + CS%beta_dx2_h(i,j) = oneOrTwo * (G%dxT(i,j)**2 + G%dyT(i,j)**2) * (sqrt(0.5 * & + ( (((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + & + ((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2) + & + (((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2 + & + ((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2) ) )) + enddo ; enddo + endif + + if (CS%calculate_cg1) then + in_use = .true. + allocate(CS%cg1(isd:ied,jsd:jed)); CS%cg1(:,:) = 0.0 + call wave_speed_init(CS%wave_speed_CSp, use_ebt_mode=CS%Resoln_use_ebt, mono_N2_depth=N2_filter_depth) + endif + + ! If nothing is being stored in this class then deallocate + if (in_use) then + CS%use_variable_mixing = .true. + else + deallocate(CS) + return + endif + end subroutine VarMix_init !> \namespace mom_lateral_mixing_coeffs diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index cbafccfb6a..0d658bbaf6 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -141,7 +141,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS use_VarMix = .false. ; Resoln_scaled = .false. ; use_stored_slopes = .false. khth_use_ebt_struct = .false. if (Associated(VarMix)) then - use_VarMix = VarMix%use_variable_mixing + use_VarMix = VarMix%use_variable_mixing .and. (CS%KHTH_Slope_Cff > 0.) Resoln_scaled = VarMix%Resoln_scaled_KhTh use_stored_slopes = VarMix%use_stored_slopes khth_use_ebt_struct = VarMix%khth_use_ebt_struct diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 179d2e2083..87b3318d1f 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -127,7 +127,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_fla Kh_v ! Tracer mixing coefficient at u-points, in m2 s-1. real :: max_CFL ! The global maximum of the diffusive CFL number. - logical :: use_VarMix, Resoln_scaled, do_online + logical :: use_VarMix, Resoln_scaled, do_online, use_Eady integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts real :: I_numitts ! The inverse of the number of iterations, num_itts. real :: scale ! The fraction of khdt_x or khdt_y that is applied in this @@ -171,10 +171,11 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_fla if (CS%debug) call MOM_tracer_chksum("Before tracer diffusion ", Reg%Tr, ntr, G) - use_VarMix = .false. ; Resoln_scaled = .false. + use_VarMix = .false. ; Resoln_scaled = .false. ; use_Eady = .false. if (Associated(VarMix)) then use_VarMix = VarMix%use_variable_mixing Resoln_scaled = VarMix%Resoln_scaled_KhTr + use_Eady = CS%KhTr_Slope_Cff > 0. endif call cpu_clock_begin(id_clock_pass) @@ -192,7 +193,8 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_fla !$OMP private(Kh_loc,Rd_dx) !$OMP do do j=js,je ; do I=is-1,ie - Kh_loc = CS%KhTr + CS%KhTr_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) + Kh_loc = CS%KhTr + if (use_Eady) Kh_loc = Kh_loc + CS%KhTr_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) if (associated(MEKE%Kh)) & Kh_Loc = Kh_Loc + MEKE%KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) @@ -208,7 +210,8 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_fla enddo ; enddo !$OMP do do J=js-1,je ; do i=is,ie - Kh_loc = CS%KhTr + CS%KhTr_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) + Kh_loc = CS%KhTr + if (use_Eady) Kh_loc = Kh_loc + CS%KhTr_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) if (associated(MEKE%Kh)) & Kh_Loc = Kh_Loc + MEKE%KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max)