source: lmdz_wrf/WRFV3/lmdz/thermcell_old.F90 @ 1

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

WRF: version v3.3
LMDZ: version v1818

More details in:

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