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