source: LMDZ4/trunk/libf/phylmd/thermcell_old.F @ 1436

Last change on this file since 1436 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files


Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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