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

    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(:, :) ! (ncum, nl) temperature (K)
    real, intent(in):: 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) potential temperature, in K
    real, intent(in):: tv(:, :) ! (klon, klev)

    real, intent(in):: lv(:, :) ! (ncum, nl)
    ! specific latent heat of vaporization of water, in J kg-1

    real, intent(in):: cpn(:, :) ! (ncum, nl)
    ! specific heat capacity at constant pressure of humid air, in J K-1 kg-1

    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(:, :)
    ! (ncum, nl) Mass flux of the unsaturated downdraft, defined
    ! positive downward, in kg m-2 s-1.  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, parameter:: tinv = 1. / 3.
    real  delti
    real afac, afac1, afac2, bfac
    real pr1, sigt, b6, c6, revap, tevap, delth
    real xf, tf, fac2, ur, sru, fac, d, af, bf
    real ampmax
    real lvcp(size(icb), nl) ! (ncum, nl) L_v / C_p, in K
    real wdtrain(size(icb)) ! (ncum)
    logical lwork(size(icb)) ! (ncum)

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

    ncum = size(icb)
    delti = 1. / delt
    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) + (rcpd * (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) + (rcpd * (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
                if (abs(mp(il, i + 1)**2 - mp(il, i)**2) > 0.1 * 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))) &
                     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) = max(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)), 0.)
                endif

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