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 * 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
                   ! * 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	193	amfac = sigd * sigd * 70. * ph(il, i) &
193			* (p(il, i - 1) - p(il, i)) &
194			* (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
195			amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
196			* mp(il, i))
197	guez	47
198	guez	193	if (amp2 > 0.1 * amfac) then
199			xf = 100. * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1))
200			tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
201			* t(il, i) / (lvcp(il, i) * sigd * th(il, i))
202			af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
203			bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
204			* mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
205			- p(il, i)) * xf * tevap
206			fac2 = 1.
207			if (bf < 0.) fac2 = - 1.
208			bf = abs(bf)
209			ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
210	guez	186
211	guez	193	if (ur >= 0.) then
212			sru = sqrt(ur)
213			fac = 1.
214			if ((0.5 * bf - sru) < 0.) fac = - 1.
215			mp(il, i) = mp(il, i + 1) * tinv + (0.5 * bf &
216			+ sru)*tinv + fac (abs(0.5 * bf - sru))**tinv
217			else
218			d = atan(2. * sqrt(- ur) / (bf + 1e-28))
219			if (fac2 < 0.)d = 3.14159 - d
220			mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) &
221			* cos(d * tinv)
222			endif
223	guez	47
224	guez	193	mp(il, i) = max(0., mp(il, i))
225	guez	47
226	guez	198	! Il y a vraisemblablement une erreur dans la ligne
227			! suivante : il faut diviser par (mp(il, i) * sigd
228			! * rg) et non par (mp(il, i) + sigd * 0.1). Et il
229			! faut bien revoir les facteurs 100.
230	guez	193	b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) &
231			* tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) &
232			- th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd &
233			* th(il, i))
234			b(il, i - 1) = max(b(il, i - 1), 0.)
235			endif
236	guez	188
237	guez	193	! Limit magnitude of mp(i) to meet CFL condition:
238			ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
239			amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
240			ampmax = min(ampmax, amp2)
241			mp(il, i) = min(mp(il, i), ampmax)
242	guez	188
243	guez	193	! Force mp to decrease linearly to zero between cloud
244			! base and the surface:
245			if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
246			* (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))
247			endif test_above_surface
248	guez	47
249	guez	193	! Find mixing ratio of precipitating downdraft
250	guez	188
251	guez	193	if (i /= inb(il)) then
252			qp(il, i) = q(il, i)
253	guez	47
254	guez	193	if (mp(il, i) > mp(il, i + 1)) then
255			qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) &
256			* (mp(il, i) - mp(il, i + 1)) + 100. * ginv * 0.5 &
257			* sigd * (ph(il, i) - ph(il, i + 1)) &
258			* (evap(il, i + 1) + evap(il, i))
259			qp(il, i) = qp(il, i) / mp(il, i)
260			up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
261			* (mp(il, i) - mp(il, i + 1))
262			up(il, i) = up(il, i) / mp(il, i)
263			vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &
264			* (mp(il, i) - mp(il, i + 1))
265			vp(il, i) = vp(il, i) / mp(il, i)
266			else
267			if (mp(il, i + 1) > 1e-16) then
268			qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd &
269			* (ph(il, i) - ph(il, i + 1)) * (evap(il, i + 1) &
270			+ evap(il, i)) / mp(il, i + 1)
271			up(il, i) = up(il, i + 1)
272			vp(il, i) = vp(il, i + 1)
273	guez	97	endif
274	guez	193	endif
275	guez	188
276	guez	193	qp(il, i) = min(qp(il, i), qs(il, i))
277			qp(il, i) = max(qp(il, i), 0.)
278	guez	97	endif
279	guez	193	endif
280	guez	195	end do loop_horizontal
281	guez	186	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