/[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 192 by guez, Thu May 12 13:00:07 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, num1
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))
52          tinv=1./3.      logical lwork(size(icb))
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)
71  !   ***  check whether ep(inb)=0, if so, skip precipitating    ***         enddo
72  !   ***             downdraft calculation                      ***      enddo
73  !  
74        ! check whether ep(inb) = 0, if so, skip precipitating
75          do il=1,ncum      ! downdraft calculation
76            lwork(il)=.TRUE.      forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4
77            if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.  
78          enddo      wdtrain = 0.
79    
80          call zilch(wdtrain,ncum)      downdraft_loop: DO i = nl - 1, 1, - 1
81           num1 = 0
82          DO 400 i=nl+1,1,-1  
83           do il = 1, ncum
84          num1=0            if (i <= inb(il) .and. lwork(il)) num1 = num1 + 1
         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  
85         enddo         enddo
86    
87         if(i.gt.1)then         if (num1 > 0) then
88          do 320 j=1,i-1            ! integrate liquid water equation to find condensed water
89           do il=1,ncum            ! and condensed water flux
90            if (i.le.inb(il) .and. lwork(il)) then  
91             awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)            ! calculate detrained precipitation
92             awat=amax1(awat,0.0)  
93             if (cvflag_grav) then            do il = 1, ncum
94              wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i)               if (i <= inb(il) .and. lwork(il)) then
95             else                  wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
96              wdtrain(il)=wdtrain(il)+10.0*awat*ment(il,j,i)               endif
97             endif            enddo
98    
99              if (i > 1) then
100                 do j = 1, i - 1
101                    do il = 1, ncum
102                       if (i <= inb(il) .and. lwork(il)) then
103                          awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
104                          awat = max(awat, 0.)
105                          wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i)
106                       endif
107                    enddo
108                 end do
109            endif            endif
          enddo  
 320     continue  
        endif  
   
 !  
 !    ***    find rain water and evaporation using provisional   ***  
 !    ***              estimates of rp(i)and rp(i-1)             ***  
 !  
   
       do 999 il=1,ncum  
   
        if (i.le.inb(il) .and. lwork(il)) then  
   
       wt(il,i)=45.0  
   
       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)  
        endif  
        mp(il,i)=amax1(0.0,mp(il,i))  
   
        if (cvflag_grav) then  
 !jyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:  
 ! il faut diviser par (mp(il,i)*sigd*grav) et non par (mp(il,i)+sigd*0.1).  
 ! Et il faut bien revoir les facteurs 100.  
         b(il,i-1)=b(il,i)+100.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)  
   
        endif  
       endif  
       rp(il,i)=amin1(rp(il,i),rs(il,i))  
       rp(il,i)=amax1(rp(il,i),0.0)  
   
       endif  
       endif  
 999   continue  
110    
111  400   continue            ! find rain water and evaporation using provisional
112              ! estimates of qp(i) and qp(i - 1)
113    
114              do il = 1, ncum
115                 if (i <= inb(il) .and. lwork(il)) then
116                    wt(il, i) = 45.
117    
118                    if (i < inb(il)) then
119                       qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) &
120                            - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)
121                       qp(il, i) = 0.5 * (qp(il, i) + q(il, i))
122                    endif
123    
124                    qp(il, i) = max(qp(il, i), 0.)
125                    qp(il, i) = min(qp(il, i), qs(il, i))
126                    qp(il, inb(il)) = q(il, inb(il))
127    
128                    if (i == 1) then
129                       afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) &
130                            / (1e4 + 2000. * p(il, 1) * qs(il, 1))
131                    else
132                       qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) &
133                            - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)
134                       qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1))
135                       qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1))
136                       qp(il, i - 1) = max(qp(il, i - 1), 0.)
137                       afac1 = p(il, i) * (qs(il, i) - qp(il, i)) &
138                            / (1e4 + 2000. * p(il, i) * qs(il, i))
139                       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                       afac = 0.5 * (afac1 + afac2)
142                    endif
143    
144                    if (i == inb(il)) afac = 0.
145                    afac = max(afac, 0.)
146                    bfac = 1. / (sigd * wt(il, i))
147    
148                    ! prise en compte de la variation progressive de sigt dans
149                    ! les couches icb et icb - 1:
150                    ! pour plcl < ph(i + 1), pr1 = 0 & pr2 = 1
151                    ! 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                          mp(il, i) = max(0., mp(il, i))
223    
224                          ! Il y a vraisemblablement une erreur dans la
225                          ! ligne suivante : il faut diviser par (mp(il,
226                          ! i) * sigd * grav) et non par (mp(il, i) + sigd
227                          ! * 0.1).  Et il faut bien revoir les facteurs
228                          ! 100.
229                          b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
230                               - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
231                               - 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
235    
236                       ! limit magnitude of mp(i) to meet cfl condition
237                       ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
238                       amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
239                       ampmax = min(ampmax, amp2)
240                       mp(il, i) = min(mp(il, i), ampmax)
241    
242                       ! force mp to decrease linearly to zero
243                       ! between cloud base and the surface
244                       if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) &
245                            * (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il)))
246                    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
273                       endif
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
279              end do
280           end if
281        end DO downdraft_loop
282    
283      end SUBROUTINE cv30_unsat
284    
285         return  end module cv30_unsat_m
        end  
Legend:
Removed from v.69  
changed lines
  Added in v.192
  ViewVC Help
Powered by ViewVC 1.1.21