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

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

  ViewVC Help
Powered by ViewVC 1.1.21