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

  ViewVC Help
Powered by ViewVC 1.1.21