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