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