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, num1
    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))
    logical lwork(size(icb))

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

    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

    wdtrain = 0.

    downdraft_loop: DO i = nl - 1, 1, - 1
       num1 = 0

       do il = 1, ncum
          if (i <= inb(il) .and. lwork(il)) num1 = num1 + 1
       enddo

       if (num1 > 0) then
          ! 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 & pr2 = 1
                ! pour plcl > ph(i), pr1 = 1 & 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 

                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 ! i == 1

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