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	0	subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
58	0	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	0	dtKd ! The layer diapycnal diffusivity times the time step [H2 ~> m2 or kg2 m-4].
96		real, dimension(SZI_(G),SZK_(GV)+1) :: &
97	0	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	0	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	0	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	0	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	0	Fprev, &! The previous estimate of F [H ~> m or kg m-2].
106	0	dFdfm, &! The partial derivative of F with respect to changes in F of the
107		! neighboring layers. [nondim]
108	0	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	0	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	0	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	0	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	0	real :: b1(SZI_(G)) ! A variable used by the tridiagonal solver [H ~> m or kg m-2]
125	0	real :: c1(SZI_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]
126
127		real, dimension(SZI_(G)) :: &
128	0	htot, & ! The total thickness above or below a layer [H ~> m or kg m-2].
129	0	Rcv, & ! Value of the coordinate variable (potential density)
130		! based on the simulated T and S and P_Ref [R ~> kg m-3].
131	0	pres, & ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].
132	0	eakb, & ! The entrainment from above by the layer below the buffer
133		! layer (i.e. layer kb) [H ~> m or kg m-2].
134	0	ea_kbp1, & ! The entrainment from above by layer kb+1 [H ~> m or kg m-2].
135	0	eb_kmb, & ! The entrainment from below by the deepest buffer layer [H ~> m or kg m-2].
136	0	dS_kb, & ! The reference potential density difference across the
137		! interface between the buffer layers and layer kb [R ~> kg m-3].
138	0	dS_anom_lim, &! The amount by which dS_kb is reduced when limits are
139		! applied [R ~> kg m-3].
140	0	I_dSkbp1, & ! The inverse of the potential density difference across the
141		! interface below layer kb [R-1 ~> m3 kg-1].
142	0	dtKd_kb, & ! The diapycnal diffusivity in layer kb times the time step
143		! [H2 ~> m2 or kg2 m-4].
144	0	maxF_correct, & ! An amount by which to correct maxF due to excessive
145		! surface heat loss [H ~> m or kg m-2].
146	0	zeros, & ! An array of all zeros. (Usually used with [H ~> m or kg m-2].)
147	0	max_eakb, & ! The maximum value of eakb that might be realized [H ~> m or kg m-2].
148	0	min_eakb, & ! The minimum value of eakb that might be realized [H ~> m or kg m-2].
149	0	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	0	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	0	err_eakb0, & ! A value of error returned by determine_Ea_kb [H2 ~> m2 or kg2 m-4].
154	0	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	0	dFdfm_kb, & ! The partial derivative of F with fm [nondim]. See dFdfm.
157	0	maxF_kb, & ! The maximum value of F_kb that might be realized [H ~> m or kg m-2].
158	0	eakb_maxF, & ! The value of eakb that gives F_kb=maxF_kb [H ~> m or kg m-2].
159	0	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	0	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	0	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	0	ds_dsp1, & ! The coordinate variable (sigma-2) difference across an
170		! interface divided by the difference across the interface
171		! below it. [nondim]
172	0	dsp1_ds, & ! The inverse coordinate variable (sigma-2) difference
173		! across an interface times the difference across the
174		! interface above it. [nondim]
175	0	I2p2dsp1_ds, & ! 1 / (2 + 2 * ds_k+1 / ds_k). [nondim]
176	0	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	0	pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].
187	0	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	0	dRho_dT, & ! The partial derivative of potential density with temperature [R C-1 ~> kg m-3 degC-1]
190	0	dRho_dS, & ! The partial derivative of potential density with salinity [R S-1 ~> kg m-3 ppt-1]
191	0	dSpV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
192	0	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	0	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	0	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	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
220	0	Angstrom = GV%Angstrom_H
221	0	h_neglect = GV%H_subroundoff
222
223	0	if (.not. CS%initialized) call MOM_error(FATAL, &
224	0	"MOM_entrain_diffusive: Module must be initialized before it is used.")
225
226	0	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	0	tolerance = CS%Tolerance_Ent
238	0	kmb = GV%nk_rho_varies
239	0	K2 = max(kmb+1,2) ; kb_min = K2
240	0	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	0	kb(:) = 0
247	0	do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ; enddo
248		endif
249
250	0	if (CS%id_diff_work > 0) allocate(diff_work(G%isd:G%ied,G%jsd:G%jed,nz+1))
251	0	if (CS%id_Kd > 0) allocate(Kd_eff(G%isd:G%ied,G%jsd:G%jed,nz))
252
253	0	if (associated(tv%eqn_of_state)) then
254	0	pres(:) = tv%P_Ref
255		else
256	0	pres(:) = 0.0
257		endif
258	0	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	0	do j=js,je
265	0	do i=is,ie ; kb(i) = 1 ; enddo
266
267	0	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	0	do k=1,nz ; do i=is,ie
284	0	dtKd(i,k) = GV%Z_to_H * (dt * Kd_lay(i,j,k))
285		enddo ; enddo
286	0	do K=1,nz+1 ; do i=is,ie
287	0	dtKd_int(i,K) = GV%Z_to_H * (dt * Kd_int(i,j,K))
288		enddo ; enddo
289		endif
290
291	0	do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
292	0	do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; enddo
293	0	do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
294
295	0	if (GV%Boussinesq .or. GV%Semi_Boussinesq) then
296	0	do k=2,nz-1 ; do i=is,ie
297	0	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	0	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	0	call set_Ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)
310
311	0	do i=is,ie
312	0	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	0	do k=2,nz-1 ; do i=is,ie
319	0	dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
320	0	I2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))
321	0	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	0	if (CS%bulkmixedlayer) then
328	0	kb_min = nz
329	0	do i=is,ie
330	0	htot(i) = h(i,j,1) - Angstrom
331		enddo
332	0	do k=2,kmb ; do i=is,ie
333	0	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
334		enddo ; enddo
335	0	do i=is,ie
336	0	max_eakb(i) = MAX(Ent_bl(i,Kmb+1) + 0.5*htot(i), htot(i))
337	0	I_dSkbp1(i) = 1.0 / (Sref(i,kmb+2) - Sref(i,kmb+1))
338	0	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	0	is, ie, G, GV, CS, maxF_kb, eakb_maxF, do_i, F_kb_maxent)
347	0	do i=is,ie ; if (kb(i) <= nz) then
348	0	maxF(i,kb(i)) = MAX(0.0, maxF_kb(i), F_kb_maxent(i))
349	0	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	0	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	0	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	0	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	0	kmb, is, ie, G, GV, CS, F_kb_maxEnt, do_i_in=do_i)
363
364	0	do i=is,ie
365	0	do_entrain_eakb = .false.
366		! If error_max_eakb0 < 0, then buffer layers are always all entrained
367	0	if (do_i(i)) then ; if (err_max_eakb0(i) < 0.0) then
368	0	do_entrain_eakb = .true.
369		endif ; endif
370
371	0	if (do_entrain_eakb) then
372	0	eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
373		else
374	0	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	0	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	0	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	0	do k=kb_min,nz-1 ; do i=is,ie
411	0	if ((k == kb(i)+1) .and. CS%bulkmixedlayer) then
412	0	maxF(i,k) = ds_dsp1(i,k)*(F_kb_maxEnt(i) + htot(i))
413	0	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
414	0	elseif (k > kb(i)) then
415	0	maxF(i,k) = ds_dsp1(i,k)*(maxF(i,k-1) + htot(i))
416	0	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
417		endif
418		enddo ; enddo
419	0	do i=is,ie
420	0	maxF(i,nz) = 0.0
421	0	if (.not.CS%bulkmixedlayer) then
422	0	maxF_correct(i) = MAX(0.0, -(maxF(i,nz-1) + htot(i)))
423		endif
424	0	htot(i) = h(i,j,nz) - Angstrom
425		enddo
426	0	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	0	do k=nz-1,kb_min,-1 ; do i=is,ie ; if (do_i(i)) then
436	0	if (k >= kb(i)) then
437	0	maxF(i,k) = MIN(maxF(i,k),dsp1_ds(i,k+1)*maxF(i,k+1) + htot(i))
438	0	htot(i) = htot(i) + (h(i,j,k) - Angstrom)
439		endif
440	0	if (k == kb(i)) then
441	0	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	0	F_kb(i) = maxF(i,k)
447	0	if ((F_kb(i) <= maxF_kb(i)) .and. (eakb_maxF(i) <= max_eakb(i))) then
448	0	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	0	G, GV, CS, eakb, Angstrom)
454	0	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	0	error=err_eakb0)
458	0	if (eakb(i) < max_eakb(i)) then
459	0	max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)
460		endif
461	0	if (eakb(i) < min_eakb(i)) then
462	0	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	0	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	0	do i=is,ie
480	0	F(i,1) = maxF(i,1)
481	0	F(i,nz) = maxF(i,nz) ; minF(i,nz) = 0.0
482		enddo
483	0	do k=nz-1,K2,-1
484	0	do i=is,ie
485	0	if ((k==kb(i)) .and. (do_i(i))) then
486	0	eakb(i) = min_eakb(i)
487	0	minF(i,k) = 0.0
488	0	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	0	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	0	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	0	fk = dtKd(i,k) * grats(i,k)
503		minF(i,k) = MIN(maxF(i,k), &
504	0	0.9(I2p2dsp1_ds(i,k) fk / (hm + sqrt(hm*hm + fk))))
505	0	if (k==kb(i)) minF(i,k) = 0.0 ! BACKWARD COMPATIBILITY - DELETE LATER?
506		else
507	0	F(i,k) = 0.0
508	0	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	0	is1 = ie+1 ; ie1 = is-1
516	0	do i=is,ie ; if (do_i(i)) then ; is1 = i ; exit ; endif ; enddo
517	0	do i=ie,is,-1 ; if (do_i(i)) then ; ie1 = i ; exit ; endif ; enddo
518
519	0	if (CS%bulkmixedlayer) then
520	0	kb_min_act = nz
521	0	do i=is,ie
522	0	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	0	do i=is1,ie1
528	0	ea_kbp1(i) = 0.0
529	0	if (do_i(i) .and. (kb(i) < nz)) &
530	0	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	0	dFdfm_kb=dFdfm_kb)
536		else
537	0	kb_min_act = kb_min
538		endif
539
540	0	do it=0,CS%max_ent_it-1
541	0	do i=is1,ie1 ; if (do_i(i)) then
542	0	if (.not.CS%bulkmixedlayer) F(i,1) = MIN(F(i,1),maxF(i,1))
543	0	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	0	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	0	if (CS%bulkmixedlayer .and. (k==kb(i))) then
551	0	F(i,k) = F_kb(i)
552	0	dFdfm(i,k) = dFdfm_kb(i)
553		else ! k > kb(i)
554	0	Fprev(i,k) = F(i,k)
555	0	fm = (F(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*F(i,k+1)
556	0	fk = grats(i,k)*dtKd(i,k)
557	0	fr = sqrt(fm*fm + fk)
558
559	0	if (fm>=0) then
560	0	F(i,k) = MIN(maxF(i,k), I2p2dsp1_ds(i,k) * (fm+fr))
561		else
562	0	F(i,k) = MIN(maxF(i,k), I2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))
563		endif
564
565	0	if ((F(i,k) >= maxF(i,k)) .or. (fr == 0.0)) then
566	0	dFdfm(i,k) = 0.0
567		else
568	0	dFdfm(i,k) = I2p2dsp1_ds(i,k) * ((fr + fm) / fr)
569		endif
570
571	0	if (k > K2) then
572		! This is part of a tridiagonal solver for the actual flux.
573	0	c1(i,k) = dFdfm(i,k-1)(dsp1_ds(i,k)b1(i))
574	0	b1(i) = 1.0 / (1.0 - c1(i,k)*dFdfm(i,k))
575	0	F(i,k) = MIN(b1(i)*(F(i,k)-Fprev(i,k)) + Fprev(i,k), maxF(i,k))
576	0	if (F(i,k) >= maxF(i,k)) dFdfm(i,k) = 0.0
577		endif
578		endif
579		endif ; enddo ; enddo
580
581	0	do k=nz-2,kb_min_act,-1 ; do i=is1,ie1
582	0	if (do_i(i) .and. (k > kb(i))) &
583	0	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	0	if (CS%bulkmixedlayer) then
587	0	do i=is1,ie1
588	0	if (do_i(i) .and. (kb(i) < nz)) then
589		! F will be increased to minF later.
590	0	ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(F(i,kb(i)+1), minF(i,kb(i)+1))
591		else
592	0	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	0	dFdfm_kb=dFdfm_kb)
599	0	do i=is1,ie1
600	0	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	0	if (it < CS%max_ent_it-1) then
606
607	0	reiterate = .false.
608	0	if (CS%bulkmixedlayer) then ; do i=is1,ie1 ; if (do_i(i)) then
609	0	eb_kmb(i) = max(2.0*Ent_bl(i,Kmb+1) - eakb(i), 0.0)
610		endif ; enddo ; endif
611	0	do i=is1,ie1
612	0	did_i(i) = do_i(i) ; do_i(i) = .false.
613		enddo
614	0	do k=kb_min_act,nz-1 ; do i=is1,ie1
615	0	if (did_i(i) .and. (k >= kb(i))) then
616	0	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	0	elseif (k > kb(i)) then
620	0	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	0	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	0	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	0	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	0	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	0	if (associated(tv%eqn_of_state)) then
700	0	do i=is,ie
701	0	h_guess(i,1) = (h(i,j,1) - Angstrom) + (eb(i,j,1) - ea(i,j,2))
702	0	h_guess(i,nz) = (h(i,j,nz) - Angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))
703	0	if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0
704	0	if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0
705		enddo
706	0	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	0	(eb(i,j,k) - ea(i,j,k+1)))
709	0	if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0
710		enddo ; enddo
711	0	if (CS%bulkmixedlayer) then
712		call determine_dSkb(h_bl, Sref, Ent_bl, eakb, is, ie, kmb, G, GV, &
713	0	.true., dS_kb, dS_anom_lim=dS_anom_lim)
714	0	do k=nz-1,kb_min,-1
715	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, EOSdom)
716	0	do i=is,ie
717	0	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	0	(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	0	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	0	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	0	0.25ds_dsp1(i,k)h_guess(i,k-1) )
737		endif
738
739	0	ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*F_cor
740	0	eb(i,j,k) = eb(i,j,k) + F_cor
741	0	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	0	dS_kb_eff = 2.0*dS_kb(i) - dS_anom_lim(i) ! Could be negative!!!
746	0	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	0	if (abs(Rho_cor) < abs(0.9eakb(i)dS_kb_eff)) then
750	0	ea_cor = -Rho_cor / dS_kb_eff
751		else
752	0	ea_cor = sign(0.9eakb(i),-Rho_cordS_kb_eff)
753		endif
754
755	0	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	0	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	0	0.25h_guess(i,k+1) / (dS_kb(i) I_dSkbp1(i)))
763		endif
764
765	0	ea(i,j,k) = ea(i,j,k) + ea_cor
766	0	eb(i,j,k) = eb(i,j,k) - (dS_kb(i) * I_dSkbp1(i)) * ea_cor
767	0	elseif (k < kb(i)) then
768		! Repetitive, unless ea(kb) has been corrected.
769	0	ea(i,j,k) = ea(i,j,k+1)
770		endif
771		enddo
772		enddo
773	0	do k=kb_min-1,K2,-1 ; do i=is,ie
774	0	ea(i,j,k) = ea(i,j,k+1)
775		enddo ; enddo
776
777		! Repetitive, unless ea(kb) has been corrected.
778	0	k=kmb
779	0	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	0	h1 = (h(i,j,k) - Angstrom) + (eb(i,j,k) - ea(i,j,k+1))
783	0	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
784		enddo
785	0	do k=kmb-1,2,-1 ; do i=is,ie
786		! Determine the entrainment from below for each buffer layer.
787	0	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	0	h1 = (h(i,j,k) - Angstrom) + (eb(i,j,k) - ea(i,j,k+1))
791	0	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
792		enddo ; enddo
793	0	do i=is,ie
794	0	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	0	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	0	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	0	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	0	if (CS%id_Kd > 0) call post_data(CS%id_Kd, Kd_eff, CS%diag)
922	0	if (CS%id_Kd > 0) deallocate(Kd_eff)
923	0	if (CS%id_diff_work > 0) call post_data(CS%id_diff_work, diff_work, CS%diag)
924	0	if (CS%id_diff_work > 0) deallocate(diff_work)
925
926	0	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	0	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	0	is = G%isc ; ie = G%iec ; nz = GV%ke
966
967	0	do i=is,ie
968	0	ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0
969		enddo
970	0	if (CS%bulkmixedlayer) then
971	0	do i=is,ie
972	0	eb(i,j,kmb) = max(2.0*Ent_bl(i,Kmb+1) - eakb(i), 0.0)
973		enddo
974	0	do k=nz-1,kmb+1,-1 ; do i=is,ie
975	0	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	0	ea(i,j,k) = dsp1_ds(i,k)*F(i,k)
979	0	eb(i,j,k) = F(i,k)
980	0	elseif (k == kb(i)) then
981	0	ea(i,j,k) = eakb(i)
982	0	eb(i,j,k) = F(i,k)
983	0	elseif (k == kb(i)-1) then
984	0	ea(i,j,k) = ea(i,j,k+1)
985	0	eb(i,j,k) = eb(i,j,kmb)
986		else
987	0	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	0	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	0	k = kmb
994	0	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	0	if (kb(i) > kmb+1) &
998	0	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	0	h1 = (h(i,j,k) - GV%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
1002	0	ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1)
1003		enddo
1004	0	do k=kmb-1,2,-1 ; do i=is,ie
1005		! Determine the entrainment from below for each buffer layer.
1006	0	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	0	h1 = (h(i,j,k) - GV%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
1010	0	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	0	do i=is,ie
1016	0	eb(i,j,1) = max(2.0*Ent_bl(i,2) - ea(i,j,2), 0.0)
1017	0	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	0	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	0	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	0	b1, d1, & ! Variables used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] and [nondim].
1081	0	Rcv, & ! Value of the coordinate variable (potential density)
1082		! based on the simulated T and S and P_Ref [R ~> kg m-3].
1083	0	pres, & ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].
1084	0	frac_rem, & ! The fraction of the diffusion remaining [nondim].
1085	0	h_interior ! The interior thickness available for entrainment [H ~> m or kg m-2].
1086		real, dimension(SZI_(G), SZK_(GV)) :: &
1087	0	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	0	is = G%isc ; ie = G%iec ; nz = GV%ke
1097
1098	0	h_neglect = GV%H_subroundoff
1099
1100	0	do i=is,ie ; pres(i) = tv%P_Ref ; enddo
1101	0	EOSdom(:) = EOS_domain(G%HI)
1102	0	do k=1,kmb
1103	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, EOSdom)
1104	0	do i=is,ie
1105	0	h_bl(i,k) = h(i,j,k) + h_neglect
1106	0	Sref(i,k) = Rcv(i) - CS%Rho_sig_off
1107		enddo
1108		enddo
1109
1110	0	do i=is,ie
1111	0	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	0	do k=2,kmb ; do i=is,ie
1116	0	if (do_i(i)) then
1117	0	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	0	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	0	do i=is,ie
1127	0	b1(i) = 1.0 / (h_bl(i,1) + Ent_bl(i,2))
1128	0	d1(i) = h_bl(i,1) * b1(i) ! = 1.0 - Ent_bl(i,2)*b1(i)
1129	0	S_est(i,1) = (h_bl(i,1)Sref(i,1)) b1(i)
1130		enddo
1131	0	do k=2,kmb-1 ; do i=is,ie
1132	0	b1(i) = 1.0 / ((h_bl(i,k) + Ent_bl(i,K+1)) + d1(i)*Ent_bl(i,K))
1133	0	d1(i) = (h_bl(i,k) + d1(i)Ent_bl(i,K)) b1(i) ! = 1.0 - Ent_bl(i,K+1)*b1(i)
1134	0	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	0	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	0	(h_bl(i,kmb) + d1(i)*Ent_bl(i,Kmb))
1139	0	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	0	do i=is,ie ; kb(i) = nz+1 ; if (do_i(i)) kb(i) = kmb+1 ; enddo
1145
1146	0	do k=kmb+1,nz ; do i=is,ie ; if (do_i(i)) then
1147	0	if ((k == kb(i)) .and. (S_est(i,kmb) > (GV%Rlay(k) - CS%Rho_sig_off))) then
1148	0	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	0	dh = max((h(i,j,k) - GV%Angstrom_H), 0.0)
1152	0	if (dh > 0.0) then
1153		frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &
1154	0	(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	0	(h_bl(i,kmb) + dh)
1157	0	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	0	(h_bl(i,kmb) + d1(i)*Ent_bl(i,Kmb))
1160		endif
1161	0	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	0	do k=nz,kmb+1,-1 ; do i=is,ie
1169	0	if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-GV%Angstrom_H)
1170	0	if (k==kb(i)) then
1171	0	h_bl(i,kmb+1) = h(i,j,k) ; Sref(i,kmb+1) = GV%Rlay(k) - CS%Rho_sig_off
1172	0	elseif (k==kb(i)+1) then
1173	0	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	0	do i=is,ie ; if (kb(i) >= nz) then
1177	0	h_bl(i,kmb+1) = h(i,j,nz)
1178	0	Sref(i,kmb+1) = GV%Rlay(nz) - CS%Rho_sig_off
1179	0	h_bl(i,kmb+2) = GV%Angstrom_H
1180	0	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	0	do i=is,ie ; if (do_i(i)) then
1187	0	Kd_x_dt = frac_rem(i) * dtKd_int(i,Kmb+1)
1188	0	if ((Kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0)) then
1189	0	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	0	Kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))
1195		endif
1196		else
1197	0	Ent_bl(i,Kmb+1) = 0.0
1198		endif ; enddo
1199
1200	0	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	0	subroutine determine_dSkb(h_bl, Sref, Ent_bl, E_kb, is, ie, kmb, G, GV, limit, &
1215	0	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	0	b1, c1, & ! b1 [H-1 ~> m-1 or m2 kg-1] and c1 [nondim] are variables used by the tridiagonal solver.
1268	0	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	0	ea, dea_dE, & ! The entrainment from above [H ~> m or kg m-2] and its derivative with E [nondim].
1270	0	eb, deb_dE ! The entrainment from below [H ~> m or kg m-2] and its derivative with E [nondim].
1271	0	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	0	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	0	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	0	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	0	h_neglect = GV%H_subroundoff
1298
1299	0	do i=is,ie
1300	0	ea(i,kmb+1) = E_kb(i) ; dea_dE(i,kmb+1) = 1.0
1301	0	S(i,kmb+1) = Sref(i,kmb+1) ; dS_dE(i,kmb+1) = 0.0
1302	0	b1(i,kmb+1) = 0.0
1303	0	d1(i) = 1.0
1304	0	do_i(i) = .true.
1305		enddo
1306	0	if (present(do_i_in)) then
1307	0	do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1308		endif
1309	0	do k=kmb,1,-1 ; do i=is,ie
1310	0	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	0	if (2.0*Ent_bl(i,K+1) > ea(i,k+1)) then
1314	0	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	0	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	0	h1 = (h_bl(i,k) - GV%Angstrom_H) + (eb(i,k) - ea(i,k+1))
1321	0	if (h1 >= 0.0) then
1322	0	ea(i,k) = Ent_bl(i,K) ; dea_dE(i,k) = 0.0
1323	0	elseif (Ent_bl(i,K) + 0.5*h1 >= 0.0) then
1324	0	ea(i,k) = Ent_bl(i,K) - 0.5*h1
1325	0	dea_dE(i,k) = 0.5*(dea_dE(i,k+1) - deb_dE(i,k))
1326		else
1327	0	ea(i,k) = -h1
1328	0	dea_dE(i,k) = dea_dE(i,k+1) - deb_dE(i,k)
1329		endif
1330		else
1331	0	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	0	h_tr = h_bl(i,k) + h_neglect
1336	0	c1(i,k) = ea(i,k+1) * b1(i,k+1)
1337	0	b_denom_1 = (h_tr + d1(i)*eb(i,k))
1338	0	b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))
1339	0	d1(i) = b_denom_1 * b1(i,k)
1340
1341	0	S(i,k) = (h_trSref(i,k) + eb(i,k)S(i,k+1)) * b1(i,k)
1342		enddo ; enddo
1343	0	do k=2,kmb ; do i=is,ie
1344	0	S(i,k) = S(i,k) + c1(i,k-1)*S(i,k-1)
1345		enddo ; enddo
1346
1347	0	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	0	do k=kmb,2,-1 ; do i=is,ie
1351	0	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	0	((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))
1357	0	else ; src = 0.0 ; endif
1358	0	dS_dE(i,k) = (src + eb(i,k)dS_dE(i,k+1)) b1(i,k)
1359		enddo ; enddo
1360	0	do i=is,ie
1361	0	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	0	(h_bl(i,1) + h_neglect + eb(i,1))
1364	0	else ; src = 0.0 ; endif
1365	0	dS_dE(i,1) = (src + eb(i,1)dS_dE(i,2)) b1(i,1)
1366		enddo
1367	0	do k=2,kmb ; do i=is,ie
1368	0	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	0	eps_dSkb = 1.0e-6 ! Should be a small, nondimensional, positive number.
1375	0	if (.not.limit) then
1376	0	do i=is,ie ; if (do_i(i)) then
1377	0	dSkb(i) = Sref(i,kmb+1) - S(i,kmb)
1378		endif ; enddo
1379	0	if (present(ddSkb_dE)) then ; do i=is,ie ; if (do_i(i)) then
1380	0	ddSkb_dE(i) = -1.0*dS_dE(i,kmb)
1381		endif ; enddo ; endif
1382
1383	0	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	0	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	0	do i=is,ie ; if (do_i(i)) then
1391		! Need to ensure that 0 < dSkb <= S_kb - Sbl
1392	0	if (Sref(i,kmb+1) - S(i,kmb) < eps_dSkb*(Sref(i,kmb+2) - Sref(i,kmb+1))) then
1393	0	dSkb(i) = eps_dSkb * (Sref(i,kmb+2) - Sref(i,kmb+1)) ; deriv_dSkb(i) = 0.0
1394		else
1395	0	dSkb(i) = Sref(i,kmb+1) - S(i,kmb) ; deriv_dSkb(i) = -1.0
1396		endif
1397	0	if (present(ddSkb_dE)) ddSkb_dE(i) = deriv_dSkb(i)*dS_dE(i,kmb)
1398		endif ; enddo
1399
1400	0	if (present(dSLay)) then
1401	0	dz_drat = 1000.0 ! The limit of large dz_drat the same as choosing a
1402		! Heaviside function.
1403	0	eps_dSLay = 1.0e-10 ! Should be ~= GV%Angstrom_H / sqrt(Kd*dt)
1404	0	do i=is,ie ; if (do_i(i)) then
1405	0	dS_kbp1 = Sref(i,kmb+2) - Sref(i,kmb+1)
1406	0	IdS_kbp1 = 1.0 / (Sref(i,kmb+2) - Sref(i,kmb+1))
1407	0	rat = (Sref(i,kmb+1) - S(i,kmb)) * IdS_kbp1
1408		! Need to ensure that 0 < dSLay <= 2*dSkb
1409	0	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	0	Inv_term = 1.0 / (1.0-rat)
1416	0	f1 = 2.0 - 0.125(Inv_term*2)
1417	0	df1_drat = - 0.25(Inv_term*3)
1418
1419		! f2 ensures that dSLay goes to 0 rapidly if rat is significantly
1420		! negative.
1421	0	z = dz_drat * rat + 4.0 ! The 4 here gives f2(0) = 0.982.
1422	0	if (z >= 18.0) then ; f2 = 1.0 ; df2_dz = 0.0
1423	0	elseif (z <= -58.0) then ; f2 = eps_dSLay ; df2_dz = 0.0
1424		else
1425	0	expz = exp(z) ; Inv_term = 1.0 / (1.0 + expz)
1426	0	f2 = (eps_dSLay + expz) * Inv_term
1427	0	df2_dz = (1.0 - eps_dSLay) * expz * Inv_term**2
1428		endif
1429
1430	0	dSLay(i) = dSkb(i) * f1 * f2
1431		deriv_dSLay = deriv_dSkb(i) * (f1 * f2) - (dSkb(i)IdS_kbp1) &
1432	0	(df1_dratf2 + f1 dz_drat * df2_dz)
1433	0	elseif (dSkb(i) <= 3.0*dS_kbp1) then
1434		! This is the "typical" case.
1435	0	dSLay(i) = 0.5 * (dSkb(i) + dS_kbp1)
1436	0	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	0	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	0	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	0	(Sref(i,kmb+1) - S(i,kmb)) )
1449		endif ; enddo ; endif
1450
1451	0	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	0	subroutine F_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &
1458	0	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	0	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	0	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	0	real :: ent_maxF(SZI_(G)) ! The value of entrainment that gives maxF [H ~> m or kg m-2]
1497	0	real :: zeros(SZI_(G)) ! An array of zero entrainments [H ~> m or kg m-2]
1498	0	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	0	dS_kbp1 = Sref(i,kmb+2) - Sref(i,kmb+1)
1504	0	max_ea = ea_kb(i) ; min_ea = 0.0
1505	0	val = dS_kbp1 * F_kb(i)
1506	0	err_min = -val
1507
1508	0	tolerance = CS%Tolerance_Ent
1509	0	if (present(tol_in)) tolerance = tol_in
1510	0	bisect_next = .true.
1511
1512		call determine_dSkb(h_bl, Sref, Ent_bl, ea_kb, i, i, kmb, G, GV, .true., &
1513	0	dS_kb, ddSkb_dE)
1514
1515	0	err = dS_kb(i) * ea_kb(i) - val
1516	0	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	0	if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea))) return
1520
1521	0	if (err == 0.0) then ; return ! The exact solution on the first guess...
1522	0	elseif (err > 0.0) then ! The root is properly bracketed.
1523	0	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	0	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	0	ea_kb(i) = ea_kb(i) - err / derr_dea
1529		else ! Use the bisection for the next guess.
1530	0	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	0	do it = 1, MAXIT
1553		call determine_dSkb(h_bl, Sref, Ent_bl, ea_kb, i, i, kmb, G, GV, .true., &
1554	0	dS_kb, ddSkb_dE)
1555
1556	0	err = dS_kb(i) * ea_kb(i) - val
1557	0	derr_dea = dS_kb(i) + ddSkb_dE(i) * ea_kb(i)
1558
1559	0	ea_prev = ea_kb(i)
1560		! Use Newton's method or the false position method to find the next value.
1561	0	Newton = .false.
1562	0	if (err > 0.0) then
1563	0	max_ea = ea_kb(i) ; err_max = err
1564	0	if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) Newton = .true.
1565		else
1566	0	min_ea = ea_kb(i) ; err_min = err
1567	0	if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) Newton = .true.
1568		endif
1569
1570	0	if (Newton) then
1571	0	ea_kb(i) = ea_kb(i) - err / derr_dea
1572	0	elseif (bisect_next) then ! Use bisection to reduce the range.
1573	0	ea_kb(i) = 0.5*(max_ea+min_ea)
1574	0	bisect_next = .false.
1575		else ! Use the false-position method for the next guess.
1576	0	ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))
1577	0	bisect_next = .true.
1578		endif
1579
1580	0	tol1 = tolerance ; if (err > 0.0) tol1 = 0.099*tolerance
1581	0	if (dS_kb(i) <= dS_kbp1) then
1582	0	if (abs(ea_kb(i) - ea_prev) <= tol1) return
1583		else
1584	0	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	0	subroutine determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &
1595	0	min_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, Ent, &
1596	0	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	0	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	0	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	0	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	0	derror_dE, & ! The derivative of err with E [H ~> m or kg m-2].
1663	0	err, & ! The "error" whose zero is being sought [H2 ~> m2 or kg2 m-4].
1664	0	E_min, E_max, & ! The minimum and maximum values of E [H ~> m or kg m-2].
1665	0	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	0	logical, dimension(SZI_(G)) :: false_position ! If true, the false position
1676		! method might be used for the next iteration.
1677	0	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	0	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	0	tolerance = CS%Tolerance_Ent
1688	0	large_err = GV%m_to_H*2 1.0e30
1689
1690	0	do i=is,ie ; redo_i(i) = do_i(i) ; enddo
1691
1692	0	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	0	E_min(i) = min_eakb(i) ; E_max(i) = max_eakb(i)
1698	0	error_minE(i) = -large_err ; error_maxE(i) = large_err
1699	0	false_position(i) = .true. ! Used to alternate between false_position and
1700		! bisection when Newton's method isn't working.
1701	0	if (present(err_min_eakb0)) error_minE(i) = err_min_eakb0(i) - E_min(i) * ea_kbp1(i)
1702	0	if (present(err_max_eakb0)) error_maxE(i) = err_max_eakb0(i) - E_max(i) * ea_kbp1(i)
1703
1704	0	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	0	if (error_maxE(i) <= 0.0) then
1707		! The errors decrease with E*ea_kbp1, so E_max is the best solution.
1708	0	Ent(i) = E_max(i) ; err(i) = error_maxE(i)
1709		else ! error_minE >= 0 is equivalent to ea_kbp1 = 0.0.
1710	0	Ent(i) = E_min(i) ; err(i) = error_minE(i)
1711		endif
1712	0	derror_dE(i) = 0.0
1713	0	redo_i(i) = .false.
1714		endif
1715		endif ; enddo ! End of i-loop
1716
1717	0	do it = 1,MAXIT
1718	0	do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo
1719	0	if (.not.do_any) exit
1720		call determine_dSkb(h_bl, Sref, Ent_bl, Ent, is, ie, kmb, G, GV, .true., dS_kb, &
1721	0	ddSkb_dE, dS_lay, ddSlay_dE, do_i_in=redo_i)
1722	0	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	0	eL = 0.0 ; if (2.0*Ent_bl(i,Kmb+1) >= Ent(i)) eL = 1.0
1726	0	fa = (1.0 + eL) + dS_kb(i)*I_dSkbp1(i)
1727	0	fk = dtKd_kb(i) * (dS_Lay(i)/dS_kb(i))
1728	0	fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + eL2.0Ent_bl(i,Kmb+1)
1729	0	if (fm > -GV%Angstrom_H) fm = fm + GV%Angstrom_H ! This could be smooth if need be.
1730	0	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	0	dtKd_kb(i) * (ddSlay_dE(i)dS_kb(i) - ddSkb_dE(i)dS_Lay(i))/(dS_kb(i)**2)
1733
1734	0	if (err(i) == 0.0) then
1735	0	redo_i(i) = .false. ; cycle
1736	0	elseif (err(i) > 0.0) then
1737	0	E_max(i) = Ent(i) ; error_maxE(i) = err(i)
1738		else
1739	0	E_min(i) = Ent(i) ; error_minE(i) = err(i)
1740		endif
1741
1742	0	E_prev = Ent(i)
1743	0	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	0	fr = sqrt(fm*2 + 4.0fa*fk)
1749	0	if (fm >= 0.0) then
1750	0	Ent(i) = (fm + fr) / (2.0 * fa)
1751		else
1752	0	Ent(i) = (2.0 * fk) / (fr - fm)
1753		endif
1754		! But make sure that the root stays bracketed, bisecting if needed.
1755	0	if ((Ent(i) > E_max(i)) .or. (Ent(i) < E_min(i))) &
1756	0	Ent(i) = 0.5*(E_max(i) + E_min(i))
1757	0	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	0	Ent(i) = Ent(i) - err(i) / derror_dE(i)
1762	0	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	0	Ent(i) = 0.5*(E_max(i) + E_min(i))
1770	0	false_position(i) = .true.
1771		endif
1772
1773	0	if (abs(E_prev - Ent(i)) < tolerance) then
1774	0	err_est = err(i) + (Ent(i) - E_prev) * derror_dE(i)
1775	0	if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &
1776	0	(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	0	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	0	dS_kb, do_i_in=do_i)
1786
1787	0	if (present(F_kb)) then ; do i=is,ie ; if (do_i(i)) then
1788	0	F_kb(i) = Ent(i) * (dS_kb(i) * I_dSkbp1(i))
1789		endif ; enddo ; endif
1790	0	if (present(error)) then ; do i=is,ie ; if (do_i(i)) then
1791	0	error(i) = err(i)
1792		endif ; enddo ; endif
1793	0	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	0	if (derror_dE(i) > 0.0) then
1798		dFdfm_kb(i) = ((dS_kb(i) + Ent(i) * ddSkb_dE(i)) * I_dSkbp1(i)) * &
1799	0	(Ent(i) / derror_dE(i))
1800		else ! Use Adcroft's division by 0 convention.
1801	0	dFdfm_kb(i) = 0.0
1802		endif
1803		endif ; enddo ; endif
1804
1805	0	end subroutine determine_Ea_kb
1806
1807		!> Maximize F = entds_kbI_dSkbp1 in the range min_ent < ent < max_ent.
1808	0	subroutine find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, &
1809	0	kmb, is, ie, G, GV, CS, maxF, ent_maxF, do_i_in, &
1810	0	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	0	ent, & ! The updated estimate of the entrainment [H ~> m or kg m-2]
1859	0	minent, maxent, ent_best, & ! Various previous estimates of the entrainment [H ~> m or kg m-2]
1860	0	F_max_ent_in, & ! The value of F that gives the input maximum value of ent [H ~> m or kg m-2]
1861	0	F_maxent, F_minent, F, F_best, & ! Various estimates of F [H ~> m or kg m-2]
1862	0	dF_dent, dF_dE_max, dF_dE_min, dF_dE_best, & ! Various derivatives of F with ent [nondim]
1863	0	dS_kb, & ! The density difference across the interface between the bottommost
1864		! buffer layer and the topmost interior layer [R ~> kg m-3]
1865	0	dS_kb_lim, dS_anom_lim, & ! Various limits on dS_kb [R ~> kg m-3]
1866	0	ddSkb_dE, & ! The partial derivative of dS_kb with ent [R H-1 ~> kg m-4 or m-1].
1867	0	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	0	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	0	tolerance = CS%Tolerance_Ent
1879
1880	0	if (present(do_i_in)) then
1881	0	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	0	dS_kb, ddSkb_dE, dS_anom_lim=dS_anom_lim)
1889	0	ie1 = is-1 ; doany = .false.
1890	0	do i=is,ie
1891	0	dS_kb_lim(i) = dS_kb(i) + dS_anom_lim(i)
1892	0	F_max_ent_in(i) = max_ent_in(i)dS_kb_lim(i)I_dSkbp1(i)
1893	0	maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)
1894	0	if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i))) then
1895	0	F_best(i) = max_ent_in(i)dS_kb(i)I_dSkbp1(i)
1896	0	ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)
1897	0	do_i(i) = .false.
1898		else
1899	0	F_maxent(i) = maxent(i) * dS_kb(i) * I_dSkbp1(i)
1900	0	dF_dE_max(i) = (dS_kb(i) + maxent(i)ddSkb_dE(i)) I_dSkbp1(i)
1901	0	doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.
1902		endif
1903		enddo
1904
1905	0	if (doany) then
1906	0	ie1 = is-1 ; do i=is,ie ; if (do_i(i)) ie1 = i ; enddo
1907	0	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	0	dS_kb, ddSkb_dE, do_i_in=do_i)
1911	0	do i=is1,ie1 ; if (do_i(i)) then
1912	0	F_minent(i) = minent(i) * dS_kb(i) * I_dSkbp1(i)
1913	0	dF_dE_min(i) = (dS_kb(i) + minent(i)ddSkb_dE(i)) I_dSkbp1(i)
1914		endif ; enddo
1915
1916	0	ratio_select_end = 0.9
1917	0	do it=1,MAXIT
1918	0	ratio_select_end = 0.5*ratio_select_end
1919	0	do i=is1,ie1 ; if (do_i(i)) then
1920	0	if (need_bracket(i)) then
1921	0	dF_dE_mean = (F_maxent(i) - F_minent(i)) / (maxent(i) - minent(i))
1922	0	maxslope = MAX(dF_dE_mean, dF_dE_min(i), dF_dE_max(i))
1923	0	minslope = MIN(dF_dE_mean, dF_dE_min(i), dF_dE_max(i))
1924	0	if (F_minent(i) >= F_maxent(i)) then
1925	0	if (dF_dE_min(i) > 0.0) then ; rat = 0.02 ! A small step should bracket the solution.
1926	0	elseif (maxslope < ratio_select_end*minslope) then
1927		! The maximum of F is at minent.
1928	0	F_best(i) = F_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0
1929	0	do_i(i) = .false.
1930	0	else ; rat = 0.382 ; endif ! Use the golden ratio
1931		else
1932	0	if (dF_dE_max(i) < 0.0) then ; rat = 0.98 ! A small step should bracket the solution.
1933	0	elseif (minslope > ratio_select_end*maxslope) then
1934		! The maximum of F is at maxent.
1935	0	F_best(i) = F_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0
1936	0	do_i(i) = .false.
1937	0	else ; rat = 0.618 ; endif ! Use the golden ratio
1938		endif
1939
1940	0	if (rat >= 0.0) ent(i) = ratmaxent(i) + (1.0-rat)minent(i)
1941	0	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	0	if (mod(it,3) == 0) then ! Re-determine the loop bounds.
1985	0	ie1 = is-1 ; do i=is1,ie ; if (do_i(i)) ie1 = i ; enddo
1986	0	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	0	dS_kb, ddSkb_dE, do_i_in=do_i)
1991	0	do i=is1,ie1 ; if (do_i(i)) then
1992	0	F(i) = ent(i)dS_kb(i)I_dSkbp1(i)
1993	0	dF_dent(i) = (dS_kb(i) + ent(i)ddSkb_dE(i)) I_dSkbp1(i)
1994		endif ; enddo
1995
1996	0	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	0	doany = .false.
2003	0	do i=is1,ie1 ; if (do_i(i)) then
2004	0	if (.not.last_it(i)) doany = .true.
2005	0	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	0	elseif (need_bracket(i)) then
2019	0	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	0	elseif ((F(i) <= F_maxent(i)) .and. (F(i) > F_minent(i))) then
2025	0	new_min_bound = .true. ! We have a new minimum bound.
2026	0	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	0	new_min_bound = .true.
2034	0	if (dF_dE_min(i) > 0.0 .or. (F_minent(i) >= F_maxent(i))) &
2035	0	new_min_bound = .false.
2036		endif
2037	0	if (need_bracket(i)) then ! Still not bracketed.
2038	0	if (new_min_bound) then
2039	0	minent(i) = ent(i) ; F_minent(i) = F(i) ; dF_dE_min(i) = dF_dent(i)
2040		else
2041	0	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	0	if (.not.doany) exit
2063		enddo
2064		endif
2065
2066	0	if (present(F_lim_maxent)) then
2067		! Return the unlimited maximum in maxF, and the limited value of F at maxent.
2068	0	do i=is,ie
2069	0	maxF(i) = F_best(i)
2070	0	F_lim_maxent(i) = F_max_ent_in(i)
2071	0	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	0	doany = .false.
2076	0	do i=is,ie
2077	0	may_use_best(i) = (ent_best(i) /= max_ent_in(i))
2078	0	if (may_use_best(i)) doany = .true.
2079		enddo
2080	0	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	0	do i=is,ie ; maxF(i) = F_max_ent_in(i) ; enddo
2098	0	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	0	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	0	'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	0	'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