/[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

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