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

Legend:
Removed from v.184  
changed lines
  Added in v.195

  ViewVC Help
Powered by ViewVC 1.1.21