phylmd/CV30_routines/cv30_unsat.f

module cv30_unsat_m

  implicit none

contains

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

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

    ! inputs:
    integer, intent(in):: ncum
    integer, intent(in):: icb(:), inb(:) ! (ncum)
    real t(klon, klev), rr(klon, klev), rs(klon, klev)
    real gz(klon, klev)
    real u(klon, klev), v(klon, klev)
    real p(klon, klev), ph(klon, klev + 1)
    real th(klon, klev)
    real tv(klon, klev)
    real lv(klon, klev)
    real cpn(klon, klev)
    real, intent(in):: ep(klon, klev), sigp(klon, klev), clw(klon, klev)
    real m(klon, klev), ment(klon, klev, klev), elij(klon, klev, klev)
    real, intent(in):: delt
    real plcl(klon)

    ! outputs:
    real mp(klon, klev), rp(klon, klev), up(klon, klev), vp(klon, klev)
    real wt(klon, klev), water(klon, klev), evap(klon, klev)
    real b(:, :) ! (ncum, klev)

    ! Local:
    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(ncum)
    logical lwork(ncum)

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

    delti = 1. / delt
    tinv = 1. / 3.
    mp = 0.

    do i = 1, nl
       do il = 1, ncum
          mp(il, i) = 0.
          rp(il, i) = rr(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.
          b(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 = amax1(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 rp(i)and rp(i - 1) 

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

                if (i < inb(il)) then
                   rp(il, i) = rp(il, i + 1) + (cpd * (t(il, i + 1) &
                        - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)
                   rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))
                endif
                rp(il, i) = amax1(rp(il, i), 0.)
                rp(il, i) = amin1(rp(il, i), rs(il, i))
                rp(il, inb(il)) = rr(il, inb(il))

                if (i == 1) then
                   afac = p(il, 1) * (rs(il, 1) - rp(il, 1)) &
                        / (1e4 + 2000. * p(il, 1) * rs(il, 1))
                else
                   rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) &
                        - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)
                   rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))
                   rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1))
                   rp(il, i - 1) = amax1(rp(il, i - 1), 0.)
                   afac1 = p(il, i) * (rs(il, i) - rp(il, i)) &
                        / (1e4 + 2000. * p(il, i) * rs(il, i))
                   afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) &
                        / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1))
                   afac = 0.5 * (afac1 + afac2)
                endif
                if (i == inb(il))afac = 0.
                afac = amax1(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 = amax1(0., evap(il, i))
                   delth = amax1(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 eqns 

                   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) = amax1(0., mp(il, i))

                      ! Il y a vraisemblablement une erreur dans la
                      ! ligne 2 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) = amax1(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 = amin1(ampmax, amp2)
                   mp(il, i) = amin1(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))) then
                      mp(il, i) = mp(il, icb(il)) * (p(il, 1) - p(il, i)) &
                           / (p(il, 1) - p(il, icb(il)))
                   endif
                endif ! i == 1

                ! find mixing ratio of precipitating downdraft 

                if (i /= inb(il)) then
                   rp(il, i) = rr(il, i)

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