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

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

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

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

Legend:
Removed from v.76  
changed lines
  Added in v.190

  ViewVC Help
Powered by ViewVC 1.1.21