← back to index

src/parameterizations/lateral/MOM_MEKE.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Implements the Mesoscale Eddy Kinetic Energy framework
6!! with topographic beta effect included in computing beta in Rhines scale
7
8module MOM_MEKE
9
10use iso_fortran_env, only : real32
11
12use MOM_coms, only : PE_here
13use MOM_database_comms, only : dbclient_type, dbcomms_CS_type
14use MOM_debugging, only : hchksum, uvchksum
15use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
16use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
17use MOM_diag_mediator, only : diag_ctrl, time_type
18use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
19use MOM_domains, only : pass_vector, pass_var
20use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, MOM_mesg, is_root_pe
21use MOM_file_parser, only : read_param, get_param, log_version, param_file_type
22use MOM_grid, only : ocean_grid_type
23use MOM_hor_index, only : hor_index_type
24use MOM_interface_heights, only : find_eta
25use MOM_interpolate, only : init_external_field, time_interp_external
26use MOM_interpolate, only : time_interp_external_init
27use MOM_interpolate, only : external_field
28use MOM_io, only : vardesc, var_desc, slasher
29use MOM_isopycnal_slopes, only : calc_isoneutral_slopes
30use MOM_restart, only : MOM_restart_CS, register_restart_field, query_initialized
31use MOM_string_functions, only : lowercase
32use MOM_time_manager, only : time_type_to_real
33use MOM_unit_scaling, only : unit_scale_type
34use MOM_variables, only : vertvisc_type, thermo_var_ptrs
35use MOM_verticalGrid, only : verticalGrid_type
36use MOM_MEKE_types, only : MEKE_type
37
38implicit none ; private
39
40#include <MOM_memory.h>
41
42public step_forward_MEKE, MEKE_init, MEKE_alloc_register_restart, MEKE_end
43
44! Constants for this module
45integer, parameter :: NUM_FEATURES = 4 !< How many features used to predict EKE
46integer, parameter :: MKE_IDX = 1 !< Index of mean kinetic energy in the feature array
47integer, parameter :: SLOPE_Z_IDX = 2 !< Index of vertically averaged isopycnal slope in the feature array
48integer, parameter :: RV_IDX = 3 !< Index of surface relative vorticity in the feature array
49integer, parameter :: RD_DX_Z_IDX = 4 !< Index of the radius of deformation over the grid size in the feature array
50
51integer, parameter :: EKE_PROG = 1 !< Use prognostic equation to calculate EKE
52integer, parameter :: EKE_FILE = 2 !< Read in EKE from a file
53integer, parameter :: EKE_DBCLIENT = 3 !< Infer EKE using a neural network
54
55!> Control structure that contains MEKE parameters and diagnostics handles
56type, public :: MEKE_CS ; private
57 logical :: initialized = .false. !< True if this control structure has been initialized.
58 ! Parameters
59 real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim]
60 real :: MEKE_bhFrCoeff!< Efficiency of conversion of ME into MEKE by the biharmonic dissipation [nondim]
61 real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim]
62 real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME [nondim]
63 real :: MEKE_damping !< Local depth-independent MEKE dissipation rate [T-1 ~> s-1].
64 real :: MEKE_Cd_scale !< The ratio of the bottom eddy velocity to the column mean
65 !! eddy velocity, i.e. sqrt(2*MEKE), [nondim]. This should be less than 1
66 !! to account for the surface intensification of MEKE.
67 real :: MEKE_Cb !< Coefficient in the \f$\gamma_{bot}\f$ expression [nondim]
68 real :: MEKE_min_gamma!< Minimum value of gamma_b^2 allowed [nondim]
69 real :: MEKE_Ct !< Coefficient in the \f$\gamma_{bt}\f$ expression [nondim]
70 logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag.
71 logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC
72 !! framework (Marshall et al., 2012)
73 real :: MEKE_GEOMETRIC_alpha !< The nondimensional coefficient governing the efficiency of the
74 !! GEOMETRIC thickness diffusion [nondim].
75 logical :: MEKE_equilibrium_alt !< If true, use an alternative calculation for the
76 !! equilibrium value of MEKE.
77 logical :: MEKE_equilibrium_restoring !< If true, restore MEKE back to its equilibrium value,
78 !! which is calculated at each time step.
79 logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather
80 !! than the streamfunction for the MEKE GM source term.
81 real :: MEKE_min_depth_tot !< The minimum total thickness over which to distribute MEKE energy
82 !! sources from GM energy conversion [H ~> m or kg m-2]. When the total
83 !! thickness is less than this, the sources are scaled away.
84 logical :: Rd_as_max_scale !< If true the length scale can not exceed the
85 !! first baroclinic deformation radius.
86 logical :: use_old_lscale !< Use the old formula for mixing length scale.
87 logical :: use_min_lscale !< Use simple minimum for mixing length scale.
88 logical :: MEKE_positive !< If true, it guarantees that MEKE will always be >= 0.
89 real :: lscale_maxval !< The ceiling on the MEKE mixing length scale when use_min_lscale is true [L ~> m].
90 real :: cdrag !< The bottom drag coefficient for MEKE, times rescaling factors [H L-1 ~> nondim or kg m-3]
91 real :: MEKE_BGsrc !< Background energy source for MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
92 real :: MEKE_dtScale !< Scale factor to accelerate time-stepping [nondim]
93 real :: MEKE_KhCoeff !< Scaling factor to convert MEKE into Kh [nondim]
94 real :: MEKE_Uscale !< MEKE velocity scale for bottom drag [L T-1 ~> m s-1]
95 real :: MEKE_KH !< Background lateral diffusion of MEKE [L2 T-1 ~> m2 s-1]
96 real :: MEKE_K4 !< Background bi-harmonic diffusivity (of MEKE) [L4 T-1 ~> m4 s-1]
97 real :: KhMEKE_Fac !< A factor relating MEKE%Kh to the diffusivity used for
98 !! MEKE itself [nondim].
99 real :: viscosity_coeff_Ku !< The scaling coefficient in the expression for
100 !! viscosity used to parameterize lateral harmonic momentum mixing
101 !! by unresolved eddies represented by MEKE [nondim].
102 real :: viscosity_coeff_Au !< The scaling coefficient in the expression for
103 !! viscosity used to parameterize lateral biharmonic momentum mixing
104 !! by unresolved eddies represented by MEKE [nondim].
105 real :: Lfixed !< Fixed mixing length scale [L ~> m].
106 real :: aDeform !< Weighting towards deformation scale of mixing length [nondim]
107 real :: aRhines !< Weighting towards Rhines scale of mixing length [nondim]
108 real :: aFrict !< Weighting towards frictional arrest scale of mixing length [nondim]
109 real :: aEady !< Weighting towards Eady scale of mixing length [nondim]
110 real :: aGrid !< Weighting towards grid scale of mixing length [nondim]
111 real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim]
112 real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered
113 !! when computing beta in Rhines scale [nondim]
114 real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its
115 !! equilibrium value [T-1 ~> s-1].
116 logical :: MEKE_advection_bug !< If true, recover a bug in the calculation of the barotropic
117 !! transport for the advection of MEKE, wherein only the transports in the
118 !! deepest layer are used.
119 logical :: fixed_total_depth !< If true, use the nominal bathymetric depth as the estimate of
120 !! the time-varying ocean depth. Otherwise base the depth on the total
121 !! ocean mass per unit area.
122 real :: rho_fixed_total_depth !< A density used to translate the nominal bathymetric depth into an
123 !! estimate of the total ocean mass per unit area when MEKE_FIXED_TOTAL_DEPTH
124 !! is true [R ~> kg m-3]
125 logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled.
126 logical :: initialize !< If True, invokes a steady state solver to calculate MEKE.
127 logical :: debug !< If true, write out checksums of data for debugging
128 integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default),
129 !! read in from a file, or inferred via a neural network
130 logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure.
131 type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output
132 !>@{ Diagnostic handles
133 integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
134 integer :: id_src_adv = -1, id_src_mom_K4 = -1, id_src_btm_drag = -1
135 integer :: id_src_GM = -1, id_src_mom_lp = -1, id_src_mom_bh = -1
136 integer :: id_Ub = -1, id_Ut = -1
137 integer :: id_GM_src = -1, id_mom_src = -1, id_mom_src_bh = -1, id_GME_snk = -1, id_decay = -1
138 integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1
139 integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1
140 integer :: id_Lrhines = -1, id_Leady = -1
141 integer :: id_MEKE_equilibrium = -1
142 !>@}
143 type(external_field) :: eke_handle !< Handle for reading in EKE from a file
144 ! Infrastructure
145 integer :: id_clock_pass !< Clock for group pass calls
146 type(group_pass_type) :: pass_MEKE !< Group halo pass handle for MEKE%MEKE and maybe MEKE%Kh_diff
147 type(group_pass_type) :: pass_Kh !< Group halo pass handle for MEKE%Kh, MEKE%Ku, and/or MEKE%Au
148
149 ! MEKE via Machine Learning
150 type(dbclient_type), pointer :: client => NULL() !< Pointer to the database client
151
152 logical :: online_analysis !< If true, post the EKE used in MOM6 at every timestep
153 character(len=5) :: model_key = 'mleke' !< Key where the ML-model is stored
154 character(len=7) :: key_suffix !< Suffix appended to every key sent to Redis
155 real :: eke_max !< The maximum value of EKE considered physically reasonable [L2 T-2 ~> m2 s-2]
156
157 ! Clock ids
158 integer :: id_client_init !< Clock id to time initialization of the client
159 integer :: id_put_tensor !< Clock id to time put_tensor routine
160 integer :: id_run_model !< Clock id to time running of the ML model
161 integer :: id_unpack_tensor !< Clock id to time retrieval of EKE prediction
162
163 ! Diagnostic ids
164 integer :: id_mke = -1 !< Diagnostic id for surface mean kinetic energy
165 integer :: id_slope_z = -1 !< Diagnostic id for vertically averaged horizontal slope magnitude
166 integer :: id_slope_x = -1 !< Diagnostic id for isopycnal slope in the x-direction
167 integer :: id_slope_y = -1 !< Diagnostic id for isopycnal slope in the y-direction
168 integer :: id_rv = -1 !< Diagnostic id for surface relative vorticity
169
170 ! Isoneutral blocking parameters
171 integer :: niblock !< The i block size used in calc_isoneutral_slopes [nondim].
172 integer :: njblock !< The j block size used in calc_isoneutral_slopes [nondim].
173 integer :: nkblock !< The k block size used in calc_isoneutral_slopes [nondim].
174
175end type MEKE_CS
176
177contains
178
179!> Integrates forward-in-time the MEKE eddy energy equation.
180!! See \ref section_MEKE_equations.
1810subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, hv, u, v, tv, Time)
182 type(MEKE_type), intent(inout) :: MEKE !< MEKE data.
183 type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
184 type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
185 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
186 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
187 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
188 real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
189 type(vertvisc_type), intent(in) :: visc !< The vertical viscosity type.
190 real, intent(in) :: dt !< Model(baroclinic) time-step [T ~> s].
191 type(MEKE_CS), intent(inout) :: CS !< MEKE control structure.
192 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: hu !< Accumulated zonal mass flux [H L2 ~> m3 or kg].
193 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: hv !< Accumulated meridional mass flux [H L2 ~> m3 or kg]
194 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
195 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
196 type(thermo_var_ptrs), intent(in) :: tv !< Type containing thermodynamic variables
197 type(time_type), intent(in) :: Time !< The time used for interpolating EKE
198
199 ! Local variables
200 real, dimension(SZI_(G),SZJ_(G)) :: &
2010 data_eke, & ! EKE from file [L2 T-2 ~> m2 s-2]
2020 mass, & ! The total mass of the water column [R Z ~> kg m-2].
2030 I_mass, & ! The inverse of mass [R-1 Z-1 ~> m2 kg-1].
2040 depth_tot, & ! The depth of the water column [H ~> m or kg m-2].
2050 src, & ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3).
2060 MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1].
2070 src_adv, & ! The MEKE source/tendency from the horizontal advection of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
2080 src_mom_K4, & ! The MEKE source/tendency from the bihamornic of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
2090 src_btm_drag, & ! The MEKE source/tendency from the bottom drag acting on MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
2100 src_GM, & ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3).
2110 src_mom_lp, & ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
2120 src_mom_bh, & ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
2130 damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].
2140 MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2].
2150 drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1]
2160 drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
2170 del2MEKE, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2].
2180 del4MEKE, & ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].
2190 LmixScale, & ! Eddy mixing length [L ~> m].
2200 barotrFac2, & ! Ratio of EKE_barotropic / EKE [nondim]
2210 bottomFac2, & ! Ratio of EKE_bottom / EKE [nondim]
2220 tmp, & ! Temporary variable for computation of diagnostic velocities [L T-1 ~> m s-1]
2230 equilibrium_value, & ! The equilibrium value of MEKE to be calculated at
224 ! each time step [L2 T-2 ~> m2 s-2]
2250 damp_rate, & ! The MEKE damping rate [T-1 ~> s-1]
2260 damping ! The net damping of a field after sdt_damp [nondim]
227
228 real, dimension(SZIB_(G),SZJ_(G)) :: &
2290 MEKE_uflux, & ! The zonal advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3].
230 ! In one place, MEKE_uflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].
2310 Kh_u, & ! The zonal diffusivity that is actually used [L2 T-1 ~> m2 s-1].
2320 baroHu, & ! Depth integrated accumulated zonal mass flux [R Z L2 ~> kg].
2330 drag_vel_u ! A piston velocity associated with bottom drag at u-points [H T-1 ~> m s-1 or kg m-2 s-1]
234 real, dimension(SZI_(G),SZJB_(G)) :: &
2350 MEKE_vflux, & ! The meridional advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3].
236 ! In one place, MEKE_vflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].
2370 Kh_v, & ! The meridional diffusivity that is actually used [L2 T-1 ~> m2 s-1].
2380 baroHv, & ! Depth integrated accumulated meridional mass flux [R Z L2 ~> kg].
2390 drag_vel_v ! A piston velocity associated with bottom drag at v-points [H T-1 ~> m s-1 or kg m-2 s-1]
240 real :: bh_coeff ! Biharmonic part of efficiency conversion in total MEKE [nondim]
241 real :: Kh_here ! The local horizontal viscosity [L2 T-1 ~> m2 s-1]
242 real :: Inv_Kh_max ! The inverse of the local horizontal viscosity [T L-2 ~> s m-2]
243 real :: K4_here ! The local horizontal biharmonic viscosity [L4 T-1 ~> m4 s-1]
244 real :: Inv_K4_max ! The inverse of the local horizontal biharmonic viscosity [T L-4 ~> s m-4]
245 real :: cdrag2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]
246 real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1]
247 real :: mass_neglect ! A negligible mass [R Z ~> kg m-2].
248 real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate)
249 real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split).
250 real :: damp_step ! Size of damping timestep relative to sdt [nondim]
251 logical :: use_drag_rate ! Flag to indicate drag_rate is finite
252 logical :: any_damping_diags_s1 ! True if any damped diagnostics are enabled in first stage
253 logical :: any_damping_diags ! True if any damped diagnostics are enabled
254 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
2550 real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features
256 ! needed for the machine learning inference, with different
257 ! units for the various subarrays [various]
258
2590 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2600 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
261
2620 if (.not.CS%initialized) call MOM_error(FATAL, &
2630 "MOM_MEKE: Module must be initialized before it is used.")
264
2650 if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then
2660 use_drag_rate = .true.
267 else
2680 use_drag_rate = .false.
269 endif
270
271 ! Only integrate the MEKE equations if MEKE is required.
2720 if (.not. allocated(MEKE%MEKE)) then
273! call MOM_error(FATAL, "MOM_MEKE: MEKE%MEKE is not associated!")
2740 return
275 endif
276
2770 select case(CS%eke_src)
278 case(EKE_PROG)
2790 if (CS%debug) then
2800 if (allocated(MEKE%mom_src)) &
2810 call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
2820 if (allocated(MEKE%mom_src_bh)) &
2830 call hchksum(MEKE%mom_src_bh, 'MEKE mom_src_bh', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
2840 if (allocated(MEKE%GME_snk)) &
2850 call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
2860 if (allocated(MEKE%GM_src)) &
2870 call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
2880 if (allocated(MEKE%MEKE)) &
2890 call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, unscale=US%L_T_to_m_s**2)
290 call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, unscale=US%s_to_T, &
2910 scalar_pair=.true.)
292 call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=0, symmetric=.true., &
2930 unscale=GV%H_to_m*US%L_to_m**2)
294 endif
295
2960 sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping
2970 mass_neglect = GV%H_to_RZ * GV%H_subroundoff
2980 cdrag2 = CS%cdrag**2
299
300 ! With a depth-dependent (and possibly strong) damping, it seems
301 ! advisable to use Strang splitting between the damping and diffusion.
3020 damp_step = 1.
3030 if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) damp_step = 0.5
3040 sdt_damp = sdt * damp_step
305
306 ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg]
3070 if (CS%MEKE_advection_factor>0.) then
3080 do j=js,je ; do I=is-1,ie
3090 baroHu(I,j) = 0.
310 enddo ; enddo
3110 do k=1,nz
3120 do j=js,je ; do I=is-1,ie
3130 baroHu(I,j) = baroHu(I,j) + hu(I,j,k) * GV%H_to_RZ
314 enddo ; enddo
315 enddo
3160 do J=js-1,je ; do i=is,ie
3170 baroHv(i,J) = 0.
318 enddo ; enddo
3190 do k=1,nz
3200 do J=js-1,je ; do i=is,ie
3210 baroHv(i,J) = baroHv(i,J) + hv(i,J,k) * GV%H_to_RZ
322 enddo ; enddo
323 enddo
3240 if (CS%MEKE_advection_bug) then
325 ! This obviously incorrect code reproduces a bug in the original implementation of
326 ! the MEKE advection.
3270 do j=js,je ; do I=is-1,ie
3280 baroHu(I,j) = hu(I,j,nz) * GV%H_to_RZ
329 enddo ; enddo
3300 do J=js-1,je ; do i=is,ie
3310 baroHv(i,J) = hv(i,J,nz) * GV%H_to_RZ
332 enddo ; enddo
333 endif
334 endif
335
336 ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow
3370 if (CS%visc_drag .and. allocated(visc%Kv_bbl_u) .and. allocated(visc%Kv_bbl_v)) then
338 !$OMP parallel do default(shared)
3390 do j=js,je ; do I=is-1,ie
3400 drag_vel_u(I,j) = 0.0
3410 if ((G%mask2dCu(I,j) > 0.0) .and. (visc%bbl_thick_u(I,j) > 0.0)) &
3420 drag_vel_u(I,j) = visc%Kv_bbl_u(I,j) / visc%bbl_thick_u(I,j)
343 enddo ; enddo
344 !$OMP parallel do default(shared)
3450 do J=js-1,je ; do i=is,ie
3460 drag_vel_v(i,J) = 0.0
3470 if ((G%mask2dCv(i,J) > 0.0) .and. (visc%bbl_thick_v(i,J) > 0.0)) &
3480 drag_vel_v(i,J) = visc%Kv_bbl_v(i,J) / visc%bbl_thick_v(i,J)
349 enddo ; enddo
350
351 !$OMP parallel do default(shared)
3520 do j=js,je ; do i=is,ie
353 drag_rate_visc(i,j) = (0.25*G%IareaT(i,j) * &
354 (((G%areaCu(I-1,j)*drag_vel_u(I-1,j)) + &
355 (G%areaCu(I,j)*drag_vel_u(I,j))) + &
356 ((G%areaCv(i,J-1)*drag_vel_v(i,J-1)) + &
3570 (G%areaCv(i,J)*drag_vel_v(i,J))) ) )
358 enddo ; enddo
359 else
360 !$OMP parallel do default(shared)
3610 do j=js,je ; do i=is,ie
3620 drag_rate_visc(i,j) = 0.
363 enddo ; enddo
364 endif
365
366 !$OMP parallel do default(shared)
3670 do j=js-1,je+1
3680 do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo
3690 do k=1,nz ; do i=is-1,ie+1
3700 mass(i,j) = mass(i,j) + G%mask2dT(i,j) * (GV%H_to_RZ * h(i,j,k)) ! [R Z ~> kg m-2]
371 enddo ; enddo
3720 do i=is-1,ie+1
3730 I_mass(i,j) = 0.0
3740 if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) ! [R-1 Z-1 ~> m2 kg-1]
375 enddo
376 enddo
377
3780 if (CS%fixed_total_depth) then
3790 if (GV%Boussinesq) then
380 !$OMP parallel do default(shared)
3810 do j=js-1,je+1 ; do i=is-1,ie+1
3820 depth_tot(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * GV%Z_to_H
383 enddo ; enddo
384 else
385 !$OMP parallel do default(shared)
3860 do j=js-1,je+1 ; do i=is-1,ie+1
3870 depth_tot(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * CS%rho_fixed_total_depth * GV%RZ_to_H
388 enddo ; enddo
389 endif
390 else
391 !$OMP parallel do default(shared)
3920 do j=js-1,je+1 ; do i=is-1,ie+1
3930 depth_tot(i,j) = mass(i,j) * GV%RZ_to_H
394 enddo ; enddo
395 endif
396
3970 if (CS%initialize) then
3980 call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot)
3990 CS%initialize = .false.
400 endif
401
402 ! Calculates bottomFac2, barotrFac2 and LmixScale
4030 call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale)
4040 if (CS%debug) then
4050 if (CS%visc_drag) &
406 call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, &
4070 unscale=GV%H_to_mks*US%s_to_T, scalar_pair=.true.)
4080 call hchksum(mass, 'MEKE mass',G%HI,haloshift=1, unscale=US%RZ_to_kg_m2)
4090 call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', G%HI, unscale=GV%H_to_mks*US%s_to_T)
4100 call hchksum(bottomFac2, 'MEKE bottomFac2', G%HI)
4110 call hchksum(barotrFac2, 'MEKE barotrFac2', G%HI)
4120 call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m)
413 endif
414
4150 if (allocated(MEKE%Le)) then
416 !$OMP parallel do default(shared)
4170 do j=js,je ; do i=is,ie
4180 MEKE%Le(i,j) = LmixScale(i,j)
419 enddo ; enddo
420 endif
421
422 ! Aggregate sources of MEKE (background, frictional and GM)
423 !$OMP parallel do default(shared)
4240 do j=js,je ; do i=is,ie
4250 src(i,j) = CS%MEKE_BGsrc
426 enddo ; enddo
427
428 ! Initialize diagnostics
4290 if (CS%id_src_adv > 0) src_adv(is:ie, js:je) = 0.
4300 if (CS%id_src_GM > 0) src_GM(is:ie, js:je) = 0.
4310 if (CS%id_src_mom_lp > 0) src_mom_lp(is:ie, js:je) = 0.
4320 if (CS%id_src_mom_bh > 0) src_mom_bh(is:ie, js:je) = 0.
4330 if (CS%id_src_mom_K4 > 0) src_mom_K4(is:ie, js:je) = 0.
4340 if (CS%id_src_btm_drag > 0) src_btm_drag(is:ie, js:je) = 0.
435
436 ! Identify any damped diagnostics in first stage of Strang splitting
437 any_damping_diags_s1 = any([ &
438 CS%id_src_GM > 0, &
439 CS%id_src_mom_lp > 0, &
440 CS%id_src_mom_bh > 0, &
441 CS%id_src_btm_drag > 0 &
4420 ])
443
444 ! Identify any damped diagnostics
445 any_damping_diags = any([ &
446 any_damping_diags_s1, &
447 CS%id_src_adv > 0, &
448 CS%id_src_mom_K4 > 0 &
4490 ])
450
4510 if (CS%MEKE_FrCoeff > 0.) then
452 !$OMP parallel do default(shared)
4530 do j=js,je ; do i=is,ie
4540 src(i,j) = src(i,j) - CS%MEKE_FrCoeff * I_mass(i,j) * MEKE%mom_src(i,j)
455 enddo ; enddo
456 endif
457
4580 if (allocated(MEKE%mom_src_bh)) then
4590 if (CS%MEKE_bhFrCoeff > 0. .and. CS%MEKE_FrCoeff > 0.) then
4600 bh_coeff = CS%MEKE_bhFrCoeff - CS%MEKE_FrCoeff
461 else
4620 bh_coeff = CS%MEKE_bhFrCoeff
463 endif
464
465 !$OMP parallel do default(shared)
4660 do j=js,je ; do i=is,ie
4670 src(i,j) = src(i,j) - bh_coeff * I_mass(i,j) * MEKE%mom_src_bh(i,j)
468 enddo ; enddo
469
4700 if (CS%id_src_mom_lp > 0) then
471 !$OMP parallel do default(shared)
4720 do j=js,je ; do i=is,ie
473 src_mom_lp(i,j) = -CS%MEKE_FrCoeff * I_mass(i,j) &
4740 * (MEKE%mom_src(i,j) - MEKE%mom_src_bh(i,j))
475 enddo ; enddo
476 endif
477
4780 if (CS%id_src_mom_bh > 0) then
479 !$OMP parallel do default(shared)
4800 do j=js,je ; do i=is,ie
4810 src_mom_bh(i,j) = -CS%MEKE_bhFrCoeff * I_mass(i,j) * MEKE%mom_src_bh(i,j)
482 enddo ; enddo
483 endif
484 endif
485
4860 if (allocated(MEKE%GME_snk)) then
487 !$OMP parallel do default(shared)
4880 do j=js,je ; do i=is,ie
4890 src(i,j) = src(i,j) - CS%MEKE_GMECoeff*I_mass(i,j)*MEKE%GME_snk(i,j)
490 enddo ; enddo
491 endif
492
4930 if (allocated(MEKE%GM_src)) then
4940 if (CS%GM_src_alt) then
495 !$OMP parallel do default(shared)
4960 do j=js,je ; do i=is,ie
497 src(i,j) = src(i,j) - CS%MEKE_GMcoeff*MEKE%GM_src(i,j) / &
4980 (GV%H_to_RZ * MAX(CS%MEKE_min_depth_tot, depth_tot(i,j)))
499 enddo ; enddo
500 else
501 !$OMP parallel do default(shared)
5020 do j=js,je ; do i=is,ie
5030 src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j)
504 enddo ; enddo
505
5060 do j=js,je ; do i=is,ie
5070 src_GM(i,j) = -CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j)
508 enddo ; enddo
509 endif
510 endif
511
5120 if (CS%MEKE_equilibrium_restoring) then
513 call MEKE_equilibrium_restoring(CS, G, GV, US, SN_u, SN_v, depth_tot, &
5140 equilibrium_value)
5150 do j=js,je ; do i=is,ie
5160 src(i,j) = src(i,j) - CS%MEKE_restoring_rate*(MEKE%MEKE(i,j) - equilibrium_value(i,j))
517 enddo ; enddo
518 endif
519
5200 if (CS%debug) then
5210 call hchksum(src, "MEKE src", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T**3)
522 endif
523
524 ! Increase EKE by a full time-steps worth of source
525 !$OMP parallel do default(shared)
5260 do j=js,je ; do i=is,ie
5270 MEKE_current(i,j) = MEKE%MEKE(i,j)
5280 MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j)
529 enddo ; enddo
530
5310 if (use_drag_rate) then
532 ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
533 !$OMP parallel do default(shared)
5340 do j=js,je ; do i=is,ie
535 drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
5360 cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
537 enddo ; enddo
538 else
539 !$OMP parallel do default(shared)
5400 do j=js,je ; do i=is,ie
5410 drag_rate(i,j) = 0.
542 enddo ; enddo
543 endif
544
545 ! First stage of Strang splitting
546
547 !$OMP parallel do default(shared)
5480 do j=js,je ; do i=is,ie
5490 damp_rate(i,j) = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
550
5510 if (MEKE%MEKE(i,j) < 0.) damp_rate(i,j) = 0.
552 ! notice that the above line ensures a damping only if MEKE is positive,
553 ! while leaving MEKE unchanged if it is negative
554 enddo ; enddo
555
556 ! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the existing
557 ! bit-reproducible solution.
558 !$OMP parallel do default(shared)
5590 do j=js,je ; do i=is,ie
5600 MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j))
561 enddo ; enddo
562
5630 if (any_damping_diags_s1) then
564 !$OMP parallel do default(shared)
5650 do j=js,je ; do i=is,ie
5660 damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j))
567 enddo ; enddo
568
5690 if (CS%id_decay > 0) then
570 !$OMP parallel do default(shared)
5710 do j=js,je ; do i=is,ie
5720 MEKE_decay(i,j) = damp_rate(i,j) * G%mask2dT(i,j)
573 enddo ; enddo
574 endif
575
5760 if (CS%id_src_GM > 0) then
577 !$OMP parallel do default(shared)
5780 do j=js,je ; do i=is,ie
5790 src_GM(i,j) = src_GM(i,j) * damping(i,j)
580 enddo ; enddo
581 endif
582
5830 if (CS%id_src_mom_lp > 0) then
584 !$OMP parallel do default(shared)
5850 do j=js,je ; do i=is,ie
5860 src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j)
587 enddo ; enddo
588 endif
589
5900 if (CS%id_src_mom_bh > 0) then
591 !$OMP parallel do default(shared)
5920 do j=js,je ; do i=is,ie
5930 src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j)
594 enddo ; enddo
595 endif
596
5970 if (CS%id_src_btm_drag > 0) then
598 !$OMP parallel do default(shared)
5990 do j=js,je ; do i=is,ie
600 src_btm_drag(i,j) = -MEKE_current(i,j) * ( &
601 damp_step * (damp_rate(i,j) * damping(i,j)) &
6020 )
603 enddo ; enddo
604
605 ! Store the effective damping rate if sdt is split
6060 if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) then
607 !$OMP parallel do default(shared)
6080 do j=js,je ; do i=is,ie
6090 damp_rate_s1(i,j) = damp_rate(i,j) * damping(i,j)
610 enddo ; enddo
611 endif
612 endif
613 endif
614
6150 if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then
616 ! Update MEKE in the halos for lateral or bi-harmonic diffusion
6170 call cpu_clock_begin(CS%id_clock_pass)
6180 call do_group_pass(CS%pass_MEKE, G%Domain)
6190 call cpu_clock_end(CS%id_clock_pass)
620 endif
621
6220 if (CS%MEKE_K4 >= 0.0) then
623 ! Calculate Laplacian of MEKE using MEKE_uflux and MEKE_vflux as temporary work space.
624 !$OMP parallel do default(shared)
6250 do j=js-1,je+1 ; do I=is-2,ie+1
626 ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
627 MEKE_uflux(I,j) = (G%dy_Cu(I,j)*G%IdxCu_OBCmask(I,j)) * &
6280 (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
629 ! This would have units of [R Z L2 T-2 ~> kg s-2]
630 ! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
631 ! ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
632 ! (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
633 enddo ; enddo
634 !$OMP parallel do default(shared)
6350 do J=js-2,je+1 ; do i=is-1,ie+1
636 ! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
637 MEKE_vflux(i,J) = (G%dx_Cv(i,J)*G%IdyCv_OBCmask(i,J)) * &
6380 (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
639 ! This would have units of [R Z L2 T-2 ~> kg s-2]
640 ! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * &
641 ! ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
642 ! (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
643 enddo ; enddo
644
645 !$OMP parallel do default(shared)
6460 do j=js-1,je+1 ; do i=is-1,ie+1 ! del2MEKE has units [T-2 ~> s-2].
647 del2MEKE(i,j) = G%IareaT(i,j) * &
6480 ((MEKE_uflux(I,j) - MEKE_uflux(I-1,j)) + (MEKE_vflux(i,J) - MEKE_vflux(i,J-1)))
649 enddo ; enddo
650
651 ! Bi-harmonic diffusion of MEKE
652 !$OMP parallel do default(shared) private(K4_here,Inv_K4_max)
6530 do j=js,je ; do I=is-1,ie
6540 K4_here = CS%MEKE_K4 ! [L4 T-1 ~> m4 s-1]
655 ! Limit Kh to avoid CFL violations.
656 Inv_K4_max = 64.0 * sdt * ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
6570 max(G%IareaT(i,j), G%IareaT(i+1,j)))**2
6580 if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max
659
660 ! Here the units of MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].
661 MEKE_uflux(I,j) = ((K4_here * (G%dy_Cu(I,j)*G%IdxCu(I,j))) * &
662 ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
6630 (del2MEKE(i+1,j) - del2MEKE(i,j))
664 enddo ; enddo
665 !$OMP parallel do default(shared) private(K4_here,Inv_K4_max)
6660 do J=js-1,je ; do i=is,ie
6670 K4_here = CS%MEKE_K4 ! [L4 T-1 ~> m4 s-1]
6680 Inv_K4_max = 64.0 * sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * max(G%IareaT(i,j), G%IareaT(i,j+1)))**2
6690 if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max
670
671 ! Here the units of MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
672 MEKE_vflux(i,J) = ((K4_here * (G%dx_Cv(i,J)*G%IdyCv(i,J))) * &
673 ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
6740 (del2MEKE(i,j+1) - del2MEKE(i,j))
675 enddo ; enddo
676 ! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2 ~> m2 s-2].
677 !$OMP parallel do default(shared)
6780 do j=js,je ; do i=is,ie
679 del4MEKE(i,j) = (sdt*(G%IareaT(i,j)*I_mass(i,j))) * &
680 ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
6810 (MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
682 src_mom_K4(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * &
683 ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
6840 (MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
685 enddo ; enddo
686 endif !
687
6880 if (CS%kh_flux_enabled) then
689 ! Lateral diffusion of MEKE
6900 Kh_here = max(0., CS%MEKE_Kh)
691 !$OMP parallel do default(shared) firstprivate(Kh_here) private(Inv_Kh_max)
6920 do j=js,je ; do I=is-1,ie
693 ! Limit Kh to avoid CFL violations.
6940 if (allocated(MEKE%Kh)) &
695 Kh_here = max(0., CS%MEKE_Kh) + &
6960 CS%KhMEKE_Fac*0.5*(MEKE%Kh(i,j)+MEKE%Kh(i+1,j))
6970 if (allocated(MEKE%Kh_diff)) &
698 Kh_here = max(0.,CS%MEKE_Kh) + &
6990 CS%KhMEKE_Fac*0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i+1,j))
700 Inv_Kh_max = 2.0*sdt * ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
7010 max(G%IareaT(i,j),G%IareaT(i+1,j)))
7020 if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max
7030 Kh_u(I,j) = Kh_here
704
705 ! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
706 MEKE_uflux(I,j) = ((Kh_here * (G%dy_Cu(I,j)*G%IdxCu(I,j))) * &
707 ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
7080 (MEKE%MEKE(i,j) - MEKE%MEKE(i+1,j))
709 enddo ; enddo
710 !$OMP parallel do default(shared) firstprivate(Kh_here) private(Inv_Kh_max)
7110 do J=js-1,je ; do i=is,ie
7120 if (allocated(MEKE%Kh)) &
7130 Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac * 0.5*(MEKE%Kh(i,j)+MEKE%Kh(i,j+1))
7140 if (allocated(MEKE%Kh_diff)) &
7150 Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac * 0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i,j+1))
7160 Inv_Kh_max = 2.0*sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * max(G%IareaT(i,j),G%IareaT(i,j+1)))
7170 if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max
7180 Kh_v(i,J) = Kh_here
719
720 ! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
721 MEKE_vflux(i,J) = ((Kh_here * (G%dx_Cv(i,J)*G%IdyCv(i,J))) * &
722 ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
7230 (MEKE%MEKE(i,j) - MEKE%MEKE(i,j+1))
724 enddo ; enddo
7250 if (CS%MEKE_advection_factor>0.) then
7260 advFac = CS%MEKE_advection_factor / sdt ! [T-1 ~> s-1]
727 !$OMP parallel do default(shared)
7280 do j=js,je ; do I=is-1,ie
729 ! Here the units of the quantities added to MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].
7300 if (baroHu(I,j)>0.) then
7310 MEKE_uflux(I,j) = MEKE_uflux(I,j) + baroHu(I,j)*MEKE%MEKE(i,j)*advFac
7320 elseif (baroHu(I,j)<0.) then
7330 MEKE_uflux(I,j) = MEKE_uflux(I,j) + baroHu(I,j)*MEKE%MEKE(i+1,j)*advFac
734 endif
735 enddo ; enddo
736 !$OMP parallel do default(shared)
7370 do J=js-1,je ; do i=is,ie
738 ! Here the units of the quantities added to MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
7390 if (baroHv(i,J)>0.) then
7400 MEKE_vflux(i,J) = MEKE_vflux(i,J) + baroHv(i,J)*MEKE%MEKE(i,j)*advFac
7410 elseif (baroHv(i,J)<0.) then
7420 MEKE_vflux(i,J) = MEKE_vflux(i,J) + baroHv(i,J)*MEKE%MEKE(i,j+1)*advFac
743 endif
744 enddo ; enddo
745 endif
746
747 !$OMP parallel do default(shared)
7480 do j=js,je ; do i=is,ie
749 MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + (sdt*(G%IareaT(i,j)*I_mass(i,j))) * &
750 ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
7510 (MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
752 enddo ; enddo
753
7540 if (CS%id_src_adv > 0) then
755 !$OMP parallel do default(shared)
7560 do j=js,je ; do i=is,ie
757 src_adv(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * &
758 ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
7590 (MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
760 enddo ; enddo
761 endif
762 endif ! MEKE_KH>0
763
764 ! Add on bi-harmonic tendency
7650 if (CS%MEKE_K4 >= 0.0) then
766 !$OMP parallel do default(shared)
7670 do j=js,je ; do i=is,ie
7680 MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + del4MEKE(i,j)
769 enddo ; enddo
770 endif
771
772 ! Second stage of Strang splitting
7730 if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.0) then
774 ! Recalculate the drag rate, since MEKE has changed.
7750 if (use_drag_rate) then
776 !$OMP parallel do default(shared)
7770 do j=js,je ; do i=is,ie
778 drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
7790 cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
780 enddo ; enddo
781 endif
782
783 !$OMP parallel do default(shared)
7840 do j=js,je ; do i=is,ie
7850 damp_rate(i,j) = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
786
7870 if (MEKE%MEKE(i,j) < 0.) damp_rate(i,j) = 0.
788 ! notice that the above line ensures a damping only if MEKE is positive,
789 ! while leaving MEKE unchanged if it is negative
790 enddo ; enddo
791
792 ! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the
793 ! existing bit-reproducible solution.
794 !$OMP parallel do default(shared)
7950 do j=js,je ; do i=is,ie
7960 MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j))
797 enddo ; enddo
798
7990 if (any_damping_diags) then
800 !$OMP parallel do default(shared)
8010 do j=js,je ; do i=is,ie
8020 damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j))
803 enddo ; enddo
804
8050 if (CS%id_decay > 0) then
806 !$OMP parallel do default(shared)
8070 do j=js,je ; do i=is,ie
8080 MEKE_decay(i,j) = damp_rate(i,j) * G%mask2dT(i,j)
809 enddo ; enddo
810 endif
811
8120 if (CS%id_src_GM > 0) then
813 !$OMP parallel do default(shared)
8140 do j=js,je ; do i=is,ie
8150 src_GM(i,j) = src_GM(i,j) * damping(i,j)
816 enddo ; enddo
817 endif
818
8190 if (CS%id_src_mom_lp > 0) then
820 !$OMP parallel do default(shared)
8210 do j=js,je ; do i=is,ie
8220 src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j)
823 enddo ; enddo
824 endif
825
8260 if (CS%id_src_mom_bh > 0) then
827 !$OMP parallel do default(shared)
8280 do j=js,je ; do i=is,ie
8290 src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j)
830 enddo ; enddo
831 endif
832
8330 if (CS%id_src_adv > 0) then
834 !$OMP parallel do default(shared)
8350 do j=js,je ; do i=is,ie
8360 src_adv(i,j) = src_adv(i,j) * damping(i,j)
837 enddo ; enddo
838 endif
839
8400 if (CS%id_src_mom_K4 > 0) then
841 !$OMP parallel do default(shared)
8420 do j=js,je ; do i=is,ie
8430 src_mom_K4(i,j) = src_mom_K4(i,j) * damping(i,j)
844 enddo ; enddo
845 endif
846
8470 if (CS%id_src_btm_drag > 0) then
848 !$OMP parallel do default(shared)
8490 do j=js,je ; do i=is,ie
850 src_btm_drag(i,j) = -MEKE_current(i,j) * (damp_step &
851 * ((damp_rate(i,j) + damp_rate_s1(i,j)) * damping(i,j)) &
8520 )
853 enddo ; enddo
854 endif
855 endif
856 endif ! MEKE_KH>=0
857
8580 if (CS%debug) then
8590 call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, unscale=US%L_T_to_m_s**2)
860 endif
861
862 case(EKE_FILE)
8630 call time_interp_external(CS%eke_handle, Time, data_eke, scale=US%m_s_to_L_T**2)
8640 do j=js,je ; do i=is,ie
8650 MEKE%MEKE(i,j) = data_eke(i,j) * G%mask2dT(i,j)
866 enddo ; enddo
8670 call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale)
868 case(EKE_DBCLIENT)
8690 call pass_vector(u, v, G%Domain)
8700 call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale)
8710 call ML_MEKE_calculate_features(G, GV, US, CS, MEKE%Rd_dx_h, u, v, tv, h, dt, features_array)
8720 call predict_MEKE(G, US, CS, SIZE(h), Time, features_array, MEKE%MEKE)
873 case default
8740 call MOM_error(FATAL,"Invalid method specified for calculating EKE")
875 end select
876
8770 if (CS%MEKE_positive) then
878 !$OMP parallel do default(shared)
8790 do j=js,je ; do i=is,ie
8800 MEKE%MEKE(i,j) = MAX(0., MEKE%MEKE(i,j))
881 enddo ; enddo
882 endif
883
8840 call cpu_clock_begin(CS%id_clock_pass)
8850 call do_group_pass(CS%pass_MEKE, G%Domain)
8860 call cpu_clock_end(CS%id_clock_pass)
887
888 ! Calculate diffusivity for main model to use
8890 if (CS%MEKE_KhCoeff>0.) then
8900 if (.not.CS%MEKE_GEOMETRIC) then
8910 if (CS%use_old_lscale) then
8920 if (CS%Rd_as_max_scale) then
893 !$OMP parallel do default(shared)
8940 do j=js,je ; do i=is,ie
895 MEKE%Kh(i,j) = (CS%MEKE_KhCoeff * &
896 sqrt(2.*max(0.,barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j)) ) * &
8970 min(MEKE%Rd_dx_h(i,j), 1.0)
898 enddo ; enddo
899 else
900 !$OMP parallel do default(shared)
9010 do j=js,je ; do i=is,ie
902 MEKE%Kh(i,j) = CS%MEKE_KhCoeff * &
9030 sqrt(2.*max(0., barotrFac2(i,j)*MEKE%MEKE(i,j))*G%areaT(i,j))
904 enddo ; enddo
905 endif
906 else
907 !$OMP parallel do default(shared)
9080 do j=js,je ; do i=is,ie
909 MEKE%Kh(i,j) = CS%MEKE_KhCoeff * &
9100 sqrt(2.*max(0., barotrFac2(i,j)*MEKE%MEKE(i,j))) * LmixScale(i,j)
911 enddo ; enddo
912 endif
913 endif
914 endif
915
916 ! Calculate viscosity for the main model to use
9170 if (CS%viscosity_coeff_Ku /=0.) then
9180 do j=js,je ; do i=is,ie
9190 MEKE%Ku(i,j) = CS%viscosity_coeff_Ku * sqrt(2.*max(0.,MEKE%MEKE(i,j))) * LmixScale(i,j)
920 enddo ; enddo
921 endif
922
9230 if (CS%viscosity_coeff_Au /=0.) then
9240 do j=js,je ; do i=is,ie
9250 MEKE%Au(i,j) = CS%viscosity_coeff_Au * sqrt(2.*max(0.,MEKE%MEKE(i,j))) * LmixScale(i,j)**3
926 enddo ; enddo
927 endif
928
929 if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
9300 .or. allocated(MEKE%Le)) then
9310 call cpu_clock_begin(CS%id_clock_pass)
9320 call do_group_pass(CS%pass_Kh, G%Domain)
9330 call cpu_clock_end(CS%id_clock_pass)
934 endif
935
936 ! Offer fields for averaging.
9370 if (any([CS%id_Ue, CS%id_Ub, CS%id_Ut] > 0)) &
9380 tmp(:,:) = 0.
9390 if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag)
9400 if (CS%id_Ue>0) then
9410 do j=js,je ; do i=is,ie
9420 tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j)))
943 enddo ; enddo
9440 call post_data(CS%id_Ue, tmp, CS%diag)
945 endif
9460 if (CS%id_Ub>0) then
9470 do j=js,je ; do i=is,ie
9480 tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * bottomFac2(i,j)))
949 enddo ; enddo
9500 call post_data(CS%id_Ub, tmp, CS%diag)
951 endif
9520 if (CS%id_Ut>0) then
9530 do j=js,je ; do i=is,ie
9540 tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * barotrFac2(i,j)))
955 enddo ; enddo
9560 call post_data(CS%id_Ut, tmp, CS%diag)
957 endif
9580 if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag)
9590 if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag)
9600 if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag)
9610 if (CS%id_KhMEKE_u>0) call post_data(CS%id_KhMEKE_u, Kh_u, CS%diag)
9620 if (CS%id_KhMEKE_v>0) call post_data(CS%id_KhMEKE_v, Kh_v, CS%diag)
9630 if (CS%id_src>0) call post_data(CS%id_src, src, CS%diag)
9640 if (CS%id_src_adv>0) call post_data(CS%id_src_adv, src_adv, CS%diag)
9650 if (CS%id_src_mom_K4>0) call post_data(CS%id_src_mom_K4, src_mom_K4, CS%diag)
9660 if (CS%id_src_btm_drag>0) call post_data(CS%id_src_btm_drag, src_btm_drag, CS%diag)
9670 if (CS%id_src_GM>0) call post_data(CS%id_src_GM, src_GM, CS%diag)
9680 if (CS%id_src_mom_lp>0) call post_data(CS%id_src_mom_lp, src_mom_lp, CS%diag)
9690 if (CS%id_src_mom_bh>0) call post_data(CS%id_src_mom_bh, src_mom_bh, CS%diag)
9700 if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag)
9710 if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag)
9720 if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag)
9730 if (CS%id_mom_src_bh>0) call post_data(CS%id_mom_src_bh, MEKE%mom_src_bh, CS%diag)
9740 if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag)
9750 if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag)
9760 if (CS%id_gamma_b>0) then
9770 do j=js,je ; do i=is,ie
9780 bottomFac2(i,j) = sqrt(bottomFac2(i,j))
979 enddo ; enddo
9800 call post_data(CS%id_gamma_b, bottomFac2, CS%diag)
981 endif
9820 if (CS%id_gamma_t>0) then
9830 do j=js,je ; do i=is,ie
9840 barotrFac2(i,j) = sqrt(barotrFac2(i,j))
985 enddo ; enddo
9860 call post_data(CS%id_gamma_t, barotrFac2, CS%diag)
987 endif
988
989end subroutine step_forward_MEKE
990
991!> Calculates the equilibrium solution where the source depends only on MEKE diffusivity
992!! and there is no lateral diffusion of MEKE.
993!! Results is in MEKE%MEKE.
9940subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot)
995 type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
996 type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
997 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
998 type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
999 type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
1000 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
1001 real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
1002 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: drag_rate_visc !< Mean flow velocity contribution
1003 !! to the MEKE drag rate [H T-1 ~> m s-1 or kg m-2 s-1]
1004 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_mass !< Inverse of column mass [R-1 Z-1 ~> m2 kg-1].
1005 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].
1006
1007 ! Local variables
1008 real :: beta ! Combined topographic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]
1009 real :: SN ! The local Eady growth rate [T-1 ~> s-1]
1010 real :: bottomFac2, barotrFac2 ! Vertical structure factors [nondim]
1011 real :: LmixScale, LRhines, LEady ! Various mixing length scales [L ~> m]
1012 real :: KhCoeff ! A copy of MEKE_KhCoeff from the control structure [nondim]
1013 real :: Kh ! A lateral diffusivity [L2 T-1 ~> m2 s-1]
1014 real :: Ubg2 ! Background (tidal?) velocity squared [L2 T-2 ~> m2 s-2]
1015 real :: cd2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]
1016 real :: drag_rate ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
1017 real :: src ! The sum of MEKE sources [L2 T-3 ~> W kg-1]
1018 real :: ldamping ! The MEKE damping rate [T-1 ~> s-1].
1019 real :: EKE, EKEmin, EKEmax, EKEerr ! [L2 T-2 ~> m2 s-2]
1020 real :: resid, ResMin, ResMax ! Residuals [L2 T-3 ~> W kg-1]
1021 real :: FatH ! Coriolis parameter at h points, used to compute topographic beta [T-1 ~> s-1]
1022 real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]
1023 real :: h_neglect ! A negligible thickness [H ~> m or kg m-2]
1024 integer :: i, j, is, ie, js, je, n1, n2
1025 real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2].
1026 logical :: useSecant, debugIteration
1027
10280 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1029
10300 debugIteration = .false.
10310 KhCoeff = CS%MEKE_KhCoeff
10320 Ubg2 = CS%MEKE_Uscale**2
10330 cd2 = CS%cdrag**2
10340 tolerance = 1.0e-12*US%m_s_to_L_T**2
10350 h_neglect = GV%H_subroundoff
1036
1037!$OMP do
10380 do j=js,je ; do i=is,ie
1039 ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
1040 ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
10410 SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1))
1042
10430 if (CS%MEKE_equilibrium_alt) then
10440 MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * depth_tot(i,j))**2 / cd2
1045 else
1046 FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
10470 (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points
1048
1049 ! Since zero-bathymetry cells are masked, this avoids calculations on land
10500 if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then
10510 beta_topo_x = 0. ; beta_topo_y = 0.
1052 else
1053 !### Consider different combinations of these estimates of topographic beta.
1054 beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1055 (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) &
1056 / max(depth_tot(i+1,j), depth_tot(i,j), h_neglect) &
1057 + (depth_tot(i,j)-depth_tot(i-1,j)) * G%IdxCu(I-1,j) &
10580 / max(depth_tot(i,j), depth_tot(i-1,j), h_neglect) )
1059 beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1060 (depth_tot(i,j+1)-depth_tot(i,j)) * G%IdyCv(i,J) &
1061 / max(depth_tot(i,j+1), depth_tot(i,j), h_neglect) + &
1062 (depth_tot(i,j)-depth_tot(i,j-1)) * G%IdyCv(i,J-1) &
10630 / max(depth_tot(i,j), depth_tot(i,j-1), h_neglect) )
1064 endif
1065 beta = sqrt(((G%dF_dx(i,j) + beta_topo_x)**2) + &
10660 ((G%dF_dy(i,j) + beta_topo_y)**2) )
1067
10680 if (KhCoeff*SN*I_mass(i,j)>0.) then
1069 ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E
10700 EKEmin = 0. ! Use the trivial root as the left bracket
10710 ResMin = 0. ! Need to detect direction of left residual
10720 EKEmax = 0.01*US%m_s_to_L_T**2 ! First guess at right bracket
10730 useSecant = .false. ! Start using a bisection method
1074
1075 ! First find right bracket for which resid<0
10760 resid = 1.0*US%m_to_L**2*US%T_to_s**3 ; n1 = 0
10770 do while (resid>0.)
10780 n1 = n1 + 1
10790 EKE = EKEmax
1080 call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, depth_tot(i,j), &
1081 MEKE%Rd_dx_h(i,j), SN, EKE, &
10820 bottomFac2, barotrFac2, LmixScale, LRhines, LEady)
1083 ! TODO: Should include resolution function in Kh
10840 Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale)
10850 src = Kh * (SN * SN)
10860 drag_rate = (GV%H_to_RZ * I_mass(i,j)) * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) )
10870 ldamping = CS%MEKE_damping + drag_rate * bottomFac2
10880 resid = src - ldamping * EKE
1089 ! if (debugIteration) then
1090 ! write(0,*) n1, 'EKE=',EKE,'resid=',resid
1091 ! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin
1092 ! write(0,*) 'src=',src,'ldamping=',ldamping
1093 ! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2
1094 ! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2
1095 ! endif
10960 if (resid>0.) then ! EKE is to the left of the root
10970 EKEmin = EKE ! so we move the left bracket here
10980 EKEmax = 10. * EKE ! and guess again for the right bracket
10990 if (resid<ResMin) useSecant = .true.
11000 ResMin = resid
11010 if (EKEmax > 2.e17*US%m_s_to_L_T**2) then
11020 if (debugIteration) stop 'Something has gone very wrong'
11030 debugIteration = .true.
11040 resid = 1. ; n1 = 0
11050 EKEmin = 0. ; ResMin = 0.
11060 EKEmax = 0.01*US%m_s_to_L_T**2
11070 useSecant = .false.
1108 endif
1109 endif
1110 enddo ! while(resid>0.) searching for right bracket
11110 ResMax = resid
1112
1113 ! Bisect the bracket
11140 n2 = 0 ; EKEerr = EKEmax - EKEmin
11150 do while (EKEerr > tolerance)
11160 n2 = n2 + 1
11170 if (useSecant) then
11180 EKE = EKEmin + (EKEmax - EKEmin) * (ResMin / (ResMin - ResMax))
1119 else
11200 EKE = 0.5 * (EKEmin + EKEmax)
1121 endif
11220 EKEerr = min( EKE-EKEmin, EKEmax-EKE )
1123 ! TODO: Should include resolution function in Kh
11240 Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale)
11250 src = Kh * (SN * SN)
11260 drag_rate = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) )
11270 ldamping = CS%MEKE_damping + drag_rate * bottomFac2
11280 resid = src - ldamping * EKE
11290 if (useSecant .and. resid>ResMin) useSecant = .false.
11300 if (resid>0.) then ! EKE is to the left of the root
11310 EKEmin = EKE ! so we move the left bracket here
11320 if (resid<ResMin) useSecant = .true.
11330 ResMin = resid ! Save this for the secant method
11340 elseif (resid<0.) then ! EKE is to the right of the root
11350 EKEmax = EKE ! so we move the right bracket here
11360 ResMax = resid ! Save this for the secant method
1137 else
11380 exit ! resid=0 => EKE is exactly at the root
1139 endif
11400 if (n2>200) stop 'Failing to converge?'
1141 enddo ! while(EKEmax-EKEmin>tolerance)
1142
1143 else
11440 EKE = 0.
1145 endif
11460 MEKE%MEKE(i,j) = EKE
1147 endif
1148 enddo ; enddo
1149
11500end subroutine MEKE_equilibrium
1151
1152
1153!< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into
1154!! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value
11550subroutine MEKE_equilibrium_restoring(CS, G, GV, US, SN_u, SN_v, depth_tot, &
11560 equilibrium_value)
1157 type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
1158 type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1159 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type.
1160 type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1161 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
1162 real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
1163 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].
1164 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: equilibrium_value
1165 !< Equilbrium value of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2]
1166
1167 ! Local variables
1168 real :: SN ! The local Eady growth rate [T-1 ~> s-1]
1169 integer :: i, j, is, ie, js, je ! local indices
1170 real :: cd2 ! The square of the drag coefficient [nondim]
1171
11720 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
11730 cd2 = CS%cdrag**2
11740 equilibrium_value(:,:) = 0.0
1175
1176!$OMP do
11770 do j=js,je ; do i=is,ie
1178 ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
1179 ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
11800 SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1))
11810 equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * depth_tot(i,j))**2 / cd2
1182 enddo ; enddo
1183
11840 if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, equilibrium_value, CS%diag)
11850end subroutine MEKE_equilibrium_restoring
1186
1187!> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$
1188!! functions that are ratios of either bottom or barotropic eddy energy to the
1189!! column eddy energy, respectively. See \ref section_MEKE_equations.
11900subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, &
11910 bottomFac2, barotrFac2, LmixScale)
1192 type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1193 type(MEKE_type), intent(in) :: MEKE !< MEKE field
1194 type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
1195 type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1196 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1197 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
1198 real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
1199 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].
1200 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].
1201 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bottomFac2 !< gamma_b^2 [nondim]
1202 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: barotrFac2 !< gamma_t^2 [nondim]
1203 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: LmixScale !< Eddy mixing length [L ~> m].
1204 ! Local variables
12050 real, dimension(SZI_(G),SZJ_(G)) :: LRhines, LEady ! Possible mixing length scales [L ~> m]
1206 real :: beta ! Combined topographic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]
1207 real :: SN ! The local Eady growth rate [T-1 ~> s-1]
1208 real :: FatH ! Coriolis parameter at h points [T-1 ~> s-1]
1209 real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]
1210 real :: h_neglect ! A negligible thickness [H ~> m or kg m-2]
1211 integer :: i, j, is, ie, js, je
1212
12130 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
12140 h_neglect = GV%H_subroundoff
1215
1216!$OMP do
12170 do j=js,je ; do i=is,ie
12180 if (.not.CS%use_old_lscale) then
12190 if (CS%aEady > 0.) then
12200 SN = 0.25 * ( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)) )
1221 else
12220 SN = 0.
1223 endif
1224 FatH = 0.25* ( ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) ) + &
12250 ( G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1) ) ) ! Coriolis parameter at h points
1226
1227 ! If depth_tot is zero, then a division by zero FPE will be raised. In this
1228 ! case, we apply Adcroft's rule of reciprocals and set the term to zero.
1229 ! Since zero-bathymetry cells are masked, this should not affect values.
12300 if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then
12310 beta_topo_x = 0. ; beta_topo_y = 0.
1232 else
1233 !### Consider different combinations of these estimates of topographic beta.
1234 beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1235 (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) &
1236 / max(depth_tot(i+1,j), depth_tot(i,j), h_neglect) &
1237 + (depth_tot(i,j)-depth_tot(i-1,j)) * G%IdxCu(I-1,j) &
12380 / max(depth_tot(i,j), depth_tot(i-1,j), h_neglect) )
1239 beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1240 (depth_tot(i,j+1)-depth_tot(i,j)) * G%IdyCv(i,J) &
1241 / max(depth_tot(i,j+1), depth_tot(i,j), h_neglect) + &
1242 (depth_tot(i,j)-depth_tot(i,j-1)) * G%IdyCv(i,J-1) &
12430 / max(depth_tot(i,j), depth_tot(i,j-1), h_neglect) )
1244 endif
1245 beta = sqrt(((G%dF_dx(i,j) + beta_topo_x)**2) + &
12460 ((G%dF_dy(i,j) + beta_topo_y)**2) )
1247
1248 else
12490 beta = 0.
1250 endif
1251 ! Returns bottomFac2, barotrFac2 and LmixScale
1252 call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, depth_tot(i,j), &
1253 MEKE%Rd_dx_h(i,j), SN, MEKE%MEKE(i,j), &
1254 bottomFac2(i,j), barotrFac2(i,j), LmixScale(i,j), &
12550 LRhines(i,j), LEady(i,j))
1256 enddo ; enddo
12570 if (CS%id_Lrhines>0) call post_data(CS%id_LRhines, LRhines, CS%diag)
12580 if (CS%id_Leady>0) call post_data(CS%id_LEady, LEady, CS%diag)
1259
12600end subroutine MEKE_lengthScales
1261
1262!> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$
1263!! functions that are ratios of either bottom or barotropic eddy energy to the
1264!! column eddy energy, respectively. See \ref section_MEKE_equations.
12650subroutine MEKE_lengthScales_0d(CS, US, area, beta, depth_tot, Rd_dx, SN, EKE, &
1266 bottomFac2, barotrFac2, LmixScale, Lrhines, Leady)
1267 type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1268 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1269 real, intent(in) :: area !< Grid cell area [L2 ~> m2]
1270 real, intent(in) :: beta !< Planetary beta = \f$ \nabla f\f$ [T-1 L-1 ~> s-1 m-1]
1271 real, intent(in) :: depth_tot !< The total thickness of the water column [H ~> m or kg m-2]
1272 real, intent(in) :: Rd_dx !< Resolution Ld/dx [nondim].
1273 real, intent(in) :: SN !< Eady growth rate [T-1 ~> s-1].
1274 real, intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].
1275 real, intent(out) :: bottomFac2 !< gamma_b^2 [nondim]
1276 real, intent(out) :: barotrFac2 !< gamma_t^2 [nondim]
1277 real, intent(out) :: LmixScale !< Eddy mixing length [L ~> m].
1278 real, intent(out) :: Lrhines !< Rhines length scale [L ~> m].
1279 real, intent(out) :: Leady !< Eady length scale [L ~> m].
1280 ! Local variables
1281 real :: Lgrid, Ldeform, Lfrict ! Length scales [L ~> m]
1282 real :: Ue ! An eddy velocity [L T-1 ~> m s-1]
1283
1284 ! Length scale for MEKE derived diffusivity
12850 Lgrid = sqrt(area) ! Grid scale
12860 Ldeform = Lgrid * Rd_dx ! Deformation scale
12870 Lfrict = depth_tot / CS%cdrag ! Frictional arrest scale
1288 ! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy
1289 ! used in calculating bottom drag
12900 bottomFac2 = CS%MEKE_CD_SCALE**2
12910 if (Lfrict*CS%MEKE_Cb>0.) bottomFac2 = bottomFac2 + 1./( 1. + CS%MEKE_Cb*(Ldeform/Lfrict) )**0.8
12920 bottomFac2 = max(bottomFac2, CS%MEKE_min_gamma)
1293 ! gamma_t^2 is the ratio of barotropic eddy energy to mean column eddy energy
1294 ! used in the velocity scale for diffusivity
12950 barotrFac2 = 1.
12960 if (Lfrict*CS%MEKE_Ct>0.) barotrFac2 = 1. / ( 1. + CS%MEKE_Ct*(Ldeform/Lfrict) )**0.25
12970 barotrFac2 = max(barotrFac2, CS%MEKE_min_gamma)
12980 if (CS%use_old_lscale) then
12990 if (CS%Rd_as_max_scale) then
13000 LmixScale = min(Ldeform, Lgrid) ! The smaller of Ld or dx
1301 else
13020 LmixScale = Lgrid
1303 endif
1304 else
13050 Ue = sqrt( 2.0 * max( 0., barotrFac2*EKE ) ) ! Barotropic eddy flow scale
13060 Lrhines = sqrt( Ue / max( beta, 1.e-30*US%T_to_s*US%L_to_m ) ) ! Rhines scale
13070 if (CS%aEady > 0.) then
13080 Leady = Ue / max( SN, 1.e-15*US%T_to_s ) ! Bound Eady time-scale < 1e15 seconds
1309 else
13100 Leady = 0.
1311 endif
13120 if (CS%use_min_lscale) then
13130 LmixScale = CS%lscale_maxval
13140 if (CS%aDeform*Ldeform > 0.) LmixScale = min(LmixScale,CS%aDeform*Ldeform)
13150 if (CS%aFrict *Lfrict > 0.) LmixScale = min(LmixScale,CS%aFrict *Lfrict)
13160 if (CS%aRhines*Lrhines > 0.) LmixScale = min(LmixScale,CS%aRhines*Lrhines)
13170 if (CS%aEady *Leady > 0.) LmixScale = min(LmixScale,CS%aEady *Leady)
13180 if (CS%aGrid *Lgrid > 0.) LmixScale = min(LmixScale,CS%aGrid *Lgrid)
13190 if (CS%Lfixed > 0.) LmixScale = min(LmixScale,CS%Lfixed)
1320 else
13210 LmixScale = 0.
13220 if (CS%aDeform*Ldeform > 0.) LmixScale = LmixScale + 1./(CS%aDeform*Ldeform)
13230 if (CS%aFrict *Lfrict > 0.) LmixScale = LmixScale + 1./(CS%aFrict *Lfrict)
13240 if (CS%aRhines*Lrhines > 0.) LmixScale = LmixScale + 1./(CS%aRhines*Lrhines)
13250 if (CS%aEady *Leady > 0.) LmixScale = LmixScale + 1./(CS%aEady *Leady)
13260 if (CS%aGrid *Lgrid > 0.) LmixScale = LmixScale + 1./(CS%aGrid *Lgrid)
13270 if (CS%Lfixed > 0.) LmixScale = LmixScale + 1./CS%Lfixed
13280 if (LmixScale > 0.) LmixScale = 1. / LmixScale
1329 endif
1330 endif
1331
13320end subroutine MEKE_lengthScales_0d
1333
1334!> Initializes the MOM_MEKE module and reads parameters.
1335!! Returns True if module is to be used, otherwise returns False.
13361logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, MEKE, restart_CS, meke_in_dynamics)
1337 type(time_type), intent(in) :: Time !< The current model time.
1338 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1339 type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1340 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1341 type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
1342 type(dbcomms_CS_type), intent(in) :: dbcomms_CS !< Database communications control structure
1343 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
1344 type(MEKE_CS), intent(inout) :: CS !< MEKE control structure.
1345 type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
1346 type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control structure
1347 logical, intent( out) :: meke_in_dynamics !< If true, MEKE is stepped forward in dynamics
1348 !! otherwise in tracer dynamics
1349
1350 ! Local variables
1351 real :: MEKE_restoring_timescale ! The timescale used to nudge MEKE toward its equilibrium value [T ~> s]
1352 real :: cdrag ! The default bottom drag coefficient [nondim].
1353 character(len=200) :: eke_filename, eke_varname, inputdir
1354 character(len=16) :: eke_source_str
1355 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1356 logical :: laplacian, biharmonic, coldStart
1357 ! This include declares and sets the variable "version".
1358# include "version_variable.h"
1359 character(len=40) :: mdl = "MOM_MEKE" ! This module's name.
1360
13611 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
13621 isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1363
1364 ! Determine whether this module will be used
13651 call get_param(param_file, mdl, "USE_MEKE", MEKE_init, default=.false., do_not_log=.true.)
13661 call log_version(param_file, mdl, version, "", all_default=.not.MEKE_init)
1367 call get_param(param_file, mdl, "USE_MEKE", MEKE_init, &
1368 "If true, turns on the MEKE scheme which calculates "// &
1369 "a sub-grid mesoscale eddy kinetic energy budget.", &
13701 default=.false.)
13711 if (.not. MEKE_init) return
13720 CS%initialized = .true.
1373 call get_param(param_file, mdl, "MEKE_IN_DYNAMICS", meke_in_dynamics, &
1374 "If true, step MEKE forward with the dynamics "// &
1375 "otherwise with the tracer timestep.", &
13760 default=.true.)
1377
1378 call get_param(param_file, mdl, "EKE_SOURCE", eke_source_str, &
1379 "Determine the where EKE comes from:\n" // &
1380 " 'prog': Calculated solving EKE equation\n"// &
1381 " 'file': Read in from a file\n" // &
13820 " 'dbclient': Retrieved from ML-database", default='prog')
1383
13840 call MOM_mesg("MEKE_init: reading parameters ", 5)
1385
13860 select case (lowercase(eke_source_str))
1387 case("file")
13880 CS%eke_src = EKE_FILE
13890 call time_interp_external_init
1390 call get_param(param_file, mdl, "EKE_FILE", eke_filename, &
1391 "A file in which to find the eddy kineteic energy variable.", &
13920 default="eke_file.nc")
1393 call get_param(param_file, mdl, "EKE_VARIABLE", eke_varname, &
1394 "The name of the eddy kinetic energy variable to read from "//&
1395 "EKE_FILE to use in MEKE.", &
13960 default="eke")
1397 call get_param(param_file, mdl, "INPUTDIR", inputdir, &
1398 "The directory in which all input files are found.", &
13990 default=".", do_not_log=.true.)
14000 inputdir = slasher(inputdir)
1401
14020 eke_filename = trim(inputdir) // trim(eke_filename)
14030 CS%eke_handle = init_external_field(eke_filename, eke_varname, domain=G%Domain%mpp_domain)
1404 case("prog")
14050 CS%eke_src = EKE_PROG
1406 ! Read all relevant parameters and write them to the model log.
1407 call get_param(param_file, mdl, "MEKE_DAMPING", CS%MEKE_damping, &
1408 "The local depth-independent MEKE dissipation rate.", &
14090 units="s-1", default=0.0, scale=US%T_to_s)
1410 call get_param(param_file, mdl, "MEKE_CD_SCALE", CS%MEKE_Cd_scale, &
1411 "The ratio of the bottom eddy velocity to the column mean "//&
1412 "eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1 "//&
1413 "to account for the surface intensification of MEKE.", &
14140 units="nondim", default=0.)
1415 call get_param(param_file, mdl, "MEKE_CB", CS%MEKE_Cb, &
1416 "A coefficient in the expression for the ratio of bottom projected "//&
1417 "eddy energy and mean column energy (see Jansen et al. 2015).",&
14180 units="nondim", default=25.)
1419 call get_param(param_file, mdl, "MEKE_MIN_GAMMA2", CS%MEKE_min_gamma, &
1420 "The minimum allowed value of gamma_b^2.",&
14210 units="nondim", default=0.0001)
1422 call get_param(param_file, mdl, "MEKE_CT", CS%MEKE_Ct, &
1423 "A coefficient in the expression for the ratio of barotropic "//&
1424 "eddy energy and mean column energy (see Jansen et al. 2015).",&
14250 units="nondim", default=50.)
1426 call get_param(param_file, mdl, "MEKE_GMCOEFF", CS%MEKE_GMcoeff, &
1427 "The efficiency of the conversion of potential energy "//&
1428 "into MEKE by the thickness mixing parameterization. "//&
1429 "If MEKE_GMCOEFF is negative, this conversion is not "//&
14300 "used or calculated.", units="nondim", default=-1.0)
1431 call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, &
1432 "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//&
14330 "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
1434 call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, &
1435 "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&
14360 "thickness diffusion.", units="nondim", default=0.05)
1437 call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", CS%MEKE_equilibrium_alt, &
1438 "If true, use an alternative formula for computing the (equilibrium) "//&
14390 "initial value of MEKE.", default=.false.)
1440 call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", CS%MEKE_equilibrium_restoring, &
1441 "If true, restore MEKE back to its equilibrium value, which is calculated at "//&
14420 "each time step.", default=.false.)
14430 if (CS%MEKE_equilibrium_restoring) then
1444 call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, &
1445 "The timescale used to nudge MEKE toward its equilibrium value.", &
14460 units="s", default=1e6, scale=US%s_to_T)
14470 CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale
1448 endif
1449
1450 call get_param(param_file, mdl, "MEKE_FRCOEFF", CS%MEKE_FrCoeff, &
1451 "The efficiency of the conversion of mean energy into "//&
1452 "MEKE. If MEKE_FRCOEFF is negative, this conversion "//&
14530 "is not used or calculated.", units="nondim", default=-1.0)
1454 call get_param(param_file, mdl, "MEKE_BHFRCOEFF", CS%MEKE_bhFrCoeff, &
1455 "The efficiency of the conversion of mean energy into "//&
1456 "MEKE by the biharmonic dissipation. If MEKE_bhFRCOEFF is negative, this conversion "//&
14570 "is not used or calculated.", units="nondim", default=-1.0)
1458 call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, &
1459 "The efficiency of the conversion of MEKE into mean energy "//&
1460 "by GME. If MEKE_GMECOEFF is negative, this conversion "//&
14610 "is not used or calculated.", units="nondim", default=-1.0)
1462 call get_param(param_file, mdl, "MEKE_BGSRC", CS%MEKE_BGsrc, &
1463 "A background energy source for MEKE.", &
14640 units="W kg-1", default=0.0, scale=US%m_to_L**2*US%T_to_s**3)
1465 call get_param(param_file, mdl, "MEKE_KH", CS%MEKE_Kh, &
1466 "A background lateral diffusivity of MEKE. "//&
1467 "Use a negative value to not apply lateral diffusion to MEKE.", &
14680 units="m2 s-1", default=-1.0, scale=US%m_to_L**2*US%T_to_s)
1469 call get_param(param_file, mdl, "MEKE_K4", CS%MEKE_K4, &
1470 "A lateral bi-harmonic diffusivity of MEKE. "//&
1471 "Use a negative value to not apply bi-harmonic diffusion to MEKE.", &
14720 units="m4 s-1", default=-1.0, scale=US%m_to_L**4*US%T_to_s)
1473 call get_param(param_file, mdl, "MEKE_DTSCALE", CS%MEKE_dtScale, &
1474 "A scaling factor to accelerate the time evolution of MEKE.", &
14750 units="nondim", default=1.0)
1476 call get_param(param_file, mdl, "MEKE_POSITIVE", CS%MEKE_positive, &
1477 "If true, it guarantees that MEKE will always be >= 0.", &
14780 default=.false.)
1479 case("dbclient")
14800 CS%eke_src = EKE_DBCLIENT
14810 call ML_MEKE_init(diag, G, US, Time, param_file, dbcomms_CS, CS)
1482 case default
14830 call MOM_error(FATAL, "Invalid method selected for calculating EKE")
1484 end select
1485 ! GMM, make sure all parameters used to calculated MEKE are within the above if
1486
1487 call get_param(param_file, mdl, "MEKE_KHCOEFF", CS%MEKE_KhCoeff, &
1488 "A scaling factor in the expression for eddy diffusivity "//&
1489 "which is otherwise proportional to the MEKE velocity- "//&
1490 "scale times an eddy mixing-length. This factor "//&
1491 "must be >0 for MEKE to contribute to the thickness/ "//&
1492 "and tracer diffusivity in the rest of the model.", &
14930 units="nondim", default=1.0)
1494 call get_param(param_file, mdl, "MEKE_USCALE", CS%MEKE_Uscale, &
1495 "The background velocity that is combined with MEKE to "//&
14960 "calculate the bottom drag.", units="m s-1", default=0.0, scale=US%m_s_to_L_T)
1497 call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, &
1498 "If true, use the GM energy conversion form S^2*N^2*kappa rather "//&
14990 "than the streamfunction for the MEKE GM source term.", default=.false.)
1500 call get_param(param_file, mdl, "MEKE_MIN_DEPTH_TOT", CS%MEKE_min_depth_tot, &
1501 "The minimum total depth over which to distribute MEKE energy sources. "//&
1502 "When the total depth is less than this, the sources are scaled away.", &
15030 units="m", default=1.0, scale=GV%m_to_H, do_not_log=.not.CS%GM_src_alt)
1504 call get_param(param_file, mdl, "MEKE_VISC_DRAG", CS%visc_drag, &
1505 "If true, use the vertvisc_type to calculate the bottom "//&
15060 "drag acting on MEKE.", default=.true.)
1507 call get_param(param_file, mdl, "MEKE_KHTH_FAC", MEKE%KhTh_fac, &
15080 "A factor that maps MEKE%Kh to KhTh.", units="nondim", default=0.0)
1509 call get_param(param_file, mdl, "MEKE_KHTR_FAC", MEKE%KhTr_fac, &
15100 "A factor that maps MEKE%Kh to KhTr.", units="nondim", default=0.0)
1511 call get_param(param_file, mdl, "MEKE_KHMEKE_FAC", CS%KhMEKE_Fac, &
1512 "A factor that maps MEKE%Kh to Kh for MEKE itself.", &
15130 units="nondim", default=0.0)
1514 call get_param(param_file, mdl, "MEKE_OLD_LSCALE", CS%use_old_lscale, &
1515 "If true, use the old formula for length scale which is "//&
1516 "a function of grid spacing and deformation radius.", &
15170 default=.false.)
1518 call get_param(param_file, mdl, "MEKE_MIN_LSCALE", CS%use_min_lscale, &
1519 "If true, use a strict minimum of provided length scales "//&
1520 "rather than harmonic mean.", &
15210 default=.false.)
1522 call get_param(param_file, mdl, "MEKE_LSCALE_MAX_VAL", CS%lscale_maxval, &
1523 "The ceiling on the value of the MEKE length scale when MEKE_MIN_LSCALE=True. "//&
1524 "The default is the distance from the equator to the pole on Earth, as "//&
1525 "estimated by enlightenment era scientists, but should probably scale with RAD_EARTH.", &
15260 units="m", default=1.0e7, scale=US%m_to_L, do_not_log=.not.CS%use_min_lscale)
1527 call get_param(param_file, mdl, "MEKE_RD_MAX_SCALE", CS%Rd_as_max_scale, &
1528 "If true, the length scale used by MEKE is the minimum of "//&
1529 "the deformation radius or grid-spacing. Only used if "//&
15300 "MEKE_OLD_LSCALE=True", default=.false.)
1531 call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_KU", CS%viscosity_coeff_Ku, &
1532 "If non-zero, is the scaling coefficient in the expression for "//&
1533 "viscosity used to parameterize harmonic lateral momentum mixing by "//&
1534 "unresolved eddies represented by MEKE. Can be negative to "//&
1535 "represent backscatter from the unresolved eddies.", &
15360 units="nondim", default=0.0)
1537 call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_AU", CS%viscosity_coeff_Au, &
1538 "If non-zero, is the scaling coefficient in the expression for "//&
1539 "viscosity used to parameterize biharmonic lateral momentum mixing by "//&
1540 "unresolved eddies represented by MEKE. Can be negative to "//&
1541 "represent backscatter from the unresolved eddies.", &
15420 units="nondim", default=0.0)
1543 call get_param(param_file, mdl, "MEKE_FIXED_MIXING_LENGTH", CS%Lfixed, &
1544 "If positive, is a fixed length contribution to the expression "//&
1545 "for mixing length used in MEKE-derived diffusivity.", &
15460 units="m", default=0.0, scale=US%m_to_L)
1547 call get_param(param_file, mdl, "MEKE_FIXED_TOTAL_DEPTH", CS%fixed_total_depth, &
1548 "If true, use the nominal bathymetric depth as the estimate of the "//&
1549 "time-varying ocean depth. Otherwise base the depth on the total ocean mass "//&
15500 "per unit area.", default=.true.)
1551 call get_param(param_file, mdl, "MEKE_TOTAL_DEPTH_RHO", CS%rho_fixed_total_depth, &
1552 "A density used to translate the nominal bathymetric depth into an estimate "//&
1553 "of the total ocean mass per unit area when MEKE_FIXED_TOTAL_DEPTH is true.", &
1554 units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R, &
15550 do_not_log=(GV%Boussinesq.or.(.not.CS%fixed_total_depth)))
1556
1557 call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", CS%aDeform, &
1558 "If positive, is a coefficient weighting the deformation scale "//&
1559 "in the expression for mixing length used in MEKE-derived diffusivity.", &
15600 units="nondim", default=0.0)
1561 call get_param(param_file, mdl, "MEKE_ALPHA_RHINES", CS%aRhines, &
1562 "If positive, is a coefficient weighting the Rhines scale "//&
1563 "in the expression for mixing length used in MEKE-derived diffusivity.", &
15640 units="nondim", default=0.0)
1565 call get_param(param_file, mdl, "MEKE_ALPHA_EADY", CS%aEady, &
1566 "If positive, is a coefficient weighting the Eady length scale "//&
1567 "in the expression for mixing length used in MEKE-derived diffusivity.", &
15680 units="nondim", default=0.0)
1569 call get_param(param_file, mdl, "MEKE_ALPHA_FRICT", CS%aFrict, &
1570 "If positive, is a coefficient weighting the frictional arrest scale "//&
1571 "in the expression for mixing length used in MEKE-derived diffusivity.", &
15720 units="nondim", default=0.0)
1573 call get_param(param_file, mdl, "MEKE_ALPHA_GRID", CS%aGrid, &
1574 "If positive, is a coefficient weighting the grid-spacing as a scale "//&
1575 "in the expression for mixing length used in MEKE-derived diffusivity.", &
15760 units="nondim", default=0.0)
1577 call get_param(param_file, mdl, "MEKE_COLD_START", coldStart, &
1578 "If true, initialize EKE to zero. Otherwise a local equilibrium solution "//&
15790 "is used as an initial condition for EKE.", default=.false.)
1580 call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", MEKE%backscatter_Ro_c, &
1581 "The coefficient in the Rossby number function for scaling the biharmonic "//&
1582 "frictional energy source. Setting to non-zero enables the Rossby number function.", &
15830 units="nondim", default=0.0)
1584 call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_POW", MEKE%backscatter_Ro_pow, &
1585 "The power in the Rossby number function for scaling the biharmonic "//&
15860 "frictional energy source.", units="nondim", default=0.0)
1587 call get_param(param_file, mdl, "MEKE_ADVECTION_FACTOR", CS%MEKE_advection_factor, &
1588 "A scale factor in front of advection of eddy energy. Zero turns advection off. "//&
1589 "Using unity would be normal but other values could accommodate a mismatch "//&
1590 "between the advecting barotropic flow and the vertical structure of MEKE.", &
15910 units="nondim", default=0.0)
1592 call get_param(param_file, mdl, "MEKE_ADVECTION_BUG", CS%MEKE_advection_bug, &
1593 "If true, recover a bug in the calculation of the barotropic transport for "//&
1594 "the advection of MEKE. With the bug, only the transports in the deepest "//&
15950 "layer are used.", default=.false., do_not_log=(CS%MEKE_advection_factor<=0.))
1596 call get_param(param_file, mdl, "MEKE_TOPOGRAPHIC_BETA", CS%MEKE_topographic_beta, &
1597 "A scale factor to determine how much topographic beta is weighed in " //&
1598 "computing beta in the expression of Rhines scale. Use 1 if full "//&
1599 "topographic beta effect is considered; use 0 if it's completely ignored.", &
16000 units="nondim", default=0.0)
1601 call get_param(param_file, mdl, "SQG_USE_MEKE", CS%sqg_use_MEKE, &
1602 "If true, the eddy scale of MEKE is used for the SQG vertical structure ",&
16030 default=.false.)
1604
1605 ! Nonlocal module parameters
1606 call get_param(param_file, mdl, "CDRAG", cdrag, &
1607 "CDRAG is the drag coefficient relating the magnitude of the velocity "//&
16080 "field to the bottom stress.", units="nondim", default=0.003)
1609 call get_param(param_file, mdl, "MEKE_CDRAG", CS%cdrag, &
1610 "Drag coefficient relating the magnitude of the velocity "//&
16110 "field to the bottom stress in MEKE.", units="nondim", default=cdrag, scale=US%L_to_m*GV%m_to_H)
16120 call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
16130 call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)
1614
16150 if (CS%viscosity_coeff_Ku/=0. .and. .not. laplacian) call MOM_error(FATAL, &
16160 "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF_KU is true.")
1617
16180 if (CS%viscosity_coeff_Au/=0. .and. .not. biharmonic) call MOM_error(FATAL, &
16190 "BIHARMONIC must be true if MEKE_VISCOSITY_COEFF_AU is true.")
1620
16210 call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
1622
1623 ! Identify if any lateral diffusive processes are active
16240 CS%kh_flux_enabled = .false.
16250 if ((CS%MEKE_KH >= 0.0) .or. (CS%KhMEKE_FAC > 0.0) .or. (CS%MEKE_advection_factor > 0.0)) &
16260 CS%kh_flux_enabled = .true.
1627
1628! Register fields for output from this module.
16290 CS%diag => diag
1630 CS%id_MEKE = register_diag_field('ocean_model', 'MEKE', diag%axesT1, Time, &
16310 'Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=US%L_T_to_m_s**2)
16320 if (.not. allocated(MEKE%MEKE)) CS%id_MEKE = -1
1633 CS%id_Kh = register_diag_field('ocean_model', 'MEKE_KH', diag%axesT1, Time, &
16340 'MEKE derived diffusivity', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
16350 if (.not. allocated(MEKE%Kh)) CS%id_Kh = -1
1636 CS%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, Time, &
16370 'MEKE derived lateral viscosity', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
16380 if (.not. allocated(MEKE%Ku)) CS%id_Ku = -1
1639 CS%id_Au = register_diag_field('ocean_model', 'MEKE_AU', diag%axesT1, Time, &
16400 'MEKE derived lateral biharmonic viscosity', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T)
16410 if (.not. allocated(MEKE%Au)) CS%id_Au = -1
1642 CS%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, Time, &
16430 'MEKE derived eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s)
16440 if (.not. allocated(MEKE%MEKE)) CS%id_Ue = -1
1645 CS%id_Ub = register_diag_field('ocean_model', 'MEKE_Ub', diag%axesT1, Time, &
16460 'MEKE derived bottom eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s)
16470 if (.not. allocated(MEKE%MEKE)) CS%id_Ub = -1
1648 CS%id_Ut = register_diag_field('ocean_model', 'MEKE_Ut', diag%axesT1, Time, &
16490 'MEKE derived barotropic eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s)
16500 if (.not. allocated(MEKE%MEKE)) CS%id_Ut = -1
1651 CS%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, Time, &
16520 'MEKE energy source', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1653
1654 CS%id_src_adv = register_diag_field('ocean_model', 'MEKE_src_adv', diag%axesT1, Time, &
16550 'MEKE energy source from the horizontal advection of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1656
1657 CS%id_src_btm_drag = register_diag_field('ocean_model', 'MEKE_src_btm_drag', diag%axesT1, Time, &
16580 'MEKE energy source from the bottom drag acting on MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1659
16600 if (CS%MEKE_K4 >= 0.) &
1661 CS%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', &
1662 diag%axesT1, Time, 'MEKE energy source from the biharmonic of MEKE', &
16630 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1664
16650 if (CS%MEKE_GMcoeff >= 0.) &
1666 CS%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', &
1667 diag%axesT1, Time, 'MEKE energy source from the thickness mixing (GM scheme)', &
16680 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1669
16700 if (CS%MEKE_FrCoeff >= 0.) &
1671 CS%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', &
1672 diag%axesT1, Time, 'MEKE energy source from the Laplacian of resolved flows', &
16730 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1674
16750 if (CS%MEKE_bhFrCoeff >= 0.) &
1676 CS%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', &
1677 diag%axesT1, Time, 'MEKE energy source from the biharmonic of resolved flows', &
16780 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T)
1679
1680 CS%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, Time, &
16810 'MEKE decay rate', 's-1', conversion=US%s_to_T)
1682 CS%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, Time, &
1683 'MEKE energy available from thickness mixing', &
16840 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
16850 if (.not. allocated(MEKE%GM_src)) CS%id_GM_src = -1
1686 CS%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, Time, &
1687 'MEKE energy available from momentum', &
16880 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
16890 if (.not. allocated(MEKE%mom_src)) CS%id_mom_src = -1
1690 CS%id_mom_src_bh = register_diag_field('ocean_model', 'MEKE_mom_src_bh',diag%axesT1, Time, &
1691 'MEKE energy available from the biharmonic dissipation of momentum', &
16920 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
16930 if (.not. allocated(MEKE%mom_src_bh)) CS%id_mom_src_bh = -1
1694 CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, &
1695 'MEKE energy lost to GME backscatter', &
16960 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
16970 if (.not. allocated(MEKE%GME_snk)) CS%id_GME_snk = -1
1698 CS%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, Time, &
16990 'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm', conversion=US%L_to_m)
1700 CS%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, Time, &
17010 'Rhines length scale used in the MEKE derived eddy diffusivity', 'm', conversion=US%L_to_m)
1702 CS%id_Leady = register_diag_field('ocean_model', 'MEKE_Leady', diag%axesT1, Time, &
17030 'Eady length scale used in the MEKE derived eddy diffusivity', 'm', conversion=US%L_to_m)
1704 CS%id_gamma_b = register_diag_field('ocean_model', 'MEKE_gamma_b', diag%axesT1, Time, &
17050 'Ratio of bottom-projected eddy velocity to column-mean eddy velocity', 'nondim')
1706 CS%id_gamma_t = register_diag_field('ocean_model', 'MEKE_gamma_t', diag%axesT1, Time, &
17070 'Ratio of barotropic eddy velocity to column-mean eddy velocity', 'nondim')
1708
17090 if (CS%kh_flux_enabled) then
1710 CS%id_KhMEKE_u = register_diag_field('ocean_model', 'KHMEKE_u', diag%axesCu1, Time, &
17110 'Zonal diffusivity of MEKE', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
1712 CS%id_KhMEKE_v = register_diag_field('ocean_model', 'KHMEKE_v', diag%axesCv1, Time, &
17130 'Meridional diffusivity of MEKE', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
1714 endif
1715
17160 if (CS%MEKE_equilibrium_restoring) then
1717 CS%id_MEKE_equilibrium = register_diag_field('ocean_model', 'MEKE_equilibrium', diag%axesT1, Time, &
17180 'Equilibrated Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=US%L_T_to_m_s**2)
1719 endif
1720
17210 CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE)
1722
1723
1724 ! Detect whether this instance of MEKE_init() is at the beginning of a run
1725 ! or after a restart. If at the beginning, we will initialize MEKE to a local
1726 ! equilibrium.
17270 CS%initialize = .not.query_initialized(MEKE%MEKE, "MEKE", restart_CS)
17280 if (coldStart) CS%initialize = .false.
17290 if (CS%initialize) call MOM_error(WARNING, &
17300 "MEKE_init: Initializing MEKE with a local equilibrium balance.")
17310 if (allocated(MEKE%Le)) then
17320 if (.not.query_initialized(MEKE%Le, "MEKE_Le", restart_CS)) then
1733 !$OMP parallel do default(shared)
17340 do j=js,je ; do i=is,ie
17350 MEKE%Le(i,j) = sqrt(G%areaT(i,j))
1736 enddo ; enddo
1737 endif
1738 endif
1739
1740 ! Set up group passes. In the case of a restart, these fields need a halo update now.
17410 if (allocated(MEKE%MEKE)) then
17420 call create_group_pass(CS%pass_MEKE, MEKE%MEKE, G%Domain)
17430 if (allocated(MEKE%Kh_diff)) call create_group_pass(CS%pass_MEKE, MEKE%Kh_diff, G%Domain)
17440 if (.not.CS%initialize) call do_group_pass(CS%pass_MEKE, G%Domain)
1745 endif
17460 if (allocated(MEKE%Kh)) call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain)
17470 if (allocated(MEKE%Ku)) call create_group_pass(CS%pass_Kh, MEKE%Ku, G%Domain)
17480 if (allocated(MEKE%Au)) call create_group_pass(CS%pass_Kh, MEKE%Au, G%Domain)
17490 if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain)
1750
1751 if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
17520 .or. allocated(MEKE%Le)) &
17530 call do_group_pass(CS%pass_Kh, G%Domain)
1754
17550end function MEKE_init
1756
1757!> Initializer for the variant of MEKE that uses ML to predict eddy kinetic energy
17580subroutine ML_MEKE_init(diag, G, US, Time, param_file, dbcomms_CS, CS)
1759 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
1760 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1761 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1762 type(time_type), intent(in) :: Time !< The current model time.
1763 type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
1764 type(dbcomms_CS_type), intent(in) :: dbcomms_CS !< Control structure for database communication
1765 type(MEKE_CS), intent(inout) :: CS !< Control structure for this module
1766
1767 character(len=200) :: inputdir, backend, model_filename
1768 integer :: db_return_code, batch_size
1769 character(len=40) :: mdl = "MOM_ML_MEKE"
1770#ifdef __NVCOMPILER_OPENMP_GPU
1771 integer, parameter :: default_niblock = 0
1772 integer, parameter :: default_njblock = 0
1773 integer, parameter :: default_nkblock = 0
1774#else
1775 integer, parameter :: default_niblock = 0
1776 integer, parameter :: default_njblock = 1
1777 integer, parameter :: default_nkblock = 1
1778#endif
1779
1780 ! Store pointers in control structure
17810 write(CS%key_suffix, '(A,I6.6)') '_', PE_here()
1782 ! Put some basic information into the database
17830 db_return_code = 0
1784 db_return_code = CS%client%put_tensor("meta"//CS%key_suffix, &
17850 REAL([G%isd_global, G%idg_offset, G%jsd_global, G%jdg_offset]),[4]) + db_return_code
17860 db_return_code = CS%client%put_tensor("geolat"//CS%key_suffix, G%geoLatT, shape(G%geoLatT)) + db_return_code
17870 db_return_code = CS%client%put_tensor("geolon"//CS%key_suffix, G%geoLonT, shape(G%geoLonT)) + db_return_code
17880 db_return_code = CS%client%put_tensor("EKE_shape"//CS%key_suffix, shape(G%geolonT), [2]) + db_return_code
1789
17900 if (CS%client%SR_error_parser(db_return_code)) call MOM_error(FATAL, "Putting metadata into the database failed")
1791
17920 call read_param(param_file, "INPUTDIR", inputdir)
17930 inputdir = slasher(inputdir)
1794
17950 call get_param(param_file, mdl, "BATCH_SIZE", batch_size, "Batch size to use for inference", default=1)
1796 call get_param(param_file, mdl, "EKE_BACKEND", backend, &
17970 "The computational backend to use for EKE inference (CPU or GPU)", default="GPU")
1798 call get_param(param_file, mdl, "EKE_MODEL", model_filename, &
17990 "Filename of the a saved pyTorch model to use", fail_if_missing = .true.)
1800 call get_param(param_file, mdl, "EKE_MAX", CS%eke_max, &
1801 "Maximum value of EKE allowed when inferring EKE", &
18020 units="m2 s-2", default=2., scale=US%m_s_to_L_T**2)
1803
1804 ! Set the machine learning model
18050 if (dbcomms_CS%colocated) then
18060 if (modulo(PE_here(),dbcomms_CS%colocated_stride) == 0) then
1807 db_return_code = CS%client%set_model_from_file(CS%model_key, trim(inputdir)//trim(model_filename), &
18080 "TORCH", backend, batch_size=batch_size)
1809 endif
1810 else
18110 if (is_root_pe()) then
1812 db_return_code = CS%client%set_model_from_file(CS%model_key, trim(inputdir)//trim(model_filename), &
18130 "TORCH", backend, batch_size=batch_size)
1814 endif
1815 endif
18160 if (CS%client%SR_error_parser(db_return_code)) then
18170 call MOM_error(FATAL, "MEKE: set_model failed")
1818 endif
1819
1820 call get_param(param_file, mdl, "ONLINE_ANALYSIS", CS%online_analysis, &
18210 "If true, post EKE used in MOM6 to the database for analysis", default=.true.)
1822
1823 ! Set various clock ids
18240 CS%id_client_init = cpu_clock_id('(ML_MEKE client init)', grain=CLOCK_ROUTINE)
18250 CS%id_put_tensor = cpu_clock_id('(ML_MEKE put tensor)', grain=CLOCK_ROUTINE)
18260 CS%id_run_model = cpu_clock_id('(ML_MEKE run model)', grain=CLOCK_ROUTINE)
18270 CS%id_unpack_tensor = cpu_clock_id('(ML_MEKE unpack tensor )', grain=CLOCK_ROUTINE)
1828
1829 ! Diagnostics for ML_MEKE
1830 CS%id_mke = register_diag_field('ocean_model', 'MEKE_MKE', diag%axesT1, Time, &
18310 'Surface mean (resolved) kinetic energy used in MEKE', 'm2 s-2', conversion=US%L_T_to_m_s**2)
1832 CS%id_slope_z= register_diag_field('ocean_model', 'MEKE_slope_z', diag%axesT1, Time, &
18330 'Vertically averaged isopyncal slope magnitude used in MEKE', 'nondim', conversion=US%Z_to_L)
1834 CS%id_slope_x= register_diag_field('ocean_model', 'MEKE_slope_x', diag%axesCui, Time, &
18350 'Isopycnal slope in the x-direction used in MEKE', 'nondim', conversion=US%Z_to_L)
1836 CS%id_slope_y= register_diag_field('ocean_model', 'MEKE_slope_y', diag%axesCvi, Time, &
18370 'Isopycnal slope in the y-direction used in MEKE', 'nondim', conversion=US%Z_to_L)
1838 CS%id_rv = register_diag_field('ocean_model', 'MEKE_RV', diag%axesT1, Time, &
18390 'Surface relative vorticity used in MEKE', 's-1', conversion=US%s_to_T)
1840
1841 ! Isoneutral blocking parameters
1842 call get_param(param_file, mdl, "ISOPYCNAL_NIBLOCK", CS%niblock, &
1843 "The i-direction block size used to calculate isopycnal slopes. "//&
1844 "If 0, or when running with OpenMP offload, "//&
1845 "the full computational domain width is used. "//&
1846 "If USE_STANLEY_ISO is true, ISOPYCNAL_NIBLOCK cannot equal 1.", &
18470 default=default_niblock, layoutParam=.true.)
1848 call get_param(param_file, mdl, "ISOPYCNAL_NJBLOCK", CS%njblock, &
1849 "The j-direction block size used to calculate isopycnal slopes. "//&
1850 "If 0, defaults to 1, except when running with OpenMP offload, "//&
1851 "in which case the full computational domain height is used. " //&
1852 "If USE_STANLEY_ISO is true, ISOPYCNAL_NJBLOCK cannot equal 1.", &
18530 default=default_njblock, layoutParam=.true.)
1854 call get_param(param_file, mdl, "ISOPYCNAL_NKBLOCK", CS%nkblock, &
1855 "The k-direction block size used to calculate isopycnal slopes. "//&
1856 "If 0, defaults to 1, except when "//&
1857 "running with OpenMP offload, in which case the full computational "//&
18580 "domain depth is used.", default=default_nkblock, layoutParam=.true.)
1859
18600 if (CS%niblock < 0) &
1861 call MOM_error(FATAL, "ISOPYCNAL_NIBLOCK must be nonnegative; "//&
18620 "use 0 to select the default block size.")
18630 if (CS%njblock < 0) &
1864 call MOM_error(FATAL, "ISOPYCNAL_NJBLOCK must be nonnegative; "//&
18650 "use 0 to select the default block size.")
18660 if (CS%nkblock < 0) &
1867 call MOM_error(FATAL, "ISOPYCNAL_NKBLOCK must be nonnegative; "//&
18680 "use 0 to select the default block size.")
1869
1870
18710end subroutine ML_MEKE_init
1872
1873!> Calculate the various features used for the machine learning prediction
18740subroutine ML_MEKE_calculate_features(G, GV, US, CS, Rd_dx_h, u, v, tv, h, dt, features_array)
1875 type(ocean_grid_type), intent(inout) :: G !< Ocean grid
1876 type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1877 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1878 type(MEKE_CS), intent(in) :: CS !< Control structure for MEKE
1879 real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: Rd_dx_h !< Rossby radius of deformation over
1880 !! the grid length scale [nondim]
1881 real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
1882 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
1883 type(thermo_var_ptrs), intent(in) :: tv !< Type containing thermodynamic variables
1884 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
1885 real, intent(in) :: dt !< Model(baroclinic) time-step [T ~> s].
1886 real(kind=real32), dimension(SIZE(h),num_features), intent( out) :: features_array
1887 !< The array of features needed for machine
1888 !! learning inference, with different units
1889 !! for the various subarrays [various]
1890
18910 real, dimension(SZI_(G),SZJ_(G)) :: mke ! Surface kinetic energy per unit mass [L2 T-2 ~> m2 s-2]
18920 real, dimension(SZI_(G),SZJ_(G)) :: slope_z ! Vertically averaged isoneutral slopes [Z L-1 ~> nondim]
18930 real, dimension(SZIB_(G),SZJB_(G)) :: rv_z ! Surface relative vorticity [T-1 ~> s-1]
18940 real, dimension(SZIB_(G),SZJB_(G)) :: rv_z_t ! Surface relative vorticity interpolated to tracer points [T-1 ~> s-1]
1895
18960 real, dimension(SZIB_(G),SZJ_(G), SZK_(G)) :: h_u ! Thickness at u point [H ~> m or kg m-2]
18970 real, dimension(SZI_(G),SZJB_(G), SZK_(G)) :: h_v ! Thickness at v point [H ~> m or kg m-2]
18980 real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: slope_x ! Isoneutral slope at U point [Z L-1 ~> nondim]
18990 real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: slope_y ! Isoneutral slope at V point [Z L-1 ~> nondim]
19000 real, dimension(SZIB_(G),SZJ_(G)) :: slope_x_vert_avg ! Isoneutral slope at U point [Z L-1 ~> nondim]
19010 real, dimension(SZI_(G),SZJB_(G)) :: slope_y_vert_avg ! Isoneutral slope at V point [Z L-1 ~> nondim]
19020 real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: e ! The interface heights relative to mean sea level [Z ~> m].
1903 real :: slope_t ! Slope interpolated to thickness points [Z L-1 ~> nondim]
1904 real :: u_t, v_t ! u and v interpolated to thickness points [L T-1 ~> m s-1]
1905 real :: dvdx, dudy ! Components of relative vorticity [T-1 ~> s-1]
1906 real :: a_e, a_w, a_n, a_s ! Fractional areas of neighboring cells for interpolating velocities [nondim]
1907 real :: Idenom ! A normalizing factor in calculating weighted averages of areas [L-2 ~> m-2]
1908 real :: sum_area ! A sum of adjacent cell areas [L2 ~> m2]
1909
1910 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1911 integer :: niblock, njblock, nkblock
1912
19130 niblock = CS%niblock
19140 njblock = CS%njblock
19150 nkblock = CS%nkblock
1916
19170 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
19180 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1919
1920 ! Calculate various features for used to infer eddy kinetic energy
1921 ! Linear interpolation to estimate thickness at a velocity points
19220 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
19230 h_u(I,j,k) = 0.5*(h(i,j,k)*G%mask2dT(i,j) + h(i+1,j,k)*G%mask2dT(i+1,j)) + GV%Angstrom_H
19240 h_v(i,J,k) = 0.5*(h(i,j,k)*G%mask2dT(i,j) + h(i,j+1,k)*G%mask2dT(i,j+1)) + GV%Angstrom_H
1925 enddo ; enddo ; enddo
1926
19270 if (niblock == 0) niblock = ie - is + 1
19280 if (njblock == 0) njblock = je - js + 1
19290 if (nkblock == 0) nkblock = nz
1930
1931 !$omp target update to(h)
1932 !$omp target enter data map(alloc: e)
19330 call find_eta(h, tv, G, GV, US, e, halo_size=2)
1934 ! Note the hard-coded dimenisional constant in the following line.
1935 ! UMW: Below is untested
1936 !$omp target enter data map(to: tv, tv%T, tv%S, slope_x, slope_y)
1937 !$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
1938 !$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
1939 call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7*GV%m2_s_to_HZ_T, .false., slope_x, slope_y, &
19400 niblock, njblock, nkblock )
1941 !$omp target exit data map(release: tv, tv%T, tv%S)
1942 !$omp target exit data map(delete: e)
1943 !$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
1944 !$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
1945 !$omp target exit data map(from: slope_x, slope_y)
19460 call pass_vector(slope_x, slope_y, G%Domain)
19470 do j=js-1,je+1 ; do i=is-1,ie+1
19480 slope_x_vert_avg(I,j) = vertical_average_interface(slope_x(i,j,:), h_u(i,j,:), GV%H_subroundoff)
19490 slope_y_vert_avg(i,J) = vertical_average_interface(slope_y(i,j,:), h_v(i,j,:), GV%H_subroundoff)
1950 enddo ; enddo
19510 slope_z(:,:) = 0.
1952
19530 call pass_vector(slope_x_vert_avg, slope_y_vert_avg, G%Domain)
19540 do j=js,je ; do i=is,ie
1955 ! Calculate weights for interpolation from velocity points to h points
19560 sum_area = G%areaCu(I-1,j) + G%areaCu(I,j)
19570 if (sum_area>0.0) then
19580 Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
19590 a_w = G%areaCu(I-1,j) * Idenom
19600 a_e = G%areaCu(I,j) * Idenom
1961 else
19620 a_w = 0.0 ; a_e = 0.0
1963 endif
1964
19650 sum_area = G%areaCv(i,J-1) + G%areaCv(i,J)
19660 if (sum_area>0.0) then
19670 Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
19680 a_s = G%areaCv(i,J-1) * Idenom
19690 a_n = G%areaCv(i,J) * Idenom
1970 else
19710 a_s = 0.0 ; a_n = 0.0
1972 endif
1973
1974 ! Calculate mean kinetic energy
19750 u_t = (a_e*u(I,j,1)) + (a_w*u(I-1,j,1))
19760 v_t = (a_n*v(i,J,1)) + (a_s*v(i,J-1,1))
19770 mke(i,j) = 0.5*( (u_t*u_t) + (v_t*v_t) )
1978
1979 ! Calculate the magnitude of the slope
19800 slope_t = slope_x_vert_avg(I,j)*a_e+slope_x_vert_avg(I-1,j)*a_w
19810 slope_z(i,j) = sqrt(slope_t*slope_t)
19820 slope_t = slope_y_vert_avg(i,J)*a_n+slope_y_vert_avg(i,J-1)*a_s
19830 slope_z(i,j) = 0.5*(slope_z(i,j) + sqrt(slope_t*slope_t))*G%mask2dT(i,j)
1984 enddo ; enddo
19850 call pass_var(slope_z, G%Domain)
1986
1987 ! Calculate relative vorticity
19880 do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
19890 dvdx = ((v(i+1,J,1)*G%dyCv(i+1,J)) - (v(i,J,1)*G%dyCv(i,J)))
19900 dudy = ((u(I,j+1,1)*G%dxCu(I,j+1)) - (u(I,j,1)*G%dxCu(I,j)))
1991 ! Assumed no slip
19920 rv_z(I,J) = (2.0-G%mask2dBu(I,J)) * (dvdx - dudy) * G%IareaBu(I,J)
1993 enddo ; enddo
1994 ! Interpolate RV to t-point, revisit this calculation to include metrics
19950 do j=js,je ; do i=is,ie
19960 rv_z_t(i,j) = 0.25*(rv_z(i-1,j) + rv_z(i,j) + rv_z(i-1,j-1) + rv_z(i,j-1))
1997 enddo ; enddo
1998
1999
2000 ! Construct the feature array
20010 features_array(:,mke_idx) = pack(mke,.true.)
20020 features_array(:,slope_z_idx) = pack(slope_z,.true.)
20030 features_array(:,rd_dx_z_idx) = pack(Rd_dx_h,.true.)
20040 features_array(:,rv_idx) = pack(rv_z_t,.true.)
2005
20060 if (CS%id_rv>0) call post_data(CS%id_rv, rv_z, CS%diag)
20070 if (CS%id_mke>0) call post_data(CS%id_mke, mke, CS%diag)
20080 if (CS%id_slope_z>0) call post_data(CS%id_slope_z, slope_z, CS%diag)
20090 if (CS%id_slope_x>0) call post_data(CS%id_slope_x, slope_x, CS%diag)
20100 if (CS%id_slope_y>0) call post_data(CS%id_slope_y, slope_y, CS%diag)
20110end subroutine ML_MEKE_calculate_features
2012
2013!> Use the machine learning interface to predict EKE
20140subroutine predict_MEKE(G, US, CS, npts, Time, features_array, MEKE)
2015 type(ocean_grid_type), intent(inout) :: G !< Ocean grid
2016 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2017 type(MEKE_CS), intent(in ) :: CS !< Control structure for MEKE
2018 integer, intent(in ) :: npts !< Number of T-grid cells on the local
2019 !! domain
2020 type(time_type), intent(in ) :: Time !< The current model time
2021 real(kind=real32), dimension(npts,num_features), intent(in ) :: features_array
2022 !< The array of features needed for machine
2023 !! learning inference, with different units
2024 !! for the various subarrays [various]
2025 real, dimension(SZI_(G),SZJ_(G)), intent( out) :: MEKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2]
2026
2027 ! Local variables
2028 integer :: db_return_code
2029 character(len=255), dimension(1) :: model_out, model_in
2030 character(len=255) :: time_suffix
20310 real(kind=real32), dimension(SIZE(MEKE)) :: MEKE_vec ! A one-dimensional array of the natural log of eddy kinetic
2032 ! energy in mks units [m2 s-2]
20330 real, dimension(size(MEKE,1),size(MEKE,2)) :: ln_MEKE ! the natural log of eddy kinetic energy
2034 ! in mks units [m2 s-2]
20350 real, dimension(size(MEKE,1),size(MEKE,2)) :: MEKE_mks ! The eddy kinetic energy in mks units [m2 s-2]
2036 integer :: i, j, is, ie, js, je
2037
20380 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
2039!> Use the database client to call a machine learning model to predict eddy kinetic energy
20400 call cpu_clock_begin(CS%id_put_tensor)
20410 db_return_code = CS%client%put_tensor("features"//CS%key_suffix, features_array, shape(features_array))
20420 call cpu_clock_end(CS%id_put_tensor)
2043
2044 ! Run the ML model to predict EKE and return the result
20450 model_out(1) = "EKE"//CS%key_suffix
20460 model_in(1) = "features"//CS%key_suffix
20470 call cpu_clock_begin(CS%id_run_model)
20480 db_return_code = CS%client%run_model(CS%model_key, model_in, model_out)
20490 call cpu_clock_end(CS%id_run_model)
20500 if (CS%client%SR_error_parser(db_return_code)) then
20510 call MOM_error(FATAL, "MEKE: run_model failed")
2052 endif
20530 call cpu_clock_begin(CS%id_unpack_tensor)
20540 db_return_code = CS%client%unpack_tensor( model_out(1), MEKE_vec, shape(MEKE_vec) )
20550 call cpu_clock_end(CS%id_unpack_tensor)
2056
20570 ln_MEKE = reshape(MEKE_vec, shape(MEKE))
2058 ! Zero out the halos. These will usually be reset by the pass_var in a few lines.
20590 MEKE_mks(:,:) = 0.0
20600 do j=js,je ; do i=is,ie
20610 MEKE_mks(i,j) = MIN(exp(ln_MEKE(i,j)), US%L_T_to_m_s**2*CS%eke_max)
2062 enddo ; enddo
20630 call pass_var(MEKE_mks, G%Domain, halo=1)
2064
20650 if (CS%online_analysis) then
20660 write(time_suffix,"(F16.0)") time_type_to_real(Time)
2067 db_return_code = CS%client%put_tensor(trim("EKE_")//trim(adjustl(time_suffix))//CS%key_suffix, &
20680 MEKE_mks, shape(MEKE))
2069 endif
2070
2071 ! Copy MEKE_mks into the argument in rescaled units.
2072 ! MEKE(:,:) = 0.0 ! This would fill in the wider halos of this intent(out) array.
20730 do j=js-1,je+1 ; do i=is-1,ie+1
20740 MEKE(i,j) = US%m_s_to_L_T**2 * MEKE_mks(i,j)
2075 enddo ; enddo
2076
20770end subroutine predict_MEKE
2078
2079!> Compute average of interface quantities weighted by the thickness of the surrounding
2080!! layers [arbitrary]
20810real function vertical_average_interface(h, w, h_min)
2082
2083 real, dimension(:), intent(in) :: h !< Layer Thicknesses [H ~> m or kg m-2]
2084 real, dimension(:), intent(in) :: w !< Quantity to average [arbitrary]
2085 real, intent(in) :: h_min !< The vanishingly small layer thickness [H ~> m or kg m-2]
2086
2087 real :: htot ! Twice the sum of the layer thicknesses interpolated to interior interfaces [H ~> m or kg m-2]
2088 real :: inv_htot ! The inverse of htot [H-1 ~> m-1 or m2 kg-1]
2089 integer :: k, nk
2090
20910 nk = size(h)
20920 htot = h_min
20930 do k=2,nk
20940 htot = htot + (h(k-1)+h(k))
2095 enddo
20960 inv_htot = 1./htot
2097
20980 vertical_average_interface = 0.
20990 do K=2,nk
21000 vertical_average_interface = vertical_average_interface + (w(k)*(h(k-1)+h(k)))*inv_htot
2101 enddo
21020end function vertical_average_interface
2103
2104!> Allocates memory and register restart fields for the MOM_MEKE module.
21051subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
2106! Arguments
2107 type(hor_index_type), intent(in) :: HI !< Horizontal index structure
2108 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2109 type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
2110 type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
2111 type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct
2112
2113 ! Local variables
2114 real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_bhFrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim]
2115 real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim]
2116 logical :: Use_KH_in_MEKE
2117 logical :: useMEKE
2118 logical :: sqg_use_MEKE
2119 integer :: isd, ied, jsd, jed
2120
2121! Determine whether this module will be used
21221 useMEKE = .false. ; call read_param(param_file,"USE_MEKE",useMEKE)
2123
2124! Read these parameters to determine what should be in the restarts
21251 MEKE_GMcoeff = -1. ; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff)
21261 MEKE_FrCoeff = -1. ; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff)
21271 MEKE_bhFrCoeff = -1. ; call read_param(param_file,"MEKE_bhFRCOEFF",MEKE_bhFrCoeff)
21281 MEKE_GMEcoeff = -1. ; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff)
21291 MEKE_KhCoeff = 1. ; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff)
21301 MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku)
21311 MEKE_viscCoeff_Au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au)
21321 Use_KH_in_MEKE = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE)
21331 sqg_use_MEKE = .false. ; call read_param(param_file,"SQG_USE_MEKE", sqg_use_MEKE)
2134
21351 if (.not. useMEKE) return
2136
2137! Allocate memory
21380 call MOM_mesg("MEKE_alloc_register_restart: allocating and registering", 5)
21390 isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed
21400 allocate(MEKE%MEKE(isd:ied,jsd:jed), source=0.0)
2141 call register_restart_field(MEKE%MEKE, "MEKE", .false., restart_CS, &
21420 longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2)
2143
21440 if (MEKE_GMcoeff>=0.) allocate(MEKE%GM_src(isd:ied,jsd:jed), source=0.0)
21450 if (MEKE_FrCoeff>=0. .or. MEKE_bhFrCoeff>=0. .or. MEKE_GMECoeff>=0.) &
21460 allocate(MEKE%mom_src(isd:ied,jsd:jed), source=0.0)
21470 if (MEKE_bhFrCoeff >= 0.) &
21480 allocate(MEKE%mom_src_bh(isd:ied,jsd:jed), source=0.0)
21490 if (MEKE_FrCoeff<0.) MEKE_FrCoeff = 0.
21500 if (MEKE_bhFrCoeff<0.) MEKE_bhFrCoeff = 0.
21510 if (MEKE_GMECoeff>=0.) allocate(MEKE%GME_snk(isd:ied,jsd:jed), source=0.0)
21520 if (MEKE_KhCoeff>=0.) then
21530 allocate(MEKE%Kh(isd:ied,jsd:jed), source=0.0)
2154 call register_restart_field(MEKE%Kh, "MEKE_Kh", .false., restart_CS, &
2155 longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy", &
21560 units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T)
2157 endif
21580 allocate(MEKE%Rd_dx_h(isd:ied,jsd:jed), source=0.0)
21590 if (MEKE_viscCoeff_Ku/=0.) then
21600 allocate(MEKE%Ku(isd:ied,jsd:jed), source=0.0)
2161 call register_restart_field(MEKE%Ku, "MEKE_Ku", .false., restart_CS, &
2162 longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", &
21630 units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T)
2164 endif
21650 if (sqg_use_MEKE) then
21660 allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0)
2167 call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, &
2168 longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", &
21690 units="m", conversion=US%L_to_m)
2170 endif
21710 if (Use_Kh_in_MEKE) then
21720 allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0)
2173 call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, &
2174 longname="Copy of thickness diffusivity for diffusing MEKE", &
21750 units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T)
2176 endif
2177
21780 if (MEKE_viscCoeff_Au/=0.) then
21790 allocate(MEKE%Au(isd:ied,jsd:jed), source=0.0)
2180 call register_restart_field(MEKE%Au, "MEKE_Au", .false., restart_CS, &
2181 longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy", &
21820 units="m4 s-1", conversion=US%L_to_m**4*US%s_to_T)
2183 endif
2184
2185end subroutine MEKE_alloc_register_restart
2186
2187!> Deallocates any variables allocated in MEKE_alloc_register_restart.
21881subroutine MEKE_end(MEKE)
2189 type(MEKE_type), intent(inout) :: MEKE !< A structure with MEKE-related fields.
2190
2191 ! NOTE: MEKE will always be allocated by MEKE_init, even if MEKE is disabled.
2192 ! So these must all be conditional, even though MEKE%MEKE and MEKE%Rd_dx_h
2193 ! are always allocated (when MEKE is enabled)
2194
21951 if (allocated(MEKE%Au)) deallocate(MEKE%Au)
21961 if (allocated(MEKE%Kh_diff)) deallocate(MEKE%Kh_diff)
21971 if (allocated(MEKE%Ku)) deallocate(MEKE%Ku)
21981 if (allocated(MEKE%Rd_dx_h)) deallocate(MEKE%Rd_dx_h)
21991 if (allocated(MEKE%Kh)) deallocate(MEKE%Kh)
22001 if (allocated(MEKE%GME_snk)) deallocate(MEKE%GME_snk)
22011 if (allocated(MEKE%mom_src)) deallocate(MEKE%mom_src)
22021 if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh)
22031 if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src)
22041 if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE)
22051 if (allocated(MEKE%Le)) deallocate(MEKE%Le)
22061end subroutine MEKE_end
2207
2208!> \namespace mom_meke
2209!!
2210!! \section section_MEKE The Mesoscale Eddy Kinetic Energy (MEKE) framework
2211!!
2212!! The MEKE framework accounts for the mean potential energy removed by
2213!! the first order closures used to parameterize mesoscale eddies.
2214!! It requires closure at the second order, namely dissipation and transport
2215!! of eddy energy.
2216!!
2217!! Monitoring the sub-grid scale eddy energy budget provides a means to predict
2218!! a sub-grid eddy-velocity scale which can be used in the lower order closures.
2219!!
2220!! \subsection section_MEKE_equations MEKE equations
2221!!
2222!! The eddy kinetic energy equation is:
2223!! \f[ \partial_{\tilde{t}} E =
2224!! \overbrace{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v
2225!! }^\text{sources}
2226!! - \overbrace{ ( \lambda + C_d | U_d | \gamma_b^2 ) E
2227!! }^\text{local dissipation}
2228!! + \overbrace{ \nabla \cdot ( ( \kappa_E + \gamma_M \kappa_M ) \nabla E
2229!! - \kappa_4 \nabla^3 E )
2230!! }^\text{smoothing}
2231!! \f]
2232!! where \f$ E \f$ is the eddy kinetic energy (variable <code>MEKE</code>) with units of
2233!! m<sup>2</sup>s<sup>-2</sup>,
2234!! and \f$\tilde{t} = a t\f$ is a scaled time. The non-dimensional factor
2235!! \f$ a\geq 1 \f$ is used to accelerate towards equilibrium.
2236!!
2237!! The MEKE equation is two-dimensional and obtained by depth averaging the
2238!! the three-dimensional eddy energy equation. In the following expressions
2239!! \f$ \left< \phi \right> = \frac{1}{H} \int^\eta_{-D} \phi \, dz \f$ maps
2240!! three dimensional terms into the two-dimensional quantities needed.
2241!!
2242!! \subsubsection section_MEKE_source_terms MEKE source terms
2243!!
2244!! The source term \f$ \dot{E}_b \f$ is a constant background source
2245!! of energy intended to avoid the limit \f$E\rightarrow 0\f$.
2246!!
2247!! The "GM" source term
2248!! \f[ \dot{E}_\eta = - \left< \overline{w^\prime b^\prime} \right>
2249!! = \left< \kappa_h N^2S^2 \right>
2250!! \approx \left< \kappa_h g\prime |\nabla_\sigma \eta|^2 \right>\f]
2251!! equals the mean potential energy removed by the Gent-McWilliams closure,
2252!! and is excluded/included in the MEKE budget by the efficiency parameter
2253!! \f$ \gamma_\eta \in [0,1] \f$.
2254!!
2255!! The "frictional" source term
2256!! \f[ \dot{E}_{v} = \left< \partial_i u_j \tau_{ij} \right> \f]
2257!! equals the mean kinetic energy removed by lateral viscous fluxes, and
2258!! is excluded/included in the MEKE budget by the efficiency parameter
2259!! \f$ \gamma_v \in [0,1] \f$.
2260!!
2261!! \subsubsection section_MEKE_dissipation_terms MEKE dissipation terms
2262!!
2263!! The local dissipation of \f$ E \f$ is parameterized through a linear
2264!! damping, \f$\lambda\f$, and bottom drag, \f$ C_d | U_d | \gamma_b^2 \f$.
2265!! The \f$ \gamma_b \f$ accounts for the weak projection of the column-mean
2266!! eddy velocity to the bottom. In other words, the bottom velocity is
2267!! estimated as \f$ \gamma_b U_e \f$.
2268!! The bottom drag coefficient, \f$ C_d \f$ is the same as that used in the bottom
2269!! friction in the mean model equations.
2270!!
2271!! The bottom drag velocity scale, \f$ U_d \f$, has contributions from the
2272!! resolved state and \f$ E \f$:
2273!! \f[ U_d = \sqrt{ U_b^2 + |u|^2_{z=-D} + |\gamma_b U_e|^2 } .\f]
2274!! where the eddy velocity scale, \f$ U_e \f$, is given by:
2275!! \f[ U_e = \sqrt{ 2 E } .\f]
2276!! \f$ U_b \f$ is a constant background bottom velocity scale and is
2277!! typically not used (i.e. set to zero).
2278!!
2279!! Following \cite jansen2015, the projection of eddy energy on to the bottom
2280!! is given by the ratio of bottom energy to column mean energy:
2281!! \f[
2282!! \gamma_b^2 = \frac{E_b}{E} = \gamma_{d0}
2283!! + \left( 1 + c_{b} \frac{L_d}{L_f} \right)^{-\frac{4}{5}}
2284!! ,
2285!! \f]
2286!! \f[
2287!! \gamma_b^2 \leftarrow \max{\left( \gamma_b^2, \gamma_{min}^2 \right)}
2288!! .
2289!! \f]
2290!!
2291!! \subsection section_MEKE_smoothing MEKE smoothing terms
2292!!
2293!! \f$ E \f$ is laterally diffused by a diffusivity \f$ \kappa_E + \gamma_M
2294!! \kappa_M \f$ where \f$ \kappa_E \f$ is a constant diffusivity and the term
2295!! \f$ \gamma_M \kappa_M \f$ is a "self diffusion" using the diffusivity
2296!! calculated in the section \ref section_MEKE_diffusivity.
2297!! \f$ \kappa_4 \f$ is a constant bi-harmonic diffusivity.
2298!!
2299!! \subsection section_MEKE_diffusivity Diffusivity derived from MEKE
2300!!
2301!! The predicted eddy velocity scale, \f$ U_e \f$, can be combined with a
2302!! mixing length scale to form a diffusivity.
2303!! The primary use of a MEKE derived diffusivity is for use in thickness
2304!! diffusion (module mom_thickness_diffuse) and optionally in along
2305!! isopycnal mixing of tracers (module mom_tracer_hor_diff).
2306!! The original form used (enabled with MEKE_OLD_LSCALE=True):
2307!!
2308!! \f[ \kappa_M = \gamma_\kappa \sqrt{ \gamma_t^2 U_e^2 A_\Delta } \f]
2309!!
2310!! where \f$ A_\Delta \f$ is the area of the grid cell.
2311!! Following \cite jansen2015, we now use
2312!!
2313!! \f[ \kappa_M = \gamma_\kappa l_M \sqrt{ \gamma_t^2 U_e^2 } \f]
2314!!
2315!! where \f$ \gamma_\kappa \in [0,1] \f$ is a non-dimensional factor and,
2316!! following \cite jansen2015, \f$\gamma_t^2\f$ is the ratio of barotropic
2317!! eddy energy to column mean eddy energy given by
2318!! \f[
2319!! \gamma_t^2 = \frac{E_t}{E} = \left( 1 + c_{t} \frac{L_d}{L_f} \right)^{-\frac{1}{4}}
2320!! ,
2321!! \f]
2322!! \f[
2323!! \gamma_t^2 \leftarrow \max{\left( \gamma_t^2, \gamma_{min}^2 \right)}
2324!! .
2325!! \f]
2326!!
2327!! The length-scale is a configurable combination of multiple length scales:
2328!!
2329!! \f[
2330!! l_M = \left(
2331!! \frac{\alpha_d}{L_d}
2332!! + \frac{\alpha_f}{L_f}
2333!! + \frac{\alpha_R}{L_R}
2334!! + \frac{\alpha_e}{L_e}
2335!! + \frac{\alpha_\Delta}{L_\Delta}
2336!! + \frac{\delta[L_c]}{L_c}
2337!! \right)^{-1}
2338!! \f]
2339!!
2340!! where
2341!!
2342!! \f{eqnarray*}{
2343!! L_d & = & \sqrt{\frac{c_g^2}{f^2+2\beta c_g}} \sim \frac{ c_g }{f} \\\\
2344!! L_R & = & \sqrt{\frac{U_e}{\beta^*}} \\\\
2345!! L_e & = & \frac{U_e}{|S| N} \\\\
2346!! L_f & = & \frac{H}{c_d} \\\\
2347!! L_\Delta & = & \sqrt{A_\Delta} .
2348!! \f}
2349!!
2350!! \f$L_c\f$ is a constant and \f$\delta[L_c]\f$ is the impulse function so that the term
2351!! \f$\frac{\delta[L_c]}{L_c}\f$ evaluates to \f$\frac{1}{L_c}\f$ when \f$L_c\f$ is non-zero
2352!! but is dropped if \f$L_c=0\f$.
2353!!
2354!! \f$\beta^*\f$ is the effective \f$\beta\f$ that combines both the planetary vorticity
2355!! gradient (i.e. \f$\beta=\nabla f\f$) and the topographic \f$\beta\f$ effect,
2356!! with the latter weighed by a weighting constant, \f$c_\beta\f$, that varies
2357!! from 0 to 1, so that \f$c_\beta=0\f$ means the topographic \f$\beta\f$ effect is ignored,
2358!! while \f$c_\beta=1\f$ means it is fully considered. The new \f$\beta^*\f$ therefore
2359!! takes the form of
2360!!
2361!! \f[
2362!! \beta^* = \sqrt{( \partial_xf - c_\beta\frac{f}{D}\partial_xD )^2 +
2363!! ( \partial_yf - c_\beta\frac{f}{D}\partial_yD )^2}
2364!! \f]
2365!! where \f$D\f$ is water column depth at T points.
2366!!
2367!! \subsection section_MEKE_viscosity Viscosity derived from MEKE
2368!!
2369!! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be
2370!! used to form a harmonic eddy viscosity,
2371!!
2372!! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } \f]
2373!!
2374!! as well as a biharmonic eddy viscosity,
2375!!
2376!! \f[ \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 } \f]
2377!!
2378!! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance
2379!!
2380!! Note that in steady-state (or when \f$ a>>1 \f$) and there is no
2381!! diffusion of \f$ E \f$ then
2382!! \f[ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
2383!! \gamma_v \dot{E}_v }{ \lambda + C_d|U_d|\gamma_b^2 } . \f]
2384!!
2385!! In the linear drag limit, where
2386!! \f$ U_e << \min(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \f$, the equilibrium becomes
2387!! \f$ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
2388!! \gamma_v \dot{E}_v }{ \lambda + C_d \sqrt{ U_b^2 + |u|^2_{z=-D} } } \f$.
2389!!
2390!! In the nonlinear drag limit, where \f$ U_e >> \max(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \f$,
2391!! the equilibrium becomes
2392!! \f$ \overline{E} \approx \left( \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
2393!! \gamma_v \dot{E}_v }{ \sqrt{2} C_d \gamma_b^3 } \right)^\frac{2}{3} \f$.
2394!!
2395!! \subsubsection section_MEKE_module_parameters MEKE module parameters
2396!!
2397!! | Symbol | Module parameter |
2398!! | ------ | --------------- |
2399!! | - | <code>USE_MEKE</code> |
2400!! | \f$ a \f$ | <code>MEKE_DTSCALE</code> |
2401!! | \f$ \dot{E}_b \f$ | <code>MEKE_BGSRC</code> |
2402!! | \f$ \gamma_\eta \f$ | <code>MEKE_GMCOEFF</code> |
2403!! | \f$ \gamma_v \f$ | <code>MEKE_FrCOEFF</code> |
2404!! | \f$ \lambda \f$ | <code>MEKE_DAMPING</code> |
2405!! | \f$ U_b \f$ | <code>MEKE_USCALE</code> |
2406!! | \f$ \gamma_{d0} \f$ | <code>MEKE_CD_SCALE</code> |
2407!! | \f$ c_{b} \f$ | <code>MEKE_CB</code> |
2408!! | \f$ c_{t} \f$ | <code>MEKE_CT</code> |
2409!! | \f$ \kappa_E \f$ | <code>MEKE_KH</code> |
2410!! | \f$ \kappa_4 \f$ | <code>MEKE_K4</code> |
2411!! | \f$ \gamma_\kappa \f$ | <code>MEKE_KHCOEFF</code> |
2412!! | \f$ \gamma_M \f$ | <code>MEKE_KHMEKE_FAC</code> |
2413!! | \f$ \gamma_u \f$ | <code>MEKE_VISCOSITY_COEFF_KU</code> |
2414!! | \f$ \gamma_4 \f$ | <code>MEKE_VISCOSITY_COEFF_AU</code> |
2415!! | \f$ \gamma_{min}^2 \f$| <code>MEKE_MIN_GAMMA2</code> |
2416!! | \f$ \alpha_d \f$ | <code>MEKE_ALPHA_DEFORM</code> |
2417!! | \f$ \alpha_f \f$ | <code>MEKE_ALPHA_FRICT</code> |
2418!! | \f$ \alpha_R \f$ | <code>MEKE_ALPHA_RHINES</code> |
2419!! | \f$ \alpha_e \f$ | <code>MEKE_ALPHA_EADY</code> |
2420!! | \f$ \alpha_\Delta \f$ | <code>MEKE_ALPHA_GRID</code> |
2421!! | \f$ L_c \f$ | <code>MEKE_FIXED_MIXING_LENGTH</code> |
2422!! | \f$ c_\beta \f$ | <code>MEKE_TOPOGRAPHIC_BETA</code> |
2423!! | - | <code>MEKE_KHTH_FAC</code> |
2424!! | - | <code>MEKE_KHTR_FAC</code> |
2425!!
2426!! | Symbol | Model parameter |
2427!! | ------ | --------------- |
2428!! | \f$ C_d \f$ | <code>CDRAG</code> |
2429!!
2430!! \subsection section_MEKE_references References
2431!!
2432!! Jansen, M. F., A. J. Adcroft, R. Hallberg, and I. M. Held, 2015: Parameterization of eddy fluxes based on a
2433!! mesoscale energy budget. Ocean Modelling, 92, 28--41, http://doi.org/10.1016/j.ocemod.2015.05.007 .
2434!!
2435!! Marshall, D. P., and A. J. Adcroft, 2010: Parameterization of ocean eddies: Potential vorticity mixing, energetics
2436!! and Arnold first stability theorem. Ocean Modelling, 32, 188--204, http://doi.org/10.1016/j.ocemod.2010.02.001 .
2437
24380end module MOM_MEKE
2439