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

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

  ViewVC Help
Powered by ViewVC 1.1.21