/[lmdze]/trunk/phylmd/CV3_routines/cv3_yield.f
ViewVC logotype

Diff of /trunk/phylmd/CV3_routines/cv3_yield.f

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

revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC revision 97 by guez, Fri Apr 25 14:58:31 2014 UTC
# Line 1  Line 1 
1    module cv3_yield_m
2    
3        SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra  &    implicit none
4                            ,icb,inb,delt &  
5                            ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th &  contains
6                            ,ep,clw,m,tp,mp,rp,up,vp,trap &  
7                            ,wt,water,evap,b &    SUBROUTINE cv3_yield(nloc,ncum,nd,na  &
8                            ,ment,qent,uent,vent,nent,elij,traent,sig &         ,icb,inb,delt &
9                            ,tv,tvp &         ,t,rr,u,v,gz,p,ph,h,hp,lv,cpn,th &
10                            ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra &         ,ep,clw,m,tp,mp,rp,up,vp &
11                            ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)         ,wt,water,evap,b &
12        use conema3_m         ,ment,qent,uent,vent,nent,elij,sig &
13              use cv3_param_m         ,tv,tvp &
14              use cvthermo         ,iflag,precip,VPrecip,ft,fr,fu,fv &
15              use cvflag         ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
16        implicit none      use conema3_m
17        use cv3_param_m
18        use cvthermo
19  ! inputs:      use cvflag
20        integer, intent(in):: ncum,nd,na,ntra,nloc  
21        integer icb(nloc), inb(nloc)      ! inputs:
22        real, intent(in):: delt      integer, intent(in):: ncum,nd,na,nloc
23        real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)      integer icb(nloc), inb(nloc)
24        real tra(nloc,nd,ntra), sig(nloc,nd)      real, intent(in):: delt
25        real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)      real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
26        real th(nloc,na), p(nloc,nd), tp(nloc,na)      real sig(nloc,nd)
27        real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)      real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
28        real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)      real th(nloc,na), p(nloc,nd), tp(nloc,na)
29        real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)      real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
30        real water(nloc,na), evap(nloc,na), b(nloc,na)      real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
31        real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)      real vp(nloc,na), wt(nloc,nd)
32  !ym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)      real water(nloc,na), evap(nloc,na), b(nloc,na)
33        real vent(nloc,na,na), elij(nloc,na,na)      real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
34        integer nent(nloc,na)      !ym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
35        real traent(nloc,na,na,ntra)      real vent(nloc,na,na), elij(nloc,na,na)
36        real tv(nloc,nd), tvp(nloc,nd)      integer nent(nloc,na)
37        real tv(nloc,nd), tvp(nloc,nd)
38  ! input/output:  
39        integer iflag(nloc)      ! input/output:
40        integer iflag(nloc)
41  ! outputs:  
42        real precip(nloc)      ! outputs:
43        real VPrecip(nloc,nd+1)      real precip(nloc)
44        real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)      real VPrecip(nloc,nd+1)
45        real ftra(nloc,nd,ntra)      real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
46        real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)      real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
47        real dnwd0(nloc,nd), mike(nloc,nd)      real dnwd0(nloc,nd), mike(nloc,nd)
48        real tls(nloc,nd), tps(nloc,nd)      real tls(nloc,nd), tps(nloc,nd)
49        real qcondc(nloc,nd)                               ! cld      real qcondc(nloc,nd)                               ! cld
50        real wd(nloc)                                      ! gust      real wd(nloc)                                      ! gust
51    
52  ! local variables:      ! local variables:
53        integer i,k,il,n,j,num1      integer i,k,il,n,j,num1
54        real rat, awat, delti      real rat, awat, delti
55        real ax, bx, cx, dx, ex      real ax, bx, cx, dx, ex
56        real cpinv, rdcp, dpinv      real cpinv, rdcp, dpinv
57        real lvcp(nloc,na), mke(nloc,na)      real lvcp(nloc,na), mke(nloc,na)
58        real am(nloc), work(nloc), ad(nloc), amp1(nloc)      real am(nloc), work(nloc), ad(nloc), amp1(nloc)
59  !!!      real up1(nloc), dn1(nloc)  !!!      real up1(nloc), dn1(nloc)
60        real up1(nloc,nd,nd), dn1(nloc,nd,nd)      real up1(nloc,nd,nd), dn1(nloc,nd,nd)
61        real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)      real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
62        real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld      real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld
63        real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld      real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
64    
65    
66  !-------------------------------------------------------------      !-------------------------------------------------------------
67    
68  ! initialization:      ! initialization:
69    
70        delti = 1.0/delt      delti = 1.0/delt
71    
72        do il=1,ncum      do il=1,ncum
73         precip(il)=0.0         precip(il)=0.0
74         wd(il)=0.0     ! gust         wd(il)=0.0     ! gust
75         VPrecip(il,nd+1)=0.         VPrecip(il,nd+1)=0.
76        enddo      enddo
77    
78        do i=1,nd      do i=1,nd
79         do il=1,ncum         do il=1,ncum
80           VPrecip(il,i)=0.0            VPrecip(il,i)=0.0
81           ft(il,i)=0.0            ft(il,i)=0.0
82           fr(il,i)=0.0            fr(il,i)=0.0
83           fu(il,i)=0.0            fu(il,i)=0.0
84           fv(il,i)=0.0            fv(il,i)=0.0
85           qcondc(il,i)=0.0                                ! cld            qcondc(il,i)=0.0                                ! cld
86           qcond(il,i)=0.0                                 ! cld            qcond(il,i)=0.0                                 ! cld
87           nqcond(il,i)=0.0                                ! cld            nqcond(il,i)=0.0                                ! cld
88         enddo         enddo
89        enddo      enddo
90    
91    
92        do i=1,nl      do i=1,nl
93         do il=1,ncum         do il=1,ncum
94           lvcp(il,i)=lv(il,i)/cpn(il,i)            lvcp(il,i)=lv(il,i)/cpn(il,i)
95         enddo         enddo
96        enddo      enddo
97    
98    
99  !      !
100  !   ***  calculate surface precipitation in mm/day     ***      !   ***  calculate surface precipitation in mm/day     ***
101  !      !
102        do il=1,ncum      do il=1,ncum
103         if(ep(il,inb(il)).ge.0.0001)then         if(ep(il,inb(il)).ge.0.0001)then
104          if (cvflag_grav) then            if (cvflag_grav) then
105           precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)               precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)
106          else            else
107           precip(il)=wt(il,1)*sigd*water(il,1)*8640.               precip(il)=wt(il,1)*sigd*water(il,1)*8640.
108          endif            endif
109         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
116           do il=1,ncum         do il=1,ncum
117            if (k.le.inb(il)) then            if (k.le.inb(il)) then
118              if (cvflag_grav) then               if (cvflag_grav) then
119               VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav                  VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav
120              else               else
121               VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10.                  VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10.
122              endif               endif
123            endif            endif
          end do  
124         end do         end do
125  !      end do
126  !      !
127  !      !
128  !   ***  calculate tendencies of lowest level potential temperature  ***      !
129  !   ***                      and mixing ratio                        ***      !   ***  calculate tendencies of lowest level potential temperature  ***
130  !      !   ***                      and mixing ratio                        ***
131        do il=1,ncum      !
132        do il=1,ncum
133         work(il)=1.0/(ph(il,1)-ph(il,2))         work(il)=1.0/(ph(il,1)-ph(il,2))
134         am(il)=0.0         am(il)=0.0
135        enddo      enddo
136    
137        do k=2,nl      do k=2,nl
138         do il=1,ncum         do il=1,ncum
139          if (k.le.inb(il)) then            if (k.le.inb(il)) then
140           am(il)=am(il)+m(il,k)               am(il)=am(il)+m(il,k)
141          endif            endif
142         enddo         enddo
143        enddo      enddo
144    
145        do il=1,ncum      do il=1,ncum
146    
147  ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4         ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
148        if (cvflag_grav) then         if (cvflag_grav) then
149        if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect            if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect
150         ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1) &            ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1) &
151                    +(gz(il,2)-gz(il,1))/cpn(il,1))                 +(gz(il,2)-gz(il,1))/cpn(il,1))
152        else         else
153         if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect            if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect
154         ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) &            ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) &
155                    +(gz(il,2)-gz(il,1))/cpn(il,1))                 +(gz(il,2)-gz(il,1))/cpn(il,1))
156        endif         endif
157    
158        ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))         ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))
159    
160        if (cvflag_grav) then         if (cvflag_grav) then
161         ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) &            ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) &
162                                     *t(il,1)*b(il,1)*work(il)                 *t(il,1)*b(il,1)*work(il)
163        else         else
164         ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)            ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)
165        endif         endif
166    
167        ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2) &         ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2) &
168        -t(il,1))*work(il)/cpn(il,1)              -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)  
169    
170         fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)         if (cvflag_grav) then
171              !jyg1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
172              ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
173              fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &
174                   +sigd*0.5*(evap(il,1)+evap(il,2))
175              !+tard     :          +sigd*evap(il,1)
176    
177         fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &            fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
178    
179              fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &
180                 +am(il)*(u(il,2)-u(il,1)))                 +am(il)*(u(il,2)-u(il,1)))
181         fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &            fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &
182                 +am(il)*(v(il,2)-v(il,1)))                 +am(il)*(v(il,2)-v(il,1)))
183        else  ! cvflag_grav         else  ! cvflag_grav
184         fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &            fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &
185                  +sigd*0.5*(evap(il,1)+evap(il,2))                 +sigd*0.5*(evap(il,1)+evap(il,2))
186         fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)            fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
187         fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &            fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &
188                 +am(il)*(u(il,2)-u(il,1)))                 +am(il)*(u(il,2)-u(il,1)))
189         fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &            fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &
190                 +am(il)*(v(il,2)-v(il,1)))                 +am(il)*(v(il,2)-v(il,1)))
191        endif ! cvflag_grav         endif ! cvflag_grav
192    
193        enddo ! il      enddo ! il
194    
195        do j=2,nl      do j=2,nl
196         do il=1,ncum         do il=1,ncum
197          if (j.le.inb(il)) then            if (j.le.inb(il)) then
198           if (cvflag_grav) then               if (cvflag_grav) then
199            fr(il,1)=fr(il,1) &                  fr(il,1)=fr(il,1) &
200               +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))                       +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))
201            fu(il,1)=fu(il,1) &                  fu(il,1)=fu(il,1) &
202               +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))                       +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
203            fv(il,1)=fv(il,1) &                  fv(il,1)=fv(il,1) &
204               +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))                       +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
205           else   ! cvflag_grav               else   ! cvflag_grav
206            fr(il,1)=fr(il,1) &                  fr(il,1)=fr(il,1) &
207                 +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))                       +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))
208            fu(il,1)=fu(il,1) &                  fu(il,1)=fu(il,1) &
209                 +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))                       +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
210            fv(il,1)=fv(il,1) &                  fv(il,1)=fv(il,1) &
211                 +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))                       +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
212           endif  ! cvflag_grav               endif  ! cvflag_grav
213          endif ! j            endif ! j
214         enddo         enddo
215        enddo      enddo
216    
217  !      !
218  !   ***  calculate tendencies of potential temperature and mixing ratio  ***      !   ***  calculate tendencies of potential temperature and mixing ratio  ***
219  !   ***               at levels above the lowest level                   ***      !   ***               at levels above the lowest level                   ***
220  !      !
221  !   ***  first find the net saturated updraft and downdraft mass fluxes  ***      !   ***  first find the net saturated updraft and downdraft mass fluxes  ***
222  !   ***                      through each level                          ***      !   ***                      through each level                          ***
223  !      !
224    
225        do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?      do  i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
226    
227         num1=0         num1=0
228         do il=1,ncum         do il=1,ncum
229          if(i.le.inb(il))num1=num1+1            if(i.le.inb(il))num1=num1+1
230         enddo         enddo
231         if(num1.le.0)go to 500         if(num1.le.0) cycle
232    
233         call zilch(amp1,ncum)         call zilch(amp1,ncum)
234         call zilch(ad,ncum)         call zilch(ad,ncum)
235    
236        do 440 k=i+1,nl+1         do  k=i+1,nl+1
237         do 441 il=1,ncum            do  il=1,ncum
238          if (i.le.inb(il) .and. k.le.(inb(il)+1)) then               if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
239           amp1(il)=amp1(il)+m(il,k)                  amp1(il)=amp1(il)+m(il,k)
240          endif               endif
241   441   continue            end do
242   440  continue         end do
243    
244        do 450 k=1,i         do  k=1,i
245         do 451 j=i+1,nl+1            do  j=i+1,nl+1
246          do 452 il=1,ncum               do  il=1,ncum
247           if (i.le.inb(il) .and. j.le.(inb(il)+1)) then                  if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
248            amp1(il)=amp1(il)+ment(il,k,j)                     amp1(il)=amp1(il)+ment(il,k,j)
249           endif                  endif
250  452     continue               end do
251  451    continue            end do
252  450   continue         end do
253    
254        do 470 k=1,i-1         do  k=1,i-1
255         do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?            do  j=i,nl+1 ! newvecto: nl au lieu nl+1?
256          do 472 il=1,ncum               do  il=1,ncum
257          if (i.le.inb(il) .and. j.le.inb(il)) then                  if (i.le.inb(il) .and. j.le.inb(il)) then
258           ad(il)=ad(il)+ment(il,j,k)                     ad(il)=ad(il)+ment(il,j,k)
259          endif                  endif
260  472     continue               end do
261  471    continue            end do
262  470   continue         end do
263    
264        do 1350 il=1,ncum         do  il=1,ncum
265        if (i.le.inb(il)) then            if (i.le.inb(il)) then
266         dpinv=1.0/(ph(il,i)-ph(il,i+1))               dpinv=1.0/(ph(il,i)-ph(il,i+1))
267         cpinv=1.0/cpn(il,i)               cpinv=1.0/cpn(il,i)
268    
269  ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4               ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
270        if (cvflag_grav) then               if (cvflag_grav) then
271         if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto                  if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
272        else               else
273         if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto                  if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
274        endif               endif
275    
276        if (cvflag_grav) then               if (cvflag_grav) then
277         ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &                  ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &
278            +(gz(il,i+1)-gz(il,i))*cpinv) &                       +(gz(il,i+1)-gz(il,i))*cpinv) &
279            -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &                       -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &
280            -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))                       -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
281         rat=cpn(il,i-1)*cpinv                  rat=cpn(il,i-1)*cpinv
282         ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &                  ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
283           -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv                       -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
284         ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &                  ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &
285            +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv                       +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
286        else  ! cvflag_grav               else  ! cvflag_grav
287         ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &                  ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &
288            +(gz(il,i+1)-gz(il,i))*cpinv) &                       +(gz(il,i+1)-gz(il,i))*cpinv) &
289            -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &                       -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &
290            -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))                       -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
291         rat=cpn(il,i-1)*cpinv                  rat=cpn(il,i-1)*cpinv
292         ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &                  ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
293           -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv                       -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
294         ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &                  ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &
295            +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv                       +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
296        endif ! cvflag_grav               endif ! cvflag_grav
297    
298    
299        ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) &               ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) &
300                   *(t(il,i+1)-t(il,i))*dpinv*cpinv                    *(t(il,i+1)-t(il,i))*dpinv*cpinv
301    
302        if (cvflag_grav) then               if (cvflag_grav) then
303         fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &                  fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &
304                   -ad(il)*(rr(il,i)-rr(il,i-1)))                       -ad(il)*(rr(il,i)-rr(il,i-1)))
305         fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &                  fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &
306                     -ad(il)*(u(il,i)-u(il,i-1)))                       -ad(il)*(u(il,i)-u(il,i-1)))
307         fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &                  fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &
308                     -ad(il)*(v(il,i)-v(il,i-1)))                       -ad(il)*(v(il,i)-v(il,i-1)))
309        else  ! cvflag_grav               else  ! cvflag_grav
310         fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &                  fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &
311                   -ad(il)*(rr(il,i)-rr(il,i-1)))                       -ad(il)*(rr(il,i)-rr(il,i-1)))
312         fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &                  fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &
313                     -ad(il)*(u(il,i)-u(il,i-1)))                       -ad(il)*(u(il,i)-u(il,i-1)))
314         fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &                  fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &
315                     -ad(il)*(v(il,i)-v(il,i-1)))                       -ad(il)*(v(il,i)-v(il,i-1)))
316        endif ! cvflag_grav               endif ! cvflag_grav
317    
318        endif ! i            endif ! i
319  1350  continue         end do
320    
321        do 480 k=1,i-1         do  k=1,i-1
322         do 1370 il=1,ncum            do  il=1,ncum
323          if (i.le.inb(il)) then               if (i.le.inb(il)) then
324           dpinv=1.0/(ph(il,i)-ph(il,i+1))                  dpinv=1.0/(ph(il,i)-ph(il,i+1))
325           cpinv=1.0/cpn(il,i)                  cpinv=1.0/cpn(il,i)
326    
327        awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)                  awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
328        awat=amax1(awat,0.0)                  awat=amax1(awat,0.0)
329    
330        if (cvflag_grav) then                  if (cvflag_grav) then
331        fr(il,i)=fr(il,i) &                     fr(il,i)=fr(il,i) &
332           +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
333        fu(il,i)=fu(il,i) &                     fu(il,i)=fu(il,i) &
334                 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
335        fv(il,i)=fv(il,i) &                     fv(il,i)=fv(il,i) &
336                 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
337        else  ! cvflag_grav                  else  ! cvflag_grav
338        fr(il,i)=fr(il,i) &                     fr(il,i)=fr(il,i) &
339           +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))                          +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
340        fu(il,i)=fu(il,i) &                     fu(il,i)=fu(il,i) &
341           +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
342        fv(il,i)=fv(il,i) &                     fv(il,i)=fv(il,i) &
343           +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))                          +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
344        endif ! cvflag_grav                  endif ! cvflag_grav
345    
346  ! (saturated updrafts resulting from mixing)        ! cld                  ! (saturated updrafts resulting from mixing)        ! cld
347          qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld                  qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
348          nqcond(il,i)=nqcond(il,i)+1.                ! cld                  nqcond(il,i)=nqcond(il,i)+1.                ! cld
349        endif ! i               endif ! i
350  1370  continue            end do
351  480   continue         end do
352    
353        do 490 k=i,nl+1         do  k=i,nl+1
354         do 1380 il=1,ncum            do  il=1,ncum
355          if (i.le.inb(il) .and. k.le.inb(il)) then               if (i.le.inb(il) .and. k.le.inb(il)) then
356           dpinv=1.0/(ph(il,i)-ph(il,i+1))                  dpinv=1.0/(ph(il,i)-ph(il,i+1))
357           cpinv=1.0/cpn(il,i)                  cpinv=1.0/cpn(il,i)
358    
359           if (cvflag_grav) then                  if (cvflag_grav) then
360           fr(il,i)=fr(il,i) &                     fr(il,i)=fr(il,i) &
361                 +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
362           fu(il,i)=fu(il,i) &                     fu(il,i)=fu(il,i) &
363                 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
364           fv(il,i)=fv(il,i) &                     fv(il,i)=fv(il,i) &
365                 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))                          +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
366           else  ! cvflag_grav                  else  ! cvflag_grav
367           fr(il,i)=fr(il,i) &                     fr(il,i)=fr(il,i) &
368                 +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))                          +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
369           fu(il,i)=fu(il,i) &                     fu(il,i)=fu(il,i) &
370                 +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))                          +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
371           fv(il,i)=fv(il,i) &                     fv(il,i)=fv(il,i) &
372                 +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))                          +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
373           endif ! cvflag_grav                  endif ! cvflag_grav
374          endif ! i and k               endif ! i and k
375  1380   continue            end do
376  490   continue         end do
377    
378        do 1400 il=1,ncum         do  il=1,ncum
379         if (i.le.inb(il)) then            if (i.le.inb(il)) then
380          dpinv=1.0/(ph(il,i)-ph(il,i+1))               dpinv=1.0/(ph(il,i)-ph(il,i+1))
381          cpinv=1.0/cpn(il,i)               cpinv=1.0/cpn(il,i)
382    
383          if (cvflag_grav) then               if (cvflag_grav) then
384  ! sb: on ne fait pas encore la correction permettant de mieux                  ! sb: on ne fait pas encore la correction permettant de mieux
385  ! conserver l'eau:                  ! conserver l'eau:
386           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)+evap(il,i+1)) &
387                +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &                       +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &
388                       *(rp(il,i)-rr(il,i-1)))*dpinv                       *(rp(il,i)-rr(il,i-1)))*dpinv
389    
390           fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &                  fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &
391                     -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv                       -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
392           fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &                  fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &
393                     -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv                       -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
394          else  ! cvflag_grav               else  ! cvflag_grav
395           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)+evap(il,i+1)) &
396                +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &                       +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &
397                       *(rp(il,i)-rr(il,i-1)))*dpinv                       *(rp(il,i)-rr(il,i-1)))*dpinv
398           fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &                  fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &
399                     -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv                       -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
400           fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &                  fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &
401                     -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv                       -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
402          endif ! cvflag_grav               endif ! cvflag_grav
403    
404        endif ! i            endif ! i
405  1400  continue         end do
406    
407  ! sb: interface with the cloud parameterization:          ! cld         ! sb: interface with the cloud parameterization:          ! cld
408    
409        do k=i+1,nl         do k=i+1,nl
410         do il=1,ncum            do il=1,ncum
411          if (k.le.inb(il) .and. i.le.inb(il)) then         ! cld               if (k.le.inb(il) .and. i.le.inb(il)) then         ! cld
412  ! (saturated downdrafts resulting from mixing)            ! cld                  ! (saturated downdrafts resulting from mixing)            ! cld
413            qcond(il,i)=qcond(il,i)+elij(il,k,i)            ! cld                  qcond(il,i)=qcond(il,i)+elij(il,k,i)            ! cld
414            nqcond(il,i)=nqcond(il,i)+1.                    ! cld                  nqcond(il,i)=nqcond(il,i)+1.                    ! cld
415          endif                                             ! cld               endif                                             ! cld
416         enddo                                              ! cld            enddo                                              ! cld
417        enddo                                               ! cld         enddo                                               ! cld
418    
419  ! (particular case: no detraining level is found)         ! cld         ! (particular case: no detraining level is found)         ! cld
420        do il=1,ncum                                        ! cld         do il=1,ncum                                        ! cld
421         if (i.le.inb(il) .and. nent(il,i).eq.0) then       ! cld            if (i.le.inb(il) .and. nent(il,i).eq.0) then       ! cld
422            qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld               qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
423            nqcond(il,i)=nqcond(il,i)+1.                    ! cld               nqcond(il,i)=nqcond(il,i)+1.                    ! cld
424         endif                                              ! cld            endif                                              ! cld
425        enddo                                               ! cld         enddo                                               ! cld
426    
427        do il=1,ncum                                        ! cld         do il=1,ncum                                        ! cld
428         if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then    ! cld            if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then    ! cld
429            qcond(il,i)=qcond(il,i)/nqcond(il,i)            ! cld               qcond(il,i)=qcond(il,i)/nqcond(il,i)            ! cld
430         endif                                              ! cld            endif                                              ! cld
431        enddo         enddo
432    
433  500   continue      end do
434    
435    
436  !   ***   move the detrainment at level inb down to level inb-1   ***      !   ***   move the detrainment at level inb down to level inb-1   ***
437  !   ***        in such a way as to preserve the vertically        ***      !   ***        in such a way as to preserve the vertically        ***
438  !   ***          integrated enthalpy and water tendencies         ***      !   ***          integrated enthalpy and water tendencies         ***
439  !      !
440        do 503 il=1,ncum      do  il=1,ncum
441    
442        ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il)) &         ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il)) &
443         +t(il,inb(il))*(cpv-cpd) &              +t(il,inb(il))*(cpv-cpd) &
444         *(rr(il,inb(il))-qent(il,inb(il),inb(il)))) &              *(rr(il,inb(il))-qent(il,inb(il),inb(il)))) &
445          /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))              /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
446        ft(il,inb(il))=ft(il,inb(il))-ax         ft(il,inb(il))=ft(il,inb(il))-ax
447        ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) &         ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) &
448            *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) &              *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) &
449            *(ph(il,inb(il)-1)-ph(il,inb(il))))              *(ph(il,inb(il)-1)-ph(il,inb(il))))
450    
451        bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il)) &         bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il)) &
452            -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))              -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
453        fr(il,inb(il))=fr(il,inb(il))-bx         fr(il,inb(il))=fr(il,inb(il))-bx
454        fr(il,inb(il)-1)=fr(il,inb(il)-1) &         fr(il,inb(il)-1)=fr(il,inb(il)-1) &
455           +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) &              +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) &
456              /(ph(il,inb(il)-1)-ph(il,inb(il)))              /(ph(il,inb(il)-1)-ph(il,inb(il)))
457    
458        cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il)) &         cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il)) &
459               -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))              -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
460        fu(il,inb(il))=fu(il,inb(il))-cx         fu(il,inb(il))=fu(il,inb(il))-cx
461        fu(il,inb(il)-1)=fu(il,inb(il)-1) &         fu(il,inb(il)-1)=fu(il,inb(il)-1) &
462             +cx*(ph(il,inb(il))-ph(il,inb(il)+1)) &              +cx*(ph(il,inb(il))-ph(il,inb(il)+1)) &
463                /(ph(il,inb(il)-1)-ph(il,inb(il)))              /(ph(il,inb(il)-1)-ph(il,inb(il)))
464    
465        dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il)) &         dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il)) &
466              -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))              -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
467        fv(il,inb(il))=fv(il,inb(il))-dx         fv(il,inb(il))=fv(il,inb(il))-dx
468        fv(il,inb(il)-1)=fv(il,inb(il)-1) &         fv(il,inb(il)-1)=fv(il,inb(il)-1) &
469            +dx*(ph(il,inb(il))-ph(il,inb(il)+1)) &              +dx*(ph(il,inb(il))-ph(il,inb(il)+1)) &
470               /(ph(il,inb(il)-1)-ph(il,inb(il)))              /(ph(il,inb(il)-1)-ph(il,inb(il)))
471    
472  503   continue      end do
473    
474  !      !
475  !   ***    homoginize tendencies below cloud base    ***      !   ***    homoginize tendencies below cloud base    ***
476  !      !
477  !      !
478        do il=1,ncum      do il=1,ncum
479         asum(il)=0.0         asum(il)=0.0
480         bsum(il)=0.0         bsum(il)=0.0
481         csum(il)=0.0         csum(il)=0.0
482         dsum(il)=0.0         dsum(il)=0.0
483        enddo      enddo
484    
485        do i=1,nl      do i=1,nl
486         do il=1,ncum         do il=1,ncum
487          if (i.le.(icb(il)-1)) then            if (i.le.(icb(il)-1)) then
488        asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))               asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
489        bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &               bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &
490                          *(ph(il,i)-ph(il,i+1))                    *(ph(il,i)-ph(il,i+1))
491        csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &               csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &
492                              *(ph(il,i)-ph(il,i+1))                    *(ph(il,i)-ph(il,i+1))
493        dsum(il)=dsum(il)+t(il,i)*(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)
494          endif            endif
495         enddo         enddo
496        enddo      enddo
497    
498  !!!!      do 700 i=1,icb(il)-1  !!!!      do 700 i=1,icb(il)-1
499        do i=1,nl      do i=1,nl
500         do il=1,ncum         do il=1,ncum
501          if (i.le.(icb(il)-1)) then            if (i.le.(icb(il)-1)) then
502           ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il))               ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il))
503           fr(il,i)=bsum(il)/csum(il)               fr(il,i)=bsum(il)/csum(il)
504          endif            endif
505         enddo         enddo
506        enddo      enddo
507    
508  !      !
509  !   ***           reset counter and return           ***      !   ***           reset counter and return           ***
510  !      !
511        do il=1,ncum      do il=1,ncum
512         sig(il,nd)=2.0         sig(il,nd)=2.0
513        enddo      enddo
514    
515    
516        do i=1,nd      do i=1,nd
517         do il=1,ncum         do il=1,ncum
518          upwd(il,i)=0.0            upwd(il,i)=0.0
519          dnwd(il,i)=0.0            dnwd(il,i)=0.0
520         enddo         enddo
521        enddo      enddo
522    
523        do i=1,nl      do i=1,nl
524         do il=1,ncum         do il=1,ncum
525          dnwd0(il,i)=-mp(il,i)            dnwd0(il,i)=-mp(il,i)
526         enddo         enddo
527        enddo      enddo
528        do i=nl+1,nd      do i=nl+1,nd
529         do il=1,ncum         do il=1,ncum
530          dnwd0(il,i)=0.            dnwd0(il,i)=0.
531         enddo         enddo
532        enddo      enddo
533    
534    
535        do i=1,nl      do i=1,nl
536         do il=1,ncum         do il=1,ncum
537          if (i.ge.icb(il) .and. i.le.inb(il)) then            if (i.ge.icb(il) .and. i.le.inb(il)) then
538            upwd(il,i)=0.0               upwd(il,i)=0.0
539            dnwd(il,i)=0.0               dnwd(il,i)=0.0
540          endif            endif
541         enddo         enddo
542        enddo      enddo
543    
544        do i=1,nl      do i=1,nl
545         do k=1,nl         do k=1,nl
546          do il=1,ncum            do il=1,ncum
547            up1(il,k,i)=0.0               up1(il,k,i)=0.0
548            dn1(il,k,i)=0.0               dn1(il,k,i)=0.0
549          enddo            enddo
550         enddo         enddo
551        enddo      enddo
552    
553        do i=1,nl      do i=1,nl
554         do k=i,nl         do k=i,nl
555          do n=1,i-1            do n=1,i-1
556           do il=1,ncum               do il=1,ncum
557            if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then                  if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
558               up1(il,k,i)=up1(il,k,i)+ment(il,n,k)                     up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
559               dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)                     dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
560            endif                  endif
561           enddo               enddo
562          enddo            enddo
563         enddo         enddo
564        enddo      enddo
565    
566        do i=2,nl      do i=2,nl
567         do k=i,nl         do k=i,nl
568          do il=1,ncum            do il=1,ncum
569  !test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then               !test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
570           if (i.le.inb(il).and.k.le.inb(il)) then               if (i.le.inb(il).and.k.le.inb(il)) then
571              upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)                  upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
572              dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)                  dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
573           endif               endif
574          enddo            enddo
575         enddo         enddo
576        enddo      enddo
577    
578    
579  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
580  !        determination de la variation de flux ascendant entre      !        determination de la variation de flux ascendant entre
581  !        deux niveau non dilue mike      !        deux niveau non dilue mike
582  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
583    
584        do i=1,nl      do i=1,nl
585         do il=1,ncum         do il=1,ncum
586          mike(il,i)=m(il,i)            mike(il,i)=m(il,i)
587         enddo         enddo
588        enddo      enddo
589    
590        do i=nl+1,nd      do i=nl+1,nd
591         do il=1,ncum         do il=1,ncum
592          mike(il,i)=0.            mike(il,i)=0.
593         enddo         enddo
594        enddo      enddo
595    
596        do i=1,nd      do i=1,nd
597         do il=1,ncum         do il=1,ncum
598          ma(il,i)=0            ma(il,i)=0
599         enddo         enddo
600        enddo      enddo
601    
602        do i=1,nl      do i=1,nl
603         do j=i,nl         do j=i,nl
604          do il=1,ncum            do il=1,ncum
605           ma(il,i)=ma(il,i)+m(il,j)               ma(il,i)=ma(il,i)+m(il,j)
606          enddo            enddo
607         enddo         enddo
608        enddo      enddo
609    
610        do i=nl+1,nd      do i=nl+1,nd
611         do il=1,ncum         do il=1,ncum
612          ma(il,i)=0.            ma(il,i)=0.
613         enddo         enddo
614        enddo      enddo
615    
616        do i=1,nl      do i=1,nl
617         do il=1,ncum         do il=1,ncum
618          if (i.le.(icb(il)-1)) then            if (i.le.(icb(il)-1)) then
619           ma(il,i)=0               ma(il,i)=0
620          endif            endif
621         enddo         enddo
622        enddo      enddo
623    
624  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
625  !        icb represente de niveau ou se trouve la      !        icb represente de niveau ou se trouve la
626  !        base du nuage , et inb le top du nuage      !        base du nuage , et inb le top du nuage
627  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
628    
629        do i=1,nd      do i=1,nd
630         do il=1,ncum         do il=1,ncum
631          mke(il,i)=upwd(il,i)+dnwd(il,i)            mke(il,i)=upwd(il,i)+dnwd(il,i)
632         enddo         enddo
633        enddo      enddo
634    
635        do i=1,nd      do i=1,nd
636         DO 999 il=1,ncum         DO  il=1,ncum
637          rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) &            rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) &
638                /(cpd*(1.-rr(il,i))+rr(il,i)*cpv)                 /(cpd*(1.-rr(il,i))+rr(il,i)*cpv)
639          tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp            tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp
640          tps(il,i)=tp(il,i)            tps(il,i)=tp(il,i)
641  999    CONTINUE         end DO
642        enddo      enddo
643    
644  !      !
645  !   *** diagnose the in-cloud mixing ratio   ***            ! cld      !   *** diagnose the in-cloud mixing ratio   ***            ! cld
646  !   ***           of condensed water         ***            ! cld      !   ***           of condensed water         ***            ! cld
647  !                                                           ! cld      !                                                           ! cld
648    
649         do i=1,nd                                            ! cld      do i=1,nd                                            ! cld
650          do il=1,ncum                                        ! cld         do il=1,ncum                                        ! cld
651           mac(il,i)=0.0                                      ! cld            mac(il,i)=0.0                                      ! cld
652           wa(il,i)=0.0                                       ! cld            wa(il,i)=0.0                                       ! cld
653           siga(il,i)=0.0                                     ! cld            siga(il,i)=0.0                                     ! cld
654           sax(il,i)=0.0                                      ! cld            sax(il,i)=0.0                                      ! cld
655          enddo                                               ! cld         enddo                                               ! cld
656         enddo                                                ! cld      enddo                                                ! cld
657    
658         do i=minorig, nl                                     ! cld      do i=minorig, nl                                     ! cld
659          do k=i+1,nl+1                                       ! cld         do k=i+1,nl+1                                       ! cld
660           do il=1,ncum                                       ! cld            do il=1,ncum                                       ! cld
661            if (i.le.inb(il) .and. k.le.(inb(il)+1)) then     ! cld               if (i.le.inb(il) .and. k.le.(inb(il)+1)) then     ! cld
662              mac(il,i)=mac(il,i)+m(il,k)                     ! cld                  mac(il,i)=mac(il,i)+m(il,k)                     ! cld
663            endif                                             ! cld               endif                                             ! cld
664           enddo                                              ! cld            enddo                                              ! cld
665          enddo                                               ! cld         enddo                                               ! cld
666         enddo                                                ! cld      enddo                                                ! cld
667    
668         do i=1,nl                                            ! cld      do i=1,nl                                            ! cld
669          do j=1,i                                            ! cld         do j=1,i                                            ! cld
670           do il=1,ncum                                       ! cld            do il=1,ncum                                       ! cld
671                 if (i.ge.icb(il) .and. i.le.(inb(il)-1) &
672                      .and. j.ge.icb(il) ) then                       ! cld
673                    sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) &
674                         *(ph(il,j)-ph(il,j+1))/p(il,j)                ! cld
675                 endif                                             ! cld
676              enddo                                              ! cld
677           enddo                                               ! cld
678        enddo                                                ! cld
679    
680        do i=1,nl                                            ! cld
681           do il=1,ncum                                        ! cld
682            if (i.ge.icb(il) .and. i.le.(inb(il)-1) &            if (i.ge.icb(il) .and. i.le.(inb(il)-1) &
683              .and. j.ge.icb(il) ) then                       ! cld                 .and. sax(il,i).gt.0.0 ) then                  ! cld
684             sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) &               wa(il,i)=sqrt(2.*sax(il,i))                      ! cld
685                *(ph(il,j)-ph(il,j+1))/p(il,j)                ! cld            endif                                              ! cld
686            endif                                             ! cld         enddo                                               ! cld
687           enddo                                              ! cld      enddo                                                ! cld
688          enddo                                               ! cld  
689         enddo                                                ! cld      do i=1,nl                                            ! cld
690           do il=1,ncum                                        ! cld
691         do i=1,nl                                            ! cld            if (wa(il,i).gt.0.0) &
692          do il=1,ncum                                        ! cld                 siga(il,i)=mac(il,i)/wa(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) &  
693                 *rrd*tvp(il,i)/p(il,i)/100./delta            ! cld                 *rrd*tvp(il,i)/p(il,i)/100./delta            ! cld
694            siga(il,i) = min(siga(il,i),1.0)                  ! cld            siga(il,i) = min(siga(il,i),1.0)                  ! cld
695  !IM cf. FH            !IM cf. FH
696           if (iflag_clw.eq.0) then            if (iflag_clw.eq.0) then
697            qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) &               qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) &
698                   + (1.-siga(il,i))*qcond(il,i)              ! cld                    + (1.-siga(il,i))*qcond(il,i)              ! cld
699           else if (iflag_clw.eq.1) then            else if (iflag_clw.eq.1) then
700            qcondc(il,i)=qcond(il,i)              ! cld               qcondc(il,i)=qcond(il,i)              ! cld
701           endif            endif
702    
703           enddo                                               ! cld
704        enddo                                                ! cld
705    
706          enddo                                               ! cld    end SUBROUTINE cv3_yield
        enddo                                                ! cld  
707    
708          return  end module cv3_yield_m
         end  

Legend:
Removed from v.82  
changed lines
  Added in v.97

  ViewVC Help
Powered by ViewVC 1.1.21