source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_old.f90 @ 5506

Last change on this file since 5506 was 5501, checked in by fhourdin, 3 days ago

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