21 use cvmix_background,
only : cvmix_init_bkgnd, cvmix_coeffs_bkgnd
23 implicit none ;
private
25 #include <MOM_memory.h>
41 real :: bryan_lewis_c1
43 real :: bryan_lewis_c2
45 real :: bryan_lewis_c3
47 real :: bryan_lewis_c4
51 real :: bckgrnd_vdc_eq
53 real :: bckgrnd_vdc_psim
55 real :: bckgrnd_vdc_banda
65 real :: kd_tanh_lat_scale
71 logical :: kd_tanh_lat_fn
75 logical :: bryan_lewis_diffusivity
77 logical :: horiz_varying_background
79 logical :: henyey_igw_background
83 logical :: henyey_igw_background_new
97 logical :: bulkmixedlayer
101 integer :: id_kd_bkgnd = -1
102 integer :: id_kv_bkgnd = -1
104 real,
allocatable,
dimension(:,:) :: kd_sfc
106 real,
allocatable,
dimension(:,:,:) :: kd_bkgnd
107 real,
allocatable,
dimension(:,:,:) :: kv_bkgnd
109 character(len=40) :: bkgnd_scheme_str =
"none"
113 character(len=40) ::
mdl =
"MOM_bkgnd_mixing"
120 type(time_type),
intent(in) :: time
125 type(
diag_ctrl),
target,
intent(inout) :: diag
131 real :: prandtl_bkgnd_comp
134 #include "version_variable.h"
136 if (
associated(cs))
then
137 call mom_error(warning,
"bkgnd_mixing_init called with an associated "// &
138 "control structure.")
145 "Adding static vertical background mixing coefficients")
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.)
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.)
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)
165 cs%bulkmixedlayer = (gv%nkml > 0)
166 if (cs%bulkmixedlayer)
then
168 call get_param(param_file,
mdl,
"KDML", cs%Kdml, default=-1.)
170 "bkgnd_mixing_init: KDML cannot be set when using"// &
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.)
186 call get_param(param_file,
mdl,
'DEBUG', cs%debug, default=.false., do_not_log=.true.)
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.)
195 if (cs%Bryan_Lewis_diffusivity)
then
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.)
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.)
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.)
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.)
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", &
221 if (cs%horiz_varying_background)
then
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)
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)
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)
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)
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)
246 if (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background)
then
248 prandtl_bkgnd_comp = cs%prandtl_bkgnd
249 if (cs%Kd /= 0.0) prandtl_bkgnd_comp = kv/cs%Kd
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")
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.)
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")
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))
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)
284 "The rotation rate of the earth.", units=
"s-1", &
285 default=7.2921e-5, scale=us%T_to_s)
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.)
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)
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.")
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.
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)
332 integer :: i, j, k, is, ie, js, je
334 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
337 deg_to_rad = atan(1.0)/45.0
341 if (.not. (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background))
then
343 do j=js,je ;
do i=is,ie
344 cs%Kd_sfc(i,j) = cs%Kd
348 if (cs%Henyey_IGW_background)
then
349 i_x30 = 2.0 /
invcosh(cs%N0_2Omega*2.0)
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) )
358 elseif (cs%Kd_tanh_lat_fn)
then
360 do j=js,je ;
do i=is,ie
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) ))
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)
380 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
382 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: n2_lay
384 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: kd_lay
386 real,
dimension(:,:,:),
pointer :: kv
388 integer,
intent(in) :: j
393 real,
dimension(SZK_(G)+1) :: depth_int
394 real,
dimension(SZK_(G)+1) :: kd_col
395 real,
dimension(SZK_(G)+1) :: kv_col
396 real,
dimension(SZI_(G)) :: depth
406 real :: bckgrnd_vdc_psin
407 real :: bckgrnd_vdc_psis
408 integer :: i, k, is, ie, js, je, nz
410 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
413 deg_to_rad = atan(1.0)/45.0
417 if (cs%Bryan_Lewis_diffusivity)
then
422 depth_int(k) = depth_int(k-1) + gv%H_to_m*h(i,j,k-1)
425 call cvmix_init_bkgnd(max_nlev=nz, &
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)
433 kd_col(:) = 0.0 ; kv_col(:) = 0.0
434 call cvmix_coeffs_bkgnd(mdiff_out=kv_col, tdiff_out=kd_col, nlev=nz, max_nlev=nz)
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)
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))
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)
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))
459 depth(i) = depth(i) + gv%H_to_Z*h(i,j,k)
462 elseif (cs%horiz_varying_background)
then
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)
468 cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_eq + bckgrnd_vdc_psin + bckgrnd_vdc_psis
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
475 cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1
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
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
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
500 cs%kv_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) * cs%prandtl_bkgnd
503 kd_lay(i,j,:) = cs%Kd_bkgnd(i,j,1)
507 elseif (cs%Henyey_IGW_background_new)
then
508 i_x30 = 2.0 /
invcosh(cs%N0_2Omega*2.0)
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)
519 do k=1,nz ;
do i=is,ie
520 kd_lay(i,j,k) = cs%Kd_sfc(i,j)
525 if (.not. (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background))
then
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
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
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)
554 default=.false., do_not_log = .true.)
562 character(len=*),
intent(in) :: str
565 if (trim(cs%bkgnd_scheme_str)==
"none")
then
566 cs%bkgnd_scheme_str = str
568 call mom_error(fatal,
"set_diffusivity_init: Cannot activate both "//trim(str)//
" and "//&
569 trim(cs%bkgnd_scheme_str)//
".")
579 if (.not.
associated(cs))
return
581 deallocate(cs%kd_bkgnd)
582 deallocate(cs%kv_bkgnd)