6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
20 implicit none ;
private
22 #include <MOM_memory.h>
42 logical :: absorb_all_sw
49 real :: bulk_ri_convective
52 real :: h_limit_fluxes
67 real :: hbuffer_rel_min
69 real :: bl_detrain_time
77 real :: bl_split_rho_tol
82 integer :: ml_presort_nz_conv_adj
88 logical :: correct_absorption
92 logical :: resolve_ekman
95 type(time_type),
pointer :: time => null()
96 logical :: tke_diagnostics = .false.
97 logical :: do_rivermix = .false.
99 real :: rivermix_depth = 0.0
102 real :: lim_det_dh_sfc
105 real :: lim_det_dh_bathy
109 logical :: use_river_heat_content
112 logical :: use_calving_heat_content
113 logical :: salt_reject_below_ml
117 real :: allowed_t_chg
119 real :: allowed_s_chg
123 real,
allocatable,
dimension(:,:) :: &
124 ml_depth, & !< The mixed layer depth [H ~> m or kg m-2].
125 diag_tke_wind, & !< The wind source of TKE.
126 diag_tke_ribulk, & !< The resolved KE source of TKE.
127 diag_tke_conv, & !< The convective source of TKE.
128 diag_tke_pen_sw, & !< The TKE sink required to mix penetrating shortwave heating.
129 diag_tke_mech_decay, & !< The decay of mechanical TKE.
130 diag_tke_conv_decay, & !< The decay of convective TKE.
131 diag_tke_mixing, & !< The work done by TKE to deepen the mixed layer.
132 diag_tke_conv_s2, & !< The convective source of TKE due to to mixing in sigma2.
133 diag_pe_detrain, & !< The spurious source of potential energy due to mixed layer
137 logical :: allow_clocks_in_omp_loops
139 type(group_pass_type) :: pass_h_sum_hmbl_prev
142 integer :: id_ml_depth = -1, id_tke_wind = -1, id_tke_mixing = -1
143 integer :: id_tke_ribulk = -1, id_tke_conv = -1, id_tke_pen_sw = -1
144 integer :: id_tke_mech_decay = -1, id_tke_conv_decay = -1, id_tke_conv_s2 = -1
145 integer :: id_pe_detrain = -1, id_pe_detrain2 = -1, id_h_mismatch = -1
146 integer :: id_hsfc_used = -1, id_hsfc_max = -1, id_hsfc_min = -1
187 subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, CS, &
188 optics, Hml, aggregate_FW_forcing, dt_diag, last_call)
192 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
193 intent(inout) :: h_3d
194 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
197 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
203 type(
forcing),
intent(inout) :: fluxes
206 real,
intent(in) :: dt
207 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
211 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
220 real,
dimension(:,:),
pointer :: hml
221 logical,
intent(in) :: aggregate_fw_forcing
225 real,
optional,
intent(in) :: dt_diag
228 logical,
optional,
intent(in) :: last_call
234 real,
dimension(SZI_(G),SZK_(GV)) :: &
243 real,
dimension(SZI_(G),SZK0_(GV)) :: &
249 real,
dimension(SZI_(G),SZK_(GV)) :: &
260 integer,
dimension(SZI_(G),SZK_(GV)) :: &
262 real,
dimension(SZI_(G),SZJ_(G)) :: &
264 real,
dimension(SZI_(G)) :: &
306 real,
dimension(max(CS%nsw,1),SZI_(G)) :: &
309 real,
dimension(max(CS%nsw,1),SZI_(G),SZK_(GV)) :: &
312 real :: cmke(2,szi_(g))
316 real :: inkml, inkmlm1
321 real,
dimension(SZI_(G)) :: &
325 real,
dimension(SZI_(G),SZK_(GV)) :: &
330 real,
dimension(SZI_(G),SZJ_(G)) :: &
340 real,
dimension(SZI_(G)) :: &
351 logical :: write_diags
352 logical :: reset_diags
353 integer :: i, j, k, is, ie, js, je, nz, nkmb, n
356 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
358 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_mixed_layer: "//&
359 "Module must be initialized before it is used.")
360 if (gv%nkml < 1)
return
362 if (.not.
associated(tv%eqn_of_state))
call mom_error(fatal, &
363 "MOM_mixed_layer: Temperature, salinity and an equation of state "//&
365 if (.NOT.
associated(fluxes%ustar))
call mom_error(fatal, &
366 "MOM_mixed_layer: No surface TKE fluxes (ustar) defined in mixedlayer!")
368 nkmb = cs%nkml+cs%nkbl
369 inkml = 1.0 / real(cs%nkml)
370 if (cs%nkml > 1) inkmlm1 = 1.0 / real(cs%nkml-1)
372 irho0 = 1.0 / (gv%Rho0)
373 dt__diag = dt ;
if (
present(dt_diag)) dt__diag = dt_diag
374 idt_diag = 1.0 / (dt__diag)
375 write_diags = .true. ;
if (
present(last_call)) write_diags = last_call
377 p_ref(:) = 0.0 ; p_ref_cv(:) = tv%P_Ref
381 if (cs%limit_det .or. (cs%id_Hsfc_min > 0))
then
383 do j=js-1,je+1 ;
do i=is-1,ie+1
384 h_sum(i,j) = 0.0 ; hmbl_prev(i,j) = 0.0
388 do k=1,nkmb ;
do i=is-1,ie+1
389 h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
390 hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k)
392 do k=nkmb+1,nz ;
do i=is-1,ie+1
393 h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
406 if (
present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
407 reset_diags = .false.
409 if (reset_diags)
then
410 if (cs%TKE_diagnostics)
then
412 do j=js,je ;
do i=is,ie
413 cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_RiBulk(i,j) = 0.0
414 cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_pen_SW(i,j) = 0.0
415 cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
416 cs%diag_TKE_conv_decay(i,j) = 0.0 ; cs%diag_TKE_conv_s2(i,j) = 0.0
419 if (
allocated(cs%diag_PE_detrain))
then
421 do j=js,je ;
do i=is,ie
422 cs%diag_PE_detrain(i,j) = 0.0
425 if (
allocated(cs%diag_PE_detrain2))
then
427 do j=js,je ;
do i=is,ie
428 cs%diag_PE_detrain2(i,j) = 0.0
433 if (cs%ML_resort)
then
434 do i=is,ie ; h_ca(i) = 0.0 ;
enddo
435 do k=1,nz ;
do i=is,ie ; dke_ca(i,k) = 0.0 ; ctke(i,k) = 0.0 ;
enddo ;
enddo
449 do k=1,nz ;
do i=is,ie
450 h(i,k) = h_3d(i,j,k) ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
451 h_orig(i,k) = h_3d(i,j,k)
452 eps(i,k) = 0.0 ;
if (k > nkmb) eps(i,k) = gv%Angstrom_H
453 t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
455 if (nsw>0)
call extract_optics_slice(optics, j, g, gv, opacity=opacity_band, opacity_scale=gv%H_to_m)
457 do k=1,nz ;
do i=is,ie
458 d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0
463 do i=is,ie ; p_ref(i) = 0.0 ;
enddo
464 do k=1,cs%nkml ;
do i=is,ie
465 p_ref(i) = p_ref(i) + 0.5*gv%H_to_Pa*h(i,k)
468 is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
470 is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
473 tv%eqn_of_state, scale=us%kg_m3_to_R)
475 ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
479 if (cs%ML_resort)
then
481 if (cs%ML_presort_nz_conv_adj > 0) &
482 call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
483 s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, us, cs, &
484 cs%ML_presort_nz_conv_adj)
486 call sort_ml(h(:,1:), r0(:,1:), eps, g, gv, cs, ksort)
489 do k=1,nz ;
do i=is,ie ; ksort(i,k) = k ;
enddo ;
enddo
495 call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
496 s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, us, cs)
497 do i=is,ie ; h_ca(i) = h(i,1) ;
enddo
502 if (
associated(fluxes%lrunoff) .and. cs%do_rivermix)
then
514 rmixconst = 0.5*cs%rivermix_depth * gv%g_Earth * irho0**2
516 tke_river(i) = max(0.0, rmixconst*dr0_ds(i)* &
517 (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * s(i,1))
520 do i=is,ie ; tke_river(i) = 0.0 ;
enddo
534 cs%H_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
535 h(:,1:), t(:,1:), netmassinout, netmassout, net_heat, net_salt, pen_sw_bnd,&
536 tv, aggregate_fw_forcing)
540 r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), &
541 r0(:,1:), rcv(:,1:), eps, &
542 dr0_dt, drcv_dt, dr0_ds, drcv_ds, &
543 netmassinout, netmassout, net_heat, net_salt, &
544 nsw, pen_sw_bnd, opacity_band, conv_en, &
545 dke_fc, j, ksort, g, gv, us, cs, tv, fluxes, dt, &
546 aggregate_fw_forcing)
559 tke, tke_river, idecay_len_tke, cmke, dt, idt_diag, &
560 j, ksort, g, gv, us, cs)
564 r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), r0(:,1:), rcv(:,1:), eps, dr0_dt, drcv_dt, &
565 cmke, idt_diag, nsw, pen_sw_bnd, opacity_band, tke, &
566 idecay_len_tke, j, ksort, g, gv, us, cs)
568 call absorbremainingsw(g, gv, us, h(:,1:), opacity_band, nsw, optics, j, dt, &
569 cs%H_limit_fluxes, cs%correct_absorption, cs%absorb_all_SW, &
570 t(:,1:), pen_sw_bnd, eps, ksort, htot, ttot)
572 if (cs%TKE_diagnostics)
then ;
do i=is,ie
573 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) - idt_diag * tke(i)
578 do i=is,ie ;
if (htot(i) > 0.0)
then
580 r0(i,0) = r0_tot(i) * ih ; rcv(i,0) = rcv_tot(i) * ih
581 t(i,0) = ttot(i) * ih ; s(i,0) = stot(i) * ih
584 t(i,0) = t(i,1) ; s(i,0) = s(i,1) ; r0(i,0) = r0(i,1) ; rcv(i,0) = rcv(i,1)
587 if (write_diags .and.
allocated(cs%ML_depth))
then ;
do i=is,ie
588 cs%ML_depth(i,j) = h(i,0) * gv%H_to_m
590 if (
associated(hml))
then ;
do i=is,ie
591 hml(i,j) = g%mask2dT(i,j) * (h(i,0) * gv%H_to_m)
604 if (cs%ML_resort)
then
606 call resort_ml(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), gv%Rlay(:), eps, &
607 d_ea, d_eb, ksort, g, gv, cs, dr0_dt, dr0_ds, drcv_dt, drcv_ds)
611 if (cs%limit_det .or. (cs%id_Hsfc_max > 0) .or. (cs%id_Hsfc_min > 0))
then
612 do i=is,ie ; hsfc(i) = h(i,0) ;
enddo
613 do k=1,nkmb ;
do i=is,ie ; hsfc(i) = hsfc(i) + h(i,k) ;
enddo ;
enddo
615 if (cs%limit_det .or. (cs%id_Hsfc_min > 0))
then
616 dhsfc = cs%lim_det_dH_sfc ; dhd = cs%lim_det_dH_bathy
618 h_nbr = min(dhsfc*max(hmbl_prev(i-1,j), hmbl_prev(i+1,j), &
619 hmbl_prev(i,j-1), hmbl_prev(i,j+1)), &
620 max(hmbl_prev(i-1,j) - dhd*min(h_sum(i,j),h_sum(i-1,j)), &
621 hmbl_prev(i+1,j) - dhd*min(h_sum(i,j),h_sum(i+1,j)), &
622 hmbl_prev(i,j-1) - dhd*min(h_sum(i,j),h_sum(i,j-1)), &
623 hmbl_prev(i,j+1) - dhd*min(h_sum(i,j),h_sum(i,j+1))) )
625 hsfc_min(i,j) = gv%H_to_Z * max(h(i,0), min(hsfc(i), h_nbr))
627 if (cs%limit_det) max_bl_det(i) = max(0.0, hsfc(i)-h_nbr)
631 if (cs%id_Hsfc_max > 0)
then ;
do i=is,ie
632 hsfc_max(i,j) = gv%H_to_Z * hsfc(i)
640 if (cs%nkbl == 1)
then
641 call mixedlayer_detrain_1(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
642 gv%Rlay(:), dt, dt__diag, d_ea, d_eb, j, g, gv, us, cs, &
643 drcv_dt, drcv_ds, max_bl_det)
644 elseif (cs%nkbl == 2)
then
645 call mixedlayer_detrain_2(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
646 gv%Rlay(:), dt, dt__diag, d_ea, j, g, gv, us, cs, &
647 dr0_dt, dr0_ds, drcv_dt, drcv_ds, max_bl_det)
650 call mom_error(fatal,
"MOM_mixed_layer: CS%nkbl must be 1 or 2 for now.")
655 if (cs%id_Hsfc_used > 0)
then
656 do i=is,ie ; hsfc_used(i,j) = gv%H_to_Z * h(i,0) ;
enddo
657 do k=cs%nkml+1,nkmb ;
do i=is,ie
658 hsfc_used(i,j) = hsfc_used(i,j) + gv%H_to_Z * h(i,k)
664 if (cs%Resolve_Ekman .and. (cs%nkml>1))
then
673 ku_star = 0.41*fluxes%ustar(i,j)
674 if (
associated(fluxes%ustar_shelf) .and. &
675 associated(fluxes%frac_shelf_h))
then
676 if (fluxes%frac_shelf_h(i,j) > 0.0) &
677 ku_star = (1.0 - fluxes%frac_shelf_h(i,j)) * ku_star + &
678 fluxes%frac_shelf_h(i,j) * (0.41*fluxes%ustar_shelf(i,j))
680 absf_x_h = 0.25 * gv%H_to_Z * h(i,0) * &
681 ((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
682 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
685 h_3d(i,j,1) = h(i,0) / (3.0 + sqrt(absf_x_h*(absf_x_h + 2.0*ku_star) / ku_star**2))
688 h_3d(i,j,k) = (h(i,0)-h_3d(i,j,1)) * inkmlm1
689 d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
690 d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
695 h_3d(i,j,1) = h(i,0) * inkml
697 do k=2,cs%nkml ;
do i=is,ie
698 h_3d(i,j,k) = h(i,0) * inkml
699 d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
700 d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
703 do i=is,ie ; h(i,0) = 0.0 ;
enddo
704 do k=1,cs%nkml ;
do i=is,ie
705 tv%T(i,j,k) = t(i,0) ; tv%S(i,j,k) = s(i,0)
721 eaml(i,nz) = (h(i,nz) - h_orig(i,nz)) - d_eb(i,nz)
723 do k=nz-1,1,-1 ;
do i=is,ie
724 ebml(i,k) = ebml(i,k+1) - d_eb(i,k+1)
725 eaml(i,k) = eaml(i,k+1) + d_ea(i,k)
727 do i=is,ie ; eaml(i,1) = netmassinout(i) ;
enddo
730 do k=cs%nkml+1,nz ;
do i=is,ie
731 h_3d(i,j,k) = h(i,k); tv%T(i,j,k) = t(i,k) ; tv%S(i,j,k) = s(i,k)
734 do k=1,nz ;
do i=is,ie
735 ea(i,j,k) = ea(i,j,k) + eaml(i,k)
736 eb(i,j,k) = eb(i,j,k) + ebml(i,k)
739 if (cs%id_h_mismatch > 0)
then
741 h_miss(i,j) = gv%H_to_Z * abs(h_3d(i,j,1) - (h_orig(i,1) + &
742 (eaml(i,1) + (ebml(i,1) - eaml(i,1+1)))))
744 do k=2,nz-1 ;
do i=is,ie
745 h_miss(i,j) = h_miss(i,j) + gv%H_to_Z * abs(h_3d(i,j,k) - (h_orig(i,k) + &
746 ((eaml(i,k) - ebml(i,k-1)) + (ebml(i,k) - eaml(i,k+1)))))
749 h_miss(i,j) = h_miss(i,j) + gv%H_to_Z * abs(h_3d(i,j,nz) - (h_orig(i,nz) + &
750 ((eaml(i,nz) - ebml(i,nz-1)) + ebml(i,nz))))
763 if (write_diags)
then
764 if (cs%id_ML_depth > 0) &
765 call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
766 if (cs%id_TKE_wind > 0) &
767 call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
768 if (cs%id_TKE_RiBulk > 0) &
769 call post_data(cs%id_TKE_RiBulk, cs%diag_TKE_RiBulk, cs%diag)
770 if (cs%id_TKE_conv > 0) &
771 call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
772 if (cs%id_TKE_pen_SW > 0) &
773 call post_data(cs%id_TKE_pen_SW, cs%diag_TKE_pen_SW, cs%diag)
774 if (cs%id_TKE_mixing > 0) &
775 call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
776 if (cs%id_TKE_mech_decay > 0) &
777 call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
778 if (cs%id_TKE_conv_decay > 0) &
779 call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
780 if (cs%id_TKE_conv_s2 > 0) &
781 call post_data(cs%id_TKE_conv_s2, cs%diag_TKE_conv_s2, cs%diag)
782 if (cs%id_PE_detrain > 0) &
783 call post_data(cs%id_PE_detrain, cs%diag_PE_detrain, cs%diag)
784 if (cs%id_PE_detrain2 > 0) &
785 call post_data(cs%id_PE_detrain2, cs%diag_PE_detrain2, cs%diag)
786 if (cs%id_h_mismatch > 0) &
787 call post_data(cs%id_h_mismatch, h_miss, cs%diag)
788 if (cs%id_Hsfc_used > 0) &
789 call post_data(cs%id_Hsfc_used, hsfc_used, cs%diag)
790 if (cs%id_Hsfc_max > 0) &
791 call post_data(cs%id_Hsfc_max, hsfc_max, cs%diag)
792 if (cs%id_Hsfc_min > 0) &
793 call post_data(cs%id_Hsfc_min, hsfc_min, cs%diag)
802 dKE_CA, cTKE, j, G, GV, US, CS, nz_conv)
805 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: h
807 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: u
809 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: v
811 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: T
812 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: S
813 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: R0
815 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: Rcv
817 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_eb
821 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: eps
823 real,
dimension(SZI_(G),SZK_(GV)),
intent(out) :: dKE_CA
826 real,
dimension(SZI_(G),SZK_(GV)),
intent(out) :: cTKE
829 integer,
intent(in) :: j
832 integer,
optional,
intent(in) :: nz_conv
841 real,
dimension(SZI_(G)) :: &
862 integer :: is, ie, nz, i, k, k1, nzc, nkmb
864 is = g%isc ; ie = g%iec ; nz = gv%ke
865 g_h2_2rho0 = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0)
866 nzc = nz ;
if (
present(nz_conv)) nzc = nz_conv
867 nkmb = cs%nkml+cs%nkbl
872 do k1=min(nzc-1,nkmb),1,-1
874 h_orig_k1(i) = h(i,k1)
875 ke_orig(i) = 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)
876 uhtot(i) = h(i,k1)*u(i,k1) ; vhtot(i) = h(i,k1)*v(i,k1)
877 r0_tot(i) = r0(i,k1) * h(i,k1)
878 ctke(i,k1) = 0.0 ; dke_ca(i,k1) = 0.0
880 rcv_tot(i) = rcv(i,k1) * h(i,k1)
881 ttot(i) = t(i,k1) * h(i,k1) ; stot(i) = s(i,k1) * h(i,k1)
885 if ((h(i,k) > eps(i,k)) .and. (r0_tot(i) > h(i,k1)*r0(i,k)))
then
886 h_ent = h(i,k)-eps(i,k)
887 ctke(i,k1) = ctke(i,k1) + h_ent * g_h2_2rho0 * &
888 (r0_tot(i) - h(i,k1)*r0(i,k)) * cs%nstar2
890 ctke(i,k1) = ctke(i,k1) + ctke(i,k)
891 dke_ca(i,k1) = dke_ca(i,k1) + dke_ca(i,k)
893 r0_tot(i) = r0_tot(i) + h_ent * r0(i,k)
894 ke_orig(i) = ke_orig(i) + 0.5*h_ent* &
895 (u(i,k)*u(i,k) + v(i,k)*v(i,k))
896 uhtot(i) = uhtot(i) + h_ent*u(i,k)
897 vhtot(i) = vhtot(i) + h_ent*v(i,k)
899 rcv_tot(i) = rcv_tot(i) + h_ent * rcv(i,k)
900 ttot(i) = ttot(i) + h_ent * t(i,k)
901 stot(i) = stot(i) + h_ent * s(i,k)
902 h(i,k1) = h(i,k1) + h_ent ; h(i,k) = eps(i,k)
904 d_eb(i,k) = d_eb(i,k) - h_ent
905 d_eb(i,k1) = d_eb(i,k1) + h_ent
911 do i=is,ie ;
if (h(i,k1) > h_orig_k1(i))
then
913 r0(i,k1) = r0_tot(i) * ih
914 u(i,k1) = uhtot(i) * ih ; v(i,k1) = vhtot(i) * ih
915 dke_ca(i,k1) = dke_ca(i,k1) + gv%H_to_Z * (cs%bulk_Ri_convective * &
916 (ke_orig(i) - 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)))
917 rcv(i,k1) = rcv_tot(i) * ih
918 t(i,k1) = ttot(i) * ih ; s(i,k1) = stot(i) * ih
923 do k=2,min(nzc,nkmb) ;
do i=is,ie ;
if (h(i,k) == 0.0)
then
925 rcv(i,k) = rcv(i,k-1) ; t(i,k) = t(i,k-1) ; s(i,k) = s(i,k-1)
926 endif ;
enddo ;
enddo
934 R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, &
935 dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, &
936 netMassInOut, netMassOut, Net_heat, Net_salt, &
937 nsw, Pen_SW_bnd, opacity_band, Conv_En, &
938 dKE_FC, j, ksort, G, GV, US, CS, tv, fluxes, dt, &
939 aggregate_FW_forcing)
942 real,
dimension(SZI_(G),SZK_(GV)), &
945 real,
dimension(SZI_(G),SZK_(GV)), &
946 intent(inout) :: d_eb
949 real,
dimension(SZI_(G)),
intent(out) :: htot
950 real,
dimension(SZI_(G)),
intent(out) :: Ttot
952 real,
dimension(SZI_(G)),
intent(out) :: Stot
954 real,
dimension(SZI_(G)),
intent(out) :: uhtot
956 real,
dimension(SZI_(G)),
intent(out) :: vhtot
958 real,
dimension(SZI_(G)),
intent(out) :: R0_tot
960 real,
dimension(SZI_(G)),
intent(out) :: Rcv_tot
962 real,
dimension(SZI_(G),SZK_(GV)), &
964 real,
dimension(SZI_(G),SZK_(GV)), &
966 real,
dimension(SZI_(G),SZK_(GV)), &
968 real,
dimension(SZI_(G),SZK_(GV)), &
970 real,
dimension(SZI_(G),SZK_(GV)), &
973 real,
dimension(SZI_(G),SZK_(GV)), &
976 real,
dimension(SZI_(G),SZK_(GV)), &
979 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
981 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
983 real,
dimension(SZI_(G)),
intent(in) :: dR0_dS
985 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
987 real,
dimension(SZI_(G)),
intent(in) :: netMassInOut
990 real,
dimension(SZI_(G)),
intent(in) :: netMassOut
992 real,
dimension(SZI_(G)),
intent(in) :: Net_heat
995 real,
dimension(SZI_(G)),
intent(in) :: Net_salt
997 integer,
intent(in) :: nsw
999 real,
dimension(max(nsw,1),SZI_(G)),
intent(inout) :: Pen_SW_bnd
1002 real,
dimension(max(nsw,1),SZI_(G),SZK_(GV)),
intent(in) :: opacity_band
1004 real,
dimension(SZI_(G)),
intent(out) :: Conv_En
1006 real,
dimension(SZI_(G)),
intent(out) :: dKE_FC
1008 integer,
intent(in) :: j
1009 integer,
dimension(SZI_(G),SZK_(GV)), &
1016 type(
forcing),
intent(inout) :: fluxes
1019 real,
intent(in) :: dt
1020 logical,
intent(in) :: aggregate_FW_forcing
1030 real,
dimension(SZI_(G)) :: &
1035 real :: Pen_absorbed
1042 real :: En_fn, Frac, x1
1044 real :: dr_ent, dr_comp
1046 real :: h_min, h_max
1061 integer :: is, ie, nz, i, k, ks, itt, n
1062 real,
dimension(max(nsw,1)) :: &
1066 angstrom = gv%Angstrom_H
1067 c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
1068 g_h2_2rho0 = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0)
1070 is = g%isc ; ie = g%iec ; nz = gv%ke
1072 do i=is,ie ;
if (ksort(i,1) > 0)
then
1075 if (aggregate_fw_forcing)
then
1077 if (netmassinout(i) < 0.0) massoutrem(i) = -netmassinout(i)
1078 netmassin(i) = netmassinout(i) + massoutrem(i)
1080 massoutrem(i) = -netmassout(i)
1081 netmassin(i) = netmassinout(i) - netmassout(i)
1085 h_ent = max(min(angstrom,h(i,k)-eps(i,k)),0.0)
1086 htot(i) = h_ent + netmassin(i)
1087 h(i,k) = h(i,k) - h_ent
1088 d_eb(i,k) = d_eb(i,k) - h_ent
1091 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1092 sw_trans = exp(-htot(i)*opacity_band(n,i,k))
1093 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0-sw_trans)
1094 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1101 ttot(i) = (net_heat(i) + (netmassin(i) * t_precip + h_ent * t(i,k))) + &
1109 stot(i) = h_ent*s(i,k) + net_salt(i)
1110 uhtot(i) = u(i,1)*netmassin(i) + u(i,k)*h_ent
1111 vhtot(i) = v(i,1)*netmassin(i) + v(i,k)*h_ent
1112 r0_tot(i) = (h_ent*r0(i,k) + netmassin(i)*r0(i,1)) + &
1114 (dr0_dt(i)*(net_heat(i) + pen_absorbed) - &
1115 dr0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1116 rcv_tot(i) = (h_ent*rcv(i,k) + netmassin(i)*rcv(i,1)) + &
1118 (drcv_dt(i)*(net_heat(i) + pen_absorbed) - &
1119 drcv_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1120 conv_en(i) = 0.0 ; dke_fc(i) = 0.0
1121 if (
associated(fluxes%heat_content_massin)) &
1122 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1123 t_precip * netmassin(i) * gv%H_to_RZ * fluxes%C_p * idt
1124 if (
associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1125 t_precip * netmassin(i) * gv%H_to_RZ
1131 do i=is,ie ;
if (ksort(i,ks) > 0)
then
1134 if ((htot(i) < angstrom) .and. (h(i,k) > eps(i,k)))
then
1137 h_ent = min(angstrom-htot(i), h(i,k)-eps(i,k))
1138 htot(i) = htot(i) + h_ent
1139 h(i,k) = h(i,k) - h_ent
1140 d_eb(i,k) = d_eb(i,k) - h_ent
1142 r0_tot(i) = r0_tot(i) + h_ent*r0(i,k)
1143 uhtot(i) = uhtot(i) + h_ent*u(i,k)
1144 vhtot(i) = vhtot(i) + h_ent*v(i,k)
1146 rcv_tot(i) = rcv_tot(i) + h_ent*rcv(i,k)
1147 ttot(i) = ttot(i) + h_ent*t(i,k)
1148 stot(i) = stot(i) + h_ent*s(i,k)
1154 if ((massoutrem(i) > 0.0) .and. (h(i,k) > eps(i,k)))
then
1155 if (massoutrem(i) > (h(i,k) - eps(i,k)))
then
1156 h_evap = h(i,k) - eps(i,k)
1158 massoutrem(i) = massoutrem(i) - h_evap
1160 h_evap = massoutrem(i)
1161 h(i,k) = h(i,k) - h_evap
1165 stot(i) = stot(i) + h_evap*s(i,k)
1166 r0_tot(i) = r0_tot(i) + dr0_ds(i)*h_evap*s(i,k)
1167 rcv_tot(i) = rcv_tot(i) + drcv_ds(i)*h_evap*s(i,k)
1168 d_eb(i,k) = d_eb(i,k) - h_evap
1174 if (
associated(fluxes%heat_content_massout)) &
1175 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) - &
1176 t(i,k)*h_evap*gv%H_to_RZ * fluxes%C_p * idt
1177 if (
associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
1178 t(i,k)*h_evap*gv%H_to_RZ
1183 h_avail = h(i,k) - eps(i,k)
1184 if (h_avail > 0.0)
then
1185 dr = r0_tot(i) - htot(i)*r0(i,k)
1189 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1190 dr0 = dr0 - (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1191 opacity_band(n,i,k)*htot(i)
1197 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1200 opacity = opacity_band(n,i,k)
1201 sw_trans = exp(-h_avail*opacity)
1202 dr_comp = dr_comp + (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1203 ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
1205 if (dr_comp >= 0.0)
then
1215 frac = dr0 / (dr0 - dr_comp)
1216 h_ent = h_avail * frac*frac
1217 h_min = 0.0 ; h_max = h_avail
1220 r_sw_top(n) = dr0_dt(i) * pen_sw_bnd(n,i)
1221 c2(n) = r_sw_top(n) * opacity_band(n,i,k)**2
1224 dr_ent = dr ; dr_dh = 0.0
1226 opacity = opacity_band(n,i,k)
1227 sw_trans = exp(-h_ent*opacity)
1228 dr_ent = dr_ent + r_sw_top(n) * ((1.0 - sw_trans) - &
1229 opacity*(htot(i)+h_ent)*sw_trans)
1230 dr_dh = dr_dh + c2(n) * (htot(i)+h_ent) * sw_trans
1233 if (dr_ent > 0.0)
then
1239 dh_newt = -dr_ent / dr_dh
1240 h_prev = h_ent ; h_ent = h_prev+dh_newt
1241 if (h_ent > h_max)
then
1242 h_ent = 0.5*(h_prev+h_max)
1243 elseif (h_ent < h_min)
then
1244 h_ent = 0.5*(h_prev+h_min)
1247 if (abs(dh_newt) < 0.2*angstrom)
exit
1254 sum_pen_en = 0.0 ; pen_absorbed = 0.0
1255 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1256 opacity = opacity_band(n,i,k)
1257 sw_trans = exp(-h_ent*opacity)
1260 if (x1 < 2.0e-5)
then
1261 en_fn = (opacity*htot(i)*(1.0 - 0.5*(x1 - c1_3*x1)) + &
1264 en_fn = ((opacity*htot(i) + 2.0) * &
1265 ((1.0-sw_trans) / x1) - 1.0 + sw_trans)
1267 sum_pen_en = sum_pen_en - (dr0_dt(i)*pen_sw_bnd(n,i)) * en_fn
1269 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1270 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1273 conv_en(i) = conv_en(i) + g_h2_2rho0 * h_ent * &
1274 ( (r0_tot(i) - r0(i,k)*htot(i)) + sum_pen_en )
1276 r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1277 stot(i) = stot(i) + h_ent * s(i,k)
1278 ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1279 rcv_tot(i) = rcv_tot(i) + (h_ent * rcv(i,k) + pen_absorbed*drcv_dt(i))
1282 if (h_ent > 0.0)
then
1283 if (htot(i) > 0.0) &
1284 dke_fc(i) = dke_fc(i) + cs%bulk_Ri_convective * 0.5 * &
1285 ((gv%H_to_Z*h_ent) / (htot(i)*(h_ent+htot(i)))) * &
1286 ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1288 htot(i) = htot(i) + h_ent
1289 h(i,k) = h(i,k) - h_ent
1290 d_eb(i,k) = d_eb(i,k) - h_ent
1291 uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
1305 subroutine find_starting_tke(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, &
1306 TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, &
1307 j, ksort, G, GV, US, CS)
1311 real,
dimension(SZI_(G)),
intent(in) :: htot
1313 real,
dimension(SZI_(G)),
intent(in) :: h_CA
1315 type(
forcing),
intent(in) :: fluxes
1318 real,
dimension(SZI_(G)),
intent(inout) :: Conv_En
1320 real,
dimension(SZI_(G)),
intent(in) :: dKE_FC
1323 real,
dimension(SZI_(G),SZK_(GV)), &
1327 real,
dimension(SZI_(G),SZK_(GV)), &
1328 intent(in) :: dKE_CA
1331 real,
dimension(SZI_(G)),
intent(out) :: TKE
1333 real,
dimension(SZI_(G)),
intent(out) :: Idecay_len_TKE
1335 real,
dimension(SZI_(G)),
intent(in) :: TKE_river
1338 real,
dimension(2,SZI_(G)),
intent(out) :: cMKE
1341 real,
intent(in) :: dt
1342 real,
intent(in) :: Idt_diag
1344 integer,
intent(in) :: j
1345 integer,
dimension(SZI_(G),SZK_(GV)), &
1368 real :: wind_TKE_src
1371 integer :: is, ie, nz, i
1373 is = g%isc ; ie = g%iec ; nz = gv%ke
1374 diag_wt = dt * idt_diag
1376 if (cs%omega_frac >= 1.0) absf = 2.0*cs%omega
1378 u_star = fluxes%ustar(i,j)
1379 if (
associated(fluxes%ustar_shelf) .and.
associated(fluxes%frac_shelf_h))
then
1380 if (fluxes%frac_shelf_h(i,j) > 0.0) &
1381 u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
1382 fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
1385 if (u_star < cs%ustar_min) u_star = cs%ustar_min
1386 if (cs%omega_frac < 1.0)
then
1387 absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
1388 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
1389 if (cs%omega_frac > 0.0) &
1390 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1392 absf_ustar = absf / u_star
1393 idecay_len_tke(i) = (absf_ustar * cs%TKE_decay) * gv%H_to_Z
1405 ih = gv%H_to_Z/(3.0*0.41*u_star*dt)
1406 cmke(1,i) = 4.0 * ih ; cmke(2,i) = (absf_ustar*gv%H_to_Z) * ih
1408 if (idecay_len_tke(i) > 0.0)
then
1409 exp_kh = exp(-htot(i)*idecay_len_tke(i))
1417 if (conv_en(i) < 0.0) conv_en(i) = 0.0
1418 if (ctke(i,1) > 0.0)
then ; tke_ca = ctke(i,1) ;
else ; tke_ca = 0.0 ;
endif
1419 if ((htot(i) >= h_ca(i)) .or. (tke_ca == 0.0))
then
1420 toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca)
1422 if (toten_z > 0.0)
then
1423 nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1424 sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
1431 if (conv_en(i) > 0.0)
then
1432 toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca * (htot(i) / h_ca(i)) )
1433 nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1434 sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
1439 toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca)
1440 if (tke_ca > 0.0)
then
1441 nstar_ca = cs%nstar * toten_z / (toten_z + 0.2 * &
1442 sqrt(0.5 * dt * (absf*(h_ca(i)*gv%H_to_Z))**3 * toten_z))
1448 if (dke_fc(i) + dke_ca(i,1) > 0.0)
then
1449 if (htot(i) >= h_ca(i))
then
1450 mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1451 mke_rate_ca = mke_rate_fc
1453 mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1454 mke_rate_ca = 1.0 / (1.0 + h_ca(i)*(cmke(1,i) + cmke(2,i)*h_ca(i)) )
1458 mke_rate_fc = 1.0 ; mke_rate_ca = 1.0
1461 dke_conv = dke_ca(i,1) * mke_rate_ca + dke_fc(i) * mke_rate_fc
1464 tke(i) = (dt*cs%mstar)*((us%Z_to_L**2*(u_star*u_star*u_star))*exp_kh) + &
1465 (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca * tke_ca)
1467 if (cs%do_rivermix)
then
1468 tke(i) = tke(i) + tke_river(i)*dt*exp_kh
1471 if (cs%TKE_diagnostics)
then
1472 wind_tke_src = cs%mstar*(us%Z_to_L**2*u_star*u_star*u_star) * diag_wt
1473 cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + &
1474 ( wind_tke_src + tke_river(i) * diag_wt )
1475 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + dke_conv*idt_diag
1476 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1477 (exp_kh-1.0)*(wind_tke_src + dke_conv*idt_diag)
1478 cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + &
1479 idt_diag * (nstar_fc*conv_en(i) + nstar_ca*tke_ca)
1480 cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + &
1481 idt_diag * ((cs%nstar-nstar_fc)*conv_en(i) + (cs%nstar-nstar_ca)*tke_ca)
1482 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
1483 idt_diag * (ctke(i,1)-tke_ca)
1491 R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, &
1492 dR0_dT, dRcv_dT, cMKE, Idt_diag, nsw, &
1493 Pen_SW_bnd, opacity_band, TKE, &
1494 Idecay_len_TKE, j, ksort, G, GV, US, CS)
1498 real,
dimension(SZI_(G),SZK_(GV)), &
1500 real,
dimension(SZI_(G),SZK_(GV)), &
1501 intent(inout) :: d_eb
1504 real,
dimension(SZI_(G)),
intent(inout) :: htot
1505 real,
dimension(SZI_(G)),
intent(inout) :: Ttot
1507 real,
dimension(SZI_(G)),
intent(inout) :: Stot
1509 real,
dimension(SZI_(G)),
intent(inout) :: uhtot
1511 real,
dimension(SZI_(G)),
intent(inout) :: vhtot
1513 real,
dimension(SZI_(G)),
intent(inout) :: R0_tot
1515 real,
dimension(SZI_(G)),
intent(inout) :: Rcv_tot
1517 real,
dimension(SZI_(G),SZK_(GV)), &
1519 real,
dimension(SZI_(G),SZK_(GV)), &
1521 real,
dimension(SZI_(G),SZK_(GV)), &
1523 real,
dimension(SZI_(G),SZK_(GV)), &
1525 real,
dimension(SZI_(G),SZK_(GV)), &
1528 real,
dimension(SZI_(G),SZK_(GV)), &
1531 real,
dimension(SZI_(G),SZK_(GV)), &
1534 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
1536 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
1538 real,
dimension(2,SZI_(G)),
intent(in) :: cMKE
1541 real,
intent(in) :: Idt_diag
1543 integer,
intent(in) :: nsw
1545 real,
dimension(max(nsw,1),SZI_(G)),
intent(inout) :: Pen_SW_bnd
1548 real,
dimension(max(nsw,1),SZI_(G),SZK_(GV)),
intent(in) :: opacity_band
1550 real,
dimension(SZI_(G)),
intent(inout) :: TKE
1553 real,
dimension(SZI_(G)),
intent(inout) :: Idecay_len_TKE
1554 integer,
intent(in) :: j
1555 integer,
dimension(SZI_(G),SZK_(GV)), &
1564 real :: Pen_absorbed
1568 real :: h_min, h_max
1578 real :: TKE_full_ent
1582 real :: Pen_En_Contrib
1591 real :: Pen_dTKE_dh_Contrib
1601 real :: f1_x1, f2_x1
1607 real :: C1_3, C1_6, C1_24
1608 integer :: is, ie, nz, i, k, ks, itt, n
1610 c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0 ; c1_24 = 1.0/24.0
1611 g_h_2rho0 = (gv%g_Earth * gv%H_to_Z) / (2.0 * gv%Rho0)
1612 hmix_min = cs%Hmix_min
1613 h_neglect = gv%H_subroundoff
1614 is = g%isc ; ie = g%iec ; nz = gv%ke
1618 do i=is,ie ;
if (ksort(i,ks) > 0)
then
1621 h_avail = h(i,k) - eps(i,k)
1622 if ((h_avail > 0.) .and. ((tke(i) > 0.) .or. (htot(i) < hmix_min)))
then
1623 drl = g_h_2rho0 * (r0(i,k)*htot(i) - r0_tot(i) )
1624 dmke = (gv%H_to_Z * cs%bulk_Ri_ML) * 0.5 * &
1625 ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1628 kh = idecay_len_tke(i)*h_avail ; exp_kh = exp(-kh)
1629 if (kh >= 2.0e-5)
then ; f1_kh = (1.0-exp_kh) / kh
1630 else ; f1_kh = (1.0 - kh*(0.5 - c1_6*kh)) ;
endif
1632 pen_en_contrib = 0.0
1633 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1634 opacity = opacity_band(n,i,k)
1638 if (idecay_len_tke(i) > opacity)
then
1639 x1 = (idecay_len_tke(i) - opacity) * h_avail
1640 if (x1 >= 2.0e-5)
then
1641 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1642 f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1644 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1645 f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1648 pen_en1 = exp(-opacity*h_avail) * &
1649 ((1.0+opacity*htot(i))*f1_x1 + opacity*h_avail*f3_x1)
1651 x1 = (opacity - idecay_len_tke(i)) * h_avail
1652 if (x1 >= 2.0e-5)
then
1653 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1654 f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1656 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1657 f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1660 pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1661 opacity*h_avail*f2_x1)
1663 pen_en_contrib = pen_en_contrib + &
1664 (g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
1667 hpe = htot(i)+h_avail
1668 mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1669 ef4_val =
ef4(htot(i)+h_neglect,h_avail,idecay_len_tke(i))
1670 tke_full_ent = (exp_kh*tke(i) - (h_avail*gv%H_to_Z)*(drl*f1_kh + pen_en_contrib)) + &
1671 mke_rate*dmke*ef4_val
1672 if ((tke_full_ent >= 0.0) .or. (h_avail+htot(i) <= hmix_min))
then
1676 if (cs%TKE_diagnostics)
then
1677 e_hxhpe = h_ent / ((htot(i)+h_neglect)*(htot(i)+h_ent+h_neglect))
1678 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1679 idt_diag * ((exp_kh-1.0)* tke(i) + (h_ent*gv%H_to_Z)*drl*(1.0-f1_kh) + &
1680 mke_rate*dmke*(ef4_val-e_hxhpe))
1681 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1682 idt_diag*(gv%H_to_Z*h_ent)*drl
1683 cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1684 idt_diag*(gv%H_to_Z*h_ent)*pen_en_contrib
1685 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1686 idt_diag*mke_rate*dmke*e_hxhpe
1689 tke(i) = tke_full_ent
1691 if (tke(i) <= 0.0) tke(i) = 1.0e-150*us%m_to_Z*us%m_s_to_L_T**2
1696 h_min = 0.0 ; h_max = h_avail
1697 if (tke(i) <= 0.0)
then
1700 h_ent = h_avail * tke(i) / (tke(i) - tke_full_ent)
1705 kh = idecay_len_tke(i)*h_ent ; exp_kh = exp(-kh)
1706 if (kh >= 2.0e-5)
then
1707 f1_kh = (1.0-exp_kh) / kh
1709 f1_kh = (1.0 - kh*(0.5 - c1_6*kh))
1713 pen_en_contrib = 0.0 ; pen_dtke_dh_contrib = 0.0
1714 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1718 opacity = opacity_band(n,i,k)
1719 sw_trans = exp(-h_ent*opacity)
1720 if (idecay_len_tke(i) > opacity)
then
1721 x1 = (idecay_len_tke(i) - opacity) * h_ent
1722 if (x1 >= 2.0e-5)
then
1723 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1724 f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1726 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1727 f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1729 pen_en1 = sw_trans * ((1.0+opacity*htot(i))*f1_x1 + &
1730 opacity*h_ent*f3_x1)
1732 x1 = (opacity - idecay_len_tke(i)) * h_ent
1733 if (x1 >= 2.0e-5)
then
1734 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1735 f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1737 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1738 f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1741 pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1742 opacity*h_ent*f2_x1)
1744 cpen1 = g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)
1745 pen_en_contrib = pen_en_contrib + cpen1*(pen_en1 - f1_kh)
1746 pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
1747 cpen1*((1.0-sw_trans) - opacity*(htot(i) + h_ent)*sw_trans)
1750 tke_ent1 = exp_kh* tke(i) - (h_ent*gv%H_to_Z)*(drl*f1_kh + pen_en_contrib)
1751 ef4_val =
ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i),def4_dh)
1753 mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1754 tke_ent = tke_ent1 + dmke*ef4_val*mke_rate
1757 dtke_dh = ((-idecay_len_tke(i)*tke_ent1 - drl*gv%H_to_Z) + &
1758 pen_dtke_dh_contrib*gv%H_to_Z) + dmke * mke_rate* &
1759 (def4_dh - ef4_val*mke_rate*(cmke(1,i)+2.0*cmke(2,i)*hpe))
1762 if (tke_ent > 0.0)
then
1763 if ((h_max-h_ent)*(-dtke_dh) > tke_ent)
then
1764 dh_newt = -tke_ent / dtke_dh
1766 dh_newt = 0.5*(h_max-h_ent)
1770 if ((h_min-h_ent)*(-dtke_dh) < tke_ent)
then
1771 dh_newt = -tke_ent / dtke_dh
1773 dh_newt = 0.5*(h_min-h_ent)
1777 h_ent = h_ent + dh_newt
1779 if (abs(dh_newt) < 0.2*gv%Angstrom_H)
exit
1783 if (h_ent < hmix_min-htot(i)) h_ent = hmix_min - htot(i)
1785 if (cs%TKE_diagnostics)
then
1787 mke_rate = 1.0/(1.0 + cmke(1,i)*hpe + cmke(2,i)*hpe**2)
1788 ef4_val =
ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i))
1790 e_hxhpe = h_ent / ((htot(i)+h_neglect)*(hpe+h_neglect))
1791 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1792 idt_diag * ((exp_kh-1.0)* tke(i) + (h_ent*gv%H_to_Z)*drl*(1.0-f1_kh) + &
1793 dmke*mke_rate*(ef4_val-e_hxhpe))
1794 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1795 idt_diag*(h_ent*gv%H_to_Z)*drl
1796 cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1797 idt_diag*(h_ent*gv%H_to_Z)*pen_en_contrib
1798 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1799 idt_diag*dmke*mke_rate*e_hxhpe
1806 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1807 sw_trans = exp(-h_ent*opacity_band(n,i,k))
1808 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1809 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1812 htot(i) = htot(i) + h_ent
1813 r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1814 h(i,k) = h(i,k) - h_ent
1815 d_eb(i,k) = d_eb(i,k) - h_ent
1817 stot(i) = stot(i) + h_ent * s(i,k)
1818 ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1819 rcv_tot(i) = rcv_tot(i) + (h_ent*rcv(i,k) + pen_absorbed*drcv_dt(i))
1821 uhtot(i) = uhtot(i) + u(i,k)*h_ent
1822 vhtot(i) = vhtot(i) + v(i,k)*h_ent
1832 subroutine sort_ml(h, R0, eps, G, GV, CS, ksort)
1835 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: h
1836 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: R0
1838 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: eps
1842 integer,
dimension(SZI_(G),SZK_(GV)),
intent(out) :: ksort
1845 real :: R0sort(SZI_(G),SZK_(GV))
1846 integer :: nsort(SZI_(G))
1847 logical :: done_sorting(SZI_(G))
1848 integer :: i, k, ks, is, ie, nz, nkmb
1850 is = g%isc ; ie = g%iec ; nz = gv%ke
1851 nkmb = cs%nkml+cs%nkbl
1861 do k=1,nz ;
do i=is,ie ; ksort(i,k) = -1 ;
enddo ;
enddo
1863 do i=is,ie ; nsort(i) = 0 ; done_sorting(i) = .false. ;
enddo
1864 do k=1,nz ;
do i=is,ie ;
if (h(i,k) > eps(i,k))
then
1865 if (done_sorting(i))
then ; ks = nsort(i) ;
else
1867 if (r0(i,k) >= r0sort(i,ks))
exit
1868 r0sort(i,ks+1) = r0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
1870 if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
1874 r0sort(i,ks+1) = r0(i,k)
1875 nsort(i) = nsort(i) + 1
1876 endif ;
enddo ;
enddo
1883 subroutine resort_ml(h, T, S, R0, Rcv, RcvTgt, eps, d_ea, d_eb, ksort, G, GV, CS, &
1884 dR0_dT, dR0_dS, dRcv_dT, dRcv_dS)
1888 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: h
1890 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: T
1891 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: S
1892 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: R0
1894 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: Rcv
1896 real,
dimension(SZK_(GV)),
intent(in) :: RcvTgt
1898 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: eps
1900 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_ea
1905 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_eb
1909 integer,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: ksort
1912 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
1916 real,
dimension(SZI_(G)),
intent(in) :: dR0_dS
1920 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
1924 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
1943 real :: h_move, h_tgt_old, I_hnew
1944 real :: dT_dS_wt2, dT_dR, dS_dR, I_denom
1946 real :: T_up, S_up, R0_up, I_hup, h_to_up
1947 real :: T_dn, S_dn, R0_dn, I_hdn, h_to_dn
1950 real :: dPE, hmin, min_dPE, min_hmin
1951 real,
dimension(SZK_(GV)) :: &
1952 h_tmp, R0_tmp, T_tmp, S_tmp, Rcv_tmp
1954 logical :: sorted, leave_in_layer
1955 integer :: ks_deep(SZI_(G)), k_count(SZI_(G)), ks2_reverse(SZI_(G), SZK_(GV))
1956 integer :: ks2(SZK_(GV))
1957 integer :: i, k, ks, is, ie, nz, k1, k2, k_tgt, k_src, k_int_top
1958 integer :: nks, nkmb, num_interior, top_interior_ks
1960 is = g%isc ; ie = g%iec ; nz = gv%ke
1961 nkmb = cs%nkml+cs%nkbl
1963 dt_ds_wt2 = cs%dT_dS_wt**2
1966 do i=is,ie ; ks_deep(i) = -1 ; k_count(i) = 0 ;
enddo
1967 do ks=nz,1,-1 ;
do i=is,ie ;
if (ksort(i,ks) > 0)
then
1970 if (h(i,k) > eps(i,k))
then
1971 if (ks_deep(i) == -1)
then
1973 ks_deep(i) = ks ; k_count(i) = k_count(i) + 1
1974 ks2_reverse(i,k_count(i)) = k
1977 k_count(i) = k_count(i) + 1
1978 ks2_reverse(i,k_count(i)) = k
1981 endif ;
enddo ;
enddo
1983 do i=is,ie ;
if (k_count(i) > 1)
then
1989 ks2(nks) = ks2_reverse(i,1)
1991 ks2(ks) = ks2_reverse(i,1+nks-ks)
1992 if (ks2(ks) > ks2(ks+1)) sorted = .false.
2000 num_interior = 0 ; top_interior_ks = 0
2001 do ks=nks,1,-1 ;
if (ks2(ks) > nkmb)
then
2002 num_interior = num_interior+1 ; top_interior_ks = ks
2005 if (num_interior >= 1)
then
2008 do k=nkmb+1,nz ;
if (rcv(i,0) < rcvtgt(k))
exit ;
enddo
2009 k_int_top = k ; rcv_int = rcvtgt(k)
2011 i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
2012 dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
2013 ds_dr = drcv_ds(i) * i_denom
2019 do ks = nks,top_interior_ks,-1
2021 leave_in_layer = .false.
2022 if ((k > nkmb) .and. (rcv(i,k) <= rcvtgt(k)))
then
2023 if (rcvtgt(k)-rcv(i,k) < cs%BL_split_rho_tol*(rcvtgt(k) - rcvtgt(k-1))) &
2024 leave_in_layer = .true.
2025 elseif (k > nkmb)
then
2026 if (rcv(i,k)-rcvtgt(k) < cs%BL_split_rho_tol*(rcvtgt(k+1) - rcvtgt(k))) &
2027 leave_in_layer = .true.
2030 if (leave_in_layer)
then
2033 elseif (rcv(i,k) < rcv_int)
then
2040 do k2=k_int_top+1,nz ;
if (rcv(i,k) < rcvtgt(k2))
exit ;
enddo
2045 dr1 = (rcvtgt(k2-1) - rcv(i,k)) ; dr2 = (rcvtgt(k2) - rcv(i,k))
2046 t_up = t(i,k) + dt_dr * dr1
2047 s_up = s(i,k) + ds_dr * dr1
2048 t_dn = t(i,k) + dt_dr * dr2
2049 s_dn = s(i,k) + ds_dr * dr2
2051 r0_up = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr1
2052 r0_dn = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr2
2055 if ((r0_up > r0(i,0)) .or. (r0_dn > r0(i,0))) &
2059 wt_dn = (rcv(i,k) - rcvtgt(k2-1)) / (rcvtgt(k2) - rcvtgt(k2-1))
2060 h_to_up = (h(i,k)-eps(i,k)) * (1.0 - wt_dn)
2061 h_to_dn = (h(i,k)-eps(i,k)) * wt_dn
2063 i_hup = 1.0 / (h(i,k2-1) + h_to_up)
2064 i_hdn = 1.0 / (h(i,k2) + h_to_dn)
2065 r0(i,k2-1) = (r0(i,k2)*h(i,k2-1) + r0_up*h_to_up) * i_hup
2066 r0(i,k2) = (r0(i,k2)*h(i,k2) + r0_dn*h_to_dn) * i_hdn
2068 t(i,k2-1) = (t(i,k2)*h(i,k2-1) + t_up*h_to_up) * i_hup
2069 t(i,k2) = (t(i,k2)*h(i,k2) + t_dn*h_to_dn) * i_hdn
2070 s(i,k2-1) = (s(i,k2)*h(i,k2-1) + s_up*h_to_up) * i_hup
2071 s(i,k2) = (s(i,k2)*h(i,k2) + s_dn*h_to_dn) * i_hdn
2072 rcv(i,k2-1) = (rcv(i,k2)*h(i,k2-1) + rcvtgt(k2-1)*h_to_up) * i_hup
2073 rcv(i,k2) = (rcv(i,k2)*h(i,k2) + rcvtgt(k2)*h_to_dn) * i_hdn
2076 h(i,k2) = h(i,k2) + h_to_dn
2077 h(i,k2-1) = h(i,k2-1) + h_to_up
2080 d_eb(i,k) = d_eb(i,k) - h_to_up
2081 d_eb(i,k2-1) = d_eb(i,k2-1) + h_to_up
2082 elseif (k < k2-1)
then
2083 d_ea(i,k) = d_ea(i,k) - h_to_up
2084 d_ea(i,k2-1) = d_ea(i,k2-1) + h_to_up
2087 d_eb(i,k) = d_eb(i,k) - h_to_dn
2088 d_eb(i,k2) = d_eb(i,k2) + h_to_dn
2089 elseif (k < k2)
then
2090 d_ea(i,k) = d_ea(i,k) - h_to_dn
2091 d_ea(i,k2) = d_ea(i,k2) + h_to_dn
2098 do while (nks > nkmb)
2107 ks_min = -1 ; min_dpe = 1.0 ; min_hmin = 0.0
2109 k1 = ks2(ks) ; k2 = ks2(ks+1)
2110 dpe = max(0.0, (r0(i,k2)-r0(i,k1)) * h(i,k1) * h(i,k2))
2111 hmin = min(h(i,k1)-eps(i,k1), h(i,k2)-eps(i,k2))
2112 if ((ks_min < 0) .or. (dpe < min_dpe) .or. &
2113 ((dpe <= 0.0) .and. (hmin < min_hmin)))
then
2114 ks_min = ks ; min_dpe = dpe ; min_hmin = hmin
2119 k1 = ks2(ks_min) ; k2 = ks2(ks_min+1)
2120 if (k1 < k2)
then ; k_tgt = k1 ; k_src = k2
2121 else ; k_tgt = k2 ; k_src = k1 ; ks2(ks_min) = k2 ;
endif
2123 h_tgt_old = h(i,k_tgt)
2124 h_move = h(i,k_src)-eps(i,k_src)
2125 h(i,k_src) = eps(i,k_src)
2126 h(i,k_tgt) = h(i,k_tgt) + h_move
2127 i_hnew = 1.0 / (h(i,k_tgt))
2128 r0(i,k_tgt) = (r0(i,k_tgt)*h_tgt_old + r0(i,k_src)*h_move) * i_hnew
2130 t(i,k_tgt) = (t(i,k_tgt)*h_tgt_old + t(i,k_src)*h_move) * i_hnew
2131 s(i,k_tgt) = (s(i,k_tgt)*h_tgt_old + s(i,k_src)*h_move) * i_hnew
2132 rcv(i,k_tgt) = (rcv(i,k_tgt)*h_tgt_old + rcv(i,k_src)*h_move) * i_hnew
2134 d_eb(i,k_src) = d_eb(i,k_src) - h_move
2135 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2138 do ks=ks_min+1,nks ; ks2(ks) = ks2(ks+1) ;
enddo
2145 do ks=1,nks-1 ;
if (ks2(ks) > ks2(ks+1)) sorted = .false. ;
enddo
2154 h_tmp(k) = h(i,k) ; r0_tmp(k) = r0(i,k)
2155 t_tmp(k) = t(i,k) ; s_tmp(k) = s(i,k) ; rcv_tmp(k) = rcv(i,k)
2161 k_tgt = nkmb - nks + ks ; k_src = ks2(ks)
2162 if (k_tgt == k_src)
then
2163 h(i,k_tgt) = h_tmp(k_tgt)
2168 if (k_src > nkmb)
then
2169 h_move = h(i,k_src)-eps(i,k_src)
2170 h(i,k_src) = eps(i,k_src)
2172 r0(i,k_tgt) = r0(i,k_src)
2174 t(i,k_tgt) = t(i,k_src) ; s(i,k_tgt) = s(i,k_src)
2175 rcv(i,k_tgt) = rcv(i,k_src)
2177 d_eb(i,k_src) = d_eb(i,k_src) - h_move
2178 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2180 h(i,k_tgt) = h_tmp(k_src)
2181 r0(i,k_tgt) = r0_tmp(k_src)
2183 t(i,k_tgt) = t_tmp(k_src) ; s(i,k_tgt) = s_tmp(k_src)
2184 rcv(i,k_tgt) = rcv_tmp(k_src)
2186 if (k_src > k_tgt)
then
2187 d_eb(i,k_src) = d_eb(i,k_src) - h_tmp(k_src)
2188 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_tmp(k_src)
2190 d_ea(i,k_src) = d_ea(i,k_src) - h_tmp(k_src)
2191 d_ea(i,k_tgt) = d_ea(i,k_tgt) + h_tmp(k_src)
2204 subroutine mixedlayer_detrain_2(h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, j, G, GV, US, CS, &
2205 dR0_dT, dR0_dS, dRcv_dT, dRcv_dS, max_BL_det)
2208 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: h
2210 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: T
2211 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: S
2212 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: R0
2214 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: Rcv
2216 real,
dimension(SZK_(GV)),
intent(in) :: RcvTgt
2218 real,
intent(in) :: dt
2219 real,
intent(in) :: dt_diag
2220 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_ea
2224 integer,
intent(in) :: j
2228 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
2232 real,
dimension(SZI_(G)),
intent(in) :: dR0_dS
2236 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
2240 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
2243 real,
dimension(SZI_(G)),
intent(in) :: max_BL_det
2271 real :: stays_min, stays_max
2274 logical :: mergeable_bl
2280 real :: stays_min_merge
2282 real :: dR0_2dz, dRcv_2dz
2291 real :: dPE_det, dPE_merge
2300 real :: h_det_to_h2, h_ml_to_h2
2301 real :: h_det_to_h1, h_ml_to_h1
2302 real :: h1_to_h2, h1_to_k0
2303 real :: h2_to_k1, h2_to_k1_rem
2308 real :: R0_det, T_det, S_det
2309 real :: Rcv_stays, R0_stays
2310 real :: T_stays, S_stays
2311 real :: dSpice_det, dSpice_stays
2315 real :: dSpice_lim, dSpice_lim2
2326 real :: dPE_time_ratio
2327 real :: dT_dS_gauge, dS_dT_gauge
2337 logical :: stable_Rcv
2346 real :: Ih, Ihdet, Ih1f, Ih2f
2347 real :: Ihk0, Ihk1, Ih12
2348 real :: dR1, dR2, dR2b, dRk1
2349 real :: dR0, dR21, dRcv
2350 real :: dRcv_stays, dRcv_det, dRcv_lim
2353 real :: h2_to_k1_lim, T_new, S_new, T_max, T_min, S_max, S_min
2354 character(len=200) :: mesg
2356 integer :: i, k, k0, k1, is, ie, nz, kb1, kb2, nkmb
2357 is = g%isc ; ie = g%iec ; nz = gv%ke
2358 kb1 = cs%nkml+1; kb2 = cs%nkml+2
2359 nkmb = cs%nkml+cs%nkbl
2360 h_neglect = gv%H_subroundoff
2361 g_2 = 0.5 * gv%g_Earth
2362 rho0xg = gv%Rho0 * gv%g_Earth
2363 idt_h2 = gv%H_to_Z**2 / dt_diag
2364 i2rho0 = 0.5 / (gv%Rho0)
2365 angstrom = gv%Angstrom_H
2368 dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 / dt_ds_gauge
2371 if (cs%nkbl /= 2)
call mom_error(fatal,
"MOM_mixed_layer"// &
2372 "CS%nkbl must be 2 in mixedlayer_detrain_2.")
2374 if (dt < cs%BL_detrain_time)
then ; dpe_time_ratio = cs%BL_detrain_time / (dt)
2375 else ; dpe_time_ratio = 1.0 ;
endif
2384 h_to_bl = 0.0 ; r0_to_bl = 0.0
2385 rcv_to_bl = 0.0 ; t_to_bl = 0.0 ; s_to_bl = 0.0
2387 do k=1,cs%nkml ;
if (h(i,k) > 0.0)
then
2388 h_to_bl = h_to_bl + h(i,k)
2389 r0_to_bl = r0_to_bl + r0(i,k)*h(i,k)
2391 rcv_to_bl = rcv_to_bl + rcv(i,k)*h(i,k)
2392 t_to_bl = t_to_bl + t(i,k)*h(i,k)
2393 s_to_bl = s_to_bl + s(i,k)*h(i,k)
2395 d_ea(i,k) = d_ea(i,k) - h(i,k)
2398 if (h_to_bl > 0.0)
then ; r0_det = r0_to_bl / h_to_bl
2399 else ; r0_det = r0(i,0) ;
endif
2421 h_min_bl = min(cs%Hbuffer_min, cs%Hbuffer_rel_min*h(i,0))
2424 if (((r0(i,kb2)-r0(i,kb1)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) &
2425 stable_rcv = .false.
2427 h1 = h(i,kb1) ; h2 = h(i,kb2)
2429 h2_to_k1_rem = (h1 + h2) + h_to_bl
2430 if ((max_bl_det(i) >= 0.0) .and. (h2_to_k1_rem > max_bl_det(i))) &
2431 h2_to_k1_rem = max_bl_det(i)
2434 if ((h2 == 0.0) .and. (h1 > 0.0))
then
2443 do k1=kb2+1,nz ;
if (h(i,k1) > 2.0*angstrom)
exit ;
enddo
2445 r0(i,kb2) = r0(i,kb1)
2447 rcv(i,kb2)=rcv(i,kb1) ; t(i,kb2)=t(i,kb1) ; s(i,kb2)=s(i,kb1)
2450 if (k1 <= nz)
then ;
if (r0(i,k1) >= r0(i,kb1))
then
2451 r0(i,kb2) = 0.5*(r0(i,kb1)+r0(i,k1))
2453 rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
2454 t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
2455 s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
2459 dpe_extrap = 0.0 ; dpe_merge = 0.0
2460 mergeable_bl = .false.
2461 if ((h1 > 0.0) .and. (h2 > 0.0) .and. (h_to_bl > 0.0) .and. &
2468 do k1=kb2+1,nz ;
if (rcvtgt(k1) >= rcv(i,kb2))
exit ;
enddo ; k0 = k1-1
2469 dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2475 h1_avail = h1 - max(0.0,h_min_bl-h_to_bl)
2476 if ((k1<=nz) .and. (h2 > h_min_bl) .and. (h1_avail > 0.0) .and. &
2477 (r0(i,kb1) < r0(i,kb2)) .and. (h_to_bl*r0(i,kb1) > r0_to_bl))
then
2478 drk1 = (rcvtgt(k1) - rcv(i,kb2)) * (r0(i,kb2) - r0(i,kb1)) / &
2479 (rcv(i,kb2) - rcv(i,kb1))
2480 b1 = drk1 / (r0(i,kb2) - r0(i,kb1))
2486 h2_to_k1 = min(h2 - h_min_bl, h2_to_k1_rem)
2489 if (h2_to_k1*(h1_avail + b1*(h1_avail + h2)) > h2*h1_avail) &
2490 h2_to_k1 = (h2*h1_avail) / (h1_avail + b1*(h1_avail + h2))
2491 if (h2_to_k1*(drk1 * h2) > (h_to_bl*r0(i,kb1) - r0_to_bl) * h1) &
2492 h2_to_k1 = (h_to_bl*r0(i,kb1) - r0_to_bl) * h1 / (drk1 * h2)
2494 if ((k1==kb2+1) .and. (cs%BL_extrap_lim > 0.))
then
2497 drcv_lim = rcv(i,kb2)-rcv(i,0)
2498 do k=1,kb2 ; drcv_lim = max(drcv_lim, rcv(i,kb2)-rcv(i,k)) ;
enddo
2499 drcv_lim = cs%BL_extrap_lim*drcv_lim
2500 if ((rcvtgt(k1) - rcv(i,kb2)) >= drcv_lim)
then
2502 elseif ((rcvtgt(k1) - rcv(i,kb2)) > 0.5*drcv_lim)
then
2503 h2_to_k1 = h2_to_k1 * (2.0 - 2.0*((rcvtgt(k1) - rcv(i,kb2)) / drcv_lim))
2507 drcv = (rcvtgt(k1) - rcv(i,kb2))
2512 dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2513 dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2514 (h2 - h2_to_k1) / (h1 + h2)
2516 if (h(i,k1) > 10.0*angstrom)
then
2517 dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2518 dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2519 if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2521 if (k1<nz)
then ;
if (h(i,k1+1) > 10.0*angstrom)
then
2522 dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2523 dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2524 if ((dspice_det*dspice_lim2 > 0.0) .and. &
2525 (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2527 if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2529 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2530 t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2531 (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2532 s_det = s(i,kb2) + i_denom * &
2533 (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2535 r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2536 (s_det-s(i,kb2)) * dr0_ds(i)
2538 if (cs%BL_extrap_lim >= 0.)
then
2541 if (h(i,k1) > 10.0*angstrom)
then
2542 t_min = min(t(i,kb1), t(i,kb2), t(i,k1)) - cs%Allowed_T_chg
2543 t_max = max(t(i,kb1), t(i,kb2), t(i,k1)) + cs%Allowed_T_chg
2544 s_min = min(s(i,kb1), s(i,kb2), s(i,k1)) - cs%Allowed_S_chg
2545 s_max = max(s(i,kb1), s(i,kb2), s(i,k1)) + cs%Allowed_S_chg
2547 t_min = min(t(i,kb1), t(i,kb2)) - cs%Allowed_T_chg
2548 t_max = max(t(i,kb1), t(i,kb2)) + cs%Allowed_T_chg
2549 s_min = min(s(i,kb1), s(i,kb2)) - cs%Allowed_S_chg
2550 s_max = max(s(i,kb1), s(i,kb2)) + cs%Allowed_S_chg
2552 ihk1 = 1.0 / (h(i,k1) + h2_to_k1)
2553 t_new = (h(i,k1)*t(i,k1) + h2_to_k1*t_det) * ihk1
2554 s_new = (h(i,k1)*s(i,k1) + h2_to_k1*s_det) * ihk1
2556 if ((t_new < t_min) .or. (t_new > t_max) .or. &
2557 (s_new < s_min) .or. (s_new > s_max)) &
2561 h1_to_h2 = b1*h2*h2_to_k1 / (h2 - (1.0+b1)*h2_to_k1)
2563 ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2564 ih2f = 1.0 / ((h(i,kb2) - h2_to_k1) + h1_to_h2)
2566 rcv(i,kb2) = ((h(i,kb2)*rcv(i,kb2) - h2_to_k1*rcvtgt(k1)) + &
2567 h1_to_h2*rcv(i,kb1))*ih2f
2568 rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2570 t(i,kb2) = ((h(i,kb2)*t(i,kb2) - h2_to_k1*t_det) + &
2571 h1_to_h2*t(i,kb1)) * ih2f
2572 t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2574 s(i,kb2) = ((h(i,kb2)*s(i,kb2) - h2_to_k1*s_det) + &
2575 h1_to_h2*s(i,kb1)) * ih2f
2576 s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2579 r0(i,kb2) = ((h(i,kb2)*r0(i,kb2) - h2_to_k1*r0_det) + &
2580 h1_to_h2*r0(i,kb1)) * ih2f
2581 r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2583 h(i,kb1) = h(i,kb1) - h1_to_h2 ; h1 = h(i,kb1)
2584 h(i,kb2) = (h(i,kb2) - h2_to_k1) + h1_to_h2 ; h2 = h(i,kb2)
2585 h(i,k1) = h(i,k1) + h2_to_k1
2587 d_ea(i,kb1) = d_ea(i,kb1) - h1_to_h2
2588 d_ea(i,kb2) = (d_ea(i,kb2) - h2_to_k1) + h1_to_h2
2589 d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2590 h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2594 if ((k1>kb2+1) .and. (rcvtgt(k1-1) >= rcv(i,kb2)))
then
2595 do k1=k1,kb2+1,-1 ;
if (rcvtgt(k1-1) < rcv(i,kb2))
exit ;
enddo
2600 dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2602 if ((k0>kb2) .and. (dr1 > 0.0) .and. (h1 > h_min_bl) .and. &
2603 (h2*dr2 < h1*dr1) .and. (r0(i,kb2) > r0(i,kb1)))
then
2608 stays_merge = 2.0*(h1+h2)*(h1*dr1 - h2*dr2) / &
2609 ((dr1+dr2)*h1 + dr1*(h1+h2) + &
2610 sqrt((dr2*h1-dr1*h2)**2 + 4*(h1+h2)*h2*(dr1+dr2)*dr2))
2612 stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
2613 h1 - (h1+h2)*(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1)))
2614 if ((stays_merge > stays_min_merge) .and. &
2615 (stays_merge + h2_to_k1_rem >= h1 + h2))
then
2616 mergeable_bl = .true.
2617 dpe_merge = g_2*(r0(i,kb2)-r0(i,kb1))*(h1-stays_merge)*(h2-stays_merge)
2621 if ((k1<=nz).and.(.not.mergeable_bl))
then
2625 dr2b = rcvtgt(k1)-rcv(i,kb2) ; dr21 = rcv(i,kb2) - rcv(i,kb1)
2626 if (dr2b*(h1+h2) < h2*dr21)
then
2628 h2_to_k1 = min(h2 - (h1+h2) * dr2b / dr21, h2_to_k1_rem)
2630 if (h2 > h2_to_k1)
then
2631 drcv = (rcvtgt(k1) - rcv(i,kb2))
2636 dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2637 dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2638 (h2 - h2_to_k1) / (h1 + h2)
2640 if (h(i,k1) > 10.0*angstrom)
then
2641 dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2642 dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2643 if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2645 if (k1<nz) then;
if (h(i,k1+1) > 10.0*angstrom)
then
2646 dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2647 dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2648 if ((dspice_det*dspice_lim2 > 0.0) .and. &
2649 (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2651 if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2653 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2654 t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2655 (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2656 s_det = s(i,kb2) + i_denom * &
2657 (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2659 r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2660 (s_det-s(i,kb2)) * dr0_ds(i)
2665 ih2f = 1.0 / (h2 - h2_to_k1)
2667 t_min = min(t(i,kb2), t(i,kb1)) - cs%Allowed_T_chg
2668 t_max = max(t(i,kb2), t(i,kb1)) + cs%Allowed_T_chg
2669 t_new = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2670 if (t_new < t_min)
then
2671 h2_to_k1_lim = h2 * (t(i,kb2) - t_min) / (t_det - t_min)
2678 h2_to_k1 = h2_to_k1_lim
2679 ih2f = 1.0 / (h2 - h2_to_k1)
2680 elseif (t_new > t_max)
then
2681 h2_to_k1_lim = h2 * (t(i,kb2) - t_max) / (t_det - t_max)
2688 h2_to_k1 = h2_to_k1_lim
2689 ih2f = 1.0 / (h2 - h2_to_k1)
2691 s_min = max(min(s(i,kb2), s(i,kb1)) - cs%Allowed_S_chg, 0.0)
2692 s_max = max(s(i,kb2), s(i,kb1)) + cs%Allowed_S_chg
2693 s_new = (h2*s(i,kb2) - h2_to_k1*s_det)*ih2f
2694 if (s_new < s_min)
then
2695 h2_to_k1_lim = h2 * (s(i,kb2) - s_min) / (s_det - s_min)
2702 h2_to_k1 = h2_to_k1_lim
2703 ih2f = 1.0 / (h2 - h2_to_k1)
2704 elseif (s_new > s_max)
then
2705 h2_to_k1_lim = h2 * (s(i,kb2) - s_max) / (s_det - s_max)
2712 h2_to_k1 = h2_to_k1_lim
2713 ih2f = 1.0 / (h2 - h2_to_k1)
2716 ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2717 rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2718 rcv(i,kb2) = rcv(i,kb2) - h2_to_k1*drcv*ih2f
2720 t(i,kb2) = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2721 t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2723 s(i,kb2) = (h2*s(i,kb2) - h2_to_k1*s_det) * ih2f
2724 s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2727 r0(i,kb2) = (h2*r0(i,kb2) - h2_to_k1*r0_det) * ih2f
2728 r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2733 ihk1 = 1.0 / (h(i,k1) + h2)
2735 rcv(i,k1) = (h(i,k1)*rcv(i,k1) + h2*rcv(i,kb2)) * ihk1
2736 t(i,k1) = (h(i,k1)*t(i,k1) + h2*t(i,kb2)) * ihk1
2737 s(i,k1) = (h(i,k1)*s(i,k1) + h2*s(i,kb2)) * ihk1
2738 r0(i,k1) = (h(i,k1)*r0(i,k1) + h2*r0(i,kb2)) * ihk1
2741 h(i,k1) = h(i,k1) + h2_to_k1
2742 h(i,kb2) = h(i,kb2) - h2_to_k1 ; h2 = h(i,kb2)
2744 dpe_extrap = i2rho0*(r0_det-r0(i,kb2))*h2_to_k1*h2
2746 d_ea(i,kb2) = d_ea(i,kb2) - h2_to_k1
2747 d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2748 h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2755 h_det_h2 = max(h_min_bl-(h1+h2), 0.0)
2756 if (h_det_h2 > 0.0)
then
2762 h_det_to_h2 = min(h_to_bl, h_det_h2)
2763 h_ml_to_h2 = h_det_h2 - h_det_to_h2
2764 h_det_to_h1 = h_to_bl - h_det_to_h2
2765 h_ml_to_h1 = max(h_min_bl-h_det_to_h1,0.0)
2768 ihdet = 0.0 ;
if (h_to_bl > 0.0) ihdet = 1.0 / h_to_bl
2769 ih1f = 1.0 / (h_det_to_h1 + h_ml_to_h1)
2771 r0(i,kb2) = ((h2*r0(i,kb2) + h1*r0(i,kb1)) + &
2772 (h_det_to_h2*r0_to_bl*ihdet + h_ml_to_h2*r0(i,0))) * ih
2773 r0(i,kb1) = (h_det_to_h1*r0_to_bl*ihdet + h_ml_to_h1*r0(i,0)) * ih1f
2775 rcv(i,kb2) = ((h2*rcv(i,kb2) + h1*rcv(i,kb1)) + &
2776 (h_det_to_h2*rcv_to_bl*ihdet + h_ml_to_h2*rcv(i,0))) * ih
2777 rcv(i,kb1) = (h_det_to_h1*rcv_to_bl*ihdet + h_ml_to_h1*rcv(i,0)) * ih1f
2779 t(i,kb2) = ((h2*t(i,kb2) + h1*t(i,kb1)) + &
2780 (h_det_to_h2*t_to_bl*ihdet + h_ml_to_h2*t(i,0))) * ih
2781 t(i,kb1) = (h_det_to_h1*t_to_bl*ihdet + h_ml_to_h1*t(i,0)) * ih1f
2783 s(i,kb2) = ((h2*s(i,kb2) + h1*s(i,kb1)) + &
2784 (h_det_to_h2*s_to_bl*ihdet + h_ml_to_h2*s(i,0))) * ih
2785 s(i,kb1) = (h_det_to_h1*s_to_bl*ihdet + h_ml_to_h1*s(i,0)) * ih1f
2788 d_ea(i,1) = d_ea(i,1) - (h_ml_to_h1 + h_ml_to_h2)
2789 d_ea(i,kb1) = d_ea(i,kb1) + ((h_det_to_h1 + h_ml_to_h1) - h1)
2790 d_ea(i,kb2) = d_ea(i,kb2) + (h_min_bl - h2)
2792 h(i,kb1) = h_det_to_h1 + h_ml_to_h1 ; h(i,kb2) = h_min_bl
2793 h(i,0) = h(i,0) - (h_ml_to_h1 + h_ml_to_h2)
2796 if (
allocated(cs%diag_PE_detrain) .or.
allocated(cs%diag_PE_detrain2))
then
2797 r0_det = r0_to_bl*ihdet
2798 s1en = g_2 * idt_h2 * ( ((r0(i,kb2)-r0(i,kb1))*h1*h2 + &
2799 h_det_to_h2*( (r0(i,kb1)-r0_det)*h1 + (r0(i,kb2)-r0_det)*h2 ) + &
2800 h_ml_to_h2*( (r0(i,kb2)-r0(i,0))*h2 + (r0(i,kb1)-r0(i,0))*h1 + &
2801 (r0_det-r0(i,0))*h_det_to_h2 ) + &
2802 h_det_to_h1*h_ml_to_h1*(r0_det-r0(i,0))) - 2.0*gv%Rho0*dpe_extrap )
2804 if (
allocated(cs%diag_PE_detrain)) &
2805 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + s1en
2807 if (
allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
2808 cs%diag_PE_detrain2(i,j) + s1en + idt_h2*rho0xg*dpe_extrap
2811 elseif ((h_to_bl > 0.0) .or. (h1 < h_min_bl) .or. (h2 < h_min_bl))
then
2816 h_from_ml = h_min_bl + max(h_min_bl-h2,0.0) - h1 - h_to_bl
2817 if (h_from_ml > 0.0)
then
2820 dpe_extrap = dpe_extrap - i2rho0*h_from_ml*(r0_to_bl - r0(i,0)*h_to_bl)
2821 r0_to_bl = r0_to_bl + h_from_ml*r0(i,0)
2822 rcv_to_bl = rcv_to_bl + h_from_ml*rcv(i,0)
2823 t_to_bl = t_to_bl + h_from_ml*t(i,0)
2824 s_to_bl = s_to_bl + h_from_ml*s(i,0)
2826 h_to_bl = h_to_bl + h_from_ml
2827 h(i,0) = h(i,0) - h_from_ml
2828 d_ea(i,1) = d_ea(i,1) - h_from_ml
2833 if (r0(i,kb2) - r0(i,kb1) > 1.0e-9*abs(r0(i,kb1) - r0_det)) &
2834 b1 = abs(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1))
2835 stays_min = max((1.0-b1)*h1 - b1*h2, 0.0, h_min_bl - h_to_bl)
2836 stays_max = h1 - max(h_min_bl-h2,0.0)
2839 if (stays_max <= stays_min)
then
2841 mergeable_bl = .false.
2842 if (stays_max < h1) scale_slope = (h1 - stays_min) / (h1 - stays_max)
2847 i_ya = (h1 + h2) / ((h1 + h2) + h_to_bl)
2849 s1 = 0.5*(h1 + (h2 - bh0) * i_ya) ; s2 = h1 - s1
2854 s3sq = i_ya*max(bh0*h1-dpe_extrap, 0.0)
2856 s3sq = i_ya*(bh0*h1-min(dpe_extrap,0.0))
2859 if (s3sq == 0.0)
then
2862 elseif (s2*s2 <= s3sq)
then
2870 if (bh0 <= 0.0)
then ; stays = h1
2871 elseif (s2 > 0.0)
then
2873 if (s1 >= stays_max)
then ; stays = stays_max
2874 elseif (s1 >= 0.0)
then ; stays = s1 + sqrt(s2*s2 - s3sq)
2875 else ; stays = (h1*(s2-s1) - s3sq) / (-s1 + sqrt(s2*s2 - s3sq))
2879 if (s1 <= stays_min)
then ; stays = stays_min
2880 else ; stays = (h1*(s1-s2) + s3sq) / (s1 + sqrt(s2*s2 - s3sq))
2889 if (stays >= stays_max)
then ; stays = stays_max
2890 elseif (stays < stays_min)
then ; stays = stays_min
2894 dpe_det = g_2*((r0(i,kb1)*h_to_bl - r0_to_bl)*stays + &
2895 (r0(i,kb2)-r0(i,kb1)) * (h1-stays) * &
2896 (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
2899 if (dpe_time_ratio*h_to_bl > h_to_bl+h(i,0))
then
2900 dpe_ratio = (h_to_bl+h(i,0)) / h_to_bl
2902 dpe_ratio = dpe_time_ratio
2905 if ((mergeable_bl) .and. (num_events*dpe_ratio*dpe_det > dpe_merge))
then
2910 h1_to_k0 = (h1-stays_merge)
2911 stays = max(h_min_bl-h_to_bl,0.0)
2912 h1_to_h2 = stays_merge - stays
2914 ihk0 = 1.0 / ((h1_to_k0 + h2) + h(i,k0))
2915 ih1f = 1.0 / (h_to_bl + stays); ih2f = 1.0 / h1_to_h2
2916 ih12 = 1.0 / (h1 + h2)
2918 drcv_2dz = (rcv(i,kb1) - rcv(i,kb2)) * ih12
2919 drcv_stays = drcv_2dz*(h1_to_k0 + h1_to_h2)
2920 drcv_det = - drcv_2dz*(stays + h1_to_h2)
2921 rcv(i,k0) = ((h1_to_k0*(rcv(i,kb1) + drcv_det) + &
2922 h2*rcv(i,kb2)) + h(i,k0)*rcv(i,k0)) * ihk0
2923 rcv(i,kb2) = rcv(i,kb1) + drcv_2dz*(h1_to_k0-stays)
2924 rcv(i,kb1) = (rcv_to_bl + stays*(rcv(i,kb1) + drcv_stays)) * ih1f
2929 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2930 dspice_2dz = (ds_dt_gauge*drcv_ds(i)*(t(i,kb1)-t(i,kb2)) - &
2931 dt_ds_gauge*drcv_dt(i)*(s(i,kb1)-s(i,kb2))) * ih12
2932 dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
2933 dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
2934 if (dspice_lim * dspice_2dz <= 0.0) dspice_2dz = 0.0
2936 if (stays > 0.0)
then
2938 if (abs(dspice_lim) < abs(dspice_2dz*(h1_to_k0 + h1_to_h2))) &
2939 dspice_2dz = dspice_lim/(h1_to_k0 + h1_to_h2)
2941 dspice_stays = dspice_2dz*(h1_to_k0 + h1_to_h2)
2942 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
2943 (dt_ds_gauge * drcv_dt(i) * drcv_stays + drcv_ds(i) * dspice_stays)
2944 s_stays = s(i,kb1) + i_denom * &
2945 (drcv_ds(i) * drcv_stays - dt_ds_gauge * drcv_dt(i) * dspice_stays)
2947 r0_stays = r0(i,kb1) + (t_stays-t(i,kb1)) * dr0_dt(i) + &
2948 (s_stays-s(i,kb1)) * dr0_ds(i)
2951 if (abs(dspice_lim) < abs(dspice_2dz*h1_to_k0)) &
2952 dspice_2dz = dspice_lim/h1_to_k0
2954 t_stays = 0.0 ; s_stays = 0.0 ; r0_stays = 0.0
2957 dspice_det = - dspice_2dz*(stays + h1_to_h2)
2958 t_det = t(i,kb1) + dt_ds_gauge * i_denom * &
2959 (dt_ds_gauge * drcv_dt(i) * drcv_det + drcv_ds(i) * dspice_det)
2960 s_det = s(i,kb1) + i_denom * &
2961 (drcv_ds(i) * drcv_det - dt_ds_gauge * drcv_dt(i) * dspice_det)
2963 r0_det = r0(i,kb1) + (t_det-t(i,kb1)) * dr0_dt(i) + &
2964 (s_det-s(i,kb1)) * dr0_ds(i)
2966 t(i,k0) = ((h1_to_k0*t_det + h2*t(i,kb2)) + h(i,k0)*t(i,k0)) * ihk0
2967 t(i,kb2) = (h1*t(i,kb1) - stays*t_stays - h1_to_k0*t_det) * ih2f
2968 t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
2970 s(i,k0) = ((h1_to_k0*s_det + h2*s(i,kb2)) + h(i,k0)*s(i,k0)) * ihk0
2971 s(i,kb2) = (h1*s(i,kb1) - stays*s_stays - h1_to_k0*s_det) * ih2f
2972 s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
2974 r0(i,k0) = ((h1_to_k0*r0_det + h2*r0(i,kb2)) + h(i,k0)*r0(i,k0)) * ihk0
2975 r0(i,kb2) = (h1*r0(i,kb1) - stays*r0_stays - h1_to_k0*r0_det) * ih2f
2976 r0(i,kb1) = (r0_to_bl + stays*r0_stays) * ih1f
2998 d_ea(i,kb1) = (d_ea(i,kb1) + h_to_bl) + (stays - h1)
2999 d_ea(i,kb2) = d_ea(i,kb2) + (h1_to_h2 - h2)
3000 d_ea(i,k0) = d_ea(i,k0) + (h1_to_k0 + h2)
3002 h(i,kb1) = stays + h_to_bl
3004 h(i,k0) = h(i,k0) + (h1_to_k0 + h2)
3005 if (
allocated(cs%diag_PE_detrain)) &
3006 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_merge
3007 if (
allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3008 cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap)
3013 h1_to_h2 = h1 - stays
3014 ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / (h2 + h1_to_h2)
3015 ih = 1.0 / (h1 + h2)
3016 dr0_2dz = (r0(i,kb1) - r0(i,kb2)) * ih
3017 r0(i,kb2) = (h2*r0(i,kb2) + h1_to_h2*(r0(i,kb1) - &
3018 scale_slope*dr0_2dz*stays)) * ih2f
3019 r0(i,kb1) = (r0_to_bl + stays*(r0(i,kb1) + &
3020 scale_slope*dr0_2dz*h1_to_h2)) * ih1f
3025 dr0 = scale_slope*dr0_2dz*h1_to_h2
3026 dspice_stays = (ds_dt_gauge*dr0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3027 dt_ds_gauge*dr0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
3028 scale_slope*h1_to_h2 * ih
3029 if (h_to_bl > 0.0)
then
3030 dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3031 dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) /&
3034 dspice_lim = ds_dt_gauge*dr0_ds(i)*(t(i,0)-t(i,kb1)) - &
3035 dt_ds_gauge*dr0_dt(i)*(s(i,0)-s(i,kb1))
3037 if (dspice_stays*dspice_lim <= 0.0)
then
3039 elseif (abs(dspice_stays) > abs(dspice_lim))
then
3040 dspice_stays = dspice_lim
3042 i_denom = 1.0 / (dr0_ds(i)**2 + (dt_ds_gauge*dr0_dt(i))**2)
3043 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3044 (dt_ds_gauge * dr0_dt(i) * dr0 + dr0_ds(i) * dspice_stays)
3045 s_stays = s(i,kb1) + i_denom * &
3046 (dr0_ds(i) * dr0 - dt_ds_gauge * dr0_dt(i) * dspice_stays)
3048 rcv_stays = rcv(i,kb1) + (t_stays-t(i,kb1)) * drcv_dt(i) + &
3049 (s_stays-s(i,kb1)) * drcv_ds(i)
3051 t(i,kb2) = (h2*t(i,kb2) + h1*t(i,kb1) - t_stays*stays) * ih2f
3052 t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3053 s(i,kb2) = (h2*s(i,kb2) + h1*s(i,kb1) - s_stays*stays) * ih2f
3054 s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3055 rcv(i,kb2) = (h2*rcv(i,kb2) + h1*rcv(i,kb1) - rcv_stays*stays) * ih2f
3056 rcv(i,kb1) = (rcv_to_bl + stays*rcv_stays) * ih1f
3075 d_ea(i,kb1) = d_ea(i,kb1) + ((stays - h1) + h_to_bl)
3076 d_ea(i,kb2) = d_ea(i,kb2) + h1_to_h2
3078 h(i,kb1) = stays + h_to_bl
3079 h(i,kb2) = h(i,kb2) + h1_to_h2
3081 if (
allocated(cs%diag_PE_detrain)) &
3082 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_det
3083 if (
allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3084 cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap)
3095 subroutine mixedlayer_detrain_1(h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, d_eb, &
3096 j, G, GV, US, CS, dRcv_dT, dRcv_dS, max_BL_det)
3099 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: h
3101 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: T
3102 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: S
3103 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: R0
3105 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: Rcv
3107 real,
dimension(SZK_(GV)),
intent(in) :: RcvTgt
3109 real,
intent(in) :: dt
3110 real,
intent(in) :: dt_diag
3112 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_ea
3116 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_eb
3120 integer,
intent(in) :: j
3124 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
3128 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
3131 real,
dimension(SZI_(G)),
intent(in) :: max_BL_det
3139 real :: max_det_rem(SZI_(G))
3140 real :: detrain(SZI_(G))
3142 real :: dT_dR, dS_dR, dRml, dR0_dRcv, dT_dS_wt2
3144 real :: Sdown, Tdown
3146 real :: g_H2_2Rho0dt
3153 logical :: splittable_BL(SZI_(G)), orthogonal_extrap
3156 integer :: i, is, ie, k, k1, nkmb, nz
3157 is = g%isc ; ie = g%iec ; nz = gv%ke
3158 nkmb = cs%nkml+cs%nkbl
3159 if (cs%nkbl /= 1)
call mom_error(fatal,
"MOM_mixed_layer: "// &
3160 "CS%nkbl must be 1 in mixedlayer_detrain_1.")
3162 dt_time = dt / cs%BL_detrain_time
3163 g_h2_2rho0dt = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0 * dt_diag)
3164 g_h2_2dt = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * dt_diag)
3168 do i=is,ie ;
if (h(i,k) > 0.0)
then
3169 ih = 1.0 / (h(i,nkmb) + h(i,k))
3170 if (cs%TKE_diagnostics) &
3171 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
3172 g_h2_2rho0dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3173 if (
allocated(cs%diag_PE_detrain)) &
3174 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + &
3175 g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3176 if (
allocated(cs%diag_PE_detrain2)) &
3177 cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + &
3178 g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3180 r0(i,nkmb) = (r0(i,nkmb)*h(i,nkmb) + r0(i,k)*h(i,k)) * ih
3181 rcv(i,nkmb) = (rcv(i,nkmb)*h(i,nkmb) + rcv(i,k)*h(i,k)) * ih
3182 t(i,nkmb) = (t(i,nkmb)*h(i,nkmb) + t(i,k)*h(i,k)) * ih
3183 s(i,nkmb) = (s(i,nkmb)*h(i,nkmb) + s(i,k)*h(i,k)) * ih
3185 d_ea(i,k) = d_ea(i,k) - h(i,k)
3186 d_ea(i,nkmb) = d_ea(i,nkmb) + h(i,k)
3187 h(i,nkmb) = h(i,nkmb) + h(i,k)
3193 max_det_rem(i) = 10.0 * h(i,nkmb)
3194 if (max_bl_det(i) >= 0.0) max_det_rem(i) = max_bl_det(i)
3209 do i=is,ie ;
if (h(i,nkmb) > 0.0)
then
3210 if ((r0(i,0)<r0(i,nz)) .and. (r0(i,nz)<r0(i,nkmb)))
then
3211 if ((r0(i,nz)-r0(i,0))*h(i,0) > &
3212 (r0(i,nkmb)-r0(i,nz))*h(i,nkmb))
then
3213 detrain(i) = (r0(i,nkmb)-r0(i,nz))*h(i,nkmb) / (r0(i,nkmb)-r0(i,0))
3215 detrain(i) = (r0(i,nz)-r0(i,0))*h(i,0) / (r0(i,nkmb)-r0(i,0))
3218 d_eb(i,cs%nkml) = d_eb(i,cs%nkml) + detrain(i)
3219 d_ea(i,cs%nkml) = d_ea(i,cs%nkml) - detrain(i)
3220 d_eb(i,nkmb) = d_eb(i,nkmb) - detrain(i)
3221 d_ea(i,nkmb) = d_ea(i,nkmb) + detrain(i)
3223 if (
allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3224 cs%diag_PE_detrain(i,j) + g_h2_2dt * detrain(i)* &
3225 (h(i,0) + h(i,nkmb)) * (r0(i,nkmb) - r0(i,0))
3227 r0(i,0) = r0(i,0) - detrain(i)*(r0(i,0)-r0(i,nkmb)) / h(i,0)
3228 r0(i,nkmb) = r0(i,nkmb) - detrain(i)*(r0(i,nkmb)-x1) / h(i,nkmb)
3230 rcv(i,0) = rcv(i,0) - detrain(i)*(rcv(i,0)-rcv(i,nkmb)) / h(i,0)
3231 rcv(i,nkmb) = rcv(i,nkmb) - detrain(i)*(rcv(i,nkmb)-x1) / h(i,nkmb)
3233 t(i,0) = t(i,0) - detrain(i)*(t(i,0)-t(i,nkmb)) / h(i,0)
3234 t(i,nkmb) = t(i,nkmb) - detrain(i)*(t(i,nkmb)-x1) / h(i,nkmb)
3236 s(i,0) = s(i,0) - detrain(i)*(s(i,0)-s(i,nkmb)) / h(i,0)
3237 s(i,nkmb) = s(i,nkmb) - detrain(i)*(s(i,nkmb)-x1) / h(i,nkmb)
3247 if (h(i,nkmb) > 0.0)
then ; splittable_bl(i) = .true.
3248 else ; splittable_bl(i) = .false. ;
endif
3251 dt_ds_wt2 = cs%dT_dS_wt**2
3253 do k=nz-1,nkmb+1,-1 ;
do i=is,ie
3254 if (splittable_bl(i))
then
3255 if (rcvtgt(k)<=rcv(i,nkmb))
then
3260 splittable_bl(i) = .false.
3262 k1 = k+1 ; orthogonal_extrap = .false.
3267 if ((h(i,k+1) < 10.0*gv%Angstrom_H) .or. &
3268 ((rcvtgt(k+1)-rcv(i,nkmb)) >= 0.9*(rcv(i,k1) - rcv(i,0))))
then
3269 if (k>=nz-1)
then ; orthogonal_extrap = .true.
3270 elseif ((h(i,k+2) <= 10.0*gv%Angstrom_H) .and. &
3271 ((rcvtgt(k+1)-rcv(i,nkmb)) < 0.9*(rcv(i,k+2)-rcv(i,0))))
then
3273 else ; orthogonal_extrap = .true. ;
endif
3276 if ((r0(i,0) >= r0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
3280 if (orthogonal_extrap)
then
3284 i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
3285 dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
3286 ds_dr = drcv_ds(i) * i_denom
3288 dt_dr = (t(i,0) - t(i,k1)) / (rcv(i,0) - rcv(i,k1))
3289 ds_dr = (s(i,0) - s(i,k1)) / (rcv(i,0) - rcv(i,k1))
3291 drml = dt_time * (r0(i,nkmb) - r0(i,0)) * &
3292 (rcv(i,0) - rcv(i,k1)) / (r0(i,0) - r0(i,k1))
3294 if (drml < 0.0) cycle
3295 dr0_drcv = (r0(i,0) - r0(i,k1)) / (rcv(i,0) - rcv(i,k1))
3297 if ((rcv(i,nkmb) - drml < rcvtgt(k)) .and. (max_det_rem(i) > h(i,nkmb)))
then
3299 detrain(i) = h(i,nkmb)*(rcv(i,nkmb) - rcvtgt(k)) / &
3300 (rcvtgt(k+1) - rcvtgt(k))
3302 if (
allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3303 cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * &
3304 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dr0_drcv
3306 tdown = detrain(i) * (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3307 t(i,k) = (h(i,k) * t(i,k) + &
3308 (h(i,nkmb) * t(i,nkmb) - tdown)) / &
3309 (h(i,k) + (h(i,nkmb) - detrain(i)))
3310 t(i,k+1) = (h(i,k+1) * t(i,k+1) + tdown)/ &
3311 (h(i,k+1) + detrain(i))
3313 sdown = detrain(i) * (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3314 s(i,k) = (h(i,k) * s(i,k) + &
3315 (h(i,nkmb) * s(i,nkmb) - sdown)) / &
3316 (h(i,k) + (h(i,nkmb) - detrain(i)))
3317 s(i,k+1) = (h(i,k+1) * s(i,k+1) + sdown)/ &
3318 (h(i,k+1) + detrain(i))
3320 rcv(i,nkmb) = rcv(i,0)
3322 d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3323 d_ea(i,k) = d_ea(i,k) + (h(i,nkmb) - detrain(i))
3324 d_ea(i,nkmb) = d_ea(i,nkmb) - h(i,nkmb)
3326 h(i,k+1) = h(i,k+1) + detrain(i)
3327 h(i,k) = h(i,k) + h(i,nkmb) - detrain(i)
3331 detrain(i) = h(i,nkmb) * drml / (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3332 if (detrain(i) > max_det_rem(i)) detrain(i) = max_det_rem(i)
3333 ih = 1.0 / (h(i,k+1) + detrain(i))
3335 tdown = (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3336 t(i,nkmb) = t(i,nkmb) - dt_dr * drml
3337 t(i,k+1) = (h(i,k+1) * t(i,k+1) + detrain(i) * tdown) * ih
3338 sdown = (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3342 s(i,nkmb) = s(i,nkmb) - ds_dr * drml
3343 s(i,k+1) = (h(i,k+1) * s(i,k+1) + detrain(i) * sdown) * ih
3345 d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3346 d_ea(i,nkmb) = d_ea(i,nkmb) - detrain(i)
3348 h(i,k+1) = h(i,k+1) + detrain(i)
3349 h(i,nkmb) = h(i,nkmb) - detrain(i)
3351 if (
allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3352 cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * dr0_drcv * &
3353 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3364 if (h(i,nkmb) < 0.1*h(i,0))
then
3365 h_ent = 0.1*h(i,0) - h(i,nkmb)
3367 t(i,nkmb) = (h(i,nkmb)*t(i,nkmb) + h_ent*t(i,0)) * ih
3368 s(i,nkmb) = (h(i,nkmb)*s(i,nkmb) + h_ent*s(i,0)) * ih
3370 d_ea(i,1) = d_ea(i,1) - h_ent
3371 d_ea(i,nkmb) = d_ea(i,nkmb) + h_ent
3373 h(i,0) = h(i,0) - h_ent
3374 h(i,nkmb) = h(i,nkmb) + h_ent
3382 type(time_type),
target,
intent(in) :: time
3388 type(
diag_ctrl),
target,
intent(inout) :: diag
3401 #include "version_variable.h"
3402 character(len=40) :: mdl =
"MOM_mixed_layer"
3403 real :: bl_detrain_time_dflt
3404 real :: omega_frac_dflt, ustar_min_dflt, hmix_min_m
3405 integer :: isd, ied, jsd, jed
3406 logical :: use_temperature, use_omega
3407 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3409 if (
associated(cs))
then
3410 call mom_error(warning,
"mixedlayer_init called with an associated"// &
3411 "associated control structure.")
3413 else ;
allocate(cs) ;
endif
3418 if (gv%nkml < 1)
return
3424 call log_param(param_file, mdl,
"NKML", cs%nkml, &
3425 "The number of sublayers within the mixed layer if "//&
3426 "BULKMIXEDLAYER is true.", units=
"nondim", default=2)
3427 cs%nkbl = gv%nk_rho_varies - gv%nkml
3428 call log_param(param_file, mdl,
"NKBL", cs%nkbl, &
3429 "The number of variable density buffer layers if "//&
3430 "BULKMIXEDLAYER is true.", units=
"nondim", default=2)
3431 call get_param(param_file, mdl,
"MSTAR", cs%mstar, &
3432 "The ratio of the friction velocity cubed to the TKE "//&
3433 "input to the mixed layer.",
"units=nondim", default=1.2)
3434 call get_param(param_file, mdl,
"NSTAR", cs%nstar, &
3435 "The portion of the buoyant potential energy imparted by "//&
3436 "surface fluxes that is available to drive entrainment "//&
3437 "at the base of mixed layer when that energy is positive.", &
3438 units=
"nondim", default=0.15)
3439 call get_param(param_file, mdl,
"BULK_RI_ML", cs%bulk_Ri_ML, &
3440 "The efficiency with which mean kinetic energy released "//&
3441 "by mechanically forced entrainment of the mixed layer "//&
3442 "is converted to turbulent kinetic energy.", units=
"nondim",&
3443 fail_if_missing=.true.)
3444 call get_param(param_file, mdl,
"ABSORB_ALL_SW", cs%absorb_all_sw, &
3445 "If true, all shortwave radiation is absorbed by the "//&
3446 "ocean, instead of passing through to the bottom mud.", &
3448 call get_param(param_file, mdl,
"TKE_DECAY", cs%TKE_decay, &
3449 "TKE_DECAY relates the vertical rate of decay of the "//&
3450 "TKE available for mechanical entrainment to the natural "//&
3451 "Ekman depth.", units=
"nondim", default=2.5)
3452 call get_param(param_file, mdl,
"NSTAR2", cs%nstar2, &
3453 "The portion of any potential energy released by "//&
3454 "convective adjustment that is available to drive "//&
3455 "entrainment at the base of mixed layer. By default "//&
3456 "NSTAR2=NSTAR.", units=
"nondim", default=cs%nstar)
3457 call get_param(param_file, mdl,
"BULK_RI_CONVECTIVE", cs%bulk_Ri_convective, &
3458 "The efficiency with which convectively released mean "//&
3459 "kinetic energy is converted to turbulent kinetic "//&
3460 "energy. By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
3461 units=
"nondim", default=cs%bulk_Ri_ML)
3462 call get_param(param_file, mdl,
"HMIX_MIN", cs%Hmix_min, &
3463 "The minimum mixed layer depth if the mixed layer depth "//&
3464 "is determined dynamically.", units=
"m", default=0.0, scale=gv%m_to_H, &
3465 unscaled=hmix_min_m)
3467 call get_param(param_file, mdl,
"LIMIT_BUFFER_DETRAIN", cs%limit_det, &
3468 "If true, limit the detrainment from the buffer layers "//&
3469 "to not be too different from the neighbors.", default=.false.)
3470 call get_param(param_file, mdl,
"ALLOWED_DETRAIN_TEMP_CHG", cs%Allowed_T_chg, &
3471 "The amount by which temperature is allowed to exceed "//&
3472 "previous values during detrainment.", units=
"K", default=0.5)
3473 call get_param(param_file, mdl,
"ALLOWED_DETRAIN_SALT_CHG", cs%Allowed_S_chg, &
3474 "The amount by which salinity is allowed to exceed "//&
3475 "previous values during detrainment.", units=
"PSU", default=0.1)
3476 call get_param(param_file, mdl,
"ML_DT_DS_WEIGHT", cs%dT_dS_wt, &
3477 "When forced to extrapolate T & S to match the layer "//&
3478 "densities, this factor (in deg C / PSU) is combined "//&
3479 "with the derivatives of density with T & S to determine "//&
3480 "what direction is orthogonal to density contours. It "//&
3481 "should be a typical value of (dR/dS) / (dR/dT) in "//&
3482 "oceanic profiles.", units=
"degC PSU-1", default=6.0)
3483 call get_param(param_file, mdl,
"BUFFER_LAYER_EXTRAP_LIMIT", cs%BL_extrap_lim, &
3484 "A limit on the density range over which extrapolation "//&
3485 "can occur when detraining from the buffer layers, "//&
3486 "relative to the density range within the mixed and "//&
3487 "buffer layers, when the detrainment is going into the "//&
3488 "lightest interior layer, nondimensional, or a negative "//&
3489 "value not to apply this limit.", units=
"nondim", default = -1.0)
3490 call get_param(param_file, mdl,
"BUFFER_LAYER_HMIN_THICK", cs%Hbuffer_min, &
3491 "The minimum buffer layer thickness when the mixed layer is very thick.", &
3492 units=
"m", default=5.0, scale=gv%m_to_H)
3493 call get_param(param_file, mdl,
"BUFFER_LAYER_HMIN_REL", cs%Hbuffer_rel_min, &
3494 "The minimum buffer layer thickness relative to the combined mixed "//&
3495 "land buffer ayer thicknesses when they are thin.", &
3496 units=
"nondim", default=0.1/cs%nkbl)
3497 bl_detrain_time_dflt = 4.0*3600.0 ;
if (cs%nkbl==1) bl_detrain_time_dflt = 86400.0*30.0
3498 call get_param(param_file, mdl,
"BUFFER_LAY_DETRAIN_TIME", cs%BL_detrain_time, &
3499 "A timescale that characterizes buffer layer detrainment events.", &
3500 units=
"s", default=bl_detrain_time_dflt, scale=us%s_to_T)
3501 call get_param(param_file, mdl,
"BUFFER_SPLIT_RHO_TOL", cs%BL_split_rho_tol, &
3502 "The fractional tolerance for matching layer target densities when splitting "//&
3503 "layers to deal with massive interior layers that are lighter than one of the "//&
3504 "mixed or buffer layers.", units=
"nondim", default=0.1)
3506 call get_param(param_file, mdl,
"DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
3507 "The surface fluxes are scaled away when the total ocean "//&
3508 "depth is less than DEPTH_LIMIT_FLUXES.", &
3509 units=
"m", default=0.1*hmix_min_m, scale=gv%m_to_H)
3510 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
3511 "The rotation rate of the earth.", &
3512 default=7.2921e-5, units=
"s-1", scale=us%T_to_s)
3513 call get_param(param_file, mdl,
"ML_USE_OMEGA", use_omega, &
3514 "If true, use the absolute rotation rate instead of the "//&
3515 "vertical component of rotation when setting the decay "//&
3516 "scale for turbulence.", default=.false., do_not_log=.true.)
3517 omega_frac_dflt = 0.0
3519 call mom_error(warning,
"ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
3520 omega_frac_dflt = 1.0
3522 call get_param(param_file, mdl,
"ML_OMEGA_FRAC", cs%omega_frac, &
3523 "When setting the decay scale for turbulence, use this "//&
3524 "fraction of the absolute rotation rate blended with the "//&
3525 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
3526 units=
"nondim", default=omega_frac_dflt)
3527 call get_param(param_file, mdl,
"ML_RESORT", cs%ML_resort, &
3528 "If true, resort the topmost layers by potential density "//&
3529 "before the mixed layer calculations.", default=.false.)
3531 call get_param(param_file, mdl,
"ML_PRESORT_NK_CONV_ADJ", cs%ML_presort_nz_conv_adj, &
3532 "Convectively mix the first ML_PRESORT_NK_CONV_ADJ "//&
3533 "layers before sorting when ML_RESORT is true.", &
3534 units=
"nondim", default=0, fail_if_missing=.true.)
3536 ustar_min_dflt = 2e-4*us%s_to_T*cs%omega*(gv%Angstrom_m + gv%H_to_m*gv%H_subroundoff)
3537 call get_param(param_file, mdl,
"BML_USTAR_MIN", cs%ustar_min, &
3538 "The minimum value of ustar that should be used by the "//&
3539 "bulk mixed layer model in setting vertical TKE decay "//&
3540 "scales. This must be greater than 0.", units=
"m s-1", &
3541 default=ustar_min_dflt, scale=us%m_to_Z*us%T_to_s)
3542 if (cs%ustar_min<=0.0)
call mom_error(fatal,
"BML_USTAR_MIN must be positive.")
3544 call get_param(param_file, mdl,
"RESOLVE_EKMAN", cs%Resolve_Ekman, &
3545 "If true, the NKML>1 layers in the mixed layer are "//&
3546 "chosen to optimally represent the impact of the Ekman "//&
3547 "transport on the mixed layer TKE budget. Otherwise, "//&
3548 "the sublayers are distributed uniformly through the "//&
3549 "mixed layer.", default=.false.)
3550 call get_param(param_file, mdl,
"CORRECT_ABSORPTION_DEPTH", cs%correct_absorption, &
3551 "If true, the average depth at which penetrating shortwave "//&
3552 "radiation is absorbed is adjusted to match the average "//&
3553 "heating depth of an exponential profile by moving some "//&
3554 "of the heating upward in the water column.", default=.false.)
3555 call get_param(param_file, mdl,
"DO_RIVERMIX", cs%do_rivermix, &
3556 "If true, apply additional mixing wherever there is "//&
3557 "runoff, so that it is mixed down to RIVERMIX_DEPTH, "//&
3558 "if the ocean is that deep.", default=.false.)
3559 if (cs%do_rivermix) &
3560 call get_param(param_file, mdl,
"RIVERMIX_DEPTH", cs%rivermix_depth, &
3561 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
3562 "defined.", units=
"m", default=0.0, scale=us%m_to_Z)
3563 call get_param(param_file, mdl,
"USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
3564 "If true, use the fluxes%runoff_Hflx field to set the "//&
3565 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
3567 call get_param(param_file, mdl,
"USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
3568 "If true, use the fluxes%calving_Hflx field to set the "//&
3569 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
3572 call get_param(param_file, mdl,
"ALLOW_CLOCKS_IN_OMP_LOOPS", &
3573 cs%allow_clocks_in_omp_loops, &
3574 "If true, clocks can be called from inside loops that can "//&
3575 "be threaded. To run with multiple threads, set to False.", &
3578 cs%id_ML_depth = register_diag_field(
'ocean_model',
'h_ML', diag%axesT1, &
3579 time,
'Surface mixed layer depth',
'm')
3580 cs%id_TKE_wind = register_diag_field(
'ocean_model',
'TKE_wind', diag%axesT1, &
3581 time,
'Wind-stirring source of mixed layer TKE', &
3582 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3583 cs%id_TKE_RiBulk = register_diag_field(
'ocean_model',
'TKE_RiBulk', diag%axesT1, &
3584 time,
'Mean kinetic energy source of mixed layer TKE', &
3585 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3586 cs%id_TKE_conv = register_diag_field(
'ocean_model',
'TKE_conv', diag%axesT1, &
3587 time,
'Convective source of mixed layer TKE',
'm3 s-3', &
3588 conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3589 cs%id_TKE_pen_SW = register_diag_field(
'ocean_model',
'TKE_pen_SW', diag%axesT1, &
3590 time,
'TKE consumed by mixing penetrative shortwave radation through the mixed layer', &
3591 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3592 cs%id_TKE_mixing = register_diag_field(
'ocean_model',
'TKE_mixing', diag%axesT1, &
3593 time,
'TKE consumed by mixing that deepens the mixed layer', &
3594 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3595 cs%id_TKE_mech_decay = register_diag_field(
'ocean_model',
'TKE_mech_decay', diag%axesT1, &
3596 time,
'Mechanical energy decay sink of mixed layer TKE', &
3597 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3598 cs%id_TKE_conv_decay = register_diag_field(
'ocean_model',
'TKE_conv_decay', diag%axesT1, &
3599 time,
'Convective energy decay sink of mixed layer TKE', &
3600 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3601 cs%id_TKE_conv_s2 = register_diag_field(
'ocean_model',
'TKE_conv_s2', diag%axesT1, &
3602 time,
'Spurious source of mixed layer TKE from sigma2', &
3603 'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3604 cs%id_PE_detrain = register_diag_field(
'ocean_model',
'PE_detrain', diag%axesT1, &
3605 time,
'Spurious source of potential energy from mixed layer detrainment', &
3606 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3607 cs%id_PE_detrain2 = register_diag_field(
'ocean_model',
'PE_detrain2', diag%axesT1, &
3608 time,
'Spurious source of potential energy from mixed layer only detrainment', &
3609 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3610 cs%id_h_mismatch = register_diag_field(
'ocean_model',
'h_miss_ML', diag%axesT1, &
3611 time,
'Summed absolute mismatch in entrainment terms',
'm', conversion=us%Z_to_m)
3612 cs%id_Hsfc_used = register_diag_field(
'ocean_model',
'Hs_used', diag%axesT1, &
3613 time,
'Surface region thickness that is used',
'm', conversion=us%Z_to_m)
3614 cs%id_Hsfc_max = register_diag_field(
'ocean_model',
'Hs_max', diag%axesT1, &
3615 time,
'Maximum surface region thickness',
'm', conversion=us%Z_to_m)
3616 cs%id_Hsfc_min = register_diag_field(
'ocean_model',
'Hs_min', diag%axesT1, &
3617 time,
'Minimum surface region thickness',
'm', conversion=us%Z_to_m)
3619 if (cs%limit_det .or. (cs%id_Hsfc_min > 0))
then
3620 call get_param(param_file, mdl,
"LIMIT_BUFFER_DET_DH_SFC", cs%lim_det_dH_sfc, &
3621 "The fractional limit in the change between grid points "//&
3622 "of the surface region (mixed & buffer layer) thickness.", &
3623 units=
"nondim", default=0.5)
3624 call get_param(param_file, mdl,
"LIMIT_BUFFER_DET_DH_BATHY", cs%lim_det_dH_bathy, &
3625 "The fraction of the total depth by which the thickness "//&
3626 "of the surface region (mixed & buffer layer) is allowed "//&
3627 "to change between grid points.", units=
"nondim", default=0.2)
3630 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", use_temperature, &
3631 "If true, temperature and salinity are used as state "//&
3632 "variables.", default=.true.)
3634 if (use_temperature)
then
3635 call get_param(param_file, mdl,
"PEN_SW_NBANDS", cs%nsw, default=1)
3639 if (max(cs%id_TKE_wind, cs%id_TKE_RiBulk, cs%id_TKE_conv, cs%id_TKE_mixing, &
3640 cs%id_TKE_pen_SW, cs%id_TKE_mech_decay, cs%id_TKE_conv_decay) > 0)
then
3641 call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
3642 call safe_alloc_alloc(cs%diag_TKE_RiBulk, isd, ied, jsd, jed)
3643 call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
3644 call safe_alloc_alloc(cs%diag_TKE_pen_SW, isd, ied, jsd, jed)
3645 call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
3646 call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
3647 call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
3648 call safe_alloc_alloc(cs%diag_TKE_conv_s2, isd, ied, jsd, jed)
3650 cs%TKE_diagnostics = .true.
3652 if (cs%id_PE_detrain > 0)
call safe_alloc_alloc(cs%diag_PE_detrain, isd, ied, jsd, jed)
3653 if (cs%id_PE_detrain2 > 0)
call safe_alloc_alloc(cs%diag_PE_detrain2, isd, ied, jsd, jed)
3654 if (cs%id_ML_depth > 0)
call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
3656 if (cs%allow_clocks_in_omp_loops)
then
3657 id_clock_detrain = cpu_clock_id(
'(Ocean mixed layer detrain)', grain=clock_routine)
3658 id_clock_mech = cpu_clock_id(
'(Ocean mixed layer mechanical entrainment)', grain=clock_routine)
3659 id_clock_conv = cpu_clock_id(
'(Ocean mixed layer convection)', grain=clock_routine)
3660 if (cs%ML_resort)
then
3661 id_clock_resort = cpu_clock_id(
'(Ocean mixed layer resorting)', grain=clock_routine)
3663 id_clock_adjustment = cpu_clock_id(
'(Ocean mixed layer convective adjustment)', grain=clock_routine)
3665 id_clock_eos = cpu_clock_id(
'(Ocean mixed layer EOS)', grain=clock_routine)
3668 if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) &
3669 id_clock_pass = cpu_clock_id(
'(Ocean mixed layer halo updates)', grain=clock_routine)
3682 function ef4(Ht, En, I_L, dR_de)
3683 real,
intent(in) :: ht
3684 real,
intent(in) :: en
3685 real,
intent(in) :: i_l
3686 real,
optional,
intent(inout) :: dr_de
3694 exp_lhpe = exp(-i_l*(en+ht))
3696 res = exp_lhpe * (en*i_hpe/ht - 0.5*i_l*log(ht*i_hpe) + 0.5*i_l*i_l*en)
3697 if (
PRESENT(dr_de)) &
3698 dr_de = -i_l*res + exp_lhpe*(i_hpe*i_hpe + 0.5*i_l*i_hpe + 0.5*i_l*i_l)