source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_old.F90 @ 5101

Last change on this file since 5101 was 5101, checked in by abarral, 4 months ago

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