Initializes the variables mixing coefficients container.
869 type(time_type),
intent(in) :: Time
870 type(ocean_grid_type),
intent(in) :: G
871 type(verticalGrid_type),
intent(in) :: GV
872 type(unit_scale_type),
intent(in) :: US
873 type(param_file_type),
intent(in) :: param_file
874 type(diag_ctrl),
target,
intent(inout) :: diag
875 type(VarMix_CS),
pointer :: CS
877 real :: KhTr_Slope_Cff, KhTh_Slope_Cff, oneOrTwo, N2_filter_depth
878 real :: KhTr_passivity_coeff
879 real :: absurdly_small_freq
881 logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
882 real :: MLE_front_length
883 real :: Leith_Lap_const
884 real :: grid_sp_u2, grid_sp_u3
885 real :: grid_sp_v2, grid_sp_v3
887 #include "version_variable.h"
888 character(len=40) :: mdl =
"MOM_lateral_mixing_coeffs"
889 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j
890 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
891 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
892 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
893 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
894 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
896 if (
associated(cs))
then
897 call mom_error(warning,
"VarMix_init called with an associated "// &
898 "control structure.")
905 cs%calculate_cg1 = .false.
906 cs%calculate_Rd_dx = .false.
907 cs%calculate_res_fns = .false.
908 cs%calculate_Eady_growth_rate = .false.
911 call log_version(param_file, mdl, version,
"")
912 call get_param(param_file, mdl,
"USE_VARIABLE_MIXING", cs%use_variable_mixing,&
913 "If true, the variable mixing code will be called. This "//&
914 "allows diagnostics to be created even if the scheme is "//&
915 "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
916 "this is set to true regardless of what is in the "//&
917 "parameter file.", default=.false.)
918 call get_param(param_file, mdl,
"USE_VISBECK", cs%use_Visbeck,&
919 "If true, use the Visbeck et al. (1997) formulation for \n"//&
920 "thickness diffusivity.", default=.false.)
921 call get_param(param_file, mdl,
"RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
922 "If true, the Laplacian lateral viscosity is scaled away "//&
923 "when the first baroclinic deformation radius is well "//&
924 "resolved.", default=.false.)
925 call get_param(param_file, mdl,
"RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
926 "If true, the interface depth diffusivity is scaled away "//&
927 "when the first baroclinic deformation radius is well "//&
928 "resolved.", default=.false.)
929 call get_param(param_file, mdl,
"RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
930 "If true, the epipycnal tracer diffusivity is scaled "//&
931 "away when the first baroclinic deformation radius is "//&
932 "well resolved.", default=.false.)
933 call get_param(param_file, mdl,
"RESOLN_USE_EBT", cs%Resoln_use_ebt, &
934 "If true, uses the equivalent barotropic wave speed instead "//&
935 "of first baroclinic wave for calculating the resolution fn.",&
937 call get_param(param_file, mdl,
"KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
938 "If true, uses the equivalent barotropic structure "//&
939 "as the vertical structure of thickness diffusivity.",&
941 call get_param(param_file, mdl,
"KHTH_SLOPE_CFF", khth_slope_cff, &
942 "The nondimensional coefficient in the Visbeck formula "//&
943 "for the interface depth diffusivity", units=
"nondim", &
945 call get_param(param_file, mdl,
"KHTR_SLOPE_CFF", khtr_slope_cff, &
946 "The nondimensional coefficient in the Visbeck formula "//&
947 "for the epipycnal tracer diffusivity", units=
"nondim", &
949 call get_param(param_file, mdl,
"USE_STORED_SLOPES", cs%use_stored_slopes,&
950 "If true, the isopycnal slopes are calculated once and "//&
951 "stored for re-use. This uses more memory but avoids calling "//&
952 "the equation of state more times than should be necessary.", &
954 call get_param(param_file, mdl,
"VERY_SMALL_FREQUENCY", absurdly_small_freq, &
955 "A miniscule frequency that is used to avoid division by 0. The default "//&
956 "value is roughly (pi / (the age of the universe)).", &
957 default=1.0e-17, units=
"s-1", scale=us%T_to_s)
958 call get_param(param_file, mdl,
"KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
959 default=.false., do_not_log=.true.)
960 cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn
961 call get_param(param_file, mdl,
"USE_MEKE", use_meke, &
962 default=.false., do_not_log=.true.)
963 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
964 cs%calculate_Eady_growth_rate = cs%calculate_Eady_growth_rate .or. use_meke
965 call get_param(param_file, mdl,
"KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
966 default=0., do_not_log=.true.)
967 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
968 call get_param(param_file, mdl,
"MLE_FRONT_LENGTH", mle_front_length, &
969 default=0., do_not_log=.true.)
970 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)
972 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false., do_not_log=.true.)
974 if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct)
then
976 call get_param(param_file, mdl,
"RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
977 "The depth below which N2 is monotonized to avoid stratification "//&
978 "artifacts from altering the equivalent barotropic mode structure.",&
979 units=
"m", default=2000.)
980 allocate(cs%ebt_struct(isd:ied,jsd:jed,g%ke)) ; cs%ebt_struct(:,:,:) = 0.0
983 if (khtr_slope_cff>0. .or. khth_slope_cff>0.)
then
984 cs%calculate_Eady_growth_rate = .true.
985 call get_param(param_file, mdl,
"VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
986 "If non-zero, is an upper bound on slopes used in the "//&
987 "Visbeck formula for diffusivity. This does not affect the "//&
988 "isopycnal slope calculation used within thickness diffusion.", &
989 units=
"nondim", default=0.0)
992 if (cs%use_stored_slopes)
then
994 allocate(cs%slope_x(isdb:iedb,jsd:jed,g%ke+1)) ; cs%slope_x(:,:,:) = 0.0
995 allocate(cs%slope_y(isd:ied,jsdb:jedb,g%ke+1)) ; cs%slope_y(:,:,:) = 0.0
996 allocate(cs%N2_u(isdb:iedb,jsd:jed,g%ke+1)) ; cs%N2_u(:,:,:) = 0.0
997 allocate(cs%N2_v(isd:ied,jsdb:jedb,g%ke+1)) ; cs%N2_v(:,:,:) = 0.0
998 call get_param(param_file, mdl,
"KD_SMOOTH", cs%kappa_smooth, &
999 "A diapycnal diffusivity that is used to interpolate "//&
1000 "more sensible values of T & S into thin layers.", &
1001 units=
"m2 s-1", default=1.0e-6, scale=us%m_to_Z**2)
1004 if (cs%calculate_Eady_growth_rate)
then
1006 allocate(cs%SN_u(isdb:iedb,jsd:jed)) ; cs%SN_u(:,:) = 0.0
1007 allocate(cs%SN_v(isd:ied,jsdb:jedb)) ; cs%SN_v(:,:) = 0.0
1008 cs%id_SN_u = register_diag_field(
'ocean_model',
'SN_u', diag%axesCu1, time, &
1009 'Inverse eddy time-scale, S*N, at u-points',
's-1')
1010 cs%id_SN_v = register_diag_field(
'ocean_model',
'SN_v', diag%axesCv1, time, &
1011 'Inverse eddy time-scale, S*N, at v-points',
's-1')
1012 call get_param(param_file, mdl,
"VARMIX_KTOP", cs%VarMix_Ktop, &
1013 "The layer number at which to start vertical integration "//&
1014 "of S*N for purposes of finding the Eady growth rate.", &
1015 units=
"nondim", default=2)
1018 if (khtr_slope_cff>0. .or. khth_slope_cff>0.)
then
1020 call get_param(param_file, mdl,
"VISBECK_L_SCALE", cs%Visbeck_L_scale, &
1021 "The fixed length scale in the Visbeck formula.", units=
"m", &
1023 allocate(cs%L2u(isdb:iedb,jsd:jed)) ; cs%L2u(:,:) = 0.0
1024 allocate(cs%L2v(isd:ied,jsdb:jedb)) ; cs%L2v(:,:) = 0.0
1025 if (cs%Visbeck_L_scale<0)
then
1026 do j=js,je ;
do i=is-1,ieq
1027 cs%L2u(i,j) = cs%Visbeck_L_scale**2*g%areaCu(i,j)
1029 do j=js-1,jeq ;
do i=is,ie
1030 cs%L2v(i,j) = cs%Visbeck_L_scale**2*g%areaCv(i,j)
1033 cs%L2u(:,:) = cs%Visbeck_L_scale**2
1034 cs%L2v(:,:) = cs%Visbeck_L_scale**2
1037 cs%id_L2u = register_diag_field(
'ocean_model',
'L2u', diag%axesCu1, time, &
1038 'Length scale squared for mixing coefficient, at u-points',
'm2')
1039 cs%id_L2v = register_diag_field(
'ocean_model',
'L2v', diag%axesCv1, time, &
1040 'Length scale squared for mixing coefficient, at v-points',
'm2')
1043 if (cs%use_stored_slopes)
then
1044 cs%id_N2_u = register_diag_field(
'ocean_model',
'N2_u', diag%axesCui, time, &
1045 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.',
's-2')
1046 cs%id_N2_v = register_diag_field(
'ocean_model',
'N2_v', diag%axesCvi, time, &
1047 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.',
's-2')
1048 cs%id_S2_u = register_diag_field(
'ocean_model',
'S2_u', diag%axesCu1, time, &
1049 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.',
's-2')
1050 cs%id_S2_v = register_diag_field(
'ocean_model',
'S2_v', diag%axesCv1, time, &
1051 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.',
's-2')
1055 if (cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. cs%Resoln_scaled_KhTr)
then
1056 cs%calculate_Rd_dx = .true.
1057 cs%calculate_res_fns = .true.
1058 allocate(cs%Res_fn_h(isd:ied,jsd:jed)) ; cs%Res_fn_h(:,:) = 0.0
1059 allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb)) ; cs%Res_fn_q(:,:) = 0.0
1060 allocate(cs%Res_fn_u(isdb:iedb,jsd:jed)) ; cs%Res_fn_u(:,:) = 0.0
1061 allocate(cs%Res_fn_v(isd:ied,jsdb:jedb)) ; cs%Res_fn_v(:,:) = 0.0
1062 allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%beta_dx2_q(:,:) = 0.0
1063 allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed)) ; cs%beta_dx2_u(:,:) = 0.0
1064 allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb)) ; cs%beta_dx2_v(:,:) = 0.0
1065 allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%f2_dx2_q(:,:) = 0.0
1066 allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed)) ; cs%f2_dx2_u(:,:) = 0.0
1067 allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb)) ; cs%f2_dx2_v(:,:) = 0.0
1069 cs%id_Res_fn = register_diag_field(
'ocean_model',
'Res_fn', diag%axesT1, time, &
1070 'Resolution function for scaling diffusivities',
'nondim')
1072 call get_param(param_file, mdl,
"KH_RES_SCALE_COEF", cs%Res_coef_khth, &
1073 "A coefficient that determines how KhTh is scaled away if "//&
1074 "RESOLN_SCALED_... is true, as "//&
1075 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
1076 units=
"nondim", default=1.0)
1077 call get_param(param_file, mdl,
"KH_RES_FN_POWER", cs%Res_fn_power_khth, &
1078 "The power of dx/Ld in the Kh resolution function. Any "//&
1079 "positive integer may be used, although even integers "//&
1080 "are more efficient to calculate. Setting this greater "//&
1081 "than 100 results in a step-function being used.", &
1082 units=
"nondim", default=2)
1083 call get_param(param_file, mdl,
"VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
1084 "A coefficient that determines how Kh is scaled away if "//&
1085 "RESOLN_SCALED_... is true, as "//&
1086 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
1087 "This function affects lateral viscosity, Kh, and not KhTh.", &
1088 units=
"nondim", default=cs%Res_coef_khth)
1089 call get_param(param_file, mdl,
"VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
1090 "The power of dx/Ld in the Kh resolution function. Any "//&
1091 "positive integer may be used, although even integers "//&
1092 "are more efficient to calculate. Setting this greater "//&
1093 "than 100 results in a step-function being used. "//&
1094 "This function affects lateral viscosity, Kh, and not KhTh.", &
1095 units=
"nondim", default=cs%Res_fn_power_khth)
1096 call get_param(param_file, mdl,
"INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
1097 "If true, interpolate the resolution function to the "//&
1098 "velocity points from the thickness points; otherwise "//&
1099 "interpolate the wave speed and calculate the resolution "//&
1100 "function independently at each point.", default=.true.)
1101 if (cs%interpolate_Res_fn)
then
1102 if (cs%Res_coef_visc /= cs%Res_coef_khth)
call mom_error(fatal, &
1103 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1104 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
1105 if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth)
call mom_error(fatal, &
1106 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1107 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
1110 call get_param(param_file, mdl,
"GILL_EQUATORIAL_LD", gill_equatorial_ld, &
1111 "If true, uses Gill's definition of the baroclinic "//&
1112 "equatorial deformation radius, otherwise, if false, use "//&
1113 "Pedlosky's definition. These definitions differ by a factor "//&
1114 "of 2 in front of the beta term in the denominator. Gill's "//&
1115 "is the more appropriate definition.", default=.false.)
1116 if (gill_equatorial_ld)
then
1120 do j=js-1,jeq ;
do i=is-1,ieq
1121 cs%f2_dx2_q(i,j) = (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * &
1122 max(g%CoriolisBu(i,j)**2, absurdly_small_freq**2)
1123 cs%beta_dx2_q(i,j) = oneortwo * (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * (sqrt(0.5 * &
1124 ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1125 ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1126 (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1127 ((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) ) ))
1130 do j=js,je ;
do i=is-1,ieq
1131 cs%f2_dx2_u(i,j) = (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * &
1132 max(0.5* (g%CoriolisBu(i,j)**2+g%CoriolisBu(i,j-1)**2), absurdly_small_freq**2)
1133 cs%beta_dx2_u(i,j) = oneortwo * (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * (sqrt( &
1134 0.25*( (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2 + &
1135 ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1136 (((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2 + &
1137 ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) ) + &
1138 ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 ))
1141 do j=js-1,jeq ;
do i=is,ie
1142 cs%f2_dx2_v(i,j) = (g%dxCv(i,j)**2 + g%dyCv(i,j)**2) * &
1143 max(0.5*(g%CoriolisBu(i,j)**2+g%CoriolisBu(i-1,j)**2), absurdly_small_freq**2)
1144 cs%beta_dx2_v(i,j) = oneortwo * (g%dxCv(i,j)**2 + g%dyCv(i,j)**2) * (sqrt( &
1145 ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1146 0.25*( (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1147 ((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2) + &
1148 (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2 + &
1149 ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1155 cs%id_Rd_dx = register_diag_field(
'ocean_model',
'Rd_dx', diag%axesT1, time, &
1156 'Ratio between deformation radius and grid spacing',
'm m-1')
1157 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)
1159 if (cs%calculate_Rd_dx)
then
1160 cs%calculate_cg1 = .true.
1161 allocate(cs%Rd_dx_h(isd:ied,jsd:jed)) ; cs%Rd_dx_h(:,:) = 0.0
1162 allocate(cs%beta_dx2_h(isd:ied,jsd:jed)); cs%beta_dx2_h(:,:) = 0.0
1163 allocate(cs%f2_dx2_h(isd:ied,jsd:jed)) ; cs%f2_dx2_h(:,:) = 0.0
1164 do j=js-1,je+1 ;
do i=is-1,ie+1
1165 cs%f2_dx2_h(i,j) = (g%dxT(i,j)**2 + g%dyT(i,j)**2) * &
1166 max(0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1167 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)), &
1168 absurdly_small_freq**2)
1169 cs%beta_dx2_h(i,j) = oneortwo * (g%dxT(i,j)**2 + g%dyT(i,j)**2) * (sqrt(0.5 * &
1170 ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1171 ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
1172 (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1173 ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1177 if (cs%calculate_cg1)
then
1179 allocate(cs%cg1(isd:ied,jsd:jed)); cs%cg1(:,:) = 0.0
1180 call wave_speed_init(cs%wave_speed_CSp, use_ebt_mode=cs%Resoln_use_ebt, mono_n2_depth=n2_filter_depth)
1184 call get_param(param_file, mdl,
"USE_QG_LEITH_GM", cs%use_QG_Leith_GM, &
1185 "If true, use the QG Leith viscosity as the GM coefficient.", &
1188 if (cs%Use_QG_Leith_GM)
then
1189 call get_param(param_file, mdl,
"LEITH_LAP_CONST", leith_lap_const, &
1190 "The nondimensional Laplacian Leith constant, \n"//&
1191 "often set to 1.0", units=
"nondim", default=0.0)
1193 call get_param(param_file, mdl,
"USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
1194 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
1197 alloc_(cs%Laplac3_const_u(isdb:iedb,jsd:jed)) ; cs%Laplac3_const_u(:,:) = 0.0
1198 alloc_(cs%Laplac3_const_v(isd:ied,jsdb:jedb)) ; cs%Laplac3_const_v(:,:) = 0.0
1199 alloc_(cs%KH_u_QG(isdb:iedb,jsd:jed,g%ke)) ; cs%KH_u_QG(:,:,:) = 0.0
1200 alloc_(cs%KH_v_QG(isd:ied,jsdb:jedb,g%ke)) ; cs%KH_v_QG(:,:,:) = 0.0
1203 cs%id_KH_u_QG = register_diag_field(
'ocean_model',
'KH_u_QG', diag%axesCuL, time, &
1204 'Horizontal viscosity from Leith QG, at u-points',
'm2 s-1')
1205 cs%id_KH_v_QG = register_diag_field(
'ocean_model',
'KH_v_QG', diag%axesCvL, time, &
1206 'Horizontal viscosity from Leith QG, at v-points',
'm2 s-1')
1208 do j=jsq,jeq+1 ;
do i=is-1,ieq
1210 grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
1211 grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
1212 cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
1214 do j=js-1,jeq ;
do i=isq,ieq+1
1216 grid_sp_v2 = g%dyCv(i,j)*g%dxCu(i,j)
1217 grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
1218 cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
1221 if (.not. cs%use_stored_slopes)
call mom_error(fatal, &
1222 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1223 "USE_STORED_SLOPES must be True when using QG Leith.")
1228 cs%use_variable_mixing = .true.