This subroutine calculates ea and eb, the rates at which a layer entrains from the layers above and below. The entrainment rates are proportional to the buoyancy flux in a layer and inversely proportional to the density differences between layers. The scheme that is used here is described in detail in Hallberg, Mon. Wea. Rev. 2000.
53 type(ocean_grid_type),
intent(in) :: G
54 type(verticalGrid_type),
intent(in) :: GV
55 type(unit_scale_type),
intent(in) :: US
56 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
58 type(thermo_var_ptrs),
intent(in) :: tv
61 type(forcing),
intent(in) :: fluxes
63 real,
intent(in) :: dt
64 type(entrain_diffusive_CS),
pointer :: CS
66 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
69 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
72 integer,
dimension(SZI_(G),SZJ_(G)), &
73 optional,
intent(inout) :: kb_out
76 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
77 optional,
intent(in) :: Kd_Lay
79 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
80 optional,
intent(in) :: Kd_int
89 real,
dimension(SZI_(G),SZK_(G)) :: &
91 real,
dimension(SZI_(G),SZK_(G)+1) :: &
93 real,
dimension(SZI_(G),SZK_(G)) :: &
106 real,
dimension(SZI_(G),SZK_(G)+1) :: &
109 real,
allocatable,
dimension(:,:,:) :: &
116 real :: hm, fm, fr, fk
119 real :: c1(SZI_(G),SZK_(G))
121 real,
dimension(SZI_(G)) :: &
152 real,
dimension(SZI_(G),SZK_(G)) :: &
160 real,
dimension(SZI_(G),SZK_(G)) :: &
175 real,
dimension(SZI_(G)) :: &
199 logical :: do_i(SZI_(G)), did_i(SZI_(G)), reiterate, correct_density
200 integer :: it, i, j, k, is, ie, js, je, nz, K2, kmb
201 integer :: kb(SZI_(G))
203 integer :: kb_min_act
205 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
206 angstrom = gv%Angstrom_H
207 h_neglect = gv%H_subroundoff
209 if (.not.
associated(cs))
call mom_error(fatal, &
210 "MOM_entrain_diffusive: Module must be initialized before it is used.")
212 if (.not.(
present(kd_lay) .or.
present(kd_int)))
call mom_error(fatal, &
213 "MOM_entrain_diffusive: Either Kd_Lay or Kd_int must be present in call.")
215 if ((.not.cs%bulkmixedlayer .and. .not.
associated(fluxes%buoy)) .and. &
216 (
associated(fluxes%lprec) .or.
associated(fluxes%evap) .or. &
217 associated(fluxes%sens) .or.
associated(fluxes%sw)))
then
218 if (is_root_pe())
call mom_error(note,
"Calculate_Entrainment: &
219 &The code to handle evaporation and precipitation without &
220 &a bulk mixed layer has not been implemented.")
221 if (is_root_pe())
call mom_error(fatal, &
222 "Either define BULKMIXEDLAYER in MOM_input or use fluxes%buoy &
223 &and a linear equation of state to drive the model.")
226 tolerance = cs%Tolerance_Ent
227 kmb = gv%nk_rho_varies
228 k2 = max(kmb+1,2) ; kb_min = k2
229 if (.not. cs%bulkmixedlayer)
then
233 do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ;
enddo
236 do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ;
enddo
239 if (cs%id_diff_work > 0)
allocate(diff_work(g%isd:g%ied,g%jsd:g%jed,nz+1))
240 if (cs%id_Kd > 0)
allocate(kd_eff(g%isd:g%ied,g%jsd:g%jed,nz))
242 correct_density = (cs%correct_density .and.
associated(tv%eqn_of_state))
243 if (correct_density)
then
266 do i=is,ie ; kb(i) = 1 ;
enddo
268 if (
present(kd_lay))
then
269 do k=1,nz ;
do i=is,ie
270 dtkd(i,k) = gv%Z_to_H**2 * (dt * kd_lay(i,j,k))
272 if (
present(kd_int))
then
273 do k=1,nz+1 ;
do i=is,ie
274 dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
277 do k=2,nz ;
do i=is,ie
278 dtkd_int(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_lay(i,j,k-1) + kd_lay(i,j,k)))
282 do k=1,nz ;
do i=is,ie
283 dtkd(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_int(i,j,k)+kd_int(i,j,k+1)))
285 dO k=1,nz+1 ;
do i=is,ie
286 dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
290 do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ;
enddo
291 do i=is,ie ; ds_dsp1(i,nz) = 0.0 ;
enddo
292 do i=is,ie ; dsp1_ds(i,nz) = 0.0 ;
enddo
294 do k=2,nz-1 ;
do i=is,ie
295 ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
298 if (cs%bulkmixedlayer)
then
302 call set_ent_bl(h, dtkd_int, tv, kb, kmb, do_i, g, gv, cs, j, ent_bl, sref, h_bl)
305 dtkd_kb(i) = 0.0 ;
if (kb(i) < nz) dtkd_kb(i) = dtkd(i,kb(i))
308 do i=is,ie ; ent_bl(i,kmb+1) = 0.0 ;
enddo
311 do k=2,nz-1 ;
do i=is,ie
312 dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
313 i2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))
314 grats(i,k) = 2.0*(2.0+(dsp1_ds(i,k)+ds_dsp1(i,k)))
320 if (cs%bulkmixedlayer)
then
323 htot(i) = h(i,j,1) - angstrom
325 do k=2,kmb ;
do i=is,ie
326 htot(i) = htot(i) + (h(i,j,k) - angstrom)
329 max_eakb(i) = max(ent_bl(i,kmb+1) + 0.5*htot(i), htot(i))
330 i_dskbp1(i) = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
338 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, max_eakb, kmb, &
339 is, ie, g, gv, cs, maxf_kb, eakb_maxf, do_i, f_kb_maxent)
340 do i=is,ie ;
if (kb(i) <= nz)
then
341 maxf(i,kb(i)) = max(0.0, maxf_kb(i), f_kb_maxent(i))
342 if ((maxf_kb(i) > f_kb_maxent(i)) .and. (eakb_maxf(i) >= htot(i))) &
343 max_eakb(i) = eakb_maxf(i)
346 do i=is,ie ; ea_kbp1(i) = 0.0 ; eakb(i) = max_eakb(i) ;
enddo
347 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &
348 max_eakb, max_eakb, kmb, is, ie, do_i, g, gv, cs, eakb, &
349 error=err_max_eakb0, f_kb=f_kb)
353 do i=is,ie ; min_eakb(i) = min(htot(i), max_eakb(i)) ;
enddo
354 call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, min_eakb, max_eakb, &
355 kmb, is, ie, g, gv, cs, f_kb_maxent, do_i_in = do_i)
358 if ((.not.do_i(i)) .or. (err_max_eakb0(i) >= 0.0))
then
359 eakb(i) = 0.0 ; min_eakb(i) = 0.0
361 eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
368 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &
369 zeros, max_eakb, kmb, is, ie, do_i, g, gv, cs, min_eakb, &
370 error=err_min_eakb0, f_kb=f_kb, err_max_eakb0=err_max_eakb0)
372 do i=is,ie ;
if ((kb(i)<nz) .and. (kb_min>kb(i))) kb_min = kb(i) ;
enddo
381 htot(i) = h(i,j,1) - angstrom
383 if (
associated(fluxes%buoy))
then ;
do i=is,ie
384 maxf(i,1) = (dt*fluxes%buoy(i,j)) / (gv%g_prime(2)*us%m_to_Z)
390 do k=kb_min,nz-1 ;
do i=is,ie
391 if ((k == kb(i)+1) .and. cs%bulkmixedlayer)
then
392 maxf(i,k) = ds_dsp1(i,k)*(f_kb_maxent(i) + htot(i))
393 htot(i) = htot(i) + (h(i,j,k) - angstrom)
394 elseif (k > kb(i))
then
395 maxf(i,k) = ds_dsp1(i,k)*(maxf(i,k-1) + htot(i))
396 htot(i) = htot(i) + (h(i,j,k) - angstrom)
401 if (.not.cs%bulkmixedlayer)
then
402 maxf_correct(i) = max(0.0, -(maxf(i,nz-1) + htot(i)))
404 htot(i) = h(i,j,nz) - angstrom
406 if (.not.cs%bulkmixedlayer)
then
407 do_any = .false. ;
do i=is,ie ;
if (maxf_correct(i) > 0.0) do_any = .true. ;
enddo
409 do k=nz-1,1,-1 ;
do i=is,ie
410 maxf(i,k) = maxf(i,k) + maxf_correct(i)
411 maxf_correct(i) = maxf_correct(i) * dsp1_ds(i,k)
415 do k=nz-1,kb_min,-1 ;
do i=is,ie ;
if (do_i(i))
then
417 maxf(i,k) = min(maxf(i,k),dsp1_ds(i,k+1)*maxf(i,k+1) + htot(i))
418 htot(i) = htot(i) + (h(i,j,k) - angstrom)
419 if ( (k == kb(i)) .and. ((maxf(i,k) < f_kb(i)) .or. &
420 (maxf(i,k) < maxf_kb(i)) .and. (eakb_maxf(i) <= max_eakb(i))) )
then
425 if ((f_kb(i) <= maxf_kb(i)) .and. (eakb_maxf(i) <= max_eakb(i)))
then
426 eakb(i) = eakb_maxf(i)
428 eakb(i) = max_eakb(i)
430 call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
431 g, gv, cs, eakb, angstrom)
432 if ((eakb(i) < max_eakb(i)) .or. (eakb(i) < min_eakb(i)))
then
433 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, zeros, &
434 eakb, eakb, kmb, i, i, do_i, g, gv, cs, eakb, &
436 if (eakb(i) < max_eakb(i))
then
437 max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)
439 if (eakb(i) < min_eakb(i))
then
440 min_eakb(i) = eakb(i) ; err_min_eakb0(i) = err_eakb0(i)
445 endif ;
enddo ;
enddo
446 if (.not.cs%bulkmixedlayer)
then
448 maxf(i,1) = min(maxf(i,1),dsp1_ds(i,2)*maxf(i,2) + htot(i))
459 f(i,nz) = maxf(i,nz) ; minf(i,nz) = 0.0
463 if ((k==kb(i)) .and. (do_i(i)))
then
464 eakb(i) = min_eakb(i)
466 elseif ((k>kb(i)) .and. (do_i(i)))
then
471 hm = h(i,j,k) + h_neglect
474 f(i,k) = min(maxf(i,k), sqrt(ds_dsp1(i,k)*dtkd(i,k)), &
475 0.5*(ds_dsp1(i,k)+1.0) * (dtkd(i,k) / hm))
480 fk = dtkd(i,k) * grats(i,k)
481 minf(i,k) = min(maxf(i,k), &
482 0.9*(i2p2dsp1_ds(i,k) * fk / (hm + sqrt(hm*hm + fk))))
483 if (k==kb(i)) minf(i,k) = 0.0
493 is1 = ie+1 ; ie1 = is-1
494 do i=is,ie ;
if (do_i(i))
then ; is1 = i ;
exit ;
endif ;
enddo
495 do i=ie,is,-1 ;
if (do_i(i))
then ; ie1 = i ;
exit ;
endif ;
enddo
497 if (cs%bulkmixedlayer)
then
500 if (do_i(i) .and. (kb(i) < kb_min_act)) kb_min_act = kb(i)
507 if (do_i(i) .and. (kb(i) < nz)) &
508 ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*f(i,kb(i)+1)
510 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
511 max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
512 err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
518 do it=0,cs%max_ent_it-1
519 do i=is1,ie1 ;
if (do_i(i))
then
520 if (.not.cs%bulkmixedlayer) f(i,1) = min(f(i,1),maxf(i,1))
526 do k=kb_min_act,nz-1 ;
do i=is1,ie1 ;
if (do_i(i) .and. (k>=kb(i)))
then
528 if (cs%bulkmixedlayer .and. (k==kb(i)))
then
530 dfdfm(i,k) = dfdfm_kb(i)
533 fm = (f(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*f(i,k+1)
534 fk = grats(i,k)*dtkd(i,k)
535 fr = sqrt(fm*fm + fk)
538 f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fm+fr))
540 f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))
543 if ((f(i,k) >= maxf(i,k)) .or. (fr == 0.0))
then
546 dfdfm(i,k) = i2p2dsp1_ds(i,k) * ((fr + fm) / fr)
551 c1(i,k) = dfdfm(i,k-1)*(dsp1_ds(i,k)*b1(i))
552 b1(i) = 1.0 / (1.0 - c1(i,k)*dfdfm(i,k))
553 f(i,k) = min(b1(i)*(f(i,k)-fprev(i,k)) + fprev(i,k), maxf(i,k))
554 if (f(i,k) >= maxf(i,k)) dfdfm(i,k) = 0.0
557 endif ;
enddo ;
enddo
559 do k=nz-2,kb_min_act,-1 ;
do i=is1,ie1
560 if (do_i(i) .and. (k > kb(i))) &
561 f(i,k) = min((f(i,k)+c1(i,k+1)*(f(i,k+1)-fprev(i,k+1))),maxf(i,k))
564 if (cs%bulkmixedlayer)
then
566 if (do_i(i) .and. (kb(i) < nz))
then
568 ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(f(i,kb(i)+1), minf(i,kb(i)+1))
573 call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
574 max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
575 err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
578 if (do_i(i) .and. (kb(i) < nz)) f(i,kb(i)) = f_kb(i)
583 if (it < cs%max_ent_it-1)
then
586 if (cs%bulkmixedlayer)
then ;
do i=is1,ie1 ;
if (do_i(i))
then
587 eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
588 endif ;
enddo ;
endif
590 did_i(i) = do_i(i) ; do_i(i) = .false.
592 do k=kb_min_act,nz-1 ;
do i=is1,ie1
593 if (did_i(i) .and. (k >= kb(i)))
then
594 if (f(i,k) < minf(i,k))
then
596 do_i(i) = .true. ; reiterate = .true.
597 elseif (k > kb(i))
then
598 if ((abs(f(i,k) - fprev(i,k)) > tolerance) .or. &
599 ((h(i,j,k) + ((1.0+dsp1_ds(i,k))*f(i,k) - &
600 (f(i,k-1) + dsp1_ds(i,k+1)*f(i,k+1)))) < 0.9*angstrom))
then
601 do_i(i) = .true. ; reiterate = .true.
607 if (h(i,j,k) + ((f(i,k) + eakb(i)) - &
608 (eb_kmb(i) + dsp1_ds(i,k+1)*f(i,k+1))) < 0.9*angstrom)
then
609 do_i(i) = .true. ; reiterate = .true.
614 if (.not.reiterate)
exit
620 if (it == (cs%max_ent_it))
then
623 do i=is1,ie1 ;
if (do_i(i))
then
624 f(i,nz-1) = max(f(i,nz-1), min(minf(i,nz-1), 0.0))
625 if (kb(i) >= nz-1)
then ; ea_kbp1(i) = 0.0 ;
endif
627 do k=nz-2,kb_min_act,-1 ;
do i=is1,ie1 ;
if (do_i(i))
then
629 f(i,k) = min(max(minf(i,k),f(i,k)), (dsp1_ds(i,k+1)*f(i,k+1) + &
630 max(((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + &
631 (h(i,j,k+1) - angstrom)), 0.5*(h(i,j,k+1)-angstrom))))
632 elseif (k==kb(i))
then
633 ea_kbp1(i) = dsp1_ds(i,k+1)*f(i,k+1)
634 h_avail = dsp1_ds(i,k+1)*f(i,k+1) + max(0.5*(h(i,j,k+1)-angstrom), &
635 ((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + (h(i,j,k+1) - angstrom)))
636 if ((f(i,k) > 0.0) .and. (f(i,k) > h_avail))
then
637 f_kb(i) = max(0.0, h_avail)
639 if ((f_kb(i) < maxf_kb(i)) .and. (eakb_maxf(i) <= eakb(i))) &
640 eakb(i) = eakb_maxf(i)
641 call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
645 endif ;
enddo ;
enddo
648 if (cs%bulkmixedlayer)
then ;
do i=is1,ie1
649 if (do_i(i) .and. (kb(i) < nz))
then
650 h_avail = eakb(i) + max(0.5*(h_bl(i,kmb+1)-angstrom), &
651 (f_kb(i)-ea_kbp1(i)) + (h_bl(i,kmb+1)-angstrom))
653 ent_bl(i,kmb+1) = min(ent_bl(i,kmb+1),0.5*(eakb(i) + h_avail))
655 eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
660 do k=kb_min_act+1,nz-1 ;
do i=is1,ie1 ;
if (do_i(i))
then
661 if ((.not.cs%bulkmixedlayer) .or. (k > kb(i)+1))
then
662 f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + &
663 dsp1_ds(i,k-1)*f(i,k-1)) - f(i,k-2)) + (h(i,j,k-1) - angstrom)))
664 f(i,k) = max(f(i,k),min(minf(i,k),0.0))
665 elseif (k == kb(i)+1)
then
666 f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + eakb(i)) - &
667 eb_kmb(i)) + (h(i,j,k-1) - angstrom)))
668 f(i,k) = max(f(i,k),min(minf(i,k),0.0))
670 endif ;
enddo ;
enddo
673 call f_to_ent(f, h, kb, kmb, j, g, gv, cs, dsp1_ds, eakb, ent_bl, ea, eb)
677 if (correct_density)
then
679 h_guess(i,1) = (h(i,j,1) - angstrom) + (eb(i,j,1) - ea(i,j,2))
680 h_guess(i,nz) = (h(i,j,nz) - angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))
681 if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0
682 if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0
684 do k=2,nz-1 ;
do i=is,ie
685 h_guess(i,k) = (h(i,j,k) - angstrom) + ((ea(i,j,k) - eb(i,j,k-1)) + &
686 (eb(i,j,k) - ea(i,j,k+1)))
687 if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0
689 if (cs%bulkmixedlayer)
then
690 call determine_dskb(h_bl, sref, ent_bl, eakb, is, ie, kmb, g, gv, &
691 .true., ds_kb, ds_anom_lim=ds_anom_lim)
693 call calculate_density(tv%T(is:ie,j,k), tv%S(is:ie,j,k), pres(is:ie), &
694 rcv(is:ie), 1, ie-is+1, tv%eqn_of_state)
696 if ((k>kb(i)) .and. (f(i,k) > 0.0))
then
703 f_cor = h(i,j,k) * min(1.0 , max(-ds_dsp1(i,k), &
704 (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )
710 if (f_cor > 0.0)
then
711 f_cor = min(f_cor, 0.9*f(i,k), ds_dsp1(i,k)*0.5*h_guess(i,k), &
714 f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
715 0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
718 ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
719 eb(i,j,k) = eb(i,j,k) + f_cor
720 elseif ((k==kb(i)) .and. (f(i,k) > 0.0))
then
724 ds_kb_eff = 2.0*ds_kb(i) - ds_anom_lim(i)
725 rho_cor = h(i,j,k) * (gv%Rlay(k)-rcv(i)) + eakb(i)*ds_anom_lim(i)
728 if (abs(rho_cor) < abs(0.9*eakb(i)*ds_kb_eff))
then
729 ea_cor = -rho_cor / ds_kb_eff
731 ea_cor = sign(0.9*eakb(i),-rho_cor*ds_kb_eff)
734 if (ea_cor > 0.0)
then
736 ea_cor = min(ea_cor, 0.5*(max_eakb(i) - eakb(i)), &
737 0.5*h_guess(i,k) / (ds_kb(i) * i_dskbp1(i)))
740 ea_cor = -min(-ea_cor, 0.5*h_guess(i,k), &
741 0.25*h_guess(i,k+1) / (ds_kb(i) * i_dskbp1(i)))
744 ea(i,j,k) = ea(i,j,k) + ea_cor
745 eb(i,j,k) = eb(i,j,k) - (ds_kb(i) * i_dskbp1(i)) * ea_cor
746 elseif (k < kb(i))
then
748 ea(i,j,k) = ea(i,j,k+1)
752 do k=kb_min-1,k2,-1 ;
do i=is,ie
753 ea(i,j,k) = ea(i,j,k+1)
761 h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
762 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
764 do k=kmb-1,2,-1 ;
do i=is,ie
766 eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)
769 h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
770 ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
773 eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
778 call calculate_density(tv%T(is:ie,j,k), tv%S(is:ie,j,k), pres(is:ie), &
779 rcv(is:ie), 1, ie-is+1, tv%eqn_of_state)
780 do i=is,ie ;
if (f(i,k) > 0.0)
then
785 f_cor = h(i,j,k) * min(dsp1_ds(i,k) , max(-1.0, &
786 (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )
792 if (f_cor >= 0.0)
then
793 f_cor = min(f_cor, 0.9*f(i,k), 0.5*dsp1_ds(i,k)*h_guess(i,k), &
796 f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
797 0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
800 ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
801 eb(i,j,k) = eb(i,j,k) + f_cor
808 if (cs%id_Kd > 0)
then
809 idt = gv%H_to_Z**2 / dt
810 do k=2,nz-1 ;
do i=is,ie
811 if (k<kb(i))
then ; kd_here = 0.0 ;
else
812 kd_here = f(i,k) * ( h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
813 (eb(i,j,k) - ea(i,j,k+1))) ) / (i2p2dsp1_ds(i,k) * grats(i,k))
816 kd_eff(i,j,k) = max(dtkd(i,k), kd_here)*idt
819 kd_eff(i,j,1) = dtkd(i,1)*idt
820 kd_eff(i,j,nz) = dtkd(i,nz)*idt
824 if (cs%id_diff_work > 0)
then
825 g_2dt = 0.5 * gv%H_to_Z**2*us%L_to_Z**2 * (gv%g_Earth / dt)
826 do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ;
enddo
827 if (
associated(tv%eqn_of_state))
then
828 if (
associated(fluxes%p_surf))
then
829 do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ;
enddo
831 do i=is,ie ; pressure(i) = 0.0 ;
enddo
834 do i=is,ie ; pressure(i) = pressure(i) + gv%H_to_Pa*h(i,j,k-1) ;
enddo
837 t_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))
838 s_eos(i) = 0.5*(tv%S(i,j,kmb) + tv%S(i,j,k))
840 t_eos(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
841 s_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
844 call calculate_density_derivs(t_eos, s_eos, pressure, &
845 drho_dt, drho_ds, is, ie-is+1, tv%eqn_of_state)
847 if ((k>kmb) .and. (k<kb(i)))
then ; diff_work(i,j,k) = 0.0
850 drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &
851 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))
853 drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
854 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))
856 diff_work(i,j,k) = g_2dt * drho * &
857 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
858 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
863 do k=2,nz ;
do i=is,ie
864 diff_work(i,j,k) = g_2dt * (gv%Rlay(k)-gv%Rlay(k-1)) * &
865 (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
866 eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
871 if (
present(kb_out))
then
872 do i=is,ie ; kb_out(i,j) = kb(i) ;
enddo
878 if (cs%id_Kd > 0)
call post_data(cs%id_Kd, kd_eff, cs%diag)
879 if (cs%id_Kd > 0)
deallocate(kd_eff)
880 if (cs%id_diff_work > 0)
call post_data(cs%id_diff_work, diff_work, cs%diag)
881 if (cs%id_diff_work > 0)
deallocate(diff_work)