MOM_set_viscosity.F90

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!> Calculates various values related to the bottom boundary layer, such as the viscosity and
6!! thickness of the BBL (set_viscous_BBL).
7module mom_set_visc
8
10use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
14use mom_debugging, only : uvchksum, hchksum
15use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
16use mom_diag_mediator, only : diag_ctrl, time_type
17use mom_domains, only : pass_var, corner
18use mom_eos, only : calculate_density, calculate_density_derivs, calculate_specific_vol_derivs
19use mom_error_handler, only : mom_error, fatal, warning
20use mom_file_parser, only : get_param, log_param, log_version, param_file_type
22use mom_forcing_type, only : forcing, mech_forcing, find_ustar
23use mom_grid, only : ocean_grid_type
25use mom_interface_heights, only : thickness_to_dz
27use mom_io, only : slasher, mom_read_data, vardesc, var_desc
29use mom_open_boundary, only : ocean_obc_type, obc_segment_type, obc_none, obc_direction_e
30use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s
31use mom_restart, only : register_restart_field, query_initialized, mom_restart_cs
32use mom_restart, only : register_restart_field_as_obsolete, register_restart_pair
33use mom_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc
37
38implicit none ; private
39
40#include <MOM_memory.h>
41
45
46! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
47! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
48! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
49! vary with the Boussinesq approximation, the Boussinesq variant is given first.
50
51!> Control structure for MOM_set_visc
52type, public :: set_visc_cs ; private
53 logical :: initialized = .false. !< True if this control structure has been initialized.
54 real :: hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2].
55 !! Runtime parameter `HBBL`.
56 real :: dz_bbl !< The static bottom boundary layer thickness in height units [Z ~> m].
57 !! Runtime parameter `HBBL`.
58 real :: cdrag !< The quadratic drag coefficient [nondim].
59 !! Runtime parameter `CDRAG`.
60 real :: c_smag !< The Laplacian Smagorinsky coefficient for
61 !! calculating the drag in channels [nondim].
62 real :: drag_bg_vel !< An assumed unresolved background velocity for
63 !! calculating the bottom drag [L T-1 ~> m s-1].
64 !! Runtime parameter `DRAG_BG_VEL`.
65 !! Should not be used if BBL_USE_TIDAL_BG is True.
66 real :: bbl_thick_min !< The minimum bottom boundary layer thickness [Z ~> m].
67 !! This might be Kv / (cdrag * drag_bg_vel) to give
68 !! Kv as the minimum near-bottom viscosity.
69 real :: htbl_shelf !< A nominal thickness of the surface boundary layer for use
70 !! in calculating the near-surface velocity [H ~> m or kg m-2].
71 real :: htbl_shelf_min !< The minimum surface boundary layer thickness [Z ~> m].
72 real :: kv_bbl_min !< The minimum viscosity in the bottom boundary layer [H Z T-1 ~> m2 s-1 or Pa s]
73 real :: kv_tbl_min !< The minimum viscosity in the top boundary layer [H Z T-1 ~> m2 s-1 or Pa s]
74 logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
75 !! drag law c_drag*|u|*u. The velocity magnitude
76 !! may be an assumed value or it may be based on the
77 !! actual velocity in the bottommost `HBBL`, depending
78 !! on whether linear_drag is true.
79 !! Runtime parameter `BOTTOMDRAGLAW`.
80 logical :: bottomdragmap !< If true, apply the spatially varying drag coefficient (cdrag_2d)
81 !! instead of the spatially uniform drag coefficient (cdrag).
82 logical :: body_force_drag !< If true, the bottom stress is imposed as an explicit body force
83 !! applied over a fixed distance from the bottom, rather than as an
84 !! implicit calculation based on an enhanced near-bottom viscosity.
85 logical :: bbl_use_eos !< If true, use the equation of state in determining
86 !! the properties of the bottom boundary layer.
87 logical :: linear_drag !< If true, the drag law is cdrag*`DRAG_BG_VEL`*u.
88 !! Runtime parameter `LINEAR_DRAG`.
89 logical :: channel_drag !< If true, the drag is exerted directly on each layer
90 !! according to what fraction of the bottom they overlie.
91 real :: chan_drag_max_vol !< The maximum bottom boundary layer volume within which the
92 !! channel drag is applied, normalized by the full cell area,
93 !! or a negative value to apply no maximum [Z ~> m].
94 real :: channel_break_depth !< When CHANNEL_DRAG is true, the bathymetric depth interpolated
95 !! to the vorticity point is a combination of the harmonic mean of the
96 !! adjacent velocity point depths below this depth [Z ~> m] and the
97 !! arithmetic mean of the adjacent depths above it, to roughly mimic a
98 !! continental shelf break profile. The internal version of this depth
99 !! uses the same offset (G%Z_ref) as the bathymetry.
100 logical :: correct_bbl_bounds !< If true, uses the correct bounds on the BBL thickness and
101 !! viscosity so that the bottom layer feels the intended drag.
102 logical :: rino_mix !< If true, use Richardson number dependent mixing.
103 logical :: dynamic_viscous_ml !< If true, use a bulk Richardson number criterion to
104 !! determine the mixed layer thickness for viscosity.
105 real :: bulk_ri_ml !< The bulk mixed layer used to determine the
106 !! thickness of the viscous mixed layer [nondim]
107 real :: omega !< The Earth's rotation rate [T-1 ~> s-1].
108 real :: ustar_min !< A minimum value of ustar to avoid numerical
109 !! problems [H T-1 ~> m s-1 or kg m-2 s-1]. If the value is
110 !! small enough, this should not affect the solution.
111 real :: tke_decay !< The ratio of the natural Ekman depth to the TKE
112 !! decay scale [nondim]
113 real :: omega_frac !< When setting the decay scale for turbulence, use this
114 !! fraction of the absolute rotation rate blended with the local
115 !! value of f, as sqrt((1-of)*f^2 + of*4*omega^2) [nondim]
116 real :: tideampfac2 !< A factor to multiply by tideamp to convert to a mean ustar,
117 !! accounts for conversion of amplitude to mean magnitude over
118 !! a time average much longer than the tidal periods and for
119 !! non-commuting conversion of mean tideamp to mean ustar**3 [nondim]
120 logical :: concave_trigonometric_l !< If true, use trigonometric expressions to determine the
121 !! fractional open interface lengths for concave topography.
122 integer :: answer_date !< The vintage of the order of arithmetic and expressions
123 !! in the set viscosity calculations. Values below
124 !! 20190101 recover the answers from the end of 2018,
125 !! while higher values use updated and more robust forms
126 !! of the same expressions. Values below 20260704 use
127 !! the non-reproducible power operator in place of
128 !! cuberoot() when finding the open distances with
129 !! CHANNEL_DRAG.
130 logical :: debug !< If true, write verbose checksums for debugging purposes.
131 logical :: bbl_use_tidal_bg !< If true, use a tidal background amplitude for the bottom velocity
132 !! when computing the bottom stress.
133 character(len=200) :: inputdir !< The directory for input files.
134 type(ocean_obc_type), pointer :: obc => null() !< Open boundaries control structure
135 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
136 !! regulate the timing of diagnostic output.
137 ! Allocatable data arrays
138 real, allocatable, dimension(:,:) :: cdrag_u !< The spatially varying quadratic drag coefficient [nondim]
139 real, allocatable, dimension(:,:) :: cdrag_v !< The spatially varying quadratic drag coefficient [nondim]
140 real, allocatable, dimension(:,:) :: tideamp !< RMS tidal amplitude at h points [L T-1 ~> m s-1]
141 ! Diagnostic arrays
142 real, allocatable, dimension(:,:) :: bbl_u !< BBL mean U current [L T-1 ~> m s-1]
143 real, allocatable, dimension(:,:) :: bbl_v !< BBL mean V current [L T-1 ~> m s-1]
144 !>@{ Diagnostics handles
145 integer :: id_bbl_thick_u = -1, id_kv_bbl_u = -1, id_bbl_u = -1
146 integer :: id_bbl_thick_v = -1, id_kv_bbl_v = -1, id_bbl_v = -1
147 integer :: id_ray_u = -1, id_ray_v = -1
148 integer :: id_nkml_visc_u = -1, id_nkml_visc_v = -1
149 !>@}
150end type set_visc_cs
151
152contains
153
154!> Calculates the thickness of the bottom boundary layer and the viscosity within that layer.
155subroutine set_viscous_bbl(u, v, h, tv, visc, G, GV, US, CS, pbv)
156 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
157 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
158 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
159 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
160 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
161 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
162 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
163 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
164 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
165 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
166 !! available thermodynamic fields. Absent fields
167 !! have NULL pointers.
168 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
169 !! related fields.
170 type(set_visc_cs), intent(inout) :: cs !< The control structure returned by a previous
171 !! call to set_visc_init.
172 type(porous_barrier_type),intent(in) :: pbv !< porous barrier fractional cell metrics
173
174 ! Local variables
175 real, dimension(SZIB_(G)) :: &
176 ustar, & ! The bottom friction velocity [H T-1 ~> m s-1 or kg m-2 s-1].
177 t_eos, & ! The temperature used to calculate the partial derivatives
178 ! of density with T and S [C ~> degC].
179 s_eos, & ! The salinity used to calculate the partial derivatives
180 ! of density with T and S [S ~> ppt].
181 dr_dt, & ! Partial derivative of the density in the bottom boundary
182 ! layer with temperature [R C-1 ~> kg m-3 degC-1].
183 dr_ds, & ! Partial derivative of the density in the bottom boundary
184 ! layer with salinity [R S-1 ~> kg m-3 ppt-1].
185 press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
186 umag_avg, & ! The average magnitude of velocities in the bottom boundary layer [L T-1 ~> m s-1].
187 h_bbl_drag, & ! The thickness over which to apply drag as a body force [H ~> m or kg m-2].
188 dz_bbl_drag ! The vertical height over which to apply drag as a body force [Z ~> m].
189 real :: htot ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
190 real :: dztot ! Distance from the bottom up to some point [Z ~> m].
191 real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
192 real :: dztot_vel ! Distance from the bottom up to some point [Z ~> m].
193
194 real :: rhtot ! Running sum of thicknesses times the layer potential
195 ! densities [H R ~> kg m-2 or kg2 m-5].
196 real, dimension(SZIB_(G),SZJ_(G)) :: &
197 d_u, & ! Bottom depth linearly interpolated to u points [Z ~> m].
198 mask_u ! A mask that disables any contributions from u points that
199 ! are land or past open boundary conditions [nondim], 0 or 1.
200 real, dimension(SZI_(G),SZJB_(G)) :: &
201 d_v, & ! Bottom depth linearly interpolated to v points [Z ~> m].
202 mask_v ! A mask that disables any contributions from v points that
203 ! are land or past open boundary conditions [nondim], 0 or 1.
204 real, dimension(SZIB_(G),SZK_(GV)) :: &
205 h_at_vel, & ! Layer thickness at a velocity point, using an upwind-biased
206 ! second order accurate estimate based on the previous velocity
207 ! direction [H ~> m or kg m-2].
208 h_vel, & ! Arithmetic mean of the layer thicknesses adjacent to a
209 ! velocity point [H ~> m or kg m-2].
210 dz_at_vel, & ! Vertical extent of a layer, using an upwind-biased
211 ! second order accurate estimate based on the previous velocity
212 ! direction [Z ~> m].
213 dz_vel, & ! Arithmetic mean of the difference in across the layers adjacent
214 ! to a velocity point [Z ~> m].
215 t_vel, & ! Arithmetic mean of the layer temperatures adjacent to a
216 ! velocity point [C ~> degC].
217 s_vel, & ! Arithmetic mean of the layer salinities adjacent to a
218 ! velocity point [S ~> ppt].
219 spv_vel, & ! Arithmetic mean of the layer averaged specific volumes adjacent to a
220 ! velocity point [R-1 ~> m3 kg-1].
221 rml_vel ! Arithmetic mean of the layer coordinate densities adjacent
222 ! to a velocity point [R ~> kg m-3].
223 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Height change across layers [Z ~> m]
224
225 real :: h_vel_pos ! The arithmetic mean thickness at a velocity point
226 ! plus H_neglect to avoid 0 values [H ~> m or kg m-2].
227 real :: ustarsq ! 400 times the square of ustar, times
228 ! Rho0 divided by G_Earth and the conversion
229 ! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
230 real :: cdrag ! The drag coefficient [nondim].
231 real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
232 real :: cdrag_sqrt_h ! Square root of the drag coefficient, times a unit conversion factor
233 ! from lateral lengths to layer thicknesses [H L-1 ~> nondim or kg m-3].
234 real :: cdrag_sqrt_h_rl ! Square root of the drag coefficient, times a unit conversion factor from
235 ! density times lateral lengths to layer thicknesses [H L-1 R-1 ~> m3 kg-1 or nondim]
236 real :: cdrag_l_to_h ! The drag coefficient times conversion factors from lateral
237 ! distance to thickness units [H L-1 ~> nondim or kg m-3]
238 real :: cdrag_rl_to_h ! The drag coefficient times conversion factors from density times lateral
239 ! distance to thickness units [H L-1 R-1 ~> m3 kg-1 or nondim]
240 real :: cdrag_conv ! The drag coefficient times a combination of static conversion factors and in
241 ! situ density or Boussinesq reference density [H L-1 ~> nondim or kg m-3]
242 real :: oldfn ! The integrated energy required to
243 ! entrain up to the bottom of the layer,
244 ! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
245 real :: dfn ! The increment in oldfn for entraining
246 ! the layer [H R ~> kg m-2 or kg2 m-5].
247 real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim]
248 real :: dh ! The increment in layer thickness from
249 ! the present layer [H ~> m or kg m-2].
250 real :: ddz ! The increment in height change from the present layer [Z ~> m].
251 real :: bbl_thick ! The thickness of the bottom boundary layer [Z ~> m].
252 real :: bbl_thick_max ! A huge upper bound on the boundary layer thickness [Z ~> m].
253 real :: kv_bbl ! The bottom boundary layer viscosity [H Z T-1 ~> m2 s-1 or Pa s]
254 real :: c2f ! C2f = 2*f at velocity points [T-1 ~> s-1].
255 real :: u2_bg(szib_(g)) ! The square of an assumed background velocity, for calculating the mean
256 ! magnitude near the bottom for use in the quadratic bottom drag [L2 T-2 ~> m2 s-2].
257 real :: hwtot ! Sum of the thicknesses used to calculate
258 ! the near-bottom velocity magnitude [H ~> m or kg m-2].
259 real :: i_hwtot ! The Adcroft reciprocal of hwtot [H-1 ~> m-1 or m2 kg-1].
260 real :: dzwtot ! The vertical extent of the region used to calculate
261 ! the near-bottom velocity magnitude [Z ~> m].
262 real :: hutot ! Running sum of thicknesses times the velocity
263 ! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
264 real :: thtot ! Running sum of thickness times temperature [C H ~> degC m or degC kg m-2].
265 real :: shtot ! Running sum of thickness times salinity [S H ~> ppt m or ppt kg m-2].
266 real :: spv_htot ! Running sum of thickness times specific volume [H R-1 ~> m4 kg-1 or m]
267 real :: hweight ! The thickness of a layer that is within Hbbl
268 ! of the bottom [H ~> m or kg m-2].
269 real :: dzweight ! The counterpart of hweight in height units [Z ~> m].
270 real :: v_at_u, u_at_v ! v at a u point or vice versa [L T-1 ~> m s-1].
271 real :: rho0x400_g ! 400*Rho0/G_Earth, times unit conversion factors
272 ! [R T2 H-1 ~> kg s2 m-4 or s2 m-1].
273 ! The 400 is a constant proposed by Killworth and Edwards, 1999.
274 real, dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
275 rml ! The mixed layer coordinate density [R ~> kg m-3].
276 real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate
277 ! density [R L2 T-2 ~> Pa] (usually set to 2e7 Pa = 2000 dbar).
278
279 real :: d_vel ! The bottom depth relative to the shelfbreak depth at a velocity point [Z ~> m].
280 real :: dp, dm ! The bottom depths at the edges of a velocity cell relative to the
281 ! shelfbreak depth [Z ~> m].
282 real :: d_vel_p, d_vel_m ! The bottom depths in adjacent velocity points relative to the
283 ! shelfbreak depth [Z ~> m].
284 real :: crv ! crv is the curvature of the bottom depth across a
285 ! cell, times the cell width squared [Z ~> m].
286 real :: slope ! The absolute value of the bottom depth slope across
287 ! a cell times the cell width [Z ~> m].
288 real :: vol_bbl_chan ! The volume of the bottom boundary layer as used in the channel
289 ! drag parameterization, normalized by the full horizontal area
290 ! of the velocity cell [Z ~> m].
291 real :: vol_below(szk_(gv)+1) ! The volume below each interface, normalized by the full
292 ! horizontal area of a velocity cell [Z ~> m].
293 real :: l(szk_(gv)+1) ! The fraction of the full cell width that is open at
294 ! the depth of each interface [nondim].
295 ! The next 9 variables are only used for debugging.
296 real :: l_trig(szk_(gv)+1) ! The fraction of the full cell width that is open at
297 ! the depth of each interface from trigonometric expressions [nondim].
298 real :: vol_err_trig(szk_(gv)+1) ! The error in the volume below based on L_trig [Z ~> m]
299 real :: vol_err_iter(szk_(gv)+1) ! The error in the volume below based on L_iter [Z ~> m]
300 real :: norm_err_trig(szk_(gv)+1) ! vol_err_trig normalized by vol_below [nondim]
301 real :: norm_err_iter(szk_(gv)+1) ! vol_err_iter normalized by vol_below [nondim]
302 real :: dl_trig_itt(szk_(gv)+1) ! The difference between estimates of the fraction of the full cell
303 ! width that is open at the depth of each interface [nondim].
304 real :: max_dl_trig_itt ! The largest difference between L and L_trig, for debugging [nondim]
305 real :: max_norm_err_trig ! The largest magnitude value of norm_err_trig in a column [nondim]
306 real :: max_norm_err_iter ! The largest magnitude value of norm_err_iter in a column [nondim]
307
308 real :: h_neglect ! A thickness that is so small it is usually lost
309 ! in roundoff and can be neglected [H ~> m or kg m-2].
310 real :: dz_neglect ! A vertical distance that is so small it is usually lost
311 ! in roundoff and can be neglected [Z ~> m].
312 real :: usth ! ustar converted to units of H T-1 [H T-1 ~> m s-1 or kg m-2 s-1].
313 real :: root ! A temporary variable [H T-1 ~> m s-1 or kg m-2 s-1].
314
315 real :: cell_width ! The transverse width of the velocity cell [L ~> m].
316 real :: rayleigh ! A factor that is multiplied by the layer's velocity magnitude
317 ! to give the Rayleigh drag velocity, times a lateral distance to
318 ! thickness conversion factor [H L-1 ~> nondim or kg m-3].
319 real :: gam ! The ratio of the change in the open interface width
320 ! to the open interface width atop a cell [nondim].
321 real :: bbl_frac ! The fraction of a layer's drag that goes into the
322 ! viscous bottom boundary layer [nondim].
323 real :: bbl_visc_frac ! The fraction of all the drag that is expressed as
324 ! a viscous bottom boundary layer [nondim].
325 real :: h_bbl_fr ! The fraction of the bottom boundary layer in a layer [nondim].
326 real :: h_sum ! The sum of the thicknesses of the layers below the one being
327 ! worked on [H ~> m or kg m-2].
328 real :: tideampfac2_x_0p5 ! tideampfac2 multiplied by the c-grid averaging factor of 0.5
329 real, parameter :: c1_3 = 1.0/3.0, c1_6 = 1.0/6.0, c1_12 = 1.0/12.0 ! Rational constants [nondim]
330 real :: tmp ! A temporary variable, sometimes in [Z ~> m]
331 logical :: use_bbl_eos, do_i(szib_(g))
332 integer, dimension(2) :: eosdom ! The computational domain for the equation of state
333 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, m, n, k2, nkmb, nkml
334 integer :: is_obc, ie_obc, js_obc, je_obc
335 type(ocean_obc_type), pointer :: obc => null()
336
337 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
338 isq = g%isc-1 ; ieq = g%IecB ; jsq = g%jsc-1 ; jeq = g%JecB
339 nkmb = gv%nk_rho_varies ; nkml = gv%nkml
340 h_neglect = gv%H_subroundoff
341 dz_neglect = gv%dZ_subroundoff
342
343 rho0x400_g = 400.0*(gv%H_to_RZ / gv%g_Earth_Z_T2)
344 tideampfac2_x_0p5 = cs%tideampfac2*0.5
345
346 if (.not.cs%initialized) call mom_error(fatal,"MOM_set_viscosity(BBL): "//&
347 "Module must be initialized before it is used.")
348
349 if (.not.cs%bottomdraglaw) return
350
351 if (cs%debug) then
352 call uvchksum("Start set_viscous_BBL [uv]", u, v, g%HI, haloshift=1, unscale=us%L_T_to_m_s)
353 call hchksum(h,"Start set_viscous_BBL h", g%HI, haloshift=1, unscale=gv%H_to_m)
354 if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", g%HI, haloshift=1, unscale=us%C_to_degC)
355 if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", g%HI, haloshift=1, unscale=us%S_to_ppt)
356 if (allocated(tv%SpV_avg)) &
357 call hchksum(tv%SpV_avg, "Start set_viscous_BBL SpV_avg", g%HI, haloshift=1, unscale=us%kg_m3_to_R)
358 if (allocated(tv%SpV_avg)) call hchksum(tv%SpV_avg, "Cornerless SpV_avg", g%HI, &
359 haloshift=1, omit_corners=.true., unscale=us%kg_m3_to_R)
360 if (associated(tv%T)) call hchksum(tv%T, "Cornerless T", g%HI, haloshift=1, &
361 omit_corners=.true., unscale=us%C_to_degC)
362 if (associated(tv%S)) call hchksum(tv%S, "Cornerless S", g%HI, haloshift=1, &
363 omit_corners=.true., unscale=us%S_to_ppt)
364 endif
365
366 use_bbl_eos = associated(tv%eqn_of_state) .and. cs%BBL_use_EOS
367 obc => cs%OBC
368
369 if (.not.cs%bottomdragmap) then
370 cdrag_sqrt = sqrt(cs%cdrag)
371 cdrag_sqrt_h = cdrag_sqrt * us%L_to_m * gv%m_to_H
372 cdrag_sqrt_h_rl = cdrag_sqrt * us%L_to_Z * gv%RZ_to_H
373 cdrag_l_to_h = cs%cdrag * us%L_to_m * gv%m_to_H
374 cdrag_rl_to_h = cs%cdrag * us%L_to_Z * gv%RZ_to_H
375 endif
376 bbl_thick_max = g%Rad_Earth_L * us%L_to_Z
377 k2 = max(nkmb+1, 2)
378
379 ! Find the vertical distances across layers.
380 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
381
382! With a linear drag law, the friction velocity is already known.
383! if (CS%linear_drag) ustar(:) = cdrag_sqrt_H*CS%drag_bg_vel
384
385 if ((nkml>0) .and. .not.use_bbl_eos) then
386 eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)
387 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
388 !$OMP parallel do default(shared)
389 do k=1,nkmb ; do j=jsq,jeq+1
390 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, rml(:,j,k), tv%eqn_of_state, &
391 eosdom)
392 enddo ; enddo
393 endif
394
395 !$OMP parallel do default(shared)
396 do j=js-1,je ; do i=is-1,ie+1
397 d_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
398 mask_v(i,j) = g%mask2dCv(i,j)
399 enddo ; enddo
400 !$OMP parallel do default(shared)
401 do j=js-1,je+1 ; do i=is-1,ie
402 d_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
403 mask_u(i,j) = g%mask2dCu(i,j)
404 enddo ; enddo
405
406 if (associated(obc) .and. cs%Channel_drag) then
407 ! Use a one-sided projection of bottom depths at OBC points.
408 if (obc%v_N_OBCs_on_PE) then
409 js_obc = max(js-1, obc%Js_v_N_obc) ; je_obc = min(je, obc%Je_v_N_obc)
410 is_obc = max(is-1, obc%is_v_N_obc) ; ie_obc = min(ie+1, obc%ie_v_N_obc)
411 !$OMP parallel do default(shared)
412 do j=js_obc,je_obc ; do i=is_obc,ie_obc
413 if (obc%segnum_v(i,j) > 0) d_v(i,j) = g%bathyT(i,j) ! OBC_DIRECTION_N
414 enddo ; enddo
415 endif
416 if (obc%v_S_OBCs_on_PE) then
417 js_obc = max(js-1, obc%Js_v_S_obc) ; je_obc = min(je, obc%Je_v_S_obc)
418 is_obc = max(is-1, obc%is_v_S_obc) ; ie_obc = min(ie+1, obc%ie_v_S_obc)
419 !$OMP parallel do default(shared)
420 do j=js_obc,je_obc ; do i=is_obc,ie_obc
421 if (obc%segnum_v(i,j) < 0) d_v(i,j) = g%bathyT(i,j+1) ! OBC_DIRECTION_S
422 enddo ; enddo
423 endif
424 if (obc%u_E_OBCs_on_PE) then
425 js_obc = max(js-1, obc%js_u_E_obc) ; je_obc = min(je+1, obc%je_u_E_obc)
426 is_obc = max(is-1, obc%Is_u_E_obc) ; ie_obc = min(ie, obc%Ie_u_E_obc)
427 !$OMP parallel do default(shared)
428 do j=js_obc,je_obc ; do i=is_obc,ie_obc
429 if (obc%segnum_u(i,j) > 0) d_u(i,j) = g%bathyT(i,j) ! OBC_DIRECTION_E
430 enddo ; enddo
431 endif
432 if (obc%u_W_OBCs_on_PE) then
433 js_obc = max(js-1, obc%js_u_W_obc) ; je_obc = min(je+1, obc%je_u_W_obc)
434 is_obc = max(is-1, obc%Is_u_W_obc) ; ie_obc = min(ie, obc%Ie_u_W_obc)
435 !$OMP parallel do default(shared)
436 do j=js_obc,je_obc ; do i=is_obc,ie_obc
437 if (obc%segnum_u(i,j) < 0) d_u(i,j) = g%bathyT(i+1,j) ! OBC_DIRECTION_W
438 enddo ; enddo
439 endif
440 endif
441
442 if (associated(obc) .and. cs%Channel_drag) then ; do n=1,obc%number_of_segments
443 ! Now project bottom depths across cell-corner points in the OBCs. The two
444 ! projections have to occur in sequence and can not be combined easily.
445 if (.not. obc%segment(n)%on_pe) cycle
446 ! Use a one-sided projection of bottom depths at OBC points.
447 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
448 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
449 do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
450 if (obc%segment(n)%direction == obc_direction_n) then
451 d_u(i,j+1) = d_u(i,j) ; mask_u(i,j+1) = 0.0
452 elseif (obc%segment(n)%direction == obc_direction_s) then
453 d_u(i,j) = d_u(i,j+1) ; mask_u(i,j) = 0.0
454 endif
455 enddo
456 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
457 do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
458 if (obc%segment(n)%direction == obc_direction_e) then
459 d_v(i+1,j) = d_v(i,j) ; mask_v(i+1,j) = 0.0
460 elseif (obc%segment(n)%direction == obc_direction_w) then
461 d_v(i,j) = d_v(i+1,j) ; mask_v(i,j) = 0.0
462 endif
463 enddo
464 endif
465 enddo ; endif
466
467 if (.not.use_bbl_eos) rml_vel(:,:) = 0.0
468
469 ! Resetting Ray_[uv] is required by body force drag.
470 if (allocated(visc%Ray_u)) visc%Ray_u(:,:,:) = 0.0
471 if (allocated(visc%Ray_v)) visc%Ray_v(:,:,:) = 0.0
472
473 !$OMP parallel do default(private) shared(u,v,h,dz,tv,visc,G,GV,US,CS,Rml,nz,nkmb,nkml,K2, &
474 !$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G, &
475 !$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, &
476 !$OMP cdrag_L_to_H,cdrag_RL_to_H,use_BBL_EOS,BBL_thick_max, &
477 !$OMP OBC,D_u,D_v,mask_u,mask_v,pbv)
478 do j=jsq,jeq ; do m=1,2
479
480 if (m==1) then
481 ! m=1 refers to u-points
482 if (j<g%Jsc) cycle
483 is = isq ; ie = ieq
484 do i=is,ie
485 do_i(i) = (g%mask2dCu(i,j) > 0.0)
486 enddo
487 else
488 ! m=2 refers to v-points
489 is = g%isc ; ie = g%iec
490 do i=is,ie
491 do_i(i) = (g%mask2dCv(i,j) > 0.0)
492 enddo
493 endif
494
495 ! Calculate thickness at velocity points (u or v depending on value of m).
496 ! Also interpolate the ML density or T/S properties.
497 if (m==1) then ! u-points
498 do k=1,nz ; do i=is,ie
499 if (do_i(i)) then
500 if (u(i,j,k) * (h(i+1,j,k) - h(i,j,k)) >= 0) then
501 ! If the flow is from thin to thick then bias towards the thinner thickness
502 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
503 (h(i,j,k) + h(i+1,j,k) + h_neglect)
504 dz_at_vel(i,k) = 2.0*dz(i,j,k)*dz(i+1,j,k) / &
505 (dz(i,j,k) + dz(i+1,j,k) + dz_neglect)
506 else
507 ! If the flow is from thick to thin then use the simple average thickness
508 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
509 dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
510 endif
511 endif
512 h_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
513 dz_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
514 enddo ; enddo
515 if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
516 ! Perhaps these should be thickness weighted.
517 t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
518 s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
519 enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
520 rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
521 enddo ; enddo ; endif
522 if (allocated(tv%SpV_avg)) then ; do k=1,nz ; do i=is,ie
523 spv_vel(i,k) = 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i+1,j,k))
524 enddo ; enddo ; endif
525 else ! v-points
526 do k=1,nz ; do i=is,ie
527 if (do_i(i)) then
528 if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
529 ! If the flow is from thin to thick then bias towards the thinner thickness
530 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
531 (h(i,j,k) + h(i,j+1,k) + h_neglect)
532 dz_at_vel(i,k) = 2.0*dz(i,j,k)*dz(i,j+1,k) / &
533 (dz(i,j,k) + dz(i,j+1,k) + dz_neglect)
534 else
535 ! If the flow is from thick to thin then use the simple average thickness
536 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
537 dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
538 endif
539 endif
540 h_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
541 dz_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
542 enddo ; enddo
543 if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
544 ! Perhaps these should be thickness weighted.
545 t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
546 s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
547 enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
548 rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
549 enddo ; enddo ; endif
550 if (allocated(tv%SpV_avg)) then ; do k=1,nz ; do i=is,ie
551 spv_vel(i,k) = 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j+1,k))
552 enddo ; enddo ; endif
553 endif
554
555 if (associated(obc)) then ; if (obc%number_of_segments > 0) then
556 ! Apply a zero gradient projection of thickness across OBC points.
557 if (m==1) then
558 do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= 0)) then
559 if (obc%segnum_u(i,j) > 0) then ! OBC_DIRECTION_E
560 do k=1,nz
561 h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
562 dz_at_vel(i,k) = dz(i,j,k) ; dz_vel(i,k) = dz(i,j,k)
563 enddo
564 if (use_bbl_eos) then
565 do k=1,nz
566 t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
567 enddo
568 else
569 do k=1,nkmb
570 rml_vel(i,k) = rml(i,j,k)
571 enddo
572 endif
573 if (allocated(tv%SpV_avg)) then ; do k=1,nz
574 spv_vel(i,k) = tv%SpV_avg(i,j,k)
575 enddo ; endif
576 elseif (obc%segnum_u(i,j) < 0) then ! OBC_DIRECTION_W
577 do k=1,nz
578 h_at_vel(i,k) = h(i+1,j,k) ; h_vel(i,k) = h(i+1,j,k)
579 dz_at_vel(i,k) = dz(i+1,j,k) ; dz_vel(i,k) = dz(i+1,j,k)
580 enddo
581 if (use_bbl_eos) then
582 do k=1,nz
583 t_vel(i,k) = tv%T(i+1,j,k) ; s_vel(i,k) = tv%S(i+1,j,k)
584 enddo
585 else
586 do k=1,nkmb
587 rml_vel(i,k) = rml(i+1,j,k)
588 enddo
589 endif
590 if (allocated(tv%SpV_avg)) then ; do k=1,nz
591 spv_vel(i,k) = tv%SpV_avg(i+1,j,k)
592 enddo ; endif
593 endif
594 endif ; enddo
595 else
596 do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= 0)) then
597 if (obc%segnum_v(i,j) > 0) then ! OBC_DIRECTION_N
598 do k=1,nz
599 h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
600 dz_at_vel(i,k) = dz(i,j,k) ; dz_vel(i,k) = dz(i,j,k)
601 enddo
602 if (use_bbl_eos) then
603 do k=1,nz
604 t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
605 enddo
606 else
607 do k=1,nkmb
608 rml_vel(i,k) = rml(i,j,k)
609 enddo
610 endif
611 if (allocated(tv%SpV_avg)) then ; do k=1,nz
612 spv_vel(i,k) = tv%SpV_avg(i,j,k)
613 enddo ; endif
614 elseif (obc%segnum_v(i,j) < 0) then ! OBC_DIRECTION_S
615 do k=1,nz
616 h_at_vel(i,k) = h(i,j+1,k) ; h_vel(i,k) = h(i,j+1,k)
617 dz_at_vel(i,k) = dz(i,j+1,k) ; dz_vel(i,k) = dz(i,j+1,k)
618 enddo
619 if (use_bbl_eos) then
620 do k=1,nz
621 t_vel(i,k) = tv%T(i,j+1,k) ; s_vel(i,k) = tv%S(i,j+1,k)
622 enddo
623 else
624 do k=1,nkmb
625 rml_vel(i,k) = rml(i,j+1,k)
626 enddo
627 endif
628 if (allocated(tv%SpV_avg)) then ; do k=1,nz
629 spv_vel(i,k) = tv%SpV_avg(i,j+1,k)
630 enddo ; endif
631 endif
632 endif ; enddo
633 endif
634 endif ; endif
635
636 ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
637 if (cs%BBL_use_tidal_bg) then
638 do i=is,ie ; if (do_i(i)) then ; if (m==1) then
639 u2_bg(i) = tideampfac2_x_0p5 * ( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
640 g%mask2dT(i+1,j)*(cs%tideamp(i+1,j)*cs%tideamp(i+1,j)) )
641 else
642 u2_bg(i) = tideampfac2_x_0p5 * ( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
643 g%mask2dT(i,j+1)*(cs%tideamp(i,j+1)*cs%tideamp(i,j+1)) )
644 endif ; endif ; enddo
645 else
646 do i=is,ie ; if (do_i(i)) then
647 u2_bg(i) = cs%drag_bg_vel * cs%drag_bg_vel
648 endif ; enddo
649 endif
650
651 if (use_bbl_eos .or. cs%body_force_drag .or. .not.cs%linear_drag) then
652 ! Calculate the mean velocity magnitude over the bottommost CS%Hbbl of
653 ! the water column for determining the quadratic bottom drag.
654 ! Used in ustar(i)
655 do i=is,ie ; if (do_i(i)) then
656 htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
657 dztot_vel = 0.0 ; dzwtot = 0.0
658 thtot = 0.0 ; shtot = 0.0 ; spv_htot = 0.0
659
660 if (cs%bottomdragmap) then
661 if (m==1) then
662 cdrag_sqrt = sqrt(cs%cdrag_u(i,j))
663 else
664 cdrag_sqrt = sqrt(cs%cdrag_v(i,j))
665 endif
666 cdrag_sqrt_h = cdrag_sqrt * us%L_to_m * gv%m_to_H
667 cdrag_sqrt_h_rl = cdrag_sqrt * us%L_to_Z * gv%RZ_to_H
668 endif
669
670 do k=nz,1,-1
671
672 if (htot_vel>=cs%Hbbl) exit ! terminate the k loop
673
674 hweight = min(cs%Hbbl - htot_vel, h_at_vel(i,k))
675 if (hweight < 1.5*gv%Angstrom_H + h_neglect) cycle
676 dzweight = min(cs%dz_bbl - dztot_vel, dz_at_vel(i,k))
677
678 htot_vel = htot_vel + h_at_vel(i,k)
679 hwtot = hwtot + hweight
680 dztot_vel = dztot_vel + dz_at_vel(i,k)
681 dzwtot = dzwtot + dzweight
682
683 if ((.not.cs%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
684 v_at_u = set_v_at_u(v, h, g, gv, i, j, k, mask_v, obc)
685 hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + v_at_u*v_at_u + u2_bg(i))
686 else
687 u_at_v = set_u_at_v(u, h, g, gv, i, j, k, mask_u, obc)
688 hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + u_at_v*u_at_v + u2_bg(i))
689 endif ; endif
690
691 if (use_bbl_eos .and. (hweight >= 0.0)) then
692 thtot = thtot + hweight * t_vel(i,k)
693 shtot = shtot + hweight * s_vel(i,k)
694 endif
695 if (allocated(tv%SpV_avg) .and. (hweight >= 0.0)) then
696 spv_htot = spv_htot + hweight * spv_vel(i,k)
697 endif
698 enddo ! end of k loop
699
700 ! Find the Adcroft reciprocal of the total thickness weights
701 i_hwtot = 0.0 ; if (hwtot > 0.0) i_hwtot = 1.0 / hwtot
702
703 ! Set u* based on u*^2 = Cdrag u_bbl^2
704 if ((hwtot <= 0.0) .or. (cs%linear_drag .and. .not.allocated(tv%SpV_avg))) then
705 ustar(i) = cdrag_sqrt_h * cs%drag_bg_vel
706 elseif (cs%linear_drag .and. allocated(tv%SpV_avg)) then
707 ustar(i) = cdrag_sqrt_h_rl * cs%drag_bg_vel * (hwtot / spv_htot)
708 elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
709 ustar(i) = cdrag_sqrt_h_rl * hutot / spv_htot
710 else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
711 ustar(i) = cdrag_sqrt_h * hutot / hwtot
712 endif
713
714 umag_avg(i) = hutot * i_hwtot
715 h_bbl_drag(i) = hwtot
716 dz_bbl_drag(i) = dzwtot
717
718 if (use_bbl_eos) then ; if (hwtot > 0.0) then
719 t_eos(i) = thtot/hwtot ; s_eos(i) = shtot/hwtot
720 else
721 t_eos(i) = 0.0 ; s_eos(i) = 0.0
722 endif ; endif
723
724 ! Diagnostic BBL flow speed at u- and v-points.
725 if (cs%id_bbl_u>0 .and. m==1) then
726 if (hwtot > 0.0) cs%bbl_u(i,j) = hutot/hwtot
727 elseif (cs%id_bbl_v>0 .and. m==2) then
728 if (hwtot > 0.0) cs%bbl_v(i,j) = hutot/hwtot
729 endif
730
731 endif ; enddo
732 else
733 do i=is,ie
734 if (cs%bottomdragmap) then
735 if (m==1) then
736 cdrag_sqrt = sqrt(cs%cdrag_u(i,j))
737 else
738 cdrag_sqrt = sqrt(cs%cdrag_v(i,j))
739 endif
740 cdrag_sqrt_h = cdrag_sqrt * us%L_to_m * gv%m_to_H
741 endif
742 ustar(i) = cdrag_sqrt_h * cs%drag_bg_vel
743 enddo
744 endif ! Not linear_drag
745
746 if (use_bbl_eos) then
747 if (associated(tv%p_surf)) then
748 if (m==1) then ; do i=is,ie ; press(i) = 0.5*(tv%p_surf(i,j) + tv%p_surf(i+1,j)) ; enddo
749 else ; do i=is,ie ; press(i) = 0.5*(tv%p_surf(i,j) + tv%p_surf(i,j+1)) ; enddo ; endif
750 else
751 do i=is,ie ; press(i) = 0.0 ; enddo
752 endif
753 do i=is,ie ; if (.not.do_i(i)) then ; t_eos(i) = 0.0 ; s_eos(i) = 0.0 ; endif ; enddo
754 do k=1,nz ; do i=is,ie
755 press(i) = press(i) + (gv%H_to_RZ*gv%g_Earth) * h_vel(i,k)
756 enddo ; enddo
757 call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, tv%eqn_of_state, &
758 (/is-g%IsdB+1,ie-g%IsdB+1/) )
759 endif
760
761 ! Find a BBL thickness given by equation 2.20 of Killworth and Edwards, 1999:
762 ! ( f h / Cn u* )^2 + ( N h / Ci u* ) = 1
763 ! where Cn=0.5 and Ci=20 (constants suggested by Zilitinkevich and Mironov, 1996).
764 ! Eq. 2.20 can be expressed in terms of boundary layer thicknesses limited by
765 ! rotation (h_f) and stratification (h_N):
766 ! ( h / h_f )^2 + ( h / h_N ) = 1
767 ! When stratification dominates h_N<<h_f, and vice versa.
768 do i=is,ie ; if (do_i(i)) then
769 ! The 400.0 in this expression is the square of a Ci introduced in KW99, eq. 2.22.
770 ustarsq = rho0x400_g * ustar(i)**2 ! Note not in units of u*^2 but [H R ~> kg m-2 or kg2 m-5]
771 htot = 0.0
772 dztot = 0.0
773
774 if (cs%bottomdragmap) then
775 if (m==1) then
776 cdrag = cs%cdrag_u(i,j)
777 else
778 cdrag = cs%cdrag_v(i,j)
779 endif
780 cdrag_l_to_h = cdrag * us%L_to_m * gv%m_to_H
781 cdrag_rl_to_h = cdrag * us%L_to_Z * gv%RZ_to_H
782 endif
783
784 ! Calculate the thickness of a stratification limited BBL ignoring rotation:
785 ! h_N = Ci u* / N (limit of KW99 eq. 2.20 for |f|->0)
786 ! For layer mode, N^2 = g'/h. Since (Ci u*)^2 = (h_N N)^2 = h_N g' then
787 ! h_N = (Ci u*)^2 / g' (KW99, eq, 2.22)
788 ! Starting from the bottom, integrate the stratification upward until h_N N balances Ci u*
789 ! or in layer mode
790 ! h_N Delta rho ~ (Ci u*)^2 rho0 / g
791 ! where the rhs is stored in variable ustarsq.
792 ! The method was described in Stephens and Hallberg 2000 (unpublished and lost manuscript).
793 if (use_bbl_eos) then
794 thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
795 do k=nz,2,-1
796 if (h_at_vel(i,k) <= 0.0) cycle
797
798 ! Delta rho * h_bbl assuming everything below is homogenized
799 oldfn = dr_dt(i)*(thtot - t_vel(i,k)*htot) + &
800 dr_ds(i)*(shtot - s_vel(i,k)*htot)
801 if (oldfn >= ustarsq) exit
802
803 ! Local Delta rho * h_bbl at interface
804 dfn = (dr_dt(i)*(t_vel(i,k) - t_vel(i,k-1)) + &
805 dr_ds(i)*(s_vel(i,k) - s_vel(i,k-1))) * &
806 (h_at_vel(i,k) + htot)
807
808 if ((oldfn + dfn) <= ustarsq) then
809 ! Use whole layer
810 dh = h_at_vel(i,k)
811 ddz = dz_at_vel(i,k)
812 else
813 ! Use only part of the layer
814 frac_used = sqrt((ustarsq-oldfn) / (dfn))
815 dh = h_at_vel(i,k) * frac_used
816 ddz = dz_at_vel(i,k) * frac_used
817 endif
818
819 ! Increment total BBL thickness and cumulative T and S
820 htot = htot + dh
821 dztot = dztot + ddz
822 thtot = thtot + t_vel(i,k)*dh ; shtot = shtot + s_vel(i,k)*dh
823 enddo
824 if ((oldfn < ustarsq) .and. h_at_vel(i,1) > 0.0) then
825 ! Layer 1 might be part of the BBL.
826 if (dr_dt(i) * (thtot - t_vel(i,1)*htot) + &
827 dr_ds(i) * (shtot - s_vel(i,1)*htot) < ustarsq) then
828 htot = htot + h_at_vel(i,1)
829 dztot = dztot + dz_at_vel(i,1)
830 endif
831 endif ! Examination of layer 1.
832 else ! Use Rlay and/or the coordinate density as density variables.
833 rhtot = 0.0
834 do k=nz,k2,-1
835 oldfn = rhtot - gv%Rlay(k)*htot
836 dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h_at_vel(i,k)+htot)
837
838 if (oldfn >= ustarsq) then
839 cycle
840 elseif ((oldfn + dfn) <= ustarsq) then
841 dh = h_at_vel(i,k)
842 ddz = dz_at_vel(i,k)
843 else
844 frac_used = sqrt((ustarsq-oldfn) / (dfn))
845 dh = h_at_vel(i,k) * frac_used
846 ddz = dz_at_vel(i,k) * frac_used
847 endif
848
849 htot = htot + dh
850 dztot = dztot + ddz
851 rhtot = rhtot + gv%Rlay(k)*dh
852 enddo
853 if (nkml>0) then
854 do k=nkmb,2,-1
855 oldfn = rhtot - rml_vel(i,k)*htot
856 dfn = (rml_vel(i,k) - rml_vel(i,k-1)) * (h_at_vel(i,k)+htot)
857
858 if (oldfn >= ustarsq) then
859 cycle
860 elseif ((oldfn + dfn) <= ustarsq) then
861 dh = h_at_vel(i,k)
862 ddz = dz_at_vel(i,k)
863 else
864 frac_used = sqrt((ustarsq-oldfn) / (dfn))
865 dh = h_at_vel(i,k) * frac_used
866 ddz = dz_at_vel(i,k) * frac_used
867 endif
868
869 htot = htot + dh
870 dztot = dztot + ddz
871 rhtot = rhtot + rml_vel(i,k)*dh
872 enddo
873 if (rhtot - rml_vel(i,1)*htot < ustarsq) then
874 htot = htot + h_at_vel(i,1)
875 dztot = dztot + dz_at_vel(i,1)
876 endif
877 else
878 if (rhtot - gv%Rlay(1)*htot < ustarsq) then
879 htot = htot + h_at_vel(i,1)
880 dztot = dztot + dz_at_vel(i,1)
881 endif
882 endif
883 endif ! use_BBL_EOS
884
885 ! Value of 2*f at u- or v-points.
886 if (m==1) then ; c2f = g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j)
887 else ; c2f = g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j) ; endif
888
889 ! The thickness of a rotation limited BBL ignoring stratification is
890 ! h_f ~ Cn u* / f (limit of KW99 eq. 2.20 for N->0).
891 ! The buoyancy limit of BBL thickness (h_N) is already in the variable htot from above.
892 ! Substituting x = h_N/h into KW99 eq. 2.20 yields the quadratic
893 ! x^2 - x = (h_N / h_f)^2
894 ! for which the positive root is
895 ! xp = 1/2 + sqrt( 1/4 + (h_N/h_f)^2 )
896 ! and thus h_bbl = h_N / xp . Since h_f = Cn u*/f and Cn=0.5
897 ! xp = 1/2 + sqrt( 1/4 + (2 f h_N/u*)^2 )
898 ! To avoid dividing by zero if u*=0 then
899 ! xp u* = 1/2 u* + sqrt( 1/4 u*^2 + (2 f h_N)^2 )
900 if (cs%cdrag * u2_bg(i) <= 0.0) then
901 ! This avoids NaNs and overflows, and could be used in all cases,
902 ! but is not bitwise identical to the current code.
903 usth = ustar(i) ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
904 if (dztot*usth <= (cs%BBL_thick_min+dz_neglect) * (0.5*usth + root)) then
905 bbl_thick = cs%BBL_thick_min
906 else
907 ! The following expression reads
908 ! h_bbl = h_N u* / ( 1/2 u* + sqrt( 1/4 u*^2 + ( 2 f h_N )^2 ) )
909 ! which is h_bbl = h_N u*/(xp u*) as described above.
910 bbl_thick = (dztot * usth) / (0.5*usth + root)
911 endif
912 else
913 ! The following expression reads
914 ! h_bbl = h_N / ( 1/2 + sqrt( 1/4 + ( 2 f h_N / u* )^2 ) )
915 ! which is h_bbl = h_N/xp as described above.
916 bbl_thick = dztot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f / (ustar(i)*ustar(i)) ) )
917
918 if (bbl_thick < cs%BBL_thick_min) bbl_thick = cs%BBL_thick_min
919 endif
920
921 ! Store the normalized bottom boundary layer volume.
922 if (cs%Channel_drag) vol_bbl_chan = bbl_thick
923
924 ! If there is Richardson number dependent mixing, that determines
925 ! the vertical extent of the bottom boundary layer, and there is no
926 ! need to set that scale here. In fact, viscously reducing the
927 ! shears over an excessively large region reduces the efficacy of
928 ! the Richardson number dependent mixing.
929 ! In other words, if using RiNo_mix then CS%dz_bbl acts as an upper bound on
930 ! bbl_thick.
931 if ((bbl_thick > 0.5*cs%dz_bbl) .and. (cs%RiNo_mix)) bbl_thick = 0.5*cs%dz_bbl
932
933 ! If drag is a body force, bbl_thick is HBBL
934 if (cs%body_force_drag) bbl_thick = dz_bbl_drag(i)
935
936 if (cs%Channel_drag) then
937
938 vol_below(nz+1) = 0.0
939 do k=nz,1,-1
940 vol_below(k) = vol_below(k+1) + dz_vel(i,k)
941 enddo
942
943 ! Find the bathymetry at adjacent points relative to the shelf break. For now this
944 ! shelf break depth is set with a global constant, but it could vary in space.
945 if (m==1) then
946 d_vel = d_u(i,j) - cs%channel_break_depth
947 d_vel_p = g%mask2dCu(i,j+1) * (d_u(i,j+1) - cs%channel_break_depth)
948 d_vel_m = g%mask2dCu(i,j-1) * (d_u(i,j-1) - cs%channel_break_depth)
949 else
950 d_vel = d_v(i,j) - cs%channel_break_depth
951 d_vel_p = g%mask2dCv(i+1,j) * (d_v(i+1,j) - cs%channel_break_depth)
952 d_vel_m = g%mask2dCv(i-1,j) * (d_v(i-1,j) - cs%channel_break_depth)
953 endif
954 ! This profile uses a harmonic mean bottom depth below some reference value to
955 ! roughly mimic the topographic shape at and beneath a continental shelf break.
956 ! Above this a simple arithmetic mean is used.
957 if ((d_vel > 0.0) .and. (d_vel_p > 0.0)) then
958 dp = 2.0 * d_vel * d_vel_p / (d_vel + d_vel_p)
959 else ! This is above the shelf-break, noting that D is positive downward.
960 dp = 0.5 * (min(d_vel, 0.0) + min(d_vel_p, 0.0))
961 endif
962 if ((d_vel > 0.0) .and. (d_vel_m > 0.0)) then
963 dm = 2.0 * d_vel * d_vel_m / (d_vel + d_vel_m)
964 else ! This is above the shelf-break, noting that D is positive downward.
965 dm = 0.5 * (min(d_vel, 0.0) + min(d_vel_m, 0.0))
966 endif
967 if (dm > dp) then ; tmp = dp ; dp = dm ; dm = tmp ; endif
968 crv = 3.0*(dp + dm - 2.0*d_vel)
969 slope = dp - dm
970
971 ! If the curvature is small enough, there is no reason not to assume
972 ! a uniformly sloping or flat bottom.
973 if (abs(crv) < 1e-2*(slope + cs%BBL_thick_min)) crv = 0.0
974
975 ! Determine the normalized open length (L) at each interface.
976 if (crv == 0.0) then
977 call find_l_open_uniform_slope(vol_below, dp, dm, l, gv)
978 elseif (crv > 0.0) then
979 if (cs%concave_trigonometric_L) then
980 call find_l_open_concave_trigonometric(vol_below, d_vel, dp, dm, l, gv)
981 else
982 call find_l_open_concave_iterative(vol_below, d_vel, dp, dm, l, gv)
983 if (cs%debug) then
984 ! The tests in this block reveal that the iterative and trigonometric solutions are
985 ! mathematically equivalent, but in some cases the iterative solution is consistent
986 ! at roundoff, but that the trigonmetric solutions have errors that can be several
987 ! orders of magnitude larger in some cases.
988 call find_l_open_concave_trigonometric(vol_below, d_vel, dp, dm, l_trig, gv)
989 call test_l_open_concave(vol_below, d_vel, dp, dm, l_trig, vol_err_trig, gv)
990 call test_l_open_concave(vol_below, d_vel, dp, dm, l, vol_err_iter, gv)
991 max_dl_trig_itt = 0.0 ; max_norm_err_trig = 0.0 ; max_norm_err_iter = 0.0
992 norm_err_trig(:) = 0.0 ; norm_err_iter(:) = 0.0
993 do k=1,nz+1
994 dl_trig_itt(k) = l_trig(k) - l(k)
995 if (abs(dl_trig_itt(k)) > abs(max_dl_trig_itt)) max_dl_trig_itt = dl_trig_itt(k)
996 norm_err_trig(k) = vol_err_trig(k) / (vol_below(k) + dz_neglect)
997 norm_err_iter(k) = vol_err_iter(k) / (vol_below(k) + dz_neglect)
998 if (abs(norm_err_trig(k)) > abs(max_norm_err_trig)) max_norm_err_trig = norm_err_trig(k)
999 if (abs(norm_err_iter(k)) > abs(max_norm_err_iter)) max_norm_err_iter = norm_err_iter(k)
1000 enddo
1001 if (abs(max_dl_trig_itt) > 1.0e-13) &
1002 k = nz+1 ! This is here only to use as a break point for a debugger.
1003 if (abs(max_norm_err_trig) > 1.0e-13) &
1004 k = nz+1 ! This is here only to use as a break point for a debugger.
1005 if (abs(max_norm_err_iter) > 1.0e-13) &
1006 k = nz+1 ! This is here only to use as a break point for a debugger.
1007 endif
1008 endif
1009 else ! crv < 0.0
1010 call find_l_open_convex(vol_below, d_vel, dp, dm, l, gv, us, cs)
1011 endif ! end of crv<0 cases.
1012
1013 ! Determine the Rayleigh drag contributions.
1014
1015 ! The drag within the bottommost Vol_bbl_chan is applied as a part of an enhanced bottom
1016 ! viscosity, while above this the drag is applied directly to the layers in question as a
1017 ! Rayleigh drag term.
1018
1019 ! Restrict the volume over which the channel drag is applied from the previously determined value.
1020 if (cs%Chan_drag_max_vol >= 0.0) vol_bbl_chan = min(vol_bbl_chan, cs%Chan_drag_max_vol)
1021
1022 bbl_visc_frac = 0.0
1023 do k=nz,1,-1
1024 !modify L(K) for porous barrier parameterization
1025 if (m==1) then ; l(k) = l(k)*pbv%por_layer_widthU(i,j,k)
1026 else ; l(k) = l(k)*pbv%por_layer_widthV(i,j,k) ; endif
1027
1028 ! Determine the drag contributing to the bottom boundary layer
1029 ! and the Rayleigh drag that acts on each layer.
1030 if (l(k) > l(k+1)) then
1031 if (vol_below(k+1) < vol_bbl_chan) then
1032 bbl_frac = (1.0-vol_below(k+1)/vol_bbl_chan)**2
1033 bbl_visc_frac = bbl_visc_frac + bbl_frac*(l(k) - l(k+1))
1034 else
1035 bbl_frac = 0.0
1036 endif
1037
1038 if (allocated(tv%SpV_avg)) then
1039 cdrag_conv = cdrag_rl_to_h / spv_vel(i,k)
1040 else
1041 cdrag_conv = cdrag_l_to_h
1042 endif
1043
1044 h_vel_pos = h_vel(i,k) + h_neglect
1045 if (m==1) then ; cell_width = g%dy_Cu(i,j)*pbv%por_face_areaU(i,j,k)
1046 else ; cell_width = g%dx_Cv(i,j)*pbv%por_face_areaV(i,j,k) ; endif
1047 gam = 1.0 - l(k+1)/l(k)
1048 rayleigh = cdrag_conv * (l(k)-l(k+1)) * (1.0-bbl_frac) * &
1049 (12.0*cs%c_Smag*h_vel_pos) / (12.0*cs%c_Smag*h_vel_pos + &
1050 cdrag_conv * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell_width)
1051 else ! This layer feels no drag.
1052 rayleigh = 0.0
1053 endif
1054
1055 if (m==1) then
1056 if (rayleigh > 0.0) then
1057 v_at_u = set_v_at_u(v, h, g, gv, i, j, k, mask_v, obc)
1058 visc%Ray_u(i,j,k) = rayleigh * sqrt(u(i,j,k)*u(i,j,k) + v_at_u*v_at_u + u2_bg(i))
1059 else ; visc%Ray_u(i,j,k) = 0.0 ; endif
1060 else
1061 if (rayleigh > 0.0) then
1062 u_at_v = set_u_at_v(u, h, g, gv, i, j, k, mask_u, obc)
1063 visc%Ray_v(i,j,k) = rayleigh * sqrt(v(i,j,k)*v(i,j,k) + u_at_v*u_at_v + u2_bg(i))
1064 else ; visc%Ray_v(i,j,k) = 0.0 ; endif
1065 endif
1066
1067 enddo ! k loop to determine visc%Ray_[uv].
1068
1069 ! Set the near-bottom viscosity to a value which will give
1070 ! the correct stress when the shear occurs over bbl_thick.
1071 ! See next block for explanation.
1072 if (cs%correct_BBL_bounds .and. &
1073 cdrag_sqrt*ustar(i)*bbl_thick*bbl_visc_frac <= cs%Kv_BBL_min) then
1074 ! If the bottom stress implies less viscosity than Kv_BBL_min then
1075 ! set kv_bbl to the bound and recompute bbl_thick to be consistent
1076 ! but with a ridiculously large upper bound on thickness (for Cd u*=0)
1077 kv_bbl = cs%Kv_BBL_min
1078 if ((cdrag_sqrt*ustar(i))*bbl_visc_frac*bbl_thick_max > kv_bbl) then
1079 bbl_thick = kv_bbl / ( (cdrag_sqrt*ustar(i)) * bbl_visc_frac )
1080 else
1081 bbl_thick = bbl_thick_max
1082 endif
1083 else
1084 kv_bbl = (cdrag_sqrt*ustar(i)) * bbl_thick*bbl_visc_frac
1085 endif
1086
1087 else ! Not Channel_drag.
1088 ! Set the near-bottom viscosity to a value which will give
1089 ! the correct stress when the shear occurs over bbl_thick.
1090 ! - The bottom stress is tau_b = Cdrag * u_bbl^2
1091 ! - u_bbl was calculated by averaging flow over CS%Hbbl
1092 ! (and includes unresolved tidal components)
1093 ! - u_bbl is embedded in u* since u*^2 = Cdrag u_bbl^2
1094 ! - The average shear in the BBL is du/dz = 2 * u_bbl / h_bbl
1095 ! (which assumes a linear profile, hence the "2")
1096 ! - bbl_thick was bounded to <= 0.5 * CS%dz_bbl
1097 ! - The viscous stress kv_bbl du/dz should balance tau_b
1098 ! Cdrag u_bbl^2 = kv_bbl du/dz
1099 ! = 2 kv_bbl u_bbl
1100 ! so
1101 ! kv_bbl = 0.5 h_bbl Cdrag u_bbl
1102 ! = 0.5 h_bbl sqrt(Cdrag) u*
1103 if (cs%correct_BBL_bounds .and. &
1104 cdrag_sqrt*ustar(i)*bbl_thick <= cs%Kv_BBL_min) then
1105 ! If the bottom stress implies less viscosity than Kv_BBL_min then
1106 ! set kv_bbl to the bound and recompute bbl_thick to be consistent
1107 ! but with a ridiculously large upper bound on thickness (for Cd u*=0)
1108 kv_bbl = cs%Kv_BBL_min
1109 if ((cdrag_sqrt*ustar(i))*bbl_thick_max > kv_bbl) then
1110 bbl_thick = kv_bbl / ( cdrag_sqrt*ustar(i) )
1111 else
1112 bbl_thick = bbl_thick_max
1113 endif
1114 else
1115 kv_bbl = (cdrag_sqrt*ustar(i)) * bbl_thick
1116 endif
1117 endif
1118
1119 if (cs%body_force_drag) then ; if (h_bbl_drag(i) > 0.0) then
1120 ! Increment the Rayleigh drag as a way introduce the bottom drag as a body force.
1121 h_sum = 0.0
1122 i_hwtot = 1.0 / h_bbl_drag(i)
1123 do k=nz,1,-1
1124 h_bbl_fr = min(h_bbl_drag(i) - h_sum, h_at_vel(i,k)) * i_hwtot
1125 if (allocated(tv%SpV_avg)) then
1126 cdrag_conv = cdrag_rl_to_h / spv_vel(i,k)
1127 else
1128 cdrag_conv = cdrag_l_to_h
1129 endif
1130 if (m==1) then
1131 visc%Ray_u(i,j,k) = visc%Ray_u(i,j,k) + (cdrag_conv * umag_avg(i)) * h_bbl_fr
1132 else
1133 visc%Ray_v(i,j,k) = visc%Ray_v(i,j,k) + (cdrag_conv * umag_avg(i)) * h_bbl_fr
1134 endif
1135 h_sum = h_sum + h_at_vel(i,k)
1136 if (h_sum >= h_bbl_drag(i)) exit ! The top of this layer is above the drag zone.
1137 enddo
1138 ! Do not enhance the near-bottom viscosity in this case.
1139 kv_bbl = cs%Kv_BBL_min
1140 endif ; endif
1141
1142 kv_bbl = max(cs%Kv_BBL_min, kv_bbl)
1143 if (m==1) then
1144 visc%bbl_thick_u(i,j) = bbl_thick
1145 if (allocated(visc%Kv_bbl_u)) visc%Kv_bbl_u(i,j) = kv_bbl
1146 else
1147 visc%bbl_thick_v(i,j) = bbl_thick
1148 if (allocated(visc%Kv_bbl_v)) visc%Kv_bbl_v(i,j) = kv_bbl
1149 endif
1150 endif ; enddo ! end of i loop
1151 enddo ; enddo ! end of m & j loops
1152
1153! Offer diagnostics for averaging
1154 if (cs%id_bbl_thick_u > 0) &
1155 call post_data(cs%id_bbl_thick_u, visc%bbl_thick_u, cs%diag)
1156 if (cs%id_kv_bbl_u > 0) &
1157 call post_data(cs%id_kv_bbl_u, visc%kv_bbl_u, cs%diag)
1158 if (cs%id_bbl_u > 0) &
1159 call post_data(cs%id_bbl_u, cs%bbl_u, cs%diag)
1160 if (cs%id_bbl_thick_v > 0) &
1161 call post_data(cs%id_bbl_thick_v, visc%bbl_thick_v, cs%diag)
1162 if (cs%id_kv_bbl_v > 0) &
1163 call post_data(cs%id_kv_bbl_v, visc%kv_bbl_v, cs%diag)
1164 if (cs%id_bbl_v > 0) &
1165 call post_data(cs%id_bbl_v, cs%bbl_v, cs%diag)
1166 if (cs%id_Ray_u > 0) &
1167 call post_data(cs%id_Ray_u, visc%Ray_u, cs%diag)
1168 if (cs%id_Ray_v > 0) &
1169 call post_data(cs%id_Ray_v, visc%Ray_v, cs%diag)
1170
1171 if (cs%debug) then
1172 if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) &
1173 call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, g%HI, haloshift=0, &
1174 unscale=gv%H_to_m*us%s_to_T, scalar_pair=.true.)
1175 if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) &
1176 call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, g%HI, &
1177 haloshift=0, unscale=gv%HZ_T_to_m2_s, scalar_pair=.true.)
1178 if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) &
1179 call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, &
1180 g%HI, haloshift=0, unscale=us%Z_to_m, scalar_pair=.true.)
1181 endif
1182
1183end subroutine set_viscous_bbl
1184
1185!> Determine the normalized open length of each interface, given the edge depths and normalized
1186!! volumes below each interface.
1187subroutine find_l_open_uniform_slope(vol_below, Dp, Dm, L, GV)
1188 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1189 real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1190 !! the full horizontal area of a velocity cell [Z ~> m]
1191 real, intent(in) :: Dp !< The larger of the two depths at the edge
1192 !! of a velocity cell [Z ~> m]
1193 real, intent(in) :: Dm !< The smaller of the two depths at the edge
1194 !! of a velocity cell [Z ~> m]
1195 real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1196 !! the depth of each interface [nondim]
1197
1198 ! Local variables
1199 real :: slope ! The absolute value of the bottom depth slope across a cell times the cell width [Z ~> m].
1200 real :: I_slope ! The inverse of the normalized slope [Z-1 ~> m-1]
1201 real :: Vol_open ! The cell volume above which it is open [Z ~> m].
1202 integer :: K, nz
1203
1204 nz = gv%ke
1205
1206 slope = abs(dp - dm)
1207 if (slope == 0.0) then
1208 l(1:nz) = 1.0 ; l(nz+1) = 0.0
1209 else
1210 vol_open = 0.5*slope
1211 i_slope = 1.0 / slope
1212
1213 l(nz+1) = 0.0
1214 do k=nz,1,-1
1215 if (vol_below(k) >= vol_open) then ; l(k) = 1.0
1216 else
1217 ! With a uniformly sloping bottom, the calculation of L(K) is the solution of a simple quadratic equation.
1218 l(k) = sqrt(2.0*vol_below(k)*i_slope)
1219 endif
1220 enddo
1221 endif
1222
1223end subroutine find_l_open_uniform_slope
1224
1225!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective)
1226!! using trigonometric expressions. In this case there can be two separate open regions.
1227subroutine find_l_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV)
1228 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1229 real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1230 !! the full horizontal area of a velocity cell [Z ~> m]
1231 real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1232 real, intent(in) :: Dp !< The larger of the two depths at the edge
1233 !! of a velocity cell [Z ~> m]
1234 real, intent(in) :: Dm !< The smaller of the two depths at the edge
1235 !! of a velocity cell [Z ~> m]
1236 real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1237 !! the depth of each interface [nondim]
1238
1239 ! Local variables
1240 real :: crv ! crv is the curvature of the bottom depth across a
1241 ! cell, times the cell width squared [Z ~> m].
1242 real :: crv_3 ! crv/3 [Z ~> m].
1243 real :: slope ! The absolute value of the bottom depth slope across
1244 ! a cell times the cell width [Z ~> m].
1245 ! The following "volumes" have units of vertical heights because they are normalized
1246 ! by the full horizontal area of a velocity cell.
1247 real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1248 real :: Vol_2_reg ! The cell volume above which there are two separate
1249 ! open areas that must be integrated [Z ~> m].
1250 real :: C24_crv ! 24/crv [Z-1 ~> m-1].
1251 real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim].
1252 real :: a2x48_apb3, Iapb ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1]
1253 real :: L0 ! A linear estimate of L appropriate for tiny volumes [nondim].
1254 real :: slope_crv ! The slope divided by the curvature [nondim]
1255 real :: tmp_val_m1_to_p1 ! A temporary variable [nondim]
1256 real, parameter :: C1_3 = 1.0/3.0, c1_12 = 1.0/12.0 ! Rational constants [nondim]
1257 real, parameter :: C2pi_3 = 8.0*atan(1.0)/3.0 ! An irrational constant, 2/3 pi. [nondim]
1258 integer :: K, nz
1259
1260 nz = gv%ke
1261
1262 ! Each cell extends from x=-1/2 to 1/2, and has a topography
1263 ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12.
1264 !crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3
1265 crv_3 = (dp + dm - (2.0*d_vel)) ; crv = 3.0*crv_3
1266 slope = dp - dm
1267
1268 ! Calculate the volume above which the entire cell is open and the volume at which the
1269 ! equation that is solved for L changes because there are two separate open regions.
1270 if (slope >= crv) then
1271 vol_open = d_vel - dm ; vol_2_reg = vol_open
1272 else
1273 slope_crv = slope / crv
1274 vol_open = 0.25*slope*slope_crv + c1_12*crv
1275 vol_2_reg = 0.5*slope_crv**2 * (crv - c1_3*slope)
1276 endif
1277 ! Define some combinations of crv & slope for later use.
1278 c24_crv = 24.0/crv ; iapb = 1.0/(crv+slope)
1279 apb_4a = (slope+crv)/(4.0*crv) ; a2x48_apb3 = (48.0*(crv*crv))*(iapb**3)
1280 ax2_3apb = 2.0*c1_3*crv*iapb
1281
1282 l(nz+1) = 0.0
1283 ! Determine the normalized open length (L) at each interface.
1284 do k=nz,1,-1
1285 if (vol_below(k) >= vol_open) then ! The whole cell is open.
1286 l(k) = 1.0
1287 elseif (vol_below(k) < vol_2_reg) then
1288 ! In this case, there is a contiguous open region and
1289 ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)).
1290 if (a2x48_apb3*vol_below(k) < 1e-8) then ! Could be 1e-7?
1291 ! There is a very good approximation here for massless layers.
1292 !L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0)
1293 l0 = sqrt(2.0*vol_below(k)*iapb) ; l(k) = l0*(1.0 + (ax2_3apb*l0))
1294 else
1295 !L(K) = apb_4a * (1.0 - &
1296 ! 2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3))
1297 l(k) = apb_4a * (1.0 - &
1298 2.0 * cos(c1_3*acos((a2x48_apb3*vol_below(k)) - 1.0) - c2pi_3))
1299 endif
1300 ! To check the answers.
1301 ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K)
1302 else ! There are two separate open regions.
1303 ! vol_below(K) = slope^2/4crv + crv/12 - (crv/12)*(1-L)^2*(1+2L)
1304 ! At the deepest volume, L = slope/crv, at the top L = 1.
1305 ! L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_crv*(Vol_open - vol_below(K))) - C2pi_3)
1306 tmp_val_m1_to_p1 = 1.0 - c24_crv*(vol_open - vol_below(k))
1307 tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
1308 l(k) = 0.5 - cos(c1_3*acos(tmp_val_m1_to_p1) - c2pi_3)
1309 ! To check the answers.
1310 ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol_below(K)
1311 endif
1312 enddo ! k loop to determine L(K) in the concave case
1313
1315
1316
1317
1318!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using
1319!! iterative methods to solve the relevant cubic equations. In this case there can be two separate open regions.
1320subroutine find_l_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV)
1321 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1322 real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1323 !! the full horizontal area of a velocity cell [Z ~> m]
1324 real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1325 real, intent(in) :: Dp !< The larger of the two depths at the edge
1326 !! of a velocity cell [Z ~> m]
1327 real, intent(in) :: Dm !< The smaller of the two depths at the edge
1328 !! of a velocity cell [Z ~> m]
1329 real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1330 !! the depth of each interface [nondim]
1331
1332 ! Local variables
1333 real :: crv ! crv is the curvature of the bottom depth across a
1334 ! cell, times the cell width squared [Z ~> m].
1335 real :: crv_3 ! crv/3 [Z ~> m].
1336 real :: slope ! The absolute value of the bottom depth slope across
1337 ! a cell times the cell width [Z ~> m].
1338
1339 ! The following "volumes" have units of vertical heights because they are normalized
1340 ! by the full horizontal area of a velocity cell.
1341 real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1342 real :: Vol_2_reg ! The cell volume above which there are two separate
1343 ! open areas that must be integrated [Z ~> m].
1344 real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim]
1345 real :: vol_inflect_1 ! The volume at which there is an inflection point in the expression
1346 ! relating L to vol_err when there is a single open region [Z ~> m]
1347 real :: vol_inflect_2 ! The volume at which there is an inflection point in the expression
1348 ! relating L to vol_err when there are two open regions [Z ~> m]
1349
1350 real :: L_inflect_1 ! The value of L that sits at an inflection point in the expression
1351 ! relating L to vol_err when there is a single open region [nondim]
1352 real :: L_inflect_2 ! The value of L that sits at an inflection point in the expression
1353 ! relating L to vol_err when there is are two open regions [nondim]
1354 real :: L_max, L_min ! Maximum and minimum bounds on the solution for L for an interface [nondim]
1355 real :: vol_err ! The difference between the volume below an interface for a given value
1356 ! of L and the target value [Z ~> m]
1357 real :: dVol_dL ! The partial derivative of the volume below with L [Z ~> m]
1358 real :: vol_err_max ! The value of vol_err when L is L_max [Z ~> m]
1359
1360 ! The following combinations of slope and crv are reused across layers, and hence are pre-calculated
1361 ! for efficiency. All are non-negative.
1362 real :: Icrvpslope ! The inverse of the sum of crv and slope [Z-1 ~> m-1]
1363 real :: slope_crv ! The slope divided by the curvature [nondim]
1364 ! These are only used if the slope exceeds or matches the curvature.
1365 real :: smc ! The slope minus the curvature [Z ~> m]
1366 real :: C3c_m_s ! 3 times the curvature minus the slope [Z ~> m]
1367 real :: I_3c_m_s ! The inverse of 3 times the curvature minus the slope [Z-1 ~> m-1]
1368 ! These are only used if the curvature exceeds the slope.
1369 real :: C4_crv ! The inverse of a quarter of the curvature [Z-1 ~> m-1]
1370 real :: sxcms_c ! The slope times the difference between the curvature and slope
1371 ! divided by the curvature [Z ~> m]
1372 real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m]
1373 real :: I_3s_m_c ! The inverse of 3 times the slope minus the curvature [Z-1 ~> m-1]
1374 real :: C3s_m_c ! 3 times the slope minus the curvature [Z ~> m]
1375
1376 real, parameter :: C1_3 = 1.0 / 3.0, c1_12 = 1.0 / 12.0 ! Rational constants [nondim]
1377 integer :: K, nz, itt
1378 integer, parameter :: max_itt = 10
1379
1380 nz = gv%ke
1381
1382 ! Each cell extends from x=-1/2 to 1/2, and has a topography
1383 ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12.
1384
1385 crv_3 = (dp + dm - 2.0*d_vel) ; crv = 3.0*crv_3
1386 slope = dp - dm
1387
1388 ! Calculate the volume above which the entire cell is open and the volume at which the
1389 ! equation that is solved for L changes because there are two separate open regions.
1390 if (slope >= crv) then
1391 vol_open = d_vel - dm ; vol_2_reg = vol_open
1392 l_2_reg = 1.0
1393 if (crv + slope >= 4.0*crv) then
1394 l_inflect_1 = 1.0 ; vol_inflect_1 = vol_open
1395 else
1396 slope_crv = slope / crv
1397 l_inflect_1 = 0.25 + 0.25*slope_crv
1398 vol_inflect_1 = 0.25*c1_12 * ((slope_crv + 1.0)**2 * (slope + crv))
1399 endif
1400 ! Precalculate some combinations of crv & slope for later use.
1401 smc = slope - crv
1402 c3c_m_s = 3.0*crv - slope
1403 if (c3c_m_s > 2.0*smc) i_3c_m_s = 1.0 / c3c_m_s
1404 else
1405 slope_crv = slope / crv
1406 vol_open = 0.25*slope*slope_crv + c1_12*crv
1407 vol_2_reg = 0.5*slope_crv**2 * (crv - c1_3*slope)
1408 l_2_reg = slope_crv
1409
1410 ! The inflection point is useful to know because below the inflection point
1411 ! Newton's method converges monotonically from above and conversely above it.
1412 ! These are the inflection point values of L and vol_below with a single open segment.
1413 vol_inflect_1 = 0.25*c1_12 * ((slope_crv + 1.0)**2 * (slope + crv))
1414 l_inflect_1 = 0.25 + 0.25*slope_crv
1415 ! These are the inflection point values of L and vol_below when there are two open segments.
1416 ! Vol_inflect_2 = Vol_open - 0.125 * crv_3, which is equivalent to:
1417 vol_inflect_2 = 0.25*slope*slope_crv + 0.125*crv_3
1418 l_inflect_2 = 0.5
1419 ! Precalculate some combinations of crv & slope for later use.
1420 c4_crv = 4.0 / crv
1421 slope2_4crv = 0.25 * slope * slope_crv
1422 sxcms_c = slope_crv*(crv - slope)
1423 c3s_m_c = 3.0*slope - crv
1424 if (c3s_m_c > 2.0*sxcms_c) i_3s_m_c = 1.0 / c3s_m_c
1425 endif
1426 ! Define some combinations of crv & slope for later use.
1427 icrvpslope = 1.0 / (crv+slope)
1428
1429 l(nz+1) = 0.0
1430 ! Determine the normalized open length (L) at each interface.
1431 do k=nz,1,-1
1432 if (vol_below(k) >= vol_open) then ! The whole cell is open.
1433 l(k) = 1.0
1434 elseif (vol_below(k) < vol_2_reg) then
1435 ! In this case, there is a single contiguous open region from x=1/2-L to 1/2.
1436 ! Changing the horizontal variable in the expression from D(x) to D(L) gives:
1437 ! x(L) = 1/2 - L
1438 ! D(L) = crv*(0.5 - L)^2 + slope*(0.5 - L) + D_vel - crv/12
1439 ! D(L) = crv*L^2 - crv*L + crv/4 + slope*(1/2 - L) + D_vel - crv/12
1440 ! D(L) = crv*L^2 - (slope+crv)*L + slope/2 + D_vel + crv/6
1441 ! D(0) = slope/2 + D_vel + crv/6 = (Dp - Dm)/2 + D_vel + (Dp + Dm - 2*D_vel)/2 = Dp
1442 ! D(1) = crv - slope - crv + slope/2 + Dvel + crv/6 = D_vel - slope/2 + crv/6 = Dm
1443 !
1444 ! vol_below = integral(y = 0 to L) D(y) dy - L * D(L)
1445 ! = crv/3*L^3 - (slope+crv)/2*L^2 + (slope/2 + D_vel + crv/6)*L -
1446 ! (crv*L^2 - (slope+crv)*L + slope/2 + D_vel + crv/6) * L
1447 ! = -2/3 * crv * L^3 + 1/2 * (slope+crv) * L^2
1448 ! vol_below(K) = 0.5*L(K)**2*(slope + crv_3*(3-4*L(K)))
1449 ! L(K) is between L(K+1) and slope_crv.
1450 l_max = min(l_2_reg, 1.0)
1451 if (vol_below(k) <= vol_inflect_1) l_max = min(l_max, l_inflect_1)
1452
1453 l_min = l(k+1)
1454 if (vol_below(k) >= vol_inflect_1) l_min = max(l_min, l_inflect_1)
1455
1456 ! Ignoring the cubic term gives an under-estimate but is very accurate for near bottom
1457 ! layers, so use this as a potential floor.
1458 if (2.0*vol_below(k)*icrvpslope > l_min**2) l_min = sqrt(2.0*vol_below(k)*icrvpslope)
1459
1460 ! Start with L_min in most cases.
1461 l(k) = l_min
1462
1463 if (vol_below(k) <= vol_inflect_1) then
1464 ! Starting with L_min below L_inflect_1, only the first overshooting iteration of Newton's
1465 ! method needs bounding.
1466 l(k) = l_min
1467 vol_err = 0.5*l(k)**2 * (slope + crv*(1.0 - 4.0*c1_3*l(k))) - vol_below(k)
1468 ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1469 if (vol_err < 0.0) then
1470 dvol_dl = l(k) * (slope + crv*(1.0 - 2.0*l(k)))
1471 if (l(k)*dvol_dl > vol_err + l_max*dvol_dl) then
1472 l(k) = l_max
1473 else
1474 l(k) = l(k) - (vol_err / dvol_dl)
1475 endif
1476
1477 ! Subsequent iterations of Newton's method do not need bounds.
1478 do itt=1,max_itt
1479 vol_err = 0.5*l(k)**2 * (slope + crv*(1.0 - 4.0*c1_3*l(k))) - vol_below(k)
1480 dvol_dl = l(k) * (slope + crv*(1.0 - 2.0*l(k)))
1481 if (abs(vol_err) < max(1.0e-15*l(k), 1.0e-25)*dvol_dl) exit
1482 l(k) = l(k) - (vol_err / dvol_dl)
1483 enddo
1484 endif
1485 else ! (vol_below(K) > vol_inflect_1)
1486 ! Iteration from below converges monotonically, but we need to deal with the case where we are
1487 ! close to the peak of the topography and Newton's method mimics the convergence of bisection.
1488
1489 ! Evaluate the error when L(K) = L_min as a possible first guess.
1490 l(k) = l_min
1491 vol_err = 0.5*l(k)**2 * (slope + crv*(1.0 - 4.0*c1_3*l(k))) - vol_below(k)
1492 ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1493 if (vol_err < 0.0) then
1494
1495 ! These two upper estimates deal with the possibility that this point may be near
1496 ! the upper extrema, where the error term might be approximately parabolic and
1497 ! Newton's method would converge slowly like simple bisection.
1498 if (slope < crv) then
1499 ! if ((L_2_reg - L_min)*(3.0*slope - crv) > 2.0*slope_crv*(crv-slope)) then
1500 if ((l_2_reg - l_min)*c3s_m_c > 2.0*sxcms_c) then
1501 ! There is a decent upper estimate of L from the approximate quadratic equation found
1502 ! by examining the error expressions at L ~= L_2_reg and ignoring the cubic term.
1503 l_max = (slope_crv*(2.0*slope) - sqrt(sxcms_c**2 + &
1504 2.0*c3s_m_c*(vol_2_reg - vol_below(k))) ) * i_3s_m_c
1505 ! The line above is equivalent to:
1506 ! L_max = (slope_crv*(2.0*slope) - sqrt(slope_crv**2*(crv-slope)**2 + &
1507 ! 2.0*(3.0*slope - crv)*(Vol_2_reg - vol_below(K))) ) / &
1508 ! (3.0*slope - crv)
1509 else
1510 l_max = slope_crv
1511 endif
1512 else ! (slope >= crv)
1513 if ((1.0 - l_min)*c3c_m_s > 2.0*smc) then
1514 ! There is a decent upper estimate of L from the approximate quadratic equation found
1515 ! by examining the error expressions at L ~= 1 and ignoring the cubic term.
1516 l_max = ( 2.0*crv - sqrt(smc**2 + 2.0*c3c_m_s * (vol_open - vol_below(k))) ) * i_3c_m_s
1517 ! The line above is equivalent to:
1518 ! L_max = ( 2.0*crv - sqrt((slope - crv)**2 + 2.0*(3.0*crv - slope) * (Vol_open - vol_below(K))) ) / &
1519 ! (3.0*crv - slope)
1520 else
1521 l_max = 1.0
1522 endif
1523 endif
1524 vol_err_max = 0.5*l_max**2 * (slope + crv*(1.0 - 4.0*c1_3*l_max)) - vol_below(k)
1525 ! if (Vol_err_max < 0.0) call MOM_error(FATAL, &
1526 ! "Vol_err_max should never be negative in find_L_open_concave_iterative.")
1527 if ((vol_err_max < abs(vol_err)) .and. (l_max < 1.0)) then
1528 ! Start with 1 bounded Newton's method step from L_max
1529 dvol_dl = l_max * (slope + crv*(1.0 - 2.0*l_max))
1530 l(k) = max(l_min, l_max - (vol_err_max / dvol_dl) )
1531 ! else ! Could use the fact that Vol_err is known to take an iteration?
1532 endif
1533
1534 ! Subsequent iterations of Newton's method do not need bounds.
1535 do itt=1,max_itt
1536 vol_err = 0.5*l(k)**2 * (slope + crv*(1.0 - 4.0*c1_3*l(k))) - vol_below(k)
1537 dvol_dl = l(k) * (slope + crv*(1.0 - 2.0*l(k)))
1538 if (abs(vol_err) < max(1.0e-15*l(k), 1.0e-25)*dvol_dl) exit
1539 l(k) = l(k) - (vol_err / dvol_dl)
1540 enddo
1541 endif
1542
1543 endif
1544
1545 ! To check the answers.
1546 ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K)
1547 else ! There are two separate open regions.
1548 ! vol_below(K) = slope^2/(4*crv) + crv/12 - (crv/12)*(1-L)^2*(1+2L)
1549 ! At the deepest volume, L = slope/crv, at the top L = 1.
1550
1551 ! To check the answers.
1552 ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol_below(K)
1553 ! or equivalently:
1554 ! Vol_err = Vol_open - 0.25*crv_3*(3.0-2.0*(1.0-L(K))) * (1.0-L(K))**2 - vol_below(K)
1555 ! ! Note that: Vol_open = 0.25*slope*slope_crv + C1_12*crv
1556 ! Vol_err = 0.25*slope*slope_crv + 0.25*crv_3*( 1.0 - (1.0 + 2.0*L(K)) * (1.0-L(K))**2 ) - vol_below(K)
1557 ! Vol_err = 0.25*crv_3*L(K)**2*( 3.0 - 2.0*L(K) ) + 0.25*slope*slope_crv - vol_below(K)
1558
1559 ! Derivation of the L_max limit below:
1560 ! Vol_open - vol_below(K) = 0.25*crv_3*(3.0-2.0*(1.0-L(K))) * (1.0-L(K))**2
1561 ! (3.0-2.0*(1.0-L(K))) * (1.0-L(K))**2 = (Vol_open - vol_below(K)) / (0.25*crv_3)
1562 ! When 1-L(K) << 1:
1563 ! 3.0 * (1.0-L_max)**2 = (Vol_open - vol_below(K)) / (0.25*crv_3)
1564 ! (1.0-L_max)**2 = (Vol_open - vol_below(K)) / (0.25*crv)
1565
1566 ! Derivation of the L_min limit below:
1567 ! Vol_err = 0.25*crv_3*L(K)**2*( 3.0 - 2.0*L(K) ) + 0.25*slope*slope_crv - vol_below(K)
1568 ! crv*L(K)**2*( 1.0 - 2.0*C1_3*L(K) ) = 4.0*vol_below(K) - slope*slope_crv
1569 ! When L(K) << 1:
1570 ! crv*L_min**2 = 4.0*vol_below(K) - slope*slope_crv
1571 ! L_min = sqrt((4.0*vol_below(K) - slope*slope_crv)/crv)
1572 ! Noting that L(K) >= slope_crv, when L(K)-slope_crv << 1:
1573 ! (crv + 2.0*C1_3*slope)*L_min**2 = 4.0*vol_below(K) - slope*slope_crv
1574 ! L_min = sqrt((4.0*vol_below(K) - slope*slope_crv)/(crv + 2.0*C1_3*slope))
1575
1576 if (vol_below(k) <= vol_inflect_2) then
1577 ! Newton's Method would converge monotonically from above, but overshoot from below.
1578 l_min = max(l(k+1), l_2_reg) ! L_2_reg = slope_crv
1579 ! This under-estimate of L(K) is accurate for L ~= slope_crv:
1580 if ((4.0*vol_below(k) - slope*slope_crv) > (crv + 2.0*c1_3*slope)*l_min**2) &
1581 l_min = max(l_min, sqrt((4.0*vol_below(k) - slope*slope_crv) / (crv + 2.0*c1_3*slope)))
1582 l_max = 0.5 ! = L_inflect_2
1583
1584 ! Starting with L_min below L_inflect_2, only the first overshooting iteration of Newton's
1585 ! method needs bounding.
1586 l(k) = l_min
1587 vol_err = crv_3*l(k)**2*( 0.75 - 0.5*l(k) ) + (slope2_4crv - vol_below(k))
1588
1589 ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1590 if (vol_err < 0.0) then
1591 dvol_dl = 0.5*crv * (l(k) * (1.0 - l(k)))
1592 if (l(k)*dvol_dl >= vol_err + l_max*dvol_dl) then
1593 l(k) = l_max
1594 else
1595 l(k) = l(k) - (vol_err / dvol_dl)
1596 endif
1597 ! Subsequent iterations of Newton's method do not need bounds.
1598 do itt=1,max_itt
1599 vol_err = crv_3 * (l(k)**2 * (0.75 - 0.5*l(k))) + (slope2_4crv - vol_below(k))
1600 dvol_dl = 0.5*crv * (l(k)*(1.0 - l(k)))
1601 if (abs(vol_err) < max(1.0e-15*l(k), 1.0e-25)*dvol_dl) exit
1602 l(k) = l(k) - (vol_err / dvol_dl)
1603 enddo
1604 endif
1605 else ! (vol_below(K) > Vol_inflect_2)
1606 ! Newton's Method would converge monotonically from below, but overshoots from above, and
1607 ! we may need to deal with the case where we are close to the peak of the topography.
1608 l_min = max(l(k+1), 0.5)
1609 l(k) = l_min
1610
1611 vol_err = crv_3 * (l(k)**2 * ( 0.75 - 0.5*l(k))) + (slope2_4crv - vol_below(k))
1612 ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L(k) is already the best solution.
1613 if (vol_err < 0.0) then
1614 ! This over-estimate of L(K) is accurate for L ~= 1:
1615 l_max = 1.0 - sqrt( (vol_open - vol_below(k)) * c4_crv )
1616 vol_err_max = crv_3 * (l_max**2 * ( 0.75 - 0.5*l_max)) + (slope2_4crv - vol_below(k))
1617 ! if (Vol_err_max < 0.0) call MOM_error(FATAL, &
1618 ! "Vol_err_max should never be negative in find_L_open_concave_iterative.")
1619 if ((vol_err_max < abs(vol_err)) .and. (l_max < 1.0)) then
1620 ! Start with 1 bounded Newton's method step from L_max
1621 dvol_dl = 0.5*crv * (l_max * (1.0 - l_max))
1622 l(k) = max(l_min, l_max - (vol_err_max / dvol_dl) )
1623 ! else ! Could use the fact that Vol_err is known to take an iteration?
1624 endif
1625
1626 ! Subsequent iterations of Newton's method do not need bounds.
1627 do itt=1,max_itt
1628 vol_err = crv_3 * (l(k)**2 * ( 0.75 - 0.5*l(k))) + (slope2_4crv - vol_below(k))
1629 dvol_dl = 0.5*crv * (l(k) * (1.0 - l(k)))
1630 if (abs(vol_err) < max(1.0e-15*l(k), 1.0e-25)*dvol_dl) exit
1631 l(k) = l(k) - (vol_err / dvol_dl)
1632 enddo
1633 endif
1634 endif
1635
1636 endif
1637 enddo ! k loop to determine L(K) in the concave case
1638
1639end subroutine find_l_open_concave_iterative
1640
1641
1642
1643!> Test the validity the normalized open lengths of each interface for concave bathymetry (from the ocean perspective)
1644!! by evaluating and returing the relevant cubic equations.
1645subroutine test_l_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err, GV)
1646 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1647 real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1648 !! the full horizontal area of a velocity cell [Z ~> m]
1649 real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1650 real, intent(in) :: Dp !< The larger of the two depths at the edge
1651 !! of a velocity cell [Z ~> m]
1652 real, intent(in) :: Dm !< The smaller of the two depths at the edge
1653 !! of a velocity cell [Z ~> m]
1654 real, dimension(SZK_(GV)+1), intent(in) :: L !< The fraction of the full cell width that is open at
1655 !! the depth of each interface [nondim]
1656 real, dimension(SZK_(GV)+1), intent(out) :: vol_err !< The difference between vol_below and the
1657 !! value obtained from using L in the cubic equation [Z ~> m]
1658
1659 ! Local variables
1660 real :: crv ! crv is the curvature of the bottom depth across a
1661 ! cell, times the cell width squared [Z ~> m].
1662 real :: crv_3 ! crv/3 [Z ~> m].
1663 real :: slope ! The absolute value of the bottom depth slope across
1664 ! a cell times the cell width [Z ~> m].
1665
1666 ! The following "volumes" have units of vertical heights because they are normalized
1667 ! by the full horizontal area of a velocity cell.
1668 real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1669 real :: Vol_2_reg ! The cell volume above which there are two separate
1670 ! open areas that must be integrated [Z ~> m].
1671 real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim]
1672
1673 ! The following combinations of slope and crv are reused across layers, and hence are pre-calculated
1674 ! for efficiency. All are non-negative.
1675 real :: slope_crv ! The slope divided by the curvature [nondim]
1676 ! These are only used if the curvature exceeds the slope.
1677 real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m]
1678
1679 real, parameter :: C1_3 = 1.0 / 3.0, c1_12 = 1.0 / 12.0 ! Rational constants [nondim]
1680 integer :: K, nz
1681
1682 nz = gv%ke
1683
1684 ! Each cell extends from x=-1/2 to 1/2, and has a topography
1685 ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12.
1686
1687 crv_3 = (dp + dm - 2.0*d_vel) ; crv = 3.0*crv_3
1688 slope = dp - dm
1689
1690 ! Calculate the volume above which the entire cell is open and the volume at which the
1691 ! equation that is solved for L changes because there are two separate open regions.
1692 if (slope >= crv) then
1693 vol_open = d_vel - dm ; vol_2_reg = vol_open
1694 l_2_reg = 1.0
1695 if (crv + slope >= 4.0*crv) then
1696 slope_crv = 1.0
1697 else
1698 slope_crv = slope / crv
1699 endif
1700 else
1701 slope_crv = slope / crv
1702 vol_open = 0.25*slope*slope_crv + c1_12*crv
1703 vol_2_reg = 0.5*slope_crv**2 * (crv - c1_3*slope)
1704 l_2_reg = slope_crv
1705 endif
1706 slope2_4crv = 0.25 * slope * slope_crv
1707
1708 ! Determine the volume error based on the normalized open length (L) at each interface.
1709 vol_err(nz+1) = 0.0
1710 do k=nz,1,-1
1711 if (l(k) >= 1.0) then
1712 vol_err(k) = max(vol_open - vol_below(k), 0.0)
1713 elseif (l(k) <= l_2_reg) then
1714 vol_err(k) = 0.5*l(k)**2 * (slope + crv*(1.0 - 4.0*c1_3*l(k))) - vol_below(k)
1715 else ! There are two separate open regions.
1716 vol_err(k) = crv_3 * (l(k)**2 * ( 0.75 - 0.5*l(k))) + (slope2_4crv - vol_below(k))
1717 endif
1718 enddo ! k loop to determine L(K) in the concave case
1719
1720end subroutine test_l_open_concave
1721
1722
1723!> Determine the normalized open length of each interface for convex bathymetry (from the ocean
1724!! perspective) using Newton's method iterations. In this case there is a single open region
1725!! with the minimum depth at one edge of the cell.
1726subroutine find_l_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS)
1727 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1728 real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1729 !! the full horizontal area of a velocity cell [Z ~> m]
1730 real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1731 real, intent(in) :: Dp !< The larger of the two depths at the edge
1732 !! of a velocity cell [Z ~> m]
1733 real, intent(in) :: Dm !< The smaller of the two depths at the edge
1734 !! of a velocity cell [Z ~> m]
1735 real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1736 !! the depth of each interface [nondim]
1737 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1738 type(set_visc_cs), intent(in) :: CS !< The control structure returned by a previous
1739 !! call to set_visc_init.
1740
1741 ! Local variables
1742 real :: crv ! crv is the curvature of the bottom depth across a
1743 ! cell, times the cell width squared [Z ~> m].
1744 real :: crv_3 ! crv/3 [Z ~> m].
1745 real :: slope ! The absolute value of the bottom depth slope across
1746 ! a cell times the cell width [Z ~> m].
1747 ! All of the following "volumes" have units of vertical heights because they are normalized
1748 ! by the full horizontal area of a velocity cell.
1749 real :: Vol_err ! The error in the volume with the latest estimate of
1750 ! L, or the error for the interface below [Z ~> m].
1751 real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m].
1752 real :: Vol_tol ! A volume error tolerance [Z ~> m].
1753 real :: Vol_open ! The cell volume above which the face is fully open [Z ~> m].
1754 real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct
1755 ! solution of a cubic equation for L.
1756 real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m]
1757 real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m]
1758 real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m].
1759 real :: dVol ! vol - Vol_0 [Z ~> m].
1760 real :: dV_dL2 ! The partial derivative of volume with L squared
1761 ! evaluated at L=L0 [Z ~> m].
1762 real :: L_direct ! The value of L above volume Vol_direct [nondim].
1763 real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim].
1764 real :: L0 ! The value of L above volume Vol_0 [nondim].
1765 real :: Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1]
1766 real :: C24_crv ! 24/crv [Z-1 ~> m-1].
1767 real :: curv_tol ! Numerator of curvature cubed, used to estimate
1768 ! accuracy of a single L(:) Newton iteration [Z5 ~> m5]
1769 real, parameter :: C1_3 = 1.0/3.0, c1_6 = 1.0/6.0 ! Rational constants [nondim]
1770 logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration
1771 integer :: K, nz, itt, maxitt=20
1772
1773 nz = gv%ke
1774
1775 ! Each cell extends from x=-1/2 to 1/2, and has a topography
1776 ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12.
1777 crv_3 = (dp + dm - 2.0*d_vel) ; crv = 3.0*crv_3
1778 slope = dp - dm
1779
1780 ! Calculate the volume above which the entire cell is open and the volume at which the
1781 ! equation that is solved for L changes because there is a direct solution.
1782 vol_open = d_vel - dm
1783 if (slope >= -crv) then
1784 iapb = 1.0e30*us%Z_to_m ; if (slope+crv /= 0.0) iapb = 1.0/(crv+slope)
1785 vol_direct = 0.0 ; l_direct = 0.0 ; c24_crv = 0.0
1786 else
1787 c24_crv = 24.0/crv ; iapb = 1.0/(crv+slope)
1788 l_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0
1789 vol_direct = -c1_6*crv*l_direct**3
1790 endif
1791 ibma_2 = 2.0 / (slope - crv)
1792
1793 if (cs%answer_date < 20190101) vol_quit = (0.9*gv%Angstrom_Z + gv%dZ_subroundoff)
1794
1795 l(nz+1) = 0.0 ; vol_err = 0.0
1796 ! Determine the normalized open length (L) at each interface.
1797 do k=nz,1,-1
1798 if (vol_below(k) >= vol_open) then
1799 l(k) = 1.0
1800 elseif (vol_below(k) <= vol_direct) then
1801 ! Both edges of the cell are bounded by walls.
1802 if (cs%answer_date < 20260704) then
1803 l(k) = (-0.25*c24_crv*vol_below(k))**c1_3
1804 else
1805 l(k) = cuberoot(-0.25*c24_crv*vol_below(k))
1806 endif
1807 else
1808 ! x_R is at 1/2 but x_L is in the interior & L is found by iteratively solving
1809 ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L))
1810
1811 ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + crv_3*(3.0-4.0*L(K+1))) - vol_below(K+1)
1812 ! Change to ...
1813 ! if (min(vol_below(K+1) + Vol_err, vol_below(K)) <= Vol_direct) then ?
1814 if (vol_below(k+1) + vol_err <= vol_direct) then
1815 l0 = l_direct ; vol_0 = vol_direct
1816 else
1817 l0 = l(k+1) ; vol_0 = vol_below(k+1) + vol_err
1818 ! Change to Vol_0 = min(vol_below(K+1) + Vol_err, vol_below(K)) ?
1819 endif
1820
1821 ! Try a relatively simple solution that usually works well
1822 ! for massless layers.
1823 dv_dl2 = 0.5*(slope+crv) - crv*l0 ; dvol = (vol_below(k)-vol_0)
1824 ! dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = max(vol_below(K)-Vol_0, 0.0)
1825
1826 use_l0 = .false.
1827 do_one_l_iter = .false.
1828 if (cs%answer_date < 20190101) then
1829 curv_tol = gv%Angstrom_Z*dv_dl2**2 &
1830 * (0.25 * dv_dl2 * gv%Angstrom_Z - crv * l0 * dvol)
1831 do_one_l_iter = (crv * crv * dvol**3) < curv_tol
1832 else
1833 ! The following code is more robust when GV%Angstrom_H=0, but
1834 ! it changes answers.
1835 use_l0 = (dvol <= 0.)
1836
1837 vol_tol = max(0.5 * gv%Angstrom_Z + gv%dZ_subroundoff, 1e-14 * vol_below(k))
1838 vol_quit = max(0.9 * gv%Angstrom_Z + gv%dZ_subroundoff, 1e-14 * vol_below(k))
1839
1840 curv_tol = vol_tol * dv_dl2**2 &
1841 * (dv_dl2 * vol_tol - 2.0 * crv * l0 * dvol)
1842 do_one_l_iter = (crv * crv * dvol**3) < curv_tol
1843 endif
1844
1845 if (use_l0) then
1846 l(k) = l0
1847 vol_err = 0.5*(l(k)*l(k))*(slope + crv_3*(3.0-4.0*l(k))) - vol_below(k)
1848 elseif (do_one_l_iter) then
1849 ! One iteration of Newton's method should give an estimate
1850 ! that is accurate to within Vol_tol.
1851 l(k) = sqrt(l0*l0 + dvol / dv_dl2)
1852 vol_err = 0.5*(l(k)*l(k))*(slope + crv_3*(3.0-4.0*l(k))) - vol_below(k)
1853 else
1854 if (dv_dl2*(1.0-l0*l0) < dvol + &
1855 dv_dl2 * (vol_open - vol_below(k))*ibma_2) then
1856 l_max = sqrt(1.0 - (vol_open - vol_below(k))*ibma_2)
1857 else
1858 l_max = sqrt(l0*l0 + dvol / dv_dl2)
1859 endif
1860 l_min = sqrt(l0*l0 + dvol / (0.5*(slope+crv) - crv*l_max))
1861
1862 vol_err_min = 0.5*(l_min**2)*(slope + crv_3*(3.0-4.0*l_min)) - vol_below(k)
1863 vol_err_max = 0.5*(l_max**2)*(slope + crv_3*(3.0-4.0*l_max)) - vol_below(k)
1864 ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
1865 if (abs(vol_err_min) <= vol_quit) then
1866 l(k) = l_min ; vol_err = vol_err_min
1867 else
1868 l(k) = sqrt((l_min**2*vol_err_max - l_max**2*vol_err_min) / &
1869 (vol_err_max - vol_err_min))
1870 do itt=1,maxitt
1871 vol_err = 0.5*(l(k)*l(k))*(slope + crv_3*(3.0-4.0*l(k))) - vol_below(k)
1872 if (abs(vol_err) <= vol_quit) exit
1873 ! Take a Newton's method iteration. This equation has proven
1874 ! robust enough not to need bracketing.
1875 l(k) = l(k) - vol_err / (l(k)* (slope + crv - 2.0*crv*l(k)))
1876 ! This would be a Newton's method iteration for L^2:
1877 ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+crv) - crv*L(K)))
1878 enddo
1879 endif ! end of iterative solver
1880 endif ! end of 1-boundary alternatives.
1881 endif ! end of 0, 1- and 2- boundary cases.
1882 enddo ! k loop to determine L(K) in the convex case
1883
1884end subroutine find_l_open_convex
1885
1886!> This subroutine finds a thickness-weighted value of v at the u-points.
1887function set_v_at_u(v, h, G, GV, i, j, k, mask2dCv, OBC)
1888 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1889 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1890 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1891 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1892 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1893 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1894 integer, intent(in) :: i !< The i-index of the u-location to work on.
1895 integer, intent(in) :: j !< The j-index of the u-location to work on.
1896 integer, intent(in) :: k !< The k-index of the u-location to work on.
1897 real, dimension(SZI_(G),SZJB_(G)),&
1898 intent(in) :: mask2dcv !< A multiplicative mask of the v-points [nondim]
1899 type(ocean_obc_type), pointer :: obc !< A pointer to an open boundary condition structure
1900 real :: set_v_at_u !< The return value of v at u points points in the
1901 !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
1902
1903 ! This subroutine finds a thickness-weighted value of v at the u-points.
1904 real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
1905 real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
1906 integer :: i0, j0, i1, j1
1907
1908 do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
1909 hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
1910 enddo ; enddo
1911
1912 if (associated(obc)) then ; if (obc%number_of_segments > 0) then
1913 do j0 = -1,0 ; do i0 = 0,1 ; if (obc%segnum_v(i+i0,j+j0) /= 0) then
1914 i1 = i+i0 ; j1 = j+j0
1915 if (obc%segnum_v(i1,j1) > 0) then ! OBC_DIRECTION_N
1916 hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
1917 elseif (obc%segnum_v(i1,j1) < 0) then ! OBC_DIRECTION_S
1918 hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
1919 endif
1920 endif ; enddo ; enddo
1921 endif ; endif
1922
1923 hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
1924 set_v_at_u = 0.0
1925 if (hwt_tot > 0.0) set_v_at_u = &
1926 (((hwt(0,0) * v(i,j,k)) + (hwt(1,-1) * v(i+1,j-1,k))) + &
1927 ((hwt(1,0) * v(i+1,j,k)) + (hwt(0,-1) * v(i,j-1,k)))) / hwt_tot
1928
1929end function set_v_at_u
1930
1931!> This subroutine finds a thickness-weighted value of u at the v-points.
1932function set_u_at_v(u, h, G, GV, i, j, k, mask2dCu, OBC)
1933 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1934 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1935 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1936 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
1937 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1938 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1939 integer, intent(in) :: i !< The i-index of the u-location to work on.
1940 integer, intent(in) :: j !< The j-index of the u-location to work on.
1941 integer, intent(in) :: k !< The k-index of the u-location to work on.
1942 real, dimension(SZIB_(G),SZJ_(G)), &
1943 intent(in) :: mask2dcu !< A multiplicative mask of the u-points [nondim]
1944 type(ocean_obc_type), pointer :: obc !< A pointer to an open boundary condition structure
1945 real :: set_u_at_v !< The return value of u at v points in the
1946 !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
1947
1948 ! This subroutine finds a thickness-weighted value of u at the v-points.
1949 real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
1950 real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
1951 integer :: i0, j0, i1, j1
1952
1953 do j0 = 0,1 ; do i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
1954 hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
1955 enddo ; enddo
1956
1957 if (associated(obc)) then ; if (obc%number_of_segments > 0) then
1958 do j0 = 0,1 ; do i0 = -1,0 ; if ((obc%segnum_u(i+i0,j+j0) /= 0)) then
1959 i1 = i+i0 ; j1 = j+j0
1960 if (obc%segnum_u(i1,j1) > 0) then ! OBC_DIRECTION_E
1961 hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
1962 elseif (obc%segnum_u(i1,j1) < 0) then ! OBC_DIRECTION_W
1963 hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
1964 endif
1965 endif ; enddo ; enddo
1966 endif ; endif
1967
1968 hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
1969 set_u_at_v = 0.0
1970 if (hwt_tot > 0.0) set_u_at_v = &
1971 (((hwt(0,0) * u(i,j,k)) + (hwt(-1,1) * u(i-1,j+1,k))) + &
1972 ((hwt(-1,0) * u(i-1,j,k)) + (hwt(0,1) * u(i,j+1,k)))) / hwt_tot
1973
1974end function set_u_at_v
1975
1976!> Calculates the thickness of the surface boundary layer for applying an elevated viscosity.
1977!!
1978!! A bulk Richardson criterion or the thickness of the topmost NKML layers (with a bulk mixed layer)
1979!! are currently used. The thicknesses are given in terms of fractional layers, so that this
1980!! thickness will move as the thickness of the topmost layers change.
1981subroutine set_viscous_ml(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
1982 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
1983 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1984 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1985 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1986 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
1987 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1988 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
1989 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1990 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1991 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
1992 !! thermodynamic fields. Absent fields have
1993 !! NULL pointers.
1994 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1995 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1996 !! related fields.
1997 real, intent(in) :: dt !< Time increment [T ~> s].
1998 type(set_visc_cs), intent(inout) :: cs !< The control structure returned by a previous
1999 !! call to set_visc_init.
2000
2001 ! Local variables
2002 real, dimension(SZIB_(G)) :: &
2003 htot, & ! The total thickness of the layers that are within the
2004 ! surface mixed layer [H ~> m or kg m-2].
2005 dztot, & ! The distance from the surface to the bottom of the layers that are
2006 ! within the surface mixed layer [Z ~> m]
2007 thtot, & ! The integrated temperature of layers that are within the
2008 ! surface mixed layer [H C ~> m degC or kg degC m-2].
2009 shtot, & ! The integrated salt of layers that are within the
2010 ! surface mixed layer [H S ~> m ppt or kg ppt m-2].
2011 spv_htot, & ! Running sum of thickness times specific volume [H R-1 ~> m4 kg-1 or m]
2012 rhtot, & ! The integrated density of layers that are within the surface mixed layer
2013 ! [H R ~> kg m-2 or kg2 m-5]. Rhtot is only used if no
2014 ! equation of state is used.
2015 uhtot, & ! The depth integrated zonal velocity within the surface
2016 ! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2017 vhtot, & ! The depth integrated meridional velocity within the surface
2018 ! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2019 idecay_len_tke, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
2020 dr_dt, & ! Partial derivative of the density at the base of layer nkml
2021 ! (roughly the base of the mixed layer) with temperature [R C-1 ~> kg m-3 degC-1].
2022 dr_ds, & ! Partial derivative of the density at the base of layer nkml
2023 ! (roughly the base of the mixed layer) with salinity [R S-1 ~> kg m-3 ppt-1].
2024 dspv_dt, & ! Partial derivative of the specific volume at the base of layer nkml
2025 ! (roughly the base of the mixed layer) with temperature [R-1 C-1 ~> m3 kg-1 degC-1].
2026 dspv_ds, & ! Partial derivative of the specific volume at the base of layer nkml
2027 ! (roughly the base of the mixed layer) with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
2028 ustar, & ! The surface friction velocity under ice shelves [H T-1 ~> m s-1 or kg m-2 s-1].
2029 press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
2030 t_eos, & ! The potential temperature at which dR_dT and dR_dS are evaluated [C ~> degC]
2031 s_eos ! The salinity at which dR_dT and dR_dS are evaluated [S ~> ppt].
2032 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Height change across layers [Z ~> m]
2033 real, dimension(SZIB_(G),SZJ_(G)) :: &
2034 mask_u ! A mask that disables any contributions from u points that
2035 ! are land or past open boundary conditions [nondim], 0 or 1.
2036 real, dimension(SZI_(G),SZJB_(G)) :: &
2037 mask_v ! A mask that disables any contributions from v points that
2038 ! are land or past open boundary conditions [nondim], 0 or 1.
2039 real :: u_star_2d(szi_(g),szj_(g)) ! The wind friction velocity in thickness-based units,
2040 ! calculated using the Boussinesq reference density or the time-evolving
2041 ! surface density in non-Boussinesq mode [H T-1 ~> m s-1 or kg m-2 s-1]
2042 real :: h_at_vel(szib_(g),szk_(gv))! Layer thickness at velocity points,
2043 ! using an upwind-biased second order accurate estimate based
2044 ! on the previous velocity direction [H ~> m or kg m-2].
2045 real :: dz_at_vel(szib_(g),szk_(gv)) ! Vertical extent of a layer at velocity points,
2046 ! using an upwind-biased second order accurate estimate based
2047 ! on the previous velocity direction [Z ~> m].
2048 integer :: k_massive(szib_(g)) ! The k-index of the deepest layer yet found
2049 ! that has more than h_tiny thickness and will be in the
2050 ! viscous mixed layer.
2051 real :: uh2 ! The squared magnitude of the difference between the velocity
2052 ! integrated through the mixed layer and the velocity of the
2053 ! interior layer layer times the depth of the mixed layer
2054 ! [H2 L2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
2055 real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
2056 real :: hwtot ! Sum of the thicknesses used to calculate
2057 ! the near-bottom velocity magnitude [H ~> m or kg m-2].
2058 real :: hutot ! Running sum of thicknesses times the velocity
2059 ! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2060 real :: hweight ! The thickness of a layer that is within Hbbl
2061 ! of the bottom [H ~> m or kg m-2].
2062 real :: tbl_thick ! The thickness of the top boundary layer [Z ~> m].
2063
2064 real :: hlay ! The layer thickness at velocity points [H ~> m or kg m-2].
2065 real :: i_2hlay ! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].
2066 real :: t_lay ! The layer temperature at velocity points [C ~> degC].
2067 real :: s_lay ! The layer salinity at velocity points [S ~> ppt].
2068 real :: rlay ! The layer potential density at velocity points [R ~> kg m-3].
2069 real :: rlb ! The potential density of the layer below [R ~> kg m-3].
2070 real :: v_at_u ! The meridional velocity at a zonal velocity point [L T-1 ~> m s-1].
2071 real :: u_at_v ! The zonal velocity at a meridional velocity point [L T-1 ~> m s-1].
2072 real :: ghprime ! The mixed-layer internal gravity wave speed squared, based
2073 ! on the mixed layer thickness and density difference across
2074 ! the base of the mixed layer [L2 T-2 ~> m2 s-2].
2075 real :: ribulk ! The bulk Richardson number below which water is in the
2076 ! viscous mixed layer, including reduction for turbulent decay [nondim]
2077 real :: dt_rho0 ! The time step divided by the conversion from the layer
2078 ! thickness to layer mass [T H Z-1 R-1 ~> s m3 kg-1 or s].
2079 real :: g_h_rho0 ! The gravitational acceleration times the conversion from H to m divided
2080 ! by the mean density [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
2081 real :: ustarsq ! 400 times the square of ustar, times
2082 ! Rho0 divided by G_Earth and the conversion
2083 ! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
2084 real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
2085 real :: cdrag_sqrt_h ! Square root of the drag coefficient, times a unit conversion
2086 ! factor from lateral lengths to layer thicknesses [H L-1 ~> nondim or kg m-3].
2087 real :: cdrag_sqrt_h_rl ! Square root of the drag coefficient, times a unit conversion factor from
2088 ! density times lateral lengths to layer thicknesses [H L-1 R-1 ~> m3 kg-1 or nondim]
2089 real :: oldfn ! The integrated energy required to
2090 ! entrain up to the bottom of the layer,
2091 ! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
2092 real :: dfn ! The increment in oldfn for entraining
2093 ! the layer [H R ~> kg m-2 or kg2 m-5].
2094 real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim]
2095 real :: dh ! The increment in layer thickness from the present layer [H ~> m or kg m-2].
2096 real :: ddz ! The increment in height change from the present layer [Z ~> m].
2097 real :: u2_bg(szib_(g)) ! The square of an assumed background velocity, for
2098 ! calculating the mean magnitude near the top for use in
2099 ! the quadratic surface drag [L2 T-2 ~> m2 s-2].
2100 real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
2101 ! h_tiny can not be the deepest in the viscous mixed layer.
2102 real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
2103 real :: u_star ! The friction velocity at velocity points [H T-1 ~> m s-1 or kg m-2 s-1].
2104 real :: h_neglect ! A thickness that is so small it is usually lost
2105 ! in roundoff and can be neglected [H ~> m or kg m-2].
2106 real :: dz_neglect ! A vertical distance that is so small it is usually lost
2107 ! in roundoff and can be neglected [Z ~> m].
2108 real :: rho0x400_g ! 400*Rho0/G_Earth, times unit conversion factors
2109 ! [R T2 H-1 ~> kg s2 m-4 or s2 m-1].
2110 ! The 400 is a constant proposed by Killworth and Edwards, 1999.
2111 real :: ustar1 ! ustar [H T-1 ~> m s-1 or kg m-2 s-1]
2112 real :: h2f2 ! (h*2*f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
2113 logical :: use_eos, do_any, do_any_shelf, do_i(szib_(g))
2114 logical :: nonbous_ml ! If true, use the non-Boussinesq form of some energy and
2115 ! stratification calculations.
2116 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, k2, nkmb, nkml, n
2117 type(ocean_obc_type), pointer :: obc => null()
2118
2119 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2120 isq = g%isc-1 ; ieq = g%IecB ; jsq = g%jsc-1 ; jeq = g%JecB
2121 nkmb = gv%nk_rho_varies ; nkml = gv%nkml
2122
2123 if (.not.cs%initialized) call mom_error(fatal,"MOM_set_viscosity(visc_ML): "//&
2124 "Module must be initialized before it is used.")
2125
2126 if (.not.(cs%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
2127 associated(forces%frac_shelf_v)) ) return
2128
2129 rho0x400_g = 400.0*(gv%H_to_RZ / gv%g_Earth_Z_T2)
2130 cdrag_sqrt = sqrt(cs%cdrag)
2131 cdrag_sqrt_h = cdrag_sqrt * us%L_to_m * gv%m_to_H
2132 cdrag_sqrt_h_rl = cdrag_sqrt * us%L_to_Z * gv%RZ_to_H
2133
2134 obc => cs%OBC
2135 use_eos = associated(tv%eqn_of_state)
2136 nonbous_ml = allocated(tv%SpV_avg)
2137 dt_rho0 = dt / gv%H_to_RZ
2138 h_neglect = gv%H_subroundoff
2139 h_tiny = 2.0*gv%Angstrom_H + h_neglect
2140 dz_neglect = gv%dZ_subroundoff
2141 g_h_rho0 = (gv%g_Earth*gv%H_to_Z) / (gv%Rho0)
2142
2143 if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) &
2144 call mom_error(fatal, "set_viscous_ML: one of forces%frac_shelf_u and "//&
2145 "forces%frac_shelf_v is associated, but the other is not.")
2146
2147 ! Extract the friction velocity from the forcing type.
2148 call find_ustar(forces, tv, u_star_2d, g, gv, us, halo=1, h_t_units=.true.)
2149
2150 if (associated(forces%frac_shelf_u)) then
2151 ! This configuration has ice shelves, and the appropriate variables need to be
2152 ! allocated. If the arrays have already been allocated, these calls do nothing.
2153 if (.not.allocated(visc%taux_shelf)) &
2154 allocate(visc%taux_shelf(g%IsdB:g%IedB, g%jsd:g%jed), source=0.0)
2155 if (.not.allocated(visc%tauy_shelf)) &
2156 allocate(visc%tauy_shelf(g%isd:g%ied, g%JsdB:g%JedB), source=0.0)
2157 if (.not.allocated(visc%tbl_thick_shelf_u)) &
2158 allocate(visc%tbl_thick_shelf_u(g%IsdB:g%IedB, g%jsd:g%jed), source=0.0)
2159 if (.not.allocated(visc%tbl_thick_shelf_v)) &
2160 allocate(visc%tbl_thick_shelf_v(g%isd:g%ied, g%JsdB:g%JedB), source=0.0)
2161 if (.not.allocated(visc%kv_tbl_shelf_u)) &
2162 allocate(visc%kv_tbl_shelf_u(g%IsdB:g%IedB, g%jsd:g%jed), source=0.0)
2163 if (.not.allocated(visc%kv_tbl_shelf_v)) &
2164 allocate(visc%kv_tbl_shelf_v(g%isd:g%ied, g%JsdB:g%JedB), source=0.0)
2165
2166 ! With a linear drag law under shelves, the friction velocity is already known.
2167! if (CS%linear_drag) ustar(:) = cdrag_sqrt_H*CS%drag_bg_vel
2168
2169 ! Find the vertical distances across layers.
2170 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
2171 endif
2172
2173 !$OMP parallel do default(shared)
2174 do j=js-1,je ; do i=is-1,ie+1
2175 mask_v(i,j) = g%mask2dCv(i,j)
2176 enddo ; enddo
2177 !$OMP parallel do default(shared)
2178 do j=js-1,je+1 ; do i=is-1,ie
2179 mask_u(i,j) = g%mask2dCu(i,j)
2180 enddo ; enddo
2181
2182 if (associated(obc)) then ; do n=1,obc%number_of_segments
2183 ! Project bottom depths across cell-corner points in the OBCs.
2184 if (.not. obc%segment(n)%on_pe) cycle
2185 ! Use a one-sided projection of bottom depths at OBC points.
2186 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
2187 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
2188 do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
2189 if (obc%segment(n)%direction == obc_direction_n) mask_u(i,j+1) = 0.0
2190 if (obc%segment(n)%direction == obc_direction_s) mask_u(i,j) = 0.0
2191 enddo
2192 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
2193 do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
2194 if (obc%segment(n)%direction == obc_direction_e) mask_v(i+1,j) = 0.0
2195 if (obc%segment(n)%direction == obc_direction_w) mask_v(i,j) = 0.0
2196 enddo
2197 endif
2198 enddo ; endif
2199
2200 !$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
2201 !$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
2202 !$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,mask_v, &
2203 !$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
2204 do j=js,je ! u-point loop
2205 if (cs%dynamic_viscous_ML) then
2206 do_any = .false.
2207 do i=isq,ieq
2208 htot(i) = 0.0
2209 if (g%mask2dCu(i,j) < 0.5) then
2210 do_i(i) = .false. ; visc%nkml_visc_u(i,j) = nkml
2211 else
2212 do_i(i) = .true. ; do_any = .true.
2213 k_massive(i) = nkml
2214 thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
2215 uhtot(i) = dt_rho0 * forces%taux(i,j)
2216 vhtot(i) = 0.25 * dt_rho0 * ((forces%tauy(i,j) + forces%tauy(i+1,j-1)) + &
2217 (forces%tauy(i,j-1) + forces%tauy(i+1,j)))
2218
2219 if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
2220 absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
2221 if (cs%omega_frac > 0.0) &
2222 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
2223 endif
2224 u_star = max(cs%ustar_min, 0.5*(u_star_2d(i,j) + u_star_2d(i+1,j)))
2225 idecay_len_tke(i) = (absf / u_star) * cs%TKE_decay
2226 endif
2227 enddo
2228
2229 if (do_any) then ; do k=1,nz
2230 if (k > nkml) then
2231 do_any = .false.
2232 if (use_eos .and. (k==nkml+1)) then
2233 ! Find dRho/dT and dRho_dS.
2234 do i=isq,ieq
2235 press(i) = (gv%H_to_RZ*gv%g_Earth) * htot(i)
2236 if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i+1,j))
2237 k2 = max(1,nkml)
2238 i_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
2239 t_eos(i) = ((h(i,j,k2)*tv%T(i,j,k2)) + (h(i+1,j,k2)*tv%T(i+1,j,k2))) * i_2hlay
2240 s_eos(i) = ((h(i,j,k2)*tv%S(i,j,k2)) + (h(i+1,j,k2)*tv%S(i+1,j,k2))) * i_2hlay
2241 enddo
2242 call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, tv%eqn_of_state, &
2243 (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
2244 if (nonbous_ml) then
2245 call calculate_specific_vol_derivs(t_eos, s_eos, press, dspv_dt, dspv_ds, tv%eqn_of_state, &
2246 (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
2247 endif
2248 endif
2249
2250 do i=isq,ieq ; if (do_i(i)) then
2251
2252 hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
2253 if (hlay > h_tiny) then ! Only consider non-vanished layers.
2254 i_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
2255 v_at_u = 0.5 * ((h(i,j,k) * (v(i,j,k) + v(i,j-1,k))) + &
2256 (h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))) * i_2hlay
2257 uh2 = (uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v_at_u)**2
2258
2259 if (use_eos) then
2260 t_lay = ((h(i,j,k)*tv%T(i,j,k)) + (h(i+1,j,k)*tv%T(i+1,j,k))) * i_2hlay
2261 s_lay = ((h(i,j,k)*tv%S(i,j,k)) + (h(i+1,j,k)*tv%S(i+1,j,k))) * i_2hlay
2262 if (nonbous_ml) then
2263 ghprime = (gv%g_Earth * gv%H_to_RZ) * (dspv_dt(i) * (thtot(i) - t_lay*htot(i)) + &
2264 dspv_ds(i) * (shtot(i) - s_lay*htot(i)))
2265 else
2266 ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
2267 dr_ds(i) * (s_lay*htot(i) - shtot(i)))
2268 endif
2269 else
2270 ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
2271 endif
2272
2273 if (ghprime > 0.0) then
2274 ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
2275 if (ribulk * uh2 <= (htot(i)**2) * ghprime) then
2276 visc%nkml_visc_u(i,j) = real(k_massive(i))
2277 do_i(i) = .false.
2278 elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
2279 visc%nkml_visc_u(i,j) = real(k-1) + &
2280 ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
2281 do_i(i) = .false.
2282 endif
2283 endif
2284 k_massive(i) = k
2285 endif ! hlay > h_tiny
2286
2287 if (do_i(i)) do_any = .true.
2288 endif ; enddo
2289
2290 if (.not.do_any) exit ! All columns are done.
2291 endif
2292
2293 do i=isq,ieq ; if (do_i(i)) then
2294 htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
2295 uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
2296 vhtot(i) = vhtot(i) + 0.25 * ((h(i,j,k) * (v(i,j,k) + v(i,j-1,k))) + &
2297 (h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))))
2298 if (use_eos) then
2299 thtot(i) = thtot(i) + 0.5 * ((h(i,j,k)*tv%T(i,j,k)) + (h(i+1,j,k)*tv%T(i+1,j,k)))
2300 shtot(i) = shtot(i) + 0.5 * ((h(i,j,k)*tv%S(i,j,k)) + (h(i+1,j,k)*tv%S(i+1,j,k)))
2301 else
2302 rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
2303 endif
2304 endif ; enddo
2305 enddo ; endif
2306
2307 if (do_any) then ; do i=isq,ieq ; if (do_i(i)) then
2308 visc%nkml_visc_u(i,j) = k_massive(i)
2309 endif ; enddo ; endif
2310 endif ! dynamic_viscous_ML
2311
2312 do_any_shelf = .false.
2313 if (associated(forces%frac_shelf_u)) then
2314 do i=isq,ieq
2315 if (forces%frac_shelf_u(i,j)*g%mask2dCu(i,j) == 0.0) then
2316 do_i(i) = .false.
2317 visc%tbl_thick_shelf_u(i,j) = 0.0 ; visc%kv_tbl_shelf_u(i,j) = 0.0
2318 else
2319 do_i(i) = .true. ; do_any_shelf = .true.
2320 endif
2321 enddo
2322 endif
2323
2324 if (do_any_shelf) then
2325 do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
2326 if (u(i,j,k) * (h(i+1,j,k) - h(i,j,k)) >= 0) then
2327 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
2328 (h(i,j,k) + h(i+1,j,k) + h_neglect)
2329 dz_at_vel(i,k) = 2.0*dz(i,j,k)*dz(i+1,j,k) / &
2330 (dz(i,j,k) + dz(i+1,j,k) + dz_neglect)
2331 else
2332 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
2333 dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
2334 endif
2335 else
2336 h_at_vel(i,k) = 0.0
2337 dz_at_vel(i,k) = 0.0
2338 ustar(i) = 0.0
2339 endif ; enddo ; enddo
2340
2341 do i=isq,ieq ; if (do_i(i)) then
2342 htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
2343 thtot(i) = 0.0 ; shtot(i) = 0.0 ; spv_htot(i) = 0.0
2344 if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
2345 if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
2346 hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
2347 if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
2348
2349 htot_vel = htot_vel + h_at_vel(i,k)
2350 hwtot = hwtot + hweight
2351
2352 if (.not.cs%linear_drag) then
2353 v_at_u = set_v_at_u(v, h, g, gv, i, j, k, mask_v, obc)
2354 ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
2355 if (cs%BBL_use_tidal_bg) then
2356 u2_bg(i) = 0.5*( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
2357 g%mask2dT(i+1,j)*(cs%tideamp(i+1,j)*cs%tideamp(i+1,j)) )
2358 else
2359 u2_bg(i) = cs%drag_bg_vel * cs%drag_bg_vel
2360 endif
2361 hutot = hutot + hweight * sqrt(u(i,j,k)**2 + v_at_u**2 + u2_bg(i))
2362 endif
2363 if (use_eos) then
2364 thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
2365 shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
2366 endif
2367 if (allocated(tv%SpV_avg)) then
2368 spv_htot(i) = spv_htot(i) + hweight * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i+1,j,k))
2369 endif
2370 enddo ; endif
2371
2372 if ((hwtot <= 0.0) .or. (cs%linear_drag .and. .not.allocated(tv%SpV_avg))) then
2373 ustar(i) = cdrag_sqrt_h * cs%drag_bg_vel
2374 elseif (cs%linear_drag .and. allocated(tv%SpV_avg)) then
2375 ustar(i) = cdrag_sqrt_h_rl * cs%drag_bg_vel * (hwtot / spv_htot(i))
2376 elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
2377 ustar(i) = cdrag_sqrt_h_rl * hutot / spv_htot(i)
2378 else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
2379 ustar(i) = cdrag_sqrt_h * hutot / hwtot
2380 endif
2381
2382 if (use_eos) then ; if (hwtot > 0.0) then
2383 t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
2384 else
2385 t_eos(i) = 0.0 ; s_eos(i) = 0.0
2386 endif ; endif
2387 ! if (allocated(tv%SpV_avg)) SpV_av(I) = SpVhtot(I) / hwtot
2388 endif ; enddo ! I-loop
2389
2390 if (use_eos) then
2391 call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), dr_dt, dr_ds, &
2392 tv%eqn_of_state, (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
2393 endif
2394
2395 do i=isq,ieq ; if (do_i(i)) then
2396 ! The 400.0 in this expression is the square of a constant proposed
2397 ! by Killworth and Edwards, 1999, in equation (2.20).
2398 ustarsq = rho0x400_g * ustar(i)**2
2399 htot(i) = 0.0 ; dztot(i) = 0.0
2400 if (use_eos) then
2401 thtot(i) = 0.0 ; shtot(i) = 0.0 ; oldfn = 0.0
2402 do k=1,nz-1
2403 if (h_at_vel(i,k) <= 0.0) cycle
2404 t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
2405 s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
2406 oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
2407 if (oldfn >= ustarsq) exit
2408
2409 dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t_lay) + &
2410 dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s_lay)) * &
2411 (h_at_vel(i,k)+htot(i))
2412 if ((oldfn + dfn) <= ustarsq) then
2413 dh = h_at_vel(i,k)
2414 ddz = dz_at_vel(i,k)
2415 else
2416 frac_used = sqrt((ustarsq-oldfn) / (dfn))
2417 dh = h_at_vel(i,k) * frac_used
2418 ddz = dz_at_vel(i,k) * frac_used
2419 endif
2420
2421 htot(i) = htot(i) + dh
2422 dztot(i) = dztot(i) + ddz
2423 thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
2424 enddo
2425 if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
2426 t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
2427 s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
2428 if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
2429 dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) then
2430 htot(i) = htot(i) + h_at_vel(i,nz)
2431 dztot(i) = dztot(i) + dz_at_vel(i,nz)
2432 endif
2433 endif ! Examination of layer nz.
2434 else ! Use Rlay as the density variable.
2435 rhtot = 0.0
2436 do k=1,nz-1
2437 rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
2438
2439 oldfn = rlay*htot(i) - rhtot(i)
2440 if (oldfn >= ustarsq) exit
2441
2442 dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
2443 if ((oldfn + dfn) <= ustarsq) then
2444 dh = h_at_vel(i,k)
2445 ddz = dz_at_vel(i,k)
2446 else
2447 frac_used = sqrt((ustarsq-oldfn) / (dfn))
2448 dh = h_at_vel(i,k) * frac_used
2449 ddz = dz_at_vel(i,k) * frac_used
2450 endif
2451
2452 htot(i) = htot(i) + dh
2453 dztot(i) = dztot(i) + ddz
2454 rhtot(i) = rhtot(i) + rlay*dh
2455 enddo
2456 if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) then
2457 htot(i) = htot(i) + h_at_vel(i,nz)
2458 dztot(i) = dztot(i) + dz_at_vel(i,nz)
2459 endif
2460 endif ! use_EOS
2461
2462 ! visc%tbl_thick_shelf_u(I,j) = max(CS%Htbl_shelf_min, &
2463 ! dztot(I) / (0.5 + sqrt(0.25 + &
2464 ! ((htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2) / &
2465 ! (ustar(i)**2) )) )
2466 ustar1 = ustar(i)
2467 h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
2468 tbl_thick = max(cs%Htbl_shelf_min, &
2469 ( dztot(i)*ustar(i) ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
2470 visc%tbl_thick_shelf_u(i,j) = tbl_thick
2471 visc%Kv_tbl_shelf_u(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar1*tbl_thick)
2472 endif ; enddo ! I-loop
2473 endif ! do_any_shelf
2474
2475 enddo ! j-loop at u-points
2476
2477 !$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
2478 !$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
2479 !$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_star_2d,mask_u, &
2480 !$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
2481 do j=jsq,jeq ! v-point loop
2482 if (cs%dynamic_viscous_ML) then
2483 do_any = .false.
2484 do i=is,ie
2485 htot(i) = 0.0
2486 if (g%mask2dCv(i,j) < 0.5) then
2487 do_i(i) = .false. ; visc%nkml_visc_v(i,j) = nkml
2488 else
2489 do_i(i) = .true. ; do_any = .true.
2490 k_massive(i) = nkml
2491 thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
2492 vhtot(i) = dt_rho0 * forces%tauy(i,j)
2493 uhtot(i) = 0.25 * dt_rho0 * ((forces%taux(i,j) + forces%taux(i-1,j+1)) + &
2494 (forces%taux(i-1,j) + forces%taux(i,j+1)))
2495
2496 if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
2497 absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
2498 if (cs%omega_frac > 0.0) &
2499 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
2500 endif
2501
2502 u_star = max(cs%ustar_min, 0.5*(u_star_2d(i,j) + u_star_2d(i,j+1)))
2503 idecay_len_tke(i) = (absf / u_star) * cs%TKE_decay
2504
2505 endif
2506 enddo
2507
2508 if (do_any) then ; do k=1,nz
2509 if (k > nkml) then
2510 do_any = .false.
2511 if (use_eos .and. (k==nkml+1)) then
2512 ! Find dRho/dT and dRho_dS.
2513 do i=is,ie
2514 press(i) = (gv%H_to_RZ * gv%g_Earth) * htot(i)
2515 if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i,j+1))
2516 k2 = max(1,nkml)
2517 i_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
2518 t_eos(i) = ((h(i,j,k2)*tv%T(i,j,k2)) + (h(i,j+1,k2)*tv%T(i,j+1,k2))) * i_2hlay
2519 s_eos(i) = ((h(i,j,k2)*tv%S(i,j,k2)) + (h(i,j+1,k2)*tv%S(i,j+1,k2))) * i_2hlay
2520 enddo
2521 call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
2522 tv%eqn_of_state, (/is-g%IsdB+1,ie-g%IsdB+1/) )
2523 if (nonbous_ml) then
2524 call calculate_specific_vol_derivs(t_eos, s_eos, press, dspv_dt, dspv_ds, tv%eqn_of_state, &
2525 (/is-g%IsdB+1,ie-g%IsdB+1/) )
2526 endif
2527 endif
2528
2529 do i=is,ie ; if (do_i(i)) then
2530
2531 hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
2532 if (hlay > h_tiny) then ! Only consider non-vanished layers.
2533 i_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
2534 u_at_v = 0.5 * ((h(i,j,k) * (u(i-1,j,k) + u(i,j,k))) + &
2535 (h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))) * i_2hlay
2536 uh2 = (vhtot(i) - htot(i)*v(i,j,k))**2 + (uhtot(i) - htot(i)*u_at_v)**2
2537
2538 if (use_eos) then
2539 t_lay = ((h(i,j,k)*tv%T(i,j,k)) + (h(i,j+1,k)*tv%T(i,j+1,k))) * i_2hlay
2540 s_lay = ((h(i,j,k)*tv%S(i,j,k)) + (h(i,j+1,k)*tv%S(i,j+1,k))) * i_2hlay
2541 if (nonbous_ml) then
2542 ghprime = (gv%g_Earth * gv%H_to_RZ) * (dspv_dt(i) * (thtot(i) - t_lay*htot(i)) + &
2543 dspv_ds(i) * (shtot(i) - s_lay*htot(i)))
2544 else
2545 ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
2546 dr_ds(i) * (s_lay*htot(i) - shtot(i)))
2547 endif
2548 else
2549 ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
2550 endif
2551
2552 if (ghprime > 0.0) then
2553 ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
2554 if (ribulk * uh2 <= htot(i)**2 * ghprime) then
2555 visc%nkml_visc_v(i,j) = real(k_massive(i))
2556 do_i(i) = .false.
2557 elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
2558 visc%nkml_visc_v(i,j) = real(k-1) + &
2559 ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
2560 do_i(i) = .false.
2561 endif
2562 endif
2563 k_massive(i) = k
2564 endif ! hlay > h_tiny
2565
2566 if (do_i(i)) do_any = .true.
2567 endif ; enddo
2568
2569 if (.not.do_any) exit ! All columns are done.
2570 endif
2571
2572 do i=is,ie ; if (do_i(i)) then
2573 htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
2574 vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
2575 uhtot(i) = uhtot(i) + 0.25 * ((h(i,j,k) * (u(i-1,j,k) + u(i,j,k))) + &
2576 (h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))))
2577 if (use_eos) then
2578 thtot(i) = thtot(i) + 0.5 * ((h(i,j,k)*tv%T(i,j,k)) + (h(i,j+1,k)*tv%T(i,j+1,k)))
2579 shtot(i) = shtot(i) + 0.5 * ((h(i,j,k)*tv%S(i,j,k)) + (h(i,j+1,k)*tv%S(i,j+1,k)))
2580 else
2581 rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
2582 endif
2583 endif ; enddo
2584 enddo ; endif
2585
2586 if (do_any) then ; do i=is,ie ; if (do_i(i)) then
2587 visc%nkml_visc_v(i,j) = k_massive(i)
2588 endif ; enddo ; endif
2589 endif ! dynamic_viscous_ML
2590
2591 do_any_shelf = .false.
2592 if (associated(forces%frac_shelf_v)) then
2593 do i=is,ie
2594 if (forces%frac_shelf_v(i,j)*g%mask2dCv(i,j) == 0.0) then
2595 do_i(i) = .false.
2596 visc%tbl_thick_shelf_v(i,j) = 0.0 ; visc%kv_tbl_shelf_v(i,j) = 0.0
2597 else
2598 do_i(i) = .true. ; do_any_shelf = .true.
2599 endif
2600 enddo
2601 endif
2602
2603 if (do_any_shelf) then
2604 do k=1,nz ; do i=is,ie ; if (do_i(i)) then
2605 if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
2606 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
2607 (h(i,j,k) + h(i,j+1,k) + h_neglect)
2608 dz_at_vel(i,k) = 2.0*dz(i,j,k)*dz(i,j+1,k) / &
2609 (dz(i,j,k) + dz(i,j+1,k) + dz_neglect)
2610 else
2611 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
2612 dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
2613 endif
2614 else
2615 h_at_vel(i,k) = 0.0
2616 dz_at_vel(i,k) = 0.0
2617 ustar(i) = 0.0
2618 endif ; enddo ; enddo
2619
2620 do i=is,ie ; if (do_i(i)) then
2621 htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
2622 thtot(i) = 0.0 ; shtot(i) = 0.0 ; spv_htot(i) = 0.0
2623 if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
2624 if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
2625 hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
2626 if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
2627
2628 htot_vel = htot_vel + h_at_vel(i,k)
2629 hwtot = hwtot + hweight
2630
2631 if (.not.cs%linear_drag) then
2632 u_at_v = set_u_at_v(u, h, g, gv, i, j, k, mask_u, obc)
2633 ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
2634 if (cs%BBL_use_tidal_bg) then
2635 u2_bg(i) = 0.5*( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
2636 g%mask2dT(i,j+1)*(cs%tideamp(i,j+1)*cs%tideamp(i,j+1)) )
2637 else
2638 u2_bg(i) = cs%drag_bg_vel * cs%drag_bg_vel
2639 endif
2640 hutot = hutot + hweight * sqrt(v(i,j,k)**2 + u_at_v**2 + u2_bg(i))
2641 endif
2642 if (use_eos) then
2643 thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
2644 shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
2645 endif
2646 if (allocated(tv%SpV_avg)) then
2647 spv_htot(i) = spv_htot(i) + hweight * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j+1,k))
2648 endif
2649 enddo ; endif
2650
2651 if ((hwtot <= 0.0) .or. (cs%linear_drag .and. .not.allocated(tv%SpV_avg))) then
2652 ustar(i) = cdrag_sqrt_h * cs%drag_bg_vel
2653 elseif (cs%linear_drag .and. allocated(tv%SpV_avg)) then
2654 ustar(i) = cdrag_sqrt_h_rl * cs%drag_bg_vel * (hwtot / spv_htot(i))
2655 elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
2656 ustar(i) = cdrag_sqrt_h_rl * hutot / spv_htot(i)
2657 else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
2658 ustar(i) = cdrag_sqrt_h * hutot / hwtot
2659 endif
2660
2661 if (use_eos) then ; if (hwtot > 0.0) then
2662 t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
2663 else
2664 t_eos(i) = 0.0 ; s_eos(i) = 0.0
2665 endif ; endif
2666 endif ; enddo ! I-loop
2667
2668 if (use_eos) then
2669 call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), dr_dt, dr_ds, &
2670 tv%eqn_of_state, (/is-g%IsdB+1,ie-g%IsdB+1/) )
2671 endif
2672
2673 do i=is,ie ; if (do_i(i)) then
2674 ! The 400.0 in this expression is the square of a constant proposed
2675 ! by Killworth and Edwards, 1999, in equation (2.20).
2676 ustarsq = rho0x400_g * ustar(i)**2
2677 htot(i) = 0.0
2678 dztot(i) = 0.0
2679 if (use_eos) then
2680 thtot(i) = 0.0 ; shtot(i) = 0.0 ; oldfn = 0.0
2681 do k=1,nz-1
2682 if (h_at_vel(i,k) <= 0.0) cycle
2683 t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
2684 s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
2685 oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
2686 if (oldfn >= ustarsq) exit
2687
2688 dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t_lay) + &
2689 dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s_lay)) * &
2690 (h_at_vel(i,k)+htot(i))
2691 if ((oldfn + dfn) <= ustarsq) then
2692 dh = h_at_vel(i,k)
2693 ddz = dz_at_vel(i,k)
2694 else
2695 frac_used = sqrt((ustarsq-oldfn) / (dfn))
2696 dh = h_at_vel(i,k) * frac_used
2697 ddz = dz_at_vel(i,k) * frac_used
2698 endif
2699
2700 htot(i) = htot(i) + dh
2701 dztot(i) = dztot(i) + ddz
2702 thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
2703 enddo
2704 if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
2705 t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
2706 s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
2707 if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
2708 dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) then
2709 htot(i) = htot(i) + h_at_vel(i,nz)
2710 dztot(i) = dztot(i) + dz_at_vel(i,nz)
2711 endif
2712 endif ! Examination of layer nz.
2713 else ! Use Rlay as the density variable.
2714 rhtot = 0.0
2715 do k=1,nz-1
2716 rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
2717
2718 oldfn = rlay*htot(i) - rhtot(i)
2719 if (oldfn >= ustarsq) exit
2720
2721 dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
2722 if ((oldfn + dfn) <= ustarsq) then
2723 dh = h_at_vel(i,k)
2724 ddz = dz_at_vel(i,k)
2725 else
2726 frac_used = sqrt((ustarsq-oldfn) / (dfn))
2727 dh = h_at_vel(i,k) * frac_used
2728 ddz = dz_at_vel(i,k) * frac_used
2729 endif
2730
2731 htot(i) = htot(i) + dh
2732 dztot(i) = dztot(i) + ddz
2733 rhtot = rhtot + rlay*dh
2734 enddo
2735 if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) then
2736 htot(i) = htot(i) + h_at_vel(i,nz)
2737 dztot(i) = dztot(i) + dz_at_vel(i,nz)
2738 endif
2739 endif ! use_EOS
2740
2741 ! visc%tbl_thick_shelf_v(i,J) = max(CS%Htbl_shelf_min, &
2742 ! dztot(i) / (0.5 + sqrt(0.25 + &
2743 ! (htot(i)*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))**2 / &
2744 ! (ustar(i))**2 )) )
2745 ustar1 = ustar(i)
2746 h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
2747 tbl_thick = max(cs%Htbl_shelf_min, &
2748 ( dztot(i)*ustar(i) ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
2749 visc%tbl_thick_shelf_v(i,j) = tbl_thick
2750 visc%Kv_tbl_shelf_v(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar1*tbl_thick)
2751
2752 endif ; enddo ! i-loop
2753 endif ! do_any_shelf
2754
2755 enddo ! J-loop at v-points
2756
2757 if (cs%debug) then
2758 if (allocated(visc%nkml_visc_u) .and. allocated(visc%nkml_visc_v)) &
2759 call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, visc%nkml_visc_v, &
2760 g%HI, haloshift=0, scalar_pair=.true.)
2761 endif
2762 if (cs%id_nkml_visc_u > 0) call post_data(cs%id_nkml_visc_u, visc%nkml_visc_u, cs%diag)
2763 if (cs%id_nkml_visc_v > 0) call post_data(cs%id_nkml_visc_v, visc%nkml_visc_v, cs%diag)
2764
2765end subroutine set_viscous_ml
2766
2767!> Register any fields associated with the vertvisc_type.
2768subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_CS, use_ice_shelf)
2769 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
2770 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2771 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2772 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2773 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2774 !! parameters.
2775 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
2776 !! viscosities and related fields.
2777 !! Allocated here.
2778 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
2779 logical, intent(in) :: use_ice_shelf !< if true, register tau_shelf restarts
2780 ! Local variables
2781 logical :: use_kappa_shear, ks_at_vertex
2782 logical :: adiabatic, usekpp, useepbl, use_ideal_age
2783 logical :: do_brine_plume, use_hor_bnd_diff, use_neutral_diffusion, use_fpmix
2784 logical :: use_cvmix_shear, mle_use_pbl_mld, mle_use_bodner, use_cvmix_conv
2785 integer :: isd, ied, jsd, jed, nz
2786 real :: hfreeze !< If hfreeze > 0 [Z ~> m], melt potential will be computed.
2787 character(len=16) :: kv_units, kd_units
2788 character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
2789 type(vardesc) :: u_desc, v_desc
2790 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
2791
2792 call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
2793 do_not_log=.true.)
2794
2795 use_kappa_shear = .false. ; ks_at_vertex = .false. ; use_cvmix_shear = .false.
2796 usekpp = .false. ; useepbl = .false. ; use_cvmix_conv = .false.
2797
2798 if (.not.adiabatic) then
2799 use_kappa_shear = kappa_shear_is_used(param_file)
2800 ks_at_vertex = kappa_shear_at_vertex(param_file)
2801 use_cvmix_shear = cvmix_shear_is_used(param_file)
2802 use_cvmix_conv = cvmix_conv_is_used(param_file)
2803 call get_param(param_file, mdl, "USE_KPP", usekpp, &
2804 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
2805 "to calculate diffusivities and non-local transport in the OBL.", &
2806 default=.false., do_not_log=.true.)
2807 call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, &
2808 "If true, use an implied energetics planetary boundary "//&
2809 "layer scheme to determine the diffusivity and viscosity "//&
2810 "in the surface boundary layer.", default=.false., do_not_log=.true.)
2811 endif
2812
2813 if (gv%Boussinesq) then
2814 kv_units = "m2 s-1" ; kd_units = "m2 s-1"
2815 else
2816 kv_units = "Pa s" ; kd_units = "kg m-1 s-1"
2817 endif
2818
2819 if (use_kappa_shear .or. usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv) then
2820 call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1)
2821 call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_cs, &
2822 "Shear-driven turbulent diffusivity at interfaces", &
2823 units=kd_units, conversion=gv%HZ_T_to_MKS, z_grid='i')
2824 endif
2825 if (usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv .or. &
2826 (use_kappa_shear .and. .not.ks_at_vertex )) then
2827 call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1)
2828 call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_cs, &
2829 "Shear-driven turbulent viscosity at interfaces", &
2830 units=kv_units, conversion=gv%HZ_T_to_MKS, z_grid='i')
2831 endif
2832 if (use_kappa_shear .and. ks_at_vertex) then
2833 call safe_alloc_ptr(visc%TKE_turb, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
2834 call safe_alloc_ptr(visc%Kv_shear_Bu, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
2835 call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_cs, &
2836 "Shear-driven turbulent viscosity at vertex interfaces", &
2837 units=kv_units, conversion=gv%HZ_T_to_MKS, hor_grid="Bu", z_grid='i')
2838 elseif (use_kappa_shear) then
2839 call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1)
2840 endif
2841
2842 if (usekpp) then
2843 ! MOM_bkgnd_mixing uses Kv_slow when KPP is defined.
2844 call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1)
2845 endif
2846
2847 ! visc%MLD and visc%h_ML are used to communicate the state of the (e)PBL or KPP to the rest of the model
2848 call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
2849 default=.false., do_not_log=.true.)
2850 ! visc%h_ML needs to be allocated when melt potential is computed (HFREEZE>0) or one of
2851 ! several other parameterizations are in use.
2852 call get_param(param_file, mdl, "HFREEZE", hfreeze, &
2853 units="m", default=-1.0, scale=us%m_to_Z, do_not_log=.true.)
2854 call get_param(param_file, mdl, "DO_BRINE_PLUME", do_brine_plume, &
2855 "If true, use a brine plume parameterization from Nguyen et al., 2009.", &
2856 default=.false., do_not_log=.true.)
2857 call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", use_hor_bnd_diff, &
2858 default=.false., do_not_log=.true.)
2859 call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", use_neutral_diffusion, &
2860 default=.false., do_not_log=.true.)
2861 if (use_neutral_diffusion) &
2862 call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", use_neutral_diffusion, &
2863 default=.false., do_not_log=.true.)
2864 call get_param(param_file, mdl, "FPMIX", use_fpmix, &
2865 default=.false., do_not_log=.true.)
2866 call get_param(param_file, mdl, "USE_IDEAL_AGE_TRACER", use_ideal_age, &
2867 default=.false., do_not_log=.true.)
2868 call openparameterblock(param_file, 'MLE', do_not_log=.true.)
2869 call get_param(param_file, mdl, "USE_BODNER23", mle_use_bodner, &
2870 default=.false., do_not_log=.true.)
2871 call closeparameterblock(param_file)
2872
2873 if (mle_use_pbl_mld .or. mle_use_bodner) then
2874 call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
2875 endif
2876 if ((hfreeze >= 0.0) .or. mle_use_pbl_mld .or. use_fpmix .or. &
2877 use_neutral_diffusion .or. use_hor_bnd_diff .or. use_ideal_age) then
2878 call safe_alloc_ptr(visc%h_ML, isd, ied, jsd, jed)
2879 endif
2880 if (do_brine_plume) then
2881 call safe_alloc_ptr(visc%h_ML_param, isd, ied, jsd, jed)
2882 call safe_alloc_ptr(visc%MLD_param, isd, ied, jsd, jed)
2883 endif
2884
2885 if (mle_use_pbl_mld) then
2886 call register_restart_field(visc%MLD, "MLD", .false., restart_cs, &
2887 "Instantaneous active mixing layer depth", units="m", conversion=us%Z_to_m)
2888 endif
2889 if (mle_use_pbl_mld .or. use_fpmix .or. use_neutral_diffusion .or. use_hor_bnd_diff) then
2890 call register_restart_field(visc%h_ML, "h_ML", .false., restart_cs, &
2891 "Instantaneous active mixing layer thickness", &
2892 units=get_thickness_units(gv), conversion=gv%H_to_mks)
2893 endif
2894 if (do_brine_plume) then
2895 call register_restart_field(visc%h_ML_param, "h_ML_param", .false., restart_cs, &
2896 "Instantaneous active mixed layer thickness", &
2897 units=get_thickness_units(gv), conversion=gv%H_to_mks)
2898 endif
2899
2900 ! visc%sfc_buoy_flx is used to communicate the state of the (e)PBL or KPP to the rest of the model
2901 if (mle_use_pbl_mld .or. mle_use_bodner) then
2902 call safe_alloc_ptr(visc%sfc_buoy_flx, isd, ied, jsd, jed)
2903 call register_restart_field(visc%sfc_buoy_flx, "SFC_BFLX", .false., restart_cs, &
2904 "Instantaneous surface buoyancy flux", "m2 s-3", &
2905 conversion=us%Z_to_m**2*us%s_to_T**3)
2906 endif
2907
2908 if (use_ice_shelf) then
2909 if (.not.allocated(visc%taux_shelf)) &
2910 allocate(visc%taux_shelf(g%IsdB:g%IedB, g%jsd:g%jed), source=0.0)
2911 if (.not.allocated(visc%tauy_shelf)) &
2912 allocate(visc%tauy_shelf(g%isd:g%ied, g%JsdB:g%JedB), source=0.0)
2913 u_desc = var_desc("u_taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", &
2914 hor_grid='Cu',z_grid='1')
2915 v_desc = var_desc("v_tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", &
2916 hor_grid='Cv',z_grid='1')
2917 call register_restart_pair(visc%taux_shelf, visc%tauy_shelf, u_desc, v_desc, &
2918 .false., restart_cs, conversion=us%RZ_T_to_kg_m2s*us%L_T_to_m_s)
2919 endif
2920
2921end subroutine set_visc_register_restarts
2922
2923!> This subroutine does remapping for the auxiliary restart variables in a vertvisc_type
2924!! that are used across timesteps
2925subroutine remap_vertvisc_aux_vars(G, GV, visc, h_old, h_new, ALE_CSp, OBC)
2926 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
2927 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
2928 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
2929 !! viscosities and related fields.
2930 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2931 intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2]
2932 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2933 intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2]
2934 type(ale_cs), pointer :: ale_csp !< ALE control structure to use when remapping
2935 type(ocean_obc_type), pointer :: obc !< Open boundary structure
2936
2937 if (associated(visc%Kd_shear)) then
2938 call ale_remap_interface_vals(ale_csp, g, gv, h_old, h_new, visc%Kd_shear)
2939 endif
2940
2941 if (associated(visc%Kv_shear)) then
2942 call ale_remap_interface_vals(ale_csp, g, gv, h_old, h_new, visc%Kv_shear)
2943 endif
2944
2945 if (associated(visc%Kv_shear_Bu)) then
2946 call ale_remap_vertex_vals(ale_csp, g, gv, h_old, h_new, visc%Kv_shear_Bu)
2947 endif
2948
2949end subroutine remap_vertvisc_aux_vars
2950
2951!> Initializes the MOM_set_visc control structure
2952subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
2953 type(time_type), target, intent(in) :: time !< The current model time.
2954 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
2955 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2956 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2957 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2958 !! parameters.
2959 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
2960 !! output.
2961 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
2962 !! related fields.
2963 type(set_visc_cs), intent(inout) :: cs !< Vertical viscosity control structure
2964 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
2965 type(ocean_obc_type), pointer :: obc !< A pointer to an open boundary condition structure
2966
2967 ! Local variables
2968 real :: csmag_chan_dflt ! The default value for SMAG_CONST_CHANNEL [nondim]
2969 real :: smag_const1 ! The default value for the Smagorinsky Laplacian coefficient [nondim]
2970 real :: tke_decay_dflt ! The default value of a coefficient scaling the vertical decay
2971 ! rate of TKE [nondim]
2972 real :: bulk_ri_ml_dflt ! The default bulk Richardson number for a bulk mixed layer [nondim]
2973 real :: kv_background ! The background kinematic viscosity in the interior [Z2 T-1 ~> m2 s-1]
2974 real :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate that
2975 ! is used in place of the absolute value of the local Coriolis
2976 ! parameter in the denominator of some expressions [nondim]
2977 real :: chan_max_thick_dflt ! The default value for CHANNEL_DRAG_MAX_THICK [Z ~> m]
2978 real :: tideamp_factor ! A factor to multiply by tideamp when converting to mean tidal magnitude [nondim]
2979 real :: shelfbreak_depth ! When CHANNEL_DRAG is true, the bathymetric depth interpolated
2980 ! to the vorticity point is a combination of the harmonic mean of the
2981 ! adjacent velocity point depths below this depth [Z ~> m] and the
2982 ! arithmetic mean of the adjacent depths above it, to roughly mimic a
2983 ! continental shelf break profile.
2984 real, allocatable, dimension(:,:) :: cdrag_h !< The spatially varying quadratic drag coefficient [nondim]
2985
2986 integer :: i, j, is, ie, js, je
2987 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
2988 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
2989 logical :: adiabatic, use_omega, mle_use_pbl_mld
2990 logical :: use_kpp
2991 logical :: use_regridding ! If true, use the ALE algorithm rather than layered
2992 ! isopycnal or stacked shallow water mode.
2993 logical :: use_temperature ! If true, temperature and salinity are used as state variables.
2994 logical :: use_eos ! If true, density calculated from T & S using an equation of state.
2995 character(len=200) :: filename, cdrag_file, tideamp_file ! Input file names or paths
2996 character(len=80) :: cdrag_var, tideamp_var ! Input file variable names
2997 ! This include declares and sets the variable "version".
2998# include "version_variable.h"
2999 character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
3000
3001 cs%initialized = .true.
3002 cs%OBC => obc
3003
3004 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3005 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
3006 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3007
3008 cs%diag => diag
3009
3010 ! Set default, read and log parameters
3011 call log_version(param_file, mdl, version, "")
3012 cs%RiNo_mix = .false.
3013 call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
3014 cs%inputdir = slasher(cs%inputdir)
3015 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
3016 "This sets the default value for the various _ANSWER_DATE parameters.", &
3017 default=99991231)
3018 call get_param(param_file, mdl, "SET_VISC_ANSWER_DATE", cs%answer_date, &
3019 "The vintage of the order of arithmetic and expressions " // &
3020 "in the set viscosity calculations. Values below " // &
3021 "20190101 recover the answers from the end of 2018, " // &
3022 "while higher values use updated and more robust forms " // &
3023 "of the same expressions. Values below 20260704 use " // &
3024 "the non-reproducible power operator in place of " // &
3025 "cuberoot() when finding the open distances with " // &
3026 "CHANNEL_DRAG.", &
3027 default=min(20260703,default_answer_date), &
3028 do_not_log=.not.gv%Boussinesq)
3029 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)
3030 call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
3031 "If true, the bottom stress is calculated with a drag "//&
3032 "law of the form c_drag*|u|*u. The velocity magnitude "//&
3033 "may be an assumed value or it may be based on the "//&
3034 "actual velocity in the bottommost HBBL, depending on "//&
3035 "LINEAR_DRAG.", default=.true.)
3036 call get_param(param_file, mdl, "DRAG_AS_BODY_FORCE", cs%body_force_drag, &
3037 "If true, the bottom stress is imposed as an explicit body force "//&
3038 "applied over a fixed distance from the bottom, rather than as an "//&
3039 "implicit calculation based on an enhanced near-bottom viscosity. "//&
3040 "The thickness of the bottom boundary layer is HBBL.", &
3041 default=.false., do_not_log=.not.cs%bottomdraglaw)
3042 call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
3043 "If true, the bottom drag is exerted directly on each "//&
3044 "layer proportional to the fraction of the bottom it overlies.", &
3045 default=.false.)
3046 call get_param(param_file, mdl, "CHANNEL_DRAG_SHELFBREAK_DEPTH", shelfbreak_depth, &
3047 "When CHANNEL_DRAG is true, the bathymetric depth interpolated to the "//&
3048 "vorticity point is a combination of the harmonic mean of the adjacent "//&
3049 "velocity point depths below this depth and the arithmetic mean of the "//&
3050 "depths above it, to roughly mimic a continental shelf break profile. "//&
3051 "Setting this to exceed MAXIMUM_DEPTH leads to linear interpolation of "//&
3052 "the topography between velocity points.", &
3053 default=0.0, units="m", scale=us%m_to_Z, do_not_log=.not.cs%Channel_drag)
3054 cs%channel_break_depth = shelfbreak_depth - g%Z_ref
3055
3056 call get_param(param_file, mdl, "LINEAR_DRAG", cs%linear_drag, &
3057 "If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
3058 "law is cdrag*DRAG_BG_VEL*u.", default=.false.)
3059 call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
3060 do_not_log=.true.)
3061 if (adiabatic) then
3062 call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
3063 "There are no diapycnal mass fluxes if ADIABATIC is true. "//&
3064 "This assumes that KD = 0.0 and that there is no buoyancy forcing, "//&
3065 "but makes the model faster by eliminating subroutine calls.", default=.false.)
3066 endif
3067
3068 if (.not.adiabatic) then
3069 cs%RiNo_mix = kappa_shear_is_used(param_file)
3070 endif
3071
3072 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
3073
3074 call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
3075 "If true, use a bulk Richardson number criterion to "//&
3076 "determine the mixed layer thickness for viscosity.", &
3077 default=.false.)
3078 if (cs%dynamic_viscous_ML) then
3079 call get_param(param_file, mdl, "BULK_RI_ML", bulk_ri_ml_dflt, units="nondim", default=0.0)
3080 call get_param(param_file, mdl, "BULK_RI_ML_VISC", cs%bulk_Ri_ML, &
3081 "The efficiency with which mean kinetic energy released by mechanically "//&
3082 "forced entrainment of the mixed layer is converted to turbulent "//&
3083 "kinetic energy. By default, BULK_RI_ML_VISC = BULK_RI_ML or 0.", &
3084 units="nondim", default=bulk_ri_ml_dflt)
3085 call get_param(param_file, mdl, "TKE_DECAY", tke_decay_dflt, units="nondim", default=0.0)
3086 call get_param(param_file, mdl, "TKE_DECAY_VISC", cs%TKE_decay, &
3087 "TKE_DECAY_VISC relates the vertical rate of decay of "//&
3088 "the TKE available for mechanical entrainment to the "//&
3089 "natural Ekman depth for use in calculating the dynamic "//&
3090 "mixed layer viscosity. By default, TKE_DECAY_VISC = TKE_DECAY or 0.", &
3091 units="nondim", default=tke_decay_dflt)
3092 call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
3093 "If true, use the absolute rotation rate instead of the "//&
3094 "vertical component of rotation when setting the decay "//&
3095 "scale for turbulence.", default=.false., do_not_log=.true.)
3096 omega_frac_dflt = 0.0
3097 if (use_omega) then
3098 call mom_error(warning, "ML_USE_OMEGA is deprecated; use ML_OMEGA_FRAC=1.0 instead.")
3099 omega_frac_dflt = 1.0
3100 endif
3101 call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
3102 "When setting the decay scale for turbulence, use this "//&
3103 "fraction of the absolute rotation rate blended with the "//&
3104 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
3105 units="nondim", default=omega_frac_dflt)
3106 call get_param(param_file, mdl, "OMEGA", cs%omega, &
3107 "The rotation rate of the earth.", &
3108 units="s-1", default=7.2921e-5, scale=us%T_to_s)
3109 ! This give a minimum decay scale that is typically much less than Angstrom.
3110 cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_H + gv%H_subroundoff)
3111 else
3112 call get_param(param_file, mdl, "OMEGA", cs%omega, &
3113 "The rotation rate of the earth.", &
3114 units="s-1", default=7.2921e-5, scale=us%T_to_s)
3115 endif
3116
3117 call get_param(param_file, mdl, "HBBL", cs%dz_bbl, &
3118 "The thickness of a bottom boundary layer with a viscosity increased by "//&
3119 "KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//&
3120 "near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//&
3121 "defined but LINEAR_DRAG is not.", &
3122 units="m", scale=us%m_to_Z, fail_if_missing=.true.) ! Rescaled later
3123 if (cs%bottomdraglaw) then
3124 call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
3125 "CDRAG is the drag coefficient relating the magnitude of "//&
3126 "the velocity field to the bottom stress. CDRAG is only "//&
3127 "used if BOTTOMDRAGLAW is defined.", units="nondim", default=0.003)
3128 call get_param(param_file, mdl, "CDRAG_MAP", cs%bottomdragmap, &
3129 "If true, apply a spatially varying scaling factor to CDRAG, "//&
3130 "specified by CDRAG_VAR in CDRAG_FILE.", default=.false.)
3131 call get_param(param_file, mdl, "CDRAG_FILE", cdrag_file, &
3132 "The name of the file with the spatially varying bottom drag "//&
3133 "scaling factor.", default="", do_not_log=.not.cs%bottomdragmap)
3134 call get_param(param_file, mdl, "CDRAG_VAR", cdrag_var, &
3135 "The name of the variable in CDRAG_FILE with the spatially "//&
3136 "varying bottom drag scaling factor at h points.", &
3137 default="", do_not_log=.not.cs%bottomdragmap)
3138 call get_param(param_file, mdl, "BBL_USE_TIDAL_BG", cs%BBL_use_tidal_bg, &
3139 "Flag to use the tidal RMS amplitude in place of constant "//&
3140 "background velocity for computing u* in the BBL. "//&
3141 "This flag is only used when BOTTOMDRAGLAW is true and "//&
3142 "LINEAR_DRAG is false.", default=.false.)
3143 if (cs%BBL_use_tidal_bg) then
3144 call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
3145 "The path to the file containing the spatially varying "//&
3146 "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
3147 call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
3148 "The name of the tidal amplitude variable in the input file.", &
3149 default="tideamp")
3150 ! This value is here only to detect whether it is inadvertently used. CS%drag_bg_vel should
3151 ! not be used if CS%BBL_use_tidal_bg is True. For this reason, we do not apply dimensions,
3152 ! nor dimensional testing in this mode. If we ever detect a dimensional sensitivity to
3153 ! this parameter, in this mode, then it means it is being used inappropriately.
3154 cs%drag_bg_vel = 1.e30
3155 call get_param(param_file, mdl, "TIDEAMP_FACTOR", tideamp_factor, &
3156 "A parameter to multiply by tideamp when converting to ustar. "//&
3157 "It accounts for converting the amplitude to a mean magintude (approx 1/sqrt(2)) "//&
3158 "and possibly also for non-commuting averaging operators when converting to ustar**3. "//&
3159 "It is ignored if negative and uncapped so it can be greater than 1 if desired.",&
3160 units="nondim", default=-1.0)
3161 if (tideamp_factor < 0.0) then
3162 cs%tideampfac2 = 1.0
3163 else
3164 cs%tideampfac2 = tideamp_factor*tideamp_factor
3165 endif
3166 else
3167 call get_param(param_file, mdl, "DRAG_BG_VEL", cs%drag_bg_vel, &
3168 "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
3169 "LINEAR_DRAG) or an unresolved velocity that is "//&
3170 "combined with the resolved velocity to estimate the "//&
3171 "velocity magnitude. DRAG_BG_VEL is only used when "//&
3172 "BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
3173 endif
3174 call get_param(param_file, mdl, "USE_REGRIDDING", use_regridding, &
3175 do_not_log=.true., default=.false. )
3176 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3177 default=.true., do_not_log=.true.)
3178 call get_param(param_file, mdl, "USE_EOS", use_eos, &
3179 default=use_temperature, do_not_log=.true.)
3180 call get_param(param_file, mdl, "BBL_USE_EOS", cs%BBL_use_EOS, &
3181 "If true, use the equation of state in determining the properties of the "//&
3182 "bottom boundary layer. Otherwise use the layer target potential densities. "//&
3183 "The default of this parameter is the value of USE_EOS.", &
3184 default=use_eos, do_not_log=.not.use_temperature)
3185 if (use_regridding .and. (.not. cs%BBL_use_EOS)) &
3186 call mom_error(fatal,"When using MOM6 in ALE mode it is required to set BBL_USE_EOS to True.")
3187 endif
3188 call get_param(param_file, mdl, "BBL_THICK_MIN", cs%BBL_thick_min, &
3189 "The minimum bottom boundary layer thickness that can be "//&
3190 "used with BOTTOMDRAGLAW. This might be "//&
3191 "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
3192 "near-bottom viscosity.", units="m", default=0.0, scale=us%m_to_Z)
3193 call get_param(param_file, mdl, "HTBL_SHELF_MIN", cs%Htbl_shelf_min, &
3194 "The minimum top boundary layer thickness that can be "//&
3195 "used with BOTTOMDRAGLAW. This might be "//&
3196 "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
3197 "near-top viscosity.", units="m", default=us%Z_to_m*cs%BBL_thick_min, scale=us%m_to_Z)
3198 call get_param(param_file, mdl, "HTBL_SHELF", cs%Htbl_shelf, &
3199 "The thickness over which near-surface velocities are "//&
3200 "averaged for the drag law under an ice shelf. By "//&
3201 "default this is the same as HBBL", &
3202 units="m", default=us%Z_to_m*cs%dz_bbl, scale=gv%m_to_H)
3203
3204 call get_param(param_file, mdl, "KV", kv_background, &
3205 "The background kinematic viscosity in the interior. "//&
3206 "The molecular value, ~1e-6 m2 s-1, may be used.", &
3207 units="m2 s-1", scale=us%m2_s_to_Z2_T, fail_if_missing=.true.)
3208
3209 call get_param(param_file, mdl, "USE_KPP", use_kpp, &
3210 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3211 "to calculate diffusivities and non-local transport in the OBL.", &
3212 do_not_log=.true., default=.false.)
3213
3214 call get_param(param_file, mdl, "KV_BBL_MIN", cs%KV_BBL_min, &
3215 "The minimum viscosities in the bottom boundary layer.", &
3216 units="m2 s-1", default=us%Z2_T_to_m2_s*kv_background, scale=gv%m2_s_to_HZ_T)
3217 call get_param(param_file, mdl, "KV_TBL_MIN", cs%KV_TBL_min, &
3218 "The minimum viscosities in the top boundary layer.", &
3219 units="m2 s-1", default=us%Z2_T_to_m2_s*kv_background, scale=gv%m2_s_to_HZ_T)
3220 call get_param(param_file, mdl, "CORRECT_BBL_BOUNDS", cs%correct_BBL_bounds, &
3221 "If true, uses the correct bounds on the BBL thickness and "//&
3222 "viscosity so that the bottom layer feels the intended drag.", &
3223 default=.false.)
3224
3225 if (cs%Channel_drag) then
3226 call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, units="nondim", default=-1.0)
3227
3228 csmag_chan_dflt = 0.15
3229 if (smag_const1 >= 0.0) csmag_chan_dflt = smag_const1
3230
3231 call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", cs%c_Smag, &
3232 "The nondimensional Laplacian Smagorinsky constant used "//&
3233 "in calculating the channel drag if it is enabled. The "//&
3234 "default is to use the same value as SMAG_LAP_CONST if "//&
3235 "it is defined, or 0.15 if it is not. The value used is "//&
3236 "also 0.15 if the specified value is negative.", &
3237 units="nondim", default=csmag_chan_dflt, do_not_log=.not.cs%Channel_drag)
3238 if (cs%c_Smag < 0.0) cs%c_Smag = 0.15
3239
3240 call get_param(param_file, mdl, "TRIG_CHANNEL_DRAG_WIDTHS", cs%concave_trigonometric_L, &
3241 "If true, use trigonometric expressions to determine the fractional open "//&
3242 "interface lengths for concave topography.", &
3243 default=.true., do_not_log=.not.cs%Channel_drag)
3244 endif
3245
3246 chan_max_thick_dflt = -1.0*us%m_to_Z
3247 if (cs%RiNo_mix) chan_max_thick_dflt = 0.5*cs%dz_bbl
3248 if (cs%body_force_drag) chan_max_thick_dflt = cs%dz_bbl
3249 call get_param(param_file, mdl, "CHANNEL_DRAG_MAX_BBL_THICK", cs%Chan_drag_max_vol, &
3250 "The maximum bottom boundary layer thickness over which the channel drag is "//&
3251 "exerted, or a negative value for no fixed limit, instead basing the BBL "//&
3252 "thickness on the bottom stress, rotation and stratification. The default is "//&
3253 "proportional to HBBL if USE_JACKSON_PARAM or DRAG_AS_BODY_FORCE is true.", &
3254 units="m", default=us%Z_to_m*chan_max_thick_dflt, scale=us%m_to_Z, &
3255 do_not_log=.not.cs%Channel_drag)
3256
3257 call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
3258 default=.false., do_not_log=.true.)
3259
3260 cs%Hbbl = cs%dz_bbl * (us%Z_to_m * gv%m_to_H) ! Rescaled for use in expressions in thickness units.
3261
3262 if (cs%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then
3263 ! This is necessary for reproducibility across restarts in non-symmetric mode.
3264 call pass_var(visc%Kv_shear_Bu, g%Domain, position=corner, complete=.true.)
3265 endif
3266
3267 if (cs%bottomdraglaw) then
3268 allocate(visc%bbl_thick_u(isdb:iedb,jsd:jed), source=0.0)
3269 allocate(visc%bbl_thick_v(isd:ied,jsdb:jedb), source=0.0)
3270 allocate(visc%kv_bbl_u(isdb:iedb,jsd:jed), source=0.0)
3271 allocate(visc%kv_bbl_v(isd:ied,jsdb:jedb), source=0.0)
3272 allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0)
3273 allocate(visc%BBL_meanKE_loss(isd:ied,jsd:jed), source=0.0)
3274 allocate(visc%BBL_meanKE_loss_sqrtCd(isd:ied,jsd:jed), source=0.0)
3275
3276 cs%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
3277 diag%axesCu1, time, 'BBL thickness at u points', 'm', conversion=us%Z_to_m)
3278 cs%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
3279 time, 'BBL viscosity at u points', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)
3280 cs%id_bbl_u = register_diag_field('ocean_model', 'bbl_u', diag%axesCu1, &
3281 time, 'BBL mean u current', 'm s-1', conversion=us%L_T_to_m_s)
3282 if (cs%id_bbl_u>0) then
3283 allocate(cs%bbl_u(isdb:iedb,jsd:jed), source=0.0)
3284 endif
3285 cs%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
3286 diag%axesCv1, time, 'BBL thickness at v points', 'm', conversion=us%Z_to_m)
3287 cs%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
3288 time, 'BBL viscosity at v points', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)
3289 cs%id_bbl_v = register_diag_field('ocean_model', 'bbl_v', diag%axesCv1, &
3290 time, 'BBL mean v current', 'm s-1', conversion=us%L_T_to_m_s)
3291 if (cs%id_bbl_v>0) then
3292 allocate(cs%bbl_v(isd:ied,jsdb:jedb), source=0.0)
3293 endif
3294 if (cs%bottomdragmap) then
3295 if (len_trim(cdrag_file)==0 .or. len_trim(cdrag_var)==0) then
3296 call mom_error(fatal,"CDRAG_FILE and CDRAG_VAR are required when using CDRAG_MAP.")
3297 endif
3298 allocate(cdrag_h(isd:ied,jsd:jed), source=0.0)
3299 allocate(cs%cdrag_u(isdb:iedb,jsd:jed), source=0.0)
3300 allocate(cs%cdrag_v(isd:ied,jsdb:jedb), source=0.0)
3301 filename = trim(cs%inputdir) // trim(cdrag_file)
3302 call log_param(param_file, mdl, "INPUTDIR/CDRAG_FILE", filename)
3303 call mom_read_data(filename, cdrag_var, cdrag_h, g%domain, scale=cs%cdrag)
3304 call pass_var(cdrag_h, g%domain)
3305 do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0) then
3306 cs%cdrag_u(i,j) = (g%mask2dT(i,j) * cdrag_h(i,j) + g%mask2dT(i+1,j) * cdrag_h(i+1,j)) / &
3307 (g%mask2dT(i,j) + g%mask2dT(i+1,j))
3308 endif ; enddo ; enddo
3309 do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0) then
3310 cs%cdrag_v(i,j) = (g%mask2dT(i,j) * cdrag_h(i,j) + g%mask2dT(i,j+1) * cdrag_h(i,j+1)) / &
3311 (g%mask2dT(i,j) + g%mask2dT(i,j+1))
3312 endif ; enddo ; enddo
3313 deallocate(cdrag_h)
3314 endif
3315 if (cs%BBL_use_tidal_bg) then
3316 allocate(cs%tideamp(isd:ied,jsd:jed), source=0.0)
3317 filename = trim(cs%inputdir) // trim(tideamp_file)
3318 call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
3319 call mom_read_data(filename, tideamp_var, cs%tideamp, g%domain, scale=us%m_s_to_L_T)
3320 call pass_var(cs%tideamp,g%domain)
3321 endif
3322 endif
3323 if (cs%Channel_drag .or. cs%body_force_drag) then
3324 allocate(visc%Ray_u(isdb:iedb,jsd:jed,nz), source=0.0)
3325 allocate(visc%Ray_v(isd:ied,jsdb:jedb,nz), source=0.0)
3326 cs%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
3327 time, 'Rayleigh drag velocity at u points', 'm s-1', conversion=gv%H_to_m*us%s_to_T)
3328 cs%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
3329 time, 'Rayleigh drag velocity at v points', 'm s-1', conversion=gv%H_to_m*us%s_to_T)
3330 endif
3331
3332
3333 if (cs%dynamic_viscous_ML) then
3334 allocate(visc%nkml_visc_u(isdb:iedb,jsd:jed), source=0.0)
3335 allocate(visc%nkml_visc_v(isd:ied,jsdb:jedb), source=0.0)
3336 cs%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
3337 diag%axesCu1, time, 'Number of layers in viscous mixed layer at u points', 'nondim')
3338 cs%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
3339 diag%axesCv1, time, 'Number of layers in viscous mixed layer at v points', 'nondim')
3340 endif
3341
3342 call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_cs)
3343 call register_restart_field_as_obsolete('Kv_turb','Kv_shear', restart_cs)
3344
3345end subroutine set_visc_init
3346
3347!> This subroutine dellocates any memory in the set_visc control structure.
3348subroutine set_visc_end(visc, CS)
3349 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
3350 !! related fields. Elements are deallocated here.
3351 type(set_visc_cs), intent(inout) :: cs !< The control structure returned by a previous
3352 !! call to set_visc_init.
3353
3354 if (allocated(visc%bbl_thick_u)) deallocate(visc%bbl_thick_u)
3355 if (allocated(visc%bbl_thick_v)) deallocate(visc%bbl_thick_v)
3356 if (allocated(visc%kv_bbl_u)) deallocate(visc%kv_bbl_u)
3357 if (allocated(visc%kv_bbl_v)) deallocate(visc%kv_bbl_v)
3358 if (allocated(cs%bbl_u)) deallocate(cs%bbl_u)
3359 if (allocated(cs%bbl_v)) deallocate(cs%bbl_v)
3360 if (allocated(visc%Ray_u)) deallocate(visc%Ray_u)
3361 if (allocated(visc%Ray_v)) deallocate(visc%Ray_v)
3362 if (allocated(visc%nkml_visc_u)) deallocate(visc%nkml_visc_u)
3363 if (allocated(visc%nkml_visc_v)) deallocate(visc%nkml_visc_v)
3364 if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear)
3365 if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow)
3366 if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
3367 if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear)
3368 if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
3369 if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
3370 if (allocated(visc%BBL_meanKE_loss)) deallocate(visc%BBL_meanKE_loss)
3371 if (allocated(visc%BBL_meanKE_loss_sqrtCd)) deallocate(visc%BBL_meanKE_loss_sqrtCd)
3372 if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf)
3373 if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
3374 if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
3375 if (allocated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
3376 if (allocated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
3377 if (allocated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)
3378end subroutine set_visc_end
3379
3380!> \namespace mom_set_visc
3381!!
3382!! This would also be the module in which other viscous quantities that are flow-independent might be set.
3383!! This information is transmitted to other modules via a vertvisc type structure.
3384!!
3385!! The same code is used for the two velocity components, by indirectly referencing the velocities and
3386!! defining a handful of direction-specific defined variables.
3387
3388end module mom_set_visc