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