/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_yield.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/CV30_routines/cv30_yield.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.69  
changed lines
  Added in v.200

  ViewVC Help
Powered by ViewVC 1.1.21