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, sigp, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, &
       evap, b)

    use cv30_param_m, only: nl, sigd
    use cv_thermo_m, only: cpd, ginv, grav
    USE dimphy, ONLY: klon, klev

    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(:, :) ! (klon, klev)
    real, intent(in):: gz(:, :) ! (klon, klev)
    real, intent(in):: u(:, :), v(:, :) ! (klon, klev)
    real, intent(in):: p(klon, klev), ph(klon, klev + 1)
    real, intent(in):: th(klon, klev)
    real, intent(in):: tv(klon, klev)
    real, intent(in):: lv(klon, klev)
    real, intent(in):: cpn(klon, klev)
    real, intent(in):: ep(:, :), sigp(:, :), 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)

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

    ! Local:
    integer ncum
    integer i, j, il, imax
    real tinv, delti
    real awat, 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(klon, klev)
    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 

       do il = 1, ncum
          if (i <= inb(il) .and. lwork(il)) then
             wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
          endif
       enddo

       if (i > 1) then
          do j = 1, i - 1
             do il = 1, ncum
                if (i <= inb(il) .and. lwork(il)) then
                   awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
                   awat = max(awat, 0.)
                   wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
                endif
             enddo
          end do
       endif

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

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