/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_unsat.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV30_routines/cv30_unsat.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.186  
changed lines
  Added in v.201

  ViewVC Help
Powered by ViewVC 1.1.21