MOM6
MOM_bkgnd_mixing.F90
Go to the documentation of this file.
1 !> Interface to background mixing schemes, including the Bryan and Lewis (1979)
2 !! which is applied via CVMix.
3 
5 
6 ! This file is part of MOM6. See LICENSE.md for the license.
7 
8 use mom_debugging, only : hchksum
9 use mom_diag_mediator, only : diag_ctrl, time_type, register_diag_field
10 use mom_diag_mediator, only : post_data
12 use mom_error_handler, only : mom_error, fatal, warning, note
15 use mom_forcing_type, only : forcing
16 use mom_grid, only : ocean_grid_type
21 use cvmix_background, only : cvmix_init_bkgnd, cvmix_coeffs_bkgnd
22 
23 implicit none ; private
24 
25 #include <MOM_memory.h>
26 
27 public bkgnd_mixing_init
28 public bkgnd_mixing_end
30 public sfc_bkgnd_mixing
31 
32 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
33 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
34 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
35 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
36 
37 !> Control structure including parameters for this module.
38 type, public :: bkgnd_mixing_cs ! TODO: private
39 
40  ! Parameters
41  real :: bryan_lewis_c1 !< The vertical diffusivity values for Bryan-Lewis profile
42  !! at |z|=D [m2 s-1]
43  real :: bryan_lewis_c2 !< The amplitude of variation in diffusivity for the
44  !! Bryan-Lewis diffusivity profile [m2 s-1]
45  real :: bryan_lewis_c3 !< The inverse length scale for transition region in the
46  !! Bryan-Lewis diffusivity profile [m-1]
47  real :: bryan_lewis_c4 !< The depth where diffusivity is Bryan_Lewis_bl1 in the
48  !! Bryan-Lewis profile [m]
49  real :: bckgrnd_vdc1 !< Background diffusivity (Ledwell) when
50  !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
51  real :: bckgrnd_vdc_eq !< Equatorial diffusivity (Gregg) when
52  !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
53  real :: bckgrnd_vdc_psim !< Max. PSI induced diffusivity (MacKinnon) when
54  !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
55  real :: bckgrnd_vdc_banda !< Banda Sea diffusivity (Gordon) when
56  !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
57  real :: kd_min !< minimum diapycnal diffusivity [Z2 T-1 ~> m2 s-1]
58  real :: kd !< interior diapycnal diffusivity [Z2 T-1 ~> m2 s-1]
59  real :: omega !< The Earth's rotation rate [T-1 ~> s-1].
60  real :: n0_2omega !< ratio of the typical Buoyancy frequency to
61  !! twice the Earth's rotation period, used with the
62  !! Henyey scaling from the mixing
63  real :: prandtl_bkgnd !< Turbulent Prandtl number used to convert
64  !! vertical background diffusivity into viscosity
65  real :: kd_tanh_lat_scale !< A nondimensional scaling for the range of
66  !! diffusivities with Kd_tanh_lat_fn. Valid values
67  !! are in the range of -2 to 2; 0.4 reproduces CM2M.
68  real :: kdml !< mixed layer diapycnal diffusivity [Z2 T-1 ~> m2 s-1]
69  !! when bulkmixedlayer==.false.
70  real :: hmix !< mixed layer thickness [Z ~> m] when bulkmixedlayer==.false.
71  logical :: kd_tanh_lat_fn !< If true, use the tanh dependence of Kd_sfc on
72  !! latitude, like GFDL CM2.1/CM2M. There is no
73  !! physical justification for this form, and it can
74  !! not be used with Henyey_IGW_background.
75  logical :: bryan_lewis_diffusivity!< If true, background vertical diffusivity
76  !! uses Bryan-Lewis (1979) like tanh profile.
77  logical :: horiz_varying_background !< If true, apply vertically uniform, latitude-dependent
78  !! background diffusivity, as described in Danabasoglu et al., 2012
79  logical :: henyey_igw_background !< If true, use a simplified variant of the
80  !! Henyey et al, JGR (1986) latitudinal scaling for the background diapycnal diffusivity,
81  !! which gives a marked decrease in the diffusivity near the equator. The simplification
82  !! here is to assume that the in-situ stratification is the same as the reference stratificaiton.
83  logical :: henyey_igw_background_new !< same as Henyey_IGW_background
84  !! but incorporate the effect of stratification on TKE dissipation,
85  !! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0
86  !! where e is the TKE dissipation, and N_0 and f_0
87  !! are the reference buoyancy frequency and inertial frequencies respectively.
88  !! e_0 is the reference dissipation at (N_0,f_0). In the previous version, N=N_0.
89  !! Additionally, the squared inverse relationship between diapycnal diffusivities
90  !! and stratification is included:
91  !!
92  !! kd = e/N^2
93  !!
94  !! where kd is the diapycnal diffusivity. This approach assumes that work done
95  !! against gravity is uniformly distributed throughout the column. Whereas, kd=kd_0*e,
96  !! as in the original version, concentrates buoyancy work in regions of strong stratification.
97  logical :: bulkmixedlayer !< If true, a refined bulk mixed layer scheme is used
98  logical :: debug !< If true, turn on debugging in this module
99  ! Daignostic handles and pointers
100  type(diag_ctrl), pointer :: diag => null() !< A structure that regulates diagnostic output
101  integer :: id_kd_bkgnd = -1 !< Diagnotic IDs
102  integer :: id_kv_bkgnd = -1 !< Diagnostic IDs
103 
104  real, allocatable, dimension(:,:) :: kd_sfc !< surface value of the diffusivity [Z2 T-1 ~> m2 s-1]
105  ! Diagnostics arrays
106  real, allocatable, dimension(:,:,:) :: kd_bkgnd !< Background diffusivity [Z2 T-1 ~> m2 s-1]
107  real, allocatable, dimension(:,:,:) :: kv_bkgnd !< Background viscosity [Z2 s-1 ~> m2 s-1]
108 
109  character(len=40) :: bkgnd_scheme_str = "none" !< Background scheme identifier
110 
111 end type bkgnd_mixing_cs
112 
113 character(len=40) :: mdl = "MOM_bkgnd_mixing" !< This module's name.
114 
115 contains
116 
117 !> Initialize the background mixing routine.
118 subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS)
119 
120  type(time_type), intent(in) :: time !< The current time.
121  type(ocean_grid_type), intent(in) :: g !< Grid structure.
122  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
123  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
124  type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
125  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
126  type(bkgnd_mixing_cs), pointer :: cs !< This module's control structure.
127 
128  ! Local variables
129  real :: kv ! The interior vertical viscosity [Z2 T-1 ~> m2 s-1] - read to set prandtl
130  ! number unless it is provided as a parameter
131  real :: prandtl_bkgnd_comp ! Kv/CS%Kd. Gets compared with user-specified prandtl_bkgnd.
132 
133 ! This include declares and sets the variable "version".
134 #include "version_variable.h"
135 
136  if (associated(cs)) then
137  call mom_error(warning, "bkgnd_mixing_init called with an associated "// &
138  "control structure.")
139  return
140  endif
141  allocate(cs)
142 
143  ! Read parameters
144  call log_version(param_file, mdl, version, &
145  "Adding static vertical background mixing coefficients")
146 
147  call get_param(param_file, mdl, "KD", cs%Kd, &
148  "The background diapycnal diffusivity of density in the "//&
149  "interior. Zero or the molecular value, ~1e-7 m2 s-1, "//&
150  "may be used.", units="m2 s-1", scale=us%m2_s_to_Z2_T, fail_if_missing=.true.)
151 
152  call get_param(param_file, mdl, "KV", kv, &
153  "The background kinematic viscosity in the interior. "//&
154  "The molecular value, ~1e-6 m2 s-1, may be used.", &
155  units="m2 s-1", scale=us%m2_s_to_Z2_T, fail_if_missing=.true.)
156 
157  call get_param(param_file, mdl, "KD_MIN", cs%Kd_min, &
158  "The minimum diapycnal diffusivity.", &
159  units="m2 s-1", default=0.01*cs%Kd*us%Z2_T_to_m2_s, scale=us%m2_s_to_Z2_T)
160 
161  ! The following is needed to set one of the choices of vertical background mixing
162 
163  ! BULKMIXEDLAYER is not always defined (e.g., CM2G63L), so the following by pass
164  ! the need to include BULKMIXEDLAYER in MOM_input
165  cs%bulkmixedlayer = (gv%nkml > 0)
166  if (cs%bulkmixedlayer) then
167  ! Check that Kdml is not set when using bulk mixed layer
168  call get_param(param_file, mdl, "KDML", cs%Kdml, default=-1.)
169  if (cs%Kdml>0.) call mom_error(fatal, &
170  "bkgnd_mixing_init: KDML cannot be set when using"// &
171  "bulk mixed layer.")
172  cs%Kdml = cs%Kd ! This is not used with a bulk mixed layer, but also
173  ! cannot be a NaN.
174  else
175  call get_param(param_file, mdl, "KDML", cs%Kdml, &
176  "If BULKMIXEDLAYER is false, KDML is the elevated "//&
177  "diapycnal diffusivity in the topmost HMIX of fluid. "//&
178  "KDML is only used if BULKMIXEDLAYER is false.", &
179  units="m2 s-1", default=cs%Kd*us%Z2_T_to_m2_s, scale=us%m2_s_to_Z2_T)
180  call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &
181  "The prescribed depth over which the near-surface "//&
182  "viscosity and diffusivity are elevated when the bulk "//&
183  "mixed layer is not used.", units="m", scale=us%m_to_Z, fail_if_missing=.true.)
184  endif
185 
186  call get_param(param_file, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
187 
188 ! call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING')
189 
190  call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", cs%Bryan_Lewis_diffusivity, &
191  "If true, use a Bryan & Lewis (JGR 1979) like tanh "//&
192  "profile of background diapycnal diffusivity with depth. "//&
193  "This is done via CVMix.", default=.false.)
194 
195  if (cs%Bryan_Lewis_diffusivity) then
196  call check_bkgnd_scheme(cs, "BRYAN_LEWIS_DIFFUSIVITY")
197 
198  call get_param(param_file, mdl, "BRYAN_LEWIS_C1", cs%Bryan_Lewis_c1, &
199  "The vertical diffusivity values for Bryan-Lewis profile at |z|=D.", &
200  units="m2 s-1", fail_if_missing=.true.)
201 
202  call get_param(param_file, mdl, "BRYAN_LEWIS_C2", cs%Bryan_Lewis_c2, &
203  "The amplitude of variation in diffusivity for the Bryan-Lewis profile", &
204  units="m2 s-1", fail_if_missing=.true.)
205 
206  call get_param(param_file, mdl, "BRYAN_LEWIS_C3", cs%Bryan_Lewis_c3, &
207  "The inverse length scale for transition region in the Bryan-Lewis profile", &
208  units="m-1", fail_if_missing=.true.)
209 
210  call get_param(param_file, mdl, "BRYAN_LEWIS_C4", cs%Bryan_Lewis_c4, &
211  "The depth where diffusivity is BRYAN_LEWIS_C1 in the Bryan-Lewis profile",&
212  units="m", fail_if_missing=.true.)
213 
214  endif ! CS%Bryan_Lewis_diffusivity
215 
216  call get_param(param_file, mdl, "HORIZ_VARYING_BACKGROUND", cs%horiz_varying_background, &
217  "If true, apply vertically uniform, latitude-dependent background "//&
218  "diffusivity, as described in Danabasoglu et al., 2012", &
219  default=.false.)
220 
221  if (cs%horiz_varying_background) then
222  call check_bkgnd_scheme(cs, "HORIZ_VARYING_BACKGROUND")
223 
224  call get_param(param_file, mdl, "BCKGRND_VDC1", cs%bckgrnd_vdc1, &
225  "Background diffusivity (Ledwell) when HORIZ_VARYING_BACKGROUND=True", &
226  units="m2 s-1",default = 0.16e-04, scale=us%m2_s_to_Z2_T)
227 
228  call get_param(param_file, mdl, "BCKGRND_VDC_EQ", cs%bckgrnd_vdc_eq, &
229  "Equatorial diffusivity (Gregg) when HORIZ_VARYING_BACKGROUND=True", &
230  units="m2 s-1",default = 0.01e-04, scale=us%m2_s_to_Z2_T)
231 
232  call get_param(param_file, mdl, "BCKGRND_VDC_PSIM", cs%bckgrnd_vdc_psim, &
233  "Max. PSI induced diffusivity (MacKinnon) when HORIZ_VARYING_BACKGROUND=True", &
234  units="m2 s-1",default = 0.13e-4, scale=us%m2_s_to_Z2_T)
235 
236  call get_param(param_file, mdl, "BCKGRND_VDC_BAN", cs%bckgrnd_vdc_Banda, &
237  "Banda Sea diffusivity (Gordon) when HORIZ_VARYING_BACKGROUND=True", &
238  units="m2 s-1",default = 1.0e-4, scale=us%m2_s_to_Z2_T)
239  endif
240 
241  call get_param(param_file, mdl, "PRANDTL_BKGND", cs%prandtl_bkgnd, &
242  "Turbulent Prandtl number used to convert vertical "//&
243  "background diffusivities into viscosities.", &
244  units="nondim", default=1.0)
245 
246  if (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background) then
247 
248  prandtl_bkgnd_comp = cs%prandtl_bkgnd
249  if (cs%Kd /= 0.0) prandtl_bkgnd_comp = kv/cs%Kd
250 
251  if ( abs(cs%prandtl_bkgnd - prandtl_bkgnd_comp)>1.e-14) then
252  call mom_error(fatal,"set_diffusivity_init: The provided KD, KV,"//&
253  "and PRANDTL_BKGND values are incompatible. The following "//&
254  "must hold: KD*PRANDTL_BKGND==KV")
255  endif
256 
257  endif
258 
259  call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", cs%Henyey_IGW_background, &
260  "If true, use a latitude-dependent scaling for the near "//&
261  "surface background diffusivity, as described in "//&
262  "Harrison & Hallberg, JPO 2008.", default=.false.)
263  if (cs%Henyey_IGW_background) call check_bkgnd_scheme(cs, "HENYEY_IGW_BACKGROUND")
264 
265 
266  call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", cs%Henyey_IGW_background_new, &
267  "If true, use a better latitude-dependent scaling for the "//&
268  "background diffusivity, as described in "//&
269  "Harrison & Hallberg, JPO 2008.", default=.false.)
270  if (cs%Henyey_IGW_background_new) call check_bkgnd_scheme(cs, "HENYEY_IGW_BACKGROUND_NEW")
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, "set_diffusivity_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.", units="s-1", &
285  default=7.2921e-5, scale=us%T_to_s)
286  endif
287 
288  call get_param(param_file, mdl, "KD_TANH_LAT_FN", &
289  cs%Kd_tanh_lat_fn, &
290  "If true, use a tanh dependence of Kd_sfc on latitude, "//&
291  "like CM2.1/CM2M. There is no physical justification "//&
292  "for this form, and it can not be used with "//&
293  "HENYEY_IGW_BACKGROUND.", default=.false.)
294 
295  if (cs%Kd_tanh_lat_fn) &
296  call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", cs%Kd_tanh_lat_scale, &
297  "A nondimensional scaling for the range ofdiffusivities "//&
298  "with KD_TANH_LAT_FN. Valid values are in the range of "//&
299  "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0)
300 
301  if (cs%Henyey_IGW_background .and. cs%Kd_tanh_lat_fn) call mom_error(fatal, &
302  "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.")
303 
304 ! call closeParameterBlock(param_file)
305 
306  ! allocate arrays and set them to zero
307  allocate(cs%Kd_bkgnd(szi_(g), szj_(g), szk_(g)+1)); cs%kd_bkgnd(:,:,:) = 0.
308  allocate(cs%kv_bkgnd(szi_(g), szj_(g), szk_(g)+1)); cs%kv_bkgnd(:,:,:) = 0.
309  allocate(cs%Kd_sfc(szi_(g), szj_(g))); cs%Kd_sfc(:,:) = 0.
310 
311  ! Register diagnostics
312  cs%diag => diag
313  cs%id_kd_bkgnd = register_diag_field('ocean_model', 'Kd_bkgnd', diag%axesTi, time, &
314  'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s', conversion=us%Z2_T_to_m2_s)
315  cs%id_kv_bkgnd = register_diag_field('ocean_model', 'Kv_bkgnd', diag%axesTi, time, &
316  'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s', conversion=us%Z2_T_to_m2_s)
317 
318 end subroutine bkgnd_mixing_init
319 
320 !> Get surface vertical background diffusivities/viscosities.
321 subroutine sfc_bkgnd_mixing(G, US, CS)
322 
323  type(ocean_grid_type), intent(in) :: g !< Grid structure.
324  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
325  type(bkgnd_mixing_cs), pointer, intent(inout) :: cs !< The control structure returned by
326  !! a previous call to bkgnd_mixing_init.
327  ! local variables
328  real :: i_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg)))
329  real :: deg_to_rad !< factor converting degrees to radians, pi/180.
330  real :: abs_sin !< absolute value of sine of latitude [nondim]
331  real :: epsilon
332  integer :: i, j, k, is, ie, js, je
333 
334  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
335 
336  ! set some parameters
337  deg_to_rad = atan(1.0)/45.0 ! = PI/180
338  epsilon = 1.e-10
339 
340 
341  if (.not. (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background)) then
342 !$OMP parallel do default(none) shared(is,ie,js,je,CS)
343  do j=js,je ; do i=is,ie
344  cs%Kd_sfc(i,j) = cs%Kd
345  enddo ; enddo
346  endif
347 
348  if (cs%Henyey_IGW_background) then
349  i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
350 !$OMP parallel do default(none) &
351 !$OMP shared(is,ie,js,je,CS,G,deg_to_rad,epsilon,I_x30) &
352 !$OMP private(abs_sin)
353  do j=js,je ; do i=is,ie
354  abs_sin = abs(sin(g%geoLatT(i,j)*deg_to_rad))
355  cs%Kd_sfc(i,j) = max(cs%Kd_min, cs%Kd_sfc(i,j) * &
356  ((abs_sin * invcosh(cs%N0_2Omega/max(epsilon,abs_sin))) * i_x30) )
357  enddo ; enddo
358  elseif (cs%Kd_tanh_lat_fn) then
359 !$OMP parallel do default(none) shared(is,ie,js,je,CS,G)
360  do j=js,je ; do i=is,ie
361  ! The transition latitude and latitude range are hard-scaled here, since
362  ! this is not really intended for wide-spread use, but rather for
363  ! comparison with CM2M / CM2.1 settings.
364  cs%Kd_sfc(i,j) = max(cs%Kd_min, cs%Kd_sfc(i,j) * (1.0 + &
365  cs%Kd_tanh_lat_scale * 0.5*tanh((abs(g%geoLatT(i,j)) - 35.0)/5.0) ))
366  enddo ; enddo
367  endif
368 
369  if (cs%debug) call hchksum(cs%Kd_sfc,"After sfc_bkgnd_mixing: Kd_sfc",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
370 
371 end subroutine sfc_bkgnd_mixing
372 
373 
374 !> Calculates the vertical background diffusivities/viscosities
375 subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kv, j, G, GV, US, CS)
376 
377  type(ocean_grid_type), intent(in) :: g !< Grid structure.
378  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
379  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
380  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
381  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
382  real, dimension(SZI_(G),SZK_(G)), intent(in) :: n2_lay !< squared buoyancy frequency associated
383  !! with layers [T-2 ~> s-2]
384  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: kd_lay !< Diapycnal diffusivity of each layer
385  !! [Z2 T-1 ~> m2 s-1].
386  real, dimension(:,:,:), pointer :: kv !< The "slow" vertical viscosity at each interface
387  !! (not layer!) [Z2 T-1 ~> m2 s-1]
388  integer, intent(in) :: j !< Meridional grid index
389  type(bkgnd_mixing_cs), pointer :: cs !< The control structure returned by
390  !! a previous call to bkgnd_mixing_init.
391 
392  ! local variables
393  real, dimension(SZK_(G)+1) :: depth_int !< distance from surface of the interfaces [m]
394  real, dimension(SZK_(G)+1) :: kd_col !< Diffusivities at the interfaces [m2 s-1]
395  real, dimension(SZK_(G)+1) :: kv_col !< Viscosities at the interfaces [m2 s-1]
396  real, dimension(SZI_(G)) :: depth !< distance from surface of an interface [Z ~> m]
397  real :: depth_c !< depth of the center of a layer [Z ~> m]
398  real :: i_hmix !< inverse of fixed mixed layer thickness [Z-1 ~> m-1]
399  real :: i_2omega !< 1/(2 Omega) [T ~> s]
400  real :: n_2omega ! The ratio of the stratification to the Earth's rotation rate [nondim]
401  real :: n02_n2 ! The ratio a reference stratification to the actual stratification [nondim]
402  real :: i_x30 !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg)))
403  real :: deg_to_rad !< factor converting degrees to radians, pi/180.
404  real :: abs_sin !< absolute value of sine of latitude [nondim]
405  real :: epsilon ! The minimum value of the sine of latitude [nondim]
406  real :: bckgrnd_vdc_psin !< PSI diffusivity in northern hemisphere [Z2 T-1 ~> m2 s-1]
407  real :: bckgrnd_vdc_psis !< PSI diffusivity in southern hemisphere [Z2 T-1 ~> m2 s-1]
408  integer :: i, k, is, ie, js, je, nz
409 
410  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
411 
412  ! set some parameters
413  deg_to_rad = atan(1.0)/45.0 ! = PI/180
414  epsilon = 1.e-10
415 
416  ! Set up the background diffusivity.
417  if (cs%Bryan_Lewis_diffusivity) then
418 
419  do i=is,ie
420  depth_int(1) = 0.0
421  do k=2,nz+1
422  depth_int(k) = depth_int(k-1) + gv%H_to_m*h(i,j,k-1)
423  enddo
424 
425  call cvmix_init_bkgnd(max_nlev=nz, &
426  zw = depth_int(:), & !< interface depths relative to the surface in m, must be positive.
427  bl1 = cs%Bryan_Lewis_c1, &
428  bl2 = cs%Bryan_Lewis_c2, &
429  bl3 = cs%Bryan_Lewis_c3, &
430  bl4 = cs%Bryan_Lewis_c4, &
431  prandtl = cs%prandtl_bkgnd)
432 
433  kd_col(:) = 0.0 ; kv_col(:) = 0.0 ! Is this line necessary?
434  call cvmix_coeffs_bkgnd(mdiff_out=kv_col, tdiff_out=kd_col, nlev=nz, max_nlev=nz)
435 
436  ! Update Kd and Kv.
437  do k=1,nz+1
438  cs%Kv_bkgnd(i,j,k) = us%m2_s_to_Z2_T*kv_col(k)
439  cs%Kd_bkgnd(i,j,k) = us%m2_s_to_Z2_T*kd_col(k)
440  enddo
441  do k=1,nz
442  kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * us%m2_s_to_Z2_T * (kd_col(k) + kd_col(k+1))
443  enddo
444  enddo ! i loop
445 
446  elseif ((.not. cs%Bryan_Lewis_diffusivity) .and. (.not.cs%bulkmixedlayer) .and. &
447  (.not. cs%horiz_varying_background) .and. (cs%Kd /= cs%Kdml)) then
448  i_hmix = 1.0 / cs%Hmix
449  do i=is,ie ; depth(i) = 0.0 ; enddo
450  do k=1,nz ; do i=is,ie
451  depth_c = depth(i) + 0.5*gv%H_to_Z*h(i,j,k)
452  if (depth_c <= cs%Hmix) then ; cs%Kd_bkgnd(i,j,k) = cs%Kdml
453  elseif (depth_c >= 2.0*cs%Hmix) then ; cs%Kd_bkgnd(i,j,k) = cs%Kd_sfc(i,j)
454  else
455  kd_lay(i,j,k) = ((cs%Kd_sfc(i,j) - cs%Kdml) * i_hmix) * depth_c + &
456  (2.0*cs%Kdml - cs%Kd_sfc(i,j))
457  endif
458 
459  depth(i) = depth(i) + gv%H_to_Z*h(i,j,k)
460  enddo ; enddo
461 
462  elseif (cs%horiz_varying_background) then
463  !### Note that there are lots of hard-coded parameters (mostly latitudes and longitudes) here.
464  do i=is,ie
465  bckgrnd_vdc_psis = cs%bckgrnd_vdc_psim * exp(-(0.4*(g%geoLatT(i,j)+28.9))**2)
466  bckgrnd_vdc_psin = cs%bckgrnd_vdc_psim * exp(-(0.4*(g%geoLatT(i,j)-28.9))**2)
467  !### Add parentheses.
468  cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_eq + bckgrnd_vdc_psin + bckgrnd_vdc_psis
469 
470  if (g%geoLatT(i,j) < -10.0) then
471  cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1
472  elseif (g%geoLatT(i,j) <= 10.0) then
473  cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1 * (g%geoLatT(i,j)/10.0)**2
474  else
475  cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1
476  endif
477 
478  ! North Banda Sea
479  if ( (g%geoLatT(i,j) < -1.0) .and. (g%geoLatT(i,j) > -4.0) .and. &
480  ( mod(g%geoLonT(i,j)+360.0,360.0) > 103.0) .and. &
481  ( mod(g%geoLonT(i,j)+360.0,360.0) < 134.0) ) then
482  cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_Banda
483  endif
484 
485  ! Middle Banda Sea
486  if ( (g%geoLatT(i,j) <= -4.0) .and. (g%geoLatT(i,j) > -7.0) .and. &
487  ( mod(g%geoLonT(i,j)+360.0,360.0) > 106.0) .and. &
488  ( mod(g%geoLonT(i,j)+360.0,360.0) < 140.0) ) then
489  cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_Banda
490  endif
491 
492  ! South Banda Sea
493  if ( (g%geoLatT(i,j) <= -7.0) .and. (g%geoLatT(i,j) > -8.3) .and. &
494  ( mod(g%geoLonT(i,j)+360.0,360.0) > 111.0) .and. &
495  ( mod(g%geoLonT(i,j)+360.0,360.0) < 142.0) ) then
496  cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_Banda
497  endif
498 
499  ! Compute kv_bkgnd
500  cs%kv_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) * cs%prandtl_bkgnd
501 
502  ! Update Kd (uniform profile; no interpolation needed)
503  kd_lay(i,j,:) = cs%Kd_bkgnd(i,j,1)
504 
505  enddo
506 
507  elseif (cs%Henyey_IGW_background_new) then
508  i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
509  i_2omega = 0.5 / cs%omega
510  do k=1,nz ; do i=is,ie
511  abs_sin = max(epsilon, abs(sin(g%geoLatT(i,j)*deg_to_rad)))
512  n_2omega = max(abs_sin, sqrt(n2_lay(i,k))*i_2omega)
513  n02_n2 = (cs%N0_2Omega/n_2omega)**2
514  kd_lay(i,j,k) = max(cs%Kd_min, cs%Kd_sfc(i,j) * &
515  ((abs_sin * invcosh(n_2omega/abs_sin)) * i_x30)*n02_n2)
516  enddo ; enddo
517 
518  else
519  do k=1,nz ; do i=is,ie
520  kd_lay(i,j,k) = cs%Kd_sfc(i,j)
521  enddo ; enddo
522  endif
523 
524  ! Update CS%kd_bkgnd and CS%kv_bkgnd for diagnostic purposes
525  if (.not. (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background)) then
526  do i=is,ie
527  cs%kd_bkgnd(i,j,1) = 0.0; cs%kv_bkgnd(i,j,1) = 0.0
528  cs%kd_bkgnd(i,j,nz+1) = 0.0; cs%kv_bkgnd(i,j,nz+1) = 0.0
529  do k=2,nz
530  cs%Kd_bkgnd(i,j,k) = 0.5*(kd_lay(i,j,k-1) + kd_lay(i,j,k))
531  cs%Kv_bkgnd(i,j,k) = cs%Kd_bkgnd(i,j,k) * cs%prandtl_bkgnd
532  enddo
533  enddo
534  endif
535 
536  ! Update Kv
537  if (associated(kv)) then
538  do k=1,nz+1 ; do i=is,ie
539  kv(i,j,k) = kv(i,j,k) + cs%Kv_bkgnd(i,j,k)
540  enddo ; enddo
541  endif
542 
543  ! TODO: In both CS%Bryan_Lewis_diffusivity and CS%horiz_varying_background, KV and KD at surface
544  ! and bottom interfaces are set to be nonzero. Make sure this is not problematic.
545 
546 end subroutine calculate_bkgnd_mixing
547 
548 !> Reads the parameter "USE_CVMix_BACKGROUND" and returns state.
549 !! This function allows other modules to know whether this parameterization will
550 !! be used without needing to duplicate the log entry.
551 logical function cvmix_bkgnd_is_used(param_file)
552  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
553  call get_param(param_file, mdl, "USE_CVMix_BACKGROUND", cvmix_bkgnd_is_used, &
554  default=.false., do_not_log = .true.)
555 
556 end function cvmix_bkgnd_is_used
557 
558 !> Sets CS%bkgnd_scheme_str to check whether multiple background diffusivity schemes are activated.
559 !! The string is also for error/log messages.
560 subroutine check_bkgnd_scheme(CS,str)
561  type(bkgnd_mixing_cs), pointer :: CS !< Control structure
562  character(len=*), intent(in) :: str !< Background scheme identifier deducted from MOM_input
563  !! parameters
564 
565  if (trim(cs%bkgnd_scheme_str)=="none") then
566  cs%bkgnd_scheme_str = str
567  else
568  call mom_error(fatal, "set_diffusivity_init: Cannot activate both "//trim(str)//" and "//&
569  trim(cs%bkgnd_scheme_str)//".")
570  endif
571 
572 end subroutine
573 
574 !> Clear pointers and dealocate memory
575 subroutine bkgnd_mixing_end(CS)
576  type(bkgnd_mixing_cs), pointer :: cs !< Control structure for this module that
577  !! will be deallocated in this subroutine
578 
579  if (.not. associated(cs)) return
580 
581  deallocate(cs%kd_bkgnd)
582  deallocate(cs%kv_bkgnd)
583  deallocate(cs)
584 
585 end subroutine bkgnd_mixing_end
586 
587 
588 end module mom_bkgnd_mixing
mom_intrinsic_functions::invcosh
real function, public invcosh(x)
Evaluate the inverse cosh, either using a math library or an equivalent expression.
Definition: MOM_intrinsic_functions.F90:16
mom_bkgnd_mixing::bkgnd_mixing_cs
Control structure including parameters for this module.
Definition: MOM_bkgnd_mixing.F90:38
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_bkgnd_mixing::calculate_bkgnd_mixing
subroutine, public calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kv, j, G, GV, US, CS)
Calculates the vertical background diffusivities/viscosities.
Definition: MOM_bkgnd_mixing.F90:376
mom_bkgnd_mixing::bkgnd_mixing_init
subroutine, public bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS)
Initialize the background mixing routine.
Definition: MOM_bkgnd_mixing.F90:119
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_bkgnd_mixing::check_bkgnd_scheme
subroutine check_bkgnd_scheme(CS, str)
Sets CSbkgnd_scheme_str to check whether multiple background diffusivity schemes are activated....
Definition: MOM_bkgnd_mixing.F90:561
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_bkgnd_mixing::bkgnd_mixing_end
subroutine, public bkgnd_mixing_end(CS)
Clear pointers and dealocate memory.
Definition: MOM_bkgnd_mixing.F90:576
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_intrinsic_functions
A module with intrinsic functions that are used by MOM but are not supported by some compilers.
Definition: MOM_intrinsic_functions.F90:3
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_variables::vertvisc_type
Vertical viscosities, drag coefficients, and related fields.
Definition: MOM_variables.F90:196
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_bkgnd_mixing
Interface to background mixing schemes, including the Bryan and Lewis (1979) which is applied via CVM...
Definition: MOM_bkgnd_mixing.F90:4
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_forcing_type::forcing
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
Definition: MOM_forcing_type.F90:50
mom_file_parser::openparameterblock
subroutine, public openparameterblock(CS, blockName, desc)
Tags blockName onto the end of the active parameter block name.
Definition: MOM_file_parser.F90:2015
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_bkgnd_mixing::mdl
character(len=40) mdl
This module's name.
Definition: MOM_bkgnd_mixing.F90:113
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_bkgnd_mixing::sfc_bkgnd_mixing
subroutine, public sfc_bkgnd_mixing(G, US, CS)
Get surface vertical background diffusivities/viscosities.
Definition: MOM_bkgnd_mixing.F90:322
mom_file_parser::closeparameterblock
subroutine, public closeparameterblock(CS)
Remove the lowest level of recursion from the active block name.
Definition: MOM_file_parser.F90:2033
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_bkgnd_mixing::cvmix_bkgnd_is_used
logical function cvmix_bkgnd_is_used(param_file)
Reads the parameter "USE_CVMix_BACKGROUND" and returns state. This function allows other modules to k...
Definition: MOM_bkgnd_mixing.F90:552
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60