source: lmdz_wrf/trunk/WRFV3/lmdz/thermcell_old.F @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 190.2 KB
Line 
1      SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep
2     s                  ,pplay,pplev,pphi
3     s                  ,pu,pv,pt,po
4     s                  ,pduadj,pdvadj,pdtadj,pdoadj
5     s                  ,fm0,entr0
6c    s                  ,pu_therm,pv_therm
7     s                  ,r_aspect,l_mix,w2di,tho)
8
9      USE dimphy
10      IMPLICIT NONE
11
12c=======================================================================
13c
14c   Calcul du transport verticale dans la couche limite en presence
15c   de "thermiques" explicitement representes
16c
17c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
18c
19c   le thermique est supposé homogène et dissipé par mélange avec
20c   son environnement. la longueur l_mix contrôle l'efficacité du
21c   mélange
22c
23c   Le calcul du transport des différentes espèces se fait en prenant
24c   en compte:
25c     1. un flux de masse montant
26c     2. un flux de masse descendant
27c     3. un entrainement
28c     4. un detrainement
29c
30c=======================================================================
31
32c-----------------------------------------------------------------------
33c   declarations:
34c   -------------
35
36#include "dimensions.h"
37cccc#include "dimphy.h"
38#include "YOMCST.h"
39
40c   arguments:
41c   ----------
42
43      INTEGER ngrid,nlay,w2di
44      REAL tho
45      real ptimestep,l_mix,r_aspect
46      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
47      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
48      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
49      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
50      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
51      real pphi(ngrid,nlay)
52
53      integer idetr
54      save idetr
55      data idetr/3/
56c$OMP THREADPRIVATE(idetr)
57
58c   local:
59c   ------
60
61      INTEGER ig,k,l,lmax(klon,klev+1),lmaxa(klon),lmix(klon)
62      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
63
64      real zlev(klon,klev+1),zlay(klon,klev)
65      REAL zh(klon,klev),zdhadj(klon,klev)
66      REAL ztv(klon,klev)
67      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
68      REAL wh(klon,klev+1)
69      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
70      real zla(klon,klev+1)
71      real zwa(klon,klev+1)
72      real zld(klon,klev+1)
73      real zwd(klon,klev+1)
74      real zsortie(klon,klev)
75      real zva(klon,klev)
76      real zua(klon,klev)
77      real zoa(klon,klev)
78
79      real zha(klon,klev)
80      real wa_moy(klon,klev+1)
81      real fraca(klon,klev+1)
82      real fracc(klon,klev+1)
83      real zf,zf2
84      real thetath2(klon,klev),wth2(klon,klev)
85!      common/comtherm/thetath2,wth2
86
87      real count_time
88      integer isplit,nsplit,ialt
89      parameter (nsplit=10)
90      data isplit/0/
91      save isplit
92c$OMP THREADPRIVATE(isplit)
93
94      logical sorties
95      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
96      real zpspsk(klon,klev)
97
98      real wmax(klon,klev),wmaxa(klon)
99
100      real wa(klon,klev,klev+1)
101      real wd(klon,klev+1)
102      real larg_part(klon,klev,klev+1)
103      real fracd(klon,klev+1)
104      real xxx(klon,klev+1)
105      real larg_cons(klon,klev+1)
106      real larg_detr(klon,klev+1)
107      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
108      real pu_therm(klon,klev),pv_therm(klon,klev)
109      real fm(klon,klev+1),entr(klon,klev)
110      real fmc(klon,klev+1)
111
112      character (len=2) :: str2
113      character (len=10) :: str10
114
115      character (len=20) :: modname='thermcell2002'
116      character (len=80) :: abort_message
117
118      LOGICAL vtest(klon),down
119
120      EXTERNAL SCOPY
121
122      integer ncorrec,ll
123      save ncorrec
124      data ncorrec/0/
125c$OMP THREADPRIVATE(ncorrec)
126
127c
128c-----------------------------------------------------------------------
129c   initialisation:
130c   ---------------
131c
132       sorties=.true.
133      IF(ngrid.NE.klon) THEN
134         PRINT*
135         PRINT*,'STOP dans convadj'
136         PRINT*,'ngrid    =',ngrid
137         PRINT*,'klon  =',klon
138      ENDIF
139c
140c-----------------------------------------------------------------------
141c   incrementation eventuelle de tendances precedentes:
142c   ---------------------------------------------------
143
144      print*,'0 OK convect8'
145
146      DO 1010 l=1,nlay
147         DO 1015 ig=1,ngrid
148            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
149            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
150            zu(ig,l)=pu(ig,l)
151            zv(ig,l)=pv(ig,l)
152            zo(ig,l)=po(ig,l)
153            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
1541015     CONTINUE
1551010  CONTINUE
156
157c     print*,'1 OK convect8'
158c                       --------------------
159c
160c
161c                       + + + + + + + + + + +
162c
163c
164c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
165c  wh,wt,wo ...
166c
167c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
168c
169c
170c                       --------------------   zlev(1)
171c                       \\\\\\\\\\\\\\\\\\\\
172c
173c
174
175c-----------------------------------------------------------------------
176c   Calcul des altitudes des couches
177c-----------------------------------------------------------------------
178
179      do l=2,nlay
180         do ig=1,ngrid
181            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
182         enddo
183      enddo
184      do ig=1,ngrid
185         zlev(ig,1)=0.
186         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
187      enddo
188      do l=1,nlay
189         do ig=1,ngrid
190            zlay(ig,l)=pphi(ig,l)/RG
191         enddo
192      enddo
193
194c     print*,'2 OK convect8'
195c-----------------------------------------------------------------------
196c   Calcul des densites
197c-----------------------------------------------------------------------
198
199      do l=1,nlay
200         do ig=1,ngrid
201            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
202         enddo
203      enddo
204
205      do l=2,nlay
206         do ig=1,ngrid
207            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
208         enddo
209      enddo
210
211      do k=1,nlay
212         do l=1,nlay+1
213            do ig=1,ngrid
214               wa(ig,k,l)=0.
215            enddo
216         enddo
217      enddo
218
219c     print*,'3 OK convect8'
220c------------------------------------------------------------------
221c   Calcul de w2, quarre de w a partir de la cape
222c   a partir de w2, on calcule wa, vitesse de l'ascendance
223c
224c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
225c   w2 est stoke dans wa
226c
227c   ATTENTION: dans convect8, on n'utilise le calcule des wa
228c   independants par couches que pour calculer l'entrainement
229c   a la base et la hauteur max de l'ascendance.
230c
231c   Indicages:
232c   l'ascendance provenant du niveau k traverse l'interface l avec
233c   une vitesse wa(k,l).
234c
235c                       --------------------
236c
237c                       + + + + + + + + + +
238c
239c  wa(k,l)   ----       --------------------    l
240c             /\
241c            /||\       + + + + + + + + + +
242c             ||
243c             ||        --------------------
244c             ||
245c             ||        + + + + + + + + + +
246c             ||
247c             ||        --------------------
248c             ||__
249c             |___      + + + + + + + + + +     k
250c
251c                       --------------------
252c
253c
254c
255c------------------------------------------------------------------
256
257
258      do k=1,nlay-1
259         do ig=1,ngrid
260            wa(ig,k,k)=0.
261            wa(ig,k,k+1)=2.*RG*(ztv(ig,k)-ztv(ig,k+1))/ztv(ig,k+1)
262     s      *(zlev(ig,k+1)-zlev(ig,k))
263         enddo
264         do l=k+1,nlay-1
265            do ig=1,ngrid
266               wa(ig,k,l+1)=wa(ig,k,l)+
267     s         2.*RG*(ztv(ig,k)-ztv(ig,l))/ztv(ig,l)
268     s         *(zlev(ig,l+1)-zlev(ig,l))
269            enddo
270         enddo
271         do ig=1,ngrid
272            wa(ig,k,nlay+1)=0.
273         enddo
274      enddo
275
276c     print*,'4 OK convect8'
277c Calcul de la couche correspondant a la hauteur du thermique
278      do k=1,nlay-1
279         do ig=1,ngrid
280            lmax(ig,k)=k
281         enddo
282         do l=nlay,k+1,-1
283            do ig=1,ngrid
284               if(wa(ig,k,l).le.1.e-10) lmax(ig,k)=l-1
285            enddo
286         enddo
287      enddo
288
289c     print*,'5 OK convect8'
290c   Calcule du w max du thermique
291      do k=1,nlay
292      do ig=1,ngrid
293         wmax(ig,k)=0.
294      enddo
295      enddo
296
297      do k=1,nlay-1
298         do l=k,nlay
299            do ig=1,ngrid
300               if (l.le.lmax(ig,k)) then
301                  wa(ig,k,l)=sqrt(wa(ig,k,l))
302                  wmax(ig,k)=max(wmax(ig,k),wa(ig,k,l))
303               else
304                  wa(ig,k,l)=0.
305               endif
306            enddo
307         enddo
308      enddo
309
310      do k=1,nlay-1
311         do ig=1,ngrid
312             pu_therm(ig,k)=sqrt(wmax(ig,k))
313             pv_therm(ig,k)=sqrt(wmax(ig,k))
314         enddo
315      enddo
316
317c     print*,'6 OK convect8'
318c   Longueur caracteristique correspondant a la hauteur des thermiques.
319      do  ig=1,ngrid
320         zmax(ig)=500.
321      enddo
322c     print*,'LMAX LMAX LMAX '
323      do k=1,nlay-1
324         do  ig=1,ngrid
325            zmax(ig)=max(zmax(ig),zlev(ig,lmax(ig,k))-zlev(ig,k))
326         enddo
327c     print*,k,lmax(1,k)
328      enddo
329c     print*,'ZMAX ZMAX ZMAX ',zmax
330c      call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX      ')
331
332c   Calcul de l'entrainement.
333c   Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
334c   de la couche d'alimentation en partant du principe que la vitesse
335c   maximum dans l'ascendance est la vitesse d'entrainement horizontale.
336      do k=1,nlay
337         do ig=1,ngrid
338            zzz=rho(ig,k)*wmax(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
339     s    /(zmax(ig)*r_aspect)
340            if(w2di.eq.2) then
341               entr(ig,k)=entr(ig,k)+
342     s         ptimestep*(zzz-entr(ig,k))/tho
343            else
344               entr(ig,k)=zzz
345            endif
346            ztva(ig,k)=ztv(ig,k)
347         enddo
348      enddo
349
350c     print*,'7 OK convect8'
351      do k=1,klev+1
352         do ig=1,ngrid
353            zw2(ig,k)=0.
354            fmc(ig,k)=0.
355            larg_cons(ig,k)=0.
356            larg_detr(ig,k)=0.
357            wa_moy(ig,k)=0.
358         enddo
359      enddo
360
361c     print*,'8 OK convect8'
362      do ig=1,ngrid
363         lmaxa(ig)=1
364         lmix(ig)=1
365         wmaxa(ig)=0.
366      enddo
367
368
369      do l=1,nlay-2
370         do ig=1,ngrid
371c           if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)) then
372c         print*,'COUCOU ',l,zw2(ig,l),ztv(ig,l),ztv(ig,l+1)
373            if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)
374     s       .and.entr(ig,l).gt.1.e-10) then
375c        print*,'COUCOU cas 1'
376c   Initialisation de l'ascendance
377c              lmix(ig)=1
378               ztva(ig,l)=ztv(ig,l)
379               fmc(ig,l)=0.
380               fmc(ig,l+1)=entr(ig,l)
381               zw2(ig,l)=0.
382c     if (.not.ztv(ig,l+1).gt.150.) then
383c     print*,'ig,l+1,ztv(ig,l+1)'
384c     print*, ig,l+1,ztv(ig,l+1)
385c     endif
386               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
387     s         *(zlev(ig,l+1)-zlev(ig,l))
388               larg_detr(ig,l)=0.
389            else if (zw2(ig,l).ge.1.e-10.and.
390     .               fmc(ig,l)+entr(ig,l).gt.1.e-10) then
391c   Incrementation...
392               fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
393c     if (.not.fmc(ig,l+1).gt.1.e-15) then
394c     print*,'ig,l+1,fmc(ig,l+1)'
395c     print*, ig,l+1,fmc(ig,l+1)
396c     print*,'Fmc ',(fmc(ig,ll),ll=1,klev+1)
397c     print*,'W2 ',(zw2(ig,ll),ll=1,klev+1)
398c     print*,'Tv ',(ztv(ig,ll),ll=1,klev)
399c     print*,'Entr ',(entr(ig,ll),ll=1,klev)
400c     endif
401               ztva(ig,l)=(fmc(ig,l)*ztva(ig,l-1)+entr(ig,l)*ztv(ig,l))
402     s          /fmc(ig,l+1)
403c  mise a jour de la vitesse ascendante (l'air entraine de la couche
404c  consideree commence avec une vitesse nulle).
405               zw2(ig,l+1)=zw2(ig,l)*(fmc(ig,l)/fmc(ig,l+1))**2+
406     s         2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
407     s         *(zlev(ig,l+1)-zlev(ig,l))
408            endif
409            if (zw2(ig,l+1).lt.0.) then
410               zw2(ig,l+1)=0.
411               lmaxa(ig)=l
412            else
413               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
414            endif
415            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
416c   lmix est le niveau de la couche ou w (wa_moy) est maximum
417               lmix(ig)=l+1
418               wmaxa(ig)=wa_moy(ig,l+1)
419            endif
420c        print*,'COUCOU cas 2 LMIX=',lmix(ig),wa_moy(ig,l+1),wmaxa(ig)
421         enddo
422      enddo
423
424c     print*,'9 OK convect8'
425c     print*,'WA1 ',wa_moy
426
427c   determination de l'indice du debut de la mixed layer ou w decroit
428
429c   calcul de la largeur de chaque ascendance dans le cas conservatif.
430c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
431c   d'une couche est égale à la hauteur de la couche alimentante.
432c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
433c   de la vitesse d'entrainement horizontal dans la couche alimentante.
434
435      do l=2,nlay
436         do ig=1,ngrid
437            if (l.le.lmaxa(ig)) then
438               zw=max(wa_moy(ig,l),1.e-10)
439               larg_cons(ig,l)=zmax(ig)*r_aspect
440     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
441            endif
442         enddo
443      enddo
444
445      do l=2,nlay
446         do ig=1,ngrid
447            if (l.le.lmaxa(ig)) then
448c              if (idetr.eq.0) then
449c  cette option est finalement en dur.
450                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
451c              else if (idetr.eq.1) then
452c                 larg_detr(ig,l)=larg_cons(ig,l)
453c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
454c              else if (idetr.eq.2) then
455c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
456c    s            *sqrt(wa_moy(ig,l))
457c              else if (idetr.eq.4) then
458c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
459c    s            *wa_moy(ig,l)
460c              endif
461            endif
462         enddo
463       enddo
464
465c     print*,'10 OK convect8'
466c     print*,'WA2 ',wa_moy
467c   calcul de la fraction de la maille concernée par l'ascendance en tenant
468c   compte de l'epluchage du thermique.
469
470      do l=2,nlay
471         do ig=1,ngrid
472            if(larg_cons(ig,l).gt.1.) then
473c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
474               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
475     s            /(r_aspect*zmax(ig))
476               if(l.gt.lmix(ig)) then
477                  xxx(ig,l)=(lmaxa(ig)+1.-l) / (lmaxa(ig)+1.-lmix(ig))
478           if (idetr.eq.0) then
479               fraca(ig,l)=fraca(ig,lmix(ig))
480           else if (idetr.eq.1) then
481               fraca(ig,l)=fraca(ig,lmix(ig))*xxx(ig,l)
482           else if (idetr.eq.2) then
483               fraca(ig,l)=fraca(ig,lmix(ig))*(1.-(1.-xxx(ig,l))**2)
484           else
485               fraca(ig,l)=fraca(ig,lmix(ig))*xxx(ig,l)**2
486           endif
487               endif
488c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
489               fraca(ig,l)=max(fraca(ig,l),0.)
490               fraca(ig,l)=min(fraca(ig,l),0.5)
491               fracd(ig,l)=1.-fraca(ig,l)
492               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
493            else
494c              wa_moy(ig,l)=0.
495               fraca(ig,l)=0.
496               fracc(ig,l)=0.
497               fracd(ig,l)=1.
498            endif
499         enddo
500      enddo
501
502c     print*,'11 OK convect8'
503c     print*,'Ea3 ',wa_moy
504c------------------------------------------------------------------
505c   Calcul de fracd, wd
506c   somme wa - wd = 0
507c------------------------------------------------------------------
508
509
510      do ig=1,ngrid
511         fm(ig,1)=0.
512         fm(ig,nlay+1)=0.
513      enddo
514
515      do l=2,nlay
516           do ig=1,ngrid
517              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
518           enddo
519         do ig=1,ngrid
520            if(fracd(ig,l).lt.0.1) then
521              abort_message = 'fracd trop petit'
522              CALL abort_gcm (modname,abort_message,1)
523           else
524c    vitesse descendante "diagnostique"
525               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
526            endif
527         enddo
528      enddo
529
530      do l=1,nlay
531         do ig=1,ngrid
532c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
533            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
534         enddo
535      enddo
536
537c     print*,'12 OK convect8'
538c     print*,'WA4 ',wa_moy
539cc------------------------------------------------------------------
540c   calcul du transport vertical
541c------------------------------------------------------------------
542
543      go to 4444
544c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
545      do l=2,nlay-1
546         do ig=1,ngrid
547            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
548     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
549c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
550c    s         ,fm(ig,l+1)*ptimestep
551c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
552             endif
553         enddo
554      enddo
555
556      do l=1,nlay
557         do ig=1,ngrid
558            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
559c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
560c    s         ,entr(ig,l)*ptimestep
561c    s         ,'   M=',masse(ig,l)
562             endif
563         enddo
564      enddo
565
566      do l=1,nlay
567         do ig=1,ngrid
568            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
569c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
570c    s         ,'   FM=',fm(ig,l)
571            endif
572            if(.not.masse(ig,l).ge.1.e-10
573     s         .or..not.masse(ig,l).le.1.e4) then
574c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
575c    s         ,'   M=',masse(ig,l)
576c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
577c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
578c     print*,'zlev(ig,l+1),zlev(ig,l)'
579c    s                ,zlev(ig,l+1),zlev(ig,l)
580c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
581c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
582            endif
583            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
584c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
585c    s         ,'   E=',entr(ig,l)
586            endif
587         enddo
588      enddo
589
5904444   continue
591
592      if (w2di.eq.1) then
593         fm0=fm0+ptimestep*(fm-fm0)/tho
594         entr0=entr0+ptimestep*(entr-entr0)/tho
595      else
596         fm0=fm
597         entr0=entr
598      endif
599
600      if (1.eq.1) then
601         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
602     .    ,zh,zdhadj,zha)
603         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
604     .    ,zo,pdoadj,zoa)
605      else
606         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
607     .    ,zh,zdhadj,zha)
608         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
609     .    ,zo,pdoadj,zoa)
610      endif
611
612      if (1.eq.0) then
613         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
614     .    ,fraca,zmax
615     .    ,zu,zv,pduadj,pdvadj,zua,zva)
616      else
617         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
618     .    ,zu,pduadj,zua)
619         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
620     .    ,zv,pdvadj,zva)
621      endif
622
623      do l=1,nlay
624         do ig=1,ngrid
625            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
626            zf2=zf/(1.-zf)
627            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
628            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
629         enddo
630      enddo
631
632
633
634c     print*,'13 OK convect8'
635c     print*,'WA5 ',wa_moy
636      do l=1,nlay
637         do ig=1,ngrid
638            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
639         enddo
640      enddo
641
642
643c     do l=1,nlay
644c        do ig=1,ngrid
645c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
646c     print*,'WARN!!! ig=',ig,'  l=',l
647c    s         ,'   pdtadj=',pdtadj(ig,l)
648c           endif
649c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
650c     print*,'WARN!!! ig=',ig,'  l=',l
651c    s         ,'   pdoadj=',pdoadj(ig,l)
652c           endif
653c        enddo
654c      enddo
655
656c     print*,'14 OK convect8'
657c------------------------------------------------------------------
658c   Calculs pour les sorties
659c------------------------------------------------------------------
660
661      if(sorties) then
662      do l=1,nlay
663         do ig=1,ngrid
664            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
665            zld(ig,l)=fracd(ig,l)*zmax(ig)
666            if(1.-fracd(ig,l).gt.1.e-10)
667     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
668         enddo
669      enddo
670
671      do l=1,nlay
672         do ig=1,ngrid
673            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
674            if (detr(ig,l).lt.0.) then
675                entr(ig,l)=entr(ig,l)-detr(ig,l)
676                detr(ig,l)=0.
677c     print*,'WARNING !!! detrainement negatif ',ig,l
678            endif
679         enddo
680      enddo
681
682c     print*,'15 OK convect8'
683
684      isplit=isplit+1
685
686
687c #define und
688        goto 123
689#ifdef und
690      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
691      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
692      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
693      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
694      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
695      CALL writeg1d(1,nlay,zla,'la      ','la      ')
696      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
697      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
698      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
699      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
700      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
701      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
702      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
703      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
704      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
705      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
706      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
707      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
708      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
709      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
710      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
711      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
712      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
713      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
714
715      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
716      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
717      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
718c   recalcul des flux en diagnostique...
719c     print*,'PAS DE TEMPS ',ptimestep
720       call dt2F(pplev,pplay,pt,pdtadj,wh)
721      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
722#endif
723123   continue
724! #define troisD
725#ifdef troisD
726c       if (sorties) then
727      print*,'Debut des wrgradsfi'
728
729c      print*,'16 OK convect8'
730         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
731         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
732         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
733         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
734         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
735         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
736c      print*,'WA6 ',wa_moy
737         call wrgradsfi(1,nlay,zla,'la        ','la        ')
738         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
739         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
740         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
741         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
742         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
743         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
744         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
745         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
746         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
747         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
748         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
749         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
750         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
751         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
752         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
753         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
754         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
755         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
756         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
757         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
758         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
759         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
760         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
761         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
762         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
763
764         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
765         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
766         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
767
768
769c      print*,'17 OK convect8'
770
771         do k=1,klev/10
772            write(str2,'(i2.2)') k
773            str10='wa'//str2
774            do l=1,nlay
775               do ig=1,ngrid
776                  zsortie(ig,l)=wa(ig,k,l)
777               enddo
778            enddo   
779            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
780            do l=1,nlay
781               do ig=1,ngrid
782                  zsortie(ig,l)=larg_part(ig,k,l)
783               enddo
784            enddo
785            str10='la'//str2
786            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
787         enddo
788
789
790c     print*,'18 OK convect8'
791c      endif
792      print*,'Fin des wrgradsfi'
793#endif
794
795      endif
796
797c     if(wa_moy(1,4).gt.1.e-10) stop
798
799c     print*,'19 OK convect8'
800      return
801      end
802
803      SUBROUTINE thermcell_cld(ngrid,nlay,ptimestep
804     s                  ,pplay,pplev,pphi,zlev,debut
805     s                  ,pu,pv,pt,po
806     s                  ,pduadj,pdvadj,pdtadj,pdoadj
807     s                  ,fm0,entr0,zqla,lmax
808     s                  ,zmax_sec,wmax_sec,zw_sec,lmix_sec
809     s                  ,ratqscth,ratqsdiff
810c    s                  ,pu_therm,pv_therm
811     s                  ,r_aspect,l_mix,w2di,tho)
812
813      USE dimphy
814      IMPLICIT NONE
815
816c=======================================================================
817c
818c   Calcul du transport verticale dans la couche limite en presence
819c   de "thermiques" explicitement representes
820c
821c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
822c
823c   le thermique est supposé homogène et dissipé par mélange avec
824c   son environnement. la longueur l_mix contrôle l'efficacité du
825c   mélange
826c
827c   Le calcul du transport des différentes espèces se fait en prenant
828c   en compte:
829c     1. un flux de masse montant
830c     2. un flux de masse descendant
831c     3. un entrainement
832c     4. un detrainement
833c
834c=======================================================================
835
836c-----------------------------------------------------------------------
837c   declarations:
838c   -------------
839
840#include "dimensions.h"
841cccc#include "dimphy.h"
842#include "YOMCST.h"
843#include "YOETHF.h"
844#include "FCTTRE.h"
845
846c   arguments:
847c   ----------
848
849      INTEGER ngrid,nlay,w2di
850      REAL tho
851      real ptimestep,l_mix,r_aspect
852      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
853      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
854      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
855      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
856      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
857      real pphi(ngrid,nlay)
858
859      integer idetr
860      save idetr
861      data idetr/3/
862c$OMP THREADPRIVATE(idetr)
863
864c   local:
865c   ------
866
867      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
868      real zsortie1d(klon)
869c CR: on remplace lmax(klon,klev+1)
870      INTEGER lmax(klon),lmin(klon),lentr(klon)
871      real linter(klon)
872      real zmix(klon), fracazmix(klon)
873      real alpha
874      save alpha
875      data alpha/1./
876c$OMP THREADPRIVATE(alpha)
877
878c RC
879      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
880      real zmax_sec(klon)
881      real zmax_sec2(klon)
882      real zw_sec(klon,klev+1)
883      INTEGER lmix_sec(klon)
884      real w_est(klon,klev+1)
885con garde le zmax du pas de temps precedent
886c      real zmax0(klon)
887c      save zmax0
888c      real zmix0(klon)
889c      save zmix0
890      REAL, SAVE, ALLOCATABLE :: zmax0(:), zmix0(:)
891c$OMP THREADPRIVATE(zmax0, zmix0)
892
893      real zlev(klon,klev+1),zlay(klon,klev)
894      real deltaz(klon,klev)
895      REAL zh(klon,klev),zdhadj(klon,klev)
896      real zthl(klon,klev),zdthladj(klon,klev)
897      REAL ztv(klon,klev)
898      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
899      real zl(klon,klev)
900      REAL wh(klon,klev+1)
901      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
902      real zla(klon,klev+1)
903      real zwa(klon,klev+1)
904      real zld(klon,klev+1)
905      real zwd(klon,klev+1)
906      real zsortie(klon,klev)
907      real zva(klon,klev)
908      real zua(klon,klev)
909      real zoa(klon,klev)
910
911      real zta(klon,klev)
912      real zha(klon,klev)
913      real wa_moy(klon,klev+1)
914      real fraca(klon,klev+1)
915      real fracc(klon,klev+1)
916      real zf,zf2
917      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
918      real q2(klon,klev)
919      real dtheta(klon,klev)
920!      common/comtherm/thetath2,wth2
921   
922      real ratqscth(klon,klev)
923      real sum
924      real sumdiff
925      real ratqsdiff(klon,klev)
926      real count_time
927      integer isplit,nsplit,ialt
928      parameter (nsplit=10)
929      data isplit/0/
930      save isplit
931c$OMP THREADPRIVATE(isplit)
932
933      logical sorties
934      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
935      real zpspsk(klon,klev)
936
937c     real wmax(klon,klev),wmaxa(klon)
938      real wmax(klon),wmaxa(klon)
939      real wmax_sec(klon)
940      real wmax_sec2(klon)
941      real wa(klon,klev,klev+1)
942      real wd(klon,klev+1)
943      real larg_part(klon,klev,klev+1)
944      real fracd(klon,klev+1)
945      real xxx(klon,klev+1)
946      real larg_cons(klon,klev+1)
947      real larg_detr(klon,klev+1)
948      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
949      real massetot(klon,klev)
950      real detr0(klon,klev)
951      real alim0(klon,klev)
952      real pu_therm(klon,klev),pv_therm(klon,klev)
953      real fm(klon,klev+1),entr(klon,klev)
954      real fmc(klon,klev+1)
955
956      real zcor,zdelta,zcvm5,qlbef
957      real Tbef(klon),qsatbef(klon)
958      real dqsat_dT,DT,num,denom
959      REAL REPS,RLvCp,DDT0
960      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
961cCR niveau de condensation
962      real nivcon(klon)
963      real zcon(klon)
964      real zqsat(klon,klev)
965      real zqsatth(klon,klev)
966      PARAMETER (DDT0=.01)
967
968
969cCR:nouvelles variables
970      real f_star(klon,klev+1),entr_star(klon,klev)
971      real detr_star(klon,klev)
972      real alim_star_tot(klon),alim_star2(klon)
973      real entr_star_tot(klon)
974      real detr_star_tot(klon)
975      real alim_star(klon,klev)
976      real alim(klon,klev)
977      real nu(klon,klev)
978      real nu_e(klon,klev)
979      real nu_min
980      real nu_max
981      real nu_r
982      real f(klon)
983c      real f(klon), f0(klon)
984c     save f0
985      REAL,SAVE, ALLOCATABLE :: f0(:)
986c$OMP THREADPRIVATE(f0)
987
988      real f_old
989      real zlevinter(klon)
990      logical, save :: first = .true.
991c$OMP THREADPRIVATE(first)
992c      data first /.false./
993c      save first
994      logical nuage
995c      save nuage
996      logical boucle
997      logical therm
998      logical debut
999      logical rale
1000      integer test(klon)
1001      integer signe_zw2
1002cRC
1003
1004      character*2 str2
1005      character*10 str10
1006
1007      character (len=20) :: modname='thermcell_cld'
1008      character (len=80) :: abort_message
1009
1010      LOGICAL vtest(klon),down
1011      LOGICAL Zsat(klon)
1012
1013      EXTERNAL SCOPY
1014
1015      integer ncorrec,ll
1016      save ncorrec
1017      data ncorrec/0/
1018c$OMP THREADPRIVATE(ncorrec)
1019
1020c
1021
1022c-----------------------------------------------------------------------
1023c   initialisation:
1024c   ---------------
1025c
1026      if (first) then
1027        allocate(zmix0(klon))
1028        allocate(zmax0(klon))
1029        allocate(f0(klon))
1030        first=.false.
1031      endif
1032
1033       sorties=.false.
1034c     print*,'NOUVEAU DETR PLUIE '
1035      IF(ngrid.NE.klon) THEN
1036         PRINT*
1037         PRINT*,'STOP dans convadj'
1038         PRINT*,'ngrid    =',ngrid
1039         PRINT*,'klon  =',klon
1040      ENDIF
1041c
1042c Initialisation
1043      RLvCp = RLVTT/RCPD
1044      REPS  = RD/RV
1045cinitialisations de zqsat
1046      DO ll=1,nlay
1047         DO ig=1,ngrid
1048            zqsat(ig,ll)=0.
1049            zqsatth(ig,ll)=0.
1050         ENDDO
1051      ENDDO
1052c
1053con met le first a true pour le premier passage de la journée
1054      do ig=1,klon
1055            test(ig)=0
1056      enddo
1057      if (debut) then
1058         do ig=1,klon
1059            test(ig)=1
1060            f0(ig)=0.
1061            zmax0(ig)=0.
1062         enddo
1063      endif
1064      do ig=1,klon     
1065         if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
1066            test(ig)=1
1067         endif
1068      enddo
1069c     do ig=1,klon
1070c        print*,'test(ig)',test(ig),zmax0(ig)
1071c     enddo
1072      nuage=.false.
1073c-----------------------------------------------------------------------
1074cAM Calcul de T,q,ql a partir de Tl et qT
1075c   ---------------------------------------------------
1076c
1077c Pr Tprec=Tl calcul de qsat
1078c Si qsat>qT T=Tl, q=qT
1079c Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
1080c On cherche DDT < DDT0
1081c
1082c defaut
1083       DO  ll=1,nlay
1084         DO ig=1,ngrid
1085            zo(ig,ll)=po(ig,ll)
1086            zl(ig,ll)=0.
1087            zh(ig,ll)=pt(ig,ll)
1088         EndDO
1089       EndDO
1090       do ig=1,ngrid
1091          Zsat(ig)=.false.
1092       enddo
1093c
1094c
1095       DO ll=1,nlay
1096c les points insatures sont definitifs
1097         DO ig=1,ngrid
1098            Tbef(ig)=pt(ig,ll)
1099            zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1100            qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
1101            qsatbef(ig)=MIN(0.5,qsatbef(ig))
1102            zcor=1./(1.-retv*qsatbef(ig))
1103            qsatbef(ig)=qsatbef(ig)*zcor
1104            Zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig)) .gt. 1.e-10)
1105         EndDO
1106
1107         DO ig=1,ngrid
1108           if (Zsat(ig).and.(1.eq.1)) then
1109            qlbef=max(0.,po(ig,ll)-qsatbef(ig))
1110c si sature: ql est surestime, d'ou la sous-relax
1111            DT = 0.5*RLvCp*qlbef
1112c            write(18,*),'DT0=',DT
1113c on pourra enchainer 2 ou 3 calculs sans Do while
1114            do while (abs(DT).gt.DDT0)
1115c il faut verifier si c,a conserve quand on repasse en insature ...
1116              Tbef(ig)=Tbef(ig)+DT
1117              zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1118              qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
1119              qsatbef(ig)=MIN(0.5,qsatbef(ig))
1120              zcor=1./(1.-retv*qsatbef(ig))
1121              qsatbef(ig)=qsatbef(ig)*zcor
1122c on veut le signe de qlbef
1123              qlbef=po(ig,ll)-qsatbef(ig)
1124              zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1125              zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
1126              zcor=1./(1.-retv*qsatbef(ig))
1127              dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
1128              num=-Tbef(ig)+pt(ig,ll)+RLvCp*qlbef
1129              denom=1.+RLvCp*dqsat_dT
1130              if (denom.lt.1.e-10) then
1131                  print*,'pb denom'
1132              endif
1133              DT=num/denom
1134            enddo
1135c on ecrit de maniere conservative (sat ou non)
1136            zl(ig,ll) = max(0.,qlbef)
1137c          T = Tl +Lv/Cp ql
1138            zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
1139            zo(ig,ll) = po(ig,ll)-zl(ig,ll)
1140           endif
1141con ecrit zqsat
1142            zqsat(ig,ll)=qsatbef(ig)     
1143         EndDO
1144       EndDO
1145cAM fin
1146c
1147c-----------------------------------------------------------------------
1148c   incrementation eventuelle de tendances precedentes:
1149c   ---------------------------------------------------
1150
1151c     print*,'0 OK convect8'
1152
1153      DO 1010 l=1,nlay
1154         DO 1015 ig=1,ngrid
1155             zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
1156c            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
1157c            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
1158            zu(ig,l)=pu(ig,l)
1159            zv(ig,l)=pv(ig,l)
1160c            zo(ig,l)=po(ig,l)
1161c            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
1162cAM attention zh est maintenant le profil de T et plus le profil de theta !
1163c
1164c   T-> Theta
1165            ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
1166cAM Theta_v
1167            ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l))
1168     s           -zl(ig,l))
1169cAM Thetal
1170            zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
1171c           
11721015     CONTINUE
11731010  CONTINUE
1174
1175c     print*,'1 OK convect8'
1176c                       --------------------
1177c
1178c
1179c                       + + + + + + + + + + +
1180c
1181c
1182c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
1183c  wh,wt,wo ...
1184c
1185c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
1186c
1187c
1188c                       --------------------   zlev(1)
1189c                       \\\\\\\\\\\\\\\\\\\\
1190c
1191c
1192
1193c-----------------------------------------------------------------------
1194c   Calcul des altitudes des couches
1195c-----------------------------------------------------------------------
1196
1197      do l=2,nlay
1198         do ig=1,ngrid
1199            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
1200         enddo
1201      enddo
1202      do ig=1,ngrid
1203         zlev(ig,1)=0.
1204         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
1205      enddo
1206      do l=1,nlay
1207         do ig=1,ngrid
1208            zlay(ig,l)=pphi(ig,l)/RG
1209         enddo
1210      enddo
1211ccalcul de deltaz
1212      do l=1,nlay
1213         do ig=1,ngrid
1214            deltaz(ig,l)=zlev(ig,l+1)-zlev(ig,l)
1215         enddo
1216      enddo
1217
1218c     print*,'2 OK convect8'
1219c-----------------------------------------------------------------------
1220c   Calcul des densites
1221c-----------------------------------------------------------------------
1222
1223      do l=1,nlay
1224         do ig=1,ngrid
1225c            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
1226             rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*ztv(ig,l))
1227         enddo
1228      enddo
1229
1230      do l=2,nlay
1231         do ig=1,ngrid
1232            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
1233         enddo
1234      enddo
1235
1236      do k=1,nlay
1237         do l=1,nlay+1
1238            do ig=1,ngrid
1239               wa(ig,k,l)=0.
1240            enddo
1241         enddo
1242      enddo
1243cCr:ajout:calcul de la masse
1244      do l=1,nlay
1245         do ig=1,ngrid
1246c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1247            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
1248         enddo
1249      enddo
1250c     print*,'3 OK convect8'
1251c------------------------------------------------------------------
1252c   Calcul de w2, quarre de w a partir de la cape
1253c   a partir de w2, on calcule wa, vitesse de l'ascendance
1254c
1255c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
1256c   w2 est stoke dans wa
1257c
1258c   ATTENTION: dans convect8, on n'utilise le calcule des wa
1259c   independants par couches que pour calculer l'entrainement
1260c   a la base et la hauteur max de l'ascendance.
1261c
1262c   Indicages:
1263c   l'ascendance provenant du niveau k traverse l'interface l avec
1264c   une vitesse wa(k,l).
1265c
1266c                       --------------------
1267c
1268c                       + + + + + + + + + +
1269c
1270c  wa(k,l)   ----       --------------------    l
1271c             /\
1272c            /||\       + + + + + + + + + +
1273c             ||
1274c             ||        --------------------
1275c             ||
1276c             ||        + + + + + + + + + +
1277c             ||
1278c             ||        --------------------
1279c             ||__
1280c             |___      + + + + + + + + + +     k
1281c
1282c                       --------------------
1283c
1284c
1285c
1286c------------------------------------------------------------------
1287
1288cCR: ponderation entrainement des couches instables
1289cdef des alim_star tels que alim=f*alim_star     
1290      do l=1,klev
1291         do ig=1,ngrid
1292            alim_star(ig,l)=0.
1293            alim(ig,l)=0.
1294         enddo
1295      enddo
1296c determination de la longueur de la couche d entrainement
1297      do ig=1,ngrid
1298         lentr(ig)=1
1299      enddo
1300
1301con ne considere que les premieres couches instables
1302      therm=.false.
1303      do k=nlay-2,1,-1
1304         do ig=1,ngrid
1305            if (ztv(ig,k).gt.ztv(ig,k+1).and.
1306     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
1307               lentr(ig)=k+1
1308               therm=.true.
1309            endif
1310          enddo
1311      enddo
1312c
1313c determination du lmin: couche d ou provient le thermique
1314      do ig=1,ngrid
1315         lmin(ig)=1
1316      enddo
1317      do ig=1,ngrid
1318         do l=nlay,2,-1
1319            if (ztv(ig,l-1).gt.ztv(ig,l)) then
1320               lmin(ig)=l-1
1321            endif
1322         enddo
1323      enddo
1324c
1325c definition de l'entrainement des couches
1326      do l=1,klev-1
1327         do ig=1,ngrid
1328            if (ztv(ig,l).gt.ztv(ig,l+1).and.
1329     s          l.ge.lmin(ig).and.l.lt.lentr(ig)) then
1330cdef possibles pour alim_star: zdthetadz, dthetadz, zdtheta
1331             alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)
1332c     s                       *(zlev(ig,l+1)-zlev(ig,l))
1333     s                       *sqrt(zlev(ig,l+1))
1334c             alim_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
1335c     s                         /zlev(ig,lentr(ig)+2)))**(3./2.)
1336            endif
1337         enddo
1338      enddo
1339     
1340c pas de thermique si couche 1 stable
1341      do ig=1,ngrid
1342c         if (lmin(ig).gt.1) then
1343cCRnouveau test
1344        if (alim_star(ig,1).lt.1.e-10) then
1345            do l=1,klev
1346                alim_star(ig,l)=0.
1347            enddo
1348         endif
1349      enddo
1350c calcul de l entrainement total
1351      do ig=1,ngrid
1352         alim_star_tot(ig)=0.
1353         entr_star_tot(ig)=0.
1354         detr_star_tot(ig)=0.
1355      enddo
1356      do ig=1,ngrid
1357         do k=1,klev
1358            alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,k)
1359         enddo
1360      enddo
1361c
1362c Calcul entrainement normalise
1363      do ig=1,ngrid
1364         if (alim_star_tot(ig).gt.1.e-10) then
1365c         do l=1,lentr(ig)
1366          do l=1,klev
1367cdef possibles pour entr_star: zdthetadz, dthetadz, zdtheta
1368            alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
1369         enddo
1370         endif
1371      enddo
1372       
1373c     print*,'fin calcul alim_star'
1374
1375cAM:initialisations
1376      do k=1,nlay
1377         do ig=1,ngrid
1378            ztva(ig,k)=ztv(ig,k)
1379            ztla(ig,k)=zthl(ig,k)
1380            zqla(ig,k)=0.
1381            zqta(ig,k)=po(ig,k)
1382            Zsat(ig) =.false.
1383         enddo
1384      enddo
1385      do k=1,klev
1386        do ig=1,ngrid
1387           detr_star(ig,k)=0.
1388           entr_star(ig,k)=0.
1389           detr(ig,k)=0.
1390           entr(ig,k)=0.
1391        enddo
1392      enddo
1393c     print*,'7 OK convect8'
1394      do k=1,klev+1
1395         do ig=1,ngrid
1396            zw2(ig,k)=0.
1397            fmc(ig,k)=0.
1398cCR
1399            f_star(ig,k)=0.
1400cRC
1401            larg_cons(ig,k)=0.
1402            larg_detr(ig,k)=0.
1403            wa_moy(ig,k)=0.
1404         enddo
1405      enddo
1406
1407cn     print*,'8 OK convect8'
1408      do ig=1,ngrid
1409         linter(ig)=1.
1410         lmaxa(ig)=1
1411         lmix(ig)=1
1412         wmaxa(ig)=0.
1413      enddo
1414
1415      nu_min=l_mix
1416      nu_max=1000.
1417c      do ig=1,ngrid
1418c      nu_max=wmax_sec(ig)
1419c      enddo
1420      do ig=1,ngrid
1421         do k=1,klev
1422            nu(ig,k)=0.
1423            nu_e(ig,k)=0.
1424         enddo
1425      enddo
1426cCalcul de l'excès de température du à la diffusion turbulente
1427      do ig=1,ngrid
1428         do l=1,klev
1429            dtheta(ig,l)=0.
1430         enddo
1431       enddo
1432      do ig=1,ngrid
1433         do l=1,lentr(ig)-1
1434      dtheta(ig,l)=sqrt(10.*0.4*zlev(ig,l+1)**2*1.
1435     s          *((ztv(ig,l+1)-ztv(ig,l))/(zlev(ig,l+1)-zlev(ig,l)))**2)
1436         enddo
1437       enddo
1438c      do l=1,nlay-2
1439      do l=1,klev-1
1440         do ig=1,ngrid
1441            if (ztv(ig,l).gt.ztv(ig,l+1)
1442     s         .and.alim_star(ig,l).gt.1.e-10
1443     s         .and.zw2(ig,l).lt.1e-10) then
1444cAM
1445ctest:on rajoute un excès de T dans couche alim
1446c               ztla(ig,l)=zthl(ig,l)+dtheta(ig,l)
1447               ztla(ig,l)=zthl(ig,l)
1448ctest: on rajoute un excès de q dans la couche alim
1449c               zqta(ig,l)=po(ig,l)+0.001
1450               zqta(ig,l)=po(ig,l)
1451               zqla(ig,l)=zl(ig,l)
1452cAM
1453               f_star(ig,l+1)=alim_star(ig,l)
1454ctest:calcul de dteta
1455               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
1456     s                     *(zlev(ig,l+1)-zlev(ig,l))
1457     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1458               w_est(ig,l+1)=zw2(ig,l+1)
1459               larg_detr(ig,l)=0.
1460c     print*,'coucou boucle 1'
1461            else if ((zw2(ig,l).ge.1e-10).and.
1462     s         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
1463c     print*,'coucou boucle 2'
1464cestimation du detrainement a partir de la geometrie du pas precedent
1465      if ((test(ig).eq.1).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
1466                  detr_star(ig,l)=0.
1467                  entr_star(ig,l)=0.
1468c     print*,'coucou test(ig)',test(ig),f0(ig),zmax0(ig)
1469             else
1470c     print*,'coucou debut detr'
1471ctests sur la definition du detr
1472        if (zqla(ig,l-1).gt.1.e-10) then
1473           nuage=.true.
1474        endif
1475
1476             w_est(ig,l+1)=zw2(ig,l)*
1477     s                   ((f_star(ig,l))**2)
1478     s                   /(f_star(ig,l)+alim_star(ig,l))**2+
1479     s                   2.*RG*(ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l)
1480     s                   *(zlev(ig,l+1)-zlev(ig,l))
1481             if (w_est(ig,l+1).lt.0.) then
1482                w_est(ig,l+1)=zw2(ig,l)
1483             endif
1484      if (l.gt.2) then
1485             if ((w_est(ig,l+1).gt.w_est(ig,l)).and.
1486     s           (zlev(ig,l+1).lt.zmax_sec(ig)).and.
1487     s            (zqla(ig,l-1).lt.1.e-10)) then
1488             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)
1489     s                *sqrt(w_est(ig,l+1))*sqrt(nu(ig,l)*zlev(ig,l+1))
1490     s       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(nu(ig,l)*zlev(ig,l)))
1491     s       /(r_aspect*zmax_sec(ig)))
1492             else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.
1493     s                (zqla(ig,l-1).lt.1.e-10)) then
1494       detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))
1495     s /(rhobarz(ig,lmix(ig))*wmaxa(ig))*
1496     s (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))
1497     s *((zmax_sec(ig)-zlev(ig,l+1))/((zmax_sec(ig)-zlev(ig,lmix(ig)))))
1498     s **2.
1499     s -rhobarz(ig,l)*sqrt(w_est(ig,l))
1500     s *((zmax_sec(ig)-zlev(ig,l))/((zmax_sec(ig)-zlev(ig,lmix(ig)))))
1501     s **2.)
1502             else
1503       detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)
1504     s                *(zlev(ig,l+1)-zlev(ig,l))
1505             
1506             endif
1507        else
1508        detr_star(ig,l)=0.
1509        endif
1510       
1511         detr_star(ig,l)=detr_star(ig,l)/f0(ig)
1512         if (nuage) then
1513         entr_star(ig,l)=0.4*detr_star(ig,l)
1514         else
1515         entr_star(ig,l)=0.4*detr_star(ig,l)
1516         endif
1517
1518             if ((detr_star(ig,l)).gt.f_star(ig,l)) then
1519              detr_star(ig,l)=f_star(ig,l)
1520c              entr_star(ig,l)=0.
1521              endif
1522
1523             if ((l.lt.lentr(ig))) then
1524                 entr_star(ig,l)=0.
1525c                 detr_star(ig,l)=0.
1526             endif 
1527
1528c           print*,'ok detr_star'
1529      endif
1530cprise en compte du detrainement dans le calcul du flux
1531             f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
1532     s                      -detr_star(ig,l)
1533ctest
1534c             if (f_star(ig,l+1).lt.0.) then
1535c                f_star(ig,l+1)=0.
1536c                entr_star(ig,l)=0.
1537c                detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)
1538c             endif
1539ctest sur le signe de f_star
1540       if (f_star(ig,l+1).gt.1.e-10) then
1541c                 then
1542ctest
1543c         if (((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10)) then
1544cAM on melange Tl et qt du thermique
1545con rajoute un excès de T dans la couche alim
1546c               if (l.lt.lentr(ig)) then
1547c           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1548c     s     (alim_star(ig,l)+entr_star(ig,l))*(zthl(ig,l)+dtheta(ig,l)))
1549c     s     /(f_star(ig,l+1)+detr_star(ig,l))
1550c               else
1551               ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1552     s                    (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))
1553     s                 /(f_star(ig,l+1)+detr_star(ig,l))
1554c     s                    /(f_star(ig,l+1))
1555c               endif
1556con rajoute un excès de q dans la couche alim
1557c               if (l.lt.lentr(ig)) then
1558c               zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1559c     s           (alim_star(ig,l)+entr_star(ig,l))*(po(ig,l)+0.001))
1560c     s                 /(f_star(ig,l+1)+detr_star(ig,l))
1561c               else
1562               zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1563     s                    (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))
1564     s                 /(f_star(ig,l+1)+detr_star(ig,l))
1565c     s                   /(f_star(ig,l+1))
1566c               endif
1567cAM on en deduit thetav et ql du thermique
1568cCR test
1569c               Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1570               Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1571               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1572               qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l) 
1573               qsatbef(ig)=MIN(0.5,qsatbef(ig))
1574               zcor=1./(1.-retv*qsatbef(ig))
1575               qsatbef(ig)=qsatbef(ig)*zcor
1576             Zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig)) .gt. 1.e-10)
1577
1578           if (Zsat(ig).and.(1.eq.1)) then
1579              qlbef=max(0.,zqta(ig,l)-qsatbef(ig))
1580              DT = 0.5*RLvCp*qlbef
1581c             write(17,*)'DT0=',DT
1582              do while (abs(DT).gt.DDT0)
1583c                 print*,'aie'
1584                 Tbef(ig)=Tbef(ig)+DT
1585                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1586                 qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
1587                 qsatbef(ig)=MIN(0.5,qsatbef(ig))
1588                 zcor=1./(1.-retv*qsatbef(ig))
1589                 qsatbef(ig)=qsatbef(ig)*zcor
1590                 qlbef=zqta(ig,l)-qsatbef(ig)
1591
1592                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1593                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
1594                 zcor=1./(1.-retv*qsatbef(ig))
1595                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
1596                 num=-Tbef(ig)+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
1597                 denom=1.+RLvCp*dqsat_dT
1598                 if (denom.lt.1.e-10) then
1599                    print*,'pb denom'
1600                 endif
1601                 DT=num/denom
1602c                 write(17,*)'DT=',DT
1603              enddo
1604              zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1605              zqla(ig,l) = max(0.,qlbef)
1606c              zqla(ig,l)=0.
1607             endif
1608c             zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1609c     
1610c on ecrit de maniere conservative (sat ou non)
1611c          T = Tl +Lv/Cp ql
1612cCR rq utilisation de humidite specifique ou rapport de melange?
1613           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1614           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1615con rajoute le calcul de zha pour diagnostiques (temp potentielle)
1616           zha(ig,l) = ztva(ig,l)
1617c           if (l.lt.lentr(ig)) then
1618c           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1619c     s              -zqla(ig,l))-zqla(ig,l)) + 0.1
1620c           else
1621           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1622     s              -zqla(ig,l))-zqla(ig,l))
1623c           endif
1624c           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1625c     s                 /(1.-retv*zqla(ig,l))
1626c           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1627c           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1628c     s                 /(1.-retv*zqta(ig,l))
1629c     s              -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1630c     s              -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1631c       write(13,*)zqla(ig,l),zqla(ig,l)/(1.-retv*zqla(ig,l))
1632con ecrit zqsat
1633           zqsatth(ig,l)=qsatbef(ig) 
1634c        enddo
1635c        DO ig=1,ngrid
1636c           if (zw2(ig,l).ge.1.e-10.and.
1637c     s               f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
1638c  mise a jour de la vitesse ascendante (l'air entraine de la couche
1639c  consideree commence avec une vitesse nulle).
1640c
1641c            if (f_star(ig,l+1).gt.1.e-10) then
1642            zw2(ig,l+1)=zw2(ig,l)*
1643c     s                  ((f_star(ig,l)-detr_star(ig,l))**2)
1644c     s                  /f_star(ig,l+1)**2+
1645     s                   ((f_star(ig,l))**2)
1646     s                   /(f_star(ig,l+1)+detr_star(ig,l))**2+
1647c     s                    /(f_star(ig,l+1))**2+           
1648     s                   2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1649     s                   *(zlev(ig,l+1)-zlev(ig,l))
1650c     s                   *(f_star(ig,l)/f_star(ig,l+1))**2
1651
1652            endif
1653        endif
1654c
1655            if (zw2(ig,l+1).lt.0.) then
1656               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
1657     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1658               zw2(ig,l+1)=0.
1659c              print*,'linter=',linter(ig)
1660c          else if ((zw2(ig,l+1).lt.1.e-10).and.(zw2(ig,l+1).ge.0.)) then
1661c               linter(ig)=l+1
1662c               print*,'linter=l',zw2(ig,l),zw2(ig,l+1)
1663            else
1664               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1665c            wa_moy(ig,l+1)=zw2(ig,l+1)
1666            endif
1667            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1668c   lmix est le niveau de la couche ou w (wa_moy) est maximum
1669               lmix(ig)=l+1
1670               wmaxa(ig)=wa_moy(ig,l+1)
1671            endif
1672         enddo
1673      enddo
1674      print*,'fin calcul zw2'
1675c
1676c Calcul de la couche correspondant a la hauteur du thermique
1677      do ig=1,ngrid
1678         lmax(ig)=lentr(ig)
1679      enddo
1680      do ig=1,ngrid
1681         do l=nlay,lentr(ig)+1,-1
1682            if (zw2(ig,l).le.1.e-10) then
1683               lmax(ig)=l-1
1684            endif
1685         enddo
1686      enddo
1687c pas de thermique si couche 1 stable
1688      do ig=1,ngrid
1689         if (lmin(ig).gt.1) then
1690            lmax(ig)=1
1691             lmin(ig)=1
1692             lentr(ig)=1
1693         endif
1694      enddo
1695c   
1696c Determination de zw2 max
1697      do ig=1,ngrid
1698         wmax(ig)=0.
1699      enddo
1700
1701      do l=1,nlay
1702         do ig=1,ngrid
1703            if (l.le.lmax(ig)) then
1704                if (zw2(ig,l).lt.0.)then
1705                  print*,'pb2 zw2<0'
1706                endif
1707                zw2(ig,l)=sqrt(zw2(ig,l))
1708                wmax(ig)=max(wmax(ig),zw2(ig,l))
1709            else
1710                 zw2(ig,l)=0.
1711            endif
1712          enddo
1713      enddo
1714
1715c   Longueur caracteristique correspondant a la hauteur des thermiques.
1716      do  ig=1,ngrid
1717         zmax(ig)=0.
1718         zlevinter(ig)=zlev(ig,1)
1719      enddo
1720      do  ig=1,ngrid
1721c calcul de zlevinter
1722          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
1723     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
1724     s    -zlev(ig,lmax(ig)))
1725cpour le cas ou on prend tjs lmin=1
1726c       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1727       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
1728       zmax0(ig)=zmax(ig)
1729       write(11,*)'ig,lmax,linter',ig,lmax(ig),linter(ig)
1730       write(12,*)'ig,zlevinter,zmax',ig,zmax(ig),zlevinter(ig)
1731      enddo
1732
1733cCalcul de zmax_sec et wmax_sec
1734      call fermeture_seche(ngrid,nlay
1735     s                  ,pplay,pplev,pphi,zlev,rhobarz,f0,zpspsk
1736     s                  ,alim,zh,zo,lentr,lmin,nu_min,nu_max,r_aspect
1737     s                  ,zmax_sec2,wmax_sec2)
1738
1739      print*,'avant fermeture'
1740c Fermeture,determination de f
1741c en lmax f=d-e
1742      do ig=1,ngrid
1743c      entr_star(ig,lmax(ig))=0.
1744c      f_star(ig,lmax(ig)+1)=0.
1745c      detr_star(ig,lmax(ig))=f_star(ig,lmax(ig))+entr_star(ig,lmax(ig))
1746c     s                       +alim_star(ig,lmax(ig))
1747      enddo
1748c
1749      do ig=1,ngrid
1750         alim_star2(ig)=0.
1751      enddo
1752ccalcul de entr_star_tot
1753      do ig=1,ngrid
1754         do k=1,lmix(ig)
1755            entr_star_tot(ig)=entr_star_tot(ig)
1756c     s                        +entr_star(ig,k)
1757     s                        +alim_star(ig,k)
1758c     s                        -detr_star(ig,k)
1759            detr_star_tot(ig)=detr_star_tot(ig)
1760c     s                        +alim_star(ig,k)
1761     s                        -detr_star(ig,k)
1762     s                        +entr_star(ig,k)
1763         enddo
1764      enddo
1765     
1766      do ig=1,ngrid
1767         if (alim_star_tot(ig).LT.1.e-10) then
1768            f(ig)=0.
1769         else   
1770c             do k=lmin(ig),lentr(ig)
1771             do k=1,lentr(ig)
1772                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2
1773     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
1774             enddo
1775             if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then
1776             f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect
1777     s             *alim_star2(ig))
1778             f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/
1779     s                     zmax_sec(ig))*wmax_sec(ig))
1780             else
1781             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect*alim_star2(ig))
1782            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/
1783     s                     zmax(ig))*wmax(ig))
1784             endif
1785         endif
1786         f0(ig)=f(ig)
1787      enddo
1788      print*,'apres fermeture'
1789c Calcul de l'entrainement
1790         do ig=1,ngrid
1791            do k=1,klev
1792            alim(ig,k)=f(ig)*alim_star(ig,k)
1793         enddo
1794      enddo
1795cCR:test pour entrainer moins que la masse
1796c       do ig=1,ngrid
1797c          do l=1,lentr(ig)
1798c             if ((alim(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
1799c                alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
1800c     s                       -0.9*masse(ig,l)/ptimestep
1801c                alim(ig,l)=0.9*masse(ig,l)/ptimestep
1802c             endif
1803c          enddo
1804c       enddo
1805c calcul du détrainement
1806         do ig=1,klon
1807             do k=1,klev
1808            detr(ig,k)=f(ig)*detr_star(ig,k)
1809            if (detr(ig,k).lt.0.) then
1810c               print*,'detr1<0!!!'
1811            endif
1812            enddo
1813            do k=1,klev
1814            entr(ig,k)=f(ig)*entr_star(ig,k)
1815            if (entr(ig,k).lt.0.) then
1816c               print*,'entr1<0!!!'
1817            endif
1818         enddo
1819      enddo
1820c
1821c       do ig=1,ngrid
1822c          do l=1,klev
1823c          if (((detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep).gt.
1824c     s          (masse(ig,l))) then 
1825c      print*,'d2+e2+a2>m2','ig=',ig,'l=',l,'lmax(ig)=',lmax(ig),'d+e+a='
1826c     s,(detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep,'m=',masse(ig,l)
1827c      endif
1828c      enddo
1829c      enddo
1830c Calcul des flux
1831
1832      do ig=1,ngrid
1833         do l=1,lmax(ig)
1834c         do l=1,klev
1835c             fmc(ig,l+1)=f(ig)*f_star(ig,l+1)
1836            fmc(ig,l+1)=fmc(ig,l)+alim(ig,l)+entr(ig,l)-detr(ig,l)
1837c        print*,'??!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1838c     s  'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1839c     s  'f+1=',fmc(ig,l+1)
1840          if (fmc(ig,l+1).lt.0.) then
1841               print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1842               fmc(ig,l+1)=fmc(ig,l)
1843               detr(ig,l)=alim(ig,l)+entr(ig,l)
1844c               fmc(ig,l+1)=0.
1845c               print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1846            endif
1847c       if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1848c          f_old=fmc(ig,l+1)
1849c          fmc(ig,l+1)=fmc(ig,l)
1850c          detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1851c       endif
1852       
1853c        if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1854c          f_old=fmc(ig,l+1)
1855c          fmc(ig,l+1)=fmc(ig,l)
1856c          detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l)
1857c       endif
1858crajout du test sur alpha croissant
1859cif test
1860c       if (1.eq.0) then
1861
1862       if (l.eq.klev) then
1863          print*,'THERMCELL PB ig=',ig,'   l=',l
1864          abort_message = 'THERMCELL PB'
1865          CALL abort_gcm (modname,abort_message,1)
1866       endif
1867!       if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1868!     s     (l.ge.lentr(ig)).and.
1869       if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1870     s         (l.ge.lentr(ig)) ) then
1871          if ( ((fmc(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.
1872     s     (fmc(ig,l)/(rhobarz(ig,l)*zw2(ig,l))))) then
1873           f_old=fmc(ig,l+1)
1874           fmc(ig,l+1)=fmc(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)
1875     s                          /(rhobarz(ig,l)*zw2(ig,l))
1876           detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1877c           detr(ig,l)=(fmc(ig,l+1)-fmc(ig,l))/(0.4-1.)
1878c           entr(ig,l)=0.4*detr(ig,l)
1879c           entr(ig,l)=fmc(ig,l+1)-fmc(ig,l)+detr(ig,l)
1880        endif
1881        endif
1882        if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1883          f_old=fmc(ig,l+1)
1884          fmc(ig,l+1)=fmc(ig,l)
1885          detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1886       endif
1887       if (detr(ig,l).gt.fmc(ig,l)) then
1888               detr(ig,l)=fmc(ig,l)
1889               entr(ig,l)=fmc(ig,l+1)-alim(ig,l)
1890        endif
1891       if (fmc(ig,l+1).lt.0.) then
1892               detr(ig,l)=detr(ig,l)+fmc(ig,l+1)
1893               fmc(ig,l+1)=0.
1894               print*,'fmc2<0',l+1,lmax(ig)
1895            endif
1896           
1897ctest pour ne pas avoir f=0 et d=e/=0
1898c       if (fmc(ig,l+1).lt.1.e-10) then
1899c          detr(ig,l+1)=0.
1900c          entr(ig,l+1)=0.
1901c          zqla(ig,l+1)=0.
1902c          zw2(ig,l+1)=0.
1903c          lmax(ig)=l+1
1904c          zmax(ig)=zlev(ig,lmax(ig))
1905c       endif
1906        if (zw2(ig,l+1).gt.1.e-10) then
1907       if ((((fmc(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.
1908     s      1.)) then
1909          f_old=fmc(ig,l+1)
1910          fmc(ig,l+1)=rhobarz(ig,l+1)*zw2(ig,l+1)
1911          zw2(ig,l+1)=0.
1912          zqla(ig,l+1)=0.
1913          detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1914         lmax(ig)=l+1
1915          zmax(ig)=zlev(ig,lmax(ig))
1916          print*,'alpha>1',l+1,lmax(ig)
1917       endif
1918        endif
1919c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
1920cendif test
1921c      endif
1922      enddo
1923      enddo
1924      do ig=1,ngrid
1925c         if (fmc(ig,lmax(ig)+1).ne.0.) then
1926         fmc(ig,lmax(ig)+1)=0.
1927         entr(ig,lmax(ig))=0.
1928         detr(ig,lmax(ig))=fmc(ig,lmax(ig))+entr(ig,lmax(ig))
1929     s                     +alim(ig,lmax(ig))
1930c         endif
1931      enddo
1932ctest sur le signe de fmc
1933       do ig=1,ngrid
1934         do l=1,klev+1
1935            if (fmc(ig,l).lt.0.) then
1936         print*,'fm1<0!!!','ig=',ig,'l=',l,'a=',alim(ig,l-1),'e='
1937     s ,entr(ig,l-1),'f=',fmc(ig,l-1),'d=',detr(ig,l-1),'f+1=',fmc(ig,l)
1938            endif
1939         enddo
1940       enddo
1941ctest de verification
1942      do ig=1,ngrid
1943       do l=1,lmax(ig)
1944       if ((abs(fmc(ig,l+1)-fmc(ig,l)-alim(ig,l)-entr(ig,l)+detr(ig,l)))
1945     s           .gt.1.e-4) then
1946c      print*,'pbcm!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1947c     s  'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1948c     s  'f+1=',fmc(ig,l+1)
1949       endif
1950       if (detr(ig,l).lt.0.) then
1951          print*,'detrdemi<0!!!'
1952       endif
1953         enddo
1954      enddo
1955c
1956cRC
1957cCR def de  zmix continu (profil parabolique des vitesses)
1958      do ig=1,ngrid
1959           if (lmix(ig).gt.1.) then
1960c test
1961              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
1962     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
1963     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
1964     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
1965     s        then
1966c             
1967            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
1968     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
1969     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
1970     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
1971     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
1972     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
1973     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
1974     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
1975              else
1976              zmix(ig)=zlev(ig,lmix(ig))
1977              print*,'pb zmix'
1978              endif
1979          else
1980              zmix(ig)=0.
1981          endif
1982ctest
1983         if ((zmax(ig)-zmix(ig)).le.0.) then
1984            zmix(ig)=0.9*zmax(ig)
1985c            print*,'pb zmix>zmax'
1986         endif
1987      enddo
1988      do ig=1,klon
1989         zmix0(ig)=zmix(ig)
1990      enddo
1991c
1992c calcul du nouveau lmix correspondant
1993      do ig=1,ngrid
1994         do l=1,klev
1995            if (zmix(ig).ge.zlev(ig,l).and.
1996     s          zmix(ig).lt.zlev(ig,l+1)) then
1997              lmix(ig)=l
1998             endif
1999          enddo
2000      enddo
2001c
2002cne devrait pas arriver!!!!!
2003      do ig=1,ngrid
2004       do l=1,klev
2005          if (detr(ig,l).gt.(fmc(ig,l)+alim(ig,l))+entr(ig,l)) then
2006             print*,'detr2>fmc2!!!','ig=',ig,'l=',l,'d=',detr(ig,l),
2007     s             'f=',fmc(ig,l),'lmax=',lmax(ig)
2008c             detr(ig,l)=fmc(ig,l)+alim(ig,l)+entr(ig,l)
2009c             entr(ig,l)=0.
2010c             fmc(ig,l+1)=0.
2011c             zw2(ig,l+1)=0.     
2012c             zqla(ig,l+1)=0.
2013           print*,'pb!fm=0 et f_star>0',l,lmax(ig)       
2014c             lmax(ig)=l
2015          endif
2016       enddo
2017      enddo
2018      do ig=1,ngrid
2019         do l=lmax(ig)+1,klev+1
2020c            fmc(ig,l)=0.
2021c            detr(ig,l)=0.
2022c            entr(ig,l)=0.
2023c            zw2(ig,l)=0.
2024c            zqla(ig,l)=0.
2025         enddo
2026      enddo
2027
2028cCalcul du detrainement lors du premier passage
2029c     print*,'9 OK convect8'
2030c     print*,'WA1 ',wa_moy
2031
2032c   determination de l'indice du debut de la mixed layer ou w decroit
2033
2034c   calcul de la largeur de chaque ascendance dans le cas conservatif.
2035c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
2036c   d'une couche est égale à la hauteur de la couche alimentante.
2037c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
2038c   de la vitesse d'entrainement horizontal dans la couche alimentante.
2039
2040      do l=2,nlay
2041         do ig=1,ngrid
2042            if (l.le.lmax(ig).and.(test(ig).eq.1)) then
2043               zw=max(wa_moy(ig,l),1.e-10)
2044               larg_cons(ig,l)=zmax(ig)*r_aspect
2045     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
2046            endif
2047         enddo
2048      enddo
2049
2050      do l=2,nlay
2051         do ig=1,ngrid
2052            if (l.le.lmax(ig).and.(test(ig).eq.1)) then
2053c              if (idetr.eq.0) then
2054c  cette option est finalement en dur.
2055                  if ((l_mix*zlev(ig,l)).lt.0.)then
2056                   print*,'pb l_mix*zlev<0'
2057                  endif
2058cCR: test: nouvelle def de lambda
2059c                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2060                  if (zw2(ig,l).gt.1.e-10) then
2061                  larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
2062                  else
2063                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2064                  endif
2065c              else if (idetr.eq.1) then
2066c                 larg_detr(ig,l)=larg_cons(ig,l)
2067c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
2068c              else if (idetr.eq.2) then
2069c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2070c    s            *sqrt(wa_moy(ig,l))
2071c              else if (idetr.eq.4) then
2072c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2073c    s            *wa_moy(ig,l)
2074c              endif
2075         endif
2076         enddo
2077       enddo
2078
2079c     print*,'10 OK convect8'
2080c     print*,'WA2 ',wa_moy
2081c   cal1cul de la fraction de la maille concernée par l'ascendance en tenant
2082c   compte de l'epluchage du thermique.
2083c
2084c
2085      do l=2,nlay
2086         do ig=1,ngrid
2087            if(larg_cons(ig,l).gt.1..and.(test(ig).eq.1)) then
2088c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
2089               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
2090     s            /(r_aspect*zmax(ig))
2091c test
2092               fraca(ig,l)=max(fraca(ig,l),0.)
2093               fraca(ig,l)=min(fraca(ig,l),0.5)
2094               fracd(ig,l)=1.-fraca(ig,l)
2095               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
2096            else
2097c              wa_moy(ig,l)=0.
2098               fraca(ig,l)=0.
2099               fracc(ig,l)=0.
2100               fracd(ig,l)=1.
2101            endif
2102         enddo
2103      enddo                 
2104cCR: calcul de fracazmix
2105       do ig=1,ngrid
2106          if (test(ig).eq.1) then
2107          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
2108     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
2109     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
2110     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
2111          endif
2112       enddo
2113c
2114       do l=2,nlay
2115          do ig=1,ngrid
2116             if(larg_cons(ig,l).gt.1..and.(test(ig).eq.1)) then
2117               if (l.gt.lmix(ig)) then
2118ctest
2119                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
2120c                   print*,'pb xxx'
2121                    xxx(ig,l)=(lmax(ig)+1.-l)/(lmax(ig)+1.-lmix(ig))
2122                 else
2123                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
2124                 endif
2125           if (idetr.eq.0) then
2126               fraca(ig,l)=fracazmix(ig)
2127           else if (idetr.eq.1) then
2128               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
2129           else if (idetr.eq.2) then
2130               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
2131           else
2132               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
2133           endif
2134c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
2135               fraca(ig,l)=max(fraca(ig,l),0.)
2136               fraca(ig,l)=min(fraca(ig,l),0.5)
2137               fracd(ig,l)=1.-fraca(ig,l)
2138               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
2139             endif
2140            endif
2141         enddo
2142      enddo
2143
2144      print*,'fin calcul fraca'
2145c     print*,'11 OK convect8'
2146c     print*,'Ea3 ',wa_moy
2147c------------------------------------------------------------------
2148c   Calcul de fracd, wd
2149c   somme wa - wd = 0
2150c------------------------------------------------------------------
2151
2152
2153      do ig=1,ngrid
2154         fm(ig,1)=0.
2155         fm(ig,nlay+1)=0.
2156      enddo
2157
2158      do l=2,nlay
2159           do ig=1,ngrid
2160              if (test(ig).eq.1) then
2161              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
2162cCR:test
2163              if (alim(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
2164     s            .and.l.gt.lmix(ig)) then
2165                 fm(ig,l)=fm(ig,l-1)
2166c                 write(1,*)'ajustement fm, l',l
2167              endif
2168c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
2169cRC
2170              endif
2171           enddo
2172         do ig=1,ngrid
2173            if(fracd(ig,l).lt.0.1.and.(test(ig).eq.1)) then
2174              abort_message = 'fracd trop petit'
2175              CALL abort_gcm (modname,abort_message,1)
2176            else
2177c    vitesse descendante "diagnostique"
2178               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
2179            endif
2180         enddo
2181      enddo
2182
2183      do l=1,nlay+1
2184         do ig=1,ngrid
2185            if (test(ig).eq.0) then
2186              fm(ig,l)=fmc(ig,l)
2187            endif
2188         enddo
2189      enddo
2190   
2191cfin du first
2192      do l=1,nlay
2193         do ig=1,ngrid
2194c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
2195            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
2196         enddo
2197      enddo
2198
2199      print*,'12 OK convect8'
2200c     print*,'WA4 ',wa_moy
2201cc------------------------------------------------------------------
2202c   calcul du transport vertical
2203c------------------------------------------------------------------
2204
2205      go to 4444
2206c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
2207      do l=2,nlay-1
2208         do ig=1,ngrid
2209            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
2210     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
2211      print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
2212     s         ,fm(ig,l+1)*ptimestep
2213     s         ,'   M=',masse(ig,l),masse(ig,l+1)
2214             endif
2215         enddo
2216      enddo
2217
2218      do l=1,nlay
2219         do ig=1,ngrid
2220            if((alim(ig,l)+entr(ig,l))*ptimestep.gt.masse(ig,l)) then
2221      print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
2222     s         ,(entr(ig,l)+alim(ig,l))*ptimestep
2223     s         ,'   M=',masse(ig,l)
2224             endif
2225         enddo
2226      enddo
2227
2228      do l=1,nlay
2229         do ig=1,ngrid
2230            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
2231c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
2232c    s         ,'   FM=',fm(ig,l)
2233            endif
2234            if(.not.masse(ig,l).ge.1.e-10
2235     s         .or..not.masse(ig,l).le.1.e4) then
2236c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
2237c    s         ,'   M=',masse(ig,l)
2238c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
2239c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
2240c     print*,'zlev(ig,l+1),zlev(ig,l)'
2241c    s                ,zlev(ig,l+1),zlev(ig,l)
2242c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
2243c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
2244            endif
2245            if(.not.alim(ig,l).ge.0..or..not.alim(ig,l).le.10.) then
2246c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
2247c    s         ,'   E=',entr(ig,l)
2248            endif
2249         enddo
2250      enddo
2251
22524444   continue
2253
2254cCR:redefinition du entr
2255cCR:test:on ne change pas la def du entr mais la def du fm
2256       do l=1,nlay
2257         do ig=1,ngrid
2258            if (test(ig).eq.1) then
2259            detr(ig,l)=fm(ig,l)+alim(ig,l)-fm(ig,l+1)
2260            if (detr(ig,l).lt.0.) then
2261c                entr(ig,l)=entr(ig,l)-detr(ig,l)
2262                fm(ig,l+1)=fm(ig,l)+alim(ig,l)
2263                detr(ig,l)=0.
2264c                write(11,*)'l,ig,entr',l,ig,entr(ig,l)
2265c     print*,'WARNING !!! detrainement negatif ',ig,l
2266            endif
2267            endif
2268         enddo
2269      enddo
2270cRC
2271
2272      if (w2di.eq.1) then
2273         fm0=fm0+ptimestep*(fm-fm0)/tho
2274         entr0=entr0+ptimestep*(alim+entr-entr0)/tho
2275      else
2276         fm0=fm
2277         entr0=alim+entr
2278         detr0=detr
2279         alim0=alim
2280c         zoa=zqta
2281c         entr0=alim
2282      endif
2283
2284      if (1.eq.1) then
2285c         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2286c     .    ,zh,zdhadj,zha)
2287c         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2288c     .    ,zo,pdoadj,zoa)
2289         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse,
2290     .                    zthl,zdthladj,zta)
2291         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse,
2292     .                   po,pdoadj,zoa)
2293      else
2294         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
2295     .    ,zh,zdhadj,zha)
2296         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
2297     .    ,zo,pdoadj,zoa)
2298      endif
2299
2300      if (1.eq.0) then
2301         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
2302     .    ,fraca,zmax
2303     .    ,zu,zv,pduadj,pdvadj,zua,zva)
2304      else
2305         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2306     .    ,zu,pduadj,zua)
2307         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2308     .    ,zv,pdvadj,zva)
2309      endif
2310
2311cCalcul des moments
2312c      do l=1,nlay
2313c         do ig=1,ngrid
2314c            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
2315c            zf2=zf/(1.-zf)
2316c            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
2317c            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
2318c         enddo
2319c      enddo
2320
2321
2322
2323
2324
2325
2326c     print*,'13 OK convect8'
2327c     print*,'WA5 ',wa_moy
2328      do l=1,nlay
2329         do ig=1,ngrid
2330c            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
2331           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
2332         enddo
2333      enddo
2334
2335
2336c     do l=1,nlay
2337c        do ig=1,ngrid
2338c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
2339c     print*,'WARN!!! ig=',ig,'  l=',l
2340c    s         ,'   pdtadj=',pdtadj(ig,l)
2341c           endif
2342c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
2343c     print*,'WARN!!! ig=',ig,'  l=',l
2344c    s         ,'   pdoadj=',pdoadj(ig,l)
2345c           endif
2346c        enddo
2347c      enddo
2348
2349      print*,'14 OK convect8'
2350c------------------------------------------------------------------
2351c   Calculs pour les sorties
2352c------------------------------------------------------------------
2353ccalcul de fraca pour les sorties
2354      do l=2,klev
2355         do ig=1,klon
2356            if (zw2(ig,l).gt.1.e-10) then
2357            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
2358            else
2359            fraca(ig,l)=0.
2360            endif
2361         enddo
2362      enddo
2363      if(sorties) then
2364      do l=1,nlay
2365         do ig=1,ngrid
2366            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
2367            zld(ig,l)=fracd(ig,l)*zmax(ig)
2368            if(1.-fracd(ig,l).gt.1.e-10)
2369     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
2370         enddo
2371      enddo
2372c CR calcul du niveau de condensation
2373c initialisation
2374      do ig=1,ngrid
2375         nivcon(ig)=0.
2376         zcon(ig)=0.
2377      enddo
2378      do k=nlay,1,-1
2379         do ig=1,ngrid
2380            if (zqla(ig,k).gt.1e-10) then
2381               nivcon(ig)=k
2382               zcon(ig)=zlev(ig,k)
2383            endif
2384c            if (zcon(ig).gt.1.e-10) then
2385c               nuage=.true.
2386c            else
2387c               nuage=.false.
2388c            endif
2389         enddo
2390      enddo
2391     
2392      do l=1,nlay
2393         do ig=1,ngrid
2394            zf=fraca(ig,l)
2395            zf2=zf/(1.-zf)
2396            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
2397            wth2(ig,l)=zf2*(zw2(ig,l))**2
2398c           print*,'wth2=',wth2(ig,l)
2399            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))
2400     s                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
2401            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
2402ctest: on calcul q2/po=ratqsc
2403c            if (nuage) then
2404            ratqscth(ig,l)=sqrt(q2(ig,l))/(po(ig,l)*1000.)
2405c            else
2406c            ratqscth(ig,l)=0.
2407c            endif
2408         enddo
2409      enddo
2410ccalcul du ratqscdiff
2411      sum=0.
2412      sumdiff=0.
2413      ratqsdiff(:,:)=0.
2414      do ig=1,ngrid
2415         do l=1,lentr(ig)
2416            sum=sum+alim_star(ig,l)*zqta(ig,l)*1000.
2417         enddo
2418      enddo
2419      do ig=1,ngrid
2420          do l=1,lentr(ig)
2421          zf=fraca(ig,l)
2422          zf2=zf/(1.-zf)
2423       sumdiff=sumdiff+alim_star(ig,l)
2424     s           *(zqta(ig,l)*1000.-sum)**2
2425c      ratqsdiff=ratqsdiff+alim_star(ig,l)*
2426c     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
2427          enddo
2428      enddo
2429      do l=1,klev
2430      do ig=1,ngrid
2431      ratqsdiff(ig,l)=sqrt(sumdiff)/(po(ig,l)*1000.)   
2432c      write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
2433      enddo
2434      enddo     
2435cdeja fait
2436c      do l=1,nlay
2437c         do ig=1,ngrid
2438c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
2439c            if (detr(ig,l).lt.0.) then
2440c                entr(ig,l)=entr(ig,l)-detr(ig,l)
2441c                detr(ig,l)=0.
2442c     print*,'WARNING !!! detrainement negatif ',ig,l
2443c            endif
2444c         enddo
2445c      enddo
2446
2447c     print*,'15 OK convect8'
2448
2449      isplit=isplit+1       
2450
2451
2452c #define und
2453        goto 123
2454#ifdef und
2455      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
2456      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
2457      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
2458      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
2459      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
2460      CALL writeg1d(1,nlay,zla,'la      ','la      ')
2461      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
2462      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
2463      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
2464      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
2465      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
2466      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
2467      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
2468      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
2469      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
2470      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
2471      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
2472      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
2473      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
2474      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
2475      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
2476      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
2477      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
2478      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
2479
2480      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
2481      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
2482      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
2483
2484c   recalcul des flux en diagnostique...
2485c     print*,'PAS DE TEMPS ',ptimestep
2486       call dt2F(pplev,pplay,pt,pdtadj,wh)
2487      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
2488#endif
2489123   continue
2490! #define troisD
2491#ifdef troisD
2492c       if (sorties) then
2493      print*,'Debut des wrgradsfi'
2494
2495c      print*,'16 OK convect8'
2496c         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
2497c         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
2498         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
2499         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
2500c         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
2501c         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
2502c      print*,'WA6 ',wa_moy
2503c         call wrgradsfi(1,nlay,zla,'la        ','la        ')
2504c         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
2505         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
2506         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
2507         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
2508         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
2509         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
2510         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
2511         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
2512         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
2513         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
2514         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
2515         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
2516         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
2517         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
2518         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
2519         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
2520         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
2521         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
2522         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
2523         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
2524         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
2525         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
2526         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
2527         call wrgradsfi(1,nlay,w_est,'w_est      ','w_est      ')
2528con sort les moments
2529         call wrgradsfi(1,nlay,thetath2,'zh2       ','zh2       ')
2530         call wrgradsfi(1,nlay,wth2,'w2       ','w2       ')
2531         call wrgradsfi(1,nlay,wth3,'w3       ','w3       ')
2532         call wrgradsfi(1,nlay,q2,'q2       ','q2       ')
2533         call wrgradsfi(1,nlay,dtheta,'dT       ','dT       ')
2534c
2535         call wrgradsfi(1,nlay,zw_sec,'zw_sec       ','zw_sec       ')
2536         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
2537         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
2538         call wrgradsfi(1,nlay,nu,'nu       ','nu       ')
2539
2540         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
2541         call wrgradsfi(1,nlay,zoa,'zoa        ','zoa        ')
2542c         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
2543c         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
2544
2545cAM:nouveaux diagnostiques
2546         call wrgradsfi(1,nlay,zthl,'zthl        ','zthl        ')
2547         call wrgradsfi(1,nlay,zta,'zta        ','zta        ')
2548         call wrgradsfi(1,nlay,zl,'zl        ','zl        ')
2549         call wrgradsfi(1,nlay,zdthladj,'zdthladj    ',
2550     s        'zdthladj    ')
2551         call wrgradsfi(1,nlay,ztla,'ztla      ','ztla      ')
2552         call wrgradsfi(1,nlay,zqta,'zqta      ','zqta      ')
2553         call wrgradsfi(1,nlay,zqla,'zqla      ','zqla      ')
2554         call wrgradsfi(1,nlay,deltaz,'deltaz      ','deltaz      ')
2555cCR:nouveaux diagnostiques
2556      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')
2557      call wrgradsfi(1,nlay,detr_star  ,'detr_star   ','detr_star   ')     
2558      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
2559      call wrgradsfi(1,nlay,zqsat    ,'zqsat   ','zqsat   ')
2560      call wrgradsfi(1,nlay,zqsatth    ,'qsath   ','qsath   ')
2561      call wrgradsfi(1,nlay,alim_star    ,'alim_star   ','alim_star   ')
2562      call wrgradsfi(1,nlay,alim    ,'alim   ','alim   ')
2563      call wrgradsfi(1,1,f,'f      ','f      ')
2564      call wrgradsfi(1,1,alim_star_tot,'a_s_t      ','a_s_t      ')
2565      call wrgradsfi(1,1,alim_star2,'a_2      ','a_2      ')
2566      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
2567      call wrgradsfi(1,1,zmax_sec,'z_sec      ','z_sec      ')
2568c      call wrgradsfi(1,1,zmax_sec2,'zz_se      ','zz_se      ')
2569      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
2570      call wrgradsfi(1,1,nivcon,'nivcon      ','nivcon      ')
2571      call wrgradsfi(1,1,zcon,'zcon      ','zcon      ')
2572      zsortie1d(:)=lmax(:)
2573      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
2574      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
2575      call wrgradsfi(1,1,wmax_sec,'w_sec      ','w_sec      ')
2576      zsortie1d(:)=lmix(:)
2577      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
2578      zsortie1d(:)=lentr(:)
2579      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
2580
2581c      print*,'17 OK convect8'
2582
2583         do k=1,klev/10
2584            write(str2,'(i2.2)') k
2585            str10='wa'//str2
2586            do l=1,nlay
2587               do ig=1,ngrid
2588                  zsortie(ig,l)=wa(ig,k,l)
2589               enddo
2590            enddo   
2591            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
2592            do l=1,nlay
2593               do ig=1,ngrid
2594                  zsortie(ig,l)=larg_part(ig,k,l)
2595               enddo
2596            enddo
2597           str10='la'//str2
2598            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
2599         enddo
2600
2601
2602c     print*,'18 OK convect8'
2603c      endif
2604      print*,'Fin des wrgradsfi'
2605#endif
2606
2607      endif
2608
2609c     if(wa_moy(1,4).gt.1.e-10) stop
2610
2611      print*,'19 OK convect8'
2612      return
2613      end
2614      SUBROUTINE thermcell_eau(ngrid,nlay,ptimestep
2615     s                  ,pplay,pplev,pphi
2616     s                  ,pu,pv,pt,po
2617     s                  ,pduadj,pdvadj,pdtadj,pdoadj
2618     s                  ,fm0,entr0
2619c    s                  ,pu_therm,pv_therm
2620     s                  ,r_aspect,l_mix,w2di,tho)
2621
2622      USE dimphy
2623      IMPLICIT NONE
2624
2625c=======================================================================
2626c
2627c   Calcul du transport verticale dans la couche limite en presence
2628c   de "thermiques" explicitement representes
2629c
2630c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
2631c
2632c   le thermique est supposé homogène et dissipé par mélange avec
2633c   son environnement. la longueur l_mix contrôle l'efficacité du
2634c   mélange
2635c
2636c   Le calcul du transport des différentes espèces se fait en prenant
2637c   en compte:
2638c     1. un flux de masse montant
2639c     2. un flux de masse descendant
2640c     3. un entrainement
2641c     4. un detrainement
2642c
2643c=======================================================================
2644
2645c-----------------------------------------------------------------------
2646c   declarations:
2647c   -------------
2648
2649#include "dimensions.h"
2650cccc#include "dimphy.h"
2651#include "YOMCST.h"
2652#include "YOETHF.h"
2653#include "FCTTRE.h"
2654
2655c   arguments:
2656c   ----------
2657
2658      INTEGER ngrid,nlay,w2di
2659      REAL tho
2660      real ptimestep,l_mix,r_aspect
2661      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
2662      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
2663      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
2664      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
2665      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
2666      real pphi(ngrid,nlay)
2667
2668      integer idetr
2669      save idetr
2670      data idetr/3/
2671c$OMP THREADPRIVATE(idetr)
2672
2673c   local:
2674c   ------
2675
2676      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
2677      real zsortie1d(klon)
2678c CR: on remplace lmax(klon,klev+1)
2679      INTEGER lmax(klon),lmin(klon),lentr(klon)
2680      real linter(klon)
2681      real zmix(klon), fracazmix(klon)
2682c RC
2683      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
2684
2685      real zlev(klon,klev+1),zlay(klon,klev)
2686      REAL zh(klon,klev),zdhadj(klon,klev)
2687      real zthl(klon,klev),zdthladj(klon,klev)
2688      REAL ztv(klon,klev)
2689      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
2690      real zl(klon,klev)
2691      REAL wh(klon,klev+1)
2692      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
2693      real zla(klon,klev+1)
2694      real zwa(klon,klev+1)
2695      real zld(klon,klev+1)
2696      real zwd(klon,klev+1)
2697      real zsortie(klon,klev)
2698      real zva(klon,klev)
2699      real zua(klon,klev)
2700      real zoa(klon,klev)
2701
2702      real zta(klon,klev)
2703      real zha(klon,klev)
2704      real wa_moy(klon,klev+1)
2705      real fraca(klon,klev+1)
2706      real fracc(klon,klev+1)
2707      real zf,zf2
2708      real thetath2(klon,klev),wth2(klon,klev)
2709!      common/comtherm/thetath2,wth2
2710
2711      real count_time
2712      integer isplit,nsplit,ialt
2713      parameter (nsplit=10)
2714      data isplit/0/
2715      save isplit
2716c$OMP THREADPRIVATE(isplit)
2717
2718      logical sorties
2719      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
2720      real zpspsk(klon,klev)
2721
2722c     real wmax(klon,klev),wmaxa(klon)
2723      real wmax(klon),wmaxa(klon)
2724      real wa(klon,klev,klev+1)
2725      real wd(klon,klev+1)
2726      real larg_part(klon,klev,klev+1)
2727      real fracd(klon,klev+1)
2728      real xxx(klon,klev+1)
2729      real larg_cons(klon,klev+1)
2730      real larg_detr(klon,klev+1)
2731      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
2732      real pu_therm(klon,klev),pv_therm(klon,klev)
2733      real fm(klon,klev+1),entr(klon,klev)
2734      real fmc(klon,klev+1)
2735
2736      real zcor,zdelta,zcvm5,qlbef
2737      real Tbef(klon),qsatbef(klon)
2738      real dqsat_dT,DT,num,denom
2739      REAL REPS,RLvCp,DDT0
2740      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
2741 
2742      PARAMETER (DDT0=.01)
2743
2744cCR:nouvelles variables
2745      real f_star(klon,klev+1),entr_star(klon,klev)
2746      real entr_star_tot(klon),entr_star2(klon)
2747      real f(klon), f0(klon)
2748      real zlevinter(klon)
2749      logical first
2750      data first /.false./
2751      save first
2752c$OMP THREADPRIVATE(first)
2753
2754cRC
2755
2756      character*2 str2
2757      character*10 str10
2758
2759      character (len=20) :: modname='thermcell_eau'
2760      character (len=80) :: abort_message
2761
2762      LOGICAL vtest(klon),down
2763      LOGICAL Zsat(klon)
2764
2765      EXTERNAL SCOPY
2766
2767      integer ncorrec,ll
2768      save ncorrec
2769      data ncorrec/0/
2770c$OMP THREADPRIVATE(ncorrec)
2771
2772c
2773
2774c-----------------------------------------------------------------------
2775c   initialisation:
2776c   ---------------
2777c
2778       sorties=.true.
2779      IF(ngrid.NE.klon) THEN
2780         PRINT*
2781         PRINT*,'STOP dans convadj'
2782         PRINT*,'ngrid    =',ngrid
2783         PRINT*,'klon  =',klon
2784      ENDIF
2785c
2786c Initialisation
2787      RLvCp = RLVTT/RCPD
2788      REPS  = RD/RV
2789c
2790c-----------------------------------------------------------------------
2791cAM Calcul de T,q,ql a partir de Tl et qT
2792c   ---------------------------------------------------
2793c
2794c Pr Tprec=Tl calcul de qsat
2795c Si qsat>qT T=Tl, q=qT
2796c Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
2797c On cherche DDT < DDT0
2798c
2799c defaut
2800       DO  ll=1,nlay
2801         DO ig=1,ngrid
2802            zo(ig,ll)=po(ig,ll)
2803            zl(ig,ll)=0.
2804            zh(ig,ll)=pt(ig,ll)
2805         EndDO
2806       EndDO
2807       do ig=1,ngrid
2808          Zsat(ig)=.false.
2809       enddo
2810c
2811c
2812       DO ll=1,nlay
2813c les points insatures sont definitifs
2814         DO ig=1,ngrid
2815            Tbef(ig)=pt(ig,ll)
2816            zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
2817            qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
2818            qsatbef(ig)=MIN(0.5,qsatbef(ig))
2819            zcor=1./(1.-retv*qsatbef(ig))
2820            qsatbef(ig)=qsatbef(ig)*zcor
2821            Zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig)) .gt. 0.00001)
2822         EndDO
2823
2824         DO ig=1,ngrid
2825           if (Zsat(ig)) then
2826            qlbef=max(0.,po(ig,ll)-qsatbef(ig))
2827c si sature: ql est surestime, d'ou la sous-relax
2828            DT = 0.5*RLvCp*qlbef
2829c on pourra enchainer 2 ou 3 calculs sans Do while
2830            do while (DT.gt.DDT0)
2831c il faut verifier si c,a conserve quand on repasse en insature ...
2832              Tbef(ig)=Tbef(ig)+DT
2833              zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
2834              qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
2835              qsatbef(ig)=MIN(0.5,qsatbef(ig))
2836              zcor=1./(1.-retv*qsatbef(ig))
2837              qsatbef(ig)=qsatbef(ig)*zcor
2838c on veut le signe de qlbef
2839              qlbef=po(ig,ll)-qsatbef(ig)
2840c          dqsat_dT
2841              zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
2842              zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
2843              zcor=1./(1.-retv*qsatbef(ig))
2844              dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
2845              num=-Tbef(ig)+pt(ig,ll)+RLvCp*qlbef
2846              denom=1.+RLvCp*dqsat_dT
2847              DT=num/denom
2848            enddo
2849c on ecrit de maniere conservative (sat ou non)
2850            zl(ig,ll) = max(0.,qlbef)
2851c          T = Tl +Lv/Cp ql
2852            zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
2853            zo(ig,ll) = po(ig,ll)-zl(ig,ll)
2854           endif
2855         EndDO
2856       EndDO
2857cAM fin
2858c
2859c-----------------------------------------------------------------------
2860c   incrementation eventuelle de tendances precedentes:
2861c   ---------------------------------------------------
2862
2863      print*,'0 OK convect8'
2864
2865      DO 1010 l=1,nlay
2866         DO 1015 ig=1,ngrid
2867            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
2868c            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
2869            zu(ig,l)=pu(ig,l)
2870            zv(ig,l)=pv(ig,l)
2871c            zo(ig,l)=po(ig,l)
2872c            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
2873cAM attention zh est maintenant le profil de T et plus le profil de theta !
2874c
2875c   T-> Theta
2876            ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
2877cAM Theta_v
2878            ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l))
2879     s           -zl(ig,l))
2880cAM Thetal
2881            zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
2882c           
28831015     CONTINUE
28841010  CONTINUE
2885
2886c     print*,'1 OK convect8'
2887c                       --------------------
2888c
2889c
2890c                       + + + + + + + + + + +
2891c
2892c
2893c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
2894c  wh,wt,wo ...
2895c
2896c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
2897c
2898c
2899c                       --------------------   zlev(1)
2900c                       \\\\\\\\\\\\\\\\\\\\
2901c
2902c
2903
2904c-----------------------------------------------------------------------
2905c   Calcul des altitudes des couches
2906c-----------------------------------------------------------------------
2907
2908      do l=2,nlay
2909         do ig=1,ngrid
2910            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
2911         enddo
2912      enddo
2913      do ig=1,ngrid
2914         zlev(ig,1)=0.
2915         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
2916      enddo
2917      do l=1,nlay
2918         do ig=1,ngrid
2919            zlay(ig,l)=pphi(ig,l)/RG
2920         enddo
2921      enddo
2922
2923c     print*,'2 OK convect8'
2924c-----------------------------------------------------------------------
2925c   Calcul des densites
2926c-----------------------------------------------------------------------
2927
2928      do l=1,nlay
2929         do ig=1,ngrid
2930c            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
2931             rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*ztv(ig,l))
2932         enddo
2933      enddo
2934
2935      do l=2,nlay
2936         do ig=1,ngrid
2937            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
2938         enddo
2939      enddo
2940
2941      do k=1,nlay
2942         do l=1,nlay+1
2943            do ig=1,ngrid
2944               wa(ig,k,l)=0.
2945            enddo
2946         enddo
2947      enddo
2948
2949c     print*,'3 OK convect8'
2950c------------------------------------------------------------------
2951c   Calcul de w2, quarre de w a partir de la cape
2952c   a partir de w2, on calcule wa, vitesse de l'ascendance
2953c
2954c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
2955c   w2 est stoke dans wa
2956c
2957c   ATTENTION: dans convect8, on n'utilise le calcule des wa
2958c   independants par couches que pour calculer l'entrainement
2959c   a la base et la hauteur max de l'ascendance.
2960c
2961c   Indicages:
2962c   l'ascendance provenant du niveau k traverse l'interface l avec
2963c   une vitesse wa(k,l).
2964c
2965c                       --------------------
2966c
2967c                       + + + + + + + + + +
2968c
2969c  wa(k,l)   ----       --------------------    l
2970c             /\
2971c            /||\       + + + + + + + + + +
2972c             ||
2973c             ||        --------------------
2974c             ||
2975c             ||        + + + + + + + + + +
2976c             ||
2977c             ||        --------------------
2978c             ||__
2979c             |___      + + + + + + + + + +     k
2980c
2981c                       --------------------
2982c
2983c
2984c
2985c------------------------------------------------------------------
2986
2987cCR: ponderation entrainement des couches instables
2988cdef des entr_star tels que entr=f*entr_star     
2989      do l=1,klev
2990         do ig=1,ngrid
2991            entr_star(ig,l)=0.
2992         enddo
2993      enddo
2994c determination de la longueur de la couche d entrainement
2995      do ig=1,ngrid
2996         lentr(ig)=1
2997      enddo
2998
2999con ne considere que les premieres couches instables
3000      do k=nlay-1,1,-1
3001         do ig=1,ngrid
3002            if (ztv(ig,k).gt.ztv(ig,k+1).and.
3003     s          ztv(ig,k+1).lt.ztv(ig,k+2)) then
3004               lentr(ig)=k
3005            endif
3006          enddo
3007      enddo
3008   
3009c determination du lmin: couche d ou provient le thermique
3010      do ig=1,ngrid
3011         lmin(ig)=1
3012      enddo
3013      do ig=1,ngrid
3014         do l=nlay,2,-1
3015            if (ztv(ig,l-1).gt.ztv(ig,l)) then
3016               lmin(ig)=l-1
3017            endif
3018         enddo
3019      enddo
3020c
3021c definition de l'entrainement des couches
3022      do l=1,klev-1
3023         do ig=1,ngrid
3024            if (ztv(ig,l).gt.ztv(ig,l+1).and.
3025     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
3026                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
3027     s                          (zlev(ig,l+1)-zlev(ig,l))
3028            endif
3029         enddo
3030      enddo
3031c pas de thermique si couche 1 stable
3032      do ig=1,ngrid
3033         if (lmin(ig).gt.1) then
3034            do l=1,klev
3035               entr_star(ig,l)=0.
3036            enddo
3037         endif
3038      enddo
3039c calcul de l entrainement total
3040      do ig=1,ngrid
3041         entr_star_tot(ig)=0.
3042      enddo
3043      do ig=1,ngrid
3044         do k=1,klev
3045            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
3046         enddo
3047      enddo
3048c
3049      do k=1,klev
3050         do ig=1,ngrid
3051            ztva(ig,k)=ztv(ig,k)
3052         enddo
3053      enddo
3054cRC
3055cAM:initialisations
3056      do k=1,nlay
3057         do ig=1,ngrid
3058            ztva(ig,k)=ztv(ig,k)
3059            ztla(ig,k)=zthl(ig,k)
3060            zqla(ig,k)=0.
3061            zqta(ig,k)=po(ig,k)
3062            Zsat(ig) =.false.
3063         enddo
3064      enddo
3065c
3066c     print*,'7 OK convect8'
3067      do k=1,klev+1
3068         do ig=1,ngrid
3069            zw2(ig,k)=0.
3070            fmc(ig,k)=0.
3071cCR
3072            f_star(ig,k)=0.
3073cRC
3074            larg_cons(ig,k)=0.
3075            larg_detr(ig,k)=0.
3076            wa_moy(ig,k)=0.
3077         enddo
3078      enddo
3079
3080c     print*,'8 OK convect8'
3081      do ig=1,ngrid
3082         linter(ig)=1.
3083         lmaxa(ig)=1
3084         lmix(ig)=1
3085         wmaxa(ig)=0.
3086      enddo
3087
3088cCR:
3089      do l=1,nlay-2
3090         do ig=1,ngrid
3091            if (ztv(ig,l).gt.ztv(ig,l+1)
3092     s         .and.entr_star(ig,l).gt.1.e-10
3093     s         .and.zw2(ig,l).lt.1e-10) then
3094cAM
3095               ztla(ig,l)=zthl(ig,l)
3096               zqta(ig,l)=po(ig,l)
3097               zqla(ig,l)=zl(ig,l)
3098cAM
3099               f_star(ig,l+1)=entr_star(ig,l)
3100ctest:calcul de dteta
3101               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
3102     s                     *(zlev(ig,l+1)-zlev(ig,l))
3103     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
3104               larg_detr(ig,l)=0.
3105            else if ((zw2(ig,l).ge.1e-10).and.
3106     s               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
3107               f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
3108c
3109cAM on melange Tl et qt du thermique
3110               ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+entr_star(ig,l)
3111     s                    *zthl(ig,l))/f_star(ig,l+1)
3112               zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+entr_star(ig,l)
3113     s                    *po(ig,l))/f_star(ig,l+1)
3114c
3115c               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
3116c     s                    *ztv(ig,l))/f_star(ig,l+1)
3117c
3118cAM on en deduit thetav et ql du thermique
3119               Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
3120               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
3121               qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
3122               qsatbef(ig)=MIN(0.5,qsatbef(ig))
3123               zcor=1./(1.-retv*qsatbef(ig))
3124               qsatbef(ig)=qsatbef(ig)*zcor
3125               Zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig)) .gt. 0.00001)
3126            endif
3127         enddo
3128         DO ig=1,ngrid
3129           if (Zsat(ig)) then
3130              qlbef=max(0.,zqta(ig,l)-qsatbef(ig))
3131              DT = 0.5*RLvCp*qlbef
3132              do while (DT.gt.DDT0)
3133                 Tbef(ig)=Tbef(ig)+DT
3134                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
3135                 qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
3136                 qsatbef(ig)=MIN(0.5,qsatbef(ig))
3137                 zcor=1./(1.-retv*qsatbef(ig))
3138                 qsatbef(ig)=qsatbef(ig)*zcor
3139                 qlbef=zqta(ig,l)-qsatbef(ig)
3140
3141                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
3142                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
3143                 zcor=1./(1.-retv*qsatbef(ig))
3144                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
3145                 num=-Tbef(ig)+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
3146                 denom=1.+RLvCp*dqsat_dT
3147                 DT=num/denom
3148              enddo
3149              zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
3150           endif
3151c on ecrit de maniere conservative (sat ou non)
3152c          T = Tl +Lv/Cp ql
3153           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
3154           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
3155           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
3156     s              -zqla(ig,l))-zqla(ig,l))
3157
3158        enddo
3159        DO ig=1,ngrid
3160           if (zw2(ig,l).ge.1.e-10.and.
3161     s               f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
3162c  mise a jour de la vitesse ascendante (l'air entraine de la couche
3163c  consideree commence avec une vitesse nulle).
3164c
3165               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
3166     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
3167     s                     *(zlev(ig,l+1)-zlev(ig,l))
3168            endif
3169c determination de zmax continu par interpolation lineaire
3170            if (zw2(ig,l+1).lt.0.) then
3171               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
3172     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
3173               zw2(ig,l+1)=0.
3174               lmaxa(ig)=l
3175            else
3176               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
3177            endif
3178            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
3179c   lmix est le niveau de la couche ou w (wa_moy) est maximum
3180               lmix(ig)=l+1
3181               wmaxa(ig)=wa_moy(ig,l+1)
3182            endif
3183         enddo
3184      enddo
3185c
3186c Calcul de la couche correspondant a la hauteur du thermique
3187      do ig=1,ngrid
3188         lmax(ig)=lentr(ig)
3189      enddo
3190      do ig=1,ngrid
3191         do l=nlay,lentr(ig)+1,-1
3192            if (zw2(ig,l).le.1.e-10) then
3193               lmax(ig)=l-1
3194            endif
3195         enddo
3196      enddo
3197c pas de thermique si couche 1 stable
3198      do ig=1,ngrid
3199         if (lmin(ig).gt.1) then
3200            lmax(ig)=1
3201            lmin(ig)=1
3202         endif
3203      enddo
3204c   
3205c Determination de zw2 max
3206      do ig=1,ngrid
3207         wmax(ig)=0.
3208      enddo
3209
3210      do l=1,nlay
3211         do ig=1,ngrid
3212            if (l.le.lmax(ig)) then
3213                zw2(ig,l)=sqrt(zw2(ig,l))
3214                wmax(ig)=max(wmax(ig),zw2(ig,l))
3215            else
3216                 zw2(ig,l)=0.
3217            endif
3218          enddo
3219      enddo
3220
3221c   Longueur caracteristique correspondant a la hauteur des thermiques.
3222      do  ig=1,ngrid
3223         zmax(ig)=500.
3224         zlevinter(ig)=zlev(ig,1)
3225      enddo
3226      do  ig=1,ngrid
3227c calcul de zlevinter
3228          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
3229     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
3230     s    -zlev(ig,lmax(ig)))
3231       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
3232      enddo
3233
3234c Fermeture,determination de f
3235      do ig=1,ngrid
3236         entr_star2(ig)=0.
3237      enddo
3238      do ig=1,ngrid
3239         if (entr_star_tot(ig).LT.1.e-10) then
3240            f(ig)=0.
3241         else
3242             do k=lmin(ig),lentr(ig)
3243                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
3244     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
3245             enddo
3246c Nouvelle fermeture
3247             f(ig)=wmax(ig)/(zmax(ig)*r_aspect*entr_star2(ig))
3248     s             *entr_star_tot(ig)
3249ctest
3250             if (first) then
3251             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
3252     s             *wmax(ig))
3253             endif
3254         endif
3255         f0(ig)=f(ig)
3256         first=.true.
3257      enddo
3258
3259c Calcul de l'entrainement
3260       do k=1,klev
3261         do ig=1,ngrid
3262            entr(ig,k)=f(ig)*entr_star(ig,k)
3263         enddo
3264      enddo
3265c Calcul des flux
3266      do ig=1,ngrid
3267         do l=1,lmax(ig)-1
3268            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
3269         enddo
3270      enddo
3271
3272cRC
3273
3274
3275c     print*,'9 OK convect8'
3276c     print*,'WA1 ',wa_moy
3277
3278c   determination de l'indice du debut de la mixed layer ou w decroit
3279
3280c   calcul de la largeur de chaque ascendance dans le cas conservatif.
3281c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
3282c   d'une couche est égale à la hauteur de la couche alimentante.
3283c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
3284c   de la vitesse d'entrainement horizontal dans la couche alimentante.
3285
3286      do l=2,nlay
3287         do ig=1,ngrid
3288            if (l.le.lmaxa(ig)) then
3289               zw=max(wa_moy(ig,l),1.e-10)
3290               larg_cons(ig,l)=zmax(ig)*r_aspect
3291     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
3292            endif
3293         enddo
3294      enddo
3295
3296      do l=2,nlay
3297         do ig=1,ngrid
3298            if (l.le.lmaxa(ig)) then
3299c              if (idetr.eq.0) then
3300c  cette option est finalement en dur.
3301                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3302c              else if (idetr.eq.1) then
3303c                 larg_detr(ig,l)=larg_cons(ig,l)
3304c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
3305c              else if (idetr.eq.2) then
3306c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3307c    s            *sqrt(wa_moy(ig,l))
3308c              else if (idetr.eq.4) then
3309c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3310c    s            *wa_moy(ig,l)
3311c              endif
3312            endif
3313         enddo
3314       enddo
3315
3316c     print*,'10 OK convect8'
3317c     print*,'WA2 ',wa_moy
3318c   calcul de la fraction de la maille concernée par l'ascendance en tenant
3319c   compte de l'epluchage du thermique.
3320c
3321cCR def de  zmix continu (profil parabolique des vitesses)
3322      do ig=1,ngrid
3323           if (lmix(ig).gt.1.) then
3324            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
3325     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
3326     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
3327     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
3328     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
3329     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
3330     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
3331     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3332         else
3333         zmix(ig)=0.
3334         endif
3335      enddo
3336c
3337c calcul du nouveau lmix correspondant
3338      do ig=1,ngrid
3339         do l=1,klev
3340            if (zmix(ig).ge.zlev(ig,l).and.
3341     s          zmix(ig).lt.zlev(ig,l+1)) then
3342              lmix(ig)=l
3343             endif
3344          enddo
3345      enddo
3346c
3347      do l=2,nlay
3348         do ig=1,ngrid
3349            if(larg_cons(ig,l).gt.1.) then
3350c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
3351               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
3352     s            /(r_aspect*zmax(ig))
3353c test
3354               fraca(ig,l)=max(fraca(ig,l),0.)
3355               fraca(ig,l)=min(fraca(ig,l),0.5)
3356               fracd(ig,l)=1.-fraca(ig,l)
3357               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
3358            else
3359c              wa_moy(ig,l)=0.
3360               fraca(ig,l)=0.
3361               fracc(ig,l)=0.
3362               fracd(ig,l)=1.
3363            endif
3364         enddo
3365      enddo                 
3366cCR: calcul de fracazmix
3367       do ig=1,ngrid
3368          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
3369     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
3370     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
3371     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3372       enddo
3373c
3374       do l=2,nlay
3375          do ig=1,ngrid
3376             if(larg_cons(ig,l).gt.1.) then
3377               if (l.gt.lmix(ig)) then
3378                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3379           if (idetr.eq.0) then
3380               fraca(ig,l)=fracazmix(ig)
3381           else if (idetr.eq.1) then
3382               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
3383           else if (idetr.eq.2) then
3384               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
3385           else
3386               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
3387           endif
3388c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3389               fraca(ig,l)=max(fraca(ig,l),0.)
3390               fraca(ig,l)=min(fraca(ig,l),0.5)
3391               fracd(ig,l)=1.-fraca(ig,l)
3392               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
3393             endif
3394            endif
3395         enddo
3396      enddo
3397
3398c     print*,'11 OK convect8'
3399c     print*,'Ea3 ',wa_moy
3400c------------------------------------------------------------------
3401c   Calcul de fracd, wd
3402c   somme wa - wd = 0
3403c------------------------------------------------------------------
3404
3405
3406      do ig=1,ngrid
3407         fm(ig,1)=0.
3408         fm(ig,nlay+1)=0.
3409      enddo
3410
3411      do l=2,nlay
3412           do ig=1,ngrid
3413              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
3414cCR:test
3415              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
3416     s            .and.l.gt.lmix(ig)) then
3417                 fm(ig,l)=fm(ig,l-1)
3418c                 write(1,*)'ajustement fm, l',l
3419              endif
3420c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3421cRC
3422           enddo
3423         do ig=1,ngrid
3424            if(fracd(ig,l).lt.0.1) then
3425              abort_message = 'fracd trop petit'
3426              CALL abort_gcm (modname,abort_message,1)
3427            else
3428c    vitesse descendante "diagnostique"
3429               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
3430            endif
3431         enddo
3432      enddo
3433
3434      do l=1,nlay
3435         do ig=1,ngrid
3436c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3437            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
3438         enddo
3439      enddo
3440
3441c     print*,'12 OK convect8'
3442c     print*,'WA4 ',wa_moy
3443cc------------------------------------------------------------------
3444c   calcul du transport vertical
3445c------------------------------------------------------------------
3446
3447      go to 4444
3448c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3449      do l=2,nlay-1
3450         do ig=1,ngrid
3451            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
3452     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
3453c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
3454c    s         ,fm(ig,l+1)*ptimestep
3455c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
3456             endif
3457         enddo
3458      enddo
3459
3460      do l=1,nlay
3461         do ig=1,ngrid
3462            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
3463c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
3464c    s         ,entr(ig,l)*ptimestep
3465c    s         ,'   M=',masse(ig,l)
3466             endif
3467         enddo
3468      enddo
3469
3470      do l=1,nlay
3471         do ig=1,ngrid
3472            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
3473c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
3474c    s         ,'   FM=',fm(ig,l)
3475            endif
3476            if(.not.masse(ig,l).ge.1.e-10
3477     s         .or..not.masse(ig,l).le.1.e4) then
3478c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
3479c    s         ,'   M=',masse(ig,l)
3480c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3481c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3482c     print*,'zlev(ig,l+1),zlev(ig,l)'
3483c    s                ,zlev(ig,l+1),zlev(ig,l)
3484c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3485c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3486            endif
3487            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
3488c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
3489c    s         ,'   E=',entr(ig,l)
3490            endif
3491         enddo
3492      enddo
3493
34944444   continue
3495
3496      if (w2di.eq.1) then
3497         fm0=fm0+ptimestep*(fm-fm0)/tho
3498         entr0=entr0+ptimestep*(entr-entr0)/tho
3499      else
3500         fm0=fm
3501         entr0=entr
3502      endif
3503
3504      if (1.eq.1) then
3505c         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3506c     .    ,zh,zdhadj,zha)
3507c         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3508c     .    ,zo,pdoadj,zoa)
3509         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3510     .    ,zthl,zdthladj,zta)
3511         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3512     .    ,po,pdoadj,zoa)
3513      else
3514         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
3515     .    ,zh,zdhadj,zha)
3516         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
3517     .    ,zo,pdoadj,zoa)
3518      endif
3519
3520      if (1.eq.0) then
3521         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
3522     .    ,fraca,zmax
3523     .    ,zu,zv,pduadj,pdvadj,zua,zva)
3524      else
3525         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3526     .    ,zu,pduadj,zua)
3527         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3528     .    ,zv,pdvadj,zva)
3529      endif
3530
3531      do l=1,nlay
3532         do ig=1,ngrid
3533            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
3534            zf2=zf/(1.-zf)
3535            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
3536            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
3537         enddo
3538      enddo
3539
3540
3541
3542c     print*,'13 OK convect8'
3543c     print*,'WA5 ',wa_moy
3544      do l=1,nlay
3545         do ig=1,ngrid
3546c            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
3547           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
3548         enddo
3549      enddo
3550
3551
3552c     do l=1,nlay
3553c        do ig=1,ngrid
3554c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
3555c     print*,'WARN!!! ig=',ig,'  l=',l
3556c    s         ,'   pdtadj=',pdtadj(ig,l)
3557c           endif
3558c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
3559c     print*,'WARN!!! ig=',ig,'  l=',l
3560c    s         ,'   pdoadj=',pdoadj(ig,l)
3561c           endif
3562c        enddo
3563c      enddo
3564
3565c     print*,'14 OK convect8'
3566c------------------------------------------------------------------
3567c   Calculs pour les sorties
3568c------------------------------------------------------------------
3569
3570      if(sorties) then
3571      do l=1,nlay
3572         do ig=1,ngrid
3573            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
3574            zld(ig,l)=fracd(ig,l)*zmax(ig)
3575            if(1.-fracd(ig,l).gt.1.e-10)
3576     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
3577         enddo
3578      enddo
3579
3580      do l=1,nlay
3581         do ig=1,ngrid
3582            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
3583            if (detr(ig,l).lt.0.) then
3584                entr(ig,l)=entr(ig,l)-detr(ig,l)
3585                detr(ig,l)=0.
3586c     print*,'WARNING !!! detrainement negatif ',ig,l
3587            endif
3588         enddo
3589      enddo
3590
3591c     print*,'15 OK convect8'
3592
3593      isplit=isplit+1
3594
3595
3596c #define und
3597        goto 123
3598#ifdef und
3599      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
3600      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
3601      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
3602      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
3603      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
3604      CALL writeg1d(1,nlay,zla,'la      ','la      ')
3605      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
3606      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
3607      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
3608      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
3609      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
3610      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
3611      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
3612      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
3613      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
3614      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
3615      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
3616      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
3617      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
3618      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
3619      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
3620      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
3621      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
3622      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
3623
3624      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
3625      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
3626      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
3627
3628c   recalcul des flux en diagnostique...
3629c     print*,'PAS DE TEMPS ',ptimestep
3630       call dt2F(pplev,pplay,pt,pdtadj,wh)
3631      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
3632#endif
3633123   continue
3634! #define troisD
3635#ifdef troisD
3636c       if (sorties) then
3637      print*,'Debut des wrgradsfi'
3638
3639c      print*,'16 OK convect8'
3640         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
3641         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
3642         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
3643         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
3644         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
3645         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
3646c      print*,'WA6 ',wa_moy
3647         call wrgradsfi(1,nlay,zla,'la        ','la        ')
3648         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
3649         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
3650         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
3651         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
3652         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
3653         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
3654         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
3655         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
3656         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
3657         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
3658         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
3659         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
3660         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
3661         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
3662         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
3663         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
3664         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
3665         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
3666         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
3667         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
3668         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
3669         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
3670         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
3671         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
3672         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
3673
3674         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
3675         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
3676         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
3677
3678cAM:nouveaux diagnostiques
3679         call wrgradsfi(1,nlay,zthl,'zthl        ','zthl        ')
3680         call wrgradsfi(1,nlay,zta,'zta        ','zta        ')
3681         call wrgradsfi(1,nlay,zl,'zl        ','zl        ')
3682         call wrgradsfi(1,nlay,zdthladj,'zdthladj    ',
3683     s        'zdthladj    ')
3684         call wrgradsfi(1,nlay,ztla,'ztla      ','ztla      ')
3685         call wrgradsfi(1,nlay,zqta,'zqta      ','zqta      ')
3686         call wrgradsfi(1,nlay,zqla,'zqla      ','zqla      ')
3687cCR:nouveaux diagnostiques
3688      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
3689      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
3690      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
3691      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
3692      zsortie1d(:)=lmax(:)
3693      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
3694      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
3695      zsortie1d(:)=lmix(:)
3696      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
3697      zsortie1d(:)=lentr(:)
3698      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
3699
3700c      print*,'17 OK convect8'
3701
3702         do k=1,klev/10
3703            write(str2,'(i2.2)') k
3704            str10='wa'//str2
3705            do l=1,nlay
3706               do ig=1,ngrid
3707                  zsortie(ig,l)=wa(ig,k,l)
3708               enddo
3709            enddo   
3710            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
3711            do l=1,nlay
3712               do ig=1,ngrid
3713                  zsortie(ig,l)=larg_part(ig,k,l)
3714               enddo
3715            enddo
3716            str10='la'//str2
3717            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
3718         enddo
3719
3720
3721c     print*,'18 OK convect8'
3722c      endif
3723      print*,'Fin des wrgradsfi'
3724#endif
3725
3726      endif
3727
3728c     if(wa_moy(1,4).gt.1.e-10) stop
3729
3730c     print*,'19 OK convect8'
3731      return
3732      end
3733
3734      SUBROUTINE thermcell(ngrid,nlay,ptimestep
3735     s                  ,pplay,pplev,pphi
3736     s                  ,pu,pv,pt,po
3737     s                  ,pduadj,pdvadj,pdtadj,pdoadj
3738     s                  ,fm0,entr0
3739c    s                  ,pu_therm,pv_therm
3740     s                  ,r_aspect,l_mix,w2di,tho)
3741
3742      USE dimphy
3743      IMPLICIT NONE
3744
3745c=======================================================================
3746c
3747c   Calcul du transport verticale dans la couche limite en presence
3748c   de "thermiques" explicitement representes
3749c
3750c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
3751c
3752c   le thermique est supposé homogène et dissipé par mélange avec
3753c   son environnement. la longueur l_mix contrôle l'efficacité du
3754c   mélange
3755c
3756c   Le calcul du transport des différentes espèces se fait en prenant
3757c   en compte:
3758c     1. un flux de masse montant
3759c     2. un flux de masse descendant
3760c     3. un entrainement
3761c     4. un detrainement
3762c
3763c=======================================================================
3764
3765c-----------------------------------------------------------------------
3766c   declarations:
3767c   -------------
3768
3769#include "dimensions.h"
3770cccc#include "dimphy.h"
3771#include "YOMCST.h"
3772
3773c   arguments:
3774c   ----------
3775
3776      INTEGER ngrid,nlay,w2di
3777      REAL tho
3778      real ptimestep,l_mix,r_aspect
3779      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
3780      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
3781      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
3782      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
3783      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
3784      real pphi(ngrid,nlay)
3785
3786      integer idetr
3787      save idetr
3788      data idetr/3/
3789c$OMP THREADPRIVATE(idetr)
3790
3791c   local:
3792c   ------
3793
3794      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
3795      real zsortie1d(klon)
3796c CR: on remplace lmax(klon,klev+1)
3797      INTEGER lmax(klon),lmin(klon),lentr(klon)
3798      real linter(klon)
3799      real zmix(klon), fracazmix(klon)
3800c RC
3801      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
3802
3803      real zlev(klon,klev+1),zlay(klon,klev)
3804      REAL zh(klon,klev),zdhadj(klon,klev)
3805      REAL ztv(klon,klev)
3806      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
3807      REAL wh(klon,klev+1)
3808      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
3809      real zla(klon,klev+1)
3810      real zwa(klon,klev+1)
3811      real zld(klon,klev+1)
3812      real zwd(klon,klev+1)
3813      real zsortie(klon,klev)
3814      real zva(klon,klev)
3815      real zua(klon,klev)
3816      real zoa(klon,klev)
3817
3818      real zha(klon,klev)
3819      real wa_moy(klon,klev+1)
3820      real fraca(klon,klev+1)
3821      real fracc(klon,klev+1)
3822      real zf,zf2
3823      real thetath2(klon,klev),wth2(klon,klev)
3824!      common/comtherm/thetath2,wth2
3825
3826      real count_time
3827      integer isplit,nsplit,ialt
3828      parameter (nsplit=10)
3829      data isplit/0/
3830      save isplit
3831c$OMP THREADPRIVATE(isplit)
3832
3833      logical sorties
3834      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
3835      real zpspsk(klon,klev)
3836
3837c     real wmax(klon,klev),wmaxa(klon)
3838      real wmax(klon),wmaxa(klon)
3839      real wa(klon,klev,klev+1)
3840      real wd(klon,klev+1)
3841      real larg_part(klon,klev,klev+1)
3842      real fracd(klon,klev+1)
3843      real xxx(klon,klev+1)
3844      real larg_cons(klon,klev+1)
3845      real larg_detr(klon,klev+1)
3846      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
3847      real pu_therm(klon,klev),pv_therm(klon,klev)
3848      real fm(klon,klev+1),entr(klon,klev)
3849      real fmc(klon,klev+1)
3850
3851cCR:nouvelles variables
3852      real f_star(klon,klev+1),entr_star(klon,klev)
3853      real entr_star_tot(klon),entr_star2(klon)
3854      real f(klon), f0(klon)
3855      real zlevinter(klon)
3856      logical first
3857      data first /.false./
3858      save first
3859c$OMP THREADPRIVATE(first)
3860cRC
3861
3862      character*2 str2
3863      character*10 str10
3864
3865      character (len=20) :: modname='thermcell'
3866      character (len=80) :: abort_message
3867
3868      LOGICAL vtest(klon),down
3869
3870      EXTERNAL SCOPY
3871
3872      integer ncorrec,ll
3873      save ncorrec
3874      data ncorrec/0/
3875c$OMP THREADPRIVATE(ncorrec)
3876     
3877c
3878c-----------------------------------------------------------------------
3879c   initialisation:
3880c   ---------------
3881c
3882       sorties=.true.
3883      IF(ngrid.NE.klon) THEN
3884         PRINT*
3885         PRINT*,'STOP dans convadj'
3886         PRINT*,'ngrid    =',ngrid
3887         PRINT*,'klon  =',klon
3888      ENDIF
3889c
3890c-----------------------------------------------------------------------
3891c   incrementation eventuelle de tendances precedentes:
3892c   ---------------------------------------------------
3893
3894       print*,'0 OK convect8'
3895
3896      DO 1010 l=1,nlay
3897         DO 1015 ig=1,ngrid
3898            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
3899            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
3900            zu(ig,l)=pu(ig,l)
3901            zv(ig,l)=pv(ig,l)
3902            zo(ig,l)=po(ig,l)
3903            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
39041015     CONTINUE
39051010  CONTINUE
3906
3907       print*,'1 OK convect8'
3908c                       --------------------
3909c
3910c
3911c                       + + + + + + + + + + +
3912c
3913c
3914c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
3915c  wh,wt,wo ...
3916c
3917c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
3918c
3919c
3920c                       --------------------   zlev(1)
3921c                       \\\\\\\\\\\\\\\\\\\\
3922c
3923c
3924
3925c-----------------------------------------------------------------------
3926c   Calcul des altitudes des couches
3927c-----------------------------------------------------------------------
3928
3929      do l=2,nlay
3930         do ig=1,ngrid
3931            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
3932         enddo
3933      enddo
3934      do ig=1,ngrid
3935         zlev(ig,1)=0.
3936         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
3937      enddo
3938      do l=1,nlay
3939         do ig=1,ngrid
3940            zlay(ig,l)=pphi(ig,l)/RG
3941         enddo
3942      enddo
3943
3944c      print*,'2 OK convect8'
3945c-----------------------------------------------------------------------
3946c   Calcul des densites
3947c-----------------------------------------------------------------------
3948
3949      do l=1,nlay
3950         do ig=1,ngrid
3951            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
3952         enddo
3953      enddo
3954
3955      do l=2,nlay
3956         do ig=1,ngrid
3957            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
3958         enddo
3959      enddo
3960
3961      do k=1,nlay
3962         do l=1,nlay+1
3963            do ig=1,ngrid
3964               wa(ig,k,l)=0.
3965            enddo
3966         enddo
3967      enddo
3968
3969c      print*,'3 OK convect8'
3970c------------------------------------------------------------------
3971c   Calcul de w2, quarre de w a partir de la cape
3972c   a partir de w2, on calcule wa, vitesse de l'ascendance
3973c
3974c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
3975c   w2 est stoke dans wa
3976c
3977c   ATTENTION: dans convect8, on n'utilise le calcule des wa
3978c   independants par couches que pour calculer l'entrainement
3979c   a la base et la hauteur max de l'ascendance.
3980c
3981c   Indicages:
3982c   l'ascendance provenant du niveau k traverse l'interface l avec
3983c   une vitesse wa(k,l).
3984c
3985c                       --------------------
3986c
3987c                       + + + + + + + + + +
3988c
3989c  wa(k,l)   ----       --------------------    l
3990c             /\
3991c            /||\       + + + + + + + + + +
3992c             ||
3993c             ||        --------------------
3994c             ||
3995c             ||        + + + + + + + + + +
3996c             ||
3997c             ||        --------------------
3998c             ||__
3999c             |___      + + + + + + + + + +     k
4000c
4001c                       --------------------
4002c
4003c
4004c
4005c------------------------------------------------------------------
4006
4007cCR: ponderation entrainement des couches instables
4008cdef des entr_star tels que entr=f*entr_star     
4009      do l=1,klev
4010         do ig=1,ngrid
4011            entr_star(ig,l)=0.
4012         enddo
4013      enddo
4014c determination de la longueur de la couche d entrainement
4015      do ig=1,ngrid
4016         lentr(ig)=1
4017      enddo
4018
4019con ne considere que les premieres couches instables
4020      do k=nlay-2,1,-1
4021         do ig=1,ngrid
4022            if (ztv(ig,k).gt.ztv(ig,k+1).and.
4023     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
4024               lentr(ig)=k
4025            endif
4026          enddo
4027      enddo
4028   
4029c determination du lmin: couche d ou provient le thermique
4030      do ig=1,ngrid
4031         lmin(ig)=1
4032      enddo
4033      do ig=1,ngrid
4034         do l=nlay,2,-1
4035            if (ztv(ig,l-1).gt.ztv(ig,l)) then
4036               lmin(ig)=l-1
4037            endif
4038         enddo
4039      enddo
4040c
4041c definition de l'entrainement des couches
4042      do l=1,klev-1
4043         do ig=1,ngrid
4044            if (ztv(ig,l).gt.ztv(ig,l+1).and.
4045     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
4046                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
4047     s                           (zlev(ig,l+1)-zlev(ig,l))
4048            endif
4049         enddo
4050      enddo
4051c pas de thermique si couches 1->5 stables
4052      do ig=1,ngrid
4053         if (lmin(ig).gt.5) then
4054            do l=1,klev
4055               entr_star(ig,l)=0.
4056            enddo
4057         endif
4058      enddo
4059c calcul de l entrainement total
4060      do ig=1,ngrid
4061         entr_star_tot(ig)=0.
4062      enddo
4063      do ig=1,ngrid
4064         do k=1,klev
4065            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
4066         enddo
4067      enddo
4068c
4069      print*,'fin calcul entr_star'
4070      do k=1,klev
4071         do ig=1,ngrid
4072            ztva(ig,k)=ztv(ig,k)
4073         enddo
4074      enddo
4075cRC
4076c      print*,'7 OK convect8'
4077      do k=1,klev+1
4078         do ig=1,ngrid
4079            zw2(ig,k)=0.
4080            fmc(ig,k)=0.
4081cCR
4082            f_star(ig,k)=0.
4083cRC
4084            larg_cons(ig,k)=0.
4085            larg_detr(ig,k)=0.
4086            wa_moy(ig,k)=0.
4087         enddo
4088      enddo
4089
4090c      print*,'8 OK convect8'
4091      do ig=1,ngrid
4092         linter(ig)=1.
4093         lmaxa(ig)=1
4094         lmix(ig)=1
4095         wmaxa(ig)=0.
4096      enddo
4097
4098cCR:
4099      do l=1,nlay-2
4100         do ig=1,ngrid
4101            if (ztv(ig,l).gt.ztv(ig,l+1)
4102     s         .and.entr_star(ig,l).gt.1.e-10
4103     s         .and.zw2(ig,l).lt.1e-10) then
4104               f_star(ig,l+1)=entr_star(ig,l)
4105ctest:calcul de dteta
4106               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
4107     s                     *(zlev(ig,l+1)-zlev(ig,l))
4108     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
4109               larg_detr(ig,l)=0.
4110            else if ((zw2(ig,l).ge.1e-10).and.
4111     s               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
4112               f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
4113               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
4114     s                    *ztv(ig,l))/f_star(ig,l+1)
4115               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
4116     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
4117     s                     *(zlev(ig,l+1)-zlev(ig,l))
4118            endif
4119c determination de zmax continu par interpolation lineaire
4120            if (zw2(ig,l+1).lt.0.) then
4121ctest
4122               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
4123                  print*,'pb linter'
4124               endif
4125               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
4126     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
4127               zw2(ig,l+1)=0.
4128               lmaxa(ig)=l
4129            else
4130               if (zw2(ig,l+1).lt.0.) then
4131                  print*,'pb1 zw2<0'
4132               endif
4133               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
4134            endif
4135            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
4136c   lmix est le niveau de la couche ou w (wa_moy) est maximum
4137               lmix(ig)=l+1
4138               wmaxa(ig)=wa_moy(ig,l+1)
4139            endif
4140         enddo
4141      enddo
4142      print*,'fin calcul zw2'
4143c
4144c Calcul de la couche correspondant a la hauteur du thermique
4145      do ig=1,ngrid
4146         lmax(ig)=lentr(ig)
4147      enddo
4148      do ig=1,ngrid
4149         do l=nlay,lentr(ig)+1,-1
4150            if (zw2(ig,l).le.1.e-10) then
4151               lmax(ig)=l-1
4152            endif
4153         enddo
4154      enddo
4155c pas de thermique si couches 1->5 stables
4156      do ig=1,ngrid
4157         if (lmin(ig).gt.5) then
4158            lmax(ig)=1
4159            lmin(ig)=1
4160         endif
4161      enddo
4162c   
4163c Determination de zw2 max
4164      do ig=1,ngrid
4165         wmax(ig)=0.
4166      enddo
4167
4168      do l=1,nlay
4169         do ig=1,ngrid
4170            if (l.le.lmax(ig)) then
4171                if (zw2(ig,l).lt.0.)then
4172                  print*,'pb2 zw2<0'
4173                endif
4174                zw2(ig,l)=sqrt(zw2(ig,l))
4175                wmax(ig)=max(wmax(ig),zw2(ig,l))
4176            else
4177                 zw2(ig,l)=0.
4178            endif
4179          enddo
4180      enddo
4181
4182c   Longueur caracteristique correspondant a la hauteur des thermiques.
4183      do  ig=1,ngrid
4184         zmax(ig)=0.
4185         zlevinter(ig)=zlev(ig,1)
4186      enddo
4187      do  ig=1,ngrid
4188c calcul de zlevinter
4189          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
4190     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
4191     s    -zlev(ig,lmax(ig)))
4192       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
4193      enddo
4194
4195      print*,'avant fermeture'
4196c Fermeture,determination de f
4197      do ig=1,ngrid
4198         entr_star2(ig)=0.
4199      enddo
4200      do ig=1,ngrid
4201         if (entr_star_tot(ig).LT.1.e-10) then
4202            f(ig)=0.
4203         else
4204             do k=lmin(ig),lentr(ig)
4205                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
4206     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
4207             enddo
4208c Nouvelle fermeture
4209             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
4210     s             *entr_star2(ig))*entr_star_tot(ig)
4211ctest
4212c             if (first) then
4213c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
4214c     s             *wmax(ig))
4215c             endif
4216         endif
4217c         f0(ig)=f(ig)
4218c         first=.true.
4219      enddo
4220      print*,'apres fermeture'
4221
4222c Calcul de l'entrainement
4223       do k=1,klev
4224         do ig=1,ngrid
4225            entr(ig,k)=f(ig)*entr_star(ig,k)
4226         enddo
4227      enddo
4228c Calcul des flux
4229      do ig=1,ngrid
4230         do l=1,lmax(ig)-1
4231            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
4232         enddo
4233      enddo
4234
4235cRC
4236
4237
4238c      print*,'9 OK convect8'
4239c     print*,'WA1 ',wa_moy
4240
4241c   determination de l'indice du debut de la mixed layer ou w decroit
4242
4243c   calcul de la largeur de chaque ascendance dans le cas conservatif.
4244c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
4245c   d'une couche est égale à la hauteur de la couche alimentante.
4246c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
4247c   de la vitesse d'entrainement horizontal dans la couche alimentante.
4248
4249      do l=2,nlay
4250         do ig=1,ngrid
4251            if (l.le.lmaxa(ig)) then
4252               zw=max(wa_moy(ig,l),1.e-10)
4253               larg_cons(ig,l)=zmax(ig)*r_aspect
4254     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
4255            endif
4256         enddo
4257      enddo
4258
4259      do l=2,nlay
4260         do ig=1,ngrid
4261            if (l.le.lmaxa(ig)) then
4262c              if (idetr.eq.0) then
4263c  cette option est finalement en dur.
4264                  if ((l_mix*zlev(ig,l)).lt.0.)then
4265                   print*,'pb l_mix*zlev<0'
4266                  endif
4267                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4268c              else if (idetr.eq.1) then
4269c                 larg_detr(ig,l)=larg_cons(ig,l)
4270c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
4271c              else if (idetr.eq.2) then
4272c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4273c    s            *sqrt(wa_moy(ig,l))
4274c              else if (idetr.eq.4) then
4275c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4276c    s            *wa_moy(ig,l)
4277c              endif
4278            endif
4279         enddo
4280       enddo
4281
4282c      print*,'10 OK convect8'
4283c     print*,'WA2 ',wa_moy
4284c   calcul de la fraction de la maille concernée par l'ascendance en tenant
4285c   compte de l'epluchage du thermique.
4286c
4287cCR def de  zmix continu (profil parabolique des vitesses)
4288      do ig=1,ngrid
4289           if (lmix(ig).gt.1.) then
4290c test
4291              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
4292     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
4293     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
4294     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
4295     s        then
4296c             
4297            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
4298     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
4299     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
4300     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
4301     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
4302     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
4303     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
4304     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
4305            else
4306            zmix(ig)=zlev(ig,lmix(ig))
4307            print*,'pb zmix'
4308            endif
4309         else
4310         zmix(ig)=0.
4311         endif
4312ctest
4313         if ((zmax(ig)-zmix(ig)).lt.0.) then
4314            zmix(ig)=0.99*zmax(ig)
4315c            print*,'pb zmix>zmax'
4316         endif
4317      enddo
4318c
4319c calcul du nouveau lmix correspondant
4320      do ig=1,ngrid
4321         do l=1,klev
4322            if (zmix(ig).ge.zlev(ig,l).and.
4323     s          zmix(ig).lt.zlev(ig,l+1)) then
4324              lmix(ig)=l
4325             endif
4326          enddo
4327      enddo
4328c
4329      do l=2,nlay
4330         do ig=1,ngrid
4331            if(larg_cons(ig,l).gt.1.) then
4332c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
4333               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
4334     s            /(r_aspect*zmax(ig))
4335c test
4336               fraca(ig,l)=max(fraca(ig,l),0.)
4337               fraca(ig,l)=min(fraca(ig,l),0.5)
4338               fracd(ig,l)=1.-fraca(ig,l)
4339               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
4340            else
4341c              wa_moy(ig,l)=0.
4342               fraca(ig,l)=0.
4343               fracc(ig,l)=0.
4344               fracd(ig,l)=1.
4345            endif
4346         enddo
4347      enddo                 
4348cCR: calcul de fracazmix
4349       do ig=1,ngrid
4350          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
4351     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
4352     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
4353     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
4354       enddo
4355c
4356       do l=2,nlay
4357          do ig=1,ngrid
4358             if(larg_cons(ig,l).gt.1.) then
4359               if (l.gt.lmix(ig)) then
4360ctest
4361                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
4362c                   print*,'pb xxx'
4363                   xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
4364                 else
4365                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
4366                 endif
4367           if (idetr.eq.0) then
4368               fraca(ig,l)=fracazmix(ig)
4369           else if (idetr.eq.1) then
4370               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
4371           else if (idetr.eq.2) then
4372               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
4373           else
4374               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
4375           endif
4376c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
4377               fraca(ig,l)=max(fraca(ig,l),0.)
4378               fraca(ig,l)=min(fraca(ig,l),0.5)
4379               fracd(ig,l)=1.-fraca(ig,l)
4380               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
4381             endif
4382            endif
4383         enddo
4384      enddo
4385     
4386      print*,'fin calcul fraca'
4387c      print*,'11 OK convect8'
4388c     print*,'Ea3 ',wa_moy
4389c------------------------------------------------------------------
4390c   Calcul de fracd, wd
4391c   somme wa - wd = 0
4392c------------------------------------------------------------------
4393
4394
4395      do ig=1,ngrid
4396         fm(ig,1)=0.
4397         fm(ig,nlay+1)=0.
4398      enddo
4399
4400      do l=2,nlay
4401           do ig=1,ngrid
4402              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
4403cCR:test
4404              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
4405     s            .and.l.gt.lmix(ig)) then
4406                 fm(ig,l)=fm(ig,l-1)
4407c                 write(1,*)'ajustement fm, l',l
4408              endif
4409c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
4410cRC
4411           enddo
4412         do ig=1,ngrid
4413            if(fracd(ig,l).lt.0.1) then
4414              abort_message = 'fracd trop petit'
4415              CALL abort_gcm (modname,abort_message,1)
4416            else
4417c    vitesse descendante "diagnostique"
4418               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
4419            endif
4420         enddo
4421      enddo
4422
4423      do l=1,nlay
4424         do ig=1,ngrid
4425c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
4426            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
4427         enddo
4428      enddo
4429
4430      print*,'12 OK convect8'
4431c     print*,'WA4 ',wa_moy
4432cc------------------------------------------------------------------
4433c   calcul du transport vertical
4434c------------------------------------------------------------------
4435
4436      go to 4444
4437c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
4438      do l=2,nlay-1
4439         do ig=1,ngrid
4440            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
4441     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
4442c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
4443c    s         ,fm(ig,l+1)*ptimestep
4444c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
4445             endif
4446         enddo
4447      enddo
4448
4449      do l=1,nlay
4450         do ig=1,ngrid
4451            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
4452c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
4453c    s         ,entr(ig,l)*ptimestep
4454c    s         ,'   M=',masse(ig,l)
4455             endif
4456         enddo
4457      enddo
4458
4459      do l=1,nlay
4460         do ig=1,ngrid
4461            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
4462c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
4463c    s         ,'   FM=',fm(ig,l)
4464            endif
4465            if(.not.masse(ig,l).ge.1.e-10
4466     s         .or..not.masse(ig,l).le.1.e4) then
4467c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
4468c    s         ,'   M=',masse(ig,l)
4469c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
4470c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
4471c     print*,'zlev(ig,l+1),zlev(ig,l)'
4472c    s                ,zlev(ig,l+1),zlev(ig,l)
4473c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
4474c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
4475            endif
4476            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
4477c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
4478c    s         ,'   E=',entr(ig,l)
4479            endif
4480         enddo
4481      enddo
4482
44834444   continue
4484
4485cCR:redefinition du entr
4486       do l=1,nlay
4487         do ig=1,ngrid
4488            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4489            if (detr(ig,l).lt.0.) then
4490                entr(ig,l)=entr(ig,l)-detr(ig,l)
4491                detr(ig,l)=0.
4492c     print*,'WARNING !!! detrainement negatif ',ig,l
4493            endif
4494         enddo
4495      enddo
4496cRC
4497      if (w2di.eq.1) then
4498         fm0=fm0+ptimestep*(fm-fm0)/tho
4499         entr0=entr0+ptimestep*(entr-entr0)/tho
4500      else
4501         fm0=fm
4502         entr0=entr
4503      endif
4504
4505      if (1.eq.1) then
4506         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4507     .    ,zh,zdhadj,zha)
4508         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4509     .    ,zo,pdoadj,zoa)
4510      else
4511         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
4512     .    ,zh,zdhadj,zha)
4513         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
4514     .    ,zo,pdoadj,zoa)
4515      endif
4516
4517      if (1.eq.0) then
4518         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
4519     .    ,fraca,zmax
4520     .    ,zu,zv,pduadj,pdvadj,zua,zva)
4521      else
4522         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4523     .    ,zu,pduadj,zua)
4524         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4525     .    ,zv,pdvadj,zva)
4526      endif
4527
4528      do l=1,nlay
4529         do ig=1,ngrid
4530            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
4531            zf2=zf/(1.-zf)
4532            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
4533            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
4534         enddo
4535      enddo
4536
4537
4538
4539c     print*,'13 OK convect8'
4540c     print*,'WA5 ',wa_moy
4541      do l=1,nlay
4542         do ig=1,ngrid
4543            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
4544         enddo
4545      enddo
4546
4547
4548c     do l=1,nlay
4549c        do ig=1,ngrid
4550c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
4551c     print*,'WARN!!! ig=',ig,'  l=',l
4552c    s         ,'   pdtadj=',pdtadj(ig,l)
4553c           endif
4554c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
4555c     print*,'WARN!!! ig=',ig,'  l=',l
4556c    s         ,'   pdoadj=',pdoadj(ig,l)
4557c           endif
4558c        enddo
4559c      enddo
4560
4561      print*,'14 OK convect8'
4562c------------------------------------------------------------------
4563c   Calculs pour les sorties
4564c------------------------------------------------------------------
4565
4566      if(sorties) then
4567      do l=1,nlay
4568         do ig=1,ngrid
4569            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
4570            zld(ig,l)=fracd(ig,l)*zmax(ig)
4571            if(1.-fracd(ig,l).gt.1.e-10)
4572     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
4573         enddo
4574      enddo
4575
4576cdeja fait
4577c      do l=1,nlay
4578c         do ig=1,ngrid
4579c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4580c            if (detr(ig,l).lt.0.) then
4581c                entr(ig,l)=entr(ig,l)-detr(ig,l)
4582c                detr(ig,l)=0.
4583c     print*,'WARNING !!! detrainement negatif ',ig,l
4584c            endif
4585c         enddo
4586c      enddo
4587
4588c     print*,'15 OK convect8'
4589
4590      isplit=isplit+1
4591
4592
4593c #define und
4594        goto 123
4595#ifdef und
4596      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
4597      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
4598      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
4599      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
4600      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
4601      CALL writeg1d(1,nlay,zla,'la      ','la      ')
4602      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
4603      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
4604      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
4605      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
4606      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
4607      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
4608      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
4609      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
4610      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
4611      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
4612      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
4613      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
4614      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
4615      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
4616      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
4617      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
4618      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
4619      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
4620
4621      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
4622      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
4623      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
4624
4625c   recalcul des flux en diagnostique...
4626c     print*,'PAS DE TEMPS ',ptimestep
4627       call dt2F(pplev,pplay,pt,pdtadj,wh)
4628      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
4629#endif
4630123   continue
4631#define troisD
4632#ifdef troisD
4633c       if (sorties) then
4634      print*,'Debut des wrgradsfi'
4635
4636c      print*,'16 OK convect8'
4637         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
4638         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
4639         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
4640         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
4641         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
4642         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
4643c      print*,'WA6 ',wa_moy
4644         call wrgradsfi(1,nlay,zla,'la        ','la        ')
4645         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
4646         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
4647         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
4648         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
4649         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
4650         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
4651         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
4652         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
4653         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
4654         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
4655         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
4656         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
4657         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
4658         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
4659         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
4660         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
4661         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
4662         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
4663         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
4664         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
4665         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
4666         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
4667         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
4668         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
4669         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
4670
4671         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
4672         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
4673         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
4674
4675cCR:nouveaux diagnostiques
4676      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
4677      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
4678      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
4679      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
4680      zsortie1d(:)=lmax(:)
4681      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
4682      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
4683      zsortie1d(:)=lmix(:)
4684      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
4685      zsortie1d(:)=lentr(:)
4686      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
4687
4688c      print*,'17 OK convect8'
4689
4690         do k=1,klev/10
4691            write(str2,'(i2.2)') k
4692            str10='wa'//str2
4693            do l=1,nlay
4694               do ig=1,ngrid
4695                  zsortie(ig,l)=wa(ig,k,l)
4696               enddo
4697            enddo   
4698            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
4699            do l=1,nlay
4700               do ig=1,ngrid
4701                  zsortie(ig,l)=larg_part(ig,k,l)
4702               enddo
4703            enddo
4704            str10='la'//str2
4705            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
4706         enddo
4707
4708
4709c     print*,'18 OK convect8'
4710c      endif
4711      print*,'Fin des wrgradsfi'
4712#endif
4713
4714      endif
4715
4716c     if(wa_moy(1,4).gt.1.e-10) stop
4717
4718      print*,'19 OK convect8'
4719      return
4720      end
4721
4722      subroutine dqthermcell(ngrid,nlay,ptimestep,fm,entr,
4723     .           masse,q,dq,qa)
4724      USE dimphy
4725      implicit none
4726
4727c=======================================================================
4728c
4729c   Calcul du transport verticale dans la couche limite en presence
4730c   de "thermiques" explicitement representes
4731c   calcul du dq/dt une fois qu'on connait les ascendances
4732c
4733c=======================================================================
4734
4735#include "dimensions.h"
4736cccc#include "dimphy.h"
4737
4738      integer ngrid,nlay
4739
4740      real ptimestep
4741      real masse(ngrid,nlay),fm(ngrid,nlay+1)
4742      real entr(ngrid,nlay)
4743      real q(ngrid,nlay)
4744      real dq(ngrid,nlay)
4745
4746      real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
4747
4748      integer ig,k
4749
4750c   calcul du detrainement
4751
4752      do k=1,nlay
4753         do ig=1,ngrid
4754            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
4755ctest
4756            if (detr(ig,k).lt.0.) then
4757               entr(ig,k)=entr(ig,k)-detr(ig,k)
4758               detr(ig,k)=0.
4759c               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
4760c     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
4761            endif
4762            if (fm(ig,k+1).lt.0.) then
4763c               print*,'fm2<0!!!'
4764            endif
4765            if (entr(ig,k).lt.0.) then
4766c               print*,'entr2<0!!!'
4767            endif
4768         enddo
4769      enddo
4770
4771c   calcul de la valeur dans les ascendances
4772      do ig=1,ngrid
4773         qa(ig,1)=q(ig,1)
4774      enddo
4775
4776      do k=2,nlay
4777         do ig=1,ngrid
4778            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
4779     s         1.e-5*masse(ig,k)) then
4780         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))
4781     s         /(fm(ig,k+1)+detr(ig,k))
4782            else
4783               qa(ig,k)=q(ig,k)
4784            endif
4785            if (qa(ig,k).lt.0.) then
4786c               print*,'qa<0!!!'
4787            endif
4788            if (q(ig,k).lt.0.) then
4789c               print*,'q<0!!!'
4790            endif
4791         enddo
4792      enddo
4793
4794      do k=2,nlay
4795         do ig=1,ngrid
4796c             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
4797            wqd(ig,k)=fm(ig,k)*q(ig,k)
4798            if (wqd(ig,k).lt.0.) then
4799c               print*,'wqd<0!!!'
4800            endif
4801         enddo
4802      enddo
4803      do ig=1,ngrid
4804         wqd(ig,1)=0.
4805         wqd(ig,nlay+1)=0.
4806      enddo
4807     
4808      do k=1,nlay
4809         do ig=1,ngrid
4810            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)
4811     s               -wqd(ig,k)+wqd(ig,k+1))
4812     s               /masse(ig,k)
4813c            if (dq(ig,k).lt.0.) then
4814c               print*,'dq<0!!!'
4815c            endif
4816         enddo
4817      enddo
4818
4819      return
4820      end
4821      subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
4822     .    ,fraca,larga
4823     .    ,u,v,du,dv,ua,va)
4824      USE dimphy
4825      implicit none
4826
4827c=======================================================================
4828c
4829c   Calcul du transport verticale dans la couche limite en presence
4830c   de "thermiques" explicitement representes
4831c   calcul du dq/dt une fois qu'on connait les ascendances
4832c
4833c=======================================================================
4834
4835#include "dimensions.h"
4836cccc#include "dimphy.h"
4837
4838      integer ngrid,nlay
4839
4840      real ptimestep
4841      real masse(ngrid,nlay),fm(ngrid,nlay+1)
4842      real fraca(ngrid,nlay+1)
4843      real larga(ngrid)
4844      real entr(ngrid,nlay)
4845      real u(ngrid,nlay)
4846      real ua(ngrid,nlay)
4847      real du(ngrid,nlay)
4848      real v(ngrid,nlay)
4849      real va(ngrid,nlay)
4850      real dv(ngrid,nlay)
4851
4852      real qa(klon,klev),detr(klon,klev)
4853      real wvd(klon,klev+1),wud(klon,klev+1)
4854      real gamma0,gamma(klon,klev+1)
4855      real dua,dva
4856      integer iter
4857
4858      integer ig,k
4859
4860c   calcul du detrainement
4861
4862      do k=1,nlay
4863         do ig=1,ngrid
4864            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
4865         enddo
4866      enddo
4867
4868c   calcul de la valeur dans les ascendances
4869      do ig=1,ngrid
4870         ua(ig,1)=u(ig,1)
4871         va(ig,1)=v(ig,1)
4872      enddo
4873
4874      do k=2,nlay
4875         do ig=1,ngrid
4876            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
4877     s         1.e-5*masse(ig,k)) then
4878c   On itère sur la valeur du coeff de freinage.
4879c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4880               gamma0=masse(ig,k)
4881     s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
4882     s         *0.5/larga(ig)
4883c              gamma0=0.
4884c   la première fois on multiplie le coefficient de freinage
4885c   par le module du vent dans la couche en dessous.
4886               dua=ua(ig,k-1)-u(ig,k-1)
4887               dva=va(ig,k-1)-v(ig,k-1)
4888               do iter=1,5
4889                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
4890                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
4891     s               +(entr(ig,k)+gamma(ig,k))*u(ig,k))
4892     s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4893                  va(ig,k)=(fm(ig,k)*va(ig,k-1)
4894     s               +(entr(ig,k)+gamma(ig,k))*v(ig,k))
4895     s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4896c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4897                  dua=ua(ig,k)-u(ig,k)
4898                  dva=va(ig,k)-v(ig,k)
4899               enddo
4900            else
4901               ua(ig,k)=u(ig,k)
4902               va(ig,k)=v(ig,k)
4903               gamma(ig,k)=0.
4904            endif
4905         enddo
4906      enddo
4907
4908      do k=2,nlay
4909         do ig=1,ngrid
4910            wud(ig,k)=fm(ig,k)*u(ig,k)
4911            wvd(ig,k)=fm(ig,k)*v(ig,k)
4912         enddo
4913      enddo
4914      do ig=1,ngrid
4915         wud(ig,1)=0.
4916         wud(ig,nlay+1)=0.
4917         wvd(ig,1)=0.
4918         wvd(ig,nlay+1)=0.
4919      enddo
4920
4921      do k=1,nlay
4922         do ig=1,ngrid
4923            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
4924     s               -(entr(ig,k)+gamma(ig,k))*u(ig,k)
4925     s               -wud(ig,k)+wud(ig,k+1))
4926     s               /masse(ig,k)
4927            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
4928     s               -(entr(ig,k)+gamma(ig,k))*v(ig,k)
4929     s               -wvd(ig,k)+wvd(ig,k+1))
4930     s               /masse(ig,k)
4931         enddo
4932      enddo
4933
4934      return
4935      end
4936      subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
4937     .    ,q,dq,qa)
4938      USE dimphy
4939      implicit none
4940
4941c=======================================================================
4942c
4943c   Calcul du transport verticale dans la couche limite en presence
4944c   de "thermiques" explicitement representes
4945c   calcul du dq/dt une fois qu'on connait les ascendances
4946c
4947c=======================================================================
4948
4949#include "dimensions.h"
4950cccc#include "dimphy.h"
4951
4952      integer ngrid,nlay
4953
4954      real ptimestep
4955      real masse(ngrid,nlay),fm(ngrid,nlay+1)
4956      real entr(ngrid,nlay),frac(ngrid,nlay)
4957      real q(ngrid,nlay)
4958      real dq(ngrid,nlay)
4959
4960      real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
4961      real qe(klon,klev),zf,zf2
4962
4963      integer ig,k
4964
4965c   calcul du detrainement
4966
4967      do k=1,nlay
4968         do ig=1,ngrid
4969            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
4970         enddo
4971      enddo
4972
4973c   calcul de la valeur dans les ascendances
4974      do ig=1,ngrid
4975         qa(ig,1)=q(ig,1)
4976         qe(ig,1)=q(ig,1)
4977      enddo
4978
4979      do k=2,nlay
4980         do ig=1,ngrid
4981            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
4982     s         1.e-5*masse(ig,k)) then
4983               zf=0.5*(frac(ig,k)+frac(ig,k+1))
4984               zf2=1./(1.-zf)
4985               qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+zf2*entr(ig,k)*q(ig,k))
4986     s         /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2)
4987               qe(ig,k)=(q(ig,k)-zf*qa(ig,k))*zf2
4988            else
4989               qa(ig,k)=q(ig,k)
4990               qe(ig,k)=q(ig,k)
4991            endif
4992         enddo
4993      enddo
4994
4995      do k=2,nlay
4996         do ig=1,ngrid
4997c             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
4998            wqd(ig,k)=fm(ig,k)*qe(ig,k)
4999         enddo
5000      enddo
5001      do ig=1,ngrid
5002         wqd(ig,1)=0.
5003         wqd(ig,nlay+1)=0.
5004      enddo
5005
5006      do k=1,nlay
5007         do ig=1,ngrid
5008            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*qe(ig,k)
5009     s               -wqd(ig,k)+wqd(ig,k+1))
5010     s               /masse(ig,k)
5011         enddo
5012      enddo
5013
5014      return
5015      end
5016      subroutine dvthermcell2(ngrid,nlay,ptimestep,fm,entr,masse
5017     .    ,fraca,larga
5018     .    ,u,v,du,dv,ua,va)
5019      use dimphy
5020      implicit none
5021
5022c=======================================================================
5023c
5024c   Calcul du transport verticale dans la couche limite en presence
5025c   de "thermiques" explicitement representes
5026c   calcul du dq/dt une fois qu'on connait les ascendances
5027c
5028c=======================================================================
5029
5030#include "dimensions.h"
5031cccc#include "dimphy.h"
5032
5033      integer ngrid,nlay
5034
5035      real ptimestep
5036      real masse(ngrid,nlay),fm(ngrid,nlay+1)
5037      real fraca(ngrid,nlay+1)
5038      real larga(ngrid)
5039      real entr(ngrid,nlay)
5040      real u(ngrid,nlay)
5041      real ua(ngrid,nlay)
5042      real du(ngrid,nlay)
5043      real v(ngrid,nlay)
5044      real va(ngrid,nlay)
5045      real dv(ngrid,nlay)
5046
5047      real qa(klon,klev),detr(klon,klev),zf,zf2
5048      real wvd(klon,klev+1),wud(klon,klev+1)
5049      real gamma0,gamma(klon,klev+1)
5050      real ue(klon,klev),ve(klon,klev)
5051      real dua,dva
5052      integer iter
5053
5054      integer ig,k
5055
5056c   calcul du detrainement
5057
5058      do k=1,nlay
5059         do ig=1,ngrid
5060            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
5061         enddo
5062      enddo
5063
5064c   calcul de la valeur dans les ascendances
5065      do ig=1,ngrid
5066         ua(ig,1)=u(ig,1)
5067         va(ig,1)=v(ig,1)
5068         ue(ig,1)=u(ig,1)
5069         ve(ig,1)=v(ig,1)
5070      enddo
5071
5072      do k=2,nlay
5073         do ig=1,ngrid
5074            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
5075     s         1.e-5*masse(ig,k)) then
5076c   On itère sur la valeur du coeff de freinage.
5077c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
5078               gamma0=masse(ig,k)
5079     s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
5080     s         *0.5/larga(ig)
5081     s         *1.
5082c    s         *0.5
5083c              gamma0=0.
5084               zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
5085               zf=0.
5086               zf2=1./(1.-zf)
5087c   la première fois on multiplie le coefficient de freinage
5088c   par le module du vent dans la couche en dessous.
5089               dua=ua(ig,k-1)-u(ig,k-1)
5090               dva=va(ig,k-1)-v(ig,k-1)
5091               do iter=1,5
5092c   On choisit une relaxation lineaire.
5093                  gamma(ig,k)=gamma0
5094c   On choisit une relaxation quadratique.
5095                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
5096                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
5097     s               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))
5098     s               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
5099     s                 +gamma(ig,k))
5100                  va(ig,k)=(fm(ig,k)*va(ig,k-1)
5101     s               +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))
5102     s               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
5103     s                 +gamma(ig,k))
5104c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
5105                  dua=ua(ig,k)-u(ig,k)
5106                  dva=va(ig,k)-v(ig,k)
5107                  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
5108                  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
5109               enddo
5110            else
5111               ua(ig,k)=u(ig,k)
5112               va(ig,k)=v(ig,k)
5113               ue(ig,k)=u(ig,k)
5114               ve(ig,k)=v(ig,k)
5115               gamma(ig,k)=0.
5116            endif
5117         enddo
5118      enddo
5119
5120      do k=2,nlay
5121         do ig=1,ngrid
5122            wud(ig,k)=fm(ig,k)*ue(ig,k)
5123            wvd(ig,k)=fm(ig,k)*ve(ig,k)
5124         enddo
5125      enddo
5126      do ig=1,ngrid
5127         wud(ig,1)=0.
5128         wud(ig,nlay+1)=0.
5129         wvd(ig,1)=0.
5130         wvd(ig,nlay+1)=0.
5131      enddo
5132
5133      do k=1,nlay
5134         do ig=1,ngrid
5135            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
5136     s               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)
5137     s               -wud(ig,k)+wud(ig,k+1))
5138     s               /masse(ig,k)
5139            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
5140     s               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)
5141     s               -wvd(ig,k)+wvd(ig,k+1))
5142     s               /masse(ig,k)
5143         enddo
5144      enddo
5145
5146      return
5147      end
5148      SUBROUTINE thermcell_sec(ngrid,nlay,ptimestep
5149     s                  ,pplay,pplev,pphi,zlev
5150     s                  ,pu,pv,pt,po
5151     s                  ,pduadj,pdvadj,pdtadj,pdoadj
5152     s                  ,fm0,entr0
5153c    s                  ,pu_therm,pv_therm
5154     s                  ,r_aspect,l_mix,w2di,tho)
5155
5156      use dimphy
5157      IMPLICIT NONE
5158
5159c=======================================================================
5160c
5161c   Calcul du transport verticale dans la couche limite en presence
5162c   de "thermiques" explicitement representes
5163c
5164c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
5165c
5166c   le thermique est supposé homogène et dissipé par mélange avec
5167c   son environnement. la longueur l_mix contrôle l'efficacité du
5168c   mélange
5169c
5170c   Le calcul du transport des différentes espèces se fait en prenant
5171c   en compte:
5172c     1. un flux de masse montant
5173c     2. un flux de masse descendant
5174c     3. un entrainement
5175c     4. un detrainement
5176c
5177c=======================================================================
5178
5179c-----------------------------------------------------------------------
5180c   declarations:
5181c   -------------
5182
5183#include "dimensions.h"
5184cccc#include "dimphy.h"
5185#include "YOMCST.h"
5186
5187c   arguments:
5188c   ----------
5189
5190      INTEGER ngrid,nlay,w2di
5191      REAL tho
5192      real ptimestep,l_mix,r_aspect
5193      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
5194      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
5195      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
5196      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
5197      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
5198      real pphi(ngrid,nlay)
5199
5200      integer idetr
5201      save idetr
5202      data idetr/3/
5203c$OMP THREADPRIVATE(idetr)
5204
5205c   local:
5206c   ------
5207
5208      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
5209      real zsortie1d(klon)
5210c CR: on remplace lmax(klon,klev+1)
5211      INTEGER lmax(klon),lmin(klon),lentr(klon)
5212      real linter(klon)
5213      real zmix(klon), fracazmix(klon)
5214c RC
5215      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
5216
5217      real zlev(klon,klev+1),zlay(klon,klev)
5218      REAL zh(klon,klev),zdhadj(klon,klev)
5219      REAL ztv(klon,klev)
5220      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
5221      REAL wh(klon,klev+1)
5222      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
5223      real zla(klon,klev+1)
5224      real zwa(klon,klev+1)
5225      real zld(klon,klev+1)
5226      real zwd(klon,klev+1)
5227      real zsortie(klon,klev)
5228      real zva(klon,klev)
5229      real zua(klon,klev)
5230      real zoa(klon,klev)
5231
5232      real zha(klon,klev)
5233      real wa_moy(klon,klev+1)
5234      real fraca(klon,klev+1)
5235      real fracc(klon,klev+1)
5236      real zf,zf2
5237      real thetath2(klon,klev),wth2(klon,klev)
5238!      common/comtherm/thetath2,wth2
5239
5240      real count_time
5241      integer isplit,nsplit,ialt
5242      parameter (nsplit=10)
5243      data isplit/0/
5244      save isplit
5245c$OMP THREADPRIVATE(isplit)
5246
5247      logical sorties
5248      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
5249      real zpspsk(klon,klev)
5250
5251c     real wmax(klon,klev),wmaxa(klon)
5252      real wmax(klon),wmaxa(klon)
5253      real wa(klon,klev,klev+1)
5254      real wd(klon,klev+1)
5255      real larg_part(klon,klev,klev+1)
5256      real fracd(klon,klev+1)
5257      real xxx(klon,klev+1)
5258      real larg_cons(klon,klev+1)
5259      real larg_detr(klon,klev+1)
5260      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
5261      real pu_therm(klon,klev),pv_therm(klon,klev)
5262      real fm(klon,klev+1),entr(klon,klev)
5263      real fmc(klon,klev+1)
5264
5265cCR:nouvelles variables
5266      real f_star(klon,klev+1),entr_star(klon,klev)
5267      real entr_star_tot(klon),entr_star2(klon)
5268      real f(klon), f0(klon)
5269      real zlevinter(klon)
5270      logical first
5271      data first /.false./
5272      save first
5273c$OMP THREADPRIVATE(first)
5274cRC
5275
5276      character*2 str2
5277      character*10 str10
5278
5279      character (len=20) :: modname='thermcell_sec'
5280      character (len=80) :: abort_message
5281
5282      LOGICAL vtest(klon),down
5283
5284      EXTERNAL SCOPY
5285
5286      integer ncorrec,ll
5287      save ncorrec
5288      data ncorrec/0/
5289c$OMP THREADPRIVATE(ncorrec)
5290     
5291c
5292c-----------------------------------------------------------------------
5293c   initialisation:
5294c   ---------------
5295c
5296       sorties=.true.
5297      IF(ngrid.NE.klon) THEN
5298         PRINT*
5299         PRINT*,'STOP dans convadj'
5300         PRINT*,'ngrid    =',ngrid
5301         PRINT*,'klon  =',klon
5302      ENDIF
5303c
5304c-----------------------------------------------------------------------
5305c   incrementation eventuelle de tendances precedentes:
5306c   ---------------------------------------------------
5307
5308c       print*,'0 OK convect8'
5309
5310      DO 1010 l=1,nlay
5311         DO 1015 ig=1,ngrid
5312            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
5313            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
5314            zu(ig,l)=pu(ig,l)
5315            zv(ig,l)=pv(ig,l)
5316            zo(ig,l)=po(ig,l)
5317            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
53181015     CONTINUE
53191010  CONTINUE
5320
5321c       print*,'1 OK convect8'
5322c                       --------------------
5323c
5324c
5325c                       + + + + + + + + + + +
5326c
5327c
5328c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
5329c  wh,wt,wo ...
5330c
5331c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
5332c
5333c
5334c                       --------------------   zlev(1)
5335c                       \\\\\\\\\\\\\\\\\\\\
5336c
5337c
5338
5339c-----------------------------------------------------------------------
5340c   Calcul des altitudes des couches
5341c-----------------------------------------------------------------------
5342
5343      do l=2,nlay
5344         do ig=1,ngrid
5345            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
5346         enddo
5347      enddo
5348      do ig=1,ngrid
5349         zlev(ig,1)=0.
5350         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
5351      enddo
5352      do l=1,nlay
5353         do ig=1,ngrid
5354            zlay(ig,l)=pphi(ig,l)/RG
5355         enddo
5356      enddo
5357
5358c      print*,'2 OK convect8'
5359c-----------------------------------------------------------------------
5360c   Calcul des densites
5361c-----------------------------------------------------------------------
5362
5363      do l=1,nlay
5364         do ig=1,ngrid
5365            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
5366         enddo
5367      enddo
5368
5369      do l=2,nlay
5370         do ig=1,ngrid
5371            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
5372         enddo
5373      enddo
5374
5375      do k=1,nlay
5376         do l=1,nlay+1
5377            do ig=1,ngrid
5378               wa(ig,k,l)=0.
5379            enddo
5380         enddo
5381      enddo
5382
5383c      print*,'3 OK convect8'
5384c------------------------------------------------------------------
5385c   Calcul de w2, quarre de w a partir de la cape
5386c   a partir de w2, on calcule wa, vitesse de l'ascendance
5387c
5388c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
5389c   w2 est stoke dans wa
5390c
5391c   ATTENTION: dans convect8, on n'utilise le calcule des wa
5392c   independants par couches que pour calculer l'entrainement
5393c   a la base et la hauteur max de l'ascendance.
5394c
5395c   Indicages:
5396c   l'ascendance provenant du niveau k traverse l'interface l avec
5397c   une vitesse wa(k,l).
5398c
5399c                       --------------------
5400c
5401c                       + + + + + + + + + +
5402c
5403c  wa(k,l)   ----       --------------------    l
5404c             /\
5405c            /||\       + + + + + + + + + +
5406c             ||
5407c             ||        --------------------
5408c             ||
5409c             ||        + + + + + + + + + +
5410c             ||
5411c             ||        --------------------
5412c             ||__
5413c             |___      + + + + + + + + + +     k
5414c
5415c                       --------------------
5416c
5417c
5418c
5419c------------------------------------------------------------------
5420
5421cCR: ponderation entrainement des couches instables
5422cdef des entr_star tels que entr=f*entr_star     
5423      do l=1,klev
5424         do ig=1,ngrid
5425            entr_star(ig,l)=0.
5426         enddo
5427      enddo
5428c determination de la longueur de la couche d entrainement
5429      do ig=1,ngrid
5430         lentr(ig)=1
5431      enddo
5432
5433con ne considere que les premieres couches instables
5434      do k=nlay-2,1,-1
5435         do ig=1,ngrid
5436            if (ztv(ig,k).gt.ztv(ig,k+1).and.
5437     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
5438               lentr(ig)=k
5439            endif
5440          enddo
5441      enddo
5442   
5443c determination du lmin: couche d ou provient le thermique
5444      do ig=1,ngrid
5445         lmin(ig)=1
5446      enddo
5447      do ig=1,ngrid
5448         do l=nlay,2,-1
5449            if (ztv(ig,l-1).gt.ztv(ig,l)) then
5450               lmin(ig)=l-1
5451            endif
5452         enddo
5453      enddo
5454c
5455c definition de l'entrainement des couches
5456      do l=1,klev-1
5457         do ig=1,ngrid
5458            if (ztv(ig,l).gt.ztv(ig,l+1).and.
5459     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
5460                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
5461c     s                           (zlev(ig,l+1)-zlev(ig,l))
5462     s                           *sqrt(zlev(ig,l+1))
5463            endif
5464         enddo
5465      enddo
5466c pas de thermique si couche 1 stable
5467      do ig=1,ngrid
5468         if (lmin(ig).gt.1) then
5469            do l=1,klev
5470               entr_star(ig,l)=0.
5471            enddo
5472         endif
5473      enddo
5474c calcul de l entrainement total
5475      do ig=1,ngrid
5476         entr_star_tot(ig)=0.
5477      enddo
5478      do ig=1,ngrid
5479         do k=1,klev
5480            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
5481         enddo
5482      enddo
5483c
5484c      print*,'fin calcul entr_star'
5485      do k=1,klev
5486         do ig=1,ngrid
5487            ztva(ig,k)=ztv(ig,k)
5488         enddo
5489      enddo
5490cRC
5491c      print*,'7 OK convect8'
5492      do k=1,klev+1
5493         do ig=1,ngrid
5494            zw2(ig,k)=0.
5495            fmc(ig,k)=0.
5496cCR
5497            f_star(ig,k)=0.
5498cRC
5499            larg_cons(ig,k)=0.
5500            larg_detr(ig,k)=0.
5501            wa_moy(ig,k)=0.
5502         enddo
5503      enddo
5504
5505c      print*,'8 OK convect8'
5506      do ig=1,ngrid
5507         linter(ig)=1.
5508         lmaxa(ig)=1
5509         lmix(ig)=1
5510         wmaxa(ig)=0.
5511      enddo
5512
5513cCR:
5514      do l=1,nlay-2
5515         do ig=1,ngrid
5516            if (ztv(ig,l).gt.ztv(ig,l+1)
5517     s         .and.entr_star(ig,l).gt.1.e-10
5518     s         .and.zw2(ig,l).lt.1e-10) then
5519               f_star(ig,l+1)=entr_star(ig,l)
5520ctest:calcul de dteta
5521               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
5522     s                     *(zlev(ig,l+1)-zlev(ig,l))
5523     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
5524               larg_detr(ig,l)=0.
5525            else if ((zw2(ig,l).ge.1e-10).and.
5526     s               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
5527               f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
5528               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
5529     s                    *ztv(ig,l))/f_star(ig,l+1)
5530               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
5531     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
5532     s                     *(zlev(ig,l+1)-zlev(ig,l))
5533            endif
5534c determination de zmax continu par interpolation lineaire
5535            if (zw2(ig,l+1).lt.0.) then
5536ctest
5537               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
5538c                  print*,'pb linter'
5539               endif
5540               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
5541     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
5542               zw2(ig,l+1)=0.
5543               lmaxa(ig)=l
5544            else
5545               if (zw2(ig,l+1).lt.0.) then
5546c                  print*,'pb1 zw2<0'
5547               endif
5548               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
5549            endif
5550            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
5551c   lmix est le niveau de la couche ou w (wa_moy) est maximum
5552               lmix(ig)=l+1
5553               wmaxa(ig)=wa_moy(ig,l+1)
5554            endif
5555         enddo
5556      enddo
5557c      print*,'fin calcul zw2'
5558c
5559c Calcul de la couche correspondant a la hauteur du thermique
5560      do ig=1,ngrid
5561         lmax(ig)=lentr(ig)
5562      enddo
5563      do ig=1,ngrid
5564         do l=nlay,lentr(ig)+1,-1
5565            if (zw2(ig,l).le.1.e-10) then
5566               lmax(ig)=l-1
5567            endif
5568         enddo
5569      enddo
5570c pas de thermique si couche 1 stable
5571      do ig=1,ngrid
5572         if (lmin(ig).gt.1) then
5573            lmax(ig)=1
5574            lmin(ig)=1
5575         endif
5576      enddo
5577c   
5578c Determination de zw2 max
5579      do ig=1,ngrid
5580         wmax(ig)=0.
5581      enddo
5582
5583      do l=1,nlay
5584         do ig=1,ngrid
5585            if (l.le.lmax(ig)) then
5586                if (zw2(ig,l).lt.0.)then
5587c                  print*,'pb2 zw2<0'
5588                endif
5589                zw2(ig,l)=sqrt(zw2(ig,l))
5590                wmax(ig)=max(wmax(ig),zw2(ig,l))
5591            else
5592                 zw2(ig,l)=0.
5593            endif
5594          enddo
5595      enddo
5596
5597c   Longueur caracteristique correspondant a la hauteur des thermiques.
5598      do  ig=1,ngrid
5599         zmax(ig)=0.
5600         zlevinter(ig)=zlev(ig,1)
5601      enddo
5602      do  ig=1,ngrid
5603c calcul de zlevinter
5604          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
5605     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
5606     s    -zlev(ig,lmax(ig)))
5607       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
5608      enddo
5609
5610c      print*,'avant fermeture'
5611c Fermeture,determination de f
5612      do ig=1,ngrid
5613         entr_star2(ig)=0.
5614      enddo
5615      do ig=1,ngrid
5616         if (entr_star_tot(ig).LT.1.e-10) then
5617            f(ig)=0.
5618         else
5619             do k=lmin(ig),lentr(ig)
5620                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
5621     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
5622             enddo
5623c Nouvelle fermeture
5624             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
5625     s             *entr_star2(ig))*entr_star_tot(ig)
5626ctest
5627c             if (first) then
5628c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
5629c     s             *wmax(ig))
5630c             endif
5631         endif
5632c         f0(ig)=f(ig)
5633c         first=.true.
5634      enddo
5635c      print*,'apres fermeture'
5636
5637c Calcul de l'entrainement
5638       do k=1,klev
5639         do ig=1,ngrid
5640            entr(ig,k)=f(ig)*entr_star(ig,k)
5641         enddo
5642      enddo
5643cCR:test pour entrainer moins que la masse
5644       do ig=1,ngrid
5645          do l=1,lentr(ig)
5646             if ((entr(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
5647                entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
5648     s                       -0.9*masse(ig,l)/ptimestep
5649                entr(ig,l)=0.9*masse(ig,l)/ptimestep
5650             endif
5651          enddo
5652       enddo
5653cCR: fin test
5654c Calcul des flux
5655      do ig=1,ngrid
5656         do l=1,lmax(ig)-1
5657            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
5658         enddo
5659      enddo
5660
5661cRC
5662
5663
5664c      print*,'9 OK convect8'
5665c     print*,'WA1 ',wa_moy
5666
5667c   determination de l'indice du debut de la mixed layer ou w decroit
5668
5669c   calcul de la largeur de chaque ascendance dans le cas conservatif.
5670c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
5671c   d'une couche est égale à la hauteur de la couche alimentante.
5672c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
5673c   de la vitesse d'entrainement horizontal dans la couche alimentante.
5674
5675      do l=2,nlay
5676         do ig=1,ngrid
5677            if (l.le.lmaxa(ig)) then
5678               zw=max(wa_moy(ig,l),1.e-10)
5679               larg_cons(ig,l)=zmax(ig)*r_aspect
5680     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
5681            endif
5682         enddo
5683      enddo
5684
5685      do l=2,nlay
5686         do ig=1,ngrid
5687            if (l.le.lmaxa(ig)) then
5688c              if (idetr.eq.0) then
5689c  cette option est finalement en dur.
5690                  if ((l_mix*zlev(ig,l)).lt.0.)then
5691c                   print*,'pb l_mix*zlev<0'
5692                  endif
5693cCR: test: nouvelle def de lambda
5694c                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5695                  if (zw2(ig,l).gt.1.e-10) then
5696                  larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
5697                  else
5698                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5699                  endif
5700cRC
5701c              else if (idetr.eq.1) then
5702c                 larg_detr(ig,l)=larg_cons(ig,l)
5703c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
5704c              else if (idetr.eq.2) then
5705c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5706c    s            *sqrt(wa_moy(ig,l))
5707c              else if (idetr.eq.4) then
5708c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5709c    s            *wa_moy(ig,l)
5710c              endif
5711            endif
5712         enddo
5713       enddo
5714
5715c      print*,'10 OK convect8'
5716c     print*,'WA2 ',wa_moy
5717c   calcul de la fraction de la maille concernée par l'ascendance en tenant
5718c   compte de l'epluchage du thermique.
5719c
5720cCR def de  zmix continu (profil parabolique des vitesses)
5721      do ig=1,ngrid
5722           if (lmix(ig).gt.1.) then
5723c test
5724              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
5725     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
5726     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
5727     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
5728     s        then
5729c             
5730            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
5731     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
5732     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
5733     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
5734     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
5735     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
5736     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
5737     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
5738            else
5739            zmix(ig)=zlev(ig,lmix(ig))
5740c            print*,'pb zmix'
5741            endif
5742         else
5743         zmix(ig)=0.
5744         endif
5745ctest
5746         if ((zmax(ig)-zmix(ig)).lt.0.) then
5747            zmix(ig)=0.99*zmax(ig)
5748c            print*,'pb zmix>zmax'
5749         endif
5750      enddo
5751c
5752c calcul du nouveau lmix correspondant
5753      do ig=1,ngrid
5754         do l=1,klev
5755            if (zmix(ig).ge.zlev(ig,l).and.
5756     s          zmix(ig).lt.zlev(ig,l+1)) then
5757              lmix(ig)=l
5758             endif
5759          enddo
5760      enddo
5761c
5762      do l=2,nlay
5763         do ig=1,ngrid
5764            if(larg_cons(ig,l).gt.1.) then
5765c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
5766               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
5767     s            /(r_aspect*zmax(ig))
5768c test
5769               fraca(ig,l)=max(fraca(ig,l),0.)
5770               fraca(ig,l)=min(fraca(ig,l),0.5)
5771               fracd(ig,l)=1.-fraca(ig,l)
5772               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
5773            else
5774c              wa_moy(ig,l)=0.
5775               fraca(ig,l)=0.
5776               fracc(ig,l)=0.
5777               fracd(ig,l)=1.
5778            endif
5779         enddo
5780      enddo                 
5781cCR: calcul de fracazmix
5782       do ig=1,ngrid
5783          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
5784     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
5785     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
5786     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
5787       enddo
5788c
5789       do l=2,nlay
5790          do ig=1,ngrid
5791             if(larg_cons(ig,l).gt.1.) then
5792               if (l.gt.lmix(ig)) then
5793ctest
5794                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
5795c                   print*,'pb xxx'
5796                   xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
5797                 else
5798                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
5799                 endif
5800           if (idetr.eq.0) then
5801               fraca(ig,l)=fracazmix(ig)
5802           else if (idetr.eq.1) then
5803               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
5804           else if (idetr.eq.2) then
5805               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
5806           else
5807               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
5808           endif
5809c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
5810               fraca(ig,l)=max(fraca(ig,l),0.)
5811               fraca(ig,l)=min(fraca(ig,l),0.5)
5812               fracd(ig,l)=1.-fraca(ig,l)
5813               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
5814             endif
5815            endif
5816         enddo
5817      enddo
5818     
5819c      print*,'fin calcul fraca'
5820c      print*,'11 OK convect8'
5821c     print*,'Ea3 ',wa_moy
5822c------------------------------------------------------------------
5823c   Calcul de fracd, wd
5824c   somme wa - wd = 0
5825c------------------------------------------------------------------
5826
5827
5828      do ig=1,ngrid
5829         fm(ig,1)=0.
5830         fm(ig,nlay+1)=0.
5831      enddo
5832
5833      do l=2,nlay
5834           do ig=1,ngrid
5835              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
5836cCR:test
5837              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
5838     s            .and.l.gt.lmix(ig)) then
5839                 fm(ig,l)=fm(ig,l-1)
5840c                 write(1,*)'ajustement fm, l',l
5841              endif
5842c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
5843cRC
5844           enddo
5845         do ig=1,ngrid
5846            if(fracd(ig,l).lt.0.1) then
5847              abort_message = 'fracd trop petit'
5848              CALL abort_gcm (modname,abort_message,1)
5849            else
5850c    vitesse descendante "diagnostique"
5851               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
5852            endif
5853         enddo
5854      enddo
5855
5856      do l=1,nlay
5857         do ig=1,ngrid
5858c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
5859            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
5860         enddo
5861      enddo
5862
5863c      print*,'12 OK convect8'
5864c     print*,'WA4 ',wa_moy
5865cc------------------------------------------------------------------
5866c   calcul du transport vertical
5867c------------------------------------------------------------------
5868
5869      go to 4444
5870c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
5871      do l=2,nlay-1
5872         do ig=1,ngrid
5873            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
5874     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
5875c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
5876c    s         ,fm(ig,l+1)*ptimestep
5877c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
5878             endif
5879         enddo
5880      enddo
5881
5882      do l=1,nlay
5883         do ig=1,ngrid
5884            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
5885c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
5886c    s         ,entr(ig,l)*ptimestep
5887c    s         ,'   M=',masse(ig,l)
5888             endif
5889         enddo
5890      enddo
5891
5892      do l=1,nlay
5893         do ig=1,ngrid
5894            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
5895c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
5896c    s         ,'   FM=',fm(ig,l)
5897            endif
5898            if(.not.masse(ig,l).ge.1.e-10
5899     s         .or..not.masse(ig,l).le.1.e4) then
5900c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
5901c    s         ,'   M=',masse(ig,l)
5902c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
5903c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
5904c     print*,'zlev(ig,l+1),zlev(ig,l)'
5905c    s                ,zlev(ig,l+1),zlev(ig,l)
5906c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
5907c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
5908            endif
5909            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
5910c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
5911c    s         ,'   E=',entr(ig,l)
5912            endif
5913         enddo
5914      enddo
5915
59164444   continue
5917
5918cCR:redefinition du entr
5919       do l=1,nlay
5920         do ig=1,ngrid
5921            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
5922            if (detr(ig,l).lt.0.) then
5923                entr(ig,l)=entr(ig,l)-detr(ig,l)
5924                detr(ig,l)=0.
5925c     print*,'WARNING !!! detrainement negatif ',ig,l
5926            endif
5927         enddo
5928      enddo
5929cRC
5930      if (w2di.eq.1) then
5931         fm0=fm0+ptimestep*(fm-fm0)/tho
5932         entr0=entr0+ptimestep*(entr-entr0)/tho
5933      else
5934         fm0=fm
5935         entr0=entr
5936      endif
5937
5938      if (1.eq.1) then
5939         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5940     .    ,zh,zdhadj,zha)
5941         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5942     .    ,zo,pdoadj,zoa)
5943      else
5944         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
5945     .    ,zh,zdhadj,zha)
5946         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
5947     .    ,zo,pdoadj,zoa)
5948      endif
5949
5950      if (1.eq.0) then
5951         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
5952     .    ,fraca,zmax
5953     .    ,zu,zv,pduadj,pdvadj,zua,zva)
5954      else
5955         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5956     .    ,zu,pduadj,zua)
5957         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5958     .    ,zv,pdvadj,zva)
5959      endif
5960
5961      do l=1,nlay
5962         do ig=1,ngrid
5963            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
5964            zf2=zf/(1.-zf)
5965            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
5966            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
5967         enddo
5968      enddo
5969
5970
5971
5972c     print*,'13 OK convect8'
5973c     print*,'WA5 ',wa_moy
5974      do l=1,nlay
5975         do ig=1,ngrid
5976            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
5977         enddo
5978      enddo
5979
5980
5981c     do l=1,nlay
5982c        do ig=1,ngrid
5983c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
5984c     print*,'WARN!!! ig=',ig,'  l=',l
5985c    s         ,'   pdtadj=',pdtadj(ig,l)
5986c           endif
5987c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
5988c     print*,'WARN!!! ig=',ig,'  l=',l
5989c    s         ,'   pdoadj=',pdoadj(ig,l)
5990c           endif
5991c        enddo
5992c      enddo
5993
5994c      print*,'14 OK convect8'
5995c------------------------------------------------------------------
5996c   Calculs pour les sorties
5997c------------------------------------------------------------------
5998
5999      if(sorties) then
6000      do l=1,nlay
6001         do ig=1,ngrid
6002            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
6003            zld(ig,l)=fracd(ig,l)*zmax(ig)
6004            if(1.-fracd(ig,l).gt.1.e-10)
6005     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
6006         enddo
6007      enddo
6008
6009cdeja fait
6010c      do l=1,nlay
6011c         do ig=1,ngrid
6012c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
6013c            if (detr(ig,l).lt.0.) then
6014c                entr(ig,l)=entr(ig,l)-detr(ig,l)
6015c                detr(ig,l)=0.
6016c     print*,'WARNING !!! detrainement negatif ',ig,l
6017c            endif
6018c         enddo
6019c      enddo
6020
6021c     print*,'15 OK convect8'
6022
6023      isplit=isplit+1
6024
6025
6026c #define und
6027        goto 123
6028#ifdef und
6029      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
6030      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
6031      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
6032      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
6033      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
6034      CALL writeg1d(1,nlay,zla,'la      ','la      ')
6035      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
6036      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
6037      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
6038      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
6039      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
6040      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
6041      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
6042      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
6043      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
6044      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
6045      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
6046      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
6047      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
6048      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
6049      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
6050      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
6051      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
6052      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
6053
6054      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
6055      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
6056      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
6057
6058c   recalcul des flux en diagnostique...
6059c     print*,'PAS DE TEMPS ',ptimestep
6060       call dt2F(pplev,pplay,pt,pdtadj,wh)
6061      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
6062#endif
6063123   continue
6064! #define troisD
6065#ifdef troisD
6066c       if (sorties) then
6067      print*,'Debut des wrgradsfi'
6068
6069c      print*,'16 OK convect8'
6070         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
6071         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
6072         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
6073         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
6074         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
6075         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
6076c      print*,'WA6 ',wa_moy
6077         call wrgradsfi(1,nlay,zla,'la        ','la        ')
6078         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
6079         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
6080         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
6081         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
6082         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
6083         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
6084         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
6085         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
6086         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
6087         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
6088         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
6089         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
6090         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
6091         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
6092         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
6093         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
6094         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
6095         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
6096         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
6097         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
6098         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
6099         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
6100         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
6101         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
6102         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
6103
6104         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
6105         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
6106         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
6107
6108cCR:nouveaux diagnostiques
6109      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
6110      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
6111      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
6112      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
6113      zsortie1d(:)=lmax(:)
6114      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
6115      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
6116      zsortie1d(:)=lmix(:)
6117      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
6118      zsortie1d(:)=lentr(:)
6119      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
6120
6121c      print*,'17 OK convect8'
6122
6123         do k=1,klev/10
6124            write(str2,'(i2.2)') k
6125            str10='wa'//str2
6126            do l=1,nlay
6127               do ig=1,ngrid
6128                  zsortie(ig,l)=wa(ig,k,l)
6129               enddo
6130            enddo   
6131            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
6132            do l=1,nlay
6133               do ig=1,ngrid
6134                  zsortie(ig,l)=larg_part(ig,k,l)
6135               enddo
6136            enddo
6137            str10='la'//str2
6138            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
6139         enddo
6140
6141
6142c     print*,'18 OK convect8'
6143c      endif
6144      print*,'Fin des wrgradsfi'
6145#endif
6146
6147      endif
6148
6149c     if(wa_moy(1,4).gt.1.e-10) stop
6150
6151c      print*,'19 OK convect8'
6152      return
6153      end
6154
Note: See TracBrowser for help on using the repository browser.