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

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

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

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