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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (hide 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 guez 54 module thermcell_m
2 guez 3
3 guez 47 IMPLICIT NONE
4 guez 3
5 guez 54 contains
6 guez 3
7 guez 54 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 guez 3
11 guez 150 ! 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 guez 157 ! environnement. La longueur "l_mix" contr\^ole l'efficacit\'e du
16 guez 150 ! m\'elange. Le calcul du transport des diff\'erentes esp\`eces se fait
17 guez 54 ! en prenant en compte :
18     ! 1. un flux de masse montant
19     ! 2. un flux de masse descendant
20 guez 150 ! 3. un entra\^inement
21     ! 4. un d\'etra\^inement
22 guez 3
23 guez 71 USE dimphy, ONLY : klev, klon
24 guez 54 USE suphec_m, ONLY : rd, rg, rkappa
25 guez 3
26 guez 54 ! arguments:
27 guez 3
28 guez 54 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 guez 91 REAL, intent(in):: pu(ngrid, nlay)
34     real pduadj(ngrid, nlay)
35     REAL, intent(in):: pv(ngrid, nlay)
36     real pdvadj(ngrid, nlay)
37 guez 54 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 guez 3
42 guez 54 integer idetr
43     save idetr
44     data idetr/3/
45 guez 3
46 guez 54 ! local:
47 guez 3
48 guez 54 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 guez 3
54 guez 150 real zmax(klon), zw, zw2(klon, klev+1), ztva(klon, klev)
55 guez 3
56 guez 54 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 guez 3
67 guez 54 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 guez 3
75 guez 105 integer isplit, nsplit
76 guez 54 parameter (nsplit=10)
77     data isplit/0/
78     save isplit
79 guez 3
80 guez 54 logical sorties
81     real rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
82     real zpspsk(klon, klev)
83 guez 3
84 guez 54 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 guez 3
95 guez 54 !CR:nouvelles variables
96     real f_star(klon, klev+1), entr_star(klon, klev)
97     real entr_star_tot(klon), entr_star2(klon)
98 guez 105 real f(klon)
99 guez 54 real zlevinter(klon)
100 guez 3
101 guez 54 EXTERNAL SCOPY
102 guez 3
103 guez 54 !-----------------------------------------------------------------------
104 guez 3
105 guez 54 ! initialisation:
106 guez 3
107 guez 54 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 guez 3
115 guez 54 ! incrementation eventuelle de tendances precedentes:
116 guez 3
117 guez 54 print *, '0 OK convect8'
118 guez 3
119 guez 54 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 guez 3
130 guez 54 print *, '1 OK convect8'
131 guez 3
132 guez 54 ! See notes, "thermcell.txt"
133     ! Calcul des altitudes des couches
134 guez 3
135 guez 54 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 guez 3
150 guez 54 ! Calcul des densites
151 guez 3
152 guez 54 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 guez 3
158 guez 54 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 guez 3
164 guez 54 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 guez 3
172 guez 54 ! Calcul de w2, quarre de w a partir de la cape
173     ! a partir de w2, on calcule wa, vitesse de l'ascendance
174 guez 3
175 guez 54 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
176     ! w2 est stoke dans wa
177 guez 3
178 guez 54 ! 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 guez 3
182 guez 54 ! 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 guez 3
187 guez 54 !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 guez 3
199 guez 54 !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 guez 3
209 guez 54 ! 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 guez 3
221 guez 54 ! 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 guez 3
249 guez 54 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 guez 3
256 guez 54 do k=1, klev+1
257     do ig=1, ngrid
258     zw2(ig, k)=0.
259     fmc(ig, k)=0.
260 guez 3
261 guez 54 f_star(ig, k)=0.
262 guez 3
263 guez 54 larg_cons(ig, k)=0.
264     larg_detr(ig, k)=0.
265     wa_moy(ig, k)=0.
266     enddo
267     enddo
268 guez 3
269 guez 54 do ig=1, ngrid
270     linter(ig)=1.
271     lmaxa(ig)=1
272     lmix(ig)=1
273     wmaxa(ig)=0.
274     enddo
275 guez 3
276 guez 54 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 guez 3
320 guez 54 ! 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 guez 3
339 guez 54 ! Determination de zw2 max
340     do ig=1, ngrid
341     wmax(ig)=0.
342     enddo
343 guez 3
344 guez 54 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 guez 3
358 guez 54 ! 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 guez 3
371 guez 54 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 guez 3
391 guez 54 ! 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 guez 3
404 guez 54 ! determination de l'indice du debut de la mixed layer ou w decroit
405 guez 3
406 guez 54 ! 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 guez 150 ! d'une couche est \'egale \`a la hauteur de la couche alimentante.
409 guez 54 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
410     ! de la vitesse d'entrainement horizontal dans la couche alimentante.
411 guez 3
412 guez 54 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 guez 3
422 guez 54 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 guez 3
433 guez 150 ! calcul de la fraction de la maille concern\'ee par l'ascendance en tenant
434 guez 54 ! compte de l'epluchage du thermique.
435 guez 3
436 guez 54 !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 guez 3
445 guez 54 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 guez 3
461 guez 54 if ((zmax(ig)-zmix(ig)).lt.0.) then
462     zmix(ig)=0.99*zmax(ig)
463     endif
464     enddo
465 guez 3
466 guez 54 ! 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 guez 3
476 guez 54 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 guez 3
500 guez 54 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 guez 3
527 guez 54 print *, 'fin calcul fraca'
528 guez 3
529 guez 54 ! Calcul de fracd, wd
530     ! somme wa - wd = 0
531 guez 3
532 guez 54 do ig=1, ngrid
533     fm(ig, 1)=0.
534     fm(ig, nlay+1)=0.
535     enddo
536 guez 3
537 guez 54 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 guez 3
555 guez 54 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 guez 3
561 guez 54 print *, '12 OK convect8'
562 guez 3
563 guez 54 ! calcul du transport vertical
564 guez 3
565 guez 54 !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 guez 3
576 guez 54 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 guez 3
584 guez 54 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 guez 3
596 guez 54 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 guez 3
607 guez 54 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 guez 3
616 guez 54 do l=1, nlay
617     do ig=1, ngrid
618     pdtadj(ig, l)=zdhadj(ig, l)*zpspsk(ig, l)
619     enddo
620     enddo
621 guez 3
622 guez 54 print *, '14 OK convect8'
623 guez 3
624 guez 54 ! Calculs pour les sorties
625 guez 3
626 guez 54 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 guez 3
636 guez 54 isplit=isplit+1
637     endif
638 guez 3
639 guez 54 print *, '19 OK convect8'
640 guez 3
641 guez 54 end SUBROUTINE thermcell
642 guez 3
643 guez 54 end module thermcell_m

Properties

Name Value
svn:executable

  ViewVC Help
Powered by ViewVC 1.1.21