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>
24 public bulkmixedlayer, bulkmixedlayer_init
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
151 integer :: id_clock_detrain=0, id_clock_mech=0, id_clock_conv=0, id_clock_adjustment=0
152 integer :: id_clock_eos=0, id_clock_resort=0, id_clock_pass=0
187 subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt_in_T, 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_in_t
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)) :: &
352 logical :: write_diags
353 logical :: reset_diags
354 integer :: i, j, k, is, ie, js, je, nz, nkmb, n
357 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
359 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_mixed_layer: "//&
360 "Module must be initialized before it is used.")
361 if (gv%nkml < 1)
return
363 if (.not.
associated(tv%eqn_of_state))
call mom_error(fatal, &
364 "MOM_mixed_layer: Temperature, salinity and an equation of state "//&
366 if (.NOT.
associated(fluxes%ustar))
call mom_error(fatal, &
367 "MOM_mixed_layer: No surface TKE fluxes (ustar) defined in mixedlayer!")
369 nkmb = cs%nkml+cs%nkbl
370 inkml = 1.0 / real(cs%nkml)
371 if (cs%nkml > 1) inkmlm1 = 1.0 / real(cs%nkml-1)
375 irho0 = 1.0 / gv%Rho0
376 dt__diag = dt_in_t ;
if (
present(dt_diag)) dt__diag = dt_diag
377 idt_diag = 1.0 / (dt__diag)
378 write_diags = .true. ;
if (
present(last_call)) write_diags = last_call
380 p_ref(:) = 0.0 ; p_ref_cv(:) = tv%P_Ref
384 if (cs%limit_det .or. (cs%id_Hsfc_min > 0))
then
386 do j=js-1,je+1 ;
do i=is-1,ie+1
387 h_sum(i,j) = 0.0 ; hmbl_prev(i,j) = 0.0
391 do k=1,nkmb ;
do i=is-1,ie+1
392 h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
393 hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k)
395 do k=nkmb+1,nz ;
do i=is-1,ie+1
396 h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
400 call cpu_clock_begin(id_clock_pass)
403 call do_group_pass(cs%pass_h_sum_hmbl_prev, g%Domain)
404 call cpu_clock_end(id_clock_pass)
409 if (
present(dt_diag) .and. write_diags .and. (dt__diag > dt_in_t)) &
410 reset_diags = .false.
412 if (reset_diags)
then
413 if (cs%TKE_diagnostics)
then
415 do j=js,je ;
do i=is,ie
416 cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_RiBulk(i,j) = 0.0
417 cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_pen_SW(i,j) = 0.0
418 cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
419 cs%diag_TKE_conv_decay(i,j) = 0.0 ; cs%diag_TKE_conv_s2(i,j) = 0.0
422 if (
allocated(cs%diag_PE_detrain))
then
424 do j=js,je ;
do i=is,ie
425 cs%diag_PE_detrain(i,j) = 0.0
428 if (
allocated(cs%diag_PE_detrain2))
then
430 do j=js,je ;
do i=is,ie
431 cs%diag_PE_detrain2(i,j) = 0.0
436 if (cs%ML_resort)
then
437 do i=is,ie ; h_ca(i) = 0.0 ;
enddo
438 do k=1,nz ;
do i=is,ie ; dke_ca(i,k) = 0.0 ; ctke(i,k) = 0.0 ;
enddo ;
enddo
452 do k=1,nz ;
do i=is,ie
453 h(i,k) = h_3d(i,j,k) ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
454 h_orig(i,k) = h_3d(i,j,k)
455 eps(i,k) = 0.0 ;
if (k > nkmb) eps(i,k) = gv%Angstrom_H
456 t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
458 if (nsw>0)
call extract_optics_slice(optics, j, g, gv, opacity=opacity_band, opacity_scale=gv%H_to_m)
460 do k=1,nz ;
do i=is,ie
461 d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0
464 if (id_clock_eos>0)
call cpu_clock_begin(id_clock_eos)
466 do i=is,ie ; p_ref(i) = 0.0 ;
enddo
467 do k=1,cs%nkml ;
do i=is,ie
468 p_ref(i) = p_ref(i) + 0.5*gv%H_to_Pa*h(i,k)
471 is, ie-is+1, tv%eqn_of_state)
473 is, ie-is+1, tv%eqn_of_state)
478 ie-is+1, tv%eqn_of_state)
480 if (id_clock_eos>0)
call cpu_clock_end(id_clock_eos)
482 if (cs%ML_resort)
then
483 if (id_clock_resort>0)
call cpu_clock_begin(id_clock_resort)
484 if (cs%ML_presort_nz_conv_adj > 0) &
485 call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
486 s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, us, cs, &
487 cs%ML_presort_nz_conv_adj)
489 call sort_ml(h(:,1:), r0(:,1:), eps, g, gv, cs, ksort)
490 if (id_clock_resort>0)
call cpu_clock_end(id_clock_resort)
492 do k=1,nz ;
do i=is,ie ; ksort(i,k) = k ;
enddo ;
enddo
494 if (id_clock_adjustment>0)
call cpu_clock_begin(id_clock_adjustment)
498 call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
499 s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, us, cs)
500 do i=is,ie ; h_ca(i) = h(i,1) ;
enddo
502 if (id_clock_adjustment>0)
call cpu_clock_end(id_clock_adjustment)
505 if (
associated(fluxes%lrunoff) .and. cs%do_rivermix)
then
517 rmixconst = 0.5*cs%rivermix_depth * (us%L_to_m**2*gv%g_Earth*us%m_to_Z) * irho0**2
519 tke_river(i) = max(0.0, rmixconst*dr0_ds(i)* &
520 us%T_to_s*(fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * s(i,1))
523 do i=is,ie ; tke_river(i) = 0.0 ;
enddo
527 if (id_clock_conv>0)
call cpu_clock_begin(id_clock_conv)
536 call extractfluxes1d(g, gv, fluxes, optics, nsw, j, us%T_to_s*dt_in_t, &
537 cs%H_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
538 h(:,1:), t(:,1:), netmassinout, netmassout, net_heat, net_salt, pen_sw_bnd,&
539 tv, aggregate_fw_forcing)
542 call mixedlayer_convection(h(:,1:), d_eb, htot, ttot, stot, uhtot, vhtot, &
543 r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), &
544 r0(:,1:), rcv(:,1:), eps, &
545 dr0_dt, drcv_dt, dr0_ds, drcv_ds, &
546 netmassinout, netmassout, net_heat, net_salt, &
547 nsw, pen_sw_bnd, opacity_band, conv_en, &
548 dke_fc, j, ksort, g, gv, us, cs, tv, fluxes, dt_in_t, &
549 aggregate_fw_forcing)
551 if (id_clock_conv>0)
call cpu_clock_end(id_clock_conv)
559 if (id_clock_mech>0)
call cpu_clock_begin(id_clock_mech)
561 call find_starting_tke(htot, h_ca, fluxes, conv_en, ctke, dke_fc, dke_ca, &
562 tke, tke_river, idecay_len_tke, cmke, dt_in_t, idt_diag, &
563 j, ksort, g, gv, us, cs)
566 call mechanical_entrainment(h(:,1:), d_eb, htot, ttot, stot, uhtot, vhtot, &
567 r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), r0(:,1:), rcv(:,1:), eps, dr0_dt, drcv_dt, &
568 cmke, idt_diag, nsw, pen_sw_bnd, opacity_band, tke, &
569 idecay_len_tke, j, ksort, g, gv, us, cs)
571 call absorbremainingsw(g, gv, us, h(:,1:), opacity_band, nsw, optics, j, dt_in_t, &
572 cs%H_limit_fluxes, cs%correct_absorption, cs%absorb_all_SW, &
573 t(:,1:), pen_sw_bnd, eps, ksort, htot, ttot)
575 if (cs%TKE_diagnostics)
then ;
do i=is,ie
576 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) - idt_diag*tke(i)
578 if (id_clock_mech>0)
call cpu_clock_end(id_clock_mech)
581 do i=is,ie ;
if (htot(i) > 0.0)
then
583 r0(i,0) = r0_tot(i) * ih ; rcv(i,0) = rcv_tot(i) * ih
584 t(i,0) = ttot(i) * ih ; s(i,0) = stot(i) * ih
587 t(i,0) = t(i,1) ; s(i,0) = s(i,1) ; r0(i,0) = r0(i,1) ; rcv(i,0) = rcv(i,1)
590 if (write_diags .and.
allocated(cs%ML_depth))
then ;
do i=is,ie
591 cs%ML_depth(i,j) = h(i,0) * gv%H_to_m
593 if (
associated(hml))
then ;
do i=is,ie
594 hml(i,j) = g%mask2dT(i,j) * (h(i,0) * gv%H_to_m)
607 if (cs%ML_resort)
then
608 if (id_clock_resort>0)
call cpu_clock_begin(id_clock_resort)
609 call resort_ml(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), gv%Rlay, eps, &
610 d_ea, d_eb, ksort, g, gv, cs, dr0_dt, dr0_ds, drcv_dt, drcv_ds)
611 if (id_clock_resort>0)
call cpu_clock_end(id_clock_resort)
614 if (cs%limit_det .or. (cs%id_Hsfc_max > 0) .or. (cs%id_Hsfc_min > 0))
then
615 do i=is,ie ; hsfc(i) = h(i,0) ;
enddo
616 do k=1,nkmb ;
do i=is,ie ; hsfc(i) = hsfc(i) + h(i,k) ;
enddo ;
enddo
618 if (cs%limit_det .or. (cs%id_Hsfc_min > 0))
then
619 dhsfc = cs%lim_det_dH_sfc ; dhd = cs%lim_det_dH_bathy
621 h_nbr = min(dhsfc*max(hmbl_prev(i-1,j), hmbl_prev(i+1,j), &
622 hmbl_prev(i,j-1), hmbl_prev(i,j+1)), &
623 max(hmbl_prev(i-1,j) - dhd*min(h_sum(i,j),h_sum(i-1,j)), &
624 hmbl_prev(i+1,j) - dhd*min(h_sum(i,j),h_sum(i+1,j)), &
625 hmbl_prev(i,j-1) - dhd*min(h_sum(i,j),h_sum(i,j-1)), &
626 hmbl_prev(i,j+1) - dhd*min(h_sum(i,j),h_sum(i,j+1))) )
628 hsfc_min(i,j) = gv%H_to_Z * max(h(i,0), min(hsfc(i), h_nbr))
630 if (cs%limit_det) max_bl_det(i) = max(0.0, hsfc(i)-h_nbr)
634 if (cs%id_Hsfc_max > 0)
then ;
do i=is,ie
635 hsfc_max(i,j) = gv%H_to_Z * hsfc(i)
642 if (id_clock_detrain>0)
call cpu_clock_begin(id_clock_detrain)
643 if (cs%nkbl == 1)
then
644 call mixedlayer_detrain_1(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
645 gv%Rlay, dt_in_t, dt__diag, d_ea, d_eb, j, g, gv, us, cs, &
646 drcv_dt, drcv_ds, max_bl_det)
647 elseif (cs%nkbl == 2)
then
648 call mixedlayer_detrain_2(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
649 gv%Rlay, dt_in_t, dt__diag, d_ea, j, g, gv, us, cs, &
650 dr0_dt, dr0_ds, drcv_dt, drcv_ds, max_bl_det)
653 call mom_error(fatal,
"MOM_mixed_layer: CS%nkbl must be 1 or 2 for now.")
655 if (id_clock_detrain>0)
call cpu_clock_end(id_clock_detrain)
658 if (cs%id_Hsfc_used > 0)
then
659 do i=is,ie ; hsfc_used(i,j) = gv%H_to_Z * h(i,0) ;
enddo
660 do k=cs%nkml+1,nkmb ;
do i=is,ie
661 hsfc_used(i,j) = hsfc_used(i,j) + gv%H_to_Z * h(i,k)
667 if (cs%Resolve_Ekman .and. (cs%nkml>1))
then
676 ku_star = 0.41*fluxes%ustar(i,j)
677 if (
associated(fluxes%ustar_shelf) .and. &
678 associated(fluxes%frac_shelf_h))
then
679 if (fluxes%frac_shelf_h(i,j) > 0.0) &
680 ku_star = (1.0 - fluxes%frac_shelf_h(i,j)) * ku_star + &
681 fluxes%frac_shelf_h(i,j) * (0.41*fluxes%ustar_shelf(i,j))
683 absf_x_h = 0.25 * gv%H_to_Z * h(i,0) * &
684 ((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
685 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
688 h_3d(i,j,1) = h(i,0) / (3.0 + sqrt(absf_x_h*(absf_x_h + 2.0*ku_star) / ku_star**2))
691 h_3d(i,j,k) = (h(i,0)-h_3d(i,j,1)) * inkmlm1
692 d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
693 d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
698 h_3d(i,j,1) = h(i,0) * inkml
700 do k=2,cs%nkml ;
do i=is,ie
701 h_3d(i,j,k) = h(i,0) * inkml
702 d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
703 d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
706 do i=is,ie ; h(i,0) = 0.0 ;
enddo
707 do k=1,cs%nkml ;
do i=is,ie
708 tv%T(i,j,k) = t(i,0) ; tv%S(i,j,k) = s(i,0)
724 eaml(i,nz) = (h(i,nz) - h_orig(i,nz)) - d_eb(i,nz)
726 do k=nz-1,1,-1 ;
do i=is,ie
727 ebml(i,k) = ebml(i,k+1) - d_eb(i,k+1)
728 eaml(i,k) = eaml(i,k+1) + d_ea(i,k)
730 do i=is,ie ; eaml(i,1) = netmassinout(i) ;
enddo
733 do k=cs%nkml+1,nz ;
do i=is,ie
734 h_3d(i,j,k) = h(i,k); tv%T(i,j,k) = t(i,k) ; tv%S(i,j,k) = s(i,k)
737 do k=1,nz ;
do i=is,ie
738 ea(i,j,k) = ea(i,j,k) + eaml(i,k)
739 eb(i,j,k) = eb(i,j,k) + ebml(i,k)
742 if (cs%id_h_mismatch > 0)
then
744 h_miss(i,j) = gv%H_to_Z * abs(h_3d(i,j,1) - (h_orig(i,1) + &
745 (eaml(i,1) + (ebml(i,1) - eaml(i,1+1)))))
747 do k=2,nz-1 ;
do i=is,ie
748 h_miss(i,j) = h_miss(i,j) + gv%H_to_Z * abs(h_3d(i,j,k) - (h_orig(i,k) + &
749 ((eaml(i,k) - ebml(i,k-1)) + (ebml(i,k) - eaml(i,k+1)))))
752 h_miss(i,j) = h_miss(i,j) + gv%H_to_Z * abs(h_3d(i,j,nz) - (h_orig(i,nz) + &
753 ((eaml(i,nz) - ebml(i,nz-1)) + ebml(i,nz))))
763 call diag_update_remap_grids(cs%diag)
766 if (write_diags)
then
767 if (cs%id_ML_depth > 0) &
768 call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
769 if (cs%id_TKE_wind > 0) &
770 call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
771 if (cs%id_TKE_RiBulk > 0) &
772 call post_data(cs%id_TKE_RiBulk, cs%diag_TKE_RiBulk, cs%diag)
773 if (cs%id_TKE_conv > 0) &
774 call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
775 if (cs%id_TKE_pen_SW > 0) &
776 call post_data(cs%id_TKE_pen_SW, cs%diag_TKE_pen_SW, cs%diag)
777 if (cs%id_TKE_mixing > 0) &
778 call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
779 if (cs%id_TKE_mech_decay > 0) &
780 call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
781 if (cs%id_TKE_conv_decay > 0) &
782 call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
783 if (cs%id_TKE_conv_s2 > 0) &
784 call post_data(cs%id_TKE_conv_s2, cs%diag_TKE_conv_s2, cs%diag)
785 if (cs%id_PE_detrain > 0) &
786 call post_data(cs%id_PE_detrain, cs%diag_PE_detrain, cs%diag)
787 if (cs%id_PE_detrain2 > 0) &
788 call post_data(cs%id_PE_detrain2, cs%diag_PE_detrain2, cs%diag)
789 if (cs%id_h_mismatch > 0) &
790 call post_data(cs%id_h_mismatch, h_miss, cs%diag)
791 if (cs%id_Hsfc_used > 0) &
792 call post_data(cs%id_Hsfc_used, hsfc_used, cs%diag)
793 if (cs%id_Hsfc_max > 0) &
794 call post_data(cs%id_Hsfc_max, hsfc_max, cs%diag)
795 if (cs%id_Hsfc_min > 0) &
796 call post_data(cs%id_Hsfc_min, hsfc_min, cs%diag)
799 end subroutine bulkmixedlayer
804 subroutine convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, &
805 dKE_CA, cTKE, j, G, GV, US, CS, nz_conv)
808 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: h
810 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: u
812 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: v
814 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: T
815 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: S
816 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: R0
818 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: Rcv
820 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_eb
824 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: eps
826 real,
dimension(SZI_(G),SZK_(GV)),
intent(out) :: dKE_CA
829 real,
dimension(SZI_(G),SZK_(GV)),
intent(out) :: cTKE
832 integer,
intent(in) :: j
835 integer,
optional,
intent(in) :: nz_conv
844 real,
dimension(SZI_(G)) :: &
865 integer :: is, ie, nz, i, k, k1, nzc, nkmb
867 is = g%isc ; ie = g%iec ; nz = gv%ke
868 g_h2_2rho0 = (us%L_to_m**2*gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0)
869 nzc = nz ;
if (
present(nz_conv)) nzc = nz_conv
870 nkmb = cs%nkml+cs%nkbl
875 do k1=min(nzc-1,nkmb),1,-1
877 h_orig_k1(i) = h(i,k1)
878 ke_orig(i) = 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)
879 uhtot(i) = h(i,k1)*u(i,k1) ; vhtot(i) = h(i,k1)*v(i,k1)
880 r0_tot(i) = r0(i,k1) * h(i,k1)
881 ctke(i,k1) = 0.0 ; dke_ca(i,k1) = 0.0
883 rcv_tot(i) = rcv(i,k1) * h(i,k1)
884 ttot(i) = t(i,k1) * h(i,k1) ; stot(i) = s(i,k1) * h(i,k1)
888 if ((h(i,k) > eps(i,k)) .and. (r0_tot(i) > h(i,k1)*r0(i,k)))
then
889 h_ent = h(i,k)-eps(i,k)
890 ctke(i,k1) = ctke(i,k1) + h_ent * g_h2_2rho0 * &
891 (r0_tot(i) - h(i,k1)*r0(i,k)) * cs%nstar2
893 ctke(i,k1) = ctke(i,k1) + ctke(i,k)
894 dke_ca(i,k1) = dke_ca(i,k1) + dke_ca(i,k)
896 r0_tot(i) = r0_tot(i) + h_ent * r0(i,k)
897 ke_orig(i) = ke_orig(i) + 0.5*h_ent* &
898 (u(i,k)*u(i,k) + v(i,k)*v(i,k))
899 uhtot(i) = uhtot(i) + h_ent*u(i,k)
900 vhtot(i) = vhtot(i) + h_ent*v(i,k)
902 rcv_tot(i) = rcv_tot(i) + h_ent * rcv(i,k)
903 ttot(i) = ttot(i) + h_ent * t(i,k)
904 stot(i) = stot(i) + h_ent * s(i,k)
905 h(i,k1) = h(i,k1) + h_ent ; h(i,k) = eps(i,k)
907 d_eb(i,k) = d_eb(i,k) - h_ent
908 d_eb(i,k1) = d_eb(i,k1) + h_ent
914 do i=is,ie ;
if (h(i,k1) > h_orig_k1(i))
then
916 r0(i,k1) = r0_tot(i) * ih
917 u(i,k1) = uhtot(i) * ih ; v(i,k1) = vhtot(i) * ih
918 dke_ca(i,k1) = dke_ca(i,k1) + gv%H_to_Z * us%T_to_s**2*(cs%bulk_Ri_convective * &
919 (ke_orig(i) - 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)))
920 rcv(i,k1) = rcv_tot(i) * ih
921 t(i,k1) = ttot(i) * ih ; s(i,k1) = stot(i) * ih
926 do k=2,min(nzc,nkmb) ;
do i=is,ie ;
if (h(i,k) == 0.0)
then
928 rcv(i,k) = rcv(i,k-1) ; t(i,k) = t(i,k-1) ; s(i,k) = s(i,k-1)
929 endif ;
enddo ;
enddo
931 end subroutine convective_adjustment
936 subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
937 R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, &
938 dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, &
939 netMassInOut, netMassOut, Net_heat, Net_salt, &
940 nsw, Pen_SW_bnd, opacity_band, Conv_en, &
941 dKE_FC, j, ksort, G, GV, US, CS, tv, fluxes, dt_in_T, &
942 aggregate_FW_forcing)
945 real,
dimension(SZI_(G),SZK_(GV)), &
948 real,
dimension(SZI_(G),SZK_(GV)), &
949 intent(inout) :: d_eb
952 real,
dimension(SZI_(G)),
intent(out) :: htot
953 real,
dimension(SZI_(G)),
intent(out) :: Ttot
955 real,
dimension(SZI_(G)),
intent(out) :: Stot
957 real,
dimension(SZI_(G)),
intent(out) :: uhtot
959 real,
dimension(SZI_(G)),
intent(out) :: vhtot
961 real,
dimension(SZI_(G)),
intent(out) :: R0_tot
963 real,
dimension(SZI_(G)),
intent(out) :: Rcv_tot
965 real,
dimension(SZI_(G),SZK_(GV)), &
967 real,
dimension(SZI_(G),SZK_(GV)), &
969 real,
dimension(SZI_(G),SZK_(GV)), &
971 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),SZK_(GV)), &
982 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
984 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
986 real,
dimension(SZI_(G)),
intent(in) :: dR0_dS
988 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
990 real,
dimension(SZI_(G)),
intent(in) :: netMassInOut
993 real,
dimension(SZI_(G)),
intent(in) :: netMassOut
995 real,
dimension(SZI_(G)),
intent(in) :: Net_heat
998 real,
dimension(SZI_(G)),
intent(in) :: Net_salt
1000 integer,
intent(in) :: nsw
1002 real,
dimension(max(nsw,1),SZI_(G)),
intent(inout) :: Pen_SW_bnd
1005 real,
dimension(max(nsw,1),SZI_(G),SZK_(GV)),
intent(in) :: opacity_band
1007 real,
dimension(SZI_(G)),
intent(out) :: Conv_en
1009 real,
dimension(SZI_(G)),
intent(out) :: dKE_FC
1011 integer,
intent(in) :: j
1012 integer,
dimension(SZI_(G),SZK_(GV)), &
1019 type(
forcing),
intent(inout) :: fluxes
1022 real,
intent(in) :: dt_in_T
1023 logical,
intent(in) :: aggregate_FW_forcing
1033 real,
dimension(SZI_(G)) :: &
1038 real :: Pen_absorbed
1045 real :: En_fn, Frac, x1
1047 real :: dr_ent, dr_comp
1049 real :: h_min, h_max
1064 integer :: is, ie, nz, i, k, ks, itt, n
1065 real,
dimension(max(nsw,1)) :: &
1069 angstrom = gv%Angstrom_H
1070 c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
1071 g_h2_2rho0 = (us%L_to_m**2*gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0)
1073 is = g%isc ; ie = g%iec ; nz = gv%ke
1075 do i=is,ie ;
if (ksort(i,1) > 0)
then
1078 if (aggregate_fw_forcing)
then
1080 if (netmassinout(i) < 0.0) massoutrem(i) = -netmassinout(i)
1081 netmassin(i) = netmassinout(i) + massoutrem(i)
1083 massoutrem(i) = -netmassout(i)
1084 netmassin(i) = netmassinout(i) - netmassout(i)
1088 h_ent = max(min(angstrom,h(i,k)-eps(i,k)),0.0)
1089 htot(i) = h_ent + netmassin(i)
1090 h(i,k) = h(i,k) - h_ent
1091 d_eb(i,k) = d_eb(i,k) - h_ent
1094 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1095 sw_trans = exp(-htot(i)*opacity_band(n,i,k))
1096 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0-sw_trans)
1097 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1104 ttot(i) = (net_heat(i) + (netmassin(i) * t_precip + h_ent * t(i,k))) + &
1112 stot(i) = h_ent*s(i,k) + net_salt(i)
1113 uhtot(i) = u(i,1)*netmassin(i) + u(i,k)*h_ent
1114 vhtot(i) = v(i,1)*netmassin(i) + v(i,k)*h_ent
1115 r0_tot(i) = (h_ent*r0(i,k) + netmassin(i)*r0(i,1)) + &
1117 (dr0_dt(i)*(net_heat(i) + pen_absorbed) - &
1118 dr0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1119 rcv_tot(i) = (h_ent*rcv(i,k) + netmassin(i)*rcv(i,1)) + &
1121 (drcv_dt(i)*(net_heat(i) + pen_absorbed) - &
1122 drcv_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1123 conv_en(i) = 0.0 ; dke_fc(i) = 0.0
1124 if (
associated(fluxes%heat_content_massin)) &
1125 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + us%s_to_T * &
1126 t_precip * netmassin(i) * gv%H_to_kg_m2 * fluxes%C_p * idt
1127 if (
associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1128 t_precip * netmassin(i) * gv%H_to_kg_m2
1134 do i=is,ie ;
if (ksort(i,ks) > 0)
then
1137 if ((htot(i) < angstrom) .and. (h(i,k) > eps(i,k)))
then
1140 h_ent = min(angstrom-htot(i), h(i,k)-eps(i,k))
1141 htot(i) = htot(i) + h_ent
1142 h(i,k) = h(i,k) - h_ent
1143 d_eb(i,k) = d_eb(i,k) - h_ent
1145 r0_tot(i) = r0_tot(i) + h_ent*r0(i,k)
1146 uhtot(i) = uhtot(i) + h_ent*u(i,k)
1147 vhtot(i) = vhtot(i) + h_ent*v(i,k)
1149 rcv_tot(i) = rcv_tot(i) + h_ent*rcv(i,k)
1150 ttot(i) = ttot(i) + h_ent*t(i,k)
1151 stot(i) = stot(i) + h_ent*s(i,k)
1157 if ((massoutrem(i) > 0.0) .and. (h(i,k) > eps(i,k)))
then
1158 if (massoutrem(i) > (h(i,k) - eps(i,k)))
then
1159 h_evap = h(i,k) - eps(i,k)
1161 massoutrem(i) = massoutrem(i) - h_evap
1163 h_evap = massoutrem(i)
1164 h(i,k) = h(i,k) - h_evap
1168 stot(i) = stot(i) + h_evap*s(i,k)
1169 r0_tot(i) = r0_tot(i) + dr0_ds(i)*h_evap*s(i,k)
1170 rcv_tot(i) = rcv_tot(i) + drcv_ds(i)*h_evap*s(i,k)
1171 d_eb(i,k) = d_eb(i,k) - h_evap
1177 if (
associated(fluxes%heat_content_massout)) &
1178 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) - us%s_to_T * &
1179 t(i,k)*h_evap*gv%H_to_kg_m2 * fluxes%C_p * idt
1180 if (
associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
1181 t(i,k)*h_evap*gv%H_to_kg_m2
1186 h_avail = h(i,k) - eps(i,k)
1187 if (h_avail > 0.0)
then
1188 dr = r0_tot(i) - htot(i)*r0(i,k)
1192 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1193 dr0 = dr0 - (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1194 opacity_band(n,i,k)*htot(i)
1200 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1203 opacity = opacity_band(n,i,k)
1204 sw_trans = exp(-h_avail*opacity)
1205 dr_comp = dr_comp + (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1206 ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
1208 if (dr_comp >= 0.0)
then
1218 frac = dr0 / (dr0 - dr_comp)
1219 h_ent = h_avail * frac*frac
1220 h_min = 0.0 ; h_max = h_avail
1223 r_sw_top(n) = dr0_dt(i) * pen_sw_bnd(n,i)
1224 c2(n) = r_sw_top(n) * opacity_band(n,i,k)**2
1227 dr_ent = dr ; dr_dh = 0.0
1229 opacity = opacity_band(n,i,k)
1230 sw_trans = exp(-h_ent*opacity)
1231 dr_ent = dr_ent + r_sw_top(n) * ((1.0 - sw_trans) - &
1232 opacity*(htot(i)+h_ent)*sw_trans)
1233 dr_dh = dr_dh + c2(n) * (htot(i)+h_ent) * sw_trans
1236 if (dr_ent > 0.0)
then
1242 dh_newt = -dr_ent / dr_dh
1243 h_prev = h_ent ; h_ent = h_prev+dh_newt
1244 if (h_ent > h_max)
then
1245 h_ent = 0.5*(h_prev+h_max)
1246 elseif (h_ent < h_min)
then
1247 h_ent = 0.5*(h_prev+h_min)
1250 if (abs(dh_newt) < 0.2*angstrom)
exit
1257 sum_pen_en = 0.0 ; pen_absorbed = 0.0
1258 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1259 opacity = opacity_band(n,i,k)
1260 sw_trans = exp(-h_ent*opacity)
1263 if (x1 < 2.0e-5)
then
1264 en_fn = (opacity*htot(i)*(1.0 - 0.5*(x1 - c1_3*x1)) + &
1267 en_fn = ((opacity*htot(i) + 2.0) * &
1268 ((1.0-sw_trans) / x1) - 1.0 + sw_trans)
1270 sum_pen_en = sum_pen_en - (dr0_dt(i)*pen_sw_bnd(n,i)) * en_fn
1272 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1273 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1276 conv_en(i) = conv_en(i) + g_h2_2rho0 * h_ent * &
1277 ( (r0_tot(i) - r0(i,k)*htot(i)) + sum_pen_en )
1279 r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1280 stot(i) = stot(i) + h_ent * s(i,k)
1281 ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1282 rcv_tot(i) = rcv_tot(i) + (h_ent * rcv(i,k) + pen_absorbed*drcv_dt(i))
1285 if (h_ent > 0.0)
then
1286 if (htot(i) > 0.0) &
1287 dke_fc(i) = dke_fc(i) + cs%bulk_Ri_convective * 0.5 * &
1288 ((gv%H_to_Z*h_ent) / (htot(i)*(h_ent+htot(i)))) * &
1289 us%T_to_s**2*((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1291 htot(i) = htot(i) + h_ent
1292 h(i,k) = h(i,k) - h_ent
1293 d_eb(i,k) = d_eb(i,k) - h_ent
1294 uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
1302 end subroutine mixedlayer_convection
1306 subroutine find_starting_tke(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, &
1307 TKE, TKE_river, Idecay_len_TKE, cMKE, dt_in_T, Idt_diag, &
1308 j, ksort, G, GV, US, CS)
1312 real,
dimension(SZI_(G)),
intent(in) :: htot
1314 real,
dimension(SZI_(G)),
intent(in) :: h_CA
1316 type(
forcing),
intent(in) :: fluxes
1319 real,
dimension(SZI_(G)),
intent(inout) :: Conv_En
1321 real,
dimension(SZI_(G)),
intent(in) :: dKE_FC
1324 real,
dimension(SZI_(G),SZK_(GV)), &
1328 real,
dimension(SZI_(G),SZK_(GV)), &
1329 intent(in) :: dKE_CA
1332 real,
dimension(SZI_(G)),
intent(out) :: TKE
1334 real,
dimension(SZI_(G)),
intent(out) :: Idecay_len_TKE
1336 real,
dimension(SZI_(G)),
intent(in) :: TKE_river
1339 real,
dimension(2,SZI_(G)),
intent(out) :: cMKE
1342 real,
intent(in) :: dt_in_T
1343 real,
intent(in) :: Idt_diag
1345 integer,
intent(in) :: j
1346 integer,
dimension(SZI_(G),SZK_(GV)), &
1369 real :: wind_TKE_src
1372 integer :: is, ie, nz, i
1374 is = g%isc ; ie = g%iec ; nz = gv%ke
1375 diag_wt = dt_in_t * idt_diag
1377 if (cs%omega_frac >= 1.0) absf = 2.0*cs%omega
1379 u_star = fluxes%ustar(i,j)
1380 if (
associated(fluxes%ustar_shelf) .and.
associated(fluxes%frac_shelf_h))
then
1381 if (fluxes%frac_shelf_h(i,j) > 0.0) &
1382 u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
1383 fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
1386 if (u_star < cs%ustar_min) u_star = cs%ustar_min
1387 if (cs%omega_frac < 1.0)
then
1388 absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
1389 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
1390 if (cs%omega_frac > 0.0) &
1391 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1393 absf_ustar = absf / u_star
1394 idecay_len_tke(i) = (absf_ustar * cs%TKE_decay) * gv%H_to_Z
1406 ih = gv%H_to_Z/(3.0*0.41*u_star*dt_in_t)
1407 cmke(1,i) = 4.0 * ih ; cmke(2,i) = (absf_ustar*gv%H_to_Z) * ih
1409 if (idecay_len_tke(i) > 0.0)
then
1410 exp_kh = exp(-htot(i)*idecay_len_tke(i))
1418 if (conv_en(i) < 0.0) conv_en(i) = 0.0
1419 if (ctke(i,1) > 0.0)
then ; tke_ca = ctke(i,1) ;
else ; tke_ca = 0.0 ;
endif
1420 if ((htot(i) >= h_ca(i)) .or. (tke_ca == 0.0))
then
1421 toten_z = us%m_to_Z**2 * (conv_en(i) + tke_ca)
1423 if (toten_z > 0.0)
then
1424 nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1425 sqrt(0.5 * dt_in_t * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
1432 if (conv_en(i) > 0.0)
then
1433 toten_z = us%m_to_Z**2 * (conv_en(i) + tke_ca * (htot(i) / h_ca(i)) )
1434 nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1435 sqrt(0.5 * dt_in_t * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
1440 toten_z = us%m_to_Z**2 * (conv_en(i) + tke_ca)
1441 if (tke_ca > 0.0)
then
1442 nstar_ca = cs%nstar * toten_z / (toten_z + 0.2 * &
1443 sqrt(0.5 * dt_in_t * (absf*(h_ca(i)*gv%H_to_Z))**3 * toten_z))
1449 if (dke_fc(i) + dke_ca(i,1) > 0.0)
then
1450 if (htot(i) >= h_ca(i))
then
1451 mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1452 mke_rate_ca = mke_rate_fc
1454 mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1455 mke_rate_ca = 1.0 / (1.0 + h_ca(i)*(cmke(1,i) + cmke(2,i)*h_ca(i)) )
1459 mke_rate_fc = 1.0 ; mke_rate_ca = 1.0
1462 dke_conv = dke_ca(i,1) * mke_rate_ca + dke_fc(i) * mke_rate_fc
1465 tke(i) = (dt_in_t*cs%mstar)*((us%Z_to_m**2*(u_star*u_star*u_star))*exp_kh) + &
1466 (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca * tke_ca)
1468 if (cs%do_rivermix)
then
1469 tke(i) = tke(i) + tke_river(i)*dt_in_t*exp_kh
1472 if (cs%TKE_diagnostics)
then
1473 wind_tke_src = cs%mstar*(us%Z_to_m**2*u_star*u_star*u_star) * diag_wt
1474 cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + &
1475 ( wind_tke_src + tke_river(i) * diag_wt )
1476 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + dke_conv*idt_diag
1477 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1478 (exp_kh-1.0)*(wind_tke_src + dke_conv*idt_diag)
1479 cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + &
1480 idt_diag * (nstar_fc*conv_en(i) + nstar_ca*tke_ca)
1481 cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + &
1482 idt_diag * ((cs%nstar-nstar_fc)*conv_en(i) + (cs%nstar-nstar_ca)*tke_ca)
1483 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
1484 idt_diag * (ctke(i,1)-tke_ca)
1488 end subroutine find_starting_tke
1491 subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
1492 R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, &
1493 dR0_dT, dRcv_dT, cMKE, Idt_diag, nsw, &
1494 Pen_SW_bnd, opacity_band, TKE, &
1495 Idecay_len_TKE, j, ksort, G, GV, US, CS)
1499 real,
dimension(SZI_(G),SZK_(GV)), &
1501 real,
dimension(SZI_(G),SZK_(GV)), &
1502 intent(inout) :: d_eb
1505 real,
dimension(SZI_(G)),
intent(inout) :: htot
1506 real,
dimension(SZI_(G)),
intent(inout) :: Ttot
1508 real,
dimension(SZI_(G)),
intent(inout) :: Stot
1510 real,
dimension(SZI_(G)),
intent(inout) :: uhtot
1512 real,
dimension(SZI_(G)),
intent(inout) :: vhtot
1514 real,
dimension(SZI_(G)),
intent(inout) :: R0_tot
1516 real,
dimension(SZI_(G)),
intent(inout) :: Rcv_tot
1518 real,
dimension(SZI_(G),SZK_(GV)), &
1520 real,
dimension(SZI_(G),SZK_(GV)), &
1522 real,
dimension(SZI_(G),SZK_(GV)), &
1524 real,
dimension(SZI_(G),SZK_(GV)), &
1526 real,
dimension(SZI_(G),SZK_(GV)), &
1529 real,
dimension(SZI_(G),SZK_(GV)), &
1532 real,
dimension(SZI_(G),SZK_(GV)), &
1535 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
1537 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
1539 real,
dimension(2,SZI_(G)),
intent(in) :: cMKE
1542 real,
intent(in) :: Idt_diag
1544 integer,
intent(in) :: nsw
1546 real,
dimension(max(nsw,1),SZI_(G)),
intent(inout) :: Pen_SW_bnd
1549 real,
dimension(max(nsw,1),SZI_(G),SZK_(GV)),
intent(in) :: opacity_band
1551 real,
dimension(SZI_(G)),
intent(inout) :: TKE
1554 real,
dimension(SZI_(G)),
intent(inout) :: Idecay_len_TKE
1555 integer,
intent(in) :: j
1556 integer,
dimension(SZI_(G),SZK_(GV)), &
1565 real :: Pen_absorbed
1569 real :: h_min, h_max
1579 real :: TKE_full_ent
1583 real :: Pen_En_Contrib
1592 real :: Pen_dTKE_dh_Contrib
1602 real :: f1_x1, f2_x1
1608 real :: C1_3, C1_6, C1_24
1609 integer :: is, ie, nz, i, k, ks, itt, n
1611 c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0 ; c1_24 = 1.0/24.0
1612 g_h_2rho0 = (us%L_to_m**2*gv%g_Earth * gv%H_to_Z) / (2.0 * gv%Rho0)
1613 hmix_min = cs%Hmix_min
1614 h_neglect = gv%H_subroundoff
1615 is = g%isc ; ie = g%iec ; nz = gv%ke
1619 do i=is,ie ;
if (ksort(i,ks) > 0)
then
1622 h_avail = h(i,k) - eps(i,k)
1623 if ((h_avail > 0.) .and. ((tke(i) > 0.) .or. (htot(i) < hmix_min)))
then
1624 drl = g_h_2rho0 * (r0(i,k)*htot(i) - r0_tot(i) )
1625 dmke = (gv%H_to_Z * cs%bulk_Ri_ML) * 0.5 * us%T_to_s**2 * &
1626 ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1629 kh = idecay_len_tke(i)*h_avail ; exp_kh = exp(-kh)
1630 if (kh >= 2.0e-5)
then ; f1_kh = (1.0-exp_kh) / kh
1631 else ; f1_kh = (1.0 - kh*(0.5 - c1_6*kh)) ;
endif
1633 pen_en_contrib = 0.0
1634 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1635 opacity = opacity_band(n,i,k)
1639 if (idecay_len_tke(i) > opacity)
then
1640 x1 = (idecay_len_tke(i) - opacity) * h_avail
1641 if (x1 >= 2.0e-5)
then
1642 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1643 f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1645 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1646 f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1649 pen_en1 = exp(-opacity*h_avail) * &
1650 ((1.0+opacity*htot(i))*f1_x1 + opacity*h_avail*f3_x1)
1652 x1 = (opacity - idecay_len_tke(i)) * h_avail
1653 if (x1 >= 2.0e-5)
then
1654 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1655 f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1657 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1658 f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1661 pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1662 opacity*h_avail*f2_x1)
1664 pen_en_contrib = pen_en_contrib + &
1665 (g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
1668 hpe = htot(i)+h_avail
1669 mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1670 ef4_val = ef4(htot(i)+h_neglect,h_avail,idecay_len_tke(i))
1671 tke_full_ent = (exp_kh*tke(i) - (h_avail*gv%H_to_Z)*(drl*f1_kh + pen_en_contrib)) + &
1672 mke_rate*dmke*ef4_val
1673 if ((tke_full_ent >= 0.0) .or. (h_avail+htot(i) <= hmix_min))
then
1677 if (cs%TKE_diagnostics)
then
1678 e_hxhpe = h_ent / ((htot(i)+h_neglect)*(htot(i)+h_ent+h_neglect))
1679 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1680 idt_diag * ((exp_kh-1.0)*tke(i) + (h_ent*gv%H_to_Z)*drl*(1.0-f1_kh) + &
1681 mke_rate*dmke*(ef4_val-e_hxhpe))
1682 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1683 idt_diag*(gv%H_to_Z*h_ent)*drl
1684 cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1685 idt_diag*(gv%H_to_Z*h_ent)*pen_en_contrib
1686 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1687 idt_diag*mke_rate*dmke*e_hxhpe
1690 tke(i) = tke_full_ent
1692 if (tke(i) <= 0.0) tke(i) = 1.0e-150*us%T_to_s**2*us%m_to_Z
1697 h_min = 0.0 ; h_max = h_avail
1698 if (tke(i) <= 0.0)
then
1701 h_ent = h_avail * tke(i) / (tke(i) - tke_full_ent)
1706 kh = idecay_len_tke(i)*h_ent ; exp_kh = exp(-kh)
1707 if (kh >= 2.0e-5)
then
1708 f1_kh = (1.0-exp_kh) / kh
1710 f1_kh = (1.0 - kh*(0.5 - c1_6*kh))
1714 pen_en_contrib = 0.0 ; pen_dtke_dh_contrib = 0.0
1715 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1719 opacity = opacity_band(n,i,k)
1720 sw_trans = exp(-h_ent*opacity)
1721 if (idecay_len_tke(i) > opacity)
then
1722 x1 = (idecay_len_tke(i) - opacity) * h_ent
1723 if (x1 >= 2.0e-5)
then
1724 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1725 f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1727 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1728 f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1730 pen_en1 = sw_trans * ((1.0+opacity*htot(i))*f1_x1 + &
1731 opacity*h_ent*f3_x1)
1733 x1 = (opacity - idecay_len_tke(i)) * h_ent
1734 if (x1 >= 2.0e-5)
then
1735 e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1736 f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1738 f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1739 f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1742 pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1743 opacity*h_ent*f2_x1)
1745 cpen1 = g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)
1746 pen_en_contrib = pen_en_contrib + cpen1*(pen_en1 - f1_kh)
1747 pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
1748 cpen1*((1.0-sw_trans) - opacity*(htot(i) + h_ent)*sw_trans)
1751 tke_ent1 = exp_kh*tke(i) - (h_ent*gv%H_to_Z)*(drl*f1_kh + pen_en_contrib)
1752 ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i),def4_dh)
1754 mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1755 tke_ent = tke_ent1 + dmke*ef4_val*mke_rate
1758 dtke_dh = ((-idecay_len_tke(i)*tke_ent1 - drl*gv%H_to_Z) + &
1759 pen_dtke_dh_contrib*gv%H_to_Z) + dmke * mke_rate* &
1760 (def4_dh - ef4_val*mke_rate*(cmke(1,i)+2.0*cmke(2,i)*hpe))
1763 if (tke_ent > 0.0)
then
1764 if ((h_max-h_ent)*(-dtke_dh) > tke_ent)
then
1765 dh_newt = -tke_ent / dtke_dh
1767 dh_newt = 0.5*(h_max-h_ent)
1771 if ((h_min-h_ent)*(-dtke_dh) < tke_ent)
then
1772 dh_newt = -tke_ent / dtke_dh
1774 dh_newt = 0.5*(h_min-h_ent)
1778 h_ent = h_ent + dh_newt
1780 if (abs(dh_newt) < 0.2*gv%Angstrom_H)
exit
1784 if (h_ent < hmix_min-htot(i)) h_ent = hmix_min - htot(i)
1786 if (cs%TKE_diagnostics)
then
1788 mke_rate = 1.0/(1.0 + cmke(1,i)*hpe + cmke(2,i)*hpe**2)
1789 ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i))
1791 e_hxhpe = h_ent / ((htot(i)+h_neglect)*(hpe+h_neglect))
1792 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1793 idt_diag * ((exp_kh-1.0)*tke(i) + (h_ent*gv%H_to_Z)*drl*(1.0-f1_kh) + &
1794 dmke*mke_rate*(ef4_val-e_hxhpe))
1795 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1796 idt_diag*(h_ent*gv%H_to_Z)*drl
1797 cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1798 idt_diag*(h_ent*gv%H_to_Z)*pen_en_contrib
1799 cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1800 idt_diag*dmke*mke_rate*e_hxhpe
1807 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
1808 sw_trans = exp(-h_ent*opacity_band(n,i,k))
1809 pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1810 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1813 htot(i) = htot(i) + h_ent
1814 r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1815 h(i,k) = h(i,k) - h_ent
1816 d_eb(i,k) = d_eb(i,k) - h_ent
1818 stot(i) = stot(i) + h_ent * s(i,k)
1819 ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1820 rcv_tot(i) = rcv_tot(i) + (h_ent*rcv(i,k) + pen_absorbed*drcv_dt(i))
1822 uhtot(i) = uhtot(i) + u(i,k)*h_ent
1823 vhtot(i) = vhtot(i) + v(i,k)*h_ent
1829 end subroutine mechanical_entrainment
1833 subroutine sort_ml(h, R0, eps, G, GV, CS, ksort)
1836 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: h
1837 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: R0
1839 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: eps
1843 integer,
dimension(SZI_(G),SZK_(GV)),
intent(out) :: ksort
1846 real :: R0sort(SZI_(G),SZK_(GV))
1847 integer :: nsort(SZI_(G))
1848 logical :: done_sorting(SZI_(G))
1849 integer :: i, k, ks, is, ie, nz, nkmb
1851 is = g%isc ; ie = g%iec ; nz = gv%ke
1852 nkmb = cs%nkml+cs%nkbl
1862 do k=1,nz ;
do i=is,ie ; ksort(i,k) = -1 ;
enddo ;
enddo
1864 do i=is,ie ; nsort(i) = 0 ; done_sorting(i) = .false. ;
enddo
1865 do k=1,nz ;
do i=is,ie ;
if (h(i,k) > eps(i,k))
then
1866 if (done_sorting(i))
then ; ks = nsort(i) ;
else
1868 if (r0(i,k) >= r0sort(i,ks))
exit
1869 r0sort(i,ks+1) = r0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
1871 if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
1875 r0sort(i,ks+1) = r0(i,k)
1876 nsort(i) = nsort(i) + 1
1877 endif ;
enddo ;
enddo
1879 end subroutine sort_ml
1884 subroutine resort_ml(h, T, S, R0, Rcv, RcvTgt, eps, d_ea, d_eb, ksort, G, GV, CS, &
1885 dR0_dT, dR0_dS, dRcv_dT, dRcv_dS)
1889 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: h
1891 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: T
1892 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: S
1893 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: R0
1895 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: Rcv
1897 real,
dimension(SZK_(GV)),
intent(in) :: RcvTgt
1899 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: eps
1901 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_ea
1906 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_eb
1910 integer,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: ksort
1913 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
1917 real,
dimension(SZI_(G)),
intent(in) :: dR0_dS
1921 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
1925 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
1944 real :: h_move, h_tgt_old, I_hnew
1945 real :: dT_dS_wt2, dT_dR, dS_dR, I_denom
1947 real :: T_up, S_up, R0_up, I_hup, h_to_up
1948 real :: T_dn, S_dn, R0_dn, I_hdn, h_to_dn
1951 real :: dPE, hmin, min_dPE, min_hmin
1952 real,
dimension(SZK_(GV)) :: &
1953 h_tmp, R0_tmp, T_tmp, S_tmp, Rcv_tmp
1955 logical :: sorted, leave_in_layer
1956 integer :: ks_deep(SZI_(G)), k_count(SZI_(G)), ks2_reverse(SZI_(G), SZK_(GV))
1957 integer :: ks2(SZK_(GV))
1958 integer :: i, k, ks, is, ie, nz, k1, k2, k_tgt, k_src, k_int_top
1959 integer :: nks, nkmb, num_interior, top_interior_ks
1961 is = g%isc ; ie = g%iec ; nz = gv%ke
1962 nkmb = cs%nkml+cs%nkbl
1964 dt_ds_wt2 = cs%dT_dS_wt**2
1967 do i=is,ie ; ks_deep(i) = -1 ; k_count(i) = 0 ;
enddo
1968 do ks=nz,1,-1 ;
do i=is,ie ;
if (ksort(i,ks) > 0)
then
1971 if (h(i,k) > eps(i,k))
then
1972 if (ks_deep(i) == -1)
then
1974 ks_deep(i) = ks ; k_count(i) = k_count(i) + 1
1975 ks2_reverse(i,k_count(i)) = k
1978 k_count(i) = k_count(i) + 1
1979 ks2_reverse(i,k_count(i)) = k
1982 endif ;
enddo ;
enddo
1984 do i=is,ie ;
if (k_count(i) > 1)
then
1990 ks2(nks) = ks2_reverse(i,1)
1992 ks2(ks) = ks2_reverse(i,1+nks-ks)
1993 if (ks2(ks) > ks2(ks+1)) sorted = .false.
2001 num_interior = 0 ; top_interior_ks = 0
2002 do ks=nks,1,-1 ;
if (ks2(ks) > nkmb)
then
2003 num_interior = num_interior+1 ; top_interior_ks = ks
2006 if (num_interior >= 1)
then
2009 do k=nkmb+1,nz ;
if (rcv(i,0) < rcvtgt(k))
exit ;
enddo
2010 k_int_top = k ; rcv_int = rcvtgt(k)
2012 i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
2013 dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
2014 ds_dr = drcv_ds(i) * i_denom
2020 do ks = nks,top_interior_ks,-1
2022 leave_in_layer = .false.
2023 if ((k > nkmb) .and. (rcv(i,k) <= rcvtgt(k)))
then
2024 if (rcvtgt(k)-rcv(i,k) < cs%BL_split_rho_tol*(rcvtgt(k) - rcvtgt(k-1))) &
2025 leave_in_layer = .true.
2026 elseif (k > nkmb)
then
2027 if (rcv(i,k)-rcvtgt(k) < cs%BL_split_rho_tol*(rcvtgt(k+1) - rcvtgt(k))) &
2028 leave_in_layer = .true.
2031 if (leave_in_layer)
then
2034 elseif (rcv(i,k) < rcv_int)
then
2041 do k2=k_int_top+1,nz ;
if (rcv(i,k) < rcvtgt(k2))
exit ;
enddo
2046 dr1 = (rcvtgt(k2-1) - rcv(i,k)) ; dr2 = (rcvtgt(k2) - rcv(i,k))
2047 t_up = t(i,k) + dt_dr * dr1
2048 s_up = s(i,k) + ds_dr * dr1
2049 t_dn = t(i,k) + dt_dr * dr2
2050 s_dn = s(i,k) + ds_dr * dr2
2052 r0_up = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr1
2053 r0_dn = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr2
2056 if ((r0_up > r0(i,0)) .or. (r0_dn > r0(i,0))) &
2060 wt_dn = (rcv(i,k) - rcvtgt(k2-1)) / (rcvtgt(k2) - rcvtgt(k2-1))
2061 h_to_up = (h(i,k)-eps(i,k)) * (1.0 - wt_dn)
2062 h_to_dn = (h(i,k)-eps(i,k)) * wt_dn
2064 i_hup = 1.0 / (h(i,k2-1) + h_to_up)
2065 i_hdn = 1.0 / (h(i,k2) + h_to_dn)
2066 r0(i,k2-1) = (r0(i,k2)*h(i,k2-1) + r0_up*h_to_up) * i_hup
2067 r0(i,k2) = (r0(i,k2)*h(i,k2) + r0_dn*h_to_dn) * i_hdn
2069 t(i,k2-1) = (t(i,k2)*h(i,k2-1) + t_up*h_to_up) * i_hup
2070 t(i,k2) = (t(i,k2)*h(i,k2) + t_dn*h_to_dn) * i_hdn
2071 s(i,k2-1) = (s(i,k2)*h(i,k2-1) + s_up*h_to_up) * i_hup
2072 s(i,k2) = (s(i,k2)*h(i,k2) + s_dn*h_to_dn) * i_hdn
2073 rcv(i,k2-1) = (rcv(i,k2)*h(i,k2-1) + rcvtgt(k2-1)*h_to_up) * i_hup
2074 rcv(i,k2) = (rcv(i,k2)*h(i,k2) + rcvtgt(k2)*h_to_dn) * i_hdn
2077 h(i,k2) = h(i,k2) + h_to_dn
2078 h(i,k2-1) = h(i,k2-1) + h_to_up
2081 d_eb(i,k) = d_eb(i,k) - h_to_up
2082 d_eb(i,k2-1) = d_eb(i,k2-1) + h_to_up
2083 elseif (k < k2-1)
then
2084 d_ea(i,k) = d_ea(i,k) - h_to_up
2085 d_ea(i,k2-1) = d_ea(i,k2-1) + h_to_up
2088 d_eb(i,k) = d_eb(i,k) - h_to_dn
2089 d_eb(i,k2) = d_eb(i,k2) + h_to_dn
2090 elseif (k < k2)
then
2091 d_ea(i,k) = d_ea(i,k) - h_to_dn
2092 d_ea(i,k2) = d_ea(i,k2) + h_to_dn
2099 do while (nks > nkmb)
2108 ks_min = -1 ; min_dpe = 1.0 ; min_hmin = 0.0
2110 k1 = ks2(ks) ; k2 = ks2(ks+1)
2111 dpe = max(0.0, (r0(i,k2)-r0(i,k1)) * h(i,k1) * h(i,k2))
2112 hmin = min(h(i,k1)-eps(i,k1), h(i,k2)-eps(i,k2))
2113 if ((ks_min < 0) .or. (dpe < min_dpe) .or. &
2114 ((dpe <= 0.0) .and. (hmin < min_hmin)))
then
2115 ks_min = ks ; min_dpe = dpe ; min_hmin = hmin
2120 k1 = ks2(ks_min) ; k2 = ks2(ks_min+1)
2121 if (k1 < k2)
then ; k_tgt = k1 ; k_src = k2
2122 else ; k_tgt = k2 ; k_src = k1 ; ks2(ks_min) = k2 ;
endif
2124 h_tgt_old = h(i,k_tgt)
2125 h_move = h(i,k_src)-eps(i,k_src)
2126 h(i,k_src) = eps(i,k_src)
2127 h(i,k_tgt) = h(i,k_tgt) + h_move
2128 i_hnew = 1.0 / (h(i,k_tgt))
2129 r0(i,k_tgt) = (r0(i,k_tgt)*h_tgt_old + r0(i,k_src)*h_move) * i_hnew
2131 t(i,k_tgt) = (t(i,k_tgt)*h_tgt_old + t(i,k_src)*h_move) * i_hnew
2132 s(i,k_tgt) = (s(i,k_tgt)*h_tgt_old + s(i,k_src)*h_move) * i_hnew
2133 rcv(i,k_tgt) = (rcv(i,k_tgt)*h_tgt_old + rcv(i,k_src)*h_move) * i_hnew
2135 d_eb(i,k_src) = d_eb(i,k_src) - h_move
2136 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2139 do ks=ks_min+1,nks ; ks2(ks) = ks2(ks+1) ;
enddo
2146 do ks=1,nks-1 ;
if (ks2(ks) > ks2(ks+1)) sorted = .false. ;
enddo
2155 h_tmp(k) = h(i,k) ; r0_tmp(k) = r0(i,k)
2156 t_tmp(k) = t(i,k) ; s_tmp(k) = s(i,k) ; rcv_tmp(k) = rcv(i,k)
2162 k_tgt = nkmb - nks + ks ; k_src = ks2(ks)
2163 if (k_tgt == k_src)
then
2164 h(i,k_tgt) = h_tmp(k_tgt)
2169 if (k_src > nkmb)
then
2170 h_move = h(i,k_src)-eps(i,k_src)
2171 h(i,k_src) = eps(i,k_src)
2173 r0(i,k_tgt) = r0(i,k_src)
2175 t(i,k_tgt) = t(i,k_src) ; s(i,k_tgt) = s(i,k_src)
2176 rcv(i,k_tgt) = rcv(i,k_src)
2178 d_eb(i,k_src) = d_eb(i,k_src) - h_move
2179 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2181 h(i,k_tgt) = h_tmp(k_src)
2182 r0(i,k_tgt) = r0_tmp(k_src)
2184 t(i,k_tgt) = t_tmp(k_src) ; s(i,k_tgt) = s_tmp(k_src)
2185 rcv(i,k_tgt) = rcv_tmp(k_src)
2187 if (k_src > k_tgt)
then
2188 d_eb(i,k_src) = d_eb(i,k_src) - h_tmp(k_src)
2189 d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_tmp(k_src)
2191 d_ea(i,k_src) = d_ea(i,k_src) - h_tmp(k_src)
2192 d_ea(i,k_tgt) = d_ea(i,k_tgt) + h_tmp(k_src)
2200 end subroutine resort_ml
2205 subroutine mixedlayer_detrain_2(h, T, S, R0, Rcv, RcvTgt, dt_in_T, dt_diag, d_ea, j, G, GV, US, CS, &
2206 dR0_dT, dR0_dS, dRcv_dT, dRcv_dS, max_BL_det)
2209 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: h
2211 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: T
2212 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: S
2213 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: R0
2215 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: Rcv
2217 real,
dimension(SZK_(GV)),
intent(in) :: RcvTgt
2219 real,
intent(in) :: dt_in_T
2220 real,
intent(in) :: dt_diag
2221 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_ea
2225 integer,
intent(in) :: j
2229 real,
dimension(SZI_(G)),
intent(in) :: dR0_dT
2233 real,
dimension(SZI_(G)),
intent(in) :: dR0_dS
2237 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
2241 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
2244 real,
dimension(SZI_(G)),
intent(in) :: max_BL_det
2272 real :: stays_min, stays_max
2275 logical :: mergeable_bl
2281 real :: stays_min_merge
2283 real :: dR0_2dz, dRcv_2dz
2292 real :: dPE_det, dPE_merge
2301 real :: h_det_to_h2, h_ml_to_h2
2302 real :: h_det_to_h1, h_ml_to_h1
2303 real :: h1_to_h2, h1_to_k0
2304 real :: h2_to_k1, h2_to_k1_rem
2309 real :: R0_det, T_det, S_det
2310 real :: Rcv_stays, R0_stays
2311 real :: T_stays, S_stays
2312 real :: dSpice_det, dSpice_stays
2316 real :: dSpice_lim, dSpice_lim2
2327 real :: dPE_time_ratio
2328 real :: dT_dS_gauge, dS_dT_gauge
2338 logical :: stable_Rcv
2347 real :: Ih, Ihdet, Ih1f, Ih2f
2348 real :: Ihk0, Ihk1, Ih12
2349 real :: dR1, dR2, dR2b, dRk1
2350 real :: dR0, dR21, dRcv
2351 real :: dRcv_stays, dRcv_det, dRcv_lim
2354 real :: h2_to_k1_lim, T_new, S_new, T_max, T_min, S_max, S_min
2355 character(len=200) :: mesg
2357 integer :: i, k, k0, k1, is, ie, nz, kb1, kb2, nkmb
2358 is = g%isc ; ie = g%iec ; nz = gv%ke
2359 kb1 = cs%nkml+1; kb2 = cs%nkml+2
2360 nkmb = cs%nkml+cs%nkbl
2361 h_neglect = gv%H_subroundoff
2362 g_2 = 0.5 * us%L_to_m**2*gv%g_Earth
2363 rho0xg = gv%Rho0 * us%L_to_m**2*gv%g_Earth
2364 idt_h2 = gv%H_to_Z**2 / dt_diag
2365 i2rho0 = 0.5 / gv%Rho0
2366 angstrom = gv%Angstrom_H
2369 dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 / dt_ds_gauge
2372 if (cs%nkbl /= 2)
call mom_error(fatal,
"MOM_mixed_layer"// &
2373 "CS%nkbl must be 2 in mixedlayer_detrain_2.")
2375 if (dt_in_t < cs%BL_detrain_time)
then ; dpe_time_ratio = cs%BL_detrain_time / (dt_in_t)
2376 else ; dpe_time_ratio = 1.0 ;
endif
2385 h_to_bl = 0.0 ; r0_to_bl = 0.0
2386 rcv_to_bl = 0.0 ; t_to_bl = 0.0 ; s_to_bl = 0.0
2388 do k=1,cs%nkml ;
if (h(i,k) > 0.0)
then
2389 h_to_bl = h_to_bl + h(i,k)
2390 r0_to_bl = r0_to_bl + r0(i,k)*h(i,k)
2392 rcv_to_bl = rcv_to_bl + rcv(i,k)*h(i,k)
2393 t_to_bl = t_to_bl + t(i,k)*h(i,k)
2394 s_to_bl = s_to_bl + s(i,k)*h(i,k)
2396 d_ea(i,k) = d_ea(i,k) - h(i,k)
2399 if (h_to_bl > 0.0)
then ; r0_det = r0_to_bl / h_to_bl
2400 else ; r0_det = r0(i,0) ;
endif
2422 h_min_bl = min(cs%Hbuffer_min, cs%Hbuffer_rel_min*h(i,0))
2425 if (((r0(i,kb2)-r0(i,kb1)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) &
2426 stable_rcv = .false.
2428 h1 = h(i,kb1) ; h2 = h(i,kb2)
2430 h2_to_k1_rem = (h1 + h2) + h_to_bl
2431 if ((max_bl_det(i) >= 0.0) .and. (h2_to_k1_rem > max_bl_det(i))) &
2432 h2_to_k1_rem = max_bl_det(i)
2435 if ((h2 == 0.0) .and. (h1 > 0.0))
then
2444 do k1=kb2+1,nz ;
if (h(i,k1) > 2.0*angstrom)
exit ;
enddo
2446 r0(i,kb2) = r0(i,kb1)
2448 rcv(i,kb2)=rcv(i,kb1) ; t(i,kb2)=t(i,kb1) ; s(i,kb2)=s(i,kb1)
2451 if (k1 <= nz)
then ;
if (r0(i,k1) >= r0(i,kb1))
then
2452 r0(i,kb2) = 0.5*(r0(i,kb1)+r0(i,k1))
2454 rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
2455 t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
2456 s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
2460 dpe_extrap = 0.0 ; dpe_merge = 0.0
2461 mergeable_bl = .false.
2462 if ((h1 > 0.0) .and. (h2 > 0.0) .and. (h_to_bl > 0.0) .and. &
2469 do k1=kb2+1,nz ;
if (rcvtgt(k1) >= rcv(i,kb2))
exit ;
enddo ; k0 = k1-1
2470 dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2476 h1_avail = h1 - max(0.0,h_min_bl-h_to_bl)
2477 if ((k1<=nz) .and. (h2 > h_min_bl) .and. (h1_avail > 0.0) .and. &
2478 (r0(i,kb1) < r0(i,kb2)) .and. (h_to_bl*r0(i,kb1) > r0_to_bl))
then
2479 drk1 = (rcvtgt(k1) - rcv(i,kb2)) * (r0(i,kb2) - r0(i,kb1)) / &
2480 (rcv(i,kb2) - rcv(i,kb1))
2481 b1 = drk1 / (r0(i,kb2) - r0(i,kb1))
2487 h2_to_k1 = min(h2 - h_min_bl, h2_to_k1_rem)
2490 if (h2_to_k1*(h1_avail + b1*(h1_avail + h2)) > h2*h1_avail) &
2491 h2_to_k1 = (h2*h1_avail) / (h1_avail + b1*(h1_avail + h2))
2492 if (h2_to_k1*(drk1 * h2) > (h_to_bl*r0(i,kb1) - r0_to_bl) * h1) &
2493 h2_to_k1 = (h_to_bl*r0(i,kb1) - r0_to_bl) * h1 / (drk1 * h2)
2495 if ((k1==kb2+1) .and. (cs%BL_extrap_lim > 0.))
then
2498 drcv_lim = rcv(i,kb2)-rcv(i,0)
2499 do k=1,kb2 ; drcv_lim = max(drcv_lim, rcv(i,kb2)-rcv(i,k)) ;
enddo
2500 drcv_lim = cs%BL_extrap_lim*drcv_lim
2501 if ((rcvtgt(k1) - rcv(i,kb2)) >= drcv_lim)
then
2503 elseif ((rcvtgt(k1) - rcv(i,kb2)) > 0.5*drcv_lim)
then
2504 h2_to_k1 = h2_to_k1 * (2.0 - 2.0*((rcvtgt(k1) - rcv(i,kb2)) / drcv_lim))
2508 drcv = (rcvtgt(k1) - rcv(i,kb2))
2513 dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2514 dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2515 (h2 - h2_to_k1) / (h1 + h2)
2517 if (h(i,k1) > 10.0*angstrom)
then
2518 dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2519 dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2520 if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2522 if (k1<nz)
then ;
if (h(i,k1+1) > 10.0*angstrom)
then
2523 dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2524 dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2525 if ((dspice_det*dspice_lim2 > 0.0) .and. &
2526 (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2528 if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2530 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2531 t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2532 (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2533 s_det = s(i,kb2) + i_denom * &
2534 (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2536 r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2537 (s_det-s(i,kb2)) * dr0_ds(i)
2539 if (cs%BL_extrap_lim >= 0.)
then
2542 if (h(i,k1) > 10.0*angstrom)
then
2543 t_min = min(t(i,kb1), t(i,kb2), t(i,k1)) - cs%Allowed_T_chg
2544 t_max = max(t(i,kb1), t(i,kb2), t(i,k1)) + cs%Allowed_T_chg
2545 s_min = min(s(i,kb1), s(i,kb2), s(i,k1)) - cs%Allowed_S_chg
2546 s_max = max(s(i,kb1), s(i,kb2), s(i,k1)) + cs%Allowed_S_chg
2548 t_min = min(t(i,kb1), t(i,kb2)) - cs%Allowed_T_chg
2549 t_max = max(t(i,kb1), t(i,kb2)) + cs%Allowed_T_chg
2550 s_min = min(s(i,kb1), s(i,kb2)) - cs%Allowed_S_chg
2551 s_max = max(s(i,kb1), s(i,kb2)) + cs%Allowed_S_chg
2553 ihk1 = 1.0 / (h(i,k1) + h2_to_k1)
2554 t_new = (h(i,k1)*t(i,k1) + h2_to_k1*t_det) * ihk1
2555 s_new = (h(i,k1)*s(i,k1) + h2_to_k1*s_det) * ihk1
2557 if ((t_new < t_min) .or. (t_new > t_max) .or. &
2558 (s_new < s_min) .or. (s_new > s_max)) &
2562 h1_to_h2 = b1*h2*h2_to_k1 / (h2 - (1.0+b1)*h2_to_k1)
2564 ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2565 ih2f = 1.0 / ((h(i,kb2) - h2_to_k1) + h1_to_h2)
2567 rcv(i,kb2) = ((h(i,kb2)*rcv(i,kb2) - h2_to_k1*rcvtgt(k1)) + &
2568 h1_to_h2*rcv(i,kb1))*ih2f
2569 rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2571 t(i,kb2) = ((h(i,kb2)*t(i,kb2) - h2_to_k1*t_det) + &
2572 h1_to_h2*t(i,kb1)) * ih2f
2573 t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2575 s(i,kb2) = ((h(i,kb2)*s(i,kb2) - h2_to_k1*s_det) + &
2576 h1_to_h2*s(i,kb1)) * ih2f
2577 s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2580 r0(i,kb2) = ((h(i,kb2)*r0(i,kb2) - h2_to_k1*r0_det) + &
2581 h1_to_h2*r0(i,kb1)) * ih2f
2582 r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2584 h(i,kb1) = h(i,kb1) - h1_to_h2 ; h1 = h(i,kb1)
2585 h(i,kb2) = (h(i,kb2) - h2_to_k1) + h1_to_h2 ; h2 = h(i,kb2)
2586 h(i,k1) = h(i,k1) + h2_to_k1
2588 d_ea(i,kb1) = d_ea(i,kb1) - h1_to_h2
2589 d_ea(i,kb2) = (d_ea(i,kb2) - h2_to_k1) + h1_to_h2
2590 d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2591 h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2595 if ((k1>kb2+1) .and. (rcvtgt(k1-1) >= rcv(i,kb2)))
then
2596 do k1=k1,kb2+1,-1 ;
if (rcvtgt(k1-1) < rcv(i,kb2))
exit ;
enddo
2601 dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2603 if ((k0>kb2) .and. (dr1 > 0.0) .and. (h1 > h_min_bl) .and. &
2604 (h2*dr2 < h1*dr1) .and. (r0(i,kb2) > r0(i,kb1)))
then
2609 stays_merge = 2.0*(h1+h2)*(h1*dr1 - h2*dr2) / &
2610 ((dr1+dr2)*h1 + dr1*(h1+h2) + &
2611 sqrt((dr2*h1-dr1*h2)**2 + 4*(h1+h2)*h2*(dr1+dr2)*dr2))
2613 stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
2614 h1 - (h1+h2)*(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1)))
2615 if ((stays_merge > stays_min_merge) .and. &
2616 (stays_merge + h2_to_k1_rem >= h1 + h2))
then
2617 mergeable_bl = .true.
2618 dpe_merge = g_2*(r0(i,kb2)-r0(i,kb1))*(h1-stays_merge)*(h2-stays_merge)
2622 if ((k1<=nz).and.(.not.mergeable_bl))
then
2626 dr2b = rcvtgt(k1)-rcv(i,kb2) ; dr21 = rcv(i,kb2) - rcv(i,kb1)
2627 if (dr2b*(h1+h2) < h2*dr21)
then
2629 h2_to_k1 = min(h2 - (h1+h2) * dr2b / dr21, h2_to_k1_rem)
2631 if (h2 > h2_to_k1)
then
2632 drcv = (rcvtgt(k1) - rcv(i,kb2))
2637 dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2638 dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2639 (h2 - h2_to_k1) / (h1 + h2)
2641 if (h(i,k1) > 10.0*angstrom)
then
2642 dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2643 dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2644 if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2646 if (k1<nz) then;
if (h(i,k1+1) > 10.0*angstrom)
then
2647 dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2648 dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2649 if ((dspice_det*dspice_lim2 > 0.0) .and. &
2650 (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2652 if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2654 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2655 t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2656 (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2657 s_det = s(i,kb2) + i_denom * &
2658 (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2660 r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2661 (s_det-s(i,kb2)) * dr0_ds(i)
2666 ih2f = 1.0 / (h2 - h2_to_k1)
2668 t_min = min(t(i,kb2), t(i,kb1)) - cs%Allowed_T_chg
2669 t_max = max(t(i,kb2), t(i,kb1)) + cs%Allowed_T_chg
2670 t_new = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2671 if (t_new < t_min)
then
2672 h2_to_k1_lim = h2 * (t(i,kb2) - t_min) / (t_det - t_min)
2679 h2_to_k1 = h2_to_k1_lim
2680 ih2f = 1.0 / (h2 - h2_to_k1)
2681 elseif (t_new > t_max)
then
2682 h2_to_k1_lim = h2 * (t(i,kb2) - t_max) / (t_det - t_max)
2689 h2_to_k1 = h2_to_k1_lim
2690 ih2f = 1.0 / (h2 - h2_to_k1)
2692 s_min = max(min(s(i,kb2), s(i,kb1)) - cs%Allowed_S_chg, 0.0)
2693 s_max = max(s(i,kb2), s(i,kb1)) + cs%Allowed_S_chg
2694 s_new = (h2*s(i,kb2) - h2_to_k1*s_det)*ih2f
2695 if (s_new < s_min)
then
2696 h2_to_k1_lim = h2 * (s(i,kb2) - s_min) / (s_det - s_min)
2703 h2_to_k1 = h2_to_k1_lim
2704 ih2f = 1.0 / (h2 - h2_to_k1)
2705 elseif (s_new > s_max)
then
2706 h2_to_k1_lim = h2 * (s(i,kb2) - s_max) / (s_det - s_max)
2713 h2_to_k1 = h2_to_k1_lim
2714 ih2f = 1.0 / (h2 - h2_to_k1)
2717 ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2718 rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2719 rcv(i,kb2) = rcv(i,kb2) - h2_to_k1*drcv*ih2f
2721 t(i,kb2) = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2722 t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2724 s(i,kb2) = (h2*s(i,kb2) - h2_to_k1*s_det) * ih2f
2725 s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2728 r0(i,kb2) = (h2*r0(i,kb2) - h2_to_k1*r0_det) * ih2f
2729 r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2734 ihk1 = 1.0 / (h(i,k1) + h2)
2736 rcv(i,k1) = (h(i,k1)*rcv(i,k1) + h2*rcv(i,kb2)) * ihk1
2737 t(i,k1) = (h(i,k1)*t(i,k1) + h2*t(i,kb2)) * ihk1
2738 s(i,k1) = (h(i,k1)*s(i,k1) + h2*s(i,kb2)) * ihk1
2739 r0(i,k1) = (h(i,k1)*r0(i,k1) + h2*r0(i,kb2)) * ihk1
2742 h(i,k1) = h(i,k1) + h2_to_k1
2743 h(i,kb2) = h(i,kb2) - h2_to_k1 ; h2 = h(i,kb2)
2745 dpe_extrap = i2rho0*(r0_det-r0(i,kb2))*h2_to_k1*h2
2747 d_ea(i,kb2) = d_ea(i,kb2) - h2_to_k1
2748 d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2749 h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2756 h_det_h2 = max(h_min_bl-(h1+h2), 0.0)
2757 if (h_det_h2 > 0.0)
then
2763 h_det_to_h2 = min(h_to_bl, h_det_h2)
2764 h_ml_to_h2 = h_det_h2 - h_det_to_h2
2765 h_det_to_h1 = h_to_bl - h_det_to_h2
2766 h_ml_to_h1 = max(h_min_bl-h_det_to_h1,0.0)
2769 ihdet = 0.0 ;
if (h_to_bl > 0.0) ihdet = 1.0 / h_to_bl
2770 ih1f = 1.0 / (h_det_to_h1 + h_ml_to_h1)
2772 r0(i,kb2) = ((h2*r0(i,kb2) + h1*r0(i,kb1)) + &
2773 (h_det_to_h2*r0_to_bl*ihdet + h_ml_to_h2*r0(i,0))) * ih
2774 r0(i,kb1) = (h_det_to_h1*r0_to_bl*ihdet + h_ml_to_h1*r0(i,0)) * ih1f
2776 rcv(i,kb2) = ((h2*rcv(i,kb2) + h1*rcv(i,kb1)) + &
2777 (h_det_to_h2*rcv_to_bl*ihdet + h_ml_to_h2*rcv(i,0))) * ih
2778 rcv(i,kb1) = (h_det_to_h1*rcv_to_bl*ihdet + h_ml_to_h1*rcv(i,0)) * ih1f
2780 t(i,kb2) = ((h2*t(i,kb2) + h1*t(i,kb1)) + &
2781 (h_det_to_h2*t_to_bl*ihdet + h_ml_to_h2*t(i,0))) * ih
2782 t(i,kb1) = (h_det_to_h1*t_to_bl*ihdet + h_ml_to_h1*t(i,0)) * ih1f
2784 s(i,kb2) = ((h2*s(i,kb2) + h1*s(i,kb1)) + &
2785 (h_det_to_h2*s_to_bl*ihdet + h_ml_to_h2*s(i,0))) * ih
2786 s(i,kb1) = (h_det_to_h1*s_to_bl*ihdet + h_ml_to_h1*s(i,0)) * ih1f
2789 d_ea(i,1) = d_ea(i,1) - (h_ml_to_h1 + h_ml_to_h2)
2790 d_ea(i,kb1) = d_ea(i,kb1) + ((h_det_to_h1 + h_ml_to_h1) - h1)
2791 d_ea(i,kb2) = d_ea(i,kb2) + (h_min_bl - h2)
2793 h(i,kb1) = h_det_to_h1 + h_ml_to_h1 ; h(i,kb2) = h_min_bl
2794 h(i,0) = h(i,0) - (h_ml_to_h1 + h_ml_to_h2)
2797 if (
allocated(cs%diag_PE_detrain) .or.
allocated(cs%diag_PE_detrain2))
then
2798 r0_det = r0_to_bl*ihdet
2799 s1en = g_2 * idt_h2 * ( ((r0(i,kb2)-r0(i,kb1))*h1*h2 + &
2800 h_det_to_h2*( (r0(i,kb1)-r0_det)*h1 + (r0(i,kb2)-r0_det)*h2 ) + &
2801 h_ml_to_h2*( (r0(i,kb2)-r0(i,0))*h2 + (r0(i,kb1)-r0(i,0))*h1 + &
2802 (r0_det-r0(i,0))*h_det_to_h2 ) + &
2803 h_det_to_h1*h_ml_to_h1*(r0_det-r0(i,0))) - 2.0*gv%Rho0*dpe_extrap )
2805 if (
allocated(cs%diag_PE_detrain)) &
2806 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + s1en
2808 if (
allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
2809 cs%diag_PE_detrain2(i,j) + s1en + idt_h2*rho0xg*dpe_extrap
2812 elseif ((h_to_bl > 0.0) .or. (h1 < h_min_bl) .or. (h2 < h_min_bl))
then
2817 h_from_ml = h_min_bl + max(h_min_bl-h2,0.0) - h1 - h_to_bl
2818 if (h_from_ml > 0.0)
then
2821 dpe_extrap = dpe_extrap - i2rho0*h_from_ml*(r0_to_bl - r0(i,0)*h_to_bl)
2822 r0_to_bl = r0_to_bl + h_from_ml*r0(i,0)
2823 rcv_to_bl = rcv_to_bl + h_from_ml*rcv(i,0)
2824 t_to_bl = t_to_bl + h_from_ml*t(i,0)
2825 s_to_bl = s_to_bl + h_from_ml*s(i,0)
2827 h_to_bl = h_to_bl + h_from_ml
2828 h(i,0) = h(i,0) - h_from_ml
2829 d_ea(i,1) = d_ea(i,1) - h_from_ml
2834 if (r0(i,kb2) - r0(i,kb1) > 1.0e-9*abs(r0(i,kb1) - r0_det)) &
2835 b1 = abs(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1))
2836 stays_min = max((1.0-b1)*h1 - b1*h2, 0.0, h_min_bl - h_to_bl)
2837 stays_max = h1 - max(h_min_bl-h2,0.0)
2840 if (stays_max <= stays_min)
then
2842 mergeable_bl = .false.
2843 if (stays_max < h1) scale_slope = (h1 - stays_min) / (h1 - stays_max)
2848 i_ya = (h1 + h2) / ((h1 + h2) + h_to_bl)
2850 s1 = 0.5*(h1 + (h2 - bh0) * i_ya) ; s2 = h1 - s1
2855 s3sq = i_ya*max(bh0*h1-dpe_extrap, 0.0)
2857 s3sq = i_ya*(bh0*h1-min(dpe_extrap,0.0))
2860 if (s3sq == 0.0)
then
2863 elseif (s2*s2 <= s3sq)
then
2871 if (bh0 <= 0.0)
then ; stays = h1
2872 elseif (s2 > 0.0)
then
2874 if (s1 >= stays_max)
then ; stays = stays_max
2875 elseif (s1 >= 0.0)
then ; stays = s1 + sqrt(s2*s2 - s3sq)
2876 else ; stays = (h1*(s2-s1) - s3sq) / (-s1 + sqrt(s2*s2 - s3sq))
2880 if (s1 <= stays_min)
then ; stays = stays_min
2881 else ; stays = (h1*(s1-s2) + s3sq) / (s1 + sqrt(s2*s2 - s3sq))
2890 if (stays >= stays_max)
then ; stays = stays_max
2891 elseif (stays < stays_min)
then ; stays = stays_min
2895 dpe_det = g_2*((r0(i,kb1)*h_to_bl - r0_to_bl)*stays + &
2896 (r0(i,kb2)-r0(i,kb1)) * (h1-stays) * &
2897 (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
2900 if (dpe_time_ratio*h_to_bl > h_to_bl+h(i,0))
then
2901 dpe_ratio = (h_to_bl+h(i,0)) / h_to_bl
2903 dpe_ratio = dpe_time_ratio
2906 if ((mergeable_bl) .and. (num_events*dpe_ratio*dpe_det > dpe_merge))
then
2911 h1_to_k0 = (h1-stays_merge)
2912 stays = max(h_min_bl-h_to_bl,0.0)
2913 h1_to_h2 = stays_merge - stays
2915 ihk0 = 1.0 / ((h1_to_k0 + h2) + h(i,k0))
2916 ih1f = 1.0 / (h_to_bl + stays); ih2f = 1.0 / h1_to_h2
2917 ih12 = 1.0 / (h1 + h2)
2919 drcv_2dz = (rcv(i,kb1) - rcv(i,kb2)) * ih12
2920 drcv_stays = drcv_2dz*(h1_to_k0 + h1_to_h2)
2921 drcv_det = - drcv_2dz*(stays + h1_to_h2)
2922 rcv(i,k0) = ((h1_to_k0*(rcv(i,kb1) + drcv_det) + &
2923 h2*rcv(i,kb2)) + h(i,k0)*rcv(i,k0)) * ihk0
2924 rcv(i,kb2) = rcv(i,kb1) + drcv_2dz*(h1_to_k0-stays)
2925 rcv(i,kb1) = (rcv_to_bl + stays*(rcv(i,kb1) + drcv_stays)) * ih1f
2930 i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2931 dspice_2dz = (ds_dt_gauge*drcv_ds(i)*(t(i,kb1)-t(i,kb2)) - &
2932 dt_ds_gauge*drcv_dt(i)*(s(i,kb1)-s(i,kb2))) * ih12
2933 dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
2934 dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
2935 if (dspice_lim * dspice_2dz <= 0.0) dspice_2dz = 0.0
2937 if (stays > 0.0)
then
2939 if (abs(dspice_lim) < abs(dspice_2dz*(h1_to_k0 + h1_to_h2))) &
2940 dspice_2dz = dspice_lim/(h1_to_k0 + h1_to_h2)
2942 dspice_stays = dspice_2dz*(h1_to_k0 + h1_to_h2)
2943 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
2944 (dt_ds_gauge * drcv_dt(i) * drcv_stays + drcv_ds(i) * dspice_stays)
2945 s_stays = s(i,kb1) + i_denom * &
2946 (drcv_ds(i) * drcv_stays - dt_ds_gauge * drcv_dt(i) * dspice_stays)
2948 r0_stays = r0(i,kb1) + (t_stays-t(i,kb1)) * dr0_dt(i) + &
2949 (s_stays-s(i,kb1)) * dr0_ds(i)
2952 if (abs(dspice_lim) < abs(dspice_2dz*h1_to_k0)) &
2953 dspice_2dz = dspice_lim/h1_to_k0
2955 t_stays = 0.0 ; s_stays = 0.0 ; r0_stays = 0.0
2958 dspice_det = - dspice_2dz*(stays + h1_to_h2)
2959 t_det = t(i,kb1) + dt_ds_gauge * i_denom * &
2960 (dt_ds_gauge * drcv_dt(i) * drcv_det + drcv_ds(i) * dspice_det)
2961 s_det = s(i,kb1) + i_denom * &
2962 (drcv_ds(i) * drcv_det - dt_ds_gauge * drcv_dt(i) * dspice_det)
2964 r0_det = r0(i,kb1) + (t_det-t(i,kb1)) * dr0_dt(i) + &
2965 (s_det-s(i,kb1)) * dr0_ds(i)
2967 t(i,k0) = ((h1_to_k0*t_det + h2*t(i,kb2)) + h(i,k0)*t(i,k0)) * ihk0
2968 t(i,kb2) = (h1*t(i,kb1) - stays*t_stays - h1_to_k0*t_det) * ih2f
2969 t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
2971 s(i,k0) = ((h1_to_k0*s_det + h2*s(i,kb2)) + h(i,k0)*s(i,k0)) * ihk0
2972 s(i,kb2) = (h1*s(i,kb1) - stays*s_stays - h1_to_k0*s_det) * ih2f
2973 s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
2975 r0(i,k0) = ((h1_to_k0*r0_det + h2*r0(i,kb2)) + h(i,k0)*r0(i,k0)) * ihk0
2976 r0(i,kb2) = (h1*r0(i,kb1) - stays*r0_stays - h1_to_k0*r0_det) * ih2f
2977 r0(i,kb1) = (r0_to_bl + stays*r0_stays) * ih1f
2999 d_ea(i,kb1) = (d_ea(i,kb1) + h_to_bl) + (stays - h1)
3000 d_ea(i,kb2) = d_ea(i,kb2) + (h1_to_h2 - h2)
3001 d_ea(i,k0) = d_ea(i,k0) + (h1_to_k0 + h2)
3003 h(i,kb1) = stays + h_to_bl
3005 h(i,k0) = h(i,k0) + (h1_to_k0 + h2)
3006 if (
allocated(cs%diag_PE_detrain)) &
3007 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_merge
3008 if (
allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3009 cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap)
3014 h1_to_h2 = h1 - stays
3015 ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / (h2 + h1_to_h2)
3016 ih = 1.0 / (h1 + h2)
3017 dr0_2dz = (r0(i,kb1) - r0(i,kb2)) * ih
3018 r0(i,kb2) = (h2*r0(i,kb2) + h1_to_h2*(r0(i,kb1) - &
3019 scale_slope*dr0_2dz*stays)) * ih2f
3020 r0(i,kb1) = (r0_to_bl + stays*(r0(i,kb1) + &
3021 scale_slope*dr0_2dz*h1_to_h2)) * ih1f
3026 dr0 = scale_slope*dr0_2dz*h1_to_h2
3027 dspice_stays = (ds_dt_gauge*dr0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3028 dt_ds_gauge*dr0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
3029 scale_slope*h1_to_h2 * ih
3030 if (h_to_bl > 0.0)
then
3031 dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3032 dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) /&
3035 dspice_lim = ds_dt_gauge*dr0_ds(i)*(t(i,0)-t(i,kb1)) - &
3036 dt_ds_gauge*dr0_dt(i)*(s(i,0)-s(i,kb1))
3038 if (dspice_stays*dspice_lim <= 0.0)
then
3040 elseif (abs(dspice_stays) > abs(dspice_lim))
then
3041 dspice_stays = dspice_lim
3043 i_denom = 1.0 / (dr0_ds(i)**2 + (dt_ds_gauge*dr0_dt(i))**2)
3044 t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3045 (dt_ds_gauge * dr0_dt(i) * dr0 + dr0_ds(i) * dspice_stays)
3046 s_stays = s(i,kb1) + i_denom * &
3047 (dr0_ds(i) * dr0 - dt_ds_gauge * dr0_dt(i) * dspice_stays)
3049 rcv_stays = rcv(i,kb1) + (t_stays-t(i,kb1)) * drcv_dt(i) + &
3050 (s_stays-s(i,kb1)) * drcv_ds(i)
3052 t(i,kb2) = (h2*t(i,kb2) + h1*t(i,kb1) - t_stays*stays) * ih2f
3053 t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3054 s(i,kb2) = (h2*s(i,kb2) + h1*s(i,kb1) - s_stays*stays) * ih2f
3055 s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3056 rcv(i,kb2) = (h2*rcv(i,kb2) + h1*rcv(i,kb1) - rcv_stays*stays) * ih2f
3057 rcv(i,kb1) = (rcv_to_bl + stays*rcv_stays) * ih1f
3076 d_ea(i,kb1) = d_ea(i,kb1) + ((stays - h1) + h_to_bl)
3077 d_ea(i,kb2) = d_ea(i,kb2) + h1_to_h2
3079 h(i,kb1) = stays + h_to_bl
3080 h(i,kb2) = h(i,kb2) + h1_to_h2
3082 if (
allocated(cs%diag_PE_detrain)) &
3083 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_det
3084 if (
allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3085 cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap)
3091 end subroutine mixedlayer_detrain_2
3096 subroutine mixedlayer_detrain_1(h, T, S, R0, Rcv, RcvTgt, dt_in_T, dt_diag, d_ea, d_eb, &
3097 j, G, GV, US, CS, dRcv_dT, dRcv_dS, max_BL_det)
3100 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: h
3102 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: T
3103 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: S
3104 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: R0
3106 real,
dimension(SZI_(G),SZK0_(GV)),
intent(inout) :: Rcv
3108 real,
dimension(SZK_(GV)),
intent(in) :: RcvTgt
3110 real,
intent(in) :: dt_in_T
3111 real,
intent(in) :: dt_diag
3113 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_ea
3117 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: d_eb
3121 integer,
intent(in) :: j
3125 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dT
3129 real,
dimension(SZI_(G)),
intent(in) :: dRcv_dS
3132 real,
dimension(SZI_(G)),
intent(in) :: max_BL_det
3140 real :: max_det_rem(SZI_(G))
3141 real :: detrain(SZI_(G))
3143 real :: dT_dR, dS_dR, dRml, dR0_dRcv, dT_dS_wt2
3145 real :: Sdown, Tdown
3147 real :: g_H2_2Rho0dt
3154 logical :: splittable_BL(SZI_(G)), orthogonal_extrap
3157 integer :: i, is, ie, k, k1, nkmb, nz
3158 is = g%isc ; ie = g%iec ; nz = gv%ke
3159 nkmb = cs%nkml+cs%nkbl
3160 if (cs%nkbl /= 1)
call mom_error(fatal,
"MOM_mixed_layer: "// &
3161 "CS%nkbl must be 1 in mixedlayer_detrain_1.")
3163 dt_time = dt_in_t / cs%BL_detrain_time
3164 g_h2_2rho0dt = (us%L_to_m**2*gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0 * dt_diag)
3165 g_h2_2dt = (us%L_to_m**2*gv%g_Earth * gv%H_to_Z**2) / (2.0 * dt_diag)
3169 do i=is,ie ;
if (h(i,k) > 0.0)
then
3170 ih = 1.0 / (h(i,nkmb) + h(i,k))
3171 if (cs%TKE_diagnostics) &
3172 cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
3173 g_h2_2rho0dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3174 if (
allocated(cs%diag_PE_detrain)) &
3175 cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + &
3176 g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3177 if (
allocated(cs%diag_PE_detrain2)) &
3178 cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + &
3179 g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3181 r0(i,nkmb) = (r0(i,nkmb)*h(i,nkmb) + r0(i,k)*h(i,k)) * ih
3182 rcv(i,nkmb) = (rcv(i,nkmb)*h(i,nkmb) + rcv(i,k)*h(i,k)) * ih
3183 t(i,nkmb) = (t(i,nkmb)*h(i,nkmb) + t(i,k)*h(i,k)) * ih
3184 s(i,nkmb) = (s(i,nkmb)*h(i,nkmb) + s(i,k)*h(i,k)) * ih
3186 d_ea(i,k) = d_ea(i,k) - h(i,k)
3187 d_ea(i,nkmb) = d_ea(i,nkmb) + h(i,k)
3188 h(i,nkmb) = h(i,nkmb) + h(i,k)
3194 max_det_rem(i) = 10.0 * h(i,nkmb)
3195 if (max_bl_det(i) >= 0.0) max_det_rem(i) = max_bl_det(i)
3210 do i=is,ie ;
if (h(i,nkmb) > 0.0)
then
3211 if ((r0(i,0)<r0(i,nz)) .and. (r0(i,nz)<r0(i,nkmb)))
then
3212 if ((r0(i,nz)-r0(i,0))*h(i,0) > &
3213 (r0(i,nkmb)-r0(i,nz))*h(i,nkmb))
then
3214 detrain(i) = (r0(i,nkmb)-r0(i,nz))*h(i,nkmb) / (r0(i,nkmb)-r0(i,0))
3216 detrain(i) = (r0(i,nz)-r0(i,0))*h(i,0) / (r0(i,nkmb)-r0(i,0))
3219 d_eb(i,cs%nkml) = d_eb(i,cs%nkml) + detrain(i)
3220 d_ea(i,cs%nkml) = d_ea(i,cs%nkml) - detrain(i)
3221 d_eb(i,nkmb) = d_eb(i,nkmb) - detrain(i)
3222 d_ea(i,nkmb) = d_ea(i,nkmb) + detrain(i)
3224 if (
allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3225 cs%diag_PE_detrain(i,j) + g_h2_2dt * detrain(i)* &
3226 (h(i,0) + h(i,nkmb)) * (r0(i,nkmb) - r0(i,0))
3228 r0(i,0) = r0(i,0) - detrain(i)*(r0(i,0)-r0(i,nkmb)) / h(i,0)
3229 r0(i,nkmb) = r0(i,nkmb) - detrain(i)*(r0(i,nkmb)-x1) / h(i,nkmb)
3231 rcv(i,0) = rcv(i,0) - detrain(i)*(rcv(i,0)-rcv(i,nkmb)) / h(i,0)
3232 rcv(i,nkmb) = rcv(i,nkmb) - detrain(i)*(rcv(i,nkmb)-x1) / h(i,nkmb)
3234 t(i,0) = t(i,0) - detrain(i)*(t(i,0)-t(i,nkmb)) / h(i,0)
3235 t(i,nkmb) = t(i,nkmb) - detrain(i)*(t(i,nkmb)-x1) / h(i,nkmb)
3237 s(i,0) = s(i,0) - detrain(i)*(s(i,0)-s(i,nkmb)) / h(i,0)
3238 s(i,nkmb) = s(i,nkmb) - detrain(i)*(s(i,nkmb)-x1) / h(i,nkmb)
3248 if (h(i,nkmb) > 0.0)
then ; splittable_bl(i) = .true.
3249 else ; splittable_bl(i) = .false. ;
endif
3252 dt_ds_wt2 = cs%dT_dS_wt**2
3254 do k=nz-1,nkmb+1,-1 ;
do i=is,ie
3255 if (splittable_bl(i))
then
3256 if (rcvtgt(k)<=rcv(i,nkmb))
then
3261 splittable_bl(i) = .false.
3263 k1 = k+1 ; orthogonal_extrap = .false.
3268 if ((h(i,k+1) < 10.0*gv%Angstrom_H) .or. &
3269 ((rcvtgt(k+1)-rcv(i,nkmb)) >= 0.9*(rcv(i,k1) - rcv(i,0))))
then
3270 if (k>=nz-1)
then ; orthogonal_extrap = .true.
3271 elseif ((h(i,k+2) <= 10.0*gv%Angstrom_H) .and. &
3272 ((rcvtgt(k+1)-rcv(i,nkmb)) < 0.9*(rcv(i,k+2)-rcv(i,0))))
then
3274 else ; orthogonal_extrap = .true. ;
endif
3277 if ((r0(i,0) >= r0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
3281 if (orthogonal_extrap)
then
3285 i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
3286 dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
3287 ds_dr = drcv_ds(i) * i_denom
3289 dt_dr = (t(i,0) - t(i,k1)) / (rcv(i,0) - rcv(i,k1))
3290 ds_dr = (s(i,0) - s(i,k1)) / (rcv(i,0) - rcv(i,k1))
3292 drml = dt_time * (r0(i,nkmb) - r0(i,0)) * &
3293 (rcv(i,0) - rcv(i,k1)) / (r0(i,0) - r0(i,k1))
3295 if (drml < 0.0) cycle
3296 dr0_drcv = (r0(i,0) - r0(i,k1)) / (rcv(i,0) - rcv(i,k1))
3298 if ((rcv(i,nkmb) - drml < rcvtgt(k)) .and. (max_det_rem(i) > h(i,nkmb)))
then
3300 detrain(i) = h(i,nkmb)*(rcv(i,nkmb) - rcvtgt(k)) / &
3301 (rcvtgt(k+1) - rcvtgt(k))
3303 if (
allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3304 cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * &
3305 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dr0_drcv
3307 tdown = detrain(i) * (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3308 t(i,k) = (h(i,k) * t(i,k) + &
3309 (h(i,nkmb) * t(i,nkmb) - tdown)) / &
3310 (h(i,k) + (h(i,nkmb) - detrain(i)))
3311 t(i,k+1) = (h(i,k+1) * t(i,k+1) + tdown)/ &
3312 (h(i,k+1) + detrain(i))
3314 sdown = detrain(i) * (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3315 s(i,k) = (h(i,k) * s(i,k) + &
3316 (h(i,nkmb) * s(i,nkmb) - sdown)) / &
3317 (h(i,k) + (h(i,nkmb) - detrain(i)))
3318 s(i,k+1) = (h(i,k+1) * s(i,k+1) + sdown)/ &
3319 (h(i,k+1) + detrain(i))
3321 rcv(i,nkmb) = rcv(i,0)
3323 d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3324 d_ea(i,k) = d_ea(i,k) + (h(i,nkmb) - detrain(i))
3325 d_ea(i,nkmb) = d_ea(i,nkmb) - h(i,nkmb)
3327 h(i,k+1) = h(i,k+1) + detrain(i)
3328 h(i,k) = h(i,k) + h(i,nkmb) - detrain(i)
3332 detrain(i) = h(i,nkmb) * drml / (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3333 if (detrain(i) > max_det_rem(i)) detrain(i) = max_det_rem(i)
3334 ih = 1.0 / (h(i,k+1) + detrain(i))
3336 tdown = (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3337 t(i,nkmb) = t(i,nkmb) - dt_dr * drml
3338 t(i,k+1) = (h(i,k+1) * t(i,k+1) + detrain(i) * tdown) * ih
3339 sdown = (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3343 s(i,nkmb) = s(i,nkmb) - ds_dr * drml
3344 s(i,k+1) = (h(i,k+1) * s(i,k+1) + detrain(i) * sdown) * ih
3346 d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3347 d_ea(i,nkmb) = d_ea(i,nkmb) - detrain(i)
3349 h(i,k+1) = h(i,k+1) + detrain(i)
3350 h(i,nkmb) = h(i,nkmb) - detrain(i)
3352 if (
allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3353 cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * dr0_drcv * &
3354 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3365 if (h(i,nkmb) < 0.1*h(i,0))
then
3366 h_ent = 0.1*h(i,0) - h(i,nkmb)
3368 t(i,nkmb) = (h(i,nkmb)*t(i,nkmb) + h_ent*t(i,0)) * ih
3369 s(i,nkmb) = (h(i,nkmb)*s(i,nkmb) + h_ent*s(i,0)) * ih
3371 d_ea(i,1) = d_ea(i,1) - h_ent
3372 d_ea(i,nkmb) = d_ea(i,nkmb) + h_ent
3374 h(i,0) = h(i,0) - h_ent
3375 h(i,nkmb) = h(i,nkmb) + h_ent
3379 end subroutine mixedlayer_detrain_1
3382 subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
3383 type(time_type),
target,
intent(in) :: time
3389 type(
diag_ctrl),
target,
intent(inout) :: diag
3402 #include "version_variable.h"
3403 character(len=40) :: mdl =
"MOM_mixed_layer"
3404 real :: bl_detrain_time_dflt
3405 real :: omega_frac_dflt, ustar_min_dflt, hmix_min_m
3406 integer :: isd, ied, jsd, jed
3407 logical :: use_temperature, use_omega
3408 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3410 if (
associated(cs))
then
3411 call mom_error(warning,
"mixedlayer_init called with an associated"// &
3412 "associated control structure.")
3414 else ;
allocate(cs) ;
endif
3419 if (gv%nkml < 1)
return
3425 call log_param(param_file, mdl,
"NKML", cs%nkml, &
3426 "The number of sublayers within the mixed layer if "//&
3427 "BULKMIXEDLAYER is true.", units=
"nondim", default=2)
3428 cs%nkbl = gv%nk_rho_varies - gv%nkml
3429 call log_param(param_file, mdl,
"NKBL", cs%nkbl, &
3430 "The number of variable density buffer layers if "//&
3431 "BULKMIXEDLAYER is true.", units=
"nondim", default=2)
3432 call get_param(param_file, mdl,
"MSTAR", cs%mstar, &
3433 "The ratio of the friction velocity cubed to the TKE "//&
3434 "input to the mixed layer.",
"units=nondim", default=1.2)
3435 call get_param(param_file, mdl,
"NSTAR", cs%nstar, &
3436 "The portion of the buoyant potential energy imparted by "//&
3437 "surface fluxes that is available to drive entrainment "//&
3438 "at the base of mixed layer when that energy is positive.", &
3439 units=
"nondim", default=0.15)
3440 call get_param(param_file, mdl,
"BULK_RI_ML", cs%bulk_Ri_ML, &
3441 "The efficiency with which mean kinetic energy released "//&
3442 "by mechanically forced entrainment of the mixed layer "//&
3443 "is converted to turbulent kinetic energy.", units=
"nondim",&
3444 fail_if_missing=.true.)
3445 call get_param(param_file, mdl,
"ABSORB_ALL_SW", cs%absorb_all_sw, &
3446 "If true, all shortwave radiation is absorbed by the "//&
3447 "ocean, instead of passing through to the bottom mud.", &
3449 call get_param(param_file, mdl,
"TKE_DECAY", cs%TKE_decay, &
3450 "TKE_DECAY relates the vertical rate of decay of the "//&
3451 "TKE available for mechanical entrainment to the natural "//&
3452 "Ekman depth.", units=
"nondim", default=2.5)
3453 call get_param(param_file, mdl,
"NSTAR2", cs%nstar2, &
3454 "The portion of any potential energy released by "//&
3455 "convective adjustment that is available to drive "//&
3456 "entrainment at the base of mixed layer. By default "//&
3457 "NSTAR2=NSTAR.", units=
"nondim", default=cs%nstar)
3458 call get_param(param_file, mdl,
"BULK_RI_CONVECTIVE", cs%bulk_Ri_convective, &
3459 "The efficiency with which convectively released mean "//&
3460 "kinetic energy is converted to turbulent kinetic "//&
3461 "energy. By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
3462 units=
"nondim", default=cs%bulk_Ri_ML)
3463 call get_param(param_file, mdl,
"HMIX_MIN", cs%Hmix_min, &
3464 "The minimum mixed layer depth if the mixed layer depth "//&
3465 "is determined dynamically.", units=
"m", default=0.0, scale=gv%m_to_H, &
3466 unscaled=hmix_min_m)
3468 call get_param(param_file, mdl,
"LIMIT_BUFFER_DETRAIN", cs%limit_det, &
3469 "If true, limit the detrainment from the buffer layers "//&
3470 "to not be too different from the neighbors.", default=.false.)
3471 call get_param(param_file, mdl,
"ALLOWED_DETRAIN_TEMP_CHG", cs%Allowed_T_chg, &
3472 "The amount by which temperature is allowed to exceed "//&
3473 "previous values during detrainment.", units=
"K", default=0.5)
3474 call get_param(param_file, mdl,
"ALLOWED_DETRAIN_SALT_CHG", cs%Allowed_S_chg, &
3475 "The amount by which salinity is allowed to exceed "//&
3476 "previous values during detrainment.", units=
"PSU", default=0.1)
3477 call get_param(param_file, mdl,
"ML_DT_DS_WEIGHT", cs%dT_dS_wt, &
3478 "When forced to extrapolate T & S to match the layer "//&
3479 "densities, this factor (in deg C / PSU) is combined "//&
3480 "with the derivatives of density with T & S to determine "//&
3481 "what direction is orthogonal to density contours. It "//&
3482 "should be a typical value of (dR/dS) / (dR/dT) in "//&
3483 "oceanic profiles.", units=
"degC PSU-1", default=6.0)
3484 call get_param(param_file, mdl,
"BUFFER_LAYER_EXTRAP_LIMIT", cs%BL_extrap_lim, &
3485 "A limit on the density range over which extrapolation "//&
3486 "can occur when detraining from the buffer layers, "//&
3487 "relative to the density range within the mixed and "//&
3488 "buffer layers, when the detrainment is going into the "//&
3489 "lightest interior layer, nondimensional, or a negative "//&
3490 "value not to apply this limit.", units=
"nondim", default = -1.0)
3491 call get_param(param_file, mdl,
"BUFFER_LAYER_HMIN_THICK", cs%Hbuffer_min, &
3492 "The minimum buffer layer thickness when the mixed layer is very thick.", &
3493 units=
"m", default=5.0, scale=gv%m_to_H)
3494 call get_param(param_file, mdl,
"BUFFER_LAYER_HMIN_REL", cs%Hbuffer_rel_min, &
3495 "The minimum buffer layer thickness relative to the combined mixed "//&
3496 "land buffer ayer thicknesses when they are thin.", &
3497 units=
"nondim", default=0.1/cs%nkbl)
3498 bl_detrain_time_dflt = 4.0*3600.0 ;
if (cs%nkbl==1) bl_detrain_time_dflt = 86400.0*30.0
3499 call get_param(param_file, mdl,
"BUFFER_LAY_DETRAIN_TIME", cs%BL_detrain_time, &
3500 "A timescale that characterizes buffer layer detrainment events.", &
3501 units=
"s", default=bl_detrain_time_dflt, scale=us%s_to_T)
3502 call get_param(param_file, mdl,
"BUFFER_SPLIT_RHO_TOL", cs%BL_split_rho_tol, &
3503 "The fractional tolerance for matching layer target densities when splitting "//&
3504 "layers to deal with massive interior layers that are lighter than one of the "//&
3505 "mixed or buffer layers.", units=
"nondim", default=0.1)
3507 call get_param(param_file, mdl,
"DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
3508 "The surface fluxes are scaled away when the total ocean "//&
3509 "depth is less than DEPTH_LIMIT_FLUXES.", &
3510 units=
"m", default=0.1*hmix_min_m, scale=gv%m_to_H)
3511 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
3512 "The rotation rate of the earth.", &
3513 default=7.2921e-5, units=
"s-1", scale=us%T_to_s)
3514 call get_param(param_file, mdl,
"ML_USE_OMEGA", use_omega, &
3515 "If true, use the absolute rotation rate instead of the "//&
3516 "vertical component of rotation when setting the decay "//&
3517 "scale for turbulence.", default=.false., do_not_log=.true.)
3518 omega_frac_dflt = 0.0
3520 call mom_error(warning,
"ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
3521 omega_frac_dflt = 1.0
3523 call get_param(param_file, mdl,
"ML_OMEGA_FRAC", cs%omega_frac, &
3524 "When setting the decay scale for turbulence, use this "//&
3525 "fraction of the absolute rotation rate blended with the "//&
3526 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
3527 units=
"nondim", default=omega_frac_dflt)
3528 call get_param(param_file, mdl,
"ML_RESORT", cs%ML_resort, &
3529 "If true, resort the topmost layers by potential density "//&
3530 "before the mixed layer calculations.", default=.false.)
3532 call get_param(param_file, mdl,
"ML_PRESORT_NK_CONV_ADJ", cs%ML_presort_nz_conv_adj, &
3533 "Convectively mix the first ML_PRESORT_NK_CONV_ADJ "//&
3534 "layers before sorting when ML_RESORT is true.", &
3535 units=
"nondim", default=0, fail_if_missing=.true.)
3537 ustar_min_dflt = 2e-4*us%s_to_T*cs%omega*(gv%Angstrom_m + gv%H_to_m*gv%H_subroundoff)
3538 call get_param(param_file, mdl,
"BML_USTAR_MIN", cs%ustar_min, &
3539 "The minimum value of ustar that should be used by the "//&
3540 "bulk mixed layer model in setting vertical TKE decay "//&
3541 "scales. This must be greater than 0.", units=
"m s-1", &
3542 default=ustar_min_dflt, scale=us%m_to_Z*us%T_to_s)
3543 if (cs%ustar_min<=0.0)
call mom_error(fatal,
"BML_USTAR_MIN must be positive.")
3545 call get_param(param_file, mdl,
"RESOLVE_EKMAN", cs%Resolve_Ekman, &
3546 "If true, the NKML>1 layers in the mixed layer are "//&
3547 "chosen to optimally represent the impact of the Ekman "//&
3548 "transport on the mixed layer TKE budget. Otherwise, "//&
3549 "the sublayers are distributed uniformly through the "//&
3550 "mixed layer.", default=.false.)
3551 call get_param(param_file, mdl,
"CORRECT_ABSORPTION_DEPTH", cs%correct_absorption, &
3552 "If true, the average depth at which penetrating shortwave "//&
3553 "radiation is absorbed is adjusted to match the average "//&
3554 "heating depth of an exponential profile by moving some "//&
3555 "of the heating upward in the water column.", default=.false.)
3556 call get_param(param_file, mdl,
"DO_RIVERMIX", cs%do_rivermix, &
3557 "If true, apply additional mixing wherever there is "//&
3558 "runoff, so that it is mixed down to RIVERMIX_DEPTH, "//&
3559 "if the ocean is that deep.", default=.false.)
3560 if (cs%do_rivermix) &
3561 call get_param(param_file, mdl,
"RIVERMIX_DEPTH", cs%rivermix_depth, &
3562 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
3563 "defined.", units=
"m", default=0.0, scale=us%m_to_Z)
3564 call get_param(param_file, mdl,
"USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
3565 "If true, use the fluxes%runoff_Hflx field to set the "//&
3566 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
3568 call get_param(param_file, mdl,
"USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
3569 "If true, use the fluxes%calving_Hflx field to set the "//&
3570 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
3573 call get_param(param_file, mdl,
"ALLOW_CLOCKS_IN_OMP_LOOPS", &
3574 cs%allow_clocks_in_omp_loops, &
3575 "If true, clocks can be called from inside loops that can "//&
3576 "be threaded. To run with multiple threads, set to False.", &
3579 cs%id_ML_depth = register_diag_field(
'ocean_model',
'h_ML', diag%axesT1, &
3580 time,
'Surface mixed layer depth',
'm')
3581 cs%id_TKE_wind = register_diag_field(
'ocean_model',
'TKE_wind', diag%axesT1, &
3582 time,
'Wind-stirring source of mixed layer TKE',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**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',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**3)
3585 cs%id_TKE_conv = register_diag_field(
'ocean_model',
'TKE_conv', diag%axesT1, &
3586 time,
'Convective source of mixed layer TKE',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**3)
3587 cs%id_TKE_pen_SW = register_diag_field(
'ocean_model',
'TKE_pen_SW', diag%axesT1, &
3588 time,
'TKE consumed by mixing penetrative shortwave radation through the mixed layer', &
3589 'm3 s-3', conversion=us%Z_to_m)
3590 cs%id_TKE_mixing = register_diag_field(
'ocean_model',
'TKE_mixing', diag%axesT1, &
3591 time,
'TKE consumed by mixing that deepens the mixed layer',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**3)
3592 cs%id_TKE_mech_decay = register_diag_field(
'ocean_model',
'TKE_mech_decay', diag%axesT1, &
3593 time,
'Mechanical energy decay sink of mixed layer TKE',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**3)
3594 cs%id_TKE_conv_decay = register_diag_field(
'ocean_model',
'TKE_conv_decay', diag%axesT1, &
3595 time,
'Convective energy decay sink of mixed layer TKE',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**3)
3596 cs%id_TKE_conv_s2 = register_diag_field(
'ocean_model',
'TKE_conv_s2', diag%axesT1, &
3597 time,
'Spurious source of mixed layer TKE from sigma2',
'm3 s-3', conversion=us%Z_to_m*us%T_to_s**3)
3598 cs%id_PE_detrain = register_diag_field(
'ocean_model',
'PE_detrain', diag%axesT1, &
3599 time,
'Spurious source of potential energy from mixed layer detrainment', &
3600 'W m-2', conversion=us%Z_to_m*us%T_to_s**3)
3601 cs%id_PE_detrain2 = register_diag_field(
'ocean_model',
'PE_detrain2', diag%axesT1, &
3602 time,
'Spurious source of potential energy from mixed layer only detrainment', &
3603 'W m-2', conversion=us%Z_to_m*us%T_to_s**3)
3604 cs%id_h_mismatch = register_diag_field(
'ocean_model',
'h_miss_ML', diag%axesT1, &
3605 time,
'Summed absolute mismatch in entrainment terms',
'm', conversion=us%Z_to_m)
3606 cs%id_Hsfc_used = register_diag_field(
'ocean_model',
'Hs_used', diag%axesT1, &
3607 time,
'Surface region thickness that is used',
'm', conversion=us%Z_to_m)
3608 cs%id_Hsfc_max = register_diag_field(
'ocean_model',
'Hs_max', diag%axesT1, &
3609 time,
'Maximum surface region thickness',
'm', conversion=us%Z_to_m)
3610 cs%id_Hsfc_min = register_diag_field(
'ocean_model',
'Hs_min', diag%axesT1, &
3611 time,
'Minimum surface region thickness',
'm', conversion=us%Z_to_m)
3613 if (cs%limit_det .or. (cs%id_Hsfc_min > 0))
then
3614 call get_param(param_file, mdl,
"LIMIT_BUFFER_DET_DH_SFC", cs%lim_det_dH_sfc, &
3615 "The fractional limit in the change between grid points "//&
3616 "of the surface region (mixed & buffer layer) thickness.", &
3617 units=
"nondim", default=0.5)
3618 call get_param(param_file, mdl,
"LIMIT_BUFFER_DET_DH_BATHY", cs%lim_det_dH_bathy, &
3619 "The fraction of the total depth by which the thickness "//&
3620 "of the surface region (mixed & buffer layer) is allowed "//&
3621 "to change between grid points.", units=
"nondim", default=0.2)
3624 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", use_temperature, &
3625 "If true, temperature and salinity are used as state "//&
3626 "variables.", default=.true.)
3628 if (use_temperature)
then
3629 call get_param(param_file, mdl,
"PEN_SW_NBANDS", cs%nsw, default=1)
3633 if (max(cs%id_TKE_wind, cs%id_TKE_RiBulk, cs%id_TKE_conv, cs%id_TKE_mixing, &
3634 cs%id_TKE_pen_SW, cs%id_TKE_mech_decay, cs%id_TKE_conv_decay) > 0)
then
3635 call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
3636 call safe_alloc_alloc(cs%diag_TKE_RiBulk, isd, ied, jsd, jed)
3637 call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
3638 call safe_alloc_alloc(cs%diag_TKE_pen_SW, isd, ied, jsd, jed)
3639 call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
3640 call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
3641 call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
3642 call safe_alloc_alloc(cs%diag_TKE_conv_s2, isd, ied, jsd, jed)
3644 cs%TKE_diagnostics = .true.
3646 if (cs%id_PE_detrain > 0)
call safe_alloc_alloc(cs%diag_PE_detrain, isd, ied, jsd, jed)
3647 if (cs%id_PE_detrain2 > 0)
call safe_alloc_alloc(cs%diag_PE_detrain2, isd, ied, jsd, jed)
3648 if (cs%id_ML_depth > 0)
call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
3650 if (cs%allow_clocks_in_omp_loops)
then
3651 id_clock_detrain = cpu_clock_id(
'(Ocean mixed layer detrain)', grain=clock_routine)
3652 id_clock_mech = cpu_clock_id(
'(Ocean mixed layer mechanical entrainment)', grain=clock_routine)
3653 id_clock_conv = cpu_clock_id(
'(Ocean mixed layer convection)', grain=clock_routine)
3654 if (cs%ML_resort)
then
3655 id_clock_resort = cpu_clock_id(
'(Ocean mixed layer resorting)', grain=clock_routine)
3657 id_clock_adjustment = cpu_clock_id(
'(Ocean mixed layer convective adjustment)', grain=clock_routine)
3659 id_clock_eos = cpu_clock_id(
'(Ocean mixed layer EOS)', grain=clock_routine)
3662 if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) &
3663 id_clock_pass = cpu_clock_id(
'(Ocean mixed layer halo updates)', grain=clock_routine)
3669 end subroutine bulkmixedlayer_init
3676 function ef4(Ht, En, I_L, dR_de)
3677 real,
intent(in) :: ht
3678 real,
intent(in) :: en
3679 real,
intent(in) :: i_l
3680 real,
optional,
intent(inout) :: dr_de
3688 exp_lhpe = exp(-i_l*(en+ht))
3690 res = exp_lhpe * (en*i_hpe/ht - 0.5*i_l*log(ht*i_hpe) + 0.5*i_l*i_l*en)
3691 if (
PRESENT(dr_de)) &
3692 dr_de = -i_l*res + exp_lhpe*(i_hpe*i_hpe + 0.5*i_l*i_hpe + 0.5*i_l*i_l)