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

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

  ViewVC Help
Powered by ViewVC 1.1.21