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