source: LMDZ4/trunk/libf/phytherm/thermcell.F @ 961

Last change on this file since 961 was 814, checked in by Laurent Fairhead, 17 years ago

Rajout de la physique utilisant les thermiques FH
LF

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