source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/thermcell_old.F90 @ 2272

Last change on this file since 2272 was 1999, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1920:1997 into testing branch

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