MOM_bulk_mixed_layer.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!> Build mixed layer parameterization
7
8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
9use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc
11use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
12use mom_eos, only : calculate_density, calculate_density_derivs, eos_domain
13use mom_eos, only : average_specific_vol, calculate_density_derivs
14use mom_eos, only : calculate_spec_vol, calculate_specific_vol_derivs
15use mom_error_handler, only : mom_error, fatal, warning
16use mom_file_parser, only : get_param, log_param, log_version, param_file_type
17use mom_forcing_type, only : extractfluxes1d, forcing, find_ustar
18use mom_grid, only : ocean_grid_type
23
24implicit none ; private
25
26#include <MOM_memory.h>
27
29
30! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34
35!> The control structure with parameters for the MOM_bulk_mixed_layer module
36type, public :: bulkmixedlayer_cs ; private
37 logical :: initialized = .false. !< True if this control structure has been initialized.
38 integer :: nkml !< The number of layers in the mixed layer.
39 integer :: nkbl !< The number of buffer layers.
40 integer :: nsw !< The number of bands of penetrating shortwave radiation.
41 real :: mstar !< The ratio of the friction velocity cubed to the
42 !! TKE input to the mixed layer [nondim].
43 real :: nstar !< The fraction of the TKE input to the mixed layer
44 !! available to drive entrainment [nondim].
45 real :: nstar2 !< The fraction of potential energy released by
46 !! convective adjustment that drives entrainment [nondim].
47 logical :: absorb_all_sw !< If true, all shortwave radiation is absorbed by the
48 !! ocean, instead of passing through to the bottom mud.
49 real :: tke_decay !< The ratio of the natural Ekman depth to the TKE
50 !! decay scale [nondim].
51 real :: bulk_ri_ml !< The efficiency with which mean kinetic energy released by
52 !! mechanically forced entrainment of the mixed layer is
53 !! converted to TKE, times conversion factors between the
54 !! natural units of mean kinetic energy and TKE [Z2 L-2 ~> nondim]
55 real :: bulk_ri_convective !< The efficiency with which convectively released mean kinetic
56 !! energy becomes TKE, times conversion factors between the natural
57 !! units of mean kinetic energy and TKE [Z2 L-2 ~> nondim]
58 real :: vonkar !< The von Karman constant as used for mixed layer viscosity [nondim]
59 real :: hmix_min !< The minimum mixed layer thickness [H ~> m or kg m-2].
60 real :: mech_tke_floor !< A tiny floor on the amount of turbulent kinetic energy that is
61 !! used when the mixed layer does not yet contain HMIX_MIN fluid
62 !! [H Z2 T-2 ~> m3 s-2 or J m-2]. The default is so small that its actual
63 !! value is irrelevant, but it is detectably greater than 0.
64 real :: h_limit_fluxes !< When the total ocean depth is less than this
65 !! value [H ~> m or kg m-2], scale away all surface forcing to
66 !! avoid boiling the ocean.
67 real :: ustar_min !< A minimum value of ustar to avoid numerical problems [Z T-1 ~> m s-1].
68 !! If the value is small enough, this should not affect the solution.
69 real :: omega !< The Earth's rotation rate [T-1 ~> s-1].
70 real :: dt_ds_wt !< When forced to extrapolate T & S to match the
71 !! layer densities, this factor [C S-1 ~> degC ppt-1] is
72 !! combined with the derivatives of density with T & S
73 !! to determines what direction is orthogonal to
74 !! density contours. It should be a typical value of
75 !! (dR/dS) / (dR/dT) in oceanic profiles.
76 !! 6 degC ppt-1 might be reasonable.
77 real :: hbuffer_min !< The minimum buffer layer thickness when the mixed layer
78 !! is very large [H ~> m or kg m-2].
79 real :: hbuffer_rel_min !< The minimum buffer layer thickness relative to the combined
80 !! mixed and buffer layer thicknesses when they are thin [nondim]
81 real :: bl_detrain_time !< A timescale that characterizes buffer layer detrainment
82 !! events [T ~> s].
83 real :: bl_extrap_lim !< A limit on the density range over which
84 !! extrapolation can occur when detraining from the
85 !! buffer layers, relative to the density range
86 !! within the mixed and buffer layers, when the
87 !! detrainment is going into the lightest interior
88 !! layer [nondim].
89 real :: bl_split_rho_tol !< The fractional tolerance for matching layer target densities
90 !! when splitting layers to deal with massive interior layers
91 !! that are lighter than one of the mixed or buffer layers [nondim].
92 logical :: ml_resort !< If true, resort the layers by density, rather than
93 !! doing convective adjustment.
94 integer :: ml_presort_nz_conv_adj !< If ML_resort is true, do convective
95 !! adjustment on this many layers (starting from the
96 !! top) before sorting the remaining layers.
97 real :: omega_frac !< When setting the decay scale for turbulence, use this fraction
98 !! of the absolute rotation rate blended with the local value of f,
99 !! as sqrt((1-of)*f^2 + of*4*omega^2) [nondim].
100 logical :: correct_absorption !< If true, the depth at which penetrating
101 !! shortwave radiation is absorbed is corrected by
102 !! moving some of the heating upward in the water
103 !! column. The default is false.
104 logical :: nonbous_energetics !< If true, use non-Boussinesq expressions for the energetic
105 !! calculations used in the bulk mixed layer calculations.
106 logical :: resolve_ekman !< If true, the nkml layers in the mixed layer are
107 !! chosen to optimally represent the impact of the
108 !! Ekman transport on the mixed layer TKE budget.
109 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
110 logical :: tke_diagnostics = .false. !< If true, calculate extensive diagnostics of the TKE budget
111 logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff
112 !! at the river mouths to rivermix_depth
113 real :: rivermix_depth = 0.0 !< The depth of mixing if do_rivermix is true [H ~> m or kg m-2].
114 logical :: limit_det !< If true, limit the extent of buffer layer
115 !! detrainment to be consistent with neighbors.
116 real :: lim_det_dh_sfc !< The fractional limit in the change between grid
117 !! points of the surface region (mixed & buffer
118 !! layer) thickness [nondim]. 0.5 by default.
119 real :: lim_det_dh_bathy !< The fraction of the total depth by which the
120 !! thickness of the surface region (mixed & buffer layers) is allowed
121 !! to change between grid points [nondim]. 0.2 by default.
122 logical :: use_river_heat_content !< If true, use the fluxes%runoff_Hflx field
123 !! to set the heat carried by runoff, instead of
124 !! using SST for temperature of liq_runoff
125 logical :: use_calving_heat_content !< Use SST for temperature of froz_runoff
126 logical :: convect_mom_bug !< If true, use code with a bug that causes a loss of momentum
127 !! conservation during mixedlayer convection.
128
129 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the
130 !! timing of diagnostic output.
131 real :: allowed_t_chg !< The amount by which temperature is allowed
132 !! to exceed previous values during detrainment [C ~> degC]
133 real :: allowed_s_chg !< The amount by which salinity is allowed
134 !! to exceed previous values during detrainment [S ~> ppt]
135
136 ! These are terms in the mixed layer TKE budget, all in [H Z2 T-3 ~> m3 s-3 or W m-2] except as noted.
137 real, allocatable, dimension(:,:) :: &
138 ml_depth, & !< The mixed layer depth [H ~> m or kg m-2].
139 diag_tke_wind, & !< The wind source of TKE [H Z2 T-3 ~> m3 s-3 or W m-2].
140 diag_tke_ribulk, & !< The resolved KE source of TKE [H Z2 T-3 ~> m3 s-3 or W m-2].
141 diag_tke_conv, & !< The convective source of TKE [H Z2 T-3 ~> m3 s-3 or W m-2].
142 diag_tke_pen_sw, & !< The TKE sink required to mix penetrating shortwave heating [H Z2 T-3 ~> m3 s-3 or W m-2].
143 diag_tke_mech_decay, & !< The decay of mechanical TKE [H Z2 T-3 ~> m3 s-3 or W m-2].
144 diag_tke_conv_decay, & !< The decay of convective TKE [H Z2 T-3 ~> m3 s-3 or W m-2].
145 diag_tke_mixing, & !< The work done by TKE to deepen the mixed layer [H Z2 T-3 ~> m3 s-3 or W m-2].
146 diag_tke_conv_s2, & !< The convective source of TKE due to to mixing in sigma2 [H Z2 T-3 ~> m3 s-3 or W m-2].
147 diag_pe_detrain, & !< The spurious source of potential energy due to mixed layer
148 !! detrainment [R Z3 T-3 ~> W m-2].
149 diag_pe_detrain2 !< The spurious source of potential energy due to mixed layer only
150 !! detrainment [R Z3 T-3 ~> W m-2].
151 type(group_pass_type) :: pass_h_sum_hmbl_prev !< For group halo pass
152
153 !>@{ Diagnostic IDs
154 integer :: id_ml_depth = -1, id_tke_wind = -1, id_tke_mixing = -1
155 integer :: id_tke_ribulk = -1, id_tke_conv = -1, id_tke_pen_sw = -1
156 integer :: id_tke_mech_decay = -1, id_tke_conv_decay = -1, id_tke_conv_s2 = -1
157 integer :: id_pe_detrain = -1, id_pe_detrain2 = -1, id_h_mismatch = -1
158 integer :: id_hsfc_used = -1, id_hsfc_max = -1, id_hsfc_min = -1
159 !>@}
160end type bulkmixedlayer_cs
161
162!>@{ CPU clock IDs
163integer :: id_clock_pass=0
164!>@}
165
166contains
167
168!> This subroutine partially steps the bulk mixed layer model.
169!! See \ref BML for more details.
170subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, CS, &
171 optics, BLD, H_ml, aggregate_FW_forcing, dt_diag, last_call)
172 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
173 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
174 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
175 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
176 intent(inout) :: h_3d !< Layer thickness [H ~> m or kg m-2].
177 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
178 intent(in) :: u_3d !< Zonal velocities interpolated to h points
179 !! [L T-1 ~> m s-1].
180 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
181 intent(in) :: v_3d !< Zonal velocities interpolated to h points
182 !! [L T-1 ~> m s-1].
183 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
184 !! available thermodynamic fields. Absent
185 !! fields have NULL pointers.
186 type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
187 !! possible forcing fields. Unused fields
188 !! have NULL pointers.
189 real, intent(in) :: dt !< Time increment [T ~> s].
190 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
191 intent(inout) :: ea !< The amount of fluid moved downward into a
192 !! layer; this should be increased due to
193 !! mixed layer detrainment [H ~> m or kg m-2].
194 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
195 intent(inout) :: eb !< The amount of fluid moved upward into a
196 !! layer; this should be increased due to
197 !! mixed layer entrainment [H ~> m or kg m-2].
198 type(bulkmixedlayer_cs), intent(inout) :: cs !< Bulk mixed layer control structure
199 type(optics_type), pointer :: optics !< The structure that can be queried for the
200 !! inverse of the vertical absorption decay
201 !! scale for penetrating shortwave radiation.
202 real, dimension(SZI_(G),SZJ_(G)), &
203 intent(inout) :: bld !< Active mixed layer depth [Z ~> m]
204 real, dimension(SZI_(G),SZJ_(G)), &
205 intent(inout) :: h_ml !< Active mixed layer thickness [H ~> m or kg m-2].
206 logical, intent(in) :: aggregate_fw_forcing !< If true, the net incoming and
207 !! outgoing surface freshwater fluxes are
208 !! combined before being applied, instead of
209 !! being applied separately.
210 real, optional, intent(in) :: dt_diag !< The diagnostic time step,
211 !! which may be less than dt if there are
212 !! two calls to mixedlayer [T ~> s].
213 logical, optional, intent(in) :: last_call !< if true, this is the last call
214 !! to mixedlayer in the current time step, so
215 !! diagnostics will be written. The default is
216 !! .true.
217
218 ! Local variables
219 real, dimension(SZI_(G),SZK_(GV)) :: &
220 eaml, & ! The amount of fluid moved downward into a layer due to mixed
221 ! layer detrainment [H ~> m or kg m-2]. (I.e. entrainment from above.)
222 ebml ! The amount of fluid moved upward into a layer due to mixed
223 ! layer detrainment [H ~> m or kg m-2]. (I.e. entrainment from below.)
224
225 ! If there is resorting, the vertical coordinate for these variables is the
226 ! new, sorted index space. Here layer 0 is an initially massless layer that
227 ! will be used to hold the new mixed layer properties.
228 real, dimension(SZI_(G),SZK0_(GV)) :: &
229 h, & ! The layer thickness [H ~> m or kg m-2].
230 t, & ! The layer temperatures [C ~> degC].
231 s, & ! The layer salinities [S ~> ppt].
232 r0, & ! The potential density referenced to the surface [R ~> kg m-3].
233 spv0, & ! The specific volume referenced to the surface [R-1 ~> m3 kg-1].
234 rcv ! The coordinate variable potential density [R ~> kg m-3].
235 real, dimension(SZI_(G),SZK_(GV)) :: &
236 u, & ! The zonal velocity [L T-1 ~> m s-1].
237 v, & ! The meridional velocity [L T-1 ~> m s-1].
238 h_orig, & ! The original thickness [H ~> m or kg m-2].
239 d_eb, & ! The downward increase across a layer in the entrainment from
240 ! below [H ~> m or kg m-2]. The sign convention is that positive values of
241 ! d_eb correspond to a gain in mass by a layer by upward motion.
242 d_ea, & ! The upward increase across a layer in the entrainment from
243 ! above [H ~> m or kg m-2]. The sign convention is that positive values of
244 ! d_ea mean a net gain in mass by a layer from downward motion.
245 eps ! The (small) thickness that must remain in a layer [H ~> m or kg m-2].
246 integer, dimension(SZI_(G),SZK_(GV)) :: &
247 ksort ! The sorted k-index that each original layer goes to.
248 real, dimension(SZI_(G),SZJ_(G)) :: &
249 h_miss ! The summed absolute mismatch [H ~> m or kg m-2].
250 real, dimension(SZI_(G),SZJ_(G)) :: &
251 u_star_2d, &! The wind friction velocity, calculated using the Boussinesq reference density or
252 ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
253 u_star_h_2d ! The wind friction velocity in thickness-based units, calculated
254 ! using the Boussinesq reference density or the time-evolving
255 ! surface density in non-Boussinesq mode [H T-1 ~> m s-1 or kg m-2 s-1]
256 real, dimension(SZI_(G)) :: &
257 tke, & ! The turbulent kinetic energy available for mixing over a
258 ! time step [H Z2 T-2 ~> m3 s-2 or J m-2].
259 conv_en, & ! The turbulent kinetic energy source due to mixing down to
260 ! the depth of free convection [H Z2 T-2 ~> m3 s-2 or J m-2].
261 htot, & ! The total depth of the layers being considered for
262 ! entrainment [H ~> m or kg m-2].
263 r0_tot, & ! The integrated potential density referenced to the surface
264 ! of the layers which are fully entrained [H R ~> kg m-2 or kg2 m-5].
265 spv0_tot, & ! The integrated specific volume referenced to the surface
266 ! of the layers which are fully entrained [H R-1 ~> m4 kg-1 or m].
267 rcv_tot, & ! The integrated coordinate value potential density of the
268 ! layers that are fully entrained [H R ~> kg m-2 or kg2 m-5].
269 ttot, & ! The integrated temperature of layers which are fully
270 ! entrained [C H ~> degC m or degC kg m-2].
271 stot, & ! The integrated salt of layers which are fully entrained
272 ! [H S ~> m ppt or ppt kg m-2].
273 uhtot, & ! The depth integrated zonal velocity in the mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1]
274 vhtot, & ! The depth integrated meridional velocity in the mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1]
275
276 netmassinout, & ! The net mass flux (if non-Boussinesq) or volume flux (if
277 ! Boussinesq - i.e. the fresh water flux (P+R-E)) into the
278 ! ocean over a time step [H ~> m or kg m-2].
279 netmassout, & ! The mass flux (if non-Boussinesq) or volume flux (if Boussinesq)
280 ! over a time step from evaporating fresh water [H ~> m or kg m-2]
281 net_heat, & ! The net heating at the surface over a time step [C H ~> degC m or degC kg m-2]
282 ! Any penetrating shortwave radiation is not included in Net_heat.
283 net_salt, & ! The surface salt flux into the ocean over a time step [S H ~> ppt m or ppt kg m-2]
284 idecay_len_tke, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
285 p_ref, & ! Reference pressure for the potential density governing mixed
286 ! layer dynamics, almost always 0 (or 1e5) [R L2 T-2 ~> Pa].
287 p_ref_cv, & ! Reference pressure for the potential density which defines
288 ! the coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
289 dr0_dt, & ! Partial derivative of the mixed layer potential density with
290 ! temperature [R C-1 ~> kg m-3 degC-1].
291 dspv0_dt, & ! Partial derivative of the mixed layer specific volume with
292 ! temperature [R-1 C-1 ~> m3 kg-1 degC-1].
293 drcv_dt, & ! Partial derivative of the coordinate variable potential
294 ! density in the mixed layer with temperature [R C-1 ~> kg m-3 degC-1].
295 dr0_ds, & ! Partial derivative of the mixed layer potential density with
296 ! salinity [R S-1 ~> kg m-3 ppt-1].
297 dspv0_ds, & ! Partial derivative of the mixed layer specific volume with
298 ! salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
299 drcv_ds, & ! Partial derivative of the coordinate variable potential
300 ! density in the mixed layer with salinity [R S-1 ~> kg m-3 ppt-1].
301 p_sfc, & ! The sea surface pressure [R L2 T-2 ~> Pa]
302 dp_ml, & ! The pressure change across the mixed layer [R L2 T-2 ~> Pa]
303 spv_ml, & ! The specific volume averaged across the mixed layer [R-1 ~> m3 kg-1]
304 tke_river ! The source of turbulent kinetic energy available for mixing
305 ! at rivermouths [H Z2 T-3 ~> m3 s-3 or W m-2].
306
307 real, dimension(max(CS%nsw,1),SZI_(G)) :: &
308 pen_sw_bnd ! The penetrating fraction of the shortwave heating integrated
309 ! over a time step in each band [C H ~> degC m or degC kg m-2].
310 real, dimension(max(CS%nsw,1),SZI_(G),SZK_(GV)) :: &
311 opacity_band ! The opacity in each band [H-1 ~> m-1 or m2 kg-1]. The indices are band, i, k.
312
313 real :: cmke(2,szi_(g)) ! Coefficients of HpE and HpE^2 used in calculating the
314 ! denominator of MKE_rate; the two elements have differing
315 ! units of [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
316 real :: irho0 ! 1.0 / rho_0 [R-1 ~> m3 kg-1]
317 real :: inkml, inkmlm1! 1.0 / REAL(nkml) and 1.0 / REAL(nkml-1) [nondim]
318 real :: ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
319 real :: idt_diag ! The inverse of the timestep used for diagnostics [T-1 ~> s-1].
320 real :: rmixconst ! A combination of constants used in the river mixing energy
321 ! calculation [H Z T-2 R-2 ~> m8 s-2 kg-2 or m5 s-2 kg-1] or
322 ! [H Z T-2 ~> m2 s-2 or kg m-1 s-2]
323 real, dimension(SZI_(G)) :: &
324 dke_fc, & ! The change in mean kinetic energy due to free convection
325 ! [H Z2 T-2 ~> m3 s-2 or J m-2].
326 h_ca ! The depth to which convective adjustment has gone [H ~> m or kg m-2].
327 real, dimension(SZI_(G),SZK_(GV)) :: &
328 dke_ca, & ! The change in mean kinetic energy due to convective
329 ! adjustment [H Z2 T-2 ~> m3 s-2 or J m-2].
330 ctke ! The turbulent kinetic energy source due to convective
331 ! adjustment [H Z2 T-2 ~> m3 s-2 or J m-2].
332 real, dimension(SZI_(G),SZJ_(G)) :: &
333 hsfc_max, & ! The thickness of the surface region (mixed and buffer layers)
334 ! after entrainment but before any buffer layer detrainment [H ~> m or kg m-2].
335 hsfc_used, & ! The thickness of the surface region after buffer layer
336 ! detrainment [H ~> m or kg m-2].
337 hsfc_min, & ! The minimum thickness of the surface region based on the
338 ! new mixed layer depth and the previous thickness of the
339 ! neighboring water columns [H ~> m or kg m-2].
340 h_sum, & ! The total thickness of the water column [H ~> m or kg m-2].
341 hmbl_prev ! The previous thickness of the mixed and buffer layers [H ~> m or kg m-2].
342 real, dimension(SZI_(G)) :: &
343 hsfc, & ! The thickness of the surface region (mixed and buffer
344 ! layers before detrainment in to the interior [H ~> m or kg m-2].
345 max_bl_det ! If non-negative, the maximum amount of entrainment from
346 ! the buffer layers that will be allowed this time step [H ~> m or kg m-2].
347 real :: dhsfc, dhd ! Local copies of nondimensional parameters [nondim]
348 real :: h_nbr ! A minimum thickness based on neighboring thicknesses [H ~> m or kg m-2].
349
350 real :: absf_x_h ! The absolute value of f times the mixed layer thickness [H T-1 ~> m s-1 or kg m-2 s-1].
351 real :: ku_star ! Ustar times the Von Karman constant [H T-1 ~> m s-1 or kg m-2 s-1].
352 real :: dt__diag ! A rescaled copy of dt_diag (if present) or dt [T ~> s].
353 logical :: write_diags ! If true, write out diagnostics with this step.
354 logical :: reset_diags ! If true, zero out the accumulated diagnostics.
355 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
356 integer :: i, j, k, is, ie, js, je, nz, nkmb
357 integer :: nsw ! The number of bands of penetrating shortwave radiation.
358
359 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
360
361 if (.not. cs%initialized) call mom_error(fatal, "MOM_bulk_mixed_layer: "//&
362 "Module must be initialized before it is used.")
363 if (gv%nkml < 1) return
364
365 if (.not. associated(tv%eqn_of_state)) call mom_error(fatal, &
366 "MOM_mixed_layer: Temperature, salinity and an equation of state "//&
367 "must now be used.")
368 if (.not. (associated(fluxes%ustar) .or. associated(fluxes%tau_mag))) call mom_error(fatal, &
369 "MOM_mixed_layer: No surface TKE fluxes (ustar or tau_mag) defined in mixedlayer!")
370
371 nkmb = cs%nkml+cs%nkbl
372 inkml = 1.0 / real(cs%nkml)
373 if (cs%nkml > 1) inkmlm1 = 1.0 / real(cs%nkml-1)
374
375 irho0 = 1.0 / gv%Rho0
376 dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
377 idt_diag = 1.0 / dt__diag
378 write_diags = .true. ; if (present(last_call)) write_diags = last_call
379
380 p_ref(:) = 0.0 ; p_ref_cv(:) = tv%P_Ref
381
382 nsw = cs%nsw
383
384 if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
385 !$OMP parallel do default(shared)
386 do j=js-1,je+1 ; do i=is-1,ie+1
387 h_sum(i,j) = 0.0 ; hmbl_prev(i,j) = 0.0
388 enddo ; enddo
389 !$OMP parallel do default(shared)
390 do j=js-1,je+1
391 do k=1,nkmb ; do i=is-1,ie+1
392 h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
393 hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k)
394 enddo ; enddo
395 do k=nkmb+1,nz ; do i=is-1,ie+1
396 h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
397 enddo ; enddo
398 enddo
399
400 call cpu_clock_begin(id_clock_pass)
401 call create_group_pass(cs%pass_h_sum_hmbl_prev, h_sum,g%Domain)
402 call create_group_pass(cs%pass_h_sum_hmbl_prev, hmbl_prev,g%Domain)
403 call do_group_pass(cs%pass_h_sum_hmbl_prev, g%Domain)
404 call cpu_clock_end(id_clock_pass)
405 endif
406
407 ! Determine whether to zero out diagnostics before accumulation.
408 reset_diags = .true.
409 if (present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
410 reset_diags = .false. ! This is the second call to mixedlayer.
411
412 if (reset_diags) then
413 if (cs%TKE_diagnostics) then
414 !$OMP parallel do default(shared)
415 do j=js,je ; do i=is,ie
416 cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_RiBulk(i,j) = 0.0
417 cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_pen_SW(i,j) = 0.0
418 cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
419 cs%diag_TKE_conv_decay(i,j) = 0.0 ; cs%diag_TKE_conv_s2(i,j) = 0.0
420 enddo ; enddo
421 endif
422 if (allocated(cs%diag_PE_detrain)) then
423 !$OMP parallel do default(shared)
424 do j=js,je ; do i=is,ie
425 cs%diag_PE_detrain(i,j) = 0.0
426 enddo ; enddo
427 endif
428 if (allocated(cs%diag_PE_detrain2)) then
429 !$OMP parallel do default(shared)
430 do j=js,je ; do i=is,ie
431 cs%diag_PE_detrain2(i,j) = 0.0
432 enddo ; enddo
433 endif
434 endif
435
436 if (cs%ML_resort) then
437 do i=is,ie ; h_ca(i) = 0.0 ; enddo
438 do k=1,nz ; do i=is,ie ; dke_ca(i,k) = 0.0 ; ctke(i,k) = 0.0 ; enddo ; enddo
439 endif
440 max_bl_det(:) = -1
441 eosdom(:) = eos_domain(g%HI)
442
443 ! Extract the friction velocity from the forcing type.
444 call find_ustar(fluxes, tv, u_star_2d, g, gv, us)
445 if (cs%Resolve_Ekman .and. (cs%nkml>1)) &
446 call find_ustar(fluxes, tv, u_star_h_2d, g, gv, us, h_t_units=.true.)
447
448 !$OMP parallel default(shared) firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) &
449 !$OMP private(h,u,v,h_orig,eps,T,S,opacity_band,d_ea,d_eb,R0,SpV0,Rcv,ksort, &
450 !$OMP dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,dSpV0_dT,dSpV0_dS,htot,Ttot,Stot,TKE,Conv_en, &
451 !$OMP RmixConst,TKE_river,Pen_SW_bnd,netMassInOut,NetMassOut, &
452 !$OMP Net_heat,Net_salt,uhtot,vhtot,R0_tot,Rcv_tot,SpV0_tot,dKE_FC, &
453 !$OMP Idecay_len_TKE,cMKE,Hsfc,dHsfc,dHD,H_nbr,kU_Star, &
454 !$OMP absf_x_H,ebml,eaml)
455 !$OMP do
456 do j=js,je
457 ! Copy the thicknesses and other fields to 2-d arrays.
458 do k=1,nz ; do i=is,ie
459 h(i,k) = h_3d(i,j,k) ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
460 h_orig(i,k) = h_3d(i,j,k)
461 eps(i,k) = 0.0 ; if (k > nkmb) eps(i,k) = gv%Angstrom_H
462 t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
463 enddo ; enddo
464 if (nsw>0) then
465 if (gv%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
466 call extract_optics_slice(optics, j, g, gv, opacity=opacity_band, opacity_scale=gv%H_to_Z)
467 else
468 call extract_optics_slice(optics, j, g, gv, opacity=opacity_band, opacity_scale=gv%H_to_RZ, &
469 spv_avg=tv%SpV_avg)
470 endif
471 endif
472
473 do k=1,nz ; do i=is,ie
474 d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0
475 enddo ; enddo
476
477 ! Calculate an estimate of the mid-mixed layer pressure [R L2 T-2 ~> Pa]
478 if (associated(tv%p_surf)) then
479 do i=is,ie ; p_ref(i) = tv%p_surf(i,j) ; enddo
480 else
481 do i=is,ie ; p_ref(i) = 0.0 ; enddo
482 endif
483 do k=1,cs%nkml ; do i=is,ie
484 p_ref(i) = p_ref(i) + 0.5*(gv%H_to_RZ*gv%g_Earth)*h(i,k)
485 enddo ; enddo
486 if (cs%nonBous_energetics) then
487 call calculate_specific_vol_derivs(t(:,1), s(:,1), p_ref, dspv0_dt, dspv0_ds, tv%eqn_of_state, eosdom)
488 do k=1,nz
489 call calculate_spec_vol(t(:,k), s(:,k), p_ref, spv0(:,k), tv%eqn_of_state, eosdom)
490 enddo
491 else
492 call calculate_density_derivs(t(:,1), s(:,1), p_ref, dr0_dt, dr0_ds, tv%eqn_of_state, eosdom)
493 do k=1,nz
494 call calculate_density(t(:,k), s(:,k), p_ref, r0(:,k), tv%eqn_of_state, eosdom)
495 enddo
496 endif
497 call calculate_density_derivs(t(:,1), s(:,1), p_ref_cv, drcv_dt, drcv_ds, tv%eqn_of_state, eosdom)
498 do k=1,nz
499 call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), tv%eqn_of_state, eosdom)
500 enddo
501
502 if (cs%ML_resort) then
503 if (cs%ML_presort_nz_conv_adj > 0) &
504 call convective_adjustment(h, u, v, r0, spv0, rcv, t, s, eps, d_eb, dke_ca, ctke, j, g, gv, &
505 us, cs, cs%ML_presort_nz_conv_adj)
506
507 call sort_ml(h, r0, spv0, eps, g, gv, cs, ksort)
508 else
509 do k=1,nz ; do i=is,ie ; ksort(i,k) = k ; enddo ; enddo
510
511 ! Undergo instantaneous entrainment into the buffer layers and mixed layers
512 ! to remove hydrostatic instabilities. Any water that is lighter than
513 ! currently in the mixed or buffer layer is entrained.
514 call convective_adjustment(h, u, v, r0, spv0, rcv, t, s, eps, d_eb, dke_ca, ctke, j, g, gv, us, cs)
515 do i=is,ie ; h_ca(i) = h(i,1) ; enddo
516
517 endif
518
519 if (associated(fluxes%lrunoff) .and. cs%do_rivermix) then
520
521 ! Here we add an additional source of TKE to the mixed layer where river
522 ! is present to simulate unresolved estuaries. The TKE input is diagnosed
523 ! as follows:
524 ! TKE_river[H Z2 T-3 ~> m3 s-3] = 0.5*rivermix_depth * g * Irho0**2 * drho_ds *
525 ! River*(Samb - Sriver) = CS%mstar*U_star^3
526 ! where River is in units of [R Z T-1 ~> kg m-2 s-1].
527 ! Samb = Ambient salinity at the mouth of the estuary
528 ! rivermix_depth = The prescribed depth over which to mix river inflow
529 ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
530 ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
531 if (cs%nonBous_energetics) then
532 rmixconst = -0.5*cs%rivermix_depth * gv%g_Earth_Z_T2
533 do i=is,ie
534 tke_river(i) = max(0.0, rmixconst * dspv0_ds(i) * &
535 ((fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) + &
536 (fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) * s(i,1))
537 enddo
538 else
539 rmixconst = 0.5*cs%rivermix_depth * gv%g_Earth_Z_T2 * irho0**2
540 do i=is,ie
541 tke_river(i) = max(0.0, rmixconst*dr0_ds(i)* &
542 ((fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) + &
543 (fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) * s(i,1))
544 enddo
545 endif
546 else
547 do i=is,ie ; tke_river(i) = 0.0 ; enddo
548 endif
549
550 ! The surface forcing is contained in the fluxes type.
551 ! We aggregate the thermodynamic forcing for a time step into the following:
552 ! netMassInOut = water [H ~> m or kg m-2] added/removed via surface fluxes
553 ! netMassOut = water [H ~> m or kg m-2] removed via evaporating surface fluxes
554 ! net_heat = heat via surface fluxes [C H ~> degC m or degC kg m-2]
555 ! net_salt = salt via surface fluxes [S H ~> ppt m or gSalt m-2]
556 ! Pen_SW_bnd = components to penetrative shortwave radiation
557 call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, &
558 cs%H_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
559 h(:,1:), t(:,1:), netmassinout, netmassout, net_heat, net_salt, pen_sw_bnd, &
560 tv, aggregate_fw_forcing)
561
562 ! This subroutine causes the mixed layer to entrain to depth of free convection.
563 call mixedlayer_convection(h, d_eb, htot, ttot, stot, uhtot, vhtot, r0_tot, spv0_tot, rcv_tot, &
564 u, v, t, s, r0, spv0, rcv, eps, dr0_dt, dspv0_dt, drcv_dt, dr0_ds, dspv0_ds, drcv_ds, &
565 netmassinout, netmassout, net_heat, net_salt, &
566 nsw, pen_sw_bnd, opacity_band, conv_en, dke_fc, &
567 j, ksort, g, gv, us, cs, tv, fluxes, dt, aggregate_fw_forcing)
568
569 ! Now the mixed layer undergoes mechanically forced entrainment.
570 ! The mixed layer may entrain down to the Monin-Obukhov depth if the
571 ! surface is becoming lighter, and is effectively detraining.
572
573 ! First the TKE at the depth of free convection that is available
574 ! to drive mixing is calculated.
575 call find_starting_tke(htot, h_ca, fluxes, u_star_2d, conv_en, ctke, dke_fc, dke_ca, &
576 tke, tke_river, idecay_len_tke, cmke, tv, dt, idt_diag, &
577 j, ksort, g, gv, us, cs)
578
579 ! Here the mechanically driven entrainment occurs.
580 call mechanical_entrainment(h, d_eb, htot, ttot, stot, uhtot, vhtot, &
581 r0_tot, spv0_tot, rcv_tot, u, v, t, s, r0, spv0, rcv, eps, &
582 dr0_dt, dspv0_dt, drcv_dt, cmke, idt_diag, nsw, pen_sw_bnd, &
583 opacity_band, tke, idecay_len_tke, j, ksort, g, gv, us, cs)
584
585 call absorbremainingsw(g, gv, us, h(:,1:), opacity_band, nsw, optics, j, dt, &
586 cs%H_limit_fluxes, cs%correct_absorption, cs%absorb_all_SW, &
587 t(:,1:), pen_sw_bnd, eps, ksort, htot, ttot)
588
589 if (cs%TKE_diagnostics) then ; do i=is,ie
590 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) - idt_diag * tke(i)
591 enddo ; endif
592
593 ! Calculate the homogeneous mixed layer properties and store them in layer 0.
594 do i=is,ie ; if (htot(i) > 0.0) then
595 ih = 1.0 / htot(i)
596 if (cs%nonBous_energetics) then
597 spv0(i,0) = spv0_tot(i) * ih
598 else
599 r0(i,0) = r0_tot(i) * ih
600 endif
601 rcv(i,0) = rcv_tot(i) * ih
602 t(i,0) = ttot(i) * ih ; s(i,0) = stot(i) * ih
603 h(i,0) = htot(i)
604 else ! This may not ever be needed?
605 t(i,0) = t(i,1) ; s(i,0) = s(i,1) ; rcv(i,0) = rcv(i,1)
606 if (cs%nonBous_energetics) then
607 spv0(i,0) = spv0(i,1)
608 else
609 r0(i,0) = r0(i,1)
610 endif
611 h(i,0) = htot(i)
612 endif ; enddo
613 if (write_diags .and. allocated(cs%ML_depth)) then ; do i=is,ie
614 cs%ML_depth(i,j) = h(i,0) ! Store the diagnostic.
615 enddo ; endif
616
617 ! Return the mixed layer depth in [Z ~> m].
618 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
619 do i=is,ie
620 bld(i,j) = g%mask2dT(i,j) * gv%H_to_Z*h(i,0)
621 enddo
622 else
623 do i=is,ie ; dp_ml(i) = gv%g_Earth * gv%H_to_RZ * h(i,0) ; enddo
624 if (associated(tv%p_surf)) then
625 do i=is,ie ; p_sfc(i) = tv%p_surf(i,j) ; enddo
626 else
627 do i=is,ie ; p_sfc(i) = 0.0 ; enddo
628 endif
629 call average_specific_vol(t(:,0), s(:,0), p_sfc, dp_ml, spv_ml, tv%eqn_of_state)
630 do i=is,ie
631 bld(i,j) = g%mask2dT(i,j) * gv%H_to_RZ * spv_ml(i) * h(i,0)
632 enddo
633 endif
634 ! Return the mixed layer thickness in [H ~> m or kg m-2].
635 do i=is,ie
636 h_ml(i,j) = g%mask2dT(i,j) * h(i,0)
637 enddo
638
639! At this point, return water to the original layers, but constrained to
640! still be sorted. After this point, all the water that is in massive
641! interior layers will be denser than water remaining in the mixed- and
642! buffer-layers. To achieve this, some of these variable density layers
643! might be split between two isopycnal layers that are denser than new
644! mixed layer or any remaining water from the old mixed- or buffer-layers.
645! Alternately, if there are fewer than nkbl of the old buffer or mixed layers
646! with any mass, relatively light interior layers might be transferred to
647! these unused layers (but not currently in the code).
648
649 if (cs%ML_resort) then
650 call resort_ml(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), spv0(:,0:), rcv(:,0:), gv%Rlay(:), eps, &
651 d_ea, d_eb, ksort, g, gv, cs, dr0_dt, dr0_ds, dspv0_dt, dspv0_ds, drcv_dt, drcv_ds)
652 endif
653
654 if (cs%limit_det .or. (cs%id_Hsfc_max > 0) .or. (cs%id_Hsfc_min > 0)) then
655 do i=is,ie ; hsfc(i) = h(i,0) ; enddo
656 do k=1,nkmb ; do i=is,ie ; hsfc(i) = hsfc(i) + h(i,k) ; enddo ; enddo
657
658 if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
659 dhsfc = cs%lim_det_dH_sfc ; dhd = cs%lim_det_dH_bathy
660 do i=is,ie
661 h_nbr = min(dhsfc*max(hmbl_prev(i-1,j), hmbl_prev(i+1,j), &
662 hmbl_prev(i,j-1), hmbl_prev(i,j+1)), &
663 max(hmbl_prev(i-1,j) - dhd*min(h_sum(i,j),h_sum(i-1,j)), &
664 hmbl_prev(i+1,j) - dhd*min(h_sum(i,j),h_sum(i+1,j)), &
665 hmbl_prev(i,j-1) - dhd*min(h_sum(i,j),h_sum(i,j-1)), &
666 hmbl_prev(i,j+1) - dhd*min(h_sum(i,j),h_sum(i,j+1))) )
667
668 hsfc_min(i,j) = max(h(i,0), min(hsfc(i), h_nbr))
669
670 if (cs%limit_det) max_bl_det(i) = max(0.0, hsfc(i)-h_nbr)
671 enddo
672 endif
673
674 if (cs%id_Hsfc_max > 0) then ; do i=is,ie
675 hsfc_max(i,j) = hsfc(i)
676 enddo ; endif
677 endif
678
679 ! Move water left in the former mixed layer into the buffer layer and
680 ! from the buffer layer into the interior. These steps might best be
681 ! treated in conjunction.
682 if (cs%nkbl == 1) then
683 call mixedlayer_detrain_1(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), spv0(:,0:), rcv(:,0:), &
684 gv%Rlay(:), dt, dt__diag, d_ea, d_eb, j, g, gv, us, cs, &
685 drcv_dt, drcv_ds, max_bl_det)
686 elseif (cs%nkbl == 2) then
687 call mixedlayer_detrain_2(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), spv0(:,0:), rcv(:,0:), &
688 gv%Rlay(:), dt, dt__diag, d_ea, j, g, gv, us, cs, &
689 dr0_dt, dr0_ds, dspv0_dt, dspv0_ds, drcv_dt, drcv_ds, max_bl_det)
690 else ! CS%nkbl not = 1 or 2
691 ! This code only works with 1 or 2 buffer layers.
692 call mom_error(fatal, "MOM_mixed_layer: CS%nkbl must be 1 or 2 for now.")
693 endif
694
695 if (cs%id_Hsfc_used > 0) then
696 do i=is,ie ; hsfc_used(i,j) = h(i,0) ; enddo
697 do k=cs%nkml+1,nkmb ; do i=is,ie
698 hsfc_used(i,j) = hsfc_used(i,j) + h(i,k)
699 enddo ; enddo
700 endif
701
702! Now set the properties of the layers in the mixed layer in the original
703! 3-d variables.
704 if (cs%Resolve_Ekman .and. (cs%nkml>1)) then
705 ! The thickness of the topmost piece of the mixed layer is given by
706 ! h_1 = H / (3 + sqrt(|f|*H^2/2*nu_max)), which asymptotes to the Ekman
707 ! layer depth and 1/3 of the mixed layer depth. This curve has been
708 ! determined to maximize the impact of the Ekman transport in the mixed
709 ! layer TKE budget with nkml=2. With nkml=3, this should also be used,
710 ! as the third piece will then optimally describe mixed layer
711 ! restratification. For nkml>=4 the whole strategy should be revisited.
712 do i=is,ie
713 ! Perhaps in the following, u* could be replaced with u*+w*?
714 ku_star = cs%vonKar * u_star_h_2d(i,j)
715 if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
716 if (fluxes%frac_shelf_h(i,j) > 0.0) then
717 if (allocated(tv%SpV_avg)) then
718 ku_star = (1.0 - fluxes%frac_shelf_h(i,j)) * ku_star + &
719 fluxes%frac_shelf_h(i,j) * ((cs%vonKar*fluxes%ustar_shelf(i,j)) / &
720 (gv%H_to_RZ * tv%SpV_avg(i,j,1)))
721 else
722 ku_star = (1.0 - fluxes%frac_shelf_h(i,j)) * ku_star + &
723 fluxes%frac_shelf_h(i,j) * (cs%vonKar*gv%Z_to_H*fluxes%ustar_shelf(i,j))
724 endif
725 endif
726 endif
727 absf_x_h = 0.25 * h(i,0) * &
728 ((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
729 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
730 ! If the mixed layer vertical viscosity specification is changed in
731 ! MOM_vert_friction.F90, this line will have to be modified accordingly.
732 h_3d(i,j,1) = h(i,0) / (3.0 + sqrt(absf_x_h*(absf_x_h + 2.0*ku_star) / ku_star**2))
733 do k=2,cs%nkml
734 ! The other layers are evenly distributed through the mixed layer.
735 h_3d(i,j,k) = (h(i,0)-h_3d(i,j,1)) * inkmlm1
736 d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
737 d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
738 enddo
739 enddo
740 else
741 do i=is,ie
742 h_3d(i,j,1) = h(i,0) * inkml
743 enddo
744 do k=2,cs%nkml ; do i=is,ie
745 h_3d(i,j,k) = h(i,0) * inkml
746 d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
747 d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
748 enddo ; enddo
749 endif
750 do i=is,ie ; h(i,0) = 0.0 ; enddo
751 do k=1,cs%nkml ; do i=is,ie
752 tv%T(i,j,k) = t(i,0) ; tv%S(i,j,k) = s(i,0)
753 enddo ; enddo
754
755 ! These sum needs to be done in the original layer space.
756
757 ! The treatment of layer 1 is atypical because evaporation shows up as
758 ! negative ea(i,1), and because all precipitation goes straight into layer 1.
759 ! The code is ordered so that any roundoff errors in ea are lost the surface.
760! do i=is,ie ; eaml(i,1) = 0.0 ; enddo
761! do k=2,nz ; do i=is,ie ; eaml(i,k) = eaml(i,k-1) - d_ea(i,k-1) ; enddo ; enddo
762! do i=is,ie ; eaml(i,1) = netMassInOut(i) ; enddo
763
764
765 do i=is,ie
766! eaml(i,nz) is derived from h(i,nz) - h_orig(i,nz) = eaml(i,nz) - ebml(i,nz-1)
767 ebml(i,nz) = 0.0
768 eaml(i,nz) = (h(i,nz) - h_orig(i,nz)) - d_eb(i,nz)
769 enddo
770 do k=nz-1,1,-1 ; do i=is,ie
771 ebml(i,k) = ebml(i,k+1) - d_eb(i,k+1)
772 eaml(i,k) = eaml(i,k+1) + d_ea(i,k)
773 enddo ; enddo
774 do i=is,ie ; eaml(i,1) = netmassinout(i) ; enddo
775
776 ! Copy the interior thicknesses and other fields back to the 3-d arrays.
777 do k=cs%nkml+1,nz ; do i=is,ie
778 h_3d(i,j,k) = h(i,k) ; tv%T(i,j,k) = t(i,k) ; tv%S(i,j,k) = s(i,k)
779 enddo ; enddo
780
781 do k=1,nz ; do i=is,ie
782 ea(i,j,k) = ea(i,j,k) + eaml(i,k)
783 eb(i,j,k) = eb(i,j,k) + ebml(i,k)
784 enddo ; enddo
785
786 if (cs%id_h_mismatch > 0) then
787 do i=is,ie
788 h_miss(i,j) = abs(h_3d(i,j,1) - (h_orig(i,1) + &
789 (eaml(i,1) + (ebml(i,1) - eaml(i,1+1)))))
790 enddo
791 do k=2,nz-1 ; do i=is,ie
792 h_miss(i,j) = h_miss(i,j) + abs(h_3d(i,j,k) - (h_orig(i,k) + &
793 ((eaml(i,k) - ebml(i,k-1)) + (ebml(i,k) - eaml(i,k+1)))))
794 enddo ; enddo
795 do i=is,ie
796 h_miss(i,j) = h_miss(i,j) + abs(h_3d(i,j,nz) - (h_orig(i,nz) + &
797 ((eaml(i,nz) - ebml(i,nz-1)) + ebml(i,nz))))
798 enddo
799 endif
800
801 enddo ! j loop
802 !$OMP end parallel
803
804 ! Whenever thickness changes let the diag manager know, target grids
805 ! for vertical remapping may need to be regenerated.
806 ! This needs to happen after the H update and before the next post_data.
807 call diag_update_remap_grids(cs%diag)
808
809
810 if (write_diags) then
811 if (cs%id_ML_depth > 0) &
812 call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
813 if (cs%id_TKE_wind > 0) &
814 call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
815 if (cs%id_TKE_RiBulk > 0) &
816 call post_data(cs%id_TKE_RiBulk, cs%diag_TKE_RiBulk, cs%diag)
817 if (cs%id_TKE_conv > 0) &
818 call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
819 if (cs%id_TKE_pen_SW > 0) &
820 call post_data(cs%id_TKE_pen_SW, cs%diag_TKE_pen_SW, cs%diag)
821 if (cs%id_TKE_mixing > 0) &
822 call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
823 if (cs%id_TKE_mech_decay > 0) &
824 call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
825 if (cs%id_TKE_conv_decay > 0) &
826 call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
827 if (cs%id_TKE_conv_s2 > 0) &
828 call post_data(cs%id_TKE_conv_s2, cs%diag_TKE_conv_s2, cs%diag)
829 if (cs%id_PE_detrain > 0) &
830 call post_data(cs%id_PE_detrain, cs%diag_PE_detrain, cs%diag)
831 if (cs%id_PE_detrain2 > 0) &
832 call post_data(cs%id_PE_detrain2, cs%diag_PE_detrain2, cs%diag)
833 if (cs%id_h_mismatch > 0) &
834 call post_data(cs%id_h_mismatch, h_miss, cs%diag)
835 if (cs%id_Hsfc_used > 0) &
836 call post_data(cs%id_Hsfc_used, hsfc_used, cs%diag)
837 if (cs%id_Hsfc_max > 0) &
838 call post_data(cs%id_Hsfc_max, hsfc_max, cs%diag)
839 if (cs%id_Hsfc_min > 0) &
840 call post_data(cs%id_Hsfc_min, hsfc_min, cs%diag)
841 endif
842
843end subroutine bulkmixedlayer
844
845!> This subroutine does instantaneous convective entrainment into the buffer
846!! layers and mixed layers to remove hydrostatic instabilities. Any water that
847!! is lighter than currently in the mixed- or buffer- layer is entrained.
848subroutine convective_adjustment(h, u, v, R0, SpV0, Rcv, T, S, eps, d_eb, &
849 dKE_CA, cTKE, j, G, GV, US, CS, nz_conv)
850 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
851 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
852 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
853 !! The units of h are referred to as H below.
854 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h
855 !! points [L T-1 ~> m s-1].
856 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h
857 !! points [L T-1 ~> m s-1].
858 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
859 !! surface pressure [R ~> kg m-3].
860 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: SpV0 !< Specific volume referenced to
861 !! surface pressure [R-1 ~> m3 kg-1].
862 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential
863 !! density [R ~> kg m-3].
864 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Layer temperatures [C ~> degC].
865 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Layer salinities [S ~> ppt].
866 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The negligibly small amount of water
867 !! that will be left in each layer [H ~> m or kg m-2].
868 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
869 !! in the entrainment from below [H ~> m or kg m-2].
870 !! Positive values go with mass gain by
871 !! a layer.
872 real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dKE_CA !< The vertically integrated change in
873 !! kinetic energy due to convective
874 !! adjustment [H Z2 T-2 ~> m3 s-2 or J m-2].
875 real, dimension(SZI_(G),SZK_(GV)), intent(out) :: cTKE !< The buoyant turbulent kinetic energy
876 !! source due to convective adjustment
877 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
878 integer, intent(in) :: j !< The j-index to work on.
879 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
880 type(bulkmixedlayer_cs), intent(in) :: CS !< Bulk mixed layer control structure
881 integer, optional, intent(in) :: nz_conv !< If present, the number of layers
882 !! over which to do convective adjustment
883 !! (perhaps CS%nkml).
884
885 ! Local variables
886 real, dimension(SZI_(G)) :: &
887 R0_tot, & ! The integrated potential density referenced to the surface
888 ! of the layers which are fully entrained [H R ~> kg m-2 or kg2 m-5].
889 spv0_tot, & ! The integrated specific volume referenced to the surface
890 ! of the layers which are fully entrained [H R-1 ~> m4 kg-1 or m].
891 rcv_tot, & ! The integrated coordinate value potential density of the
892 ! layers that are fully entrained [H R ~> kg m-2 or kg2 m-5].
893 ttot, & ! The integrated temperature of layers which are fully
894 ! entrained [C H ~> degC m or degC kg m-2].
895 stot, & ! The integrated salt of layers which are fully entrained
896 ! [H S ~> m ppt or ppt kg m-2].
897 uhtot, & ! The depth integrated zonal velocities in the mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1]
898 vhtot, & ! The depth integrated meridional velocities in the mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1]
899 ke_orig, & ! The total mean kinetic energy per unit area in the mixed layer before
900 ! convection, [H L2 T-2 ~> m3 s-2 or kg s-2].
901 h_orig_k1 ! The depth of layer k1 before convective adjustment [H ~> m or kg m-2].
902 real :: h_ent ! The thickness from a layer that is entrained [H ~> m or kg m-2].
903 real :: Ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
904 real :: g_H_2Rho0 ! Half the gravitational acceleration times
905 ! the conversion from H to Z divided by the mean density,
906 ! in [Z2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
907 logical :: unstable
908 integer :: is, ie, nz, i, k, k1, nzc, nkmb
909
910 is = g%isc ; ie = g%iec ; nz = gv%ke
911 g_h_2rho0 = (gv%g_Earth_Z_T2 * gv%H_to_Z) / (2.0 * gv%Rho0)
912 nzc = nz ; if (present(nz_conv)) nzc = nz_conv
913 nkmb = cs%nkml+cs%nkbl
914
915! Undergo instantaneous entrainment into the buffer layers and mixed layers
916! to remove hydrostatic instabilities. Any water that is lighter than currently
917! in the layer is entrained.
918 do k1=min(nzc-1,nkmb),1,-1
919 do i=is,ie
920 h_orig_k1(i) = h(i,k1)
921 ke_orig(i) = 0.5*h(i,k1)*((u(i,k1)**2) + (v(i,k1)**2))
922 uhtot(i) = h(i,k1)*u(i,k1) ; vhtot(i) = h(i,k1)*v(i,k1)
923 if (cs%nonBous_energetics) then
924 spv0_tot(i) = spv0(i,k1) * h(i,k1)
925 else
926 r0_tot(i) = r0(i,k1) * h(i,k1)
927 endif
928 ctke(i,k1) = 0.0 ; dke_ca(i,k1) = 0.0
929
930 rcv_tot(i) = rcv(i,k1) * h(i,k1)
931 ttot(i) = t(i,k1) * h(i,k1) ; stot(i) = s(i,k1) * h(i,k1)
932 enddo
933 do k=k1+1,nzc
934 do i=is,ie
935 if (cs%nonBous_energetics) then
936 unstable = (spv0_tot(i) < h(i,k1)*spv0(i,k))
937 else
938 unstable = (r0_tot(i) > h(i,k1)*r0(i,k))
939 endif
940 if ((h(i,k) > eps(i,k)) .and. unstable) then
941 h_ent = h(i,k)-eps(i,k)
942 if (cs%nonBous_energetics) then
943 ! This and the other energy calculations assume that specific volume is
944 ! conserved during mixing, which ignores certain thermobaric contributions.
945 ctke(i,k1) = ctke(i,k1) + 0.5 * h_ent * (gv%g_Earth_Z_T2 * gv%H_to_RZ) * &
946 (h(i,k1)*spv0(i,k) - spv0_tot(i)) * cs%nstar2
947 spv0_tot(i) = spv0_tot(i) + h_ent * spv0(i,k)
948 else
949 ctke(i,k1) = ctke(i,k1) + h_ent * g_h_2rho0 * &
950 (r0_tot(i) - h(i,k1)*r0(i,k)) * cs%nstar2
951 r0_tot(i) = r0_tot(i) + h_ent * r0(i,k)
952 endif
953 if (k < nkmb) then
954 ctke(i,k1) = ctke(i,k1) + ctke(i,k)
955 dke_ca(i,k1) = dke_ca(i,k1) + dke_ca(i,k)
956 endif
957 ke_orig(i) = ke_orig(i) + 0.5*h_ent* &
958 ((u(i,k)*u(i,k)) + (v(i,k)*v(i,k)))
959 uhtot(i) = uhtot(i) + h_ent*u(i,k)
960 vhtot(i) = vhtot(i) + h_ent*v(i,k)
961
962 rcv_tot(i) = rcv_tot(i) + h_ent * rcv(i,k)
963 ttot(i) = ttot(i) + h_ent * t(i,k)
964 stot(i) = stot(i) + h_ent * s(i,k)
965 h(i,k1) = h(i,k1) + h_ent ; h(i,k) = eps(i,k)
966
967 d_eb(i,k) = d_eb(i,k) - h_ent
968 d_eb(i,k1) = d_eb(i,k1) + h_ent
969 endif
970 enddo
971 enddo
972! Determine the temperature, salinity, and velocities of the mixed or buffer
973! layer in question, if it has entrained.
974 do i=is,ie ; if (h(i,k1) > h_orig_k1(i)) then
975 ih = 1.0 / h(i,k1)
976 if (cs%nonBous_energetics) then
977 spv0(i,k1) = spv0_tot(i) * ih
978 else
979 r0(i,k1) = r0_tot(i) * ih
980 endif
981 u(i,k1) = uhtot(i) * ih ; v(i,k1) = vhtot(i) * ih
982 dke_ca(i,k1) = dke_ca(i,k1) + cs%bulk_Ri_convective * &
983 (ke_orig(i) - 0.5*h(i,k1)*((u(i,k1)**2) + (v(i,k1)**2)))
984 rcv(i,k1) = rcv_tot(i) * ih
985 t(i,k1) = ttot(i) * ih ; s(i,k1) = stot(i) * ih
986 endif ; enddo
987 enddo
988! If lower mixed or buffer layers are massless, give them the properties of the
989! layer above.
990 do k=2,min(nzc,nkmb) ; do i=is,ie ; if (h(i,k) == 0.0) then
991 if (cs%nonBous_energetics) then
992 spv0(i,k) = spv0(i,k-1)
993 else
994 r0(i,k) = r0(i,k-1)
995 endif
996 rcv(i,k) = rcv(i,k-1) ; t(i,k) = t(i,k-1) ; s(i,k) = s(i,k-1)
997 endif ; enddo ; enddo
998
999end subroutine convective_adjustment
1000
1001!> This subroutine causes the mixed layer to entrain to the depth of free
1002!! convection. The depth of free convection is the shallowest depth at which the
1003!! fluid is denser than the average of the fluid above.
1004subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
1005 R0_tot, SpV0_tot, Rcv_tot, u, v, T, S, R0, SpV0, Rcv, eps, &
1006 dR0_dT, dSpV0_dT, dRcv_dT, dR0_dS, dSpV0_dS, dRcv_dS, &
1007 netMassInOut, netMassOut, Net_heat, Net_salt, &
1008 nsw, Pen_SW_bnd, opacity_band, Conv_En, &
1009 dKE_FC, j, ksort, G, GV, US, CS, tv, fluxes, dt, &
1010 aggregate_FW_forcing)
1011 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1012 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1013 real, dimension(SZI_(G),SZK0_(GV)), &
1014 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
1015 !! The units of h are referred to as H below.
1016 real, dimension(SZI_(G),SZK_(GV)), &
1017 intent(inout) :: d_eb !< The downward increase across a layer in the
1018 !! layer in the entrainment from below [H ~> m or kg m-2].
1019 !! Positive values go with mass gain by a layer.
1020 real, dimension(SZI_(G)), intent(out) :: htot !< The accumulated mixed layer thickness [H ~> m or kg m-2].
1021 real, dimension(SZI_(G)), intent(out) :: Ttot !< The depth integrated mixed layer temperature
1022 !! [C H ~> degC m or degC kg m-2].
1023 real, dimension(SZI_(G)), intent(out) :: Stot !< The depth integrated mixed layer salinity
1024 !! [S H ~> ppt m or ppt kg m-2].
1025 real, dimension(SZI_(G)), intent(out) :: uhtot !< The depth integrated mixed layer zonal
1026 !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1027 real, dimension(SZI_(G)), intent(out) :: vhtot !< The integrated mixed layer meridional
1028 !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1029 real, dimension(SZI_(G)), intent(out) :: R0_tot !< The integrated mixed layer potential density referenced
1030 !! to 0 pressure [H R ~> kg m-2 or kg2 m-5].
1031 real, dimension(SZI_(G)), intent(out) :: SpV0_tot !< The integrated mixed layer specific volume referenced
1032 !! to 0 pressure [H R-1 ~> m4 kg-1 or m].
1033 real, dimension(SZI_(G)), intent(out) :: Rcv_tot !< The integrated mixed layer coordinate
1034 !! variable potential density [H R ~> kg m-2 or kg2 m-5].
1035 real, dimension(SZI_(G),SZK_(GV)), &
1036 intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
1037 real, dimension(SZI_(G),SZK_(GV)), &
1038 intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
1039 real, dimension(SZI_(G),SZK0_(GV)), &
1040 intent(in) :: T !< Layer temperatures [C ~> degC].
1041 real, dimension(SZI_(G),SZK0_(GV)), &
1042 intent(in) :: S !< Layer salinities [S ~> ppt].
1043 real, dimension(SZI_(G),SZK0_(GV)), &
1044 intent(in) :: R0 !< Potential density referenced to
1045 !! surface pressure [R ~> kg m-3].
1046 real, dimension(SZI_(G),SZK0_(GV)), &
1047 intent(in) :: SpV0 !< Specific volume referenced to
1048 !! surface pressure [R-1 ~> m3 kg-1].
1049 real, dimension(SZI_(G),SZK0_(GV)), &
1050 intent(in) :: Rcv !< The coordinate defining potential
1051 !! density [R ~> kg m-3].
1052 real, dimension(SZI_(G),SZK_(GV)), &
1053 intent(in) :: eps !< The negligibly small amount of water
1054 !! that will be left in each layer [H ~> m or kg m-2].
1055 real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of R0 with respect to
1056 !! temperature [R C-1 ~> kg m-3 degC-1].
1057 real, dimension(SZI_(G)), intent(in) :: dSpV0_dT !< The partial derivative of SpV0 with respect to
1058 !! temperature [R-1 C-1 ~> m3 kg-1 degC-1].
1059 real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of Rcv with respect to
1060 !! temperature [R C-1 ~> kg m-3 degC-1].
1061 real, dimension(SZI_(G)), intent(in) :: dR0_dS !< The partial derivative of R0 with respect to
1062 !! salinity [R S-1 ~> kg m-3 ppt-1].
1063 real, dimension(SZI_(G)), intent(in) :: dSpV0_dS !< The partial derivative of SpV0 with respect to
1064 !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
1065 real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of Rcv with respect to
1066 !! salinity [R S-1 ~> kg m-3 ppt-1].
1067 real, dimension(SZI_(G)), intent(in) :: netMassInOut !< The net mass flux (if non-Boussinesq)
1068 !! or volume flux (if Boussinesq) into the ocean
1069 !! within a time step [H ~> m or kg m-2]. (I.e. P+R-E.)
1070 real, dimension(SZI_(G)), intent(in) :: netMassOut !< The mass or volume flux out of the ocean
1071 !! within a time step [H ~> m or kg m-2].
1072 real, dimension(SZI_(G)), intent(in) :: Net_heat !< The net heating at the surface over a time
1073 !! step [C H ~> degC m or degC kg m-2]. Any penetrating
1074 !! shortwave radiation is not included in Net_heat.
1075 real, dimension(SZI_(G)), intent(in) :: Net_salt !< The net surface salt flux into the ocean
1076 !! over a time step [S H ~> ppt m or ppt kg m-2].
1077 integer, intent(in) :: nsw !< The number of bands of penetrating
1078 !! shortwave radiation.
1079 real, dimension(max(nsw,1),SZI_(G)), intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
1080 !! heating at the sea surface in each penetrating
1081 !! band [C H ~> degC m or degC kg m-2].
1082 real, dimension(max(nsw,1),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< The opacity in each band of
1083 !! penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
1084 real, dimension(SZI_(G)), intent(out) :: Conv_En !< The buoyant turbulent kinetic energy source
1085 !! due to free convection [H Z2 T-2 ~> m3 s-2 or J m-2].
1086 real, dimension(SZI_(G)), intent(out) :: dKE_FC !< The vertically integrated change in kinetic
1087 !! energy due to free convection [H Z2 T-2 ~> m3 s-2 or J m-2].
1088 integer, intent(in) :: j !< The j-index to work on.
1089 integer, dimension(SZI_(G),SZK_(GV)), &
1090 intent(in) :: ksort !< The density-sorted k-indices.
1091 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1092 type(bulkmixedlayer_cs), intent(in) :: CS !< Bulk mixed layer control structure
1093 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
1094 !! available thermodynamic fields. Absent
1095 !! fields have NULL pointers.
1096 type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
1097 !! possible forcing fields. Unused fields
1098 !! have NULL pointers.
1099 real, intent(in) :: dt !< Time increment [T ~> s].
1100 logical, intent(in) :: aggregate_FW_forcing !< If true, the net incoming and
1101 !! outgoing surface freshwater fluxes are
1102 !! combined before being applied, instead of
1103 !! being applied separately.
1104
1105! This subroutine causes the mixed layer to entrain to the depth of free
1106! convection. The depth of free convection is the shallowest depth at which the
1107! fluid is denser than the average of the fluid above.
1108
1109 ! Local variables
1110 real, dimension(SZI_(G)) :: &
1111 massOutRem, & ! Evaporation that remains to be supplied [H ~> m or kg m-2].
1112 netMassIn ! mass entering through ocean surface [H ~> m or kg m-2]
1113 real :: SW_trans ! The fraction of shortwave radiation
1114 ! that is not absorbed in a layer [nondim].
1115 real :: Pen_absorbed ! The amount of penetrative shortwave radiation
1116 ! that is absorbed in a layer [C H ~> degC m or degC kg m-2].
1117 real :: h_avail ! The thickness in a layer available for
1118 ! entrainment [H ~> m or kg m-2].
1119 real :: h_ent ! The thickness from a layer that is entrained [H ~> m or kg m-2].
1120 real :: T_precip ! The temperature of the precipitation [C ~> degC].
1121 real :: C1_3, C1_6 ! 1/3 and 1/6 [nondim]
1122 real :: En_fn, Frac, x1 ! Nondimensional temporary variables [nondim].
1123 real :: dr, dr0 ! Temporary variables [R H ~> kg m-2 or kg2 m-5] or [H R-1 ~> m4 kg-1 or m].
1124 real :: dr_ent, dr_comp ! Temporary variables [R H ~> kg m-2 or kg2 m-5].
1125 real :: dr_dh ! The partial derivative of dr_ent with h_ent [R ~> kg m-3].
1126 real :: h_min, h_max ! The minimum and maximum estimates for h_ent [H ~> m or kg m-2]
1127 real :: h_prev ! The previous estimate for h_ent [H ~> m or kg m-2]
1128 real :: h_evap ! The thickness that is evaporated [H ~> m or kg m-2].
1129 real :: dh_Newt ! The Newton's method estimate of the change in
1130 ! h_ent between iterations [H ~> m or kg m-2].
1131 real :: g_H_2Rho0 ! Half the gravitational acceleration times
1132 ! the conversion from H to Z divided by the mean density,
1133 ! [Z2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
1134 real :: Angstrom ! The minimum layer thickness [H ~> m or kg m-2].
1135 real :: opacity ! The opacity converted to inverse thickness units [H-1 ~> m-1 or m2 kg-1]
1136 real :: sum_Pen_En ! The potential energy change due to penetrating
1137 ! shortwave radiation, integrated over a layer
1138 ! [H R ~> kg m-2 or kg2 m-5].
1139 real :: Idt ! 1.0/dt [T-1 ~> s-1]
1140 integer :: is, ie, nz, i, k, ks, itt, n
1141 real, dimension(max(nsw,1)) :: &
1142 C2, & ! Temporary variable [R H-1 ~> kg m-4 or m-1].
1143 r_SW_top ! Temporary variables [H R ~> kg m-2 or kg2 m-5].
1144
1145 angstrom = gv%Angstrom_H
1146 c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
1147 g_h_2rho0 = (gv%g_Earth_Z_T2 * gv%H_to_Z) / (2.0 * gv%Rho0)
1148 idt = 1.0 / dt
1149 is = g%isc ; ie = g%iec ; nz = gv%ke
1150
1151 do i=is,ie ; if (ksort(i,1) > 0) then
1152 k = ksort(i,1)
1153
1154 if (aggregate_fw_forcing) then
1155 massoutrem(i) = 0.0
1156 if (netmassinout(i) < 0.0) massoutrem(i) = -netmassinout(i)
1157 netmassin(i) = netmassinout(i) + massoutrem(i)
1158 else
1159 massoutrem(i) = -netmassout(i)
1160 netmassin(i) = netmassinout(i) - netmassout(i)
1161 endif
1162
1163 ! htot is an Angstrom (taken from layer 1) plus any net precipitation.
1164 h_ent = max(min(angstrom,h(i,k)-eps(i,k)),0.0)
1165 htot(i) = h_ent + netmassin(i)
1166 h(i,k) = h(i,k) - h_ent
1167 d_eb(i,k) = d_eb(i,k) - h_ent
1168
1169 pen_absorbed = 0.0
1170 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1171 sw_trans = exp(-htot(i)*opacity_band(n,i,k))
1172 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0-sw_trans)
1173 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1174 endif ; enddo
1175
1176 ! Precipitation is assumed to have the same temperature and velocity
1177 ! as layer 1. Because layer 1 might not be the topmost layer, this
1178 ! involves multiple terms.
1179 t_precip = t(i,1)
1180 ttot(i) = (net_heat(i) + (netmassin(i) * t_precip + h_ent * t(i,k))) + &
1181 pen_absorbed
1182 ! Net_heat contains both heat fluxes and the heat content of mass fluxes.
1183 !! Ttot(i) = netMassIn(i) * T_precip + h_ent * T(i,k)
1184 !! Ttot(i) = Net_heat(i) + Ttot(i)
1185 !! Ttot(i) = Ttot(i) + Pen_absorbed
1186 ! smg:
1187 ! Ttot(i) = (Net_heat(i) + (h_ent * T(i,k))) + Pen_absorbed
1188 stot(i) = h_ent*s(i,k) + net_salt(i)
1189 uhtot(i) = u(i,1)*netmassin(i) + u(i,k)*h_ent
1190 vhtot(i) = v(i,1)*netmassin(i) + v(i,k)*h_ent
1191 if (cs%nonBous_energetics) then
1192 spv0_tot(i) = (h_ent*spv0(i,k) + netmassin(i)*spv0(i,1)) + &
1193! dSpV0_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1194 (dspv0_dt(i)*(net_heat(i) + pen_absorbed) - &
1195 dspv0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1196 else
1197 r0_tot(i) = (h_ent*r0(i,k) + netmassin(i)*r0(i,1)) + &
1198! dR0_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1199 (dr0_dt(i)*(net_heat(i) + pen_absorbed) - &
1200 dr0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1201 endif
1202 rcv_tot(i) = (h_ent*rcv(i,k) + netmassin(i)*rcv(i,1)) + &
1203! dRcv_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1204 (drcv_dt(i)*(net_heat(i) + pen_absorbed) - &
1205 drcv_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1206 conv_en(i) = 0.0 ; dke_fc(i) = 0.0
1207 if (associated(fluxes%heat_content_massin)) &
1208 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1209 t_precip * netmassin(i) * gv%H_to_RZ * tv%C_p * idt
1210 if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1211 t_precip * netmassin(i) * gv%H_to_RZ
1212 else ! This is a massless column, but zero out the summed variables anyway for safety.
1213 htot(i) = 0.0 ; ttot(i) = 0.0 ; stot(i) = 0.0 ; rcv_tot = 0.0
1214 r0_tot(i) = 0.0 ; spv0_tot(i) = 0.0
1215 uhtot(i) = 0.0 ; vhtot(i) = 0.0 ; conv_en(i) = 0.0 ; dke_fc(i) = 0.0
1216 endif ; enddo
1217
1218 ! Now do netMassOut case in this block.
1219 ! At this point htot contains an Angstrom of fluid from layer 0 plus netMassIn.
1220 do ks=1,nz
1221 do i=is,ie ; if (ksort(i,ks) > 0) then
1222 k = ksort(i,ks)
1223
1224 if ((htot(i) < angstrom) .and. (h(i,k) > eps(i,k))) then
1225 ! If less than an Angstrom was available from the layers above plus
1226 ! any precipitation, add more fluid from this layer.
1227 h_ent = min(angstrom-htot(i), h(i,k)-eps(i,k))
1228 htot(i) = htot(i) + h_ent
1229 h(i,k) = h(i,k) - h_ent
1230 d_eb(i,k) = d_eb(i,k) - h_ent
1231
1232 if (cs%nonBous_energetics) then
1233 spv0_tot(i) = spv0_tot(i) + h_ent*spv0(i,k)
1234 else
1235 r0_tot(i) = r0_tot(i) + h_ent*r0(i,k)
1236 endif
1237 uhtot(i) = uhtot(i) + h_ent*u(i,k)
1238 vhtot(i) = vhtot(i) + h_ent*v(i,k)
1239
1240 rcv_tot(i) = rcv_tot(i) + h_ent*rcv(i,k)
1241 ttot(i) = ttot(i) + h_ent*t(i,k)
1242 stot(i) = stot(i) + h_ent*s(i,k)
1243 endif
1244
1245 ! Water is removed from the topmost layers with any mass.
1246 ! We may lose layers if they are thin enough.
1247 ! The salt that is left behind goes into Stot.
1248 if ((massoutrem(i) > 0.0) .and. (h(i,k) > eps(i,k))) then
1249 if (massoutrem(i) > (h(i,k) - eps(i,k))) then
1250 h_evap = h(i,k) - eps(i,k)
1251 h(i,k) = eps(i,k)
1252 massoutrem(i) = massoutrem(i) - h_evap
1253 else
1254 h_evap = massoutrem(i)
1255 h(i,k) = h(i,k) - h_evap
1256 massoutrem(i) = 0.0
1257 endif
1258
1259 stot(i) = stot(i) + h_evap*s(i,k)
1260 if (cs%nonBous_energetics) then
1261 spv0_tot(i) = spv0_tot(i) + dspv0_ds(i)*h_evap*s(i,k)
1262 else
1263 r0_tot(i) = r0_tot(i) + dr0_ds(i)*h_evap*s(i,k)
1264 endif
1265 rcv_tot(i) = rcv_tot(i) + drcv_ds(i)*h_evap*s(i,k)
1266 d_eb(i,k) = d_eb(i,k) - h_evap
1267
1268 ! smg: when resolve the A=B code, we will set
1269 ! heat_content_massout = heat_content_massout - T(i,k)*h_evap*GV%H_to_RZ*tv%C_p*Idt
1270 ! by uncommenting the lines here.
1271 ! we will also then completely remove TempXpme from the model.
1272 if (associated(fluxes%heat_content_massout)) &
1273 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) - &
1274 t(i,k)*h_evap*gv%H_to_RZ * tv%C_p * idt
1275 if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
1276 t(i,k)*h_evap*gv%H_to_RZ
1277
1278 endif
1279
1280 ! The following section calculates how much fluid will be entrained.
1281 h_avail = h(i,k) - eps(i,k)
1282 if (h_avail > 0.0) then
1283 h_ent = 0.0
1284
1285 if (cs%nonBous_energetics) then
1286 dr = htot(i)*spv0(i,k) - spv0_tot(i)
1287
1288 dr0 = dr
1289 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1290 dr0 = dr0 + (dspv0_dt(i)*pen_sw_bnd(n,i)) * &
1291 opacity_band(n,i,k)*htot(i)
1292 endif ; enddo
1293 else
1294 dr = r0_tot(i) - htot(i)*r0(i,k)
1295
1296 dr0 = dr
1297 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1298 dr0 = dr0 - (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1299 opacity_band(n,i,k)*htot(i)
1300 endif ; enddo
1301 endif
1302
1303 ! Some entrainment will occur from this layer.
1304 if (dr0 > 0.0) then
1305 dr_comp = dr
1306 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1307 ! Compare the density at the bottom of a layer with the
1308 ! density averaged over the mixed layer and that layer.
1309 opacity = opacity_band(n,i,k)
1310 sw_trans = exp(-h_avail*opacity)
1311 if (cs%nonBous_energetics) then
1312 dr_comp = dr_comp - (dspv0_dt(i)*pen_sw_bnd(n,i)) * &
1313 ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
1314 else
1315 dr_comp = dr_comp + (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1316 ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
1317 endif
1318 endif ; enddo
1319 if (dr_comp >= 0.0) then
1320 ! The entire layer is entrained.
1321 h_ent = h_avail
1322 else
1323 ! The layer is partially entrained. Iterate to determine how much
1324 ! entrainment occurs. Solve for the h_ent at which dr_ent = 0.
1325
1326 ! Instead of assuming that the curve is linear between the two end
1327 ! points, assume that the change is concentrated near small values
1328 ! of entrainment. On average, this saves about 1 iteration.
1329 frac = dr0 / (dr0 - dr_comp)
1330 h_ent = h_avail * frac*frac
1331 h_min = 0.0 ; h_max = h_avail
1332
1333 do n=1,nsw
1334 if (cs%nonBous_energetics) then
1335 r_sw_top(n) = -dspv0_dt(i) * pen_sw_bnd(n,i)
1336 else
1337 r_sw_top(n) = dr0_dt(i) * pen_sw_bnd(n,i)
1338 endif
1339 c2(n) = r_sw_top(n) * opacity_band(n,i,k)**2
1340 enddo
1341 do itt=1,10
1342 dr_ent = dr ; dr_dh = 0.0
1343 do n=1,nsw
1344 opacity = opacity_band(n,i,k)
1345 sw_trans = exp(-h_ent*opacity)
1346 dr_ent = dr_ent + r_sw_top(n) * ((1.0 - sw_trans) - &
1347 opacity*(htot(i)+h_ent)*sw_trans)
1348 dr_dh = dr_dh + c2(n) * (htot(i)+h_ent) * sw_trans
1349 enddo
1350
1351 if (dr_ent > 0.0) then
1352 h_min = h_ent
1353 else
1354 h_max = h_ent
1355 endif
1356
1357 dh_newt = -dr_ent / dr_dh
1358 h_prev = h_ent ; h_ent = h_prev+dh_newt
1359 if (h_ent > h_max) then
1360 h_ent = 0.5*(h_prev+h_max)
1361 elseif (h_ent < h_min) then
1362 h_ent = 0.5*(h_prev+h_min)
1363 endif
1364
1365 if (abs(dh_newt) < 0.2*angstrom) exit
1366 enddo
1367
1368 endif
1369
1370 ! Now that the amount of entrainment (h_ent) has been determined,
1371 ! calculate changes in various terms.
1372 sum_pen_en = 0.0 ; pen_absorbed = 0.0
1373 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1374 opacity = opacity_band(n,i,k)
1375 sw_trans = exp(-h_ent*opacity)
1376
1377 x1 = h_ent*opacity
1378 if (x1 < 2.0e-5) then
1379 en_fn = (opacity*htot(i)*(1.0 - 0.5*(x1 - c1_3*x1)) + &
1380 x1*x1*c1_6)
1381 else
1382 en_fn = ((opacity*htot(i) + 2.0) * &
1383 ((1.0-sw_trans) / x1) - 1.0 + sw_trans)
1384 endif
1385 if (cs%nonBous_energetics) then
1386 sum_pen_en = sum_pen_en + (dspv0_dt(i)*pen_sw_bnd(n,i)) * en_fn
1387 else
1388 sum_pen_en = sum_pen_en - (dr0_dt(i)*pen_sw_bnd(n,i)) * en_fn
1389 endif
1390
1391 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1392 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1393 endif ; enddo
1394
1395 if (cs%nonBous_energetics) then
1396 ! This and the other energy calculations assume that specific volume is
1397 ! conserved during mixing, which ignores certain thermobaric contributions.
1398 conv_en(i) = conv_en(i) + 0.5 * (gv%g_Earth_Z_T2 * gv%H_to_RZ) * h_ent * &
1399 ( (spv0(i,k)*htot(i) - spv0_tot(i)) + sum_pen_en )
1400 spv0_tot(i) = spv0_tot(i) + (h_ent * spv0(i,k) + pen_absorbed*dspv0_dt(i))
1401 else
1402 conv_en(i) = conv_en(i) + g_h_2rho0 * h_ent * &
1403 ( (r0_tot(i) - r0(i,k)*htot(i)) + sum_pen_en )
1404 r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1405 endif
1406
1407 stot(i) = stot(i) + h_ent * s(i,k)
1408 ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1409 rcv_tot(i) = rcv_tot(i) + (h_ent * rcv(i,k) + pen_absorbed*drcv_dt(i))
1410 endif ! dr0 > 0.0
1411
1412
1413 if ((h_ent > 0.0) .and. (htot(i) > 0.0)) &
1414 dke_fc(i) = dke_fc(i) + cs%bulk_Ri_convective * 0.5 * &
1415 ((h_ent) / (htot(i)*(h_ent+htot(i)))) * &
1416 (((uhtot(i)-u(i,k)*htot(i))**2) + ((vhtot(i)-v(i,k)*htot(i))**2))
1417
1418 if (h_ent > 0.0) then
1419 htot(i) = htot(i) + h_ent
1420 h(i,k) = h(i,k) - h_ent
1421 d_eb(i,k) = d_eb(i,k) - h_ent
1422 if (cs%convect_mom_bug) then
1423 uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
1424 else
1425 uhtot(i) = uhtot(i) + h_ent*u(i,k) ; vhtot(i) = vhtot(i) + h_ent*v(i,k)
1426 endif
1427 endif
1428
1429 endif ! h_avail>0
1430 endif ; enddo ! i loop
1431 enddo ! k loop
1432
1433end subroutine mixedlayer_convection
1434
1435!> This subroutine determines the TKE available at the depth of free
1436!! convection to drive mechanical entrainment.
1437subroutine find_starting_tke(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, &
1438 TKE, TKE_river, Idecay_len_TKE, cMKE, tv, dt, Idt_diag, &
1439 j, ksort, G, GV, US, CS)
1440 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1441 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1442 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1443 real, dimension(SZI_(G)), intent(in) :: htot !< The accumulated mixed layer thickness
1444 !! [H ~> m or kg m-2]
1445 real, dimension(SZI_(G)), intent(in) :: h_CA !< The mixed layer depth after convective
1446 !! adjustment [H ~> m or kg m-2].
1447 type(forcing), intent(in) :: fluxes !< A structure containing pointers to any
1448 !! possible forcing fields. Unused fields
1449 !! have NULL pointers.
1450 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: U_star_2d !< The wind friction velocity, calculated
1451 !! using the Boussinesq reference density or
1452 !! the time-evolving surface density in
1453 !! non-Boussinesq mode [Z T-1 ~> m s-1]
1454 real, dimension(SZI_(G)), intent(inout) :: Conv_En !< The buoyant turbulent kinetic energy source
1455 !! due to free convection [H Z2 T-2 ~> m3 s-2 or J m-2].
1456 real, dimension(SZI_(G)), intent(in) :: dKE_FC !< The vertically integrated change in
1457 !! kinetic energy due to free convection
1458 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
1459 real, dimension(SZI_(G),SZK_(GV)), &
1460 intent(in) :: cTKE !< The buoyant turbulent kinetic energy
1461 !! source due to convective adjustment
1462 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
1463 real, dimension(SZI_(G),SZK_(GV)), &
1464 intent(in) :: dKE_CA !< The vertically integrated change in
1465 !! kinetic energy due to convective
1466 !! adjustment [H Z2 T-2 ~> m3 s-2 or J m-2].
1467 real, dimension(SZI_(G)), intent(out) :: TKE !< The turbulent kinetic energy available for
1468 !! mixing over a time step [H Z2 T-2 ~> m3 s-2 or J m-2]
1469 real, dimension(SZI_(G)), intent(out) :: Idecay_len_TKE !< The inverse of the vertical decay
1470 !! scale for TKE [H-1 ~> m-1 or m2 kg-1].
1471 real, dimension(SZI_(G)), intent(in) :: TKE_river !< The source of turbulent kinetic energy
1472 !! available for driving mixing at river mouths
1473 !! [H Z2 T-3 ~> m3 s-3 or W m-2].
1474 real, dimension(2,SZI_(G)), intent(out) :: cMKE !< Coefficients of HpE and HpE^2 in
1475 !! calculating the denominator of MKE_rate,
1476 !! [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
1477 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
1478 !! available thermodynamic fields.
1479 real, intent(in) :: dt !< The time step [T ~> s].
1480 real, intent(in) :: Idt_diag !< The inverse of the accumulated diagnostic
1481 !! time interval [T-1 ~> s-1].
1482 integer, intent(in) :: j !< The j-index to work on.
1483 integer, dimension(SZI_(G),SZK_(GV)), &
1484 intent(in) :: ksort !< The density-sorted k-indices.
1485 type(bulkmixedlayer_cs), intent(inout) :: CS !< Bulk mixed layer control structure
1486
1487! This subroutine determines the TKE available at the depth of free
1488! convection to drive mechanical entrainment.
1489
1490 ! Local variables
1491 real :: dKE_conv ! The change in mean kinetic energy due to all convection [H Z2 T-2 ~> m3 s-2 or J m-2].
1492 real :: nstar_FC ! The effective efficiency with which the energy released by
1493 ! free convection is converted to TKE, often ~0.2 [nondim].
1494 real :: nstar_CA ! The effective efficiency with which the energy released by
1495 ! convective adjustment is converted to TKE, often ~0.2 [nondim].
1496 real :: TKE_CA ! The potential energy released by convective adjustment if
1497 ! that release is positive [H Z2 T-2 ~> m3 s-2 or J m-2].
1498 real :: MKE_rate_CA ! MKE_rate for convective adjustment [nondim], 0 to 1.
1499 real :: MKE_rate_FC ! MKE_rate for free convection [nondim], 0 to 1.
1500 real :: totEn_Z ! The total potential energy released by convection, [H Z2 T-2 ~> m3 s-2 or J m-2].
1501 real :: Ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
1502 real :: exp_kh ! The nondimensional decay of TKE across a layer [nondim].
1503 real :: absf ! The absolute value of f averaged to thickness points [T-1 ~> s-1].
1504 real :: U_star ! The friction velocity [Z T-1 ~> m s-1].
1505 real :: absf_Ustar ! The absolute value of f divided by U_star converted to thickness units [H-1 ~> m-1 or m2 kg-1]
1506 real :: wind_TKE_src ! The surface wind source of TKE [H Z2 T-3 ~> m3 s-3 or W m-2].
1507 real :: diag_wt ! The ratio of the current timestep to the diagnostic
1508 ! timestep (which may include 2 calls) [nondim].
1509 real :: H_to_Z ! The thickness to depth conversion factor, which in non-Boussinesq mode is
1510 ! based on the layer-averaged specific volume [Z H-1 ~> nondim or m3 kg-1]
1511 integer :: is, ie, nz, i
1512
1513 is = g%isc ; ie = g%iec ; nz = gv%ke
1514 diag_wt = dt * idt_diag
1515
1516 if (cs%omega_frac >= 1.0) absf = 2.0*cs%omega
1517 do i=is,ie
1518 u_star = u_star_2d(i,j)
1519
1520 if (gv%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
1521 h_to_z = gv%H_to_Z
1522 else
1523 h_to_z = gv%H_to_RZ * tv%SpV_avg(i,j,1)
1524 endif
1525
1526 if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
1527 if (fluxes%frac_shelf_h(i,j) > 0.0) &
1528 u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
1529 fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
1530 endif
1531
1532 if (u_star < cs%ustar_min) u_star = cs%ustar_min
1533
1534 if (cs%omega_frac < 1.0) then
1535 absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
1536 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
1537 if (cs%omega_frac > 0.0) &
1538 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1539 endif
1540 absf_ustar = h_to_z * absf / u_star
1541 idecay_len_tke(i) = absf_ustar * cs%TKE_decay
1542
1543! The first number in the denominator could be anywhere up to 16. The
1544! value of 3 was chosen to minimize the time-step dependence of the amount
1545! of shear-driven mixing in 10 days of a 1-degree global model, emphasizing
1546! the equatorial areas. Although it is not cast as a parameter, it should
1547! be considered an empirical parameter, and it might depend strongly on the
1548! number of sublayers in the mixed layer and their locations.
1549! This equation assumes that small & large scales contribute to mixed layer
1550! deepening at similar rates, even though small scales are dissipated more
1551! rapidly (implying they are less efficient).
1552! Ih = H_to_Z / (16.0*CS%vonKar*U_star*dt)
1553 ih = h_to_z / (3.0*cs%vonKar*u_star*dt)
1554 cmke(1,i) = 4.0 * ih ; cmke(2,i) = absf_ustar * ih
1555
1556 if (idecay_len_tke(i) > 0.0) then
1557 exp_kh = exp(-htot(i)*idecay_len_tke(i))
1558 else
1559 exp_kh = 1.0
1560 endif
1561
1562! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based
1563! on a curve fit from the data of Wang (GRL, 2003).
1564! Note: Ro = 1.0/sqrt(0.5 * dt * (absf*htot(i))**3 / totEn)
1565 if (conv_en(i) < 0.0) conv_en(i) = 0.0
1566 if (ctke(i,1) > 0.0) then ; tke_ca = ctke(i,1) ; else ; tke_ca = 0.0 ; endif
1567 if ((htot(i) >= h_ca(i)) .or. (tke_ca == 0.0)) then
1568 toten_z = (conv_en(i) + tke_ca)
1569
1570 if (toten_z > 0.0) then
1571 nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1572 sqrt(0.5 * dt * (h_to_z**2*(absf*htot(i))**3) * toten_z))
1573 else
1574 nstar_fc = cs%nstar
1575 endif
1576 nstar_ca = nstar_fc
1577 else
1578 ! This reconstructs the Buoyancy flux within the topmost htot of water.
1579 if (conv_en(i) > 0.0) then
1580 toten_z = (conv_en(i) + tke_ca * (htot(i) / h_ca(i)) )
1581 nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1582 sqrt(0.5 * dt * (h_to_z**2*(absf*htot(i))**3) * toten_z))
1583 else
1584 nstar_fc = cs%nstar
1585 endif
1586
1587 toten_z = (conv_en(i) + tke_ca)
1588 if (tke_ca > 0.0) then
1589 nstar_ca = cs%nstar * toten_z / (toten_z + 0.2 * &
1590 sqrt(0.5 * dt * (h_to_z**2*(absf*h_ca(i))**3) * toten_z))
1591 else
1592 nstar_ca = cs%nstar
1593 endif
1594 endif
1595
1596 if (dke_fc(i) + dke_ca(i,1) > 0.0) then
1597 if (htot(i) >= h_ca(i)) then
1598 mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1599 mke_rate_ca = mke_rate_fc
1600 else
1601 mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1602 mke_rate_ca = 1.0 / (1.0 + h_ca(i)*(cmke(1,i) + cmke(2,i)*h_ca(i)) )
1603 endif
1604 else
1605 ! This branch just saves unnecessary calculations.
1606 mke_rate_fc = 1.0 ; mke_rate_ca = 1.0
1607 endif
1608
1609 dke_conv = dke_ca(i,1) * mke_rate_ca + dke_fc(i) * mke_rate_fc
1610! At this point, it is assumed that cTKE is positive and stored in TKE_CA!
1611! Note: Removed factor of 2 in u*^3 terms.
1612 if (gv%Boussinesq .or. gv%semi_Boussinesq .or. .not.(associated(fluxes%tau_mag))) then
1613 tke(i) = (dt*cs%mstar)*((gv%Z_to_H*(u_star*u_star*u_star))*exp_kh) + &
1614 (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca * tke_ca)
1615 else
1616 ! Note that GV%Z_to_H*U_star**3 = GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star
1617 tke(i) = (dt*cs%mstar) * ((gv%RZ_to_H * fluxes%tau_mag(i,j) * u_star)*exp_kh) + &
1618 (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca * tke_ca)
1619 endif
1620
1621 if (cs%do_rivermix) then ! Add additional TKE at river mouths
1622 tke(i) = tke(i) + tke_river(i)*dt*exp_kh
1623 endif
1624
1625 if (cs%TKE_diagnostics) then
1626 if (gv%Boussinesq .or. gv%semi_Boussinesq .or. .not.(associated(fluxes%tau_mag))) then
1627 wind_tke_src = cs%mstar*(gv%Z_to_H*u_star*u_star*u_star) * diag_wt
1628 else
1629 wind_tke_src = cs%mstar*(gv%RZ_to_H * fluxes%tau_mag(i,j) * u_star) * diag_wt
1630 endif
1631 cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + &
1632 ( wind_tke_src + tke_river(i) * diag_wt )
1633 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + dke_conv*idt_diag
1634 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1635 (exp_kh-1.0)*(wind_tke_src + dke_conv*idt_diag)
1636 cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + &
1637 idt_diag * (nstar_fc*conv_en(i) + nstar_ca*tke_ca)
1638 cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + &
1639 idt_diag * ((cs%nstar-nstar_fc)*conv_en(i) + (cs%nstar-nstar_ca)*tke_ca)
1640 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
1641 idt_diag * (ctke(i,1)-tke_ca)
1642 endif
1643 enddo
1644
1645end subroutine find_starting_tke
1646
1647!> This subroutine calculates mechanically driven entrainment.
1648subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
1649 R0_tot, SpV0_tot, Rcv_tot, u, v, T, S, R0, SpV0, Rcv, eps, &
1650 dR0_dT, dSpV0_dT, dRcv_dT, cMKE, Idt_diag, nsw, &
1651 Pen_SW_bnd, opacity_band, TKE, &
1652 Idecay_len_TKE, j, ksort, G, GV, US, CS)
1653 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1654 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1655 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1656 real, dimension(SZI_(G),SZK0_(GV)), &
1657 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
1658 real, dimension(SZI_(G),SZK_(GV)), &
1659 intent(inout) :: d_eb !< The downward increase across a layer in the
1660 !! layer in the entrainment from below [H ~> m or kg m-2].
1661 !! Positive values go with mass gain by a layer.
1662 real, dimension(SZI_(G)), intent(inout) :: htot !< The accumulated mixed layer thickness [H ~> m or kg m-2].
1663 real, dimension(SZI_(G)), intent(inout) :: Ttot !< The depth integrated mixed layer temperature
1664 !! [C H ~> degC m or degC kg m-2].
1665 real, dimension(SZI_(G)), intent(inout) :: Stot !< The depth integrated mixed layer salinity
1666 !! [S H ~> ppt m or ppt kg m-2].
1667 real, dimension(SZI_(G)), intent(inout) :: uhtot !< The depth integrated mixed layer zonal
1668 !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1669 real, dimension(SZI_(G)), intent(inout) :: vhtot !< The integrated mixed layer meridional
1670 !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1671 real, dimension(SZI_(G)), intent(inout) :: R0_tot !< The integrated mixed layer potential density
1672 !! referenced to 0 pressure [H R ~> kg m-2 or kg2 m-5].
1673 real, dimension(SZI_(G)), intent(inout) :: SpV0_tot !< The integrated mixed layer specific volume referenced
1674 !! to 0 pressure [H R-1 ~> m4 kg-1 or m].
1675 real, dimension(SZI_(G)), intent(inout) :: Rcv_tot !< The integrated mixed layer coordinate variable
1676 !! potential density [H R ~> kg m-2 or kg2 m-5].
1677 real, dimension(SZI_(G),SZK_(GV)), &
1678 intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
1679 real, dimension(SZI_(G),SZK_(GV)), &
1680 intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
1681 real, dimension(SZI_(G),SZK0_(GV)), &
1682 intent(in) :: T !< Layer temperatures [C ~> degC].
1683 real, dimension(SZI_(G),SZK0_(GV)), &
1684 intent(in) :: S !< Layer salinities [S ~> ppt].
1685 real, dimension(SZI_(G),SZK0_(GV)), &
1686 intent(in) :: R0 !< Potential density referenced to
1687 !! surface pressure [R ~> kg m-3].
1688 real, dimension(SZI_(G),SZK0_(GV)), &
1689 intent(in) :: SpV0 !< Specific volume referenced to
1690 !! surface pressure [R-1 ~> m3 kg-1].
1691 real, dimension(SZI_(G),SZK0_(GV)), &
1692 intent(in) :: Rcv !< The coordinate defining potential
1693 !! density [R ~> kg m-3].
1694 real, dimension(SZI_(G),SZK_(GV)), &
1695 intent(in) :: eps !< The negligibly small amount of water
1696 !! that will be left in each layer [H ~> m or kg m-2].
1697 real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of R0 with respect to
1698 !! temperature [R C-1 ~> kg m-3 degC-1].
1699 real, dimension(SZI_(G)), intent(in) :: dSpV0_dT !< The partial derivative of SpV0 with respect to
1700 !! temperature [R-1 C-1 ~> m3 kg-1 degC-1].
1701 real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of Rcv with respect to
1702 !! temperature [R C-1 ~> kg m-3 degC-1].
1703 real, dimension(2,SZI_(G)), intent(in) :: cMKE !< Coefficients of HpE and HpE^2 used in calculating the
1704 !! denominator of MKE_rate; the two elements have differing
1705 !! units of [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
1706 real, intent(in) :: Idt_diag !< The inverse of the accumulated diagnostic
1707 !! time interval [T-1 ~> s-1].
1708 integer, intent(in) :: nsw !< The number of bands of penetrating
1709 !! shortwave radiation.
1710 real, dimension(max(nsw,1),SZI_(G)), intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
1711 !! heating at the sea surface in each penetrating
1712 !! band [C H ~> degC m or degC kg m-2].
1713 real, dimension(max(nsw,1),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< The opacity in each band of
1714 !! penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
1715 real, dimension(SZI_(G)), intent(inout) :: TKE !< The turbulent kinetic energy
1716 !! available for mixing over a time
1717 !! step [H Z2 T-2 ~> m3 s-2 or J m-2].
1718 real, dimension(SZI_(G)), intent(inout) :: Idecay_len_TKE !< The vertical TKE decay rate [H-1 ~> m-1 or m2 kg-1].
1719 integer, intent(in) :: j !< The j-index to work on.
1720 integer, dimension(SZI_(G),SZK_(GV)), &
1721 intent(in) :: ksort !< The density-sorted k-indices.
1722 type(bulkmixedlayer_cs), intent(inout) :: CS !< Bulk mixed layer control structure
1723
1724! This subroutine calculates mechanically driven entrainment.
1725
1726 ! Local variables
1727 real :: SW_trans ! The fraction of shortwave radiation that is not
1728 ! absorbed in a layer [nondim].
1729 real :: Pen_absorbed ! The amount of penetrative shortwave radiation
1730 ! that is absorbed in a layer [C H ~> degC m or degC kg m-2].
1731 real :: h_avail ! The thickness in a layer available for entrainment [H ~> m or kg m-2].
1732 real :: h_ent ! The thickness from a layer that is entrained [H ~> m or kg m-2].
1733 real :: h_min, h_max ! Limits on the solution for h_ent [H ~> m or kg m-2].
1734 real :: dh_Newt ! The Newton's method estimate of the change in
1735 ! h_ent between iterations [H ~> m or kg m-2].
1736 real :: MKE_rate ! The fraction of the energy in resolved shears
1737 ! within the mixed layer that will be eliminated
1738 ! within a timestep [nondim], 0 to 1.
1739 real :: HpE ! The current thickness plus entrainment [H ~> m or kg m-2].
1740 real :: g_H_2Rho0 ! Half the gravitational acceleration times the
1741 ! conversion from H to m divided by the mean density,
1742 ! in [Z2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
1743 real :: TKE_full_ent ! The TKE remaining if a layer is fully entrained
1744 ! [H Z2 T-2 ~> m3 s-2 or J m-2].
1745 real :: dRL ! Work required to mix water from the next layer
1746 ! across the mixed layer [Z2 T-2 ~> m2 s-2].
1747 real :: Pen_En_Contrib ! Penetrating SW contributions to the changes in
1748 ! TKE, divided by layer thickness in m [Z2 T-2 ~> m2 s-2].
1749 real :: Cpen1 ! A temporary variable [Z2 T-2 ~> m2 s-2].
1750 real :: dMKE ! A temporary variable related to the release of mean
1751 ! kinetic energy [H2 Z2 T-2 ~> m4 s-2 or kg2 m-2 s-2]
1752 real :: TKE_ent ! The TKE that remains if h_ent were entrained [H Z2 T-2 ~> m3 s-2 or J m-2]
1753 real :: TKE_ent1 ! The TKE that would remain, without considering the
1754 ! release of mean kinetic energy [H Z2 T-2 ~> m3 s-2 or J m-2]
1755 real :: dTKE_dh ! The partial derivative of TKE with h_ent [Z2 T-2 ~> m2 s-2]
1756 real :: Pen_dTKE_dh_Contrib ! The penetrating shortwave contribution to
1757 ! dTKE_dh [Z2 T-2 ~> m2 s-2].
1758 real :: EF4_val ! The result of EF4() (see later) [H-1 ~> m-1 or m2 kg-1].
1759 real :: h_neglect ! A thickness that is so small it is usually lost
1760 ! in roundoff and can be neglected [H ~> m or kg m-2].
1761 real :: dEF4_dh ! The partial derivative of EF4 with h [H-2 ~> m-2 or m4 kg-2].
1762 real :: Pen_En1 ! A nondimensional temporary variable [nondim].
1763 real :: kh, exp_kh, f1_kh ! Nondimensional temporary variables related to the
1764 ! fractional decay of TKE across a layer [nondim].
1765 real :: x1, e_x1 ! Nondimensional temporary variables related to the relative decay
1766 ! of TKE and SW radiation across a layer [nondim]
1767 real :: f1_x1, f2_x1, f3_x1 ! Exponential-related functions of x1 [nondim].
1768 real :: E_HxHpE ! Entrainment divided by the product of the new and old
1769 ! thicknesses [H-1 ~> m-1 or m2 kg-1].
1770 real :: Hmix_min ! The minimum mixed layer depth [H ~> m or kg m-2].
1771 real :: opacity ! The opacity of a layer in a band of shortwave radiation [H-1 ~> m-1 or m2 kg-1]
1772 real :: C1_3, C1_6, C1_24 ! 1/3, 1/6, and 1/24. [nondim]
1773 integer :: is, ie, nz, i, k, ks, itt, n
1774
1775 c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0 ; c1_24 = 1.0/24.0
1776 g_h_2rho0 = (gv%g_Earth_Z_T2 * gv%H_to_Z) / (2.0 * gv%Rho0)
1777 hmix_min = cs%Hmix_min
1778 h_neglect = gv%H_subroundoff
1779 is = g%isc ; ie = g%iec ; nz = gv%ke
1780
1781 do ks=1,nz
1782
1783 do i=is,ie ; if (ksort(i,ks) > 0) then
1784 k = ksort(i,ks)
1785
1786 h_avail = h(i,k) - eps(i,k)
1787 if ((h_avail > 0.) .and. ((tke(i) > 0.) .or. (htot(i) < hmix_min))) then
1788 if (cs%nonBous_energetics) then
1789 drl = 0.5 * (gv%g_Earth_Z_T2 * gv%H_to_RZ) * (spv0_tot(i) - spv0(i,k)*htot(i))
1790 else
1791 drl = g_h_2rho0 * (r0(i,k)*htot(i) - r0_tot(i) )
1792 endif
1793 dmke = cs%bulk_Ri_ML * 0.5 * &
1794 (((uhtot(i)-u(i,k)*htot(i))**2) + ((vhtot(i)-v(i,k)*htot(i))**2))
1795
1796! Find the TKE that would remain if the entire layer were entrained.
1797 kh = idecay_len_tke(i)*h_avail ; exp_kh = exp(-kh)
1798 if (kh >= 2.0e-5) then ; f1_kh = (1.0-exp_kh) / kh
1799 else ; f1_kh = (1.0 - kh*(0.5 - c1_6*kh)) ; endif
1800
1801 pen_en_contrib = 0.0
1802 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1803 opacity = opacity_band(n,i,k)
1804! Two different forms are used here to make sure that only negative
1805! values are taken into exponentials to avoid excessively large
1806! numbers. They are, of course, mathematically identical.
1807 if (idecay_len_tke(i) > opacity) then
1808 x1 = (idecay_len_tke(i) - opacity) * h_avail
1809 if (x1 >= 2.0e-5) then
1810 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1811 f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1812 else
1813 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1814 f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1815 endif
1816
1817 pen_en1 = exp(-opacity*h_avail) * &
1818 ((1.0+opacity*htot(i))*f1_x1 + opacity*h_avail*f3_x1)
1819 else
1820 x1 = (opacity - idecay_len_tke(i)) * h_avail
1821 if (x1 >= 2.0e-5) then
1822 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1823 f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1824 else
1825 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1826 f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1827 endif
1828
1829 pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1830 opacity*h_avail*f2_x1)
1831 endif
1832 if (cs%nonBous_energetics) then
1833 pen_en_contrib = pen_en_contrib - &
1834 (0.5 * (gv%g_Earth_Z_T2 * gv%H_to_RZ) * dspv0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
1835 else
1836 pen_en_contrib = pen_en_contrib + &
1837 (g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
1838 endif
1839 endif ; enddo
1840
1841 hpe = htot(i)+h_avail
1842 mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1843 ef4_val = ef4(htot(i)+h_neglect,h_avail,idecay_len_tke(i))
1844 tke_full_ent = (exp_kh*tke(i) - h_avail*(drl*f1_kh + pen_en_contrib)) + &
1845 mke_rate*dmke*ef4_val
1846 if ((tke_full_ent >= 0.0) .or. (h_avail+htot(i) <= hmix_min)) then
1847 ! The layer will be fully entrained.
1848 h_ent = h_avail
1849
1850 if (cs%TKE_diagnostics) then
1851 e_hxhpe = h_ent / ((htot(i)+h_neglect)*(htot(i)+h_ent+h_neglect))
1852 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1853 idt_diag * ((exp_kh-1.0)* tke(i) + h_ent*drl*(1.0-f1_kh) + &
1854 mke_rate*dmke*(ef4_val-e_hxhpe))
1855 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - idt_diag*h_ent*drl
1856 cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1857 idt_diag*h_ent*pen_en_contrib
1858 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1859 idt_diag*mke_rate*dmke*e_hxhpe
1860 endif
1861
1862 tke(i) = tke_full_ent
1863
1864 if (tke(i) <= 0.0) tke(i) = cs%mech_TKE_floor
1865 else
1866! The layer is only partially entrained. The amount that will be
1867! entrained is determined iteratively. No further layers will be
1868! entrained.
1869 h_min = 0.0 ; h_max = h_avail
1870 if (tke(i) <= 0.0) then
1871 h_ent = 0.0
1872 else
1873 h_ent = h_avail * tke(i) / (tke(i) - tke_full_ent)
1874
1875 do itt=1,15
1876 ! Evaluate the TKE that would remain if h_ent were entrained.
1877
1878 kh = idecay_len_tke(i)*h_ent ; exp_kh = exp(-kh)
1879 if (kh >= 2.0e-5) then
1880 f1_kh = (1.0-exp_kh) / kh
1881 else
1882 f1_kh = (1.0 - kh*(0.5 - c1_6*kh))
1883 endif
1884
1885
1886 pen_en_contrib = 0.0 ; pen_dtke_dh_contrib = 0.0
1887 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1888 ! Two different forms are used here to make sure that only negative
1889 ! values are taken into exponentials to avoid excessively large
1890 ! numbers. They are, of course, mathematically identical.
1891 opacity = opacity_band(n,i,k)
1892 sw_trans = exp(-h_ent*opacity)
1893 if (idecay_len_tke(i) > opacity) then
1894 x1 = (idecay_len_tke(i) - opacity) * h_ent
1895 if (x1 >= 2.0e-5) then
1896 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1897 f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1898 else
1899 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1900 f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1901 endif
1902 pen_en1 = sw_trans * ((1.0+opacity*htot(i))*f1_x1 + &
1903 opacity*h_ent*f3_x1)
1904 else
1905 x1 = (opacity - idecay_len_tke(i)) * h_ent
1906 if (x1 >= 2.0e-5) then
1907 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1908 f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1909 else
1910 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1911 f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1912 endif
1913
1914 pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1915 opacity*h_ent*f2_x1)
1916 endif
1917 if (cs%nonBous_energetics) then
1918 cpen1 = -0.5 * (gv%g_Earth_Z_T2 * gv%H_to_RZ) * dspv0_dt(i) * pen_sw_bnd(n,i)
1919 else
1920 cpen1 = g_h_2rho0 * dr0_dt(i) * pen_sw_bnd(n,i)
1921 endif
1922 pen_en_contrib = pen_en_contrib + cpen1*(pen_en1 - f1_kh)
1923 pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
1924 cpen1*((1.0-sw_trans) - opacity*(htot(i) + h_ent)*sw_trans)
1925 endif ; enddo ! (Pen_SW_bnd(n,i) > 0.0)
1926
1927 tke_ent1 = exp_kh* tke(i) - h_ent*(drl*f1_kh + pen_en_contrib)
1928 ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i),def4_dh)
1929 hpe = htot(i)+h_ent
1930 mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1931 tke_ent = tke_ent1 + dmke*ef4_val*mke_rate
1932 ! TKE_ent is the TKE that would remain if h_ent were entrained.
1933
1934 dtke_dh = ((-idecay_len_tke(i)*tke_ent1 - drl) + &
1935 pen_dtke_dh_contrib) + dmke * mke_rate* &
1936 (def4_dh - ef4_val*mke_rate*(cmke(1,i)+2.0*cmke(2,i)*hpe))
1937 ! dh_Newt = -TKE_ent / dTKE_dh
1938 ! Bisect if the Newton's method prediction is outside of the bounded range.
1939 if (tke_ent > 0.0) then
1940 if ((h_max-h_ent)*(-dtke_dh) > tke_ent) then
1941 dh_newt = -tke_ent / dtke_dh
1942 else
1943 dh_newt = 0.5*(h_max-h_ent)
1944 endif
1945 h_min = h_ent
1946 else
1947 if ((h_min-h_ent)*(-dtke_dh) < tke_ent) then
1948 dh_newt = -tke_ent / dtke_dh
1949 else
1950 dh_newt = 0.5*(h_min-h_ent)
1951 endif
1952 h_max = h_ent
1953 endif
1954 h_ent = h_ent + dh_newt
1955
1956 if (abs(dh_newt) < 0.2*gv%Angstrom_H) exit
1957 enddo
1958 endif
1959
1960 if (h_ent < hmix_min-htot(i)) h_ent = hmix_min - htot(i)
1961
1962 if (cs%TKE_diagnostics) then
1963 hpe = htot(i)+h_ent
1964 mke_rate = 1.0/(1.0 + cmke(1,i)*hpe + cmke(2,i)*hpe**2)
1965 ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i))
1966
1967 e_hxhpe = h_ent / ((htot(i)+h_neglect)*(hpe+h_neglect))
1968 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1969 idt_diag * ((exp_kh-1.0)* tke(i) + h_ent*drl*(1.0-f1_kh) + &
1970 dmke*mke_rate*(ef4_val-e_hxhpe))
1971 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - idt_diag*h_ent*drl
1972 cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - idt_diag*h_ent*pen_en_contrib
1973 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + idt_diag*dmke*mke_rate*e_hxhpe
1974 endif
1975
1976 tke(i) = 0.0
1977 endif ! TKE_full_ent > 0.0
1978
1979 pen_absorbed = 0.0
1980 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1981 sw_trans = exp(-h_ent*opacity_band(n,i,k))
1982 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1983 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1984 endif ; enddo
1985
1986 htot(i) = htot(i) + h_ent
1987 if (cs%nonBous_energetics) then
1988 spv0_tot(i) = spv0_tot(i) + (h_ent * spv0(i,k) + pen_absorbed*dspv0_dt(i))
1989 else
1990 r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1991 endif
1992 h(i,k) = h(i,k) - h_ent
1993 d_eb(i,k) = d_eb(i,k) - h_ent
1994
1995 stot(i) = stot(i) + h_ent * s(i,k)
1996 ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1997 rcv_tot(i) = rcv_tot(i) + (h_ent*rcv(i,k) + pen_absorbed*drcv_dt(i))
1998
1999 uhtot(i) = uhtot(i) + u(i,k)*h_ent
2000 vhtot(i) = vhtot(i) + v(i,k)*h_ent
2001 endif ! h_avail > 0.0 .and. TKE(i) > 0.0
2002
2003 endif ; enddo ! i loop
2004 enddo ! k loop
2005
2006end subroutine mechanical_entrainment
2007
2008!> This subroutine generates an array of indices that are sorted by layer
2009!! density.
2010subroutine sort_ml(h, R0, SpV0, eps, G, GV, CS, ksort)
2011 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2012 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2013 real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
2014 real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: R0 !< The potential density used to sort
2015 !! the layers [R ~> kg m-3].
2016 real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: SpV0 !< Specific volume referenced to
2017 !! surface pressure [R-1 ~> m3 kg-1]
2018 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The (small) thickness that must
2019 !! remain in each layer [H ~> m or kg m-2].
2020 type(bulkmixedlayer_cs), intent(in) :: CS !< Bulk mixed layer control structure
2021 integer, dimension(SZI_(G),SZK_(GV)), intent(out) :: ksort !< The k-index to use in the sort.
2022
2023 ! Local variables
2024 real :: R0sort(SZI_(G),SZK_(GV)) ! The sorted potential density [R ~> kg m-3]
2025 real :: SpV0sort(SZI_(G),SZK_(GV)) ! The sorted specific volume [R-1 ~> m3 kg-1]
2026 integer :: nsort(SZI_(G)) ! The number of layers left to sort
2027 logical :: done_sorting(SZI_(G))
2028 integer :: i, k, ks, is, ie, nz, nkmb
2029
2030 is = g%isc ; ie = g%iec ; nz = gv%ke
2031 nkmb = cs%nkml+cs%nkbl
2032
2033 ! Come up with a sorted index list of layers with increasing R0.
2034 ! Assume that the layers below nkmb are already stably stratified.
2035 ! Only layers that are thicker than eps are in the list. Extra elements
2036 ! have an index of -1.
2037
2038 ! This is coded using straight insertion, on the assumption that the
2039 ! layers are usually in the right order (or massless) anyway.
2040
2041 do k=1,nz ; do i=is,ie ; ksort(i,k) = -1 ; enddo ; enddo
2042
2043 do i=is,ie ; nsort(i) = 0 ; done_sorting(i) = .false. ; enddo
2044
2045 if (cs%nonBous_energetics) then
2046 do k=1,nz ; do i=is,ie ; if (h(i,k) > eps(i,k)) then
2047 if (done_sorting(i)) then ; ks = nsort(i) ; else
2048 do ks=nsort(i),1,-1
2049 if (spv0(i,k) <= spv0sort(i,ks)) exit
2050 spv0sort(i,ks+1) = spv0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
2051 enddo
2052 if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
2053 endif
2054
2055 ksort(i,ks+1) = k
2056 spv0sort(i,ks+1) = spv0(i,k)
2057 nsort(i) = nsort(i) + 1
2058 endif ; enddo ; enddo
2059 else
2060 do k=1,nz ; do i=is,ie ; if (h(i,k) > eps(i,k)) then
2061 if (done_sorting(i)) then ; ks = nsort(i) ; else
2062 do ks=nsort(i),1,-1
2063 if (r0(i,k) >= r0sort(i,ks)) exit
2064 r0sort(i,ks+1) = r0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
2065 enddo
2066 if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
2067 endif
2068
2069 ksort(i,ks+1) = k
2070 r0sort(i,ks+1) = r0(i,k)
2071 nsort(i) = nsort(i) + 1
2072 endif ; enddo ; enddo
2073 endif
2074
2075end subroutine sort_ml
2076
2077!> This subroutine actually moves properties between layers to achieve a
2078!! resorted state, with all of the resorted water either moved into the correct
2079!! interior layers or in the top nkmb layers.
2080subroutine resort_ml(h, T, S, R0, SpV0, Rcv, RcvTgt, eps, d_ea, d_eb, ksort, G, GV, CS, &
2081 dR0_dT, dR0_dS, dSpV0_dT, dSpV0_dS, dRcv_dT, dRcv_dS)
2082 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2083 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid
2084 !! structure.
2085 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
2086 !! Layer 0 is the new mixed layer.
2087 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Layer temperatures [C ~> degC].
2088 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Layer salinities [S ~> ppt].
2089 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
2090 !! surface pressure [R ~> kg m-3].
2091 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: SpV0 !< Specific volume referenced to
2092 !! surface pressure [R-1 ~> m3 kg-1]
2093 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining
2094 !! potential density [R ~> kg m-3].
2095 real, dimension(SZK_(GV)), intent(in) :: RcvTgt !< The target value of Rcv for each
2096 !! layer [R ~> kg m-3].
2097 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: eps !< The (small) thickness that must
2098 !! remain in each layer [H ~> m or kg m-2].
2099 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a
2100 !! layer in the entrainment from
2101 !! above [H ~> m or kg m-2].
2102 !! Positive d_ea goes with layer
2103 !! thickness increases.
2104 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a
2105 !! layer in the entrainment from
2106 !! below [H ~> m or kg m-2]. Positive values go
2107 !! with mass gain by a layer.
2108 integer, dimension(SZI_(G),SZK_(GV)), intent(in) :: ksort !< The density-sorted k-indices.
2109 type(bulkmixedlayer_cs), intent(in) :: CS !< Bulk mixed layer control structure
2110 real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of
2111 !! potential density referenced
2112 !! to the surface with potential
2113 !! temperature [R C-1 ~> kg m-3 degC-1].
2114 real, dimension(SZI_(G)), intent(in) :: dR0_dS !< The partial derivative of
2115 !! potential density referenced
2116 !! to the surface with salinity,
2117 !! [R S-1 ~> kg m-3 ppt-1].
2118 real, dimension(SZI_(G)), intent(in) :: dSpV0_dT !< The partial derivative of SpV0 with respect
2119 !! to temperature [R-1 C-1 ~> m3 kg-1 degC-1]
2120 real, dimension(SZI_(G)), intent(in) :: dSpV0_dS !< The partial derivative of SpV0 with respect
2121 !! to salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
2122 real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of
2123 !! coordinate defining potential
2124 !! density with potential
2125 !! temperature [R C-1 ~> kg m-3 degC-1].
2126 real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of
2127 !! coordinate defining potential
2128 !! density with salinity,
2129 !! [R S-1 ~> kg m-3 ppt-1].
2130
2131! If there are no massive light layers above the deepest of the mixed- and
2132! buffer layers, do nothing (except perhaps to reshuffle these layers).
2133! If there are nkbl or fewer layers above the deepest mixed- or buffer-
2134! layers, move them (in sorted order) into the buffer layers, even if they
2135! were previously interior layers.
2136! If there are interior layers that are intermediate in density (both in-situ
2137! and the coordinate density (sigma-2)) between the newly forming mixed layer
2138! and a residual buffer- or mixed layer, and the number of massive layers above
2139! the deepest massive buffer or mixed layer is greater than nkbl, then split
2140! those buffer layers into pieces that match the target density of the two
2141! nearest interior layers.
2142! Otherwise, if there are more than nkbl+1 remaining massive layers
2143
2144 ! Local variables
2145 real :: h_move ! The thickness of water being moved between layers [H ~> m or kg m-2]
2146 real :: h_tgt_old ! The previous thickness of the recipient layer [H ~> m or kg m-2]
2147 real :: I_hnew ! The inverse of a new layer thickness [H-1 ~> m-1 or m2 kg-1]
2148 real :: dT_dS_wt2 ! The square of the relative weighting of temperature and salinity changes
2149 ! when extrapolating to match a target density [C2 S-2 ~> degC2 ppt-2]
2150 real :: dT_dR ! The ratio of temperature changes to density changes when
2151 ! extrapolating [C R-1 ~> degC m3 kg-1]
2152 real :: dS_dR ! The ratio of salinity changes to density changes when
2153 ! extrapolating [S R-1 ~> ppt m3 kg-1]
2154 real :: I_denom ! A work variable with units of [S2 R-2 ~> ppt2 m6 kg-2].
2155 real :: Rcv_int ! The target coordinate density of an interior layer [R ~> kg m-3]
2156 real :: T_up, T_dn ! Temperatures projected to match the target densities of two layers [C ~> degC]
2157 real :: S_up, S_dn ! Salinities projected to match the target densities of two layers [S ~> ppt]
2158 real :: R0_up, R0_dn ! Potential densities projected to match the target coordinate
2159 ! densities of two layers [R ~> kg m-3]
2160 real :: SpV0_up, SpV0_dn ! Specific volumes projected to be consistent with the target coordinate
2161 ! densities of two layers [R-1 ~> m3 kg-1]
2162 real :: I_hup, I_hdn ! Inverse of the new thicknesses of the two layers [H-1 ~> m-1 or m2 kg-1]
2163 real :: h_to_up, h_to_dn ! Thickness transferred to two layers [H ~> m or kg m-2]
2164 real :: wt_dn ! Fraction of the thickness transferred to the deeper layer [nondim]
2165 real :: dR1, dR2 ! Density difference with the target densities of two layers [R ~> kg m-3]
2166 real :: dPE, min_dPE ! Values proportional to the potential energy change due to the merging of a
2167 ! pair of layers [R H2 ~> kg m-1 or kg3 m-7] or [R-1 H2 ~> m5 kg-1 or kg m-1]
2168 real :: hmin, min_hmin ! The thickness of the thinnest layer [H ~> m or kg m-2]
2169 real :: h_tmp(SZK_(GV)) ! A copy of the original layer thicknesses [H ~> m or kg m-2]
2170 real :: R0_tmp(SZK_(GV)) ! A copy of the original layer potential densities [R ~> kg m-3]
2171 real :: SpV0_tmp(SZK_(GV)) ! A copy of the original layer specific volumes [R ~> kg m-3]
2172 real :: T_tmp(SZK_(GV)) ! A copy of the original layer temperatures [C ~> degC]
2173 real :: S_tmp(SZK_(GV)) ! A copy of the original layer salinities [S ~> ppt]
2174 real :: Rcv_tmp(SZK_(GV)) ! A copy of the original layer coordinate densities [R ~> kg m-3]
2175 integer :: ks_min
2176 logical :: sorted, leave_in_layer
2177 integer :: ks_deep(SZI_(G)), k_count(SZI_(G)), ks2_reverse(SZI_(G), SZK_(GV))
2178 integer :: ks2(SZK_(GV))
2179 integer :: i, k, ks, is, ie, nz, k1, k2, k_tgt, k_src, k_int_top
2180 integer :: nks, nkmb, num_interior, top_interior_ks
2181
2182 is = g%isc ; ie = g%iec ; nz = gv%ke
2183 nkmb = cs%nkml+cs%nkbl
2184
2185 dt_ds_wt2 = cs%dT_dS_wt**2
2186
2187! Find out how many massive layers are above the deepest buffer or mixed layer.
2188 do i=is,ie ; ks_deep(i) = -1 ; k_count(i) = 0 ; enddo
2189 do ks=nz,1,-1 ; do i=is,ie ; if (ksort(i,ks) > 0) then
2190 k = ksort(i,ks)
2191
2192 if (h(i,k) > eps(i,k)) then
2193 if (ks_deep(i) == -1) then
2194 if (k <= nkmb) then
2195 ks_deep(i) = ks ; k_count(i) = k_count(i) + 1
2196 ks2_reverse(i,k_count(i)) = k
2197 endif
2198 else
2199 k_count(i) = k_count(i) + 1
2200 ks2_reverse(i,k_count(i)) = k
2201 endif
2202 endif
2203 endif ; enddo ; enddo
2204
2205 do i=is,ie ; if (k_count(i) > 1) then
2206 ! This column might need to be reshuffled.
2207 nks = k_count(i)
2208
2209 ! Put ks2 in the right order and check whether reshuffling is needed.
2210 sorted = .true.
2211 ks2(nks) = ks2_reverse(i,1)
2212 do ks=nks-1,1,-1
2213 ks2(ks) = ks2_reverse(i,1+nks-ks)
2214 if (ks2(ks) > ks2(ks+1)) sorted = .false.
2215 enddo
2216
2217 ! Go to the next column of no reshuffling is needed.
2218 if (sorted) cycle
2219
2220 ! Find out how many interior layers are being reshuffled. If none,
2221 ! then this is a simple swapping procedure.
2222 num_interior = 0 ; top_interior_ks = 0
2223 do ks=nks,1,-1 ; if (ks2(ks) > nkmb) then
2224 num_interior = num_interior+1 ; top_interior_ks = ks
2225 endif ; enddo
2226
2227 if (num_interior >= 1) then
2228 ! Find the lightest interior layer with a target coordinate density
2229 ! greater than the newly forming mixed layer.
2230 do k=nkmb+1,nz ; if (rcv(i,0) < rcvtgt(k)) exit ; enddo
2231 k_int_top = k ; rcv_int = rcvtgt(k)
2232
2233 i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
2234 dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
2235 ds_dr = drcv_ds(i) * i_denom
2236
2237
2238 ! Examine whether layers can be split out of existence. Stop when there
2239 ! is a layer that cannot be handled this way, or when the topmost formerly
2240 ! interior layer has been dealt with.
2241 do ks = nks,top_interior_ks,-1
2242 k = ks2(ks)
2243 leave_in_layer = .false.
2244 if ((k > nkmb) .and. (rcv(i,k) <= rcvtgt(k))) then
2245 if (rcvtgt(k)-rcv(i,k) < cs%BL_split_rho_tol*(rcvtgt(k) - rcvtgt(k-1))) &
2246 leave_in_layer = .true.
2247 elseif (k > nkmb) then
2248 if (rcv(i,k)-rcvtgt(k) < cs%BL_split_rho_tol*(rcvtgt(k+1) - rcvtgt(k))) &
2249 leave_in_layer = .true.
2250 endif
2251
2252 if (leave_in_layer) then
2253 ! Just drop this layer from the sorted list.
2254 nks = nks-1
2255 elseif (rcv(i,k) < rcv_int) then
2256 ! There is no interior layer with a target density that is intermediate
2257 ! between this layer and the mixed layer.
2258 exit
2259 else
2260 ! Try splitting the layer between two interior isopycnal layers.
2261 ! Find the target densities that bracket this layer.
2262 do k2=k_int_top+1,nz ; if (rcv(i,k) < rcvtgt(k2)) exit ; enddo
2263 if (k2>nz) exit
2264
2265 ! This layer is bracketed in density between layers k2-1 and k2.
2266
2267 dr1 = (rcvtgt(k2-1) - rcv(i,k)) ; dr2 = (rcvtgt(k2) - rcv(i,k))
2268 t_up = t(i,k) + dt_dr * dr1
2269 s_up = s(i,k) + ds_dr * dr1
2270 t_dn = t(i,k) + dt_dr * dr2
2271 s_dn = s(i,k) + ds_dr * dr2
2272
2273 if (cs%nonBous_energetics) then
2274 spv0_up = spv0(i,k) + (dt_dr*dspv0_dt(i) + ds_dr*dspv0_ds(i)) * dr1
2275 spv0_dn = spv0(i,k) + (dt_dr*dspv0_dt(i) + ds_dr*dspv0_ds(i)) * dr2
2276
2277 ! Make sure the new properties are acceptable, and avoid creating obviously unstable profiles.
2278 if ((spv0_up < spv0(i,0)) .or. (spv0_dn < spv0(i,0))) exit
2279 else
2280 r0_up = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr1
2281 r0_dn = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr2
2282
2283 ! Make sure the new properties are acceptable, and avoid creating obviously unstable profiles.
2284 if ((r0_up > r0(i,0)) .or. (r0_dn > r0(i,0))) exit
2285 endif
2286
2287 wt_dn = (rcv(i,k) - rcvtgt(k2-1)) / (rcvtgt(k2) - rcvtgt(k2-1))
2288 h_to_up = (h(i,k)-eps(i,k)) * (1.0 - wt_dn)
2289 h_to_dn = (h(i,k)-eps(i,k)) * wt_dn
2290
2291 i_hup = 1.0 / (h(i,k2-1) + h_to_up)
2292 i_hdn = 1.0 / (h(i,k2) + h_to_dn)
2293 if (cs%nonBous_energetics) then
2294 spv0(i,k2-1) = (spv0(i,k2)*h(i,k2-1) + spv0_up*h_to_up) * i_hup
2295 spv0(i,k2) = (spv0(i,k2)*h(i,k2) + spv0_dn*h_to_dn) * i_hdn
2296 else
2297 r0(i,k2-1) = (r0(i,k2)*h(i,k2-1) + r0_up*h_to_up) * i_hup
2298 r0(i,k2) = (r0(i,k2)*h(i,k2) + r0_dn*h_to_dn) * i_hdn
2299 endif
2300
2301 t(i,k2-1) = (t(i,k2)*h(i,k2-1) + t_up*h_to_up) * i_hup
2302 t(i,k2) = (t(i,k2)*h(i,k2) + t_dn*h_to_dn) * i_hdn
2303 s(i,k2-1) = (s(i,k2)*h(i,k2-1) + s_up*h_to_up) * i_hup
2304 s(i,k2) = (s(i,k2)*h(i,k2) + s_dn*h_to_dn) * i_hdn
2305 rcv(i,k2-1) = (rcv(i,k2)*h(i,k2-1) + rcvtgt(k2-1)*h_to_up) * i_hup
2306 rcv(i,k2) = (rcv(i,k2)*h(i,k2) + rcvtgt(k2)*h_to_dn) * i_hdn
2307
2308 h(i,k) = eps(i,k)
2309 h(i,k2) = h(i,k2) + h_to_dn
2310 h(i,k2-1) = h(i,k2-1) + h_to_up
2311
2312 if (k > k2-1) then
2313 d_eb(i,k) = d_eb(i,k) - h_to_up
2314 d_eb(i,k2-1) = d_eb(i,k2-1) + h_to_up
2315 elseif (k < k2-1) then
2316 d_ea(i,k) = d_ea(i,k) - h_to_up
2317 d_ea(i,k2-1) = d_ea(i,k2-1) + h_to_up
2318 endif
2319 if (k > k2) then
2320 d_eb(i,k) = d_eb(i,k) - h_to_dn
2321 d_eb(i,k2) = d_eb(i,k2) + h_to_dn
2322 elseif (k < k2) then
2323 d_ea(i,k) = d_ea(i,k) - h_to_dn
2324 d_ea(i,k2) = d_ea(i,k2) + h_to_dn
2325 endif
2326 nks = nks-1
2327 endif
2328 enddo
2329 endif
2330
2331 do while (nks > nkmb)
2332 ! Having already tried to move surface layers into the interior, there
2333 ! are still too many layers, and layers must be merged until nks=nkmb.
2334 ! Examine every merger of a layer with its neighbors, and merge the ones
2335 ! that increase the potential energy the least. If there are layers
2336 ! with (apparently?) negative potential energy change, choose the one
2337 ! with the smallest total thickness. Repeat until nkmb layers remain.
2338 ! Choose the smaller value for the remaining index for convenience.
2339
2340 ks_min = -1 ; min_dpe = 1.0 ; min_hmin = 0.0
2341 do ks=1,nks-1
2342 k1 = ks2(ks) ; k2 = ks2(ks+1)
2343 if (cs%nonBous_energetics) then
2344 dpe = max(0.0, (spv0(i,k1) - spv0(i,k2)) * (h(i,k1) * h(i,k2)))
2345 else
2346 dpe = max(0.0, (r0(i,k2) - r0(i,k1)) * h(i,k1) * h(i,k2))
2347 endif
2348 hmin = min(h(i,k1)-eps(i,k1), h(i,k2)-eps(i,k2))
2349 if ((ks_min < 0) .or. (dpe < min_dpe) .or. &
2350 ((dpe <= 0.0) .and. (hmin < min_hmin))) then
2351 ks_min = ks ; min_dpe = dpe ; min_hmin = hmin
2352 endif
2353 enddo
2354
2355 ! Now merge the two layers that do the least damage.
2356 k1 = ks2(ks_min) ; k2 = ks2(ks_min+1)
2357 if (k1 < k2) then ; k_tgt = k1 ; k_src = k2
2358 else ; k_tgt = k2 ; k_src = k1 ; ks2(ks_min) = k2 ; endif
2359
2360 h_tgt_old = h(i,k_tgt)
2361 h_move = h(i,k_src)-eps(i,k_src)
2362 h(i,k_src) = eps(i,k_src)
2363 h(i,k_tgt) = h(i,k_tgt) + h_move
2364 i_hnew = 1.0 / (h(i,k_tgt))
2365 if (cs%nonBous_energetics) then
2366 spv0(i,k_tgt) = (spv0(i,k_tgt)*h_tgt_old + spv0(i,k_src)*h_move) * i_hnew
2367 else
2368 r0(i,k_tgt) = (r0(i,k_tgt)*h_tgt_old + r0(i,k_src)*h_move) * i_hnew
2369 endif
2370
2371 t(i,k_tgt) = (t(i,k_tgt)*h_tgt_old + t(i,k_src)*h_move) * i_hnew
2372 s(i,k_tgt) = (s(i,k_tgt)*h_tgt_old + s(i,k_src)*h_move) * i_hnew
2373 rcv(i,k_tgt) = (rcv(i,k_tgt)*h_tgt_old + rcv(i,k_src)*h_move) * i_hnew
2374
2375 d_eb(i,k_src) = d_eb(i,k_src) - h_move
2376 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2377
2378 ! Remove the newly missing layer from the sorted list.
2379 do ks=ks_min+1,nks ; ks2(ks) = ks2(ks+1) ; enddo
2380 nks = nks-1
2381 enddo
2382
2383 ! Check again whether the layers are sorted, and go on to the next column
2384 ! if they are.
2385 sorted = .true.
2386 do ks=1,nks-1 ; if (ks2(ks) > ks2(ks+1)) sorted = .false. ; enddo
2387 if (sorted) cycle
2388
2389 if (nks > 1) then
2390 ! At this point, actually swap the properties of the layers, and place
2391 ! the remaining layers in order starting with nkmb.
2392
2393 ! Save all the properties of the nkmb layers that might be replaced.
2394 do k=1,nkmb
2395 h_tmp(k) = h(i,k)
2396 if (cs%nonBous_energetics) then
2397 spv0_tmp(k) = spv0(i,k)
2398 else
2399 r0_tmp(k) = r0(i,k)
2400 endif
2401 t_tmp(k) = t(i,k) ; s_tmp(k) = s(i,k) ; rcv_tmp(k) = rcv(i,k)
2402
2403 h(i,k) = 0.0
2404 enddo
2405
2406 do ks=nks,1,-1
2407 k_tgt = nkmb - nks + ks ; k_src = ks2(ks)
2408 if (k_tgt == k_src) then
2409 h(i,k_tgt) = h_tmp(k_tgt) ! This layer doesn't move, so put the water back.
2410 cycle
2411 endif
2412
2413 ! Note below that eps=0 for k<=nkmb.
2414 if (k_src > nkmb) then
2415 h_move = h(i,k_src)-eps(i,k_src)
2416 h(i,k_src) = eps(i,k_src)
2417 h(i,k_tgt) = h_move
2418 if (cs%nonBous_energetics) then
2419 spv0(i,k_tgt) = spv0(i,k_src)
2420 else
2421 r0(i,k_tgt) = r0(i,k_src)
2422 endif
2423
2424 t(i,k_tgt) = t(i,k_src) ; s(i,k_tgt) = s(i,k_src)
2425 rcv(i,k_tgt) = rcv(i,k_src)
2426
2427 d_eb(i,k_src) = d_eb(i,k_src) - h_move
2428 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2429 else
2430 h(i,k_tgt) = h_tmp(k_src)
2431 if (cs%nonBous_energetics) then
2432 spv0(i,k_tgt) = spv0_tmp(k_src)
2433 else
2434 r0(i,k_tgt) = r0_tmp(k_src)
2435 endif
2436
2437 t(i,k_tgt) = t_tmp(k_src) ; s(i,k_tgt) = s_tmp(k_src)
2438 rcv(i,k_tgt) = rcv_tmp(k_src)
2439
2440 if (k_src > k_tgt) then
2441 d_eb(i,k_src) = d_eb(i,k_src) - h_tmp(k_src)
2442 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_tmp(k_src)
2443 else
2444 d_ea(i,k_src) = d_ea(i,k_src) - h_tmp(k_src)
2445 d_ea(i,k_tgt) = d_ea(i,k_tgt) + h_tmp(k_src)
2446 endif
2447 endif
2448 enddo
2449 endif
2450
2451 endif ; enddo
2452
2453end subroutine resort_ml
2454
2455!> This subroutine moves any water left in the former mixed layers into the
2456!! two buffer layers and may also move buffer layer water into the interior
2457!! isopycnal layers.
2458subroutine mixedlayer_detrain_2(h, T, S, R0, Spv0, Rcv, RcvTgt, dt, dt_diag, d_ea, j, G, GV, US, CS, &
2459 dR0_dT, dR0_dS, dSpV0_dT, dSpV0_dS, dRcv_dT, dRcv_dS, max_BL_det)
2460 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2461 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2462 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
2463 !! Layer 0 is the new mixed layer.
2464 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Potential temperature [C ~> degC].
2465 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Salinity [S ~> ppt].
2466 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
2467 !! surface pressure [R ~> kg m-3].
2468 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: SpV0 !< Specific volume referenced to
2469 !! surface pressure [R-1 ~> m3 kg-1]
2470 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential
2471 !! density [R ~> kg m-3].
2472 real, dimension(SZK_(GV)), intent(in) :: RcvTgt !< The target value of Rcv for each
2473 !! layer [R ~> kg m-3].
2474 real, intent(in) :: dt !< Time increment [T ~> s].
2475 real, intent(in) :: dt_diag !< The diagnostic time step [T ~> s].
2476 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a layer in
2477 !! the entrainment from above
2478 !! [H ~> m or kg m-2]. Positive d_ea
2479 !! goes with layer thickness increases.
2480 integer, intent(in) :: j !< The meridional row to work on.
2481 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2482 type(bulkmixedlayer_cs), intent(inout) :: CS !< Bulk mixed layer control structure
2483 real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of
2484 !! potential density referenced to the
2485 !! surface with potential temperature,
2486 !! [R C-1 ~> kg m-3 degC-1].
2487 real, dimension(SZI_(G)), intent(in) :: dR0_dS !< The partial derivative of
2488 !! potential density referenced to the
2489 !! surface with salinity
2490 !! [R S-1 ~> kg m-3 ppt-1].
2491 real, dimension(SZI_(G)), intent(in) :: dSpV0_dT !< The partial derivative of specific
2492 !! volume with respect to temeprature
2493 !! [R-1 C-1 ~> m3 kg-1 degC-1]
2494 real, dimension(SZI_(G)), intent(in) :: dSpV0_dS !< The partial derivative of specific
2495 !! volume with respect to salinity
2496 !! [R-1 S-1 ~> m3 kg-1 ppt-1]
2497 real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of
2498 !! coordinate defining potential density
2499 !! with potential temperature,
2500 !! [R C-1 ~> kg m-3 degC-1].
2501 real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of
2502 !! coordinate defining potential density
2503 !! with salinity [R S-1 ~> kg m-3 ppt-1].
2504 real, dimension(SZI_(G)), intent(in) :: max_BL_det !< If non-negative, the maximum
2505 !! detrainment permitted from the buffer
2506 !! layers [H ~> m or kg m-2].
2507
2508! This subroutine moves any water left in the former mixed layers into the
2509! two buffer layers and may also move buffer layer water into the interior
2510! isopycnal layers.
2511
2512 ! Local variables
2513 real :: h_to_bl ! The total thickness detrained to the buffer
2514 ! layers [H ~> m or kg m-2].
2515 real :: R0_to_bl ! The depth integrated amount of R0 that is detrained to the
2516 ! buffer layer [H R ~> kg m-2 or kg2 m-5]
2517 real :: SpV0_to_bl ! The depth integrated amount of SpV0 that is detrained to the
2518 ! buffer layer [H R-1 ~> m4 kg-1 or m]
2519 real :: Rcv_to_bl ! The depth integrated amount of Rcv that is detrained to the
2520 ! buffer layer [H R ~> kg m-2 or kg2 m-5]
2521 real :: T_to_bl ! The depth integrated amount of T that is detrained to the
2522 ! buffer layer [C H ~> degC m or degC kg m-2]
2523 real :: S_to_bl ! The depth integrated amount of S that is detrained to the
2524 ! buffer layer [S H ~> ppt m or ppt kg m-2]
2525 real :: h_min_bl ! The minimum buffer layer thickness [H ~> m or kg m-2].
2526
2527 real :: h1, h2 ! Scalar variables holding the values of
2528 ! h(i,CS%nkml+1) and h(i,CS%nkml+2) [H ~> m or kg m-2].
2529 real :: h1_avail ! The thickness of the upper buffer layer
2530 ! available to move into the lower buffer
2531 ! layer [H ~> m or kg m-2].
2532 real :: stays ! stays is the thickness of the upper buffer
2533 ! layer that remains there [H ~> m or kg m-2].
2534 real :: stays_min, stays_max ! The minimum and maximum permitted values of
2535 ! stays [H ~> m or kg m-2].
2536
2537 logical :: intermediate ! True if the water in layer kb1 is intermediate in density
2538 ! between the water in kb2 and the water being detrained.
2539 logical :: mergeable_bl ! If true, it is an option to combine the two
2540 ! buffer layers and create water that matches
2541 ! the target density of an interior layer.
2542 logical :: better_to_merge ! True if it is energetically favorable to merge layers
2543 real :: stays_merge ! If the two buffer layers can be combined
2544 ! stays_merge is the thickness of the upper
2545 ! layer that remains [H ~> m or kg m-2].
2546 real :: stays_min_merge ! The minimum allowed value of stays_merge [H ~> m or kg m-2].
2547
2548 real :: dR0_2dz, dRcv_2dz ! Half the vertical gradients of R0 and Rcv [R H-1 ~> kg m-4 or m-1]
2549 real :: dSpV0_2dz ! Half the vertical gradients of SpV0 and Rcv [R-1 H-1 ~> m2 kg-1 or m5 kg-2]
2550! real :: dT_2dz ! Half the vertical gradient of T [C H-1 ~> degC m-1 or degC m2 kg-1]
2551! real :: dS_2dz ! Half the vertical gradient of S [S H-1 ~> ppt m-1 or ppt m2 kg-1]
2552 real :: scale_slope ! A nondimensional number < 1 used to scale down
2553 ! the slope within the upper buffer layer when
2554 ! water MUST be detrained to the lower layer [nondim].
2555
2556 real :: dPE_extrap_rhoG ! The potential energy change due to dispersive
2557 ! advection or mixing layers, divided by
2558 ! rho_0*g [H2 ~> m2 or kg2 m-4].
2559 real :: dPE_extrapolate ! The potential energy change due to dispersive advection or
2560 ! mixing layers [R Z3 T-2 ~> J m-2].
2561 real :: dPE_det, dPE_merge ! The energy required to mix the detrained water
2562 ! into the buffer layer or the merge the two
2563 ! buffer layers [R H2 Z T-2 ~> J m-2 or J kg2 m-8].
2564 real :: dPE_det_nB, dPE_merge_nB ! The energy required to mix the detrained water
2565 ! into the buffer layer or the merge the two
2566 ! buffer layers [R Z3 T-2 ~> J m-2].
2567
2568 real :: h_from_ml ! The amount of additional water that must be
2569 ! drawn from the mixed layer [H ~> m or kg m-2].
2570 real :: h_det_h2 ! The amount of detrained water and mixed layer
2571 ! water that will go directly into the lower
2572 ! buffer layer [H ~> m or kg m-2].
2573
2574 real :: h_det_to_h2, h_ml_to_h2 ! The fluxes of detrained and mixed layer water to
2575 ! the lower buffer layer [H ~> m or kg m-2].
2576 real :: h_det_to_h1, h_ml_to_h1 ! The fluxes of detrained and mixed layer water to
2577 ! the upper buffer layer [H ~> m or kg m-2].
2578 real :: h1_to_h2, h1_to_k0 ! The fluxes of upper buffer layer water to the lower buffer layer
2579 ! and to an interior layer that is just denser than the lower
2580 ! buffer layer [H ~> m or kg m-2].
2581 real :: h2_to_k1, h2_to_k1_rem ! Fluxes of lower buffer layer water to the interior layer that
2582 ! is just denser than the lower buffer layer [H ~> m or kg m-2].
2583
2584 real :: R0_det ! Detrained value of potential density referenced to the surface [R ~> kg m-3]
2585 real :: SpV0_det ! Detrained value of specific volume referenced to the surface [R-1 ~> m3 kg-1]
2586 real :: T_det, S_det ! Detrained values of temperature [C ~> degC] and salinity [S ~> ppt]
2587 real :: Rcv_stays, R0_stays ! Values of Rcv and R0 that stay in a layer [R ~> kg m-3]
2588 real :: SpV0_stays ! Values of SpV0 that stay in a layer [R-1 ~> m3 kg-1]
2589 real :: T_stays, S_stays ! Values of T and S that stay in a layer, [C ~> degC] and S [S ~> ppt]
2590 real :: dSpice_det, dSpice_stays! The spiciness difference between an original
2591 ! buffer layer and the water that moves into
2592 ! an interior layer or that stays in that
2593 ! layer [R ~> kg m-3].
2594 real :: dSpice_lim, dSpice_lim2 ! Limits to the spiciness difference between
2595 ! the lower buffer layer and the water that
2596 ! moves into an interior layer [R ~> kg m-3].
2597 real :: dSpice_2dz ! The vertical gradient of spiciness used for
2598 ! advection [R H-1 ~> kg m-4 or m-1].
2599 real :: dSpiceSpV_stays ! The specific volume based spiciness difference between an original
2600 ! buffer layer and the water that stays in that layer [R-1 ~> m3 kg-1]
2601 real :: dSpiceSpV_lim ! A limit on the specific volume based spiciness difference
2602 ! between the lower buffer layer and the water that
2603 ! moves into an interior layer [R-1 ~> m3 kg-1]
2604 real :: dPE_ratio ! Multiplier of dPE_det at which merging is
2605 ! permitted - here (detrainment_per_day/dt)*30
2606 ! days? [nondim]
2607 real :: num_events ! The number of detrainment events over which
2608 ! to prefer merging the buffer layers [nondim].
2609 real :: dPE_time_ratio ! Larger of 1 and the detrainment timescale over dt [nondim].
2610 real :: dT_dS_gauge, dS_dT_gauge ! The relative scales of temperature and
2611 ! salinity changes in defining spiciness, in
2612 ! [C S-1 ~> degC ppt-1] and [S C-1 ~> ppt degC-1].
2613 real :: I_denom ! A work variable with units of [S2 R-2 ~> ppt2 m6 kg-2] or [R2 S2 ~> ppt2 kg2 m-6].
2614
2615 real :: g_2 ! 1/2 g_Earth [Z T-2 ~> m s-2].
2616 real :: Rho0xG ! Rho0 times G_Earth [R Z T-2 ~> kg m-2 s-2].
2617 real :: I2Rho0 ! 1 / (2 Rho0) [R-1 ~> m3 kg-1].
2618 real :: Idt_diag ! The inverse of the timestep used for diagnostics [T-1 ~> s-1].
2619 real :: Idt_H2 ! The square of the conversion from thickness to Z
2620 ! divided by the time step [Z2 H-2 T-1 ~> s-1 or m6 kg-2 s-1].
2621 logical :: stable_Rcv ! If true, the buffer layers are stable with
2622 ! respect to the coordinate potential density.
2623 real :: h_neglect ! A thickness that is so small it is usually lost
2624 ! in roundoff and can be neglected [H ~> m or kg m-2].
2625
2626 real :: s1en ! A work variable [R Z3 T-3 ~> W m-2]
2627 real :: s1, s2, bh0 ! Work variables [H ~> m or kg m-2].
2628 real :: s3sq ! A work variable [H2 ~> m2 or kg2 m-4].
2629 real :: I_ya, b1 ! Nondimensional work variables [nondim]
2630 real :: Ih, Ihdet, Ih1f, Ih2f ! Assorted inverse thickness work variables [H-1 ~> m-1 or m2 kg-1]
2631 real :: Ihk0, Ihk1, Ih12 ! Assorted inverse thickness work variables [H-1 ~> m-1 or m2 kg-1]
2632 real :: dR1, dR2, dR2b, dRk1 ! Assorted density difference work variables [R ~> kg m-3]
2633 real :: dR0, dR21, dRcv ! Assorted density difference work variables [R ~> kg m-3]
2634 real :: dSpV0, dSpVk1 ! Assorted specific volume difference work variables [R-1 ~> m3 kg-1]
2635 real :: dRcv_stays, dRcv_det, dRcv_lim ! Assorted densities [R ~> kg m-3]
2636 real :: Angstrom ! The minimum layer thickness [H ~> m or kg m-2].
2637
2638 real :: h2_to_k1_lim ! A limit on the thickness that can be detrained to layer k1 [H ~> m or kg m-2]
2639 real :: T_new, T_max, T_min ! Temperature of the detrained water and limits on it [C ~> degC]
2640 real :: S_new, S_max, S_min ! Salinity of the detrained water and limits on it [S ~> ppt]
2641 logical :: stable
2642 integer :: i, k, k0, k1, is, ie, nz, kb1, kb2, nkmb
2643
2644 is = g%isc ; ie = g%iec ; nz = gv%ke
2645 kb1 = cs%nkml+1 ; kb2 = cs%nkml+2
2646 nkmb = cs%nkml+cs%nkbl
2647 h_neglect = gv%H_subroundoff
2648 g_2 = 0.5 * gv%g_Earth_Z_T2
2649 rho0xg = gv%Rho0 * gv%g_Earth_Z_T2
2650 idt_diag = 1.0 / dt_diag
2651 idt_h2 = gv%H_to_Z**2 / dt_diag
2652 i2rho0 = 0.5 / gv%Rho0
2653 angstrom = gv%Angstrom_H
2654
2655 ! This is hard coding of arbitrary and dimensional numbers.
2656 dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 / dt_ds_gauge
2657 num_events = 10.0
2658
2659 if (cs%nkbl /= 2) call mom_error(fatal, "MOM_mixed_layer: "// &
2660 "CS%nkbl must be 2 in mixedlayer_detrain_2.")
2661
2662 if (dt < cs%BL_detrain_time) then ; dpe_time_ratio = cs%BL_detrain_time / (dt)
2663 else ; dpe_time_ratio = 1.0 ; endif
2664
2665 do i=is,ie
2666
2667 ! Determine all of the properties being detrained from the mixed layer.
2668
2669 ! As coded this has the k and i loop orders switched, but k is CS%nkml is
2670 ! often just 1 or 2, so this seems like it should not be a problem, especially
2671 ! since it means that a number of variables can now be scalars, not arrays.
2672 h_to_bl = 0.0 ; r0_to_bl = 0.0 ; spv0_to_bl = 0.0
2673 rcv_to_bl = 0.0 ; t_to_bl = 0.0 ; s_to_bl = 0.0
2674
2675 do k=1,cs%nkml ; if (h(i,k) > 0.0) then
2676 h_to_bl = h_to_bl + h(i,k)
2677 if (cs%nonBous_energetics) then
2678 spv0_to_bl = spv0_to_bl + spv0(i,k)*h(i,k)
2679 else
2680 r0_to_bl = r0_to_bl + r0(i,k)*h(i,k)
2681 endif
2682
2683 rcv_to_bl = rcv_to_bl + rcv(i,k)*h(i,k)
2684 t_to_bl = t_to_bl + t(i,k)*h(i,k)
2685 s_to_bl = s_to_bl + s(i,k)*h(i,k)
2686
2687 d_ea(i,k) = d_ea(i,k) - h(i,k)
2688 h(i,k) = 0.0
2689 endif ; enddo
2690
2691 if (cs%nonBous_energetics) then
2692 if (h_to_bl > 0.0) then ; spv0_det = spv0_to_bl / h_to_bl
2693 else ; spv0_det = spv0(i,0) ; endif
2694 else
2695 if (h_to_bl > 0.0) then ; r0_det = r0_to_bl / h_to_bl
2696 else ; r0_det = r0(i,0) ; endif
2697 endif
2698
2699 ! This code does both downward detrainment from both the mixed layer and the
2700 ! buffer layers.
2701 ! Several considerations apply in detraining water into the interior:
2702 ! (1) Water only moves into the interior from the deeper buffer layer,
2703 ! so the deeper buffer layer must have some mass.
2704 ! (2) The upper buffer layer must have some mass so the extrapolation of
2705 ! density is meaningful (i.e. there is not detrainment from the buffer
2706 ! layers when there is strong mixed layer entrainment).
2707 ! (3) The lower buffer layer density extrapolated to its base with a
2708 ! linear fit between the two layers must exceed the density of the
2709 ! next denser interior layer.
2710 ! (4) The average extrapolated coordinate density that is moved into the
2711 ! isopycnal interior matches the target value for that layer.
2712 ! (5) The potential energy change is calculated and might be used later
2713 ! to allow the upper buffer layer to mix more into the lower buffer
2714 ! layer.
2715
2716 ! Determine whether more must be detrained from the mixed layer to keep a
2717 ! minimal amount of mass in the buffer layers. In this case the 5% of the
2718 ! mixed layer thickness is hard-coded, but probably shouldn't be!
2719 h_min_bl = min(cs%Hbuffer_min, cs%Hbuffer_rel_min*h(i,0))
2720
2721 stable_rcv = .true.
2722 if (cs%nonBous_energetics) then
2723 if (((spv0(i,kb1)-spv0(i,kb2)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) stable_rcv = .false.
2724 else
2725 if (((r0(i,kb2)-r0(i,kb1)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) stable_rcv = .false.
2726 endif
2727
2728 h1 = h(i,kb1) ; h2 = h(i,kb2)
2729
2730 h2_to_k1_rem = (h1 + h2) + h_to_bl
2731 if ((max_bl_det(i) >= 0.0) .and. (h2_to_k1_rem > max_bl_det(i))) &
2732 h2_to_k1_rem = max_bl_det(i)
2733
2734
2735 if ((h2 == 0.0) .and. (h1 > 0.0)) then
2736 ! The lower buffer layer has been eliminated either by convective
2737 ! adjustment or entrainment from the interior, and its current properties
2738 ! are not meaningful, but may later be used to determine the properties of
2739 ! waters moving into the lower buffer layer. So the properties of the
2740 ! lower buffer layer are set to be between those of the upper buffer layer
2741 ! and the next denser interior layer, measured by R0 or SpV0. This probably does
2742 ! not happen very often, so I am not too worried about the inefficiency of
2743 ! the following loop.
2744 do k1=kb2+1,nz ; if (h(i,k1) > 2.0*angstrom) exit ; enddo
2745
2746 rcv(i,kb2) = rcv(i,kb1) ; t(i,kb2) = t(i,kb1) ; s(i,kb2) = s(i,kb1)
2747
2748 if (cs%nonBous_energetics) then
2749 spv0(i,kb2) = spv0(i,kb1)
2750 if (k1 <= nz) then ; if (spv0(i,k1) <= spv0(i,kb1)) then
2751 spv0(i,kb2) = 0.5*(spv0(i,kb1)+spv0(i,k1))
2752
2753 rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
2754 t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
2755 s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
2756 endif ; endif
2757 else
2758 r0(i,kb2) = r0(i,kb1)
2759
2760 if (k1 <= nz) then ; if (r0(i,k1) >= r0(i,kb1)) then
2761 r0(i,kb2) = 0.5*(r0(i,kb1)+r0(i,k1))
2762
2763 rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
2764 t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
2765 s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
2766 endif ; endif
2767 endif
2768 endif ! (h2 = 0 && h1 > 0)
2769
2770 dpe_extrap_rhog = 0.0 ; dpe_extrapolate = 0.0 ; dpe_merge = 0.0 ; dpe_merge_nb = 0.0
2771 mergeable_bl = .false.
2772 if ((h1 > 0.0) .and. (h2 > 0.0) .and. (h_to_bl > 0.0) .and. &
2773 (stable_rcv)) then
2774 ! Check whether it is permissible for the buffer layers to detrain
2775 ! into the interior isopycnal layers.
2776
2777 ! Determine the layer that has the lightest target density that is
2778 ! denser than the lowermost buffer layer.
2779 do k1=kb2+1,nz ; if (rcvtgt(k1) >= rcv(i,kb2)) exit ; enddo ; k0 = k1-1
2780 dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2781
2782 ! Use an energy-balanced combination of downwind advection into the next
2783 ! denser interior layer and upwind advection from the upper buffer layer
2784 ! into the lower one, each with an energy change that equals that required
2785 ! to mix the detrained water with the upper buffer layer.
2786 h1_avail = h1 - max(0.0,h_min_bl-h_to_bl)
2787 if (cs%nonBous_energetics) then
2788 intermediate = (spv0(i,kb1) > spv0(i,kb2)) .and. (h_to_bl*spv0(i,kb1) < spv0_to_bl)
2789 else
2790 intermediate = (r0(i,kb1) < r0(i,kb2)) .and. (h_to_bl*r0(i,kb1) > r0_to_bl)
2791 endif
2792
2793 if ((k1<=nz) .and. (h2 > h_min_bl) .and. (h1_avail > 0.0) .and. intermediate) then
2794 if (cs%nonBous_energetics) then
2795 dspvk1 = (rcvtgt(k1) - rcv(i,kb2)) * (spv0(i,kb2) - spv0(i,kb1)) / &
2796 (rcv(i,kb2) - rcv(i,kb1))
2797 b1 = (rcvtgt(k1) - rcv(i,kb2)) / (rcv(i,kb2) - rcv(i,kb1))
2798 else
2799 drk1 = (rcvtgt(k1) - rcv(i,kb2)) * (r0(i,kb2) - r0(i,kb1)) / &
2800 (rcv(i,kb2) - rcv(i,kb1))
2801 b1 = drk1 / (r0(i,kb2) - r0(i,kb1))
2802 ! b1 = RcvTgt(k1) - Rcv(i,kb2)) / (Rcv(i,kb2) - Rcv(i,kb1))
2803 endif
2804
2805 ! Apply several limits to the detrainment.
2806 ! Entrain less than the mass in h2, and keep the base of the buffer
2807 ! layers from becoming shallower than any neighbors.
2808 h2_to_k1 = min(h2 - h_min_bl, h2_to_k1_rem)
2809 ! Balance downwind advection of density into the layer below the
2810 ! buffer layers with upwind advection from the layer above.
2811 if (h2_to_k1*(h1_avail + b1*(h1_avail + h2)) > h2*h1_avail) &
2812 h2_to_k1 = (h2*h1_avail) / (h1_avail + b1*(h1_avail + h2))
2813
2814 if (cs%nonBous_energetics) then
2815 if (h2_to_k1*(dspvk1 * h2) < (h_to_bl*spv0(i,kb1) - spv0_to_bl) * h1) &
2816 h2_to_k1 = (h_to_bl*spv0(i,kb1) - spv0_to_bl) * h1 / (dspvk1 * h2)
2817 else
2818 if (h2_to_k1*(drk1 * h2) > (h_to_bl*r0(i,kb1) - r0_to_bl) * h1) &
2819 h2_to_k1 = (h_to_bl*r0(i,kb1) - r0_to_bl) * h1 / (drk1 * h2)
2820 endif
2821
2822 if ((k1==kb2+1) .and. (cs%BL_extrap_lim > 0.)) then
2823 ! Simply do not detrain very light water into the lightest isopycnal
2824 ! coordinate layers if the density jump is too large.
2825 drcv_lim = rcv(i,kb2)-rcv(i,0)
2826 do k=1,kb2 ; drcv_lim = max(drcv_lim, rcv(i,kb2)-rcv(i,k)) ; enddo
2827 drcv_lim = cs%BL_extrap_lim*drcv_lim
2828 if ((rcvtgt(k1) - rcv(i,kb2)) >= drcv_lim) then
2829 h2_to_k1 = 0.0
2830 elseif ((rcvtgt(k1) - rcv(i,kb2)) > 0.5*drcv_lim) then
2831 h2_to_k1 = h2_to_k1 * (2.0 - 2.0*((rcvtgt(k1) - rcv(i,kb2)) / drcv_lim))
2832 endif
2833 endif
2834
2835 drcv = (rcvtgt(k1) - rcv(i,kb2))
2836
2837 ! Use 2nd order upwind advection of spiciness, limited by the values
2838 ! in deeper thick layers to determine the detrained temperature and
2839 ! salinity.
2840 dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2841 dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2842 (h2 - h2_to_k1) / (h1 + h2)
2843 dspice_lim = 0.0
2844 if (h(i,k1) > 10.0*angstrom) then
2845 dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2846 dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2847 if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2848 endif
2849 if (k1<nz) then ; if (h(i,k1+1) > 10.0*angstrom) then
2850 dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2851 dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2852 if ((dspice_det*dspice_lim2 > 0.0) .and. &
2853 (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2854 endif ; endif
2855 if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2856
2857 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2858 t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2859 (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2860 s_det = s(i,kb2) + i_denom * &
2861 (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2862
2863 ! The detrained values of R0 or SpV0 are based on changes in T and S.
2864 if (cs%nonBous_energetics) then
2865 spv0_det = spv0(i,kb2) + (t_det-t(i,kb2)) * dspv0_dt(i) + &
2866 (s_det-s(i,kb2)) * dspv0_ds(i)
2867 else
2868 r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2869 (s_det-s(i,kb2)) * dr0_ds(i)
2870 endif
2871
2872 if (cs%BL_extrap_lim >= 0.) then
2873 ! Only do this detrainment if the new layer's temperature and salinity
2874 ! are not too far outside of the range of previous values.
2875 if (h(i,k1) > 10.0*angstrom) then
2876 t_min = min(t(i,kb1), t(i,kb2), t(i,k1)) - cs%Allowed_T_chg
2877 t_max = max(t(i,kb1), t(i,kb2), t(i,k1)) + cs%Allowed_T_chg
2878 s_min = min(s(i,kb1), s(i,kb2), s(i,k1)) - cs%Allowed_S_chg
2879 s_max = max(s(i,kb1), s(i,kb2), s(i,k1)) + cs%Allowed_S_chg
2880 else
2881 t_min = min(t(i,kb1), t(i,kb2)) - cs%Allowed_T_chg
2882 t_max = max(t(i,kb1), t(i,kb2)) + cs%Allowed_T_chg
2883 s_min = min(s(i,kb1), s(i,kb2)) - cs%Allowed_S_chg
2884 s_max = max(s(i,kb1), s(i,kb2)) + cs%Allowed_S_chg
2885 endif
2886 ihk1 = 1.0 / (h(i,k1) + h2_to_k1)
2887 t_new = (h(i,k1)*t(i,k1) + h2_to_k1*t_det) * ihk1
2888 s_new = (h(i,k1)*s(i,k1) + h2_to_k1*s_det) * ihk1
2889 ! A less restrictive limit might be used here.
2890 if ((t_new < t_min) .or. (t_new > t_max) .or. &
2891 (s_new < s_min) .or. (s_new > s_max)) &
2892 h2_to_k1 = 0.0
2893 endif
2894
2895 h1_to_h2 = b1*h2*h2_to_k1 / (h2 - (1.0+b1)*h2_to_k1)
2896
2897 ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2898 ih2f = 1.0 / ((h(i,kb2) - h2_to_k1) + h1_to_h2)
2899
2900 rcv(i,kb2) = ((h(i,kb2)*rcv(i,kb2) - h2_to_k1*rcvtgt(k1)) + &
2901 h1_to_h2*rcv(i,kb1))*ih2f
2902 rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2903
2904 t(i,kb2) = ((h(i,kb2)*t(i,kb2) - h2_to_k1*t_det) + &
2905 h1_to_h2*t(i,kb1)) * ih2f
2906 t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2907
2908 s(i,kb2) = ((h(i,kb2)*s(i,kb2) - h2_to_k1*s_det) + &
2909 h1_to_h2*s(i,kb1)) * ih2f
2910 s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2911
2912 ! Changes in R0 or SpV0 are based on changes in T and S.
2913 if (cs%nonBous_energetics) then
2914 spv0(i,kb2) = ((h(i,kb2)*spv0(i,kb2) - h2_to_k1*spv0_det) + h1_to_h2*spv0(i,kb1)) * ih2f
2915 spv0(i,k1) = ((h(i,k1)+h_neglect)*spv0(i,k1) + h2_to_k1*spv0_det) * ihk1
2916 else
2917 r0(i,kb2) = ((h(i,kb2)*r0(i,kb2) - h2_to_k1*r0_det) + h1_to_h2*r0(i,kb1)) * ih2f
2918 r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2919 endif
2920
2921 h(i,kb1) = h(i,kb1) - h1_to_h2 ; h1 = h(i,kb1)
2922 h(i,kb2) = (h(i,kb2) - h2_to_k1) + h1_to_h2 ; h2 = h(i,kb2)
2923 h(i,k1) = h(i,k1) + h2_to_k1
2924
2925 d_ea(i,kb1) = d_ea(i,kb1) - h1_to_h2
2926 d_ea(i,kb2) = (d_ea(i,kb2) - h2_to_k1) + h1_to_h2
2927 d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2928 h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2929
2930 ! The lower buffer layer has become lighter - it may be necessary to
2931 ! adjust k1 lighter.
2932 if ((k1>kb2+1) .and. (rcvtgt(k1-1) >= rcv(i,kb2))) then
2933 do k1=k1,kb2+1,-1 ; if (rcvtgt(k1-1) < rcv(i,kb2)) exit ; enddo
2934 endif
2935 endif
2936
2937 k0 = k1-1
2938 dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2939
2940 if (cs%nonBous_energetics) then
2941 stable = (spv0(i,kb2) < spv0(i,kb1))
2942 else
2943 stable = (r0(i,kb2) > r0(i,kb1))
2944 endif
2945
2946 if ((k0>kb2) .and. (dr1 > 0.0) .and. (h1 > h_min_bl) .and. (h2*dr2 < h1*dr1) .and. stable) then
2947 ! An interior isopycnal layer (k0) is intermediate in density between
2948 ! the two buffer layers, and there can be detrainment. The entire
2949 ! lower buffer layer is combined with a portion of the upper buffer
2950 ! layer to match the target density of layer k0.
2951 stays_merge = 2.0*(h1+h2)*(h1*dr1 - h2*dr2) / &
2952 ((dr1+dr2)*h1 + dr1*(h1+h2) + &
2953 sqrt((dr2*h1-dr1*h2)**2 + 4*(h1+h2)*h2*(dr1+dr2)*dr2))
2954
2955 if (cs%nonBous_energetics) then
2956 stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
2957 h1 - (h1+h2)*(spv0(i,kb1) - spv0_det) / (spv0(i,kb2) - spv0(i,kb1)))
2958 if ((stays_merge > stays_min_merge) .and. (stays_merge + h2_to_k1_rem >= h1 + h2)) then
2959 mergeable_bl = .true.
2960 dpe_merge_nb = g_2*gv%H_to_RZ**2*(spv0(i,kb1)-spv0(i,kb2)) * ((h1-stays_merge)*(h2-stays_merge))
2961 endif
2962 else
2963 stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
2964 h1 - (h1+h2)*(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1)))
2965 if ((stays_merge > stays_min_merge) .and. (stays_merge + h2_to_k1_rem >= h1 + h2)) then
2966 mergeable_bl = .true.
2967 dpe_merge = g_2*(r0(i,kb2)-r0(i,kb1)) * (h1-stays_merge)*(h2-stays_merge)
2968 endif
2969 endif
2970 endif
2971
2972 if ((k1<=nz).and.(.not.mergeable_bl)) then
2973 ! Check whether linear extrapolation of density (i.e. 2nd order upwind
2974 ! advection) will allow some of the lower buffer layer to detrain into
2975 ! the next denser interior layer (k1).
2976 dr2b = rcvtgt(k1)-rcv(i,kb2) ; dr21 = rcv(i,kb2) - rcv(i,kb1)
2977 if (dr2b*(h1+h2) < h2*dr21) then
2978 ! Some of layer kb2 is denser than k1.
2979 h2_to_k1 = min(h2 - (h1+h2) * dr2b / dr21, h2_to_k1_rem)
2980
2981 if (h2 > h2_to_k1) then
2982 drcv = (rcvtgt(k1) - rcv(i,kb2))
2983
2984 ! Use 2nd order upwind advection of spiciness, limited by the values
2985 ! in deeper thick layers to determine the detrained temperature and
2986 ! salinity.
2987 dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2988 dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2989 (h2 - h2_to_k1) / (h1 + h2)
2990 dspice_lim = 0.0
2991 if (h(i,k1) > 10.0*angstrom) then
2992 dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2993 dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2994 if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2995 endif
2996 if (k1<nz) then ; if (h(i,k1+1) > 10.0*angstrom) then
2997 dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2998 dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2999 if ((dspice_det*dspice_lim2 > 0.0) .and. &
3000 (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
3001 endif ; endif
3002 if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
3003
3004 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
3005 t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
3006 (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
3007 s_det = s(i,kb2) + i_denom * &
3008 (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
3009 ! The detrained values of R0 or SpV0 are based on changes in T and S.
3010 if (cs%nonBous_energetics) then
3011 spv0_det = spv0(i,kb2) + (t_det-t(i,kb2)) * dspv0_dt(i) + &
3012 (s_det-s(i,kb2)) * dspv0_ds(i)
3013 else
3014 r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
3015 (s_det-s(i,kb2)) * dr0_ds(i)
3016 endif
3017
3018 ! Now that the properties of the detrained water are known,
3019 ! potentially limit the amount of water that is detrained to
3020 ! avoid creating unphysical properties in the remaining water.
3021 ih2f = 1.0 / (h2 - h2_to_k1)
3022
3023 t_min = min(t(i,kb2), t(i,kb1)) - cs%Allowed_T_chg
3024 t_max = max(t(i,kb2), t(i,kb1)) + cs%Allowed_T_chg
3025 t_new = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
3026 if (t_new < t_min) then
3027 h2_to_k1_lim = h2 * (t(i,kb2) - t_min) / (t_det - t_min)
3028! write(mesg,'("Low temperature limits det to ", &
3029! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
3030! & 5(1pe12.5))') &
3031! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
3032! T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_min
3033! call MOM_error(WARNING, mesg)
3034 h2_to_k1 = h2_to_k1_lim
3035 ih2f = 1.0 / (h2 - h2_to_k1)
3036 elseif (t_new > t_max) then
3037 h2_to_k1_lim = h2 * (t(i,kb2) - t_max) / (t_det - t_max)
3038! write(mesg,'("High temperature limits det to ", &
3039! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
3040! & 5(1pe12.5))') &
3041! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
3042! T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_max
3043! call MOM_error(WARNING, mesg)
3044 h2_to_k1 = h2_to_k1_lim
3045 ih2f = 1.0 / (h2 - h2_to_k1)
3046 endif
3047 s_min = max(min(s(i,kb2), s(i,kb1)) - cs%Allowed_S_chg, 0.0)
3048 s_max = max(s(i,kb2), s(i,kb1)) + cs%Allowed_S_chg
3049 s_new = (h2*s(i,kb2) - h2_to_k1*s_det)*ih2f
3050 if (s_new < s_min) then
3051 h2_to_k1_lim = h2 * (s(i,kb2) - s_min) / (s_det - s_min)
3052! write(mesg,'("Low salinity limits det to ", &
3053! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
3054! & 5(1pe12.5))') &
3055! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
3056! S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_min
3057! call MOM_error(WARNING, mesg)
3058 h2_to_k1 = h2_to_k1_lim
3059 ih2f = 1.0 / (h2 - h2_to_k1)
3060 elseif (s_new > s_max) then
3061 h2_to_k1_lim = h2 * (s(i,kb2) - s_max) / (s_det - s_max)
3062! write(mesg,'("High salinity limits det to ", &
3063! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
3064! & 5(1pe12.5))') &
3065! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
3066! S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_max
3067! call MOM_error(WARNING, mesg)
3068 h2_to_k1 = h2_to_k1_lim
3069 ih2f = 1.0 / (h2 - h2_to_k1)
3070 endif
3071
3072 ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
3073 rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
3074 rcv(i,kb2) = rcv(i,kb2) - h2_to_k1*drcv*ih2f
3075
3076 t(i,kb2) = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
3077 t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
3078
3079 s(i,kb2) = (h2*s(i,kb2) - h2_to_k1*s_det) * ih2f
3080 s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
3081
3082 ! Changes in R0 or SpV0 are based on changes in T and S.
3083 if (cs%nonBous_energetics) then
3084 spv0(i,kb2) = (h2*spv0(i,kb2) - h2_to_k1*spv0_det) * ih2f
3085 spv0(i,k1) = ((h(i,k1)+h_neglect)*spv0(i,k1) + h2_to_k1*spv0_det) * ihk1
3086 else
3087 r0(i,kb2) = (h2*r0(i,kb2) - h2_to_k1*r0_det) * ih2f
3088 r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
3089 endif
3090 else
3091 ! h2==h2_to_k1 can happen if dR2b = 0 exactly, but this is very
3092 ! unlikely. In this case the entirety of layer kb2 is detrained.
3093 h2_to_k1 = h2 ! These 2 lines are probably unnecessary.
3094 ihk1 = 1.0 / (h(i,k1) + h2)
3095
3096 rcv(i,k1) = (h(i,k1)*rcv(i,k1) + h2*rcv(i,kb2)) * ihk1
3097 t(i,k1) = (h(i,k1)*t(i,k1) + h2*t(i,kb2)) * ihk1
3098 s(i,k1) = (h(i,k1)*s(i,k1) + h2*s(i,kb2)) * ihk1
3099 if (cs%nonBous_energetics) then
3100 spv0(i,k1) = (h(i,k1)*spv0(i,k1) + h2*spv0(i,kb2)) * ihk1
3101 else
3102 r0(i,k1) = (h(i,k1)*r0(i,k1) + h2*r0(i,kb2)) * ihk1
3103 endif
3104 endif
3105
3106 h(i,k1) = h(i,k1) + h2_to_k1
3107 h(i,kb2) = h(i,kb2) - h2_to_k1 ; h2 = h(i,kb2)
3108 ! dPE_extrap_rhoG should be positive here.
3109 if (cs%nonBous_energetics) then
3110 dpe_extrap_rhog = 0.5*(spv0(i,kb2)-spv0_det) * (h2_to_k1*h2) / spv0(i,k1)
3111 dpe_extrapolate = 0.5*gv%g_Earth_Z_T2*gv%H_to_RZ**2*(spv0(i,kb2)-spv0_det) * (h2_to_k1*h2)
3112 else
3113 dpe_extrap_rhog = i2rho0*(r0_det-r0(i,kb2))*h2_to_k1*h2
3114 endif
3115
3116 d_ea(i,kb2) = d_ea(i,kb2) - h2_to_k1
3117 d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
3118 h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
3119 endif
3120 endif ! Detrainment by extrapolation.
3121
3122 endif ! Detrainment to the interior at all.
3123
3124 ! Does some of the detrained water go into the lower buffer layer?
3125 h_det_h2 = max(h_min_bl-(h1+h2), 0.0)
3126 if (h_det_h2 > 0.0) then
3127 ! Detrained water will go into both upper and lower buffer layers.
3128 ! h(kb2) will be h_min_bl, but h(kb1) may be larger if there was already
3129 ! ample detrainment; all water in layer kb1 moves into layer kb2.
3130
3131 ! Determine the fluxes between the various layers.
3132 h_det_to_h2 = min(h_to_bl, h_det_h2)
3133 h_ml_to_h2 = h_det_h2 - h_det_to_h2
3134 h_det_to_h1 = h_to_bl - h_det_to_h2
3135 h_ml_to_h1 = max(h_min_bl-h_det_to_h1,0.0)
3136
3137 ih = 1.0/h_min_bl
3138 ihdet = 0.0 ; if (h_to_bl > 0.0) ihdet = 1.0 / h_to_bl
3139 ih1f = 1.0 / (h_det_to_h1 + h_ml_to_h1)
3140
3141 if (cs%nonBous_energetics) then
3142 spv0(i,kb2) = ((h2*spv0(i,kb2) + h1*spv0(i,kb1)) + &
3143 (h_det_to_h2*spv0_to_bl*ihdet + h_ml_to_h2*spv0(i,0))) * ih
3144 spv0(i,kb1) = (h_det_to_h1*spv0_to_bl*ihdet + h_ml_to_h1*spv0(i,0)) * ih1f
3145 else
3146 r0(i,kb2) = ((h2*r0(i,kb2) + h1*r0(i,kb1)) + &
3147 (h_det_to_h2*r0_to_bl*ihdet + h_ml_to_h2*r0(i,0))) * ih
3148 r0(i,kb1) = (h_det_to_h1*r0_to_bl*ihdet + h_ml_to_h1*r0(i,0)) * ih1f
3149 endif
3150
3151 rcv(i,kb2) = ((h2*rcv(i,kb2) + h1*rcv(i,kb1)) + &
3152 (h_det_to_h2*rcv_to_bl*ihdet + h_ml_to_h2*rcv(i,0))) * ih
3153 rcv(i,kb1) = (h_det_to_h1*rcv_to_bl*ihdet + h_ml_to_h1*rcv(i,0)) * ih1f
3154
3155 t(i,kb2) = ((h2*t(i,kb2) + h1*t(i,kb1)) + &
3156 (h_det_to_h2*t_to_bl*ihdet + h_ml_to_h2*t(i,0))) * ih
3157 t(i,kb1) = (h_det_to_h1*t_to_bl*ihdet + h_ml_to_h1*t(i,0)) * ih1f
3158
3159 s(i,kb2) = ((h2*s(i,kb2) + h1*s(i,kb1)) + &
3160 (h_det_to_h2*s_to_bl*ihdet + h_ml_to_h2*s(i,0))) * ih
3161 s(i,kb1) = (h_det_to_h1*s_to_bl*ihdet + h_ml_to_h1*s(i,0)) * ih1f
3162
3163 ! Recall that h1 = h(i,kb1) & h2 = h(i,kb2).
3164 d_ea(i,1) = d_ea(i,1) - (h_ml_to_h1 + h_ml_to_h2)
3165 d_ea(i,kb1) = d_ea(i,kb1) + ((h_det_to_h1 + h_ml_to_h1) - h1)
3166 d_ea(i,kb2) = d_ea(i,kb2) + (h_min_bl - h2)
3167
3168 h(i,kb1) = h_det_to_h1 + h_ml_to_h1 ; h(i,kb2) = h_min_bl
3169 h(i,0) = h(i,0) - (h_ml_to_h1 + h_ml_to_h2)
3170
3171
3172 if (allocated(cs%diag_PE_detrain) .or. allocated(cs%diag_PE_detrain2)) then
3173 if (cs%nonBous_energetics) then
3174 spv0_det = spv0_to_bl*ihdet
3175 s1en = idt_diag * ( -gv%H_to_RZ**2 * g_2 * ((spv0(i,kb2)-spv0(i,kb1))*h1*h2 + &
3176 h_det_to_h2*( (spv0(i,kb1)-spv0_det)*h1 + (spv0(i,kb2)-spv0_det)*h2 ) + &
3177 h_ml_to_h2*( (spv0(i,kb2)-spv0(i,0))*h2 + (spv0(i,kb1)-spv0(i,0))*h1 + &
3178 (spv0_det-spv0(i,0))*h_det_to_h2 ) + &
3179 h_det_to_h1*h_ml_to_h1*(spv0_det-spv0(i,0))) - dpe_extrapolate )
3180
3181 if (allocated(cs%diag_PE_detrain2)) &
3182 cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + s1en + idt_diag*dpe_extrapolate
3183 else
3184 r0_det = r0_to_bl*ihdet
3185 s1en = g_2 * idt_h2 * ( ((r0(i,kb2)-r0(i,kb1))*h1*h2 + &
3186 h_det_to_h2*( (r0(i,kb1)-r0_det)*h1 + (r0(i,kb2)-r0_det)*h2 ) + &
3187 h_ml_to_h2*( (r0(i,kb2)-r0(i,0))*h2 + (r0(i,kb1)-r0(i,0))*h1 + &
3188 (r0_det-r0(i,0))*h_det_to_h2 ) + &
3189 h_det_to_h1*h_ml_to_h1*(r0_det-r0(i,0))) - 2.0*gv%Rho0*dpe_extrap_rhog )
3190
3191 if (allocated(cs%diag_PE_detrain2)) &
3192 cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + s1en + idt_h2*rho0xg*dpe_extrap_rhog
3193 endif
3194
3195 if (allocated(cs%diag_PE_detrain)) &
3196 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + s1en
3197 endif
3198
3199 elseif ((h_to_bl > 0.0) .or. (h1 < h_min_bl) .or. (h2 < h_min_bl)) then
3200 ! Determine how much of the upper buffer layer will be moved into
3201 ! the lower buffer layer and the properties with which it is moving.
3202 ! This implementation assumes a 2nd-order upwind advection of density
3203 ! from the uppermost buffer layer into the next one down.
3204 h_from_ml = h_min_bl + max(h_min_bl-h2,0.0) - h1 - h_to_bl
3205 if (h_from_ml > 0.0) then
3206 ! Some water needs to be moved from the mixed layer so that the upper
3207 ! (and perhaps lower) buffer layers exceed their minimum thicknesses.
3208 if (cs%nonBous_energetics) then
3209 ! The choice of which specific volume to use in the denominator could be revisited.
3210 ! dPE_extrap_rhoG = dPE_extrap_rhoG + 0.5*h_from_ml*(SpV0_to_bl - SpV0(i,0)*h_to_bl) / SpV0(i,0)
3211 dpe_extrap_rhog = dpe_extrap_rhog + 0.5*h_from_ml*(spv0_to_bl - spv0(i,0)*h_to_bl) * &
3212 ( (h_to_bl + h_from_ml) / (spv0_to_bl + h_from_ml*spv0(i,0)) )
3213 dpe_extrapolate = dpe_extrapolate + 0.5*gv%g_Earth_Z_T2*gv%H_to_RZ**2 * &
3214 h_from_ml*(spv0_to_bl - spv0(i,0)*h_to_bl)
3215 spv0_to_bl = spv0_to_bl + h_from_ml*spv0(i,0)
3216 else
3217 dpe_extrap_rhog = dpe_extrap_rhog - i2rho0*h_from_ml*(r0_to_bl - r0(i,0)*h_to_bl)
3218 r0_to_bl = r0_to_bl + h_from_ml*r0(i,0)
3219 endif
3220 rcv_to_bl = rcv_to_bl + h_from_ml*rcv(i,0)
3221 t_to_bl = t_to_bl + h_from_ml*t(i,0)
3222 s_to_bl = s_to_bl + h_from_ml*s(i,0)
3223
3224 h_to_bl = h_to_bl + h_from_ml
3225 h(i,0) = h(i,0) - h_from_ml
3226 d_ea(i,1) = d_ea(i,1) - h_from_ml
3227 endif
3228
3229 ! The absolute value should be unnecessary and 1e9 is just a large number.
3230 b1 = 1.0e9
3231 if (cs%nonBous_energetics) then
3232 if (spv0(i,kb1) - spv0(i,kb2) > 1.0e-9*abs(spv0_det - spv0(i,kb1))) &
3233 b1 = abs(spv0_det - spv0(i,kb1)) / (spv0(i,kb1) - spv0(i,kb2))
3234 else
3235 if (r0(i,kb2) - r0(i,kb1) > 1.0e-9*abs(r0(i,kb1) - r0_det)) &
3236 b1 = abs(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1))
3237 endif
3238 stays_min = max((1.0-b1)*h1 - b1*h2, 0.0, h_min_bl - h_to_bl)
3239 stays_max = h1 - max(h_min_bl-h2,0.0)
3240
3241 scale_slope = 1.0
3242 if (stays_max <= stays_min) then
3243 stays = stays_max
3244 mergeable_bl = .false.
3245 if (stays_max < h1) scale_slope = (h1 - stays_min) / (h1 - stays_max)
3246 else
3247 ! There are numerous temporary variables used here that should not be
3248 ! used outside of this "else" branch: s1, s2, s3sq, I_ya, bh0
3249 bh0 = b1*h_to_bl
3250 i_ya = (h1 + h2) / ((h1 + h2) + h_to_bl)
3251 ! s1 is the amount staying that minimizes the PE increase.
3252 s1 = 0.5*(h1 + (h2 - bh0) * i_ya) ; s2 = h1 - s1
3253
3254 if (s2 < 0.0) then
3255 ! The energy released by detrainment from the lower buffer layer can be
3256 ! used to mix water from the upper buffer layer into the lower one.
3257 s3sq = i_ya*max(bh0*h1-dpe_extrap_rhog, 0.0)
3258 else
3259 s3sq = i_ya*(bh0*h1-min(dpe_extrap_rhog,0.0))
3260 endif
3261
3262 if (s3sq == 0.0) then
3263 ! There is a simple, exact solution to the quadratic equation, namely:
3264 stays = h1 ! This will revert to stays_max later.
3265 elseif (s2*s2 <= s3sq) then
3266 ! There is no solution with 0 PE change - use the minimum energy input.
3267 stays = s1
3268 else
3269 ! The following choose the solutions that are continuous with all water
3270 ! staying in the upper buffer layer when there is no detrainment,
3271 ! namely the + root when s2>0 and the - root otherwise. They also
3272 ! carefully avoid differencing large numbers, using s2 = (h1-s).
3273 if (bh0 <= 0.0) then ; stays = h1
3274 elseif (s2 > 0.0) then
3275 ! stays = s + sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s
3276 if (s1 >= stays_max) then ; stays = stays_max
3277 elseif (s1 >= 0.0) then ; stays = s1 + sqrt(s2*s2 - s3sq)
3278 else ; stays = (h1*(s2-s1) - s3sq) / (-s1 + sqrt(s2*s2 - s3sq))
3279 endif
3280 else
3281 ! stays = s - sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s & stays_min >= 0
3282 if (s1 <= stays_min) then ; stays = stays_min
3283 else ; stays = (h1*(s1-s2) + s3sq) / (s1 + sqrt(s2*s2 - s3sq))
3284 endif
3285 endif
3286 endif
3287
3288 ! Limit the amount that stays so that the motion of water is from the
3289 ! upper buffer layer into the lower, but no more than is in the upper
3290 ! layer, and the water left in the upper layer is no lighter than the
3291 ! detrained water.
3292 if (stays >= stays_max) then ; stays = stays_max
3293 elseif (stays < stays_min) then ; stays = stays_min
3294 endif
3295 endif
3296
3297 if (cs%nonBous_energetics) then
3298 dpe_det_nb = -g_2*gv%H_to_RZ**2*((spv0(i,kb1)*h_to_bl - spv0_to_bl)*stays + &
3299 (spv0(i,kb2)-spv0(i,kb1)) * (h1-stays) * &
3300 (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
3301 dpe_extrapolate
3302 else
3303 dpe_det = g_2*((r0(i,kb1)*h_to_bl - r0_to_bl)*stays + &
3304 (r0(i,kb2)-r0(i,kb1)) * (h1-stays) * &
3305 (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
3306 rho0xg*dpe_extrap_rhog
3307 endif
3308
3309 if (dpe_time_ratio*h_to_bl > h_to_bl+h(i,0)) then
3310 dpe_ratio = (h_to_bl+h(i,0)) / h_to_bl
3311 else
3312 dpe_ratio = dpe_time_ratio
3313 endif
3314
3315 if (cs%nonBous_energetics) then
3316 better_to_merge = (num_events*dpe_ratio*dpe_det_nb > dpe_merge_nb)
3317 else
3318 better_to_merge = (num_events*dpe_ratio*dpe_det > dpe_merge)
3319 endif
3320
3321 if (mergeable_bl .and. better_to_merge) then
3322 ! It is energetically preferable to merge the two buffer layers, detrain
3323 ! them into interior layer (k0), move the remaining upper buffer layer
3324 ! water into the lower buffer layer, and detrain undiluted into the
3325 ! upper buffer layer.
3326 h1_to_k0 = (h1-stays_merge)
3327 stays = max(h_min_bl-h_to_bl,0.0)
3328 h1_to_h2 = stays_merge - stays
3329
3330 ihk0 = 1.0 / ((h1_to_k0 + h2) + h(i,k0))
3331 ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / h1_to_h2
3332 ih12 = 1.0 / (h1 + h2)
3333
3334 drcv_2dz = (rcv(i,kb1) - rcv(i,kb2)) * ih12
3335 drcv_stays = drcv_2dz*(h1_to_k0 + h1_to_h2)
3336 drcv_det = - drcv_2dz*(stays + h1_to_h2)
3337 rcv(i,k0) = ((h1_to_k0*(rcv(i,kb1) + drcv_det) + &
3338 h2*rcv(i,kb2)) + h(i,k0)*rcv(i,k0)) * ihk0
3339 rcv(i,kb2) = rcv(i,kb1) + drcv_2dz*(h1_to_k0-stays)
3340 rcv(i,kb1) = (rcv_to_bl + stays*(rcv(i,kb1) + drcv_stays)) * ih1f
3341
3342 ! Use 2nd order upwind advection of spiciness, limited by the value in
3343 ! the water from the mixed layer to determine the temperature and
3344 ! salinity of the water that stays in the buffer layers.
3345 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
3346 dspice_2dz = (ds_dt_gauge*drcv_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3347 dt_ds_gauge*drcv_dt(i)*(s(i,kb1)-s(i,kb2))) * ih12
3348 if (cs%nonBous_energetics) then
3349 ! Use the specific volume differences to limit the coordinate density change.
3350 dspice_lim = -rcv(i,kb1) * (ds_dt_gauge*dspv0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3351 dt_ds_gauge*dspv0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / (spv0(i,kb1) * h_to_bl)
3352 else
3353 dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3354 dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
3355 endif
3356 if (dspice_lim * dspice_2dz <= 0.0) dspice_2dz = 0.0
3357
3358 if (stays > 0.0) then
3359 ! Limit the spiciness of the water that stays in the upper buffer layer.
3360 if (abs(dspice_lim) < abs(dspice_2dz*(h1_to_k0 + h1_to_h2))) &
3361 dspice_2dz = dspice_lim/(h1_to_k0 + h1_to_h2)
3362
3363 dspice_stays = dspice_2dz*(h1_to_k0 + h1_to_h2)
3364 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3365 (dt_ds_gauge * drcv_dt(i) * drcv_stays + drcv_ds(i) * dspice_stays)
3366 s_stays = s(i,kb1) + i_denom * &
3367 (drcv_ds(i) * drcv_stays - dt_ds_gauge * drcv_dt(i) * dspice_stays)
3368 ! The values of R0 or SpV0 are based on changes in T and S.
3369 if (cs%nonBous_energetics) then
3370 spv0_stays = spv0(i,kb1) + (t_stays-t(i,kb1)) * dspv0_dt(i) + &
3371 (s_stays-s(i,kb1)) * dspv0_ds(i)
3372 else
3373 r0_stays = r0(i,kb1) + (t_stays-t(i,kb1)) * dr0_dt(i) + &
3374 (s_stays-s(i,kb1)) * dr0_ds(i)
3375 endif
3376 else
3377 ! Limit the spiciness of the water that moves into the lower buffer layer.
3378 if (abs(dspice_lim) < abs(dspice_2dz*h1_to_k0)) &
3379 dspice_2dz = dspice_lim/h1_to_k0
3380 ! These will be multiplied by 0 later.
3381 t_stays = 0.0 ; s_stays = 0.0 ; r0_stays = 0.0 ; spv0_stays = 0.0
3382 endif
3383
3384 dspice_det = - dspice_2dz*(stays + h1_to_h2)
3385 t_det = t(i,kb1) + dt_ds_gauge * i_denom * &
3386 (dt_ds_gauge * drcv_dt(i) * drcv_det + drcv_ds(i) * dspice_det)
3387 s_det = s(i,kb1) + i_denom * &
3388 (drcv_ds(i) * drcv_det - dt_ds_gauge * drcv_dt(i) * dspice_det)
3389 ! The values of R0 or SpV0 are based on changes in T and S.
3390 if (cs%nonBous_energetics) then
3391 spv0_det = spv0(i,kb1) + (t_det-t(i,kb1)) * dspv0_dt(i) + &
3392 (s_det-s(i,kb1)) * dspv0_ds(i)
3393 else
3394 r0_det = r0(i,kb1) + (t_det-t(i,kb1)) * dr0_dt(i) + &
3395 (s_det-s(i,kb1)) * dr0_ds(i)
3396 endif
3397
3398 t(i,k0) = ((h1_to_k0*t_det + h2*t(i,kb2)) + h(i,k0)*t(i,k0)) * ihk0
3399 t(i,kb2) = (h1*t(i,kb1) - stays*t_stays - h1_to_k0*t_det) * ih2f
3400 t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3401
3402 s(i,k0) = ((h1_to_k0*s_det + h2*s(i,kb2)) + h(i,k0)*s(i,k0)) * ihk0
3403 s(i,kb2) = (h1*s(i,kb1) - stays*s_stays - h1_to_k0*s_det) * ih2f
3404 s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3405
3406 if (cs%nonBous_energetics) then
3407 spv0(i,k0) = ((h1_to_k0*spv0_det + h2*spv0(i,kb2)) + h(i,k0)*spv0(i,k0)) * ihk0
3408 spv0(i,kb2) = (h1*spv0(i,kb1) - stays*spv0_stays - h1_to_k0*spv0_det) * ih2f
3409 spv0(i,kb1) = (spv0_to_bl + stays*spv0_stays) * ih1f
3410 else
3411 r0(i,k0) = ((h1_to_k0*r0_det + h2*r0(i,kb2)) + h(i,k0)*r0(i,k0)) * ihk0
3412 r0(i,kb2) = (h1*r0(i,kb1) - stays*r0_stays - h1_to_k0*r0_det) * ih2f
3413 r0(i,kb1) = (r0_to_bl + stays*r0_stays) * ih1f
3414 endif
3415
3416! ! The following is 2nd-order upwind advection without limiters.
3417! dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih12
3418! T(i,k0) = (h1_to_k0*(T(i,kb1) - dT_2dz*(stays+h1_to_h2)) + &
3419! h2*T(i,kb2) + h(i,k0)*T(i,k0)) * Ihk0
3420! T(i,kb2) = T(i,kb1) + dT_2dz*(h1_to_k0-stays)
3421! T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + dT_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3422! dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih12
3423! S(i,k0) = (h1_to_k0*(S(i,kb1) - dS_2dz*(stays+h1_to_h2)) + &
3424! h2*S(i,kb2) + h(i,k0)*S(i,k0)) * Ihk0
3425! S(i,kb2) = S(i,kb1) + dS_2dz*(h1_to_k0-stays)
3426! S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + dS_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3427! if (CS%nonBous_energetics) then
3428! dSpV0_2dz = (SpV0(i,kb1) - SpV0(i,kb2)) * Ih12
3429! SpV0(i,k0) = (h1_to_k0*(SpV0(i,kb1) - dSpV0_2dz*(stays+h1_to_h2)) + &
3430! h2*SpV0(i,kb2) + h(i,k0)*SpV0(i,k0)) * Ihk0
3431! SpV0(i,kb2) = SpV0(i,kb1) + dSpV0_2dz*(h1_to_k0-stays)
3432! SpV0(i,kb1) = (SpV0_to_bl + stays*(SpV0(i,kb1) + dSpV0_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3433! else
3434! dR0_2dz = (R0(i,kb1) - R0(i,kb2)) * Ih12
3435! R0(i,k0) = (h1_to_k0*(R0(i,kb1) - dR0_2dz*(stays+h1_to_h2)) + &
3436! h2*R0(i,kb2) + h(i,k0)*R0(i,k0)) * Ihk0
3437! R0(i,kb2) = R0(i,kb1) + dR0_2dz*(h1_to_k0-stays)
3438! R0(i,kb1) = (R0_to_bl + stays*(R0(i,kb1) + dR0_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3439! endif
3440
3441 d_ea(i,kb1) = (d_ea(i,kb1) + h_to_bl) + (stays - h1)
3442 d_ea(i,kb2) = d_ea(i,kb2) + (h1_to_h2 - h2)
3443 d_ea(i,k0) = d_ea(i,k0) + (h1_to_k0 + h2)
3444
3445 h(i,kb1) = stays + h_to_bl
3446 h(i,kb2) = h1_to_h2
3447 h(i,k0) = h(i,k0) + (h1_to_k0 + h2)
3448 if (cs%nonBous_energetics) then
3449 if (allocated(cs%diag_PE_detrain)) &
3450 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_diag*dpe_merge_nb
3451 if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3452 cs%diag_PE_detrain2(i,j) + idt_diag*(dpe_det_nb + dpe_extrapolate)
3453 else
3454 if (allocated(cs%diag_PE_detrain)) &
3455 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_merge
3456 if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3457 cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap_rhog)
3458 endif
3459 else ! Not mergeable_bl.
3460 ! There is no further detrainment from the buffer layers, and the
3461 ! upper buffer layer water is distributed optimally between the
3462 ! upper and lower buffer layer.
3463 h1_to_h2 = h1 - stays
3464 ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / (h2 + h1_to_h2)
3465 ih = 1.0 / (h1 + h2)
3466 if (cs%nonBous_energetics) then
3467 dspv0_2dz = (spv0(i,kb1) - spv0(i,kb2)) * ih
3468 spv0(i,kb2) = (h2*spv0(i,kb2) + h1_to_h2*(spv0(i,kb1) - scale_slope*dspv0_2dz*stays)) * ih2f
3469 spv0(i,kb1) = (spv0_to_bl + stays*(spv0(i,kb1) + scale_slope*dspv0_2dz*h1_to_h2)) * ih1f
3470 else
3471 dr0_2dz = (r0(i,kb1) - r0(i,kb2)) * ih
3472 r0(i,kb2) = (h2*r0(i,kb2) + h1_to_h2*(r0(i,kb1) - scale_slope*dr0_2dz*stays)) * ih2f
3473 r0(i,kb1) = (r0_to_bl + stays*(r0(i,kb1) + scale_slope*dr0_2dz*h1_to_h2)) * ih1f
3474 endif
3475
3476 ! Use 2nd order upwind advection of spiciness, limited by the value in the
3477 ! detrained water to determine the detrained temperature and salinity.
3478 if (cs%nonBous_energetics) then
3479 dspv0 = scale_slope*dspv0_2dz*h1_to_h2
3480 dspicespv_stays = (ds_dt_gauge*dspv0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3481 dt_ds_gauge*dspv0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
3482 scale_slope*h1_to_h2 * ih
3483 if (h_to_bl > 0.0) then
3484 dspicespv_lim = (ds_dt_gauge*dspv0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3485 dt_ds_gauge*dspv0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
3486 else
3487 dspicespv_lim = ds_dt_gauge*dspv0_ds(i)*(t(i,0)-t(i,kb1)) - &
3488 dt_ds_gauge*dspv0_dt(i)*(s(i,0)-s(i,kb1))
3489 endif
3490 if (dspicespv_stays*dspicespv_lim <= 0.0) then
3491 dspicespv_stays = 0.0
3492 elseif (abs(dspicespv_stays) > abs(dspicespv_lim)) then
3493 dspicespv_stays = dspicespv_lim
3494 endif
3495 i_denom = 1.0 / (dspv0_ds(i)**2 + (dt_ds_gauge*dspv0_dt(i))**2)
3496 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3497 (dt_ds_gauge * dspv0_dt(i) * dspv0 + dspv0_ds(i) * dspicespv_stays)
3498 s_stays = s(i,kb1) + i_denom * &
3499 (dspv0_ds(i) * dspv0 - dt_ds_gauge * dspv0_dt(i) * dspicespv_stays)
3500 else
3501 dr0 = scale_slope*dr0_2dz*h1_to_h2
3502 dspice_stays = (ds_dt_gauge*dr0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3503 dt_ds_gauge*dr0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
3504 scale_slope*h1_to_h2 * ih
3505 if (h_to_bl > 0.0) then
3506 dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3507 dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
3508 else
3509 dspice_lim = ds_dt_gauge*dr0_ds(i)*(t(i,0)-t(i,kb1)) - &
3510 dt_ds_gauge*dr0_dt(i)*(s(i,0)-s(i,kb1))
3511 endif
3512 if (dspice_stays*dspice_lim <= 0.0) then
3513 dspice_stays = 0.0
3514 elseif (abs(dspice_stays) > abs(dspice_lim)) then
3515 dspice_stays = dspice_lim
3516 endif
3517 i_denom = 1.0 / (dr0_ds(i)**2 + (dt_ds_gauge*dr0_dt(i))**2)
3518 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3519 (dt_ds_gauge * dr0_dt(i) * dr0 + dr0_ds(i) * dspice_stays)
3520 s_stays = s(i,kb1) + i_denom * &
3521 (dr0_ds(i) * dr0 - dt_ds_gauge * dr0_dt(i) * dspice_stays)
3522 endif
3523
3524 ! The detrained values of Rcv are based on changes in T and S.
3525 rcv_stays = rcv(i,kb1) + (t_stays-t(i,kb1)) * drcv_dt(i) + &
3526 (s_stays-s(i,kb1)) * drcv_ds(i)
3527
3528 t(i,kb2) = (h2*t(i,kb2) + h1*t(i,kb1) - t_stays*stays) * ih2f
3529 t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3530 s(i,kb2) = (h2*s(i,kb2) + h1*s(i,kb1) - s_stays*stays) * ih2f
3531 s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3532 rcv(i,kb2) = (h2*rcv(i,kb2) + h1*rcv(i,kb1) - rcv_stays*stays) * ih2f
3533 rcv(i,kb1) = (rcv_to_bl + stays*rcv_stays) * ih1f
3534
3535! ! The following is 2nd-order upwind advection without limiters.
3536! dRcv_2dz = (Rcv(i,kb1) - Rcv(i,kb2)) * Ih
3537! dRcv = scale_slope*dRcv_2dz*h1_to_h2
3538! Rcv(i,kb2) = (h2*Rcv(i,kb2) + h1_to_h2*(Rcv(i,kb1) - &
3539! scale_slope*dRcv_2dz*stays)) * Ih2f
3540! Rcv(i,kb1) = (Rcv_to_bl + stays*(Rcv(i,kb1) + dRcv)) * Ih1f
3541! dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih
3542! T(i,kb2) = (h2*T(i,kb2) + h1_to_h2*(T(i,kb1) - &
3543! scale_slope*dT_2dz*stays)) * Ih2f
3544! T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
3545! scale_slope*dT_2dz*h1_to_h2)) * Ih1f
3546! dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih
3547! S(i,kb2) = (h2*S(i,kb2) + h1_to_h2*(S(i,kb1) - &
3548! scale_slope*dS_2dz*stays)) * Ih2f
3549! S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
3550! scale_slope*dS_2dz*h1_to_h2)) * Ih1f
3551
3552 d_ea(i,kb1) = d_ea(i,kb1) + ((stays - h1) + h_to_bl)
3553 d_ea(i,kb2) = d_ea(i,kb2) + h1_to_h2
3554
3555 h(i,kb1) = stays + h_to_bl
3556 h(i,kb2) = h(i,kb2) + h1_to_h2
3557
3558 if (cs%nonBous_energetics) then
3559 if (allocated(cs%diag_PE_detrain)) &
3560 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_diag*dpe_det_nb
3561 if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3562 cs%diag_PE_detrain2(i,j) + idt_diag*(dpe_det_nb + dpe_extrapolate)
3563 else
3564 ! Recasting dPE_det into the same units as dPE_det_nB changes these diagnostics slightly
3565 ! in some cases for reasons that are not understood.
3566 if (allocated(cs%diag_PE_detrain)) &
3567 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_det
3568 if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3569 cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap_rhog)
3570 endif
3571 endif
3572 endif ! End of detrainment...
3573
3574 enddo ! i loop
3575
3576end subroutine mixedlayer_detrain_2
3577
3578!> This subroutine moves any water left in the former mixed layers into the
3579!! single buffer layers and may also move buffer layer water into the interior
3580!! isopycnal layers.
3581subroutine mixedlayer_detrain_1(h, T, S, R0, SpV0, Rcv, RcvTgt, dt, dt_diag, d_ea, d_eb, &
3582 j, G, GV, US, CS, dRcv_dT, dRcv_dS, max_BL_det)
3583 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3584 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3585 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
3586 !! Layer 0 is the new mixed layer.
3587 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Potential temperature [C ~> degC].
3588 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Salinity [S ~> ppt].
3589 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
3590 !! surface pressure [R ~> kg m-3].
3591 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: SpV0 !< Specific volume referenced to
3592 !! surface pressure [R-1 ~> m3 kg-1]
3593 real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential
3594 !! density [R ~> kg m-3].
3595 real, dimension(SZK_(GV)), intent(in) :: RcvTgt !< The target value of Rcv for each
3596 !! layer [R ~> kg m-3].
3597 real, intent(in) :: dt !< Time increment [T ~> s].
3598 real, intent(in) :: dt_diag !< The accumulated time interval for
3599 !! diagnostics [T ~> s].
3600 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a layer in
3601 !! the entrainment from above
3602 !! [H ~> m or kg m-2]. Positive d_ea
3603 !! goes with layer thickness increases.
3604 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
3605 !! in the entrainment from below [H ~> m or kg m-2].
3606 !! Positive values go with mass gain by
3607 !! a layer.
3608 integer, intent(in) :: j !< The meridional row to work on.
3609 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3610 type(bulkmixedlayer_cs), intent(inout) :: CS !< Bulk mixed layer control structure
3611 real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of
3612 !! coordinate defining potential density
3613 !! with potential temperature
3614 !! [R C-1 ~> kg m-3 degC-1].
3615 real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of
3616 !! coordinate defining potential density
3617 !! with salinity [R S-1 ~> kg m-3 ppt-1].
3618 real, dimension(SZI_(G)), intent(in) :: max_BL_det !< If non-negative, the maximum
3619 !! detrainment permitted from the buffer
3620 !! layers [H ~> m or kg m-2].
3621
3622 ! Local variables
3623 real :: Ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
3624 real :: h_ent ! The thickness from a layer that is
3625 ! entrained [H ~> m or kg m-2].
3626 real :: max_det_rem(SZI_(G)) ! Remaining permitted detrainment [H ~> m or kg m-2].
3627 real :: detrain(SZI_(G)) ! The thickness of fluid to detrain
3628 ! from the mixed layer [H ~> m or kg m-2].
3629 real :: dT_dS_wt2 ! The square of the relative weighting of temperature and salinity changes
3630 ! when extraploating to match a target density [C2 S-2 ~> degC2 ppt-2]
3631 real :: dT_dR ! The ratio of temperature changes to density changes when
3632 ! extrapolating [C R-1 ~> degC m3 kg-1]
3633 real :: dS_dR ! The ratio of salinity changes to density changes when
3634 ! extrapolating [S R-1 ~> ppt m3 kg-1]
3635 real :: dRml ! The density range within the extent of the mixed layers [R ~> kg m-3]
3636 real :: dR0_dRcv ! The relative changes in the potential density and the coordinate density [nondim]
3637 real :: dSpV0_dRcv ! The relative changes in the specific volume and the coordinate density [R-2 ~> m6 kg-2]
3638 real :: I_denom ! A work variable [S2 R-2 ~> ppt2 m6 kg-2].
3639 real :: Sdown ! The salinity of the detrained water [S ~> ppt]
3640 real :: Tdown ! The temperature of the detrained water [C ~> degC]
3641 real :: dt_Time ! The timestep divided by the detrainment timescale [nondim].
3642 real :: g_H_2Rho0dt ! Half the gravitational acceleration times the
3643 ! conversion from H to m divided by the mean density times the time
3644 ! step [Z2 T-3 H-1 R-1 ~> m4 s-3 kg-1 or m7 s-3 kg-2].
3645 real :: g_H2_2dt ! Half the gravitational acceleration times the square of the
3646 ! conversion from H to Z divided by the diagnostic time step
3647 ! [Z3 H-2 T-3 ~> m s-3 or m7 kg-2 s-3].
3648 real :: nB_g_H_2dt ! Half the gravitational acceleration times the conversion from
3649 ! H to RZ divided by the diagnostic time step
3650 ! [R Z2 H-1 T-3 ~> kg m-2 s-3 or m s-3].
3651 real :: nB_gRZ_H2_2dt ! Half the gravitational acceleration times the conversion from
3652 ! H to RZ squared divided by the diagnostic time step
3653 ! [R2 Z3 H-2 T-3 ~> kg2 m-5 s-3 or m s-3]
3654 real :: x1 ! A temporary work variable [various]
3655 logical :: splittable_BL(SZI_(G)), orthogonal_extrap
3656 logical :: must_unmix
3657 integer :: i, is, ie, k, k1, nkmb, nz
3658
3659 is = g%isc ; ie = g%iec ; nz = gv%ke
3660 nkmb = cs%nkml+cs%nkbl
3661 if (cs%nkbl /= 1) call mom_error(fatal,"MOM_mixed_layer: "// &
3662 "CS%nkbl must be 1 in mixedlayer_detrain_1.")
3663
3664 dt_time = dt / cs%BL_detrain_time
3665
3666 if (cs%nonBous_energetics) then
3667 nb_g_h_2dt = (gv%g_Earth_Z_T2 * gv%H_to_RZ) / (2.0 * dt_diag)
3668 nb_grz_h2_2dt = gv%H_to_RZ * nb_g_h_2dt
3669 else
3670 g_h2_2dt = (gv%g_Earth_Z_T2 * gv%H_to_Z**2) / (2.0 * dt_diag)
3671 g_h_2rho0dt = g_h2_2dt * gv%RZ_to_H
3672 endif
3673
3674 ! Move detrained water into the buffer layer.
3675 do k=1,cs%nkml
3676 do i=is,ie ; if (h(i,k) > 0.0) then
3677 ih = 1.0 / (h(i,nkmb) + h(i,k))
3678
3679 if (cs%nonBous_energetics) then
3680 if (cs%TKE_diagnostics) &
3681 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) - &
3682 nb_g_h_2dt * (h(i,k) * h(i,nkmb)) * (spv0(i,nkmb) - spv0(i,k))
3683 if (allocated(cs%diag_PE_detrain)) &
3684 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) - &
3685 nb_grz_h2_2dt * (h(i,k) * h(i,nkmb)) * (spv0(i,nkmb) - spv0(i,k))
3686 if (allocated(cs%diag_PE_detrain2)) &
3687 cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) - &
3688 nb_grz_h2_2dt * (h(i,k) * h(i,nkmb)) * (spv0(i,nkmb) - spv0(i,k))
3689
3690 spv0(i,nkmb) = (spv0(i,nkmb)*h(i,nkmb) + spv0(i,k)*h(i,k)) * ih
3691 else
3692 if (cs%TKE_diagnostics) &
3693 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
3694 g_h_2rho0dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3695 if (allocated(cs%diag_PE_detrain)) &
3696 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + &
3697 g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3698 if (allocated(cs%diag_PE_detrain2)) &
3699 cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + &
3700 g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3701
3702 r0(i,nkmb) = (r0(i,nkmb)*h(i,nkmb) + r0(i,k)*h(i,k)) * ih
3703 endif
3704 rcv(i,nkmb) = (rcv(i,nkmb)*h(i,nkmb) + rcv(i,k)*h(i,k)) * ih
3705 t(i,nkmb) = (t(i,nkmb)*h(i,nkmb) + t(i,k)*h(i,k)) * ih
3706 s(i,nkmb) = (s(i,nkmb)*h(i,nkmb) + s(i,k)*h(i,k)) * ih
3707
3708 d_ea(i,k) = d_ea(i,k) - h(i,k)
3709 d_ea(i,nkmb) = d_ea(i,nkmb) + h(i,k)
3710 h(i,nkmb) = h(i,nkmb) + h(i,k)
3711 h(i,k) = 0.0
3712 endif ; enddo
3713 enddo
3714
3715 do i=is,ie
3716 max_det_rem(i) = 10.0 * h(i,nkmb)
3717 if (max_bl_det(i) >= 0.0) max_det_rem(i) = max_bl_det(i)
3718 enddo
3719
3720! If the mixed layer was denser than the densest interior layer,
3721! but is now lighter than this layer, leaving a buffer layer that
3722! is denser than this layer, there are problems. This should prob-
3723! ably be considered a case of an inadequate choice of resolution in
3724! density space and should be avoided. To make the model run sens-
3725! ibly in this case, it will make the mixed layer denser while making
3726! the buffer layer the density of the densest interior layer (pro-
3727! vided that the this will not make the mixed layer denser than the
3728! interior layer). Otherwise, make the mixed layer the same density
3729! as the densest interior layer and lighten the buffer layer with
3730! the released buoyancy. With multiple buffer layers, much more
3731! graceful options are available.
3732 do i=is,ie ; if (h(i,nkmb) > 0.0) then
3733 if (cs%nonBous_energetics) then
3734 must_unmix = (spv0(i,0) > spv0(i,nz)) .and. (spv0(i,nz) > spv0(i,nkmb))
3735 else
3736 must_unmix = (r0(i,0) < r0(i,nz)) .and. (r0(i,nz) < r0(i,nkmb))
3737 endif
3738 if (must_unmix) then
3739 if (cs%nonBous_energetics) then
3740 if ((spv0(i,0)-spv0(i,nz))*h(i,0) > (spv0(i,nz)-spv0(i,nkmb))*h(i,nkmb)) then
3741 detrain(i) = (spv0(i,nz)-spv0(i,nkmb))*h(i,nkmb) / (spv0(i,0)-spv0(i,nkmb))
3742 else
3743 detrain(i) = (spv0(i,0)-spv0(i,nz))*h(i,0) / (spv0(i,0)-spv0(i,nkmb))
3744 endif
3745 else
3746 if ((r0(i,nz)-r0(i,0))*h(i,0) > (r0(i,nkmb)-r0(i,nz))*h(i,nkmb)) then
3747 detrain(i) = (r0(i,nkmb)-r0(i,nz))*h(i,nkmb) / (r0(i,nkmb)-r0(i,0))
3748 else
3749 detrain(i) = (r0(i,nz)-r0(i,0))*h(i,0) / (r0(i,nkmb)-r0(i,0))
3750 endif
3751 endif
3752
3753 d_eb(i,cs%nkml) = d_eb(i,cs%nkml) + detrain(i)
3754 d_ea(i,cs%nkml) = d_ea(i,cs%nkml) - detrain(i)
3755 d_eb(i,nkmb) = d_eb(i,nkmb) - detrain(i)
3756 d_ea(i,nkmb) = d_ea(i,nkmb) + detrain(i)
3757
3758 if (cs%nonBous_energetics) then
3759 if (allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3760 cs%diag_PE_detrain(i,j) - nb_grz_h2_2dt * detrain(i)* &
3761 (h(i,0) + h(i,nkmb)) * (spv0(i,nkmb) - spv0(i,0))
3762 x1 = spv0(i,0)
3763 spv0(i,0) = spv0(i,0) - detrain(i)*(spv0(i,0)-spv0(i,nkmb)) / h(i,0)
3764 spv0(i,nkmb) = spv0(i,nkmb) - detrain(i)*(spv0(i,nkmb)-x1) / h(i,nkmb)
3765 else
3766 if (allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3767 cs%diag_PE_detrain(i,j) + g_h2_2dt * detrain(i)* &
3768 (h(i,0) + h(i,nkmb)) * (r0(i,nkmb) - r0(i,0))
3769 x1 = r0(i,0)
3770 r0(i,0) = r0(i,0) - detrain(i)*(r0(i,0)-r0(i,nkmb)) / h(i,0)
3771 r0(i,nkmb) = r0(i,nkmb) - detrain(i)*(r0(i,nkmb)-x1) / h(i,nkmb)
3772 endif
3773
3774 x1 = rcv(i,0)
3775 rcv(i,0) = rcv(i,0) - detrain(i)*(rcv(i,0)-rcv(i,nkmb)) / h(i,0)
3776 rcv(i,nkmb) = rcv(i,nkmb) - detrain(i)*(rcv(i,nkmb)-x1) / h(i,nkmb)
3777 x1 = t(i,0)
3778 t(i,0) = t(i,0) - detrain(i)*(t(i,0)-t(i,nkmb)) / h(i,0)
3779 t(i,nkmb) = t(i,nkmb) - detrain(i)*(t(i,nkmb)-x1) / h(i,nkmb)
3780 x1 = s(i,0)
3781 s(i,0) = s(i,0) - detrain(i)*(s(i,0)-s(i,nkmb)) / h(i,0)
3782 s(i,nkmb) = s(i,nkmb) - detrain(i)*(s(i,nkmb)-x1) / h(i,nkmb)
3783
3784 endif
3785 endif ; enddo
3786
3787 ! Move water out of the buffer layer, if convenient.
3788! Split the buffer layer if possible, and replace the buffer layer
3789! with a small amount of fluid from the mixed layer.
3790! This is the exponential-in-time splitting, circa 2005.
3791 do i=is,ie
3792 if (h(i,nkmb) > 0.0) then ; splittable_bl(i) = .true.
3793 else ; splittable_bl(i) = .false. ; endif
3794 enddo
3795
3796 dt_ds_wt2 = cs%dT_dS_wt**2
3797
3798 do k=nz-1,nkmb+1,-1 ; do i=is,ie
3799 if (splittable_bl(i)) then
3800 if (rcvtgt(k) <= rcv(i,nkmb)) then
3801! Estimate dR/drho, dTheta/dR, and dS/dR, where R is the coordinate variable
3802! and rho is in-situ (or surface) potential density.
3803! There is no "right" way to do this, so this keeps things reasonable, if
3804! slightly arbitrary.
3805 splittable_bl(i) = .false.
3806
3807 k1 = k+1 ; orthogonal_extrap = .false.
3808 ! Here we try to find a massive layer to use for interpolating the
3809 ! temperature and salinity. If none is available a pseudo-orthogonal
3810 ! extrapolation is used. The 10.0 and 0.9 in the following are
3811 ! arbitrary but probably about right.
3812 if ((h(i,k+1) < 10.0*gv%Angstrom_H) .or. &
3813 ((rcvtgt(k+1)-rcv(i,nkmb)) >= 0.9*(rcv(i,k1) - rcv(i,0)))) then
3814 if (k>=nz-1) then ; orthogonal_extrap = .true.
3815 elseif ((h(i,k+2) <= 10.0*gv%Angstrom_H) .and. &
3816 ((rcvtgt(k+1)-rcv(i,nkmb)) < 0.9*(rcv(i,k+2)-rcv(i,0)))) then
3817 k1 = k+2
3818 else ; orthogonal_extrap = .true. ; endif
3819 endif
3820
3821 ! Check for the case when there is an inversion of in-situ density relative to
3822 ! the coordinate variable. Do not detrain from the buffer layer in this case.
3823 if (cs%nonBous_energetics) then
3824 if ((spv0(i,0) <= spv0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
3825 else
3826 if ((r0(i,0) >= r0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
3827 endif
3828
3829 if (orthogonal_extrap) then
3830 ! 36 here is a typical oceanic value of (dR/dS) / (dR/dT) - it says
3831 ! that the relative weights of T & S changes is a plausible 6:1.
3832 ! Also, this was coded on Athena's 6th birthday!
3833 i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
3834 dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
3835 ds_dr = drcv_ds(i) * i_denom
3836 else
3837 dt_dr = (t(i,0) - t(i,k1)) / (rcv(i,0) - rcv(i,k1))
3838 ds_dr = (s(i,0) - s(i,k1)) / (rcv(i,0) - rcv(i,k1))
3839 endif
3840
3841 if (cs%nonBous_energetics) then
3842 drml = dt_time * (spv0(i,0) - spv0(i,nkmb)) * &
3843 (rcv(i,0) - rcv(i,k1)) / (spv0(i,k1) - spv0(i,0))
3844 if (drml < 0.0) cycle ! Once again, there is an apparent density inversion in Rcv.
3845 dspv0_drcv = (spv0(i,0) - spv0(i,k1)) / (rcv(i,0) - rcv(i,k1))
3846 else
3847 drml = dt_time * (r0(i,nkmb) - r0(i,0)) * &
3848 (rcv(i,0) - rcv(i,k1)) / (r0(i,0) - r0(i,k1))
3849 if (drml < 0.0) cycle ! Once again, there is an apparent density inversion in Rcv.
3850 dr0_drcv = (r0(i,0) - r0(i,k1)) / (rcv(i,0) - rcv(i,k1))
3851 endif
3852
3853 if ((rcv(i,nkmb) - drml < rcvtgt(k)) .and. (max_det_rem(i) > h(i,nkmb))) then
3854 ! In this case, the buffer layer is split into two isopycnal layers.
3855 detrain(i) = h(i,nkmb) * (rcv(i,nkmb) - rcvtgt(k)) / &
3856 (rcvtgt(k+1) - rcvtgt(k))
3857
3858 if (allocated(cs%diag_PE_detrain)) then
3859 if (cs%nonBous_energetics) then
3860 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + nb_grz_h2_2dt * detrain(i) * &
3861 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dspv0_drcv
3862 else
3863 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * &
3864 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dr0_drcv
3865 endif
3866 endif
3867
3868 tdown = detrain(i) * (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3869 t(i,k) = (h(i,k) * t(i,k) + &
3870 (h(i,nkmb) * t(i,nkmb) - tdown)) / &
3871 (h(i,k) + (h(i,nkmb) - detrain(i)))
3872 t(i,k+1) = (h(i,k+1) * t(i,k+1) + tdown)/ &
3873 (h(i,k+1) + detrain(i))
3874 t(i,nkmb) = t(i,0)
3875 sdown = detrain(i) * (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3876 s(i,k) = (h(i,k) * s(i,k) + &
3877 (h(i,nkmb) * s(i,nkmb) - sdown)) / &
3878 (h(i,k) + (h(i,nkmb) - detrain(i)))
3879 s(i,k+1) = (h(i,k+1) * s(i,k+1) + sdown)/ &
3880 (h(i,k+1) + detrain(i))
3881 s(i,nkmb) = s(i,0)
3882 rcv(i,nkmb) = rcv(i,0)
3883
3884 d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3885 d_ea(i,k) = d_ea(i,k) + (h(i,nkmb) - detrain(i))
3886 d_ea(i,nkmb) = d_ea(i,nkmb) - h(i,nkmb)
3887
3888 h(i,k+1) = h(i,k+1) + detrain(i)
3889 h(i,k) = h(i,k) + h(i,nkmb) - detrain(i)
3890 h(i,nkmb) = 0.0
3891 else
3892 ! Here only part of the buffer layer is moved into the interior.
3893 detrain(i) = h(i,nkmb) * drml / (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3894 if (detrain(i) > max_det_rem(i)) detrain(i) = max_det_rem(i)
3895 ih = 1.0 / (h(i,k+1) + detrain(i))
3896
3897 tdown = (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3898 t(i,nkmb) = t(i,nkmb) - dt_dr * drml
3899 t(i,k+1) = (h(i,k+1) * t(i,k+1) + detrain(i) * tdown) * ih
3900 sdown = (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3901! The following two expressions updating S(nkmb) are mathematically identical.
3902! S(i,nkmb) = (h(i,nkmb) * S(i,nkmb) - detrain(i) * Sdown) / &
3903! (h(i,nkmb) - detrain(i))
3904 s(i,nkmb) = s(i,nkmb) - ds_dr * drml
3905 s(i,k+1) = (h(i,k+1) * s(i,k+1) + detrain(i) * sdown) * ih
3906
3907 d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3908 d_ea(i,nkmb) = d_ea(i,nkmb) - detrain(i)
3909
3910 h(i,k+1) = h(i,k+1) + detrain(i)
3911 h(i,nkmb) = h(i,nkmb) - detrain(i)
3912
3913 if (allocated(cs%diag_PE_detrain)) then
3914 if (cs%nonBous_energetics) then
3915 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + nb_grz_h2_2dt * detrain(i) * dspv0_drcv * &
3916 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3917 else
3918 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * dr0_drcv * &
3919 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3920 endif
3921 endif
3922 endif
3923 endif ! (RcvTgt(k) <= Rcv(i,nkmb))
3924 endif ! splittable_BL
3925 enddo ; enddo ! i & k loops
3926
3927! The numerical behavior of the buffer layer is dramatically improved
3928! if it is always at least a small fraction (say 10%) of the thickness
3929! of the mixed layer. As the physical distinction between the mixed
3930! and buffer layers is vague anyway, this seems hard to argue against.
3931 do i=is,ie
3932 if (h(i,nkmb) < 0.1*h(i,0)) then
3933 h_ent = 0.1*h(i,0) - h(i,nkmb)
3934 ih = 10.0/h(i,0)
3935 t(i,nkmb) = (h(i,nkmb)*t(i,nkmb) + h_ent*t(i,0)) * ih
3936 s(i,nkmb) = (h(i,nkmb)*s(i,nkmb) + h_ent*s(i,0)) * ih
3937
3938 d_ea(i,1) = d_ea(i,1) - h_ent
3939 d_ea(i,nkmb) = d_ea(i,nkmb) + h_ent
3940
3941 h(i,0) = h(i,0) - h_ent
3942 h(i,nkmb) = h(i,nkmb) + h_ent
3943 endif
3944 enddo
3945
3946end subroutine mixedlayer_detrain_1
3947
3948!> This subroutine initializes the MOM bulk mixed layer module.
3949subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
3950 type(time_type), target, intent(in) :: time !< The model's clock with the current time.
3951 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
3952 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
3953 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
3954 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
3955 !! parameters.
3956 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
3957 !! output.
3958 type(bulkmixedlayer_cs), intent(inout) :: cs !< Bulk mixed layer control structure
3959
3960 ! This include declares and sets the variable "version".
3961# include "version_variable.h"
3962 character(len=40) :: mdl = "MOM_mixed_layer" ! This module's name.
3963 real :: omega_frac_dflt ! The default value for ML_OMEGA_FRAC [nondim]
3964 real :: ustar_min_dflt ! The default value for BML_USTAR_MIN [Z T-1 ~> m s-1]
3965 real :: hmix_min_z ! HMIX_MIN in units of vertical extent [Z ~> m], used to set other defaults
3966 integer :: isd, ied, jsd, jed
3967 logical :: use_temperature, use_omega
3968 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3969
3970 cs%initialized = .true.
3971 cs%diag => diag
3972 cs%Time => time
3973
3974 if (gv%nkml < 1) return
3975
3976! Set default, read and log parameters
3977 call log_version(param_file, mdl, version, "")
3978
3979 cs%nkml = gv%nkml
3980 call log_param(param_file, mdl, "NKML", cs%nkml, &
3981 "The number of sublayers within the mixed layer if "//&
3982 "BULKMIXEDLAYER is true.", units="nondim", default=2)
3983 cs%nkbl = gv%nk_rho_varies - gv%nkml
3984 call log_param(param_file, mdl, "NKBL", cs%nkbl, &
3985 "The number of variable density buffer layers if "//&
3986 "BULKMIXEDLAYER is true.", units="nondim", default=2)
3987 call get_param(param_file, mdl, "MSTAR", cs%mstar, &
3988 "The ratio of the friction velocity cubed to the TKE "//&
3989 "input to the mixed layer.", units="nondim", default=1.2)
3990 call get_param(param_file, mdl, "NSTAR", cs%nstar, &
3991 "The portion of the buoyant potential energy imparted by "//&
3992 "surface fluxes that is available to drive entrainment "//&
3993 "at the base of mixed layer when that energy is positive.", &
3994 units="nondim", default=0.15)
3995 call get_param(param_file, mdl, "BULK_RI_ML", cs%bulk_Ri_ML, &
3996 "The efficiency with which mean kinetic energy released "//&
3997 "by mechanically forced entrainment of the mixed layer "//&
3998 "is converted to turbulent kinetic energy.", &
3999 units="nondim", fail_if_missing=.true., scale=us%L_to_Z**2)
4000 call get_param(param_file, mdl, "ABSORB_ALL_SW", cs%absorb_all_sw, &
4001 "If true, all shortwave radiation is absorbed by the "//&
4002 "ocean, instead of passing through to the bottom mud.", &
4003 default=.false.)
4004 call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
4005 "TKE_DECAY relates the vertical rate of decay of the "//&
4006 "TKE available for mechanical entrainment to the natural "//&
4007 "Ekman depth.", units="nondim", default=2.5)
4008 call get_param(param_file, mdl, "NSTAR2", cs%nstar2, &
4009 "The portion of any potential energy released by "//&
4010 "convective adjustment that is available to drive "//&
4011 "entrainment at the base of mixed layer. By default NSTAR2=NSTAR.", &
4012 units="nondim", default=cs%nstar)
4013 call get_param(param_file, mdl, "BULK_RI_CONVECTIVE", cs%bulk_Ri_convective, &
4014 "The efficiency with which convectively released mean "//&
4015 "kinetic energy is converted to turbulent kinetic "//&
4016 "energy. By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
4017 units="nondim", default=us%Z_to_L**2*cs%bulk_Ri_ML, scale=us%L_to_Z**2)
4018 call get_param(param_file, mdl, 'VON_KARMAN_CONST', cs%vonKar, &
4019 'The value the von Karman constant as used for mixed layer viscosity.', &
4020 units='nondim', default=0.41)
4021 call get_param(param_file, mdl, "HMIX_MIN", hmix_min_z, &
4022 "The minimum mixed layer depth if the mixed layer depth "//&
4023 "is determined dynamically.", units="m", default=0.0, scale=us%m_to_Z)
4024 cs%Hmix_min = gv%m_to_H * (us%Z_to_m * hmix_min_z)
4025 call get_param(param_file, mdl, "MECH_TKE_FLOOR", cs%mech_TKE_floor, &
4026 "A tiny floor on the amount of turbulent kinetic energy that is used when "//&
4027 "the mixed layer does not yet contain HMIX_MIN fluid. The default is so "//&
4028 "small that its actual value is irrelevant, so long as it is greater than 0.", &
4029 units="m3 s-2", default=1.0e-150, scale=gv%m_to_H*us%m_s_to_L_T**2*us%L_to_Z**2, &
4030 do_not_log=(hmix_min_z<=0.0))
4031
4032 call get_param(param_file, mdl, "LIMIT_BUFFER_DETRAIN", cs%limit_det, &
4033 "If true, limit the detrainment from the buffer layers "//&
4034 "to not be too different from the neighbors.", default=.false.)
4035 call get_param(param_file, mdl, "ALLOWED_DETRAIN_TEMP_CHG", cs%Allowed_T_chg, &
4036 "The amount by which temperature is allowed to exceed previous values "//&
4037 "during detrainment.", units="K", default=0.5, scale=us%degC_to_C)
4038 call get_param(param_file, mdl, "ALLOWED_DETRAIN_SALT_CHG", cs%Allowed_S_chg, &
4039 "The amount by which salinity is allowed to exceed previous values "//&
4040 "during detrainment.", units="ppt", default=0.1, scale=us%ppt_to_S)
4041 call get_param(param_file, mdl, "ML_DT_DS_WEIGHT", cs%dT_dS_wt, &
4042 "When forced to extrapolate T & S to match the layer "//&
4043 "densities, this factor (in deg C / PSU) is combined "//&
4044 "with the derivatives of density with T & S to determine "//&
4045 "what direction is orthogonal to density contours. It "//&
4046 "should be a typical value of (dR/dS) / (dR/dT) in oceanic profiles.", &
4047 units="degC ppt-1", default=6.0, scale=us%degC_to_C*us%S_to_ppt)
4048 call get_param(param_file, mdl, "BUFFER_LAYER_EXTRAP_LIMIT", cs%BL_extrap_lim, &
4049 "A limit on the density range over which extrapolation "//&
4050 "can occur when detraining from the buffer layers, "//&
4051 "relative to the density range within the mixed and "//&
4052 "buffer layers, when the detrainment is going into the "//&
4053 "lightest interior layer, nondimensional, or a negative "//&
4054 "value not to apply this limit.", units="nondim", default=-1.0)
4055 call get_param(param_file, mdl, "BUFFER_LAYER_HMIN_THICK", cs%Hbuffer_min, &
4056 "The minimum buffer layer thickness when the mixed layer is very thick.", &
4057 units="m", default=5.0, scale=gv%m_to_H)
4058 call get_param(param_file, mdl, "BUFFER_LAYER_HMIN_REL", cs%Hbuffer_rel_min, &
4059 "The minimum buffer layer thickness relative to the combined mixed "//&
4060 "land buffer ayer thicknesses when they are thin.", &
4061 units="nondim", default=0.1/cs%nkbl)
4062 if (cs%nkbl==1) then
4063 call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", cs%BL_detrain_time, &
4064 "A timescale that characterizes buffer layer detrainment events.", &
4065 units="s", default=86400.0*30.0, scale=us%s_to_T)
4066 else
4067 call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", cs%BL_detrain_time, &
4068 "A timescale that characterizes buffer layer detrainment events.", &
4069 units="s", default=4.0*3600.0, scale=us%s_to_T)
4070 endif
4071 call get_param(param_file, mdl, "BUFFER_SPLIT_RHO_TOL", cs%BL_split_rho_tol, &
4072 "The fractional tolerance for matching layer target densities when splitting "//&
4073 "layers to deal with massive interior layers that are lighter than one of the "//&
4074 "mixed or buffer layers.", units="nondim", default=0.1)
4075
4076 call get_param(param_file, mdl, "DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
4077 "The surface fluxes are scaled away when the total ocean "//&
4078 "depth is less than DEPTH_LIMIT_FLUXES.", &
4079 units="m", default=0.1*us%Z_to_m*hmix_min_z, scale=gv%m_to_H)
4080 call get_param(param_file, mdl, "OMEGA", cs%omega, &
4081 "The rotation rate of the earth.", &
4082 default=7.2921e-5, units="s-1", scale=us%T_to_s)
4083 call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
4084 "If true, use the absolute rotation rate instead of the "//&
4085 "vertical component of rotation when setting the decay "//&
4086 "scale for turbulence.", default=.false., do_not_log=.true.)
4087 omega_frac_dflt = 0.0
4088 if (use_omega) then
4089 call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
4090 omega_frac_dflt = 1.0
4091 endif
4092 call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
4093 "When setting the decay scale for turbulence, use this "//&
4094 "fraction of the absolute rotation rate blended with the "//&
4095 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
4096 units="nondim", default=omega_frac_dflt)
4097 call get_param(param_file, mdl, "ML_RESORT", cs%ML_resort, &
4098 "If true, resort the topmost layers by potential density "//&
4099 "before the mixed layer calculations.", default=.false.)
4100 if (cs%ML_resort) &
4101 call get_param(param_file, mdl, "ML_PRESORT_NK_CONV_ADJ", cs%ML_presort_nz_conv_adj, &
4102 "Convectively mix the first ML_PRESORT_NK_CONV_ADJ "//&
4103 "layers before sorting when ML_RESORT is true.", &
4104 units="nondim", default=0, fail_if_missing=.true.) ! Fail added by AJA.
4105 ! This gives a minimum decay scale that is typically much less than Angstrom.
4106 ustar_min_dflt = 2e-4*cs%omega*(gv%Angstrom_Z + gv%dZ_subroundoff)
4107 call get_param(param_file, mdl, "BML_USTAR_MIN", cs%ustar_min, &
4108 "The minimum value of ustar that should be used by the "//&
4109 "bulk mixed layer model in setting vertical TKE decay "//&
4110 "scales. This must be greater than 0.", &
4111 units="m s-1", default=us%Z_to_m*us%s_to_T*ustar_min_dflt, scale=us%m_to_Z*us%T_to_s)
4112 if (cs%ustar_min<=0.0) call mom_error(fatal, "BML_USTAR_MIN must be positive.")
4113
4114 call get_param(param_file, mdl, "BML_NONBOUSINESQ", cs%nonBous_energetics, &
4115 "If true, use non-Boussinesq expressions for the energetic calculations "//&
4116 "used in the bulk mixed layer calculations.", &
4117 default=.not.(gv%Boussinesq.or.gv%semi_Boussinesq))
4118
4119 call get_param(param_file, mdl, "RESOLVE_EKMAN", cs%Resolve_Ekman, &
4120 "If true, the NKML>1 layers in the mixed layer are "//&
4121 "chosen to optimally represent the impact of the Ekman "//&
4122 "transport on the mixed layer TKE budget. Otherwise, "//&
4123 "the sublayers are distributed uniformly through the "//&
4124 "mixed layer.", default=.false.)
4125 call get_param(param_file, mdl, "CORRECT_ABSORPTION_DEPTH", cs%correct_absorption, &
4126 "If true, the average depth at which penetrating shortwave "//&
4127 "radiation is absorbed is adjusted to match the average "//&
4128 "heating depth of an exponential profile by moving some "//&
4129 "of the heating upward in the water column.", default=.false.)
4130 call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
4131 "If true, apply additional mixing wherever there is "//&
4132 "runoff, so that it is mixed down to RIVERMIX_DEPTH, "//&
4133 "if the ocean is that deep.", default=.false.)
4134 if (cs%do_rivermix) &
4135 call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
4136 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
4137 "defined.", units="m", default=0.0, scale=gv%m_to_H)
4138 call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
4139 "If true, use the fluxes%runoff_Hflx field to set the "//&
4140 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
4141 default=.false.)
4142 call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
4143 "If true, use the fluxes%calving_Hflx field to set the "//&
4144 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
4145 default=.false.)
4146 call get_param(param_file, mdl, "BULKML_CONV_MOMENTUM_BUG", cs%convect_mom_bug, &
4147 "If true, use code with a bug that causes a loss of momentum conservation "//&
4148 "during mixedlayer convection.", default=.false.)
4149
4150 cs%id_ML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, &
4151 time, 'Surface mixed layer depth', 'm', conversion=gv%H_to_m)
4152 cs%id_TKE_wind = register_diag_field('ocean_model', 'TKE_wind', diag%axesT1, &
4153 time, 'Wind-stirring source of mixed layer TKE', &
4154 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4155 cs%id_TKE_RiBulk = register_diag_field('ocean_model', 'TKE_RiBulk', diag%axesT1, &
4156 time, 'Mean kinetic energy source of mixed layer TKE', &
4157 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4158 cs%id_TKE_conv = register_diag_field('ocean_model', 'TKE_conv', diag%axesT1, &
4159 time, 'Convective source of mixed layer TKE', &
4160 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4161 cs%id_TKE_pen_SW = register_diag_field('ocean_model', 'TKE_pen_SW', diag%axesT1, &
4162 time, 'TKE consumed by mixing penetrative shortwave radation through the mixed layer', &
4163 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4164 cs%id_TKE_mixing = register_diag_field('ocean_model', 'TKE_mixing', diag%axesT1, &
4165 time, 'TKE consumed by mixing that deepens the mixed layer', &
4166 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4167 cs%id_TKE_mech_decay = register_diag_field('ocean_model', 'TKE_mech_decay', diag%axesT1, &
4168 time, 'Mechanical energy decay sink of mixed layer TKE', &
4169 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4170 cs%id_TKE_conv_decay = register_diag_field('ocean_model', 'TKE_conv_decay', diag%axesT1, &
4171 time, 'Convective energy decay sink of mixed layer TKE', &
4172 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4173 cs%id_TKE_conv_s2 = register_diag_field('ocean_model', 'TKE_conv_s2', diag%axesT1, &
4174 time, 'Spurious source of mixed layer TKE from sigma2', &
4175 'm3 s-3', conversion=gv%H_to_m*(us%Z_to_m**2)*(us%s_to_T**3))
4176 cs%id_PE_detrain = register_diag_field('ocean_model', 'PE_detrain', diag%axesT1, &
4177 time, 'Spurious source of potential energy from mixed layer detrainment', &
4178 'W m-2', conversion=us%RZ3_T3_to_W_m2)
4179 cs%id_PE_detrain2 = register_diag_field('ocean_model', 'PE_detrain2', diag%axesT1, &
4180 time, 'Spurious source of potential energy from mixed layer only detrainment', &
4181 'W m-2', conversion=us%RZ3_T3_to_W_m2)
4182 cs%id_h_mismatch = register_diag_field('ocean_model', 'h_miss_ML', diag%axesT1, &
4183 time, 'Summed absolute mismatch in entrainment terms', 'm', conversion=gv%H_to_m)
4184 cs%id_Hsfc_used = register_diag_field('ocean_model', 'Hs_used', diag%axesT1, &
4185 time, 'Surface region thickness that is used', 'm', conversion=gv%H_to_m)
4186 cs%id_Hsfc_max = register_diag_field('ocean_model', 'Hs_max', diag%axesT1, &
4187 time, 'Maximum surface region thickness', 'm', conversion=gv%H_to_m)
4188 cs%id_Hsfc_min = register_diag_field('ocean_model', 'Hs_min', diag%axesT1, &
4189 time, 'Minimum surface region thickness', 'm', conversion=gv%H_to_m)
4190 !CS%lim_det_dH_sfc = 0.5 ; CS%lim_det_dH_bathy = 0.2 ! Technically these should not get used if limit_det is false?
4191 if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
4192 call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_SFC", cs%lim_det_dH_sfc, &
4193 "The fractional limit in the change between grid points "//&
4194 "of the surface region (mixed & buffer layer) thickness.", &
4195 units="nondim", default=0.5)
4196 call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_BATHY", cs%lim_det_dH_bathy, &
4197 "The fraction of the total depth by which the thickness "//&
4198 "of the surface region (mixed & buffer layer) is allowed "//&
4199 "to change between grid points.", units="nondim", default=0.2)
4200 endif
4201
4202 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
4203 "If true, temperature and salinity are used as state "//&
4204 "variables.", default=.true.)
4205 cs%nsw = 0
4206 if (use_temperature) then
4207 call get_param(param_file, mdl, "PEN_SW_NBANDS", cs%nsw, default=1)
4208 endif
4209
4210
4211 if (max(cs%id_TKE_wind, cs%id_TKE_RiBulk, cs%id_TKE_conv, cs%id_TKE_mixing, &
4212 cs%id_TKE_pen_SW, cs%id_TKE_mech_decay, cs%id_TKE_conv_decay) > 0) then
4213 call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
4214 call safe_alloc_alloc(cs%diag_TKE_RiBulk, isd, ied, jsd, jed)
4215 call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
4216 call safe_alloc_alloc(cs%diag_TKE_pen_SW, isd, ied, jsd, jed)
4217 call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
4218 call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
4219 call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
4220 call safe_alloc_alloc(cs%diag_TKE_conv_s2, isd, ied, jsd, jed)
4221
4222 cs%TKE_diagnostics = .true.
4223 endif
4224 if (cs%id_PE_detrain > 0) call safe_alloc_alloc(cs%diag_PE_detrain, isd, ied, jsd, jed)
4225 if (cs%id_PE_detrain2 > 0) call safe_alloc_alloc(cs%diag_PE_detrain2, isd, ied, jsd, jed)
4226 if (cs%id_ML_depth > 0) call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
4227
4228 if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) &
4229 id_clock_pass = cpu_clock_id('(Ocean mixed layer halo updates)', grain=clock_routine)
4230
4231end subroutine bulkmixedlayer_init
4232
4233!> This subroutine returns an approximation to the integral
4234!! R = exp(-L*(H+E)) integral(LH to L(H+E)) L/(1-(1+x)exp(-x)) dx.
4235!! The approximation to the integrand is good to within -2% at x~.3
4236!! and +25% at x~3.5, but the exponential deemphasizes the importance of
4237!! large x. When L=0, EF4 returns E/((Ht+E)*Ht).
4238function ef4(Ht, En, I_L, dR_de)
4239 real, intent(in) :: ht !< Total thickness [H ~> m or kg m-2].
4240 real, intent(in) :: en !< Entrainment [H ~> m or kg m-2].
4241 real, intent(in) :: i_l !< The e-folding scale [H-1 ~> m-1 or m2 kg-1]
4242 real, optional, intent(inout) :: dr_de !< The partial derivative of the result R with E [H-2 ~> m-2 or m4 kg-2].
4243 real :: ef4 !< The integral [H-1 ~> m-1 or m2 kg-1].
4244
4245 ! Local variables
4246 real :: exp_lhpe ! A nondimensional exponential decay [nondim].
4247 real :: i_hpe ! An inverse thickness plus entrainment [H-1 ~> m-1 or m2 kg-1].
4248 real :: res ! The result of the integral above [H-1 ~> m-1 or m2 kg-1].
4249
4250 exp_lhpe = exp(-i_l*(en+ht))
4251 i_hpe = 1.0/(ht+en)
4252 res = exp_lhpe * (en*i_hpe/ht - 0.5*i_l*log(ht*i_hpe) + 0.5*i_l*i_l*en)
4253 if (PRESENT(dr_de)) &
4254 dr_de = -i_l*res + exp_lhpe*(i_hpe*i_hpe + 0.5*i_l*i_hpe + 0.5*i_l*i_l)
4255 ef4 = res
4256
4257end function ef4
4258
4259!> \namespace mom_bulk_mixed_layer
4260!!
4261!! By Robert Hallberg, 1997 - 2005.
4262!!
4263!! This file contains the subroutine (bulkmixedlayer) that
4264!! implements a Kraus-Turner-like bulk mixed layer, based on the work
4265!! of various people, as described in the review paper by \cite niiler1977,
4266!! with particular attention to the form proposed by \cite Oberhuber1993a,
4267!! with an extension to a refined bulk mixed layer as described in
4268!! Hallberg (\cite muller2003). The physical processes portrayed in
4269!! this subroutine include convective adjustment and mixed layer entrainment
4270!! and detrainment. Penetrating shortwave radiation and an exponential decay
4271!! of TKE fluxes are also supported by this subroutine. Several constants
4272!! can alternately be set to give a traditional Kraus-Turner mixed
4273!! layer scheme, although that is not the preferred option. The
4274!! physical processes and arguments are described in detail in \ref BML.
4275
4276end module mom_bulk_mixed_layer