17 implicit none ;
private
19 #include <MOM_memory.h>
30 logical :: bulkmixedlayer
32 logical :: correct_density
42 integer :: id_diff_work = -1
52 subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
53 kb_out, Kd_Lay, Kd_int)
57 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
62 type(
forcing),
intent(in) :: fluxes
64 real,
intent(in) :: dt
67 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
70 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
73 integer,
dimension(SZI_(G),SZJ_(G)), &
74 optional,
intent(inout) :: kb_out
77 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
78 optional,
intent(in) :: kd_lay
80 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
81 optional,
intent(in) :: kd_int
90 real,
dimension(SZI_(G),SZK_(G)) :: &
92 real,
dimension(SZI_(G),SZK_(G)+1) :: &
94 real,
dimension(SZI_(G),SZK_(G)) :: &
107 real,
dimension(SZI_(G),SZK_(G)+1) :: &
110 real,
allocatable,
dimension(:,:,:) :: &
117 real :: hm, fm, fr, fk
120 real :: c1(szi_(g),szk_(g))
122 real,
dimension(SZI_(G)) :: &
153 real,
dimension(SZI_(G),SZK_(G)) :: &
161 real,
dimension(SZI_(G),SZK_(G)) :: &
176 real,
dimension(SZI_(G)) :: &
200 logical :: do_entrain_eakb
201 logical :: do_i(szi_(g)), did_i(szi_(g)), reiterate, correct_density
202 integer :: it, i, j, k, is, ie, js, je, nz, k2, kmb
203 integer :: kb(szi_(g))
205 integer :: kb_min_act
207 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
208 angstrom = gv%Angstrom_H
209 h_neglect = gv%H_subroundoff
211 if (.not.
associated(cs))
call mom_error(fatal, &
212 "MOM_entrain_diffusive: Module must be initialized before it is used.")
214 if (.not.(
present(kd_lay) .or.
present(kd_int)))
call mom_error(fatal, &
215 "MOM_entrain_diffusive: Either Kd_Lay or Kd_int must be present in call.")
217 if ((.not.cs%bulkmixedlayer .and. .not.
associated(fluxes%buoy)) .and. &
218 (
associated(fluxes%lprec) .or.
associated(fluxes%evap) .or. &
219 associated(fluxes%sens) .or.
associated(fluxes%sw)))
then
221 &The code to handle evaporation and precipitation without &
222 &a bulk mixed layer has not been implemented.")
224 "Either define BULKMIXEDLAYER in MOM_input or use fluxes%buoy &
225 &and a linear equation of state to drive the model.")
228 tolerance = cs%Tolerance_Ent
229 kmb = gv%nk_rho_varies
230 k2 = max(kmb+1,2) ; kb_min = k2
231 if (.not. cs%bulkmixedlayer)
then
235 do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ;
enddo
238 do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ;
enddo
241 if (cs%id_diff_work > 0)
allocate(diff_work(g%isd:g%ied,g%jsd:g%jed,nz+1))
242 if (cs%id_Kd > 0)
allocate(kd_eff(g%isd:g%ied,g%jsd:g%jed,nz))
244 correct_density = (cs%correct_density .and.
associated(tv%eqn_of_state))
245 if (correct_density)
then
268 do i=is,ie ; kb(i) = 1 ;
enddo
270 if (
present(kd_lay))
then
271 do k=1,nz ;
do i=is,ie
272 dtkd(i,k) = gv%Z_to_H**2 * (dt * kd_lay(i,j,k))
274 if (
present(kd_int))
then
275 do k=1,nz+1 ;
do i=is,ie
276 dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
279 do k=2,nz ;
do i=is,ie
280 dtkd_int(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_lay(i,j,k-1) + kd_lay(i,j,k)))
284 do k=1,nz ;
do i=is,ie
285 dtkd(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_int(i,j,k)+kd_int(i,j,k+1)))
287 dO k=1,nz+1 ;
do i=is,ie
288 dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
292 do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ;
enddo
293 do i=is,ie ; ds_dsp1(i,nz) = 0.0 ;
enddo
294 do i=is,ie ; dsp1_ds(i,nz) = 0.0 ;
enddo
296 do k=2,nz-1 ;
do i=is,ie
297 ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
300 if (cs%bulkmixedlayer)
then
304 call set_ent_bl(h, dtkd_int, tv, kb, kmb, do_i, g, gv, us, cs, j, ent_bl, sref, h_bl)
307 dtkd_kb(i) = 0.0 ;
if (kb(i) < nz) dtkd_kb(i) = dtkd(i,kb(i))
310 do i=is,ie ; ent_bl(i,kmb+1) = 0.0 ;
enddo
313 do k=2,nz-1 ;
do i=is,ie
314 dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
315 i2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))
316 grats(i,k) = 2.0*(2.0+(dsp1_ds(i,k)+ds_dsp1(i,k)))
322 if (cs%bulkmixedlayer)
then
325 htot(i) = h(i,j,1) - angstrom
327 do k=2,kmb ;
do i=is,ie
328 htot(i) = htot(i) + (h(i,j,k) - angstrom)
331 max_eakb(i) = max(ent_bl(i,kmb+1) + 0.5*htot(i), htot(i))
332 i_dskbp1(i) = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
340 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, max_eakb, kmb, &
341 is, ie, g, gv, cs, maxf_kb, eakb_maxf, do_i, f_kb_maxent)
342 do i=is,ie ;
if (kb(i) <= nz)
then
343 maxf(i,kb(i)) = max(0.0, maxf_kb(i), f_kb_maxent(i))
344 if ((maxf_kb(i) > f_kb_maxent(i)) .and. (eakb_maxf(i) >= htot(i))) &
345 max_eakb(i) = eakb_maxf(i)
348 do i=is,ie ; ea_kbp1(i) = 0.0 ; eakb(i) = max_eakb(i) ;
enddo
350 max_eakb, max_eakb, kmb, is, ie, do_i, g, gv, cs, eakb, &
351 error=err_max_eakb0, f_kb=f_kb)
355 do i=is,ie ; min_eakb(i) = min(htot(i), max_eakb(i)) ;
enddo
356 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, min_eakb, max_eakb, &
357 kmb, is, ie, g, gv, cs, f_kb_maxent, do_i_in = do_i)
360 do_entrain_eakb = .false.
362 if (do_i(i))
then ;
if (err_max_eakb0(i) < 0.0)
then
363 do_entrain_eakb = .true.
366 if (do_entrain_eakb)
then
367 eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
369 eakb(i) = 0.0 ; min_eakb(i) = 0.0
377 zeros, max_eakb, kmb, is, ie, do_i, g, gv, cs, min_eakb, &
378 error=err_min_eakb0, f_kb=f_kb, err_max_eakb0=err_max_eakb0)
380 do i=is,ie ;
if ((kb(i)<nz) .and. (kb_min>kb(i))) kb_min = kb(i) ;
enddo
389 htot(i) = h(i,j,1) - angstrom
391 if (
associated(fluxes%buoy))
then ;
do i=is,ie
392 maxf(i,1) = gv%Z_to_H * (dt*fluxes%buoy(i,j)) / gv%g_prime(2)
398 do k=kb_min,nz-1 ;
do i=is,ie
399 if ((k == kb(i)+1) .and. cs%bulkmixedlayer)
then
400 maxf(i,k) = ds_dsp1(i,k)*(f_kb_maxent(i) + htot(i))
401 htot(i) = htot(i) + (h(i,j,k) - angstrom)
402 elseif (k > kb(i))
then
403 maxf(i,k) = ds_dsp1(i,k)*(maxf(i,k-1) + htot(i))
404 htot(i) = htot(i) + (h(i,j,k) - angstrom)
409 if (.not.cs%bulkmixedlayer)
then
410 maxf_correct(i) = max(0.0, -(maxf(i,nz-1) + htot(i)))
412 htot(i) = h(i,j,nz) - angstrom
414 if (.not.cs%bulkmixedlayer)
then
415 do_any = .false. ;
do i=is,ie ;
if (maxf_correct(i) > 0.0) do_any = .true. ;
enddo
417 do k=nz-1,1,-1 ;
do i=is,ie
418 maxf(i,k) = maxf(i,k) + maxf_correct(i)
419 maxf_correct(i) = maxf_correct(i) * dsp1_ds(i,k)
423 do k=nz-1,kb_min,-1 ;
do i=is,ie ;
if (do_i(i))
then
425 maxf(i,k) = min(maxf(i,k),dsp1_ds(i,k+1)*maxf(i,k+1) + htot(i))
426 htot(i) = htot(i) + (h(i,j,k) - angstrom)
429 if ((maxf(i,k) < f_kb(i)) .or. (maxf(i,k) < maxf_kb(i)) &
430 .and. (eakb_maxf(i) <= max_eakb(i)))
then
435 if ((f_kb(i) <= maxf_kb(i)) .and. (eakb_maxf(i) <= max_eakb(i)))
then
436 eakb(i) = eakb_maxf(i)
438 eakb(i) = max_eakb(i)
440 call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
441 g, gv, cs, eakb, angstrom)
442 if ((eakb(i) < max_eakb(i)) .or. (eakb(i) < min_eakb(i)))
then
444 eakb, eakb, kmb, i, i, do_i, g, gv, cs, eakb, &
446 if (eakb(i) < max_eakb(i))
then
447 max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)
449 if (eakb(i) < min_eakb(i))
then
450 min_eakb(i) = eakb(i) ; err_min_eakb0(i) = err_eakb0(i)
455 endif ;
enddo ;
enddo
456 if (.not.cs%bulkmixedlayer)
then
458 maxf(i,1) = min(maxf(i,1),dsp1_ds(i,2)*maxf(i,2) + htot(i))
469 f(i,nz) = maxf(i,nz) ; minf(i,nz) = 0.0
473 if ((k==kb(i)) .and. (do_i(i)))
then
474 eakb(i) = min_eakb(i)
476 elseif ((k>kb(i)) .and. (do_i(i)))
then
481 hm = h(i,j,k) + h_neglect
484 f(i,k) = min(maxf(i,k), sqrt(ds_dsp1(i,k)*dtkd(i,k)), &
485 0.5*(ds_dsp1(i,k)+1.0) * (dtkd(i,k) / hm))
490 fk = dtkd(i,k) * grats(i,k)
491 minf(i,k) = min(maxf(i,k), &
492 0.9*(i2p2dsp1_ds(i,k) * fk / (hm + sqrt(hm*hm + fk))))
493 if (k==kb(i)) minf(i,k) = 0.0
503 is1 = ie+1 ; ie1 = is-1
504 do i=is,ie ;
if (do_i(i))
then ; is1 = i ;
exit ;
endif ;
enddo
505 do i=ie,is,-1 ;
if (do_i(i))
then ; ie1 = i ;
exit ;
endif ;
enddo
507 if (cs%bulkmixedlayer)
then
510 if (do_i(i) .and. (kb(i) < kb_min_act)) kb_min_act = kb(i)
517 if (do_i(i) .and. (kb(i) < nz)) &
518 ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*f(i,kb(i)+1)
520 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
521 max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
522 err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
528 do it=0,cs%max_ent_it-1
529 do i=is1,ie1 ;
if (do_i(i))
then
530 if (.not.cs%bulkmixedlayer) f(i,1) = min(f(i,1),maxf(i,1))
536 do k=kb_min_act,nz-1 ;
do i=is1,ie1 ;
if (do_i(i) .and. (k>=kb(i)))
then
538 if (cs%bulkmixedlayer .and. (k==kb(i)))
then
540 dfdfm(i,k) = dfdfm_kb(i)
543 fm = (f(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*f(i,k+1)
544 fk = grats(i,k)*dtkd(i,k)
545 fr = sqrt(fm*fm + fk)
548 f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fm+fr))
550 f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))
553 if ((f(i,k) >= maxf(i,k)) .or. (fr == 0.0))
then
556 dfdfm(i,k) = i2p2dsp1_ds(i,k) * ((fr + fm) / fr)
561 c1(i,k) = dfdfm(i,k-1)*(dsp1_ds(i,k)*b1(i))
562 b1(i) = 1.0 / (1.0 - c1(i,k)*dfdfm(i,k))
563 f(i,k) = min(b1(i)*(f(i,k)-fprev(i,k)) + fprev(i,k), maxf(i,k))
564 if (f(i,k) >= maxf(i,k)) dfdfm(i,k) = 0.0
567 endif ;
enddo ;
enddo
569 do k=nz-2,kb_min_act,-1 ;
do i=is1,ie1
570 if (do_i(i) .and. (k > kb(i))) &
571 f(i,k) = min((f(i,k)+c1(i,k+1)*(f(i,k+1)-fprev(i,k+1))),maxf(i,k))
574 if (cs%bulkmixedlayer)
then
576 if (do_i(i) .and. (kb(i) < nz))
then
578 ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(f(i,kb(i)+1), minf(i,kb(i)+1))
583 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
584 max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
585 err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
588 if (do_i(i) .and. (kb(i) < nz)) f(i,kb(i)) = f_kb(i)
593 if (it < cs%max_ent_it-1)
then
596 if (cs%bulkmixedlayer)
then ;
do i=is1,ie1 ;
if (do_i(i))
then
597 eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
598 endif ;
enddo ;
endif
600 did_i(i) = do_i(i) ; do_i(i) = .false.
602 do k=kb_min_act,nz-1 ;
do i=is1,ie1
603 if (did_i(i) .and. (k >= kb(i)))
then
604 if (f(i,k) < minf(i,k))
then
606 do_i(i) = .true. ; reiterate = .true.
607 elseif (k > kb(i))
then
608 if ((abs(f(i,k) - fprev(i,k)) > tolerance) .or. &
609 ((h(i,j,k) + ((1.0+dsp1_ds(i,k))*f(i,k) - &
610 (f(i,k-1) + dsp1_ds(i,k+1)*f(i,k+1)))) < 0.9*angstrom))
then
611 do_i(i) = .true. ; reiterate = .true.
617 if (h(i,j,k) + ((f(i,k) + eakb(i)) - &
618 (eb_kmb(i) + dsp1_ds(i,k+1)*f(i,k+1))) < 0.9*angstrom)
then
619 do_i(i) = .true. ; reiterate = .true.
624 if (.not.reiterate)
exit
630 if (it == (cs%max_ent_it))
then
633 do i=is1,ie1 ;
if (do_i(i))
then
634 f(i,nz-1) = max(f(i,nz-1), min(minf(i,nz-1), 0.0))
635 if (kb(i) >= nz-1)
then ; ea_kbp1(i) = 0.0 ;
endif
637 do k=nz-2,kb_min_act,-1 ;
do i=is1,ie1 ;
if (do_i(i))
then
639 f(i,k) = min(max(minf(i,k),f(i,k)), (dsp1_ds(i,k+1)*f(i,k+1) + &
640 max(((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + &
641 (h(i,j,k+1) - angstrom)), 0.5*(h(i,j,k+1)-angstrom))))
642 elseif (k==kb(i))
then
643 ea_kbp1(i) = dsp1_ds(i,k+1)*f(i,k+1)
644 h_avail = dsp1_ds(i,k+1)*f(i,k+1) + max(0.5*(h(i,j,k+1)-angstrom), &
645 ((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + (h(i,j,k+1) - angstrom)))
646 if ((f(i,k) > 0.0) .and. (f(i,k) > h_avail))
then
647 f_kb(i) = max(0.0, h_avail)
649 if ((f_kb(i) < maxf_kb(i)) .and. (eakb_maxf(i) <= eakb(i))) &
650 eakb(i) = eakb_maxf(i)
651 call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
655 endif ;
enddo ;
enddo
658 if (cs%bulkmixedlayer)
then ;
do i=is1,ie1
659 if (do_i(i) .and. (kb(i) < nz))
then
660 h_avail = eakb(i) + max(0.5*(h_bl(i,kmb+1)-angstrom), &
661 (f_kb(i)-ea_kbp1(i)) + (h_bl(i,kmb+1)-angstrom))
663 ent_bl(i,kmb+1) = min(ent_bl(i,kmb+1),0.5*(eakb(i) + h_avail))
665 eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
670 do k=kb_min_act+1,nz-1 ;
do i=is1,ie1 ;
if (do_i(i))
then
671 if ((.not.cs%bulkmixedlayer) .or. (k > kb(i)+1))
then
672 f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + &
673 dsp1_ds(i,k-1)*f(i,k-1)) - f(i,k-2)) + (h(i,j,k-1) - angstrom)))
674 f(i,k) = max(f(i,k),min(minf(i,k),0.0))
675 elseif (k == kb(i)+1)
then
676 f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + eakb(i)) - &
677 eb_kmb(i)) + (h(i,j,k-1) - angstrom)))
678 f(i,k) = max(f(i,k),min(minf(i,k),0.0))
680 endif ;
enddo ;
enddo
683 call f_to_ent(f, h, kb, kmb, j, g, gv, cs, dsp1_ds, eakb, ent_bl, ea, eb)
687 if (correct_density)
then
689 h_guess(i,1) = (h(i,j,1) - angstrom) + (eb(i,j,1) - ea(i,j,2))
690 h_guess(i,nz) = (h(i,j,nz) - angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))
691 if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0
692 if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0
694 do k=2,nz-1 ;
do i=is,ie
695 h_guess(i,k) = (h(i,j,k) - angstrom) + ((ea(i,j,k) - eb(i,j,k-1)) + &
696 (eb(i,j,k) - ea(i,j,k+1)))
697 if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0
699 if (cs%bulkmixedlayer)
then
700 call determine_dskb(h_bl, sref, ent_bl, eakb, is, ie, kmb, g, gv, &
701 .true., ds_kb, ds_anom_lim=ds_anom_lim)
704 rcv(is:ie), 1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
706 if ((k>kb(i)) .and. (f(i,k) > 0.0))
then
713 f_cor = h(i,j,k) * min(1.0 , max(-ds_dsp1(i,k), &
714 (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )
720 if (f_cor > 0.0)
then
721 f_cor = min(f_cor, 0.9*f(i,k), ds_dsp1(i,k)*0.5*h_guess(i,k), &
724 f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
725 0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
728 ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
729 eb(i,j,k) = eb(i,j,k) + f_cor
730 elseif ((k==kb(i)) .and. (f(i,k) > 0.0))
then
734 ds_kb_eff = 2.0*ds_kb(i) - ds_anom_lim(i)
735 rho_cor = h(i,j,k) * (gv%Rlay(k)-rcv(i)) + eakb(i)*ds_anom_lim(i)
738 if (abs(rho_cor) < abs(0.9*eakb(i)*ds_kb_eff))
then
739 ea_cor = -rho_cor / ds_kb_eff
741 ea_cor = sign(0.9*eakb(i),-rho_cor*ds_kb_eff)
744 if (ea_cor > 0.0)
then
746 ea_cor = min(ea_cor, 0.5*(max_eakb(i) - eakb(i)), &
747 0.5*h_guess(i,k) / (ds_kb(i) * i_dskbp1(i)))
750 ea_cor = -min(-ea_cor, 0.5*h_guess(i,k), &
751 0.25*h_guess(i,k+1) / (ds_kb(i) * i_dskbp1(i)))
754 ea(i,j,k) = ea(i,j,k) + ea_cor
755 eb(i,j,k) = eb(i,j,k) - (ds_kb(i) * i_dskbp1(i)) * ea_cor
756 elseif (k < kb(i))
then
758 ea(i,j,k) = ea(i,j,k+1)
762 do k=kb_min-1,k2,-1 ;
do i=is,ie
763 ea(i,j,k) = ea(i,j,k+1)
771 h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
772 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
774 do k=kmb-1,2,-1 ;
do i=is,ie
776 eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)
779 h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
780 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
783 eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
789 rcv(is:ie), 1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
790 do i=is,ie ;
if (f(i,k) > 0.0)
then
795 f_cor = h(i,j,k) * min(dsp1_ds(i,k) , max(-1.0, &
796 (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )
802 if (f_cor >= 0.0)
then
803 f_cor = min(f_cor, 0.9*f(i,k), 0.5*dsp1_ds(i,k)*h_guess(i,k), &
806 f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
807 0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
810 ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
811 eb(i,j,k) = eb(i,j,k) + f_cor
818 if (cs%id_Kd > 0)
then
819 idt = gv%H_to_Z**2 / dt
820 do k=2,nz-1 ;
do i=is,ie
821 if (k<kb(i))
then ; kd_here = 0.0 ;
else
822 kd_here = f(i,k) * ( h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
823 (eb(i,j,k) - ea(i,j,k+1))) ) / (i2p2dsp1_ds(i,k) * grats(i,k))
826 kd_eff(i,j,k) = max(dtkd(i,k), kd_here)*idt
829 kd_eff(i,j,1) = dtkd(i,1)*idt
830 kd_eff(i,j,nz) = dtkd(i,nz)*idt
834 if (cs%id_diff_work > 0)
then
835 g_2dt = 0.5 * gv%H_to_Z**2*us%L_to_Z**2 * (gv%g_Earth / dt)
836 do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ;
enddo
837 if (
associated(tv%eqn_of_state))
then
838 if (
associated(fluxes%p_surf))
then
839 do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ;
enddo
841 do i=is,ie ; pressure(i) = 0.0 ;
enddo
844 do i=is,ie ; pressure(i) = pressure(i) + gv%H_to_Pa*h(i,j,k-1) ;
enddo
847 t_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))
848 s_eos(i) = 0.5*(tv%S(i,j,kmb) + tv%S(i,j,k))
850 t_eos(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
851 s_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
855 drho_dt, drho_ds, is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
857 if ((k>kmb) .and. (k<kb(i)))
then ; diff_work(i,j,k) = 0.0
860 drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &
861 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))
863 drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
864 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))
866 diff_work(i,j,k) = g_2dt * drho * &
867 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
868 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
873 do k=2,nz ;
do i=is,ie
874 diff_work(i,j,k) = g_2dt * (gv%Rlay(k)-gv%Rlay(k-1)) * &
875 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
876 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
881 if (
present(kb_out))
then
882 do i=is,ie ; kb_out(i,j) = kb(i) ;
enddo
888 if (cs%id_Kd > 0)
call post_data(cs%id_Kd, kd_eff, cs%diag)
889 if (cs%id_Kd > 0)
deallocate(kd_eff)
890 if (cs%id_diff_work > 0)
call post_data(cs%id_diff_work, diff_work, cs%diag)
891 if (cs%id_diff_work > 0)
deallocate(diff_work)
898 subroutine f_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, do_i_in)
901 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: F
905 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
907 integer,
dimension(SZI_(G)),
intent(in) :: kb
909 integer,
intent(in) :: kmb
910 integer,
intent(in) :: j
912 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: dsp1_ds
916 real,
dimension(SZI_(G)),
intent(in) :: eakb
918 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: Ent_bl
921 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
924 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
927 logical,
dimension(SZI_(G)), &
928 optional,
intent(in) :: do_i_in
935 logical :: do_i(SZI_(G))
936 integer :: i, k, is, ie, nz
938 is = g%isc ; ie = g%iec ; nz = g%ke
940 if (
present(do_i_in))
then
941 do i=is,ie ; do_i(i) = do_i_in(i) ;
enddo
942 do i=g%isc,g%iec ;
if (do_i(i))
then
945 do i=g%iec,g%isc,-1 ;
if (do_i(i))
then
949 do i=is,ie ; do_i(i) = .true. ;
enddo
953 ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0
955 if (cs%bulkmixedlayer)
then
957 eb(i,j,kmb) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
959 do k=nz-1,kmb+1,-1 ;
do i=is,ie ;
if (do_i(i))
then
963 ea(i,j,k) = dsp1_ds(i,k)*f(i,k)
965 elseif (k == kb(i))
then
968 elseif (k == kb(i)-1)
then
969 ea(i,j,k) = ea(i,j,k+1)
970 eb(i,j,k) = eb(i,j,kmb)
972 ea(i,j,k) = ea(i,j,k+1)
975 eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)
977 endif ;
enddo ;
enddo
979 do i=is,ie ;
if (do_i(i))
then
983 eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)
986 h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
987 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
989 do k=kmb-1,2,-1 ;
do i=is,ie ;
if (do_i(i))
then
991 eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)
994 h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
995 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
999 endif ;
enddo ;
enddo
1000 do i=is,ie ;
if (do_i(i))
then
1001 eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
1009 ea(i,j,1) = 0.0 ; eb(i,j,1) = max(f(i,1),0.0)
1010 ea(i,j,2) = dsp1_ds(i,2)*f(i,2) - min(f(i,1),0.0)
1013 do k=2,nz-1 ;
do i=is,ie
1014 eb(i,j,k) = max(f(i,k),0.0)
1015 ea(i,j,k+1) = dsp1_ds(i,k+1)*f(i,k+1) - (f(i,k)-eb(i,j,k))
1016 if (ea(i,j,k+1) < 0.0)
then
1017 eb(i,j,k) = eb(i,j,k) - ea(i,j,k+1)
1028 subroutine set_ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)
1031 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1033 real,
dimension(SZI_(G),SZK_(G)+1), &
1034 intent(in) :: dtKd_int
1040 integer,
dimension(SZI_(G)),
intent(inout) :: kb
1043 integer,
intent(in) :: kmb
1044 logical,
dimension(SZI_(G)),
intent(in) :: do_i
1048 integer,
intent(in) :: j
1049 real,
dimension(SZI_(G),SZK_(G)+1), &
1050 intent(out) :: Ent_bl
1053 real,
dimension(SZI_(G),SZK_(G)),
intent(out) :: Sref
1055 real,
dimension(SZI_(G),SZK_(G)),
intent(out) :: h_bl
1064 real,
dimension(SZI_(G)) :: &
1071 real,
dimension(SZI_(G), SZK_(G)) :: &
1080 integer :: i, k, is, ie, nz
1081 is = g%isc ; ie = g%iec ; nz = g%ke
1084 max_ent = 1.0e4*gv%m_to_H
1085 h_neglect = gv%H_subroundoff
1087 do i=is,ie ; pres(i) = tv%P_Ref ;
enddo
1090 rcv(is:ie), 1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1092 h_bl(i,k) = h(i,j,k) + h_neglect
1093 sref(i,k) = rcv(i) - cs%Rho_sig_off
1098 h_interior(i) = 0.0 ; ent_bl(i,1) = 0.0
1102 do k=2,kmb ;
do i=is,ie
1104 ent_bl(i,k) = min(2.0 * dtkd_int(i,k) / (h(i,j,k-1) + h(i,j,k) + h_neglect), &
1106 else ; ent_bl(i,k) = 0.0 ;
endif
1115 b1(i) = 1.0 / (h_bl(i,1) + ent_bl(i,2))
1116 d1(i) = h_bl(i,1) * b1(i)
1117 s_est(i,1) = (h_bl(i,1)*sref(i,1)) * b1(i)
1119 do k=2,kmb-1 ;
do i=is,ie
1120 b1(i) = 1.0 / ((h_bl(i,k) + ent_bl(i,k+1)) + d1(i)*ent_bl(i,k))
1121 d1(i) = (h_bl(i,k) + d1(i)*ent_bl(i,k)) * b1(i)
1122 s_est(i,k) = (h_bl(i,k)*sref(i,k) + ent_bl(i,k)*s_est(i,k-1)) * b1(i)
1125 s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
1126 (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
1132 do i=is,ie ; kb(i) = nz+1 ;
if (do_i(i)) kb(i) = kmb+1 ;
enddo
1134 do k=kmb+1,nz ;
do i=is,ie ;
if (do_i(i))
then
1135 if ((k == kb(i)) .and. (s_est(i,kmb) > (gv%Rlay(k) - cs%Rho_sig_off)))
then
1136 if (4.0*dtkd_int(i,kmb+1)*frac_rem(i) > &
1137 (h_bl(i,kmb) + h(i,j,k)) * (h(i,j,k) - gv%Angstrom_H))
then
1139 dh = max((h(i,j,k) - gv%Angstrom_H), 0.0)
1141 frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &
1142 (4.0*dtkd_int(i,kmb+1))
1143 sref(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + dh*(gv%Rlay(k)-cs%Rho_sig_off)) / &
1145 h_bl(i,kmb) = h_bl(i,kmb) + dh
1146 s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
1147 (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
1152 endif ;
enddo ;
enddo
1156 do k=nz,kmb+1,-1 ;
do i=is,ie
1157 if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-gv%Angstrom_H)
1159 h_bl(i,kmb+1) = h(i,j,k) ; sref(i,kmb+1) = gv%Rlay(k) - cs%Rho_sig_off
1160 elseif (k==kb(i)+1)
then
1161 h_bl(i,kmb+2) = h(i,j,k) ; sref(i,kmb+2) = gv%Rlay(k) - cs%Rho_sig_off
1164 do i=is,ie ;
if (kb(i) >= nz)
then
1165 h_bl(i,kmb+1) = h(i,j,nz)
1166 sref(i,kmb+1) = gv%Rlay(nz) - cs%Rho_sig_off
1167 h_bl(i,kmb+2) = gv%Angstrom_H
1168 sref(i,kmb+2) = sref(i,kmb+1) + (gv%Rlay(nz) - gv%Rlay(nz-1))
1174 do i=is,ie ;
if (do_i(i))
then
1175 kd_x_dt = frac_rem(i) * dtkd_int(i,kmb+1)
1176 if ((kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0))
then
1177 ent_bl(i,kmb+1) = 0.0
1181 ent_bl(i,kmb+1) = min(0.5*h_interior(i), sqrt(kd_x_dt), &
1182 kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))
1185 ent_bl(i,kmb+1) = 0.0
1202 subroutine determine_dskb(h_bl, Sref, Ent_bl, E_kb, is, ie, kmb, G, GV, limit, &
1203 dSkb, ddSkb_dE, dSlay, ddSlay_dE, dS_anom_lim, do_i_in)
1207 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: h_bl
1208 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: Sref
1209 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: Ent_bl
1212 real,
dimension(SZI_(G)),
intent(in) :: E_kb
1214 integer,
intent(in) :: is
1215 integer,
intent(in) :: ie
1216 integer,
intent(in) :: kmb
1217 logical,
intent(in) :: limit
1219 real,
dimension(SZI_(G)),
intent(inout) :: dSkb
1224 real,
dimension(SZI_(G)),
optional,
intent(inout) :: ddSkb_dE
1226 real,
dimension(SZI_(G)),
optional,
intent(inout) :: dSlay
1229 real,
dimension(SZI_(G)),
optional,
intent(inout) :: ddSlay_dE
1231 real,
dimension(SZI_(G)),
optional,
intent(inout) :: dS_anom_lim
1234 logical,
dimension(SZI_(G)),
optional,
intent(in) :: do_i_in
1254 real,
dimension(SZI_(G),SZK_(G)) :: &
1259 real :: deriv_dSkb(SZI_(G))
1264 logical,
dimension(SZI_(G)) :: do_i
1271 real :: dS_kbp1, IdS_kbp1
1274 real :: f1, df1_drat
1275 real :: z, dz_drat, f2, df2_dz, expz
1276 real :: eps_dSLay, eps_dSkb
1279 if (
present(ddslay_de) .and. .not.
present(dslay))
call mom_error(fatal, &
1280 "In deterimine_dSkb, ddSLay_dE may only be present if dSlay is.")
1282 h_neglect = gv%H_subroundoff
1285 ea(i,kmb+1) = e_kb(i) ; dea_de(i,kmb+1) = 1.0
1286 s(i,kmb+1) = sref(i,kmb+1) ; ds_de(i,kmb+1) = 0.0
1291 if (
present(do_i_in))
then
1292 do i=is,ie ; do_i(i) = do_i_in(i) ;
enddo
1294 do k=kmb,1,-1 ;
do i=is,ie
1298 if (2.0*ent_bl(i,k+1) > ea(i,k+1))
then
1299 eb(i,k) = 2.0*ent_bl(i,k+1) - ea(i,k+1) ; deb_de(i,k) = -dea_de(i,k+1)
1301 eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
1305 h1 = (h_bl(i,k) - gv%Angstrom_H) + (eb(i,k) - ea(i,k+1))
1307 ea(i,k) = ent_bl(i,k) ; dea_de(i,k) = 0.0
1308 elseif (ent_bl(i,k) + 0.5*h1 >= 0.0)
then
1309 ea(i,k) = ent_bl(i,k) - 0.5*h1
1310 dea_de(i,k) = 0.5*(dea_de(i,k+1) - deb_de(i,k))
1313 dea_de(i,k) = dea_de(i,k+1) - deb_de(i,k)
1316 ea(i,k) = 0.0 ; dea_de(i,k) = 0.0 ; eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
1320 h_tr = h_bl(i,k) + h_neglect
1321 c1(i,k) = ea(i,k+1) * b1(i,k+1)
1322 b_denom_1 = (h_tr + d1(i)*eb(i,k))
1323 b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))
1324 d1(i) = b_denom_1 * b1(i,k)
1326 s(i,k) = (h_tr*sref(i,k) + eb(i,k)*s(i,k+1)) * b1(i,k)
1328 do k=2,kmb ;
do i=is,ie
1329 s(i,k) = s(i,k) + c1(i,k-1)*s(i,k-1)
1332 if (
present(ddskb_de) .or.
present(ddslay_de))
then
1335 do k=kmb,2,-1 ;
do i=is,ie
1336 if (do_i(i) .and. (dea_de(i,k) - deb_de(i,k) > 0.0))
then
1337 src = (((s(i,k+1) - sref(i,k)) * (h_bl(i,k) + h_neglect) + &
1338 (s(i,k+1) - s(i,k-1)) * ea(i,k)) * deb_de(i,k) - &
1339 ((sref(i,k) - s(i,k-1)) * h_bl(i,k) + &
1340 (s(i,k+1) - s(i,k-1)) * eb(i,k)) * dea_de(i,k)) / &
1341 ((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))
1342 else ; src = 0.0 ;
endif
1343 ds_de(i,k) = (src + eb(i,k)*ds_de(i,k+1)) * b1(i,k)
1346 if (do_i(i) .and. (deb_de(i,1) < 0.0))
then
1347 src = (((s(i,2) - sref(i,1)) * (h_bl(i,1) + h_neglect)) * deb_de(i,1)) / &
1348 (h_bl(i,1) + h_neglect + eb(i,1))
1349 else ; src = 0.0 ;
endif
1350 ds_de(i,1) = (src + eb(i,1)*ds_de(i,2)) * b1(i,1)
1352 do k=2,kmb ;
do i=is,ie
1353 ds_de(i,k) = ds_de(i,k) + c1(i,k-1)*ds_de(i,k-1)
1360 if (.not.limit)
then
1361 do i=is,ie ;
if (do_i(i))
then
1362 dskb(i) = sref(i,kmb+1) - s(i,kmb)
1364 if (
present(ddskb_de))
then ;
do i=is,ie ;
if (do_i(i))
then
1365 ddskb_de(i) = -1.0*ds_de(i,kmb)
1366 endif ;
enddo ;
endif
1368 if (
present(dslay))
then ;
do i=is,ie ;
if (do_i(i))
then
1369 dslay(i) = 0.5 * (sref(i,kmb+2) - s(i,kmb))
1370 endif ;
enddo ;
endif
1371 if (
present(ddslay_de))
then ;
do i=is,ie ;
if (do_i(i))
then
1372 ddslay_de(i) = -0.5*ds_de(i,kmb)
1373 endif ;
enddo ;
endif
1375 do i=is,ie ;
if (do_i(i))
then
1377 if (sref(i,kmb+1) - s(i,kmb) < eps_dskb*(sref(i,kmb+2) - sref(i,kmb+1)))
then
1378 dskb(i) = eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) ; deriv_dskb(i) = 0.0
1380 dskb(i) = sref(i,kmb+1) - s(i,kmb) ; deriv_dskb(i) = -1.0
1382 if (
present(ddskb_de)) ddskb_de(i) = deriv_dskb(i)*ds_de(i,kmb)
1385 if (
present(dslay))
then
1389 do i=is,ie ;
if (do_i(i))
then
1390 ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
1391 ids_kbp1 = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
1392 rat = (sref(i,kmb+1) - s(i,kmb)) * ids_kbp1
1400 inv_term = 1.0 / (1.0-rat)
1401 f1 = 2.0 - 0.125*(inv_term**2)
1402 df1_drat = - 0.25*(inv_term**3)
1406 z = dz_drat * rat + 4.0
1407 if (z >= 18.0)
then ; f2 = 1.0 ; df2_dz = 0.0
1408 elseif (z <= -58.0)
then ; f2 = eps_dslay ; df2_dz = 0.0
1410 expz = exp(z) ; inv_term = 1.0 / (1.0 + expz)
1411 f2 = (eps_dslay + expz) * inv_term
1412 df2_dz = (1.0 - eps_dslay) * expz * inv_term**2
1415 dslay(i) = dskb(i) * f1 * f2
1416 deriv_dslay = deriv_dskb(i) * (f1 * f2) - (dskb(i)*ids_kbp1) * &
1417 (df1_drat*f2 + f1 * dz_drat * df2_dz)
1418 elseif (dskb(i) <= 3.0*ds_kbp1)
then
1420 dslay(i) = 0.5 * (dskb(i) + ds_kbp1)
1421 deriv_dslay = 0.5 * deriv_dskb(i)
1423 dslay(i) = 2.0*ds_kbp1
1426 if (
present(ddslay_de)) ddslay_de(i) = deriv_dslay*ds_de(i,kmb)
1431 if (
present(ds_anom_lim))
then ;
do i=is,ie ;
if (do_i(i))
then
1432 ds_anom_lim(i) = max(0.0, eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) - &
1433 (sref(i,kmb+1) - s(i,kmb)) )
1434 endif ;
enddo ;
endif
1442 subroutine f_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &
1443 G, GV, CS, ea_kb, tol_in)
1446 real,
dimension(SZI_(G),SZK_(G)), &
1449 real,
dimension(SZI_(G),SZK_(G)), &
1453 real,
dimension(SZI_(G),SZK_(G)), &
1454 intent(in) :: Ent_bl
1457 real,
dimension(SZI_(G)),
intent(in) :: I_dSkbp1
1460 real,
dimension(SZI_(G)),
intent(in) :: F_kb
1462 integer,
intent(in) :: kmb
1463 integer,
intent(in) :: i
1465 real,
dimension(SZI_(G)),
intent(inout) :: ea_kb
1467 real,
optional,
intent(in) :: tol_in
1470 real :: max_ea, min_ea
1471 real :: err, err_min, err_max
1473 real :: val, tolerance, tol1
1476 logical :: bisect_next, Newton
1477 real,
dimension(SZI_(G)) :: dS_kb
1478 real,
dimension(SZI_(G)) :: maxF, ent_maxF, zeros
1479 real,
dimension(SZI_(G)) :: ddSkb_dE
1481 integer,
parameter :: MAXIT = 30
1483 ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
1484 max_ea = ea_kb(i) ; min_ea = 0.0
1485 val = ds_kbp1 * f_kb(i)
1488 tolerance = cs%Tolerance_Ent
1489 if (
present(tol_in)) tolerance = tol_in
1490 bisect_next = .true.
1492 call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
1495 err = ds_kb(i) * ea_kb(i) - val
1496 derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
1499 if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea)))
return
1501 if (err == 0.0)
then ;
return
1502 elseif (err > 0.0)
then
1503 max_ea = ea_kb(i) ; err_max = err
1506 if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i) - min_ea) > err) .and. &
1507 (derr_dea*(max_ea - ea_kb(i)) > -1.0*err))
then
1508 ea_kb(i) = ea_kb(i) - err / derr_dea
1510 ea_kb(i) = 0.5*(max_ea+min_ea)
1516 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, ea_kb, &
1517 kmb, i, i, g, gv, cs, maxf, ent_maxf, f_thresh = f_kb)
1518 err_max = ds_kbp1 * maxf(i) - val
1521 if (err_max <= 0.0)
then
1522 ea_kb(i) = ent_maxf(i) ;
return
1524 max_ea = ent_maxf(i)
1525 ea_kb(i) = 0.5*(max_ea+min_ea)
1533 call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
1536 err = ds_kb(i) * ea_kb(i) - val
1537 derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
1543 max_ea = ea_kb(i) ; err_max = err
1544 if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) newton = .true.
1546 min_ea = ea_kb(i) ; err_min = err
1547 if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) newton = .true.
1551 ea_kb(i) = ea_kb(i) - err / derr_dea
1552 elseif (bisect_next)
then
1553 ea_kb(i) = 0.5*(max_ea+min_ea)
1554 bisect_next = .false.
1556 ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))
1557 bisect_next = .true.
1560 tol1 = tolerance ;
if (err > 0.0) tol1 = 0.099*tolerance
1561 if (ds_kb(i) <= ds_kbp1)
then
1562 if (abs(ea_kb(i) - ea_prev) <= tol1)
return
1564 if (ds_kbp1*abs(ea_kb(i) - ea_prev) <= ds_kb(i)*tol1)
return
1574 subroutine determine_ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &
1575 min_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, Ent, &
1576 error, err_min_eakb0, err_max_eakb0, F_kb, dFdfm_kb)
1579 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: h_bl
1581 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: Sref
1585 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: Ent_bl
1588 real,
dimension(SZI_(G)),
intent(in) :: I_dSkbp1
1592 real,
dimension(SZI_(G)),
intent(in) :: dtKd_kb
1595 real,
dimension(SZI_(G)),
intent(in) :: ea_kbp1
1597 real,
dimension(SZI_(G)),
intent(in) :: min_eakb
1599 real,
dimension(SZI_(G)),
intent(in) :: max_eakb
1601 integer,
intent(in) :: kmb
1602 integer,
intent(in) :: is
1603 integer,
intent(in) :: ie
1604 logical,
dimension(SZI_(G)),
intent(in) :: do_i
1607 real,
dimension(SZI_(G)),
intent(inout) :: Ent
1610 real,
dimension(SZI_(G)),
optional,
intent(out) :: error
1613 real,
dimension(SZI_(G)),
optional,
intent(in) :: err_min_eakb0
1616 real,
dimension(SZI_(G)),
optional,
intent(in) :: err_max_eakb0
1619 real,
dimension(SZI_(G)),
optional,
intent(out) :: F_kb
1623 real,
dimension(SZI_(G)),
optional,
intent(out) :: dFdfm_kb
1631 real,
dimension(SZI_(G)) :: &
1638 ddskb_de, ddslay_de, &
1643 error_mine, error_maxe
1653 logical,
dimension(SZI_(G)) :: false_position
1655 logical,
dimension(SZI_(G)) :: redo_i
1659 integer,
parameter :: MAXIT = 30
1661 if (.not.cs%bulkmixedlayer)
then
1662 call mom_error(fatal,
"determine_Ea_kb should not be called "//&
1663 "unless BULKMIXEDLAYER is defined.")
1665 tolerance = cs%Tolerance_Ent
1666 large_err = gv%m_to_H**2 * 1.0e30
1668 do i=is,ie ; redo_i(i) = do_i(i) ;
enddo
1670 do i=is,ie ;
if (do_i(i))
then
1675 e_min(i) = min_eakb(i) ; e_max(i) = max_eakb(i)
1676 error_mine(i) = -large_err ; error_maxe(i) = large_err
1677 false_position(i) = .true.
1679 if (
present(err_min_eakb0)) error_mine(i) = err_min_eakb0(i) - e_min(i) * ea_kbp1(i)
1680 if (
present(err_max_eakb0)) error_maxe(i) = err_max_eakb0(i) - e_max(i) * ea_kbp1(i)
1682 if ((error_maxe(i) <= 0.0) .or. (error_mine(i) >= 0.0))
then
1684 if (error_maxe(i) <= 0.0)
then
1686 ent(i) = e_max(i) ; err(i) = error_maxe(i)
1688 ent(i) = e_min(i) ; err(i) = error_mine(i)
1696 do_any = .false. ;
do i=is,ie ;
if (redo_i(i)) do_any = .true. ;
enddo
1697 if (.not.do_any)
exit
1698 call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., ds_kb, &
1699 ddskb_de, ds_lay, ddslay_de, do_i_in = redo_i)
1700 do i=is,ie ;
if (redo_i(i))
then
1703 el = 0.0 ;
if (2.0*ent_bl(i,kmb+1) >= ent(i)) el = 1.0
1704 fa = (1.0 + el) + ds_kb(i)*i_dskbp1(i)
1705 fk = dtkd_kb(i) * (ds_lay(i)/ds_kb(i))
1706 fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + el*2.0*ent_bl(i,kmb+1)
1707 if (fm > -gv%Angstrom_H) fm = fm + gv%Angstrom_H
1708 err(i) = (fa * ent(i)**2 - fm * ent(i)) - fk
1709 derror_de(i) = ((2.0*fa + (ddskb_de(i)*i_dskbp1(i))*ent(i))*ent(i) - fm) - &
1710 dtkd_kb(i) * (ddslay_de(i)*ds_kb(i) - ddskb_de(i)*ds_lay(i))/(ds_kb(i)**2)
1712 if (err(i) == 0.0)
then
1713 redo_i(i) = .false. ; cycle
1714 elseif (err(i) > 0.0)
then
1715 e_max(i) = ent(i) ; error_maxe(i) = err(i)
1717 e_min(i) = ent(i) ; error_mine(i) = err(i)
1721 if ((it == 1) .or. (derror_de(i) <= 0.0))
then
1726 fr = sqrt(fm**2 + 4.0*fa*fk)
1728 ent(i) = (fm + fr) / (2.0 * fa)
1730 ent(i) = (2.0 * fk) / (fr - fm)
1733 if ((ent(i) > e_max(i)) .or. (ent(i) < e_min(i))) &
1734 ent(i) = 0.5*(e_max(i) + e_min(i))
1735 elseif (((e_max(i)-ent(i))*derror_de(i) > -err(i)) .and. &
1736 ((ent(i)-e_min(i))*derror_de(i) > err(i)) )
then
1739 ent(i) = ent(i) - err(i) / derror_de(i)
1740 elseif (false_position(i) .and. &
1741 (error_maxe(i) - error_mine(i) < 0.9*large_err))
then
1743 ent(i) = e_min(i) + (e_max(i)-e_min(i)) * &
1744 (-error_mine(i)/(error_maxe(i) - error_mine(i)))
1745 false_position(i) = .false.
1747 ent(i) = 0.5*(e_max(i) + e_min(i))
1748 false_position(i) = .true.
1751 if (abs(e_prev - ent(i)) < tolerance)
then
1752 err_est = err(i) + (ent(i) - e_prev) * derror_de(i)
1753 if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &
1754 (abs(err_est) < abs(tolerance*derror_de(i)))) redo_i(i) = .false.
1761 if (
present(f_kb) .or.
present(dfdfm_kb)) &
1762 call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., &
1763 ds_kb, do_i_in = do_i)
1765 if (
present(f_kb))
then ;
do i=is,ie ;
if (do_i(i))
then
1766 f_kb(i) = ent(i) * (ds_kb(i) * i_dskbp1(i))
1767 endif ;
enddo ;
endif
1768 if (
present(error))
then ;
do i=is,ie ;
if (do_i(i))
then
1770 endif ;
enddo ;
endif
1771 if (
present(dfdfm_kb))
then ;
do i=is,ie ;
if (do_i(i))
then
1775 if (derror_de(i) > 0.0)
then
1776 dfdfm_kb(i) = ((ds_kb(i) + ent(i) * ddskb_de(i)) * i_dskbp1(i)) * &
1777 (ent(i) / derror_de(i))
1781 endif ;
enddo ;
endif
1786 subroutine find_maxf_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, &
1787 kmb, is, ie, G, GV, CS, maxF, ent_maxF, do_i_in, &
1788 F_lim_maxent, F_thresh)
1791 real,
dimension(SZI_(G),SZK_(G)), &
1793 real,
dimension(SZI_(G),SZK_(G)), &
1795 real,
dimension(SZI_(G),SZK_(G)), &
1796 intent(in) :: Ent_bl
1799 real,
dimension(SZI_(G)),
intent(in) :: I_dSkbp1
1803 real,
dimension(SZI_(G)),
intent(in) :: min_ent_in
1805 real,
dimension(SZI_(G)),
intent(in) :: max_ent_in
1807 integer,
intent(in) :: kmb
1808 integer,
intent(in) :: is
1809 integer,
intent(in) :: ie
1811 real,
dimension(SZI_(G)),
intent(out) :: maxF
1814 real,
dimension(SZI_(G)), &
1815 optional,
intent(out) :: ent_maxF
1816 logical,
dimension(SZI_(G)), &
1817 optional,
intent(in) :: do_i_in
1819 real,
dimension(SZI_(G)), &
1820 optional,
intent(out) :: F_lim_maxent
1824 real,
dimension(SZI_(G)), &
1825 optional,
intent(in) :: F_thresh
1835 real,
dimension(SZI_(G)) :: &
1837 minent, maxent, ent_best, &
1839 F_maxent, F_minent, F, F_best, &
1840 dF_dent, dF_dE_max, dF_dE_min, dF_dE_best, &
1841 dS_kb, dS_kb_lim, ddSkb_dE, dS_anom_lim, &
1842 chg_prev, chg_pre_prev
1843 real :: dF_dE_mean, maxslope, minslope
1845 real :: ratio_select_end
1846 real :: rat, max_chg, min_chg, chg1, chg2, chg
1847 logical,
dimension(SZI_(G)) :: do_i, last_it, need_bracket, may_use_best
1848 logical :: doany, OK1, OK2, bisect, new_min_bound
1849 integer :: i, it, is1, ie1
1850 integer,
parameter :: MAXIT = 20
1852 tolerance = cs%Tolerance_Ent
1854 if (
present(do_i_in))
then
1855 do i=is,ie ; do_i(i) = do_i_in(i) ;
enddo
1857 do i=is,ie ; do_i(i) = .true. ;
enddo
1861 call determine_dskb(h_bl, sref, ent_bl, max_ent_in, is, ie, kmb, g, gv, .false., &
1862 ds_kb, ddskb_de, ds_anom_lim=ds_anom_lim)
1863 ie1 = is-1 ; doany = .false.
1865 ds_kb_lim(i) = ds_kb(i) + ds_anom_lim(i)
1866 f_max_ent_in(i) = max_ent_in(i)*ds_kb_lim(i)*i_dskbp1(i)
1867 maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)
1868 if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i)))
then
1869 f_best(i) = max_ent_in(i)*ds_kb(i)*i_dskbp1(i)
1870 ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)
1873 f_maxent(i) = maxent(i) * ds_kb(i) * i_dskbp1(i)
1874 df_de_max(i) = (ds_kb(i) + maxent(i)*ddskb_de(i)) * i_dskbp1(i)
1875 doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.
1880 ie1 = is-1 ;
do i=is,ie ;
if (do_i(i)) ie1 = i ;
enddo
1881 do i=ie1,is,-1 ;
if (do_i(i)) is1 = i ;
enddo
1883 call determine_dskb(h_bl, sref, ent_bl, minent, is1, ie1, kmb, g, gv, .false., &
1884 ds_kb, ddskb_de, do_i_in = do_i)
1885 do i=is1,ie1 ;
if (do_i(i))
then
1886 f_minent(i) = minent(i) * ds_kb(i) * i_dskbp1(i)
1887 df_de_min(i) = (ds_kb(i) + minent(i)*ddskb_de(i)) * i_dskbp1(i)
1890 ratio_select_end = 0.9
1892 ratio_select_end = 0.5*ratio_select_end
1893 do i=is1,ie1 ;
if (do_i(i))
then
1894 if (need_bracket(i))
then
1895 df_de_mean = (f_maxent(i) - f_minent(i)) / (maxent(i) - minent(i))
1896 maxslope = max(df_de_mean, df_de_min(i), df_de_max(i))
1897 minslope = min(df_de_mean, df_de_min(i), df_de_max(i))
1898 if (f_minent(i) >= f_maxent(i))
then
1899 if (df_de_min(i) > 0.0)
then ; rat = 0.02
1900 elseif (maxslope < ratio_select_end*minslope)
then
1902 f_best(i) = f_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0
1904 else ; rat = 0.382 ;
endif
1906 if (df_de_max(i) < 0.0)
then ; rat = 0.98
1907 elseif (minslope > ratio_select_end*maxslope)
then
1909 f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0
1911 else ; rat = 0.618 ;
endif
1914 if (rat >= 0.0) ent(i) = rat*maxent(i) + (1.0-rat)*minent(i)
1915 if (((maxent(i) - minent(i)) < tolerance) .or. (it==maxit)) &
1918 chg1 = 2.0*(maxent(i) - minent(i)) ; chg2 = chg1
1919 if (df_de_best(i) > 0)
then
1920 max_chg = maxent(i) - ent_best(i) ; min_chg = 0.0
1922 max_chg = 0.0 ; min_chg = minent(i) - ent_best(i)
1924 if (max_chg - min_chg < 2.0*tolerance) last_it(i) = .true.
1925 if (df_de_max(i) /= df_de_best(i)) &
1926 chg1 = (maxent(i) - ent_best(i))*df_de_best(i) / &
1927 (df_de_best(i) - df_de_max(i))
1928 if (df_de_min(i) /= df_de_best(i)) &
1929 chg2 = (minent(i) - ent_best(i))*df_de_best(i) / &
1930 (df_de_best(i) - df_de_min(i))
1931 ok1 = ((chg1 < max_chg) .and. (chg1 > min_chg))
1932 ok2 = ((chg2 < max_chg) .and. (chg2 > min_chg))
1933 if (.not.(ok1 .or. ok2))
then ; bisect = .true. ;
else
1934 if (ok1 .and. ok2)
then
1935 chg = chg1 ;
if (abs(chg2) < abs(chg1)) chg = chg2
1936 elseif (ok1)
then ; chg = chg1
1937 else ; chg = chg2 ;
endif
1938 if (abs(chg) > 0.5*abs(chg_pre_prev(i)))
then ; bisect = .true.
1939 else ; bisect = .false. ;
endif
1941 chg_pre_prev(i) = chg_prev(i)
1943 if (df_de_best(i) > 0.0)
then
1944 ent(i) = 0.5*(maxent(i) + ent_best(i))
1945 chg_prev(i) = 0.5*(maxent(i) - ent_best(i))
1947 ent(i) = 0.5*(minent(i) + ent_best(i))
1948 chg_prev(i) = 0.5*(minent(i) - ent_best(i))
1951 if (abs(chg) < tolerance) chg = sign(tolerance,chg)
1952 ent(i) = ent_best(i) + chg
1958 if (mod(it,3) == 0)
then
1959 ie1 = is-1 ;
do i=is1,ie ;
if (do_i(i)) ie1 = i ;
enddo
1960 do i=ie1,is,-1 ;
if (do_i(i)) is1 = i ;
enddo
1963 call determine_dskb(h_bl, sref, ent_bl, ent, is1, ie1, kmb, g, gv, .false., &
1964 ds_kb, ddskb_de, do_i_in = do_i)
1965 do i=is1,ie1 ;
if (do_i(i))
then
1966 f(i) = ent(i)*ds_kb(i)*i_dskbp1(i)
1967 df_dent(i) = (ds_kb(i) + ent(i)*ddskb_de(i)) * i_dskbp1(i)
1970 if (
present(f_thresh))
then ;
do i=is1,ie1 ;
if (do_i(i))
then
1971 if (f(i) >= f_thresh(i))
then
1972 f_best(i) = f(i) ; ent_best(i) = ent(i) ; do_i(i) = .false.
1974 endif ;
enddo ;
endif
1977 do i=is1,ie1 ;
if (do_i(i))
then
1978 if (.not.last_it(i)) doany = .true.
1979 if (last_it(i))
then
1980 if (need_bracket(i))
then
1981 if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i)))
then
1982 f_best(i) = f(i) ; ent_best(i) = ent(i)
1983 elseif (f_maxent(i) > f_minent(i))
then
1984 f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i)
1986 f_best(i) = f_minent(i) ; ent_best(i) = minent(i)
1988 elseif (f(i) > f_best(i))
then
1989 f_best(i) = f(i) ; ent_best(i) = ent(i)
1992 elseif (need_bracket(i))
then
1993 if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i)))
then
1994 need_bracket(i) = .false.
1995 chg_prev(i) = (maxent(i) - minent(i))
1996 chg_pre_prev(i) = 2.0*chg_prev(i)
1997 ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
1998 elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i)))
then
1999 new_min_bound = .true.
2000 elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i)))
then
2001 new_min_bound = .false.
2007 new_min_bound = .true.
2008 if (df_de_min(i) > 0.0 .or. (f_minent(i) >= f_maxent(i))) &
2009 new_min_bound = .false.
2011 if (need_bracket(i))
then
2012 if (new_min_bound)
then
2013 minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
2015 maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
2019 if (f(i) >= f_best(i))
then
2020 if (ent(i) > ent_best(i))
then
2021 minent(i) = ent_best(i) ; f_minent(i) = f_best(i) ; df_de_min(i) = df_de_best(i)
2023 maxent(i) = ent_best(i) ; f_maxent(i) = f_best(i) ; df_de_max(i) = df_de_best(i)
2025 ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
2027 if (ent(i) < ent_best(i))
then
2028 minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
2030 maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
2033 if ((maxent(i) - minent(i)) <= tolerance) do_i(i) = .false.
2036 if (.not.doany)
exit
2040 if (
present(f_lim_maxent))
then
2044 f_lim_maxent(i) = f_max_ent_in(i)
2045 if (
present(ent_maxf)) ent_maxf(i) = ent_best(i)
2051 may_use_best(i) = (ent_best(i) /= max_ent_in(i))
2052 if (may_use_best(i)) doany = .true.
2056 call determine_dskb(h_bl, sref, ent_bl, ent_best, is, ie, kmb, g, gv, .true., &
2059 f_best(i) = ent_best(i)*ds_kb_lim(i)*i_dskbp1(i)
2062 if ((f_best(i) > f_max_ent_in(i)) .and. (may_use_best(i)))
then
2064 if (
present(ent_maxf)) ent_maxf(i) = ent_best(i)
2066 maxf(i) = f_max_ent_in(i)
2067 if (
present(ent_maxf)) ent_maxf(i) = max_ent_in(i)
2072 do i=is,ie ; maxf(i) = f_max_ent_in(i) ;
enddo
2073 if (
present(ent_maxf))
then
2074 do i=is,ie ; ent_maxf(i) = max_ent_in(i) ;
enddo
2084 type(time_type),
intent(in) :: time
2090 type(
diag_ctrl),
target,
intent(inout) :: diag
2103 real :: decay_length, dt, kd
2105 #include "version_variable.h"
2106 character(len=40) :: mdl =
"MOM_entrain_diffusive"
2108 if (
associated(cs))
then
2109 call mom_error(warning,
"entrain_diffusive_init called with an associated "// &
2110 "control structure.")
2117 cs%bulkmixedlayer = (gv%nkml > 0)
2121 call get_param(param_file, mdl,
"CORRECT_DENSITY", cs%correct_density, &
2122 "If true, and USE_EOS is true, the layer densities are "//&
2123 "restored toward their target values by the diapycnal "//&
2124 "mixing, as described in Hallberg (MWR, 2000).", &
2126 call get_param(param_file, mdl,
"MAX_ENT_IT", cs%max_ent_it, &
2127 "The maximum number of iterations that may be used to "//&
2128 "calculate the interior diapycnal entrainment.", default=5)
2130 call get_param(param_file, mdl,
"KD", kd, fail_if_missing=.true.)
2131 call get_param(param_file, mdl,
"DT", dt, &
2132 "The (baroclinic) dynamics time step.", units =
"s", &
2133 fail_if_missing=.true.)
2135 call get_param(param_file, mdl,
"TOLERANCE_ENT", cs%Tolerance_Ent, &
2136 "The tolerance with which to solve for entrainment values.", &
2137 units=
"m", default=max(100.0*gv%Angstrom_m,1.0e-4*sqrt(dt*kd)), scale=gv%m_to_H)
2139 cs%Rho_sig_off = 1000.0*us%kg_m3_to_R
2141 cs%id_Kd = register_diag_field(
'ocean_model',
'Kd_effective', diag%axesTL, time, &
2142 'Diapycnal diffusivity as applied',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
2143 cs%id_diff_work = register_diag_field(
'ocean_model',
'diff_work', diag%axesTi, time, &
2144 'Work actually done by diapycnal diffusion across each interface',
'W m-2', &
2145 conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2154 if (
associated(cs))
deallocate(cs)