This subroutine moves any water left in the former mixed layers into the two buffer layers and may also move buffer layer water into the interior isopycnal layers.
2207 type(ocean_grid_type),
intent(in) :: G
2208 type(verticalGrid_type),
intent(in) :: GV
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
2226 type(unit_scale_type),
intent(in) :: US
2227 type(bulkmixedlayer_CS),
pointer :: CS
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)