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

Legend:
Removed from v.178  
changed lines
  Added in v.190

  ViewVC Help
Powered by ViewVC 1.1.21