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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 69 - (hide annotations)
Mon Feb 18 16:33:12 2013 UTC (11 years, 3 months ago) by guez
Original Path: trunk/libf/phylmd/CV3_routines/cv3_yield.f90
File size: 23628 byte(s)
Deleted files cvparam3.f90 and nuagecom.f90. Moved variables from
module cvparam3 to module cv3_param_m. Moved variables rad_chau1 and
rad_chau2 from module nuagecom to module conf_phys_m.

Read clesphys2_nml from conf_phys instead of gcm.

Removed argument iflag_con from several procedures. Access module
variable instead.

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

  ViewVC Help
Powered by ViewVC 1.1.21