source: LMDZ6/branches/blowing_snow/libf/phylmd/thermcell_old.F90 @ 5229

Last change on this file since 5229 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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