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

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

  ViewVC Help
Powered by ViewVC 1.1.21