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