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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (hide annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/CV3_routines/cv3_yield.f90
File size: 27971 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

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

  ViewVC Help
Powered by ViewVC 1.1.21