portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Interface to background mixing schemes, including the Bryan and Lewis (1979) | |
| 6 | !! which is applied via CVMix. | |
| 7 | ||
| 8 | module MOM_bkgnd_mixing | |
| 9 | ||
| 10 | use MOM_debugging, only : hchksum | |
| 11 | use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field | |
| 12 | use MOM_diag_mediator, only : post_data | |
| 13 | use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE | |
| 14 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 15 | use MOM_file_parser, only : openParameterBlock, closeParameterBlock | |
| 16 | use MOM_forcing_type, only : forcing | |
| 17 | use MOM_grid, only : ocean_grid_type | |
| 18 | use MOM_interface_heights, only : thickness_to_dz | |
| 19 | use MOM_unit_scaling, only : unit_scale_type | |
| 20 | use MOM_verticalGrid, only : verticalGrid_type | |
| 21 | use MOM_variables, only : thermo_var_ptrs, vertvisc_type | |
| 22 | use MOM_intrinsic_functions, only : invcosh | |
| 23 | use CVMix_background, only : CVMix_init_bkgnd, CVMix_coeffs_bkgnd | |
| 24 | ||
| 25 | implicit none ; private | |
| 26 | ||
| 27 | #include <MOM_memory.h> | |
| 28 | ||
| 29 | public bkgnd_mixing_init | |
| 30 | public bkgnd_mixing_end | |
| 31 | public calculate_bkgnd_mixing | |
| 32 | ||
| 33 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 34 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 35 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 36 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 37 | ||
| 38 | !> Control structure including parameters for this module. | |
| 39 | type, public :: bkgnd_mixing_cs ; private | |
| 40 | ||
| 41 | ! Parameters | |
| 42 | real :: Bryan_Lewis_c1 !< The vertical diffusivity values for Bryan-Lewis profile | |
| 43 | !! at |z|=D [Z2 T-1 ~> m2 s-1] | |
| 44 | real :: Bryan_Lewis_c2 !< The amplitude of variation in diffusivity for the | |
| 45 | !! Bryan-Lewis diffusivity profile [Z2 T-1 ~> m2 s-1] | |
| 46 | real :: Bryan_Lewis_c3 !< The inverse length scale for transition region in the | |
| 47 | !! Bryan-Lewis diffusivity profile [Z-1 ~> m-1] | |
| 48 | real :: Bryan_Lewis_c4 !< The depth where diffusivity is Bryan_Lewis_bl1 in the | |
| 49 | !! Bryan-Lewis profile [Z ~> m] | |
| 50 | real :: bckgrnd_vdc1 !< Background diffusivity (Ledwell) when | |
| 51 | !! horiz_varying_background=.true. [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 52 | real :: bckgrnd_vdc_eq !< Equatorial diffusivity (Gregg) when | |
| 53 | !! horiz_varying_background=.true. [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 54 | real :: bckgrnd_vdc_psim !< Max. PSI induced diffusivity (MacKinnon) when | |
| 55 | !! horiz_varying_background=.true. [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 56 | real :: bckgrnd_vdc_Banda !< Banda Sea diffusivity (Gordon) when | |
| 57 | !! horiz_varying_background=.true. [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 58 | real :: Kd_min !< minimum diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 59 | real :: Kd !< interior diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 60 | real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. | |
| 61 | real :: N0_2Omega !< ratio of the typical Buoyancy frequency to | |
| 62 | !! twice the Earth's rotation period, used with the | |
| 63 | !! Henyey scaling from the mixing [nondim] | |
| 64 | real :: Henyey_max_lat !< A latitude poleward of which the Henyey profile | |
| 65 | !! is returned to the minimum diffusivity [degrees_N] | |
| 66 | real :: prandtl_bkgnd !< Turbulent Prandtl number used to convert | |
| 67 | !! vertical background diffusivity into viscosity [nondim] | |
| 68 | real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of | |
| 69 | !! diffusivities with Kd_tanh_lat_fn [nondim]. Valid values | |
| 70 | !! are in the range of -2 to 2; 0.4 reproduces CM2M. | |
| 71 | real :: Kd_tot_ml !< The mixed layer diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 72 | !! when no other physically based mixed layer turbulence | |
| 73 | !! parameterization is being used. | |
| 74 | real :: Hmix !< mixed layer thickness [H ~> m or kg m-2] when no physically based | |
| 75 | !! ocean surface boundary layer parameterization is used. | |
| 76 | logical :: Kd_tanh_lat_fn !< If true, use the tanh dependence of Kd_sfc on | |
| 77 | !! latitude, like GFDL CM2.1/CM2M. There is no | |
| 78 | !! physical justification for this form, and it can | |
| 79 | !! not be used with Henyey_IGW_background. | |
| 80 | logical :: Bryan_Lewis_diffusivity!< If true, background vertical diffusivity | |
| 81 | !! uses Bryan-Lewis (1979) like tanh profile. | |
| 82 | logical :: horiz_varying_background !< If true, apply vertically uniform, latitude-dependent | |
| 83 | !! background diffusivity, as described in Danabasoglu et al., 2012 | |
| 84 | logical :: Henyey_IGW_background !< If true, use a simplified variant of the | |
| 85 | !! Henyey et al, JGR (1986) latitudinal scaling for the background diapycnal diffusivity, | |
| 86 | !! which gives a marked decrease in the diffusivity near the equator. The simplification | |
| 87 | !! here is to assume that the in-situ stratification is the same as the reference stratificaiton. | |
| 88 | logical :: physical_OBL_scheme !< If true, a physically-based scheme is used to determine mixing in the | |
| 89 | !! ocean's surface boundary layer, such as ePBL, KPP, or a refined bulk mixed layer scheme. | |
| 90 | logical :: debug !< If true, turn on debugging in this module | |
| 91 | ! Diagnostic handles and pointers | |
| 92 | type(diag_ctrl), pointer :: diag => NULL() !< A structure that regulates diagnostic output | |
| 93 | ||
| 94 | character(len=40) :: bkgnd_scheme_str = "none" !< Background scheme identifier | |
| 95 | ||
| 96 | end type bkgnd_mixing_cs | |
| 97 | ||
| 98 | character(len=40) :: mdl = "MOM_bkgnd_mixing" !< This module's name. | |
| 99 | ||
| 100 | contains | |
| 101 | ||
| 102 | !> Initialize the background mixing routine. | |
| 103 | 1 | subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL_scheme) |
| 104 | ||
| 105 | type(time_type), intent(in) :: Time !< The current time. | |
| 106 | type(ocean_grid_type), intent(in) :: G !< Grid structure. | |
| 107 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 108 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 109 | type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle | |
| 110 | type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. | |
| 111 | type(bkgnd_mixing_cs), pointer :: CS !< This module's control structure. | |
| 112 | logical, intent(in) :: physical_OBL_scheme !< If true, a physically based | |
| 113 | !! parameterization (like KPP or ePBL or a bulk mixed | |
| 114 | !! layer) is used outside of set_diffusivity to | |
| 115 | !! specify the mixing that occurs in the ocean's | |
| 116 | !! surface boundary layer. | |
| 117 | ||
| 118 | ! Local variables | |
| 119 | real :: Kv ! The interior vertical viscosity [H Z T-1 ~> m2 s-1 or Pa s] - read to set Prandtl | |
| 120 | ! number unless it is provided as a parameter | |
| 121 | real :: Kd_z ! The background diapycnal diffusivity in [Z2 T-1 ~> m2 s-1] for use | |
| 122 | ! in setting the default for other diffusivities. | |
| 123 | real :: prandtl_bkgnd_comp ! Kv/CS%Kd [nondim]. Gets compared with user-specified prandtl_bkgnd. | |
| 124 | ||
| 125 | ! This include declares and sets the variable "version". | |
| 126 | # include "version_variable.h" | |
| 127 | ||
| 128 | 1 | if (associated(CS)) then |
| 129 | call MOM_error(WARNING, "bkgnd_mixing_init called with an associated "// & | |
| 130 | 0 | "control structure.") |
| 131 | 0 | return |
| 132 | endif | |
| 133 | 1 | allocate(CS) |
| 134 | ||
| 135 | ! Read parameters | |
| 136 | call log_version(param_file, mdl, version, & | |
| 137 | 1 | "Adding static vertical background mixing coefficients") |
| 138 | ||
| 139 | call get_param(param_file, mdl, "KD", Kd_z, & | |
| 140 | "The background diapycnal diffusivity of density in the "//& | |
| 141 | "interior. Zero or the molecular value, ~1e-7 m2 s-1, "//& | |
| 142 | 1 | "may be used.", default=0.0, units="m2 s-1", scale=US%m2_s_to_Z2_T) |
| 143 | 1 | CS%Kd = (GV%m2_s_to_HZ_T*US%Z2_T_to_m2_s) * Kd_z |
| 144 | ||
| 145 | call get_param(param_file, mdl, "KV", Kv, & | |
| 146 | "The background kinematic viscosity in the interior. "//& | |
| 147 | "The molecular value, ~1e-6 m2 s-1, may be used.", & | |
| 148 | 1 | units="m2 s-1", scale=GV%m2_s_to_HZ_T, fail_if_missing=.true.) |
| 149 | ||
| 150 | call get_param(param_file, mdl, "KD_MIN", CS%Kd_min, & | |
| 151 | "The minimum diapycnal diffusivity.", & | |
| 152 | 1 | units="m2 s-1", default=0.01*Kd_z*US%Z2_T_to_m2_s, scale=GV%m2_s_to_HZ_T) |
| 153 | ||
| 154 | ! The following is needed to set one of the choices of vertical background mixing | |
| 155 | ||
| 156 | 1 | CS%physical_OBL_scheme = physical_OBL_scheme |
| 157 | 1 | if (CS%physical_OBL_scheme) then |
| 158 | ! Check that Kdml is not set when using bulk mixed layer | |
| 159 | call get_param(param_file, mdl, "KDML", CS%Kd_tot_ml, & | |
| 160 | 1 | units="m2 s-1", default=-1., scale=GV%m2_s_to_HZ_T, do_not_log=.true.) |
| 161 | 1 | if (CS%Kd_tot_ml>0.) call MOM_error(FATAL, & |
| 162 | 0 | "bkgnd_mixing_init: KDML is a depricated parameter that should not be used.") |
| 163 | call get_param(param_file, mdl, "KD_ML_TOT", CS%Kd_tot_ml, & | |
| 164 | 1 | units="m2 s-1", default=-1.0, scale=GV%m2_s_to_HZ_T, do_not_log=.true.) |
| 165 | 1 | if (CS%Kd_tot_ml>0.) call MOM_error(FATAL, & |
| 166 | "bkgnd_mixing_init: KD_ML_TOT cannot be set when using a physically based ocean "//& | |
| 167 | 0 | "boundary layer mixing parameterization.") |
| 168 | 1 | CS%Kd_tot_ml = CS%Kd ! This is not used with a bulk mixed layer, but also cannot be a NaN. |
| 169 | else | |
| 170 | call get_param(param_file, mdl, "KD_ML_TOT", CS%Kd_tot_ml, & | |
| 171 | "The total diapcynal diffusivity in the surface mixed layer when there is "//& | |
| 172 | "not a physically based parameterization of mixing in the mixed layer, such "//& | |
| 173 | "as bulk mixed layer or KPP or ePBL.", & | |
| 174 | 0 | units="m2 s-1", default=Kd_z*US%Z2_T_to_m2_s, scale=GV%m2_s_to_HZ_T, do_not_log=.true.) |
| 175 | 0 | if (abs(CS%Kd_tot_ml - CS%Kd) <= 1.0e-15*abs(CS%Kd)) then |
| 176 | call get_param(param_file, mdl, "KDML", CS%Kd_tot_ml, & | |
| 177 | "If BULKMIXEDLAYER is false, KDML is the elevated "//& | |
| 178 | "diapycnal diffusivity in the topmost HMIX of fluid. "//& | |
| 179 | "KDML is only used if BULKMIXEDLAYER is false.", & | |
| 180 | 0 | units="m2 s-1", default=Kd_z*US%Z2_T_to_m2_s, scale=GV%m2_s_to_HZ_T, do_not_log=.true.) |
| 181 | 0 | if (abs(CS%Kd_tot_ml - CS%Kd) > 1.0e-15*abs(CS%Kd)) & |
| 182 | 0 | call MOM_error(WARNING, "KDML is a depricated parameter. Use KD_ML_TOT instead.") |
| 183 | endif | |
| 184 | call log_param(param_file, mdl, "KD_ML_TOT", CS%Kd_tot_ml, & | |
| 185 | "The total diapcynal diffusivity in the surface mixed layer when there is "//& | |
| 186 | "not a physically based parameterization of mixing in the mixed layer, such "//& | |
| 187 | "as bulk mixed layer or KPP or ePBL.", & | |
| 188 | 0 | units="m2 s-1", default=Kd_z*US%Z2_T_to_m2_s, unscale=GV%HZ_T_to_m2_s) |
| 189 | ||
| 190 | call get_param(param_file, mdl, "HMIX_FIXED", CS%Hmix, & | |
| 191 | "The prescribed depth over which the near-surface "//& | |
| 192 | "viscosity and diffusivity are elevated when the bulk "//& | |
| 193 | 0 | "mixed layer is not used.", units="m", scale=GV%m_to_H, fail_if_missing=.true.) |
| 194 | endif | |
| 195 | ||
| 196 | 1 | call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) |
| 197 | ||
| 198 | ! call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING') | |
| 199 | ||
| 200 | call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", CS%Bryan_Lewis_diffusivity, & | |
| 201 | "If true, use a Bryan & Lewis (JGR 1979) like tanh "//& | |
| 202 | "profile of background diapycnal diffusivity with depth. "//& | |
| 203 | 1 | "This is done via CVMix.", default=.false.) |
| 204 | ||
| 205 | 1 | if (CS%Bryan_Lewis_diffusivity) then |
| 206 | 0 | call check_bkgnd_scheme(CS, "BRYAN_LEWIS_DIFFUSIVITY") |
| 207 | ||
| 208 | call get_param(param_file, mdl, "BRYAN_LEWIS_C1", CS%Bryan_Lewis_c1, & | |
| 209 | "The vertical diffusivity values for Bryan-Lewis profile at |z|=D.", & | |
| 210 | 0 | units="m2 s-1", scale=US%m2_s_to_Z2_T, fail_if_missing=.true.) |
| 211 | ||
| 212 | call get_param(param_file, mdl, "BRYAN_LEWIS_C2", CS%Bryan_Lewis_c2, & | |
| 213 | "The amplitude of variation in diffusivity for the Bryan-Lewis profile", & | |
| 214 | 0 | units="m2 s-1", scale=US%m2_s_to_Z2_T, fail_if_missing=.true.) |
| 215 | ||
| 216 | call get_param(param_file, mdl, "BRYAN_LEWIS_C3", CS%Bryan_Lewis_c3, & | |
| 217 | "The inverse length scale for transition region in the Bryan-Lewis profile", & | |
| 218 | 0 | units="m-1", scale=US%Z_to_m, fail_if_missing=.true.) |
| 219 | ||
| 220 | call get_param(param_file, mdl, "BRYAN_LEWIS_C4", CS%Bryan_Lewis_c4, & | |
| 221 | "The depth where diffusivity is BRYAN_LEWIS_C1 in the Bryan-Lewis profile",& | |
| 222 | 0 | units="m", scale=US%m_to_Z, fail_if_missing=.true.) |
| 223 | ||
| 224 | endif ! CS%Bryan_Lewis_diffusivity | |
| 225 | ||
| 226 | call get_param(param_file, mdl, "HORIZ_VARYING_BACKGROUND", CS%horiz_varying_background, & | |
| 227 | "If true, apply vertically uniform, latitude-dependent background "//& | |
| 228 | "diffusivity, as described in Danabasoglu et al., 2012", & | |
| 229 | 1 | default=.false.) |
| 230 | ||
| 231 | 1 | if (CS%horiz_varying_background) then |
| 232 | 0 | call check_bkgnd_scheme(CS, "HORIZ_VARYING_BACKGROUND") |
| 233 | ||
| 234 | call get_param(param_file, mdl, "BCKGRND_VDC1", CS%bckgrnd_vdc1, & | |
| 235 | "Background diffusivity (Ledwell) when HORIZ_VARYING_BACKGROUND=True", & | |
| 236 | 0 | units="m2 s-1", default=0.16e-04, scale=GV%m2_s_to_HZ_T) |
| 237 | ||
| 238 | call get_param(param_file, mdl, "BCKGRND_VDC_EQ", CS%bckgrnd_vdc_eq, & | |
| 239 | "Equatorial diffusivity (Gregg) when HORIZ_VARYING_BACKGROUND=True", & | |
| 240 | 0 | units="m2 s-1", default=0.01e-04, scale=GV%m2_s_to_HZ_T) |
| 241 | ||
| 242 | call get_param(param_file, mdl, "BCKGRND_VDC_PSIM", CS%bckgrnd_vdc_psim, & | |
| 243 | "Max. PSI induced diffusivity (MacKinnon) when HORIZ_VARYING_BACKGROUND=True", & | |
| 244 | 0 | units="m2 s-1", default=0.13e-4, scale=GV%m2_s_to_HZ_T) |
| 245 | ||
| 246 | call get_param(param_file, mdl, "BCKGRND_VDC_BAN", CS%bckgrnd_vdc_Banda, & | |
| 247 | "Banda Sea diffusivity (Gordon) when HORIZ_VARYING_BACKGROUND=True", & | |
| 248 | 0 | units="m2 s-1", default=1.0e-4, scale=GV%m2_s_to_HZ_T) |
| 249 | endif | |
| 250 | ||
| 251 | call get_param(param_file, mdl, "PRANDTL_BKGND", CS%prandtl_bkgnd, & | |
| 252 | "Turbulent Prandtl number used to convert vertical "//& | |
| 253 | "background diffusivities into viscosities.", & | |
| 254 | 1 | units="nondim", default=1.0) |
| 255 | ||
| 256 | 1 | if (CS%Bryan_Lewis_diffusivity .or. CS%horiz_varying_background) then |
| 257 | 0 | prandtl_bkgnd_comp = CS%prandtl_bkgnd |
| 258 | 0 | if (CS%Kd /= 0.0) prandtl_bkgnd_comp = Kv / CS%Kd |
| 259 | ||
| 260 | 0 | if ( abs(CS%prandtl_bkgnd - prandtl_bkgnd_comp)>1.e-14) then |
| 261 | call MOM_error(FATAL, "bkgnd_mixing_init: The provided KD, KV and PRANDTL_BKGND values "//& | |
| 262 | 0 | "are incompatible. The following must hold: KD*PRANDTL_BKGND==KV") |
| 263 | endif | |
| 264 | endif | |
| 265 | ||
| 266 | call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", CS%Henyey_IGW_background, & | |
| 267 | "If true, use a latitude-dependent scaling for the near "//& | |
| 268 | "surface background diffusivity, as described in "//& | |
| 269 | 1 | "Harrison & Hallberg, JPO 2008.", default=.false.) |
| 270 | 1 | if (CS%Henyey_IGW_background) call check_bkgnd_scheme(CS, "HENYEY_IGW_BACKGROUND") |
| 271 | ||
| 272 | 1 | if (CS%Kd>0.0 .and. (trim(CS%bkgnd_scheme_str)=="BRYAN_LEWIS_DIFFUSIVITY" .or.& |
| 273 | trim(CS%bkgnd_scheme_str)=="HORIZ_VARYING_BACKGROUND" )) then | |
| 274 | call MOM_error(WARNING, "bkgnd_mixing_init: a nonzero constant background "//& | |
| 275 | 0 | "diffusivity (KD) is specified along with "//trim(CS%bkgnd_scheme_str)) |
| 276 | endif | |
| 277 | ||
| 278 | 1 | if (CS%Henyey_IGW_background) then |
| 279 | call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", CS%N0_2Omega, & | |
| 280 | "The ratio of the typical Buoyancy frequency to twice "//& | |
| 281 | "the Earth's rotation period, used with the Henyey "//& | |
| 282 | 1 | "scaling from the mixing.", units="nondim", default=20.0) |
| 283 | call get_param(param_file, mdl, "OMEGA", CS%omega, & | |
| 284 | "The rotation rate of the earth.", & | |
| 285 | 1 | units="s-1", default=7.2921e-5, scale=US%T_to_s) |
| 286 | call get_param(param_file, mdl, "HENYEY_MAX_LAT", CS%Henyey_max_lat, & | |
| 287 | "A latitude poleward of which the Henyey profile "//& | |
| 288 | "is returned to the minimum diffusivity", & | |
| 289 | 1 | units="degN", default=95.0) |
| 290 | endif | |
| 291 | ||
| 292 | call get_param(param_file, mdl, "KD_TANH_LAT_FN", CS%Kd_tanh_lat_fn, & | |
| 293 | "If true, use a tanh dependence of Kd_sfc on latitude, "//& | |
| 294 | "like CM2.1/CM2M. There is no physical justification "//& | |
| 295 | "for this form, and it can not be used with "//& | |
| 296 | 1 | "HENYEY_IGW_BACKGROUND.", default=.false.) |
| 297 | ||
| 298 | 1 | if (CS%Kd_tanh_lat_fn) & |
| 299 | call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", CS%Kd_tanh_lat_scale, & | |
| 300 | "A nondimensional scaling for the range ofdiffusivities "//& | |
| 301 | "with KD_TANH_LAT_FN. Valid values are in the range of "//& | |
| 302 | 0 | "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0) |
| 303 | ||
| 304 | 1 | if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & |
| 305 | 0 | "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") |
| 306 | ||
| 307 | ! call closeParameterBlock(param_file) | |
| 308 | ||
| 309 | end subroutine bkgnd_mixing_init | |
| 310 | ||
| 311 | !> Calculates the vertical background diffusivities/viscosities | |
| 312 | 720 | subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, GV, US, CS) |
| 313 | ||
| 314 | type(ocean_grid_type), intent(in) :: G !< Grid structure. | |
| 315 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 316 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. | |
| 317 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. | |
| 318 | real, dimension(SZI_(G),SZK_(GV)), intent(in) :: N2_lay !< squared buoyancy frequency associated | |
| 319 | !! with layers [T-2 ~> s-2] | |
| 320 | real, dimension(SZI_(G),SZK_(GV)), intent(out) :: Kd_lay !< The background diapycnal diffusivity of each | |
| 321 | !! layer [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 322 | real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: Kd_int !< The background diapycnal diffusivity of each | |
| 323 | !! interface [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 324 | real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: Kv_bkgnd !< The background vertical viscosity at | |
| 325 | !! each interface [H Z T-1 ~> m2 s-1 or Pa s] | |
| 326 | integer, intent(in) :: j !< Meridional grid index | |
| 327 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 328 | type(bkgnd_mixing_cs), pointer :: CS !< The control structure returned by | |
| 329 | !! a previous call to bkgnd_mixing_init. | |
| 330 | ||
| 331 | ! local variables | |
| 332 | 1440 | real, dimension(SZK_(GV)+1) :: depth_int !< Distance from surface of the interfaces [m] |
| 333 | 1440 | real, dimension(SZK_(GV)+1) :: Kd_col !< Diffusivities at the interfaces [m2 s-1] |
| 334 | 1440 | real, dimension(SZK_(GV)+1) :: Kv_col !< Viscosities at the interfaces [m2 s-1] |
| 335 | 1440 | real, dimension(SZI_(G)) :: Kd_sfc !< Surface value of the diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] |
| 336 | 1440 | real, dimension(SZI_(G)) :: depth !< Distance from surface of an interface [H ~> m or kg m-2] |
| 337 | 1440 | real, dimension(SZI_(G),SZK_(GV)) :: dz !< Height change across layers [Z ~> m] |
| 338 | real :: depth_c !< depth of the center of a layer [H ~> m or kg m-2] | |
| 339 | real :: I_Hmix !< inverse of fixed mixed layer thickness [H-1 ~> m-1 or m2 kg-1] | |
| 340 | real :: I_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg))) [nondim] | |
| 341 | real :: deg_to_rad !< factor converting degrees to radians [radians degree-1], pi/180. | |
| 342 | real :: abs_sinlat !< absolute value of sine of latitude [nondim] | |
| 343 | real :: min_sinlat ! The minimum value of the sine of latitude [nondim] | |
| 344 | real :: bckgrnd_vdc_psin !< PSI diffusivity in northern hemisphere [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 345 | real :: bckgrnd_vdc_psis !< PSI diffusivity in southern hemisphere [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 346 | integer :: i, k, is, ie, js, je, nz | |
| 347 | ||
| 348 | 720 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 349 | ||
| 350 | ! set some parameters | |
| 351 | 720 | deg_to_rad = atan(1.0)/45.0 ! = PI/180 |
| 352 | 720 | min_sinlat = 1.e-10 |
| 353 | ||
| 354 | ! Start with a constant value that may be replaced below. | |
| 355 | 6966720 | Kd_lay(:,:) = CS%Kd |
| 356 | 7059600 | Kv_bkgnd(:,:) = 0.0 |
| 357 | ||
| 358 | ! Set up the background diffusivity. | |
| 359 | 720 | if (CS%Bryan_Lewis_diffusivity) then |
| 360 | ||
| 361 | 0 | call thickness_to_dz(h, tv, dz, j, G, GV) |
| 362 | ||
| 363 | 0 | do i=is,ie |
| 364 | 0 | depth_int(1) = 0.0 |
| 365 | 0 | do k=2,nz+1 |
| 366 | 0 | depth_int(k) = depth_int(k-1) + US%Z_to_m*dz(i,k-1) |
| 367 | enddo | |
| 368 | ||
| 369 | call CVMix_init_bkgnd(max_nlev=nz, & | |
| 370 | zw = depth_int(:), & !< interface depths relative to the surface in m, must be positive. | |
| 371 | bl1 = US%Z2_T_to_m2_s*CS%Bryan_Lewis_c1, & | |
| 372 | bl2 = US%Z2_T_to_m2_s*CS%Bryan_Lewis_c2, & | |
| 373 | bl3 = US%m_to_Z*CS%Bryan_Lewis_c3, & | |
| 374 | bl4 = US%Z_to_m*CS%Bryan_Lewis_c4, & | |
| 375 | 0 | prandtl = CS%prandtl_bkgnd) |
| 376 | ||
| 377 | 0 | Kd_col(:) = 0.0 ; Kv_col(:) = 0.0 ! Is this line necessary? |
| 378 | 0 | call CVMix_coeffs_bkgnd(Mdiff_out=Kv_col, Tdiff_out=Kd_col, nlev=nz, max_nlev=nz) |
| 379 | ||
| 380 | ! Update Kd and Kv. | |
| 381 | 0 | do K=1,nz+1 |
| 382 | 0 | Kv_bkgnd(i,K) = GV%m2_s_to_HZ_T * Kv_col(K) |
| 383 | 0 | Kd_int(i,K) = GV%m2_s_to_HZ_T*Kd_col(K) |
| 384 | enddo | |
| 385 | 0 | do k=1,nz |
| 386 | 0 | Kd_lay(i,k) = Kd_lay(i,k) + 0.5 * GV%m2_s_to_HZ_T * (Kd_col(K) + Kd_col(K+1)) |
| 387 | enddo | |
| 388 | enddo ! i loop | |
| 389 | ||
| 390 | 720 | elseif (CS%horiz_varying_background) then |
| 391 | !### Note that there are lots of hard-coded parameters (mostly latitudes and longitudes) here. | |
| 392 | 0 | do i=is,ie |
| 393 | 0 | bckgrnd_vdc_psis = CS%bckgrnd_vdc_psim * exp(-(0.4*(G%geoLatT(i,j)+28.9))**2) |
| 394 | 0 | bckgrnd_vdc_psin = CS%bckgrnd_vdc_psim * exp(-(0.4*(G%geoLatT(i,j)-28.9))**2) |
| 395 | 0 | Kd_int(i,1) = (CS%bckgrnd_vdc_eq + bckgrnd_vdc_psin) + bckgrnd_vdc_psis |
| 396 | ||
| 397 | 0 | if (G%geoLatT(i,j) < -10.0) then |
| 398 | 0 | Kd_int(i,1) = Kd_int(i,1) + CS%bckgrnd_vdc1 |
| 399 | 0 | elseif (G%geoLatT(i,j) <= 10.0) then |
| 400 | 0 | Kd_int(i,1) = Kd_int(i,1) + CS%bckgrnd_vdc1 * (G%geoLatT(i,j)/10.0)**2 |
| 401 | else | |
| 402 | 0 | Kd_int(i,1) = Kd_int(i,1) + CS%bckgrnd_vdc1 |
| 403 | endif | |
| 404 | ||
| 405 | ! North Banda Sea | |
| 406 | if ( (G%geoLatT(i,j) < -1.0) .and. (G%geoLatT(i,j) > -4.0) .and. & | |
| 407 | 0 | ( mod(G%geoLonT(i,j)+360.0,360.0) > 103.0) .and. & |
| 408 | ( mod(G%geoLonT(i,j)+360.0,360.0) < 134.0) ) then | |
| 409 | 0 | Kd_int(i,1) = CS%bckgrnd_vdc_Banda |
| 410 | endif | |
| 411 | ||
| 412 | ! Middle Banda Sea | |
| 413 | if ( (G%geoLatT(i,j) <= -4.0) .and. (G%geoLatT(i,j) > -7.0) .and. & | |
| 414 | 0 | ( mod(G%geoLonT(i,j)+360.0,360.0) > 106.0) .and. & |
| 415 | ( mod(G%geoLonT(i,j)+360.0,360.0) < 140.0) ) then | |
| 416 | 0 | Kd_int(i,1) = CS%bckgrnd_vdc_Banda |
| 417 | endif | |
| 418 | ||
| 419 | ! South Banda Sea | |
| 420 | if ( (G%geoLatT(i,j) <= -7.0) .and. (G%geoLatT(i,j) > -8.3) .and. & | |
| 421 | 0 | ( mod(G%geoLonT(i,j)+360.0,360.0) > 111.0) .and. & |
| 422 | 0 | ( mod(G%geoLonT(i,j)+360.0,360.0) < 142.0) ) then |
| 423 | 0 | Kd_int(i,1) = CS%bckgrnd_vdc_Banda |
| 424 | endif | |
| 425 | ||
| 426 | enddo | |
| 427 | ! Update interior values of Kd and Kv (uniform profile; no interpolation needed) | |
| 428 | 0 | do K=1,nz+1 ; do i=is,ie |
| 429 | 0 | Kd_int(i,K) = Kd_int(i,1) |
| 430 | 0 | Kv_bkgnd(i,K) = Kd_int(i,1) * CS%prandtl_bkgnd |
| 431 | enddo ; enddo | |
| 432 | 0 | do k=1,nz ; do i=is,ie |
| 433 | 0 | Kd_lay(i,k) = Kd_int(i,1) |
| 434 | enddo ; enddo | |
| 435 | ||
| 436 | else | |
| 437 | ! Set a potentially spatially varying surface value of diffusivity. | |
| 438 | 720 | if (CS%Henyey_IGW_background) then |
| 439 | 720 | I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. |
| 440 | 87120 | do i=is,ie |
| 441 | 86400 | abs_sinlat = abs(sin(G%geoLatT(i,j)*deg_to_rad)) |
| 442 | 86400 | if (abs(G%geoLatT(i,j))>CS%Henyey_max_lat) abs_sinlat = min_sinlat |
| 443 | Kd_sfc(i) = max(CS%Kd_min, CS%Kd * & | |
| 444 | 87120 | ((abs_sinlat * invcosh(CS%N0_2Omega / max(min_sinlat, abs_sinlat))) * I_x30) ) |
| 445 | enddo | |
| 446 | 0 | elseif (CS%Kd_tanh_lat_fn) then |
| 447 | 0 | do i=is,ie |
| 448 | ! The transition latitude and latitude range are hard-scaled here, since | |
| 449 | ! this is not really intended for wide-spread use, but rather for | |
| 450 | ! comparison with CM2M / CM2.1 settings. | |
| 451 | Kd_sfc(i) = max(CS%Kd_min, CS%Kd * (1.0 + & | |
| 452 | 0 | CS%Kd_tanh_lat_scale * 0.5*tanh((abs(G%geoLatT(i,j)) - 35.0)/5.0) )) |
| 453 | enddo | |
| 454 | else ! Use a spatially constant surface value. | |
| 455 | 0 | do i=is,ie |
| 456 | 0 | Kd_sfc(i) = CS%Kd |
| 457 | enddo | |
| 458 | endif | |
| 459 | ||
| 460 | ! Now set background diffusivities based on these surface values, possibly with vertical structure. | |
| 461 | 720 | if ((.not.CS%physical_OBL_scheme) .and. (CS%Kd /= CS%Kd_tot_ml)) then |
| 462 | ! This is a crude way to put in a diffusive boundary layer without an explicit boundary | |
| 463 | ! layer turbulence scheme. It should not be used for any realistic ocean models. | |
| 464 | 0 | I_Hmix = 1.0 / (CS%Hmix + GV%H_subroundoff) |
| 465 | 0 | do i=is,ie ; depth(i) = 0.0 ; enddo |
| 466 | 0 | do k=1,nz ; do i=is,ie |
| 467 | 0 | depth_c = depth(i) + 0.5*h(i,j,k) |
| 468 | 0 | if (depth_c <= CS%Hmix) then ; Kd_lay(i,k) = CS%Kd_tot_ml |
| 469 | 0 | elseif (depth_c >= 2.0*CS%Hmix) then ; Kd_lay(i,k) = Kd_sfc(i) |
| 470 | else | |
| 471 | 0 | Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kd_tot_ml) * I_Hmix) * depth_c + (2.0*CS%Kd_tot_ml - Kd_sfc(i)) |
| 472 | endif | |
| 473 | ||
| 474 | 0 | depth(i) = depth(i) + h(i,j,k) |
| 475 | enddo ; enddo | |
| 476 | ||
| 477 | else ! There is no vertical structure to the background diffusivity. | |
| 478 | 6534720 | do k=1,nz ; do i=is,ie |
| 479 | 6534000 | Kd_lay(i,k) = Kd_sfc(i) |
| 480 | enddo ; enddo | |
| 481 | endif | |
| 482 | ||
| 483 | ! Update Kd_int and Kv_bkgnd, based on Kd_lay. These might be just used for diagnostic purposes. | |
| 484 | 87120 | do i=is,ie |
| 485 | 86400 | Kd_int(i,1) = 0.0 ; Kv_bkgnd(i,1) = 0.0 |
| 486 | 87120 | Kd_int(i,nz+1) = 0.0 ; Kv_bkgnd(i,nz+1) = 0.0 |
| 487 | enddo | |
| 488 | 6447600 | do K=2,nz ; do i=is,ie |
| 489 | 6393600 | Kd_int(i,K) = 0.5*(Kd_lay(i,k-1) + Kd_lay(i,k)) |
| 490 | 6446880 | Kv_bkgnd(i,K) = Kd_int(i,K) * CS%prandtl_bkgnd |
| 491 | enddo ; enddo | |
| 492 | endif | |
| 493 | ||
| 494 | 720 | end subroutine calculate_bkgnd_mixing |
| 495 | ||
| 496 | !> Reads the parameter "USE_CVMix_BACKGROUND" and returns state. | |
| 497 | !! This function allows other modules to know whether this parameterization will | |
| 498 | !! be used without needing to duplicate the log entry. | |
| 499 | 0 | logical function CVMix_bkgnd_is_used(param_file) |
| 500 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 501 | call get_param(param_file, mdl, "USE_CVMix_BACKGROUND", CVMix_bkgnd_is_used, & | |
| 502 | 0 | default=.false., do_not_log=.true.) |
| 503 | ||
| 504 | 0 | end function CVMix_bkgnd_is_used |
| 505 | ||
| 506 | !> Sets CS%bkgnd_scheme_str to check whether multiple background diffusivity schemes are activated. | |
| 507 | !! The string is also for error/log messages. | |
| 508 | 1 | subroutine check_bkgnd_scheme(CS, str) |
| 509 | type(bkgnd_mixing_cs), pointer :: CS !< Control structure | |
| 510 | character(len=*), intent(in) :: str !< Background scheme identifier deducted from MOM_input | |
| 511 | !! parameters | |
| 512 | ||
| 513 | 1 | if (trim(CS%bkgnd_scheme_str)=="none") then |
| 514 | 1 | CS%bkgnd_scheme_str = str |
| 515 | else | |
| 516 | call MOM_error(FATAL, "bkgnd_mixing_init: Cannot activate both "//trim(str)//" and "//& | |
| 517 | 0 | trim(CS%bkgnd_scheme_str)//".") |
| 518 | endif | |
| 519 | ||
| 520 | 1 | end subroutine |
| 521 | ||
| 522 | !> Clear pointers and deallocate memory | |
| 523 | 1 | subroutine bkgnd_mixing_end(CS) |
| 524 | type(bkgnd_mixing_cs), pointer :: CS !< Control structure for this module that | |
| 525 | !! will be deallocated in this subroutine | |
| 526 | ||
| 527 | 1 | if (.not. associated(CS)) return |
| 528 | 1 | deallocate(CS) |
| 529 | ||
| 530 | end subroutine bkgnd_mixing_end | |
| 531 | ||
| 532 | ||
| 533 | 0 | end module MOM_bkgnd_mixing |