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)

    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(:, :) ! (klon, 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(:) ! (klon)

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