/[lmdze]/trunk/phylmd/CV30_routines/cv30_yield.f
ViewVC logotype

Diff of /trunk/phylmd/CV30_routines/cv30_yield.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/phylmd/CV3_routines/cv3_yield.f revision 105 by guez, Thu Sep 4 10:40:24 2014 UTC trunk/Sources/phylmd/CV30_routines/cv30_yield.f revision 188 by guez, Tue Mar 22 16:31:39 2016 UTC
# Line 1  Line 1 
1  module cv3_yield_m  module cv30_yield_m
2    
3    implicit none    implicit none
4    
5  contains  contains
6    
7    SUBROUTINE cv3_yield(nloc,ncum,nd,na  &    SUBROUTINE cv30_yield(icb, inb, delt, t, rr, u, v, gz, p, ph, h, hp, lv, &
8         ,icb,inb,delt &         cpn, th, ep, clw, m, tp, mp, rp, up, vp, wt, water, evap, b, ment, &
9         ,t,rr,u,v,gz,p,ph,h,hp,lv,cpn,th &         qent, uent, vent, nent, elij, sig, tv, tvp, iflag, precip, VPrecip, &
10         ,ep,clw,m,tp,mp,rp,up,vp &         ft, fr, fu, fv, upwd, dnwd, dnwd0, ma, mike, tls, tps, qcondc, wd)
11         ,wt,water,evap,b &  
12         ,ment,qent,uent,vent,nent,elij,sig &      use conema3_m, only: iflag_clw
13         ,tv,tvp &      use cv30_param_m, only: delta, minorig, nl, sigd
14         ,iflag,precip,VPrecip,ft,fr,fu,fv &      use cvthermo, only: cl, cpd, cpv, grav, rowl, rrd, rrv
15         ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)      USE dimphy, ONLY: klev, klon
     use conema3_m  
     use cv3_param_m  
     use cvthermo  
     use cvflag  
16    
17      ! inputs:      ! inputs:
18      integer, intent(in):: ncum,nd,na,nloc      integer, intent(in):: icb(:), inb(:) ! (ncum)
     integer icb(nloc), inb(nloc)  
19      real, intent(in):: delt      real, intent(in):: delt
20      real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)      real t(klon, klev), rr(klon, klev), u(klon, klev), v(klon, klev)
21      real sig(nloc,nd)      real gz(klon, klev)
22      real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)      real p(klon, klev)
23      real th(nloc,na), p(nloc,nd), tp(nloc,na)      real ph(klon, klev + 1), h(klon, klev), hp(klon, klev)
24      real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)      real lv(klon, klev), cpn(klon, klev)
25      real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)      real th(klon, klev)
26      real vp(nloc,na), wt(nloc,nd)      real ep(klon, klev), clw(klon, klev)
27      real water(nloc,na), evap(nloc,na), b(nloc,na)      real m(klon, klev)
28      real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)      real tp(klon, klev)
29      !ym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)      real mp(klon, klev), rp(klon, klev), up(klon, klev)
30      real vent(nloc,na,na), elij(nloc,na,na)      real vp(klon, klev), wt(klon, klev)
31      integer nent(nloc,na)      real, intent(in):: water(:, :), evap(:, :), b(:, :) ! (ncum, nl)
32      real tv(nloc,nd), tvp(nloc,nd)      real ment(klon, klev, klev), qent(klon, klev, klev), uent(klon, klev, klev)
33        real vent(klon, klev, klev)
34        integer nent(klon, klev)
35        real elij(klon, klev, klev)
36        real sig(klon, klev)
37        real tv(klon, klev), tvp(klon, klev)
38    
39      ! input/output:      ! input / output:
40      integer iflag(nloc)      integer iflag(klon)
41    
42      ! outputs:      ! outputs:
43      real precip(nloc)      real precip(klon)
44      real VPrecip(nloc,nd+1)      real VPrecip(klon, klev + 1)
45      real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)      real ft(klon, klev), fr(klon, klev), fu(klon, klev), fv(klon, klev)
46      real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)      real upwd(klon, klev), dnwd(klon, klev)
47      real dnwd0(nloc,nd), mike(nloc,nd)      real dnwd0(klon, klev)
48      real tls(nloc,nd), tps(nloc,nd)      real ma(klon, klev)
49      real qcondc(nloc,nd)                               ! cld      real mike(klon, klev)
50      real wd(nloc)                                      ! gust      real tls(klon, klev), tps(klon, klev)
51        real qcondc(klon, klev)
52      ! local variables:      real wd(klon) ! gust
53      integer i,k,il,n,j,num1  
54        ! Local:
55        integer ncum
56        integer i, k, il, n, j, num1
57      real rat, awat, delti      real rat, awat, delti
58      real ax, bx, cx, dx      real ax, bx, cx, dx
59      real cpinv, rdcp, dpinv      real cpinv, rdcp, dpinv
60      real lvcp(nloc,na), mke(nloc,na)      real lvcp(klon, klev)
61      real am(nloc), work(nloc), ad(nloc), amp1(nloc)      real am(klon), work(klon), ad(klon), amp1(klon)
62  !!!      real up1(nloc), dn1(nloc)      real up1(klon, klev, klev), dn1(klon, klev, klev)
63      real up1(nloc,nd,nd), dn1(nloc,nd,nd)      real asum(klon), bsum(klon), csum(klon), dsum(klon)
64      real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)      real qcond(klon, klev), nqcond(klon, klev), wa(klon, klev)
65      real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld      real siga(klon, klev), sax(klon, klev), mac(klon, klev)
     real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld  
   
66    
67      !-------------------------------------------------------------      !-------------------------------------------------------------
68    
69        ncum = size(icb)
70    
71      ! initialization:      ! initialization:
72    
73      delti = 1.0/delt      delti = 1.0 / delt
74    
75      do il=1,ncum      do il = 1, ncum
76         precip(il)=0.0         precip(il) = 0.0
77         wd(il)=0.0     ! gust         wd(il) = 0.0 ! gust
78         VPrecip(il,nd+1)=0.         VPrecip(il, klev + 1) = 0.
79      enddo      enddo
80    
81      do i=1,nd      do i = 1, klev
82         do il=1,ncum         do il = 1, ncum
83            VPrecip(il,i)=0.0            VPrecip(il, i) = 0.0
84            ft(il,i)=0.0            ft(il, i) = 0.0
85            fr(il,i)=0.0            fr(il, i) = 0.0
86            fu(il,i)=0.0            fu(il, i) = 0.0
87            fv(il,i)=0.0            fv(il, i) = 0.0
88            qcondc(il,i)=0.0                                ! cld            qcondc(il, i) = 0.0
89            qcond(il,i)=0.0                                 ! cld            qcond(il, i) = 0.0
90            nqcond(il,i)=0.0                                ! cld            nqcond(il, i) = 0.0
91         enddo         enddo
92      enddo      enddo
93    
94        do i = 1, nl
95      do i=1,nl         do il = 1, ncum
96         do il=1,ncum            lvcp(il, i) = lv(il, i) / cpn(il, i)
           lvcp(il,i)=lv(il,i)/cpn(il,i)  
97         enddo         enddo
98      enddo      enddo
99    
100        ! calculate surface precipitation in mm / day
101    
102      !      do il = 1, ncum
103      !   ***  calculate surface precipitation in mm/day     ***         if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd &
104      !              * water(il, 1) * 86400. * 1000. / (rowl * grav)
     do il=1,ncum  
        if(ep(il,inb(il)).ge.0.0001)then  
           if (cvflag_grav) then  
              precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)  
           else  
              precip(il)=wt(il,1)*sigd*water(il,1)*8640.  
           endif  
        endif  
105      enddo      enddo
106    
107      !   ***  CALCULATE VERTICAL PROFILE OF  PRECIPITATIONs IN kg/m2/s  ===      ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s ===
108      !  
109      ! MAF rajout pour lessivage      ! MAF rajout pour lessivage
110      do k=1,nl      do k = 1, nl - 1
111         do il=1,ncum         do il = 1, ncum
112            if (k.le.inb(il)) then            if (k <= inb(il)) VPrecip(il, k) = wt(il, k) * sigd * water(il, k) &
113               if (cvflag_grav) then                 / grav
                 VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav  
              else  
                 VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10.  
              endif  
           endif  
114         end do         end do
115      end do      end do
116      !  
117      !      ! calculate tendencies of lowest level potential temperature
118      !      ! and mixing ratio
119      !   ***  calculate tendencies of lowest level potential temperature  ***  
120      !   ***                      and mixing ratio                        ***      do il = 1, ncum
121      !         work(il) = 1.0 / (ph(il, 1) - ph(il, 2))
122      do il=1,ncum         am(il) = 0.0
        work(il)=1.0/(ph(il,1)-ph(il,2))  
        am(il)=0.0  
123      enddo      enddo
124    
125      do k=2,nl      do k = 2, nl
126         do il=1,ncum         do il = 1, ncum
127            if (k.le.inb(il)) then            if (k <= inb(il)) am(il) = am(il) + m(il, k)
              am(il)=am(il)+m(il,k)  
           endif  
128         enddo         enddo
129      enddo      enddo
130    
131      do il=1,ncum      do il = 1, ncum
132           ! Consist vect:
133           if (0.01 * grav * work(il) * am(il) >= delti) iflag(il) = 1
134    
135           ft(il, 1) = 0.01 * grav * work(il) * am(il) * (t(il, 2) - t(il, 1) &
136                + (gz(il, 2) - gz(il, 1)) / cpn(il, 1))
137    
138           ft(il, 1) = ft(il, 1) - 0.5 * lvcp(il, 1) * sigd * (evap(il, 1) &
139                + evap(il, 2))
140    
141           ft(il, 1) = ft(il, 1) - 0.009 * grav * sigd * mp(il, 2) &
142                * t(il, 1) * b(il, 1) * work(il)
143    
144           ft(il, 1) = ft(il, 1) + 0.01 * sigd * wt(il, 1) * (cl - cpd) &
145                * water(il, 2) * (t(il, 2) - t(il, 1)) * work(il) / cpn(il, 1)
146    
147           !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
148           ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
149           fr(il, 1) = 0.01 * grav * mp(il, 2) * (rp(il, 2) - rr(il, 1)) &
150                * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2))
151           ! + tard : + sigd * evap(il, 1)
152    
153         ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4         fr(il, 1) = fr(il, 1) + 0.01 * grav * am(il) * (rr(il, 2) - rr(il, 1)) &
154         if (cvflag_grav) then              * work(il)
           if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect  
           ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1) &  
                +(gz(il,2)-gz(il,1))/cpn(il,1))  
        else  
           if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect  
           ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) &  
                +(gz(il,2)-gz(il,1))/cpn(il,1))  
        endif  
   
        ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))  
   
        if (cvflag_grav) then  
           ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) &  
                *t(il,1)*b(il,1)*work(il)  
        else  
           ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)  
        endif  
   
        ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2) &  
             -t(il,1))*work(il)/cpn(il,1)  
   
        if (cvflag_grav) then  
           !jyg1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)  
           ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)  
           fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &  
                +sigd*0.5*(evap(il,1)+evap(il,2))  
           !+tard     :          +sigd*evap(il,1)  
   
           fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)  
   
           fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &  
                +am(il)*(u(il,2)-u(il,1)))  
           fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &  
                +am(il)*(v(il,2)-v(il,1)))  
        else  ! cvflag_grav  
           fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &  
                +sigd*0.5*(evap(il,1)+evap(il,2))  
           fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)  
           fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &  
                +am(il)*(u(il,2)-u(il,1)))  
           fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &  
                +am(il)*(v(il,2)-v(il,1)))  
        endif ! cvflag_grav  
155    
156           fu(il, 1) = fu(il, 1) + 0.01 * grav * work(il) * (mp(il, 2) &
157                * (up(il, 2) - u(il, 1)) + am(il) * (u(il, 2) - u(il, 1)))
158           fv(il, 1) = fv(il, 1) + 0.01 * grav * work(il) * (mp(il, 2) &
159                * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1)))
160      enddo ! il      enddo ! il
161    
162      do j=2,nl      do j = 2, nl
163         do il=1,ncum         do il = 1, ncum
164            if (j.le.inb(il)) then            if (j <= inb(il)) then
165               if (cvflag_grav) then               fr(il, 1) = fr(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) &
166                  fr(il,1)=fr(il,1) &                    * (qent(il, j, 1) - rr(il, 1))
167                       +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))               fu(il, 1) = fu(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) &
168                  fu(il,1)=fu(il,1) &                    * (uent(il, j, 1) - u(il, 1))
169                       +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))               fv(il, 1) = fv(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) &
170                  fv(il,1)=fv(il,1) &                    * (vent(il, j, 1) - v(il, 1))
171                       +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))            endif
              else   ! cvflag_grav  
                 fr(il,1)=fr(il,1) &  
                      +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))  
                 fu(il,1)=fu(il,1) &  
                      +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))  
                 fv(il,1)=fv(il,1) &  
                      +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))  
              endif  ! cvflag_grav  
           endif ! j  
172         enddo         enddo
173      enddo      enddo
174    
175      !      ! calculate tendencies of potential temperature and mixing ratio
176      !   ***  calculate tendencies of potential temperature and mixing ratio  ***      ! at levels above the lowest level
177      !   ***               at levels above the lowest level                   ***  
178      !      ! first find the net saturated updraft and downdraft mass fluxes
179      !   ***  first find the net saturated updraft and downdraft mass fluxes  ***      ! through each level
     !   ***                      through each level                          ***  
     !  
180    
181      do  i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?      loop_i: do i = 2, nl - 1
182           num1 = 0
183    
184         num1=0         do il = 1, ncum
185         do il=1,ncum            if (i <= inb(il)) num1 = num1 + 1
           if(i.le.inb(il))num1=num1+1  
186         enddo         enddo
        if(num1.le.0) cycle  
187    
188         call zilch(amp1,ncum)         if (num1 > 0) then
189         call zilch(ad,ncum)            amp1(:ncum) = 0.
190              ad(:ncum) = 0.
191    
192         do  k=i+1,nl+1            do k = i + 1, nl + 1
193            do  il=1,ncum               do il = 1, ncum
194               if (i.le.inb(il) .and. k.le.(inb(il)+1)) then                  if (i <= inb(il) .and. k <= (inb(il) + 1)) then
195                  amp1(il)=amp1(il)+m(il,k)                     amp1(il) = amp1(il) + m(il, k)
196               endif                  endif
197                 end do
198            end do            end do
        end do  
199    
200         do  k=1,i            do k = 1, i
201            do  j=i+1,nl+1               do j = i + 1, nl + 1
202               do  il=1,ncum                  do il = 1, ncum
203                  if (i.le.inb(il) .and. j.le.(inb(il)+1)) then                     if (i <= inb(il) .and. j <= (inb(il) + 1)) then
204                     amp1(il)=amp1(il)+ment(il,k,j)                        amp1(il) = amp1(il) + ment(il, k, j)
205                  endif                     endif
206                    end do
207               end do               end do
208            end do            end do
        end do  
209    
210         do  k=1,i-1            do k = 1, i - 1
211            do  j=i,nl+1 ! newvecto: nl au lieu nl+1?               do j = i, nl + 1 ! newvecto: nl au lieu nl + 1?
212               do  il=1,ncum                  do il = 1, ncum
213                  if (i.le.inb(il) .and. j.le.inb(il)) then                     if (i <= inb(il) .and. j <= inb(il)) then
214                     ad(il)=ad(il)+ment(il,j,k)                        ad(il) = ad(il) + ment(il, j, k)
215                  endif                     endif
216                    end do
217               end do               end do
218            end do            end do
        end do  
219    
220         do  il=1,ncum            do il = 1, ncum
221            if (i.le.inb(il)) then               if (i <= inb(il)) then
222               dpinv=1.0/(ph(il,i)-ph(il,i+1))                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
223               cpinv=1.0/cpn(il,i)                  cpinv = 1.0 / cpn(il, i)
224    
225               ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4                  ! Vecto:
226               if (cvflag_grav) then                  if (0.01 * grav * dpinv * amp1(il) >= delti) iflag(il) = 1
227                  if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto  
228               else                  ft(il, i) = 0.01 * grav * dpinv * (amp1(il) * (t(il, i + 1) &
229                  if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto                       - t(il, i) + (gz(il, i + 1) - gz(il, i)) * cpinv) &
230                         - ad(il) * (t(il, i) - t(il, i - 1) + (gz(il, i) &
231                         - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) &
232                         * (evap(il, i) + evap(il, i + 1))
233                    rat = cpn(il, i - 1) * cpinv
234                    ft(il, i) = ft(il, i) - 0.009 * grav * sigd * (mp(il, i + 1) &
235                         * t(il, i) * b(il, i) - mp(il, i) * t(il, i - 1) * rat &
236                         * b(il, i - 1)) * dpinv
237                    ft(il, i) = ft(il, i) + 0.01 * grav * dpinv * ment(il, i, i) &
238                         * (hp(il, i) - h(il, i) + t(il, i) * (cpv - cpd) &
239                         * (rr(il, i) - qent(il, i, i))) * cpinv
240    
241                    ft(il, i) = ft(il, i) + 0.01 * sigd * wt(il, i) * (cl - cpd) &
242                         * water(il, i + 1) * (t(il, i + 1) - t(il, i)) * dpinv &
243                         * cpinv
244    
245                    fr(il, i) = 0.01 * grav * dpinv * (amp1(il) * (rr(il, i + 1) &
246                         - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
247                    fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * (amp1(il) &
248                         * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) &
249                         - u(il, i - 1)))
250                    fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * (amp1(il) &
251                         * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) &
252                         - v(il, i - 1)))
253               endif               endif
254              end do
255    
256               if (cvflag_grav) then            do k = 1, i - 1
257                  ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &               do il = 1, ncum
258                       +(gz(il,i+1)-gz(il,i))*cpinv) &                  if (i <= inb(il)) then
259                       -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &                     dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
260                       -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))                     cpinv = 1.0 / cpn(il, i)
261                  rat=cpn(il,i-1)*cpinv  
262                  ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &                     awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
263                       -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv                     awat = amax1(awat, 0.0)
264                  ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &  
265                       +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv                     fr(il, i) = fr(il, i) + 0.01 * grav * dpinv &
266               else  ! cvflag_grav                          * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i))
267                  ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &                     fu(il, i) = fu(il, i) + 0.01 * grav * dpinv &
268                       +(gz(il,i+1)-gz(il,i))*cpinv) &                          * ment(il, k, i) * (uent(il, k, i) - u(il, i))
269                       -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &                     fv(il, i) = fv(il, i) + 0.01 * grav * dpinv &
270                       -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))                          * ment(il, k, i) * (vent(il, k, i) - v(il, i))
271                  rat=cpn(il,i-1)*cpinv  
272                  ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &                     ! (saturated updrafts resulting from mixing)
273                       -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv                     qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat)
274                  ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &                     nqcond(il, i) = nqcond(il, i) + 1.
275                       +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv                  endif ! i
276               endif ! cvflag_grav               end do
   
   
              ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) &  
                   *(t(il,i+1)-t(il,i))*dpinv*cpinv  
   
              if (cvflag_grav) then  
                 fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &  
                      -ad(il)*(rr(il,i)-rr(il,i-1)))  
                 fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &  
                      -ad(il)*(u(il,i)-u(il,i-1)))  
                 fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &  
                      -ad(il)*(v(il,i)-v(il,i-1)))  
              else  ! cvflag_grav  
                 fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &  
                      -ad(il)*(rr(il,i)-rr(il,i-1)))  
                 fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &  
                      -ad(il)*(u(il,i)-u(il,i-1)))  
                 fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &  
                      -ad(il)*(v(il,i)-v(il,i-1)))  
              endif ! cvflag_grav  
   
           endif ! i  
        end do  
   
        do  k=1,i-1  
           do  il=1,ncum  
              if (i.le.inb(il)) then  
                 dpinv=1.0/(ph(il,i)-ph(il,i+1))  
                 cpinv=1.0/cpn(il,i)  
   
                 awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)  
                 awat=amax1(awat,0.0)  
   
                 if (cvflag_grav) then  
                    fr(il,i)=fr(il,i) &  
                         +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))  
                    fu(il,i)=fu(il,i) &  
                         +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))  
                    fv(il,i)=fv(il,i) &  
                         +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))  
                 else  ! cvflag_grav  
                    fr(il,i)=fr(il,i) &  
                         +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))  
                    fu(il,i)=fu(il,i) &  
                         +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))  
                    fv(il,i)=fv(il,i) &  
                         +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))  
                 endif ! cvflag_grav  
   
                 ! (saturated updrafts resulting from mixing)        ! cld  
                 qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld  
                 nqcond(il,i)=nqcond(il,i)+1.                ! cld  
              endif ! i  
277            end do            end do
        end do  
278    
279         do  k=i,nl+1            do k = i, nl + 1
280            do  il=1,ncum               do il = 1, ncum
281               if (i.le.inb(il) .and. k.le.inb(il)) then                  if (i <= inb(il) .and. k <= inb(il)) then
282                  dpinv=1.0/(ph(il,i)-ph(il,i+1))                     dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
283                  cpinv=1.0/cpn(il,i)                     cpinv = 1.0 / cpn(il, i)
284    
285                  if (cvflag_grav) then                     fr(il, i) = fr(il, i) + 0.01 * grav * dpinv &
286                     fr(il,i)=fr(il,i) &                          * ment(il, k, i) * (qent(il, k, i) - rr(il, i))
287                          +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))                     fu(il, i) = fu(il, i) + 0.01 * grav * dpinv &
288                     fu(il,i)=fu(il,i) &                          * ment(il, k, i) * (uent(il, k, i) - u(il, i))
289                          +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                     fv(il, i) = fv(il, i) + 0.01 * grav * dpinv &
290                     fv(il,i)=fv(il,i) &                          * ment(il, k, i) * (vent(il, k, i) - v(il, i))
291                          +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))                  endif
292                  else  ! cvflag_grav               end do
                    fr(il,i)=fr(il,i) &  
                         +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))  
                    fu(il,i)=fu(il,i) &  
                         +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))  
                    fv(il,i)=fv(il,i) &  
                         +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))  
                 endif ! cvflag_grav  
              endif ! i and k  
293            end do            end do
        end do  
294    
295         do  il=1,ncum            do il = 1, ncum
296            if (i.le.inb(il)) then               if (i <= inb(il)) then
297               dpinv=1.0/(ph(il,i)-ph(il,i+1))                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
298               cpinv=1.0/cpn(il,i)                  cpinv = 1.0 / cpn(il, i)
299    
              if (cvflag_grav) then  
300                  ! sb: on ne fait pas encore la correction permettant de mieux                  ! sb: on ne fait pas encore la correction permettant de mieux
301                  ! conserver l'eau:                  ! conserver l'eau:
302                  fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &                  fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) &
303                       +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &                       + evap(il, i + 1)) + 0.01 * grav * (mp(il, i + 1) &
304                       *(rp(il,i)-rr(il,i-1)))*dpinv                       * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) &
305                         - rr(il, i - 1))) * dpinv
306                  fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &  
307                       -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv                  fu(il, i) = fu(il, i) + 0.01 * grav * (mp(il, i + 1) &
308                  fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &                       * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) &
309                       -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv                       - u(il, i - 1))) * dpinv
310               else  ! cvflag_grav                  fv(il, i) = fv(il, i) + 0.01 * grav * (mp(il, i + 1) &
311                  fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &                       * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) &
312                       +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &                       - v(il, i - 1))) * dpinv
313                       *(rp(il,i)-rr(il,i-1)))*dpinv               endif
314                  fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &            end do
                      -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv  
                 fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &  
                      -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv  
              endif ! cvflag_grav  
   
           endif ! i  
        end do  
   
        ! sb: interface with the cloud parameterization:          ! cld  
315    
316         do k=i+1,nl            ! sb: interface with the cloud parameterization:
           do il=1,ncum  
              if (k.le.inb(il) .and. i.le.inb(il)) then         ! cld  
                 ! (saturated downdrafts resulting from mixing)            ! cld  
                 qcond(il,i)=qcond(il,i)+elij(il,k,i)            ! cld  
                 nqcond(il,i)=nqcond(il,i)+1.                    ! cld  
              endif                                             ! cld  
           enddo                                              ! cld  
        enddo                                               ! cld  
   
        ! (particular case: no detraining level is found)         ! cld  
        do il=1,ncum                                        ! cld  
           if (i.le.inb(il) .and. nent(il,i).eq.0) then       ! cld  
              qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld  
              nqcond(il,i)=nqcond(il,i)+1.                    ! cld  
           endif                                              ! cld  
        enddo                                               ! cld  
   
        do il=1,ncum                                        ! cld  
           if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then    ! cld  
              qcond(il,i)=qcond(il,i)/nqcond(il,i)            ! cld  
           endif                                              ! cld  
        enddo  
317    
318      end do            do k = i + 1, nl
319                 do il = 1, ncum
320                    if (k <= inb(il) .and. i <= inb(il)) then
321                       ! (saturated downdrafts resulting from mixing)
322                       qcond(il, i) = qcond(il, i) + elij(il, k, i)
323                       nqcond(il, i) = nqcond(il, i) + 1.
324                    endif
325                 enddo
326              enddo
327    
328              ! (particular case: no detraining level is found)
329              do il = 1, ncum
330                 if (i <= inb(il) .and. nent(il, i) == 0) then
331                    qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i)
332                    nqcond(il, i) = nqcond(il, i) + 1.
333                 endif
334              enddo
335    
336      !   ***   move the detrainment at level inb down to level inb-1   ***            do il = 1, ncum
337      !   ***        in such a way as to preserve the vertically        ***               if (i <= inb(il) .and. nqcond(il, i) /= 0.) then
338      !   ***          integrated enthalpy and water tendencies         ***                  qcond(il, i) = qcond(il, i) / nqcond(il, i)
339      !               endif
340      do  il=1,ncum            enddo
341           end if
342        end do loop_i
343    
344         ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il)) &      ! move the detrainment at level inb down to level inb - 1
345              +t(il,inb(il))*(cpv-cpd) &      ! in such a way as to preserve the vertically
346              *(rr(il,inb(il))-qent(il,inb(il),inb(il)))) &      ! integrated enthalpy and water tendencies
347              /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))  
348         ft(il,inb(il))=ft(il,inb(il))-ax      do il = 1, ncum
349         ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) &         ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) &
350              *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) &              - h(il, inb(il)) + t(il, inb(il)) * (cpv - cpd) &
351              *(ph(il,inb(il)-1)-ph(il,inb(il))))              * (rr(il, inb(il)) - qent(il, inb(il), inb(il)))) &
352                / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1)))
353         bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il)) &         ft(il, inb(il)) = ft(il, inb(il)) - ax
354              -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))         ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) &
355         fr(il,inb(il))=fr(il,inb(il))-bx              * (ph(il, inb(il)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) &
356         fr(il,inb(il)-1)=fr(il,inb(il)-1) &              * (ph(il, inb(il) - 1) - ph(il, inb(il))))
357              +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) &  
358              /(ph(il,inb(il)-1)-ph(il,inb(il)))         bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) &
359                - rr(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
360         cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il)) &         fr(il, inb(il)) = fr(il, inb(il)) - bx
361              -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))         fr(il, inb(il) - 1) = fr(il, inb(il) - 1) &
362         fu(il,inb(il))=fu(il,inb(il))-cx              + bx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
363         fu(il,inb(il)-1)=fu(il,inb(il)-1) &              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
364              +cx*(ph(il,inb(il))-ph(il,inb(il)+1)) &  
365              /(ph(il,inb(il)-1)-ph(il,inb(il)))         cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) &
366                - u(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
367         dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il)) &         fu(il, inb(il)) = fu(il, inb(il)) - cx
368              -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))         fu(il, inb(il) - 1) = fu(il, inb(il) - 1) &
369         fv(il,inb(il))=fv(il,inb(il))-dx              + cx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
370         fv(il,inb(il)-1)=fv(il,inb(il)-1) &              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
371              +dx*(ph(il,inb(il))-ph(il,inb(il)+1)) &  
372              /(ph(il,inb(il)-1)-ph(il,inb(il)))         dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) &
373                - v(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
374           fv(il, inb(il)) = fv(il, inb(il)) - dx
375           fv(il, inb(il) - 1) = fv(il, inb(il) - 1) &
376                + dx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
377                / (ph(il, inb(il) - 1) - ph(il, inb(il)))
378    
379      end do      end do
380    
381      !      ! homoginize tendencies below cloud base
382      !   ***    homoginize tendencies below cloud base    ***  
383      !      do il = 1, ncum
384      !         asum(il) = 0.0
385      do il=1,ncum         bsum(il) = 0.0
386         asum(il)=0.0         csum(il) = 0.0
387         bsum(il)=0.0         dsum(il) = 0.0
388         csum(il)=0.0      enddo
389         dsum(il)=0.0  
390      enddo      do i = 1, nl
391           do il = 1, ncum
392      do i=1,nl            if (i <= (icb(il) - 1)) then
393         do il=1,ncum               asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
394            if (i.le.(icb(il)-1)) then               bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (cl - cpd) &
395               asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))                    * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
396               bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &               csum(il) = csum(il) + (lv(il, i) + (cl - cpd) * (t(il, i) &
397                    *(ph(il,i)-ph(il,i+1))                    - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
398               csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &               dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) &
399                    *(ph(il,i)-ph(il,i+1))                    / th(il, i)
              dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)  
400            endif            endif
401         enddo         enddo
402      enddo      enddo
403    
404  !!!!      do 700 i=1,icb(il)-1      do i = 1, nl
405      do i=1,nl         do il = 1, ncum
406         do il=1,ncum            if (i <= (icb(il) - 1)) then
407            if (i.le.(icb(il)-1)) then               ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
408               ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il))               fr(il, i) = bsum(il) / csum(il)
              fr(il,i)=bsum(il)/csum(il)  
409            endif            endif
410         enddo         enddo
411      enddo      enddo
412    
413      !      ! reset counter and return
     !   ***           reset counter and return           ***  
     !  
     do il=1,ncum  
        sig(il,nd)=2.0  
     enddo  
414    
415        do il = 1, ncum
416           sig(il, klev) = 2.0
417        enddo
418    
419      do i=1,nd      do i = 1, klev
420         do il=1,ncum         do il = 1, ncum
421            upwd(il,i)=0.0            upwd(il, i) = 0.0
422            dnwd(il,i)=0.0            dnwd(il, i) = 0.0
423         enddo         enddo
424      enddo      enddo
425    
426      do i=1,nl      do i = 1, nl
427         do il=1,ncum         do il = 1, ncum
428            dnwd0(il,i)=-mp(il,i)            dnwd0(il, i) = - mp(il, i)
429         enddo         enddo
430      enddo      enddo
431      do i=nl+1,nd      do i = nl + 1, klev
432         do il=1,ncum         do il = 1, ncum
433            dnwd0(il,i)=0.            dnwd0(il, i) = 0.
434         enddo         enddo
435      enddo      enddo
436    
437        do i = 1, nl
438      do i=1,nl         do il = 1, ncum
439         do il=1,ncum            if (i >= icb(il) .and. i <= inb(il)) then
440            if (i.ge.icb(il) .and. i.le.inb(il)) then               upwd(il, i) = 0.0
441               upwd(il,i)=0.0               dnwd(il, i) = 0.0
              dnwd(il,i)=0.0  
442            endif            endif
443         enddo         enddo
444      enddo      enddo
445    
446      do i=1,nl      do i = 1, nl
447         do k=1,nl         do k = 1, nl
448            do il=1,ncum            do il = 1, ncum
449               up1(il,k,i)=0.0               up1(il, k, i) = 0.0
450               dn1(il,k,i)=0.0               dn1(il, k, i) = 0.0
451            enddo            enddo
452         enddo         enddo
453      enddo      enddo
454    
455      do i=1,nl      do i = 1, nl
456         do k=i,nl         do k = i, nl
457            do n=1,i-1            do n = 1, i - 1
458               do il=1,ncum               do il = 1, ncum
459                  if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then                  if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then
460                     up1(il,k,i)=up1(il,k,i)+ment(il,n,k)                     up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
461                     dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)                     dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
462                  endif                  endif
463               enddo               enddo
464            enddo            enddo
465         enddo         enddo
466      enddo      enddo
467    
468      do i=2,nl      do i = 2, nl
469         do k=i,nl         do k = i, nl
470            do il=1,ncum            do il = 1, ncum
471               !test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then               if (i <= inb(il).and.k <= inb(il)) then
472               if (i.le.inb(il).and.k.le.inb(il)) then                  upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
473                  upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)                  dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
                 dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)  
474               endif               endif
475            enddo            enddo
476         enddo         enddo
477      enddo      enddo
478    
   
479      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
480      !        determination de la variation de flux ascendant entre      ! determination de la variation de flux ascendant entre
481      !        deux niveau non dilue mike      ! deux niveau non dilue mike
482      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
483    
484      do i=1,nl      do i = 1, nl
485         do il=1,ncum         do il = 1, ncum
486            mike(il,i)=m(il,i)            mike(il, i) = m(il, i)
487         enddo         enddo
488      enddo      enddo
489    
490      do i=nl+1,nd      do i = nl + 1, klev
491         do il=1,ncum         do il = 1, ncum
492            mike(il,i)=0.            mike(il, i) = 0.
493         enddo         enddo
494      enddo      enddo
495    
496      do i=1,nd      do i = 1, klev
497         do il=1,ncum         do il = 1, ncum
498            ma(il,i)=0            ma(il, i) = 0
499         enddo         enddo
500      enddo      enddo
501    
502      do i=1,nl      do i = 1, nl
503         do j=i,nl         do j = i, nl
504            do il=1,ncum            do il = 1, ncum
505               ma(il,i)=ma(il,i)+m(il,j)               ma(il, i) = ma(il, i) + m(il, j)
506            enddo            enddo
507         enddo         enddo
508      enddo      enddo
509    
510      do i=nl+1,nd      do i = nl + 1, klev
511         do il=1,ncum         do il = 1, ncum
512            ma(il,i)=0.            ma(il, i) = 0.
513         enddo         enddo
514      enddo      enddo
515    
516      do i=1,nl      do i = 1, nl
517         do il=1,ncum         do il = 1, ncum
518            if (i.le.(icb(il)-1)) then            if (i <= (icb(il) - 1)) then
519               ma(il,i)=0               ma(il, i) = 0
520            endif            endif
521         enddo         enddo
522      enddo      enddo
523    
524      !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
525      !        icb represente de niveau ou se trouve la      ! icb represente de niveau ou se trouve la
526      !        base du nuage , et inb le top du nuage      ! base du nuage, et inb le top du nuage
527      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
528    
529      do i=1,nd      do i = 1, klev
530         do il=1,ncum         DO il = 1, ncum
531            mke(il,i)=upwd(il,i)+dnwd(il,i)            rdcp = (rrd * (1. - rr(il, i)) - rr(il, i) * rrv) &
532                   / (cpd * (1. - rr(il, i)) + rr(il, i) * cpv)
533              tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp
534              tps(il, i) = tp(il, i)
535           end DO
536        enddo
537    
538        ! diagnose the in-cloud mixing ratio
539        ! of condensed water
540        !
541    
542        do i = 1, klev
543           do il = 1, ncum
544              mac(il, i) = 0.0
545              wa(il, i) = 0.0
546              siga(il, i) = 0.0
547              sax(il, i) = 0.0
548         enddo         enddo
549      enddo      enddo
550    
551      do i=1,nd      do i = minorig, nl
552         DO  il=1,ncum         do k = i + 1, nl + 1
553            rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) &            do il = 1, ncum
554                 /(cpd*(1.-rr(il,i))+rr(il,i)*cpv)               if (i <= inb(il) .and. k <= (inb(il) + 1)) then
555            tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp                  mac(il, i) = mac(il, i) + m(il, k)
556            tps(il,i)=tp(il,i)               endif
557         end DO            enddo
558           enddo
559      enddo      enddo
560    
561      !      do i = 1, nl
562      !   *** diagnose the in-cloud mixing ratio   ***            ! cld         do j = 1, i
563      !   ***           of condensed water         ***            ! cld            do il = 1, ncum
564      !                                                           ! cld               if (i >= icb(il) .and. i <= (inb(il) - 1) &
565                      .and. j >= icb(il)) then
566      do i=1,nd                                            ! cld                  sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) &
567         do il=1,ncum                                        ! cld                       * (ph(il, j) - ph(il, j + 1)) / p(il, j)
568            mac(il,i)=0.0                                      ! cld               endif
569            wa(il,i)=0.0                                       ! cld            enddo
570            siga(il,i)=0.0                                     ! cld         enddo
571            sax(il,i)=0.0                                      ! cld      enddo
572         enddo                                               ! cld  
573      enddo                                                ! cld      do i = 1, nl
574           do il = 1, ncum
575      do i=minorig, nl                                     ! cld            if (i >= icb(il) .and. i <= (inb(il) - 1) &
576         do k=i+1,nl+1                                       ! cld                 .and. sax(il, i) > 0.0) then
577            do il=1,ncum                                       ! cld               wa(il, i) = sqrt(2. * sax(il, i))
              if (i.le.inb(il) .and. k.le.(inb(il)+1)) then     ! cld  
                 mac(il,i)=mac(il,i)+m(il,k)                     ! cld  
              endif                                             ! cld  
           enddo                                              ! cld  
        enddo                                               ! cld  
     enddo                                                ! cld  
   
     do i=1,nl                                            ! cld  
        do j=1,i                                            ! cld  
           do il=1,ncum                                       ! cld  
              if (i.ge.icb(il) .and. i.le.(inb(il)-1) &  
                   .and. j.ge.icb(il) ) then                       ! cld  
                 sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) &  
                      *(ph(il,j)-ph(il,j+1))/p(il,j)                ! cld  
              endif                                             ! cld  
           enddo                                              ! cld  
        enddo                                               ! cld  
     enddo                                                ! cld  
   
     do i=1,nl                                            ! cld  
        do il=1,ncum                                        ! cld  
           if (i.ge.icb(il) .and. i.le.(inb(il)-1) &  
                .and. sax(il,i).gt.0.0 ) then                  ! cld  
              wa(il,i)=sqrt(2.*sax(il,i))                      ! cld  
           endif                                              ! cld  
        enddo                                               ! cld  
     enddo                                                ! cld  
   
     do i=1,nl                                            ! cld  
        do il=1,ncum                                        ! cld  
           if (wa(il,i).gt.0.0) &  
                siga(il,i)=mac(il,i)/wa(il,i) &  
                *rrd*tvp(il,i)/p(il,i)/100./delta            ! cld  
           siga(il,i) = min(siga(il,i),1.0)                  ! cld  
           !IM cf. FH  
           if (iflag_clw.eq.0) then  
              qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) &  
                   + (1.-siga(il,i))*qcond(il,i)              ! cld  
           else if (iflag_clw.eq.1) then  
              qcondc(il,i)=qcond(il,i)              ! cld  
578            endif            endif
579           enddo
580        enddo
581    
582         enddo                                               ! cld      do i = 1, nl
583      enddo                                                ! cld         do il = 1, ncum
584              if (wa(il, i) > 0.0) siga(il, i) = mac(il, i) / wa(il, i) * rrd &
585                   * tvp(il, i) / p(il, i) / 100. / delta
586              siga(il, i) = min(siga(il, i), 1.0)
587    
588              if (iflag_clw == 0) then
589                 qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) &
590                      + (1. - siga(il, i)) * qcond(il, i)
591              else if (iflag_clw == 1) then
592                 qcondc(il, i) = qcond(il, i)
593              endif
594           enddo
595        enddo
596    
597    end SUBROUTINE cv3_yield    end SUBROUTINE cv30_yield
598    
599  end module cv3_yield_m  end module cv30_yield_m

Legend:
Removed from v.105  
changed lines
  Added in v.188

  ViewVC Help
Powered by ViewVC 1.1.21