/[lmdze]/trunk/libf/phylmd/Thermcell/thermcell.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/Thermcell/thermcell.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (hide annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years ago) by guez
Original Path: trunk/libf/phylmd/thermcell.f
File size: 35964 byte(s) 
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".
1 guez 3 SUBROUTINE thermcell(ngrid,nlay,ptimestep
2     s ,pplay,pplev,pphi
3     s ,pu,pv,pt,po
4     s ,pduadj,pdvadj,pdtadj,pdoadj
5     s ,fm0,entr0
6     c s ,pu_therm,pv_therm
7     s ,r_aspect,l_mix,w2di,tho)
8    
9     use dimens_m
10     use dimphy
11     use YOMCST
12     IMPLICIT NONE
13    
14     c=======================================================================
15     c
16     c Calcul du transport verticale dans la couche limite en presence
17     c de "thermiques" explicitement representes
18     c
19     c Réécriture à partir d'un listing papier à Habas, le 14/02/00
20     c
21     c le thermique est supposé homogène et dissipé par mélange avec
22     c son environnement. la longueur l_mix contrôle l'efficacité du
23     c mélange
24     c
25     c Le calcul du transport des différentes espèces se fait en prenant
26     c en compte:
27     c 1. un flux de masse montant
28     c 2. un flux de masse descendant
29     c 3. un entrainement
30     c 4. un detrainement
31     c
32     c=======================================================================
33    
34     c-----------------------------------------------------------------------
35     c declarations:
36     c -------------
37    
38    
39     c arguments:
40     c ----------
41    
42     INTEGER ngrid,nlay,w2di,tho
43     real ptimestep,l_mix,r_aspect
44     REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
45     REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
46     REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
47     REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
48 guez 10 REAL, intent(in):: pplay(ngrid,nlay)
49 guez 3 real, intent(in):: pplev(ngrid,nlay+1)
50     real pphi(ngrid,nlay)
51    
52     integer idetr
53     save idetr
54     data idetr/3/
55    
56     c local:
57     c ------
58    
59     INTEGER ig,k,l,lmaxa(klon),lmix(klon)
60     real zsortie1d(klon)
61     c CR: on remplace lmax(klon,klev+1)
62     INTEGER lmax(klon),lmin(klon),lentr(klon)
63     real linter(klon)
64     real zmix(klon), fracazmix(klon)
65     c RC
66     real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
67    
68     real zlev(klon,klev+1),zlay(klon,klev)
69     REAL zh(klon,klev),zdhadj(klon,klev)
70     REAL ztv(klon,klev)
71     real zu(klon,klev),zv(klon,klev),zo(klon,klev)
72     REAL wh(klon,klev+1)
73     real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
74     real zla(klon,klev+1)
75     real zwa(klon,klev+1)
76     real zld(klon,klev+1)
77     real zwd(klon,klev+1)
78     real zsortie(klon,klev)
79     real zva(klon,klev)
80     real zua(klon,klev)
81     real zoa(klon,klev)
82    
83     real zha(klon,klev)
84     real wa_moy(klon,klev+1)
85     real fraca(klon,klev+1)
86     real fracc(klon,klev+1)
87     real zf,zf2
88     real thetath2(klon,klev),wth2(klon,klev)
89     common/comtherm/thetath2,wth2
90    
91     real count_time
92     integer isplit,nsplit,ialt
93     parameter (nsplit=10)
94     data isplit/0/
95     save isplit
96    
97     logical sorties
98     real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
99     real zpspsk(klon,klev)
100    
101     c real wmax(klon,klev),wmaxa(klon)
102     real wmax(klon),wmaxa(klon)
103     real wa(klon,klev,klev+1)
104     real wd(klon,klev+1)
105     real larg_part(klon,klev,klev+1)
106     real fracd(klon,klev+1)
107     real xxx(klon,klev+1)
108     real larg_cons(klon,klev+1)
109     real larg_detr(klon,klev+1)
110     real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
111     real pu_therm(klon,klev),pv_therm(klon,klev)
112     real fm(klon,klev+1),entr(klon,klev)
113     real fmc(klon,klev+1)
114    
115     cCR:nouvelles variables
116     real f_star(klon,klev+1),entr_star(klon,klev)
117     real entr_star_tot(klon),entr_star2(klon)
118     real f(klon), f0(klon)
119     real zlevinter(klon)
120     logical first
121     data first /.false./
122     save first
123     cRC
124    
125     character*2 str2
126     character*10 str10
127    
128     LOGICAL vtest(klon),down
129    
130     EXTERNAL SCOPY
131    
132     integer ncorrec,ll
133     save ncorrec
134     data ncorrec/0/
135    
136     c
137     c-----------------------------------------------------------------------
138     c initialisation:
139     c ---------------
140     c
141     sorties=.true.
142     IF(ngrid.NE.klon) THEN
143     PRINT*
144     PRINT*,'STOP dans convadj'
145     PRINT*,'ngrid =',ngrid
146     PRINT*,'klon =',klon
147     ENDIF
148     c
149     c-----------------------------------------------------------------------
150     c incrementation eventuelle de tendances precedentes:
151     c ---------------------------------------------------
152    
153     print*,'0 OK convect8'
154    
155     DO 1010 l=1,nlay
156     DO 1015 ig=1,ngrid
157     zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
158     zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
159     zu(ig,l)=pu(ig,l)
160     zv(ig,l)=pv(ig,l)
161     zo(ig,l)=po(ig,l)
162     ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
163     1015 CONTINUE
164     1010 CONTINUE
165    
166     print*,'1 OK convect8'
167     c --------------------
168     c
169     c
170     c + + + + + + + + + + +
171     c
172     c
173     c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
174     c wh,wt,wo ...
175     c
176     c + + + + + + + + + + + zh,zu,zv,zo,rho
177     c
178     c
179     c -------------------- zlev(1)
180     c \\\\\\\\\\\\\\\\\\\\
181     c
182     c
183    
184     c-----------------------------------------------------------------------
185     c Calcul des altitudes des couches
186     c-----------------------------------------------------------------------
187    
188     do l=2,nlay
189     do ig=1,ngrid
190     zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
191     enddo
192     enddo
193     do ig=1,ngrid
194     zlev(ig,1)=0.
195     zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
196     enddo
197     do l=1,nlay
198     do ig=1,ngrid
199     zlay(ig,l)=pphi(ig,l)/RG
200     enddo
201     enddo
202    
203     c print*,'2 OK convect8'
204     c-----------------------------------------------------------------------
205     c Calcul des densites
206     c-----------------------------------------------------------------------
207    
208     do l=1,nlay
209     do ig=1,ngrid
210     rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
211     enddo
212     enddo
213    
214     do l=2,nlay
215     do ig=1,ngrid
216     rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
217     enddo
218     enddo
219    
220     do k=1,nlay
221     do l=1,nlay+1
222     do ig=1,ngrid
223     wa(ig,k,l)=0.
224     enddo
225     enddo
226     enddo
227    
228     c print*,'3 OK convect8'
229     c------------------------------------------------------------------
230     c Calcul de w2, quarre de w a partir de la cape
231     c a partir de w2, on calcule wa, vitesse de l'ascendance
232     c
233     c ATTENTION: Dans cette version, pour cause d'economie de memoire,
234     c w2 est stoke dans wa
235     c
236     c ATTENTION: dans convect8, on n'utilise le calcule des wa
237     c independants par couches que pour calculer l'entrainement
238     c a la base et la hauteur max de l'ascendance.
239     c
240     c Indicages:
241     c l'ascendance provenant du niveau k traverse l'interface l avec
242     c une vitesse wa(k,l).
243     c
244     c --------------------
245     c
246     c + + + + + + + + + +
247     c
248     c wa(k,l) ---- -------------------- l
249     c /\
250     c /||\ + + + + + + + + + +
251     c ||
252     c || --------------------
253     c ||
254     c || + + + + + + + + + +
255     c ||
256     c || --------------------
257     c ||__
258     c |___ + + + + + + + + + + k
259     c
260     c --------------------
261     c
262     c
263     c
264     c------------------------------------------------------------------
265    
266     cCR: ponderation entrainement des couches instables
267     cdef des entr_star tels que entr=f*entr_star
268     do l=1,klev
269     do ig=1,ngrid
270     entr_star(ig,l)=0.
271     enddo
272     enddo
273     c determination de la longueur de la couche d entrainement
274     do ig=1,ngrid
275     lentr(ig)=1
276     enddo
277    
278     con ne considere que les premieres couches instables
279     do k=nlay-2,1,-1
280     do ig=1,ngrid
281     if (ztv(ig,k).gt.ztv(ig,k+1).and.
282     s ztv(ig,k+1).le.ztv(ig,k+2)) then
283     lentr(ig)=k
284     endif
285     enddo
286     enddo
287    
288     c determination du lmin: couche d ou provient le thermique
289     do ig=1,ngrid
290     lmin(ig)=1
291     enddo
292     do ig=1,ngrid
293     do l=nlay,2,-1
294     if (ztv(ig,l-1).gt.ztv(ig,l)) then
295     lmin(ig)=l-1
296     endif
297     enddo
298     enddo
299     c
300     c definition de l'entrainement des couches
301     do l=1,klev-1
302     do ig=1,ngrid
303     if (ztv(ig,l).gt.ztv(ig,l+1).and.
304     s l.ge.lmin(ig).and.l.le.lentr(ig)) then
305     entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
306     s (zlev(ig,l+1)-zlev(ig,l))
307     endif
308     enddo
309     enddo
310     c pas de thermique si couches 1->5 stables
311     do ig=1,ngrid
312     if (lmin(ig).gt.5) then
313     do l=1,klev
314     entr_star(ig,l)=0.
315     enddo
316     endif
317     enddo
318     c calcul de l entrainement total
319     do ig=1,ngrid
320     entr_star_tot(ig)=0.
321     enddo
322     do ig=1,ngrid
323     do k=1,klev
324     entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
325     enddo
326     enddo
327     c
328     print*,'fin calcul entr_star'
329     do k=1,klev
330     do ig=1,ngrid
331     ztva(ig,k)=ztv(ig,k)
332     enddo
333     enddo
334     cRC
335     c print*,'7 OK convect8'
336     do k=1,klev+1
337     do ig=1,ngrid
338     zw2(ig,k)=0.
339     fmc(ig,k)=0.
340     cCR
341     f_star(ig,k)=0.
342     cRC
343     larg_cons(ig,k)=0.
344     larg_detr(ig,k)=0.
345     wa_moy(ig,k)=0.
346     enddo
347     enddo
348    
349     c print*,'8 OK convect8'
350     do ig=1,ngrid
351     linter(ig)=1.
352     lmaxa(ig)=1
353     lmix(ig)=1
354     wmaxa(ig)=0.
355     enddo
356    
357     cCR:
358     do l=1,nlay-2
359     do ig=1,ngrid
360     if (ztv(ig,l).gt.ztv(ig,l+1)
361     s .and.entr_star(ig,l).gt.1.e-10
362     s .and.zw2(ig,l).lt.1e-10) then
363     f_star(ig,l+1)=entr_star(ig,l)
364     ctest:calcul de dteta
365     zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
366     s *(zlev(ig,l+1)-zlev(ig,l))
367     s *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
368     larg_detr(ig,l)=0.
369     else if ((zw2(ig,l).ge.1e-10).and.
370     s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
371     f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
372     ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
373     s *ztv(ig,l))/f_star(ig,l+1)
374     zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
375     s 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
376     s *(zlev(ig,l+1)-zlev(ig,l))
377     endif
378     c determination de zmax continu par interpolation lineaire
379     if (zw2(ig,l+1).lt.0.) then
380     ctest
381     if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
382     print*,'pb linter'
383     endif
384     linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
385     s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
386     zw2(ig,l+1)=0.
387     lmaxa(ig)=l
388     else
389     if (zw2(ig,l+1).lt.0.) then
390     print*,'pb1 zw2<0'
391     endif
392     wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
393     endif
394     if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
395     c lmix est le niveau de la couche ou w (wa_moy) est maximum
396     lmix(ig)=l+1
397     wmaxa(ig)=wa_moy(ig,l+1)
398     endif
399     enddo
400     enddo
401     print*,'fin calcul zw2'
402     c
403     c Calcul de la couche correspondant a la hauteur du thermique
404     do ig=1,ngrid
405     lmax(ig)=lentr(ig)
406     enddo
407     do ig=1,ngrid
408     do l=nlay,lentr(ig)+1,-1
409     if (zw2(ig,l).le.1.e-10) then
410     lmax(ig)=l-1
411     endif
412     enddo
413     enddo
414     c pas de thermique si couches 1->5 stables
415     do ig=1,ngrid
416     if (lmin(ig).gt.5) then
417     lmax(ig)=1
418     lmin(ig)=1
419     endif
420     enddo
421     c
422     c Determination de zw2 max
423     do ig=1,ngrid
424     wmax(ig)=0.
425     enddo
426    
427     do l=1,nlay
428     do ig=1,ngrid
429     if (l.le.lmax(ig)) then
430     if (zw2(ig,l).lt.0.)then
431     print*,'pb2 zw2<0'
432     endif
433     zw2(ig,l)=sqrt(zw2(ig,l))
434     wmax(ig)=max(wmax(ig),zw2(ig,l))
435     else
436     zw2(ig,l)=0.
437     endif
438     enddo
439     enddo
440    
441     c Longueur caracteristique correspondant a la hauteur des thermiques.
442     do ig=1,ngrid
443     zmax(ig)=0.
444     zlevinter(ig)=zlev(ig,1)
445     enddo
446     do ig=1,ngrid
447     c calcul de zlevinter
448     zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
449     s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
450     s -zlev(ig,lmax(ig)))
451     zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
452     enddo
453    
454     print*,'avant fermeture'
455     c Fermeture,determination de f
456     do ig=1,ngrid
457     entr_star2(ig)=0.
458     enddo
459     do ig=1,ngrid
460     if (entr_star_tot(ig).LT.1.e-10) then
461     f(ig)=0.
462     else
463     do k=lmin(ig),lentr(ig)
464     entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
465     s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
466     enddo
467     c Nouvelle fermeture
468     f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
469     s *entr_star2(ig))*entr_star_tot(ig)
470     ctest
471     c if (first) then
472     c f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
473     c s *wmax(ig))
474     c endif
475     endif
476     c f0(ig)=f(ig)
477     c first=.true.
478     enddo
479     print*,'apres fermeture'
480    
481     c Calcul de l'entrainement
482     do k=1,klev
483     do ig=1,ngrid
484     entr(ig,k)=f(ig)*entr_star(ig,k)
485     enddo
486     enddo
487     c Calcul des flux
488     do ig=1,ngrid
489     do l=1,lmax(ig)-1
490     fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
491     enddo
492     enddo
493    
494     cRC
495    
496    
497     c print*,'9 OK convect8'
498     c print*,'WA1 ',wa_moy
499    
500     c determination de l'indice du debut de la mixed layer ou w decroit
501    
502     c calcul de la largeur de chaque ascendance dans le cas conservatif.
503     c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
504     c d'une couche est égale à la hauteur de la couche alimentante.
505     c La vitesse maximale dans l'ascendance est aussi prise comme estimation
506     c de la vitesse d'entrainement horizontal dans la couche alimentante.
507    
508     do l=2,nlay
509     do ig=1,ngrid
510     if (l.le.lmaxa(ig)) then
511     zw=max(wa_moy(ig,l),1.e-10)
512     larg_cons(ig,l)=zmax(ig)*r_aspect
513     s *fmc(ig,l)/(rhobarz(ig,l)*zw)
514     endif
515     enddo
516     enddo
517    
518     do l=2,nlay
519     do ig=1,ngrid
520     if (l.le.lmaxa(ig)) then
521     c if (idetr.eq.0) then
522     c cette option est finalement en dur.
523     if ((l_mix*zlev(ig,l)).lt.0.)then
524     print*,'pb l_mix*zlev<0'
525     endif
526     larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
527     c else if (idetr.eq.1) then
528     c larg_detr(ig,l)=larg_cons(ig,l)
529     c s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
530     c else if (idetr.eq.2) then
531     c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
532     c s *sqrt(wa_moy(ig,l))
533     c else if (idetr.eq.4) then
534     c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
535     c s *wa_moy(ig,l)
536     c endif
537     endif
538     enddo
539     enddo
540    
541     c print*,'10 OK convect8'
542     c print*,'WA2 ',wa_moy
543     c calcul de la fraction de la maille concernée par l'ascendance en tenant
544     c compte de l'epluchage du thermique.
545     c
546     cCR def de zmix continu (profil parabolique des vitesses)
547     do ig=1,ngrid
548     if (lmix(ig).gt.1.) then
549     c test
550     if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
551     s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
552     s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
553     s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
554     s then
555     c
556     zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
557     s *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
558     s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
559     s *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
560     s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
561     s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
562     s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
563     s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
564     else
565     zmix(ig)=zlev(ig,lmix(ig))
566     print*,'pb zmix'
567     endif
568     else
569     zmix(ig)=0.
570     endif
571     ctest
572     if ((zmax(ig)-zmix(ig)).lt.0.) then
573     zmix(ig)=0.99*zmax(ig)
574     c print*,'pb zmix>zmax'
575     endif
576     enddo
577     c
578     c calcul du nouveau lmix correspondant
579     do ig=1,ngrid
580     do l=1,klev
581     if (zmix(ig).ge.zlev(ig,l).and.
582     s zmix(ig).lt.zlev(ig,l+1)) then
583     lmix(ig)=l
584     endif
585     enddo
586     enddo
587     c
588     do l=2,nlay
589     do ig=1,ngrid
590     if(larg_cons(ig,l).gt.1.) then
591     c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
592     fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
593     s /(r_aspect*zmax(ig))
594     c test
595     fraca(ig,l)=max(fraca(ig,l),0.)
596     fraca(ig,l)=min(fraca(ig,l),0.5)
597     fracd(ig,l)=1.-fraca(ig,l)
598     fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
599     else
600     c wa_moy(ig,l)=0.
601     fraca(ig,l)=0.
602     fracc(ig,l)=0.
603     fracd(ig,l)=1.
604     endif
605     enddo
606     enddo
607     cCR: calcul de fracazmix
608     do ig=1,ngrid
609     fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
610     s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
611     s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
612     s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
613     enddo
614     c
615     do l=2,nlay
616     do ig=1,ngrid
617     if(larg_cons(ig,l).gt.1.) then
618     if (l.gt.lmix(ig)) then
619     ctest
620     if (zmax(ig)-zmix(ig).lt.1.e-10) then
621     c print*,'pb xxx'
622     xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
623     else
624     xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
625     endif
626     if (idetr.eq.0) then
627     fraca(ig,l)=fracazmix(ig)
628     else if (idetr.eq.1) then
629     fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
630     else if (idetr.eq.2) then
631     fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
632     else
633     fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
634     endif
635     c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
636     fraca(ig,l)=max(fraca(ig,l),0.)
637     fraca(ig,l)=min(fraca(ig,l),0.5)
638     fracd(ig,l)=1.-fraca(ig,l)
639     fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
640     endif
641     endif
642     enddo
643     enddo
644    
645     print*,'fin calcul fraca'
646     c print*,'11 OK convect8'
647     c print*,'Ea3 ',wa_moy
648     c------------------------------------------------------------------
649     c Calcul de fracd, wd
650     c somme wa - wd = 0
651     c------------------------------------------------------------------
652    
653    
654     do ig=1,ngrid
655     fm(ig,1)=0.
656     fm(ig,nlay+1)=0.
657     enddo
658    
659     do l=2,nlay
660     do ig=1,ngrid
661     fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
662     cCR:test
663     if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
664     s .and.l.gt.lmix(ig)) then
665     fm(ig,l)=fm(ig,l-1)
666     c write(1,*)'ajustement fm, l',l
667     endif
668     c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
669     cRC
670     enddo
671     do ig=1,ngrid
672     if(fracd(ig,l).lt.0.1) then
673     stop'fracd trop petit'
674     else
675     c vitesse descendante "diagnostique"
676     wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
677     endif
678     enddo
679     enddo
680    
681     do l=1,nlay
682     do ig=1,ngrid
683     c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
684     masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
685     enddo
686     enddo
687    
688     print*,'12 OK convect8'
689     c print*,'WA4 ',wa_moy
690     cc------------------------------------------------------------------
691     c calcul du transport vertical
692     c------------------------------------------------------------------
693    
694     go to 4444
695     c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
696     do l=2,nlay-1
697     do ig=1,ngrid
698     if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
699     s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
700     c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
701     c s ,fm(ig,l+1)*ptimestep
702     c s ,' M=',masse(ig,l),masse(ig,l+1)
703     endif
704     enddo
705     enddo
706    
707     do l=1,nlay
708     do ig=1,ngrid
709     if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
710     c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
711     c s ,entr(ig,l)*ptimestep
712     c s ,' M=',masse(ig,l)
713     endif
714     enddo
715     enddo
716    
717     do l=1,nlay
718     do ig=1,ngrid
719     if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
720     c print*,'WARN!!! fm exagere ig=',ig,' l=',l
721     c s ,' FM=',fm(ig,l)
722     endif
723     if(.not.masse(ig,l).ge.1.e-10
724     s .or..not.masse(ig,l).le.1.e4) then
725     endif
726     if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
727     c print*,'WARN!!! entr exagere ig=',ig,' l=',l
728     c s ,' E=',entr(ig,l)
729     endif
730     enddo
731     enddo
732    
733     4444 continue
734    
735     cCR:redefinition du entr
736     do l=1,nlay
737     do ig=1,ngrid
738     detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
739     if (detr(ig,l).lt.0.) then
740     entr(ig,l)=entr(ig,l)-detr(ig,l)
741     detr(ig,l)=0.
742     c print*,'WARNING !!! detrainement negatif ',ig,l
743     endif
744     enddo
745     enddo
746     cRC
747     if (w2di.eq.1) then
748     fm0=fm0+ptimestep*(fm-fm0)/float(tho)
749     entr0=entr0+ptimestep*(entr-entr0)/float(tho)
750     else
751     fm0=fm
752     entr0=entr
753     endif
754    
755     if (1.eq.1) then
756     call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
757     . ,zh,zdhadj,zha)
758     call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
759     . ,zo,pdoadj,zoa)
760     else
761     call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
762     . ,zh,zdhadj,zha)
763     call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
764     . ,zo,pdoadj,zoa)
765     endif
766    
767     if (1.eq.0) then
768     call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
769     . ,fraca,zmax
770     . ,zu,zv,pduadj,pdvadj,zua,zva)
771     else
772     call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
773     . ,zu,pduadj,zua)
774     call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
775     . ,zv,pdvadj,zva)
776     endif
777    
778     do l=1,nlay
779     do ig=1,ngrid
780     zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
781     zf2=zf/(1.-zf)
782     thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
783     wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
784     enddo
785     enddo
786    
787    
788    
789     c print*,'13 OK convect8'
790     c print*,'WA5 ',wa_moy
791     do l=1,nlay
792     do ig=1,ngrid
793     pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
794     enddo
795     enddo
796    
797    
798     c do l=1,nlay
799     c do ig=1,ngrid
800     c if(abs(pdtadj(ig,l))*86400..gt.500.) then
801     c print*,'WARN!!! ig=',ig,' l=',l
802     c s ,' pdtadj=',pdtadj(ig,l)
803     c endif
804     c if(abs(pdoadj(ig,l))*86400..gt.1.) then
805     c print*,'WARN!!! ig=',ig,' l=',l
806     c s ,' pdoadj=',pdoadj(ig,l)
807     c endif
808     c enddo
809     c enddo
810    
811     print*,'14 OK convect8'
812     c------------------------------------------------------------------
813     c Calculs pour les sorties
814     c------------------------------------------------------------------
815    
816     if(sorties) then
817     do l=1,nlay
818     do ig=1,ngrid
819     zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
820     zld(ig,l)=fracd(ig,l)*zmax(ig)
821     if(1.-fracd(ig,l).gt.1.e-10)
822     s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
823     enddo
824     enddo
825    
826     cdeja fait
827     c do l=1,nlay
828     c do ig=1,ngrid
829     c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
830     c if (detr(ig,l).lt.0.) then
831     c entr(ig,l)=entr(ig,l)-detr(ig,l)
832     c detr(ig,l)=0.
833     c print*,'WARNING !!! detrainement negatif ',ig,l
834     c endif
835     c enddo
836     c enddo
837    
838     c print*,'15 OK convect8'
839    
840     isplit=isplit+1
841    
842    
843     goto 123
844     123 continue
845    
846     endif
847    
848     c if(wa_moy(1,4).gt.1.e-10) stop
849    
850     print*,'19 OK convect8'
851     return
852     end
853    
854     subroutine dqthermcell(ngrid,nlay,ptimestep,fm,entr,masse
855     . ,q,dq,qa)
856     use dimens_m
857     use dimphy
858     implicit none
859    
860     c=======================================================================
861     c
862     c Calcul du transport verticale dans la couche limite en presence
863     c de "thermiques" explicitement representes
864     c calcul du dq/dt une fois qu'on connait les ascendances
865     c
866     c=======================================================================
867    
868    
869     integer ngrid,nlay
870    
871     real ptimestep
872     real, intent(in):: masse(ngrid,nlay)
873     real fm(ngrid,nlay+1)
874     real entr(ngrid,nlay)
875     real q(ngrid,nlay)
876     real dq(ngrid,nlay)
877    
878     real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
879    
880     integer ig,k
881    
882     c calcul du detrainement
883    
884     do k=1,nlay
885     do ig=1,ngrid
886     detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
887     enddo
888     enddo
889    
890     c calcul de la valeur dans les ascendances
891     do ig=1,ngrid
892     qa(ig,1)=q(ig,1)
893     enddo
894    
895     do k=2,nlay
896     do ig=1,ngrid
897     if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
898     s 1.e-5*masse(ig,k)) then
899     qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))
900     s /(fm(ig,k+1)+detr(ig,k))
901     else
902     qa(ig,k)=q(ig,k)
903     endif
904     enddo
905     enddo
906    
907     do k=2,nlay
908     do ig=1,ngrid
909     c wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
910     wqd(ig,k)=fm(ig,k)*q(ig,k)
911     enddo
912     enddo
913     do ig=1,ngrid
914     wqd(ig,1)=0.
915     wqd(ig,nlay+1)=0.
916     enddo
917    
918     do k=1,nlay
919     do ig=1,ngrid
920     dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)
921     s -wqd(ig,k)+wqd(ig,k+1))
922     s /masse(ig,k)
923     enddo
924     enddo
925    
926     return
927     end
928     subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
929     . ,fraca,larga
930     . ,u,v,du,dv,ua,va)
931     use dimens_m
932     use dimphy
933     implicit none
934    
935     c=======================================================================
936     c
937     c Calcul du transport verticale dans la couche limite en presence
938     c de "thermiques" explicitement representes
939     c calcul du dq/dt une fois qu'on connait les ascendances
940     c
941     c=======================================================================
942    
943    
944     integer ngrid,nlay
945    
946     real ptimestep
947     real masse(ngrid,nlay),fm(ngrid,nlay+1)
948     real fraca(ngrid,nlay+1)
949     real larga(ngrid)
950     real entr(ngrid,nlay)
951     real u(ngrid,nlay)
952     real ua(ngrid,nlay)
953     real du(ngrid,nlay)
954     real v(ngrid,nlay)
955     real va(ngrid,nlay)
956     real dv(ngrid,nlay)
957    
958     real qa(klon,klev),detr(klon,klev)
959     real wvd(klon,klev+1),wud(klon,klev+1)
960     real gamma0,gamma(klon,klev+1)
961     real dua,dva
962     integer iter
963    
964     integer ig,k
965    
966     c calcul du detrainement
967    
968     do k=1,nlay
969     do ig=1,ngrid
970     detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
971     enddo
972     enddo
973    
974     c calcul de la valeur dans les ascendances
975     do ig=1,ngrid
976     ua(ig,1)=u(ig,1)
977     va(ig,1)=v(ig,1)
978     enddo
979    
980     do k=2,nlay
981     do ig=1,ngrid
982     if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
983     s 1.e-5*masse(ig,k)) then
984     c On itère sur la valeur du coeff de freinage.
985     c gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
986     gamma0=masse(ig,k)
987     s *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
988     s *0.5/larga(ig)
989     c gamma0=0.
990     c la première fois on multiplie le coefficient de freinage
991     c par le module du vent dans la couche en dessous.
992     dua=ua(ig,k-1)-u(ig,k-1)
993     dva=va(ig,k-1)-v(ig,k-1)
994     do iter=1,5
995     gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
996     ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
997     s +(entr(ig,k)+gamma(ig,k))*u(ig,k))
998     s /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
999     va(ig,k)=(fm(ig,k)*va(ig,k-1)
1000     s +(entr(ig,k)+gamma(ig,k))*v(ig,k))
1001     s /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
1002     c print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
1003     dua=ua(ig,k)-u(ig,k)
1004     dva=va(ig,k)-v(ig,k)
1005     enddo
1006     else
1007     ua(ig,k)=u(ig,k)
1008     va(ig,k)=v(ig,k)
1009     gamma(ig,k)=0.
1010     endif
1011     enddo
1012     enddo
1013    
1014     do k=2,nlay
1015     do ig=1,ngrid
1016     wud(ig,k)=fm(ig,k)*u(ig,k)
1017     wvd(ig,k)=fm(ig,k)*v(ig,k)
1018     enddo
1019     enddo
1020     do ig=1,ngrid
1021     wud(ig,1)=0.
1022     wud(ig,nlay+1)=0.
1023     wvd(ig,1)=0.
1024     wvd(ig,nlay+1)=0.
1025     enddo
1026    
1027     do k=1,nlay
1028     do ig=1,ngrid
1029     du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
1030     s -(entr(ig,k)+gamma(ig,k))*u(ig,k)
1031     s -wud(ig,k)+wud(ig,k+1))
1032     s /masse(ig,k)
1033     dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
1034     s -(entr(ig,k)+gamma(ig,k))*v(ig,k)
1035     s -wvd(ig,k)+wvd(ig,k+1))
1036     s /masse(ig,k)
1037     enddo
1038     enddo
1039    
1040     return
1041     end
1042     subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
1043     . ,q,dq,qa)
1044     use dimens_m
1045     use dimphy
1046     implicit none
1047    
1048     c=======================================================================
1049     c
1050     c Calcul du transport verticale dans la couche limite en presence
1051     c de "thermiques" explicitement representes
1052     c calcul du dq/dt une fois qu'on connait les ascendances
1053     c
1054     c=======================================================================
1055    
1056    
1057     integer ngrid,nlay
1058    
1059     real ptimestep
1060     real masse(ngrid,nlay),fm(ngrid,nlay+1)
1061     real entr(ngrid,nlay),frac(ngrid,nlay)
1062     real q(ngrid,nlay)
1063     real dq(ngrid,nlay)
1064    
1065     real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
1066     real qe(klon,klev),zf,zf2
1067    
1068     integer ig,k
1069    
1070     c calcul du detrainement
1071    
1072     do k=1,nlay
1073     do ig=1,ngrid
1074     detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1075     enddo
1076     enddo
1077    
1078     c calcul de la valeur dans les ascendances
1079     do ig=1,ngrid
1080     qa(ig,1)=q(ig,1)
1081     qe(ig,1)=q(ig,1)
1082     enddo
1083    
1084     do k=2,nlay
1085     do ig=1,ngrid
1086     if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
1087     s 1.e-5*masse(ig,k)) then
1088     zf=0.5*(frac(ig,k)+frac(ig,k+1))
1089     zf2=1./(1.-zf)
1090     qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+zf2*entr(ig,k)*q(ig,k))
1091     s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2)
1092     qe(ig,k)=(q(ig,k)-zf*qa(ig,k))*zf2
1093     else
1094     qa(ig,k)=q(ig,k)
1095     qe(ig,k)=q(ig,k)
1096     endif
1097     enddo
1098     enddo
1099    
1100     do k=2,nlay
1101     do ig=1,ngrid
1102     c wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
1103     wqd(ig,k)=fm(ig,k)*qe(ig,k)
1104     enddo
1105     enddo
1106     do ig=1,ngrid
1107     wqd(ig,1)=0.
1108     wqd(ig,nlay+1)=0.
1109     enddo
1110    
1111     do k=1,nlay
1112     do ig=1,ngrid
1113     dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*qe(ig,k)
1114     s -wqd(ig,k)+wqd(ig,k+1))
1115     s /masse(ig,k)
1116     enddo
1117     enddo
1118    
1119     return
1120     end
1121     subroutine dvthermcell2(ngrid,nlay,ptimestep,fm,entr,masse
1122     . ,fraca,larga
1123     . ,u,v,du,dv,ua,va)
1124     use dimens_m
1125     use dimphy
1126     implicit none
1127    
1128     c=======================================================================
1129     c
1130     c Calcul du transport verticale dans la couche limite en presence
1131     c de "thermiques" explicitement representes
1132     c calcul du dq/dt une fois qu'on connait les ascendances
1133     c
1134     c=======================================================================
1135    
1136    
1137     integer ngrid,nlay
1138    
1139     real ptimestep
1140     real masse(ngrid,nlay),fm(ngrid,nlay+1)
1141     real fraca(ngrid,nlay+1)
1142     real larga(ngrid)
1143     real entr(ngrid,nlay)
1144     real u(ngrid,nlay)
1145     real ua(ngrid,nlay)
1146     real du(ngrid,nlay)
1147     real v(ngrid,nlay)
1148     real va(ngrid,nlay)
1149     real dv(ngrid,nlay)
1150    
1151     real qa(klon,klev),detr(klon,klev),zf,zf2
1152     real wvd(klon,klev+1),wud(klon,klev+1)
1153     real gamma0,gamma(klon,klev+1)
1154     real ue(klon,klev),ve(klon,klev)
1155     real dua,dva
1156     integer iter
1157    
1158     integer ig,k
1159    
1160     c calcul du detrainement
1161    
1162     do k=1,nlay
1163     do ig=1,ngrid
1164     detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1165     enddo
1166     enddo
1167    
1168     c calcul de la valeur dans les ascendances
1169     do ig=1,ngrid
1170     ua(ig,1)=u(ig,1)
1171     va(ig,1)=v(ig,1)
1172     ue(ig,1)=u(ig,1)
1173     ve(ig,1)=v(ig,1)
1174     enddo
1175    
1176     do k=2,nlay
1177     do ig=1,ngrid
1178     if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
1179     s 1.e-5*masse(ig,k)) then
1180     c On itère sur la valeur du coeff de freinage.
1181     c gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
1182     gamma0=masse(ig,k)
1183     s *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
1184     s *0.5/larga(ig)
1185     s *1.
1186     c s *0.5
1187     c gamma0=0.
1188     zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
1189     zf=0.
1190     zf2=1./(1.-zf)
1191     c la première fois on multiplie le coefficient de freinage
1192     c par le module du vent dans la couche en dessous.
1193     dua=ua(ig,k-1)-u(ig,k-1)
1194     dva=va(ig,k-1)-v(ig,k-1)
1195     do iter=1,5
1196     c On choisit une relaxation lineaire.
1197     gamma(ig,k)=gamma0
1198     c On choisit une relaxation quadratique.
1199     gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
1200     ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
1201     s +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))
1202     s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
1203     s +gamma(ig,k))
1204     va(ig,k)=(fm(ig,k)*va(ig,k-1)
1205     s +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))
1206     s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
1207     s +gamma(ig,k))
1208     c print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
1209     dua=ua(ig,k)-u(ig,k)
1210     dva=va(ig,k)-v(ig,k)
1211     ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
1212     ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
1213     enddo
1214     else
1215     ua(ig,k)=u(ig,k)
1216     va(ig,k)=v(ig,k)
1217     ue(ig,k)=u(ig,k)
1218     ve(ig,k)=v(ig,k)
1219     gamma(ig,k)=0.
1220     endif
1221     enddo
1222     enddo
1223    
1224     do k=2,nlay
1225     do ig=1,ngrid
1226     wud(ig,k)=fm(ig,k)*ue(ig,k)
1227     wvd(ig,k)=fm(ig,k)*ve(ig,k)
1228     enddo
1229     enddo
1230     do ig=1,ngrid
1231     wud(ig,1)=0.
1232     wud(ig,nlay+1)=0.
1233     wvd(ig,1)=0.
1234     wvd(ig,nlay+1)=0.
1235     enddo
1236    
1237     do k=1,nlay
1238     do ig=1,ngrid
1239     du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
1240     s -(entr(ig,k)+gamma(ig,k))*ue(ig,k)
1241     s -wud(ig,k)+wud(ig,k+1))
1242     s /masse(ig,k)
1243     dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
1244     s -(entr(ig,k)+gamma(ig,k))*ve(ig,k)
1245     s -wvd(ig,k)+wvd(ig,k+1))
1246     s /masse(ig,k)
1247     enddo
1248     enddo
1249    
1250     return
1251     end

Properties

Name Value
svn:executable
  ViewVC Help
Powered by ViewVC 1.1.21