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