/[lmdze]/trunk/Sources/phylmd/Thermcell/thermcell.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/Thermcell/thermcell.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
File size: 18477 byte(s)
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.

1 module thermcell_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE thermcell(ngrid, nlay, ptimestep, pplay, pplev, pphi, pu, pv, pt, &
8 po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0, r_aspect, l_mix, w2di, &
9 tho)
10
11 ! Calcul du transport vertical dans la couche limite en pr\'esence
12 ! de "thermiques" explicitement repr\'esent\'es. R\'ecriture \`a partir
13 ! d'un listing papier \`a Habas, le 14/02/00. Le thermique est
14 ! suppos\'e homog\`ene et dissip\'e par m\'elange avec son
15 ! environnement. La longueur "l_mix" contr\^ole l'efficacit\'e du
16 ! m\'elange. Le calcul du transport des diff\'erentes esp\`eces se fait
17 ! en prenant en compte :
18 ! 1. un flux de masse montant
19 ! 2. un flux de masse descendant
20 ! 3. un entra\^inement
21 ! 4. un d\'etra\^inement
22
23 USE dimphy, ONLY : klev, klon
24 USE suphec_m, ONLY : rd, rg, rkappa
25
26 ! arguments:
27
28 INTEGER ngrid, nlay, w2di
29 real tho
30 real ptimestep, l_mix, r_aspect
31 REAL, intent(in):: pt(ngrid, nlay)
32 real pdtadj(ngrid, nlay)
33 REAL, intent(in):: pu(ngrid, nlay)
34 real pduadj(ngrid, nlay)
35 REAL, intent(in):: pv(ngrid, nlay)
36 real pdvadj(ngrid, nlay)
37 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
38 REAL, intent(in):: pplay(ngrid, nlay)
39 real, intent(in):: pplev(ngrid, nlay+1)
40 real, intent(in):: pphi(ngrid, nlay)
41
42 integer idetr
43 save idetr
44 data idetr/3/
45
46 ! local:
47
48 INTEGER ig, k, l, lmaxa(klon), lmix(klon)
49 ! CR: on remplace lmax(klon, klev+1)
50 INTEGER lmax(klon), lmin(klon), lentr(klon)
51 real linter(klon)
52 real zmix(klon), fracazmix(klon)
53
54 real zmax(klon), zw, zw2(klon, klev+1), ztva(klon, klev)
55
56 real zlev(klon, klev+1), zlay(klon, klev)
57 REAL zh(klon, klev), zdhadj(klon, klev)
58 REAL ztv(klon, klev)
59 real zu(klon, klev), zv(klon, klev), zo(klon, klev)
60 real zla(klon, klev+1)
61 real zwa(klon, klev+1)
62 real zld(klon, klev+1)
63 real zva(klon, klev)
64 real zua(klon, klev)
65 real zoa(klon, klev)
66
67 real zha(klon, klev)
68 real wa_moy(klon, klev+1)
69 real fraca(klon, klev+1)
70 real fracc(klon, klev+1)
71 real zf, zf2
72 real thetath2(klon, klev), wth2(klon, klev)
73 common/comtherm/thetath2, wth2
74
75 integer isplit, nsplit
76 parameter (nsplit=10)
77 data isplit/0/
78 save isplit
79
80 logical sorties
81 real rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
82 real zpspsk(klon, klev)
83
84 real wmax(klon), wmaxa(klon)
85 real wa(klon, klev, klev+1)
86 real wd(klon, klev+1)
87 real fracd(klon, klev+1)
88 real xxx(klon, klev+1)
89 real larg_cons(klon, klev+1)
90 real larg_detr(klon, klev+1)
91 real fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
92 real fm(klon, klev+1), entr(klon, klev)
93 real fmc(klon, klev+1)
94
95 !CR:nouvelles variables
96 real f_star(klon, klev+1), entr_star(klon, klev)
97 real entr_star_tot(klon), entr_star2(klon)
98 real f(klon)
99 real zlevinter(klon)
100
101 EXTERNAL SCOPY
102
103 !-----------------------------------------------------------------------
104
105 ! initialisation:
106
107 sorties=.true.
108 IF(ngrid.NE.klon) THEN
109 PRINT *
110 PRINT *, 'STOP dans convadj'
111 PRINT *, 'ngrid =', ngrid
112 PRINT *, 'klon =', klon
113 ENDIF
114
115 ! incrementation eventuelle de tendances precedentes:
116
117 print *, '0 OK convect8'
118
119 DO l=1, nlay
120 DO ig=1, ngrid
121 zpspsk(ig, l)=(pplay(ig, l)/pplev(ig, 1))**RKAPPA
122 zh(ig, l)=pt(ig, l)/zpspsk(ig, l)
123 zu(ig, l)=pu(ig, l)
124 zv(ig, l)=pv(ig, l)
125 zo(ig, l)=po(ig, l)
126 ztv(ig, l)=zh(ig, l)*(1.+0.61*zo(ig, l))
127 end DO
128 end DO
129
130 print *, '1 OK convect8'
131
132 ! See notes, "thermcell.txt"
133 ! Calcul des altitudes des couches
134
135 do l=2, nlay
136 do ig=1, ngrid
137 zlev(ig, l)=0.5*(pphi(ig, l)+pphi(ig, l-1))/RG
138 enddo
139 enddo
140 do ig=1, ngrid
141 zlev(ig, 1)=0.
142 zlev(ig, nlay+1)=(2.*pphi(ig, klev)-pphi(ig, klev-1))/RG
143 enddo
144 do l=1, nlay
145 do ig=1, ngrid
146 zlay(ig, l)=pphi(ig, l)/RG
147 enddo
148 enddo
149
150 ! Calcul des densites
151
152 do l=1, nlay
153 do ig=1, ngrid
154 rho(ig, l)=pplay(ig, l)/(zpspsk(ig, l)*RD*zh(ig, l))
155 enddo
156 enddo
157
158 do l=2, nlay
159 do ig=1, ngrid
160 rhobarz(ig, l)=0.5*(rho(ig, l)+rho(ig, l-1))
161 enddo
162 enddo
163
164 do k=1, nlay
165 do l=1, nlay+1
166 do ig=1, ngrid
167 wa(ig, k, l)=0.
168 enddo
169 enddo
170 enddo
171
172 ! Calcul de w2, quarre de w a partir de la cape
173 ! a partir de w2, on calcule wa, vitesse de l'ascendance
174
175 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
176 ! w2 est stoke dans wa
177
178 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
179 ! independants par couches que pour calculer l'entrainement
180 ! a la base et la hauteur max de l'ascendance.
181
182 ! Indicages:
183 ! l'ascendance provenant du niveau k traverse l'interface l avec
184 ! une vitesse wa(k, l).
185 ! See notes, "thermcell.txt".
186
187 !CR: ponderation entrainement des couches instables
188 !def des entr_star tels que entr=f*entr_star
189 do l=1, klev
190 do ig=1, ngrid
191 entr_star(ig, l)=0.
192 enddo
193 enddo
194 ! determination de la longueur de la couche d entrainement
195 do ig=1, ngrid
196 lentr(ig)=1
197 enddo
198
199 !on ne considere que les premieres couches instables
200 do k=nlay-2, 1, -1
201 do ig=1, ngrid
202 if (ztv(ig, k).gt.ztv(ig, k+1).and. &
203 ztv(ig, k+1).le.ztv(ig, k+2)) then
204 lentr(ig)=k
205 endif
206 enddo
207 enddo
208
209 ! determination du lmin: couche d ou provient le thermique
210 do ig=1, ngrid
211 lmin(ig)=1
212 enddo
213 do ig=1, ngrid
214 do l=nlay, 2, -1
215 if (ztv(ig, l-1).gt.ztv(ig, l)) then
216 lmin(ig)=l-1
217 endif
218 enddo
219 enddo
220
221 ! definition de l'entrainement des couches
222 do l=1, klev-1
223 do ig=1, ngrid
224 if (ztv(ig, l).gt.ztv(ig, l+1).and. &
225 l.ge.lmin(ig).and.l.le.lentr(ig)) then
226 entr_star(ig, l)=(ztv(ig, l)-ztv(ig, l+1))* &
227 (zlev(ig, l+1)-zlev(ig, l))
228 endif
229 enddo
230 enddo
231 ! pas de thermique si couches 1->5 stables
232 do ig=1, ngrid
233 if (lmin(ig).gt.5) then
234 do l=1, klev
235 entr_star(ig, l)=0.
236 enddo
237 endif
238 enddo
239 ! calcul de l entrainement total
240 do ig=1, ngrid
241 entr_star_tot(ig)=0.
242 enddo
243 do ig=1, ngrid
244 do k=1, klev
245 entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig, k)
246 enddo
247 enddo
248
249 print *, 'fin calcul entr_star'
250 do k=1, klev
251 do ig=1, ngrid
252 ztva(ig, k)=ztv(ig, k)
253 enddo
254 enddo
255
256 do k=1, klev+1
257 do ig=1, ngrid
258 zw2(ig, k)=0.
259 fmc(ig, k)=0.
260
261 f_star(ig, k)=0.
262
263 larg_cons(ig, k)=0.
264 larg_detr(ig, k)=0.
265 wa_moy(ig, k)=0.
266 enddo
267 enddo
268
269 do ig=1, ngrid
270 linter(ig)=1.
271 lmaxa(ig)=1
272 lmix(ig)=1
273 wmaxa(ig)=0.
274 enddo
275
276 do l=1, nlay-2
277 do ig=1, ngrid
278 if (ztv(ig, l).gt.ztv(ig, l+1) &
279 .and.entr_star(ig, l).gt.1.e-10 &
280 .and.zw2(ig, l).lt.1e-10) then
281 f_star(ig, l+1)=entr_star(ig, l)
282 !test:calcul de dteta
283 zw2(ig, l+1)=2.*RG*(ztv(ig, l)-ztv(ig, l+1))/ztv(ig, l+1) &
284 *(zlev(ig, l+1)-zlev(ig, l)) &
285 *0.4*pphi(ig, l)/(pphi(ig, l+1)-pphi(ig, l))
286 larg_detr(ig, l)=0.
287 else if ((zw2(ig, l).ge.1e-10).and. &
288 (f_star(ig, l)+entr_star(ig, l).gt.1.e-10)) then
289 f_star(ig, l+1)=f_star(ig, l)+entr_star(ig, l)
290 ztva(ig, l)=(f_star(ig, l)*ztva(ig, l-1)+entr_star(ig, l) &
291 *ztv(ig, l))/f_star(ig, l+1)
292 zw2(ig, l+1)=zw2(ig, l)*(f_star(ig, l)/f_star(ig, l+1))**2+ &
293 2.*RG*(ztva(ig, l)-ztv(ig, l))/ztv(ig, l) &
294 *(zlev(ig, l+1)-zlev(ig, l))
295 endif
296 ! determination de zmax continu par interpolation lineaire
297 if (zw2(ig, l+1).lt.0.) then
298 if (abs(zw2(ig, l+1)-zw2(ig, l)).lt.1e-10) then
299 print *, 'pb linter'
300 endif
301 linter(ig)=(l*(zw2(ig, l+1)-zw2(ig, l)) &
302 -zw2(ig, l))/(zw2(ig, l+1)-zw2(ig, l))
303 zw2(ig, l+1)=0.
304 lmaxa(ig)=l
305 else
306 if (zw2(ig, l+1).lt.0.) then
307 print *, 'pb1 zw2<0'
308 endif
309 wa_moy(ig, l+1)=sqrt(zw2(ig, l+1))
310 endif
311 if (wa_moy(ig, l+1).gt.wmaxa(ig)) then
312 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
313 lmix(ig)=l+1
314 wmaxa(ig)=wa_moy(ig, l+1)
315 endif
316 enddo
317 enddo
318 print *, 'fin calcul zw2'
319
320 ! Calcul de la couche correspondant a la hauteur du thermique
321 do ig=1, ngrid
322 lmax(ig)=lentr(ig)
323 enddo
324 do ig=1, ngrid
325 do l=nlay, lentr(ig)+1, -1
326 if (zw2(ig, l).le.1.e-10) then
327 lmax(ig)=l-1
328 endif
329 enddo
330 enddo
331 ! pas de thermique si couches 1->5 stables
332 do ig=1, ngrid
333 if (lmin(ig).gt.5) then
334 lmax(ig)=1
335 lmin(ig)=1
336 endif
337 enddo
338
339 ! Determination de zw2 max
340 do ig=1, ngrid
341 wmax(ig)=0.
342 enddo
343
344 do l=1, nlay
345 do ig=1, ngrid
346 if (l.le.lmax(ig)) then
347 if (zw2(ig, l).lt.0.)then
348 print *, 'pb2 zw2<0'
349 endif
350 zw2(ig, l)=sqrt(zw2(ig, l))
351 wmax(ig)=max(wmax(ig), zw2(ig, l))
352 else
353 zw2(ig, l)=0.
354 endif
355 enddo
356 enddo
357
358 ! Longueur caracteristique correspondant a la hauteur des thermiques.
359 do ig=1, ngrid
360 zmax(ig)=0.
361 zlevinter(ig)=zlev(ig, 1)
362 enddo
363 do ig=1, ngrid
364 ! calcul de zlevinter
365 zlevinter(ig)=(zlev(ig, lmax(ig)+1)-zlev(ig, lmax(ig)))* &
366 linter(ig)+zlev(ig, lmax(ig))-lmax(ig)*(zlev(ig, lmax(ig)+1) &
367 -zlev(ig, lmax(ig)))
368 zmax(ig)=max(zmax(ig), zlevinter(ig)-zlev(ig, lmin(ig)))
369 enddo
370
371 print *, 'avant fermeture'
372 ! Fermeture, determination de f
373 do ig=1, ngrid
374 entr_star2(ig)=0.
375 enddo
376 do ig=1, ngrid
377 if (entr_star_tot(ig).LT.1.e-10) then
378 f(ig)=0.
379 else
380 do k=lmin(ig), lentr(ig)
381 entr_star2(ig)=entr_star2(ig)+entr_star(ig, k)**2 &
382 /(rho(ig, k)*(zlev(ig, k+1)-zlev(ig, k)))
383 enddo
384 ! Nouvelle fermeture
385 f(ig)=wmax(ig)/(max(500., zmax(ig))*r_aspect &
386 *entr_star2(ig))*entr_star_tot(ig)
387 endif
388 enddo
389 print *, 'apres fermeture'
390
391 ! Calcul de l'entrainement
392 do k=1, klev
393 do ig=1, ngrid
394 entr(ig, k)=f(ig)*entr_star(ig, k)
395 enddo
396 enddo
397 ! Calcul des flux
398 do ig=1, ngrid
399 do l=1, lmax(ig)-1
400 fmc(ig, l+1)=fmc(ig, l)+entr(ig, l)
401 enddo
402 enddo
403
404 ! determination de l'indice du debut de la mixed layer ou w decroit
405
406 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
407 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
408 ! d'une couche est \'egale \`a la hauteur de la couche alimentante.
409 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
410 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
411
412 do l=2, nlay
413 do ig=1, ngrid
414 if (l.le.lmaxa(ig)) then
415 zw=max(wa_moy(ig, l), 1.e-10)
416 larg_cons(ig, l)=zmax(ig)*r_aspect &
417 *fmc(ig, l)/(rhobarz(ig, l)*zw)
418 endif
419 enddo
420 enddo
421
422 do l=2, nlay
423 do ig=1, ngrid
424 if (l.le.lmaxa(ig)) then
425 if ((l_mix*zlev(ig, l)).lt.0.)then
426 print *, 'pb l_mix*zlev<0'
427 endif
428 larg_detr(ig, l)=sqrt(l_mix*zlev(ig, l))
429 endif
430 enddo
431 enddo
432
433 ! calcul de la fraction de la maille concern\'ee par l'ascendance en tenant
434 ! compte de l'epluchage du thermique.
435
436 !CR def de zmix continu (profil parabolique des vitesses)
437 do ig=1, ngrid
438 if (lmix(ig).gt.1.) then
439 if (((zw2(ig, lmix(ig)-1)-zw2(ig, lmix(ig))) &
440 *((zlev(ig, lmix(ig)))-(zlev(ig, lmix(ig)+1))) &
441 -(zw2(ig, lmix(ig))-zw2(ig, lmix(ig)+1)) &
442 *((zlev(ig, lmix(ig)-1))-(zlev(ig, lmix(ig))))).gt.1e-10) &
443 then
444
445 zmix(ig)=((zw2(ig, lmix(ig)-1)-zw2(ig, lmix(ig))) &
446 *((zlev(ig, lmix(ig)))**2-(zlev(ig, lmix(ig)+1))**2) &
447 -(zw2(ig, lmix(ig))-zw2(ig, lmix(ig)+1)) &
448 *((zlev(ig, lmix(ig)-1))**2-(zlev(ig, lmix(ig)))**2)) &
449 /(2.*((zw2(ig, lmix(ig)-1)-zw2(ig, lmix(ig))) &
450 *((zlev(ig, lmix(ig)))-(zlev(ig, lmix(ig)+1))) &
451 -(zw2(ig, lmix(ig))-zw2(ig, lmix(ig)+1)) &
452 *((zlev(ig, lmix(ig)-1))-(zlev(ig, lmix(ig))))))
453 else
454 zmix(ig)=zlev(ig, lmix(ig))
455 print *, 'pb zmix'
456 endif
457 else
458 zmix(ig)=0.
459 endif
460
461 if ((zmax(ig)-zmix(ig)).lt.0.) then
462 zmix(ig)=0.99*zmax(ig)
463 endif
464 enddo
465
466 ! calcul du nouveau lmix correspondant
467 do ig=1, ngrid
468 do l=1, klev
469 if (zmix(ig).ge.zlev(ig, l).and. &
470 zmix(ig).lt.zlev(ig, l+1)) then
471 lmix(ig)=l
472 endif
473 enddo
474 enddo
475
476 do l=2, nlay
477 do ig=1, ngrid
478 if(larg_cons(ig, l).gt.1.) then
479 fraca(ig, l)=(larg_cons(ig, l)-larg_detr(ig, l)) &
480 /(r_aspect*zmax(ig))
481 fraca(ig, l)=max(fraca(ig, l), 0.)
482 fraca(ig, l)=min(fraca(ig, l), 0.5)
483 fracd(ig, l)=1.-fraca(ig, l)
484 fracc(ig, l)=larg_cons(ig, l)/(r_aspect*zmax(ig))
485 else
486 fraca(ig, l)=0.
487 fracc(ig, l)=0.
488 fracd(ig, l)=1.
489 endif
490 enddo
491 enddo
492 !CR: calcul de fracazmix
493 do ig=1, ngrid
494 fracazmix(ig)=(fraca(ig, lmix(ig)+1)-fraca(ig, lmix(ig)))/ &
495 (zlev(ig, lmix(ig)+1)-zlev(ig, lmix(ig)))*zmix(ig) &
496 +fraca(ig, lmix(ig))-zlev(ig, lmix(ig))*(fraca(ig, lmix(ig)+1) &
497 -fraca(ig, lmix(ig)))/(zlev(ig, lmix(ig)+1)-zlev(ig, lmix(ig)))
498 enddo
499
500 do l=2, nlay
501 do ig=1, ngrid
502 if(larg_cons(ig, l).gt.1.) then
503 if (l.gt.lmix(ig)) then
504 if (zmax(ig)-zmix(ig).lt.1.e-10) then
505 xxx(ig, l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
506 else
507 xxx(ig, l)=(zmax(ig)-zlev(ig, l))/(zmax(ig)-zmix(ig))
508 endif
509 if (idetr.eq.0) then
510 fraca(ig, l)=fracazmix(ig)
511 else if (idetr.eq.1) then
512 fraca(ig, l)=fracazmix(ig)*xxx(ig, l)
513 else if (idetr.eq.2) then
514 fraca(ig, l)=fracazmix(ig)*(1.-(1.-xxx(ig, l))**2)
515 else
516 fraca(ig, l)=fracazmix(ig)*xxx(ig, l)**2
517 endif
518 fraca(ig, l)=max(fraca(ig, l), 0.)
519 fraca(ig, l)=min(fraca(ig, l), 0.5)
520 fracd(ig, l)=1.-fraca(ig, l)
521 fracc(ig, l)=larg_cons(ig, l)/(r_aspect*zmax(ig))
522 endif
523 endif
524 enddo
525 enddo
526
527 print *, 'fin calcul fraca'
528
529 ! Calcul de fracd, wd
530 ! somme wa - wd = 0
531
532 do ig=1, ngrid
533 fm(ig, 1)=0.
534 fm(ig, nlay+1)=0.
535 enddo
536
537 do l=2, nlay
538 do ig=1, ngrid
539 fm(ig, l)=fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
540 if (entr(ig, l-1).lt.1e-10.and.fm(ig, l).gt.fm(ig, l-1) &
541 .and.l.gt.lmix(ig)) then
542 fm(ig, l)=fm(ig, l-1)
543 endif
544 enddo
545 do ig=1, ngrid
546 if(fracd(ig, l).lt.0.1) then
547 stop'fracd trop petit'
548 else
549 ! vitesse descendante "diagnostique"
550 wd(ig, l)=fm(ig, l)/(fracd(ig, l)*rhobarz(ig, l))
551 endif
552 enddo
553 enddo
554
555 do l=1, nlay
556 do ig=1, ngrid
557 masse(ig, l)=(pplev(ig, l)-pplev(ig, l+1))/RG
558 enddo
559 enddo
560
561 print *, '12 OK convect8'
562
563 ! calcul du transport vertical
564
565 !CR:redefinition du entr
566 do l=1, nlay
567 do ig=1, ngrid
568 detr(ig, l)=fm(ig, l)+entr(ig, l)-fm(ig, l+1)
569 if (detr(ig, l).lt.0.) then
570 entr(ig, l)=entr(ig, l)-detr(ig, l)
571 detr(ig, l)=0.
572 endif
573 enddo
574 enddo
575
576 if (w2di.eq.1) then
577 fm0=fm0+ptimestep*(fm-fm0)/tho
578 entr0=entr0+ptimestep*(entr-entr0)/tho
579 else
580 fm0=fm
581 entr0=entr
582 endif
583
584 if (1.eq.1) then
585 call dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse &
586 , zh, zdhadj, zha)
587 call dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse &
588 , zo, pdoadj, zoa)
589 else
590 call dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca &
591 , zh, zdhadj, zha)
592 call dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca &
593 , zo, pdoadj, zoa)
594 endif
595
596 if (1.eq.0) then
597 call dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse &
598 , fraca, zmax &
599 , zu, zv, pduadj, pdvadj, zua, zva)
600 else
601 call dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse &
602 , zu, pduadj, zua)
603 call dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse &
604 , zv, pdvadj, zva)
605 endif
606
607 do l=1, nlay
608 do ig=1, ngrid
609 zf=0.5*(fracc(ig, l)+fracc(ig, l+1))
610 zf2=zf/(1.-zf)
611 thetath2(ig, l)=zf2*(zha(ig, l)-zh(ig, l))**2
612 wth2(ig, l)=zf2*(0.5*(wa_moy(ig, l)+wa_moy(ig, l+1)))**2
613 enddo
614 enddo
615
616 do l=1, nlay
617 do ig=1, ngrid
618 pdtadj(ig, l)=zdhadj(ig, l)*zpspsk(ig, l)
619 enddo
620 enddo
621
622 print *, '14 OK convect8'
623
624 ! Calculs pour les sorties
625
626 if(sorties) then
627 do l=1, nlay
628 do ig=1, ngrid
629 zla(ig, l)=(1.-fracd(ig, l))*zmax(ig)
630 zld(ig, l)=fracd(ig, l)*zmax(ig)
631 if(1.-fracd(ig, l).gt.1.e-10) &
632 zwa(ig, l)=wd(ig, l)*fracd(ig, l)/(1.-fracd(ig, l))
633 enddo
634 enddo
635
636 isplit=isplit+1
637 endif
638
639 print *, '19 OK convect8'
640
641 end SUBROUTINE thermcell
642
643 end module thermcell_m

Properties

Name Value
svn:executable

  ViewVC Help
Powered by ViewVC 1.1.21