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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 76 - (show annotations)
Fri Nov 15 18:45:49 2013 UTC (10 years, 6 months ago) by guez
Original Path: trunk/phylmd/CV3_routines/cv3_yield.f90
File size: 23628 byte(s)
Moved everything out of libf.
1
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 cv3_param_m
13 use cvthermo
14 use cvflag
15 implicit none
16
17
18 ! inputs:
19 integer, intent(in):: 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
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