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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21