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
    use SUPHEC_M, only: rg

    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(:, :) ! (ncum, nl)
    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(:, :) ! (ncum, nl)
    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)
    ! mass flux of the unsaturated downdraft, defined positive downward
    ! M_p in Emanuel (1991 928)

    real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl)
    real, intent(out):: wt(:, :) ! (ncum, nl)

    real, intent(out):: water(:, :) ! (ncum, nl)
    ! precipitation mixing ratio, l_p in Emanuel (1991 928)

    real, intent(out):: evap(:, :) ! (ncum, nl)
    ! sigt * rate of evaporation of precipitation, in s-1
    ! \sigma_s E in Emanuel (1991 928)

    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) = rg &
            * (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**2 * 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)**2 - mp(il, i)**2)

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