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