MOM6
Data Types | Functions/Subroutines
mom_kappa_shear Module Reference

Detailed Description

Shear-dependent mixing following Jackson et al. 2008.

By Laura Jackson and Robert Hallberg, 2006-2008.

This file contains the subroutines that determine the diapycnal diffusivity driven by resolved shears, as specified by the parameterizations described in Jackson and Hallberg (JPO, 2008).

The technique by which the 6 equations (for kappa, TKE, u, v, T, and S) are solved simultaneously has been dramatically revised from the previous version. The previous version was not converging in some cases, especially near the surface mixed layer, while the revised version does. The revised version solves for kappa and TKE with shear and stratification fixed, then marches the density and velocities forward with an adaptive (and aggressive) time step in a predictor-corrector-corrector emulation of a trapezoidal scheme. Run-time-settable parameters determine the tolerence to which the kappa and TKE equations are solved and the minimum time step that can be taken.

Data Types

type  kappa_shear_cs
 This control structure holds the parameters that regulate shear mixing. More...
 

Functions/Subroutines

subroutine, public calculate_kappa_shear (u_in, v_in, h, tv, p_surf, kappa_io, tke_io, kv_io, dt, G, GV, US, CS, initialize_all)
 Subroutine for calculating shear-driven diffusivity and TKE in tracer columns. More...
 
subroutine, public calc_kappa_shear_vertex (u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_io, tke_io, kv_io, dt, G, GV, US, CS, initialize_all)
 Subroutine for calculating shear-driven diffusivity and TKE in corner columns. More...
 
subroutine kappa_shear_column (kappa, tke, dt, nzc, f2, surface_pres, dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, tv, CS, GV, US, I_Ld2_1d, dz_Int_1d)
 This subroutine calculates shear-driven diffusivity and TKE in a single column. More...
 
subroutine calculate_projected_state (kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int, dbuoy_dT, dbuoy_dS, u, v, T, Sal, GV, US, N2, S2, ks_int, ke_int, vel_underflow)
 This subroutine calculates the velocities, temperature and salinity that the water column will have after mixing for dt with diffusivities kappa. It may also calculate the projected buoyancy frequency and shear. More...
 
subroutine find_kappa_tke (N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, nz, CS, GV, US, K_Q, tke, kappa, kappa_src, local_src)
 This subroutine calculates new, consistent estimates of TKE and kappa. More...
 
logical function, public kappa_shear_init (Time, G, GV, US, param_file, diag, CS)
 This subroutine initializes the parameters that regulate shear-driven mixing. More...
 
logical function, public kappa_shear_is_used (param_file)
 This function indicates to other modules whether the Jackson et al shear mixing parameterization will be used without needing to duplicate the log entry. More...
 
logical function, public kappa_shear_at_vertex (param_file)
 This function indicates to other modules whether the Jackson et al shear mixing parameterization will be used at the vertices without needing to duplicate the log entry. It returns false if the Jackson et al scheme is not used or if it is used via calculations at the tracer points. More...
 

Function/Subroutine Documentation

◆ calc_kappa_shear_vertex()

subroutine, public mom_kappa_shear::calc_kappa_shear_vertex ( real, dimension(szib_(g),szj_(g),szk_(gv)), intent(in)  u_in,
real, dimension(szi_(g),szjb_(g),szk_(gv)), intent(in)  v_in,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  T_in,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  S_in,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension(:,:), pointer  p_surf,
real, dimension(szi_(g),szj_(g),szk_(gv)+1), intent(out)  kappa_io,
real, dimension(szib_(g),szjb_(g),szk_(gv)+1), intent(out)  tke_io,
real, dimension(szib_(g),szjb_(g),szk_(gv)+1), intent(inout)  kv_io,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(kappa_shear_cs), pointer  CS,
logical, intent(in), optional  initialize_all 
)

Subroutine for calculating shear-driven diffusivity and TKE in corner columns.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]u_inInitial zonal velocity [L T-1 ~> m s-1].
[in]v_inInitial meridional velocity [L T-1 ~> m s-1].
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]t_inLayer potential temperatures [degC]
[in]s_inLayer salinities in ppt.
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
p_surfThe pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
[out]kappa_ioThe diapycnal diffusivity at each interface
[out]tke_ioThe turbulent kinetic energy per unit mass at
[in,out]kv_ioThe vertical viscosity at each interface [Z2 T-1 ~> m2 s-1].
[in]dtTime increment [T ~> s].
csThe control structure returned by a previous call to kappa_shear_init.
[in]initialize_allIf present and false, the previous value of kappa is used to start the iterations

Definition at line 342 of file MOM_kappa_shear.F90.

344  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
345  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
346  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
347  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
348  intent(in) :: u_in !< Initial zonal velocity [L T-1 ~> m s-1].
349  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
350  intent(in) :: v_in !< Initial meridional velocity [L T-1 ~> m s-1].
351  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
352  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
353  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
354  intent(in) :: T_in !< Layer potential temperatures [degC]
355  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
356  intent(in) :: S_in !< Layer salinities in ppt.
357  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
358  !! available thermodynamic fields. Absent fields
359  !! have NULL ptrs.
360  real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
361  !! (or NULL).
362  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
363  intent(out) :: kappa_io !< The diapycnal diffusivity at each interface
364  !! (not layer!) [Z2 T-1 ~> m2 s-1].
365  real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
366  intent(out) :: tke_io !< The turbulent kinetic energy per unit mass at
367  !! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
368  real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
369  intent(inout) :: kv_io !< The vertical viscosity at each interface [Z2 T-1 ~> m2 s-1].
370  !! The previous value is used to initialize kappa
371  !! in the vertex columes as Kappa = Kv/Prandtl
372  !! to accelerate the iteration toward covergence.
373  real, intent(in) :: dt !< Time increment [T ~> s].
374  type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
375  !! call to kappa_shear_init.
376  logical, optional, intent(in) :: initialize_all !< If present and false, the previous
377  !! value of kappa is used to start the iterations
378 
379  ! Local variables
380  real, dimension(SZIB_(G),SZK_(GV)) :: &
381  h_2d, & ! A 2-D version of h, but converted to [Z ~> m].
382  u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
383  T_2d, S_2d, rho_2d ! 2-D versions of T [degC], S [ppt], and rho [R ~> kg m-3].
384  real, dimension(SZIB_(G),SZK_(GV)+1,2) :: &
385  kappa_2d ! Quasi 2-D versions of kappa_io [Z2 T-1 ~> m2 s-1].
386  real, dimension(SZIB_(G),SZK_(GV)+1) :: &
387  tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
388  real, dimension(SZK_(GV)) :: &
389  Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
390  dz, & ! The layer thickness [Z ~> m].
391  u0xdz, & ! The initial zonal velocity times dz [L Z T-1 ~> m2 s-1].
392  v0xdz, & ! The initial meridional velocity times dz [L Z T-1 ~> m2 s-1].
393  T0xdz, & ! The initial temperature times dz [degC Z ~> degC m].
394  S0xdz ! The initial salinity times dz [ppt Z ~> ppt m].
395  real, dimension(SZK_(GV)+1) :: &
396  kappa, & ! The shear-driven diapycnal diffusivity at an interface [Z2 T-1 ~> m2 s-1].
397  tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
398  kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
399  tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
400  real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2].
401  real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa].
402 
403  real :: dz_in_lay ! The running sum of the thickness in a layer [Z ~> m].
404  real :: k0dt ! The background diffusivity times the timestep [Z2 ~> m2].
405  real :: dz_massless ! A layer thickness that is considered massless [Z ~> m].
406  real :: I_hwt ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1].
407  real :: I_Prandtl ! The inverse of the turbulent Prandtl number [nondim].
408  logical :: use_temperature ! If true, temperature and salinity have been
409  ! allocated and are being used as state variables.
410  logical :: new_kappa = .true. ! If true, ignore the value of kappa from the
411  ! last call to this subroutine.
412  logical :: do_i ! If true, work on this column.
413 
414  integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
415  ! interfaces and the interfaces with massless layers
416  ! merged into nearby massive layers.
417  real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
418  ! interpolating back to the original index space [nondim].
419  integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc, J2, J2m1
420 
421  ! Diagnostics that should be deleted?
422  isb = g%isc-1 ; ieb = g%iecB ; jsb = g%jsc-1 ; jeb = g%jecB ; nz = gv%ke
423 
424  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
425  new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all
426 
427  k0dt = dt*cs%kappa_0
428  dz_massless = 0.1*sqrt(k0dt)
429  i_prandtl = 0.0 ; if (cs%Prandtl_turb > 0.0) i_prandtl = 1.0 / cs%Prandtl_turb
430 
431  !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,new_kappa, &
432  !$OMP tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl)
433  do j=jsb,jeb
434  j2 = mod(j,2)+1 ; j2m1 = 3-j2 ! = mod(J-1,2)+1
435 
436  ! Interpolate the various quantities to the corners, using masks.
437  do k=1,nz ; do i=isb,ieb
438  u_2d(i,k) = (u_in(i,j,k) * (g%mask2dCu(i,j) * (h(i,j,k) + h(i+1,j,k))) + &
439  u_in(i,j+1,k) * (g%mask2dCu(i,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) ) / &
440  ((g%mask2dCu(i,j) * (h(i,j,k) + h(i+1,j,k)) + &
441  g%mask2dCu(i,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) + gv%H_subroundoff)
442  v_2d(i,k) = (v_in(i,j,k) * (g%mask2dCv(i,j) * (h(i,j,k) + h(i,j+1,k))) + &
443  v_in(i+1,j,k) * (g%mask2dCv(i+1,j) * (h(i+1,j,k) + h(i+1,j+1,k))) ) / &
444  ((g%mask2dCv(i,j) * (h(i,j,k) + h(i,j+1,k)) + &
445  g%mask2dCv(i+1,j) * (h(i+1,j,k) + h(i+1,j+1,k))) + gv%H_subroundoff)
446  i_hwt = 1.0 / (((g%mask2dT(i,j) * h(i,j,k) + g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
447  (g%mask2dT(i+1,j) * h(i+1,j,k) + g%mask2dT(i,j+1) * h(i,j+1,k))) + &
448  gv%H_subroundoff)
449  if (use_temperature) then
450  t_2d(i,k) = ( ((g%mask2dT(i,j) * h(i,j,k)) * t_in(i,j,k) + &
451  (g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * t_in(i+1,j+1,k)) + &
452  ((g%mask2dT(i+1,j) * h(i+1,j,k)) * t_in(i+1,j,k) + &
453  (g%mask2dT(i,j+1) * h(i,j+1,k)) * t_in(i,j+1,k)) ) * i_hwt
454  s_2d(i,k) = ( ((g%mask2dT(i,j) * h(i,j,k)) * s_in(i,j,k) + &
455  (g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * s_in(i+1,j+1,k)) + &
456  ((g%mask2dT(i+1,j) * h(i+1,j,k)) * s_in(i+1,j,k) + &
457  (g%mask2dT(i,j+1) * h(i,j+1,k)) * s_in(i,j+1,k)) ) * i_hwt
458  endif
459  h_2d(i,k) = gv%H_to_Z * ((g%mask2dT(i,j) * h(i,j,k) + g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
460  (g%mask2dT(i+1,j) * h(i+1,j,k) + g%mask2dT(i,j+1) * h(i,j+1,k)) ) / &
461  ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + &
462  (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)) + 1.0e-36 )
463 ! h_2d(I,k) = 0.25*((h(i,j,k) + h(i+1,j+1,k)) + (h(i+1,j,k) + h(i,j+1,k)))*GV%H_to_Z
464 ! h_2d(I,k) = ((h(i,j,k)**2 + h(i+1,j+1,k)**2) + &
465 ! (h(i+1,j,k)**2 + h(i,j+1,k)**2))*GV%H_to_Z * I_hwt
466  enddo ; enddo
467  if (.not.use_temperature) then ; do k=1,nz ; do i=isb,ieb
468  rho_2d(i,k) = gv%Rlay(k)
469  enddo ; enddo ; endif
470  if (.not.new_kappa) then ; do k=1,nz+1 ; do i=isb,ieb
471  kappa_2d(i,k,j2) = kv_io(i,j,k) * i_prandtl
472  enddo ; enddo ; endif
473 
474 !---------------------------------------
475 ! Work on each column.
476 !---------------------------------------
477  do i=isb,ieb ; if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
478  (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0) then
479  ! call cpu_clock_begin(Id_clock_setup)
480  ! Store a transposed version of the initial arrays.
481  ! Any elimination of massless layers would occur here.
482  if (cs%eliminate_massless) then
483  nzc = 1
484  do k=1,nz
485  ! Zero out the thicknesses of all layers, even if they are unused.
486  dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
487  t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
488 
489  ! Add a new layer if this one has mass.
490 ! if ((dz(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless)) nzc = nzc+1
491  if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
492  (h_2d(i,k) > dz_massless)) nzc = nzc+1
493 
494  ! Only merge clusters of massless layers.
495 ! if ((dz(nzc) > dz_massless) .or. &
496 ! ((dz(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless))) nzc = nzc+1
497 
498  kc(k) = nzc
499  dz(nzc) = dz(nzc) + h_2d(i,k)
500  u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
501  v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
502  if (use_temperature) then
503  t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
504  s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
505  else
506  t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
507  s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
508  endif
509  enddo
510  kc(nz+1) = nzc+1
511 
512  ! Set up Idz as the inverse of layer thicknesses.
513  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
514 
515  ! Now determine kf, the fractional weight of interface kc when
516  ! interpolating between interfaces kc and kc+1.
517  kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
518  do k=2,nz
519  if (kc(k) > kc(k-1)) then
520  kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
521  else
522  kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
523  endif
524  enddo
525  kf(nz+1) = 0.0
526  else
527  do k=1,nz
528  dz(k) = h_2d(i,k)
529  u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
530  enddo
531  if (use_temperature) then
532  do k=1,nz
533  t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
534  enddo
535  else
536  do k=1,nz
537  t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
538  enddo
539  endif
540  nzc = nz
541  do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
542  endif
543  f2 = g%CoriolisBu(i,j)**2
544  surface_pres = 0.0
545  if (associated(p_surf)) then
546  if (cs%psurf_bug) then
547  ! This is wrong because it is averaging values from land in some places.
548  surface_pres = 0.25 * ((p_surf(i,j) + p_surf(i+1,j+1)) + &
549  (p_surf(i+1,j) + p_surf(i,j+1)))
550  else
551  surface_pres = ((g%mask2dT(i,j) * p_surf(i,j) + g%mask2dT(i+1,j+1) * p_surf(i+1,j+1)) + &
552  (g%mask2dT(i+1,j) * p_surf(i+1,j) + g%mask2dT(i,j+1) * p_surf(i,j+1)) ) / &
553  ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + &
554  (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)) + 1.0e-36 )
555  endif
556  endif
557 
558  ! ----------------------------------------------------
559  ! Set the initial guess for kappa, here defined at interfaces.
560  ! ----------------------------------------------------
561  if (new_kappa) then
562  do k=1,nzc+1 ; kappa(k) = us%m2_s_to_Z2_T*1.0 ; enddo
563  else
564  do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k,j2) ; enddo
565  endif
566 
567  call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
568  dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
569  tke_avg, tv, cs, gv, us)
570  ! call cpu_clock_begin(Id_clock_setup)
571  ! Extrapolate from the vertically reduced grid back to the original layers.
572  if (nz == nzc) then
573  do k=1,nz+1
574  kappa_2d(i,k,j2) = kappa_avg(k)
575  if (cs%all_layer_TKE_bug) then
576  tke_2d(i,k) = tke(k)
577  else
578  tke_2d(i,k) = tke_avg(k)
579  endif
580  enddo
581  else
582  do k=1,nz+1
583  if (kf(k) == 0.0) then
584  kappa_2d(i,k,j2) = kappa_avg(kc(k))
585  tke_2d(i,k) = tke_avg(kc(k))
586  else
587  kappa_2d(i,k,j2) = (1.0-kf(k)) * kappa_avg(kc(k)) + kf(k) * kappa_avg(kc(k)+1)
588  tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + kf(k) * tke_avg(kc(k)+1)
589  endif
590  enddo
591  endif
592  ! call cpu_clock_end(Id_clock_setup)
593  else ! Land points, still inside the i-loop.
594  do k=1,nz+1
595  kappa_2d(i,k,j2) = 0.0 ; tke_2d(i,k) = 0.0
596  enddo
597  endif ; enddo ! i-loop
598 
599  do k=1,nz+1 ; do i=isb,ieb
600  tke_io(i,j,k) = g%mask2dBu(i,j) * tke_2d(i,k)
601  kv_io(i,j,k) = ( g%mask2dBu(i,j) * kappa_2d(i,k,j2) ) * cs%Prandtl_turb
602  enddo ; enddo
603  if (j>=g%jsc) then ; do k=1,nz+1 ; do i=g%isc,g%iec
604  ! Set the diffusivities in tracer columns from the values at vertices.
605  kappa_io(i,j,k) = g%mask2dT(i,j) * 0.25 * &
606  ((kappa_2d(i-1,k,j2m1) + kappa_2d(i,k,j2)) + &
607  (kappa_2d(i-1,k,j2) + kappa_2d(i,k,j2m1)))
608  enddo ; enddo ; endif
609 
610  enddo ! end of J-loop
611 
612  if (cs%debug) then
613  call hchksum(kappa_io, "kappa", g%HI, scale=us%Z2_T_to_m2_s)
614  call bchksum(tke_io, "tke", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
615  endif
616 
617  if (cs%id_Kd_shear > 0) call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
618  if (cs%id_TKE > 0) call post_data(cs%id_TKE, tke_io, cs%diag)
619 

◆ calculate_kappa_shear()

subroutine, public mom_kappa_shear::calculate_kappa_shear ( real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  u_in,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  v_in,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension(:,:), pointer  p_surf,
real, dimension(szi_(g),szj_(g),szk_(gv)+1), intent(inout)  kappa_io,
real, dimension(szi_(g),szj_(g),szk_(gv)+1), intent(out)  tke_io,
real, dimension(szi_(g),szj_(g),szk_(gv)+1), intent(inout)  kv_io,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(kappa_shear_cs), pointer  CS,
logical, intent(in), optional  initialize_all 
)

Subroutine for calculating shear-driven diffusivity and TKE in tracer columns.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]u_inInitial zonal velocity [L T-1 ~> m s-1].
[in]v_inInitial meridional velocity [L T-1 ~> m s-1].
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
p_surfThe pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
[in,out]kappa_ioThe diapycnal diffusivity at each interface
[out]tke_ioThe turbulent kinetic energy per unit mass at
[in,out]kv_ioThe vertical viscosity at each interface
[in]dtTime increment [T ~> s].
csThe control structure returned by a previous call to kappa_shear_init.
[in]initialize_allIf present and false, the previous value of kappa is used to start the iterations

Definition at line 109 of file MOM_kappa_shear.F90.

111  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
112  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
113  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
114  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
115  intent(in) :: u_in !< Initial zonal velocity [L T-1 ~> m s-1].
116  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
117  intent(in) :: v_in !< Initial meridional velocity [L T-1 ~> m s-1].
118  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
119  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
120  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
121  !! available thermodynamic fields. Absent fields
122  !! have NULL ptrs.
123  real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
124  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
125  intent(inout) :: kappa_io !< The diapycnal diffusivity at each interface
126  !! (not layer!) [Z2 T-1 ~> m2 s-1]. Initially this is the
127  !! value from the previous timestep, which may
128  !! accelerate the iteration toward convergence.
129  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
130  intent(out) :: tke_io !< The turbulent kinetic energy per unit mass at
131  !! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
132  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
133  intent(inout) :: kv_io !< The vertical viscosity at each interface
134  !! (not layer!) [Z2 T-1 ~> m2 s-1]. This discards any
135  !! previous value (i.e. it is intent out) and
136  !! simply sets Kv = Prandtl * Kd_shear
137  real, intent(in) :: dt !< Time increment [T ~> s].
138  type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
139  !! call to kappa_shear_init.
140  logical, optional, intent(in) :: initialize_all !< If present and false, the previous
141  !! value of kappa is used to start the iterations
142 
143  ! Local variables
144  real, dimension(SZI_(G),SZK_(GV)) :: &
145  h_2d, & ! A 2-D version of h, but converted to [Z ~> m].
146  u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
147  T_2d, S_2d, rho_2d ! 2-D versions of T [degC], S [ppt], and rho [R ~> kg m-3].
148  real, dimension(SZI_(G),SZK_(GV)+1) :: &
149  kappa_2d, & ! 2-D version of kappa_io [Z2 T-1 ~> m2 s-1].
150  tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
151  real, dimension(SZK_(GV)) :: &
152  Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
153  dz, & ! The layer thickness [Z ~> m].
154  u0xdz, & ! The initial zonal velocity times dz [Z L T-1 ~> m2 s-1].
155  v0xdz, & ! The initial meridional velocity times dz [Z L T-1 ~> m2 s-1].
156  T0xdz, & ! The initial temperature times dz [degC Z ~> degC m].
157  S0xdz ! The initial salinity times dz [ppt Z ~> ppt m].
158  real, dimension(SZK_(GV)+1) :: &
159  kappa, & ! The shear-driven diapycnal diffusivity at an interface [Z2 T-1 ~> m2 s-1].
160  tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
161  kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
162  tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
163  real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2].
164  real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa].
165 
166  real :: dz_in_lay ! The running sum of the thickness in a layer [Z ~> m].
167  real :: k0dt ! The background diffusivity times the timestep [Z2 ~> m2].
168  real :: dz_massless ! A layer thickness that is considered massless [Z ~> m].
169  logical :: use_temperature ! If true, temperature and salinity have been
170  ! allocated and are being used as state variables.
171  logical :: new_kappa = .true. ! If true, ignore the value of kappa from the
172  ! last call to this subroutine.
173 
174  integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
175  ! interfaces and the interfaces with massless layers
176  ! merged into nearby massive layers.
177  real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
178  ! interpolating back to the original index space [nondim].
179  integer :: is, ie, js, je, i, j, k, nz, nzc
180 
181  is = g%isc ; ie = g%iec; js = g%jsc ; je = g%jec ; nz = gv%ke
182 
183  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
184  new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all
185 
186  k0dt = dt*cs%kappa_0
187  dz_massless = 0.1*sqrt(k0dt)
188 
189  !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,new_kappa, &
190  !$OMP tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io)
191  do j=js,je
192  do k=1,nz ; do i=is,ie
193  h_2d(i,k) = h(i,j,k)*gv%H_to_Z
194  u_2d(i,k) = u_in(i,j,k) ; v_2d(i,k) = v_in(i,j,k)
195  enddo ; enddo
196  if (use_temperature) then ; do k=1,nz ; do i=is,ie
197  t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
198  enddo ; enddo ; else ; do k=1,nz ; do i=is,ie
199  rho_2d(i,k) = gv%Rlay(k) ! Could be tv%Rho(i,j,k) ?
200  enddo ; enddo ; endif
201  if (.not.new_kappa) then ; do k=1,nz+1 ; do i=is,ie
202  kappa_2d(i,k) = kappa_io(i,j,k)
203  enddo ; enddo ; endif
204 
205 !---------------------------------------
206 ! Work on each column.
207 !---------------------------------------
208  do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
209  ! call cpu_clock_begin(id_clock_setup)
210  ! Store a transposed version of the initial arrays.
211  ! Any elimination of massless layers would occur here.
212  if (cs%eliminate_massless) then
213  nzc = 1
214  do k=1,nz
215  ! Zero out the thicknesses of all layers, even if they are unused.
216  dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
217  t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
218 
219  ! Add a new layer if this one has mass.
220 ! if ((dz(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless)) nzc = nzc+1
221  if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
222  (h_2d(i,k) > dz_massless)) nzc = nzc+1
223 
224  ! Only merge clusters of massless layers.
225 ! if ((dz(nzc) > dz_massless) .or. &
226 ! ((dz(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless))) nzc = nzc+1
227 
228  kc(k) = nzc
229  dz(nzc) = dz(nzc) + h_2d(i,k)
230  u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
231  v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
232  if (use_temperature) then
233  t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
234  s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
235  else
236  t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
237  s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
238  endif
239  enddo
240  kc(nz+1) = nzc+1
241 
242  ! Set up Idz as the inverse of layer thicknesses.
243  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
244 
245  ! Now determine kf, the fractional weight of interface kc when
246  ! interpolating between interfaces kc and kc+1.
247  kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
248  do k=2,nz
249  if (kc(k) > kc(k-1)) then
250  kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
251  else
252  kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
253  endif
254  enddo
255  kf(nz+1) = 0.0
256  else
257  do k=1,nz
258  dz(k) = h_2d(i,k)
259  u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
260  enddo
261  if (use_temperature) then
262  do k=1,nz
263  t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
264  enddo
265  else
266  do k=1,nz
267  t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
268  enddo
269  endif
270  nzc = nz
271  do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
272  endif
273  f2 = 0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
274  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
275  surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j)
276 
277  ! ---------------------------------------------------- I_Ld2_1d, dz_Int_1d
278 
279  ! Set the initial guess for kappa, here defined at interfaces.
280  ! ----------------------------------------------------
281  if (new_kappa) then
282  do k=1,nzc+1 ; kappa(k) = us%m2_s_to_Z2_T*1.0 ; enddo
283  else
284  do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k) ; enddo
285  endif
286 
287  call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
288  dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
289  tke_avg, tv, cs, gv, us)
290 
291  ! call cpu_clock_begin(id_clock_setup)
292  ! Extrapolate from the vertically reduced grid back to the original layers.
293  if (nz == nzc) then
294  do k=1,nz+1
295  kappa_2d(i,k) = kappa_avg(k)
296  if (cs%all_layer_TKE_bug) then
297  tke_2d(i,k) = tke(k)
298  else
299  tke_2d(i,k) = tke_avg(k)
300  endif
301  enddo
302  else
303  do k=1,nz+1
304  if (kf(k) == 0.0) then
305  kappa_2d(i,k) = kappa_avg(kc(k))
306  tke_2d(i,k) = tke_avg(kc(k))
307  else
308  kappa_2d(i,k) = (1.0-kf(k)) * kappa_avg(kc(k)) + &
309  kf(k) * kappa_avg(kc(k)+1)
310  tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + &
311  kf(k) * tke_avg(kc(k)+1)
312  endif
313  enddo
314  endif
315  ! call cpu_clock_end(id_clock_setup)
316  else ! Land points, still inside the i-loop.
317  do k=1,nz+1
318  kappa_2d(i,k) = 0.0 ; tke_2d(i,k) = 0.0
319  enddo
320  endif ; enddo ! i-loop
321 
322  do k=1,nz+1 ; do i=is,ie
323  kappa_io(i,j,k) = g%mask2dT(i,j) * kappa_2d(i,k)
324  tke_io(i,j,k) = g%mask2dT(i,j) * tke_2d(i,k)
325  kv_io(i,j,k) = ( g%mask2dT(i,j) * kappa_2d(i,k) ) * cs%Prandtl_turb
326  enddo ; enddo
327 
328  enddo ! end of j-loop
329 
330  if (cs%debug) then
331  call hchksum(kappa_io, "kappa", g%HI, scale=us%Z2_T_to_m2_s)
332  call hchksum(tke_io, "tke", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
333  endif
334 
335  if (cs%id_Kd_shear > 0) call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
336  if (cs%id_TKE > 0) call post_data(cs%id_TKE, tke_io, cs%diag)
337 

◆ calculate_projected_state()

subroutine mom_kappa_shear::calculate_projected_state ( real, dimension(nz+1), intent(in)  kappa,
real, dimension(nz), intent(in)  u0,
real, dimension(nz), intent(in)  v0,
real, dimension(nz), intent(in)  T0,
real, dimension(nz), intent(in)  S0,
real, intent(in)  dt,
integer, intent(in)  nz,
real, dimension(nz), intent(in)  dz,
real, dimension(nz+1), intent(in)  I_dz_int,
real, dimension(nz+1), intent(in)  dbuoy_dT,
real, dimension(nz+1), intent(in)  dbuoy_dS,
real, dimension(nz), intent(inout)  u,
real, dimension(nz), intent(inout)  v,
real, dimension(nz), intent(inout)  T,
real, dimension(nz), intent(inout)  Sal,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension(nz+1), intent(inout), optional  N2,
real, dimension(nz+1), intent(inout), optional  S2,
integer, intent(in), optional  ks_int,
integer, intent(in), optional  ke_int,
real, intent(in), optional  vel_underflow 
)
private

This subroutine calculates the velocities, temperature and salinity that the water column will have after mixing for dt with diffusivities kappa. It may also calculate the projected buoyancy frequency and shear.

Parameters
[in]nzThe number of layers (after eliminating massless layers?).
[in]kappaThe diapycnal diffusivity at interfaces, [Z2 T-1 ~> m2 s-1].
[in]u0The initial zonal velocity [L T-1 ~> m s-1].
[in]v0The initial meridional velocity [L T-1 ~> m s-1].
[in]t0The initial temperature [degC].
[in]s0The initial salinity [ppt].
[in]dzThe grid spacing of layers [Z ~> m].
[in]i_dz_intThe inverse of the layer's thicknesses [Z-1 ~> m-1].
[in]dbuoy_dtThe partial derivative of buoyancy with temperature [Z T-2 degC-1 ~> m s-2 degC-1].
[in]dbuoy_dsThe partial derivative of buoyancy with salinity [Z T-2 ppt-1 ~> m s-2 ppt-1].
[in]dtThe time step [T ~> s].
[in,out]uThe zonal velocity after dt [L T-1 ~> m s-1].
[in,out]vThe meridional velocity after dt [L T-1 ~> m s-1].
[in,out]tThe temperature after dt [degC].
[in,out]salThe salinity after dt [ppt].
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in,out]n2The buoyancy frequency squared at interfaces [T-2 ~> s-2].
[in,out]s2The squared shear at interfaces [T-2 ~> s-2].
[in]ks_intThe topmost k-index with a non-zero diffusivity.
[in]ke_intThe bottommost k-index with a non-zero diffusivity.
[in]vel_underflowIf present and true, any velocities that are smaller in magnitude than this value are set to 0 [L T-1 ~> m s-1].

Definition at line 1079 of file MOM_kappa_shear.F90.

1082  integer, intent(in) :: nz !< The number of layers (after eliminating massless
1083  !! layers?).
1084  real, dimension(nz+1), intent(in) :: kappa !< The diapycnal diffusivity at interfaces,
1085  !! [Z2 T-1 ~> m2 s-1].
1086  real, dimension(nz), intent(in) :: u0 !< The initial zonal velocity [L T-1 ~> m s-1].
1087  real, dimension(nz), intent(in) :: v0 !< The initial meridional velocity [L T-1 ~> m s-1].
1088  real, dimension(nz), intent(in) :: T0 !< The initial temperature [degC].
1089  real, dimension(nz), intent(in) :: S0 !< The initial salinity [ppt].
1090  real, dimension(nz), intent(in) :: dz !< The grid spacing of layers [Z ~> m].
1091  real, dimension(nz+1), intent(in) :: I_dz_int !< The inverse of the layer's thicknesses
1092  !! [Z-1 ~> m-1].
1093  real, dimension(nz+1), intent(in) :: dbuoy_dT !< The partial derivative of buoyancy with
1094  !! temperature [Z T-2 degC-1 ~> m s-2 degC-1].
1095  real, dimension(nz+1), intent(in) :: dbuoy_dS !< The partial derivative of buoyancy with
1096  !! salinity [Z T-2 ppt-1 ~> m s-2 ppt-1].
1097  real, intent(in) :: dt !< The time step [T ~> s].
1098  real, dimension(nz), intent(inout) :: u !< The zonal velocity after dt [L T-1 ~> m s-1].
1099  real, dimension(nz), intent(inout) :: v !< The meridional velocity after dt [L T-1 ~> m s-1].
1100  real, dimension(nz), intent(inout) :: T !< The temperature after dt [degC].
1101  real, dimension(nz), intent(inout) :: Sal !< The salinity after dt [ppt].
1102  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1103  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1104  real, dimension(nz+1), optional, &
1105  intent(inout) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
1106  real, dimension(nz+1), optional, &
1107  intent(inout) :: S2 !< The squared shear at interfaces [T-2 ~> s-2].
1108  integer, optional, intent(in) :: ks_int !< The topmost k-index with a non-zero diffusivity.
1109  integer, optional, intent(in) :: ke_int !< The bottommost k-index with a non-zero
1110  !! diffusivity.
1111  real, optional, intent(in) :: vel_underflow !< If present and true, any velocities that
1112  !! are smaller in magnitude than this value are
1113  !! set to 0 [L T-1 ~> m s-1].
1114 
1115  ! Local variables
1116  real, dimension(nz+1) :: c1
1117  real :: L2_to_Z2 ! A conversion factor from horizontal length units to vertical depth
1118  ! units squared [Z2 s2 T-2 m-2 ~> 1].
1119  real :: underflow_vel ! Velocities smaller in magnitude than underflow_vel are set to 0 [L T-1 ~> m s-1].
1120  real :: a_a, a_b, b1, d1, bd1, b1nz_0
1121  integer :: k, ks, ke
1122 
1123  ks = 1 ; ke = nz
1124  if (present(ks_int)) ks = max(ks_int-1,1)
1125  if (present(ke_int)) ke = min(ke_int,nz)
1126  underflow_vel = 0.0 ; if (present(vel_underflow)) underflow_vel = vel_underflow
1127 
1128  if (ks > ke) return
1129 
1130  if (dt > 0.0) then
1131  a_b = dt*(kappa(ks+1)*i_dz_int(ks+1))
1132  b1 = 1.0 / (dz(ks) + a_b)
1133  c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1 ! = 1 - c1
1134 
1135  u(ks) = (b1 * dz(ks))*u0(ks) ; v(ks) = (b1 * dz(ks))*v0(ks)
1136  t(ks) = (b1 * dz(ks))*t0(ks) ; sal(ks) = (b1 * dz(ks))*s0(ks)
1137  do k=ks+1,ke-1
1138  a_a = a_b
1139  a_b = dt*(kappa(k+1)*i_dz_int(k+1))
1140  bd1 = dz(k) + d1*a_a
1141  b1 = 1.0 / (bd1 + a_b)
1142  c1(k+1) = a_b * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
1143 
1144  u(k) = b1 * (dz(k)*u0(k) + a_a*u(k-1))
1145  v(k) = b1 * (dz(k)*v0(k) + a_a*v(k-1))
1146  t(k) = b1 * (dz(k)*t0(k) + a_a*t(k-1))
1147  sal(k) = b1 * (dz(k)*s0(k) + a_a*sal(k-1))
1148  enddo
1149  ! T and S have insulating boundary conditions, u & v use no-slip
1150  ! bottom boundary conditions at the solid bottom.
1151 
1152  ! For insulating boundary conditions or mixing simply stopping, use...
1153  a_a = a_b
1154  b1 = 1.0 / (dz(ke) + d1*a_a)
1155  t(ke) = b1 * (dz(ke)*t0(ke) + a_a*t(ke-1))
1156  sal(ke) = b1 * (dz(ke)*s0(ke) + a_a*sal(ke-1))
1157 
1158  ! There is no distinction between the effective boundary conditions for
1159  ! tracers and velocities if the mixing is separated from the bottom, but if
1160  ! the mixing goes all the way to the bottom, use no-slip BCs for velocities.
1161  if (ke == nz) then
1162  a_b = dt*(kappa(nz+1)*i_dz_int(nz+1))
1163  b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
1164  else
1165  b1nz_0 = b1
1166  endif
1167  u(ke) = b1nz_0 * (dz(ke)*u0(ke) + a_a*u(ke-1))
1168  v(ke) = b1nz_0 * (dz(ke)*v0(ke) + a_a*v(ke-1))
1169  if (abs(u(ke)) < underflow_vel) u(ke) = 0.0
1170  if (abs(v(ke)) < underflow_vel) v(ke) = 0.0
1171 
1172  do k=ke-1,ks,-1
1173  u(k) = u(k) + c1(k+1)*u(k+1)
1174  v(k) = v(k) + c1(k+1)*v(k+1)
1175  if (abs(u(k)) < underflow_vel) u(k) = 0.0
1176  if (abs(v(k)) < underflow_vel) v(k) = 0.0
1177  t(k) = t(k) + c1(k+1)*t(k+1)
1178  sal(k) = sal(k) + c1(k+1)*sal(k+1)
1179  enddo
1180  else ! dt <= 0.0
1181  do k=1,nz
1182  u(k) = u0(k) ; v(k) = v0(k) ; t(k) = t0(k) ; sal(k) = s0(k)
1183  if (abs(u(k)) < underflow_vel) u(k) = 0.0
1184  if (abs(v(k)) < underflow_vel) v(k) = 0.0
1185  enddo
1186  endif
1187 
1188  if (present(s2)) then
1189  ! L2_to_Z2 = US%m_to_Z**2 * US%T_to_s**2
1190  l2_to_z2 = us%L_to_Z**2
1191  s2(1) = 0.0 ; s2(nz+1) = 0.0
1192  if (ks > 1) &
1193  s2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (l2_to_z2*i_dz_int(ks)**2)
1194  do k=ks+1,ke
1195  s2(k) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (l2_to_z2*i_dz_int(k)**2)
1196  enddo
1197  if (ke<nz) &
1198  s2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * (l2_to_z2*i_dz_int(ke+1)**2)
1199  endif
1200 
1201  if (present(n2)) then
1202  n2(1) = 0.0 ; n2(nz+1) = 0.0
1203  if (ks > 1) &
1204  n2(ks) = max(0.0, i_dz_int(ks) * &
1205  (dbuoy_dt(ks) * (t0(ks-1)-t(ks)) + dbuoy_ds(ks) * (s0(ks-1)-sal(ks))))
1206  do k=ks+1,ke
1207  n2(k) = max(0.0, i_dz_int(k) * &
1208  (dbuoy_dt(k) * (t(k-1)-t(k)) + dbuoy_ds(k) * (sal(k-1)-sal(k))))
1209  enddo
1210  if (ke<nz) &
1211  n2(ke+1) = max(0.0, i_dz_int(ke+1) * &
1212  (dbuoy_dt(ke+1) * (t(ke)-t0(ke+1)) + dbuoy_ds(ke+1) * (sal(ke)-s0(ke+1))))
1213  endif
1214 

◆ find_kappa_tke()

subroutine mom_kappa_shear::find_kappa_tke ( real, dimension(nz+1), intent(in)  N2,
real, dimension(nz+1), intent(in)  S2,
real, dimension(nz+1), intent(in)  kappa_in,
real, dimension(nz), intent(in)  Idz,
real, dimension(nz+1), intent(in)  dz_Int,
real, dimension(nz+1), intent(in)  I_L2_bdry,
real, intent(in)  f2,
integer, intent(in)  nz,
type(kappa_shear_cs), pointer  CS,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension(nz+1), intent(inout)  K_Q,
real, dimension(nz+1), intent(out)  tke,
real, dimension(nz+1), intent(out)  kappa,
real, dimension(nz+1), intent(out), optional  kappa_src,
real, dimension(nz+1), intent(out), optional  local_src 
)
private

This subroutine calculates new, consistent estimates of TKE and kappa.

Parameters
[in]nzThe number of layers to work on.
[in]n2The buoyancy frequency squared at interfaces [T-2 ~> s-2].
[in]s2The squared shear at interfaces [T-2 ~> s-2].
[in]kappa_inThe initial guess at the diffusivity [Z2 T-1 ~> m2 s-1].
[in]dz_intThe thicknesses associated with interfaces [Z-1 ~> m-1].
[in]i_l2_bdryThe inverse of the squared distance to boundaries [Z-2 ~> m-2].
[in]idzThe inverse grid spacing of layers [Z-1 ~> m-1].
[in]f2The squared Coriolis parameter [T-2 ~> s-2].
csA pointer to this module's control structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in,out]k_qThe shear-driven diapycnal diffusivity divided by the turbulent kinetic energy per unit mass at interfaces [T ~> s].
[out]tkeThe turbulent kinetic energy per unit mass at interfaces [Z2 T-2 ~> m2 s-2].
[out]kappaThe diapycnal diffusivity at interfaces [Z2 T-1 ~> m2 s-1].
[out]kappa_srcThe source term for kappa [T-1 ~> s-1].
[out]local_srcThe sum of all local sources for kappa,

Definition at line 1218 of file MOM_kappa_shear.F90.

1220  integer, intent(in) :: nz !< The number of layers to work on.
1221  real, dimension(nz+1), intent(in) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
1222  real, dimension(nz+1), intent(in) :: S2 !< The squared shear at interfaces [T-2 ~> s-2].
1223  real, dimension(nz+1), intent(in) :: kappa_in !< The initial guess at the diffusivity
1224  !! [Z2 T-1 ~> m2 s-1].
1225  real, dimension(nz+1), intent(in) :: dz_Int !< The thicknesses associated with interfaces
1226  !! [Z-1 ~> m-1].
1227  real, dimension(nz+1), intent(in) :: I_L2_bdry !< The inverse of the squared distance to
1228  !! boundaries [Z-2 ~> m-2].
1229  real, dimension(nz), intent(in) :: Idz !< The inverse grid spacing of layers [Z-1 ~> m-1].
1230  real, intent(in) :: f2 !< The squared Coriolis parameter [T-2 ~> s-2].
1231  type(Kappa_shear_CS), pointer :: CS !< A pointer to this module's control structure.
1232  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1233  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1234  real, dimension(nz+1), intent(inout) :: K_Q !< The shear-driven diapycnal diffusivity divided by
1235  !! the turbulent kinetic energy per unit mass at
1236  !! interfaces [T ~> s].
1237  real, dimension(nz+1), intent(out) :: tke !< The turbulent kinetic energy per unit mass at
1238  !! interfaces [Z2 T-2 ~> m2 s-2].
1239  real, dimension(nz+1), intent(out) :: kappa !< The diapycnal diffusivity at interfaces
1240  !! [Z2 T-1 ~> m2 s-1].
1241  real, dimension(nz+1), optional, &
1242  intent(out) :: kappa_src !< The source term for kappa [T-1 ~> s-1].
1243  real, dimension(nz+1), optional, &
1244  intent(out) :: local_src !< The sum of all local sources for kappa,
1245  !! [T-1 ~> s-1].
1246  ! This subroutine calculates new, consistent estimates of TKE and kappa.
1247 
1248  ! Local variables
1249  real, dimension(nz) :: &
1250  aQ, & ! aQ is the coupling between adjacent interfaces in the TKE equations [Z T-1 ~> m s-1].
1251  dQdz ! Half the partial derivative of TKE with depth [Z T-2 ~> m s-2].
1252  real, dimension(nz+1) :: &
1253  dK, & ! The change in kappa [Z2 T-1 ~> m2 s-1].
1254  dQ, & ! The change in TKE [Z2 T-2 ~> m2 s-2].
1255  cQ, cK, & ! cQ and cK are the upward influences in the tridiagonal and
1256  ! hexadiagonal solvers for the TKE and kappa equations [nondim].
1257  i_ld2, & ! 1/Ld^2, where Ld is the effective decay length scale for kappa [Z-2 ~> m-2].
1258  tke_decay, & ! The local TKE decay rate [T-1 ~> s-1].
1259  k_src, & ! The source term in the kappa equation [T-1 ~> s-1].
1260  dqmdk, & ! With Newton's method the change in dQ(k-1) due to dK(k) [T ~> s].
1261  dkdq, & ! With Newton's method the change in dK(k) due to dQ(k) [T-1 ~> s-1].
1262  e1 ! The fractional change in a layer TKE due to a change in the
1263  ! TKE of the layer above when all the kappas below are 0.
1264  ! e1 is nondimensional, and 0 < e1 < 1.
1265  real :: tke_src ! The net source of TKE due to mixing against the shear
1266  ! and stratification [Z2 T-3 ~> m2 s-3]. (For convenience,
1267  ! a term involving the non-dissipation of q0 is also
1268  ! included here.)
1269  real :: bQ ! The inverse of the pivot in the tridiagonal equations [T Z-1 ~> s m-1].
1270  real :: bK ! The inverse of the pivot in the tridiagonal equations [Z-1 ~> m-1].
1271  real :: bQd1 ! A term in the denominator of bQ [Z T-1 ~> m s-1].
1272  real :: bKd1 ! A term in the denominator of bK [Z ~> m].
1273  real :: cQcomp, cKcomp ! 1 - cQ or 1 - cK in the tridiagonal equations.
1274  real :: c_s2 ! The coefficient for the decay of TKE due to
1275  ! shear (i.e. proportional to |S|*tke), nondimensional.
1276  real :: c_n2 ! The coefficient for the decay of TKE due to
1277  ! stratification (i.e. proportional to N*tke) [nondim].
1278  real :: Ri_crit ! The critical shear Richardson number for shear-
1279  ! driven mixing. The theoretical value is 0.25.
1280  real :: q0 ! The background level of TKE [Z2 T-2 ~> m2 s-2].
1281  real :: Ilambda2 ! 1.0 / CS%lambda**2 [nondim]
1282  real :: TKE_min ! The minimum value of shear-driven TKE that can be
1283  ! solved for [Z2 T-2 ~> m2 s-2].
1284  real :: kappa0 ! The background diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
1285  real :: kappa_trunc ! Diffusivities smaller than this are rounded to 0 [Z2 T-1 ~> m2 s-1].
1286 
1287  real :: eden1, eden2, I_eden, ome ! Variables used in calculating e1.
1288  real :: diffusive_src ! The diffusive source in the kappa equation [Z T-1 ~> m s-1].
1289  real :: chg_by_k0 ! The value of k_src that leads to an increase of
1290  ! kappa_0 if only the diffusive term is a sink [T-1 ~> s-1].
1291 
1292  real :: kappa_mean ! A mean value of kappa [Z2 T-1 ~> m2 s-1].
1293  real :: Newton_test ! The value of relative error that will cause the next
1294  ! iteration to use Newton's method.
1295  ! Temporary variables used in the Newton's method iterations.
1296  real :: decay_term_k ! The decay term in the diffusivity equation
1297  real :: decay_term_Q ! The decay term in the TKE equation - proportional to [T-1 ~> s-1]
1298  real :: I_Q ! The inverse of TKE [T2 Z-2 ~> s2 m-2]
1299  real :: kap_src
1300  real :: v1 ! A temporary variable proportional to [T-1 ~> s-1]
1301  real :: v2
1302  real :: tol_err ! The tolerance for max_err that determines when to
1303  ! stop iterating.
1304  real :: Newton_err ! The tolerance for max_err that determines when to
1305  ! start using Newton's method. Empirically, an initial
1306  ! value of about 0.2 seems to be most efficient.
1307  real, parameter :: roundoff = 1.0e-16 ! A negligible fractional change in TKE.
1308  ! This could be larger but performance gains are small.
1309 
1310  logical :: tke_noflux_bottom_BC = .false. ! Specify the boundary conditions
1311  logical :: tke_noflux_top_BC = .false. ! that are applied to the TKE eqns.
1312  logical :: do_Newton ! If .true., use Newton's method for the next iteration.
1313  logical :: abort_Newton ! If .true., an Newton's method has encountered a 0
1314  ! pivot, and should not have been used.
1315  logical :: was_Newton ! The value of do_Newton before checking convergence.
1316  logical :: within_tolerance ! If .true., all points are within tolerance to
1317  ! enable this subroutine to return.
1318  integer :: ks_src, ke_src ! The range indices that have nonzero k_src.
1319  integer :: ks_kappa, ke_kappa, ke_tke ! The ranges of k-indices that are or
1320  integer :: ks_kappa_prev, ke_kappa_prev ! were being worked on.
1321  integer :: itt, k, k2
1322 
1323  ! These variables are used only for debugging.
1324  logical, parameter :: debug_soln = .false.
1325  real :: K_err_lin, Q_err_lin, TKE_src_norm
1326  real, dimension(nz+1) :: &
1327  I_Ld2_debug, & ! A separate version of I_Ld2 for debugging [Z-2 ~> m-2].
1328  kappa_prev, & ! The value of kappa at the start of the current iteration [Z2 T-1 ~> m2 s-1].
1329  TKE_prev ! The value of TKE at the start of the current iteration [Z2 T-2 ~> m2 s-2].
1330 
1331  c_n2 = cs%C_N**2 ; c_s2 = cs%C_S**2
1332  q0 = cs%TKE_bg ; kappa0 = cs%kappa_0
1333  tke_min = max(cs%TKE_bg, 1.0e-20*us%m_to_Z**2*us%T_to_s**2)
1334  ri_crit = cs%Rino_crit
1335  ilambda2 = 1.0 / cs%lambda**2
1336  kappa_trunc = cs%kappa_trunc
1337  do_newton = .false. ; abort_newton = .false.
1338  tol_err = cs%kappa_tol_err
1339  newton_err = 0.2 ! This initial value may be automatically reduced later.
1340 
1341  ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
1342 
1343  ke_src = 0 ; ks_src = nz+1
1344  do k=2,nz
1345  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
1346 ! Ri = N2(K) / S2(K)
1347 ! k_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1348 ! ((Ri_crit - Ri) / (Ri_crit + CS%FRi_curvature*Ri))
1349  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1350  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1351  ke_src = k
1352  if (ks_src > k) ks_src = k
1353  else
1354  k_src(k) = 0.0
1355  endif
1356  enddo
1357 
1358  ! If there is no source anywhere, return kappa(K) = 0.
1359  if (ks_src > ke_src) then
1360  do k=1,nz+1
1361  kappa(k) = 0.0 ; k_q(k) = 0.0 ; tke(k) = tke_min
1362  enddo
1363  if (present(kappa_src)) then ; do k=1,nz+1 ; kappa_src(k) = 0.0 ; enddo ; endif
1364  if (present(local_src)) then ; do k=1,nz+1 ; local_src(k) = 0.0 ; enddo ; endif
1365  return
1366  endif
1367 
1368  do k=1,nz+1
1369  kappa(k) = kappa_in(k)
1370 ! TKE_decay(K) = c_n*sqrt(N2(K)) + c_s*sqrt(S2(K)) ! The expression in JHL.
1371  tke_decay(k) = sqrt(c_n2*n2(k) + c_s2*s2(k))
1372  if ((kappa(k) > 0.0) .and. (k_q(k) > 0.0)) then
1373  tke(k) = kappa(k) / k_q(k) ! Perhaps take the max with TKE_min
1374  else
1375  tke(k) = tke_min
1376  endif
1377  enddo
1378  ! Apply boundary conditions to kappa.
1379  kappa(1) = 0.0 ; kappa(nz+1) = 0.0
1380 
1381  ! Calculate the term (e1) that allows changes in TKE to be calculated quickly
1382  ! below the deepest nonzero value of kappa. If kappa = 0, below interface
1383  ! k-1, the final changes in TKE are related by dQ(K+1) = e1(K+1)*dQ(K).
1384  eden2 = kappa0 * idz(nz)
1385  if (tke_noflux_bottom_bc) then
1386  eden1 = dz_int(nz+1)*tke_decay(nz+1)
1387  i_eden = 1.0 / (eden2 + eden1)
1388  e1(nz+1) = eden2 * i_eden ; ome = eden1 * i_eden
1389  else
1390  e1(nz+1) = 0.0 ; ome = 1.0
1391  endif
1392  do k=nz,2,-1
1393  eden1 = dz_int(k)*tke_decay(k) + ome * eden2
1394  eden2 = kappa0 * idz(k-1)
1395  i_eden = 1.0 / (eden2 + eden1)
1396  e1(k) = eden2 * i_eden ; ome = eden1 * i_eden ! = 1-e1
1397  enddo
1398  e1(1) = 0.0
1399 
1400 
1401  ! Iterate here to convergence to within some tolerance of order tol_err.
1402  do itt=1,cs%max_RiNo_it
1403 
1404  ! ----------------------------------------------------
1405  ! Calculate TKE
1406  ! ----------------------------------------------------
1407 
1408  if (debug_soln) then ; do k=1,nz+1 ; kappa_prev(k) = kappa(k) ; tke_prev(k) = tke(k) ; enddo ; endif
1409 
1410  if (.not.do_newton) then
1411  ! Use separate steps of the TKE and kappa equations, that are
1412  ! explicit in the nonlinear source terms, implicit in a linearized
1413  ! version of the nonlinear sink terms, and implicit in the linear
1414  ! terms.
1415 
1416  ke_tke = max(ke_kappa,ke_kappa_prev)+1
1417  ! aQ is the coupling between adjacent interfaces [Z T-1 ~> m s-1].
1418  do k=1,min(ke_tke,nz)
1419  aq(k) = (0.5*(kappa(k)+kappa(k+1)) + kappa0) * idz(k)
1420  enddo
1421  dq(1) = -tke(1)
1422  if (tke_noflux_top_bc) then
1423  tke_src = kappa0*s2(1) + q0 * tke_decay(1) ! Uses that kappa(1) = 0
1424  bqd1 = dz_int(1) * tke_decay(1)
1425  bq = 1.0 / (bqd1 + aq(1))
1426  tke(1) = bq * (dz_int(1)*tke_src)
1427  cq(2) = aq(1) * bq ; cqcomp = bqd1 * bq ! = 1 - cQ
1428  else
1429  tke(1) = q0 ; cq(2) = 0.0 ; cqcomp = 1.0
1430  endif
1431  do k=2,ke_tke-1
1432  dq(k) = -tke(k)
1433  tke_src = (kappa(k) + kappa0)*s2(k) + q0*tke_decay(k)
1434  bqd1 = dz_int(k)*(tke_decay(k) + n2(k)*k_q(k)) + cqcomp*aq(k-1)
1435  bq = 1.0 / (bqd1 + aq(k))
1436  tke(k) = bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1))
1437  cq(k+1) = aq(k) * bq ; cqcomp = bqd1 * bq ! = 1 - cQ
1438  enddo
1439  if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_bc)) then
1440  tke(nz+1) = tke_min
1441  dq(nz+1) = 0.0
1442  else
1443  k = ke_tke
1444  tke_src = kappa0*s2(k) + q0*tke_decay(k) ! Uses that kappa(ke_tke) = 0
1445  if (k == nz+1) then
1446  dq(k) = -tke(k)
1447  bq = 1.0 / (dz_int(k)*tke_decay(k) + cqcomp*aq(k-1))
1448  tke(k) = max(tke_min, bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)))
1449  dq(k) = tke(k) + dq(k)
1450  else
1451  bq = 1.0 / ((dz_int(k)*tke_decay(k) + cqcomp*aq(k-1)) + aq(k))
1452  cq(k+1) = aq(k) * bq
1453  ! Account for all changes deeper in the water column.
1454  dq(k) = -tke(k)
1455  tke(k) = max((bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)) + &
1456  cq(k+1)*(tke(k+1) - e1(k+1)*tke(k))) / (1.0 - cq(k+1)*e1(k+1)), tke_min)
1457  dq(k) = tke(k) + dq(k)
1458 
1459  ! Adjust TKE deeper in the water column in case ke_tke increases.
1460  ! This might not be strictly necessary?
1461  do k=ke_tke+1,nz+1
1462  dq(k) = e1(k)*dq(k-1)
1463  tke(k) = max(tke(k) + dq(k), tke_min)
1464  if (abs(dq(k)) < roundoff*tke(k)) exit
1465  enddo
1466  do k2=k+1,nz
1467  if (dq(k2) == 0.0) exit
1468  dq(k2) = 0.0
1469  enddo
1470  endif
1471  endif
1472  do k=ke_tke-1,1,-1
1473  tke(k) = max(tke(k) + cq(k+1)*tke(k+1), tke_min)
1474  dq(k) = tke(k) + dq(k)
1475  enddo
1476 
1477  ! ----------------------------------------------------
1478  ! Calculate kappa, here defined at interfaces.
1479  ! ----------------------------------------------------
1480 
1481  ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
1482 
1483  dk(1) = 0.0 ! kappa takes boundary values of 0.
1484  ck(2) = 0.0 ; ckcomp = 1.0
1485  if (itt == 1) then ; dO k=2,nz
1486  i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1487  enddo ; endif
1488  do k=2,nz
1489  dk(k) = -kappa(k)
1490  if (itt>1) &
1491  i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1492  bkd1 = dz_int(k)*i_ld2(k) + ckcomp*idz(k-1)
1493  bk = 1.0 / (bkd1 + idz(k))
1494 
1495  kappa(k) = bk * (idz(k-1)*kappa(k-1) + dz_int(k) * k_src(k))
1496  ck(k+1) = idz(k) * bk ; ckcomp = bkd1 * bk ! = 1 - cK(K+1)
1497 
1498  ! Neglect values that are smaller than kappa_trunc.
1499  if (kappa(k) < ckcomp*kappa_trunc) then
1500  kappa(k) = 0.0
1501  if (k > ke_src) then ; ke_kappa = k-1 ; k_q(k) = 0.0 ; exit ; endif
1502  elseif (kappa(k) < 2.0*ckcomp*kappa_trunc) then
1503  kappa(k) = 2.0 * (kappa(k) - ckcomp*kappa_trunc)
1504  endif
1505  enddo
1506  k_q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
1507  dk(ke_kappa) = dk(ke_kappa) + kappa(ke_kappa)
1508  do k=ke_kappa+2,ke_kappa_prev
1509  dk(k) = -kappa(k) ; kappa(k) = 0.0 ; k_q(k) = 0.0
1510  enddo
1511  do k=ke_kappa-1,2,-1
1512  kappa(k) = kappa(k) + ck(k+1)*kappa(k+1)
1513  ! Neglect values that are smaller than kappa_trunc.
1514  if (kappa(k) <= kappa_trunc) then
1515  kappa(k) = 0.0
1516  if (k < ks_src) then ; ks_kappa = k+1 ; k_q(k) = 0.0 ; exit ; endif
1517  elseif (kappa(k) < 2.0*kappa_trunc) then
1518  kappa(k) = 2.0 * (kappa(k) - kappa_trunc)
1519  endif
1520 
1521  dk(k) = dk(k) + kappa(k)
1522  k_q(k) = kappa(k) / tke(k)
1523  enddo
1524  do k=ks_kappa_prev,ks_kappa-2 ; kappa(k) = 0.0 ; k_q(k) = 0.0 ; enddo
1525 
1526  else ! do_Newton is .true.
1527 ! Once the solutions are close enough, use a Newton's method solver of the
1528 ! whole system to accelerate convergence.
1529  ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
1530  ks_kappa = 2
1531  dk(1) = 0.0 ; ck(2) = 0.0 ; ckcomp = 1.0 ; dkdq(1) = 0.0
1532  aq(1) = (0.5*(kappa(1)+kappa(2))+kappa0) * idz(1)
1533  dqdz(1) = 0.5*(tke(1) - tke(2))*idz(1)
1534  if (tke_noflux_top_bc) then
1535  tke_src = dz_int(1) * (kappa0*s2(1) - (tke(1) - q0)*tke_decay(1)) - &
1536  aq(1) * (tke(1) - tke(2))
1537 
1538  bq = 1.0 / (aq(1) + dz_int(1)*tke_decay(1))
1539  cq(2) = aq(1) * bq
1540  cqcomp = (dz_int(1)*tke_decay(1)) * bq ! = 1 - cQ(2)
1541  dqmdk(2) = -dqdz(1) * bq
1542  dq(1) = bq * tke_src
1543  else
1544  dq(1) = 0.0 ; cq(2) = 0.0 ; cqcomp = 1.0 ; dqmdk(2) = 0.0
1545  endif
1546  do k=2,nz
1547  i_q = 1.0 / tke(k)
1548  i_ld2(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1549 
1550  kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa(k)) + &
1551  idz(k-1)*(kappa(k-1)-kappa(k)) - idz(k)*(kappa(k)-kappa(k+1))
1552 
1553  ! Ensure that the pivot is always positive, and that 0 <= cK <= 1.
1554  ! Otherwise do not use Newton's method.
1555  decay_term_k = -idz(k-1)*dqmdk(k)*dkdq(k-1) + dz_int(k)*i_ld2(k)
1556  if (decay_term_k < 0.0) then ; abort_newton = .true. ; exit ; endif
1557  bk = 1.0 / (idz(k) + idz(k-1)*ckcomp + decay_term_k)
1558 
1559  ck(k+1) = bk * idz(k)
1560  ckcomp = bk * (idz(k-1)*ckcomp + decay_term_k) ! = 1-cK(K+1)
1561  if (cs%dKdQ_iteration_bug) then
1562  dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1563  us%m_to_Z*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
1564  else
1565  dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1566  dz_int(k)*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
1567  endif
1568  dk(k) = bk * (kap_src + idz(k-1)*dk(k-1) + idz(k-1)*dkdq(k-1)*dq(k-1))
1569 
1570  ! Truncate away negligibly small values of kappa.
1571  if (dk(k) <= ckcomp*(kappa_trunc - kappa(k))) then
1572  dk(k) = -ckcomp*kappa(k)
1573 ! if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1574  elseif (dk(k) < ckcomp*(2.0*kappa_trunc - kappa(k))) then
1575  dk(k) = 2.0 * dk(k) - ckcomp*(2.0*kappa_trunc - kappa(k))
1576  endif
1577 
1578  ! Solve for dQ(K)...
1579  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1580  dqdz(k) = 0.5*(tke(k) - tke(k+1))*idz(k)
1581  tke_src = dz_int(k) * (((kappa(k) + kappa0)*s2(k) - kappa(k)*n2(k)) - &
1582  (tke(k) - q0)*tke_decay(k)) - &
1583  (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k)))
1584  v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1585  v2 = (v1*dqmdk(k) + dqdz(k-1)*ck(k)) + &
1586  ((dqdz(k-1) - dqdz(k)) + dz_int(k)*(s2(k) - n2(k)))
1587 
1588  ! Ensure that the pivot is always positive, and that 0 <= cQ <= 1.
1589  ! Otherwise do not use Newton's method.
1590  decay_term_q = dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k) - v2*dkdq(k)
1591  if (decay_term_q < 0.0) then ; abort_newton = .true. ; exit ; endif
1592  bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
1593 
1594  cq(k+1) = aq(k) * bq
1595  cqcomp = (cqcomp*aq(k-1) + decay_term_q) * bq
1596  dqmdk(k+1) = (v2 * ck(k+1) - dqdz(k)) * bq
1597 
1598  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1599  dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + &
1600  (v2 * dk(k) + tke_src)), cqcomp*(-0.5*tke(k)))
1601 
1602  ! Check whether the next layer will be affected by any nonzero kappas.
1603  if ((itt > 1) .and. (k > ke_src) .and. (dk(k) == 0.0) .and. &
1604  ((kappa(k) + kappa(k+1)) == 0.0)) then
1605  ! Could also do .and. (bQ*abs(tke_src) < roundoff*TKE(K)) then
1606  ke_kappa = k-1 ; exit
1607  endif
1608  enddo
1609  if ((ke_kappa == nz) .and. (.not. abort_newton)) then
1610  dk(nz+1) = 0.0 ; dkdq(nz+1) = 0.0
1611  if (tke_noflux_bottom_bc) then
1612  k = nz+1
1613  tke_src = dz_int(k) * (kappa0*s2(k) - (tke(k) - q0)*tke_decay(k)) + &
1614  aq(k-1) * (tke(k-1) - tke(k))
1615 
1616  v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1617  decay_term_q = max(0.0, dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k))
1618  if (decay_term_q < 0.0) then
1619  abort_newton = .true.
1620  else
1621  bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
1622  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1623  dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + tke_src), -0.5*tke(k))
1624  tke(k) = max(tke(k) + dq(k), tke_min)
1625  endif
1626  else
1627  dq(nz+1) = 0.0
1628  endif
1629  elseif (.not. abort_newton) then
1630  ! Alter the first-guess determination of dQ(K).
1631  dq(ke_kappa+1) = dq(ke_kappa+1) / (1.0 - cq(ke_kappa+2)*e1(ke_kappa+2))
1632  tke(ke_kappa+1) = max(tke(ke_kappa+1) + dq(ke_kappa+1), tke_min)
1633  do k=ke_kappa+2,nz+1
1634  if (debug_soln .and. (k < nz+1)) then
1635  ! Ignore this source?
1636  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1637  ! tke_src_norm = (dz_Int(K) * (kappa0*S2(K) - (TKE(K)-q0)*TKE_decay(K)) - &
1638  ! (aQ(k) * (TKE(K) - TKE(K+1)) - aQ(k-1) * (TKE(K-1) - TKE(K))) ) / &
1639  ! (aQ(k) + (aQ(k-1) + dz_Int(K)*TKE_decay(K)))
1640  endif
1641  dk(k) = 0.0
1642  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1643  dq(k) = max(e1(k)*dq(k-1),-0.5*tke(k))
1644  tke(k) = max(tke(k) + dq(k), tke_min)
1645  if (abs(dq(k)) < roundoff*tke(k)) exit
1646  enddo
1647  if (debug_soln) then ; do k2=k+1,nz+1 ; dq(k2) = 0.0 ; dk(k2) = 0.0 ; enddo ; endif
1648  endif
1649  if (.not. abort_newton) then
1650  do k=ke_kappa,2,-1
1651  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1652  dq(k) = max(dq(k) + (cq(k+1)*dq(k+1) + dqmdk(k+1) * dk(k+1)), -0.5*tke(k))
1653  tke(k) = max(tke(k) + dq(k), tke_min)
1654  dk(k) = dk(k) + (ck(k+1)*dk(k+1) + dkdq(k) * dq(k))
1655  ! Truncate away negligibly small values of kappa.
1656  if (dk(k) <= kappa_trunc - kappa(k)) then
1657  dk(k) = -kappa(k)
1658  kappa(k) = 0.0
1659  if ((k < ks_src) .and. (k+1 > ks_kappa)) ks_kappa = k+1
1660  elseif (dk(k) < 2.0*kappa_trunc - kappa(k)) then
1661  dk(k) = 2.0*dk(k) - (2.0*kappa_trunc - kappa(k))
1662  kappa(k) = max(kappa(k) + dk(k), 0.0) ! The max is for paranoia.
1663  if (k<=ks_kappa) ks_kappa = 2
1664  else
1665  kappa(k) = kappa(k) + dk(k)
1666  if (k<=ks_kappa) ks_kappa = 2
1667  endif
1668  enddo
1669  dq(1) = max(dq(1) + cq(2)*dq(2) + dqmdk(2) * dk(2), tke_min - tke(1))
1670  tke(1) = max(tke(1) + dq(1), tke_min)
1671  dk(1) = 0.0
1672  endif
1673 
1674  ! Check these solutions for consistency.
1675  ! The unit conversions here have not been carefully tested.
1676  if (debug_soln) then ; do k=2,nz
1677  ! In these equations, K_err_lin and Q_err_lin should be at round-off levels
1678  ! compared with the dominant terms, perhaps, dz_Int*I_Ld2*kappa and
1679  ! dz_Int*TKE_decay*TKE. The exception is where, either 1) the decay term has been
1680  ! been increased to ensure a positive pivot, or 2) negative TKEs have been
1681  ! truncated, or 3) small or negative kappas have been rounded toward 0.
1682  i_q = 1.0 / tke(k)
1683  i_ld2_debug(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1684 
1685  kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa_prev(k)) + &
1686  (idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - &
1687  idz(k)*(kappa_prev(k)-kappa_prev(k+1)))
1688  k_err_lin = -idz(k-1)*(dk(k-1)-dk(k)) + idz(k)*(dk(k)-dk(k+1)) + &
1689  dz_int(k)*i_ld2_debug(k)*dk(k) - kap_src - &
1690  dz_int(k)*(n2(k)*ilambda2 + f2)*i_q**2*kappa_prev(k) * dq(k)
1691 
1692  tke_src = dz_int(k) * ((kappa_prev(k) + kappa0)*s2(k) - &
1693  kappa_prev(k)*n2(k) - (tke_prev(k) - q0)*tke_decay(k)) - &
1694  (aq(k) * (tke_prev(k) - tke_prev(k+1)) - aq(k-1) * (tke_prev(k-1) - tke_prev(k)))
1695  q_err_lin = tke_src + (aq(k-1) * (dq(k-1)-dq(k)) - aq(k) * (dq(k)-dq(k+1))) - &
1696  0.5*(tke_prev(k)-tke_prev(k+1))*idz(k) * (dk(k) + dk(k+1)) - &
1697  0.5*(tke_prev(k)-tke_prev(k-1))*idz(k-1)* (dk(k-1) + dk(k)) + &
1698  dz_int(k) * (dk(k) * (s2(k) - n2(k)) - dq(k)*tke_decay(k))
1699  enddo ; endif
1700 
1701  endif ! End of the Newton's method solver.
1702 
1703  ! Test kappa for convergence...
1704  if ((tol_err < newton_err) .and. (.not.abort_newton)) then
1705  ! A lower tolerance is used to switch to Newton's method than to switch back.
1706  newton_test = newton_err ; if (do_newton) newton_test = 2.0*newton_err
1707  was_newton = do_newton
1708  within_tolerance = .true. ; do_newton = .true.
1709  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1710  kappa_mean = kappa0 + (kappa(k) - 0.5*dk(k))
1711  if (abs(dk(k)) > newton_test * kappa_mean) then
1712  if (do_newton) abort_newton = .true.
1713  within_tolerance = .false. ; do_newton = .false. ; exit
1714  elseif (abs(dk(k)) > tol_err * kappa_mean) then
1715  within_tolerance = .false. ; if (.not.do_newton) exit
1716  endif
1717  if (abs(dq(k)) > newton_test*(tke(k) - 0.5*dq(k))) then
1718  if (do_newton) abort_newton = .true.
1719  do_newton = .false. ; if (.not.within_tolerance) exit
1720  endif
1721  enddo
1722 
1723  else ! Newton's method will not be used again, so no need to check.
1724  within_tolerance = .true.
1725  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1726  if (abs(dk(k)) > tol_err * (kappa0 + (kappa(k) - 0.5*dk(k)))) then
1727  within_tolerance = .false. ; exit
1728  endif
1729  enddo
1730  endif
1731 
1732  if (abort_newton) then
1733  do_newton = .false. ; abort_newton = .false.
1734  ! We went to Newton too quickly last time, so restrict the tolerance.
1735  newton_err = 0.5*newton_err
1736  ke_kappa_prev = nz
1737  do k=2,nz ; k_q(k) = kappa(k) / max(tke(k), tke_min) ; enddo
1738  endif
1739 
1740  if (within_tolerance) exit
1741 
1742  enddo
1743 
1744  if (do_newton) then ! K_Q needs to be calculated.
1745  do k=1,ks_kappa-1 ; k_q(k) = 0.0 ; enddo
1746  do k=ks_kappa,ke_kappa ; k_q(k) = kappa(k) / tke(k) ; enddo
1747  do k=ke_kappa+1,nz+1 ; k_q(k) = 0.0 ; enddo
1748  endif
1749 
1750  if (present(local_src)) then
1751  local_src(1) = 0.0 ; local_src(nz+1) = 0.0
1752  do k=2,nz
1753  diffusive_src = idz(k-1)*(kappa(k-1)-kappa(k)) + idz(k)*(kappa(k+1)-kappa(k))
1754  chg_by_k0 = kappa0 * ((idz(k-1)+idz(k)) / dz_int(k) + i_ld2(k))
1755  if (diffusive_src <= 0.0) then
1756  local_src(k) = k_src(k) + chg_by_k0
1757  else
1758  local_src(k) = (k_src(k) + chg_by_k0) + diffusive_src / dz_int(k)
1759  endif
1760  enddo
1761  endif
1762  if (present(kappa_src)) then
1763  kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
1764  do k=2,nz
1765  kappa_src(k) = k_src(k)
1766  enddo
1767  endif
1768 

◆ kappa_shear_at_vertex()

logical function, public mom_kappa_shear::kappa_shear_at_vertex ( type(param_file_type), intent(in)  param_file)

This function indicates to other modules whether the Jackson et al shear mixing parameterization will be used at the vertices without needing to duplicate the log entry. It returns false if the Jackson et al scheme is not used or if it is used via calculations at the tracer points.

Parameters
[in]param_fileA structure to parse for run-time parameters

Definition at line 1974 of file MOM_kappa_shear.F90.

1975  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1976 
1977  ! Local variables
1978  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
1979  logical :: do_kappa_shear
1980  ! This function returns true only if the parameters "USE_JACKSON_PARAM" and "VERTEX_SHEAR" are both true.
1981 
1982  kappa_shear_at_vertex = .false.
1983 
1984  call get_param(param_file, mdl, "USE_JACKSON_PARAM", do_kappa_shear, &
1985  default=.false., do_not_log=.true.)
1986  if (do_kappa_shear) &
1987  call get_param(param_file, mdl, "VERTEX_SHEAR", kappa_shear_at_vertex, &
1988  "If true, do the calculations of the shear-driven mixing "//&
1989  "at the cell vertices (i.e., the vorticity points).", &
1990  default=.false., do_not_log=.true.)
1991 

◆ kappa_shear_column()

subroutine mom_kappa_shear::kappa_shear_column ( real, dimension( gv %ke+1), intent(inout)  kappa,
real, dimension( gv %ke+1), intent(out)  tke,
real, intent(in)  dt,
integer, intent(in)  nzc,
real, intent(in)  f2,
real, intent(in)  surface_pres,
real, dimension( gv %ke), intent(in)  dz,
real, dimension( gv %ke), intent(in)  u0xdz,
real, dimension( gv %ke), intent(in)  v0xdz,
real, dimension( gv %ke), intent(in)  T0xdz,
real, dimension( gv %ke), intent(in)  S0xdz,
real, dimension( gv %ke+1), intent(out)  kappa_avg,
real, dimension( gv %ke+1), intent(out)  tke_avg,
type(thermo_var_ptrs), intent(in)  tv,
type(kappa_shear_cs), pointer  CS,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( gv %ke+1), intent(out), optional  I_Ld2_1d,
real, dimension( gv %ke+1), intent(out), optional  dz_Int_1d 
)
private

This subroutine calculates shear-driven diffusivity and TKE in a single column.

Parameters
[in]gvThe ocean's vertical grid structure.
[in,out]kappaThe time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
[out]tkeThe Turbulent Kinetic Energy per unit mass at
[in]nzcThe number of active layers in the column.
[in]f2The square of the Coriolis parameter [T-2 ~> s-2].
[in]surface_presThe surface pressure [R L2 T-2 ~> Pa].
[in]dzThe layer thickness [Z ~> m].
[in]u0xdzThe initial zonal velocity times dz [Z L T-1 ~> m2 s-1].
[in]v0xdzThe initial meridional velocity times dz [Z L T-1 ~> m2 s-1].
[in]t0xdzThe initial temperature times dz [degC Z ~> degC m].
[in]s0xdzThe initial salinity times dz [ppt Z ~> ppt m].
[out]kappa_avgThe time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
[out]tke_avgThe time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
[in]dtTime increment [T ~> s].
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
csThe control structure returned by a previous call to kappa_shear_init.
[in]usA dimensional unit scaling type
[out]i_ld2_1dThe inverse of the squared mixing length [Z-2 ~> m-2].
[out]dz_int_1dThe extent of a finite-volume space surrounding an interface,

Definition at line 624 of file MOM_kappa_shear.F90.

627  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
628  real, dimension(SZK_(GV)+1), &
629  intent(inout) :: kappa !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
630  real, dimension(SZK_(GV)+1), &
631  intent(out) :: tke !< The Turbulent Kinetic Energy per unit mass at
632  !! an interface [Z2 T-2 ~> m2 s-2].
633  integer, intent(in) :: nzc !< The number of active layers in the column.
634  real, intent(in) :: f2 !< The square of the Coriolis parameter [T-2 ~> s-2].
635  real, intent(in) :: surface_pres !< The surface pressure [R L2 T-2 ~> Pa].
636  real, dimension(SZK_(GV)), &
637  intent(in) :: dz !< The layer thickness [Z ~> m].
638  real, dimension(SZK_(GV)), &
639  intent(in) :: u0xdz !< The initial zonal velocity times dz [Z L T-1 ~> m2 s-1].
640  real, dimension(SZK_(GV)), &
641  intent(in) :: v0xdz !< The initial meridional velocity times dz [Z L T-1 ~> m2 s-1].
642  real, dimension(SZK_(GV)), &
643  intent(in) :: T0xdz !< The initial temperature times dz [degC Z ~> degC m].
644  real, dimension(SZK_(GV)), &
645  intent(in) :: S0xdz !< The initial salinity times dz [ppt Z ~> ppt m].
646  real, dimension(SZK_(GV)+1), &
647  intent(out) :: kappa_avg !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
648  real, dimension(SZK_(GV)+1), &
649  intent(out) :: tke_avg !< The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
650  real, intent(in) :: dt !< Time increment [T ~> s].
651  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
652  !! available thermodynamic fields. Absent fields
653  !! have NULL ptrs.
654  type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
655  !! call to kappa_shear_init.
656  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
657  real, dimension(SZK_(GV)+1), &
658  optional, intent(out) :: I_Ld2_1d !< The inverse of the squared mixing length [Z-2 ~> m-2].
659  real, dimension(SZK_(GV)+1), &
660  optional, intent(out) :: dz_Int_1d !< The extent of a finite-volume space surrounding an interface,
661  !! as used in calculating kappa and TKE [Z ~> m].
662 
663  ! Local variables
664  real, dimension(nzc) :: &
665  u, & ! The zonal velocity after a timestep of mixing [L T-1 ~> m s-1].
666  v, & ! The meridional velocity after a timestep of mixing [L T-1 ~> m s-1].
667  Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
668  T, & ! The potential temperature after a timestep of mixing [degC].
669  Sal, & ! The salinity after a timestep of mixing [ppt].
670  u_test, v_test, & ! Temporary velocities [L T-1 ~> m s-1].
671  T_test, S_test ! Temporary temperatures [degC] and salinities [ppt].
672 
673  real, dimension(nzc+1) :: &
674  N2, & ! The squared buoyancy frequency at an interface [T-2 ~> s-2].
675  dz_Int, & ! The extent of a finite-volume space surrounding an interface,
676  ! as used in calculating kappa and TKE [Z ~> m].
677  i_dz_int, & ! The inverse of the distance between velocity & density points
678  ! above and below an interface [Z-1 ~> m-1]. This is used to
679  ! calculate N2, shear, and fluxes, and it might differ from
680  ! 1/dz_Int, as they have different uses.
681  s2, & ! The squared shear at an interface [T-2 ~> s-2].
682  a1, & ! a1 is the coupling between adjacent interfaces in the TKE,
683  ! velocity, and density equations [Z s-1 ~> m s-1] or [Z ~> m]
684  c1, & ! c1 is used in the tridiagonal (and similar) solvers.
685  k_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
686  kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
687  kappa_out, & ! The kappa that results from the kappa equation [Z2 T-1 ~> m2 s-1].
688  kappa_mid, & ! The average of the initial and predictor estimates of kappa [Z2 T-1 ~> m2 s-1].
689  tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2].
690  kappa_pred, & ! The value of kappa from a predictor step [Z2 T-1 ~> m2 s-1].
691  pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].
692  t_int, & ! The temperature interpolated to an interface [degC].
693  sal_int, & ! The salinity interpolated to an interface [ppt].
694  dbuoy_dt, & ! The partial derivatives of buoyancy with changes in temperature
695  dbuoy_ds, & ! and salinity, [Z T-2 degC-1 ~> m s-2 degC-1] and [Z T-2 ppt-1 ~> m s-2 ppt-1].
696  i_l2_bdry, & ! The inverse of the square of twice the harmonic mean
697  ! distance to the top and bottom boundaries [Z-2 ~> m-2].
698  k_q, & ! Diffusivity divided by TKE [T ~> s].
699  k_q_tmp, & ! A temporary copy of diffusivity divided by TKE [T ~> s].
700  local_src_avg, & ! The time-integral of the local source [nondim].
701  tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1].
702  tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1].
703  tol_chg, & ! The tolerated kappa change integrated over a timestep [nondim].
704  dist_from_top, & ! The distance from the top surface [Z ~> m].
705  local_src ! The sum of all sources of kappa, including kappa_src and
706  ! sources from the elliptic term [T-1 ~> s-1].
707 
708  real :: dist_from_bot ! The distance from the bottom surface [Z ~> m].
709  real :: b1 ! The inverse of the pivot in the tridiagonal equations.
710  real :: bd1 ! A term in the denominator of b1.
711  real :: d1 ! 1 - c1 in the tridiagonal equations.
712  real :: gR0 ! A conversion factor from Z to pressure, given by Rho_0 times g
713  ! [R L2 T-2 Z-1 ~> kg m-2 s-2].
714  real :: g_R0 ! g_R0 is a rescaled version of g/Rho [Z R-1 T-2 ~> m4 kg-1 s-2].
715  real :: Norm ! A factor that normalizes two weights to 1 [Z-2 ~> m-2].
716  real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration
717  ! relative to the local source [nondim]. This must be greater than 1.
718  real :: tol2 ! The tolerance for the change in the kappa source within an iteration
719  ! relative to the average local source over previous iterations [nondim].
720  real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc
721  ! within an iteration [nondim]. 0 < tol_dksrc_low < 1.
722  real :: Ri_crit ! The critical shear Richardson number for shear-
723  ! driven mixing [nondim]. The theoretical value is 0.25.
724  real :: dt_rem ! The remaining time to advance the solution [T ~> s].
725  real :: dt_now ! The time step used in the current iteration [T ~> s].
726  real :: dt_wt ! The fractional weight of the current iteration [nondim].
727  real :: dt_test ! A time-step that is being tested for whether it
728  ! gives acceptably small changes in k_src [T ~> s].
729  real :: Idtt ! Idtt = 1 / dt_test [T-1 ~> s-1].
730  real :: dt_inc ! An increment to dt_test that is being tested [T ~> s].
731 
732  real :: k0dt ! The background diffusivity times the timestep [Z2 ~> m2].
733  logical :: valid_dt ! If true, all levels so far exhibit acceptably small changes in k_src.
734  logical :: use_temperature ! If true, temperature and salinity have been
735  ! allocated and are being used as state variables.
736  integer :: ks_kappa, ke_kappa ! The k-range with nonzero kappas.
737  integer :: dt_halvings ! The number of times that the time-step is halved
738  ! in seeking an acceptable timestep. If none is
739  ! found, dt_rem*0.5^dt_halvings is used.
740  integer :: dt_refinements ! The number of 2-fold refinements that will be used
741  ! to estimate the maximum permitted time step. I.e.,
742  ! the resolution is 1/2^dt_refinements.
743  integer :: k, itt, itt_dt
744 
745  ! This calculation of N2 is for debugging only.
746  ! real, dimension(SZK_(GV)+1) :: &
747  ! N2_debug, & ! A version of N2 for debugging [T-2 ~> s-2]
748 
749  ri_crit = cs%Rino_crit
750  gr0 = gv%Rho0 * gv%g_Earth
751  g_r0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
752  k0dt = dt*cs%kappa_0
753 
754  tol_dksrc = cs%kappa_src_max_chg
755  if (tol_dksrc == 10.0) then
756  ! This is equivalent to the expression below, but avoids changes at roundoff for the default value.
757  tol_dksrc_low = 0.95
758  else
759  tol_dksrc_low = (tol_dksrc - 0.5)/tol_dksrc
760  endif
761  tol2 = 2.0*cs%kappa_tol_err
762  dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
763  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
764 
765 
766  ! Set up Idz as the inverse of layer thicknesses.
767  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
768  ! Set up I_dz_int as the inverse of the distance between
769  ! adjacent layer centers.
770  i_dz_int(1) = 2.0 / dz(1)
771  dist_from_top(1) = 0.0
772  do k=2,nzc
773  i_dz_int(k) = 2.0 / (dz(k-1) + dz(k))
774  dist_from_top(k) = dist_from_top(k-1) + dz(k-1)
775  enddo
776  i_dz_int(nzc+1) = 2.0 / dz(nzc)
777 
778  ! Determine the velocities and thicknesses after eliminating massless
779  ! layers and applying a time-step of background diffusion.
780  if (nzc > 1) then
781  a1(2) = k0dt*i_dz_int(2)
782  b1 = 1.0 / (dz(1) + a1(2))
783  u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
784  t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
785  c1(2) = a1(2) * b1 ; d1 = dz(1) * b1 ! = 1 - c1
786  do k=2,nzc-1
787  bd1 = dz(k) + d1*a1(k)
788  a1(k+1) = k0dt*i_dz_int(k+1)
789  b1 = 1.0 / (bd1 + a1(k+1))
790  u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
791  v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
792  t(k) = b1 * (t0xdz(k) + a1(k)*t(k-1))
793  sal(k) = b1 * (s0xdz(k) + a1(k)*sal(k-1))
794  c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
795  enddo
796  ! rho or T and S have insulating boundary conditions, u & v use no-slip
797  ! bottom boundary conditions (if kappa0 > 0).
798  ! For no-slip bottom boundary conditions
799  b1 = 1.0 / ((dz(nzc) + d1*a1(nzc)) + k0dt*i_dz_int(nzc+1))
800  u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
801  v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
802  ! For insulating boundary conditions
803  b1 = 1.0 / (dz(nzc) + d1*a1(nzc))
804  t(nzc) = b1 * (t0xdz(nzc) + a1(nzc)*t(nzc-1))
805  sal(nzc) = b1 * (s0xdz(nzc) + a1(nzc)*sal(nzc-1))
806  do k=nzc-1,1,-1
807  u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
808  t(k) = t(k) + c1(k+1)*t(k+1) ; sal(k) = sal(k) + c1(k+1)*sal(k+1)
809  enddo
810  else
811  ! This is correct, but probably unnecessary.
812  b1 = 1.0 / (dz(1) + k0dt*i_dz_int(2))
813  u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
814  b1 = 1.0 / dz(1)
815  t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
816  endif
817 
818  ! This uses half the harmonic mean of thicknesses to provide two estimates
819  ! of the boundary between cells, and the inverse of the harmonic mean to
820  ! weight the two estimates. The net effect is that interfaces around thin
821  ! layers have thin cells, and the total thickness adds up properly.
822  ! The top- and bottom- interfaces have zero thickness, consistent with
823  ! adding additional zero thickness layers.
824  dz_int(1) = 0.0 ; dz_int(2) = dz(1)
825  do k=2,nzc-1
826  norm = 1.0 / (dz(k)*(dz(k-1)+dz(k+1)) + 2.0*dz(k-1)*dz(k+1))
827  dz_int(k) = dz_int(k) + dz(k) * ( ((dz(k)+dz(k+1)) * dz(k-1)) * norm)
828  dz_int(k+1) = dz(k) * ( ((dz(k-1)+dz(k)) * dz(k+1)) * norm)
829  enddo
830  dz_int(nzc) = dz_int(nzc) + dz(nzc) ; dz_int(nzc+1) = 0.0
831 
832  dist_from_bot = 0.0
833  do k=nzc,2,-1
834  dist_from_bot = dist_from_bot + dz(k)
835  i_l2_bdry(k) = (dist_from_top(k) + dist_from_bot)**2 / &
836  (dist_from_top(k) * dist_from_bot)**2
837  enddo
838 
839  ! Calculate thermodynamic coefficients and an initial estimate of N2.
840  if (use_temperature) then
841  pressure(1) = surface_pres
842  do k=2,nzc
843  pressure(k) = pressure(k-1) + gr0*dz(k-1)
844  t_int(k) = 0.5*(t(k-1) + t(k))
845  sal_int(k) = 0.5*(sal(k-1) + sal(k))
846  enddo
847  call calculate_density_derivs(t_int, sal_int, pressure, dbuoy_dt, dbuoy_ds, &
848  tv%eqn_of_state, (/2,nzc/), scale=-g_r0 )
849  else
850  do k=1,nzc+1 ; dbuoy_dt(k) = -g_r0 ; dbuoy_ds(k) = 0.0 ; enddo
851  endif
852 
853  ! N2_debug(1) = 0.0 ; N2_debug(nzc+1) = 0.0
854  ! do K=2,nzc
855  ! N2_debug(K) = max((dbuoy_dT(K) * (T0xdz(k-1)*Idz(k-1) - T0xdz(k)*Idz(k)) + &
856  ! dbuoy_dS(K) * (S0xdz(k-1)*Idz(k-1) - S0xdz(k)*Idz(k))) * &
857  ! I_dz_int(K), 0.0)
858  ! enddo
859 
860  ! This call just calculates N2 and S2.
861  call calculate_projected_state(kappa, u, v, t, sal, 0.0, nzc, dz, i_dz_int, &
862  dbuoy_dt, dbuoy_ds, u, v, t, sal, gv, us, &
863  n2=n2, s2=s2, vel_underflow=cs%vel_underflow)
864 ! ----------------------------------------------------
865 ! Iterate
866 ! ----------------------------------------------------
867  dt_rem = dt
868  do k=1,nzc+1
869  k_q(k) = 0.0
870  kappa_avg(k) = 0.0 ; tke_avg(k) = 0.0
871  local_src_avg(k) = 0.0
872  ! Use the grid spacings to scale errors in the source.
873  if ( dz_int(k) > 0.0 ) &
874  local_src_avg(k) = 0.1 * k0dt * i_dz_int(k) / dz_int(k)
875  enddo
876 
877 ! call cpu_clock_end(id_clock_setup)
878 
879 ! do itt=1,CS%max_RiNo_it
880  do itt=1,cs%max_KS_it
881 
882 ! ----------------------------------------------------
883 ! Calculate new values of u, v, rho, N^2 and S.
884 ! ----------------------------------------------------
885 
886  ! call cpu_clock_begin(id_clock_KQ)
887  call find_kappa_tke(n2, s2, kappa, idz, dz_int, i_l2_bdry, f2, &
888  nzc, cs, gv, us, k_q, tke, kappa_out, kappa_src, local_src)
889  ! call cpu_clock_end(id_clock_KQ)
890 
891  ! call cpu_clock_begin(id_clock_avg)
892  ! Determine the range of non-zero values of kappa_out.
893  ks_kappa = gv%ke+1 ; ke_kappa = 0
894  do k=2,nzc ; if (kappa_out(k) > 0.0) then
895  ks_kappa = k ; exit
896  endif ; enddo
897  do k=nzc,ks_kappa,-1 ; if (kappa_out(k) > 0.0) then
898  ke_kappa = k ; exit
899  endif ; enddo
900  if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
901  ! call cpu_clock_end(id_clock_avg)
902 
903  ! Determine how long to use this value of kappa (dt_now).
904 
905  ! call cpu_clock_begin(id_clock_project)
906  if ((ke_kappa < ks_kappa) .or. (itt==cs%max_KS_it)) then
907  dt_now = dt_rem
908  else
909  ! Limit dt_now so that |k_src(k)-kappa_src(k)| < tol * local_src(k)
910  dt_test = dt_rem
911  do k=2,nzc
912  tol_max(k) = kappa_src(k) + tol_dksrc * local_src(k)
913  tol_min(k) = kappa_src(k) - tol_dksrc_low * local_src(k)
914  tol_chg(k) = tol2 * local_src_avg(k)
915  enddo
916 
917  do itt_dt=1,(cs%max_KS_it+1-itt)/2
918  ! The maximum number of times that the time-step is halved in
919  ! seeking an acceptable timestep is reduced with each iteration,
920  ! so that as the maximum number of iterations is approached, the
921  ! whole remaining timestep is used. Typically, an acceptable
922  ! timestep is found long before the minimum is reached, so the
923  ! value of max_KS_it may be unimportant, especially if it is large
924  ! enough.
925  call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*dt_test, nzc, dz, i_dz_int, &
926  dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
927  gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, &
928  vel_underflow=cs%vel_underflow)
929  valid_dt = .true.
930  idtt = 1.0 / dt_test
931  do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
932  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
933  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
934  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
935  if (cs%restrictive_tolerance_check) then
936  if ((k_src(k) > min(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
937  (k_src(k) < max(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
938  valid_dt = .false. ; exit
939  endif
940  else
941  if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
942  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
943  valid_dt = .false. ; exit
944  endif
945  endif
946  else
947  if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
948  valid_dt = .false. ; k_src(k) = 0.0 ; exit
949  endif
950  endif
951  enddo
952 
953  if (valid_dt) exit
954  dt_test = 0.5*dt_test
955  enddo
956  if ((dt_test < dt_rem) .and. valid_dt) then
957  dt_inc = 0.5*dt_test
958  do itt_dt=1,dt_refinements
959  call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*(dt_test+dt_inc), &
960  nzc, dz, i_dz_int, dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
961  gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, vel_underflow=cs%vel_underflow)
962  valid_dt = .true.
963  idtt = 1.0 / (dt_test+dt_inc)
964  do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
965  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
966  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
967  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
968  if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
969  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
970  valid_dt = .false. ; exit
971  endif
972  else
973  if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
974  valid_dt = .false. ; k_src(k) = 0.0 ; exit
975  endif
976  endif
977  enddo
978 
979  if (valid_dt) dt_test = dt_test + dt_inc
980  dt_inc = 0.5*dt_inc
981  enddo
982  else
983  dt_inc = 0.0
984  endif
985 
986  dt_now = min(dt_test*(1.0+cs%kappa_tol_err)+dt_inc, dt_rem)
987  do k=2,nzc
988  local_src_avg(k) = local_src_avg(k) + dt_now * local_src(k)
989  enddo
990  endif ! Are all the values of kappa_out 0?
991  ! call cpu_clock_end(id_clock_project)
992 
993  ! The state has already been projected forward. Now find new values of kappa.
994 
995  if (ke_kappa < ks_kappa) then
996  ! There is no mixing now, and will not be again.
997  ! call cpu_clock_begin(id_clock_avg)
998  dt_wt = dt_rem / dt ; dt_rem = 0.0
999  do k=1,nzc+1
1000  kappa_mid(k) = 0.0
1001  ! This would be here but does nothing.
1002  ! kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
1003  tke_avg(k) = tke_avg(k) + dt_wt*tke(k)
1004  enddo
1005  ! call cpu_clock_end(id_clock_avg)
1006  else
1007  ! call cpu_clock_begin(id_clock_project)
1008  call calculate_projected_state(kappa_out, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
1009  dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1010  gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1011  vel_underflow=cs%vel_underflow)
1012  ! call cpu_clock_end(id_clock_project)
1013 
1014  ! call cpu_clock_begin(id_clock_KQ)
1015  do k=1,nzc+1 ; k_q_tmp(k) = k_q(k) ; enddo
1016  call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1017  nzc, cs, gv, us, k_q_tmp, tke_pred, kappa_pred)
1018  ! call cpu_clock_end(id_clock_KQ)
1019 
1020  ks_kappa = gv%ke+1 ; ke_kappa = 0
1021  do k=1,nzc+1
1022  kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1023  if ((kappa_mid(k) > 0.0) .and. (k<ks_kappa)) ks_kappa = k
1024  if (kappa_mid(k) > 0.0) ke_kappa = k
1025  enddo
1026 
1027  ! call cpu_clock_begin(id_clock_project)
1028  call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
1029  dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1030  gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1031  vel_underflow=cs%vel_underflow)
1032  ! call cpu_clock_end(id_clock_project)
1033 
1034  ! call cpu_clock_begin(id_clock_KQ)
1035  call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1036  nzc, cs, gv, us, k_q, tke_pred, kappa_pred)
1037  ! call cpu_clock_end(id_clock_KQ)
1038 
1039  ! call cpu_clock_begin(id_clock_avg)
1040  dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
1041  do k=1,nzc+1
1042  kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1043  kappa_avg(k) = kappa_avg(k) + kappa_mid(k)*dt_wt
1044  tke_avg(k) = tke_avg(k) + dt_wt*0.5*(tke_pred(k) + tke(k))
1045  kappa(k) = kappa_pred(k) ! First guess for the next iteration.
1046  enddo
1047  ! call cpu_clock_end(id_clock_avg)
1048  endif
1049 
1050  if (dt_rem > 0.0) then
1051  ! Update the values of u, v, T, Sal, N2, and S2 for the next iteration.
1052  ! call cpu_clock_begin(id_clock_project)
1053  call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, &
1054  dz, i_dz_int, dbuoy_dt, dbuoy_ds, u, v, t, sal, &
1055  gv, us, n2, s2, vel_underflow=cs%vel_underflow)
1056  ! call cpu_clock_end(id_clock_project)
1057  endif
1058 
1059  if (dt_rem <= 0.0) exit
1060 
1061  enddo ! end itt loop
1062 
1063  if (present(i_ld2_1d)) then
1064  do k=1,gv%ke+1 ; i_ld2_1d(k) = 0.0 ; enddo
1065  do k=2,nzc ; if (tke(k) > 0.0) &
1066  i_ld2_1d(k) = i_l2_bdry(k) + (n2(k) / cs%lambda**2 + f2) / tke(k)
1067  enddo
1068  endif
1069  if (present(dz_int_1d)) then
1070  do k=1,nzc+1 ; dz_int_1d(k) = dz_int(k) ; enddo
1071  do k=nzc+2,gv%ke ; dz_int_1d(k) = 0.0 ; enddo
1072  endif
1073 

◆ kappa_shear_init()

logical function, public mom_kappa_shear::kappa_shear_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(kappa_shear_cs), pointer  CS 
)

This subroutine initializes the parameters that regulate shear-driven mixing.

Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagA structure that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure for this module
Returns
True if module is to be used, False otherwise

Definition at line 1772 of file MOM_kappa_shear.F90.

1773  type(time_type), intent(in) :: Time !< The current model time.
1774  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1775  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1776  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1777  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1778  !! parameters.
1779  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1780  !! output.
1781  type(Kappa_shear_CS), pointer :: CS !< A pointer that is set to point to the control
1782  !! structure for this module
1783  logical :: kappa_shear_init !< True if module is to be used, False otherwise
1784 
1785  ! Local variables
1786  logical :: merge_mixedlayer
1787  logical :: debug_shear
1788  logical :: just_read ! If true, this module is not used, so only read the parameters.
1789  ! This include declares and sets the variable "version".
1790 # include "version_variable.h"
1791  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
1792  real :: kappa_0_unscaled ! The value of kappa_0 in MKS units [m2 s-1]
1793  real :: KD_normal ! The KD of the main model, read here only as a parameter
1794  ! for setting the default of KD_SMOOTH in MKS units [m2 s-1]
1795 
1796  if (associated(cs)) then
1797  call mom_error(warning, "kappa_shear_init called with an associated "// &
1798  "control structure.")
1799  return
1800  endif
1801  allocate(cs)
1802 
1803  ! The Jackson-Hallberg-Legg shear mixing parameterization uses the following
1804  ! 6 nondimensional coefficients. That paper gives 3 best fit parameter sets.
1805  ! Ri_Crit Rate FRi_Curv K_buoy TKE_N TKE_Shear
1806  ! p1: 0.25 0.089 -0.97 0.82 0.24 0.14
1807  ! p2: 0.30 0.085 -0.94 0.86 0.26 0.13
1808  ! p3: 0.35 0.088 -0.89 0.81 0.28 0.12
1809  ! Future research will reveal how these should be modified to take
1810  ! subgridscale inhomogeneity into account.
1811 
1812 ! Set default, read and log parameters
1813  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, default=.false., do_not_log=.true.)
1814  call log_version(param_file, mdl, version, &
1815  "Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008", &
1816  log_to_all=.true., debugging=kappa_shear_init, all_default=.not.kappa_shear_init)
1817  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, &
1818  "If true, use the Jackson-Hallberg-Legg (JPO 2008) "//&
1819  "shear mixing parameterization.", default=.false.)
1820  just_read = .not.kappa_shear_init
1821  call get_param(param_file, mdl, "VERTEX_SHEAR", cs%KS_at_vertex, &
1822  "If true, do the calculations of the shear-driven mixing "//&
1823  "at the cell vertices (i.e., the vorticity points).", &
1824  default=.false., do_not_log=just_read)
1825  call get_param(param_file, mdl, "RINO_CRIT", cs%RiNo_crit, &
1826  "The critical Richardson number for shear mixing.", &
1827  units="nondim", default=0.25, do_not_log=just_read)
1828  call get_param(param_file, mdl, "SHEARMIX_RATE", cs%Shearmix_rate, &
1829  "A nondimensional rate scale for shear-driven entrainment. "//&
1830  "Jackson et al find values in the range of 0.085-0.089.", &
1831  units="nondim", default=0.089, do_not_log=just_read)
1832  call get_param(param_file, mdl, "MAX_RINO_IT", cs%max_RiNo_it, &
1833  "The maximum number of iterations that may be used to "//&
1834  "estimate the Richardson number driven mixing.", &
1835  units="nondim", default=50, do_not_log=just_read)
1836  call get_param(param_file, mdl, "KD", kd_normal, default=0.0, do_not_log=.true.)
1837  call get_param(param_file, mdl, "KD_KAPPA_SHEAR_0", cs%kappa_0, &
1838  "The background diffusivity that is used to smooth the "//&
1839  "density and shear profiles before solving for the "//&
1840  "diffusivities. The default is the greater of KD and 1e-7 m2 s-1.", &
1841  units="m2 s-1", default=max(kd_normal, 1.0e-7), scale=us%m2_s_to_Z2_T, &
1842  unscaled=kappa_0_unscaled, do_not_log=just_read)
1843  call get_param(param_file, mdl, "KD_TRUNC_KAPPA_SHEAR", cs%kappa_trunc, &
1844  "The value of shear-driven diffusivity that is considered negligible "//&
1845  "and is rounded down to 0. The default is 1% of KD_KAPPA_SHEAR_0.", &
1846  units="m2 s-1", default=0.01*kappa_0_unscaled, scale=us%m2_s_to_Z2_T, do_not_log=just_read)
1847  call get_param(param_file, mdl, "FRI_CURVATURE", cs%FRi_curvature, &
1848  "The nondimensional curvature of the function of the "//&
1849  "Richardson number in the kappa source term in the "//&
1850  "Jackson et al. scheme.", units="nondim", default=-0.97, do_not_log=just_read)
1851  call get_param(param_file, mdl, "TKE_N_DECAY_CONST", cs%C_N, &
1852  "The coefficient for the decay of TKE due to "//&
1853  "stratification (i.e. proportional to N*tke). "//&
1854  "The values found by Jackson et al. are 0.24-0.28.", &
1855  units="nondim", default=0.24, do_not_log=just_read)
1856 ! call get_param(param_file, mdl, "LAYER_KAPPA_STAGGER", CS%layer_stagger, &
1857 ! default=.false., do_not_log=just_read)
1858  call get_param(param_file, mdl, "TKE_SHEAR_DECAY_CONST", cs%C_S, &
1859  "The coefficient for the decay of TKE due to shear (i.e. "//&
1860  "proportional to |S|*tke). The values found by Jackson "//&
1861  "et al. are 0.14-0.12.", units="nondim", default=0.14, do_not_log=just_read)
1862  call get_param(param_file, mdl, "KAPPA_BUOY_SCALE_COEF", cs%lambda, &
1863  "The coefficient for the buoyancy length scale in the "//&
1864  "kappa equation. The values found by Jackson et al. are "//&
1865  "in the range of 0.81-0.86.", units="nondim", default=0.82, do_not_log=just_read)
1866  call get_param(param_file, mdl, "KAPPA_N_OVER_S_SCALE_COEF2", cs%lambda2_N_S, &
1867  "The square of the ratio of the coefficients of the "//&
1868  "buoyancy and shear scales in the diffusivity equation, "//&
1869  "Set this to 0 (the default) to eliminate the shear scale. "//&
1870  "This is only used if USE_JACKSON_PARAM is true.", &
1871  units="nondim", default=0.0, do_not_log=just_read)
1872  call get_param(param_file, mdl, "KAPPA_SHEAR_TOL_ERR", cs%kappa_tol_err, &
1873  "The fractional error in kappa that is tolerated. "//&
1874  "Iteration stops when changes between subsequent "//&
1875  "iterations are smaller than this everywhere in a "//&
1876  "column. The peak diffusivities usually converge most "//&
1877  "rapidly, and have much smaller errors than this.", &
1878  units="nondim", default=0.1, do_not_log=just_read)
1879  call get_param(param_file, mdl, "TKE_BACKGROUND", cs%TKE_bg, &
1880  "A background level of TKE used in the first iteration "//&
1881  "of the kappa equation. TKE_BACKGROUND could be 0.", &
1882  units="m2 s-2", default=0.0, scale=us%m_to_Z**2*us%T_to_s**2)
1883  call get_param(param_file, mdl, "KAPPA_SHEAR_ELIM_MASSLESS", cs%eliminate_massless, &
1884  "If true, massless layers are merged with neighboring "//&
1885  "massive layers in this calculation. The default is "//&
1886  "true and I can think of no good reason why it should "//&
1887  "be false. This is only used if USE_JACKSON_PARAM is true.", &
1888  default=.true., do_not_log=just_read)
1889  call get_param(param_file, mdl, "MAX_KAPPA_SHEAR_IT", cs%max_KS_it, &
1890  "The maximum number of iterations that may be used to "//&
1891  "estimate the time-averaged diffusivity.", units="nondim", &
1892  default=13, do_not_log=just_read)
1893  call get_param(param_file, mdl, "PRANDTL_TURB", cs%Prandtl_turb, &
1894  "The turbulent Prandtl number applied to shear instability.", &
1895  units="nondim", default=1.0, do_not_log=.true.)
1896  call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &
1897  "A negligibly small velocity magnitude below which velocity components are set "//&
1898  "to 0. A reasonable value might be 1e-30 m/s, which is less than an "//&
1899  "Angstrom divided by the age of the universe.", &
1900  units="m s-1", default=0.0, scale=us%m_s_to_L_T, do_not_log=just_read)
1901  call get_param(param_file, mdl, "KAPPA_SHEAR_MAX_KAP_SRC_CHG", cs%kappa_src_max_chg, &
1902  "The maximum permitted increase in the kappa source within an iteration relative "//&
1903  "to the local source; this must be greater than 1. The lower limit for the "//&
1904  "permitted fractional decrease is (1 - 0.5/kappa_src_max_chg). These limits "//&
1905  "could perhaps be made dynamic with an improved iterative solver.", &
1906  default=10.0, units="nondim", do_not_log=just_read)
1907 
1908  call get_param(param_file, mdl, "DEBUG", cs%debug, &
1909  "If true, write out verbose debugging data.", &
1910  default=.false., debuggingparam=.true., do_not_log=just_read)
1911  call get_param(param_file, mdl, "DEBUG_KAPPA_SHEAR", debug_shear, &
1912  "If true, write debugging data for the kappa-shear code.", &
1913  default=.false., debuggingparam=.true., do_not_log=.true.)
1914  if (debug_shear) cs%debug = .true.
1915  call get_param(param_file, mdl, "KAPPA_SHEAR_VERTEX_PSURF_BUG", cs%psurf_bug, &
1916  "If true, do a simple average of the cell surface pressures to get a pressure "//&
1917  "at the corner if VERTEX_SHEAR=True. Otherwise mask out any land points in "//&
1918  "the average.", default=.true., do_not_log=(just_read .or. (.not.cs%KS_at_vertex)))
1919 
1920  call get_param(param_file, mdl, "KAPPA_SHEAR_ITER_BUG", cs%dKdQ_iteration_bug, &
1921  "If true, use an older, dimensionally inconsistent estimate of the "//&
1922  "derivative of diffusivity with energy in the Newton's method iteration. "//&
1923  "The bug causes undercorrections when dz > 1 m.", default=.false., do_not_log=just_read)
1924  call get_param(param_file, mdl, "KAPPA_SHEAR_ALL_LAYER_TKE_BUG", cs%all_layer_TKE_bug, &
1925  "If true, report back the latest estimate of TKE instead of the time average "//&
1926  "TKE when there is mass in all layers. Otherwise always report the time "//&
1927  "averaged TKE, as is currently done when there are some massless layers.", &
1928  default=.false., do_not_log=just_read)
1929  call get_param(param_file, mdl, "USE_RESTRICTIVE_TOLERANCE_CHECK", cs%restrictive_tolerance_check, &
1930  "If true, uses the more restrictive tolerance check to determine if a timestep "//&
1931  "is acceptable for the KS_it outer iteration loop. False uses the original less "//&
1932  "restrictive check.", default=.false., do_not_log=just_read)
1933 ! id_clock_KQ = cpu_clock_id('Ocean KS kappa_shear', grain=CLOCK_ROUTINE)
1934 ! id_clock_avg = cpu_clock_id('Ocean KS avg', grain=CLOCK_ROUTINE)
1935 ! id_clock_project = cpu_clock_id('Ocean KS project', grain=CLOCK_ROUTINE)
1936 ! id_clock_setup = cpu_clock_id('Ocean KS setup', grain=CLOCK_ROUTINE)
1937 
1938  cs%nkml = 1
1939  if (gv%nkml>0) then
1940  call get_param(param_file, mdl, "KAPPA_SHEAR_MERGE_ML",merge_mixedlayer, &
1941  "If true, combine the mixed layers together before solving the "//&
1942  "kappa-shear equations.", default=.true., do_not_log=just_read)
1943  if (merge_mixedlayer) cs%nkml = gv%nkml
1944  endif
1945 
1946 ! Forego remainder of initialization if not using this scheme
1947  if (.not. kappa_shear_init) return
1948 
1949  cs%diag => diag
1950 
1951  cs%id_Kd_shear = register_diag_field('ocean_model','Kd_shear', diag%axesTi, time, &
1952  'Shear-driven Diapycnal Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1953  cs%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesTi, time, &
1954  'Shear-driven Turbulent Kinetic Energy', 'm2 s-2', conversion=us%Z_to_m**2*us%s_to_T**2)
1955 

◆ kappa_shear_is_used()

logical function, public mom_kappa_shear::kappa_shear_is_used ( type(param_file_type), intent(in)  param_file)

This function indicates to other modules whether the Jackson et al shear mixing parameterization will be used without needing to duplicate the log entry.

Parameters
[in]param_fileA structure to parse for run-time parameters

Definition at line 1960 of file MOM_kappa_shear.F90.

1961  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1962 
1963  ! Local variables
1964  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
1965  ! This function reads the parameter "USE_JACKSON_PARAM" and returns its value.
1966 
1967  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_is_used, &
1968  default=.false., do_not_log=.true.)