MOM6
mom_meke Module Reference

Detailed Description

Implements the Mesoscale Eddy Kinetic Energy framework.

The Mesoscale Eddy Kinetic Energy (MEKE) framework

The MEKE framework accounts for the mean potential energy removed by the first order closures used to parameterize mesoscale eddies. It requires closure at the second order, namely dissipation and transport of eddy energy.

Monitoring the sub-grid scale eddy energy budget provides a means to predict a sub-grid eddy-velocity scale which can be used in the lower order closures.

MEKE equations

The eddy kinetic energy equation is:

\[ \partial_\tilde{t} E = \overbrace{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v }^\text{sources} - \overbrace{ ( \lambda + C_d | U_d | \gamma_b^2 ) E }^\text{local dissipation} + \overbrace{ \nabla \cdot ( ( \kappa_E + \gamma_M \kappa_M ) \nabla E - \kappa_4 \nabla^3 E ) }^\text{smoothing} \]

where \( E \) is the eddy kinetic energy (variable MEKE) with units of m2s-2, and \(\tilde{t} = a t\) is a scaled time. The non-dimensional factor \( a\geq 1 \) is used to accelerate towards equilibrium.

The MEKE equation is two-dimensional and obtained by depth averaging the the three-dimensional eddy energy equation. In the following expressions \( \left< \phi \right> = \frac{1}{H} \int^\eta_{-D} \phi \, dz \) maps three dimensional terms into the two-dimensional quantities needed.

MEKE source terms

The source term \( \dot{E}_b \) is a constant background source of energy intended to avoid the limit \(E\rightarrow 0\).

The "GM" source term

\[ \dot{E}_\eta = - \left< \overline{w^\prime b^\prime} \right> = \left< \kappa_h N^2S^2 \right> \approx \left< \kappa_h g\prime |\nabla_\sigma \eta|^2 \right>\]

equals the mean potential energy removed by the Gent-McWilliams closure, and is excluded/included in the MEKE budget by the efficiency parameter \( \gamma_\eta \in [0,1] \).

The "frictional" source term

\[ \dot{E}_{v} = \left< u \cdot \tau_h \right> \]

equals the mean kinetic energy removed by lateral viscous fluxes, and is excluded/included in the MEKE budget by the efficiency parameter \( \gamma_v \in [0,1] \).

MEKE dissipation terms

The local dissipation of \( E \) is parameterized through a linear damping, \(\lambda\), and bottom drag, \( C_d | U_d | \gamma_b^2 \). The \( \gamma_b \) accounts for the weak projection of the column-mean eddy velocty to the bottom. In other words, the bottom velocity is estimated as \( \gamma_b U_e \). The bottom drag coefficient, \( C_d \) is the same as that used in the bottom friction in the mean model equations.

The bottom drag velocity scale, \( U_d \), has contributions from the resolved state and \( E \):

\[ U_d = \sqrt{ U_b^2 + |u|^2_{z=-D} + |\gamma_b U_e|^2 } .\]

where the eddy velocity scale, \( U_e \), is given by:

\[ U_e = \sqrt{ 2 E } .\]

\( U_b \) is a constant background bottom velocity scale and is typically not used (i.e. set to zero).

Following Jansen et al., 2015, the projection of eddy energy on to the bottom is given by the ratio of bottom energy to column mean energy:

\[ \gamma_b^2 = \frac{E_b}{E} = \gamma_{d0} + \left( 1 + c_{b} \frac{L_d}{L_f} \right)^{-\frac{4}{5}} , \]

\[ \gamma_b^2 \leftarrow \max{\left( \gamma_b^2, \gamma_{min}^2 \right)} . \]

MEKE smoothing terms

\( E \) is laterally diffused by a diffusivity \( \kappa_E + \gamma_M \kappa_M \) where \( \kappa_E \) is a constant diffusivity and the term \( \gamma_M \kappa_M \) is a "self diffusion" using the diffusivity calculated in the section Diffusivity derived from MEKE. \( \kappa_4 \) is a constant bi-harmonic diffusivity.

Diffusivity derived from MEKE

The predicted eddy velocity scale, \( U_e \), can be combined with a mixing length scale to form a diffusivity. The primary use of a MEKE derived diffusivity is for use in thickness diffusion (module mom_thickness_diffuse) and optionally in along isopycnal mixing of tracers (module mom_tracer_hor_diff). The original form used (enabled with MEKE_OLD_LSCALE=True):

\[ \kappa_M = \gamma_\kappa \sqrt{ \gamma_t^2 U_e^2 A_\Delta } \]

where \( A_\Delta \) is the area of the grid cell. Following Jansen et al., 2015, we now use

\[ \kappa_M = \gamma_\kappa l_M \sqrt{ \gamma_t^2 U_e^2 } \]

where \( \gamma_\kappa \in [0,1] \) is a non-dimensional factor and, following Jansen et al., 2015, \(\gamma_t^2\) is the ratio of barotropic eddy energy to column mean eddy energy given by

\[ \gamma_t^2 = \frac{E_t}{E} = \left( 1 + c_{t} \frac{L_d}{L_f} \right)^{-\frac{1}{4}} , \]

\[ \gamma_t^2 \leftarrow \max{\left( \gamma_t^2, \gamma_{min}^2 \right)} . \]

The length-scale is a configurable combination of multiple length scales:

\[ l_M = \left( \frac{\alpha_d}{L_d} + \frac{\alpha_f}{L_f} + \frac{\alpha_R}{L_R} + \frac{\alpha_e}{L_e} + \frac{\alpha_\Delta}{L_\Delta} + \frac{\delta[L_c]}{L_c} \right)^{-1} \]

where

\begin{eqnarray*} L_d & = & \sqrt{\frac{c_g^2}{f^2+2\beta c_g}} \sim \frac{ c_g }{f} \\\ L_R & = & \sqrt{\frac{U_e}{\beta}} \\\ L_e & = & \frac{U_e}{|S| N} \\\ L_f & = & \frac{H}{c_d} \\\ L_\Delta & = & \sqrt{A_\Delta} . \end{eqnarray*}

\(L_c\) is a constant and \(\delta[L_c]\) is the impulse function so that the term \(\frac{\delta[L_c]}{L_c}\) evaluates to \(\frac{1}{L_c}\) when \(L_c\) is non-zero but is dropped if \(L_c=0\).

Viscosity derived from MEKE

As for \( \kappa_M \), the predicted eddy velocity scale can be used to form an eddy viscosity:

\[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } . \]

Limit cases for local source-dissipative balance

Note that in steady-state (or when \( a>>1 \)) and there is no diffusion of \( E \) then

\[ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v }{ \lambda + C_d|U_d|\gamma_b^2 } . \]

In the linear drag limit, where \( U_e << \min(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \), the equilibrium becomes \( \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v }{ \lambda + C_d \sqrt{ U_b^2 + |u|^2_{z=-D} } } \).

In the nonlinear drag limit, where \( U_e >> \max(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \), the equilibrium becomes \( \overline{E} \approx \left( \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v }{ \sqrt{2} C_d \gamma_b^3 } \right)^\frac{2}{3} \).

MEKE module parameters

Symbol Module parameter
- USE_MEKE
\( a \) MEKE_DTSCALE
\( \dot{E}_b \) MEKE_BGSRC
\( \gamma_\eta \) MEKE_GMCOEFF
\( \gamma_v \) MEKE_FrCOEFF
\( \lambda \) MEKE_DAMPING
\( U_b \) MEKE_USCALE
\( \gamma_{d0} \) MEKE_CD_SCALE
\( c_{b} \) MEKE_CB
\( c_{t} \) MEKE_CT
\( \kappa_E \) MEKE_KH
\( \kappa_4 \) MEKE_K4
\( \gamma_\kappa \) MEKE_KHCOEFF
\( \gamma_M \) MEKE_KHMEKE_FAC
\( \gamma_u \) MEKE_VISCOSITY_COEFF
\( \gamma_{min}^2 \)MEKE_MIN_GAMMA2
\( \alpha_d \) MEKE_ALPHA_DEFORM
\( \alpha_f \) MEKE_ALPHA_FRICT
\( \alpha_R \) MEKE_ALPHA_RHINES
\( \alpha_e \) MEKE_ALPHA_EADY
\( \alpha_\Delta \) MEKE_ALPHA_GRID
\( L_c \) MEKE_FIXED_MIXING_LENGTH
- MEKE_KHTH_FAC
- MEKE_KHTR_FAC
Symbol Model parameter
\( C_d \) CDRAG

References

Jansen, M. F., A. J. Adcroft, R. Hallberg, and I. M. Held, 2015: Parameterization of eddy fluxes based on a mesoscale energy budget. Ocean Modelling, 92, 28–41, http://doi.org/10.1016/j.ocemod.2015.05.007 .

Marshall, D. P., and A. J. Adcroft, 2010: Parameterization of ocean eddies: Potential vorticity mixing, energetics and Arnold first stability theorem. Ocean Modelling, 32, 188–204, http://doi.org/10.1016/j.ocemod.2010.02.001 .

Data Types

type  meke_cs
 Control structure that contains MEKE parameters and diagnostics handles. More...
 

Functions/Subroutines

subroutine, public step_forward_meke (MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
 Integrates forward-in-time the MEKE eddy energy equation. See MEKE equations. More...
 
subroutine meke_equilibrium (CS, MEKE, G, GV, SN_u, SN_v, drag_rate_visc, I_mass)
 Calculates the equilibrium solutino where the source depends only on MEKE diffusivity and there is no lateral diffusion of MEKE. Results is in MEKEMEKE. More...
 
subroutine meke_lengthscales (CS, MEKE, G, SN_u, SN_v, EKE, bottomFac2, barotrFac2, LmixScale)
 Calculates the eddy mixing length scale and \(\gamma_b\) and \(\gamma_t\) functions that are ratios of either bottom or barotropic eddy energy to the column eddy energy, respectively. See MEKE equations. More...
 
subroutine meke_lengthscales_0d (CS, area, beta, depth, Rd_dx, SN, EKE, bottomFac2, barotrFac2, LmixScale, Lrhines, Leady)
 Calculates the eddy mixing length scale and \(\gamma_b\) and \(\gamma_t\) functions that are ratios of either bottom or barotropic eddy energy to the column eddy energy, respectively. See MEKE equations. More...
 
logical function, public meke_init (Time, G, param_file, diag, CS, MEKE, restart_CS)
 Initializes the MOM_MEKE module and reads parameters. Returns True if module is to be used, otherwise returns False. More...
 
subroutine, public meke_alloc_register_restart (HI, param_file, MEKE, restart_CS)
 Allocates memory and register restart fields for the MOM_MEKE module. More...
 
subroutine, public meke_end (MEKE, CS)
 Deallocates any variables allocated in MEKE_init or MEKE_alloc_register_restart. More...
 

Function/Subroutine Documentation

◆ meke_alloc_register_restart()

subroutine, public mom_meke::meke_alloc_register_restart ( type(hor_index_type), intent(in)  HI,
type(param_file_type), intent(in)  param_file,
type(meke_type), pointer  MEKE,
type(mom_restart_cs), pointer  restart_CS 
)

Allocates memory and register restart fields for the MOM_MEKE module.

Parameters
[in]hiHorizontal index structure
[in]param_fileParameter file parser structure.
mekeA structure with MEKE-related fields.
restart_csRestart control structure for MOM_MEKE.

Definition at line 1043 of file MOM_MEKE.F90.

References mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), and mom_io::var_desc().

Referenced by mom::initialize_mom().

1043 ! Arguments
1044  type(hor_index_type), intent(in) :: hi !< Horizontal index structure
1045  type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
1046  type(meke_type), pointer :: meke !< A structure with MEKE-related fields.
1047  type(mom_restart_cs), pointer :: restart_cs !< Restart control structure for MOM_MEKE.
1048 ! Local variables
1049  type(vardesc) :: vd
1050  real :: meke_gmcoeff, meke_frcoeff, meke_khcoeff, meke_visccoeff
1051  logical :: usemeke
1052  integer :: isd, ied, jsd, jed
1053 
1054 ! Determine whether this module will be used
1055  usemeke = .false.; call read_param(param_file,"USE_MEKE",usemeke)
1056 
1057 ! Read these parameters to determine what should be in the restarts
1058  meke_gmcoeff =-1.; call read_param(param_file,"MEKE_GMCOEFF",meke_gmcoeff)
1059  meke_frcoeff =-1.; call read_param(param_file,"MEKE_FRCOEFF",meke_frcoeff)
1060  meke_khcoeff =1.; call read_param(param_file,"MEKE_KHCOEFF",meke_khcoeff)
1061  meke_visccoeff =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF",meke_visccoeff)
1062 
1063 ! Allocate control structure
1064  if (associated(meke)) then
1065  call mom_error(warning, "MEKE_alloc_register_restart called with an associated "// &
1066  "MEKE type.")
1067  return
1068  else; allocate(meke); endif
1069 
1070  if (.not. usemeke) return
1071 
1072 ! Allocate memory
1073  call mom_mesg("MEKE_alloc_register_restart: allocating and registering", 5)
1074  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
1075  allocate(meke%MEKE(isd:ied,jsd:jed)) ; meke%MEKE(:,:) = 0.0
1076  vd = var_desc("MEKE", "m2 s-2", hor_grid='h', z_grid='1', &
1077  longname="Mesoscale Eddy Kinetic Energy")
1078  call register_restart_field(meke%MEKE, vd, .false., restart_cs)
1079  if (meke_gmcoeff>=0.) then
1080  allocate(meke%GM_src(isd:ied,jsd:jed)) ; meke%GM_src(:,:) = 0.0
1081  endif
1082  if (meke_frcoeff>=0.) then
1083  allocate(meke%mom_src(isd:ied,jsd:jed)) ; meke%mom_src(:,:) = 0.0
1084  endif
1085  if (meke_khcoeff>=0.) then
1086  allocate(meke%Kh(isd:ied,jsd:jed)) ; meke%Kh(:,:) = 0.0
1087  vd = var_desc("MEKE_Kh", "m2 s-1",hor_grid='h',z_grid='1', &
1088  longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy")
1089  call register_restart_field(meke%Kh, vd, .false., restart_cs)
1090  endif
1091  allocate(meke%Rd_dx_h(isd:ied,jsd:jed)) ; meke%Rd_dx_h(:,:) = 0.0
1092  if (meke_visccoeff/=0.) then
1093  allocate(meke%Ku(isd:ied,jsd:jed)) ; meke%Ku(:,:) = 0.0
1094  vd = var_desc("MEKE_Ah", "m2 s-1", hor_grid='h', z_grid='1', &
1095  longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy")
1096  call register_restart_field(meke%Ku, vd, .false., restart_cs)
1097  endif
1098 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ meke_end()

subroutine, public mom_meke::meke_end ( type(meke_type), pointer  MEKE,
type(meke_cs), pointer  CS 
)

Deallocates any variables allocated in MEKE_init or MEKE_alloc_register_restart.

Parameters
mekeA structure with MEKE-related fields.
csThe control structure for MOM_MEKE.

Definition at line 1104 of file MOM_MEKE.F90.

1104  type(meke_type), pointer :: meke !< A structure with MEKE-related fields.
1105  type(meke_cs), pointer :: cs !< The control structure for MOM_MEKE.
1106 
1107  if (associated(cs)) deallocate(cs)
1108 
1109  if (.not.associated(meke)) return
1110 
1111  if (associated(meke%MEKE)) deallocate(meke%MEKE)
1112  if (associated(meke%GM_src)) deallocate(meke%GM_src)
1113  if (associated(meke%mom_src)) deallocate(meke%mom_src)
1114  if (associated(meke%Kh)) deallocate(meke%Kh)
1115  if (associated(meke%Ku)) deallocate(meke%Ku)
1116  if (allocated(cs%del2MEKE)) deallocate(cs%del2MEKE)
1117  deallocate(meke)
1118 

◆ meke_equilibrium()

subroutine mom_meke::meke_equilibrium ( type(meke_cs), pointer  CS,
type(meke_type), pointer  MEKE,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in)  SN_u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in)  SN_v,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  drag_rate_visc,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  I_mass 
)
private

Calculates the equilibrium solutino where the source depends only on MEKE diffusivity and there is no lateral diffusion of MEKE. Results is in MEKEMEKE.

Parameters
[in,out]gOcean grid.
[in]gvOcean vertical grid structure.
csMEKE control structure.
mekeMEKE data.
[in]sn_uEady growth rate at u-points (s-1).
[in]sn_vEady growth rate at u-points (s-1).
[in]drag_rate_viscMean flow contrib. to drag rate
[in]i_massInverse of column mass.

Definition at line 568 of file MOM_MEKE.F90.

References meke_lengthscales_0d().

Referenced by step_forward_meke().

568  type(ocean_grid_type), intent(inout) :: g !< Ocean grid.
569  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure.
570  type(meke_cs), pointer :: cs !< MEKE control structure.
571  type(meke_type), pointer :: meke !< MEKE data.
572  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: sn_u !< Eady growth rate at u-points (s-1).
573  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: sn_v !< Eady growth rate at u-points (s-1).
574  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: drag_rate_visc !< Mean flow contrib. to drag rate
575  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: i_mass !< Inverse of column mass.
576  ! Local variables
577  real :: beta, sn, bottomfac2, barotrfac2, lmixscale, lrhines, leady
578  real :: i_h, khcoeff, kh, ubg2, cd2, drag_rate, ldamping, src
579  real :: eke, ekemin, ekemax, resid, resmin, resmax, ekeerr
580  integer :: i, j, is, ie, js, je, n1, n2
581  real, parameter :: tolerance = 1.e-12 ! Width of EKE bracket in m^2 s^-2.
582  logical :: usesecant, debugiteration
583 
584  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
585 
586  debugiteration = .false.
587  khcoeff = cs%MEKE_KhCoeff
588  ubg2 = cs%MEKE_Uscale**2
589  cd2 = cs%cdrag**2
590 
591 !$OMP do
592  do j=js,je ; do i=is,ie
593  !SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
594  ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
595  sn = min( min(sn_u(i,j) , sn_u(i-1,j)) , min(sn_v(i,j), sn_v(i,j-1)) )
596  beta = sqrt( g%dF_dx(i,j)**2 + g%dF_dy(i,j)**2 )
597  i_h = gv%Rho0 * i_mass(i,j)
598 
599  if (khcoeff*sn*i_h>0.) then
600  ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E
601  ekemin = 0. ! Use the trivial root as the left bracket
602  resmin = 0. ! Need to detect direction of left residual
603  ekemax = 0.01 ! First guess at right bracket
604  usesecant = .false. ! Start using a bisection method
605 
606  ! First find right bracket for which resid<0
607  resid = 1. ; n1 = 0
608  do while (resid>0.)
609  n1 = n1 + 1
610  eke = ekemax
611  call meke_lengthscales_0d(cs, g%areaT(i,j), beta, g%bathyT(i,j), &
612  meke%Rd_dx_h(i,j), sn, eke, &
613  bottomfac2, barotrfac2, lmixscale, &
614  lrhines, leady)
615  ! TODO: Should include resolution function in Kh
616  kh = (khcoeff * sqrt(2.*barotrfac2*eke) * lmixscale)
617  src = kh * (sn * sn)
618  drag_rate = i_h * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomfac2*eke + ubg2 ) )
619  ldamping = cs%MEKE_damping + drag_rate * bottomfac2
620  resid = src - ldamping * eke
621  if (debugiteration) then
622  write(0,*) n1, 'EKE=',eke,'resid=',resid
623  write(0,*) 'EKEmin=',ekemin,'ResMin=',resmin
624  write(0,*) 'src=',src,'ldamping=',ldamping
625  write(0,*) 'gamma-b=',bottomfac2,'gamma-t=',barotrfac2
626  write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',ubg2
627  endif
628  if (resid>0.) then ! EKE is to the left of the root
629  ekemin = eke ! so we move the left bracket here
630  ekemax = 10. * eke ! and guess again for the right bracket
631  if (resid<resmin) usesecant = .true.
632  resmin = resid
633  if (ekemax > 2.e17) then
634  if (debugiteration) stop 'Something has gone very wrong'
635  debugiteration = .true.
636  resid = 1. ; n1 = 0
637  ekemin = 0. ; resmin = 0.
638  ekemax = 0.01
639  usesecant = .false.
640  endif
641  endif
642  enddo ! while(resid>0.) searching for right bracket
643  resmax = resid
644 
645  ! Bisect the bracket
646  n2 = 0 ; ekeerr = ekemax - ekemin
647  do while (ekeerr>tolerance)
648  n2 = n2 + 1
649  if (usesecant) then
650  eke = ekemin + (ekemax - ekemin) * (resmin / (resmin - resmax))
651  else
652  eke = 0.5 * (ekemin + ekemax)
653  endif
654  ekeerr = min( eke-ekemin, ekemax-eke )
655  ! TODO: Should include resolution function in Kh
656  kh = (khcoeff * sqrt(2.*barotrfac2*eke) * lmixscale)
657  src = kh * (sn * sn)
658  drag_rate = i_h * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomfac2*eke + ubg2 ) )
659  ldamping = cs%MEKE_damping + drag_rate * bottomfac2
660  resid = src - ldamping * eke
661  if (usesecant.and.resid>resmin) usesecant = .false.
662  if (resid>0.) then ! EKE is to the left of the root
663  ekemin = eke ! so we move the left bracket here
664  if (resid<resmin) usesecant = .true.
665  resmin = resid ! Save this for the secant method
666  elseif (resid<0.) then ! EKE is to the right of the root
667  ekemax = eke ! so we move the right bracket here
668  resmax = resid ! Save this for the secant method
669  else
670  exit ! resid=0 => EKE is exactly at the root
671  endif
672  if (n2>200) stop 'Failing to converge?'
673  enddo ! while(EKEmax-EKEmin>tolerance)
674 
675  else
676  eke = 0.
677  endif
678  meke%MEKE(i,j) = eke
679  enddo ; enddo
680 
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ meke_init()

logical function, public mom_meke::meke_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(meke_cs), pointer  CS,
type(meke_type), pointer  MEKE,
type(mom_restart_cs), pointer  restart_CS 
)

Initializes the MOM_MEKE module and reads parameters. Returns True if module is to be used, otherwise returns False.

Parameters
[in]timeThe current model time.
[in,out]gThe ocean's grid structure.
[in]param_fileParameter file parser structure.
[in,out]diagDiagnostics structure.
csMEKE control structure.
mekeMEKE-related fields.
restart_csRestart control structure for MOM_MEKE.

Definition at line 789 of file MOM_MEKE.F90.

References mom_error_handler::mom_error(), and mom_error_handler::mom_mesg().

Referenced by mom::initialize_mom().

789  type(time_type), intent(in) :: time !< The current model time.
790  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
791  type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
792  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
793  type(meke_cs), pointer :: cs !< MEKE control structure.
794  type(meke_type), pointer :: meke !< MEKE-related fields.
795  type(mom_restart_cs), pointer :: restart_cs !< Restart control structure for MOM_MEKE.
796 ! Local variables
797  integer :: is, ie, js, je, isd, ied, jsd, jed, nz
798  logical :: laplacian, usevarmix, coldstart
799 ! This include declares and sets the variable "version".
800 #include "version_variable.h"
801  character(len=40) :: mdl = "MOM_MEKE" ! This module's name.
802 
803  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
804  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
805 
806  ! Determine whether this module will be used
807  call log_version(param_file, mdl, version, "")
808  call get_param(param_file, mdl, "USE_MEKE", meke_init, &
809  "If true, turns on the MEKE scheme which calculates\n"// &
810  "a sub-grid mesoscale eddy kinetic energy budget.", &
811  default=.false.)
812  if (.not. meke_init) return
813 
814  if (.not. associated(meke)) then
815  ! The MEKE structure should have been allocated in MEKE_alloc_register_restart()
816  call mom_error(warning, "MEKE_init called with NO associated "// &
817  "MEKE-type structure.")
818  return
819  endif
820  if (associated(cs)) then
821  call mom_error(warning, &
822  "MEKE_init called with an associated control structure.")
823  return
824  else ; allocate(cs) ; endif
825 
826  call mom_mesg("MEKE_init: reading parameters ", 5)
827 
828  ! Read all relevant parameters and write them to the model log.
829  call get_param(param_file, mdl, "MEKE_DAMPING", cs%MEKE_damping, &
830  "The local depth-indepented MEKE dissipation rate.", &
831  units="s-1", default=0.0)
832  call get_param(param_file, mdl, "MEKE_CD_SCALE", cs%MEKE_Cd_scale, &
833  "The ratio of the bottom eddy velocity to the column mean\n"//&
834  "eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1\n"//&
835  "to account for the surface intensification of MEKE.", &
836  units="nondim", default=0.)
837  call get_param(param_file, mdl, "MEKE_CB", cs%MEKE_Cb, &
838  "A coefficient in the expression for the ratio of bottom projected\n"//&
839  "eddy energy and mean column energy (see Jansen et al. 2015).",&
840  units="nondim", default=25.)
841  call get_param(param_file, mdl, "MEKE_MIN_GAMMA2", cs%MEKE_min_gamma, &
842  "The minimum allowed value of gamma_b^2.",&
843  units="nondim", default=0.0001)
844  call get_param(param_file, mdl, "MEKE_CT", cs%MEKE_Ct, &
845  "A coefficient in the expression for the ratio of barotropic\n"//&
846  "eddy energy and mean column energy (see Jansen et al. 2015).",&
847  units="nondim", default=50.)
848  call get_param(param_file, mdl, "MEKE_GMCOEFF", cs%MEKE_GMcoeff, &
849  "The efficiency of the conversion of potential energy \n"//&
850  "into MEKE by the thickness mixing parameterization. \n"//&
851  "If MEKE_GMCOEFF is negative, this conversion is not \n"//&
852  "used or calculated.", units="nondim", default=-1.0)
853  call get_param(param_file, mdl, "MEKE_FRCOEFF", cs%MEKE_FrCoeff, &
854  "The efficiency of the conversion of mean energy into \n"//&
855  "MEKE. If MEKE_FRCOEFF is negative, this conversion \n"//&
856  "is not used or calculated.", units="nondim", default=-1.0)
857  call get_param(param_file, mdl, "MEKE_BGSRC", cs%MEKE_BGsrc, &
858  "A background energy source for MEKE.", units="W kg-1", &
859  default=0.0)
860  call get_param(param_file, mdl, "MEKE_KH", cs%MEKE_Kh, &
861  "A background lateral diffusivity of MEKE.\n"//&
862  "Use a negative value to not apply lateral diffusion to MEKE.", &
863  units="m2 s-1", default=-1.0)
864  call get_param(param_file, mdl, "MEKE_K4", cs%MEKE_K4, &
865  "A lateral bi-harmonic diffusivity of MEKE.\n"//&
866  "Use a negative value to not apply bi-harmonic diffusion to MEKE.", &
867  units="m4 s-1", default=-1.0)
868  call get_param(param_file, mdl, "MEKE_DTSCALE", cs%MEKE_dtScale, &
869  "A scaling factor to accelerate the time evolution of MEKE.", &
870  units="nondim", default=1.0)
871  call get_param(param_file, mdl, "MEKE_KHCOEFF", cs%MEKE_KhCoeff, &
872  "A scaling factor in the expression for eddy diffusivity\n"//&
873  "which is otherwise proportional to the MEKE velocity-\n"//&
874  "scale times an eddy mixing-length. This factor\n"//&
875  "must be >0 for MEKE to contribute to the thickness/\n"//&
876  "and tracer diffusivity in the rest of the model.", &
877  units="nondim", default=1.0)
878  call get_param(param_file, mdl, "MEKE_USCALE", cs%MEKE_Uscale, &
879  "The background velocity that is combined with MEKE to \n"//&
880  "calculate the bottom drag.", units="m s-1", default=0.0)
881  call get_param(param_file, mdl, "MEKE_VISC_DRAG", cs%visc_drag, &
882  "If true, use the vertvisc_type to calculate the bottom \n"//&
883  "drag acting on MEKE.", default=.true.)
884  call get_param(param_file, mdl, "MEKE_KHTH_FAC", meke%KhTh_fac, &
885  "A factor that maps MEKE%Kh to KhTh.", units="nondim", &
886  default=0.0)
887  call get_param(param_file, mdl, "MEKE_KHTR_FAC", meke%KhTr_fac, &
888  "A factor that maps MEKE%Kh to KhTr.", units="nondim", &
889  default=0.0)
890  call get_param(param_file, mdl, "MEKE_KHMEKE_FAC", cs%KhMEKE_Fac, &
891  "A factor that maps MEKE%Kh to Kh for MEKE itself.", &
892  units="nondim", default=0.0)
893  call get_param(param_file, mdl, "MEKE_OLD_LSCALE", cs%use_old_lscale, &
894  "If true, use the old formula for length scale which is\n"//&
895  "a function of grid spacing and deformation radius.", &
896  default=.false.)
897  call get_param(param_file, mdl, "MEKE_RD_MAX_SCALE", cs%Rd_as_max_scale, &
898  "If true, the length scale used by MEKE is the minimum of\n"//&
899  "the deformation radius or grid-spacing. Only used if\n"//&
900  "MEKE_OLD_LSCALE=True", units="nondim", default=.false.)
901  call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF", cs%viscosity_coeff, &
902  "If non-zero, is the scaling coefficient in the expression for\n"//&
903  "viscosity used to parameterize lateral momentum mixing by\n"//&
904  "unresolved eddies represented by MEKE. Can be negative to\n"//&
905  "represent backscatter from the unresolved eddies.", &
906  units="nondim", default=0.0)
907  call get_param(param_file, mdl, "MEKE_FIXED_MIXING_LENGTH", cs%Lfixed, &
908  "If positive, is a fixed length contribution to the expression\n"//&
909  "for mixing length used in MEKE-derived diffusiviity.", &
910  units="m", default=0.0)
911  call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", cs%aDeform, &
912  "If positive, is a coefficient weighting the deformation scale\n"//&
913  "in the expression for mixing length used in MEKE-derived diffusiviity.", &
914  units="nondim", default=0.0)
915  call get_param(param_file, mdl, "MEKE_ALPHA_RHINES", cs%aRhines, &
916  "If positive, is a coefficient weighting the Rhines scale\n"//&
917  "in the expression for mixing length used in MEKE-derived diffusiviity.", &
918  units="nondim", default=0.05)
919  call get_param(param_file, mdl, "MEKE_ALPHA_EADY", cs%aEady, &
920  "If positive, is a coefficient weighting the Eady length scale\n"//&
921  "in the expression for mixing length used in MEKE-derived diffusiviity.", &
922  units="nondim", default=0.05)
923  call get_param(param_file, mdl, "MEKE_ALPHA_FRICT", cs%aFrict, &
924  "If positive, is a coefficient weighting the frictional arrest scale\n"//&
925  "in the expression for mixing length used in MEKE-derived diffusiviity.", &
926  units="nondim", default=0.0)
927  call get_param(param_file, mdl, "MEKE_ALPHA_GRID", cs%aGrid, &
928  "If positive, is a coefficient weighting the grid-spacing as a scale\n"//&
929  "in the expression for mixing length used in MEKE-derived diffusiviity.", &
930  units="nondim", default=0.0)
931  call get_param(param_file, mdl, "MEKE_COLD_START", coldstart, &
932  "If true, initialize EKE to zero. Otherwise a local equilibrium solution\n"//&
933  "is used as an initial condition for EKE.", default=.false.)
934  call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", meke%backscatter_Ro_c, &
935  "The coefficient in the Rossby number function for scaling the buharmonic\n"//&
936  "frictional energy source. Setting to non-zero enables the Rossby number function.", &
937  units="nondim", default=0.0)
938  call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_POW", meke%backscatter_Ro_pow, &
939  "The power in the Rossby number function for scaling the biharmomnic\n"//&
940  "frictional energy source.", units="nondim", default=0.0)
941  call get_param(param_file, mdl, "MEKE_ADVECTION_FACTOR", cs%MEKE_advection_factor, &
942  "A scale factor in front of advection of eddy energy. Zero turns advection off.\n"//&
943  "Using unity would be normal but other values could accomodate a mismatch\n"//&
944  "between the advecting barotropic flow and the vertical structure of MEKE.", &
945  units="nondim", default=0.0)
946 
947  ! Nonlocal module parameters
948  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
949  "CDRAG is the drag coefficient relating the magnitude of \n"//&
950  "the velocity field to the bottom stress.", units="nondim", &
951  default=0.003)
952  call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
953  if (cs%viscosity_coeff/=0. .and. .not. laplacian) call mom_error(fatal, &
954  "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF is true.")
955 
956  call get_param(param_file, mdl, "USE_VARIABLE_MIXING", usevarmix, default=.false., do_not_log=.true.)
957  if (.not. usevarmix .and. cs%aEady>0.) call mom_error(fatal, &
958  "USE_VARIABLE_MIXING must be true if USE_MEKE is true and MEKE_ALPHA_EADY>0.")
959  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
960 
961  ! Allocation of storage NOT shared with other modules
962  if (cs%MEKE_K4>=0.) then
963  allocate(cs%del2MEKE(isd:ied,jsd:jed)) ; cs%del2MEKE(:,:) = 0.0
964  endif
965 
966 ! In the case of a restart, these fields need a halo update
967  if (associated(meke%MEKE)) then
968  call create_group_pass(cs%pass_MEKE, meke%MEKE, g%Domain)
969  call do_group_pass(cs%pass_MEKE, g%Domain)
970  endif
971  if (associated(meke%Kh)) then
972  call create_group_pass(cs%pass_Kh, meke%Kh, g%Domain)
973  call do_group_pass(cs%pass_Kh, g%Domain)
974  endif
975  if (associated(meke%Ku)) then
976  call create_group_pass(cs%pass_Ku, meke%Ku, g%Domain)
977  call do_group_pass(cs%pass_Ku, g%Domain)
978  endif
979  if (allocated(cs%del2MEKE)) then
980  call create_group_pass(cs%pass_del2MEKE, cs%del2MEKE, g%Domain)
981  call do_group_pass(cs%pass_del2MEKE, g%Domain)
982  endif
983 
984 ! Register fields for output from this module.
985  cs%diag => diag
986  cs%id_MEKE = register_diag_field('ocean_model', 'MEKE', diag%axesT1, time, &
987  'Mesoscale Eddy Kinetic Energy', 'meter2 second-2')
988  if (.not. associated(meke%MEKE)) cs%id_MEKE = -1
989  cs%id_Kh = register_diag_field('ocean_model', 'MEKE_KH', diag%axesT1, time, &
990  'MEKE derived diffusivity', 'meter2 second-1')
991  if (.not. associated(meke%Kh)) cs%id_Kh = -1
992  cs%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, time, &
993  'MEKE derived lateral viscosity', 'meter2 second-1')
994  if (.not. associated(meke%Ku)) cs%id_Ku = -1
995  cs%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, time, &
996  'MEKE derived eddy-velocity scale', 'meter second-1')
997  if (.not. associated(meke%MEKE)) cs%id_Ue = -1
998  cs%id_Ub = register_diag_field('ocean_model', 'MEKE_Ub', diag%axesT1, time, &
999  'MEKE derived bottom eddy-velocity scale', 'meter second-1')
1000  if (.not. associated(meke%MEKE)) cs%id_Ub = -1
1001  cs%id_Ut = register_diag_field('ocean_model', 'MEKE_Ut', diag%axesT1, time, &
1002  'MEKE derived barotropic eddy-velocity scale', 'meter second-1')
1003  if (.not. associated(meke%MEKE)) cs%id_Ut = -1
1004  cs%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, time, &
1005  'MEKE energy source', 'meter2 second-3')
1006  cs%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, time, &
1007  'MEKE decay rate', 'second-1')
1008  cs%id_KhMEKE_u = register_diag_field('ocean_model', 'KHMEKE_u', diag%axesCu1, time, &
1009  'Zonal diffusivity of MEKE', 'meter2 second-1')
1010  cs%id_KhMEKE_v = register_diag_field('ocean_model', 'KHMEKE_v', diag%axesCv1, time, &
1011  'Meridional diffusivity of MEKE', 'meter2 second-1')
1012  cs%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, time, &
1013  'MEKE energy available from thickness mixing', 'Watt meter-2')
1014  if (.not. associated(meke%GM_src)) cs%id_GM_src = -1
1015  cs%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, time, &
1016  'MEKE energy available from momentum', 'Watt meter-2')
1017  if (.not. associated(meke%mom_src)) cs%id_mom_src = -1
1018  cs%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, time, &
1019  'Eddy mixing length used in the MEKE derived eddy diffusivity', 'meter')
1020  cs%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, time, &
1021  'Rhines length scale used in the MEKE derived eddy diffusivity', 'meter')
1022  cs%id_Leady = register_diag_field('ocean_model', 'MEKE_Leady', diag%axesT1, time, &
1023  'Eady length scale used in the MEKE derived eddy diffusivity', 'meter')
1024  cs%id_gamma_b = register_diag_field('ocean_model', 'MEKE_gamma_b', diag%axesT1, time, &
1025  'Ratio of bottom-projected eddy velocity to column-mean eddy velocity', 'nondim')
1026  cs%id_gamma_t = register_diag_field('ocean_model', 'MEKE_gamma_t', diag%axesT1, time, &
1027  'Ratio of barotropic eddy velocity to column-mean eddy velocity', 'nondim')
1028 
1029  cs%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=clock_routine)
1030 
1031  ! Detect whether this instant of MEKE_init() is at the beginning of a run
1032  ! or after a restart. If at the beginning, we will initialize MEKE to a local
1033  ! equilibrium.
1034  cs%initialize = .not.query_initialized(meke%MEKE,"MEKE",restart_cs)
1035  if (coldstart) cs%initialize = .false.
1036  if (cs%initialize) call mom_error(warning, &
1037  "MEKE_init: Initializing MEKE with a local equilibrium balance.")
1038 
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ meke_lengthscales()

subroutine mom_meke::meke_lengthscales ( type(meke_cs), pointer  CS,
type(meke_type), pointer  MEKE,
type(ocean_grid_type), intent(inout)  G,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in)  SN_u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in)  SN_v,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  EKE,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  bottomFac2,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  barotrFac2,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  LmixScale 
)
private

Calculates the eddy mixing length scale and \(\gamma_b\) and \(\gamma_t\) functions that are ratios of either bottom or barotropic eddy energy to the column eddy energy, respectively. See MEKE equations.

Parameters
csMEKE control structure.
mekeMEKE data.
[in,out]gOcean grid.
[in]sn_uEady growth rate at u-points (s-1).
[in]sn_vEady growth rate at u-points (s-1).
[in]ekeEddy kinetic energy (m2/s2).
[out]bottomfac2gamma_b^2
[out]barotrfac2gamma_t^2
[out]lmixscaleEddy mixing length (m).

Definition at line 689 of file MOM_MEKE.F90.

References meke_lengthscales_0d().

Referenced by step_forward_meke().

689  type(meke_cs), pointer :: cs !< MEKE control structure.
690  type(meke_type), pointer :: meke !< MEKE data.
691  type(ocean_grid_type), intent(inout) :: g !< Ocean grid.
692  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: sn_u !< Eady growth rate at u-points (s-1).
693  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: sn_v !< Eady growth rate at u-points (s-1).
694  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eke !< Eddy kinetic energy (m2/s2).
695  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bottomfac2 !< gamma_b^2
696  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: barotrfac2 !< gamma_t^2
697  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: lmixscale !< Eddy mixing length (m).
698  ! Local variables
699  real, dimension(SZI_(G),SZJ_(G)) :: lrhines, leady
700  real :: beta, sn
701  integer :: i, j, is, ie, js, je
702 
703  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
704 
705 !$OMP do
706  do j=js,je ; do i=is,ie
707  if (.not.cs%use_old_lscale) then
708  if (cs%aEady > 0.) then
709  sn = 0.25*( (sn_u(i,j) + sn_u(i-1,j)) + (sn_v(i,j) + sn_v(i,j-1)) )
710  else
711  sn = 0.
712  endif
713  beta = sqrt( g%dF_dx(i,j)**2 + g%dF_dy(i,j)**2 )
714  endif
715  ! Returns bottomFac2, barotrFac2 and LmixScale
716  call meke_lengthscales_0d(cs, g%areaT(i,j), beta, g%bathyT(i,j), &
717  meke%Rd_dx_h(i,j), sn, meke%MEKE(i,j), &
718  bottomfac2(i,j), barotrfac2(i,j), lmixscale(i,j), &
719  lrhines(i,j), leady(i,j))
720  enddo ; enddo
721  if (cs%id_Lrhines>0) call post_data(cs%id_Lrhines, lrhines, cs%diag)
722  if (cs%id_Leady>0) call post_data(cs%id_Leady, leady, cs%diag)
723 
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ meke_lengthscales_0d()

subroutine mom_meke::meke_lengthscales_0d ( type(meke_cs), pointer  CS,
real, intent(in)  area,
real, intent(in)  beta,
real, intent(in)  depth,
real, intent(in)  Rd_dx,
real, intent(in)  SN,
real, intent(in)  EKE,
real, intent(out)  bottomFac2,
real, intent(out)  barotrFac2,
real, intent(out)  LmixScale,
real, intent(out)  Lrhines,
real, intent(out)  Leady 
)
private

Calculates the eddy mixing length scale and \(\gamma_b\) and \(\gamma_t\) functions that are ratios of either bottom or barotropic eddy energy to the column eddy energy, respectively. See MEKE equations.

Parameters
csMEKE control structure.
[in]areaGrid cell area (m2)
[in]betaPlanetary beta = |grad F| (s-1 m-1)
[in]depthOcean depth (m)
[in]rd_dxResolution Ld/dx (nondim).
[in]snEady growth rate (s-1).
[in]ekeEddy kinetic energy (m s-1).
[out]bottomfac2gamma_b^2
[out]barotrfac2gamma_t^2
[out]lmixscaleEddy mixing length (m).
[out]lrhinesRhines length scale (m).
[out]leadyEady length scale (m).

Definition at line 731 of file MOM_MEKE.F90.

Referenced by meke_equilibrium(), and meke_lengthscales().

731  type(meke_cs), pointer :: cs !< MEKE control structure.
732  real, intent(in) :: area !< Grid cell area (m2)
733  real, intent(in) :: beta !< Planetary beta = |grad F| (s-1 m-1)
734  real, intent(in) :: depth !< Ocean depth (m)
735  real, intent(in) :: rd_dx !< Resolution Ld/dx (nondim).
736  real, intent(in) :: sn !< Eady growth rate (s-1).
737  real, intent(in) :: eke !< Eddy kinetic energy (m s-1).
738  real, intent(out) :: bottomfac2 !< gamma_b^2
739  real, intent(out) :: barotrfac2 !< gamma_t^2
740  real, intent(out) :: lmixscale !< Eddy mixing length (m).
741  real, intent(out) :: lrhines !< Rhines length scale (m).
742  real, intent(out) :: leady !< Eady length scale (m).
743  ! Local variables
744  real :: lgrid, ldeform, ldeformlim, ue, lfrict
745 
746  ! Length scale for MEKE derived diffusivity
747  lgrid = sqrt(area) ! Grid scale
748  ldeform = lgrid * rd_dx ! Deformation scale
749  lfrict = depth / cs%cdrag ! Frictional arrest scale
750  ! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy
751  ! used in calculating bottom drag
752  bottomfac2 = cs%MEKE_CD_SCALE**2
753  if (lfrict*cs%MEKE_Cb>0.) bottomfac2 = bottomfac2 + 1./( 1. + cs%MEKE_Cb*(ldeform/lfrict) )**0.8
754  bottomfac2 = max(bottomfac2, cs%MEKE_min_gamma)
755  ! gamma_t^2 is the ratio of barotropic eddy energy to mean column eddy energy
756  ! used in the velocity scale for diffusivity
757  barotrfac2 = 1.
758  if (lfrict*cs%MEKE_Ct>0.) barotrfac2 = 1./( 1. + cs%MEKE_Ct*(ldeform/lfrict) )**0.25
759  barotrfac2 = max(barotrfac2, cs%MEKE_min_gamma)
760  if (cs%use_old_lscale) then
761  if (cs%Rd_as_max_scale) then
762  lmixscale = min(ldeform, lgrid) ! The smaller of Ld or dx
763  else
764  lmixscale = lgrid
765  endif
766  else
767  ue = sqrt( 2.0 * max( 0., barotrfac2*eke ) ) ! Barotropic eddy flow scale
768  lrhines = sqrt( ue / max( beta, 1.e-30 ) ) ! Rhines scale
769  if (cs%aEady > 0.) then
770  leady = ue / max( sn, 1.e-15 ) ! Bound Eady time-scale < 1e15 seconds
771  else
772  leady = 0.
773  endif
774  lmixscale = 0.
775  if (cs%aDeform*ldeform > 0.) lmixscale = lmixscale + 1./(cs%aDeform*ldeform)
776  if (cs%aFrict *lfrict > 0.) lmixscale = lmixscale + 1./(cs%aFrict *lfrict)
777  if (cs%aRhines*lrhines > 0.) lmixscale = lmixscale + 1./(cs%aRhines*lrhines)
778  if (cs%aEady *leady > 0.) lmixscale = lmixscale + 1./(cs%aEady *leady)
779  if (cs%aGrid *lgrid > 0.) lmixscale = lmixscale + 1./(cs%aGrid *lgrid)
780  if (cs%Lfixed > 0.) lmixscale = lmixscale + 1./cs%Lfixed
781  if (lmixscale > 0.) lmixscale = 1. / lmixscale
782  endif
783 
Here is the caller graph for this function:

◆ step_forward_meke()

subroutine, public mom_meke::step_forward_meke ( type(meke_type), pointer  MEKE,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in)  SN_u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in)  SN_v,
type(vertvisc_type), intent(in)  visc,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(meke_cs), pointer  CS,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  hu,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  hv 
)

Integrates forward-in-time the MEKE eddy energy equation. See MEKE equations.

Parameters
mekeMEKE data.
[in,out]gOcean grid.
[in]gvOcean vertical grid structure.
[in]hLayer thickness (m or kg m-2).
[in]sn_uEady growth rate at u-points (s-1).
[in]sn_vEady growth rate at u-points (s-1).
[in]viscThe vertical viscosity type.
[in]dtModel(baroclinic) time-step (s).
csMEKE control structure.
[in]huZonal flux flux (m3).
[in]hvMeridional mass flux (m3).

Definition at line 88 of file MOM_MEKE.F90.

References meke_equilibrium(), meke_lengthscales(), and mom_error_handler::mom_error().

Referenced by mom::step_mom().

88  type(meke_type), pointer :: meke !< MEKE data.
89  type(ocean_grid_type), intent(inout) :: g !< Ocean grid.
90  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure.
91  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg m-2).
92  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: sn_u !< Eady growth rate at u-points (s-1).
93  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: sn_v !< Eady growth rate at u-points (s-1).
94  type(vertvisc_type), intent(in) :: visc !< The vertical viscosity type.
95  real, intent(in) :: dt !< Model(baroclinic) time-step (s).
96  type(meke_cs), pointer :: cs !< MEKE control structure.
97  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: hu !< Zonal flux flux (m3).
98  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: hv !< Meridional mass flux (m3).
99  ! Local variables
100  real, dimension(SZI_(G),SZJ_(G)) :: &
101  mass, & ! The total mass of the water column, in kg m-2.
102  i_mass, & ! The inverse of mass, in m2 kg-1.
103  src, & ! The sum of all MEKE sources, in m2 s-3.
104  meke_decay, & ! The MEKE decay timescale, in s-1.
105  meke_gm_src, & ! The MEKE source from thickness mixing, in m2 s-3.
106  meke_mom_src, & ! The MEKE source from momentum, in m2 s-3.
107  drag_rate_visc, &
108  drag_rate, & ! The MEKE spindown timescale due to bottom drag, in s-1.
109  lmixscale, & ! Square of eddy mixing length, in m2.
110  barotrfac2, & ! Ratio of EKE_barotropic / EKE (nondim)/
111  bottomfac2 ! Ratio of EKE_bottom / EKE (nondim)/
112  real, dimension(SZIB_(G),SZJ_(G)) :: &
113  meke_uflux, & ! The zonal diffusive flux of MEKE, in kg m2 s-3.
114  kh_u, & ! The zonal diffusivity that is actually used, in m2 s-1.
115  barohu, & ! Depth integrated zonal mass flux (m3).
116  drag_vel_u ! A (vertical) viscosity associated with bottom drag at
117  ! u-points, in m s-1.
118  real, dimension(SZI_(G),SZJB_(G)) :: &
119  meke_vflux, & ! The meridional diffusive flux of MEKE, in kg m2 s-3.
120  kh_v, & ! The meridional diffusivity that is actually used, in m2 s-1.
121  barohv, & ! Depth integrated meridional mass flux (m3).
122  drag_vel_v ! A (vertical) viscosity associated with bottom drag at
123  ! v-points, in m s-1.
124  real :: kh_here, inv_kh_max, k4_here
125  real :: cdrag2
126  real :: advfac
127  real :: mass_neglect ! A negligible mass, in kg m-2.
128  real :: ldamping ! The MEKE damping rate in s-1.
129  real :: rho0 ! A density used to convert mass to distance, in kg m-3.
130  real :: sdt ! dt to use locally (could be scaled to accelerate)
131  real :: sdt_damp ! dt for damping (sdt could be split).
132  logical :: use_drag_rate ! Flag to indicate drag_rate is finite
133  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
134 
135  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
136  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
137 
138  if (.not.associated(cs)) call mom_error(fatal, &
139  "MOM_MEKE: Module must be initialized before it is used.")
140  if (.not.associated(meke)) call mom_error(fatal, &
141  "MOM_MEKE: MEKE must be initialized before it is used.")
142 
143  rho0 = gv%H_to_kg_m2 * gv%m_to_H
144  mass_neglect = gv%H_to_kg_m2 * gv%H_subroundoff
145  sdt = dt*cs%MEKE_dtScale ! Scaled dt to use for time-stepping
146  if (cs%MEKE_damping + cs%MEKE_Cd_scale > 0.0 .or. cs%MEKE_Cb>0. &
147  .or. cs%visc_drag) then
148  use_drag_rate = .true.
149  else
150  use_drag_rate = .false.
151  endif
152 
153  ! Only integrate the MEKE equations if MEKE is required.
154  if (associated(meke%MEKE)) then
155 
156  if (cs%debug) then
157  if (associated(meke%mom_src)) call hchksum(meke%mom_src, 'MEKE mom_src',g%HI)
158  if (associated(meke%GM_src)) call hchksum(meke%GM_src, 'MEKE GM_src',g%HI)
159  if (associated(meke%MEKE)) call hchksum(meke%MEKE, 'MEKE MEKE',g%HI)
160  call uvchksum("MEKE SN_[uv]", sn_u, sn_v, g%HI)
161  call uvchksum("MEKE h[uv]", hu, hv, g%HI, haloshift=1)
162  endif
163 
164  ! Why are these 3 lines repeated from above?
165  sdt = dt*cs%MEKE_dtScale ! Scaled dt to use for time-stepping
166  rho0 = gv%H_to_kg_m2 * gv%m_to_H
167  mass_neglect = gv%H_to_kg_m2 * gv%H_subroundoff
168  cdrag2 = cs%cdrag**2
169 
170  ! With a depth-dependent (and possibly strong) damping, it seems
171  ! advisable to use Strang splitting between the damping and diffusion.
172  sdt_damp = sdt ; if (cs%MEKE_KH >= 0.0 .or. cs%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt
173 
174  ! Calculate depth integrated mass flux if doing advection
175  if (cs%MEKE_advection_factor>0.) then
176  do j=js,je ; do i=is-1,ie
177  barohu(i,j) = 0.
178  enddo ; enddo
179  do k=1,nz
180  do j=js,je ; do i=is-1,ie
181  barohu(i,j) = hu(i,j,k)
182  enddo ; enddo
183  enddo
184  do j=js-1,je ; do i=is,ie
185  barohv(i,j) = 0.
186  enddo ; enddo
187  do k=1,nz
188  do j=js-1,je ; do i=is,ie
189  barohv(i,j) = hv(i,j,k)
190  enddo ; enddo
191  enddo
192  endif
193 
194 !$OMP parallel default(none) shared(MEKE,CS,is,ie,js,je,nz,src,mass,G,GV,h,I_mass, &
195 !$OMP sdt,drag_vel_u,visc,drag_vel_v,drag_rate_visc, &
196 !$OMP drag_rate,Rho0,MEKE_decay,sdt_damp,cdrag2, &
197 !$OMP bottomFac2) &
198 !$OMP private(ldamping)
199 
200  if (cs%MEKE_Cd_scale == 0.0 .and. .not. cs%visc_drag) then
201 !$OMP do
202  do j=js,je ; do i=is,ie
203  drag_rate(i,j) = 0.
204  enddo ; enddo
205  endif
206 
207  ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow
208  if (cs%visc_drag) then
209 !$OMP do
210  do j=js,je ; do i=is-1,ie
211  drag_vel_u(i,j) = 0.0
212  if ((g%mask2dCu(i,j) > 0.0) .and. (visc%bbl_thick_u(i,j) > 0.0)) &
213  drag_vel_u(i,j) = visc%kv_bbl_u(i,j) / visc%bbl_thick_u(i,j)
214  enddo ; enddo
215 !$OMP do
216  do j=js-1,je ; do i=is,ie
217  drag_vel_v(i,j) = 0.0
218  if ((g%mask2dCv(i,j) > 0.0) .and. (visc%bbl_thick_v(i,j) > 0.0)) &
219  drag_vel_v(i,j) = visc%kv_bbl_v(i,j) / visc%bbl_thick_v(i,j)
220  enddo ; enddo
221 
222 !$OMP do
223  do j=js,je ; do i=is,ie
224  drag_rate_visc(i,j) = (0.25*g%IareaT(i,j) * &
225  ((g%areaCu(i-1,j)*drag_vel_u(i-1,j) + &
226  g%areaCu(i,j)*drag_vel_u(i,j)) + &
227  (g%areaCv(i,j-1)*drag_vel_v(i,j-1) + &
228  g%areaCv(i,j)*drag_vel_v(i,j)) ) )
229  enddo ; enddo
230  else
231 !$OMP do
232  do j=js,je ; do i=is,ie
233  drag_rate_visc(i,j) = 0.
234  enddo ; enddo
235  endif
236 
237 !$OMP do
238  do j=js-1,je+1
239  do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo
240  do k=1,nz ; do i=is-1,ie+1
241  mass(i,j) = mass(i,j) + g%mask2dT(i,j) * (gv%H_to_kg_m2 * h(i,j,k))
242  enddo ; enddo
243  do i=is-1,ie+1
244  i_mass(i,j) = 0.0
245  if (mass(i,j) > 0.0) i_mass(i,j) = 1.0 / mass(i,j)
246  enddo
247  enddo
248 !$OMP end parallel
249 
250  if (cs%initialize) then
251  call meke_equilibrium(cs, meke, g, gv, sn_u, sn_v, drag_rate_visc, i_mass)
252  cs%initialize = .false.
253  endif
254 
255  ! Calculates bottomFac2, barotrFac2 and LmixScale
256  call meke_lengthscales(cs, meke, g, sn_u, sn_v, meke%MEKE, bottomfac2, barotrfac2, lmixscale)
257  if (cs%debug) then
258  call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, g%HI)
259  call hchksum(mass, 'MEKE mass',g%HI,haloshift=1)
260  call hchksum(drag_rate_visc, 'MEKE drag_rate_visc',g%HI)
261  call hchksum(bottomfac2, 'MEKE bottomFac2',g%HI)
262  call hchksum(barotrfac2, 'MEKE barotrFac2',g%HI)
263  call hchksum(lmixscale, 'MEKE LmixScale',g%HI)
264  endif
265 
266 !$OMP parallel default(none) shared(MEKE,CS,is,ie,js,je,nz,src,mass,G,h,I_mass, &
267 !$OMP sdt,drag_vel_u,visc,drag_vel_v,drag_rate_visc, &
268 !$OMP drag_rate,Rho0,MEKE_decay,sdt_damp,cdrag2, &
269 !$OMP bottomFac2,barotrFac2,use_drag_rate) &
270 !$OMP private(ldamping)
271 
272  ! Aggregate sources of MEKE (background, frictional and GM)
273 !$OMP do
274  do j=js,je ; do i=is,ie
275  src(i,j) = cs%MEKE_BGsrc
276  enddo ; enddo
277 
278  if (associated(meke%mom_src)) then
279 !$OMP do
280  do j=js,je ; do i=is,ie
281  src(i,j) = src(i,j) - cs%MEKE_FrCoeff*i_mass(i,j)*meke%mom_src(i,j)
282  enddo ; enddo
283  endif
284 
285  if (associated(meke%GM_src)) then
286 !$OMP do
287  do j=js,je ; do i=is,ie
288  src(i,j) = src(i,j) - cs%MEKE_GMcoeff*i_mass(i,j)*meke%GM_src(i,j)
289  enddo ; enddo
290  endif
291 
292  ! Increase EKE by a full time-steps worth of source
293 !$OMP do
294  do j=js,je ; do i=is,ie
295  meke%MEKE(i,j) = (meke%MEKE(i,j) + sdt*src(i,j) )*g%mask2dT(i,j)
296  enddo ; enddo
297 
298  if (use_drag_rate) then
299  ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
300 !$OMP do
301  do j=js,je ; do i=is,ie
302  drag_rate(i,j) = (rho0 * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
303  + cdrag2 * ( max(0.0, 2.0*bottomfac2(i,j)*meke%MEKE(i,j)) + cs%MEKE_Uscale**2 ) )
304  enddo ; enddo
305  endif
306 
307  ! First stage of Strang splitting
308 !$OMP do
309  do j=js,je ; do i=is,ie
310  ldamping = cs%MEKE_damping + drag_rate(i,j) * bottomfac2(i,j)
311  if (meke%MEKE(i,j)<0.) ldamping = 0.
312  ! notice that the above line ensures a damping only if MEKE is positive,
313  ! while leaving MEKE unchanged if it is negative
314  meke%MEKE(i,j) = meke%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
315  meke_decay(i,j) = ldamping*g%mask2dT(i,j)
316  enddo ; enddo
317 !$OMP end parallel
318 
319  if (cs%MEKE_KH >= 0.0 .or. cs%KhMEKE_FAC > 0.0 .or. cs%MEKE_K4 >= 0.0) then
320  ! Update halos for lateral or bi-harmonic diffusion
321  call cpu_clock_begin(cs%id_clock_pass)
322  call do_group_pass(cs%pass_MEKE, g%Domain)
323  call cpu_clock_end(cs%id_clock_pass)
324  endif
325 
326  if (cs%MEKE_K4 >= 0.0) then
327  ! Calculate Laplacian of MEKE
328 !$OMP parallel default(none) shared(is,ie,js,je,MEKE_uflux,G,MEKE,MEKE_vflux,CS)
329 !$OMP do
330  do j=js,je ; do i=is-1,ie
331  meke_uflux(i,j) = ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * g%mask2dCu(i,j)) * &
332  (meke%MEKE(i+1,j) - meke%MEKE(i,j))
333  ! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
334  ! ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
335  ! (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
336  enddo ; enddo
337 !$OMP do
338  do j=js-1,je ; do i=is,ie
339  meke_vflux(i,j) = ((g%dx_Cv(i,j)*g%IdyCv(i,j)) * g%mask2dCv(i,j)) * &
340  (meke%MEKE(i,j+1) - meke%MEKE(i,j))
341  ! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * &
342  ! ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
343  ! (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
344  enddo ; enddo
345 !$OMP do
346  do j=js,je ; do i=is,ie
347  cs%del2MEKE(i,j) = g%IareaT(i,j) * &
348  ((meke_uflux(i,j) - meke_uflux(i-1,j)) + (meke_vflux(i,j) - meke_vflux(i,j-1)))
349  ! CS%del2MEKE(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * &
350  ! ((MEKE_uflux(I,j) - MEKE_uflux(I-1,j)) + (MEKE_vflux(i,J) - MEKE_vflux(i,J-1)))
351  enddo ; enddo
352 !$OMP end parallel
353  call cpu_clock_begin(cs%id_clock_pass)
354  call do_group_pass(cs%pass_del2MEKE, g%Domain)
355  call cpu_clock_end(cs%id_clock_pass)
356 
357  ! Bi-harmonic diffusion of MEKE
358 !$OMP parallel default(none) shared(is,ie,js,je,MEKE_uflux,G,CS,sdt,mass, &
359 !$OMP mass_neglect,MEKE_vflux,I_mass) &
360 !$OMP private(K4_here,Inv_Kh_max)
361 !$OMP do
362  do j=js,je ; do i=is-1,ie
363  k4_here = cs%MEKE_K4
364  ! Limit Kh to avoid CFL violations.
365  inv_kh_max = 64.0*sdt * (((g%dy_Cu(i,j)*g%IdxCu(i,j)) * &
366  max(g%IareaT(i,j),g%IareaT(i+1,j))))**2.0
367  if (k4_here*inv_kh_max > 0.3) k4_here = 0.3 / inv_kh_max
368 
369  meke_uflux(i,j) = ((k4_here * (g%dy_Cu(i,j)*g%IdxCu(i,j))) * &
370  ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
371  (cs%del2MEKE(i+1,j) - cs%del2MEKE(i,j))
372  enddo ; enddo
373 !$OMP do
374  do j=js-1,je ; do i=is,ie
375  k4_here = cs%MEKE_K4
376  inv_kh_max = 64.0*sdt * (((g%dx_Cv(i,j)*g%IdyCv(i,j)) * &
377  max(g%IareaT(i,j),g%IareaT(i,j+1))))**2.0
378  if (k4_here*inv_kh_max > 0.3) k4_here = 0.3 / inv_kh_max
379 
380  meke_vflux(i,j) = ((k4_here * (g%dx_Cv(i,j)*g%IdyCv(i,j))) * &
381  ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
382  (cs%del2MEKE(i,j+1) - cs%del2MEKE(i,j))
383  enddo ; enddo
384 !$OMP do
385  ! Store tendency of bi-harmonic in del2MEKE
386  do j=js,je ; do i=is,ie
387  cs%del2MEKE(i,j) = (sdt*(g%IareaT(i,j)*i_mass(i,j))) * &
388  ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &
389  (meke_vflux(i,j-1) - meke_vflux(i,j)))
390  enddo ; enddo
391 !$OMP end parallel
392  endif !
393 
394 !$OMP parallel default(none) shared(is,ie,js,je,MEKE,CS,sdt,G,Kh_u,MEKE_uflux, &
395 !$OMP mass,mass_neglect,Kh_v,MEKE_vflux,I_mass, &
396 !$OMP sdt_damp,drag_rate,Rho0,drag_rate_visc, &
397 !$OMP cdrag2,bottomFac2,MEKE_decay,barotrFac2, &
398 !$OMP use_drag_rate,dt,baroHu,baroHv) &
399 !$OMP private(Kh_here,Inv_Kh_max,ldamping,advFac)
400  if (cs%MEKE_KH >= 0.0 .or. cs%KhMEKE_FAC > 0.0 .or. cs%MEKE_advection_factor >0.0) then
401  ! Lateral diffusion of MEKE
402  kh_here = max(0.,cs%MEKE_Kh)
403 !$OMP do
404  do j=js,je ; do i=is-1,ie
405  ! Limit Kh to avoid CFL violations.
406  if (associated(meke%Kh)) &
407  kh_here = max(0.,cs%MEKE_Kh) + cs%KhMEKE_Fac*0.5*(meke%Kh(i,j)+meke%Kh(i+1,j))
408  inv_kh_max = 2.0*sdt * ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * &
409  max(g%IareaT(i,j),g%IareaT(i+1,j)))
410  if (kh_here*inv_kh_max > 0.25) kh_here = 0.25 / inv_kh_max
411  kh_u(i,j) = kh_here
412 
413  meke_uflux(i,j) = ((kh_here * (g%dy_Cu(i,j)*g%IdxCu(i,j))) * &
414  ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
415  (meke%MEKE(i,j) - meke%MEKE(i+1,j))
416  enddo ; enddo
417 !$OMP do
418  do j=js-1,je ; do i=is,ie
419  if (associated(meke%Kh)) &
420  kh_here = max(0.,cs%MEKE_Kh) + cs%KhMEKE_Fac*0.5*(meke%Kh(i,j)+meke%Kh(i,j+1))
421  inv_kh_max = 2.0*sdt * ((g%dx_Cv(i,j)*g%IdyCv(i,j)) * &
422  max(g%IareaT(i,j),g%IareaT(i,j+1)))
423  if (kh_here*inv_kh_max > 0.25) kh_here = 0.25 / inv_kh_max
424  kh_v(i,j) = kh_here
425 
426  meke_vflux(i,j) = ((kh_here * (g%dx_Cv(i,j)*g%IdyCv(i,j))) * &
427  ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
428  (meke%MEKE(i,j) - meke%MEKE(i,j+1))
429  enddo ; enddo
430  if (cs%MEKE_advection_factor>0.) then
431  advfac = cs%MEKE_advection_factor / dt
432 !$OMP do
433  do j=js,je ; do i=is-1,ie
434  if (barohu(i,j)>0.) then
435  meke_uflux(i,j) = meke_uflux(i,j) + barohu(i,j)*meke%MEKE(i,j)*advfac
436  elseif (barohu(i,j)<0.) then
437  meke_uflux(i,j) = meke_uflux(i,j) + barohu(i,j)*meke%MEKE(i+1,j)*advfac
438  endif
439  enddo ; enddo
440 !$OMP do
441  do j=js-1,je ; do i=is,ie
442  if (barohv(i,j)>0.) then
443  meke_vflux(i,j) = meke_vflux(i,j) + barohv(i,j)*meke%MEKE(i,j)*advfac
444  elseif (barohv(i,j)<0.) then
445  meke_vflux(i,j) = meke_vflux(i,j) + barohv(i,j)*meke%MEKE(i,j+1)*advfac
446  endif
447  enddo ; enddo
448  endif
449 !$OMP do
450  do j=js,je ; do i=is,ie
451  meke%MEKE(i,j) = meke%MEKE(i,j) + (sdt*(g%IareaT(i,j)*i_mass(i,j))) * &
452  ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &
453  (meke_vflux(i,j-1) - meke_vflux(i,j)))
454  enddo ; enddo
455  endif ! MEKE_KH>0
456 
457  ! Add on bi-harmonic tendency
458  if (cs%MEKE_K4 >= 0.0) then
459 !$OMP do
460  do j=js,je ; do i=is,ie
461  meke%MEKE(i,j) = meke%MEKE(i,j) + cs%del2MEKE(i,j)
462  enddo ; enddo
463  endif
464 
465  ! Second stage of Strang splitting
466  if (cs%MEKE_KH >= 0.0 .or. cs%MEKE_K4 >= 0.0) then
467  if (sdt>sdt_damp) then
468  ! Recalculate the drag rate, since MEKE has changed.
469  if (use_drag_rate) then
470 !$OMP do
471  do j=js,je ; do i=is,ie
472  drag_rate(i,j) = (rho0 * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 &
473  + cdrag2 * ( max(0.0, 2.0*bottomfac2(i,j)*meke%MEKE(i,j)) + cs%MEKE_Uscale**2 ) )
474  enddo ; enddo
475  endif
476 !$OMP do
477  do j=js,je ; do i=is,ie
478  ldamping = cs%MEKE_damping + drag_rate(i,j) * bottomfac2(i,j)
479  if (meke%MEKE(i,j)<0.) ldamping = 0.
480  ! notice that the above line ensures a damping only if MEKE is positive,
481  ! while leaving MEKE unchanged if it is negative
482  meke%MEKE(i,j) = meke%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
483  meke_decay(i,j) = 0.5 * g%mask2dT(i,j) * (meke_decay(i,j) + ldamping)
484  enddo ; enddo
485  endif
486  endif ! MEKE_KH>=0
487 !$OMP end parallel
488 
489  call cpu_clock_begin(cs%id_clock_pass)
490  call do_group_pass(cs%pass_MEKE, g%Domain)
491  call cpu_clock_end(cs%id_clock_pass)
492 
493  ! Calculate diffusivity for main model to use
494  if (cs%MEKE_KhCoeff>0.) then
495  if (cs%use_old_lscale) then
496  if (cs%Rd_as_max_scale) then
497 !$OMP parallel do default(none) shared(is,ie,js,je,MEKE,CS,G,barotrFac2)
498  do j=js,je ; do i=is,ie
499  meke%Kh(i,j) = (cs%MEKE_KhCoeff &
500  * sqrt(2.*max(0.,barotrfac2(i,j)*meke%MEKE(i,j))*g%areaT(i,j))) &
501  * min(meke%Rd_dx_h(i,j), 1.0)
502  enddo ; enddo
503  else
504 !$OMP parallel do default(none) shared(is,ie,js,je,MEKE,CS,G,barotrFac2)
505  do j=js,je ; do i=is,ie
506  meke%Kh(i,j) = cs%MEKE_KhCoeff*sqrt(2.*max(0.,barotrfac2(i,j)*meke%MEKE(i,j))*g%areaT(i,j))
507  enddo ; enddo
508  endif
509  else
510 !$OMP parallel do default(none) shared(is,ie,js,je,MEKE,LmixScale,CS,G,barotrFac2)
511  do j=js,je ; do i=is,ie
512  meke%Kh(i,j) = (cs%MEKE_KhCoeff*sqrt(2.*max(0.,barotrfac2(i,j)*meke%MEKE(i,j)))*lmixscale(i,j))
513  enddo ; enddo
514  endif
515  call cpu_clock_begin(cs%id_clock_pass)
516  call do_group_pass(cs%pass_Kh, g%Domain)
517  call cpu_clock_end(cs%id_clock_pass)
518  endif
519 
520  ! Calculate viscosity for the main model to use
521  if (cs%viscosity_coeff/=0.) then
522 !aja: should make range jsq:jeq, isq:ieq
523  do j=js-1,je+1 ; do i=is-1,ie+1
524  meke%Ku(i,j) = cs%viscosity_coeff*sqrt(2.*max(0.,meke%MEKE(i,j)))*lmixscale(i,j)
525  enddo ; enddo
526  call cpu_clock_begin(cs%id_clock_pass)
527  call do_group_pass(cs%pass_Ku, g%Domain)
528  call cpu_clock_end(cs%id_clock_pass)
529  endif
530 
531 ! Offer fields for averaging.
532  if (cs%id_MEKE>0) call post_data(cs%id_MEKE, meke%MEKE, cs%diag)
533  if (cs%id_Ue>0) call post_data(cs%id_Ue, sqrt(max(0.,2.0*meke%MEKE)), cs%diag)
534  if (cs%id_Ub>0) call post_data(cs%id_Ub, sqrt(max(0.,2.0*meke%MEKE*bottomfac2)), cs%diag)
535  if (cs%id_Ut>0) call post_data(cs%id_Ut, sqrt(max(0.,2.0*meke%MEKE*barotrfac2)), cs%diag)
536  if (cs%id_Kh>0) call post_data(cs%id_Kh, meke%Kh, cs%diag)
537  if (cs%id_Ku>0) call post_data(cs%id_Ku, meke%Ku, cs%diag)
538  if (cs%id_KhMEKE_u>0) call post_data(cs%id_KhMEKE_u, kh_u, cs%diag)
539  if (cs%id_KhMEKE_v>0) call post_data(cs%id_KhMEKE_v, kh_v, cs%diag)
540  if (cs%id_src>0) call post_data(cs%id_src, src, cs%diag)
541  if (cs%id_decay>0) call post_data(cs%id_decay, meke_decay, cs%diag)
542  if (cs%id_GM_src>0) call post_data(cs%id_GM_src, meke%GM_src, cs%diag)
543  if (cs%id_mom_src>0) call post_data(cs%id_mom_src, meke%mom_src, cs%diag)
544  if (cs%id_Le>0) call post_data(cs%id_Le, lmixscale, cs%diag)
545  if (cs%id_gamma_b>0) then
546  do j=js,je ; do i=is,ie
547  bottomfac2(i,j) = sqrt(bottomfac2(i,j))
548  enddo ; enddo
549  call post_data(cs%id_gamma_b, bottomfac2, cs%diag)
550  endif
551  if (cs%id_gamma_t>0) then
552  do j=js,je ; do i=is,ie
553  barotrfac2(i,j) = sqrt(barotrfac2(i,j))
554  enddo ; enddo
555  call post_data(cs%id_gamma_t, barotrfac2, cs%diag)
556  endif
557 
558 ! else ! if MEKE%MEKE
559 ! call MOM_error(FATAL, "MOM_MEKE: MEKE%MEKE is not associated!")
560  endif
561 
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function: