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

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

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

trunk/libf/phylmd/CV3_routines/cv3_yield.f90 revision 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/Sources/phylmd/CV30_routines/cv30_yield.f revision 200 by guez, Thu Jun 2 15:40:30 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, intent(in):: t(klon, klev), rr(klon, klev)
25        real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)      real, intent(in):: u(klon, klev), v(klon, klev)
26        real th(nloc,na), p(nloc,nd), tp(nloc,na)      real gz(klon, klev)
27        real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)      real p(klon, klev)
28        real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)      real ph(klon, klev + 1), h(klon, klev), hp(klon, klev)
29        real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)      real lv(klon, klev), cpn(klon, klev)
30        real water(nloc,na), evap(nloc,na), b(nloc,na)      real th(klon, klev)
31        real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)      real ep(klon, klev), clw(klon, klev)
32  !ym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)      real m(klon, klev)
33        real vent(nloc,na,na), elij(nloc,na,na)      real tp(klon, klev)
34        integer nent(nloc,na)      real mp(klon, klev), rp(klon, klev), up(klon, klev)
35        real traent(nloc,na,na,ntra)      real, intent(in):: vp(:, 2:) ! (ncum, 2:nl)
36        real tv(nloc,nd), tvp(nloc,nd)      real, intent(in):: wt(:, :) ! (ncum, nl - 1)
37        real, intent(in):: water(:, :), evap(:, :) ! (ncum, nl)
38  ! input/output:      real, intent(in):: b(:, :) ! (ncum, nl - 1)
39        integer iflag(nloc)      real ment(klon, klev, klev), qent(klon, klev, klev), uent(klon, klev, klev)
40        real vent(klon, klev, klev)
41  ! outputs:      integer nent(klon, klev)
42        real precip(nloc)      real elij(klon, klev, klev)
43        real VPrecip(nloc,nd+1)      real sig(klon, klev)
44        real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)      real tv(klon, klev), tvp(klon, klev)
45        real ftra(nloc,nd,ntra)  
46        real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)      ! outputs:
47        real dnwd0(nloc,nd), mike(nloc,nd)      integer, intent(out):: iflag(:) ! (ncum)
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
62        real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld      real 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         if (0.01 * rg * work(il) * am(il) >= delti) iflag(il) = 1
138    
139        do k=2,nl         ft(il, 1) = 0.01 * rg * work(il) * am(il) * (t(il, 2) - t(il, 1) &
140         do il=1,ncum              + (gz(il, 2) - gz(il, 1)) / cpn(il, 1)) - 0.5 * lvcp(il, 1) &
141          if (k.le.inb(il)) then              * sigd * (evap(il, 1) + evap(il, 2)) - 0.009 * rg * sigd &
142           am(il)=am(il)+m(il,k)              * mp(il, 2) * t(il, 1) * b(il, 1) * work(il) + 0.01 * sigd &
143          endif              * wt(il, 1) * (cl - cpd) * water(il, 2) * (t(il, 2) - t(il, 1)) &
144         enddo              * work(il) / cpn(il, 1)
145        enddo  
146           !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
147        do il=1,ncum         ! (sb: pour l'instant, on ne fait que le chgt concernant rg, pas evap)
148           fr(il, 1) = 0.01 * rg * mp(il, 2) * (rp(il, 2) - rr(il, 1)) &
149  ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4              * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2))
150        if (cvflag_grav) then         ! + tard : + sigd * evap(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) &         fr(il, 1) = fr(il, 1) + 0.01 * rg * am(il) * (rr(il, 2) - rr(il, 1)) &
153                    +(gz(il,2)-gz(il,1))/cpn(il,1))              * work(il)
154        else  
155         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) &
156         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)))
157                    +(gz(il,2)-gz(il,1))/cpn(il,1))         fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
158        endif              * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1)))
159        enddo
160        ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))  
161        do j = 2, nl
162        if (cvflag_grav) then         do il = 1, ncum
163         ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) &            if (j <= inb(il)) then
164                                     *t(il,1)*b(il,1)*work(il)               fr(il, 1) = fr(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
165        else                    * (qent(il, j, 1) - rr(il, 1))
166         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) &
167        endif                    * (uent(il, j, 1) - u(il, 1))
168                 fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
169        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))
170        -t(il,1))*work(il)/cpn(il,1)            endif
171           enddo
172        if (cvflag_grav) then      enddo
173  !jyg1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)  
174  ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)      ! calculate tendencies of potential temperature and mixing ratio
175         fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &      ! at levels above the lowest level
176                  +sigd*0.5*(evap(il,1)+evap(il,2))  
177  !+tard     :          +sigd*evap(il,1)      ! first find the net saturated updraft and downdraft mass fluxes
178        ! through each level
179         fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)  
180        loop_i: do i = 2, nl - 1
181         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
182                 +am(il)*(u(il,2)-u(il,1)))            amp1(:ncum) = 0.
183         fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &            ad(:ncum) = 0.
184                 +am(il)*(v(il,2)-v(il,1)))  
185        else  ! cvflag_grav            do k = i + 1, nl + 1
186         fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &               do il = 1, ncum
187                  +sigd*0.5*(evap(il,1)+evap(il,2))                  if (i <= inb(il) .and. k <= (inb(il) + 1)) then
188         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)
189         fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &                  endif
190                 +am(il)*(u(il,2)-u(il,1)))               end do
191         fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &            end do
192                 +am(il)*(v(il,2)-v(il,1)))  
193        endif ! cvflag_grav            do k = 1, i
194                 do j = i + 1, nl + 1
195        enddo ! il                  do il = 1, ncum
196                       if (i <= inb(il) .and. j <= (inb(il) + 1)) then
197  !      do j=1,ntra                        amp1(il) = amp1(il) + ment(il, k, j)
198  !       do il=1,ncum                     endif
199  !        if (cvflag_grav) then                  end do
200  !         ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)               end do
201  !     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))            end do
202  !     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))  
203  !        else            do k = 1, i - 1
204  !         ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)               do j = i, nl + 1 ! newvecto: nl au lieu nl + 1?
205  !     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))                  do il = 1, ncum
206  !     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))                     if (i <= inb(il) .and. j <= inb(il)) then
207  !        endif                        ad(il) = ad(il) + ment(il, j, k)
208  !       enddo                     endif
209  !      enddo                  end do
210                 end do
211        do j=2,nl            end do
212         do il=1,ncum  
213          if (j.le.inb(il)) then            do il = 1, ncum
214           if (cvflag_grav) then               if (i <= inb(il)) then
215            fr(il,1)=fr(il,1) &                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
216               +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))                  cpinv = 1.0 / cpn(il, i)
217            fu(il,1)=fu(il,1) &  
218               +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
219            fv(il,1)=fv(il,1) &  
220               +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) &
221           else   ! cvflag_grav                       - t(il, i) + (gz(il, i + 1) - gz(il, i)) * cpinv) &
222            fr(il,1)=fr(il,1) &                       - ad(il) * (t(il, i) - t(il, i - 1) + (gz(il, i) &
223                 +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) &
224            fu(il,1)=fu(il,1) &                       * (evap(il, i) + evap(il, i + 1)) - 0.009 * rg * sigd &
225                 +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) &
226            fv(il,1)=fv(il,1) &                       * t(il, i - 1) * cpn(il, i - 1) * cpinv * b(il, i - 1)) &
227                 +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))                       * dpinv + 0.01 * rg * dpinv * ment(il, i, i) &
228           endif  ! cvflag_grav                       * (hp(il, i) - h(il, i) + t(il, i) * (cpv - cpd) &
229          endif ! j                       * (rr(il, i) - qent(il, i, i))) * cpinv + 0.01 * sigd &
230         enddo                       * wt(il, i) * (cl - cpd) * water(il, i + 1) &
231        enddo                       * (t(il, i + 1) - t(il, i)) * dpinv * cpinv
232                    fr(il, i) = 0.01 * rg * dpinv * (amp1(il) * (rr(il, i + 1) &
233  !      do k=1,ntra                       - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
234  !       do j=2,nl                  fu(il, i) = fu(il, i) + 0.01 * rg * dpinv * (amp1(il) &
235  !        do il=1,ncum                       * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) &
236  !         if (j.le.inb(il)) then                       - u(il, i - 1)))
237                    fv(il, i) = fv(il, i) + 0.01 * rg * dpinv * (amp1(il) &
238  !          if (cvflag_grav) then                       * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) &
239  !           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)                       - v(il, i - 1)))
240  !     :                *(traent(il,j,1,k)-tra(il,1,k))               endif
241  !          else            end do
242  !           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)  
243  !     :                *(traent(il,j,1,k)-tra(il,1,k))            do k = 1, i - 1
244  !          endif               do il = 1, ncum
245                    if (i <= inb(il)) then
246  !         endif                     dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
247  !        enddo                     cpinv = 1.0 / cpn(il, i)
248  !       enddo  
249  !      enddo                     awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
250                       awat = amax1(awat, 0.0)
251  !  
252  !   ***  calculate tendencies of potential temperature and mixing ratio  ***                     fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
253  !   ***               at levels above the lowest level                   ***                          * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i))
254  !                     fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
255  !   ***  first find the net saturated updraft and downdraft mass fluxes  ***                          * ment(il, k, i) * (uent(il, k, i) - u(il, i))
256  !   ***                      through each level                          ***                     fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
257  !                          * ment(il, k, i) * (vent(il, k, i) - v(il, i))
258    
259        do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?                     ! (saturated updrafts resulting from mixing)
260                       qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat)
261         num1=0                     nqcond(il, i) = nqcond(il, i) + 1.
262         do il=1,ncum                  endif ! i
263          if(i.le.inb(il))num1=num1+1               end do
264         enddo            end do
265         if(num1.le.0)go to 500  
266              do k = i, nl + 1
267         call zilch(amp1,ncum)               do il = 1, ncum
268         call zilch(ad,ncum)                  if (i <= inb(il) .and. k <= inb(il)) then
269                       dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
270        do 440 k=i+1,nl+1                     cpinv = 1.0 / cpn(il, i)
271         do 441 il=1,ncum  
272          if (i.le.inb(il) .and. k.le.(inb(il)+1)) then                     fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
273           amp1(il)=amp1(il)+m(il,k)                          * ment(il, k, i) * (qent(il, k, i) - rr(il, i))
274          endif                     fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
275   441   continue                          * ment(il, k, i) * (uent(il, k, i) - u(il, i))
276   440  continue                     fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
277                            * ment(il, k, i) * (vent(il, k, i) - v(il, i))
278        do 450 k=1,i                  endif
279         do 451 j=i+1,nl+1               end do
280          do 452 il=1,ncum            end do
281           if (i.le.inb(il) .and. j.le.(inb(il)+1)) then  
282            amp1(il)=amp1(il)+ment(il,k,j)            do il = 1, ncum
283           endif               if (i <= inb(il)) then
284  452     continue                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
285  451    continue                  cpinv = 1.0 / cpn(il, i)
286  450   continue  
287                    ! sb: on ne fait pas encore la correction permettant de mieux
288        do 470 k=1,i-1                  ! conserver l'eau:
289         do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?                  fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) &
290          do 472 il=1,ncum                       + evap(il, i + 1)) + 0.01 * rg * (mp(il, i + 1) &
291          if (i.le.inb(il) .and. j.le.inb(il)) then                       * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) &
292           ad(il)=ad(il)+ment(il,j,k)                       - rr(il, i - 1))) * dpinv
293          endif  
294  472     continue                  fu(il, i) = fu(il, i) + 0.01 * rg * (mp(il, i + 1) &
295  471    continue                       * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) &
296  470   continue                       - u(il, i - 1))) * dpinv
297                    fv(il, i) = fv(il, i) + 0.01 * rg * (mp(il, i + 1) &
298        do 1350 il=1,ncum                       * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) &
299        if (i.le.inb(il)) then                       - v(il, i - 1))) * dpinv
300         dpinv=1.0/(ph(il,i)-ph(il,i+1))               endif
301         cpinv=1.0/cpn(il,i)            end do
302    
303  ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4            ! sb: interface with the cloud parameterization:
304        if (cvflag_grav) then  
305         if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto            do k = i + 1, nl
306        else               do il = 1, ncum
307         if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto                  if (k <= inb(il) .and. i <= inb(il)) then
308        endif                     ! (saturated downdrafts resulting from mixing)
309                       qcond(il, i) = qcond(il, i) + elij(il, k, i)
310        if (cvflag_grav) then                     nqcond(il, i) = nqcond(il, i) + 1.
311         ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &                  endif
312            +(gz(il,i+1)-gz(il,i))*cpinv) &               enddo
313            -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &            enddo
314            -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))  
315         rat=cpn(il,i-1)*cpinv            ! (particular case: no detraining level is found)
316         ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &            do il = 1, ncum
317           -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv               if (i <= inb(il) .and. nent(il, i) == 0) then
318         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)
319            +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv                  nqcond(il, i) = nqcond(il, i) + 1.
320        else  ! cvflag_grav               endif
321         ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &            enddo
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 il = 1, ncum
324            -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))               if (i <= inb(il) .and. nqcond(il, i) /= 0.) then
325         rat=cpn(il,i-1)*cpinv                  qcond(il, i) = qcond(il, i) / nqcond(il, i)
326         ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &               endif
327           -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv            enddo
328         ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &         end if
329            +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv      end do loop_i
330        endif ! cvflag_grav  
331        ! move the detrainment at level inb down to level inb - 1
332        ! in such a way as to preserve the vertically
333        ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) &      ! integrated enthalpy and water tendencies
334                   *(t(il,i+1)-t(il,i))*dpinv*cpinv  
335        do il = 1, ncum
336        if (cvflag_grav) then         ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) &
337         fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &              - h(il, inb(il)) + t(il, inb(il)) * (cpv - cpd) &
338                   -ad(il)*(rr(il,i)-rr(il,i-1)))              * (rr(il, inb(il)) - qent(il, inb(il), inb(il)))) &
339         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)))
340                     -ad(il)*(u(il,i)-u(il,i-1)))         ft(il, inb(il)) = ft(il, inb(il)) - ax
341         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)) &
342                     -ad(il)*(v(il,i)-v(il,i-1)))              * (ph(il, inb(il)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) &
343        else  ! cvflag_grav              * (ph(il, inb(il) - 1) - ph(il, inb(il))))
344         fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &  
345                   -ad(il)*(rr(il,i)-rr(il,i-1)))         bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) &
346         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))
347                     -ad(il)*(u(il,i)-u(il,i-1)))         fr(il, inb(il)) = fr(il, inb(il)) - bx
348         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) &
349                     -ad(il)*(v(il,i)-v(il,i-1)))              + bx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
350        endif ! cvflag_grav              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
351    
352        endif ! i         cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) &
353  1350  continue              - u(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
354           fu(il, inb(il)) = fu(il, inb(il)) - cx
355  !      do k=1,ntra         fu(il, inb(il) - 1) = fu(il, inb(il) - 1) &
356  !       do il=1,ncum              + cx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
357  !        if (i.le.inb(il)) then              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
358  !         dpinv=1.0/(ph(il,i)-ph(il,i+1))  
359  !         cpinv=1.0/cpn(il,i)         dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) &
360  !         if (cvflag_grav) then              - v(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
361  !           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv         fv(il, inb(il)) = fv(il, inb(il)) - dx
362  !     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))         fv(il, inb(il) - 1) = fv(il, inb(il) - 1) &
363  !     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))              + dx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
364  !         else              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
365  !           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv  
366  !     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))      end do
367  !     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))  
368  !         endif      ! homoginize tendencies below cloud base
369  !        endif  
370  !       enddo      do il = 1, ncum
371  !      enddo         asum(il) = 0.0
372           bsum(il) = 0.0
373        do 480 k=1,i-1         csum(il) = 0.0
374         do 1370 il=1,ncum         dsum(il) = 0.0
375          if (i.le.inb(il)) then      enddo
376           dpinv=1.0/(ph(il,i)-ph(il,i+1))  
377           cpinv=1.0/cpn(il,i)      do i = 1, nl
378           do il = 1, ncum
379        awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)            if (i <= (icb(il) - 1)) then
380        awat=amax1(awat,0.0)               asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
381                 bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (cl - cpd) &
382        if (cvflag_grav) then                    * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
383        fr(il,i)=fr(il,i) &               csum(il) = csum(il) + (lv(il, i) + (cl - cpd) * (t(il, i) &
384           +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))
385        fu(il,i)=fu(il,i) &               dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) &
386                 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                    / th(il, i)
387        fv(il,i)=fv(il,i) &            endif
388                 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))         enddo
389        else  ! cvflag_grav      enddo
390        fr(il,i)=fr(il,i) &  
391           +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))      do i = 1, nl
392        fu(il,i)=fu(il,i) &         do il = 1, ncum
393           +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))            if (i <= (icb(il) - 1)) then
394        fv(il,i)=fv(il,i) &               ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
395           +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))               fr(il, i) = bsum(il) / csum(il)
396        endif ! cvflag_grav            endif
397           enddo
398  ! (saturated updrafts resulting from mixing)        ! cld      enddo
399          qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld  
400          nqcond(il,i)=nqcond(il,i)+1.                ! cld      ! reset counter and return
401        endif ! i  
402  1370  continue      do il = 1, ncum
403  480   continue         sig(il, klev) = 2.0
404        enddo
405  !      do j=1,ntra  
406  !       do k=1,i-1      do i = 1, klev
407  !        do il=1,ncum         do il = 1, ncum
408  !         if (i.le.inb(il)) then            upwd(il, i) = 0.0
409  !          dpinv=1.0/(ph(il,i)-ph(il,i+1))            dnwd(il, i) = 0.0
410  !          cpinv=1.0/cpn(il,i)         enddo
411  !          if (cvflag_grav) then      enddo
412  !           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)  
413  !     :        *(traent(il,k,i,j)-tra(il,i,j))      do i = 1, nl
414  !          else         do il = 1, ncum
415  !           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)            dnwd0(il, i) = - mp(il, i)
416  !     :        *(traent(il,k,i,j)-tra(il,i,j))         enddo
417  !          endif      enddo
418  !         endif      do i = nl + 1, klev
419  !        enddo         do il = 1, ncum
420  !       enddo            dnwd0(il, i) = 0.
421  !      enddo         enddo
422        enddo
423        do 490 k=i,nl+1  
424         do 1380 il=1,ncum      do i = 1, nl
425          if (i.le.inb(il) .and. k.le.inb(il)) then         do il = 1, ncum
426           dpinv=1.0/(ph(il,i)-ph(il,i+1))            if (i >= icb(il) .and. i <= inb(il)) then
427           cpinv=1.0/cpn(il,i)               upwd(il, i) = 0.0
428                 dnwd(il, i) = 0.0
429           if (cvflag_grav) then            endif
430           fr(il,i)=fr(il,i) &         enddo
431                 +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))      enddo
432           fu(il,i)=fu(il,i) &  
433                 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))      do i = 1, nl
434           fv(il,i)=fv(il,i) &         do k = 1, nl
435                 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))            do il = 1, ncum
436           else  ! cvflag_grav               up1(il, k, i) = 0.0
437           fr(il,i)=fr(il,i) &               dn1(il, k, i) = 0.0
438                 +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))            enddo
439           fu(il,i)=fu(il,i) &         enddo
440                 +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))      enddo
441           fv(il,i)=fv(il,i) &  
442                 +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))      do i = 1, nl
443           endif ! cvflag_grav         do k = i, nl
444          endif ! i and k            do n = 1, i - 1
445  1380   continue               do il = 1, ncum
446  490   continue                  if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then
447                       up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
448  !      do j=1,ntra                     dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
449  !       do k=i,nl+1                  endif
450  !        do il=1,ncum               enddo
451  !         if (i.le.inb(il) .and. k.le.inb(il)) then            enddo
452  !          dpinv=1.0/(ph(il,i)-ph(il,i+1))         enddo
453  !          cpinv=1.0/cpn(il,i)      enddo
454  !          if (cvflag_grav) then  
455  !           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)      do i = 2, nl
456  !     :         *(traent(il,k,i,j)-tra(il,i,j))         do k = i, nl
457  !          else            do il = 1, ncum
458  !           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)               if (i <= inb(il).and.k <= inb(il)) then
459  !     :             *(traent(il,k,i,j)-tra(il,i,j))                  upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
460  !          endif                  dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
461  !         endif ! i and k               endif
462  !        enddo            enddo
463  !       enddo         enddo
464  !      enddo      enddo
465    
466        do 1400 il=1,ncum      ! D\'etermination de la variation de flux ascendant entre
467         if (i.le.inb(il)) then      ! deux niveaux non dilu\'es mike
468          dpinv=1.0/(ph(il,i)-ph(il,i+1))  
469          cpinv=1.0/cpn(il,i)      do i = 1, nl
470           do il = 1, ncum
471          if (cvflag_grav) then            mike(il, i) = m(il, i)
472  ! sb: on ne fait pas encore la correction permettant de mieux         enddo
473  ! conserver l'eau:      enddo
474           fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &  
475                +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &      do i = nl + 1, klev
476                       *(rp(il,i)-rr(il,i-1)))*dpinv         do il = 1, ncum
477              mike(il, i) = 0.
478           fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &         enddo
479                     -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv      enddo
480           fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &  
481                     -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv      do i = 1, klev
482          else  ! cvflag_grav         do il = 1, ncum
483           fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &            ma(il, i) = 0
484                +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &         enddo
485                       *(rp(il,i)-rr(il,i-1)))*dpinv      enddo
486           fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &  
487                     -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv      do i = 1, nl
488           fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &         do j = i, nl
489                     -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv            do il = 1, ncum
490          endif ! cvflag_grav               ma(il, i) = ma(il, i) + m(il, j)
491              enddo
492        endif ! i         enddo
493  1400  continue      enddo
494    
495  ! sb: interface with the cloud parameterization:          ! cld      do i = nl + 1, klev
496           do il = 1, ncum
497        do k=i+1,nl            ma(il, i) = 0.
498         do il=1,ncum         enddo
499          if (k.le.inb(il) .and. i.le.inb(il)) then         ! cld      enddo
500  ! (saturated downdrafts resulting from mixing)            ! cld  
501            qcond(il,i)=qcond(il,i)+elij(il,k,i)            ! cld      do i = 1, nl
502            nqcond(il,i)=nqcond(il,i)+1.                    ! cld         do il = 1, ncum
503          endif                                             ! cld            if (i <= (icb(il) - 1)) then
504         enddo                                              ! cld               ma(il, i) = 0
505        enddo                                               ! cld            endif
506           enddo
507  ! (particular case: no detraining level is found)         ! cld      enddo
508        do il=1,ncum                                        ! cld  
509         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
510            qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld      ! inb le sommet du nuage
511            nqcond(il,i)=nqcond(il,i)+1.                    ! cld  
512         endif                                              ! cld      do i = 1, klev
513        enddo                                               ! cld         DO il = 1, ncum
514              rdcp = (rrd * (1. - rr(il, i)) - rr(il, i) * rrv) &
515        do il=1,ncum                                        ! cld                 / (cpd * (1. - rr(il, i)) + rr(il, i) * cpv)
516         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
517            qcond(il,i)=qcond(il,i)/nqcond(il,i)            ! cld            tps(il, i) = tp(il, i)
518         endif                                              ! cld         end DO
519        enddo      enddo
520    
521  !      do j=1,ntra      ! Diagnose the in-cloud mixing ratio of condensed water
522  !       do il=1,ncum  
523  !        if (i.le.inb(il)) then      do i = 1, klev
524  !         dpinv=1.0/(ph(il,i)-ph(il,i+1))         do il = 1, ncum
525  !         cpinv=1.0/cpn(il,i)            mac(il, i) = 0.0
526              wa(il, i) = 0.0
527  !         if (cvflag_grav) then            siga(il, i) = 0.0
528  !          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv            sax(il, i) = 0.0
529  !     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))         enddo
530  !     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))      enddo
531  !         else  
532  !          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv      do i = minorig, nl
533  !     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))         do k = i + 1, nl + 1
534  !     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))            do il = 1, ncum
535  !         endif               if (i <= inb(il) .and. k <= (inb(il) + 1)) then
536  !        endif ! i                  mac(il, i) = mac(il, i) + m(il, k)
537  !       enddo               endif
538  !      enddo            enddo
539           enddo
540  500   continue      enddo
541    
542        do i = 1, nl
543  !   ***   move the detrainment at level inb down to level inb-1   ***         do j = 1, i
544  !   ***        in such a way as to preserve the vertically        ***            do il = 1, ncum
545  !   ***          integrated enthalpy and water tendencies         ***               if (i >= icb(il) .and. i <= (inb(il) - 1) &
546  !                    .and. j >= icb(il)) then
547        do 503 il=1,ncum                  sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) &
548                         * (ph(il, j) - ph(il, j + 1)) / p(il, j)
549        ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il)) &               endif
550         +t(il,inb(il))*(cpv-cpd) &            enddo
551         *(rr(il,inb(il))-qent(il,inb(il),inb(il)))) &         enddo
552          /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))      enddo
553        ft(il,inb(il))=ft(il,inb(il))-ax  
554        ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) &      do i = 1, nl
555            *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) &         do il = 1, ncum
556            *(ph(il,inb(il)-1)-ph(il,inb(il))))            if (i >= icb(il) .and. i <= (inb(il) - 1) &
557                   .and. sax(il, i) > 0.0) then
558        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)  
559            endif            endif
          enddo  
         enddo  
560         enddo         enddo
561        enddo      enddo
562    
563        do i=2,nl      do i = 1, nl
564         do k=i,nl         do il = 1, ncum
565          do il=1,ncum            if (wa(il, i) > 0.0) siga(il, i) = mac(il, i) / wa(il, i) * rrd &
566  !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
567           if (i.le.inb(il).and.k.le.inb(il)) then            siga(il, i) = min(siga(il, i), 1.0)
568              upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)  
569              dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)            if (iflag_clw == 0) then
570           endif               qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) &
571          enddo                    + (1. - siga(il, i)) * qcond(il, i)
572         enddo            else if (iflag_clw == 1) then
573        enddo               qcondc(il, i) = qcond(il, i)
574              endif
575           enddo
576  !!!!      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  
577    
578          enddo                                               ! cld    end SUBROUTINE cv30_yield
        enddo                                                ! cld  
579    
580          return  end module cv30_yield_m
         end  

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

  ViewVC Help
Powered by ViewVC 1.1.21