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

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

  ViewVC Help
Powered by ViewVC 1.1.21