src/parameterizations/vertical/MOM_entrain_diffusive.F90

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

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Diapycnal mixing and advection in isopycnal mode
6		module MOM_entrain_diffusive
7
8		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
9		use MOM_diag_mediator, only : diag_ctrl, time_type
10		use MOM_EOS, only : calculate_density, calculate_density_derivs
11		use MOM_EOS, only : calculate_specific_vol_derivs, EOS_domain
12		use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
13		use MOM_file_parser, only : get_param, log_version, param_file_type
14		use MOM_forcing_type, only : forcing
15		use MOM_grid, only : ocean_grid_type
16		use MOM_unit_scaling, only : unit_scale_type
17		use MOM_variables, only : thermo_var_ptrs
18		use MOM_verticalGrid, only : verticalGrid_type
19
20		implicit none ; private
21
22		#include <MOM_memory.h>
23
24		public entrainment_diffusive, entrain_diffusive_init
25
26		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
27		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
28		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
29		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
30
31		!> The control structure holding parametes for the MOM_entrain_diffusive module
32		type, public :: entrain_diffusive_CS ; private
33		logical :: initialized = .false. !< True if this control structure has been initialized.
34		logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with
35		!! GV%nk_rho_varies variable density mixed & buffer layers.
36		integer :: max_ent_it !< The maximum number of iterations that may be used to
37		!! calculate the diapycnal entrainment.
38		real :: Tolerance_Ent !< The tolerance with which to solve for entrainment values
39		!! [H ~> m or kg m-2].
40		real :: max_Ent !< A large ceiling on the maximum permitted amount of entrainment
41		!! across each interface between the mixed and buffer layers within
42		!! a timestep [H ~> m or kg m-2].
43		real :: Rho_sig_off !< The offset between potential density and a sigma value [R ~> kg m-3]
44		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
45		!! regulate the timing of diagnostic output.
46		integer :: id_Kd = -1 !< Diagnostic ID for diffusivity
47		integer :: id_diff_work = -1 !< Diagnostic ID for mixing work
48		end type entrain_diffusive_CS
49
50		contains
51
52		!> This subroutine calculates ea and eb, the rates at which a layer entrains
53		!! from the layers above and below. The entrainment rates are proportional to
54		!! the buoyancy flux in a layer and inversely proportional to the density
55		!! differences between layers. The scheme that is used here is described in
56		!! detail in Hallberg, Mon. Wea. Rev. 2000.
57	1	subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
58	1	kb_out, Kd_Lay, Kd_int)
59		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
60		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
61		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
62		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
63		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
64		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
65		!! thermodynamic fields. Absent fields have NULL
66		!! ptrs.
67		type(forcing), intent(in) :: fluxes !< A structure of surface fluxes that may
68		!! be used.
69		real, intent(in) :: dt !< The time increment [T ~> s].
70		type(entrain_diffusive_CS), intent(in) :: CS !< The control structure returned by a previous
71		!! call to entrain_diffusive_init.
72		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
73		intent(out) :: ea !< The amount of fluid entrained from the layer
74		!! above within this time step [H ~> m or kg m-2].
75		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
76		intent(out) :: eb !< The amount of fluid entrained from the layer
77		!! below within this time step [H ~> m or kg m-2].
78		integer, dimension(SZI_(G),SZJ_(G)), &
79		intent(inout) :: kb_out !< The index of the lightest layer denser than
80		!! the buffer layer.
81		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
82		intent(in) :: Kd_Lay !< The diapycnal diffusivity of layers
83		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
84		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
85		intent(in) :: Kd_int !< The diapycnal diffusivity of interfaces
86		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
87
88		! This subroutine calculates ea and eb, the rates at which a layer entrains
89		! from the layers above and below. The entrainment rates are proportional to
90		! the buoyancy flux in a layer and inversely proportional to the density
91		! differences between layers. The scheme that is used here is described in
92		! detail in Hallberg, Mon. Wea. Rev. 2000.
93
94		real, dimension(SZI_(G),SZK_(GV)) :: &
95	2	dtKd ! The layer diapycnal diffusivity times the time step [H2 ~> m2 or kg2 m-4].
96		real, dimension(SZI_(G),SZK_(GV)+1) :: &
97	2	dtKd_int ! The diapycnal diffusivity at the interfaces times the time step [H2 ~> m2 or kg2 m-4]
98		real, dimension(SZI_(G),SZK_(GV)) :: &
99	2	F, & ! The density flux through a layer within a time step divided by the
100		! density difference across the interface below the layer [H ~> m or kg m-2].
101	2	maxF, & ! maxF is the maximum value of F that will not deplete all of the
102		! layers above or below a layer within a timestep [H ~> m or kg m-2].
103	2	minF, & ! minF is the minimum flux that should be expected in the absence of
104		! interactions between layers [H ~> m or kg m-2].
105	2	Fprev, &! The previous estimate of F [H ~> m or kg m-2].
106	2	dFdfm, &! The partial derivative of F with respect to changes in F of the
107		! neighboring layers. [nondim]
108	2	h_guess ! An estimate of the layer thicknesses after entrainment, but
109		! before the entrainments are adjusted to drive the layer
110		! densities toward their target values [H ~> m or kg m-2].
111		real, dimension(SZI_(G),SZK_(GV)+1) :: &
112	2	Ent_bl ! The average entrainment upward and downward across
113		! each interface around the buffer layers [H ~> m or kg m-2].
114		real, allocatable, dimension(:,:,:) :: &
115	1	Kd_eff, & ! The effective diffusivity that actually applies to each
116		! layer after the effects of boundary conditions are
117		! considered [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
118	1	diff_work ! The work actually done by diffusion across each
119		! interface [R Z3 T-3 ~> W m-2]. Sum vertically for the total work.
120
121		real :: hm, fm, fr ! Work variables with units of [H ~> m or kg m-2].
122		real :: fk ! A Work variable with units of [H2 ~> m2 or kg2 m-4]
123
124	2	real :: b1(SZI_(G)) ! A variable used by the tridiagonal solver [H ~> m or kg m-2]
125	2	real :: c1(SZI_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]
126
127		real, dimension(SZI_(G)) :: &
128	2	htot, & ! The total thickness above or below a layer [H ~> m or kg m-2].
129	2	Rcv, & ! Value of the coordinate variable (potential density)
130		! based on the simulated T and S and P_Ref [R ~> kg m-3].
131	2	pres, & ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].
132	2	eakb, & ! The entrainment from above by the layer below the buffer
133		! layer (i.e. layer kb) [H ~> m or kg m-2].
134	2	ea_kbp1, & ! The entrainment from above by layer kb+1 [H ~> m or kg m-2].
135	2	eb_kmb, & ! The entrainment from below by the deepest buffer layer [H ~> m or kg m-2].
136	2	dS_kb, & ! The reference potential density difference across the
137		! interface between the buffer layers and layer kb [R ~> kg m-3].
138	2	dS_anom_lim, &! The amount by which dS_kb is reduced when limits are
139		! applied [R ~> kg m-3].
140	2	I_dSkbp1, & ! The inverse of the potential density difference across the
141		! interface below layer kb [R-1 ~> m3 kg-1].
142	2	dtKd_kb, & ! The diapycnal diffusivity in layer kb times the time step
143		! [H2 ~> m2 or kg2 m-4].
144	2	maxF_correct, & ! An amount by which to correct maxF due to excessive
145		! surface heat loss [H ~> m or kg m-2].
146	2	zeros, & ! An array of all zeros. (Usually used with [H ~> m or kg m-2].)
147	2	max_eakb, & ! The maximum value of eakb that might be realized [H ~> m or kg m-2].
148	2	min_eakb, & ! The minimum value of eakb that might be realized [H ~> m or kg m-2].
149	2	err_max_eakb0, & ! The value of error returned by determine_Ea_kb when eakb = max_eakb
150		! and ea_kbp1 = 0 [H2 ~> m2 or kg2 m-4].
151	2	err_min_eakb0, & ! The value of error returned by determine_Ea_kb when eakb = min_eakb
152		! and ea_kbp1 = 0 [H2 ~> m2 or kg2 m-4].
153	2	err_eakb0, & ! A value of error returned by determine_Ea_kb [H2 ~> m2 or kg2 m-4].
154	2	F_kb, & ! The value of F in layer kb, or equivalently the entrainment
155		! from below by layer kb [H ~> m or kg m-2].
156	2	dFdfm_kb, & ! The partial derivative of F with fm [nondim]. See dFdfm.
157	2	maxF_kb, & ! The maximum value of F_kb that might be realized [H ~> m or kg m-2].
158	2	eakb_maxF, & ! The value of eakb that gives F_kb=maxF_kb [H ~> m or kg m-2].
159	2	F_kb_maxEnt ! The value of F_kb when eakb = max_eakb [H ~> m or kg m-2].
160		real, dimension(SZI_(G),SZK_(GV)) :: &
161	2	Sref, & ! The reference potential density of the mixed and buffer layers,
162		! and of the two lightest interior layers (kb and kb+1) copied
163		! into layers kmb+1 and kmb+2 [R ~> kg m-3].
164	2	h_bl ! The thicknesses of the mixed and buffer layers, and of the two
165		! lightest interior layers (kb and kb+1) copied into layers kmb+1
166		! and kmb+2 [H ~> m or kg m-2].
167
168		real, dimension(SZI_(G),SZK_(GV)) :: &
169	2	ds_dsp1, & ! The coordinate variable (sigma-2) difference across an
170		! interface divided by the difference across the interface
171		! below it. [nondim]
172	2	dsp1_ds, & ! The inverse coordinate variable (sigma-2) difference
173		! across an interface times the difference across the
174		! interface above it. [nondim]
175	2	I2p2dsp1_ds, & ! 1 / (2 + 2 * ds_k+1 / ds_k). [nondim]
176	2	grats ! 2*(2 + ds_k+1 / ds_k + ds_k / ds_k+1) =
177		! 4ds_Lay(1/ds_k + 1/ds_k+1). [nondim]
178
179		real :: dRho ! The change in locally referenced potential density between
180		! the layers above and below an interface [R ~> kg m-3]
181		real :: dSpV ! The change in locally referenced specific volume between
182		! the layers above and below an interface [R-1 ~> m3 kg-1]
183		real :: g_2dt ! 0.5 * G_Earth / dt, times unit conversion factors
184		! [Z3 H-2 T-3 or R2 Z3 H-2 T-3 ~> m s-3].
185		real, dimension(SZI_(G)) :: &
186	2	pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].
187	2	T_eos, S_eos, & ! The potential temperature and salinity at which to
188		! evaluate dRho_dT and dRho_dS [C ~> degC] and [S ~> ppt].
189	2	dRho_dT, & ! The partial derivative of potential density with temperature [R C-1 ~> kg m-3 degC-1]
190	2	dRho_dS, & ! The partial derivative of potential density with salinity [R S-1 ~> kg m-3 ppt-1]
191	2	dSpV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
192	2	dSpV_dS ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
193
194		real :: tolerance ! The tolerance within which E must be converged [H ~> m or kg m-2].
195		real :: Angstrom ! The minimum layer thickness [H ~> m or kg m-2].
196		real :: h_neglect ! A thickness that is so small it is usually lost
197		! in roundoff and can be neglected [H ~> m or kg m-2].
198		real :: F_cor ! A correction to the amount of F that is used to
199		! entrain from the layer above [H ~> m or kg m-2].
200		real :: Kd_here ! The effective diapycnal diffusivity times the timestep [H2 ~> m2 or kg2 m-4].
201		real :: h_avail ! The thickness that is available for entrainment [H ~> m or kg m-2].
202		real :: dS_kb_eff ! The value of dS_kb after limiting is taken into account [R ~> kg m-3].
203		real :: Rho_cor ! The depth-integrated potential density anomaly that
204		! needs to be corrected for [H R ~> kg m-2 or kg2 m-5].
205		real :: ea_cor ! The corrective adjustment to eakb [H ~> m or kg m-2].
206		real :: h1 ! The layer thickness after entrainment through the
207		! interface below is taken into account [H ~> m or kg m-2].
208		real :: Idt ! The inverse of the time step [Z H-1 T-1 ~> s-1 or m3 kg-1 s-1].
209
210		logical :: do_any
211		logical :: do_entrain_eakb ! True if buffer layer is entrained
212	3	logical :: do_i(SZI_(G)), did_i(SZI_(G)), reiterate
213		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
214		integer :: it, i, j, k, is, ie, js, je, nz, K2, kmb
215	2	integer :: kb(SZI_(G)) ! The value of kb in row j.
216		integer :: kb_min ! The minimum value of kb in the current j-row.
217		integer :: kb_min_act ! The minimum active value of kb in the current j-row.
218		integer :: is1, ie1 ! The minimum and maximum active values of i in the current j-row.
219	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
220	1	Angstrom = GV%Angstrom_H
221	1	h_neglect = GV%H_subroundoff
222
223	1	if (.not. CS%initialized) call MOM_error(FATAL, &
224	0	"MOM_entrain_diffusive: Module must be initialized before it is used.")
225
226	1	if ((.not.CS%bulkmixedlayer .and. .not.associated(fluxes%buoy)) .and. &
227		(associated(fluxes%lprec) .or. associated(fluxes%evap) .or. &
228		associated(fluxes%sens) .or. associated(fluxes%sw))) then
229	0	if (is_root_pe()) call MOM_error(NOTE, "Calculate_Entrainment: &
230		&The code to handle evaporation and precipitation without &
231	0	&a bulk mixed layer has not been implemented.")
232	0	if (is_root_pe()) call MOM_error(FATAL, &
233		"Either define BULKMIXEDLAYER in MOM_input or use fluxes%buoy &
234	0	&and a linear equation of state to drive the model.")
235		endif
236
237	1	tolerance = CS%Tolerance_Ent
238	1	kmb = GV%nk_rho_varies
239	1	K2 = max(kmb+1,2) ; kb_min = K2
240	1	if (.not. CS%bulkmixedlayer) then
241	0	kb(:) = 1
242		! These lines fill in values that are arbitrary, but needed because
243		! they are used to normalize the buoyancy flux in layer nz.
244	0	do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ; enddo
245		else
246	369	kb(:) = 0
247	361	do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ; enddo
248		endif
249
250	1	if (CS%id_diff_work > 0) allocate(diff_work(G%isd:G%ied,G%jsd:G%jed,nz+1))
251	1	if (CS%id_Kd > 0) allocate(Kd_eff(G%isd:G%ied,G%jsd:G%jed,nz))
252
253	1	if (associated(tv%eqn_of_state)) then
254	369	pres(:) = tv%P_Ref
255		else
256	0	pres(:) = 0.0
257		endif
258	1	EOSdom(:) = EOS_domain(G%HI)
259
260		!$OMP parallel do default(private) shared(is,ie,js,je,nz,Kd_Lay,G,GV,US,dt,CS,h,tv, &
261		!$OMP kmb,Angstrom,fluxes,K2,h_neglect,tolerance, &
262		!$OMP ea,eb,Kd_int,Kd_eff,EOSdom,diff_work,g_2dt, kb_out) &
263		!$OMP firstprivate(kb,ds_dsp1,dsp1_ds,pres,kb_min)
264	181	do j=js,je
265	64980	do i=is,ie ; kb(i) = 1 ; enddo
266
267	180	if (allocated(tv%SpV_avg)) then
268	0	do k=1,nz ; do i=is,ie
269	0	dtKd(i,k) = GV%RZ_to_H * (dt * Kd_lay(i,j,k)) / tv%SpV_avg(i,j,k)
270		enddo ; enddo
271	0	do i=is,ie
272	0	dtKd_int(i,1) = GV%RZ_to_H * (dt * Kd_int(i,j,1)) / tv%SpV_avg(i,j,1)
273	0	dtKd_int(i,nz+1) = GV%RZ_to_H * (dt * Kd_int(i,j,nz+1)) / tv%SpV_avg(i,j,nz)
274		enddo
275		! Use the mass-weighted average specific volume to translate thicknesses to verti distances.
276	0	do K=2,nz ; do i=is,ie
277		dtKd_int(i,K) = GV%RZ_to_H * (dt * Kd_int(i,j,K)) * &
278		( (h(i,j,k-1) + h(i,j,k) + 2.0*h_neglect) / &
279		((h(i,j,k-1)+h_neglect) * tv%SpV_avg(i,j,k-1) + &
280	0	(h(i,j,k)+h_neglect) * tv%SpV_avg(i,j,k)) )
281		enddo ; enddo
282		else
283	1429740	do k=1,nz ; do i=is,ie
284	1429560	dtKd(i,k) = GV%Z_to_H * (dt * Kd_lay(i,j,k))
285		enddo ; enddo
286	1494720	do K=1,nz+1 ; do i=is,ie
287	1494540	dtKd_int(i,K) = GV%Z_to_H * (dt * Kd_int(i,j,K))
288		enddo ; enddo
289		endif
290
291	64980	do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
292	64980	do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; enddo
293	64980	do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
294
295	180	if (GV%Boussinesq .or. GV%Semi_Boussinesq) then
296	1299780	do k=2,nz-1 ; do i=is,ie
297	1299600	ds_dsp1(i,k) = GV%g_prime(k) / GV%g_prime(k+1)
298		enddo ; enddo
299		else ! Use a mathematically equivalent form that avoids any dependency on RHO_0.
300	0	do k=2,nz-1 ; do i=is,ie
301	0	ds_dsp1(i,k) = (GV%Rlay(k) - GV%Rlay(k-1)) / (GV%Rlay(k+1) - GV%Rlay(k))
302		enddo ; enddo
303		endif
304
305	180	if (CS%bulkmixedlayer) then
306		! This subroutine determines the averaged entrainment across each
307		! interface and causes thin and relatively light interior layers to be
308		! entrained by the deepest buffer layer. This also determines kb.
309	180	call set_Ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)
310
311	64980	do i=is,ie
312	64980	dtKd_kb(i) = 0.0 ; if (kb(i) < nz) dtKd_kb(i) = dtKd(i,kb(i))
313		enddo
314		else
315	0	do i=is,ie ; Ent_bl(i,Kmb+1) = 0.0 ; enddo
316		endif
317
318	1299780	do k=2,nz-1 ; do i=is,ie
319	1296000	dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
320	1296000	I2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))
321	1299600	grats(i,k) = 2.0*(2.0+(dsp1_ds(i,k)+ds_dsp1(i,k)))
322		enddo ; enddo
323
324		! Determine the maximum flux, maxF, for each of the isopycnal layers.
325		! Also determine when the fluxes start entraining
326		! from various buffer or mixed layers, where appropriate.
327	180	if (CS%bulkmixedlayer) then
328	180	kb_min = nz
329	64980	do i=is,ie
330	64980	htot(i) = h(i,j,1) - Angstrom
331		enddo
332	195120	do k=2,kmb ; do i=is,ie
333	194940	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
334		enddo ; enddo
335	64980	do i=is,ie
336	64800	max_eakb(i) = MAX(Ent_bl(i,Kmb+1) + 0.5*htot(i), htot(i))
337	64800	I_dSkbp1(i) = 1.0 / (Sref(i,kmb+2) - Sref(i,kmb+1))
338	64980	zeros(i) = 0.0
339		enddo
340
341		! Find the maximum amount of entrainment from below that the top
342		! interior layer could exhibit (maxF(i,kb)), approximating that
343		! entrainment by (eakb*max(dS_kb/dSkbp1,0)). eakb is in the range
344		! from 0 to max_eakb.
345		call find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, zeros, max_eakb, kmb, &
346	180	is, ie, G, GV, CS, maxF_kb, eakb_maxF, do_i, F_kb_maxent)
347	64980	do i=is,ie ; if (kb(i) <= nz) then
348	46250	maxF(i,kb(i)) = MAX(0.0, maxF_kb(i), F_kb_maxent(i))
349	46250	if ((maxF_kb(i) > F_kb_maxent(i)) .and. (eakb_maxF(i) >= htot(i))) &
350	0	max_eakb(i) = eakb_maxF(i)
351		endif ; enddo
352
353	64980	do i=is,ie ; ea_kbp1(i) = 0.0 ; eakb(i) = max_eakb(i) ; enddo
354		call determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &
355		max_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, eakb, &
356	180	error=err_max_eakb0, F_kb=F_kb)
357
358		! The maximum value of F(kb) between htot and max_eakb determines
359		! what maxF(kb+1) should be.
360	64980	do i=is,ie ; min_eakb(i) = MIN(htot(i), max_eakb(i)) ; enddo
361		call find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_eakb, max_eakb, &
362	180	kmb, is, ie, G, GV, CS, F_kb_maxEnt, do_i_in=do_i)
363
364	64980	do i=is,ie
365	64800	do_entrain_eakb = .false.
366		! If error_max_eakb0 < 0, then buffer layers are always all entrained
367	64800	if (do_i(i)) then ; if (err_max_eakb0(i) < 0.0) then
368	0	do_entrain_eakb = .true.
369		endif ; endif
370
371	64980	if (do_entrain_eakb) then
372	0	eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
373		else
374	64800	eakb(i) = 0.0 ; min_eakb(i) = 0.0
375		endif
376		enddo
377
378		! Find the amount of entrainment of the buffer layers that would occur
379		! if there were no entrainment by the deeper interior layers. Also find
380		! how much entrainment of the deeper layers would occur.
381		call determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &
382		zeros, max_eakb, kmb, is, ie, do_i, G, GV, CS, min_eakb, &
383	180	error=err_min_eakb0, F_kb=F_kb, err_max_eakb0=err_max_eakb0)
384		! Error_min_eakb0 should be ~0 unless error_max_eakb0 < 0.
385	64980	do i=is,ie ; if ((kb(i)<nz) .and. (kb_min>kb(i))) kb_min = kb(i) ; enddo
386		else
387		! Without a bulk mixed layer, surface fluxes are applied in this
388		! subroutine. (Otherwise, they are handled in mixedlayer.)
389		! Initially the maximum flux in layer zero is given by the surface
390		! buoyancy flux. It will later be limited if the surface flux is
391		! too large. Here buoy is the surface buoyancy flux.
392	0	do i=is,ie
393	0	maxF(i,1) = 0.0
394	0	htot(i) = h(i,j,1) - Angstrom
395		enddo
396	0	if (associated(fluxes%buoy) .and. GV%Boussinesq) then
397	0	do i=is,ie
398	0	maxF(i,1) = GV%Z_to_H * (dt*fluxes%buoy(i,j)) / GV%g_prime(2)
399		enddo
400	0	elseif (associated(fluxes%buoy)) then
401	0	do i=is,ie
402		maxF(i,1) = (GV%RZ_to_H * 0.5(GV%Rlay(1) + GV%Rlay(2)) (dt*fluxes%buoy(i,j))) / &
403	0	GV%g_prime(2)
404		enddo
405		endif
406		endif
407
408		! The following code calculates the maximum flux, maxF, for the interior
409		! layers.
410	612436	do k=kb_min,nz-1 ; do i=is,ie
411	612256	if ((k == kb(i)+1) .and. CS%bulkmixedlayer) then
412	46250	maxF(i,k) = ds_dsp1(i,k)*(F_kb_maxEnt(i) + htot(i))
413	46250	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
414	564310	elseif (k > kb(i)) then
415	272689	maxF(i,k) = ds_dsp1(i,k)*(maxF(i,k-1) + htot(i))
416	272689	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
417		endif
418		enddo ; enddo
419	64980	do i=is,ie
420	64800	maxF(i,nz) = 0.0
421	64800	if (.not.CS%bulkmixedlayer) then
422	0	maxF_correct(i) = MAX(0.0, -(maxF(i,nz-1) + htot(i)))
423		endif
424	64980	htot(i) = h(i,j,nz) - Angstrom
425		enddo
426	180	if (.not.CS%bulkmixedlayer) then
427	0	do_any = .false. ; do i=is,ie ; if (maxF_correct(i) > 0.0) do_any = .true. ; enddo
428	0	if (do_any) then
429	0	do k=nz-1,1,-1 ; do i=is,ie
430	0	maxF(i,k) = maxF(i,k) + maxF_correct(i)
431	0	maxF_correct(i) = maxF_correct(i) * dsp1_ds(i,k)
432		enddo ; enddo
433		endif
434		endif
435	612436	do k=nz-1,kb_min,-1 ; do i=is,ie ; if (do_i(i)) then
436	365388	if (k >= kb(i)) then
437	365189	maxF(i,k) = MIN(maxF(i,k),dsp1_ds(i,k+1)*maxF(i,k+1) + htot(i))
438	365189	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
439		endif
440	365388	if (k == kb(i)) then
441	46250	if ((maxF(i,k) < F_kb(i)) .or. (maxF(i,k) < maxF_kb(i)) &
442		.and. (eakb_maxF(i) <= max_eakb(i))) then
443		! In this case, too much was being entrained by the topmost interior
444		! layer, even with the minimum initial estimate. The buffer layer
445		! will always entrain the maximum amount.
446	948	F_kb(i) = maxF(i,k)
447	948	if ((F_kb(i) <= maxF_kb(i)) .and. (eakb_maxF(i) <= max_eakb(i))) then
448	948	eakb(i) = eakb_maxF(i)
449		else
450	0	eakb(i) = max_eakb(i)
451		endif
452		call F_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &
453	948	G, GV, CS, eakb, Angstrom)
454	948	if ((eakb(i) < max_eakb(i)) .or. (eakb(i) < min_eakb(i))) then
455		call determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, zeros, &
456		eakb, eakb, kmb, i, i, do_i, G, GV, CS, eakb, &
457	948	error=err_eakb0)
458	948	if (eakb(i) < max_eakb(i)) then
459	948	max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)
460		endif
461	948	if (eakb(i) < min_eakb(i)) then
462	765	min_eakb(i) = eakb(i) ; err_min_eakb0(i) = err_eakb0(i)
463		endif
464		endif
465		endif
466		endif
467		endif ; enddo ; enddo
468	180	if (.not.CS%bulkmixedlayer) then
469	0	do i=is,ie
470	0	maxF(i,1) = MIN(maxF(i,1),dsp1_ds(i,2)*maxF(i,2) + htot(i))
471		enddo
472		endif
473
474		! The following code provides an initial estimate of the flux in
475		! each layer, F. The initial guess for the layer diffusive flux is
476		! the smaller of a forward discretization or the maximum diffusive
477		! flux starting from zero thickness in one time step without
478		! considering adjacent layers.
479	64980	do i=is,ie
480	64800	F(i,1) = maxF(i,1)
481	64980	F(i,nz) = maxF(i,nz) ; minF(i,nz) = 0.0
482		enddo
483	3240	do k=nz-1,K2,-1
484	1104840	do i=is,ie
485	1104660	if ((k==kb(i)) .and. (do_i(i))) then
486	46250	eakb(i) = min_eakb(i)
487	46250	minF(i,k) = 0.0
488	1055350	elseif ((k>kb(i)) .and. (do_i(i))) then
489		! Here the layer flux is estimated, assuming no entrainment from
490		! the surrounding layers. The estimate is a forward (steady) flux,
491		! limited by the maximum flux for a layer starting with zero
492		! thickness. This is often a good guess and leads to few iterations.
493	318939	hm = h(i,j,k) + h_neglect
494		! Note: Tried sqrt((0.5ds_dsp1(i,k))dtKd(i,k)) for the second limit,
495		! but it usually doesn't work as well.
496		F(i,k) = MIN(maxF(i,k), sqrt(ds_dsp1(i,k)*dtKd(i,k)), &
497	318939	0.5(ds_dsp1(i,k)+1.0) (dtKd(i,k) / hm))
498
499		! Calculate the minimum flux that can be expected if there is no entrainment
500		! from the neighboring layers. The 0.9 is used to give used to give a 10%
501		! known error tolerance.
502	318939	fk = dtKd(i,k) * grats(i,k)
503		minF(i,k) = MIN(maxF(i,k), &
504	318939	0.9(I2p2dsp1_ds(i,k) fk / (hm + sqrt(hm*hm + fk))))
505	318939	if (k==kb(i)) minF(i,k) = 0.0 ! BACKWARD COMPATIBILITY - DELETE LATER?
506		else
507	736411	F(i,k) = 0.0
508	736411	minF(i,k) = 0.0
509		endif
510		enddo ! end of i loop
511		enddo ! end of k loop
512
513		! This is where the fluxes are actually calculated.
514
515	180	is1 = ie+1 ; ie1 = is-1
516	4578	do i=is,ie ; if (do_i(i)) then ; is1 = i ; exit ; endif ; enddo
517	4578	do i=ie,is,-1 ; if (do_i(i)) then ; ie1 = i ; exit ; endif ; enddo
518
519	180	if (CS%bulkmixedlayer) then
520	180	kb_min_act = nz
521	64980	do i=is,ie
522	64980	if (do_i(i) .and. (kb(i) < kb_min_act)) kb_min_act = kb(i)
523		enddo
524		! Solve for the entrainment rate from above in the topmost interior
525		! layer, eakb, such that
526		! eakbdS_implicit = dtKd*dS_layer_implicit / h_implicit.
527	56184	do i=is1,ie1
528	56004	ea_kbp1(i) = 0.0
529	56004	if (do_i(i) .and. (kb(i) < nz)) &
530	46430	ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*F(i,kb(i)+1)
531		enddo
532		call determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, min_eakb, &
533		max_eakb, kmb, is1, ie1, do_i, G, GV, CS, eakb, F_kb=F_kb, &
534		err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
535	180	dFdfm_kb=dFdfm_kb)
536		else
537	0	kb_min_act = kb_min
538		endif
539
540	180	do it=0,CS%max_ent_it-1
541	56184	do i=is1,ie1 ; if (do_i(i)) then
542	46250	if (.not.CS%bulkmixedlayer) F(i,1) = MIN(F(i,1),maxF(i,1))
543	46250	b1(i) = 1.0
544		endif ; enddo ! end of i loop
545
546		! F_kb has already been found for this iteration, either above or at
547		! the end of the code for the previous iteration.
548	493018	do k=kb_min_act,nz-1 ; do i=is1,ie1 ; if (do_i(i) .and. (k>=kb(i))) then
549		! Calculate the flux in layer k.
550	365189	if (CS%bulkmixedlayer .and. (k==kb(i))) then
551	46250	F(i,k) = F_kb(i)
552	46250	dFdfm(i,k) = dFdfm_kb(i)
553		else ! k > kb(i)
554	318939	Fprev(i,k) = F(i,k)
555	318939	fm = (F(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*F(i,k+1)
556	318939	fk = grats(i,k)*dtKd(i,k)
557	318939	fr = sqrt(fm*fm + fk)
558
559	318939	if (fm>=0) then
560	0	F(i,k) = MIN(maxF(i,k), I2p2dsp1_ds(i,k) * (fm+fr))
561		else
562	318939	F(i,k) = MIN(maxF(i,k), I2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))
563		endif
564
565	318939	if ((F(i,k) >= maxF(i,k)) .or. (fr == 0.0)) then
566	59487	dFdfm(i,k) = 0.0
567		else
568	259452	dFdfm(i,k) = I2p2dsp1_ds(i,k) * ((fr + fm) / fr)
569		endif
570
571	318939	if (k > K2) then
572		! This is part of a tridiagonal solver for the actual flux.
573	318939	c1(i,k) = dFdfm(i,k-1)(dsp1_ds(i,k)b1(i))
574	318939	b1(i) = 1.0 / (1.0 - c1(i,k)*dFdfm(i,k))
575	318939	F(i,k) = MIN(b1(i)*(F(i,k)-Fprev(i,k)) + Fprev(i,k), maxF(i,k))
576	318939	if (F(i,k) >= maxF(i,k)) dFdfm(i,k) = 0.0
577		endif
578		endif
579		endif ; enddo ; enddo
580
581	436834	do k=nz-2,kb_min_act,-1 ; do i=is1,ie1
582	435138	if (do_i(i) .and. (k > kb(i))) &
583	274205	F(i,k) = MIN((F(i,k)+c1(i,k+1)*(F(i,k+1)-Fprev(i,k+1))),maxF(i,k))
584		enddo ; enddo
585
586	180	if (CS%bulkmixedlayer) then
587	56184	do i=is1,ie1
588	56184	if (do_i(i) .and. (kb(i) < nz)) then
589		! F will be increased to minF later.
590	46250	ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(F(i,kb(i)+1), minF(i,kb(i)+1))
591		else
592	9754	ea_kbp1(i) = 0.0
593		endif
594		enddo
595		call determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, min_eakb, &
596		max_eakb, kmb, is1, ie1, do_i, G, GV, CS, eakb, F_kb=F_kb, &
597		err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
598	180	dFdfm_kb=dFdfm_kb)
599	56184	do i=is1,ie1
600	56184	if (do_i(i) .and. (kb(i) < nz)) F(i,kb(i)) = F_kb(i)
601		enddo
602		endif
603
604		! Determine whether to do another iteration.
605	180	if (it < CS%max_ent_it-1) then
606
607	180	reiterate = .false.
608	56184	if (CS%bulkmixedlayer) then ; do i=is1,ie1 ; if (do_i(i)) then
609	46250	eb_kmb(i) = max(2.0*Ent_bl(i,Kmb+1) - eakb(i), 0.0)
610		endif ; enddo ; endif
611	56184	do i=is1,ie1
612	56184	did_i(i) = do_i(i) ; do_i(i) = .false.
613		enddo
614	493018	do k=kb_min_act,nz-1 ; do i=is1,ie1
615	492838	if (did_i(i) .and. (k >= kb(i))) then
616	365189	if (F(i,k) < minF(i,k)) then
617	0	F(i,k) = minF(i,k)
618	0	do_i(i) = .true. ; reiterate = .true.
619	365189	elseif (k > kb(i)) then
620	318939	if ((abs(F(i,k) - Fprev(i,k)) > tolerance) .or. &
621		((h(i,j,k) + ((1.0+dsp1_ds(i,k))*F(i,k) - &
622		(F(i,k-1) + dsp1_ds(i,k+1)F(i,k+1)))) < 0.9Angstrom)) then
623	0	do_i(i) = .true. ; reiterate = .true.
624		endif
625		else ! (k == kb(i))
626		! A more complicated test is required for the layer beneath the buffer layer,
627		! since its flux may be partially used to entrain directly from the mixed layer.
628		! Negative fluxes should not occur with the bulk mixed layer.
629	46250	if (h(i,j,k) + ((F(i,k) + eakb(i)) - &
630		(eb_kmb(i) + dsp1_ds(i,k+1)F(i,k+1))) < 0.9Angstrom) then
631	0	do_i(i) = .true. ; reiterate = .true.
632		endif
633		endif
634		endif
635		enddo ; enddo
636	180	if (.not.reiterate) exit
637		endif ! end of if (it < CS%max_ent_it-1)
638		enddo ! end of it loop
639		! This is the end of the section that might be iterated.
640
641
642	180	if (it == (CS%max_ent_it)) then
643		! Limit the flux so that the layer below is not depleted.
644		! This should only be applied to the last iteration.
645	0	do i=is1,ie1 ; if (do_i(i)) then
646	0	F(i,nz-1) = MAX(F(i,nz-1), MIN(minF(i,nz-1), 0.0))
647	0	if (kb(i) >= nz-1) then ; ea_kbp1(i) = 0.0 ; endif
648		endif ; enddo
649	0	do k=nz-2,kb_min_act,-1 ; do i=is1,ie1 ; if (do_i(i)) then
650	0	if (k>kb(i)) then
651		F(i,k) = MIN(MAX(minF(i,k),F(i,k)), (dsp1_ds(i,k+1)*F(i,k+1) + &
652		MAX(((F(i,k+1)-dsp1_ds(i,k+2)*F(i,k+2)) + &
653	0	(h(i,j,k+1) - Angstrom)), 0.5*(h(i,j,k+1)-Angstrom))))
654	0	elseif (k==kb(i)) then
655	0	ea_kbp1(i) = dsp1_ds(i,k+1)*F(i,k+1)
656		h_avail = dsp1_ds(i,k+1)F(i,k+1) + MAX(0.5(h(i,j,k+1)-Angstrom), &
657	0	((F(i,k+1)-dsp1_ds(i,k+2)*F(i,k+2)) + (h(i,j,k+1) - Angstrom)))
658	0	if ((F(i,k) > 0.0) .and. (F(i,k) > h_avail)) then
659	0	F_kb(i) = MAX(0.0, h_avail)
660	0	F(i,k) = F_kb(i)
661	0	if ((F_kb(i) < maxF_kb(i)) .and. (eakb_maxF(i) <= eakb(i))) &
662	0	eakb(i) = eakb_maxF(i)
663		call F_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &
664	0	G, GV, CS, eakb)
665		endif
666		endif
667		endif ; enddo ; enddo
668
669
670	0	if (CS%bulkmixedlayer) then ; do i=is1,ie1
671	0	if (do_i(i) .and. (kb(i) < nz)) then
672		h_avail = eakb(i) + MAX(0.5*(h_bl(i,kmb+1)-Angstrom), &
673	0	(F_kb(i)-ea_kbp1(i)) + (h_bl(i,kmb+1)-Angstrom))
674		! Ensure that 0 < eb_kmb < h_avail.
675	0	Ent_bl(i,Kmb+1) = MIN(Ent_bl(i,Kmb+1),0.5*(eakb(i) + h_avail))
676
677	0	eb_kmb(i) = max(2.0*Ent_bl(i,Kmb+1) - eakb(i), 0.0)
678		endif
679		enddo ; endif
680
681		! Limit the flux so that the layer above is not depleted.
682	0	do k=kb_min_act+1,nz-1 ; do i=is1,ie1 ; if (do_i(i)) then
683	0	if ((.not.CS%bulkmixedlayer) .or. (k > kb(i)+1)) then
684		F(i,k) = MIN(F(i,k), ds_dsp1(i,k)*( ((F(i,k-1) + &
685	0	dsp1_ds(i,k-1)*F(i,k-1)) - F(i,k-2)) + (h(i,j,k-1) - Angstrom)))
686	0	F(i,k) = MAX(F(i,k),MIN(minF(i,k),0.0))
687	0	elseif (k == kb(i)+1) then
688		F(i,k) = MIN(F(i,k), ds_dsp1(i,k)*( ((F(i,k-1) + eakb(i)) - &
689	0	eb_kmb(i)) + (h(i,j,k-1) - Angstrom)))
690	0	F(i,k) = MAX(F(i,k),MIN(minF(i,k),0.0))
691		endif
692		endif ; enddo ; enddo
693		endif ! (it == (CS%max_ent_it))
694
695	180	call F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb)
696
697		! Calculate the layer thicknesses after the entrainment to constrain the
698		! corrective fluxes.
699	180	if (associated(tv%eqn_of_state)) then
700	64980	do i=is,ie
701	64800	h_guess(i,1) = (h(i,j,1) - Angstrom) + (eb(i,j,1) - ea(i,j,2))
702	64800	h_guess(i,nz) = (h(i,j,nz) - Angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))
703	64800	if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0
704	64980	if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0
705		enddo
706	1299780	do k=2,nz-1 ; do i=is,ie
707		h_guess(i,k) = (h(i,j,k) - Angstrom) + ((ea(i,j,k) - eb(i,j,k-1)) + &
708	1296000	(eb(i,j,k) - ea(i,j,k+1)))
709	1299600	if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0
710		enddo ; enddo
711	180	if (CS%bulkmixedlayer) then
712		call determine_dSkb(h_bl, Sref, Ent_bl, eakb, is, ie, kmb, G, GV, &
713	180	.true., dS_kb, dS_anom_lim=dS_anom_lim)
714	1876	do k=nz-1,kb_min,-1
715	1696	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, EOSdom)
716	612436	do i=is,ie
717	612256	if ((k>kb(i)) .and. (F(i,k) > 0.0)) then
718		! Within a time step, a layer may entrain no more than its
719		! thickness for correction. This limitation should apply
720		! extremely rarely, but precludes undesirable behavior.
721		! Note: Corrected a sign/logic error & factor of 2 error, and
722		! the layers tracked the target density better, mostly due to
723		! the factor of 2 error.
724		F_cor = h(i,j,k) * MIN(1.0 , MAX(-ds_dsp1(i,k), &
725	314014	(GV%Rlay(k) - Rcv(i)) / (GV%Rlay(k+1)-GV%Rlay(k))) )
726
727		! Ensure that (1) Entrainments are positive, (2) Corrections in
728		! a layer cannot deplete the layer itself (very generously), and
729		! (3) a layer can take no more than a quarter the mass of its
730		! neighbor.
731	314014	if (F_cor > 0.0) then
732		F_cor = MIN(F_cor, 0.9F(i,k), ds_dsp1(i,k)0.5*h_guess(i,k), &
733	32876	0.25*h_guess(i,k+1))
734		else
735		F_cor = -MIN(-F_cor, 0.9F(i,k), 0.5h_guess(i,k), &
736	281138	0.25ds_dsp1(i,k)h_guess(i,k-1) )
737		endif
738
739	314014	ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*F_cor
740	314014	eb(i,j,k) = eb(i,j,k) + F_cor
741	296546	elseif ((k==kb(i)) .and. (F(i,k) > 0.0)) then
742		! Rho_cor is the density anomaly that needs to be corrected,
743		! taking into account that the true potential density of the
744		! deepest buffer layer is not exactly what is returned as dS_kb.
745	46208	dS_kb_eff = 2.0*dS_kb(i) - dS_anom_lim(i) ! Could be negative!!!
746	46208	Rho_cor = h(i,j,k) * (GV%Rlay(k)-Rcv(i)) + eakb(i)*dS_anom_lim(i)
747
748		! Ensure that -.9eakb < ea_cor < .9eakb
749	46208	if (abs(Rho_cor) < abs(0.9eakb(i)dS_kb_eff)) then
750	43198	ea_cor = -Rho_cor / dS_kb_eff
751		else
752	3010	ea_cor = sign(0.9eakb(i),-Rho_cordS_kb_eff)
753		endif
754
755	46208	if (ea_cor > 0.0) then
756		! Ensure that -F_cor < 0.5*h_guess
757		ea_cor = MIN(ea_cor, 0.5*(max_eakb(i) - eakb(i)), &
758	2707	0.5h_guess(i,k) / (dS_kb(i) I_dSkbp1(i)))
759		else
760		! Ensure that -ea_cor < 0.5h_guess & F_cor < 0.25h_guess(k+1)
761		ea_cor = -MIN(-ea_cor, 0.5*h_guess(i,k), &
762	43501	0.25h_guess(i,k+1) / (dS_kb(i) I_dSkbp1(i)))
763		endif
764
765	46208	ea(i,j,k) = ea(i,j,k) + ea_cor
766	46208	eb(i,j,k) = eb(i,j,k) - (dS_kb(i) * I_dSkbp1(i)) * ea_cor
767	250338	elseif (k < kb(i)) then
768		! Repetitive, unless ea(kb) has been corrected.
769	245371	ea(i,j,k) = ea(i,j,k+1)
770		endif
771		enddo
772		enddo
773	492584	do k=kb_min-1,K2,-1 ; do i=is,ie
774	492404	ea(i,j,k) = ea(i,j,k+1)
775		enddo ; enddo
776
777		! Repetitive, unless ea(kb) has been corrected.
778	180	k=kmb
779	64980	do i=is,ie
780		! Do not adjust eb through the base of the buffer layers, but it
781		! may be necessary to change entrainment from above.
782	64800	h1 = (h(i,j,k) - Angstrom) + (eb(i,j,k) - ea(i,j,k+1))
783	64980	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
784		enddo
785	130140	do k=kmb-1,2,-1 ; do i=is,ie
786		! Determine the entrainment from below for each buffer layer.
787	129600	eb(i,j,k) = max(2.0*Ent_bl(i,K+1) - ea(i,j,k+1), 0.0)
788
789		! Determine the entrainment from above for each buffer layer.
790	129600	h1 = (h(i,j,k) - Angstrom) + (eb(i,j,k) - ea(i,j,k+1))
791	129960	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
792		enddo ; enddo
793	64980	do i=is,ie
794	64980	eb(i,j,1) = max(2.0*Ent_bl(i,2) - ea(i,j,2), 0.0)
795		enddo
796
797		else ! not bulkmixedlayer
798	0	do k=K2,nz-1
799	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, EOSdom)
800	0	do i=is,ie ; if (F(i,k) > 0.0) then
801		! Within a time step, a layer may entrain no more than
802		! its thickness for correction. This limitation should
803		! apply extremely rarely, but precludes undesirable
804		! behavior.
805		F_cor = h(i,j,k) * MIN(dsp1_ds(i,k) , MAX(-1.0, &
806	0	(GV%Rlay(k) - Rcv(i)) / (GV%Rlay(k+1)-GV%Rlay(k))) )
807
808		! Ensure that (1) Entrainments are positive, (2) Corrections in
809		! a layer cannot deplete the layer itself (very generously), and
810		! (3) a layer can take no more than a quarter the mass of its
811		! neighbor.
812	0	if (F_cor >= 0.0) then
813		F_cor = MIN(F_cor, 0.9F(i,k), 0.5dsp1_ds(i,k)*h_guess(i,k), &
814	0	0.25*h_guess(i,k+1))
815		else
816		F_cor = -MIN(-F_cor, 0.9F(i,k), 0.5h_guess(i,k), &
817	0	0.25ds_dsp1(i,k)h_guess(i,k-1) )
818		endif
819
820	0	ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*F_cor
821	0	eb(i,j,k) = eb(i,j,k) + F_cor
822		endif ; enddo
823		enddo
824		endif
825
826		endif ! associated(tv%eqn_of_state))
827
828	180	if (CS%id_Kd > 0) then
829	0	Idt = (GV%H_to_m*US%m_to_Z) / dt
830	0	do k=2,nz-1 ; do i=is,ie
831	0	if (k<kb(i)) then ; Kd_here = 0.0 ; else
832		Kd_here = F(i,k) * ( h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
833	0	(eb(i,j,k) - ea(i,j,k+1))) ) / (I2p2dsp1_ds(i,k) * grats(i,k))
834		endif
835
836	0	Kd_eff(i,j,k) = MAX(dtKd(i,k), Kd_here)*Idt
837		enddo ; enddo
838	0	do i=is,ie
839	0	Kd_eff(i,j,1) = dtKd(i,1)*Idt
840	0	Kd_eff(i,j,nz) = dtKd(i,nz)*Idt
841		enddo
842		endif
843
844	180	if (CS%id_diff_work > 0) then
845	0	if (GV%Boussinesq .or. .not.associated(tv%eqn_of_state)) then
846	0	g_2dt = 0.5 * GV%H_to_Z*2 (GV%g_Earth_Z_T2 / dt)
847		else
848	0	g_2dt = 0.5 * GV%H_to_RZ*2 (GV%g_Earth_Z_T2 / dt)
849		endif
850	0	do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo
851	0	if (associated(tv%eqn_of_state)) then
852	0	if (associated(fluxes%p_surf)) then
853	0	do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ; enddo
854		else
855	0	do i=is,ie ; pressure(i) = 0.0 ; enddo
856		endif
857	0	do K=2,nz
858	0	do i=is,ie ; pressure(i) = pressure(i) + (GV%g_EarthGV%H_to_RZ)h(i,j,k-1) ; enddo
859	0	do i=is,ie
860	0	if (k==kb(i)) then
861	0	T_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))
862	0	S_eos(i) = 0.5*(tv%S(i,j,kmb) + tv%S(i,j,k))
863		else
864	0	T_eos(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
865	0	S_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
866		endif
867		enddo
868	0	if (GV%Boussinesq) then
869		call calculate_density_derivs(T_EOS, S_EOS, pressure, dRho_dT, dRho_dS, &
870	0	tv%eqn_of_state, EOSdom)
871	0	do i=is,ie
872	0	if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,K) = 0.0
873		else
874	0	if (k==kb(i)) then
875		dRho = dRho_dT(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &
876	0	dRho_dS(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))
877		else
878		dRho = dRho_dT(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
879	0	dRho_dS(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))
880		endif
881		diff_work(i,j,K) = g_2dt * dRho * &
882		(ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
883	0	eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
884		endif
885		enddo
886		else
887		call calculate_specific_vol_derivs(T_EOS, S_EOS, pressure, dSpV_dT, dSpV_dS, &
888	0	tv%eqn_of_state, EOSdom)
889
890	0	do i=is,ie
891	0	if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,K) = 0.0
892		else
893	0	if (k==kb(i)) then
894		dSpV = dSpV_dT(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &
895	0	dSpV_dS(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))
896		else
897		dSpV = dSpV_dT(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
898	0	dSpV_dS(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))
899		endif
900		diff_work(i,j,K) = -g_2dt * dSpV * &
901		(ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
902	0	eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
903		endif
904		enddo
905		endif
906		enddo
907		else
908	0	do K=2,nz ; do i=is,ie
909		diff_work(i,j,K) = g_2dt * (GV%Rlay(k)-GV%Rlay(k-1)) * &
910		(ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
911	0	eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
912		enddo ; enddo
913		endif
914		endif
915
916	64981	do i=is,ie ; kb_out(i,j) = kb(i) ; enddo
917
918		enddo ! end of j loop
919
920		! Offer diagnostic fields for averaging.
921	1	if (CS%id_Kd > 0) call post_data(CS%id_Kd, Kd_eff, CS%diag)
922	1	if (CS%id_Kd > 0) deallocate(Kd_eff)
923	1	if (CS%id_diff_work > 0) call post_data(CS%id_diff_work, diff_work, CS%diag)
924	1	if (CS%id_diff_work > 0) deallocate(diff_work)
925
926	1	end subroutine entrainment_diffusive
927
928		!> This subroutine calculates the actual entrainments (ea and eb) and the
929		!! amount of surface forcing that is applied to each layer if there is no bulk
930		!! mixed layer.
931	180	subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb)
932		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
933		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
934		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: F !< The density flux through a layer within
935		!! a time step divided by the density
936		!! difference across the interface below
937		!! the layer [H ~> m or kg m-2].
938		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
939		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
940		integer, dimension(SZI_(G)), intent(in) :: kb !< The index of the lightest layer denser than
941		!! the deepest buffer layer.
942		integer, intent(in) :: kmb !< The number of mixed and buffer layers.
943		integer, intent(in) :: j !< The meridional index upon which to work.
944		type(entrain_diffusive_CS), intent(in) :: CS !< This module's control structure.
945		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dsp1_ds !< The ratio of coordinate variable
946		!! differences across the interfaces below
947		!! a layer over the difference across the
948		!! interface above the layer [nondim].
949		real, dimension(SZI_(G)), intent(in) :: eakb !< The entrainment from above by the layer
950		!! below the buffer layer [H ~> m or kg m-2].
951		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Ent_bl !< The average entrainment upward and
952		!! downward across each interface around
953		!! the buffer layers [H ~> m or kg m-2].
954		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
955		intent(inout) :: ea !< The amount of fluid entrained from the layer
956		!! above within this time step [H ~> m or kg m-2].
957		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
958		intent(inout) :: eb !< The amount of fluid entrained from the layer
959		!! below within this time step [H ~> m or kg m-2].
960
961		real :: h1 ! The thickness in excess of the minimum that will remain
962		! after exchange with the layer below [H ~> m or kg m-2].
963		integer :: i, k, is, ie, nz
964
965	180	is = G%isc ; ie = G%iec ; nz = GV%ke
966
967	64980	do i=is,ie
968	64980	ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0
969		enddo
970	180	if (CS%bulkmixedlayer) then
971	64980	do i=is,ie
972	64980	eb(i,j,kmb) = max(2.0*Ent_bl(i,Kmb+1) - eakb(i), 0.0)
973		enddo
974	1104840	do k=nz-1,kmb+1,-1 ; do i=is,ie
975	1104660	if (k > kb(i)) then
976		! With a bulk mixed layer, surface buoyancy fluxes are applied
977		! elsewhere, so F should always be nonnegative.
978	318939	ea(i,j,k) = dsp1_ds(i,k)*F(i,k)
979	318939	eb(i,j,k) = F(i,k)
980	782661	elseif (k == kb(i)) then
981	46250	ea(i,j,k) = eakb(i)
982	46250	eb(i,j,k) = F(i,k)
983	736411	elseif (k == kb(i)-1) then
984	46250	ea(i,j,k) = ea(i,j,k+1)
985	46250	eb(i,j,k) = eb(i,j,kmb)
986		else
987	690161	ea(i,j,k) = ea(i,j,k+1)
988		! Add the entrainment of the thin interior layers to eb going
989		! up into the buffer layer.
990	690161	eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - GV%Angstrom_H)
991		endif
992		enddo ; enddo
993	180	k = kmb
994	64980	do i=is,ie
995		! Adjust the previously calculated entrainment from below by the deepest
996		! buffer layer to account for entrainment of thin interior layers .
997	64800	if (kb(i) > kmb+1) &
998	64800	eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - GV%Angstrom_H)
999
1000		! Determine the entrainment from above for each buffer layer.
1001	64800	h1 = (h(i,j,k) - GV%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
1002	64980	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
1003		enddo
1004	130140	do k=kmb-1,2,-1 ; do i=is,ie
1005		! Determine the entrainment from below for each buffer layer.
1006	129600	eb(i,j,k) = max(2.0*Ent_bl(i,K+1) - ea(i,j,k+1), 0.0)
1007
1008		! Determine the entrainment from above for each buffer layer.
1009	129600	h1 = (h(i,j,k) - GV%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
1010	129960	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
1011		! if (h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)
1012		! elseif (Ent_bl(i,K)+0.5h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)-0.5h1
1013		! else ; ea(i,j,k) = -h1 ; endif
1014		enddo ; enddo
1015	64980	do i=is,ie
1016	64800	eb(i,j,1) = max(2.0*Ent_bl(i,2) - ea(i,j,2), 0.0)
1017	64980	ea(i,j,1) = 0.0
1018		enddo
1019		else ! not BULKMIXEDLAYER
1020		! Calculate the entrainment by each layer from above and below.
1021		! Entrainment is always positive, but F may be negative due to
1022		! surface buoyancy fluxes.
1023	0	do i=is,ie
1024	0	ea(i,j,1) = 0.0 ; eb(i,j,1) = MAX(F(i,1),0.0)
1025	0	ea(i,j,2) = dsp1_ds(i,2)*F(i,2) - MIN(F(i,1),0.0)
1026		enddo
1027
1028	0	do k=2,nz-1 ; do i=is,ie
1029	0	eb(i,j,k) = MAX(F(i,k),0.0)
1030	0	ea(i,j,k+1) = dsp1_ds(i,k+1)*F(i,k+1) - (F(i,k)-eb(i,j,k))
1031	0	if (ea(i,j,k+1) < 0.0) then
1032	0	eb(i,j,k) = eb(i,j,k) - ea(i,j,k+1)
1033	0	ea(i,j,k+1) = 0.0
1034		endif
1035		enddo ; enddo
1036		endif ! end BULKMIXEDLAYER
1037	180	end subroutine F_to_ent
1038
1039		!> This subroutine sets the average entrainment across each of the interfaces
1040		!! between buffer layers within a timestep. It also causes thin and relatively
1041		!! light interior layers to be entrained by the deepest buffer layer.
1042		!! Also find the initial coordinate potential densities (Sref) of each layer.
1043	180	subroutine set_Ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)
1044		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1045		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1046		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1047		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1048		real, dimension(SZI_(G),SZK_(GV)+1), &
1049		intent(in) :: dtKd_int !< The diapycnal diffusivity across
1050		!! each interface times the time step
1051		!! [H2 ~> m2 or kg2 m-4].
1052		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
1053		!! available thermodynamic fields. Absent
1054		!! fields have NULL ptrs.
1055		integer, dimension(SZI_(G)), intent(inout) :: kb !< The index of the lightest layer denser
1056		!! than the buffer layer or 1 if there is
1057		!! no buffer layer.
1058		integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1059		logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which
1060		!! i-points to work on.
1061		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1062		type(entrain_diffusive_CS), intent(in) :: CS !< This module's control structure.
1063		integer, intent(in) :: j !< The meridional index upon which to work.
1064		real, dimension(SZI_(G),SZK_(GV)+1), &
1065		intent(out) :: Ent_bl !< The average entrainment upward and
1066		!! downward across each interface around
1067		!! the buffer layers [H ~> m or kg m-2].
1068		real, dimension(SZI_(G),SZK_(GV)), intent(out) :: Sref !< The coordinate potential density minus
1069		!! 1000 for each layer [R ~> kg m-3].
1070		real, dimension(SZI_(G),SZK_(GV)), intent(out) :: h_bl !< The thickness of each layer [H ~> m or kg m-2].
1071
1072		! This subroutine sets the average entrainment across each of the interfaces
1073		! between buffer layers within a timestep. It also causes thin and relatively
1074		! light interior layers to be entrained by the deepest buffer layer.
1075		! Also find the initial coordinate potential densities (Sref) of each layer.
1076		! Does there need to be limiting when the layers below are all thin?
1077
1078		! Local variables
1079		real, dimension(SZI_(G)) :: &
1080	360	b1, d1, & ! Variables used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] and [nondim].
1081	360	Rcv, & ! Value of the coordinate variable (potential density)
1082		! based on the simulated T and S and P_Ref [R ~> kg m-3].
1083	360	pres, & ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].
1084	360	frac_rem, & ! The fraction of the diffusion remaining [nondim].
1085	360	h_interior ! The interior thickness available for entrainment [H ~> m or kg m-2].
1086		real, dimension(SZI_(G), SZK_(GV)) :: &
1087	360	S_est ! An estimate of the coordinate potential density - 1000 after
1088		! entrainment for each layer [R ~> kg m-3].
1089		real :: dh ! An available thickness [H ~> m or kg m-2].
1090		real :: Kd_x_dt ! The diffusion that remains after thin layers are
1091		! entrained [H2 ~> m2 or kg2 m-4].
1092		real :: h_neglect ! A thickness that is so small it is usually lost
1093		! in roundoff and can be neglected [H ~> m or kg m-2].
1094		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1095		integer :: i, k, is, ie, nz
1096	180	is = G%isc ; ie = G%iec ; nz = GV%ke
1097
1098	180	h_neglect = GV%H_subroundoff
1099
1100	64980	do i=is,ie ; pres(i) = tv%P_Ref ; enddo
1101	180	EOSdom(:) = EOS_domain(G%HI)
1102	900	do k=1,kmb
1103	720	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, EOSdom)
1104	260100	do i=is,ie
1105	259200	h_bl(i,k) = h(i,j,k) + h_neglect
1106	259920	Sref(i,k) = Rcv(i) - CS%Rho_sig_off
1107		enddo
1108		enddo
1109
1110	64980	do i=is,ie
1111	64980	h_interior(i) = 0.0 ; Ent_bl(i,1) = 0.0
1112		! if (kb(i) > nz) Ent_bl(i,Kmb+1) = 0.0
1113		enddo
1114
1115	195120	do k=2,kmb ; do i=is,ie
1116	194940	if (do_i(i)) then
1117	138750	Ent_bl(i,K) = min(2.0 * dtKd_int(i,K) / (h(i,j,k-1) + h(i,j,k) + h_neglect), CS%max_Ent)
1118	55650	else ; Ent_bl(i,K) = 0.0 ; endif
1119		enddo ; enddo
1120
1121		! Determine the coordinate density of the bottommost buffer layer if there
1122		! is no entrainment from the layers below. This is a partial solver, based
1123		! on the first pass of a tridiagonal solver, as the values in the upper buffer
1124		! layers are not needed.
1125
1126	64980	do i=is,ie
1127	64800	b1(i) = 1.0 / (h_bl(i,1) + Ent_bl(i,2))
1128	64800	d1(i) = h_bl(i,1) * b1(i) ! = 1.0 - Ent_bl(i,2)*b1(i)
1129	64980	S_est(i,1) = (h_bl(i,1)Sref(i,1)) b1(i)
1130		enddo
1131	130140	do k=2,kmb-1 ; do i=is,ie
1132	129600	b1(i) = 1.0 / ((h_bl(i,k) + Ent_bl(i,K+1)) + d1(i)*Ent_bl(i,K))
1133	129600	d1(i) = (h_bl(i,k) + d1(i)Ent_bl(i,K)) b1(i) ! = 1.0 - Ent_bl(i,K+1)*b1(i)
1134	129960	S_est(i,k) = (h_bl(i,k)Sref(i,k) + Ent_bl(i,K)S_est(i,k-1)) * b1(i)
1135		enddo ; enddo
1136	64980	do i=is,ie
1137		S_est(i,kmb) = (h_bl(i,kmb)Sref(i,kmb) + Ent_bl(i,Kmb)S_est(i,kmb-1)) / &
1138	64800	(h_bl(i,kmb) + d1(i)*Ent_bl(i,Kmb))
1139	64980	frac_rem(i) = 1.0
1140		enddo
1141
1142		! Entrain any thin interior layers that are lighter (in the coordinate
1143		! potential density) than the deepest buffer layer will be, and adjust kb.
1144	64980	do i=is,ie ; kb(i) = nz+1 ; if (do_i(i)) kb(i) = kmb+1 ; enddo
1145
1146	1169820	do k=kmb+1,nz ; do i=is,ie ; if (do_i(i)) then
1147	832500	if ((k == kb(i)) .and. (S_est(i,kmb) > (GV%Rlay(k) - CS%Rho_sig_off))) then
1148	423742	if (4.0dtKd_int(i,Kmb+1)frac_rem(i) > &
1149		(h_bl(i,kmb) + h(i,j,k)) * (h(i,j,k) - GV%Angstrom_H)) then
1150		! Entrain this layer into the buffer layer and move kb down.
1151	421061	dh = max((h(i,j,k) - GV%Angstrom_H), 0.0)
1152	421061	if (dh > 0.0) then
1153		frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &
1154	232188	(4.0*dtKd_int(i,Kmb+1))
1155		Sref(i,kmb) = (h_bl(i,kmb)Sref(i,kmb) + dh(GV%Rlay(k)-CS%Rho_sig_off)) / &
1156	232188	(h_bl(i,kmb) + dh)
1157	232188	h_bl(i,kmb) = h_bl(i,kmb) + dh
1158		S_est(i,kmb) = (h_bl(i,kmb)Sref(i,kmb) + Ent_bl(i,Kmb)S_est(i,kmb-1)) / &
1159	232188	(h_bl(i,kmb) + d1(i)*Ent_bl(i,Kmb))
1160		endif
1161	421061	kb(i) = kb(i) + 1
1162		endif
1163		endif
1164		endif ; enddo ; enddo
1165
1166		! This is where variables are be set up with a different vertical grid
1167		! in which the (newly?) massless layers are taken out.
1168	1169820	do k=nz,kmb+1,-1 ; do i=is,ie
1169	1166400	if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-GV%Angstrom_H)
1170	1169640	if (k==kb(i)) then
1171	46250	h_bl(i,kmb+1) = h(i,j,k) ; Sref(i,kmb+1) = GV%Rlay(k) - CS%Rho_sig_off
1172	1120150	elseif (k==kb(i)+1) then
1173	46250	h_bl(i,kmb+2) = h(i,j,k) ; Sref(i,kmb+2) = GV%Rlay(k) - CS%Rho_sig_off
1174		endif
1175		enddo ; enddo
1176	64980	do i=is,ie ; if (kb(i) >= nz) then
1177	18550	h_bl(i,kmb+1) = h(i,j,nz)
1178	18550	Sref(i,kmb+1) = GV%Rlay(nz) - CS%Rho_sig_off
1179	18550	h_bl(i,kmb+2) = GV%Angstrom_H
1180	18550	Sref(i,kmb+2) = Sref(i,kmb+1) + (GV%Rlay(nz) - GV%Rlay(nz-1))
1181		endif ; enddo
1182
1183		! Perhaps we should revisit the way that the average entrainment between the
1184		! buffer layer and the interior is calculated so that it is not unduly
1185		! limited when the layers are less than sqrt(Kd * dt) thick?
1186	64980	do i=is,ie ; if (do_i(i)) then
1187	46250	Kd_x_dt = frac_rem(i) * dtKd_int(i,Kmb+1)
1188	46250	if ((Kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0)) then
1189	21	Ent_bl(i,Kmb+1) = 0.0
1190		else
1191		! If the combined layers are exceptionally thin, use sqrt(Kd*dt) as the
1192		! estimate of the thickness in the denominator of the thickness diffusion.
1193		Ent_bl(i,Kmb+1) = MIN(0.5*h_interior(i), sqrt(Kd_x_dt), &
1194	46229	Kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))
1195		endif
1196		else
1197	18550	Ent_bl(i,Kmb+1) = 0.0
1198		endif ; enddo
1199
1200	180	end subroutine set_Ent_bl
1201
1202		!> This subroutine determines the reference density difference between the
1203		!! bottommost buffer layer and the first interior after the mixing between mixed
1204		!! and buffer layers and mixing with the layer below. Within the mixed and buffer
1205		!! layers, entrainment from the layer above is increased when it is necessary to
1206		!! keep the layers from developing a negative thickness; otherwise it equals
1207		!! Ent_bl. At each interface, the upward and downward fluxes average out to
1208		!! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.
1209		!! The density difference across the first interior layer may also be returned.
1210		!! It could also be limited to avoid negative values or values that greatly
1211		!! exceed the density differences across an interface.
1212		!! Additionally, the partial derivatives of dSkb and dSlay with E_kb could
1213		!! also be returned.
1214	24600	subroutine determine_dSkb(h_bl, Sref, Ent_bl, E_kb, is, ie, kmb, G, GV, limit, &
1215	18719	dSkb, ddSkb_dE, dSlay, ddSlay_dE, dS_anom_lim, do_i_in)
1216		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1217		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid
1218		!! structure.
1219		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h_bl !< Layer thickness [H ~> m or kg m-2]
1220		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Sref !< Reference potential density [R ~> kg m-3]
1221		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Ent_bl !< The average entrainment upward and
1222		!! downward across each interface
1223		!! around the buffer layers [H ~> m or kg m-2].
1224		real, dimension(SZI_(G)), intent(in) :: E_kb !< The entrainment by the top interior
1225		!! layer [H ~> m or kg m-2].
1226		integer, intent(in) :: is !< The start of the i-index range to work on.
1227		integer, intent(in) :: ie !< The end of the i-index range to work on.
1228		integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1229		logical, intent(in) :: limit !< If true, limit dSkb and dSlay to
1230		!! avoid negative values.
1231		real, dimension(SZI_(G)), intent(inout) :: dSkb !< The limited potential density
1232		!! difference across the interface
1233		!! between the bottommost buffer layer
1234		!! and the topmost interior layer. [R ~> kg m-3]
1235		!! dSkb > 0.
1236		real, dimension(SZI_(G)), optional, intent(inout) :: ddSkb_dE !< The partial derivative of dSkb
1237		!! with E [R H-1 ~> kg m-4 or m-1].
1238		real, dimension(SZI_(G)), optional, intent(inout) :: dSlay !< The limited potential density
1239		!! difference across the topmost
1240		!! interior layer. 0 < dSkb [R ~> kg m-3]
1241		real, dimension(SZI_(G)), optional, intent(inout) :: ddSlay_dE !< The partial derivative of dSlay
1242		!! with E [R H-1 ~> kg m-4 or m-1].
1243		real, dimension(SZI_(G)), optional, intent(inout) :: dS_anom_lim !< A limiting value to use for
1244		!! the density anomalies below the
1245		!! buffer layer [R ~> kg m-3].
1246		logical, dimension(SZI_(G)), optional, intent(in) :: do_i_in !< If present, determines which
1247		!! columns are worked on.
1248
1249		! Note that dSkb, ddSkb_dE, dSlay, ddSlay_dE, and dS_anom_lim are declared
1250		! intent inout because they should not change where do_i_in is false.
1251
1252		! This subroutine determines the reference density difference between the
1253		! bottommost buffer layer and the first interior after the mixing between mixed
1254		! and buffer layers and mixing with the layer below. Within the mixed and buffer
1255		! layers, entrainment from the layer above is increased when it is necessary to
1256		! keep the layers from developing a negative thickness; otherwise it equals
1257		! Ent_bl. At each interface, the upward and downward fluxes average out to
1258		! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.
1259		! The density difference across the first interior layer may also be returned.
1260		! It could also be limited to avoid negative values or values that greatly
1261		! exceed the density differences across an interface.
1262		! Additionally, the partial derivatives of dSkb and dSlay with E_kb could
1263		! also be returned.
1264
1265		! Local variables
1266		real, dimension(SZI_(G),SZK_(GV)) :: &
1267	16400	b1, c1, & ! b1 [H-1 ~> m-1 or m2 kg-1] and c1 [nondim] are variables used by the tridiagonal solver.
1268	16400	S, dS_dE, & ! The coordinate density [R ~> kg m-3] and its derivative with E [R H-1 ~> kg m-4 or m-1].
1269	16400	ea, dea_dE, & ! The entrainment from above [H ~> m or kg m-2] and its derivative with E [nondim].
1270	16400	eb, deb_dE ! The entrainment from below [H ~> m or kg m-2] and its derivative with E [nondim].
1271	16400	real :: deriv_dSkb(SZI_(G)) ! The limited derivative of the new density difference across the base of
1272		! the buffer layers with the new density of the bottommost buffer layer [nondim]
1273	16400	real :: d1(SZI_(G)) ! d1 = 1.0-c1 is also used by the tridiagonal solver [nondim].
1274		real :: src ! A source term for dS_dR [R ~> kg m-3].
1275		real :: h1 ! The thickness in excess of the minimum that will remain
1276		! after exchange with the layer below [H ~> m or kg m-2].
1277	16400	logical, dimension(SZI_(G)) :: do_i
1278		real :: h_neglect ! A thickness that is so small it is usually lost
1279		! in roundoff and can be neglected [H ~> m or kg m-2].
1280		real :: h_tr ! h_tr is h at tracer points with a tiny thickness
1281		! added to ensure positive definiteness [H ~> m or kg m-2].
1282		real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
1283		real :: rat ! A ratio of density differences [nondim]
1284		real :: dS_kbp1 ! The density difference between the top two interior layers [R ~> kg m-3].
1285		real :: IdS_kbp1 ! The inverse of dS_kbp1 [R-1 ~> m3 kg-1]
1286		real :: deriv_dSLay ! The derivative of the projected density difference across the topmost interior
1287		! layer with the density difference across the interface above it [nondim]
1288		real :: Inv_term ! The inverse of a nondimensional expression [nondim]
1289		real :: f1, df1_drat ! Temporary variables [nondim].
1290		real :: z, dz_drat, f2, df2_dz, expz ! Temporary variables [nondim].
1291		real :: eps_dSLay, eps_dSkb ! Small nondimensional constants [nondim].
1292		integer :: i, k
1293
1294	8200	if (present(ddSlay_dE) .and. .not.present(dSlay)) call MOM_error(FATAL, &
1295	0	"In deterimine_dSkb, ddSLay_dE may only be present if dSlay is.")
1296
1297	8200	h_neglect = GV%H_subroundoff
1298
1299	1062957	do i=is,ie
1300	1054757	ea(i,kmb+1) = E_kb(i) ; dea_dE(i,kmb+1) = 1.0
1301	1054757	S(i,kmb+1) = Sref(i,kmb+1) ; dS_dE(i,kmb+1) = 0.0
1302	1054757	b1(i,kmb+1) = 0.0
1303	1054757	d1(i) = 1.0
1304	1062957	do_i(i) = .true.
1305		enddo
1306	8200	if (present(do_i_in)) then
1307	861599	do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1308		endif
1309	4260028	do k=kmb,1,-1 ; do i=is,ie
1310	4219028	if (do_i(i)) then
1311		! The do_i test here is only for efficiency.
1312		! Determine the entrainment from below for each buffer layer.
1313	3002508	if (2.0*Ent_bl(i,K+1) > ea(i,k+1)) then
1314	2193836	eb(i,k) = 2.0*Ent_bl(i,K+1) - ea(i,k+1) ; deb_dE(i,k) = -dea_dE(i,k+1)
1315		else
1316	808672	eb(i,k) = 0.0 ; deb_dE(i,k) = 0.0
1317		endif
1318
1319		! Determine the entrainment from above for each buffer layer.
1320	3002508	h1 = (h_bl(i,k) - GV%Angstrom_H) + (eb(i,k) - ea(i,k+1))
1321	3002508	if (h1 >= 0.0) then
1322	1835991	ea(i,k) = Ent_bl(i,K) ; dea_dE(i,k) = 0.0
1323	1166517	elseif (Ent_bl(i,K) + 0.5*h1 >= 0.0) then
1324	637103	ea(i,k) = Ent_bl(i,K) - 0.5*h1
1325	637103	dea_dE(i,k) = 0.5*(dea_dE(i,k+1) - deb_dE(i,k))
1326		else
1327	529414	ea(i,k) = -h1
1328	529414	dea_dE(i,k) = dea_dE(i,k+1) - deb_dE(i,k)
1329		endif
1330		else
1331	1216520	ea(i,k) = 0.0 ; dea_dE(i,k) = 0.0 ; eb(i,k) = 0.0 ; deb_dE(i,k) = 0.0
1332		endif
1333
1334		! This is the first-pass of a tridiagonal solver for S.
1335	4219028	h_tr = h_bl(i,k) + h_neglect
1336	4219028	c1(i,k) = ea(i,k+1) * b1(i,k+1)
1337	4219028	b_denom_1 = (h_tr + d1(i)*eb(i,k))
1338	4219028	b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))
1339	4219028	d1(i) = b_denom_1 * b1(i,k)
1340
1341	4251828	S(i,k) = (h_trSref(i,k) + eb(i,k)S(i,k+1)) * b1(i,k)
1342		enddo ; enddo
1343	3197071	do k=2,kmb ; do i=is,ie
1344	3188871	S(i,k) = S(i,k) + c1(i,k-1)*S(i,k-1)
1345		enddo ; enddo
1346
1347	8200	if (present(ddSkb_dE) .or. present(ddSlay_dE)) then
1348		! These two tridiagonal solvers cannot be combined because the solutions for
1349		! S are required as a source for dS_dE.
1350	2274247	do k=kmb,2,-1 ; do i=is,ie
1351	2245047	if (do_i(i) .and. (dea_dE(i,k) - deb_dE(i,k) > 0.0)) then
1352		src = (((S(i,k+1) - Sref(i,k)) * (h_bl(i,k) + h_neglect) + &
1353		(S(i,k+1) - S(i,k-1)) * ea(i,k)) * deb_dE(i,k) - &
1354		((Sref(i,k) - S(i,k-1)) * h_bl(i,k) + &
1355		(S(i,k+1) - S(i,k-1)) * eb(i,k)) * dea_dE(i,k)) / &
1356	1080757	((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))
1357	1164290	else ; src = 0.0 ; endif
1358	2266947	dS_dE(i,k) = (src + eb(i,k)dS_dE(i,k+1)) b1(i,k)
1359		enddo ; enddo
1360	755649	do i=is,ie
1361	748349	if (do_i(i) .and. (deb_dE(i,1) < 0.0)) then
1362		src = (((S(i,2) - Sref(i,1)) * (h_bl(i,1) + h_neglect)) * deb_dE(i,1)) / &
1363	286	(h_bl(i,1) + h_neglect + eb(i,1))
1364	748063	else ; src = 0.0 ; endif
1365	755649	dS_dE(i,1) = (src + eb(i,1)dS_dE(i,2)) b1(i,1)
1366		enddo
1367	2274247	do k=2,kmb ; do i=is,ie
1368	2266947	dS_dE(i,k) = dS_dE(i,k) + c1(i,k-1)*dS_dE(i,k-1)
1369		enddo ; enddo
1370		endif
1371
1372		! Now, apply any limiting and return the requested variables.
1373
1374	8200	eps_dSkb = 1.0e-6 ! Should be a small, nondimensional, positive number.
1375	8200	if (.not.limit) then
1376	288839	do i=is,ie ; if (do_i(i)) then
1377	209152	dSkb(i) = Sref(i,kmb+1) - S(i,kmb)
1378		endif ; enddo
1379	288839	if (present(ddSkb_dE)) then ; do i=is,ie ; if (do_i(i)) then
1380	209152	ddSkb_dE(i) = -1.0*dS_dE(i,kmb)
1381		endif ; enddo ; endif
1382
1383	877	if (present(dSlay)) then ; do i=is,ie ; if (do_i(i)) then
1384	0	dSlay(i) = 0.5 * (Sref(i,kmb+2) - S(i,kmb))
1385		endif ; enddo ; endif
1386	877	if (present(ddSlay_dE)) then ; do i=is,ie ; if (do_i(i)) then
1387	0	ddSlay_dE(i) = -0.5*dS_dE(i,kmb)
1388		endif ; enddo ; endif
1389		else
1390	774118	do i=is,ie ; if (do_i(i)) then
1391		! Need to ensure that 0 < dSkb <= S_kb - Sbl
1392	541475	if (Sref(i,kmb+1) - S(i,kmb) < eps_dSkb*(Sref(i,kmb+2) - Sref(i,kmb+1))) then
1393	27056	dSkb(i) = eps_dSkb * (Sref(i,kmb+2) - Sref(i,kmb+1)) ; deriv_dSkb(i) = 0.0
1394		else
1395	514419	dSkb(i) = Sref(i,kmb+1) - S(i,kmb) ; deriv_dSkb(i) = -1.0
1396		endif
1397	541475	if (present(ddSkb_dE)) ddSkb_dE(i) = deriv_dSkb(i)*dS_dE(i,kmb)
1398		endif ; enddo
1399
1400	7323	if (present(dSLay)) then
1401	3214	dz_drat = 1000.0 ! The limit of large dz_drat the same as choosing a
1402		! Heaviside function.
1403	3214	eps_dSLay = 1.0e-10 ! Should be ~= GV%Angstrom_H / sqrt(Kd*dt)
1404	460392	do i=is,ie ; if (do_i(i)) then
1405	288466	dS_kbp1 = Sref(i,kmb+2) - Sref(i,kmb+1)
1406	288466	IdS_kbp1 = 1.0 / (Sref(i,kmb+2) - Sref(i,kmb+1))
1407	288466	rat = (Sref(i,kmb+1) - S(i,kmb)) * IdS_kbp1
1408		! Need to ensure that 0 < dSLay <= 2*dSkb
1409	288466	if (rat < 0.5) then
1410		! The coefficients here are chosen so that at rat = 0.5, the value (1.5)
1411		! and first derivative (-0.5) match with the "typical" case (next).
1412		! The functional form here is arbitrary.
1413		! f1 provides a reasonable profile that matches the value and derivative
1414		! of the "typical" case at rat = 0.5, and has a maximum of less than 2.
1415	217481	Inv_term = 1.0 / (1.0-rat)
1416	217481	f1 = 2.0 - 0.125(Inv_term*2)
1417	217481	df1_drat = - 0.25(Inv_term*3)
1418
1419		! f2 ensures that dSLay goes to 0 rapidly if rat is significantly
1420		! negative.
1421	217481	z = dz_drat * rat + 4.0 ! The 4 here gives f2(0) = 0.982.
1422	217481	if (z >= 18.0) then ; f2 = 1.0 ; df2_dz = 0.0
1423	31787	elseif (z <= -58.0) then ; f2 = eps_dSLay ; df2_dz = 0.0
1424		else
1425	27865	expz = exp(z) ; Inv_term = 1.0 / (1.0 + expz)
1426	27865	f2 = (eps_dSLay + expz) * Inv_term
1427	27865	df2_dz = (1.0 - eps_dSLay) * expz * Inv_term**2
1428		endif
1429
1430	217481	dSLay(i) = dSkb(i) * f1 * f2
1431		deriv_dSLay = deriv_dSkb(i) * (f1 * f2) - (dSkb(i)IdS_kbp1) &
1432	217481	(df1_dratf2 + f1 dz_drat * df2_dz)
1433	70985	elseif (dSkb(i) <= 3.0*dS_kbp1) then
1434		! This is the "typical" case.
1435	70985	dSLay(i) = 0.5 * (dSkb(i) + dS_kbp1)
1436	70985	deriv_dSLay = 0.5 * deriv_dSkb(i) ! = -0.5
1437		else
1438	0	dSLay(i) = 2.0*dS_kbp1
1439	0	deriv_dSLay = 0.0
1440		endif
1441	288466	if (present(ddSlay_dE)) ddSlay_dE(i) = deriv_dSLay*dS_dE(i,kmb)
1442		endif ; enddo
1443		endif ! present(dSlay)
1444		endif ! Not limited.
1445
1446	202600	if (present(dS_anom_lim)) then ; do i=is,ie ; if (do_i(i)) then
1447		dS_anom_lim(i) = max(0.0, eps_dSkb * (Sref(i,kmb+2) - Sref(i,kmb+1)) - &
1448	194400	(Sref(i,kmb+1) - S(i,kmb)) )
1449		endif ; enddo ; endif
1450
1451	26379	end subroutine determine_dSkb
1452
1453		!> Given an entrainment from below for layer kb, determine a consistent
1454		!! entrainment from above, such that dSkb * ea_kb = dSkbp1 * F_kb. The input
1455		!! value of ea_kb is both the maximum value that can be obtained and the first
1456		!! guess of the iterations. Ideally ea_kb should be an under-estimate
1457	948	subroutine F_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &
1458	948	G, GV, CS, ea_kb, tol_in)
1459		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1460		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1461		real, dimension(SZI_(G),SZK_(GV)), &
1462		intent(in) :: h_bl !< Layer thickness, with the top interior
1463		!! layer at k-index kmb+1 [H ~> m or kg m-2].
1464		real, dimension(SZI_(G),SZK_(GV)), &
1465		intent(in) :: Sref !< The coordinate reference potential density,
1466		!! with the value of the topmost interior layer
1467		!! at index kmb+1 [R ~> kg m-3].
1468		real, dimension(SZI_(G),SZK_(GV)), &
1469		intent(in) :: Ent_bl !< The average entrainment upward and downward
1470		!! across each interface around the buffer layers,
1471		!! [H ~> m or kg m-2].
1472		real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in reference
1473		!! potential density across the base of the
1474		!! uppermost interior layer [R-1 ~> m3 kg-1].
1475		real, dimension(SZI_(G)), intent(in) :: F_kb !< The entrainment from below by the
1476		!! uppermost interior layer [H ~> m or kg m-2]
1477		integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1478		integer, intent(in) :: i !< The i-index to work on
1479		type(entrain_diffusive_CS), intent(in) :: CS !< This module's control structure.
1480		real, dimension(SZI_(G)), intent(inout) :: ea_kb !< The entrainment from above by the layer below
1481		!! the buffer layer (i.e. layer kb) [H ~> m or kg m-2].
1482		real, optional, intent(in) :: tol_in !< A tolerance for the iterative determination
1483		!! of the entrainment [H ~> m or kg m-2].
1484
1485		real :: max_ea, min_ea ! Bounds on the estimated entraiment [H ~> m or kg m-2]
1486		real :: err, err_min, err_max ! Errors in the mass flux balance [H R ~> kg m-2 or kg2 m-5]
1487		real :: derr_dea ! The change in error with the change in ea [R ~> kg m-3]
1488		real :: val ! An estimate mass flux [H R ~> kg m-2 or kg2 m-5]
1489		real :: tolerance, tol1 ! Tolerances for the determination of the entrainment [H ~> m or kg m-2]
1490		real :: ea_prev ! A previous estimate of ea_kb [H ~> m or kg m-2]
1491		real :: dS_kbp1 ! The density difference between two interior layers [R ~> kg m-3]
1492	1896	real :: dS_kb(SZI_(G)) ! The limited potential density difference across the interface
1493		! between the bottommost buffer layer and the topmost interior layer [R ~> kg m-3]
1494	1896	real :: maxF(SZI_(G)) ! The maximum value of F (the density flux divided by density
1495		! differences) found in the range min_ent < ent < max_ent [H ~> m or kg m-2].
1496	1896	real :: ent_maxF(SZI_(G)) ! The value of entrainment that gives maxF [H ~> m or kg m-2]
1497	1896	real :: zeros(SZI_(G)) ! An array of zero entrainments [H ~> m or kg m-2]
1498	948	real :: ddSkb_dE(SZI_(G)) ! The partial derivative of dS_kb with ea_kb [R H-1 ~> kg m-4 or m-1]
1499		logical :: bisect_next, Newton ! These indicate what method the next iteration should use
1500		integer :: it
1501		integer, parameter :: MAXIT = 30
1502
1503	948	dS_kbp1 = Sref(i,kmb+2) - Sref(i,kmb+1)
1504	948	max_ea = ea_kb(i) ; min_ea = 0.0
1505	948	val = dS_kbp1 * F_kb(i)
1506	948	err_min = -val
1507
1508	948	tolerance = CS%Tolerance_Ent
1509	948	if (present(tol_in)) tolerance = tol_in
1510	948	bisect_next = .true.
1511
1512		call determine_dSkb(h_bl, Sref, Ent_bl, ea_kb, i, i, kmb, G, GV, .true., &
1513	948	dS_kb, ddSkb_dE)
1514
1515	948	err = dS_kb(i) * ea_kb(i) - val
1516	948	derr_dea = dS_kb(i) + ddSkb_dE(i) * ea_kb(i)
1517		! Return if Newton's method on the first guess would give a tolerably small
1518		! change in the value of ea_kb.
1519	948	if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea))) return
1520
1521	948	if (err == 0.0) then ; return ! The exact solution on the first guess...
1522	948	elseif (err > 0.0) then ! The root is properly bracketed.
1523	948	max_ea = ea_kb(i) ; err_max = err
1524		! Use Newton's method (if it stays bounded) or the false position method
1525		! to find the next value.
1526	948	if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i) - min_ea) > err) .and. &
1527		(derr_dea(max_ea - ea_kb(i)) > -1.0err)) then
1528	430	ea_kb(i) = ea_kb(i) - err / derr_dea
1529		else ! Use the bisection for the next guess.
1530	518	ea_kb(i) = 0.5*(max_ea+min_ea)
1531		endif
1532		else
1533		! Try to bracket the root first. If unable to bracket the root, return
1534		! the maximum.
1535	0	zeros(i) = 0.0
1536		call find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, zeros, ea_kb, &
1537	0	kmb, i, i, G, GV, CS, maxF, ent_maxF, F_thresh=F_kb)
1538	0	err_max = dS_kbp1 * maxF(i) - val
1539		! If err_max is negative, there is no good solution, so use the maximum
1540		! value of F in the valid range.
1541	0	if (err_max <= 0.0) then
1542	0	ea_kb(i) = ent_maxF(i) ; return
1543		else
1544	0	max_ea = ent_maxF(i)
1545	0	ea_kb(i) = 0.5*(max_ea+min_ea) ! Use bisection for the next guess.
1546		endif
1547		endif
1548
1549		! Exit if the range between max_ea and min_ea already acceptable.
1550		! if (abs(max_ea - min_ea) < 0.1*tolerance) return
1551
1552	2261	do it = 1, MAXIT
1553		call determine_dSkb(h_bl, Sref, Ent_bl, ea_kb, i, i, kmb, G, GV, .true., &
1554	2261	dS_kb, ddSkb_dE)
1555
1556	2261	err = dS_kb(i) * ea_kb(i) - val
1557	2261	derr_dea = dS_kb(i) + ddSkb_dE(i) * ea_kb(i)
1558
1559	2261	ea_prev = ea_kb(i)
1560		! Use Newton's method or the false position method to find the next value.
1561	2261	Newton = .false.
1562	2261	if (err > 0.0) then
1563	1460	max_ea = ea_kb(i) ; err_max = err
1564	1460	if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) Newton = .true.
1565		else
1566	801	min_ea = ea_kb(i) ; err_min = err
1567	801	if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) Newton = .true.
1568		endif
1569
1570	2261	if (Newton) then
1571	1541	ea_kb(i) = ea_kb(i) - err / derr_dea
1572	720	elseif (bisect_next) then ! Use bisection to reduce the range.
1573	467	ea_kb(i) = 0.5*(max_ea+min_ea)
1574	467	bisect_next = .false.
1575		else ! Use the false-position method for the next guess.
1576	253	ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))
1577	253	bisect_next = .true.
1578		endif
1579
1580	2261	tol1 = tolerance ; if (err > 0.0) tol1 = 0.099*tolerance
1581	2261	if (dS_kb(i) <= dS_kbp1) then
1582	2189	if (abs(ea_kb(i) - ea_prev) <= tol1) return
1583		else
1584	72	if (dS_kbp1abs(ea_kb(i) - ea_prev) <= dS_kb(i)tol1) return
1585		endif
1586		enddo
1587
1588		end subroutine F_kb_to_ea_kb
1589
1590
1591		!> This subroutine determines the entrainment from above by the top interior
1592		!! layer (labeled kb elsewhere) given an entrainment by the layer below it,
1593		!! constrained to be within the provided bounds.
1594	5004	subroutine determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &
1595	3336	min_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, Ent, &
1596	3288	error, err_min_eakb0, err_max_eakb0, F_kb, dFdfm_kb)
1597		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1598		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1599		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h_bl !< Layer thickness, with the top interior
1600		!! layer at k-index kmb+1 [H ~> m or kg m-2].
1601		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Sref !< The coordinate reference potential
1602		!! density, with the value of the
1603		!! topmost interior layer at layer
1604		!! kmb+1 [R ~> kg m-3].
1605		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Ent_bl !< The average entrainment upward and
1606		!! downward across each interface around
1607		!! the buffer layers [H ~> m or kg m-2].
1608		real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in
1609		!! reference potential density across
1610		!! the base of the uppermost interior
1611		!! layer [R-1 ~> m3 kg-1].
1612		real, dimension(SZI_(G)), intent(in) :: dtKd_kb !< The diapycnal diffusivity in the top
1613		!! interior layer times the time step
1614		!! [H2 ~> m2 or kg2 m-4].
1615		real, dimension(SZI_(G)), intent(in) :: ea_kbp1 !< The entrainment from above by layer
1616		!! kb+1 [H ~> m or kg m-2].
1617		real, dimension(SZI_(G)), intent(in) :: min_eakb !< The minimum permissible rate of
1618		!! entrainment [H ~> m or kg m-2].
1619		real, dimension(SZI_(G)), intent(in) :: max_eakb !< The maximum permissible rate of
1620		!! entrainment [H ~> m or kg m-2].
1621		integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1622		integer, intent(in) :: is !< The start of the i-index range to work on.
1623		integer, intent(in) :: ie !< The end of the i-index range to work on.
1624		logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which
1625		!! i-points to work on.
1626		type(entrain_diffusive_CS), intent(in) :: CS !< This module's control structure.
1627		real, dimension(SZI_(G)), intent(inout) :: Ent !< The entrainment rate of the uppermost
1628		!! interior layer [H ~> m or kg m-2].
1629		!! The input value is the first guess.
1630		real, dimension(SZI_(G)), optional, intent(out) :: error !< The error (locally defined in this
1631		!! routine) associated with the returned
1632		!! solution [H2 ~> m2 or kg2 m-4]
1633		real, dimension(SZI_(G)), optional, intent(in) :: err_min_eakb0 !< The errors (locally defined)
1634		!! associated with min_eakb when ea_kbp1 = 0,
1635		!! returned from a previous call to this
1636		!! subroutine [H2 ~> m2 or kg2 m-4].
1637		real, dimension(SZI_(G)), optional, intent(in) :: err_max_eakb0 !< The errors (locally defined)
1638		!! associated with min_eakb when ea_kbp1 = 0,
1639		!! returned from a previous call to this
1640		!! subroutine [H2 ~> m2 or kg2 m-4].
1641		real, dimension(SZI_(G)), optional, intent(out) :: F_kb !< The entrainment from below by the
1642		!! uppermost interior layer
1643		!! corresponding to the returned
1644		!! value of Ent [H ~> m or kg m-2].
1645		real, dimension(SZI_(G)), optional, intent(out) :: dFdfm_kb !< The partial derivative of F_kb with
1646		!! ea_kbp1 [nondim].
1647
1648		! This subroutine determines the entrainment from above by the top interior
1649		! layer (labeled kb elsewhere) given an entrainment by the layer below it,
1650		! constrained to be within the provided bounds.
1651
1652		! Local variables
1653		real, dimension(SZI_(G)) :: &
1654	3336	dS_kb, & ! The coordinate-density difference between the
1655		! layer kb and deepest buffer layer, limited to
1656		! ensure that it is positive [R ~> kg m-3].
1657	3336	dS_Lay, & ! The coordinate-density difference across layer
1658		! kb, limited to ensure that it is positive and not
1659		! too much bigger than dS_kb or dS_kbp1 [R ~> kg m-3].
1660	3336	ddSkb_dE, ddSlay_dE, & ! The derivatives of dS_kb and dS_Lay with E
1661		! [R H-1 ~> kg m-4 or m-1].
1662	3336	derror_dE, & ! The derivative of err with E [H ~> m or kg m-2].
1663	3336	err, & ! The "error" whose zero is being sought [H2 ~> m2 or kg2 m-4].
1664	3336	E_min, E_max, & ! The minimum and maximum values of E [H ~> m or kg m-2].
1665	3336	error_minE, error_maxE ! err when E = E_min or E = E_max [H2 ~> m2 or kg2 m-4].
1666		real :: err_est ! An estimate of what err will be [H2 ~> m2 or kg2 m-4].
1667		real :: eL ! 1 or 0, depending on whether increases in E lead
1668		! to decreases in the entrainment from below by the
1669		! deepest buffer layer [nondim].
1670		real :: fa ! Temporary variable used to calculate err [nondim].
1671		real :: fk ! Temporary variable used to calculate err [H2 ~> m2 or kg2 m-4].
1672		real :: fm, fr ! Temporary variables used to calculate err [H ~> m or kg m-2].
1673		real :: tolerance ! The tolerance within which E must be converged [H ~> m or kg m-2].
1674		real :: E_prev ! The previous value of E [H ~> m or kg m-2].
1675	3336	logical, dimension(SZI_(G)) :: false_position ! If true, the false position
1676		! method might be used for the next iteration.
1677	3336	logical, dimension(SZI_(G)) :: redo_i ! If true, more work is needed on this column.
1678		logical :: do_any
1679		real :: large_err ! A large error measure [H2 ~> m2 or kg2 m-4].
1680		integer :: i, it
1681		integer, parameter :: MAXIT = 30
1682
1683	1668	if (.not.CS%bulkmixedlayer) then
1684		call MOM_error(FATAL, "determine_Ea_kb should not be called "//&
1685	0	"unless BULKMIXEDLAYER is defined.")
1686		endif
1687	1668	tolerance = CS%Tolerance_Ent
1688	1668	large_err = GV%m_to_H*2 1.0e30
1689
1690	244224	do i=is,ie ; redo_i(i) = do_i(i) ; enddo
1691
1692	244224	do i=is,ie ; if (do_i(i)) then
1693		! The first guess of Ent was the value from the previous iteration.
1694
1695		! These were previously calculated and provide good limits and estimates
1696		! of the errors there. By construction the errors increase with R*ea_kbp1.
1697	185948	E_min(i) = min_eakb(i) ; E_max(i) = max_eakb(i)
1698	185948	error_minE(i) = -large_err ; error_maxE(i) = large_err
1699	185948	false_position(i) = .true. ! Used to alternate between false_position and
1700		! bisection when Newton's method isn't working.
1701	185948	if (present(err_min_eakb0)) error_minE(i) = err_min_eakb0(i) - E_min(i) * ea_kbp1(i)
1702	185948	if (present(err_max_eakb0)) error_maxE(i) = err_max_eakb0(i) - E_max(i) * ea_kbp1(i)
1703
1704	185948	if ((error_maxE(i) <= 0.0) .or. (error_minE(i) >= 0.0)) then
1705		! The root is not bracketed and one of the limiting values should be used.
1706	39890	if (error_maxE(i) <= 0.0) then
1707		! The errors decrease with E*ea_kbp1, so E_max is the best solution.
1708	1530	Ent(i) = E_max(i) ; err(i) = error_maxE(i)
1709		else ! error_minE >= 0 is equivalent to ea_kbp1 = 0.0.
1710	38360	Ent(i) = E_min(i) ; err(i) = error_minE(i)
1711		endif
1712	39890	derror_dE(i) = 0.0
1713	39890	redo_i(i) = .false.
1714		endif
1715		endif ; enddo ! End of i-loop
1716
1717	4882	do it = 1,MAXIT
1718	704616	do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo
1719	4882	if (.not.do_any) exit
1720		call determine_dSkb(h_bl, Sref, Ent_bl, Ent, is, ie, kmb, G, GV, .true., dS_kb, &
1721	3214	ddSkb_dE, dS_lay, ddSlay_dE, do_i_in=redo_i)
1722	460392	do i=is,ie ; if (redo_i(i)) then
1723		! The correct root is bracketed between E_min and E_max.
1724		! Note the following limits: Ent >= 0 ; fa > 1 ; fk > 0
1725	288466	eL = 0.0 ; if (2.0*Ent_bl(i,Kmb+1) >= Ent(i)) eL = 1.0
1726	288466	fa = (1.0 + eL) + dS_kb(i)*I_dSkbp1(i)
1727	288466	fk = dtKd_kb(i) * (dS_Lay(i)/dS_kb(i))
1728	288466	fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + eL2.0Ent_bl(i,Kmb+1)
1729	288466	if (fm > -GV%Angstrom_H) fm = fm + GV%Angstrom_H ! This could be smooth if need be.
1730	288466	err(i) = (fa * Ent(i)*2 - fm Ent(i)) - fk
1731		derror_dE(i) = ((2.0fa + (ddSkb_dE(i)I_dSkbp1(i))Ent(i))Ent(i) - fm) - &
1732	288466	dtKd_kb(i) * (ddSlay_dE(i)dS_kb(i) - ddSkb_dE(i)dS_Lay(i))/(dS_kb(i)**2)
1733
1734	288466	if (err(i) == 0.0) then
1735	0	redo_i(i) = .false. ; cycle
1736	288466	elseif (err(i) > 0.0) then
1737	207547	E_max(i) = Ent(i) ; error_maxE(i) = err(i)
1738		else
1739	80919	E_min(i) = Ent(i) ; error_minE(i) = err(i)
1740		endif
1741
1742	288466	E_prev = Ent(i)
1743	288466	if ((it == 1) .or. (derror_dE(i) <= 0.0)) then
1744		! Assuming that the coefficients of the quadratic equation are correct
1745		! will usually give a very good first guess. Also, if derror_dE < 0.0,
1746		! R is on the wrong side of the approximate parabola. In either case,
1747		! try assuming that the error is approximately a parabola and solve.
1748	146059	fr = sqrt(fm*2 + 4.0fa*fk)
1749	146059	if (fm >= 0.0) then
1750	14598	Ent(i) = (fm + fr) / (2.0 * fa)
1751		else
1752	131461	Ent(i) = (2.0 * fk) / (fr - fm)
1753		endif
1754		! But make sure that the root stays bracketed, bisecting if needed.
1755	146059	if ((Ent(i) > E_max(i)) .or. (Ent(i) < E_min(i))) &
1756	47199	Ent(i) = 0.5*(E_max(i) + E_min(i))
1757	142407	elseif (((E_max(i)-Ent(i))*derror_dE(i) > -err(i)) .and. &
1758		((Ent(i)-E_min(i))*derror_dE(i) > err(i)) ) then
1759		! Use Newton's method for the next estimate, provided it will
1760		! remain bracketed between Rmin and Rmax.
1761	95210	Ent(i) = Ent(i) - err(i) / derror_dE(i)
1762	47197	elseif (false_position(i) .and. &
1763		(error_maxE(i) - error_minE(i) < 0.9*large_err)) then
1764		! Use the false position method if there are decent error estimates.
1765		Ent(i) = E_min(i) + (E_max(i)-E_min(i)) * &
1766	0	(-error_minE(i)/(error_maxE(i) - error_minE(i)))
1767	0	false_position(i) = .false.
1768		else ! Bisect as a last resort or if the false position method was used last.
1769	47197	Ent(i) = 0.5*(E_max(i) + E_min(i))
1770	47197	false_position(i) = .true.
1771		endif
1772
1773	288466	if (abs(E_prev - Ent(i)) < tolerance) then
1774	193256	err_est = err(i) + (Ent(i) - E_prev) * derror_dE(i)
1775	193256	if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &
1776	146058	(abs(err_est) < abs(tolerance*derror_dE(i)))) redo_i(i) = .false.
1777		endif
1778
1779		endif ; enddo ! End of i-loop
1780		enddo ! End of iterations to determine Ent(i).
1781
1782		! Update the value of dS_kb for consistency with Ent.
1783	1668	if (present(F_kb) .or. present(dFdfm_kb)) &
1784		call determine_dSkb(h_bl, Sref, Ent_bl, Ent, is, ie, kmb, G, GV, .true., &
1785	720	dS_kb, do_i_in=do_i)
1786
1787	243276	if (present(F_kb)) then ; do i=is,ie ; if (do_i(i)) then
1788	185000	F_kb(i) = Ent(i) * (dS_kb(i) * I_dSkbp1(i))
1789		endif ; enddo ; endif
1790	132216	if (present(error)) then ; do i=is,ie ; if (do_i(i)) then
1791	93448	error(i) = err(i)
1792		endif ; enddo ; endif
1793	113676	if (present(dFdfm_kb)) then ; do i=is,ie ; if (do_i(i)) then
1794		! derror_dE and ddSkb_dE are _not_ recalculated here, since dFdfm_kb is
1795		! only used in Newton's method, and slightly increasing the accuracy of the
1796		! estimate is unlikely to speed convergence.
1797	92500	if (derror_dE(i) > 0.0) then
1798		dFdfm_kb(i) = ((dS_kb(i) + Ent(i) * ddSkb_dE(i)) * I_dSkbp1(i)) * &
1799	52610	(Ent(i) / derror_dE(i))
1800		else ! Use Adcroft's division by 0 convention.
1801	39890	dFdfm_kb(i) = 0.0
1802		endif
1803		endif ; enddo ; endif
1804
1805	3648	end subroutine determine_Ea_kb
1806
1807		!> Maximize F = entds_kbI_dSkbp1 in the range min_ent < ent < max_ent.
1808	720	subroutine find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, &
1809	540	kmb, is, ie, G, GV, CS, maxF, ent_maxF, do_i_in, &
1810	180	F_lim_maxent, F_thresh)
1811		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1812		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1813		real, dimension(SZI_(G),SZK_(GV)), &
1814		intent(in) :: h_bl !< Layer thickness [H ~> m or kg m-2]
1815		real, dimension(SZI_(G),SZK_(GV)), &
1816		intent(in) :: Sref !< Reference potential density [R ~> kg m-3].
1817		real, dimension(SZI_(G),SZK_(GV)), &
1818		intent(in) :: Ent_bl !< The average entrainment upward and
1819		!! downward across each interface around
1820		!! the buffer layers [H ~> m or kg m-2].
1821		real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in
1822		!! reference potential density across the
1823		!! base of the uppermost interior layer
1824		!! [R-1 ~> m3 kg-1].
1825		real, dimension(SZI_(G)), intent(in) :: min_ent_in !< The minimum value of ent to search,
1826		!! [H ~> m or kg m-2].
1827		real, dimension(SZI_(G)), intent(in) :: max_ent_in !< The maximum value of ent to search,
1828		!! [H ~> m or kg m-2].
1829		integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1830		integer, intent(in) :: is !< The start of the i-index range to work on.
1831		integer, intent(in) :: ie !< The end of the i-index range to work on.
1832		type(entrain_diffusive_CS), intent(in) :: CS !< This module's control structure.
1833		real, dimension(SZI_(G)), intent(out) :: maxF !< The maximum value of F
1834		!! = entds_kbI_dSkbp1 found in the range
1835		!! min_ent < ent < max_ent [H ~> m or kg m-2].
1836		real, dimension(SZI_(G)), &
1837		optional, intent(out) :: ent_maxF !< The value of ent at that maximum [H ~> m or kg m-2].
1838		logical, dimension(SZI_(G)), &
1839		optional, intent(in) :: do_i_in !< A logical array indicating which columns
1840		!! to work on.
1841		real, dimension(SZI_(G)), &
1842		optional, intent(out) :: F_lim_maxent !< If present, do not apply the limit in
1843		!! finding the maximum value, but return the
1844		!! limited value at ent=max_ent_in in this
1845		!! array [H ~> m or kg m-2].
1846		real, dimension(SZI_(G)), &
1847		optional, intent(in) :: F_thresh !< If F_thresh is present, return the first value
1848		!! found that has F > F_thresh [H ~> m or kg m-2], or
1849		!! the maximum root if it is absent.
1850
1851		! Maximize F = entds_kbI_dSkbp1 in the range min_ent < ent < max_ent.
1852		! ds_kb may itself be limited to positive values in determine_dSkb, which gives
1853		! the prospect of two local maxima in the range - one at max_ent_in with that
1854		! minimum value of ds_kb, and the other due to the unlimited (potentially
1855		! negative) value. It is faster to find the true maximum by first finding the
1856		! unlimited maximum and comparing it to the limited value at max_ent_in.
1857		real, dimension(SZI_(G)) :: &
1858	720	ent, & ! The updated estimate of the entrainment [H ~> m or kg m-2]
1859	720	minent, maxent, ent_best, & ! Various previous estimates of the entrainment [H ~> m or kg m-2]
1860	720	F_max_ent_in, & ! The value of F that gives the input maximum value of ent [H ~> m or kg m-2]
1861	720	F_maxent, F_minent, F, F_best, & ! Various estimates of F [H ~> m or kg m-2]
1862	720	dF_dent, dF_dE_max, dF_dE_min, dF_dE_best, & ! Various derivatives of F with ent [nondim]
1863	720	dS_kb, & ! The density difference across the interface between the bottommost
1864		! buffer layer and the topmost interior layer [R ~> kg m-3]
1865	720	dS_kb_lim, dS_anom_lim, & ! Various limits on dS_kb [R ~> kg m-3]
1866	720	ddSkb_dE, & ! The partial derivative of dS_kb with ent [R H-1 ~> kg m-4 or m-1].
1867	720	chg_prev, chg_pre_prev ! Changes in estimates of the entrainment from previous iterations [H ~> m or kg m-2]
1868		real :: dF_dE_mean, maxslope, minslope ! Various derivatives of F with ent [nondim]
1869		real :: tolerance ! The tolerance within which ent must be converged [H ~> m or kg m-2]
1870		real :: ratio_select_end, rat ! Fractional changes in the value of ent to use for the next iteration
1871		! relative to its bounded range [nondim]
1872		real :: max_chg, min_chg, chg1, chg2, chg ! Changes in entrainment estimates [H ~> m or kg m-2]
1873	720	logical, dimension(SZI_(G)) :: do_i, last_it, need_bracket, may_use_best
1874		logical :: doany, OK1, OK2, bisect, new_min_bound
1875		integer :: i, it, is1, ie1
1876		integer, parameter :: MAXIT = 20
1877
1878	360	tolerance = CS%Tolerance_Ent
1879
1880	360	if (present(do_i_in)) then
1881	129960	do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1882		else
1883	0	do i=is,ie ; do_i(i) = .true. ; enddo
1884		endif
1885
1886		! The most likely value is at max_ent.
1887		call determine_dSkb(h_bl, Sref, Ent_bl, max_ent_in, is, ie, kmb, G, GV, .false., &
1888	360	dS_kb, ddSkb_dE, dS_anom_lim=dS_anom_lim)
1889	360	ie1 = is-1 ; doany = .false.
1890	129960	do i=is,ie
1891	129600	dS_kb_lim(i) = dS_kb(i) + dS_anom_lim(i)
1892	129600	F_max_ent_in(i) = max_ent_in(i)dS_kb_lim(i)I_dSkbp1(i)
1893	129600	maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)
1894	129960	if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i))) then
1895	83350	F_best(i) = max_ent_in(i)dS_kb(i)I_dSkbp1(i)
1896	83350	ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)
1897	83350	do_i(i) = .false.
1898		else
1899	46250	F_maxent(i) = maxent(i) * dS_kb(i) * I_dSkbp1(i)
1900	46250	dF_dE_max(i) = (dS_kb(i) + maxent(i)ddSkb_dE(i)) I_dSkbp1(i)
1901	46250	doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.
1902		endif
1903		enddo
1904
1905	360	if (doany) then
1906	64980	ie1 = is-1 ; do i=is,ie ; if (do_i(i)) ie1 = i ; enddo
1907	60582	do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
1908		! Find the value of F and its derivative at min_ent.
1909		call determine_dSkb(h_bl, Sref, Ent_bl, minent, is1, ie1, kmb, G, GV, .false., &
1910	180	dS_kb, ddSkb_dE, do_i_in=do_i)
1911	56184	do i=is1,ie1 ; if (do_i(i)) then
1912	46250	F_minent(i) = minent(i) * dS_kb(i) * I_dSkbp1(i)
1913	46250	dF_dE_min(i) = (dS_kb(i) + minent(i)ddSkb_dE(i)) I_dSkbp1(i)
1914		endif ; enddo
1915
1916	180	ratio_select_end = 0.9
1917	337	do it=1,MAXIT
1918	337	ratio_select_end = 0.5*ratio_select_end
1919	103775	do i=is1,ie1 ; if (do_i(i)) then
1920	79552	if (need_bracket(i)) then
1921	79552	dF_dE_mean = (F_maxent(i) - F_minent(i)) / (maxent(i) - minent(i))
1922	79552	maxslope = MAX(dF_dE_mean, dF_dE_min(i), dF_dE_max(i))
1923	79552	minslope = MIN(dF_dE_mean, dF_dE_min(i), dF_dE_max(i))
1924	79552	if (F_minent(i) >= F_maxent(i)) then
1925	8526	if (dF_dE_min(i) > 0.0) then ; rat = 0.02 ! A small step should bracket the solution.
1926	8526	elseif (maxslope < ratio_select_end*minslope) then
1927		! The maximum of F is at minent.
1928	2681	F_best(i) = F_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0
1929	2681	do_i(i) = .false.
1930	5845	else ; rat = 0.382 ; endif ! Use the golden ratio
1931		else
1932	71026	if (dF_dE_max(i) < 0.0) then ; rat = 0.98 ! A small step should bracket the solution.
1933	71026	elseif (minslope > ratio_select_end*maxslope) then
1934		! The maximum of F is at maxent.
1935	43569	F_best(i) = F_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0
1936	43569	do_i(i) = .false.
1937	27457	else ; rat = 0.618 ; endif ! Use the golden ratio
1938		endif
1939
1940	79552	if (rat >= 0.0) ent(i) = ratmaxent(i) + (1.0-rat)minent(i)
1941	79552	if (((maxent(i) - minent(i)) < tolerance) .or. (it==MAXIT)) &
1942	0	last_it(i) = .true.
1943		else ! The maximum is bracketed by minent, ent_best, and maxent.
1944	0	chg1 = 2.0*(maxent(i) - minent(i)) ; chg2 = chg1
1945	0	if (dF_dE_best(i) > 0) then
1946	0	max_chg = maxent(i) - ent_best(i) ; min_chg = 0.0
1947		else
1948	0	max_chg = 0.0 ; min_chg = minent(i) - ent_best(i) ! < 0
1949		endif
1950	0	if (max_chg - min_chg < 2.0*tolerance) last_it(i) = .true.
1951	0	if (dF_dE_max(i) /= dF_dE_best(i)) &
1952		chg1 = (maxent(i) - ent_best(i))*dF_dE_best(i) / &
1953	0	(dF_dE_best(i) - dF_dE_max(i))
1954	0	if (dF_dE_min(i) /= dF_dE_best(i)) &
1955		chg2 = (minent(i) - ent_best(i))*dF_dE_best(i) / &
1956	0	(dF_dE_best(i) - dF_dE_min(i))
1957	0	OK1 = ((chg1 < max_chg) .and. (chg1 > min_chg))
1958	0	OK2 = ((chg2 < max_chg) .and. (chg2 > min_chg))
1959	0	if (.not.(OK1 .or. OK2)) then ; bisect = .true. ; else
1960	0	if (OK1 .and. OK2) then ! Take the acceptable smaller change.
1961	0	chg = chg1 ; if (abs(chg2) < abs(chg1)) chg = chg2
1962	0	elseif (OK1) then ; chg = chg1
1963	0	else ; chg = chg2 ; endif
1964	0	if (abs(chg) > 0.5*abs(chg_pre_prev(i))) then ; bisect = .true.
1965	0	else ; bisect = .false. ; endif
1966		endif
1967	0	chg_pre_prev(i) = chg_prev(i)
1968	0	if (bisect) then
1969	0	if (dF_dE_best(i) > 0.0) then
1970	0	ent(i) = 0.5*(maxent(i) + ent_best(i))
1971	0	chg_prev(i) = 0.5*(maxent(i) - ent_best(i))
1972		else
1973	0	ent(i) = 0.5*(minent(i) + ent_best(i))
1974	0	chg_prev(i) = 0.5*(minent(i) - ent_best(i))
1975		endif
1976		else
1977	0	if (abs(chg) < tolerance) chg = SIGN(tolerance,chg)
1978	0	ent(i) = ent_best(i) + chg
1979	0	chg_prev(i) = chg
1980		endif
1981		endif
1982		endif ; enddo
1983
1984	337	if (mod(it,3) == 0) then ! Re-determine the loop bounds.
1985	2116	ie1 = is-1 ; do i=is1,ie ; if (do_i(i)) ie1 = i ; enddo
1986	1036	do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
1987		endif
1988
1989		call determine_dSkb(h_bl, Sref, Ent_bl, ent, is1, ie1, kmb, G, GV, .false., &
1990	337	dS_kb, ddSkb_dE, do_i_in=do_i)
1991	102695	do i=is1,ie1 ; if (do_i(i)) then
1992	33302	F(i) = ent(i)dS_kb(i)I_dSkbp1(i)
1993	33302	dF_dent(i) = (dS_kb(i) + ent(i)ddSkb_dE(i)) I_dSkbp1(i)
1994		endif ; enddo
1995
1996	337	if (present(F_thresh)) then ; do i=is1,ie1 ; if (do_i(i)) then
1997	0	if (F(i) >= F_thresh(i)) then
1998	0	F_best(i) = F(i) ; ent_best(i) = ent(i) ; do_i(i) = .false.
1999		endif
2000		endif ; enddo ; endif
2001
2002	337	doany = .false.
2003	102695	do i=is1,ie1 ; if (do_i(i)) then
2004	33302	if (.not.last_it(i)) doany = .true.
2005	33302	if (last_it(i)) then
2006	0	if (need_bracket(i)) then
2007	0	if ((F(i) > F_maxent(i)) .and. (F(i) > F_minent(i))) then
2008	0	F_best(i) = F(i) ; ent_best(i) = ent(i)
2009	0	elseif (F_maxent(i) > F_minent(i)) then
2010	0	F_best(i) = F_maxent(i) ; ent_best(i) = maxent(i)
2011		else
2012	0	F_best(i) = F_minent(i) ; ent_best(i) = minent(i)
2013		endif
2014	0	elseif (F(i) > F_best(i)) then
2015	0	F_best(i) = F(i) ; ent_best(i) = ent(i)
2016		endif
2017	0	do_i(i) = .false.
2018	33302	elseif (need_bracket(i)) then
2019	33302	if ((F(i) > F_maxent(i)) .and. (F(i) > F_minent(i))) then
2020	0	need_bracket(i) = .false. ! The maximum is now bracketed.
2021	0	chg_prev(i) = (maxent(i) - minent(i))
2022	0	chg_pre_prev(i) = 2.0*chg_prev(i)
2023	0	ent_best(i) = ent(i) ; F_best(i) = F(i) ; dF_dE_best(i) = dF_dent(i)
2024	33302	elseif ((F(i) <= F_maxent(i)) .and. (F(i) > F_minent(i))) then
2025	27457	new_min_bound = .true. ! We have a new minimum bound.
2026	5845	elseif ((F(i) <= F_maxent(i)) .and. (F(i) > F_minent(i))) then
2027	0	new_min_bound = .false. ! We have a new maximum bound.
2028		else ! This case would bracket a minimum. Weird.
2029		! Unless the derivative indicates that there is a maximum near the
2030		! lower bound, try keeping the end with the larger value of F
2031		! in a tie keep the minimum as the answer here will be compared
2032		! with the maximum input value later.
2033	5845	new_min_bound = .true.
2034	5845	if (dF_dE_min(i) > 0.0 .or. (F_minent(i) >= F_maxent(i))) &
2035	5845	new_min_bound = .false.
2036		endif
2037	33302	if (need_bracket(i)) then ! Still not bracketed.
2038	33302	if (new_min_bound) then
2039	27457	minent(i) = ent(i) ; F_minent(i) = F(i) ; dF_dE_min(i) = dF_dent(i)
2040		else
2041	5845	maxent(i) = ent(i) ; F_maxent(i) = F(i) ; dF_dE_max(i) = dF_dent(i)
2042		endif
2043		endif
2044		else ! The root was previously bracketed.
2045	0	if (F(i) >= F_best(i)) then ! There is a new maximum.
2046	0	if (ent(i) > ent_best(i)) then ! Replace minent with ent_prev.
2047	0	minent(i) = ent_best(i) ; F_minent(i) = F_best(i) ; dF_dE_min(i) = dF_dE_best(i)
2048		else ! Replace maxent with ent_best.
2049	0	maxent(i) = ent_best(i) ; F_maxent(i) = F_best(i) ; dF_dE_max(i) = dF_dE_best(i)
2050		endif
2051	0	ent_best(i) = ent(i) ; F_best(i) = F(i) ; dF_dE_best(i) = dF_dent(i)
2052		else
2053	0	if (ent(i) < ent_best(i)) then ! Replace the minent with ent.
2054	0	minent(i) = ent(i) ; F_minent(i) = F(i) ; dF_dE_min(i) = dF_dent(i)
2055		else ! Replace maxent with ent_prev.
2056	0	maxent(i) = ent(i) ; F_maxent(i) = F(i) ; dF_dE_max(i) = dF_dent(i)
2057		endif
2058		endif
2059	0	if ((maxent(i) - minent(i)) <= tolerance) do_i(i) = .false. ! Done.
2060		endif ! need_bracket.
2061		endif ; enddo
2062	337	if (.not.doany) exit
2063		enddo
2064		endif
2065
2066	360	if (present(F_lim_maxent)) then
2067		! Return the unlimited maximum in maxF, and the limited value of F at maxent.
2068	64980	do i=is,ie
2069	64800	maxF(i) = F_best(i)
2070	64800	F_lim_maxent(i) = F_max_ent_in(i)
2071	64980	if (present(ent_maxF)) ent_maxF(i) = ent_best(i)
2072		enddo
2073		else
2074		! Now compare the two? potential maxima using the limited value of dF_kb.
2075	180	doany = .false.
2076	64980	do i=is,ie
2077	64800	may_use_best(i) = (ent_best(i) /= max_ent_in(i))
2078	64980	if (may_use_best(i)) doany = .true.
2079		enddo
2080	180	if (doany) then
2081		! For efficiency, could save previous value of dS_anom_lim_best?
2082	0	call determine_dSkb(h_bl, Sref, Ent_bl, ent_best, is, ie, kmb, G, GV, .true., dS_kb_lim)
2083	0	do i=is,ie
2084	0	F_best(i) = ent_best(i)dS_kb_lim(i)I_dSkbp1(i)
2085		! The second test seems necessary because of roundoff differences that
2086		! can arise during compilation.
2087	0	if ((F_best(i) > F_max_ent_in(i)) .and. (may_use_best(i))) then
2088	0	maxF(i) = F_best(i)
2089	0	if (present(ent_maxF)) ent_maxF(i) = ent_best(i)
2090		else
2091	0	maxF(i) = F_max_ent_in(i)
2092	0	if (present(ent_maxF)) ent_maxF(i) = max_ent_in(i)
2093		endif
2094		enddo
2095		else
2096		! All of the maxima are at the maximum entrainment.
2097	64980	do i=is,ie ; maxF(i) = F_max_ent_in(i) ; enddo
2098	180	if (present(ent_maxF)) then
2099	0	do i=is,ie ; ent_maxF(i) = max_ent_in(i) ; enddo
2100		endif
2101		endif
2102		endif
2103
2104	900	end subroutine find_maxF_kb
2105
2106		!> This subroutine initializes the parameters and memory associated with the
2107		!! entrain_diffusive module.
2108	1	subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_read_params)
2109		type(time_type), intent(in) :: Time !< The current model time.
2110		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2111		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2112		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2113		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2114		!! parameters.
2115		type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
2116		!! output.
2117		type(entrain_diffusive_CS), intent(inout) :: CS !< Entrainment diffusion control structure
2118		logical, intent(in) :: just_read_params !< If true, this call will only read
2119		!! and log parameters without registering
2120		!! any diagnostics
2121
2122		! Local variables
2123		real :: dt ! The dynamics timestep, used here in the default for TOLERANCE_ENT [T ~> s]
2124		real :: Kd ! A diffusivity used in the default for TOLERANCE_ENT [Z2 T-1 ~> m2 s-1]
2125		! This include declares and sets the variable "version".
2126		# include "version_variable.h"
2127		character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name.
2128
2129	1	CS%initialized = .true.
2130	1	CS%diag => diag
2131
2132	1	CS%bulkmixedlayer = (GV%nkml > 0)
2133
2134		! Set default, read and log parameters
2135	1	if (.not.just_read_params) call log_version(param_file, mdl, version, "")
2136		call get_param(param_file, mdl, "MAX_ENT_IT", CS%max_ent_it, &
2137		"The maximum number of iterations that may be used to "//&
2138	1	"calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read_params)
2139		! In this module, KD is only used to set the default for TOLERANCE_ENT. [Z2 T-1 ~> m2 s-1]
2140	1	call get_param(param_file, mdl, "KD", Kd, units="m2 s-1", default=0.0, scale=US%m2_s_to_Z2_T)
2141		call get_param(param_file, mdl, "DT", dt, &
2142		"The (baroclinic) dynamics time step.", &
2143	1	units="s", scale=US%s_to_T, fail_if_missing=.true., do_not_log=just_read_params)
2144		call get_param(param_file, mdl, "TOLERANCE_ENT", CS%Tolerance_Ent, &
2145		"The tolerance with which to solve for entrainment values.", &
2146		units="m", default=US%Z_to_mMAX(100.0GV%Angstrom_Z,1.0e-4sqrt(dtKd)), scale=GV%m_to_H, &
2147	1	do_not_log=just_read_params)
2148		call get_param(param_file, mdl, "ENTRAIN_DIFFUSIVE_MAX_ENT", CS%max_Ent, &
2149		"A large ceiling on the maximum permitted amount of entrainment across each "//&
2150		"interface between the mixed and buffer layers within a timestep.", &
2151	1	units="m", default=1.0e4, scale=GV%m_to_H, do_not_log=.not.CS%bulkmixedlayer)
2152
2153	1	CS%Rho_sig_off = 1000.0*US%kg_m3_to_R
2154
2155	1	if (.not.just_read_params) then
2156		CS%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, Time, &
2157	1	'Diapycnal diffusivity as applied', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2158		CS%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, Time, &
2159		'Work actually done by diapycnal diffusion across each interface', &
2160	1	'W m-2', conversion=US%RZ3_T3_to_W_m2)
2161		endif
2162	1	end subroutine entrain_diffusive_init
2163
2164		!> \namespace mom_entrain_diffusive
2165		!!
2166		!! By Robert Hallberg, September 1997 - July 2000
2167		!!
2168		!! This file contains the subroutines that implement diapycnal
2169		!! mixing and advection in isopycnal layers. The main subroutine,
2170		!! calculate_entrainment, returns the entrainment by each layer
2171		!! across the interfaces above and below it. These are calculated
2172		!! subject to the constraints that no layers can be driven to negative
2173		!! thickness and that the each layer maintains its target density,
2174		!! using the scheme described in Hallberg (MWR 2000). There may or
2175		!! may not be a bulk mixed layer above the isopycnal layers.
2176		!! The solution is iterated until the change in the entrainment
2177		!! between successive iterations is less than some small tolerance.
2178		!!
2179		!! The dual-stream entrainment scheme of MacDougall and Dewar
2180		!! (JPO 1997) is used for combined diapycnal advection and diffusion,
2181		!! modified as described in Hallberg (MWR 2000) to be solved
2182		!! implicitly in time. Any profile of diffusivities may be used.
2183		!! Diapycnal advection is fundamentally the residual of diapycnal
2184		!! diffusion, so the fully implicit upwind differencing scheme that
2185		!! is used is entirely appropriate. The downward buoyancy flux in
2186		!! each layer is determined from an implicit calculation based on
2187		!! the previously calculated flux of the layer above and an estimated
2188		!! flux in the layer below. This flux is subject to the following
2189		!! conditions: (1) the flux in the top and bottom layers are
2190		!! set by the boundary conditions, and (2) no layer may be driven
2191		!! below an Angstrom thickness. If there is a bulk mixed layer, the
2192		!! mixed and buffer layers are treated as Eulerian layers, whose
2193		!! thicknesses only change due to entrainment by the interior layers.
2194
2195	0	end module MOM_entrain_diffusive