/[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

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