source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_old.F90 @ 5428

Last change on this file since 5428 was 5285, checked in by abarral, 2 months ago

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