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

Legend:
Removed from v.76  
changed lines
  Added in v.205

  ViewVC Help
Powered by ViewVC 1.1.21