/[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 185 by guez, Wed Mar 16 15:04:46 2016 UTC revision 186 by guez, Mon Mar 21 15:36:26 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 &    SUBROUTINE cv30_unsat(nloc, ncum, nd, na, icb, inb, t, rr, rs, gz, u, v, p, &
8         ,t,rr,rs,gz,u,v,p,ph &         ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, &
9         ,th,tv,lv,cpn,ep,sigp,clw &         up, vp, wt, water, evap, b)
10         ,m,ment,elij,delt,plcl &  
11         ,mp,rp,up,vp,wt,water,evap,b)      use cv30_param_m, only: nl, sigd
12      use cv30_param_m      use cvflag, only: cvflag_grav
13      use cvthermo      use cvthermo, only: cpd, ginv, grav
     use cvflag  
   
14    
15      ! inputs:      ! inputs:
16      integer, intent(in):: ncum, nd, na, nloc      integer, intent(in):: nloc, ncum, nd, na
17      integer icb(nloc), inb(nloc)      integer, intent(in):: icb(:), inb(:) ! (ncum)
18        real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
19        real gz(nloc, na)
20        real u(nloc, nd), v(nloc, nd)
21        real p(nloc, nd), ph(nloc, nd + 1)
22        real th(nloc, na)
23        real tv(nloc, na)
24        real lv(nloc, na)
25        real cpn(nloc, na)
26        real ep(nloc, na), sigp(nloc, na), clw(nloc, na)
27        real m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
28      real, intent(in):: delt      real, intent(in):: delt
29      real plcl(nloc)      real plcl(nloc)
     real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)  
     real u(nloc,nd), v(nloc,nd)  
     real p(nloc,nd), ph(nloc,nd+1)  
     real th(nloc,na), gz(nloc,na)  
     real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)  
     real cpn(nloc,na), tv(nloc,na)  
     real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)  
30    
31      ! outputs:      ! outputs:
32      real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)      real mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
33      real water(nloc,na), evap(nloc,na), wt(nloc,na)      real wt(nloc, na), water(nloc, na), evap(nloc, na)
34      real b(nloc,na)      real b(:, :) ! (nloc, na)
35    
36      ! local variables      ! Local:
37      integer i,j,il,num1      integer i, j, il, num1
38      real tinv, delti      real tinv, delti
39      real awat, afac, afac1, afac2, bfac      real awat, afac, afac1, afac2, bfac
40      real pr1, pr2, sigt, b6, c6, revap, tevap, delth      real pr1, pr2, sigt, b6, c6, revap, tevap, delth
41      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
42      real ampmax      real ampmax
43      real lvcp(nloc,na)      real lvcp(nloc, na)
44      real wdtrain(nloc)      real wdtrain(nloc)
45      logical lwork(nloc)      logical lwork(nloc)
46    
   
47      !------------------------------------------------------      !------------------------------------------------------
48    
49      delti = 1./delt      delti = 1. / delt
50      tinv=1./3.      tinv = 1. / 3.
51        mp = 0.
52      mp(:,:)=0.  
53        do i = 1, nl
54      do i=1,nl         do il = 1, ncum
55         do il=1,ncum            mp(il, i) = 0.
56            mp(il,i)=0.0            rp(il, i) = rr(il, i)
57            rp(il,i)=rr(il,i)            up(il, i) = u(il, i)
58            up(il,i)=u(il,i)            vp(il, i) = v(il, i)
59            vp(il,i)=v(il,i)            wt(il, i) = 0.001
60            wt(il,i)=0.001            water(il, i) = 0.
61            water(il,i)=0.0            evap(il, i) = 0.
62            evap(il,i)=0.0            b(il, i) = 0.
63            b(il,i)=0.0            lvcp(il, i) = lv(il, i) / cpn(il, i)
           lvcp(il,i)=lv(il,i)/cpn(il,i)  
64         enddo         enddo
65      enddo      enddo
66    
67      !      ! check whether ep(inb) = 0, if so, skip precipitating
68      !   ***  check whether ep(inb)=0, if so, skip precipitating    ***      ! downdraft calculation
69      !   ***             downdraft calculation                      ***  
70      !      do il = 1, ncum
71           lwork(il) = .TRUE.
72      do il=1,ncum         if (ep(il, inb(il)) < 0.0001) lwork(il) = .FALSE.
        lwork(il)=.TRUE.  
        if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.  
73      enddo      enddo
74    
75      call zilch(wdtrain,ncum)      wdtrain(:ncum) = 0.
76    
77      DO  i=nl+1,1,-1      downdraft_loop: DO i = nl + 1, 1, - 1
78           num1 = 0
79    
80         num1=0         do il = 1, ncum
81         do il=1,ncum            if (i <= inb(il) .and. lwork(il)) num1 = num1 + 1
           if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1  
82         enddo         enddo
        if (num1.le.0) cycle  
   
        !  
        !   ***  integrate liquid water equation to find condensed water   ***  
        !   ***                and condensed water flux                    ***  
        !  
   
        !  
        !    ***                    begin downdraft loop                    ***  
        !  
   
        !  
        !    ***              calculate detrained precipitation             ***  
        !  
        do il=1,ncum  
           if (i.le.inb(il) .and. lwork(il)) then  
              if (cvflag_grav) then  
                 wdtrain(il)=grav*ep(il,i)*m(il,i)*clw(il,i)  
              else  
                 wdtrain(il)=10.0*ep(il,i)*m(il,i)*clw(il,i)  
              endif  
           endif  
        enddo  
   
        if(i.gt.1)then  
           do j=1,i-1  
              do il=1,ncum  
                 if (i.le.inb(il) .and. lwork(il)) then  
                    awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)  
                    awat=amax1(awat,0.0)  
                    if (cvflag_grav) then  
                       wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i)  
                    else  
                       wdtrain(il)=wdtrain(il)+10.0*awat*ment(il,j,i)  
                    endif  
                 endif  
              enddo  
           end do  
        endif  
   
        !  
        !    ***    find rain water and evaporation using provisional   ***  
        !    ***              estimates of rp(i)and rp(i-1)             ***  
        !  
   
        do  il=1,ncum  
83    
84            if (i.le.inb(il) .and. lwork(il)) then         if (num1 > 0) then
85              ! integrate liquid water equation to find condensed water
86              ! and condensed water flux
87    
88               wt(il,i)=45.0            ! calculate detrained precipitation
   
              if(i.lt.inb(il))then  
                 rp(il,i)=rp(il,i+1) &  
                      +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)  
                 rp(il,i)=0.5*(rp(il,i)+rr(il,i))  
              endif  
              rp(il,i)=amax1(rp(il,i),0.0)  
              rp(il,i)=amin1(rp(il,i),rs(il,i))  
              rp(il,inb(il))=rr(il,inb(il))  
   
              if(i.eq.1)then  
                 afac=p(il,1)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1))  
              else  
                 rp(il,i-1)=rp(il,i) &  
                      +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)  
                 rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))  
                 rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))  
                 rp(il,i-1)=amax1(rp(il,i-1),0.0)  
                 afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i))  
                 afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1)) &  
                      /(1.0e4+2000.0*p(il,i-1)*rs(il,i-1))  
                 afac=0.5*(afac1+afac2)  
              endif  
              if(i.eq.inb(il))afac=0.0  
              afac=amax1(afac,0.0)  
              bfac=1./(sigd*wt(il,i))  
              !  
              !jyg1  
              !cc        sigt=1.0  
              !cc        if(i.ge.icb)sigt=sigp(i)  
              ! prise en compte de la variation progressive de sigt dans  
              ! les couches icb et icb-1:  
              !     pour plcl<ph(i+1), pr1=0 & pr2=1  
              !     pour plcl>ph(i),   pr1=1 & pr2=0  
              !     pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval  
              !    sur le nuage, et pr2 est la proportion sous la base du  
              !    nuage.  
              pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))  
              pr1=max(0.,min(1.,pr1))  
              pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))  
              pr2=max(0.,min(1.,pr2))  
              sigt=sigp(il,i)*pr1+pr2  
              !jyg2  
              !  
              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.gt.0.0)then  
                 revap=0.5*(-b6+sqrt(b6*b6+4.*c6))  
                 evap(il,i)=sigt*afac*revap  
                 water(il,i)=revap*revap  
              else  
                 evap(il,i)=-evap(il,i+1) &  
                      +0.02*(wdtrain(il)+sigd*wt(il,i)*water(il,i+1)) &  
                      /(sigd*(ph(il,i)-ph(il,i+1)))  
              end if  
              !  
              !    ***  calculate precipitating downdraft mass flux under     ***  
              !    ***              hydrostatic approximation                 ***  
              !  
              if (i.ne.1) then  
89    
90                  tevap=amax1(0.0,evap(il,i))            do il = 1, ncum
91                  delth=amax1(0.001,(th(il,i)-th(il,i-1)))               if (i <= inb(il) .and. lwork(il)) then
92                  if (cvflag_grav) then                  if (cvflag_grav) then
93                     mp(il,i)=100.*ginv*lvcp(il,i)*sigd*tevap &                     wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
                         *(p(il,i-1)-p(il,i))/delth  
94                  else                  else
95                     mp(il,i)=10.*lvcp(il,i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth                     wdtrain(il) = 10. * ep(il, i) * m(il, i) * clw(il, i)
96                  endif                  endif
97                  !               endif
98                  !    ***           if hydrostatic assumption fails,             ***            enddo
99                  !    ***   solve cubic difference equation for downdraft theta  ***  
100                  !    ***  and mass flux from two simultaneous differential eqns ***            if (i > 1) then
101                  !               do j = 1, i - 1
102                  amfac=sigd*sigd*70.0*ph(il,i)*(p(il,i-1)-p(il,i)) &                  do il = 1, ncum
103                       *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))                     if (i <= inb(il) .and. lwork(il)) then
104                  amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))                        awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
105                  if(amp2.gt.(0.1*amfac))then                        awat = amax1(awat, 0.)
106                     xf=100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1))                        if (cvflag_grav) then
107                     tf=b(il,i)-5.0*(th(il,i)-th(il,i-1))*t(il,i) &                           wdtrain(il) = wdtrain(il) + grav * awat &
108                          /(lvcp(il,i)*sigd*th(il,i))                                * ment(il, j, i)
109                     af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv                        else
110                     bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf &                           wdtrain(il) = wdtrain(il) + 10. * awat * ment(il, j, i)
111                          +50.*(p(il,i-1)-p(il,i))*xf*tevap                        endif
                    fac2=1.0  
                    if(bf.lt.0.0)fac2=-1.0  
                    bf=abs(bf)  
                    ur=0.25*bf*bf-af*af*af*tinv*tinv*tinv  
                    if(ur.ge.0.0)then  
                       sru=sqrt(ur)  
                       fac=1.0  
                       if((0.5*bf-sru).lt.0.0)fac=-1.0  
                       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+1.0e-28))  
                       if(fac2.lt.0.0)d=3.14159-d  
                       mp(il,i)=mp(il,i+1)*tinv+2.*sqrt(af*tinv)*cos(d*tinv)  
112                     endif                     endif
113                     mp(il,i)=amax1(0.0,mp(il,i))                  enddo
114                 end do
115              endif
116    
117              ! find rain water and evaporation using provisional
118              ! estimates of rp(i)and rp(i - 1)
119    
120              do il = 1, ncum
121                 if (i <= inb(il) .and. lwork(il)) then
122                    wt(il, i) = 45.
123    
124                    if (i < inb(il)) then
125                       rp(il, i) = rp(il, i + 1) + (cpd * (t(il, i + 1) &
126                            - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)
127                       rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))
128                    endif
129                    rp(il, i) = amax1(rp(il, i), 0.)
130                    rp(il, i) = amin1(rp(il, i), rs(il, i))
131                    rp(il, inb(il)) = rr(il, inb(il))
132    
133                    if (i == 1) then
134                       afac = p(il, 1) * (rs(il, 1) - rp(il, 1)) &
135                            / (1e4 + 2000. * p(il, 1) * rs(il, 1))
136                    else
137                       rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) &
138                            - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)
139                       rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))
140                       rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1))
141                       rp(il, i - 1) = amax1(rp(il, i - 1), 0.)
142                       afac1 = p(il, i) * (rs(il, i) - rp(il, i)) &
143                            / (1e4 + 2000. * p(il, i) * rs(il, i))
144                       afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) &
145                            / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1))
146                       afac = 0.5 * (afac1 + afac2)
147                    endif
148                    if (i == inb(il))afac = 0.
149                    afac = amax1(afac, 0.)
150                    bfac = 1. / (sigd * wt(il, i))
151    
152                    ! prise en compte de la variation progressive de sigt dans
153                    ! les couches icb et icb - 1:
154                    ! pour plcl < ph(i + 1), pr1 = 0 & pr2 = 1
155                    ! pour plcl > ph(i), pr1 = 1 & pr2 = 0
156                    ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion a cheval
157                    ! sur le nuage, et pr2 est la proportion sous la base du
158                    ! nuage.
159                    pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))
160                    pr1 = max(0., min(1., pr1))
161                    pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))
162                    pr2 = max(0., min(1., pr2))
163                    sigt = sigp(il, i) * pr1 + pr2
164    
165                    b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt &
166                         * afac
167                    c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
168                         * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)
169                    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) &
175                            + sigd * wt(il, i) * water(il, i + 1)) &
176                            / (sigd * (ph(il, i) - ph(il, i + 1)))
177                    end if
178    
179                    ! calculate precipitating downdraft mass flux under
180                    ! hydrostatic approximation
181    
182                    if (i /= 1) then
183                       tevap = amax1(0., evap(il, i))
184                       delth = amax1(0.001, (th(il, i) - th(il, i - 1)))
185                     if (cvflag_grav) then                     if (cvflag_grav) then
186                        !jyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:                        mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
187                        ! il faut diviser par (mp(il,i)*sigd*grav) et non par (mp(il,i)+sigd*0.1).                             * (p(il, i - 1) - p(il, i)) / delth
                       ! Et il faut bien revoir les facteurs 100.  
                       b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap &  
                            /(mp(il,i)+sigd*0.1) &  
                            -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i))  
188                     else                     else
189                        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap &                        mp(il, i) = 10. * lvcp(il, i) * sigd * tevap &
190                             /(mp(il,i)+sigd*0.1) &                             * (p(il, i - 1) - p(il, i)) / delth
                            -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i))  
191                     endif                     endif
                    b(il,i-1)=amax1(b(il,i-1),0.0)  
                 endif  
                 !  
                 !   ***         limit magnitude of mp(i) to meet cfl condition      ***  
                 !  
                 ampmax=2.0*(ph(il,i)-ph(il,i+1))*delti  
                 amp2=2.0*(ph(il,i-1)-ph(il,i))*delti  
                 ampmax=amin1(ampmax,amp2)  
                 mp(il,i)=amin1(mp(il,i),ampmax)  
                 !  
                 !    ***      force mp to decrease linearly to zero                 ***  
                 !    ***       between cloud base and the surface                   ***  
                 !  
                 if(p(il,i).gt.p(il,icb(il)))then  
                    mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))  
                 endif  
   
              endif ! i.eq.1  
              !  
              !    ***       find mixing ratio of precipitating downdraft     ***  
              !  
192    
193               if (i.ne.inb(il)) then                     ! if hydrostatic assumption fails,
194                       ! solve cubic difference equation for downdraft theta
195                       ! and mass flux from two simultaneous differential eqns
196    
197                       amfac = sigd * sigd * 70. * ph(il, i) &
198                            * (p(il, i - 1) - p(il, i)) &
199                            * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
200                       amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
201                            * mp(il, i))
202                       if (amp2 > (0.1 * amfac)) then
203                          xf = 100. * sigd * sigd * sigd * (ph(il, i) &
204                               - ph(il, i + 1))
205                          tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
206                               * t(il, i) / (lvcp(il, i) * sigd * th(il, i))
207                          af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
208                          bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
209                               * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
210                               - p(il, i)) * xf * tevap
211                          fac2 = 1.
212                          if (bf < 0.)fac2 = - 1.
213                          bf = abs(bf)
214                          ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
215                          if (ur >= 0.) then
216                             sru = sqrt(ur)
217                             fac = 1.
218                             if ((0.5 * bf - sru) < 0.)fac = - 1.
219                             mp(il, i) = mp(il, i + 1) * tinv &
220                                  + (0.5 * bf + sru)**tinv &
221                                  + fac * (abs(0.5 * bf - sru))**tinv
222                          else
223                             d = atan(2. * sqrt(- ur) / (bf + 1e-28))
224                             if (fac2 < 0.)d = 3.14159 - d
225                             mp(il, i) = mp(il, i + 1) * tinv + 2. &
226                                  * sqrt(af * tinv) * cos(d * tinv)
227                          endif
228                          mp(il, i) = amax1(0., mp(il, i))
229    
230                  rp(il,i)=rr(il,i)                        if (cvflag_grav) then
231                             ! Il y a vraisemblablement une erreur dans la
232                             ! ligne 2 suivante: il faut diviser par (mp(il,
233                             ! i) * sigd * grav) et non par (mp(il, i) + sigd
234                             ! * 0.1).  Et il faut bien revoir les facteurs
235                             ! 100.
236                             b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
237                                  - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
238                                  - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
239                                  / (lvcp(il, i) * sigd * th(il, i))
240                          else
241                             b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
242                                  - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
243                                  - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
244                                  / (lvcp(il, i) * sigd * th(il, i))
245                          endif
246                          b(il, i - 1) = amax1(b(il, i - 1), 0.)
247                       endif
248    
249                  if(mp(il,i).gt.mp(il,i+1))then                     ! limit magnitude of mp(i) to meet cfl condition
250    
251                     if (cvflag_grav) then                     ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
252                        rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1)) &                     amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
253                             +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1)) &                     ampmax = amin1(ampmax, amp2)
254                             *(evap(il,i+1)+evap(il,i))                     mp(il, i) = amin1(mp(il, i), ampmax)
255                     else  
256                        rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1)) &                     ! force mp to decrease linearly to zero
257                             +5.*sigd*(ph(il,i)-ph(il,i+1)) &                     ! between cloud base and the surface
258                             *(evap(il,i+1)+evap(il,i))  
259                       if (p(il, i) > p(il, icb(il))) then
260                          mp(il, i) = mp(il, icb(il)) * (p(il, 1) - p(il, i)) &
261                               / (p(il, 1) - p(il, icb(il)))
262                     endif                     endif
263                     rp(il,i)=rp(il,i)/mp(il,i)                  endif ! i == 1
                    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)  
264    
265                  else                  ! find mixing ratio of precipitating downdraft
266    
267                    if (i /= inb(il)) then
268                       rp(il, i) = rr(il, i)
269    
270                     if(mp(il,i+1).gt.1.0e-16)then                     if (mp(il, i) > mp(il, i + 1)) then
271                        if (cvflag_grav) then                        if (cvflag_grav) then
272                           rp(il,i)=rp(il,i+1) &                           rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &
273                                +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1)) &                                * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &
274                                *(evap(il,i+1)+evap(il,i))/mp(il,i+1)                                * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &
275                                  * (evap(il, i + 1) + evap(il, i))
276                        else                        else
277                           rp(il,i)=rp(il,i+1) &                           rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &
278                                +5.*sigd*(ph(il,i)-ph(il,i+1)) &                                * (mp(il, i) - mp(il, i + 1)) + 5. * sigd &
279                                *(evap(il,i+1)+evap(il,i))/mp(il,i+1)                                * (ph(il, i) - ph(il, i + 1)) &
280                                  * (evap(il, i + 1) + evap(il, i))
281                          endif
282                          rp(il, i) = rp(il, i) / mp(il, i)
283                          up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
284                               * (mp(il, i) - mp(il, i + 1))
285                          up(il, i) = up(il, i) / mp(il, i)
286                          vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &
287                               * (mp(il, i) - mp(il, i + 1))
288                          vp(il, i) = vp(il, i) / mp(il, i)
289                       else
290                          if (mp(il, i + 1) > 1e-16) then
291                             if (cvflag_grav) then
292                                rp(il, i) = rp(il, i + 1) &
293                                     + 100. * ginv * 0.5 * sigd * (ph(il, i) &
294                                     - ph(il, i + 1)) &
295                                     * (evap(il, i + 1) + evap(il, i)) &
296                                     / mp(il, i + 1)
297                             else
298                                rp(il, i) = rp(il, i + 1) &
299                                     + 5. * sigd * (ph(il, i) - ph(il, i + 1)) &
300                                     * (evap(il, i + 1) + evap(il, i)) &
301                                     / mp(il, i + 1)
302                             endif
303                             up(il, i) = up(il, i + 1)
304                             vp(il, i) = vp(il, i + 1)
305                        endif                        endif
                       up(il,i)=up(il,i+1)  
                       vp(il,i)=vp(il,i+1)  
   
306                     endif                     endif
307                       rp(il, i) = amin1(rp(il, i), rs(il, i))
308                       rp(il, i) = amax1(rp(il, i), 0.)
309                  endif                  endif
                 rp(il,i)=amin1(rp(il,i),rs(il,i))  
                 rp(il,i)=amax1(rp(il,i),0.0)  
   
310               endif               endif
311            endif            end do
312         end do         end if
313        end DO downdraft_loop
     end DO  
314    
315    end SUBROUTINE cv30_unsat    end SUBROUTINE cv30_unsat
316    
Legend:
Removed from v.185  
changed lines
  Added in v.186
  ViewVC Help
Powered by ViewVC 1.1.21