source: LMDZ6/branches/contrails/libf/phylmd/lmdz_thermcell_old.f90 @ 5456

Last change on this file since 5456 was 5450, checked in by aborella, 2 days ago

Merge with trunk

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