MOM6
Data Types
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 subroutineinitializesthe 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 without needing to duplicate the log entry. 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 [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 364 of file MOM_kappa_shear.F90.

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

References kappa_shear_column().

Referenced by mom_set_diffusivity::set_diffusivity().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 [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 101 of file MOM_kappa_shear.F90.

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

References kappa_shear_column().

Here is the call graph for this function:

◆ 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 1216 of file MOM_kappa_shear.F90.

1216  integer, intent(in) :: nz !< The number of layers (after eliminating massless
1217  !! layers?).
1218  real, dimension(nz+1), intent(in) :: kappa !< The diapycnal diffusivity at interfaces,
1219  !! [Z2 T-1 ~> m2 s-1].
1220  real, dimension(nz), intent(in) :: u0 !< The initial zonal velocity [L T-1 ~> m s-1].
1221  real, dimension(nz), intent(in) :: v0 !< The initial meridional velocity [L T-1 ~> m s-1].
1222  real, dimension(nz), intent(in) :: T0 !< The initial temperature [degC].
1223  real, dimension(nz), intent(in) :: S0 !< The initial salinity [ppt].
1224  real, dimension(nz), intent(in) :: dz !< The grid spacing of layers [Z ~> m].
1225  real, dimension(nz+1), intent(in) :: I_dz_int !< The inverse of the layer's thicknesses
1226  !! [Z-1 ~> m-1].
1227  real, dimension(nz+1), intent(in) :: dbuoy_dT !< The partial derivative of buoyancy with
1228  !! temperature [Z T-2 degC-1 ~> m s-2 degC-1].
1229  real, dimension(nz+1), intent(in) :: dbuoy_dS !< The partial derivative of buoyancy with
1230  !! salinity [Z T-2 ppt-1 ~> m s-2 ppt-1].
1231  real, intent(in) :: dt !< The time step [T ~> s].
1232  real, dimension(nz), intent(inout) :: u !< The zonal velocity after dt [L T-1 ~> m s-1].
1233  real, dimension(nz), intent(inout) :: v !< The meridional velocity after dt [L T-1 ~> m s-1].
1234  real, dimension(nz), intent(inout) :: T !< The temperature after dt [degC].
1235  real, dimension(nz), intent(inout) :: Sal !< The salinity after dt [ppt].
1236  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1237  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1238  real, dimension(nz+1), optional, &
1239  intent(inout) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
1240  real, dimension(nz+1), optional, &
1241  intent(inout) :: S2 !< The squared shear at interfaces [T-2 ~> s-2].
1242  integer, optional, intent(in) :: ks_int !< The topmost k-index with a non-zero diffusivity.
1243  integer, optional, intent(in) :: ke_int !< The bottommost k-index with a non-zero
1244  !! diffusivity.
1245  real, optional, intent(in) :: vel_underflow !< If present and true, any velocities that
1246  !! are smaller in magnitude than this value are
1247  !! set to 0 [L T-1 ~> m s-1].
1248 
1249  ! Local variables
1250  real, dimension(nz+1) :: c1
1251  real :: L2_to_Z2 ! A conversion factor from horizontal length units to vertical depth
1252  ! units squared [Z2 s2 T-2 m-2 ~> 1].
1253  real :: underflow_vel ! Velocities smaller in magnitude than underflow_vel are set to 0 [L T-1 ~> m s-1].
1254  real :: a_a, a_b, b1, d1, bd1, b1nz_0
1255  integer :: k, ks, ke
1256 
1257  ks = 1 ; ke = nz
1258  if (present(ks_int)) ks = max(ks_int-1,1)
1259  if (present(ke_int)) ke = min(ke_int,nz)
1260  underflow_vel = 0.0 ; if (present(vel_underflow)) underflow_vel = vel_underflow
1261 
1262  if (ks > ke) return
1263 
1264  if (dt > 0.0) then
1265  a_b = dt*(kappa(ks+1)*i_dz_int(ks+1))
1266  b1 = 1.0 / (dz(ks) + a_b)
1267  c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1 ! = 1 - c1
1268 
1269  u(ks) = (b1 * dz(ks))*u0(ks) ; v(ks) = (b1 * dz(ks))*v0(ks)
1270  t(ks) = (b1 * dz(ks))*t0(ks) ; sal(ks) = (b1 * dz(ks))*s0(ks)
1271  do k=ks+1,ke-1
1272  a_a = a_b
1273  a_b = dt*(kappa(k+1)*i_dz_int(k+1))
1274  bd1 = dz(k) + d1*a_a
1275  b1 = 1.0 / (bd1 + a_b)
1276  c1(k+1) = a_b * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
1277 
1278  u(k) = b1 * (dz(k)*u0(k) + a_a*u(k-1))
1279  v(k) = b1 * (dz(k)*v0(k) + a_a*v(k-1))
1280  t(k) = b1 * (dz(k)*t0(k) + a_a*t(k-1))
1281  sal(k) = b1 * (dz(k)*s0(k) + a_a*sal(k-1))
1282  enddo
1283  ! T and S have insulating boundary conditions, u & v use no-slip
1284  ! bottom boundary conditions at the solid bottom.
1285 
1286  ! For insulating boundary conditions or mixing simply stopping, use...
1287  a_a = a_b
1288  b1 = 1.0 / (dz(ke) + d1*a_a)
1289  t(ke) = b1 * (dz(ke)*t0(ke) + a_a*t(ke-1))
1290  sal(ke) = b1 * (dz(ke)*s0(ke) + a_a*sal(ke-1))
1291 
1292  ! There is no distinction between the effective boundary conditions for
1293  ! tracers and velocities if the mixing is separated from the bottom, but if
1294  ! the mixing goes all the way to the bottom, use no-slip BCs for velocities.
1295  if (ke == nz) then
1296  a_b = dt*(kappa(nz+1)*i_dz_int(nz+1))
1297  b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
1298  else
1299  b1nz_0 = b1
1300  endif
1301  u(ke) = b1nz_0 * (dz(ke)*u0(ke) + a_a*u(ke-1))
1302  v(ke) = b1nz_0 * (dz(ke)*v0(ke) + a_a*v(ke-1))
1303  if (abs(u(ke)) < underflow_vel) u(ke) = 0.0
1304  if (abs(v(ke)) < underflow_vel) v(ke) = 0.0
1305 
1306  do k=ke-1,ks,-1
1307  u(k) = u(k) + c1(k+1)*u(k+1)
1308  v(k) = v(k) + c1(k+1)*v(k+1)
1309  if (abs(u(k)) < underflow_vel) u(k) = 0.0
1310  if (abs(v(k)) < underflow_vel) v(k) = 0.0
1311  t(k) = t(k) + c1(k+1)*t(k+1)
1312  sal(k) = sal(k) + c1(k+1)*sal(k+1)
1313  enddo
1314  else ! dt <= 0.0
1315  do k=1,nz
1316  u(k) = u0(k) ; v(k) = v0(k) ; t(k) = t0(k) ; sal(k) = s0(k)
1317  if (abs(u(k)) < underflow_vel) u(k) = 0.0
1318  if (abs(v(k)) < underflow_vel) v(k) = 0.0
1319  enddo
1320  endif
1321 
1322  if (present(s2)) then
1323  ! L2_to_Z2 = US%m_to_Z**2 * US%T_to_s**2
1324  l2_to_z2 = us%L_to_Z**2
1325  s2(1) = 0.0 ; s2(nz+1) = 0.0
1326  if (ks > 1) &
1327  s2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (l2_to_z2*i_dz_int(ks)**2)
1328  do k=ks+1,ke
1329  s2(k) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (l2_to_z2*i_dz_int(k)**2)
1330  enddo
1331  if (ke<nz) &
1332  s2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * (l2_to_z2*i_dz_int(ke+1)**2)
1333  endif
1334 
1335  if (present(n2)) then
1336  n2(1) = 0.0 ; n2(nz+1) = 0.0
1337  if (ks > 1) &
1338  n2(ks) = max(0.0, i_dz_int(ks) * &
1339  (dbuoy_dt(ks) * (t0(ks-1)-t(ks)) + dbuoy_ds(ks) * (s0(ks-1)-sal(ks))))
1340  do k=ks+1,ke
1341  n2(k) = max(0.0, i_dz_int(k) * &
1342  (dbuoy_dt(k) * (t(k-1)-t(k)) + dbuoy_ds(k) * (sal(k-1)-sal(k))))
1343  enddo
1344  if (ke<nz) &
1345  n2(ke+1) = max(0.0, i_dz_int(ke+1) * &
1346  (dbuoy_dt(ke+1) * (t(ke)-t0(ke+1)) + dbuoy_ds(ke+1) * (sal(ke)-s0(ke+1))))
1347  endif
1348 

Referenced by kappa_shear_column().

Here is the caller graph for this function:

◆ 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 [Z2 m-2 s2 T-1 ~> 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 1354 of file MOM_kappa_shear.F90.

1354  integer, intent(in) :: nz !< The number of layers to work on.
1355  real, dimension(nz+1), intent(in) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
1356  real, dimension(nz+1), intent(in) :: S2 !< The squared shear at interfaces [T-2 ~> s-2].
1357  real, dimension(nz+1), intent(in) :: kappa_in !< The initial guess at the diffusivity
1358  !! [Z2 T-1 ~> m2 s-1].
1359  real, dimension(nz+1), intent(in) :: dz_Int !< The thicknesses associated with interfaces
1360  !! [Z-1 ~> m-1].
1361  real, dimension(nz+1), intent(in) :: I_L2_bdry !< The inverse of the squared distance to
1362  !! boundaries [Z-2 ~> m-2].
1363  real, dimension(nz), intent(in) :: Idz !< The inverse grid spacing of layers [Z-1 ~> m-1].
1364  real, intent(in) :: f2 !< The squared Coriolis parameter [T-2 ~> s-2].
1365  type(Kappa_shear_CS), pointer :: CS !< A pointer to this module's control structure.
1366  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1367  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1368  real, dimension(nz+1), intent(inout) :: K_Q !< The shear-driven diapycnal diffusivity divided by
1369  !! the turbulent kinetic energy per unit mass at
1370  !! interfaces [Z2 m-2 s2 T-1 ~> s].
1371  real, dimension(nz+1), intent(out) :: tke !< The turbulent kinetic energy per unit mass at
1372  !! interfaces [Z2 T-2 ~> m2 s-2].
1373  real, dimension(nz+1), intent(out) :: kappa !< The diapycnal diffusivity at interfaces
1374  !! [Z2 T-1 ~> m2 s-1].
1375  real, dimension(nz+1), optional, &
1376  intent(out) :: kappa_src !< The source term for kappa [T-1 ~> s-1].
1377  real, dimension(nz+1), optional, &
1378  intent(out) :: local_src !< The sum of all local sources for kappa,
1379  !! [T-1 ~> s-1].
1380 ! This subroutine calculates new, consistent estimates of TKE and kappa.
1381 
1382  ! Local variables
1383  real, dimension(nz) :: &
1384  aQ, & ! aQ is the coupling between adjacent interfaces in the TKE equations [Z T-1 ~> m s-1].
1385  dQdz ! Half the partial derivative of TKE with depth [Z T-2 ~> m s-2].
1386  real, dimension(nz+1) :: &
1387  dK, & ! The change in kappa [Z2 T-1 ~> m2 s-1].
1388  dQ, & ! The change in TKE [Z2 T-2 ~> m2 s-2].
1389  cQ, cK, & ! cQ and cK are the upward influences in the tridiagonal and
1390  ! hexadiagonal solvers for the TKE and kappa equations [nondim].
1391  i_ld2, & ! 1/Ld^2, where Ld is the effective decay length scale for kappa [Z-2 ~> m-2].
1392  tke_decay, & ! The local TKE decay rate [T-1 ~> s-1].
1393  k_src, & ! The source term in the kappa equation [T-1 ~> s-1].
1394  dqmdk, & ! With Newton's method the change in dQ(k-1) due to dK(k) [T ~> s].
1395  dkdq, & ! With Newton's method the change in dK(k) due to dQ(k) [T-1 ~> s-1].
1396  e1 ! The fractional change in a layer TKE due to a change in the
1397  ! TKE of the layer above when all the kappas below are 0.
1398  ! e1 is nondimensional, and 0 < e1 < 1.
1399  real :: tke_src ! The net source of TKE due to mixing against the shear
1400  ! and stratification [Z2 T-3 ~> m2 s-3]. (For convenience,
1401  ! a term involving the non-dissipation of q0 is also
1402  ! included here.)
1403  real :: bQ ! The inverse of the pivot in the tridiagonal equations [T Z-1 ~> s m-1].
1404  real :: bK ! The inverse of the pivot in the tridiagonal equations [Z-1 ~> m-1].
1405  real :: bQd1 ! A term in the denominator of bQ [Z T-1 ~> m s-1].
1406  real :: bKd1 ! A term in the denominator of bK [Z ~> m].
1407  real :: cQcomp, cKcomp ! 1 - cQ or 1 - cK in the tridiagonal equations.
1408  real :: c_s2 ! The coefficient for the decay of TKE due to
1409  ! shear (i.e. proportional to |S|*tke), nondimensional.
1410  real :: c_n2 ! The coefficient for the decay of TKE due to
1411  ! stratification (i.e. proportional to N*tke) [nondim].
1412  real :: Ri_crit ! The critical shear Richardson number for shear-
1413  ! driven mixing. The theoretical value is 0.25.
1414  real :: q0 ! The background level of TKE [Z2 T-2 ~> m2 s-2].
1415  real :: Ilambda2 ! 1.0 / CS%lambda**2 [nondim]
1416  real :: TKE_min ! The minimum value of shear-driven TKE that can be
1417  ! solved for [Z2 T-2 ~> m2 s-2].
1418  real :: kappa0 ! The background diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
1419  real :: kappa_trunc ! Diffusivities smaller than this are rounded to 0 [Z2 T-1 ~> m2 s-1].
1420 
1421  real :: eden1, eden2, I_eden, ome ! Variables used in calculating e1.
1422  real :: diffusive_src ! The diffusive source in the kappa equation [Z T-1 ~> m s-1].
1423  real :: chg_by_k0 ! The value of k_src that leads to an increase of
1424  ! kappa_0 if only the diffusive term is a sink [T-1 ~> s-1].
1425 
1426  real :: kappa_mean ! A mean value of kappa [Z2 T-1 ~> m2 s-1].
1427  real :: Newton_test ! The value of relative error that will cause the next
1428  ! iteration to use Newton's method.
1429  ! Temporary variables used in the Newton's method iterations.
1430  real :: decay_term_k ! The decay term in the diffusivity equation
1431  real :: decay_term_Q ! The decay term in the TKE equation - proportional to [T-1 ~> s-1]
1432  real :: I_Q ! The inverse of TKE [T2 Z-2 ~> s2 m-2]
1433  real :: kap_src
1434  real :: v1 ! A temporary variable proportional to [T-1 ~> s-1]
1435  real :: v2
1436  real :: tol_err ! The tolerance for max_err that determines when to
1437  ! stop iterating.
1438  real :: Newton_err ! The tolerance for max_err that determines when to
1439  ! start using Newton's method. Empirically, an initial
1440  ! value of about 0.2 seems to be most efficient.
1441  real, parameter :: roundoff = 1.0e-16 ! A negligible fractional change in TKE.
1442  ! This could be larger but performance gains are small.
1443 
1444  logical :: tke_noflux_bottom_BC = .false. ! Specify the boundary conditions
1445  logical :: tke_noflux_top_BC = .false. ! that are applied to the TKE eqns.
1446  logical :: do_Newton ! If .true., use Newton's method for the next iteration.
1447  logical :: abort_Newton ! If .true., an Newton's method has encountered a 0
1448  ! pivot, and should not have been used.
1449  logical :: was_Newton ! The value of do_Newton before checking convergence.
1450  logical :: within_tolerance ! If .true., all points are within tolerance to
1451  ! enable this subroutine to return.
1452  integer :: ks_src, ke_src ! The range indices that have nonzero k_src.
1453  integer :: ks_kappa, ke_kappa, ke_tke ! The ranges of k-indices that are or
1454  integer :: ks_kappa_prev, ke_kappa_prev ! were being worked on.
1455  integer :: itt, k, k2
1456 #ifdef DEBUG
1457  integer :: max_debug_itt ; parameter(max_debug_itt=20)
1458  real :: K_err_lin, Q_err_lin, TKE_src_norm
1459  real, dimension(nz+1) :: &
1460  I_Ld2_debug, & ! A separate version of I_Ld2 for debugging [Z-2 ~> m-2].
1461  kappa_prev, & ! The value of kappa at the start of the current iteration [Z2 T-1 ~> m2 s-1].
1462  TKE_prev ! The value of TKE at the start of the current iteration [Z2 T-2 ~> m2 s-2].
1463  real, dimension(nz+1,1:max_debug_itt) :: &
1464  tke_it1, kappa_it1, kprev_it1, & ! Various values from each iteration.
1465  dkappa_it1, K_Q_it1, d_dkappa_it1, dkappa_norm_it1
1466  integer :: it2
1467 #endif
1468 
1469  c_n2 = cs%C_N**2 ; c_s2 = cs%C_S**2
1470  q0 = cs%TKE_bg ; kappa0 = cs%kappa_0
1471  tke_min = max(cs%TKE_bg, 1.0e-20*us%m_to_Z**2*us%T_to_s**2)
1472  ri_crit = cs%Rino_crit
1473  ilambda2 = 1.0 / cs%lambda**2
1474  kappa_trunc = cs%kappa_trunc
1475  do_newton = .false. ; abort_newton = .false.
1476  tol_err = cs%kappa_tol_err
1477  newton_err = 0.2 ! This initial value may be automatically reduced later.
1478 
1479  ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
1480 
1481  ke_src = 0 ; ks_src = nz+1
1482  do k=2,nz
1483  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
1484 ! Ri = N2(K) / S2(K)
1485 ! k_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1486 ! ((Ri_crit - Ri) / (Ri_crit + CS%FRi_curvature*Ri))
1487  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1488  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1489  ke_src = k
1490  if (ks_src > k) ks_src = k
1491  else
1492  k_src(k) = 0.0
1493  endif
1494  enddo
1495 
1496  ! If there is no source anywhere, return kappa(K) = 0.
1497  if (ks_src > ke_src) then
1498  do k=1,nz+1
1499  kappa(k) = 0.0 ; k_q(k) = 0.0 ; tke(k) = tke_min
1500  enddo
1501  if (present(kappa_src)) then ; do k=1,nz+1 ; kappa_src(k) = 0.0 ; enddo ; endif
1502  if (present(local_src)) then ; do k=1,nz+1 ; local_src(k) = 0.0 ; enddo ; endif
1503  return
1504  endif
1505 
1506  do k=1,nz+1
1507  kappa(k) = kappa_in(k)
1508 ! TKE_decay(K) = c_n*sqrt(N2(K)) + c_s*sqrt(S2(K)) ! The expression in JHL.
1509  tke_decay(k) = sqrt(c_n2*n2(k) + c_s2*s2(k))
1510  if ((kappa(k) > 0.0) .and. (k_q(k) > 0.0)) then
1511  tke(k) = kappa(k) / k_q(k)
1512  else
1513  tke(k) = tke_min
1514  endif
1515  enddo
1516  ! Apply boundary conditions to kappa.
1517  kappa(1) = 0.0 ; kappa(nz+1) = 0.0
1518 
1519  ! Calculate the term (e1) that allows changes in TKE to be calculated quickly
1520  ! below the deepest nonzero value of kappa. If kappa = 0, below interface
1521  ! k-1, the final changes in TKE are related by dQ(K+1) = e1(K+1)*dQ(K).
1522  eden2 = kappa0 * idz(nz)
1523  if (tke_noflux_bottom_bc) then
1524  eden1 = dz_int(nz+1)*tke_decay(nz+1)
1525  i_eden = 1.0 / (eden2 + eden1)
1526  e1(nz+1) = eden2 * i_eden ; ome = eden1 * i_eden
1527  else
1528  e1(nz+1) = 0.0 ; ome = 1.0
1529  endif
1530  do k=nz,2,-1
1531  eden1 = dz_int(k)*tke_decay(k) + ome * eden2
1532  eden2 = kappa0 * idz(k-1)
1533  i_eden = 1.0 / (eden2 + eden1)
1534  e1(k) = eden2 * i_eden ; ome = eden1 * i_eden ! = 1-e1
1535  enddo
1536  e1(1) = 0.0
1537 
1538 
1539  ! Iterate here to convergence to within some tolerance of order tol_err.
1540  do itt=1,cs%max_RiNo_it
1541 
1542  ! ----------------------------------------------------
1543  ! Calculate TKE
1544  ! ----------------------------------------------------
1545 
1546 #ifdef DEBUG
1547  do k=1,nz+1 ; kappa_prev(k) = kappa(k) ; tke_prev(k) = tke(k) ; enddo
1548 #endif
1549 
1550  if (.not.do_newton) then
1551  ! Use separate steps of the TKE and kappa equations, that are
1552  ! explicit in the nonlinear source terms, implicit in a linearized
1553  ! version of the nonlinear sink terms, and implicit in the linear
1554  ! terms.
1555 
1556  ke_tke = max(ke_kappa,ke_kappa_prev)+1
1557  ! aQ is the coupling between adjacent interfaces [Z T-1 ~> m s-1].
1558  do k=1,min(ke_tke,nz)
1559  aq(k) = (0.5*(kappa(k)+kappa(k+1)) + kappa0) * idz(k)
1560  enddo
1561  dq(1) = -tke(1)
1562  if (tke_noflux_top_bc) then
1563  tke_src = kappa0*s2(1) + q0 * tke_decay(1) ! Uses that kappa(1) = 0
1564  bqd1 = dz_int(1) * tke_decay(1)
1565  bq = 1.0 / (bqd1 + aq(1))
1566  tke(1) = bq * (dz_int(1)*tke_src)
1567  cq(2) = aq(1) * bq ; cqcomp = bqd1 * bq ! = 1 - cQ
1568  else
1569  tke(1) = q0 ; cq(2) = 0.0 ; cqcomp = 1.0
1570  endif
1571  do k=2,ke_tke-1
1572  dq(k) = -tke(k)
1573  tke_src = (kappa(k) + kappa0)*s2(k) + q0*tke_decay(k)
1574  bqd1 = dz_int(k)*(tke_decay(k) + n2(k)*k_q(k)) + cqcomp*aq(k-1)
1575  bq = 1.0 / (bqd1 + aq(k))
1576  tke(k) = bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1))
1577  cq(k+1) = aq(k) * bq ; cqcomp = bqd1 * bq ! = 1 - cQ
1578  enddo
1579  if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_bc)) then
1580  tke(nz+1) = tke_min
1581  dq(nz+1) = 0.0
1582  else
1583  k = ke_tke
1584  tke_src = kappa0*s2(k) + q0*tke_decay(k) ! Uses that kappa(ke_tke) = 0
1585  if (k == nz+1) then
1586  dq(k) = -tke(k)
1587  bq = 1.0 / (dz_int(k)*tke_decay(k) + cqcomp*aq(k-1))
1588  tke(k) = max(tke_min, bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)))
1589  dq(k) = tke(k) + dq(k)
1590  else
1591  bq = 1.0 / ((dz_int(k)*tke_decay(k) + cqcomp*aq(k-1)) + aq(k))
1592  cq(k+1) = aq(k) * bq
1593  ! Account for all changes deeper in the water column.
1594  dq(k) = -tke(k)
1595  tke(k) = max((bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)) + &
1596  cq(k+1)*(tke(k+1) - e1(k+1)*tke(k))) / (1.0 - cq(k+1)*e1(k+1)), tke_min)
1597  dq(k) = tke(k) + dq(k)
1598 
1599  ! Adjust TKE deeper in the water column in case ke_tke increases.
1600  ! This might not be strictly necessary?
1601  do k=ke_tke+1,nz+1
1602  dq(k) = e1(k)*dq(k-1)
1603  tke(k) = max(tke(k) + dq(k), tke_min)
1604  if (abs(dq(k)) < roundoff*tke(k)) exit
1605  enddo
1606  do k2=k+1,nz
1607  if (dq(k2) == 0.0) exit
1608  dq(k2) = 0.0
1609  enddo
1610  endif
1611  endif
1612  do k=ke_tke-1,1,-1
1613  tke(k) = max(tke(k) + cq(k+1)*tke(k+1), tke_min)
1614  dq(k) = tke(k) + dq(k)
1615  enddo
1616 
1617  ! ----------------------------------------------------
1618  ! Calculate kappa, here defined at interfaces.
1619  ! ----------------------------------------------------
1620 
1621  ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
1622 
1623  dk(1) = 0.0 ! kappa takes boundary values of 0.
1624  ck(2) = 0.0 ; ckcomp = 1.0
1625  if (itt == 1) then ; dO k=2,nz
1626  i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1627  enddo ; endif
1628  do k=2,nz
1629  dk(k) = -kappa(k)
1630  if (itt>1) &
1631  i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1632  bkd1 = dz_int(k)*i_ld2(k) + ckcomp*idz(k-1)
1633  bk = 1.0 / (bkd1 + idz(k))
1634 
1635  kappa(k) = bk * (idz(k-1)*kappa(k-1) + dz_int(k) * k_src(k))
1636  ck(k+1) = idz(k) * bk ; ckcomp = bkd1 * bk ! = 1 - cK(K+1)
1637 
1638  ! Neglect values that are smaller than kappa_trunc.
1639  if (kappa(k) < ckcomp*kappa_trunc) then
1640  kappa(k) = 0.0
1641  if (k > ke_src) then ; ke_kappa = k-1 ; k_q(k) = 0.0 ; exit ; endif
1642  elseif (kappa(k) < 2.0*ckcomp*kappa_trunc) then
1643  kappa(k) = 2.0 * (kappa(k) - ckcomp*kappa_trunc)
1644  endif
1645  enddo
1646  k_q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
1647  dk(ke_kappa) = dk(ke_kappa) + kappa(ke_kappa)
1648  do k=ke_kappa+2,ke_kappa_prev
1649  dk(k) = -kappa(k) ; kappa(k) = 0.0 ; k_q(k) = 0.0
1650  enddo
1651  do k=ke_kappa-1,2,-1
1652  kappa(k) = kappa(k) + ck(k+1)*kappa(k+1)
1653  ! Neglect values that are smaller than kappa_trunc.
1654  if (kappa(k) <= kappa_trunc) then
1655  kappa(k) = 0.0
1656  if (k < ks_src) then ; ks_kappa = k+1 ; k_q(k) = 0.0 ; exit ; endif
1657  elseif (kappa(k) < 2.0*kappa_trunc) then
1658  kappa(k) = 2.0 * (kappa(k) - kappa_trunc)
1659  endif
1660 
1661  dk(k) = dk(k) + kappa(k)
1662  k_q(k) = kappa(k) / tke(k)
1663  enddo
1664  do k=ks_kappa_prev,ks_kappa-2 ; kappa(k) = 0.0 ; k_q(k) = 0.0 ; enddo
1665 
1666  else ! do_Newton is .true.
1667 ! Once the solutions are close enough, use a Newton's method solver of the
1668 ! whole system to accelerate convergence.
1669  ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
1670  ks_kappa = 2
1671  dk(1) = 0.0 ; ck(2) = 0.0 ; ckcomp = 1.0 ; dkdq(1) = 0.0
1672  aq(1) = (0.5*(kappa(1)+kappa(2))+kappa0) * idz(1)
1673  dqdz(1) = 0.5*(tke(1) - tke(2))*idz(1)
1674  if (tke_noflux_top_bc) then
1675  tke_src = dz_int(1) * (kappa0*s2(1) - (tke(1) - q0)*tke_decay(1)) - &
1676  aq(1) * (tke(1) - tke(2))
1677 
1678  bq = 1.0 / (aq(1) + dz_int(1)*tke_decay(1))
1679  cq(2) = aq(1) * bq
1680  cqcomp = (dz_int(1)*tke_decay(1)) * bq ! = 1 - cQ(2)
1681  dqmdk(2) = -dqdz(1) * bq
1682  dq(1) = bq * tke_src
1683  else
1684  dq(1) = 0.0 ; cq(2) = 0.0 ; cqcomp = 1.0 ; dqmdk(2) = 0.0
1685  endif
1686  do k=2,nz
1687  i_q = 1.0 / tke(k)
1688  i_ld2(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1689 
1690  kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa(k)) + &
1691  idz(k-1)*(kappa(k-1)-kappa(k)) - idz(k)*(kappa(k)-kappa(k+1))
1692 
1693  ! Ensure that the pivot is always positive, and that 0 <= cK <= 1.
1694  ! Otherwise do not use Newton's method.
1695  decay_term_k = -idz(k-1)*dqmdk(k)*dkdq(k-1) + dz_int(k)*i_ld2(k)
1696  if (decay_term_k < 0.0) then ; abort_newton = .true. ; exit ; endif
1697  bk = 1.0 / (idz(k) + idz(k-1)*ckcomp + decay_term_k)
1698 
1699  ck(k+1) = bk * idz(k)
1700  ckcomp = bk * (idz(k-1)*ckcomp + decay_term_k) ! = 1-cK(K+1)
1701  if (cs%dKdQ_iteration_bug) then
1702  dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1703  us%m_to_Z*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
1704  else
1705  dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1706  dz_int(k)*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
1707  endif
1708  dk(k) = bk * (kap_src + idz(k-1)*dk(k-1) + idz(k-1)*dkdq(k-1)*dq(k-1))
1709 
1710  ! Truncate away negligibly small values of kappa.
1711  if (dk(k) <= ckcomp*(kappa_trunc - kappa(k))) then
1712  dk(k) = -ckcomp*kappa(k)
1713 ! if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1714  elseif (dk(k) < ckcomp*(2.0*kappa_trunc - kappa(k))) then
1715  dk(k) = 2.0 * dk(k) - ckcomp*(2.0*kappa_trunc - kappa(k))
1716  endif
1717 
1718  ! Solve for dQ(K)...
1719  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1720  dqdz(k) = 0.5*(tke(k) - tke(k+1))*idz(k)
1721  tke_src = dz_int(k) * (((kappa(k) + kappa0)*s2(k) - kappa(k)*n2(k)) - &
1722  (tke(k) - q0)*tke_decay(k)) - &
1723  (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k)))
1724  v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1725  v2 = (v1*dqmdk(k) + dqdz(k-1)*ck(k)) + &
1726  ((dqdz(k-1) - dqdz(k)) + dz_int(k)*(s2(k) - n2(k)))
1727 
1728  ! Ensure that the pivot is always positive, and that 0 <= cQ <= 1.
1729  ! Otherwise do not use Newton's method.
1730  decay_term_q = dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k) - v2*dkdq(k)
1731  if (decay_term_q < 0.0) then ; abort_newton = .true. ; exit ; endif
1732  bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
1733 
1734  cq(k+1) = aq(k) * bq
1735  cqcomp = (cqcomp*aq(k-1) + decay_term_q) * bq
1736  dqmdk(k+1) = (v2 * ck(k+1) - dqdz(k)) * bq
1737 
1738  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1739  dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + &
1740  (v2 * dk(k) + tke_src)), cqcomp*(-0.5*tke(k)))
1741 
1742  ! Check whether the next layer will be affected by any nonzero kappas.
1743  if ((itt > 1) .and. (k > ke_src) .and. (dk(k) == 0.0) .and. &
1744  ((kappa(k) + kappa(k+1)) == 0.0)) then
1745  ! Could also do .and. (bQ*abs(tke_src) < roundoff*TKE(K)) then
1746  ke_kappa = k-1 ; exit
1747  endif
1748  enddo
1749  if ((ke_kappa == nz) .and. (.not. abort_newton)) then
1750  dk(nz+1) = 0.0 ; dkdq(nz+1) = 0.0
1751  if (tke_noflux_bottom_bc) then
1752  k = nz+1
1753  tke_src = dz_int(k) * (kappa0*s2(k) - (tke(k) - q0)*tke_decay(k)) + &
1754  aq(k-1) * (tke(k-1) - tke(k))
1755 
1756  v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1757  decay_term_q = max(0.0, dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k))
1758  if (decay_term_q < 0.0) then
1759  abort_newton = .true.
1760  else
1761  bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
1762  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1763  dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + tke_src), -0.5*tke(k))
1764  tke(k) = max(tke(k) + dq(k), tke_min)
1765  endif
1766  else
1767  dq(nz+1) = 0.0
1768  endif
1769  elseif (.not. abort_newton) then
1770  ! Alter the first-guess determination of dQ(K).
1771  dq(ke_kappa+1) = dq(ke_kappa+1) / (1.0 - cq(ke_kappa+2)*e1(ke_kappa+2))
1772  tke(ke_kappa+1) = max(tke(ke_kappa+1) + dq(ke_kappa+1), tke_min)
1773  do k=ke_kappa+2,nz+1
1774 #ifdef DEBUG
1775  if (k < nz+1) then
1776  ! Ignore this source?
1777  aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1778  tke_src_norm = (dz_int(k) * (kappa0*s2(k) - (tke(k)-q0)*tke_decay(k)) - &
1779  (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k))) ) / &
1780  (aq(k) + (aq(k-1) + dz_int(k)*tke_decay(k)))
1781  endif
1782 #endif
1783  dk(k) = 0.0
1784  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1785  dq(k) = max(e1(k)*dq(k-1),-0.5*tke(k))
1786  tke(k) = max(tke(k) + dq(k), tke_min)
1787  if (abs(dq(k)) < roundoff*tke(k)) exit
1788  enddo
1789 #ifdef DEBUG
1790  do k2=k+1,ke_kappa_prev+1 ; dq(k2) = 0.0 ; dk(k2) = 0.0 ; enddo
1791  do k=k2,nz+1 ; if (dq(k) == 0.0) exit ; dq(k) = 0.0 ; dk(k) = 0.0 ; enddo
1792 #endif
1793  endif
1794  if (.not. abort_newton) then
1795  do k=ke_kappa,2,-1
1796  ! Ensure that TKE+dQ will not drop below 0.5*TKE.
1797  dq(k) = max(dq(k) + (cq(k+1)*dq(k+1) + dqmdk(k+1) * dk(k+1)), -0.5*tke(k))
1798  tke(k) = max(tke(k) + dq(k), tke_min)
1799  dk(k) = dk(k) + (ck(k+1)*dk(k+1) + dkdq(k) * dq(k))
1800  ! Truncate away negligibly small values of kappa.
1801  if (dk(k) <= kappa_trunc - kappa(k)) then
1802  dk(k) = -kappa(k)
1803  kappa(k) = 0.0
1804  if ((k < ks_src) .and. (k+1 > ks_kappa)) ks_kappa = k+1
1805  elseif (dk(k) < 2.0*kappa_trunc - kappa(k)) then
1806  dk(k) = 2.0*dk(k) - (2.0*kappa_trunc - kappa(k))
1807  kappa(k) = max(kappa(k) + dk(k), 0.0) ! The max is for paranoia.
1808  if (k<=ks_kappa) ks_kappa = 2
1809  else
1810  kappa(k) = kappa(k) + dk(k)
1811  if (k<=ks_kappa) ks_kappa = 2
1812  endif
1813  enddo
1814  dq(1) = max(dq(1) + cq(2)*dq(2) + dqmdk(2) * dk(2), tke_min - tke(1))
1815  tke(1) = max(tke(1) + dq(1), tke_min)
1816  dk(1) = 0.0
1817  endif
1818 
1819 #ifdef DEBUG
1820  ! Check these solutions for consistency.
1821  ! The unit conversions here have not been carefully tested.
1822  do k=2,nz
1823  ! In these equations, K_err_lin and Q_err_lin should be at round-off levels
1824  ! compared with the dominant terms, perhaps, dz_Int*I_Ld2*kappa and
1825  ! dz_Int*TKE_decay*TKE. The exception is where, either 1) the decay term has been
1826  ! been increased to ensure a positive pivot, or 2) negative TKEs have been
1827  ! truncated, or 3) small or negative kappas have been rounded toward 0.
1828  i_q = 1.0 / tke(k)
1829  i_ld2_debug(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1830 
1831  kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa_prev(k)) + &
1832  (idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - &
1833  idz(k)*(kappa_prev(k)-kappa_prev(k+1)))
1834  k_err_lin = -idz(k-1)*(dk(k-1)-dk(k)) + idz(k)*(dk(k)-dk(k+1)) + &
1835  dz_int(k)*i_ld2_debug(k)*dk(k) - kap_src - &
1836  dz_int(k)*(n2(k)*ilambda2 + f2)*i_q**2*kappa_prev(k) * dq(k)
1837 
1838  tke_src = dz_int(k) * ((kappa_prev(k) + kappa0)*s2(k) - &
1839  kappa_prev(k)*n2(k) - (tke_prev(k) - q0)*tke_decay(k)) - &
1840  (aq(k) * (tke_prev(k) - tke_prev(k+1)) - aq(k-1) * (tke_prev(k-1) - tke_prev(k)))
1841  q_err_lin = tke_src + (aq(k-1) * (dq(k-1)-dq(k)) - aq(k) * (dq(k)-dq(k+1))) - &
1842  0.5*(tke_prev(k)-tke_prev(k+1))*idz(k) * (dk(k) + dk(k+1)) - &
1843  0.5*(tke_prev(k)-tke_prev(k-1))*idz(k-1)* (dk(k-1) + dk(k)) + &
1844  dz_int(k) * (dk(k) * (s2(k) - n2(k)) - dq(k)*tke_decay(k))
1845  enddo
1846 #endif
1847  endif ! End of the Newton's method solver.
1848 
1849  ! Test kappa for convergence...
1850  if ((tol_err < newton_err) .and. (.not.abort_newton)) then
1851  ! A lower tolerance is used to switch to Newton's method than to switch back.
1852  newton_test = newton_err ; if (do_newton) newton_test = 2.0*newton_err
1853  was_newton = do_newton
1854  within_tolerance = .true. ; do_newton = .true.
1855  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1856  kappa_mean = kappa0 + (kappa(k) - 0.5*dk(k))
1857  if (abs(dk(k)) > newton_test * kappa_mean) then
1858  if (do_newton) abort_newton = .true.
1859  within_tolerance = .false. ; do_newton = .false. ; exit
1860  elseif (abs(dk(k)) > tol_err * kappa_mean) then
1861  within_tolerance = .false. ; if (.not.do_newton) exit
1862  endif
1863  if (abs(dq(k)) > newton_test*(tke(k) - 0.5*dq(k))) then
1864  if (do_newton) abort_newton = .true.
1865  do_newton = .false. ; if (.not.within_tolerance) exit
1866  endif
1867  enddo
1868 
1869  else ! Newton's method will not be used again, so no need to check.
1870  within_tolerance = .true.
1871  do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1872  if (abs(dk(k)) > tol_err * (kappa0 + (kappa(k) - 0.5*dk(k)))) then
1873  within_tolerance = .false. ; exit
1874  endif
1875  enddo
1876  endif
1877 
1878  if (abort_newton) then
1879  do_newton = .false. ; abort_newton = .false.
1880  ! We went to Newton too quickly last time, so restrict the tolerance.
1881  newton_err = 0.5*newton_err
1882  ke_kappa_prev = nz
1883  do k=2,nz ; k_q(k) = kappa(k) / max(tke(k), tke_min) ; enddo
1884  endif
1885 
1886 #ifdef DEBUG
1887  if (itt <= max_debug_itt) then
1888  do k=1,nz+1
1889  kprev_it1(k,itt) = kappa_prev(k)
1890  kappa_it1(k,itt) = kappa(k) ; tke_it1(k,itt) = tke(k)
1891  dkappa_it1(k,itt) = kappa(k) - kappa_prev(k)
1892  dkappa_norm_it1(k,itt) = (kappa(k) - kappa_prev(k)) / &
1893  (kappa0 + 0.5*(kappa(k) + kappa_prev(k)))
1894  k_q_it1(k,itt) = kappa(k) / max(tke(k),tke_min)
1895  d_dkappa_it1(k,itt) = 0.0
1896  if (itt > 1) then ; if (abs(dkappa_it1(k,itt-1)) > 1e-20*us%m2_s_to_Z2_T) &
1897  d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
1898  endif
1899  enddo
1900  endif
1901 #endif
1902 
1903  if (within_tolerance) exit
1904 
1905  enddo
1906 
1907 #ifdef DEBUG
1908  do it2=itt+1,max_debug_itt ; do k=1,nz+1
1909  kprev_it1(k,it2) = 0.0 ; kappa_it1(k,it2) = 0.0 ; tke_it1(k,it2) = 0.0
1910  dkappa_it1(k,it2) = 0.0 ; k_q_it1(k,it2) = 0.0 ; d_dkappa_it1(k,it2) = 0.0
1911  enddo ; enddo
1912 #endif
1913 
1914  if (do_newton) then ! K_Q needs to be calculated.
1915  do k=1,ks_kappa-1 ; k_q(k) = 0.0 ; enddo
1916  do k=ks_kappa,ke_kappa ; k_q(k) = kappa(k) / tke(k) ; enddo
1917  do k=ke_kappa+1,nz+1 ; k_q(k) = 0.0 ; enddo
1918  endif
1919 
1920  if (present(local_src)) then
1921  local_src(1) = 0.0 ; local_src(nz+1) = 0.0
1922  do k=2,nz
1923  diffusive_src = idz(k-1)*(kappa(k-1)-kappa(k)) + idz(k)*(kappa(k+1)-kappa(k))
1924  chg_by_k0 = kappa0 * ((idz(k-1)+idz(k)) / dz_int(k) + i_ld2(k))
1925  if (diffusive_src <= 0.0) then
1926  local_src(k) = k_src(k) + chg_by_k0
1927  else
1928  local_src(k) = (k_src(k) + chg_by_k0) + diffusive_src / dz_int(k)
1929  endif
1930  enddo
1931  endif
1932  if (present(kappa_src)) then
1933  kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
1934  do k=2,nz
1935  kappa_src(k) = k_src(k)
1936  enddo
1937  endif
1938 

Referenced by kappa_shear_column().

Here is the caller graph for this function:

◆ 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 without needing to duplicate the log entry.

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

Definition at line 2120 of file MOM_kappa_shear.F90.

2120  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
2121 ! Reads the parameter "USE_JACKSON_PARAM" and returns state.
2122  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
2123 
2124  logical :: do_kappa_shear
2125 
2126  call get_param(param_file, mdl, "USE_JACKSON_PARAM", do_kappa_shear, &
2127  default=.false., do_not_log=.true.)
2128  kappa_shear_at_vertex = .false.
2129  if (do_kappa_shear) &
2130  call get_param(param_file, mdl, "VERTEX_SHEAR", kappa_shear_at_vertex, &
2131  "If true, do the calculations of the shear-driven mixing "//&
2132  "at the cell vertices (i.e., the vorticity points).", &
2133  default=.false., do_not_log=.true.)
2134 

Referenced by mom_set_diffusivity::set_diffusivity_init(), mom_set_visc::set_visc_init(), and mom_set_visc::set_visc_register_restarts().

Here is the caller graph for this function:

◆ 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]usA dimensional unit scaling type
[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 [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.
[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 666 of file MOM_kappa_shear.F90.

666  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
667  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
668  real, dimension(SZK_(GV)+1), &
669  intent(inout) :: kappa !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
670  real, dimension(SZK_(GV)+1), &
671  intent(out) :: tke !< The Turbulent Kinetic Energy per unit mass at
672  !! an interface [Z2 T-2 ~> m2 s-2].
673  integer, intent(in) :: nzc !< The number of active layers in the column.
674  real, intent(in) :: f2 !< The square of the Coriolis parameter [T-2 ~> s-2].
675  real, intent(in) :: surface_pres !< The surface pressure [Pa].
676  real, dimension(SZK_(GV)), &
677  intent(in) :: dz !< The layer thickness [Z ~> m].
678  real, dimension(SZK_(GV)), &
679  intent(in) :: u0xdz !< The initial zonal velocity times dz [Z L T-1 ~> m2 s-1].
680  real, dimension(SZK_(GV)), &
681  intent(in) :: v0xdz !< The initial meridional velocity times dz [Z L T-1 ~> m2 s-1].
682  real, dimension(SZK_(GV)), &
683  intent(in) :: T0xdz !< The initial temperature times dz [degC Z ~> degC m].
684  real, dimension(SZK_(GV)), &
685  intent(in) :: S0xdz !< The initial salinity times dz [ppt Z ~> ppt m].
686  real, dimension(SZK_(GV)+1), &
687  intent(out) :: kappa_avg !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
688  real, dimension(SZK_(GV)+1), &
689  intent(out) :: tke_avg !< The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
690  real, intent(in) :: dt !< Time increment [T ~> s].
691  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
692  !! available thermodynamic fields. Absent fields
693  !! have NULL ptrs.
694  type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
695  !! call to kappa_shear_init.
696  real, dimension(SZK_(GV)+1), &
697  optional, intent(out) :: I_Ld2_1d !< The inverse of the squared mixing length [Z-2 ~> m-2].
698  real, dimension(SZK_(GV)+1), &
699  optional, intent(out) :: dz_Int_1d !< The extent of a finite-volume space surrounding an interface,
700  !! as used in calculating kappa and TKE [Z ~> m].
701 
702  real, dimension(nzc) :: &
703  u, & ! The zonal velocity after a timestep of mixing [L T-1 ~> m s-1].
704  v, & ! The meridional velocity after a timestep of mixing [L T-1 ~> m s-1].
705  Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
706  T, & ! The potential temperature after a timestep of mixing [degC].
707  Sal, & ! The salinity after a timestep of mixing [ppt].
708  u_test, v_test, & ! Temporary velocities [L T-1 ~> m s-1].
709  T_test, S_test ! Temporary temperatures [degC] and salinities [ppt].
710 
711  real, dimension(nzc+1) :: &
712  N2, & ! The squared buoyancy frequency at an interface [T-2 ~> s-2].
713  dz_Int, & ! The extent of a finite-volume space surrounding an interface,
714  ! as used in calculating kappa and TKE [Z ~> m].
715  i_dz_int, & ! The inverse of the distance between velocity & density points
716  ! above and below an interface [Z-1 ~> m-1]. This is used to
717  ! calculate N2, shear, and fluxes, and it might differ from
718  ! 1/dz_Int, as they have different uses.
719  s2, & ! The squared shear at an interface [T-2 ~> s-2].
720  a1, & ! a1 is the coupling between adjacent interfaces in the TKE,
721  ! velocity, and density equations [Z s-1 ~> m s-1] or [Z ~> m]
722  c1, & ! c1 is used in the tridiagonal (and similar) solvers.
723  k_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
724  kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
725  kappa_out, & ! The kappa that results from the kappa equation [Z2 T-1 ~> m2 s-1].
726  kappa_mid, & ! The average of the initial and predictor estimates of kappa [Z2 T-1 ~> m2 s-1].
727  tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2].
728  kappa_pred, & ! The value of kappa from a predictor step [Z2 T-1 ~> m2 s-1].
729  pressure, & ! The pressure at an interface [Pa].
730  t_int, & ! The temperature interpolated to an interface [degC].
731  sal_int, & ! The salinity interpolated to an interface [ppt].
732  dbuoy_dt, & ! The partial derivatives of buoyancy with changes in temperature
733  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].
734  i_l2_bdry, & ! The inverse of the square of twice the harmonic mean
735  ! distance to the top and bottom boundaries [Z-2 ~> m-2].
736  k_q, & ! Diffusivity divided by TKE [Z2 m-2 s2 T-1 ~> s].
737  k_q_tmp, & ! A temporary copy of diffusivity divided by TKE [Z2 m-2 s2 T-1 ~> s].
738  local_src_avg, & ! The time-integral of the local source [nondim].
739  tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1].
740  tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1].
741  tol_chg, & ! The tolerated kappa change integrated over a timestep [nondim].
742  dist_from_top, & ! The distance from the top surface [Z ~> m].
743  local_src ! The sum of all sources of kappa, including kappa_src and
744  ! sources from the elliptic term [T-1 ~> s-1].
745 
746  real :: dist_from_bot ! The distance from the bottom surface [Z ~> m].
747  real :: b1 ! The inverse of the pivot in the tridiagonal equations.
748  real :: bd1 ! A term in the denominator of b1.
749  real :: d1 ! 1 - c1 in the tridiagonal equations.
750  real :: gR0 ! A conversion factor from Z to Pa equal to Rho_0 times g
751  ! [Pa Z-1 = kg m-1 s-2 Z-1 ~> kg m-2 s-2].
752  real :: g_R0 ! g_R0 is a rescaled version of g/Rho [Z R-1 T-2 ~> m4 kg-1 s-2].
753  real :: Norm ! A factor that normalizes two weights to 1 [Z-2 ~> m-2].
754  real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration
755  ! relative to the local source [nondim].
756  real :: tol2 ! The tolerance for the change in the kappa source within an iteration
757  ! relative to the average local source over previous iterations [nondim].
758  real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc
759  ! within an iteration [nondim]. 0 < tol_dksrc_low < 1.
760  real :: Ri_crit ! The critical shear Richardson number for shear-
761  ! driven mixing [nondim]. The theoretical value is 0.25.
762  real :: dt_rem ! The remaining time to advance the solution [T ~> s].
763  real :: dt_now ! The time step used in the current iteration [T ~> s].
764  real :: dt_wt ! The fractional weight of the current iteration [nondim].
765  real :: dt_test ! A time-step that is being tested for whether it
766  ! gives acceptably small changes in k_src [T ~> s].
767  real :: Idtt ! Idtt = 1 / dt_test [T-1 ~> s-1].
768  real :: dt_inc ! An increment to dt_test that is being tested [T ~> s].
769 
770  real :: k0dt ! The background diffusivity times the timestep [Z2 ~> m2].
771  logical :: valid_dt ! If true, all levels so far exhibit acceptably small changes in k_src.
772  logical :: use_temperature ! If true, temperature and salinity have been
773  ! allocated and are being used as state variables.
774  integer :: ks_kappa, ke_kappa ! The k-range with nonzero kappas.
775  integer :: dt_halvings ! The number of times that the time-step is halved
776  ! in seeking an acceptable timestep. If none is
777  ! found, dt_rem*0.5^dt_halvings is used.
778  integer :: dt_refinements ! The number of 2-fold refinements that will be used
779  ! to estimate the maximum permitted time step. I.e.,
780  ! the resolution is 1/2^dt_refinements.
781  integer :: k, itt, itt_dt
782 #ifdef DEBUG
783  integer :: max_debug_itt ; parameter(max_debug_itt=20)
784  real :: wt(SZK_(GV)+1), wt_tot, I_wt_tot, wt_itt
785  real, dimension(SZK_(GV)+1) :: &
786  Ri_k, tke_prev, dtke, dkappa, dtke_norm, &
787  N2_debug, & ! A version of N2 for debugging [T-2 ~> s-2]
788  ksrc_av ! The average through the iterations of k_src [T-1 ~> s-1].
789  real, dimension(SZK_(GV)+1,0:max_debug_itt) :: &
790  tke_it1, N2_it1, Sh2_it1, ksrc_it1, kappa_it1, kprev_it1
791  real, dimension(SZK_(GV)+1,1:max_debug_itt) :: &
792  dkappa_it1, wt_it1, K_Q_it1, d_dkappa_it1, dkappa_norm
793  real, dimension(SZK_(GV),0:max_debug_itt) :: &
794  u_it1, v_it1, rho_it1, T_it1, S_it1
795  real, dimension(0:max_debug_itt) :: &
796  dk_wt_it1, dkpos_wt_it1, dkneg_wt_it1, k_mag
797  real, dimension(max_debug_itt) :: dt_it1
798 #endif
799 
800  ri_crit = cs%Rino_crit
801  gr0 = gv%z_to_H*gv%H_to_Pa
802  g_r0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
803  k0dt = dt*cs%kappa_0
804  !### These 3 tolerances are hard-coded and fixed for now. Perhaps these could be made dynamic later?
805  ! tol_dksrc = 0.5*tol_ksrc_chg ; tol_dksrc_low = 1.0 - 1.0/tol_ksrc_chg ?
806  tol_dksrc = 10.0 ; tol_dksrc_low = 0.95 ; tol2 = 2.0*cs%kappa_tol_err
807  dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
808  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
809 
810 
811  ! Set up Idz as the inverse of layer thicknesses.
812  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
813  ! Set up I_dz_int as the inverse of the distance between
814  ! adjacent layer centers.
815  i_dz_int(1) = 2.0 / dz(1)
816  dist_from_top(1) = 0.0
817  do k=2,nzc
818  i_dz_int(k) = 2.0 / (dz(k-1) + dz(k))
819  dist_from_top(k) = dist_from_top(k-1) + dz(k-1)
820  enddo
821  i_dz_int(nzc+1) = 2.0 / dz(nzc)
822 
823  ! Determine the velocities and thicknesses after eliminating massless
824  ! layers and applying a time-step of background diffusion.
825  if (nzc > 1) then
826  a1(2) = k0dt*i_dz_int(2)
827  b1 = 1.0 / (dz(1) + a1(2))
828  u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
829  t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
830  c1(2) = a1(2) * b1 ; d1 = dz(1) * b1 ! = 1 - c1
831  do k=2,nzc-1
832  bd1 = dz(k) + d1*a1(k)
833  a1(k+1) = k0dt*i_dz_int(k+1)
834  b1 = 1.0 / (bd1 + a1(k+1))
835  u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
836  v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
837  t(k) = b1 * (t0xdz(k) + a1(k)*t(k-1))
838  sal(k) = b1 * (s0xdz(k) + a1(k)*sal(k-1))
839  c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
840  enddo
841  ! rho or T and S have insulating boundary conditions, u & v use no-slip
842  ! bottom boundary conditions (if kappa0 > 0).
843  ! For no-slip bottom boundary conditions
844  b1 = 1.0 / ((dz(nzc) + d1*a1(nzc)) + k0dt*i_dz_int(nzc+1))
845  u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
846  v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
847  ! For insulating boundary conditions
848  b1 = 1.0 / (dz(nzc) + d1*a1(nzc))
849  t(nzc) = b1 * (t0xdz(nzc) + a1(nzc)*t(nzc-1))
850  sal(nzc) = b1 * (s0xdz(nzc) + a1(nzc)*sal(nzc-1))
851  do k=nzc-1,1,-1
852  u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
853  t(k) = t(k) + c1(k+1)*t(k+1) ; sal(k) = sal(k) + c1(k+1)*sal(k+1)
854  enddo
855  else
856  ! This is correct, but probably unnecessary.
857  b1 = 1.0 / (dz(1) + k0dt*i_dz_int(2))
858  u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
859  b1 = 1.0 / dz(1)
860  t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
861  endif
862 
863  ! This uses half the harmonic mean of thicknesses to provide two estimates
864  ! of the boundary between cells, and the inverse of the harmonic mean to
865  ! weight the two estimates. The net effect is that interfaces around thin
866  ! layers have thin cells, and the total thickness adds up properly.
867  ! The top- and bottom- interfaces have zero thickness, consistent with
868  ! adding additional zero thickness layers.
869  dz_int(1) = 0.0 ; dz_int(2) = dz(1)
870  do k=2,nzc-1
871  norm = 1.0 / (dz(k)*(dz(k-1)+dz(k+1)) + 2.0*dz(k-1)*dz(k+1))
872  dz_int(k) = dz_int(k) + dz(k) * ( ((dz(k)+dz(k+1)) * dz(k-1)) * norm)
873  dz_int(k+1) = dz(k) * ( ((dz(k-1)+dz(k)) * dz(k+1)) * norm)
874  enddo
875  dz_int(nzc) = dz_int(nzc) + dz(nzc) ; dz_int(nzc+1) = 0.0
876 
877  dist_from_bot = 0.0
878  do k=nzc,2,-1
879  dist_from_bot = dist_from_bot + dz(k)
880  i_l2_bdry(k) = (dist_from_top(k) + dist_from_bot)**2 / &
881  (dist_from_top(k) * dist_from_bot)**2
882  enddo
883 
884  ! Calculate thermodynamic coefficients and an initial estimate of N2.
885  if (use_temperature) then
886  pressure(1) = surface_pres
887  do k=2,nzc
888  pressure(k) = pressure(k-1) + gr0*dz(k-1)
889  t_int(k) = 0.5*(t(k-1) + t(k))
890  sal_int(k) = 0.5*(sal(k-1) + sal(k))
891  enddo
892  call calculate_density_derivs(t_int, sal_int, pressure, dbuoy_dt, &
893  dbuoy_ds, 2, nzc-1, tv%eqn_of_state, scale=-g_r0*us%kg_m3_to_R)
894  else
895  do k=1,nzc+1 ; dbuoy_dt(k) = -g_r0 ; dbuoy_ds(k) = 0.0 ; enddo
896  endif
897 
898 #ifdef DEBUG
899  n2_debug(1) = 0.0 ; n2_debug(nzc+1) = 0.0
900  do k=2,nzc
901  n2_debug(k) = max((dbuoy_dt(k) * (t0xdz(k-1)*idz(k-1) - t0xdz(k)*idz(k)) + &
902  dbuoy_ds(k) * (s0xdz(k-1)*idz(k-1) - s0xdz(k)*idz(k))) * &
903  i_dz_int(k), 0.0)
904  enddo
905  do k=1,nzc
906  u_it1(k,0) = u0xdz(k)*idz(k) ; v_it1(k,0) = v0xdz(k)*idz(k)
907  t_it1(k,0) = t0xdz(k)*idz(k) ; s_it1(k,0) = s0xdz(k)*idz(k)
908  enddo
909  do k=1,nzc+1
910  kprev_it1(k,0) = kappa(k) ; kappa_it1(k,0) = kappa(k)
911  tke_it1(k,0) = 0.0
912  n2_it1(k,0) = n2_debug(k) ; sh2_it1(k,0) = s2(k) ; ksrc_it1(k,0) = k_src(k)
913  enddo
914  do k=nzc+1,gv%ke
915  u_it1(k,0) = 0.0 ; v_it1(k,0) = 0.0
916  t_it1(k,0) = 0.0 ; s_it1(k,0) = 0.0
917  kprev_it1(k+1,0) = 0.0 ; kappa_it1(k+1,0) = 0.0 ; tke_it1(k+1,0) = 0.0
918  n2_it1(k+1,0) = 0.0 ; sh2_it1(k+1,0) = 0.0 ; ksrc_it1(k+1,0) = 0.0
919  enddo
920  do itt=1,max_debug_itt
921  dt_it1(itt) = 0.0
922  do k=1,gv%ke
923  u_it1(k,itt) = 0.0 ; v_it1(k,itt) = 0.0
924  t_it1(k,itt) = 0.0 ; s_it1(k,itt) = 0.0
925  rho_it1(k,itt) = 0.0
926  enddo
927  do k=1,gv%ke+1
928  kprev_it1(k,itt) = 0.0 ; kappa_it1(k,itt) = 0.0 ; tke_it1(k,itt) = 0.0
929  n2_it1(k,itt) = 0.0 ; sh2_it1(k,itt) = 0.0
930  ksrc_it1(k,itt) = 0.0
931  dkappa_it1(k,itt) = 0.0 ; wt_it1(k,itt) = 0.0
932  k_q_it1(k,itt) = 0.0 ; d_dkappa_it1(k,itt) = 0.0
933  enddo
934  enddo
935  do k=1,gv%ke+1 ; ksrc_av(k) = 0.0 ; enddo
936 #endif
937 
938  ! This call just calculates N2 and S2.
939  call calculate_projected_state(kappa, u, v, t, sal, 0.0, nzc, dz, i_dz_int, &
940  dbuoy_dt, dbuoy_ds, u, v, t, sal, gv, us, &
941  n2=n2, s2=s2, vel_underflow=cs%vel_underflow)
942 ! ----------------------------------------------------
943 ! Iterate
944 ! ----------------------------------------------------
945  dt_rem = dt
946  do k=1,nzc+1
947  k_q(k) = 0.0
948  kappa_avg(k) = 0.0 ; tke_avg(k) = 0.0
949  local_src_avg(k) = 0.0
950  ! Use the grid spacings to scale errors in the source.
951  if ( dz_int(k) > 0.0 ) &
952  local_src_avg(k) = 0.1 * k0dt * i_dz_int(k) / dz_int(k)
953  enddo
954 
955 ! call cpu_clock_end(id_clock_setup)
956 
957 ! do itt=1,CS%max_RiNo_it
958  do itt=1,cs%max_KS_it
959 
960 ! ----------------------------------------------------
961 ! Calculate new values of u, v, rho, N^2 and S.
962 ! ----------------------------------------------------
963 #ifdef DEBUG
964  do k=1,nzc+1
965  ri_k(k) = 1e3 ; if (s2(k) > 1e-3*n2(k)) ri_k(k) = n2(k) / s2(k)
966  if (itt > 1) then ; tke_prev(k) = tke(k) ; else ; tke_prev(k) = 0.0 ; endif
967  enddo
968 #endif
969 
970  ! call cpu_clock_begin(id_clock_KQ)
971  call find_kappa_tke(n2, s2, kappa, idz, dz_int, i_l2_bdry, f2, &
972  nzc, cs, gv, us, k_q, tke, kappa_out, kappa_src, local_src)
973  ! call cpu_clock_end(id_clock_KQ)
974 
975  ! call cpu_clock_begin(id_clock_avg)
976  ! Determine the range of non-zero values of kappa_out.
977  ks_kappa = gv%ke+1 ; ke_kappa = 0
978  do k=2,nzc ; if (kappa_out(k) > 0.0) then
979  ks_kappa = k ; exit
980  endif ; enddo
981  do k=nzc,ks_kappa,-1 ; if (kappa_out(k) > 0.0) then
982  ke_kappa = k ; exit
983  endif ; enddo
984  if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
985  ! call cpu_clock_end(id_clock_avg)
986 
987  ! Determine how long to use this value of kappa (dt_now).
988 
989  ! call cpu_clock_begin(id_clock_project)
990  if ((ke_kappa < ks_kappa) .or. (itt==cs%max_RiNo_it)) then
991  dt_now = dt_rem
992  else
993  ! Limit dt_now so that |k_src(k)-kappa_src(k)| < tol * local_src(k)
994  dt_test = dt_rem
995  do k=2,nzc
996  tol_max(k) = kappa_src(k) + tol_dksrc * local_src(k)
997  tol_min(k) = kappa_src(k) - tol_dksrc_low * local_src(k)
998  tol_chg(k) = tol2 * local_src_avg(k)
999  enddo
1000 
1001  do itt_dt=1,(cs%max_KS_it+1-itt)/2
1002  ! The maximum number of times that the time-step is halved in
1003  ! seeking an acceptable timestep is reduced with each iteration,
1004  ! so that as the maximum number of iterations is approached, the
1005  ! whole remaining timestep is used. Typically, an acceptable
1006  ! timestep is found long before the minimum is reached, so the
1007  ! value of max_KS_it may be unimportant, especially if it is large
1008  ! enough.
1009  call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*dt_test, nzc, dz, i_dz_int, &
1010  dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1011  gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1012  vel_underflow=cs%vel_underflow)
1013  valid_dt = .true.
1014  idtt = 1.0 / dt_test
1015  do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1016  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
1017  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1018  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1019  if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
1020  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
1021  valid_dt = .false. ; exit
1022  endif
1023  else
1024  if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
1025  valid_dt = .false. ; k_src(k) = 0.0 ; exit
1026  endif
1027  endif
1028  enddo
1029 
1030  if (valid_dt) exit
1031  dt_test = 0.5*dt_test
1032  enddo
1033  if ((dt_test < dt_rem) .and. valid_dt) then
1034  dt_inc = 0.5*dt_test
1035  do itt_dt=1,dt_refinements
1036  call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*(dt_test+dt_inc), &
1037  nzc, dz, i_dz_int, dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1038  gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, vel_underflow=cs%vel_underflow)
1039  valid_dt = .true.
1040  idtt = 1.0 / (dt_test+dt_inc)
1041  do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1042  if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
1043  k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1044  ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1045  if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
1046  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
1047  valid_dt = .false. ; exit
1048  endif
1049  else
1050  if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
1051  valid_dt = .false. ; k_src(k) = 0.0 ; exit
1052  endif
1053  endif
1054  enddo
1055 
1056  if (valid_dt) dt_test = dt_test + dt_inc
1057  dt_inc = 0.5*dt_inc
1058  enddo
1059  else
1060  dt_inc = 0.0
1061  endif
1062 
1063  dt_now = min(dt_test*(1.0+cs%kappa_tol_err)+dt_inc, dt_rem)
1064  do k=2,nzc
1065  local_src_avg(k) = local_src_avg(k) + dt_now * local_src(k)
1066  enddo
1067  endif ! Are all the values of kappa_out 0?
1068  ! call cpu_clock_end(id_clock_project)
1069 
1070  ! The state has already been projected forward. Now find new values of kappa.
1071 
1072  if (ke_kappa < ks_kappa) then
1073  ! There is no mixing now, and will not be again.
1074  ! call cpu_clock_begin(id_clock_avg)
1075  dt_wt = dt_rem / dt ; dt_rem = 0.0
1076  do k=1,nzc+1
1077  kappa_mid(k) = 0.0
1078  ! This would be here but does nothing.
1079  ! kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
1080  tke_avg(k) = tke_avg(k) + dt_wt*tke(k)
1081 #ifdef DEBUG
1082  tke_pred(k) = tke(k) ; kappa_pred(k) = 0.0 ; kappa(k) = 0.0
1083 #endif
1084  enddo
1085  ! call cpu_clock_end(id_clock_avg)
1086  else
1087  ! call cpu_clock_begin(id_clock_project)
1088  call calculate_projected_state(kappa_out, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
1089  dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1090  gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1091  vel_underflow=cs%vel_underflow)
1092  ! call cpu_clock_end(id_clock_project)
1093 
1094  ! call cpu_clock_begin(id_clock_KQ)
1095  do k=1,nzc+1 ; k_q_tmp(k) = k_q(k) ; enddo
1096  call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1097  nzc, cs, gv, us, k_q_tmp, tke_pred, kappa_pred)
1098  ! call cpu_clock_end(id_clock_KQ)
1099 
1100  ks_kappa = gv%ke+1 ; ke_kappa = 0
1101  do k=1,nzc+1
1102  kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1103  if ((kappa_mid(k) > 0.0) .and. (k<ks_kappa)) ks_kappa = k
1104  if (kappa_mid(k) > 0.0) ke_kappa = k
1105  enddo
1106 
1107  ! call cpu_clock_begin(id_clock_project)
1108  call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
1109  dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1110  gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1111  vel_underflow=cs%vel_underflow)
1112  ! call cpu_clock_end(id_clock_project)
1113 
1114  ! call cpu_clock_begin(id_clock_KQ)
1115  call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1116  nzc, cs, gv, us, k_q, tke_pred, kappa_pred)
1117  ! call cpu_clock_end(id_clock_KQ)
1118 
1119  ! call cpu_clock_begin(id_clock_avg)
1120  dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
1121  do k=1,nzc+1
1122  kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1123  kappa_avg(k) = kappa_avg(k) + kappa_mid(k)*dt_wt
1124  tke_avg(k) = tke_avg(k) + dt_wt*0.5*(tke_pred(k) + tke(k))
1125  kappa(k) = kappa_pred(k) ! First guess for the next iteration.
1126  enddo
1127  ! call cpu_clock_end(id_clock_avg)
1128  endif
1129 
1130  if (dt_rem > 0.0) then
1131  ! Update the values of u, v, T, Sal, N2, and S2 for the next iteration.
1132  ! call cpu_clock_begin(id_clock_project)
1133  call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, &
1134  dz, i_dz_int, dbuoy_dt, dbuoy_ds, u, v, t, sal, &
1135  gv, us, n2, s2, vel_underflow=cs%vel_underflow)
1136  ! call cpu_clock_end(id_clock_project)
1137  endif
1138 
1139 #ifdef DEBUG
1140  if (itt <= max_debug_itt) then
1141  dt_it1(itt) = dt_now
1142  dk_wt_it1(itt) = 0.0 ; dkpos_wt_it1(itt) = 0.0 ; dkneg_wt_it1(itt) = 0.0
1143  k_mag(itt) = 0.0
1144  wt_itt = 1.0/real(itt) ; wt_tot = 0.0
1145  do k=1,nzc+1
1146  ksrc_av(k) = (1.0-wt_itt)*ksrc_av(k) + wt_itt*k_src(k)
1147  wt_tot = wt_tot + dz_int(k) * ksrc_av(k)
1148  enddo
1149  ! Use the 1/0=0 convention.
1150  i_wt_tot = 0.0 ; if (wt_tot > 0.0) i_wt_tot = 1.0/wt_tot
1151 
1152  do k=1,nzc+1
1153  wt(k) = (dz_int(k)*ksrc_av(k)) * i_wt_tot
1154  k_mag(itt) = k_mag(itt) + wt(k)*kappa_mid(k)
1155  dkappa_it1(k,itt) = kappa_pred(k) - kappa_out(k)
1156  dk_wt_it1(itt) = dk_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1157  if (dkappa_it1(k,itt) > 0.0) then
1158  dkpos_wt_it1(itt) = dkpos_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1159  else
1160  dkneg_wt_it1(itt) = dkneg_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1161  endif
1162  wt_it1(k,itt) = wt(k)
1163  enddo
1164  endif
1165  do k=1,nzc+1
1166  ri_k(k) = 1e3 ; if (n2(k) < 1e3 * s2(k)) ri_k(k) = n2(k) / s2(k)
1167  dtke(k) = tke_pred(k) - tke(k)
1168  dtke_norm(k) = dtke(k) / (0.5*(tke(k) + tke_pred(k)))
1169  dkappa(k) = kappa_pred(k) - kappa_out(k)
1170  enddo
1171  if (itt <= max_debug_itt) then
1172  do k=1,nzc
1173  u_it1(k,itt) = u(k) ; v_it1(k,itt) = v(k)
1174  t_it1(k,itt) = t(k) ; s_it1(k,itt) = sal(k)
1175  enddo
1176  do k=1,nzc+1
1177  kprev_it1(k,itt) = kappa_out(k)
1178  kappa_it1(k,itt) = kappa_mid(k) ; tke_it1(k,itt) = 0.5*(tke(k)+tke_pred(k))
1179  n2_it1(k,itt)=n2(k) ; sh2_it1(k,itt)=s2(k)
1180  ksrc_it1(k,itt) = kappa_src(k)
1181  k_q_it1(k,itt) = kappa_out(k) / (tke(k))
1182  if (itt > 1) then
1183  if (abs(dkappa_it1(k,itt-1)) > 1e-20) &
1184  d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
1185  endif
1186  dkappa_norm(k,itt) = dkappa(k) / max(0.5*(kappa_pred(k) + kappa_out(k)), us%m2_s_to_Z2_T*1e-100)
1187  enddo
1188  endif
1189 #endif
1190 
1191  if (dt_rem <= 0.0) exit
1192 
1193  enddo ! end itt loop
1194 
1195 #ifdef ADD_DIAGNOSTICS
1196  if (present(i_ld2_1d)) then
1197  do k=1,gv%ke+1 ; i_ld2_1d(k) = 0.0 ; enddo
1198  do k=2,nzc ; if (tke(k) > 0.0) &
1199  i_ld2_1d(k) = i_l2_bdry(k) + (n2(k) / cs%lambda**2 + f2) / tke(k)
1200  enddo
1201  endif
1202  if (present(dz_int_1d)) then
1203  do k=1,nzc+1 ; dz_int_1d(k) = dz_int(k) ; enddo
1204  do k=nzc+2,gv%ke ; dz_int_1d(k) = 0.0 ; enddo
1205  endif
1206 #endif
1207 

References calculate_projected_state(), and find_kappa_tke().

Referenced by calc_kappa_shear_vertex(), and calculate_kappa_shear().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 subroutineinitializesthe 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 1943 of file MOM_kappa_shear.F90.

1943  type(time_type), intent(in) :: Time !< The current model time.
1944  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1945  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1946  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1947  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1948  !! parameters.
1949  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1950  !! output.
1951  type(Kappa_shear_CS), pointer :: CS !< A pointer that is set to point to the control
1952  !! structure for this module
1953  logical :: kappa_shear_init !< True if module is to be used, False otherwise
1954 
1955  ! Local variables
1956  logical :: merge_mixedlayer
1957 ! This include declares and sets the variable "version".
1958 #include "version_variable.h"
1959  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
1960  real :: kappa_0_unscaled ! The value of kappa_0 in MKS units [m2 s-1]
1961  real :: KD_normal ! The KD of the main model, read here only as a parameter
1962  ! for setting the default of KD_SMOOTH in MKS units [m2 s-1]
1963 
1964  if (associated(cs)) then
1965  call mom_error(warning, "kappa_shear_init called with an associated "// &
1966  "control structure.")
1967  return
1968  endif
1969  allocate(cs)
1970 
1971  ! The Jackson-Hallberg-Legg shear mixing parameterization uses the following
1972  ! 6 nondimensional coefficients. That paper gives 3 best fit parameter sets.
1973  ! Ri_Crit Rate FRi_Curv K_buoy TKE_N TKE_Shear
1974  ! p1: 0.25 0.089 -0.97 0.82 0.24 0.14
1975  ! p2: 0.30 0.085 -0.94 0.86 0.26 0.13
1976  ! p3: 0.35 0.088 -0.89 0.81 0.28 0.12
1977  ! Future research will reveal how these should be modified to take
1978  ! subgridscale inhomogeneity into account.
1979 
1980 ! Set default, read and log parameters
1981  call log_version(param_file, mdl, version, &
1982  "Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008")
1983  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, &
1984  "If true, use the Jackson-Hallberg-Legg (JPO 2008) "//&
1985  "shear mixing parameterization.", default=.false.)
1986  call get_param(param_file, mdl, "VERTEX_SHEAR", cs%KS_at_vertex, &
1987  "If true, do the calculations of the shear-driven mixing "//&
1988  "at the cell vertices (i.e., the vorticity points).", &
1989  default=.false.)
1990  call get_param(param_file, mdl, "RINO_CRIT", cs%RiNo_crit, &
1991  "The critical Richardson number for shear mixing.", &
1992  units="nondim", default=0.25)
1993  call get_param(param_file, mdl, "SHEARMIX_RATE", cs%Shearmix_rate, &
1994  "A nondimensional rate scale for shear-driven entrainment. "//&
1995  "Jackson et al find values in the range of 0.085-0.089.", &
1996  units="nondim", default=0.089)
1997  call get_param(param_file, mdl, "MAX_RINO_IT", cs%max_RiNo_it, &
1998  "The maximum number of iterations that may be used to "//&
1999  "estimate the Richardson number driven mixing.", &
2000  units="nondim", default=50)
2001  call get_param(param_file, mdl, "KD", kd_normal, default=1.0e-7, do_not_log=.true.)
2002  call get_param(param_file, mdl, "KD_KAPPA_SHEAR_0", cs%kappa_0, &
2003  "The background diffusivity that is used to smooth the "//&
2004  "density and shear profiles before solving for the "//&
2005  "diffusivities. Defaults to value of KD.", &
2006  units="m2 s-1", default=kd_normal, scale=us%m2_s_to_Z2_T, unscaled=kappa_0_unscaled)
2007  call get_param(param_file, mdl, "KD_TRUNC_KAPPA_SHEAR", cs%kappa_trunc, &
2008  "The value of shear-driven diffusivity that is considered negligible "//&
2009  "and is rounded down to 0. The default is 1% of KD_KAPPA_SHEAR_0.", &
2010  units="m2 s-1", default=0.01*kappa_0_unscaled, scale=us%m2_s_to_Z2_T)
2011  call get_param(param_file, mdl, "FRI_CURVATURE", cs%FRi_curvature, &
2012  "The nondimensional curvature of the function of the "//&
2013  "Richardson number in the kappa source term in the "//&
2014  "Jackson et al. scheme.", units="nondim", default=-0.97)
2015  call get_param(param_file, mdl, "TKE_N_DECAY_CONST", cs%C_N, &
2016  "The coefficient for the decay of TKE due to "//&
2017  "stratification (i.e. proportional to N*tke). "//&
2018  "The values found by Jackson et al. are 0.24-0.28.", &
2019  units="nondim", default=0.24)
2020 ! call get_param(param_file, mdl, "LAYER_KAPPA_STAGGER", CS%layer_stagger, &
2021 ! default=.false.)
2022  call get_param(param_file, mdl, "TKE_SHEAR_DECAY_CONST", cs%C_S, &
2023  "The coefficient for the decay of TKE due to shear (i.e. "//&
2024  "proportional to |S|*tke). The values found by Jackson "//&
2025  "et al. are 0.14-0.12.", units="nondim", default=0.14)
2026  call get_param(param_file, mdl, "KAPPA_BUOY_SCALE_COEF", cs%lambda, &
2027  "The coefficient for the buoyancy length scale in the "//&
2028  "kappa equation. The values found by Jackson et al. are "//&
2029  "in the range of 0.81-0.86.", units="nondim", default=0.82)
2030  call get_param(param_file, mdl, "KAPPA_N_OVER_S_SCALE_COEF2", cs%lambda2_N_S, &
2031  "The square of the ratio of the coefficients of the "//&
2032  "buoyancy and shear scales in the diffusivity equation, "//&
2033  "Set this to 0 (the default) to eliminate the shear scale. "//&
2034  "This is only used if USE_JACKSON_PARAM is true.", &
2035  units="nondim", default=0.0)
2036  call get_param(param_file, mdl, "KAPPA_SHEAR_TOL_ERR", cs%kappa_tol_err, &
2037  "The fractional error in kappa that is tolerated. "//&
2038  "Iteration stops when changes between subsequent "//&
2039  "iterations are smaller than this everywhere in a "//&
2040  "column. The peak diffusivities usually converge most "//&
2041  "rapidly, and have much smaller errors than this.", &
2042  units="nondim", default=0.1)
2043  call get_param(param_file, mdl, "TKE_BACKGROUND", cs%TKE_bg, &
2044  "A background level of TKE used in the first iteration "//&
2045  "of the kappa equation. TKE_BACKGROUND could be 0.", &
2046  units="m2 s-2", default=0.0, scale=us%m_to_Z**2*us%T_to_s**2)
2047  call get_param(param_file, mdl, "KAPPA_SHEAR_ELIM_MASSLESS", cs%eliminate_massless, &
2048  "If true, massless layers are merged with neighboring "//&
2049  "massive layers in this calculation. The default is "//&
2050  "true and I can think of no good reason why it should "//&
2051  "be false. This is only used if USE_JACKSON_PARAM is true.", &
2052  default=.true.)
2053  call get_param(param_file, mdl, "MAX_KAPPA_SHEAR_IT", cs%max_KS_it, &
2054  "The maximum number of iterations that may be used to "//&
2055  "estimate the time-averaged diffusivity.", units="nondim", &
2056  default=13)
2057  call get_param(param_file, mdl, "PRANDTL_TURB", cs%Prandtl_turb, &
2058  "The turbulent Prandtl number applied to shear "//&
2059  "instability.", units="nondim", default=1.0, do_not_log=.true.)
2060  call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &
2061  "A negligibly small velocity magnitude below which velocity "//&
2062  "components are set to 0. A reasonable value might be "//&
2063  "1e-30 m/s, which is less than an Angstrom divided by "//&
2064  "the age of the universe.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
2065 
2066  call get_param(param_file, mdl, "DEBUG_KAPPA_SHEAR", cs%debug, &
2067  "If true, write debugging data for the kappa-shear code. \n"//&
2068  "Caution: this option is _very_ verbose and should only "//&
2069  "be used in single-column mode!", &
2070  default=.false., debuggingparam=.true.)
2071  call get_param(param_file, mdl, "KAPPA_SHEAR_ITER_BUG", cs%dKdQ_iteration_bug, &
2072  "If true. use an older, dimensionally inconsistent estimate of the "//&
2073  "derivative of diffusivity with energy in the Newton's method iteration. "//&
2074  "The bug causes undercorrections when dz > 1m.", default=.true.)
2075 ! id_clock_KQ = cpu_clock_id('Ocean KS kappa_shear', grain=CLOCK_ROUTINE)
2076 ! id_clock_avg = cpu_clock_id('Ocean KS avg', grain=CLOCK_ROUTINE)
2077 ! id_clock_project = cpu_clock_id('Ocean KS project', grain=CLOCK_ROUTINE)
2078 ! id_clock_setup = cpu_clock_id('Ocean KS setup', grain=CLOCK_ROUTINE)
2079 
2080  cs%nkml = 1
2081  if (gv%nkml>0) then
2082  call get_param(param_file, mdl, "KAPPA_SHEAR_MERGE_ML",merge_mixedlayer, &
2083  "If true, combine the mixed layers together before "//&
2084  "solving the kappa-shear equations.", default=.true.)
2085  if (merge_mixedlayer) cs%nkml = gv%nkml
2086  endif
2087 
2088 ! Forego remainder of initialization if not using this scheme
2089  if (.not. kappa_shear_init) return
2090 
2091  cs%diag => diag
2092 
2093  cs%id_Kd_shear = register_diag_field('ocean_model','Kd_shear',diag%axesTi,time, &
2094  'Shear-driven Diapycnal Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2095  cs%id_TKE = register_diag_field('ocean_model','TKE_shear',diag%axesTi,time, &
2096  'Shear-driven Turbulent Kinetic Energy', 'm2 s-2', conversion=us%Z_to_m**2*us%s_to_T**2)
2097 #ifdef ADD_DIAGNOSTICS
2098  cs%id_ILd2 = register_diag_field('ocean_model','ILd2_shear',diag%axesTi,time, &
2099  'Inverse kappa decay scale at interfaces', 'm-2', conversion=us%m_to_Z**2)
2100  cs%id_dz_Int = register_diag_field('ocean_model','dz_Int_shear',diag%axesTi,time, &
2101  'Finite volume thickness of interfaces', 'm', conversion=us%Z_to_m)
2102 #endif
2103 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ 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 2109 of file MOM_kappa_shear.F90.

2109  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
2110 ! Reads the parameter "USE_JACKSON_PARAM" and returns state.
2111  character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
2112 
2113  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_is_used, &
2114  default=.false., do_not_log=.true.)

Referenced by mom_cvmix_shear::cvmix_shear_init(), mom_diabatic_driver::diabatic_driver_init(), mom_set_visc::set_visc_init(), and mom_set_visc::set_visc_register_restarts().

Here is the caller graph for this function: