source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/thermcell.F @ 1440

Last change on this file since 1440 was 987, checked in by Laurent Fairhead, 16 years ago

Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9
et ORCHIDEE YM
LF

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