source: trunk/LMDZ.TERRE/libf/phylmd/thermcell.F @ 134

Last change on this file since 134 was 1, checked in by emillour, 14 years ago

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