MOM_entrain_diffusive.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.

2! See the LICENSE file for licensing information.

3! SPDX-License-Identifier: Apache-2.0

4

5!> Diapycnal mixing and advection in isopycnal mode

6module mom_entrain_diffusive

7

8use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr

9use mom_diag_mediator, only : diag_ctrl, time_type

10use mom_eos, only : calculate_density, calculate_density_derivs

11use mom_eos, only : calculate_specific_vol_derivs, eos_domain

12use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note

13use mom_file_parser, only : get_param, log_version, param_file_type

14use mom_forcing_type, only : forcing

15use mom_grid, only : ocean_grid_type

16use mom_unit_scaling, only : unit_scale_type

17use mom_variables, only : thermo_var_ptrs

18use mom_verticalgrid, only : verticalgrid_type

19

20implicit none ; private

21

22#include <MOM_memory.h>

23

24public entrainment_diffusive, entrain_diffusive_init

25

26! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional

27! consistency testing. These are noted in comments with units like Z, H, L, and T, along with

28! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units

29! vary with the Boussinesq approximation, the Boussinesq variant is given first.

30

31!> The control structure holding parametes for the MOM_entrain_diffusive module

32type, public :: entrain_diffusive_cs ; private

33 logical :: initialized = .false. !< True if this control structure has been initialized.

34 logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with

35 !! GV%nk_rho_varies variable density mixed & buffer layers.

36 integer :: max_ent_it !< The maximum number of iterations that may be used to

37 !! calculate the diapycnal entrainment.

38 real :: tolerance_ent !< The tolerance with which to solve for entrainment values

39 !! [H ~> m or kg m-2].

40 real :: max_ent !< A large ceiling on the maximum permitted amount of entrainment

41 !! across each interface between the mixed and buffer layers within

42 !! a timestep [H ~> m or kg m-2].

43 real :: rho_sig_off !< The offset between potential density and a sigma value [R ~> kg m-3]

44 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to

45 !! regulate the timing of diagnostic output.

46 integer :: id_kd = -1 !< Diagnostic ID for diffusivity

47 integer :: id_diff_work = -1 !< Diagnostic ID for mixing work

48end type entrain_diffusive_cs

49

50contains

51

52!> This subroutine calculates ea and eb, the rates at which a layer entrains

53!! from the layers above and below. The entrainment rates are proportional to

54!! the buoyancy flux in a layer and inversely proportional to the density

55!! differences between layers. The scheme that is used here is described in

56!! detail in Hallberg, Mon. Wea. Rev. 2000.

57subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &

58 kb_out, Kd_Lay, Kd_int)

59 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

60 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.

61 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

62 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

63 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].

64 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available

65 !! thermodynamic fields. Absent fields have NULL

66 !! ptrs.

67 type(forcing), intent(in) :: fluxes !< A structure of surface fluxes that may

68 !! be used.

69 real, intent(in) :: dt !< The time increment [T ~> s].

70 type(entrain_diffusive_cs), intent(in) :: cs !< The control structure returned by a previous

71 !! call to entrain_diffusive_init.

72 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

73 intent(out) :: ea !< The amount of fluid entrained from the layer

74 !! above within this time step [H ~> m or kg m-2].

75 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

76 intent(out) :: eb !< The amount of fluid entrained from the layer

77 !! below within this time step [H ~> m or kg m-2].

78 integer, dimension(SZI_(G),SZJ_(G)), &

79 intent(inout) :: kb_out !< The index of the lightest layer denser than

80 !! the buffer layer.

81 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

82 intent(in) :: kd_lay !< The diapycnal diffusivity of layers

83 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

84 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &

85 intent(in) :: kd_int !< The diapycnal diffusivity of interfaces

86 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

87

88! This subroutine calculates ea and eb, the rates at which a layer entrains

89! from the layers above and below. The entrainment rates are proportional to

90! the buoyancy flux in a layer and inversely proportional to the density

91! differences between layers. The scheme that is used here is described in

92! detail in Hallberg, Mon. Wea. Rev. 2000.

93

94 real, dimension(SZI_(G),SZK_(GV)) :: &

95 dtkd ! The layer diapycnal diffusivity times the time step [H2 ~> m2 or kg2 m-4].

96 real, dimension(SZI_(G),SZK_(GV)+1) :: &

97 dtkd_int ! The diapycnal diffusivity at the interfaces times the time step [H2 ~> m2 or kg2 m-4]

98 real, dimension(SZI_(G),SZK_(GV)) :: &

99 f, & ! The density flux through a layer within a time step divided by the

100 ! density difference across the interface below the layer [H ~> m or kg m-2].

101 maxf, & ! maxF is the maximum value of F that will not deplete all of the

102 ! layers above or below a layer within a timestep [H ~> m or kg m-2].

103 minf, & ! minF is the minimum flux that should be expected in the absence of

104 ! interactions between layers [H ~> m or kg m-2].

105 fprev, &! The previous estimate of F [H ~> m or kg m-2].

106 dfdfm, &! The partial derivative of F with respect to changes in F of the

107 ! neighboring layers. [nondim]

108 h_guess ! An estimate of the layer thicknesses after entrainment, but

109 ! before the entrainments are adjusted to drive the layer

110 ! densities toward their target values [H ~> m or kg m-2].

111 real, dimension(SZI_(G),SZK_(GV)+1) :: &

112 ent_bl ! The average entrainment upward and downward across

113 ! each interface around the buffer layers [H ~> m or kg m-2].

114 real, allocatable, dimension(:,:,:) :: &

115 kd_eff, & ! The effective diffusivity that actually applies to each

116 ! layer after the effects of boundary conditions are

117 ! considered [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

118 diff_work ! The work actually done by diffusion across each

119 ! interface [R Z3 T-3 ~> W m-2]. Sum vertically for the total work.

120

121 real :: hm, fm, fr ! Work variables with units of [H ~> m or kg m-2].

122 real :: fk ! A Work variable with units of [H2 ~> m2 or kg2 m-4]

123

124 real :: b1(szi_(g)) ! A variable used by the tridiagonal solver [H ~> m or kg m-2]

125 real :: c1(szi_(g),szk_(gv)) ! A variable used by the tridiagonal solver [nondim]

126

127 real, dimension(SZI_(G)) :: &

128 htot, & ! The total thickness above or below a layer [H ~> m or kg m-2].

129 rcv, & ! Value of the coordinate variable (potential density)

130 ! based on the simulated T and S and P_Ref [R ~> kg m-3].

131 pres, & ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].

132 eakb, & ! The entrainment from above by the layer below the buffer

133 ! layer (i.e. layer kb) [H ~> m or kg m-2].

134 ea_kbp1, & ! The entrainment from above by layer kb+1 [H ~> m or kg m-2].

135 eb_kmb, & ! The entrainment from below by the deepest buffer layer [H ~> m or kg m-2].

136 ds_kb, & ! The reference potential density difference across the

137 ! interface between the buffer layers and layer kb [R ~> kg m-3].

138 ds_anom_lim, &! The amount by which dS_kb is reduced when limits are

139 ! applied [R ~> kg m-3].

140 i_dskbp1, & ! The inverse of the potential density difference across the

141 ! interface below layer kb [R-1 ~> m3 kg-1].

142 dtkd_kb, & ! The diapycnal diffusivity in layer kb times the time step

143 ! [H2 ~> m2 or kg2 m-4].

144 maxf_correct, & ! An amount by which to correct maxF due to excessive

145 ! surface heat loss [H ~> m or kg m-2].

146 zeros, & ! An array of all zeros. (Usually used with [H ~> m or kg m-2].)

147 max_eakb, & ! The maximum value of eakb that might be realized [H ~> m or kg m-2].

148 min_eakb, & ! The minimum value of eakb that might be realized [H ~> m or kg m-2].

149 err_max_eakb0, & ! The value of error returned by determine_Ea_kb when eakb = max_eakb

150 ! and ea_kbp1 = 0 [H2 ~> m2 or kg2 m-4].

151 err_min_eakb0, & ! The value of error returned by determine_Ea_kb when eakb = min_eakb

152 ! and ea_kbp1 = 0 [H2 ~> m2 or kg2 m-4].

153 err_eakb0, & ! A value of error returned by determine_Ea_kb [H2 ~> m2 or kg2 m-4].

154 f_kb, & ! The value of F in layer kb, or equivalently the entrainment

155 ! from below by layer kb [H ~> m or kg m-2].

156 dfdfm_kb, & ! The partial derivative of F with fm [nondim]. See dFdfm.

157 maxf_kb, & ! The maximum value of F_kb that might be realized [H ~> m or kg m-2].

158 eakb_maxf, & ! The value of eakb that gives F_kb=maxF_kb [H ~> m or kg m-2].

159 f_kb_maxent ! The value of F_kb when eakb = max_eakb [H ~> m or kg m-2].

160 real, dimension(SZI_(G),SZK_(GV)) :: &

161 sref, & ! The reference potential density of the mixed and buffer layers,

162 ! and of the two lightest interior layers (kb and kb+1) copied

163 ! into layers kmb+1 and kmb+2 [R ~> kg m-3].

164 h_bl ! The thicknesses of the mixed and buffer layers, and of the two

165 ! lightest interior layers (kb and kb+1) copied into layers kmb+1

166 ! and kmb+2 [H ~> m or kg m-2].

167

168 real, dimension(SZI_(G),SZK_(GV)) :: &

169 ds_dsp1, & ! The coordinate variable (sigma-2) difference across an

170 ! interface divided by the difference across the interface

171 ! below it. [nondim]

172 dsp1_ds, & ! The inverse coordinate variable (sigma-2) difference

173 ! across an interface times the difference across the

174 ! interface above it. [nondim]

175 i2p2dsp1_ds, & ! 1 / (2 + 2 * ds_k+1 / ds_k). [nondim]

176 grats ! 2*(2 + ds_k+1 / ds_k + ds_k / ds_k+1) =

177 ! 4*ds_Lay*(1/ds_k + 1/ds_k+1). [nondim]

178

179 real :: drho ! The change in locally referenced potential density between

180 ! the layers above and below an interface [R ~> kg m-3]

181 real :: dspv ! The change in locally referenced specific volume between

182 ! the layers above and below an interface [R-1 ~> m3 kg-1]

183 real :: g_2dt ! 0.5 * G_Earth / dt, times unit conversion factors

184 ! [Z3 H-2 T-3 or R2 Z3 H-2 T-3 ~> m s-3].

185 real, dimension(SZI_(G)) :: &

186 pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].

187 t_eos, s_eos, & ! The potential temperature and salinity at which to

188 ! evaluate dRho_dT and dRho_dS [C ~> degC] and [S ~> ppt].

189 drho_dt, & ! The partial derivative of potential density with temperature [R C-1 ~> kg m-3 degC-1]

190 drho_ds, & ! The partial derivative of potential density with salinity [R S-1 ~> kg m-3 ppt-1]

191 dspv_dt, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]

192 dspv_ds ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]

193

194 real :: tolerance ! The tolerance within which E must be converged [H ~> m or kg m-2].

195 real :: angstrom ! The minimum layer thickness [H ~> m or kg m-2].

196 real :: h_neglect ! A thickness that is so small it is usually lost

197 ! in roundoff and can be neglected [H ~> m or kg m-2].

198 real :: f_cor ! A correction to the amount of F that is used to

199 ! entrain from the layer above [H ~> m or kg m-2].

200 real :: kd_here ! The effective diapycnal diffusivity times the timestep [H2 ~> m2 or kg2 m-4].

201 real :: h_avail ! The thickness that is available for entrainment [H ~> m or kg m-2].

202 real :: ds_kb_eff ! The value of dS_kb after limiting is taken into account [R ~> kg m-3].

203 real :: rho_cor ! The depth-integrated potential density anomaly that

204 ! needs to be corrected for [H R ~> kg m-2 or kg2 m-5].

205 real :: ea_cor ! The corrective adjustment to eakb [H ~> m or kg m-2].

206 real :: h1 ! The layer thickness after entrainment through the

207 ! interface below is taken into account [H ~> m or kg m-2].

208 real :: idt ! The inverse of the time step [Z H-1 T-1 ~> s-1 or m3 kg-1 s-1].

209

210 logical :: do_any

211 logical :: do_entrain_eakb ! True if buffer layer is entrained

212 logical :: do_i(szi_(g)), did_i(szi_(g)), reiterate

213 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state

214 integer :: it, i, j, k, is, ie, js, je, nz, k2, kmb

215 integer :: kb(szi_(g)) ! The value of kb in row j.

216 integer :: kb_min ! The minimum value of kb in the current j-row.

217 integer :: kb_min_act ! The minimum active value of kb in the current j-row.

218 integer :: is1, ie1 ! The minimum and maximum active values of i in the current j-row.

219 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

220 angstrom = gv%Angstrom_H

221 h_neglect = gv%H_subroundoff

222

223 if (.not. cs%initialized) call mom_error(fatal, &

224 "MOM_entrain_diffusive: Module must be initialized before it is used.")

225

226 if ((.not.cs%bulkmixedlayer .and. .not.associated(fluxes%buoy)) .and. &

227 (associated(fluxes%lprec) .or. associated(fluxes%evap) .or. &

228 associated(fluxes%sens) .or. associated(fluxes%sw))) then

229 if (is_root_pe()) call mom_error(note, "Calculate_Entrainment: &

230 &The code to handle evaporation and precipitation without &

231 &a bulk mixed layer has not been implemented.")

232 if (is_root_pe()) call mom_error(fatal, &

233 "Either define BULKMIXEDLAYER in MOM_input or use fluxes%buoy &

234 &and a linear equation of state to drive the model.")

235 endif

236

237 tolerance = cs%Tolerance_Ent

238 kmb = gv%nk_rho_varies

239 k2 = max(kmb+1,2) ; kb_min = k2

240 if (.not. cs%bulkmixedlayer) then

241 kb(:) = 1

242 ! These lines fill in values that are arbitrary, but needed because

243 ! they are used to normalize the buoyancy flux in layer nz.

244 do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ; enddo

245 else

246 kb(:) = 0

247 do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ; enddo

248 endif

249

250 if (cs%id_diff_work > 0) allocate(diff_work(g%isd:g%ied,g%jsd:g%jed,nz+1))

251 if (cs%id_Kd > 0) allocate(kd_eff(g%isd:g%ied,g%jsd:g%jed,nz))

252

253 if (associated(tv%eqn_of_state)) then

254 pres(:) = tv%P_Ref

255 else

256 pres(:) = 0.0

257 endif

258 eosdom(:) = eos_domain(g%HI)

259

260 !$OMP parallel do default(private) shared(is,ie,js,je,nz,Kd_Lay,G,GV,US,dt,CS,h,tv, &

261 !$OMP kmb,Angstrom,fluxes,K2,h_neglect,tolerance, &

262 !$OMP ea,eb,Kd_int,Kd_eff,EOSdom,diff_work,g_2dt, kb_out) &

263 !$OMP firstprivate(kb,ds_dsp1,dsp1_ds,pres,kb_min)

264 do j=js,je

265 do i=is,ie ; kb(i) = 1 ; enddo

266

267 if (allocated(tv%SpV_avg)) then

268 do k=1,nz ; do i=is,ie

269 dtkd(i,k) = gv%RZ_to_H * (dt * kd_lay(i,j,k)) / tv%SpV_avg(i,j,k)

270 enddo ; enddo

271 do i=is,ie

272 dtkd_int(i,1) = gv%RZ_to_H * (dt * kd_int(i,j,1)) / tv%SpV_avg(i,j,1)

273 dtkd_int(i,nz+1) = gv%RZ_to_H * (dt * kd_int(i,j,nz+1)) / tv%SpV_avg(i,j,nz)

274 enddo

275 ! Use the mass-weighted average specific volume to translate thicknesses to verti distances.

276 do k=2,nz ; do i=is,ie

277 dtkd_int(i,k) = gv%RZ_to_H * (dt * kd_int(i,j,k)) * &

278 ( (h(i,j,k-1) + h(i,j,k) + 2.0*h_neglect) / &

279 ((h(i,j,k-1)+h_neglect) * tv%SpV_avg(i,j,k-1) + &

280 (h(i,j,k)+h_neglect) * tv%SpV_avg(i,j,k)) )

281 enddo ; enddo

282 else

283 do k=1,nz ; do i=is,ie

284 dtkd(i,k) = gv%Z_to_H * (dt * kd_lay(i,j,k))

285 enddo ; enddo

286 do k=1,nz+1 ; do i=is,ie

287 dtkd_int(i,k) = gv%Z_to_H * (dt * kd_int(i,j,k))

288 enddo ; enddo

289 endif

290

291 do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.0) ; enddo

292 do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; enddo

293 do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo

294

295 if (gv%Boussinesq .or. gv%Semi_Boussinesq) then

296 do k=2,nz-1 ; do i=is,ie

297 ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)

298 enddo ; enddo

299 else ! Use a mathematically equivalent form that avoids any dependency on RHO_0.

300 do k=2,nz-1 ; do i=is,ie

301 ds_dsp1(i,k) = (gv%Rlay(k) - gv%Rlay(k-1)) / (gv%Rlay(k+1) - gv%Rlay(k))

302 enddo ; enddo

303 endif

304

305 if (cs%bulkmixedlayer) then

306 ! This subroutine determines the averaged entrainment across each

307 ! interface and causes thin and relatively light interior layers to be

308 ! entrained by the deepest buffer layer. This also determines kb.

309 call set_ent_bl(h, dtkd_int, tv, kb, kmb, do_i, g, gv, us, cs, j, ent_bl, sref, h_bl)

310

311 do i=is,ie

312 dtkd_kb(i) = 0.0 ; if (kb(i) < nz) dtkd_kb(i) = dtkd(i,kb(i))

313 enddo

314 else

315 do i=is,ie ; ent_bl(i,kmb+1) = 0.0 ; enddo

316 endif

317

318 do k=2,nz-1 ; do i=is,ie

319 dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)

320 i2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))

321 grats(i,k) = 2.0*(2.0+(dsp1_ds(i,k)+ds_dsp1(i,k)))

322 enddo ; enddo

323

324! Determine the maximum flux, maxF, for each of the isopycnal layers.

325! Also determine when the fluxes start entraining

326! from various buffer or mixed layers, where appropriate.

327 if (cs%bulkmixedlayer) then

328 kb_min = nz

329 do i=is,ie

330 htot(i) = h(i,j,1) - angstrom

331 enddo

332 do k=2,kmb ; do i=is,ie

333 htot(i) = htot(i) + (h(i,j,k) - angstrom)

334 enddo ; enddo

335 do i=is,ie

336 max_eakb(i) = max(ent_bl(i,kmb+1) + 0.5*htot(i), htot(i))

337 i_dskbp1(i) = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))

338 zeros(i) = 0.0

339 enddo

340

341 ! Find the maximum amount of entrainment from below that the top

342 ! interior layer could exhibit (maxF(i,kb)), approximating that

343 ! entrainment by (eakb*max(dS_kb/dSkbp1,0)). eakb is in the range

344 ! from 0 to max_eakb.

345 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, max_eakb, kmb, &

346 is, ie, g, gv, cs, maxf_kb, eakb_maxf, do_i, f_kb_maxent)

347 do i=is,ie ; if (kb(i) <= nz) then

348 maxf(i,kb(i)) = max(0.0, maxf_kb(i), f_kb_maxent(i))

349 if ((maxf_kb(i) > f_kb_maxent(i)) .and. (eakb_maxf(i) >= htot(i))) &

350 max_eakb(i) = eakb_maxf(i)

351 endif ; enddo

352

353 do i=is,ie ; ea_kbp1(i) = 0.0 ; eakb(i) = max_eakb(i) ; enddo

354 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &

355 max_eakb, max_eakb, kmb, is, ie, do_i, g, gv, cs, eakb, &

356 error=err_max_eakb0, f_kb=f_kb)

357

358 ! The maximum value of F(kb) between htot and max_eakb determines

359 ! what maxF(kb+1) should be.

360 do i=is,ie ; min_eakb(i) = min(htot(i), max_eakb(i)) ; enddo

361 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, min_eakb, max_eakb, &

362 kmb, is, ie, g, gv, cs, f_kb_maxent, do_i_in=do_i)

363

364 do i=is,ie

365 do_entrain_eakb = .false.

366 ! If error_max_eakb0 < 0, then buffer layers are always all entrained

367 if (do_i(i)) then ; if (err_max_eakb0(i) < 0.0) then

368 do_entrain_eakb = .true.

369 endif ; endif

370

371 if (do_entrain_eakb) then

372 eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)

373 else

374 eakb(i) = 0.0 ; min_eakb(i) = 0.0

375 endif

376 enddo

377

378 ! Find the amount of entrainment of the buffer layers that would occur

379 ! if there were no entrainment by the deeper interior layers. Also find

380 ! how much entrainment of the deeper layers would occur.

381 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &

382 zeros, max_eakb, kmb, is, ie, do_i, g, gv, cs, min_eakb, &

383 error=err_min_eakb0, f_kb=f_kb, err_max_eakb0=err_max_eakb0)

384 ! Error_min_eakb0 should be ~0 unless error_max_eakb0 < 0.

385 do i=is,ie ; if ((kb(i)<nz) .and. (kb_min>kb(i))) kb_min = kb(i) ; enddo

386 else

387 ! Without a bulk mixed layer, surface fluxes are applied in this

388 ! subroutine. (Otherwise, they are handled in mixedlayer.)

389 ! Initially the maximum flux in layer zero is given by the surface

390 ! buoyancy flux. It will later be limited if the surface flux is

391 ! too large. Here buoy is the surface buoyancy flux.

392 do i=is,ie

393 maxf(i,1) = 0.0

394 htot(i) = h(i,j,1) - angstrom

395 enddo

396 if (associated(fluxes%buoy) .and. gv%Boussinesq) then

397 do i=is,ie

398 maxf(i,1) = gv%Z_to_H * (dt*fluxes%buoy(i,j)) / gv%g_prime(2)

399 enddo

400 elseif (associated(fluxes%buoy)) then

401 do i=is,ie

402 maxf(i,1) = (gv%RZ_to_H * 0.5*(gv%Rlay(1) + gv%Rlay(2)) * (dt*fluxes%buoy(i,j))) / &

403 gv%g_prime(2)

404 enddo

405 endif

406 endif

407

408! The following code calculates the maximum flux, maxF, for the interior

409! layers.

410 do k=kb_min,nz-1 ; do i=is,ie

411 if ((k == kb(i)+1) .and. cs%bulkmixedlayer) then

412 maxf(i,k) = ds_dsp1(i,k)*(f_kb_maxent(i) + htot(i))

413 htot(i) = htot(i) + (h(i,j,k) - angstrom)

414 elseif (k > kb(i)) then

415 maxf(i,k) = ds_dsp1(i,k)*(maxf(i,k-1) + htot(i))

416 htot(i) = htot(i) + (h(i,j,k) - angstrom)

417 endif

418 enddo ; enddo

419 do i=is,ie

420 maxf(i,nz) = 0.0

421 if (.not.cs%bulkmixedlayer) then

422 maxf_correct(i) = max(0.0, -(maxf(i,nz-1) + htot(i)))

423 endif

424 htot(i) = h(i,j,nz) - angstrom

425 enddo

426 if (.not.cs%bulkmixedlayer) then

427 do_any = .false. ; do i=is,ie ; if (maxf_correct(i) > 0.0) do_any = .true. ; enddo

428 if (do_any) then

429 do k=nz-1,1,-1 ; do i=is,ie

430 maxf(i,k) = maxf(i,k) + maxf_correct(i)

431 maxf_correct(i) = maxf_correct(i) * dsp1_ds(i,k)

432 enddo ; enddo

433 endif

434 endif

435 do k=nz-1,kb_min,-1 ; do i=is,ie ; if (do_i(i)) then

436 if (k >= kb(i)) then

437 maxf(i,k) = min(maxf(i,k),dsp1_ds(i,k+1)*maxf(i,k+1) + htot(i))

438 htot(i) = htot(i) + (h(i,j,k) - angstrom)

439 endif

440 if (k == kb(i)) then

441 if ((maxf(i,k) < f_kb(i)) .or. (maxf(i,k) < maxf_kb(i)) &

442 .and. (eakb_maxf(i) <= max_eakb(i))) then

443 ! In this case, too much was being entrained by the topmost interior

444 ! layer, even with the minimum initial estimate. The buffer layer

445 ! will always entrain the maximum amount.

446 f_kb(i) = maxf(i,k)

447 if ((f_kb(i) <= maxf_kb(i)) .and. (eakb_maxf(i) <= max_eakb(i))) then

448 eakb(i) = eakb_maxf(i)

449 else

450 eakb(i) = max_eakb(i)

451 endif

452 call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &

453 g, gv, cs, eakb, angstrom)

454 if ((eakb(i) < max_eakb(i)) .or. (eakb(i) < min_eakb(i))) then

455 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, zeros, &

456 eakb, eakb, kmb, i, i, do_i, g, gv, cs, eakb, &

457 error=err_eakb0)

458 if (eakb(i) < max_eakb(i)) then

459 max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)

460 endif

461 if (eakb(i) < min_eakb(i)) then

462 min_eakb(i) = eakb(i) ; err_min_eakb0(i) = err_eakb0(i)

463 endif

464 endif

465 endif

466 endif

467 endif ; enddo ; enddo

468 if (.not.cs%bulkmixedlayer) then

469 do i=is,ie

470 maxf(i,1) = min(maxf(i,1),dsp1_ds(i,2)*maxf(i,2) + htot(i))

471 enddo

472 endif

473

474! The following code provides an initial estimate of the flux in

475! each layer, F. The initial guess for the layer diffusive flux is

476! the smaller of a forward discretization or the maximum diffusive

477! flux starting from zero thickness in one time step without

478! considering adjacent layers.

479 do i=is,ie

480 f(i,1) = maxf(i,1)

481 f(i,nz) = maxf(i,nz) ; minf(i,nz) = 0.0

482 enddo

483 do k=nz-1,k2,-1

484 do i=is,ie

485 if ((k==kb(i)) .and. (do_i(i))) then

486 eakb(i) = min_eakb(i)

487 minf(i,k) = 0.0

488 elseif ((k>kb(i)) .and. (do_i(i))) then

489! Here the layer flux is estimated, assuming no entrainment from

490! the surrounding layers. The estimate is a forward (steady) flux,

491! limited by the maximum flux for a layer starting with zero

492! thickness. This is often a good guess and leads to few iterations.

493 hm = h(i,j,k) + h_neglect

494 ! Note: Tried sqrt((0.5*ds_dsp1(i,k))*dtKd(i,k)) for the second limit,

495 ! but it usually doesn't work as well.

496 f(i,k) = min(maxf(i,k), sqrt(ds_dsp1(i,k)*dtkd(i,k)), &

497 0.5*(ds_dsp1(i,k)+1.0) * (dtkd(i,k) / hm))

498

499! Calculate the minimum flux that can be expected if there is no entrainment

500! from the neighboring layers. The 0.9 is used to give used to give a 10%

501! known error tolerance.

502 fk = dtkd(i,k) * grats(i,k)

503 minf(i,k) = min(maxf(i,k), &

504 0.9*(i2p2dsp1_ds(i,k) * fk / (hm + sqrt(hm*hm + fk))))

505 if (k==kb(i)) minf(i,k) = 0.0 ! BACKWARD COMPATIBILITY - DELETE LATER?

506 else

507 f(i,k) = 0.0

508 minf(i,k) = 0.0

509 endif

510 enddo ! end of i loop

511 enddo ! end of k loop

512

513 ! This is where the fluxes are actually calculated.

514

515 is1 = ie+1 ; ie1 = is-1

516 do i=is,ie ; if (do_i(i)) then ; is1 = i ; exit ; endif ; enddo

517 do i=ie,is,-1 ; if (do_i(i)) then ; ie1 = i ; exit ; endif ; enddo

518

519 if (cs%bulkmixedlayer) then

520 kb_min_act = nz

521 do i=is,ie

522 if (do_i(i) .and. (kb(i) < kb_min_act)) kb_min_act = kb(i)

523 enddo

524 ! Solve for the entrainment rate from above in the topmost interior

525 ! layer, eakb, such that

526 ! eakb*dS_implicit = dt*Kd*dS_layer_implicit / h_implicit.

527 do i=is1,ie1

528 ea_kbp1(i) = 0.0

529 if (do_i(i) .and. (kb(i) < nz)) &

530 ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*f(i,kb(i)+1)

531 enddo

532 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &

533 max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &

534 err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &

535 dfdfm_kb=dfdfm_kb)

536 else

537 kb_min_act = kb_min

538 endif

539

540 do it=0,cs%max_ent_it-1

541 do i=is1,ie1 ; if (do_i(i)) then

542 if (.not.cs%bulkmixedlayer) f(i,1) = min(f(i,1),maxf(i,1))

543 b1(i) = 1.0

544 endif ; enddo ! end of i loop

545

546 ! F_kb has already been found for this iteration, either above or at

547 ! the end of the code for the previous iteration.

548 do k=kb_min_act,nz-1 ; do i=is1,ie1 ; if (do_i(i) .and. (k>=kb(i))) then

549 ! Calculate the flux in layer k.

550 if (cs%bulkmixedlayer .and. (k==kb(i))) then

551 f(i,k) = f_kb(i)

552 dfdfm(i,k) = dfdfm_kb(i)

553 else ! k > kb(i)

554 fprev(i,k) = f(i,k)

555 fm = (f(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*f(i,k+1)

556 fk = grats(i,k)*dtkd(i,k)

557 fr = sqrt(fm*fm + fk)

558

559 if (fm>=0) then

560 f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fm+fr))

561 else

562 f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))

563 endif

564

565 if ((f(i,k) >= maxf(i,k)) .or. (fr == 0.0)) then

566 dfdfm(i,k) = 0.0

567 else

568 dfdfm(i,k) = i2p2dsp1_ds(i,k) * ((fr + fm) / fr)

569 endif

570

571 if (k > k2) then

572 ! This is part of a tridiagonal solver for the actual flux.

573 c1(i,k) = dfdfm(i,k-1)*(dsp1_ds(i,k)*b1(i))

574 b1(i) = 1.0 / (1.0 - c1(i,k)*dfdfm(i,k))

575 f(i,k) = min(b1(i)*(f(i,k)-fprev(i,k)) + fprev(i,k), maxf(i,k))

576 if (f(i,k) >= maxf(i,k)) dfdfm(i,k) = 0.0

577 endif

578 endif

579 endif ; enddo ; enddo

580

581 do k=nz-2,kb_min_act,-1 ; do i=is1,ie1

582 if (do_i(i) .and. (k > kb(i))) &

583 f(i,k) = min((f(i,k)+c1(i,k+1)*(f(i,k+1)-fprev(i,k+1))),maxf(i,k))

584 enddo ; enddo

585

586 if (cs%bulkmixedlayer) then

587 do i=is1,ie1

588 if (do_i(i) .and. (kb(i) < nz)) then

589 ! F will be increased to minF later.

590 ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(f(i,kb(i)+1), minf(i,kb(i)+1))

591 else

592 ea_kbp1(i) = 0.0

593 endif

594 enddo

595 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &

596 max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &

597 err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &

598 dfdfm_kb=dfdfm_kb)

599 do i=is1,ie1

600 if (do_i(i) .and. (kb(i) < nz)) f(i,kb(i)) = f_kb(i)

601 enddo

602 endif

603

604! Determine whether to do another iteration.

605 if (it < cs%max_ent_it-1) then

606

607 reiterate = .false.

608 if (cs%bulkmixedlayer) then ; do i=is1,ie1 ; if (do_i(i)) then

609 eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)

610 endif ; enddo ; endif

611 do i=is1,ie1

612 did_i(i) = do_i(i) ; do_i(i) = .false.

613 enddo

614 do k=kb_min_act,nz-1 ; do i=is1,ie1

615 if (did_i(i) .and. (k >= kb(i))) then

616 if (f(i,k) < minf(i,k)) then

617 f(i,k) = minf(i,k)

618 do_i(i) = .true. ; reiterate = .true.

619 elseif (k > kb(i)) then

620 if ((abs(f(i,k) - fprev(i,k)) > tolerance) .or. &

621 ((h(i,j,k) + ((1.0+dsp1_ds(i,k))*f(i,k) - &

622 (f(i,k-1) + dsp1_ds(i,k+1)*f(i,k+1)))) < 0.9*angstrom)) then

623 do_i(i) = .true. ; reiterate = .true.

624 endif

625 else ! (k == kb(i))

626! A more complicated test is required for the layer beneath the buffer layer,

627! since its flux may be partially used to entrain directly from the mixed layer.

628! Negative fluxes should not occur with the bulk mixed layer.

629 if (h(i,j,k) + ((f(i,k) + eakb(i)) - &

630 (eb_kmb(i) + dsp1_ds(i,k+1)*f(i,k+1))) < 0.9*angstrom) then

631 do_i(i) = .true. ; reiterate = .true.

632 endif

633 endif

634 endif

635 enddo ; enddo

636 if (.not.reiterate) exit

637 endif ! end of if (it < CS%max_ent_it-1)

638 enddo ! end of it loop

639! This is the end of the section that might be iterated.

640

641

642 if (it == (cs%max_ent_it)) then

643 ! Limit the flux so that the layer below is not depleted.

644 ! This should only be applied to the last iteration.

645 do i=is1,ie1 ; if (do_i(i)) then

646 f(i,nz-1) = max(f(i,nz-1), min(minf(i,nz-1), 0.0))

647 if (kb(i) >= nz-1) then ; ea_kbp1(i) = 0.0 ; endif

648 endif ; enddo

649 do k=nz-2,kb_min_act,-1 ; do i=is1,ie1 ; if (do_i(i)) then

650 if (k>kb(i)) then

651 f(i,k) = min(max(minf(i,k),f(i,k)), (dsp1_ds(i,k+1)*f(i,k+1) + &

652 max(((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + &

653 (h(i,j,k+1) - angstrom)), 0.5*(h(i,j,k+1)-angstrom))))

654 elseif (k==kb(i)) then

655 ea_kbp1(i) = dsp1_ds(i,k+1)*f(i,k+1)

656 h_avail = dsp1_ds(i,k+1)*f(i,k+1) + max(0.5*(h(i,j,k+1)-angstrom), &

657 ((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + (h(i,j,k+1) - angstrom)))

658 if ((f(i,k) > 0.0) .and. (f(i,k) > h_avail)) then

659 f_kb(i) = max(0.0, h_avail)

660 f(i,k) = f_kb(i)

661 if ((f_kb(i) < maxf_kb(i)) .and. (eakb_maxf(i) <= eakb(i))) &

662 eakb(i) = eakb_maxf(i)

663 call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &

664 g, gv, cs, eakb)

665 endif

666 endif

667 endif ; enddo ; enddo

668

669

670 if (cs%bulkmixedlayer) then ; do i=is1,ie1

671 if (do_i(i) .and. (kb(i) < nz)) then

672 h_avail = eakb(i) + max(0.5*(h_bl(i,kmb+1)-angstrom), &

673 (f_kb(i)-ea_kbp1(i)) + (h_bl(i,kmb+1)-angstrom))

674 ! Ensure that 0 < eb_kmb < h_avail.

675 ent_bl(i,kmb+1) = min(ent_bl(i,kmb+1),0.5*(eakb(i) + h_avail))

676

677 eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)

678 endif

679 enddo ; endif

680

681 ! Limit the flux so that the layer above is not depleted.

682 do k=kb_min_act+1,nz-1 ; do i=is1,ie1 ; if (do_i(i)) then

683 if ((.not.cs%bulkmixedlayer) .or. (k > kb(i)+1)) then

684 f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + &

685 dsp1_ds(i,k-1)*f(i,k-1)) - f(i,k-2)) + (h(i,j,k-1) - angstrom)))

686 f(i,k) = max(f(i,k),min(minf(i,k),0.0))

687 elseif (k == kb(i)+1) then

688 f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + eakb(i)) - &

689 eb_kmb(i)) + (h(i,j,k-1) - angstrom)))

690 f(i,k) = max(f(i,k),min(minf(i,k),0.0))

691 endif

692 endif ; enddo ; enddo

693 endif ! (it == (CS%max_ent_it))

694

695 call f_to_ent(f, h, kb, kmb, j, g, gv, cs, dsp1_ds, eakb, ent_bl, ea, eb)

696

697 ! Calculate the layer thicknesses after the entrainment to constrain the

698 ! corrective fluxes.

699 if (associated(tv%eqn_of_state)) then

700 do i=is,ie

701 h_guess(i,1) = (h(i,j,1) - angstrom) + (eb(i,j,1) - ea(i,j,2))

702 h_guess(i,nz) = (h(i,j,nz) - angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))

703 if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0

704 if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0

705 enddo

706 do k=2,nz-1 ; do i=is,ie

707 h_guess(i,k) = (h(i,j,k) - angstrom) + ((ea(i,j,k) - eb(i,j,k-1)) + &

708 (eb(i,j,k) - ea(i,j,k+1)))

709 if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0

710 enddo ; enddo

711 if (cs%bulkmixedlayer) then

712 call determine_dskb(h_bl, sref, ent_bl, eakb, is, ie, kmb, g, gv, &

713 .true., ds_kb, ds_anom_lim=ds_anom_lim)

714 do k=nz-1,kb_min,-1

715 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, rcv, tv%eqn_of_state, eosdom)

716 do i=is,ie

717 if ((k>kb(i)) .and. (f(i,k) > 0.0)) then

718 ! Within a time step, a layer may entrain no more than its

719 ! thickness for correction. This limitation should apply

720 ! extremely rarely, but precludes undesirable behavior.

721 ! Note: Corrected a sign/logic error & factor of 2 error, and

722 ! the layers tracked the target density better, mostly due to

723 ! the factor of 2 error.

724 f_cor = h(i,j,k) * min(1.0 , max(-ds_dsp1(i,k), &

725 (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )

726

727 ! Ensure that (1) Entrainments are positive, (2) Corrections in

728 ! a layer cannot deplete the layer itself (very generously), and

729 ! (3) a layer can take no more than a quarter the mass of its

730 ! neighbor.

731 if (f_cor > 0.0) then

732 f_cor = min(f_cor, 0.9*f(i,k), ds_dsp1(i,k)*0.5*h_guess(i,k), &

733 0.25*h_guess(i,k+1))

734 else

735 f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &

736 0.25*ds_dsp1(i,k)*h_guess(i,k-1) )

737 endif

738

739 ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor

740 eb(i,j,k) = eb(i,j,k) + f_cor

741 elseif ((k==kb(i)) .and. (f(i,k) > 0.0)) then

742 ! Rho_cor is the density anomaly that needs to be corrected,

743 ! taking into account that the true potential density of the

744 ! deepest buffer layer is not exactly what is returned as dS_kb.

745 ds_kb_eff = 2.0*ds_kb(i) - ds_anom_lim(i) ! Could be negative!!!

746 rho_cor = h(i,j,k) * (gv%Rlay(k)-rcv(i)) + eakb(i)*ds_anom_lim(i)

747

748 ! Ensure that -.9*eakb < ea_cor < .9*eakb

749 if (abs(rho_cor) < abs(0.9*eakb(i)*ds_kb_eff)) then

750 ea_cor = -rho_cor / ds_kb_eff

751 else

752 ea_cor = sign(0.9*eakb(i),-rho_cor*ds_kb_eff)

753 endif

754

755 if (ea_cor > 0.0) then

756 ! Ensure that -F_cor < 0.5*h_guess

757 ea_cor = min(ea_cor, 0.5*(max_eakb(i) - eakb(i)), &

758 0.5*h_guess(i,k) / (ds_kb(i) * i_dskbp1(i)))

759 else

760 ! Ensure that -ea_cor < 0.5*h_guess & F_cor < 0.25*h_guess(k+1)

761 ea_cor = -min(-ea_cor, 0.5*h_guess(i,k), &

762 0.25*h_guess(i,k+1) / (ds_kb(i) * i_dskbp1(i)))

763 endif

764

765 ea(i,j,k) = ea(i,j,k) + ea_cor

766 eb(i,j,k) = eb(i,j,k) - (ds_kb(i) * i_dskbp1(i)) * ea_cor

767 elseif (k < kb(i)) then

768 ! Repetitive, unless ea(kb) has been corrected.

769 ea(i,j,k) = ea(i,j,k+1)

770 endif

771 enddo

772 enddo

773 do k=kb_min-1,k2,-1 ; do i=is,ie

774 ea(i,j,k) = ea(i,j,k+1)

775 enddo ; enddo

776

777 ! Repetitive, unless ea(kb) has been corrected.

778 k=kmb

779 do i=is,ie

780 ! Do not adjust eb through the base of the buffer layers, but it

781 ! may be necessary to change entrainment from above.

782 h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))

783 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)

784 enddo

785 do k=kmb-1,2,-1 ; do i=is,ie

786 ! Determine the entrainment from below for each buffer layer.

787 eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)

788

789 ! Determine the entrainment from above for each buffer layer.

790 h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))

791 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)

792 enddo ; enddo

793 do i=is,ie

794 eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)

795 enddo

796

797 else ! not bulkmixedlayer

798 do k=k2,nz-1

799 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, rcv, tv%eqn_of_state, eosdom)

800 do i=is,ie ; if (f(i,k) > 0.0) then

801 ! Within a time step, a layer may entrain no more than

802 ! its thickness for correction. This limitation should

803 ! apply extremely rarely, but precludes undesirable

804 ! behavior.

805 f_cor = h(i,j,k) * min(dsp1_ds(i,k) , max(-1.0, &

806 (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )

807

808 ! Ensure that (1) Entrainments are positive, (2) Corrections in

809 ! a layer cannot deplete the layer itself (very generously), and

810 ! (3) a layer can take no more than a quarter the mass of its

811 ! neighbor.

812 if (f_cor >= 0.0) then

813 f_cor = min(f_cor, 0.9*f(i,k), 0.5*dsp1_ds(i,k)*h_guess(i,k), &

814 0.25*h_guess(i,k+1))

815 else

816 f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &

817 0.25*ds_dsp1(i,k)*h_guess(i,k-1) )

818 endif

819

820 ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor

821 eb(i,j,k) = eb(i,j,k) + f_cor

822 endif ; enddo

823 enddo

824 endif

825

826 endif ! associated(tv%eqn_of_state))

827

828 if (cs%id_Kd > 0) then

829 idt = (gv%H_to_m*us%m_to_Z) / dt

830 do k=2,nz-1 ; do i=is,ie

831 if (k<kb(i)) then ; kd_here = 0.0 ; else

832 kd_here = f(i,k) * ( h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &

833 (eb(i,j,k) - ea(i,j,k+1))) ) / (i2p2dsp1_ds(i,k) * grats(i,k))

834 endif

835

836 kd_eff(i,j,k) = max(dtkd(i,k), kd_here)*idt

837 enddo ; enddo

838 do i=is,ie

839 kd_eff(i,j,1) = dtkd(i,1)*idt

840 kd_eff(i,j,nz) = dtkd(i,nz)*idt

841 enddo

842 endif

843

844 if (cs%id_diff_work > 0) then

845 if (gv%Boussinesq .or. .not.associated(tv%eqn_of_state)) then

846 g_2dt = 0.5 * gv%H_to_Z**2 * (gv%g_Earth_Z_T2 / dt)

847 else

848 g_2dt = 0.5 * gv%H_to_RZ**2 * (gv%g_Earth_Z_T2 / dt)

849 endif

850 do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo

851 if (associated(tv%eqn_of_state)) then

852 if (associated(fluxes%p_surf)) then

853 do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ; enddo

854 else

855 do i=is,ie ; pressure(i) = 0.0 ; enddo

856 endif

857 do k=2,nz

858 do i=is,ie ; pressure(i) = pressure(i) + (gv%g_Earth*gv%H_to_RZ)*h(i,j,k-1) ; enddo

859 do i=is,ie

860 if (k==kb(i)) then

861 t_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))

862 s_eos(i) = 0.5*(tv%S(i,j,kmb) + tv%S(i,j,k))

863 else

864 t_eos(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))

865 s_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))

866 endif

867 enddo

868 if (gv%Boussinesq) then

869 call calculate_density_derivs(t_eos, s_eos, pressure, drho_dt, drho_ds, &

870 tv%eqn_of_state, eosdom)

871 do i=is,ie

872 if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,k) = 0.0

873 else

874 if (k==kb(i)) then

875 drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &

876 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))

877 else

878 drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &

879 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))

880 endif

881 diff_work(i,j,k) = g_2dt * drho * &

882 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &

883 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))

884 endif

885 enddo

886 else

887 call calculate_specific_vol_derivs(t_eos, s_eos, pressure, dspv_dt, dspv_ds, &

888 tv%eqn_of_state, eosdom)

889

890 do i=is,ie

891 if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,k) = 0.0

892 else

893 if (k==kb(i)) then

894 dspv = dspv_dt(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &

895 dspv_ds(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))

896 else

897 dspv = dspv_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &

898 dspv_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))

899 endif

900 diff_work(i,j,k) = -g_2dt * dspv * &

901 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &

902 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))

903 endif

904 enddo

905 endif

906 enddo

907 else

908 do k=2,nz ; do i=is,ie

909 diff_work(i,j,k) = g_2dt * (gv%Rlay(k)-gv%Rlay(k-1)) * &

910 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &

911 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))

912 enddo ; enddo

913 endif

914 endif

915

916 do i=is,ie ; kb_out(i,j) = kb(i) ; enddo

917

918 enddo ! end of j loop

919

920! Offer diagnostic fields for averaging.

921 if (cs%id_Kd > 0) call post_data(cs%id_Kd, kd_eff, cs%diag)

922 if (cs%id_Kd > 0) deallocate(kd_eff)

923 if (cs%id_diff_work > 0) call post_data(cs%id_diff_work, diff_work, cs%diag)

924 if (cs%id_diff_work > 0) deallocate(diff_work)

925

926end subroutine entrainment_diffusive

927

928!> This subroutine calculates the actual entrainments (ea and eb) and the

929!! amount of surface forcing that is applied to each layer if there is no bulk

930!! mixed layer.

931subroutine f_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb)

932 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure

933 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure

934 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: F !< The density flux through a layer within

935 !! a time step divided by the density

936 !! difference across the interface below

937 !! the layer [H ~> m or kg m-2].

938 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

939 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]

940 integer, dimension(SZI_(G)), intent(in) :: kb !< The index of the lightest layer denser than

941 !! the deepest buffer layer.

942 integer, intent(in) :: kmb !< The number of mixed and buffer layers.

943 integer, intent(in) :: j !< The meridional index upon which to work.

944 type(entrain_diffusive_cs), intent(in) :: CS !< This module's control structure.

945 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dsp1_ds !< The ratio of coordinate variable

946 !! differences across the interfaces below

947 !! a layer over the difference across the

948 !! interface above the layer [nondim].

949 real, dimension(SZI_(G)), intent(in) :: eakb !< The entrainment from above by the layer

950 !! below the buffer layer [H ~> m or kg m-2].

951 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Ent_bl !< The average entrainment upward and

952 !! downward across each interface around

953 !! the buffer layers [H ~> m or kg m-2].

954 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

955 intent(inout) :: ea !< The amount of fluid entrained from the layer

956 !! above within this time step [H ~> m or kg m-2].

957 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

958 intent(inout) :: eb !< The amount of fluid entrained from the layer

959 !! below within this time step [H ~> m or kg m-2].

960

961 real :: h1 ! The thickness in excess of the minimum that will remain

962 ! after exchange with the layer below [H ~> m or kg m-2].

963 integer :: i, k, is, ie, nz

964

965 is = g%isc ; ie = g%iec ; nz = gv%ke

966

967 do i=is,ie

968 ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0

969 enddo

970 if (cs%bulkmixedlayer) then

971 do i=is,ie

972 eb(i,j,kmb) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)

973 enddo

974 do k=nz-1,kmb+1,-1 ; do i=is,ie

975 if (k > kb(i)) then

976 ! With a bulk mixed layer, surface buoyancy fluxes are applied

977 ! elsewhere, so F should always be nonnegative.

978 ea(i,j,k) = dsp1_ds(i,k)*f(i,k)

979 eb(i,j,k) = f(i,k)

980 elseif (k == kb(i)) then

981 ea(i,j,k) = eakb(i)

982 eb(i,j,k) = f(i,k)

983 elseif (k == kb(i)-1) then

984 ea(i,j,k) = ea(i,j,k+1)

985 eb(i,j,k) = eb(i,j,kmb)

986 else

987 ea(i,j,k) = ea(i,j,k+1)

988 ! Add the entrainment of the thin interior layers to eb going

989 ! up into the buffer layer.

990 eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)

991 endif

992 enddo ; enddo

993 k = kmb

994 do i=is,ie

995 ! Adjust the previously calculated entrainment from below by the deepest

996 ! buffer layer to account for entrainment of thin interior layers .

997 if (kb(i) > kmb+1) &

998 eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)

999

1000 ! Determine the entrainment from above for each buffer layer.

1001 h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))

1002 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)

1003 enddo

1004 do k=kmb-1,2,-1 ; do i=is,ie

1005 ! Determine the entrainment from below for each buffer layer.

1006 eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)

1007

1008 ! Determine the entrainment from above for each buffer layer.

1009 h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))

1010 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)

1011! if (h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)

1012! elseif (Ent_bl(i,K)+0.5*h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)-0.5*h1

1013! else ; ea(i,j,k) = -h1 ; endif

1014 enddo ; enddo

1015 do i=is,ie

1016 eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)

1017 ea(i,j,1) = 0.0

1018 enddo

1019 else ! not BULKMIXEDLAYER

1020 ! Calculate the entrainment by each layer from above and below.

1021 ! Entrainment is always positive, but F may be negative due to

1022 ! surface buoyancy fluxes.

1023 do i=is,ie

1024 ea(i,j,1) = 0.0 ; eb(i,j,1) = max(f(i,1),0.0)

1025 ea(i,j,2) = dsp1_ds(i,2)*f(i,2) - min(f(i,1),0.0)

1026 enddo

1027

1028 do k=2,nz-1 ; do i=is,ie

1029 eb(i,j,k) = max(f(i,k),0.0)

1030 ea(i,j,k+1) = dsp1_ds(i,k+1)*f(i,k+1) - (f(i,k)-eb(i,j,k))

1031 if (ea(i,j,k+1) < 0.0) then

1032 eb(i,j,k) = eb(i,j,k) - ea(i,j,k+1)

1033 ea(i,j,k+1) = 0.0

1034 endif

1035 enddo ; enddo

1036 endif ! end BULKMIXEDLAYER

1037end subroutine f_to_ent

1038

1039!> This subroutine sets the average entrainment across each of the interfaces

1040!! between buffer layers within a timestep. It also causes thin and relatively

1041!! light interior layers to be entrained by the deepest buffer layer.

1042!! Also find the initial coordinate potential densities (Sref) of each layer.

1043subroutine set_ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)

1044 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.

1045 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.

1046 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1047 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]

1048 real, dimension(SZI_(G),SZK_(GV)+1), &

1049 intent(in) :: dtKd_int !< The diapycnal diffusivity across

1050 !! each interface times the time step

1051 !! [H2 ~> m2 or kg2 m-4].

1052 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any

1053 !! available thermodynamic fields. Absent

1054 !! fields have NULL ptrs.

1055 integer, dimension(SZI_(G)), intent(inout) :: kb !< The index of the lightest layer denser

1056 !! than the buffer layer or 1 if there is

1057 !! no buffer layer.

1058 integer, intent(in) :: kmb !< The number of mixed and buffer layers.

1059 logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which

1060 !! i-points to work on.

1061 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1062 type(entrain_diffusive_cs), intent(in) :: CS !< This module's control structure.

1063 integer, intent(in) :: j !< The meridional index upon which to work.

1064 real, dimension(SZI_(G),SZK_(GV)+1), &

1065 intent(out) :: Ent_bl !< The average entrainment upward and

1066 !! downward across each interface around

1067 !! the buffer layers [H ~> m or kg m-2].

1068 real, dimension(SZI_(G),SZK_(GV)), intent(out) :: Sref !< The coordinate potential density minus

1069 !! 1000 for each layer [R ~> kg m-3].

1070 real, dimension(SZI_(G),SZK_(GV)), intent(out) :: h_bl !< The thickness of each layer [H ~> m or kg m-2].

1071

1072! This subroutine sets the average entrainment across each of the interfaces

1073! between buffer layers within a timestep. It also causes thin and relatively

1074! light interior layers to be entrained by the deepest buffer layer.

1075! Also find the initial coordinate potential densities (Sref) of each layer.

1076! Does there need to be limiting when the layers below are all thin?

1077

1078 ! Local variables

1079 real, dimension(SZI_(G)) :: &

1080 b1, d1, & ! Variables used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] and [nondim].

1081 Rcv, & ! Value of the coordinate variable (potential density)

1082 ! based on the simulated T and S and P_Ref [R ~> kg m-3].

1083 pres, & ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].

1084 frac_rem, & ! The fraction of the diffusion remaining [nondim].

1085 h_interior ! The interior thickness available for entrainment [H ~> m or kg m-2].

1086 real, dimension(SZI_(G), SZK_(GV)) :: &

1087 S_est ! An estimate of the coordinate potential density - 1000 after

1088 ! entrainment for each layer [R ~> kg m-3].

1089 real :: dh ! An available thickness [H ~> m or kg m-2].

1090 real :: Kd_x_dt ! The diffusion that remains after thin layers are

1091 ! entrained [H2 ~> m2 or kg2 m-4].

1092 real :: h_neglect ! A thickness that is so small it is usually lost

1093 ! in roundoff and can be neglected [H ~> m or kg m-2].

1094 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state

1095 integer :: i, k, is, ie, nz

1096 is = g%isc ; ie = g%iec ; nz = gv%ke

1097

1098 h_neglect = gv%H_subroundoff

1099

1100 do i=is,ie ; pres(i) = tv%P_Ref ; enddo

1101 eosdom(:) = eos_domain(g%HI)

1102 do k=1,kmb

1103 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, rcv, tv%eqn_of_state, eosdom)

1104 do i=is,ie

1105 h_bl(i,k) = h(i,j,k) + h_neglect

1106 sref(i,k) = rcv(i) - cs%Rho_sig_off

1107 enddo

1108 enddo

1109

1110 do i=is,ie

1111 h_interior(i) = 0.0 ; ent_bl(i,1) = 0.0

1112! if (kb(i) > nz) Ent_bl(i,Kmb+1) = 0.0

1113 enddo

1114

1115 do k=2,kmb ; do i=is,ie

1116 if (do_i(i)) then

1117 ent_bl(i,k) = min(2.0 * dtkd_int(i,k) / (h(i,j,k-1) + h(i,j,k) + h_neglect), cs%max_Ent)

1118 else ; ent_bl(i,k) = 0.0 ; endif

1119 enddo ; enddo

1120

1121 ! Determine the coordinate density of the bottommost buffer layer if there

1122 ! is no entrainment from the layers below. This is a partial solver, based

1123 ! on the first pass of a tridiagonal solver, as the values in the upper buffer

1124 ! layers are not needed.

1125

1126 do i=is,ie

1127 b1(i) = 1.0 / (h_bl(i,1) + ent_bl(i,2))

1128 d1(i) = h_bl(i,1) * b1(i) ! = 1.0 - Ent_bl(i,2)*b1(i)

1129 s_est(i,1) = (h_bl(i,1)*sref(i,1)) * b1(i)

1130 enddo

1131 do k=2,kmb-1 ; do i=is,ie

1132 b1(i) = 1.0 / ((h_bl(i,k) + ent_bl(i,k+1)) + d1(i)*ent_bl(i,k))

1133 d1(i) = (h_bl(i,k) + d1(i)*ent_bl(i,k)) * b1(i) ! = 1.0 - Ent_bl(i,K+1)*b1(i)

1134 s_est(i,k) = (h_bl(i,k)*sref(i,k) + ent_bl(i,k)*s_est(i,k-1)) * b1(i)

1135 enddo ; enddo

1136 do i=is,ie

1137 s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &

1138 (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))

1139 frac_rem(i) = 1.0

1140 enddo

1141

1142 ! Entrain any thin interior layers that are lighter (in the coordinate

1143 ! potential density) than the deepest buffer layer will be, and adjust kb.

1144 do i=is,ie ; kb(i) = nz+1 ; if (do_i(i)) kb(i) = kmb+1 ; enddo

1145

1146 do k=kmb+1,nz ; do i=is,ie ; if (do_i(i)) then

1147 if ((k == kb(i)) .and. (s_est(i,kmb) > (gv%Rlay(k) - cs%Rho_sig_off))) then

1148 if (4.0*dtkd_int(i,kmb+1)*frac_rem(i) > &

1149 (h_bl(i,kmb) + h(i,j,k)) * (h(i,j,k) - gv%Angstrom_H)) then

1150 ! Entrain this layer into the buffer layer and move kb down.

1151 dh = max((h(i,j,k) - gv%Angstrom_H), 0.0)

1152 if (dh > 0.0) then

1153 frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &

1154 (4.0*dtkd_int(i,kmb+1))

1155 sref(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + dh*(gv%Rlay(k)-cs%Rho_sig_off)) / &

1156 (h_bl(i,kmb) + dh)

1157 h_bl(i,kmb) = h_bl(i,kmb) + dh

1158 s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &

1159 (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))

1160 endif

1161 kb(i) = kb(i) + 1

1162 endif

1163 endif

1164 endif ; enddo ; enddo

1165

1166 ! This is where variables are be set up with a different vertical grid

1167 ! in which the (newly?) massless layers are taken out.

1168 do k=nz,kmb+1,-1 ; do i=is,ie

1169 if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-gv%Angstrom_H)

1170 if (k==kb(i)) then

1171 h_bl(i,kmb+1) = h(i,j,k) ; sref(i,kmb+1) = gv%Rlay(k) - cs%Rho_sig_off

1172 elseif (k==kb(i)+1) then

1173 h_bl(i,kmb+2) = h(i,j,k) ; sref(i,kmb+2) = gv%Rlay(k) - cs%Rho_sig_off

1174 endif

1175 enddo ; enddo

1176 do i=is,ie ; if (kb(i) >= nz) then

1177 h_bl(i,kmb+1) = h(i,j,nz)

1178 sref(i,kmb+1) = gv%Rlay(nz) - cs%Rho_sig_off

1179 h_bl(i,kmb+2) = gv%Angstrom_H

1180 sref(i,kmb+2) = sref(i,kmb+1) + (gv%Rlay(nz) - gv%Rlay(nz-1))

1181 endif ; enddo

1182

1183 ! Perhaps we should revisit the way that the average entrainment between the

1184 ! buffer layer and the interior is calculated so that it is not unduly

1185 ! limited when the layers are less than sqrt(Kd * dt) thick?

1186 do i=is,ie ; if (do_i(i)) then

1187 kd_x_dt = frac_rem(i) * dtkd_int(i,kmb+1)

1188 if ((kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0)) then

1189 ent_bl(i,kmb+1) = 0.0

1190 else

1191 ! If the combined layers are exceptionally thin, use sqrt(Kd*dt) as the

1192 ! estimate of the thickness in the denominator of the thickness diffusion.

1193 ent_bl(i,kmb+1) = min(0.5*h_interior(i), sqrt(kd_x_dt), &

1194 kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))

1195 endif

1196 else

1197 ent_bl(i,kmb+1) = 0.0

1198 endif ; enddo

1199

1200end subroutine set_ent_bl

1201

1202!> This subroutine determines the reference density difference between the

1203!! bottommost buffer layer and the first interior after the mixing between mixed

1204!! and buffer layers and mixing with the layer below. Within the mixed and buffer

1205!! layers, entrainment from the layer above is increased when it is necessary to

1206!! keep the layers from developing a negative thickness; otherwise it equals

1207!! Ent_bl. At each interface, the upward and downward fluxes average out to

1208!! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.

1209!! The density difference across the first interior layer may also be returned.

1210!! It could also be limited to avoid negative values or values that greatly

1211!! exceed the density differences across an interface.

1212!! Additionally, the partial derivatives of dSkb and dSlay with E_kb could

1213!! also be returned.

1214subroutine determine_dskb(h_bl, Sref, Ent_bl, E_kb, is, ie, kmb, G, GV, limit, &

1215 dSkb, ddSkb_dE, dSlay, ddSlay_dE, dS_anom_lim, do_i_in)

1216 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.

1217 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid

1218 !! structure.

1219 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h_bl !< Layer thickness [H ~> m or kg m-2]

1220 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Sref !< Reference potential density [R ~> kg m-3]

1221 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Ent_bl !< The average entrainment upward and

1222 !! downward across each interface

1223 !! around the buffer layers [H ~> m or kg m-2].

1224 real, dimension(SZI_(G)), intent(in) :: E_kb !< The entrainment by the top interior

1225 !! layer [H ~> m or kg m-2].

1226 integer, intent(in) :: is !< The start of the i-index range to work on.

1227 integer, intent(in) :: ie !< The end of the i-index range to work on.

1228 integer, intent(in) :: kmb !< The number of mixed and buffer layers.

1229 logical, intent(in) :: limit !< If true, limit dSkb and dSlay to

1230 !! avoid negative values.

1231 real, dimension(SZI_(G)), intent(inout) :: dSkb !< The limited potential density

1232 !! difference across the interface

1233 !! between the bottommost buffer layer

1234 !! and the topmost interior layer. [R ~> kg m-3]

1235 !! dSkb > 0.

1236 real, dimension(SZI_(G)), optional, intent(inout) :: ddSkb_dE !< The partial derivative of dSkb

1237 !! with E [R H-1 ~> kg m-4 or m-1].

1238 real, dimension(SZI_(G)), optional, intent(inout) :: dSlay !< The limited potential density

1239 !! difference across the topmost

1240 !! interior layer. 0 < dSkb [R ~> kg m-3]

1241 real, dimension(SZI_(G)), optional, intent(inout) :: ddSlay_dE !< The partial derivative of dSlay

1242 !! with E [R H-1 ~> kg m-4 or m-1].

1243 real, dimension(SZI_(G)), optional, intent(inout) :: dS_anom_lim !< A limiting value to use for

1244 !! the density anomalies below the

1245 !! buffer layer [R ~> kg m-3].

1246 logical, dimension(SZI_(G)), optional, intent(in) :: do_i_in !< If present, determines which

1247 !! columns are worked on.

1248

1249! Note that dSkb, ddSkb_dE, dSlay, ddSlay_dE, and dS_anom_lim are declared

1250! intent inout because they should not change where do_i_in is false.

1251

1252! This subroutine determines the reference density difference between the

1253! bottommost buffer layer and the first interior after the mixing between mixed

1254! and buffer layers and mixing with the layer below. Within the mixed and buffer

1255! layers, entrainment from the layer above is increased when it is necessary to

1256! keep the layers from developing a negative thickness; otherwise it equals

1257! Ent_bl. At each interface, the upward and downward fluxes average out to

1258! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.

1259! The density difference across the first interior layer may also be returned.

1260! It could also be limited to avoid negative values or values that greatly

1261! exceed the density differences across an interface.

1262! Additionally, the partial derivatives of dSkb and dSlay with E_kb could

1263! also be returned.

1264

1265 ! Local variables

1266 real, dimension(SZI_(G),SZK_(GV)) :: &

1267 b1, c1, & ! b1 [H-1 ~> m-1 or m2 kg-1] and c1 [nondim] are variables used by the tridiagonal solver.

1268 S, dS_dE, & ! The coordinate density [R ~> kg m-3] and its derivative with E [R H-1 ~> kg m-4 or m-1].

1269 ea, dea_dE, & ! The entrainment from above [H ~> m or kg m-2] and its derivative with E [nondim].

1270 eb, deb_dE ! The entrainment from below [H ~> m or kg m-2] and its derivative with E [nondim].

1271 real :: deriv_dSkb(SZI_(G)) ! The limited derivative of the new density difference across the base of

1272 ! the buffer layers with the new density of the bottommost buffer layer [nondim]

1273 real :: d1(SZI_(G)) ! d1 = 1.0-c1 is also used by the tridiagonal solver [nondim].

1274 real :: src ! A source term for dS_dR [R ~> kg m-3].

1275 real :: h1 ! The thickness in excess of the minimum that will remain

1276 ! after exchange with the layer below [H ~> m or kg m-2].

1277 logical, dimension(SZI_(G)) :: do_i

1278 real :: h_neglect ! A thickness that is so small it is usually lost

1279 ! in roundoff and can be neglected [H ~> m or kg m-2].

1280 real :: h_tr ! h_tr is h at tracer points with a tiny thickness

1281 ! added to ensure positive definiteness [H ~> m or kg m-2].

1282 real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].

1283 real :: rat ! A ratio of density differences [nondim]

1284 real :: dS_kbp1 ! The density difference between the top two interior layers [R ~> kg m-3].

1285 real :: IdS_kbp1 ! The inverse of dS_kbp1 [R-1 ~> m3 kg-1]

1286 real :: deriv_dSLay ! The derivative of the projected density difference across the topmost interior

1287 ! layer with the density difference across the interface above it [nondim]

1288 real :: Inv_term ! The inverse of a nondimensional expression [nondim]

1289 real :: f1, df1_drat ! Temporary variables [nondim].

1290 real :: z, dz_drat, f2, df2_dz, expz ! Temporary variables [nondim].

1291 real :: eps_dSLay, eps_dSkb ! Small nondimensional constants [nondim].

1292 integer :: i, k

1293

1294 if (present(ddslay_de) .and. .not.present(dslay)) call mom_error(fatal, &

1295 "In deterimine_dSkb, ddSLay_dE may only be present if dSlay is.")

1296

1297 h_neglect = gv%H_subroundoff

1298

1299 do i=is,ie

1300 ea(i,kmb+1) = e_kb(i) ; dea_de(i,kmb+1) = 1.0

1301 s(i,kmb+1) = sref(i,kmb+1) ; ds_de(i,kmb+1) = 0.0

1302 b1(i,kmb+1) = 0.0

1303 d1(i) = 1.0

1304 do_i(i) = .true.

1305 enddo

1306 if (present(do_i_in)) then

1307 do i=is,ie ; do_i(i) = do_i_in(i) ; enddo

1308 endif

1309 do k=kmb,1,-1 ; do i=is,ie

1310 if (do_i(i)) then

1311 ! The do_i test here is only for efficiency.

1312 ! Determine the entrainment from below for each buffer layer.

1313 if (2.0*ent_bl(i,k+1) > ea(i,k+1)) then

1314 eb(i,k) = 2.0*ent_bl(i,k+1) - ea(i,k+1) ; deb_de(i,k) = -dea_de(i,k+1)

1315 else

1316 eb(i,k) = 0.0 ; deb_de(i,k) = 0.0

1317 endif

1318

1319 ! Determine the entrainment from above for each buffer layer.

1320 h1 = (h_bl(i,k) - gv%Angstrom_H) + (eb(i,k) - ea(i,k+1))

1321 if (h1 >= 0.0) then

1322 ea(i,k) = ent_bl(i,k) ; dea_de(i,k) = 0.0

1323 elseif (ent_bl(i,k) + 0.5*h1 >= 0.0) then

1324 ea(i,k) = ent_bl(i,k) - 0.5*h1

1325 dea_de(i,k) = 0.5*(dea_de(i,k+1) - deb_de(i,k))

1326 else

1327 ea(i,k) = -h1

1328 dea_de(i,k) = dea_de(i,k+1) - deb_de(i,k)

1329 endif

1330 else

1331 ea(i,k) = 0.0 ; dea_de(i,k) = 0.0 ; eb(i,k) = 0.0 ; deb_de(i,k) = 0.0

1332 endif

1333

1334 ! This is the first-pass of a tridiagonal solver for S.

1335 h_tr = h_bl(i,k) + h_neglect

1336 c1(i,k) = ea(i,k+1) * b1(i,k+1)

1337 b_denom_1 = (h_tr + d1(i)*eb(i,k))

1338 b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))

1339 d1(i) = b_denom_1 * b1(i,k)

1340

1341 s(i,k) = (h_tr*sref(i,k) + eb(i,k)*s(i,k+1)) * b1(i,k)

1342 enddo ; enddo

1343 do k=2,kmb ; do i=is,ie

1344 s(i,k) = s(i,k) + c1(i,k-1)*s(i,k-1)

1345 enddo ; enddo

1346

1347 if (present(ddskb_de) .or. present(ddslay_de)) then

1348 ! These two tridiagonal solvers cannot be combined because the solutions for

1349 ! S are required as a source for dS_dE.

1350 do k=kmb,2,-1 ; do i=is,ie

1351 if (do_i(i) .and. (dea_de(i,k) - deb_de(i,k) > 0.0)) then

1352 src = (((s(i,k+1) - sref(i,k)) * (h_bl(i,k) + h_neglect) + &

1353 (s(i,k+1) - s(i,k-1)) * ea(i,k)) * deb_de(i,k) - &

1354 ((sref(i,k) - s(i,k-1)) * h_bl(i,k) + &

1355 (s(i,k+1) - s(i,k-1)) * eb(i,k)) * dea_de(i,k)) / &

1356 ((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))

1357 else ; src = 0.0 ; endif

1358 ds_de(i,k) = (src + eb(i,k)*ds_de(i,k+1)) * b1(i,k)

1359 enddo ; enddo

1360 do i=is,ie

1361 if (do_i(i) .and. (deb_de(i,1) < 0.0)) then

1362 src = (((s(i,2) - sref(i,1)) * (h_bl(i,1) + h_neglect)) * deb_de(i,1)) / &

1363 (h_bl(i,1) + h_neglect + eb(i,1))

1364 else ; src = 0.0 ; endif

1365 ds_de(i,1) = (src + eb(i,1)*ds_de(i,2)) * b1(i,1)

1366 enddo

1367 do k=2,kmb ; do i=is,ie

1368 ds_de(i,k) = ds_de(i,k) + c1(i,k-1)*ds_de(i,k-1)

1369 enddo ; enddo

1370 endif

1371

1372 ! Now, apply any limiting and return the requested variables.

1373

1374 eps_dskb = 1.0e-6 ! Should be a small, nondimensional, positive number.

1375 if (.not.limit) then

1376 do i=is,ie ; if (do_i(i)) then

1377 dskb(i) = sref(i,kmb+1) - s(i,kmb)

1378 endif ; enddo

1379 if (present(ddskb_de)) then ; do i=is,ie ; if (do_i(i)) then

1380 ddskb_de(i) = -1.0*ds_de(i,kmb)

1381 endif ; enddo ; endif

1382

1383 if (present(dslay)) then ; do i=is,ie ; if (do_i(i)) then

1384 dslay(i) = 0.5 * (sref(i,kmb+2) - s(i,kmb))

1385 endif ; enddo ; endif

1386 if (present(ddslay_de)) then ; do i=is,ie ; if (do_i(i)) then

1387 ddslay_de(i) = -0.5*ds_de(i,kmb)

1388 endif ; enddo ; endif

1389 else

1390 do i=is,ie ; if (do_i(i)) then

1391 ! Need to ensure that 0 < dSkb <= S_kb - Sbl

1392 if (sref(i,kmb+1) - s(i,kmb) < eps_dskb*(sref(i,kmb+2) - sref(i,kmb+1))) then

1393 dskb(i) = eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) ; deriv_dskb(i) = 0.0

1394 else

1395 dskb(i) = sref(i,kmb+1) - s(i,kmb) ; deriv_dskb(i) = -1.0

1396 endif

1397 if (present(ddskb_de)) ddskb_de(i) = deriv_dskb(i)*ds_de(i,kmb)

1398 endif ; enddo

1399

1400 if (present(dslay)) then

1401 dz_drat = 1000.0 ! The limit of large dz_drat the same as choosing a

1402 ! Heaviside function.

1403 eps_dslay = 1.0e-10 ! Should be ~= GV%Angstrom_H / sqrt(Kd*dt)

1404 do i=is,ie ; if (do_i(i)) then

1405 ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)

1406 ids_kbp1 = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))

1407 rat = (sref(i,kmb+1) - s(i,kmb)) * ids_kbp1

1408 ! Need to ensure that 0 < dSLay <= 2*dSkb

1409 if (rat < 0.5) then

1410 ! The coefficients here are chosen so that at rat = 0.5, the value (1.5)

1411 ! and first derivative (-0.5) match with the "typical" case (next).

1412 ! The functional form here is arbitrary.

1413 ! f1 provides a reasonable profile that matches the value and derivative

1414 ! of the "typical" case at rat = 0.5, and has a maximum of less than 2.

1415 inv_term = 1.0 / (1.0-rat)

1416 f1 = 2.0 - 0.125*(inv_term**2)

1417 df1_drat = - 0.25*(inv_term**3)

1418

1419 ! f2 ensures that dSLay goes to 0 rapidly if rat is significantly

1420 ! negative.

1421 z = dz_drat * rat + 4.0 ! The 4 here gives f2(0) = 0.982.

1422 if (z >= 18.0) then ; f2 = 1.0 ; df2_dz = 0.0

1423 elseif (z <= -58.0) then ; f2 = eps_dslay ; df2_dz = 0.0

1424 else

1425 expz = exp(z) ; inv_term = 1.0 / (1.0 + expz)

1426 f2 = (eps_dslay + expz) * inv_term

1427 df2_dz = (1.0 - eps_dslay) * expz * inv_term**2

1428 endif

1429

1430 dslay(i) = dskb(i) * f1 * f2

1431 deriv_dslay = deriv_dskb(i) * (f1 * f2) - (dskb(i)*ids_kbp1) * &

1432 (df1_drat*f2 + f1 * dz_drat * df2_dz)

1433 elseif (dskb(i) <= 3.0*ds_kbp1) then

1434 ! This is the "typical" case.

1435 dslay(i) = 0.5 * (dskb(i) + ds_kbp1)

1436 deriv_dslay = 0.5 * deriv_dskb(i) ! = -0.5

1437 else

1438 dslay(i) = 2.0*ds_kbp1

1439 deriv_dslay = 0.0

1440 endif

1441 if (present(ddslay_de)) ddslay_de(i) = deriv_dslay*ds_de(i,kmb)

1442 endif ; enddo

1443 endif ! present(dSlay)

1444 endif ! Not limited.

1445

1446 if (present(ds_anom_lim)) then ; do i=is,ie ; if (do_i(i)) then

1447 ds_anom_lim(i) = max(0.0, eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) - &

1448 (sref(i,kmb+1) - s(i,kmb)) )

1449 endif ; enddo ; endif

1450

1451end subroutine determine_dskb

1452

1453!> Given an entrainment from below for layer kb, determine a consistent

1454!! entrainment from above, such that dSkb * ea_kb = dSkbp1 * F_kb. The input

1455!! value of ea_kb is both the maximum value that can be obtained and the first

1456!! guess of the iterations. Ideally ea_kb should be an under-estimate

1457subroutine f_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &

1458 G, GV, CS, ea_kb, tol_in)

1459 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure

1460 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure

1461 real, dimension(SZI_(G),SZK_(GV)), &

1462 intent(in) :: h_bl !< Layer thickness, with the top interior

1463 !! layer at k-index kmb+1 [H ~> m or kg m-2].

1464 real, dimension(SZI_(G),SZK_(GV)), &

1465 intent(in) :: Sref !< The coordinate reference potential density,

1466 !! with the value of the topmost interior layer

1467 !! at index kmb+1 [R ~> kg m-3].

1468 real, dimension(SZI_(G),SZK_(GV)), &

1469 intent(in) :: Ent_bl !< The average entrainment upward and downward

1470 !! across each interface around the buffer layers,

1471 !! [H ~> m or kg m-2].

1472 real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in reference

1473 !! potential density across the base of the

1474 !! uppermost interior layer [R-1 ~> m3 kg-1].

1475 real, dimension(SZI_(G)), intent(in) :: F_kb !< The entrainment from below by the

1476 !! uppermost interior layer [H ~> m or kg m-2]

1477 integer, intent(in) :: kmb !< The number of mixed and buffer layers.

1478 integer, intent(in) :: i !< The i-index to work on

1479 type(entrain_diffusive_cs), intent(in) :: CS !< This module's control structure.

1480 real, dimension(SZI_(G)), intent(inout) :: ea_kb !< The entrainment from above by the layer below

1481 !! the buffer layer (i.e. layer kb) [H ~> m or kg m-2].

1482 real, optional, intent(in) :: tol_in !< A tolerance for the iterative determination

1483 !! of the entrainment [H ~> m or kg m-2].

1484

1485 real :: max_ea, min_ea ! Bounds on the estimated entraiment [H ~> m or kg m-2]

1486 real :: err, err_min, err_max ! Errors in the mass flux balance [H R ~> kg m-2 or kg2 m-5]

1487 real :: derr_dea ! The change in error with the change in ea [R ~> kg m-3]

1488 real :: val ! An estimate mass flux [H R ~> kg m-2 or kg2 m-5]

1489 real :: tolerance, tol1 ! Tolerances for the determination of the entrainment [H ~> m or kg m-2]

1490 real :: ea_prev ! A previous estimate of ea_kb [H ~> m or kg m-2]

1491 real :: dS_kbp1 ! The density difference between two interior layers [R ~> kg m-3]

1492 real :: dS_kb(SZI_(G)) ! The limited potential density difference across the interface

1493 ! between the bottommost buffer layer and the topmost interior layer [R ~> kg m-3]

1494 real :: maxF(SZI_(G)) ! The maximum value of F (the density flux divided by density

1495 ! differences) found in the range min_ent < ent < max_ent [H ~> m or kg m-2].

1496 real :: ent_maxF(SZI_(G)) ! The value of entrainment that gives maxF [H ~> m or kg m-2]

1497 real :: zeros(SZI_(G)) ! An array of zero entrainments [H ~> m or kg m-2]

1498 real :: ddSkb_dE(SZI_(G)) ! The partial derivative of dS_kb with ea_kb [R H-1 ~> kg m-4 or m-1]

1499 logical :: bisect_next, Newton ! These indicate what method the next iteration should use

1500 integer :: it

1501 integer, parameter :: MAXIT = 30

1502

1503 ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)

1504 max_ea = ea_kb(i) ; min_ea = 0.0

1505 val = ds_kbp1 * f_kb(i)

1506 err_min = -val

1507

1508 tolerance = cs%Tolerance_Ent

1509 if (present(tol_in)) tolerance = tol_in

1510 bisect_next = .true.

1511

1512 call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &

1513 ds_kb, ddskb_de)

1514

1515 err = ds_kb(i) * ea_kb(i) - val

1516 derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)

1517 ! Return if Newton's method on the first guess would give a tolerably small

1518 ! change in the value of ea_kb.

1519 if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea))) return

1520

1521 if (err == 0.0) then ; return ! The exact solution on the first guess...

1522 elseif (err > 0.0) then ! The root is properly bracketed.

1523 max_ea = ea_kb(i) ; err_max = err

1524 ! Use Newton's method (if it stays bounded) or the false position method

1525 ! to find the next value.

1526 if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i) - min_ea) > err) .and. &

1527 (derr_dea*(max_ea - ea_kb(i)) > -1.0*err)) then

1528 ea_kb(i) = ea_kb(i) - err / derr_dea

1529 else ! Use the bisection for the next guess.

1530 ea_kb(i) = 0.5*(max_ea+min_ea)

1531 endif

1532 else

1533 ! Try to bracket the root first. If unable to bracket the root, return

1534 ! the maximum.

1535 zeros(i) = 0.0

1536 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, ea_kb, &

1537 kmb, i, i, g, gv, cs, maxf, ent_maxf, f_thresh=f_kb)

1538 err_max = ds_kbp1 * maxf(i) - val

1539 ! If err_max is negative, there is no good solution, so use the maximum

1540 ! value of F in the valid range.

1541 if (err_max <= 0.0) then

1542 ea_kb(i) = ent_maxf(i) ; return

1543 else

1544 max_ea = ent_maxf(i)

1545 ea_kb(i) = 0.5*(max_ea+min_ea) ! Use bisection for the next guess.

1546 endif

1547 endif

1548

1549 ! Exit if the range between max_ea and min_ea already acceptable.

1550 ! if (abs(max_ea - min_ea) < 0.1*tolerance) return

1551

1552 do it = 1, maxit

1553 call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &

1554 ds_kb, ddskb_de)

1555

1556 err = ds_kb(i) * ea_kb(i) - val

1557 derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)

1558

1559 ea_prev = ea_kb(i)

1560 ! Use Newton's method or the false position method to find the next value.

1561 newton = .false.

1562 if (err > 0.0) then

1563 max_ea = ea_kb(i) ; err_max = err

1564 if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) newton = .true.

1565 else

1566 min_ea = ea_kb(i) ; err_min = err

1567 if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) newton = .true.

1568 endif

1569

1570 if (newton) then

1571 ea_kb(i) = ea_kb(i) - err / derr_dea

1572 elseif (bisect_next) then ! Use bisection to reduce the range.

1573 ea_kb(i) = 0.5*(max_ea+min_ea)

1574 bisect_next = .false.

1575 else ! Use the false-position method for the next guess.

1576 ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))

1577 bisect_next = .true.

1578 endif

1579

1580 tol1 = tolerance ; if (err > 0.0) tol1 = 0.099*tolerance

1581 if (ds_kb(i) <= ds_kbp1) then

1582 if (abs(ea_kb(i) - ea_prev) <= tol1) return

1583 else

1584 if (ds_kbp1*abs(ea_kb(i) - ea_prev) <= ds_kb(i)*tol1) return

1585 endif

1586 enddo

1587

1588end subroutine f_kb_to_ea_kb

1589

1590

1591!> This subroutine determines the entrainment from above by the top interior

1592!! layer (labeled kb elsewhere) given an entrainment by the layer below it,

1593!! constrained to be within the provided bounds.

1594subroutine determine_ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &

1595 min_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, Ent, &

1596 error, err_min_eakb0, err_max_eakb0, F_kb, dFdfm_kb)

1597 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.

1598 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.

1599 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h_bl !< Layer thickness, with the top interior

1600 !! layer at k-index kmb+1 [H ~> m or kg m-2].

1601 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Sref !< The coordinate reference potential

1602 !! density, with the value of the

1603 !! topmost interior layer at layer

1604 !! kmb+1 [R ~> kg m-3].

1605 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: Ent_bl !< The average entrainment upward and

1606 !! downward across each interface around

1607 !! the buffer layers [H ~> m or kg m-2].

1608 real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in

1609 !! reference potential density across

1610 !! the base of the uppermost interior

1611 !! layer [R-1 ~> m3 kg-1].

1612 real, dimension(SZI_(G)), intent(in) :: dtKd_kb !< The diapycnal diffusivity in the top

1613 !! interior layer times the time step

1614 !! [H2 ~> m2 or kg2 m-4].

1615 real, dimension(SZI_(G)), intent(in) :: ea_kbp1 !< The entrainment from above by layer

1616 !! kb+1 [H ~> m or kg m-2].

1617 real, dimension(SZI_(G)), intent(in) :: min_eakb !< The minimum permissible rate of

1618 !! entrainment [H ~> m or kg m-2].

1619 real, dimension(SZI_(G)), intent(in) :: max_eakb !< The maximum permissible rate of

1620 !! entrainment [H ~> m or kg m-2].

1621 integer, intent(in) :: kmb !< The number of mixed and buffer layers.

1622 integer, intent(in) :: is !< The start of the i-index range to work on.

1623 integer, intent(in) :: ie !< The end of the i-index range to work on.

1624 logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which

1625 !! i-points to work on.

1626 type(entrain_diffusive_cs), intent(in) :: CS !< This module's control structure.

1627 real, dimension(SZI_(G)), intent(inout) :: Ent !< The entrainment rate of the uppermost

1628 !! interior layer [H ~> m or kg m-2].

1629 !! The input value is the first guess.

1630 real, dimension(SZI_(G)), optional, intent(out) :: error !< The error (locally defined in this

1631 !! routine) associated with the returned

1632 !! solution [H2 ~> m2 or kg2 m-4]

1633 real, dimension(SZI_(G)), optional, intent(in) :: err_min_eakb0 !< The errors (locally defined)

1634 !! associated with min_eakb when ea_kbp1 = 0,

1635 !! returned from a previous call to this

1636 !! subroutine [H2 ~> m2 or kg2 m-4].

1637 real, dimension(SZI_(G)), optional, intent(in) :: err_max_eakb0 !< The errors (locally defined)

1638 !! associated with min_eakb when ea_kbp1 = 0,

1639 !! returned from a previous call to this

1640 !! subroutine [H2 ~> m2 or kg2 m-4].

1641 real, dimension(SZI_(G)), optional, intent(out) :: F_kb !< The entrainment from below by the

1642 !! uppermost interior layer

1643 !! corresponding to the returned

1644 !! value of Ent [H ~> m or kg m-2].

1645 real, dimension(SZI_(G)), optional, intent(out) :: dFdfm_kb !< The partial derivative of F_kb with

1646 !! ea_kbp1 [nondim].

1647

1648! This subroutine determines the entrainment from above by the top interior

1649! layer (labeled kb elsewhere) given an entrainment by the layer below it,

1650! constrained to be within the provided bounds.

1651

1652 ! Local variables

1653 real, dimension(SZI_(G)) :: &

1654 dS_kb, & ! The coordinate-density difference between the

1655 ! layer kb and deepest buffer layer, limited to

1656 ! ensure that it is positive [R ~> kg m-3].

1657 ds_lay, & ! The coordinate-density difference across layer

1658 ! kb, limited to ensure that it is positive and not

1659 ! too much bigger than dS_kb or dS_kbp1 [R ~> kg m-3].

1660 ddskb_de, ddslay_de, & ! The derivatives of dS_kb and dS_Lay with E

1661 ! [R H-1 ~> kg m-4 or m-1].

1662 derror_de, & ! The derivative of err with E [H ~> m or kg m-2].

1663 err, & ! The "error" whose zero is being sought [H2 ~> m2 or kg2 m-4].

1664 e_min, e_max, & ! The minimum and maximum values of E [H ~> m or kg m-2].

1665 error_mine, error_maxe ! err when E = E_min or E = E_max [H2 ~> m2 or kg2 m-4].

1666 real :: err_est ! An estimate of what err will be [H2 ~> m2 or kg2 m-4].

1667 real :: eL ! 1 or 0, depending on whether increases in E lead

1668 ! to decreases in the entrainment from below by the

1669 ! deepest buffer layer [nondim].

1670 real :: fa ! Temporary variable used to calculate err [nondim].

1671 real :: fk ! Temporary variable used to calculate err [H2 ~> m2 or kg2 m-4].

1672 real :: fm, fr ! Temporary variables used to calculate err [H ~> m or kg m-2].

1673 real :: tolerance ! The tolerance within which E must be converged [H ~> m or kg m-2].

1674 real :: E_prev ! The previous value of E [H ~> m or kg m-2].

1675 logical, dimension(SZI_(G)) :: false_position ! If true, the false position

1676 ! method might be used for the next iteration.

1677 logical, dimension(SZI_(G)) :: redo_i ! If true, more work is needed on this column.

1678 logical :: do_any

1679 real :: large_err ! A large error measure [H2 ~> m2 or kg2 m-4].

1680 integer :: i, it

1681 integer, parameter :: MAXIT = 30

1682

1683 if (.not.cs%bulkmixedlayer) then

1684 call mom_error(fatal, "determine_Ea_kb should not be called "//&

1685 "unless BULKMIXEDLAYER is defined.")

1686 endif

1687 tolerance = cs%Tolerance_Ent

1688 large_err = gv%m_to_H**2 * 1.0e30

1689

1690 do i=is,ie ; redo_i(i) = do_i(i) ; enddo

1691

1692 do i=is,ie ; if (do_i(i)) then

1693 ! The first guess of Ent was the value from the previous iteration.

1694

1695 ! These were previously calculated and provide good limits and estimates

1696 ! of the errors there. By construction the errors increase with R*ea_kbp1.

1697 e_min(i) = min_eakb(i) ; e_max(i) = max_eakb(i)

1698 error_mine(i) = -large_err ; error_maxe(i) = large_err

1699 false_position(i) = .true. ! Used to alternate between false_position and

1700 ! bisection when Newton's method isn't working.

1701 if (present(err_min_eakb0)) error_mine(i) = err_min_eakb0(i) - e_min(i) * ea_kbp1(i)

1702 if (present(err_max_eakb0)) error_maxe(i) = err_max_eakb0(i) - e_max(i) * ea_kbp1(i)

1703

1704 if ((error_maxe(i) <= 0.0) .or. (error_mine(i) >= 0.0)) then

1705 ! The root is not bracketed and one of the limiting values should be used.

1706 if (error_maxe(i) <= 0.0) then

1707 ! The errors decrease with E*ea_kbp1, so E_max is the best solution.

1708 ent(i) = e_max(i) ; err(i) = error_maxe(i)

1709 else ! error_minE >= 0 is equivalent to ea_kbp1 = 0.0.

1710 ent(i) = e_min(i) ; err(i) = error_mine(i)

1711 endif

1712 derror_de(i) = 0.0

1713 redo_i(i) = .false.

1714 endif

1715 endif ; enddo ! End of i-loop

1716

1717 do it = 1,maxit

1718 do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo

1719 if (.not.do_any) exit

1720 call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., ds_kb, &

1721 ddskb_de, ds_lay, ddslay_de, do_i_in=redo_i)

1722 do i=is,ie ; if (redo_i(i)) then

1723 ! The correct root is bracketed between E_min and E_max.

1724 ! Note the following limits: Ent >= 0 ; fa > 1 ; fk > 0

1725 el = 0.0 ; if (2.0*ent_bl(i,kmb+1) >= ent(i)) el = 1.0

1726 fa = (1.0 + el) + ds_kb(i)*i_dskbp1(i)

1727 fk = dtkd_kb(i) * (ds_lay(i)/ds_kb(i))

1728 fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + el*2.0*ent_bl(i,kmb+1)

1729 if (fm > -gv%Angstrom_H) fm = fm + gv%Angstrom_H ! This could be smooth if need be.

1730 err(i) = (fa * ent(i)**2 - fm * ent(i)) - fk

1731 derror_de(i) = ((2.0*fa + (ddskb_de(i)*i_dskbp1(i))*ent(i))*ent(i) - fm) - &

1732 dtkd_kb(i) * (ddslay_de(i)*ds_kb(i) - ddskb_de(i)*ds_lay(i))/(ds_kb(i)**2)

1733

1734 if (err(i) == 0.0) then

1735 redo_i(i) = .false. ; cycle

1736 elseif (err(i) > 0.0) then

1737 e_max(i) = ent(i) ; error_maxe(i) = err(i)

1738 else

1739 e_min(i) = ent(i) ; error_mine(i) = err(i)

1740 endif

1741

1742 e_prev = ent(i)

1743 if ((it == 1) .or. (derror_de(i) <= 0.0)) then

1744 ! Assuming that the coefficients of the quadratic equation are correct

1745 ! will usually give a very good first guess. Also, if derror_dE < 0.0,

1746 ! R is on the wrong side of the approximate parabola. In either case,

1747 ! try assuming that the error is approximately a parabola and solve.

1748 fr = sqrt(fm**2 + 4.0*fa*fk)

1749 if (fm >= 0.0) then

1750 ent(i) = (fm + fr) / (2.0 * fa)

1751 else

1752 ent(i) = (2.0 * fk) / (fr - fm)

1753 endif

1754 ! But make sure that the root stays bracketed, bisecting if needed.

1755 if ((ent(i) > e_max(i)) .or. (ent(i) < e_min(i))) &

1756 ent(i) = 0.5*(e_max(i) + e_min(i))

1757 elseif (((e_max(i)-ent(i))*derror_de(i) > -err(i)) .and. &

1758 ((ent(i)-e_min(i))*derror_de(i) > err(i)) ) then

1759 ! Use Newton's method for the next estimate, provided it will

1760 ! remain bracketed between Rmin and Rmax.

1761 ent(i) = ent(i) - err(i) / derror_de(i)

1762 elseif (false_position(i) .and. &

1763 (error_maxe(i) - error_mine(i) < 0.9*large_err)) then

1764 ! Use the false position method if there are decent error estimates.

1765 ent(i) = e_min(i) + (e_max(i)-e_min(i)) * &

1766 (-error_mine(i)/(error_maxe(i) - error_mine(i)))

1767 false_position(i) = .false.

1768 else ! Bisect as a last resort or if the false position method was used last.

1769 ent(i) = 0.5*(e_max(i) + e_min(i))

1770 false_position(i) = .true.

1771 endif

1772

1773 if (abs(e_prev - ent(i)) < tolerance) then

1774 err_est = err(i) + (ent(i) - e_prev) * derror_de(i)

1775 if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &

1776 (abs(err_est) < abs(tolerance*derror_de(i)))) redo_i(i) = .false.

1777 endif

1778

1779 endif ; enddo ! End of i-loop

1780 enddo ! End of iterations to determine Ent(i).

1781

1782 ! Update the value of dS_kb for consistency with Ent.

1783 if (present(f_kb) .or. present(dfdfm_kb)) &

1784 call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., &

1785 ds_kb, do_i_in=do_i)

1786

1787 if (present(f_kb)) then ; do i=is,ie ; if (do_i(i)) then

1788 f_kb(i) = ent(i) * (ds_kb(i) * i_dskbp1(i))

1789 endif ; enddo ; endif

1790 if (present(error)) then ; do i=is,ie ; if (do_i(i)) then

1791 error(i) = err(i)

1792 endif ; enddo ; endif

1793 if (present(dfdfm_kb)) then ; do i=is,ie ; if (do_i(i)) then

1794 ! derror_dE and ddSkb_dE are _not_ recalculated here, since dFdfm_kb is

1795 ! only used in Newton's method, and slightly increasing the accuracy of the

1796 ! estimate is unlikely to speed convergence.

1797 if (derror_de(i) > 0.0) then

1798 dfdfm_kb(i) = ((ds_kb(i) + ent(i) * ddskb_de(i)) * i_dskbp1(i)) * &

1799 (ent(i) / derror_de(i))

1800 else ! Use Adcroft's division by 0 convention.

1801 dfdfm_kb(i) = 0.0

1802 endif

1803 endif ; enddo ; endif

1804

1805end subroutine determine_ea_kb

1806

1807!> Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.

1808subroutine find_maxf_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, &

1809 kmb, is, ie, G, GV, CS, maxF, ent_maxF, do_i_in, &

1810 F_lim_maxent, F_thresh)

1811 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.

1812 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.

1813 real, dimension(SZI_(G),SZK_(GV)), &

1814 intent(in) :: h_bl !< Layer thickness [H ~> m or kg m-2]

1815 real, dimension(SZI_(G),SZK_(GV)), &

1816 intent(in) :: Sref !< Reference potential density [R ~> kg m-3].

1817 real, dimension(SZI_(G),SZK_(GV)), &

1818 intent(in) :: Ent_bl !< The average entrainment upward and

1819 !! downward across each interface around

1820 !! the buffer layers [H ~> m or kg m-2].

1821 real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in

1822 !! reference potential density across the

1823 !! base of the uppermost interior layer

1824 !! [R-1 ~> m3 kg-1].

1825 real, dimension(SZI_(G)), intent(in) :: min_ent_in !< The minimum value of ent to search,

1826 !! [H ~> m or kg m-2].

1827 real, dimension(SZI_(G)), intent(in) :: max_ent_in !< The maximum value of ent to search,

1828 !! [H ~> m or kg m-2].

1829 integer, intent(in) :: kmb !< The number of mixed and buffer layers.

1830 integer, intent(in) :: is !< The start of the i-index range to work on.

1831 integer, intent(in) :: ie !< The end of the i-index range to work on.

1832 type(entrain_diffusive_cs), intent(in) :: CS !< This module's control structure.

1833 real, dimension(SZI_(G)), intent(out) :: maxF !< The maximum value of F

1834 !! = ent*ds_kb*I_dSkbp1 found in the range

1835 !! min_ent < ent < max_ent [H ~> m or kg m-2].

1836 real, dimension(SZI_(G)), &

1837 optional, intent(out) :: ent_maxF !< The value of ent at that maximum [H ~> m or kg m-2].

1838 logical, dimension(SZI_(G)), &

1839 optional, intent(in) :: do_i_in !< A logical array indicating which columns

1840 !! to work on.

1841 real, dimension(SZI_(G)), &

1842 optional, intent(out) :: F_lim_maxent !< If present, do not apply the limit in

1843 !! finding the maximum value, but return the

1844 !! limited value at ent=max_ent_in in this

1845 !! array [H ~> m or kg m-2].

1846 real, dimension(SZI_(G)), &

1847 optional, intent(in) :: F_thresh !< If F_thresh is present, return the first value

1848 !! found that has F > F_thresh [H ~> m or kg m-2], or

1849 !! the maximum root if it is absent.

1850

1851! Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.

1852! ds_kb may itself be limited to positive values in determine_dSkb, which gives

1853! the prospect of two local maxima in the range - one at max_ent_in with that

1854! minimum value of ds_kb, and the other due to the unlimited (potentially

1855! negative) value. It is faster to find the true maximum by first finding the

1856! unlimited maximum and comparing it to the limited value at max_ent_in.

1857 real, dimension(SZI_(G)) :: &

1858 ent, & ! The updated estimate of the entrainment [H ~> m or kg m-2]

1859 minent, maxent, ent_best, & ! Various previous estimates of the entrainment [H ~> m or kg m-2]

1860 F_max_ent_in, & ! The value of F that gives the input maximum value of ent [H ~> m or kg m-2]

1861 F_maxent, F_minent, F, F_best, & ! Various estimates of F [H ~> m or kg m-2]

1862 dF_dent, dF_dE_max, dF_dE_min, dF_dE_best, & ! Various derivatives of F with ent [nondim]

1863 dS_kb, & ! The density difference across the interface between the bottommost

1864 ! buffer layer and the topmost interior layer [R ~> kg m-3]

1865 ds_kb_lim, ds_anom_lim, & ! Various limits on dS_kb [R ~> kg m-3]

1866 ddskb_de, & ! The partial derivative of dS_kb with ent [R H-1 ~> kg m-4 or m-1].

1867 chg_prev, chg_pre_prev ! Changes in estimates of the entrainment from previous iterations [H ~> m or kg m-2]

1868 real :: dF_dE_mean, maxslope, minslope ! Various derivatives of F with ent [nondim]

1869 real :: tolerance ! The tolerance within which ent must be converged [H ~> m or kg m-2]

1870 real :: ratio_select_end, rat ! Fractional changes in the value of ent to use for the next iteration

1871 ! relative to its bounded range [nondim]

1872 real :: max_chg, min_chg, chg1, chg2, chg ! Changes in entrainment estimates [H ~> m or kg m-2]

1873 logical, dimension(SZI_(G)) :: do_i, last_it, need_bracket, may_use_best

1874 logical :: doany, OK1, OK2, bisect, new_min_bound

1875 integer :: i, it, is1, ie1

1876 integer, parameter :: MAXIT = 20

1877

1878 tolerance = cs%Tolerance_Ent

1879

1880 if (present(do_i_in)) then

1881 do i=is,ie ; do_i(i) = do_i_in(i) ; enddo

1882 else

1883 do i=is,ie ; do_i(i) = .true. ; enddo

1884 endif

1885

1886 ! The most likely value is at max_ent.

1887 call determine_dskb(h_bl, sref, ent_bl, max_ent_in, is, ie, kmb, g, gv, .false., &

1888 ds_kb, ddskb_de, ds_anom_lim=ds_anom_lim)

1889 ie1 = is-1 ; doany = .false.

1890 do i=is,ie

1891 ds_kb_lim(i) = ds_kb(i) + ds_anom_lim(i)

1892 f_max_ent_in(i) = max_ent_in(i)*ds_kb_lim(i)*i_dskbp1(i)

1893 maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)

1894 if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i))) then

1895 f_best(i) = max_ent_in(i)*ds_kb(i)*i_dskbp1(i)

1896 ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)

1897 do_i(i) = .false.

1898 else

1899 f_maxent(i) = maxent(i) * ds_kb(i) * i_dskbp1(i)

1900 df_de_max(i) = (ds_kb(i) + maxent(i)*ddskb_de(i)) * i_dskbp1(i)

1901 doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.

1902 endif

1903 enddo

1904

1905 if (doany) then

1906 ie1 = is-1 ; do i=is,ie ; if (do_i(i)) ie1 = i ; enddo

1907 do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo

1908 ! Find the value of F and its derivative at min_ent.

1909 call determine_dskb(h_bl, sref, ent_bl, minent, is1, ie1, kmb, g, gv, .false., &

1910 ds_kb, ddskb_de, do_i_in=do_i)

1911 do i=is1,ie1 ; if (do_i(i)) then

1912 f_minent(i) = minent(i) * ds_kb(i) * i_dskbp1(i)

1913 df_de_min(i) = (ds_kb(i) + minent(i)*ddskb_de(i)) * i_dskbp1(i)

1914 endif ; enddo

1915

1916 ratio_select_end = 0.9

1917 do it=1,maxit

1918 ratio_select_end = 0.5*ratio_select_end

1919 do i=is1,ie1 ; if (do_i(i)) then

1920 if (need_bracket(i)) then

1921 df_de_mean = (f_maxent(i) - f_minent(i)) / (maxent(i) - minent(i))

1922 maxslope = max(df_de_mean, df_de_min(i), df_de_max(i))

1923 minslope = min(df_de_mean, df_de_min(i), df_de_max(i))

1924 if (f_minent(i) >= f_maxent(i)) then

1925 if (df_de_min(i) > 0.0) then ; rat = 0.02 ! A small step should bracket the solution.

1926 elseif (maxslope < ratio_select_end*minslope) then

1927 ! The maximum of F is at minent.

1928 f_best(i) = f_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0

1929 do_i(i) = .false.

1930 else ; rat = 0.382 ; endif ! Use the golden ratio

1931 else

1932 if (df_de_max(i) < 0.0) then ; rat = 0.98 ! A small step should bracket the solution.

1933 elseif (minslope > ratio_select_end*maxslope) then

1934 ! The maximum of F is at maxent.

1935 f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0

1936 do_i(i) = .false.

1937 else ; rat = 0.618 ; endif ! Use the golden ratio

1938 endif

1939

1940 if (rat >= 0.0) ent(i) = rat*maxent(i) + (1.0-rat)*minent(i)

1941 if (((maxent(i) - minent(i)) < tolerance) .or. (it==maxit)) &

1942 last_it(i) = .true.

1943 else ! The maximum is bracketed by minent, ent_best, and maxent.

1944 chg1 = 2.0*(maxent(i) - minent(i)) ; chg2 = chg1

1945 if (df_de_best(i) > 0) then

1946 max_chg = maxent(i) - ent_best(i) ; min_chg = 0.0

1947 else

1948 max_chg = 0.0 ; min_chg = minent(i) - ent_best(i) ! < 0

1949 endif

1950 if (max_chg - min_chg < 2.0*tolerance) last_it(i) = .true.

1951 if (df_de_max(i) /= df_de_best(i)) &

1952 chg1 = (maxent(i) - ent_best(i))*df_de_best(i) / &

1953 (df_de_best(i) - df_de_max(i))

1954 if (df_de_min(i) /= df_de_best(i)) &

1955 chg2 = (minent(i) - ent_best(i))*df_de_best(i) / &

1956 (df_de_best(i) - df_de_min(i))

1957 ok1 = ((chg1 < max_chg) .and. (chg1 > min_chg))

1958 ok2 = ((chg2 < max_chg) .and. (chg2 > min_chg))

1959 if (.not.(ok1 .or. ok2)) then ; bisect = .true. ; else

1960 if (ok1 .and. ok2) then ! Take the acceptable smaller change.

1961 chg = chg1 ; if (abs(chg2) < abs(chg1)) chg = chg2

1962 elseif (ok1) then ; chg = chg1

1963 else ; chg = chg2 ; endif

1964 if (abs(chg) > 0.5*abs(chg_pre_prev(i))) then ; bisect = .true.

1965 else ; bisect = .false. ; endif

1966 endif

1967 chg_pre_prev(i) = chg_prev(i)

1968 if (bisect) then

1969 if (df_de_best(i) > 0.0) then

1970 ent(i) = 0.5*(maxent(i) + ent_best(i))

1971 chg_prev(i) = 0.5*(maxent(i) - ent_best(i))

1972 else

1973 ent(i) = 0.5*(minent(i) + ent_best(i))

1974 chg_prev(i) = 0.5*(minent(i) - ent_best(i))

1975 endif

1976 else

1977 if (abs(chg) < tolerance) chg = sign(tolerance,chg)

1978 ent(i) = ent_best(i) + chg

1979 chg_prev(i) = chg

1980 endif

1981 endif

1982 endif ; enddo

1983

1984 if (mod(it,3) == 0) then ! Re-determine the loop bounds.

1985 ie1 = is-1 ; do i=is1,ie ; if (do_i(i)) ie1 = i ; enddo

1986 do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo

1987 endif

1988

1989 call determine_dskb(h_bl, sref, ent_bl, ent, is1, ie1, kmb, g, gv, .false., &

1990 ds_kb, ddskb_de, do_i_in=do_i)

1991 do i=is1,ie1 ; if (do_i(i)) then

1992 f(i) = ent(i)*ds_kb(i)*i_dskbp1(i)

1993 df_dent(i) = (ds_kb(i) + ent(i)*ddskb_de(i)) * i_dskbp1(i)

1994 endif ; enddo

1995

1996 if (present(f_thresh)) then ; do i=is1,ie1 ; if (do_i(i)) then

1997 if (f(i) >= f_thresh(i)) then

1998 f_best(i) = f(i) ; ent_best(i) = ent(i) ; do_i(i) = .false.

1999 endif

2000 endif ; enddo ; endif

2001

2002 doany = .false.

2003 do i=is1,ie1 ; if (do_i(i)) then

2004 if (.not.last_it(i)) doany = .true.

2005 if (last_it(i)) then

2006 if (need_bracket(i)) then

2007 if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then

2008 f_best(i) = f(i) ; ent_best(i) = ent(i)

2009 elseif (f_maxent(i) > f_minent(i)) then

2010 f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i)

2011 else

2012 f_best(i) = f_minent(i) ; ent_best(i) = minent(i)

2013 endif

2014 elseif (f(i) > f_best(i)) then

2015 f_best(i) = f(i) ; ent_best(i) = ent(i)

2016 endif

2017 do_i(i) = .false.

2018 elseif (need_bracket(i)) then

2019 if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then

2020 need_bracket(i) = .false. ! The maximum is now bracketed.

2021 chg_prev(i) = (maxent(i) - minent(i))

2022 chg_pre_prev(i) = 2.0*chg_prev(i)

2023 ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)

2024 elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then

2025 new_min_bound = .true. ! We have a new minimum bound.

2026 elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then

2027 new_min_bound = .false. ! We have a new maximum bound.

2028 else ! This case would bracket a minimum. Weird.

2029 ! Unless the derivative indicates that there is a maximum near the

2030 ! lower bound, try keeping the end with the larger value of F

2031 ! in a tie keep the minimum as the answer here will be compared

2032 ! with the maximum input value later.

2033 new_min_bound = .true.

2034 if (df_de_min(i) > 0.0 .or. (f_minent(i) >= f_maxent(i))) &

2035 new_min_bound = .false.

2036 endif

2037 if (need_bracket(i)) then ! Still not bracketed.

2038 if (new_min_bound) then

2039 minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)

2040 else

2041 maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)

2042 endif

2043 endif

2044 else ! The root was previously bracketed.

2045 if (f(i) >= f_best(i)) then ! There is a new maximum.

2046 if (ent(i) > ent_best(i)) then ! Replace minent with ent_prev.

2047 minent(i) = ent_best(i) ; f_minent(i) = f_best(i) ; df_de_min(i) = df_de_best(i)

2048 else ! Replace maxent with ent_best.

2049 maxent(i) = ent_best(i) ; f_maxent(i) = f_best(i) ; df_de_max(i) = df_de_best(i)

2050 endif

2051 ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)

2052 else

2053 if (ent(i) < ent_best(i)) then ! Replace the minent with ent.

2054 minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)

2055 else ! Replace maxent with ent_prev.

2056 maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)

2057 endif

2058 endif

2059 if ((maxent(i) - minent(i)) <= tolerance) do_i(i) = .false. ! Done.

2060 endif ! need_bracket.

2061 endif ; enddo

2062 if (.not.doany) exit

2063 enddo

2064 endif

2065

2066 if (present(f_lim_maxent)) then

2067 ! Return the unlimited maximum in maxF, and the limited value of F at maxent.

2068 do i=is,ie

2069 maxf(i) = f_best(i)

2070 f_lim_maxent(i) = f_max_ent_in(i)

2071 if (present(ent_maxf)) ent_maxf(i) = ent_best(i)

2072 enddo

2073 else

2074 ! Now compare the two? potential maxima using the limited value of dF_kb.

2075 doany = .false.

2076 do i=is,ie

2077 may_use_best(i) = (ent_best(i) /= max_ent_in(i))

2078 if (may_use_best(i)) doany = .true.

2079 enddo

2080 if (doany) then

2081 ! For efficiency, could save previous value of dS_anom_lim_best?

2082 call determine_dskb(h_bl, sref, ent_bl, ent_best, is, ie, kmb, g, gv, .true., ds_kb_lim)

2083 do i=is,ie

2084 f_best(i) = ent_best(i)*ds_kb_lim(i)*i_dskbp1(i)

2085 ! The second test seems necessary because of roundoff differences that

2086 ! can arise during compilation.

2087 if ((f_best(i) > f_max_ent_in(i)) .and. (may_use_best(i))) then

2088 maxf(i) = f_best(i)

2089 if (present(ent_maxf)) ent_maxf(i) = ent_best(i)

2090 else

2091 maxf(i) = f_max_ent_in(i)

2092 if (present(ent_maxf)) ent_maxf(i) = max_ent_in(i)

2093 endif

2094 enddo

2095 else

2096 ! All of the maxima are at the maximum entrainment.

2097 do i=is,ie ; maxf(i) = f_max_ent_in(i) ; enddo

2098 if (present(ent_maxf)) then

2099 do i=is,ie ; ent_maxf(i) = max_ent_in(i) ; enddo

2100 endif

2101 endif

2102 endif

2103

2104end subroutine find_maxf_kb

2105

2106!> This subroutine initializes the parameters and memory associated with the

2107!! entrain_diffusive module.

2108subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_read_params)

2109 type(time_type), intent(in) :: time !< The current model time.

2110 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

2111 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.

2112 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2113 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time

2114 !! parameters.

2115 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic

2116 !! output.

2117 type(entrain_diffusive_cs), intent(inout) :: cs !< Entrainment diffusion control structure

2118 logical, intent(in) :: just_read_params !< If true, this call will only read

2119 !! and log parameters without registering

2120 !! any diagnostics

2121

2122 ! Local variables

2123 real :: dt ! The dynamics timestep, used here in the default for TOLERANCE_ENT [T ~> s]

2124 real :: kd ! A diffusivity used in the default for TOLERANCE_ENT [Z2 T-1 ~> m2 s-1]

2125 ! This include declares and sets the variable "version".

2126# include "version_variable.h"

2127 character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name.

2128

2129 cs%initialized = .true.

2130 cs%diag => diag

2131

2132 cs%bulkmixedlayer = (gv%nkml > 0)

2133

2134 ! Set default, read and log parameters

2135 if (.not.just_read_params) call log_version(param_file, mdl, version, "")

2136 call get_param(param_file, mdl, "MAX_ENT_IT", cs%max_ent_it, &

2137 "The maximum number of iterations that may be used to "//&

2138 "calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read_params)

2139 ! In this module, KD is only used to set the default for TOLERANCE_ENT. [Z2 T-1 ~> m2 s-1]

2140 call get_param(param_file, mdl, "KD", kd, units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)

2141 call get_param(param_file, mdl, "DT", dt, &

2142 "The (baroclinic) dynamics time step.", &

2143 units="s", scale=us%s_to_T, fail_if_missing=.true., do_not_log=just_read_params)

2144 call get_param(param_file, mdl, "TOLERANCE_ENT", cs%Tolerance_Ent, &

2145 "The tolerance with which to solve for entrainment values.", &

2146 units="m", default=us%Z_to_m*max(100.0*gv%Angstrom_Z,1.0e-4*sqrt(dt*kd)), scale=gv%m_to_H, &

2147 do_not_log=just_read_params)

2148 call get_param(param_file, mdl, "ENTRAIN_DIFFUSIVE_MAX_ENT", cs%max_Ent, &

2149 "A large ceiling on the maximum permitted amount of entrainment across each "//&

2150 "interface between the mixed and buffer layers within a timestep.", &

2151 units="m", default=1.0e4, scale=gv%m_to_H, do_not_log=.not.cs%bulkmixedlayer)

2152

2153 cs%Rho_sig_off = 1000.0*us%kg_m3_to_R

2154

2155 if (.not.just_read_params) then

2156 cs%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, time, &

2157 'Diapycnal diffusivity as applied', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

2158 cs%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, time, &

2159 'Work actually done by diapycnal diffusion across each interface', &

2160 'W m-2', conversion=us%RZ3_T3_to_W_m2)

2161 endif

2162end subroutine entrain_diffusive_init

2163

2164!> \namespace mom_entrain_diffusive

2165!!

2166!! By Robert Hallberg, September 1997 - July 2000

2167!!

2168!! This file contains the subroutines that implement diapycnal

2169!! mixing and advection in isopycnal layers. The main subroutine,

2170!! calculate_entrainment, returns the entrainment by each layer

2171!! across the interfaces above and below it. These are calculated

2172!! subject to the constraints that no layers can be driven to negative

2173!! thickness and that the each layer maintains its target density,

2174!! using the scheme described in Hallberg (MWR 2000). There may or

2175!! may not be a bulk mixed layer above the isopycnal layers.

2176!! The solution is iterated until the change in the entrainment

2177!! between successive iterations is less than some small tolerance.

2178!!

2179!! The dual-stream entrainment scheme of MacDougall and Dewar

2180!! (JPO 1997) is used for combined diapycnal advection and diffusion,

2181!! modified as described in Hallberg (MWR 2000) to be solved

2182!! implicitly in time. Any profile of diffusivities may be used.

2183!! Diapycnal advection is fundamentally the residual of diapycnal

2184!! diffusion, so the fully implicit upwind differencing scheme that

2185!! is used is entirely appropriate. The downward buoyancy flux in

2186!! each layer is determined from an implicit calculation based on

2187!! the previously calculated flux of the layer above and an estimated

2188!! flux in the layer below. This flux is subject to the following

2189!! conditions: (1) the flux in the top and bottom layers are

2190!! set by the boundary conditions, and (2) no layer may be driven

2191!! below an Angstrom thickness. If there is a bulk mixed layer, the

2192!! mixed and buffer layers are treated as Eulerian layers, whose

2193!! thicknesses only change due to entrainment by the interior layers.

2194

2195end module mom_entrain_diffusive