phylmd/CV30_routines/cv30_unsat.f

module cv30_unsat_m

  implicit none

contains

  SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, &
       ep, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, evap, b)

    ! Unsaturated (precipitating) downdrafts

    use cv30_param_m, only: nl, sigd
    use cv_thermo_m, only: cpd, ginv, grav

    integer, intent(in):: icb(:) ! (ncum)
    ! {2 <= icb <= nl - 3}
    ! {ph(i, icb(i) + 1) < plcl(i) <= ph(i, icb(i))}

    integer, intent(in):: inb(:) ! (ncum)
    ! first model level above the level of neutral buoyancy of the
    ! parcel (1 <= inb <= nl - 1)

    real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (ncum, nl)
    real, intent(in):: gz(:, :) ! (klon, klev)
    real, intent(in):: u(:, :), v(:, :) ! (klon, klev)
    real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa
    real, intent(in):: ph(:, :) ! (ncum, klev + 1)
    real, intent(in):: th(:, :) ! (ncum, nl - 1)
    real, intent(in):: tv(:, :) ! (klon, klev)
    real, intent(in):: lv(:, :) ! (klon, klev)
    real, intent(in):: cpn(:, :) ! (klon, klev)
    real, intent(in):: ep(:, :) ! (ncum, klev)
    real, intent(in):: clw(:, :) ! (ncum, klev)
    real, intent(in):: m(:, :) ! (ncum, klev)
    real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)
    real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)
    real, intent(in):: delt
    real, intent(in):: plcl(:) ! (ncum)

    real, intent(out):: mp(:, :) ! (klon, klev)
    real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
    real, intent(out):: wt(:, :) ! (ncum, nl)
    real, intent(out):: water(:, :), evap(:, :) ! (ncum, nl)
    real, intent(out):: b(:, :) ! (ncum, nl - 1)

    ! Local:

    real, parameter:: sigp = 0.15
    ! fraction of precipitation falling outside of cloud, \sig_s in
    ! Emanuel (1991 928)

    integer ncum
    integer i, il, imax
    real tinv, delti
    real afac, afac1, afac2, bfac
    real pr1, pr2, sigt, b6, c6, revap, tevap, delth
    real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
    real ampmax
    real lvcp(size(icb), nl) ! (ncum, nl)
    real wdtrain(size(icb)) ! (ncum)
    logical lwork(size(icb)) ! (ncum)

    !------------------------------------------------------

    ncum = size(icb)
    delti = 1. / delt
    tinv = 1. / 3.
    mp = 0.
    b = 0.

    do i = 1, nl
       do il = 1, ncum
          qp(il, i) = q(il, i)
          up(il, i) = u(il, i)
          vp(il, i) = v(il, i)
          wt(il, i) = 0.001
          water(il, i) = 0.
          evap(il, i) = 0.
          lvcp(il, i) = lv(il, i) / cpn(il, i)
       enddo
    enddo

    ! Check whether ep(inb) = 0. If so, skip precipitating downdraft
    ! calculation.

    forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4

    imax = nl - 1
    do while (.not. any(inb >= imax .and. lwork) .and. imax >= 1)
       imax = imax - 1
    end do

    downdraft_loop: DO i = imax, 1, - 1
       ! Integrate liquid water equation to find condensed water 
       ! and condensed water flux 

       ! Calculate detrained precipitation 
       forall (il = 1:ncum, inb(il) >= i .and. lwork(il)) wdtrain(il) = grav &
            * (ep(il, i) * m(il, i) * clw(il, i) &
            + sum(max(elij(il, :i - 1, i) - (1. - ep(il, i)) * clw(il, i), 0.) &
            * ment(il, :i - 1, i)))

       ! Find rain water and evaporation using provisional 
       ! estimates of qp(i) and qp(i - 1) 

       loop_horizontal: do il = 1, ncum
          if (i <= inb(il) .and. lwork(il)) then
             wt(il, i) = 45.

             if (i < inb(il)) then
                qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) - t(il, i)) &
                     + gz(il, i + 1) - gz(il, i)) / lv(il, i)
                qp(il, i) = 0.5 * (qp(il, i) + q(il, i))
             endif

             qp(il, i) = max(qp(il, i), 0.)
             qp(il, i) = min(qp(il, i), qs(il, i))
             qp(il, inb(il)) = q(il, inb(il))

             if (i == 1) then
                afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &
                     / (1e4 + 2000. * p(il, 1) * qs(il, 1))
             else
                qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) - t(il, i - 1)) &
                     + gz(il, i) - gz(il, i - 1)) / lv(il, i)
                qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))
                qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))
                qp(il, i - 1) = max(qp(il, i - 1), 0.)
                afac1 = p(il, i) * (qs(il, i) - qp(il, i)) &
                     / (1e4 + 2000. * p(il, i) * qs(il, i))
                afac2 = p(il, i - 1) * (qs(il, i - 1) - qp(il, i - 1)) &
                     / (1e4 + 2000. * p(il, i - 1) * qs(il, i - 1))
                afac = 0.5 * (afac1 + afac2)
             endif

             if (i == inb(il)) afac = 0.
             afac = max(afac, 0.)
             bfac = 1. / (sigd * wt(il, i))

             ! Prise en compte de la variation progressive de sigt dans
             ! les couches icb et icb - 1:
             ! pour plcl <= ph(i + 1), pr1 = 0 et pr2 = 1
             ! pour plcl >= ph(i), pr1 = 1 et pr2 = 0
             ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion \`a cheval
             ! sur le nuage, et pr2 est la proportion sous la base du
             ! nuage.
             pr1 = max(0., min(1., &
                  (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))))
             pr2 = max(0., min(1., &
                  (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))))
             sigt = sigp * pr1 + pr2

             b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt * afac
             c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
                  * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)

             if (c6 > 0.) then
                revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))
                evap(il, i) = sigt * afac * revap
                water(il, i) = revap * revap
             else
                evap(il, i) = - evap(il, i + 1) + 0.02 * (wdtrain(il) + sigd &
                     * wt(il, i) * water(il, i + 1)) / (sigd * (ph(il, i) &
                     - ph(il, i + 1)))
             end if

             ! Calculate precipitating downdraft mass flux under 
             ! hydrostatic approximation 

             test_above_surface: if (i /= 1) then
                tevap = max(0., evap(il, i))
                delth = max(0.001, (th(il, i) - th(il, i - 1)))
                mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
                     * (p(il, i - 1) - p(il, i)) / delth

                ! If hydrostatic assumption fails, solve cubic
                ! difference equation for downdraft theta and mass
                ! flux from two simultaneous differential equations

                amfac = sigd * sigd * 70. * ph(il, i) &
                     * (p(il, i - 1) - p(il, i)) &
                     * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
                amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
                     * mp(il, i))

                if (amp2 > 0.1 * amfac) then
                   xf = 100. * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))
                   tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
                        * t(il, i) / (lvcp(il, i) * sigd * th(il, i))
                   af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
                   bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
                        * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
                        - p(il, i)) * xf * tevap
                   fac2 = 1.
                   if (bf < 0.) fac2 = - 1.
                   bf = abs(bf)
                   ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv

                   if (ur >= 0.) then
                      sru = sqrt(ur)
                      fac = 1.
                      if ((0.5 * bf - sru) < 0.) fac = - 1.
                      mp(il, i) = mp(il, i + 1) * tinv + (0.5 * bf &
                           + sru)**tinv + fac * (abs(0.5 * bf - sru))**tinv
                   else
                      d = atan(2. * sqrt(- ur) / (bf + 1e-28))
                      if (fac2 < 0.)d = 3.14159 - d
                      mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) &
                           * cos(d * tinv)
                   endif

                   mp(il, i) = max(0., mp(il, i))

                   ! Il y a vraisemblablement une erreur dans la
                   ! ligne suivante : il faut diviser par (mp(il,
                   ! i) * sigd * grav) et non par (mp(il, i) + sigd
                   ! * 0.1).  Et il faut bien revoir les facteurs
                   ! 100.
                   b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &
                        * tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &
                        - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &
                        * th(il, i))
                   b(il, i - 1) = max(b(il, i - 1), 0.)
                endif

                ! Limit magnitude of mp(i) to meet CFL condition:
                ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
                amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
                ampmax = min(ampmax, amp2)
                mp(il, i) = min(mp(il, i), ampmax)

                ! Force mp to decrease linearly to zero between cloud
                ! base and the surface:
                if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
                     * (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))
             endif test_above_surface

             ! Find mixing ratio of precipitating downdraft 

             if (i /= inb(il)) then
                qp(il, i) = q(il, i)

                if (mp(il, i) > mp(il, i + 1)) then
                   qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &
                        * (mp(il, i) - mp(il, i + 1)) + 100. * ginv * 0.5 &
                        * sigd * (ph(il, i) - ph(il, i + 1)) &
                        * (evap(il, i + 1) + evap(il, i))
                   qp(il, i) = qp(il, i) / mp(il, i)
                   up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
                        * (mp(il, i) - mp(il, i + 1))
                   up(il, i) = up(il, i) / mp(il, i)
                   vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &
                        * (mp(il, i) - mp(il, i + 1))
                   vp(il, i) = vp(il, i) / mp(il, i)
                else
                   if (mp(il, i + 1) > 1e-16) then
                      qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &
                           * (ph(il, i) - ph(il, i + 1)) * (evap(il, i + 1) &
                           + evap(il, i)) / mp(il, i + 1)
                      up(il, i) = up(il, i + 1)
                      vp(il, i) = vp(il, i + 1)
                   endif
                endif

                qp(il, i) = min(qp(il, i), qs(il, i))
                qp(il, i) = max(qp(il, i), 0.)
             endif
          endif
       end do loop_horizontal
    end DO downdraft_loop

  end SUBROUTINE cv30_unsat

end module cv30_unsat_m
1	guez	185	module cv30_unsat_m
2	guez	47
3	guez	97	implicit none
4	guez	47
5	guez	97	contains
6	guez	47
7	guez	189	SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, &
8	guez	195	ep, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, evap, b)
9	guez	47
10	guez	195	! Unsaturated (precipitating) downdrafts
11
12	guez	186	use cv30_param_m, only: nl, sigd
13	guez	190	use cv_thermo_m, only: cpd, ginv, grav
14	guez	47
15	guez	192	integer, intent(in):: icb(:) ! (ncum)
16	guez	196	! {2 <= icb <= nl - 3}
17			! {ph(i, icb(i) + 1) < plcl(i) <= ph(i, icb(i))}
18	guez	192
19			integer, intent(in):: inb(:) ! (ncum)
20			! first model level above the level of neutral buoyancy of the
21			! parcel (1 <= inb <= nl - 1)
22
23	guez	195	real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (ncum, nl)
24	guez	189	real, intent(in):: gz(:, :) ! (klon, klev)
25			real, intent(in):: u(:, :), v(:, :) ! (klon, klev)
26	guez	195	real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa
27	guez	196	real, intent(in):: ph(:, :) ! (ncum, klev + 1)
28	guez	195	real, intent(in):: th(:, :) ! (ncum, nl - 1)
29			real, intent(in):: tv(:, :) ! (klon, klev)
30			real, intent(in):: lv(:, :) ! (klon, klev)
31			real, intent(in):: cpn(:, :) ! (klon, klev)
32			real, intent(in):: ep(:, :) ! (ncum, klev)
33			real, intent(in):: clw(:, :) ! (ncum, klev)
34	guez	192	real, intent(in):: m(:, :) ! (ncum, klev)
35			real, intent(in):: ment(:, :, :) ! (ncum, klev, klev)
36			real, intent(in):: elij(:, :, :) ! (ncum, klev, klev)
37	guez	97	real, intent(in):: delt
38	guez	196	real, intent(in):: plcl(:) ! (ncum)
39	guez	47
40	guez	195	real, intent(out):: mp(:, :) ! (klon, klev)
41	guez	189	real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
42	guez	195	real, intent(out):: wt(:, :) ! (ncum, nl)
43			real, intent(out):: water(:, :), evap(:, :) ! (ncum, nl)
44	guez	189	real, intent(out):: b(:, :) ! (ncum, nl - 1)
45	guez	47
46	guez	186	! Local:
47	guez	195
48			real, parameter:: sigp = 0.15
49			! fraction of precipitation falling outside of cloud, \sig_s in
50			! Emanuel (1991 928)
51
52	guez	188	integer ncum
53	guez	194	integer i, il, imax
54	guez	97	real tinv, delti
55	guez	194	real afac, afac1, afac2, bfac
56	guez	97	real pr1, pr2, sigt, b6, c6, revap, tevap, delth
57			real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
58			real ampmax
59	guez	195	real lvcp(size(icb), nl) ! (ncum, nl)
60	guez	193	real wdtrain(size(icb)) ! (ncum)
61			logical lwork(size(icb)) ! (ncum)
62	guez	47
63	guez	97	!------------------------------------------------------
64	guez	47
65	guez	188	ncum = size(icb)
66	guez	186	delti = 1. / delt
67			tinv = 1. / 3.
68			mp = 0.
69	guez	189	b = 0.
70	guez	47
71	guez	186	do i = 1, nl
72			do il = 1, ncum
73	guez	189	qp(il, i) = q(il, i)
74	guez	186	up(il, i) = u(il, i)
75			vp(il, i) = v(il, i)
76			wt(il, i) = 0.001
77			water(il, i) = 0.
78			evap(il, i) = 0.
79			lvcp(il, i) = lv(il, i) / cpn(il, i)
80	guez	97	enddo
81			enddo
82	guez	47
83	guez	193	! Check whether ep(inb) = 0. If so, skip precipitating downdraft
84			! calculation.
85
86	guez	187	forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4
87	guez	47
88	guez	193	imax = nl - 1
89			do while (.not. any(inb >= imax .and. lwork) .and. imax >= 1)
90			imax = imax - 1
91			end do
92	guez	47
93	guez	193	downdraft_loop: DO i = imax, 1, - 1
94			! Integrate liquid water equation to find condensed water
95			! and condensed water flux
96	guez	47
97	guez	193	! Calculate detrained precipitation
98	guez	194	forall (il = 1:ncum, inb(il) >= i .and. lwork(il)) wdtrain(il) = grav &
99			* (ep(il, i) * m(il, i) * clw(il, i) &
100			+ sum(max(elij(il, :i - 1, i) - (1. - ep(il, i)) * clw(il, i), 0.) &
101			* ment(il, :i - 1, i)))
102	guez	193
103			! Find rain water and evaporation using provisional
104			! estimates of qp(i) and qp(i - 1)
105	guez	47
106	guez	195	loop_horizontal: do il = 1, ncum
107	guez	193	if (i <= inb(il) .and. lwork(il)) then
108			wt(il, i) = 45.
109
110			if (i < inb(il)) then
111			qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) - t(il, i)) &
112			+ gz(il, i + 1) - gz(il, i)) / lv(il, i)
113			qp(il, i) = 0.5 * (qp(il, i) + q(il, i))
114	guez	97	endif
115	guez	47
116	guez	193	qp(il, i) = max(qp(il, i), 0.)
117			qp(il, i) = min(qp(il, i), qs(il, i))
118			qp(il, inb(il)) = q(il, inb(il))
119	guez	47
120	guez	193	if (i == 1) then
121			afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &
122			/ (1e4 + 2000. * p(il, 1) * qs(il, 1))
123			else
124			qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) - t(il, i - 1)) &
125			+ gz(il, i) - gz(il, i - 1)) / lv(il, i)
126			qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))
127			qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))
128			qp(il, i - 1) = max(qp(il, i - 1), 0.)
129			afac1 = p(il, i) * (qs(il, i) - qp(il, i)) &
130			/ (1e4 + 2000. * p(il, i) * qs(il, i))
131			afac2 = p(il, i - 1) * (qs(il, i - 1) - qp(il, i - 1)) &
132			/ (1e4 + 2000. * p(il, i - 1) * qs(il, i - 1))
133			afac = 0.5 * (afac1 + afac2)
134			endif
135	guez	47
136	guez	193	if (i == inb(il)) afac = 0.
137			afac = max(afac, 0.)
138			bfac = 1. / (sigd * wt(il, i))
139	guez	47
140	guez	193	! Prise en compte de la variation progressive de sigt dans
141			! les couches icb et icb - 1:
142	guez	195	! pour plcl <= ph(i + 1), pr1 = 0 et pr2 = 1
143			! pour plcl >= ph(i), pr1 = 1 et pr2 = 0
144	guez	193	! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion \`a cheval
145			! sur le nuage, et pr2 est la proportion sous la base du
146			! nuage.
147	guez	195	pr1 = max(0., min(1., &
148			(plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))))
149			pr2 = max(0., min(1., &
150			(ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))))
151			sigt = sigp * pr1 + pr2
152	guez	188
153	guez	193	b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt * afac
154			c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
155			* (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)
156	guez	195
157	guez	193	if (c6 > 0.) then
158			revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))
159			evap(il, i) = sigt * afac * revap
160			water(il, i) = revap * revap
161			else
162			evap(il, i) = - evap(il, i + 1) + 0.02 * (wdtrain(il) + sigd &
163			* wt(il, i) * water(il, i + 1)) / (sigd * (ph(il, i) &
164			- ph(il, i + 1)))
165			end if
166	guez	47
167	guez	193	! Calculate precipitating downdraft mass flux under
168			! hydrostatic approximation
169	guez	188
170	guez	193	test_above_surface: if (i /= 1) then
171			tevap = max(0., evap(il, i))
172			delth = max(0.001, (th(il, i) - th(il, i - 1)))
173			mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
174			* (p(il, i - 1) - p(il, i)) / delth
175	guez	47
176	guez	193	! If hydrostatic assumption fails, solve cubic
177			! difference equation for downdraft theta and mass
178			! flux from two simultaneous differential equations
179	guez	47
180	guez	193	amfac = sigd * sigd * 70. * ph(il, i) &
181			* (p(il, i - 1) - p(il, i)) &
182			* (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
183			amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
184			* mp(il, i))
185	guez	47
186	guez	193	if (amp2 > 0.1 * amfac) then
187			xf = 100. * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))
188			tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
189			* t(il, i) / (lvcp(il, i) * sigd * th(il, i))
190			af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
191			bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
192			* mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
193			- p(il, i)) * xf * tevap
194			fac2 = 1.
195			if (bf < 0.) fac2 = - 1.
196			bf = abs(bf)
197			ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
198	guez	186
199	guez	193	if (ur >= 0.) then
200			sru = sqrt(ur)
201			fac = 1.
202			if ((0.5 * bf - sru) < 0.) fac = - 1.
203			mp(il, i) = mp(il, i + 1) * tinv + (0.5 * bf &
204			+ sru)*tinv + fac (abs(0.5 * bf - sru))**tinv
205			else
206			d = atan(2. * sqrt(- ur) / (bf + 1e-28))
207			if (fac2 < 0.)d = 3.14159 - d
208			mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) &
209			* cos(d * tinv)
210			endif
211	guez	47
212	guez	193	mp(il, i) = max(0., mp(il, i))
213	guez	47
214	guez	193	! Il y a vraisemblablement une erreur dans la
215			! ligne suivante : il faut diviser par (mp(il,
216			! i) * sigd * grav) et non par (mp(il, i) + sigd
217			! * 0.1). Et il faut bien revoir les facteurs
218			! 100.
219			b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &
220			* tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &
221			- th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &
222			* th(il, i))
223			b(il, i - 1) = max(b(il, i - 1), 0.)
224			endif
225	guez	188
226	guez	193	! Limit magnitude of mp(i) to meet CFL condition:
227			ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
228			amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
229			ampmax = min(ampmax, amp2)
230			mp(il, i) = min(mp(il, i), ampmax)
231	guez	188
232	guez	193	! Force mp to decrease linearly to zero between cloud
233			! base and the surface:
234			if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
235			* (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))
236			endif test_above_surface
237	guez	47
238	guez	193	! Find mixing ratio of precipitating downdraft
239	guez	188
240	guez	193	if (i /= inb(il)) then
241			qp(il, i) = q(il, i)
242	guez	47
243	guez	193	if (mp(il, i) > mp(il, i + 1)) then
244			qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &
245			* (mp(il, i) - mp(il, i + 1)) + 100. * ginv * 0.5 &
246			* sigd * (ph(il, i) - ph(il, i + 1)) &
247			* (evap(il, i + 1) + evap(il, i))
248			qp(il, i) = qp(il, i) / mp(il, i)
249			up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
250			* (mp(il, i) - mp(il, i + 1))
251			up(il, i) = up(il, i) / mp(il, i)
252			vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &
253			* (mp(il, i) - mp(il, i + 1))
254			vp(il, i) = vp(il, i) / mp(il, i)
255			else
256			if (mp(il, i + 1) > 1e-16) then
257			qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &
258			* (ph(il, i) - ph(il, i + 1)) * (evap(il, i + 1) &
259			+ evap(il, i)) / mp(il, i + 1)
260			up(il, i) = up(il, i + 1)
261			vp(il, i) = vp(il, i + 1)
262	guez	97	endif
263	guez	193	endif
264	guez	188
265	guez	193	qp(il, i) = min(qp(il, i), qs(il, i))
266			qp(il, i) = max(qp(il, i), 0.)
267	guez	97	endif
268	guez	193	endif
269	guez	195	end do loop_horizontal
270	guez	186	end DO downdraft_loop
271	guez	47
272	guez	185	end SUBROUTINE cv30_unsat
273	guez	97
274	guez	185	end module cv30_unsat_m