/[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/libf/phylmd/CV3_routines/cv3_yield.f90 revision 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/phylmd/CV30_routines/cv30_yield.f revision 254 by guez, Mon Feb 5 10:39:38 2018 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 cvparam3      ! processes, etc.
14              use cvthermo  
15              use cvflag      use conf_phys_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 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  !      do j=1,ntra      ! initialization:
92  !       do i=1,nd  
93  !        do il=1,ncum      delti = 1.0 / delt
94  !          ftra(il,i,j)=0.0  
95  !        enddo      do il = 1, ncum
96  !       enddo         precip(il) = 0.0
97  !      enddo         VPrecip(il, klev + 1) = 0.
98        enddo
99        do i=1,nl  
100         do il=1,ncum      do i = 1, klev
101           lvcp(il,i)=lv(il,i)/cpn(il,i)         do il = 1, ncum
102         enddo            VPrecip(il, i) = 0.0
103        enddo            ft(il, i) = 0.0
104              fr(il, i) = 0.0
105              fu(il, i) = 0.0
106  !            fv(il, i) = 0.0
107  !   ***  calculate surface precipitation in mm/day     ***            qcondc(il, i) = 0.0
108  !            qcond(il, i) = 0.0
109        do il=1,ncum            nqcond(il, i) = 0.0
110         if(ep(il,inb(il)).ge.0.0001)then         enddo
111          if (cvflag_grav) then      enddo
112           precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)  
113          else      do i = 1, nl
114           precip(il)=wt(il,1)*sigd*water(il,1)*8640.         do il = 1, ncum
115          endif            lvcp(il, i) = lv(il, i) / cpn(il, i)
116         endif         enddo
117        enddo      enddo
118    
119  !   ***  CALCULATE VERTICAL PROFILE OF  PRECIPITATIONs IN kg/m2/s  ===      ! calculate surface precipitation in mm / day
120  !  
121  ! MAF rajout pour lessivage      do il = 1, ncum
122         do k=1,nl         if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd &
123           do il=1,ncum              * water(il, 1) * 86400. * 1000. / (rowl * rg)
124            if (k.le.inb(il)) then      enddo
125              if (cvflag_grav) then  
126               VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav      ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s ===
127              else  
128               VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10.      ! MAF rajout pour lessivage
129              endif      do k = 1, nl - 1
130            endif         do il = 1, ncum
131           end do            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 downdraft velocity scale    ***      ! calculate tendencies of lowest level potential temperature
137  !   ***  NE PAS UTILISER POUR L'INSTANT ***      ! and mixing ratio
138  !  
139  !!      do il=1,ncum      do il = 1, ncum
140  !!        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))         work(il) = 1.0 / (ph(il, 1) - ph(il, 2))
141  !!     :                                  /(sigd*p(il,icb(il)))         am(il) = 0.0
142  !!      enddo      enddo
143    
144  !      do k = 2, nl
145  !   ***  calculate tendencies of lowest level potential temperature  ***         do il = 1, ncum
146  !   ***                      and mixing ratio                        ***            if (k <= inb(il)) am(il) = am(il) + m(il, k)
147  !         enddo
148        do il=1,ncum      enddo
149         work(il)=1.0/(ph(il,1)-ph(il,2))  
150         am(il)=0.0      do il = 1, ncum
151        enddo         if (0.01 * rg * work(il) * am(il) >= delti) iflag(il) = 1
152    
153        do k=2,nl         ft(il, 1) = 0.01 * rg * work(il) * am(il) * (t(il, 2) - t(il, 1) &
154         do il=1,ncum              + (gz(il, 2) - gz(il, 1)) / cpn(il, 1)) - 0.5 * lvcp(il, 1) &
155          if (k.le.inb(il)) then              * sigd * (evap(il, 1) + evap(il, 2)) - 0.009 * rg * sigd &
156           am(il)=am(il)+m(il,k)              * mp(il, 2) * t(il, 1) * b(il, 1) * work(il) + 0.01 * sigd &
157          endif              * wt(il, 1) * (rcw - rcpd) * water(il, 2) * (t(il, 2) - t(il, 1)) &
158         enddo              * work(il) / cpn(il, 1)
159        enddo  
160           ! jyg Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
161        do il=1,ncum         ! (sb: pour l'instant, on ne fait que le chgt concernant rg, pas evap)
162           fr(il, 1) = 0.01 * rg * mp(il, 2) * (qp(il, 2) - rr(il, 1)) &
163  ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4              * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2))
164        if (cvflag_grav) then         ! + tard : + sigd * evap(il, 1)
165        if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect  
166         ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1) &         fr(il, 1) = fr(il, 1) + 0.01 * rg * am(il) * (rr(il, 2) - rr(il, 1)) &
167                    +(gz(il,2)-gz(il,1))/cpn(il,1))              * work(il)
168        else  
169         if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect         fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
170         ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) &              * (up(il, 2) - u(il, 1)) + am(il) * (u(il, 2) - u(il, 1)))
171                    +(gz(il,2)-gz(il,1))/cpn(il,1))         fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
172        endif              * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1)))
173        enddo
174        ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))  
175        do j = 2, nl
176        if (cvflag_grav) then         do il = 1, ncum
177         ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) &            if (j <= inb(il)) then
178                                     *t(il,1)*b(il,1)*work(il)               fr(il, 1) = fr(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
179        else                    * (qent(il, j, 1) - rr(il, 1))
180         ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)               fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
181        endif                    * (uent(il, j, 1) - u(il, 1))
182                 fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
183        ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2) &                    * (vent(il, j, 1) - v(il, 1))
184        -t(il,1))*work(il)/cpn(il,1)            endif
185           enddo
186        if (cvflag_grav) then      enddo
187  !jyg1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)  
188  ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)      ! calculate tendencies of potential temperature and mixing ratio
189         fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &      ! at levels above the lowest level
190                  +sigd*0.5*(evap(il,1)+evap(il,2))  
191  !+tard     :          +sigd*evap(il,1)      ! first find the net saturated updraft and downdraft mass fluxes
192        ! through each level
193         fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)  
194        loop_i: do i = 2, nl - 1
195         fu(il,1)=fu(il,1)+0.01*grav*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.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &            ad(:ncum) = 0.
198                 +am(il)*(v(il,2)-v(il,1)))  
199        else  ! cvflag_grav            do k = i + 1, nl + 1
200         fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &               do il = 1, ncum
201                  +sigd*0.5*(evap(il,1)+evap(il,2))                  if (i <= inb(il) .and. k <= (inb(il) + 1)) then
202         fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)                     amp1(il) = amp1(il) + m(il, k)
203         fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &                  endif
204                 +am(il)*(u(il,2)-u(il,1)))               end do
205         fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &            end do
206                 +am(il)*(v(il,2)-v(il,1)))  
207        endif ! cvflag_grav            do k = 1, i
208                 do j = i + 1, nl + 1
209        enddo ! il                  do il = 1, ncum
210                       if (i <= inb(il) .and. j <= (inb(il) + 1)) then
211  !      do j=1,ntra                        amp1(il) = amp1(il) + ment(il, k, j)
212  !       do il=1,ncum                     endif
213  !        if (cvflag_grav) then                  end do
214  !         ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)               end do
215  !     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))            end do
216  !     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))  
217  !        else            do k = 1, i - 1
218  !         ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)               do j = i, nl + 1 ! newvecto: nl au lieu nl + 1?
219  !     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))                  do il = 1, ncum
220  !     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))                     if (i <= inb(il) .and. j <= inb(il)) then
221  !        endif                        ad(il) = ad(il) + ment(il, j, k)
222  !       enddo                     endif
223  !      enddo                  end do
224                 end do
225        do j=2,nl            end do
226         do il=1,ncum  
227          if (j.le.inb(il)) then            do il = 1, ncum
228           if (cvflag_grav) then               if (i <= inb(il)) then
229            fr(il,1)=fr(il,1) &                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
230               +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))                  cpinv = 1.0 / cpn(il, i)
231            fu(il,1)=fu(il,1) &  
232               +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))                  if (0.01 * rg * dpinv * amp1(il) >= delti) iflag(il) = 1
233            fv(il,1)=fv(il,1) &  
234               +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))                  ft(il, i) = 0.01 * rg * dpinv * (amp1(il) * (t(il, i + 1) &
235           else   ! cvflag_grav                       - t(il, i) + (gz(il, i + 1) - gz(il, i)) * cpinv) &
236            fr(il,1)=fr(il,1) &                       - ad(il) * (t(il, i) - t(il, i - 1) + (gz(il, i) &
237                 +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))                       - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) &
238            fu(il,1)=fu(il,1) &                       * (evap(il, i) + evap(il, i + 1)) - 0.009 * rg * sigd &
239                 +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))                       * (mp(il, i + 1) * t(il, i) * b(il, i) - mp(il, i) &
240            fv(il,1)=fv(il,1) &                       * t(il, i - 1) * cpn(il, i - 1) * cpinv * b(il, i - 1)) &
241                 +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))                       * dpinv + 0.01 * rg * dpinv * ment(il, i, i) &
242           endif  ! cvflag_grav                       * (hp(il, i) - h(il, i) + t(il, i) * (rcpv - rcpd) &
243          endif ! j                       * (rr(il, i) - qent(il, i, i))) * cpinv + 0.01 * sigd &
244         enddo                       * wt(il, i) * (rcw - rcpd) * water(il, i + 1) &
245        enddo                       * (t(il, i + 1) - t(il, i)) * dpinv * cpinv
246                    fr(il, i) = 0.01 * rg * dpinv * (amp1(il) * (rr(il, i + 1) &
247  !      do k=1,ntra                       - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
248  !       do j=2,nl                  fu(il, i) = fu(il, i) + 0.01 * rg * dpinv * (amp1(il) &
249  !        do il=1,ncum                       * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) &
250  !         if (j.le.inb(il)) then                       - u(il, i - 1)))
251                    fv(il, i) = fv(il, i) + 0.01 * rg * dpinv * (amp1(il) &
252  !          if (cvflag_grav) then                       * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) &
253  !           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)                       - v(il, i - 1)))
254  !     :                *(traent(il,j,1,k)-tra(il,1,k))               endif
255  !          else            end do
256  !           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)  
257  !     :                *(traent(il,j,1,k)-tra(il,1,k))            do k = 1, i - 1
258  !          endif               do il = 1, ncum
259                    if (i <= inb(il)) then
260  !         endif                     dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
261  !        enddo                     cpinv = 1.0 / cpn(il, i)
262  !       enddo  
263  !      enddo                     awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
264                       awat = amax1(awat, 0.0)
265  !  
266  !   ***  calculate tendencies of potential temperature and mixing ratio  ***                     fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
267  !   ***               at levels above the lowest level                   ***                          * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i))
268  !                     fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
269  !   ***  first find the net saturated updraft and downdraft mass fluxes  ***                          * ment(il, k, i) * (uent(il, k, i) - u(il, i))
270  !   ***                      through each level                          ***                     fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
271  !                          * ment(il, k, i) * (vent(il, k, i) - v(il, i))
272    
273        do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?                     ! (saturated updrafts resulting from mixing)
274                       qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat)
275         num1=0                     nqcond(il, i) = nqcond(il, i) + 1.
276         do il=1,ncum                  endif ! i
277          if(i.le.inb(il))num1=num1+1               end do
278         enddo            end do
279         if(num1.le.0)go to 500  
280              do k = i, nl + 1
281         call zilch(amp1,ncum)               do il = 1, ncum
282         call zilch(ad,ncum)                  if (i <= inb(il) .and. k <= inb(il)) then
283                       dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
284        do 440 k=i+1,nl+1                     cpinv = 1.0 / cpn(il, i)
285         do 441 il=1,ncum  
286          if (i.le.inb(il) .and. k.le.(inb(il)+1)) then                     fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
287           amp1(il)=amp1(il)+m(il,k)                          * ment(il, k, i) * (qent(il, k, i) - rr(il, i))
288          endif                     fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
289   441   continue                          * ment(il, k, i) * (uent(il, k, i) - u(il, i))
290   440  continue                     fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
291                            * ment(il, k, i) * (vent(il, k, i) - v(il, i))
292        do 450 k=1,i                  endif
293         do 451 j=i+1,nl+1               end do
294          do 452 il=1,ncum            end do
295           if (i.le.inb(il) .and. j.le.(inb(il)+1)) then  
296            amp1(il)=amp1(il)+ment(il,k,j)            do il = 1, ncum
297           endif               if (i <= inb(il)) then
298  452     continue                  dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
299  451    continue                  cpinv = 1.0 / cpn(il, i)
300  450   continue  
301                    ! sb: on ne fait pas encore la correction permettant de mieux
302        do 470 k=1,i-1                  ! conserver l'eau:
303         do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?                  fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) &
304          do 472 il=1,ncum                       + evap(il, i + 1)) + 0.01 * rg * (mp(il, i + 1) &
305          if (i.le.inb(il) .and. j.le.inb(il)) then                       * (qp(il, i + 1) - rr(il, i)) - mp(il, i) * (qp(il, i) &
306           ad(il)=ad(il)+ment(il,j,k)                       - rr(il, i - 1))) * dpinv
307          endif  
308  472     continue                  fu(il, i) = fu(il, i) + 0.01 * rg * (mp(il, i + 1) &
309  471    continue                       * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) &
310  470   continue                       - u(il, i - 1))) * dpinv
311                    fv(il, i) = fv(il, i) + 0.01 * rg * (mp(il, i + 1) &
312        do 1350 il=1,ncum                       * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) &
313        if (i.le.inb(il)) then                       - v(il, i - 1))) * dpinv
314         dpinv=1.0/(ph(il,i)-ph(il,i+1))               endif
315         cpinv=1.0/cpn(il,i)            end do
316    
317  ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4            ! sb: interface with the cloud parameterization:
318        if (cvflag_grav) then  
319         if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto            do k = i + 1, nl
320        else               do il = 1, ncum
321         if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto                  if (k <= inb(il) .and. i <= inb(il)) then
322        endif                     ! (saturated downdrafts resulting from mixing)
323                       qcond(il, i) = qcond(il, i) + elij(il, k, i)
324        if (cvflag_grav) then                     nqcond(il, i) = nqcond(il, i) + 1.
325         ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &                  endif
326            +(gz(il,i+1)-gz(il,i))*cpinv) &               enddo
327            -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &            enddo
328            -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))  
329         rat=cpn(il,i-1)*cpinv            ! (particular case: no detraining level is found)
330         ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &            do il = 1, ncum
331           -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv               if (i <= inb(il) .and. nent(il, i) == 0) then
332         ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &                  qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i)
333            +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv                  nqcond(il, i) = nqcond(il, i) + 1.
334        else  ! cvflag_grav               endif
335         ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &            enddo
336            +(gz(il,i+1)-gz(il,i))*cpinv) &  
337            -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &            do il = 1, ncum
338            -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))               if (i <= inb(il) .and. nqcond(il, i) /= 0.) then
339         rat=cpn(il,i-1)*cpinv                  qcond(il, i) = qcond(il, i) / nqcond(il, i)
340         ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &               endif
341           -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv            enddo
342         ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &         end if
343            +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv      end do loop_i
344        endif ! cvflag_grav  
345        ! move the detrainment at level inb down to level inb - 1
346        ! in such a way as to preserve the vertically
347        ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) &      ! integrated enthalpy and water tendencies
348                   *(t(il,i+1)-t(il,i))*dpinv*cpinv  
349        do il = 1, ncum
350        if (cvflag_grav) then         ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) &
351         fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &              - h(il, inb(il)) + t(il, inb(il)) * (rcpv - rcpd) &
352                   -ad(il)*(rr(il,i)-rr(il,i-1)))              * (rr(il, inb(il)) - qent(il, inb(il), inb(il)))) &
353         fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &              / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1)))
354                     -ad(il)*(u(il,i)-u(il,i-1)))         ft(il, inb(il)) = ft(il, inb(il)) - ax
355         fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &         ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) &
356                     -ad(il)*(v(il,i)-v(il,i-1)))              * (ph(il, inb(il)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) &
357        else  ! cvflag_grav              * (ph(il, inb(il) - 1) - ph(il, inb(il))))
358         fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &  
359                   -ad(il)*(rr(il,i)-rr(il,i-1)))         bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) &
360         fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &              - rr(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
361                     -ad(il)*(u(il,i)-u(il,i-1)))         fr(il, inb(il)) = fr(il, inb(il)) - bx
362         fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &         fr(il, inb(il) - 1) = fr(il, inb(il) - 1) &
363                     -ad(il)*(v(il,i)-v(il,i-1)))              + bx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
364        endif ! cvflag_grav              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
365    
366        endif ! i         cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) &
367  1350  continue              - u(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
368           fu(il, inb(il)) = fu(il, inb(il)) - cx
369  !      do k=1,ntra         fu(il, inb(il) - 1) = fu(il, inb(il) - 1) &
370  !       do il=1,ncum              + cx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
371  !        if (i.le.inb(il)) then              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
372  !         dpinv=1.0/(ph(il,i)-ph(il,i+1))  
373  !         cpinv=1.0/cpn(il,i)         dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) &
374  !         if (cvflag_grav) then              - v(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
375  !           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv         fv(il, inb(il)) = fv(il, inb(il)) - dx
376  !     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))         fv(il, inb(il) - 1) = fv(il, inb(il) - 1) &
377  !     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))              + dx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
378  !         else              / (ph(il, inb(il) - 1) - ph(il, inb(il)))
379  !           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv  
380  !     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))      end do
381  !     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))  
382  !         endif      ! homoginize tendencies below cloud base
383  !        endif  
384  !       enddo      do il = 1, ncum
385  !      enddo         asum(il) = 0.0
386           bsum(il) = 0.0
387        do 480 k=1,i-1         csum(il) = 0.0
388         do 1370 il=1,ncum         dsum(il) = 0.0
389          if (i.le.inb(il)) then      enddo
390           dpinv=1.0/(ph(il,i)-ph(il,i+1))  
391           cpinv=1.0/cpn(il,i)      do i = 1, nl
392           do il = 1, ncum
393        awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)            if (i <= (icb(il) - 1)) then
394        awat=amax1(awat,0.0)               asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
395                 bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (rcw - rcpd) &
396        if (cvflag_grav) then                    * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
397        fr(il,i)=fr(il,i) &               csum(il) = csum(il) + (lv(il, i) + (rcw - rcpd) * (t(il, i) &
398           +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))                    - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
399        fu(il,i)=fu(il,i) &               dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) &
400                 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                    / th(il, i)
401        fv(il,i)=fv(il,i) &            endif
402                 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))         enddo
403        else  ! cvflag_grav      enddo
404        fr(il,i)=fr(il,i) &  
405           +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))      do i = 1, nl
406        fu(il,i)=fu(il,i) &         do il = 1, ncum
407           +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))            if (i <= (icb(il) - 1)) then
408        fv(il,i)=fv(il,i) &               ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
409           +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))               fr(il, i) = bsum(il) / csum(il)
       endif ! cvflag_grav  
   
 ! (saturated updrafts resulting from mixing)        ! cld  
         qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld  
         nqcond(il,i)=nqcond(il,i)+1.                ! cld  
       endif ! i  
 1370  continue  
 480   continue  
   
 !      do j=1,ntra  
 !       do k=1,i-1  
 !        do il=1,ncum  
 !         if (i.le.inb(il)) then  
 !          dpinv=1.0/(ph(il,i)-ph(il,i+1))  
 !          cpinv=1.0/cpn(il,i)  
 !          if (cvflag_grav) then  
 !           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)  
 !     :        *(traent(il,k,i,j)-tra(il,i,j))  
 !          else  
 !           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)  
 !     :        *(traent(il,k,i,j)-tra(il,i,j))  
 !          endif  
 !         endif  
 !        enddo  
 !       enddo  
 !      enddo  
   
       do 490 k=i,nl+1  
        do 1380 il=1,ncum  
         if (i.le.inb(il) .and. k.le.inb(il)) then  
          dpinv=1.0/(ph(il,i)-ph(il,i+1))  
          cpinv=1.0/cpn(il,i)  
   
          if (cvflag_grav) then  
          fr(il,i)=fr(il,i) &  
                +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))  
          fu(il,i)=fu(il,i) &  
                +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))  
          fv(il,i)=fv(il,i) &  
                +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))  
          else  ! cvflag_grav  
          fr(il,i)=fr(il,i) &  
                +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))  
          fu(il,i)=fu(il,i) &  
                +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))  
          fv(il,i)=fv(il,i) &  
                +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))  
          endif ! cvflag_grav  
         endif ! i and k  
 1380   continue  
 490   continue  
   
 !      do j=1,ntra  
 !       do k=i,nl+1  
 !        do il=1,ncum  
 !         if (i.le.inb(il) .and. k.le.inb(il)) then  
 !          dpinv=1.0/(ph(il,i)-ph(il,i+1))  
 !          cpinv=1.0/cpn(il,i)  
 !          if (cvflag_grav) then  
 !           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)  
 !     :         *(traent(il,k,i,j)-tra(il,i,j))  
 !          else  
 !           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)  
 !     :             *(traent(il,k,i,j)-tra(il,i,j))  
 !          endif  
 !         endif ! i and k  
 !        enddo  
 !       enddo  
 !      enddo  
   
       do 1400 il=1,ncum  
        if (i.le.inb(il)) then  
         dpinv=1.0/(ph(il,i)-ph(il,i+1))  
         cpinv=1.0/cpn(il,i)  
   
         if (cvflag_grav) then  
 ! sb: on ne fait pas encore la correction permettant de mieux  
 ! conserver l'eau:  
          fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &  
               +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &  
                      *(rp(il,i)-rr(il,i-1)))*dpinv  
   
          fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &  
                    -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv  
          fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &  
                    -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv  
         else  ! cvflag_grav  
          fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &  
               +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &  
                      *(rp(il,i)-rr(il,i-1)))*dpinv  
          fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &  
                    -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv  
          fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &  
                    -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv  
         endif ! cvflag_grav  
   
       endif ! i  
 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  
   
 !      do j=1,ntra  
 !       do il=1,ncum  
 !        if (i.le.inb(il)) then  
 !         dpinv=1.0/(ph(il,i)-ph(il,i+1))  
 !         cpinv=1.0/cpn(il,i)  
   
 !         if (cvflag_grav) then  
 !          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv  
 !     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))  
 !     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))  
 !         else  
 !          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv  
 !     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))  
 !     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))  
 !         endif  
 !        endif ! i  
 !       enddo  
 !      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  
   
 !      do j=1,ntra  
 !       do il=1,ncum  
 !        ex=0.1*ment(il,inb(il),inb(il))  
 !     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))  
 !     :      /(ph(il,inb(il))-ph(il,inb(il)+1))  
 !        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex  
 !        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)  
 !     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))  
 !     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))  
 !       enddo  
 !      enddo  
   
 !  
 !   ***    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  !!!!      DO il=1,ncum      do i = 1, nl
428  !!!!      do i=icb(il),inb(il)         do il = 1, ncum
429  !!!!            if (i >= icb(il) .and. i <= inb(il)) then
430  !!!!      upwd(il,i)=0.0               upwd(il, i) = 0.0
431  !!!!      dnwd(il,i)=0.0               dnwd(il, i) = 0.0
432  !!!!      do k=i,inb(il)            endif
433  !!!!      up1=0.0         enddo
434  !!!!      dn1=0.0      enddo
435  !!!!      do n=1,i-1  
436  !!!!      up1=up1+ment(il,n,k)      do i = 1, nl
437  !!!!      dn1=dn1-ment(il,k,n)         do k = 1, nl
438  !!!!      enddo            do il = 1, ncum
439  !!!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1               up1(il, k, i) = 0.0
440  !!!!      dnwd(il,i)=dnwd(il,i)+dn1               dn1(il, k, i) = 0.0
441  !!!!      enddo            enddo
442  !!!!      enddo         enddo
443  !!!!      enddo
444  !!!!      ENDDO  
445        do i = 1, nl
446  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc         do k = i, nl
447  !        determination de la variation de flux ascendant entre            do n = 1, i - 1
448  !        deux niveau non dilue mike               do il = 1, ncum
449  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc                  if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then
450                       up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
451        do i=1,nl                     dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
452         do il=1,ncum                  endif
453          mike(il,i)=m(il,i)               enddo
454         enddo            enddo
455        enddo         enddo
456        enddo
457        do i=nl+1,nd  
458         do il=1,ncum      do i = 2, nl
459          mike(il,i)=0.         do k = i, nl
460         enddo            do il = 1, ncum
461        enddo               if (i <= inb(il).and.k <= inb(il)) then
462                    upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
463        do i=1,nd                  dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
464         do il=1,ncum               endif
465          ma(il,i)=0            enddo
466         enddo         enddo
467        enddo      enddo
468    
469        do i=1,nl      ! D\'etermination de la variation de flux ascendant entre
470         do j=i,nl      ! deux niveaux non dilu\'es mike
471          do il=1,ncum  
472           ma(il,i)=ma(il,i)+m(il,j)      do i = 1, nl
473          enddo         do il = 1, ncum
474         enddo            mike(il, i) = m(il, i)
475        enddo         enddo
476        enddo
477        do i=nl+1,nd  
478         do il=1,ncum      do i = nl + 1, klev
479          ma(il,i)=0.         do il = 1, ncum
480         enddo            mike(il, i) = 0.
481        enddo         enddo
482        enddo
483        do i=1,nl  
484         do il=1,ncum      do i = 1, klev
485          if (i.le.(icb(il)-1)) then         do il = 1, ncum
486           ma(il,i)=0            ma(il, i) = 0
487          endif         enddo
488         enddo      enddo
489        enddo  
490        do i = 1, nl
491  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc         do j = i, nl
492  !        icb represente de niveau ou se trouve la            do il = 1, ncum
493  !        base du nuage , et inb le top du nuage               ma(il, i) = ma(il, i) + m(il, j)
494  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc            enddo
495           enddo
496        do i=1,nd      enddo
497         do il=1,ncum  
498          mke(il,i)=upwd(il,i)+dnwd(il,i)      do i = nl + 1, klev
499         enddo         do il = 1, ncum
500        enddo            ma(il, i) = 0.
501           enddo
502        do i=1,nd      enddo
503         DO 999 il=1,ncum  
504          rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) &      do i = 1, nl
505                /(cpd*(1.-rr(il,i))+rr(il,i)*cpv)         do il = 1, ncum
506          tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp            if (i <= (icb(il) - 1)) then
507          tps(il,i)=tp(il,i)               ma(il, i) = 0
508  999    CONTINUE            endif
509        enddo         enddo
510        enddo
511  !  
512  !   *** diagnose the in-cloud mixing ratio   ***            ! cld      ! icb repr\'esente le niveau o\`u se trouve la base du nuage, et
513  !   ***           of condensed water         ***            ! cld      ! inb le sommet du nuage
514  !                                                           ! cld  
515        do i = 1, klev
516         do i=1,nd                                            ! cld         DO il = 1, ncum
517          do il=1,ncum                                        ! cld            rdcp = (rd * (1. - rr(il, i)) - rr(il, i) * rv) &
518           mac(il,i)=0.0                                      ! cld                 / (rcpd * (1. - rr(il, i)) + rr(il, i) * rcpv)
519           wa(il,i)=0.0                                       ! cld            tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp
520           siga(il,i)=0.0                                     ! cld            tps(il, i) = tp(il, i)
521           sax(il,i)=0.0                                      ! cld         end DO
522          enddo                                               ! cld      enddo
523         enddo                                                ! cld  
524        ! Diagnose the in-cloud mixing ratio of condensed water
525         do i=minorig, nl                                     ! cld  
526          do k=i+1,nl+1                                       ! cld      do i = 1, klev
527           do il=1,ncum                                       ! cld         do il = 1, ncum
528            if (i.le.inb(il) .and. k.le.(inb(il)+1)) then     ! cld            mac(il, i) = 0.0
529              mac(il,i)=mac(il,i)+m(il,k)                     ! cld            wa(il, i) = 0.0
530            endif                                             ! cld            siga(il, i) = 0.0
531           enddo                                              ! cld            sax(il, i) = 0.0
532          enddo                                               ! cld         enddo
533         enddo                                                ! cld      enddo
534    
535         do i=1,nl                                            ! cld      do i = minorig, nl
536          do j=1,i                                            ! cld         do k = i + 1, nl + 1
537           do il=1,ncum                                       ! cld            do il = 1, ncum
538            if (i.ge.icb(il) .and. i.le.(inb(il)-1) &               if (i <= inb(il) .and. k <= (inb(il) + 1)) then
539              .and. j.ge.icb(il) ) then                       ! cld                  mac(il, i) = mac(il, i) + m(il, k)
540             sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) &               endif
541                *(ph(il,j)-ph(il,j+1))/p(il,j)                ! cld            enddo
542            endif                                             ! cld         enddo
543           enddo                                              ! cld      enddo
544          enddo                                               ! cld  
545         enddo                                                ! cld      do i = 1, nl
546           do j = 1, i
547         do i=1,nl                                            ! cld            do il = 1, ncum
548          do il=1,ncum                                        ! cld               if (i >= icb(il) .and. i <= (inb(il) - 1) &
549           if (i.ge.icb(il) .and. i.le.(inb(il)-1) &                    .and. j >= icb(il)) then
550               .and. sax(il,i).gt.0.0 ) then                  ! cld                  sax(il, i) = sax(il, i) + rd * (tvp(il, j) - tv(il, j)) &
551             wa(il,i)=sqrt(2.*sax(il,i))                      ! cld                       * (ph(il, j) - ph(il, j + 1)) / p(il, j)
552           endif                                              ! cld               endif
553          enddo                                               ! cld            enddo
554         enddo                                                ! cld         enddo
555        enddo
556         do i=1,nl                                            ! cld  
557          do il=1,ncum                                        ! cld      do i = 1, nl
558           if (wa(il,i).gt.0.0) &         do il = 1, ncum
559             siga(il,i)=mac(il,i)/wa(il,i) &            if (i >= icb(il) .and. i <= (inb(il) - 1) &
560                 *rrd*tvp(il,i)/p(il,i)/100./delta            ! cld                 .and. sax(il, i) > 0.0) then
561            siga(il,i) = min(siga(il,i),1.0)                  ! cld               wa(il, i) = sqrt(2. * sax(il, i))
562  !IM cf. FH            endif
563           if (iflag_clw.eq.0) then         enddo
564            qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) &      enddo
565                   + (1.-siga(il,i))*qcond(il,i)              ! cld  
566           else if (iflag_clw.eq.1) then      do i = 1, nl
567            qcondc(il,i)=qcond(il,i)              ! cld         do il = 1, ncum
568           endif            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.47  
changed lines
  Added in v.254

  ViewVC Help
Powered by ViewVC 1.1.21