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}

    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, 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))

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

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