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