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

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

  ViewVC Help
Powered by ViewVC 1.1.21