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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21