MOM_bkgnd_mixing.F90

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
9
10use mom_debugging, only : hchksum
12use mom_diag_mediator, only : post_data
13use mom_error_handler, only : mom_error, fatal, warning, note
14use mom_file_parser, only : get_param, log_param, log_version, param_file_type
16use mom_forcing_type, only : forcing
17use mom_grid, only : ocean_grid_type
18use mom_interface_heights, only : thickness_to_dz
23use cvmix_background, only : cvmix_init_bkgnd, cvmix_coeffs_bkgnd
24
25implicit none ; private
26
27#include <MOM_memory.h>
28
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.
39type, 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
96end type bkgnd_mixing_cs
97
98character(len=40) :: mdl = "MOM_bkgnd_mixing" !< This module's name.
99
100contains
101
102!> Initialize the background mixing routine.
103subroutine 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 if (associated(cs)) then
129 call mom_error(warning, "bkgnd_mixing_init called with an associated "// &
130 "control structure.")
131 return
132 endif
133 allocate(cs)
134
135 ! Read parameters
136 call log_version(param_file, mdl, version, &
137 "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 "may be used.", default=0.0, units="m2 s-1", scale=us%m2_s_to_Z2_T)
143 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 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 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 cs%physical_OBL_scheme = physical_obl_scheme
157 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 units="m2 s-1", default=-1., scale=gv%m2_s_to_HZ_T, do_not_log=.true.)
161 if (cs%Kd_tot_ml>0.) call mom_error(fatal, &
162 "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 units="m2 s-1", default=-1.0, scale=gv%m2_s_to_HZ_T, do_not_log=.true.)
165 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 "boundary layer mixing parameterization.")
168 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 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 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 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 if (abs(cs%Kd_tot_ml - cs%Kd) > 1.0e-15*abs(cs%Kd)) &
182 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 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 "mixed layer is not used.", units="m", scale=gv%m_to_H, fail_if_missing=.true.)
194 endif
195
196 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 "This is done via CVMix.", default=.false.)
204
205 if (cs%Bryan_Lewis_diffusivity) then
206 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 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 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 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 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 default=.false.)
230
231 if (cs%horiz_varying_background) then
232 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 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 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 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 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 units="nondim", default=1.0)
255
256 if (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background) then
257 prandtl_bkgnd_comp = cs%prandtl_bkgnd
258 if (cs%Kd /= 0.0) prandtl_bkgnd_comp = kv / cs%Kd
259
260 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 "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 "Harrison & Hallberg, JPO 2008.", default=.false.)
270 if (cs%Henyey_IGW_background) call check_bkgnd_scheme(cs, "HENYEY_IGW_BACKGROUND")
271
272 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 "diffusivity (KD) is specified along with "//trim(cs%bkgnd_scheme_str))
276 endif
277
278 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 "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 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 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 "HENYEY_IGW_BACKGROUND.", default=.false.)
297
298 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 "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0)
303
304 if (cs%Henyey_IGW_background .and. cs%Kd_tanh_lat_fn) call mom_error(fatal, &
305 "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.")
306
307! call closeParameterBlock(param_file)
308
309end subroutine bkgnd_mixing_init
310
311!> Calculates the vertical background diffusivities/viscosities
312subroutine 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 real, dimension(SZK_(GV)+1) :: depth_int !< Distance from surface of the interfaces [m]
333 real, dimension(SZK_(GV)+1) :: kd_col !< Diffusivities at the interfaces [m2 s-1]
334 real, dimension(SZK_(GV)+1) :: kv_col !< Viscosities at the interfaces [m2 s-1]
335 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 real, dimension(SZI_(G)) :: depth !< Distance from surface of an interface [H ~> m or kg m-2]
337 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
349
350 ! set some parameters
351 deg_to_rad = atan(1.0)/45.0 ! = PI/180
352 min_sinlat = 1.e-10
353
354 ! Start with a constant value that may be replaced below.
355 kd_lay(:,:) = cs%Kd
356 kv_bkgnd(:,:) = 0.0
357
358 ! Set up the background diffusivity.
359 if (cs%Bryan_Lewis_diffusivity) then
360
361 call thickness_to_dz(h, tv, dz, j, g, gv)
362
363 do i=is,ie
364 depth_int(1) = 0.0
365 do k=2,nz+1
366 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 prandtl = cs%prandtl_bkgnd)
376
377 kd_col(:) = 0.0 ; kv_col(:) = 0.0 ! Is this line necessary?
378 call cvmix_coeffs_bkgnd(mdiff_out=kv_col, tdiff_out=kd_col, nlev=nz, max_nlev=nz)
379
380 ! Update Kd and Kv.
381 do k=1,nz+1
382 kv_bkgnd(i,k) = gv%m2_s_to_HZ_T * kv_col(k)
383 kd_int(i,k) = gv%m2_s_to_HZ_T*kd_col(k)
384 enddo
385 do k=1,nz
386 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 elseif (cs%horiz_varying_background) then
391 !### Note that there are lots of hard-coded parameters (mostly latitudes and longitudes) here.
392 do i=is,ie
393 bckgrnd_vdc_psis = cs%bckgrnd_vdc_psim * exp(-(0.4*(g%geoLatT(i,j)+28.9))**2)
394 bckgrnd_vdc_psin = cs%bckgrnd_vdc_psim * exp(-(0.4*(g%geoLatT(i,j)-28.9))**2)
395 kd_int(i,1) = (cs%bckgrnd_vdc_eq + bckgrnd_vdc_psin) + bckgrnd_vdc_psis
396
397 if (g%geoLatT(i,j) < -10.0) then
398 kd_int(i,1) = kd_int(i,1) + cs%bckgrnd_vdc1
399 elseif (g%geoLatT(i,j) <= 10.0) then
400 kd_int(i,1) = kd_int(i,1) + cs%bckgrnd_vdc1 * (g%geoLatT(i,j)/10.0)**2
401 else
402 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 ( 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 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 ( 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 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 ( mod(g%geoLonT(i,j)+360.0,360.0) > 111.0) .and. &
422 ( mod(g%geoLonT(i,j)+360.0,360.0) < 142.0) ) then
423 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 do k=1,nz+1 ; do i=is,ie
429 kd_int(i,k) = kd_int(i,1)
430 kv_bkgnd(i,k) = kd_int(i,1) * cs%prandtl_bkgnd
431 enddo ; enddo
432 do k=1,nz ; do i=is,ie
433 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 if (cs%Henyey_IGW_background) then
439 i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
440 do i=is,ie
441 abs_sinlat = abs(sin(g%geoLatT(i,j)*deg_to_rad))
442 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 ((abs_sinlat * invcosh(cs%N0_2Omega / max(min_sinlat, abs_sinlat))) * i_x30) )
445 enddo
446 elseif (cs%Kd_tanh_lat_fn) then
447 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 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 do i=is,ie
456 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 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 i_hmix = 1.0 / (cs%Hmix + gv%H_subroundoff)
465 do i=is,ie ; depth(i) = 0.0 ; enddo
466 do k=1,nz ; do i=is,ie
467 depth_c = depth(i) + 0.5*h(i,j,k)
468 if (depth_c <= cs%Hmix) then ; kd_lay(i,k) = cs%Kd_tot_ml
469 elseif (depth_c >= 2.0*cs%Hmix) then ; kd_lay(i,k) = kd_sfc(i)
470 else
471 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 depth(i) = depth(i) + h(i,j,k)
475 enddo ; enddo
476
477 else ! There is no vertical structure to the background diffusivity.
478 do k=1,nz ; do i=is,ie
479 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 do i=is,ie
485 kd_int(i,1) = 0.0 ; kv_bkgnd(i,1) = 0.0
486 kd_int(i,nz+1) = 0.0 ; kv_bkgnd(i,nz+1) = 0.0
487 enddo
488 do k=2,nz ; do i=is,ie
489 kd_int(i,k) = 0.5*(kd_lay(i,k-1) + kd_lay(i,k))
490 kv_bkgnd(i,k) = kd_int(i,k) * cs%prandtl_bkgnd
491 enddo ; enddo
492 endif
493
494end 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.
499logical 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 default=.false., do_not_log=.true.)
503
504end 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.
508subroutine 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 if (trim(cs%bkgnd_scheme_str)=="none") then
514 cs%bkgnd_scheme_str = str
515 else
516 call mom_error(fatal, "bkgnd_mixing_init: Cannot activate both "//trim(str)//" and "//&
517 trim(cs%bkgnd_scheme_str)//".")
518 endif
519
520end subroutine
521
522!> Clear pointers and deallocate memory
523subroutine 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 if (.not. associated(cs)) return
528 deallocate(cs)
529
530end subroutine bkgnd_mixing_end
531
532
533end module mom_bkgnd_mixing