My Project
 All Classes Files Functions Variables Macros
thermcell_main.F90
Go to the documentation of this file.
1 !
2 ! $Id: thermcell_main.F90 1738 2013-04-02 09:48:33Z idelkadi $
3 !
4  SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep &
5  & ,pplay,pplev,pphi,debut &
6  & ,pu,pv,pt,po &
7  & ,pduadj,pdvadj,pdtadj,pdoadj &
8  & ,fm0,entr0,detr0,zqta,zqla,lmax &
9  & ,ratqscth,ratqsdiff,zqsatth &
10  & ,ale_bl,alp_bl,lalim_conv,wght_th &
11  & ,zmax0, f0,zw2,fraca,ztv &
12  & ,zpspsk,ztla,zthl &
13 !!! nrlmd le 10/04/2012
14  & ,pbl_tke,pctsrf,omega,airephy &
15  & ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
16  & ,n2,s2,ale_bl_stat &
17  & ,therm_tke_max,env_tke_max &
18  & ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
19  & ,alp_bl_conv,alp_bl_stat &
20 !!! fin nrlmd le 10/04/2012
21  & )
22 
23  USE dimphy
24  USE ioipsl
25  USE comgeomphy , ONLY:rlond,rlatd
26  IMPLICIT NONE
27 
28 !=======================================================================
29 ! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
30 ! Version du 09.02.07
31 ! Calcul du transport vertical dans la couche limite en presence
32 ! de "thermiques" explicitement representes avec processus nuageux
33 !
34 ! Reecriture a partir d'un listing papier a Habas, le 14/02/00
35 !
36 ! le thermique est suppose homogene et dissipe par melange avec
37 ! son environnement. la longueur l_mix controle l'efficacite du
38 ! melange
39 !
40 ! Le calcul du transport des differentes especes se fait en prenant
41 ! en compte:
42 ! 1. un flux de masse montant
43 ! 2. un flux de masse descendant
44 ! 3. un entrainement
45 ! 4. un detrainement
46 !
47 ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
48 ! Introduction of an implicit computation of vertical advection in
49 ! the environment of thermal plumes in thermcell_dq
50 ! impl = 0 : explicit, 1 : implicit, -1 : old version
51 ! controled by iflag_thermals =
52 ! 15, 16 run with impl=-1 : numerical convergence with NPv3
53 ! 17, 18 run with impl=1 : more stable
54 ! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
55 !
56 !=======================================================================
57 
58 
59 !-----------------------------------------------------------------------
60 ! declarations:
61 ! -------------
62 
63 #include "dimensions.h"
64 #include "YOMCST.h"
65 #include "YOETHF.h"
66 #include "FCTTRE.h"
67 #include "iniprint.h"
68 #include "thermcell.h"
69 !!! nrlmd le 10/04/2012
70 #include "indicesol.h"
71 !!! fin nrlmd le 10/04/2012
72 
73 ! arguments:
74 ! ----------
75 
76 !IM 140508
77  INTEGER itap
78 
79  INTEGER ngrid,nlay
80  real ptimestep
81  REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
82  REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
83  REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
84  REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
85  REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
86  real pphi(ngrid,nlay)
87 
88 ! local:
89 ! ------
90 
91  integer icount
92 
93  integer, save :: dvdq=1,dqimpl=-1
94 !$OMP THREADPRIVATE(dvdq,dqimpl)
95  data icount/0/
96  save icount
97 !$OMP THREADPRIVATE(icount)
98 
99  integer,save :: igout=1
100 !$OMP THREADPRIVATE(igout)
101  integer,save :: lunout1=6
102 !$OMP THREADPRIVATE(lunout1)
103  integer,save :: lev_out=10
104 !$OMP THREADPRIVATE(lev_out)
105 
106  REAL susqr2pi, reuler
107 
108  INTEGER ig,k,l,ll,ierr
109  real zsortie1d(klon)
110  INTEGER lmax(klon),lmin(klon),lalim(klon)
111  INTEGER lmix(klon)
112  INTEGER lmix_bis(klon)
113  real linter(klon)
114  real zmix(klon)
115  real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
116 ! real fraca(klon,klev)
117 
118  real zmax_sec(klon)
119 !on garde le zmax du pas de temps precedent
120  real zmax0(klon)
121 !FH/IM save zmax0
122 
123  real lambda
124 
125  real zlev(klon,klev+1),zlay(klon,klev)
126  real deltaz(klon,klev)
127  REAL zh(klon,klev)
128  real zthl(klon,klev),zdthladj(klon,klev)
129  REAL ztv(klon,klev)
130  real zu(klon,klev),zv(klon,klev),zo(klon,klev)
131  real zl(klon,klev)
132  real zsortie(klon,klev)
133  real zva(klon,klev)
134  real zua(klon,klev)
135  real zoa(klon,klev)
136 
137  real zta(klon,klev)
138  real zha(klon,klev)
139  real fraca(klon,klev+1)
140  real zf,zf2
141  real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
142  real q2(klon,klev)
143 ! FH probleme de dimensionnement avec l'allocation dynamique
144 ! common/comtherm/thetath2,wth2
145  real wq(klon,klev)
146  real wthl(klon,klev)
147  real wthv(klon,klev)
148 
149  real ratqscth(klon,klev)
150  real var
151  real vardiff
152  real ratqsdiff(klon,klev)
153 
154  logical sorties
155  real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
156  real zpspsk(klon,klev)
157 
158  real wmax(klon)
159  real wmax_tmp(klon)
160  real wmax_sec(klon)
161  real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
162  real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
163 
164  real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
165 !niveau de condensation
166  integer nivcon(klon)
167  real zcon(klon)
168  REAL chi
169  real zcon2(klon)
170  real pcon(klon)
171  real zqsat(klon,klev)
172  real zqsatth(klon,klev)
173 
174  real f_star(klon,klev+1),entr_star(klon,klev)
175  real detr_star(klon,klev)
176  real alim_star_tot(klon)
177  real alim_star(klon,klev)
178  real alim_star_clos(klon,klev)
179  real f(klon), f0(klon)
180 !FH/IM save f0
181  real zlevinter(klon)
182  logical debut
183  real seuil
184  real csc(klon,klev)
185 
186 !!! nrlmd le 10/04/2012
187 
188 !------Entrées
189  real pbl_tke(klon,klev+1,nbsrf)
190  real pctsrf(klon,nbsrf)
191  real omega(klon,klev)
192  real airephy(klon)
193 !------Sorties
194  real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
195  real therm_tke_max0(klon),env_tke_max0(klon)
196  real n2(klon),s2(klon)
197  real ale_bl_stat(klon)
198  real therm_tke_max(klon,klev),env_tke_max(klon,klev)
199  real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
200 !------Local
201  integer nsrf
202  real rhobarz0(klon) ! Densité au LCL
203  logical ok_lcl(klon) ! Existence du LCL des thermiques
204  integer klcl(klon) ! Niveau du LCL
205  real interp(klon) ! Coef d'interpolation pour le LCL
206 !--Triggering
207  real su ! Surface unité: celle d'un updraft élémentaire
208  parameter(su=4e4)
209  real hcoef ! Coefficient directeur pour le calcul de s2
210  parameter(hcoef=1)
211  real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
212  parameter(hmincoef=0.3)
213  real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
214  parameter(eps1=0.3)
215  real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2
216  real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
217  real zmax_moy_coef
218  parameter(zmax_moy_coef=0.33)
219  real depth(klon) ! Epaisseur moyenne du cumulus
220  real w_max(klon) ! Vitesse max statistique
221  real s_max(klon)
222 !--Closure
223  real pbl_tke_max(klon,klev) ! Profil de TKE moyenne
224  real pbl_tke_max0(klon) ! TKE moyenne au LCL
225  real w_ls(klon,klev) ! Vitesse verticale grande échelle (m/s)
226  real coef_m ! On considère un rendement pour alp_bl_fluct_m
227  parameter(coef_m=1.)
228  real coef_tke ! On considère un rendement pour alp_bl_fluct_tke
229  parameter(coef_tke=1.)
230 
231 !!! fin nrlmd le 10/04/2012
232 
233 !
234  !nouvelles variables pour la convection
235  real ale_bl(klon)
236  real alp_bl(klon)
237  real alp_int(klon),dp_int(klon),zdp
238  real ale_int(klon)
239  integer n_int(klon)
240  real fm_tot(klon)
241  real wght_th(klon,klev)
242  integer lalim_conv(klon)
243 !v1d logical therm
244 !v1d save therm
245 
246  character*2 str2
247  character*10 str10
248 
249  character (len=20) :: modname='thermcell_main'
250  character (len=80) :: abort_message
251 
252  EXTERNAL scopy
253 !
254 
255 !-----------------------------------------------------------------------
256 ! initialisation:
257 ! ---------------
258 !
259 
260  seuil=0.25
261 
262  if (debut) then
263 ! call getin('dvdq',dvdq)
264 ! call getin('dqimpl',dqimpl)
265 
266  if (iflag_thermals==15.or.iflag_thermals==16) then
267  dvdq=0
268  dqimpl=-1
269  else
270  dvdq=1
271  dqimpl=1
272  endif
273 
274  fm0=0.
275  entr0=0.
276  detr0=0.
277 
278 
279 #undef wrgrads_thermcell
280 #ifdef wrgrads_thermcell
281 ! Initialisation des sorties grads pour les thermiques.
282 ! Pour l'instant en 1D sur le point igout.
283 ! Utilise par thermcell_out3d.h
284  str10='therm'
285  call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
286  & rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
287  & ptimestep,str10,'therm ')
288 #endif
289 
290 
291 
292  endif
293 
294  fm=0. ; entr=0. ; detr=0.
295 
296 
297  icount=icount+1
298 
299 !IM 090508 beg
300 !print*,'====================================================================='
301 !print*,'====================================================================='
302 !print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
303 !print*,'====================================================================='
304 !print*,'====================================================================='
305 !IM 090508 end
306 
307  if (prt_level.ge.1) print*,'thermcell_main V4'
308 
309  sorties=.true.
310  IF(ngrid.NE.klon) THEN
311  print*
312  print*,'STOP dans convadj'
313  print*,'ngrid =',ngrid
314  print*,'klon =',klon
315  ENDIF
316 !
317 ! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
318  do ig=1,klon
319  f0(ig)=max(f0(ig),1.e-2)
320  zmax0(ig)=max(zmax0(ig),40.)
321 !IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
322  enddo
323 
324  if (prt_level.ge.20) then
325  do ig=1,ngrid
326  print*,'th_main ig f0',ig,f0(ig)
327  enddo
328  endif
329 !-----------------------------------------------------------------------
330 ! Calcul de T,q,ql a partir de Tl et qT dans l environnement
331 ! --------------------------------------------------------------------
332 !
333  CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
334  & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
335 
336  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
337 
338 !------------------------------------------------------------------------
339 ! --------------------
340 !
341 !
342 ! + + + + + + + + + + +
343 !
344 !
345 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
346 ! wh,wt,wo ...
347 !
348 ! + + + + + + + + + + + zh,zu,zv,zo,rho
349 !
350 !
351 ! -------------------- zlev(1)
352 ! \\\\\\\\\\\\\\\\\\\\
353 !
354 !
355 
356 !-----------------------------------------------------------------------
357 ! Calcul des altitudes des couches
358 !-----------------------------------------------------------------------
359 
360  do l=2,nlay
361  zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/rg
362  enddo
363  zlev(:,1)=0.
364  zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/rg
365  do l=1,nlay
366  zlay(:,l)=pphi(:,l)/rg
367  enddo
368 !calcul de l epaisseur des couches
369  do l=1,nlay
370  deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
371  enddo
372 
373 ! print*,'2 OK convect8'
374 !-----------------------------------------------------------------------
375 ! Calcul des densites
376 !-----------------------------------------------------------------------
377 
378  rho(:,:)=pplay(:,:)/(zpspsk(:,:)*rd*ztv(:,:))
379 
380  if (prt_level.ge.10)write(lunout,*) &
381  & 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
382  rhobarz(:,1)=rho(:,1)
383 
384  do l=2,nlay
385  rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
386  enddo
387 
388 !calcul de la masse
389  do l=1,nlay
390  masse(:,l)=(pplev(:,l)-pplev(:,l+1))/rg
391  enddo
392 
393  if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
394 
395 !------------------------------------------------------------------
396 !
397 ! /|\
398 ! -------- | F_k+1 -------
399 ! ----> D_k
400 ! /|\ <---- E_k , A_k
401 ! -------- | F_k ---------
402 ! ----> D_k-1
403 ! <---- E_k-1 , A_k-1
404 !
405 !
406 !
407 !
408 !
409 ! ---------------------------
410 !
411 ! ----- F_lmax+1=0 ---------- \
412 ! lmax (zmax) |
413 ! --------------------------- |
414 ! |
415 ! --------------------------- |
416 ! |
417 ! --------------------------- |
418 ! |
419 ! --------------------------- |
420 ! |
421 ! --------------------------- |
422 ! | E
423 ! --------------------------- | D
424 ! |
425 ! --------------------------- |
426 ! |
427 ! --------------------------- \ |
428 ! lalim | |
429 ! --------------------------- | |
430 ! | |
431 ! --------------------------- | |
432 ! | A |
433 ! --------------------------- | |
434 ! | |
435 ! --------------------------- | |
436 ! lmin (=1 pour le moment) | |
437 ! ----- F_lmin=0 ------------ / /
438 !
439 ! ---------------------------
440 ! //////////////////////////
441 !
442 !
443 !=============================================================================
444 ! Calculs initiaux ne faisant pas intervenir les changements de phase
445 !=============================================================================
446 
447 !------------------------------------------------------------------
448 ! 1. alim_star est le profil vertical de l'alimentation a la base du
449 ! panache thermique, calcule a partir de la flotabilite de l'air sec
450 ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
451 !------------------------------------------------------------------
452 !
453  entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
454  lmin=1
455 
456 !-----------------------------------------------------------------------------
457 ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
458 ! panache sec conservatif (e=d=0) alimente selon alim_star
459 ! Il s'agit d'un calcul de type CAPE
460 ! zmax_sec est utilise pour determiner la geometrie du thermique.
461 !------------------------------------------------------------------------------
462 !---------------------------------------------------------------------------------
463 !calcul du melange et des variables dans le thermique
464 !--------------------------------------------------------------------------------
465 !
466  if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
467 !IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
468 
469 ! Gestion temporaire de plusieurs appels à thermcell_plume au travers
470 ! de la variable iflag_thermals
471 
472 ! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
473  if (iflag_thermals_ed<=9) then
474 ! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
475  CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
476  & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
477  & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
478  & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
479  & ,lev_out,lunout1,igout)
480 
481  elseif (iflag_thermals_ed>9) then
482 ! print*,'THERM RIO et al 2010, version d Arnaud'
483  CALL thermcellv1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
484  & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
485  & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
486  & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
487  & ,lev_out,lunout1,igout)
488 
489  endif
490 
491  if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
492 
493  call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
494  call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
495 
496  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
497  if (prt_level.ge.10) then
498  write(lunout1,*) 'Dans thermcell_main 2'
499  write(lunout1,*) 'lmin ',lmin(igout)
500  write(lunout1,*) 'lalim ',lalim(igout)
501  write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
502  write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
503  & ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
504  endif
505 
506 !-------------------------------------------------------------------------------
507 ! Calcul des caracteristiques du thermique:zmax,zmix,wmax
508 !-------------------------------------------------------------------------------
509 !
510  CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
511  & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
512 ! Attention, w2 est transforme en sa racine carree dans cette routine
513 ! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
514  wmax_tmp=0.
515  do l=1,nlay
516  wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
517  enddo
518 ! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
519 
520 
521 
522  call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
523  call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
524  call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
525  call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
526 
527  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
528 
529 !-------------------------------------------------------------------------------
530 ! Fermeture,determination de f
531 !-------------------------------------------------------------------------------
532 !
533 !
534 !! write(lunout,*)'THERM NOUVEAU XXXXX'
535  CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
536  & lalim,lmin,zmax_sec,wmax_sec,lev_out)
537 
538 call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
539 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
540 
541  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
542  if (prt_level.ge.10) then
543  write(lunout1,*) 'Dans thermcell_main 1b'
544  write(lunout1,*) 'lmin ',lmin(igout)
545  write(lunout1,*) 'lalim ',lalim(igout)
546  write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
547  write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
548  & ,l=1,lalim(igout)+4)
549  endif
550 
551 
552 
553 
554 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
555 ! Apparemment sans importance
556  alim_star_clos(:,:)=alim_star(:,:)
557  alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
558 
559 ! Appel avec la version seche
560  CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
561  & zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
562 
563 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
564 ! Appel avec les zmax et wmax tenant compte de la condensation
565 ! Semble moins bien marcher
566 ! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
567 ! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
568 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
569 
570  if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
571 
572  if (tau_thermals>1.) then
573  lambda=exp(-ptimestep/tau_thermals)
574  f0=(1.-lambda)*f+lambda*f0
575  else
576  f0=f
577  endif
578 
579 ! Test valable seulement en 1D mais pas genant
580  if (.not. (f0(1).ge.0.) ) then
581  abort_message = .not..ge.' (f0(1)0.)'
582  CALL abort_gcm(modname,abort_message,1)
583  endif
584 
585 !-------------------------------------------------------------------------------
586 !deduction des flux
587 !-------------------------------------------------------------------------------
588 
589  CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
590  & lalim,lmax,alim_star, &
591  & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
592  & detr,zqla,lev_out,lunout1,igout)
593 !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
594 
595  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
596  call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
597  call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
598 
599 !------------------------------------------------------------------
600 ! On ne prend pas directement les profils issus des calculs precedents
601 ! mais on s'autorise genereusement une relaxation vers ceci avec
602 ! une constante de temps tau_thermals (typiquement 1800s).
603 !------------------------------------------------------------------
604 
605  if (tau_thermals>1.) then
606  lambda=exp(-ptimestep/tau_thermals)
607  fm0=(1.-lambda)*fm+lambda*fm0
608  entr0=(1.-lambda)*entr+lambda*entr0
609  detr0=(1.-lambda)*detr+lambda*detr0
610  else
611  fm0=fm
612  entr0=entr
613  detr0=detr
614  endif
615 
616 !c------------------------------------------------------------------
617 ! calcul du transport vertical
618 !------------------------------------------------------------------
619 
620  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
621  & zthl,zdthladj,zta,lev_out)
622  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
623  & po,pdoadj,zoa,lev_out)
624 
625 !------------------------------------------------------------------
626 ! Calcul de la fraction de l'ascendance
627 !------------------------------------------------------------------
628  do ig=1,klon
629  fraca(ig,1)=0.
630  fraca(ig,nlay+1)=0.
631  enddo
632  do l=2,nlay
633  do ig=1,klon
634  if (zw2(ig,l).gt.1.e-10) then
635  fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
636  else
637  fraca(ig,l)=0.
638  endif
639  enddo
640  enddo
641 
642 !------------------------------------------------------------------
643 ! calcul du transport vertical du moment horizontal
644 !------------------------------------------------------------------
645 
646 !IM 090508
647  if (dvdq == 0 ) then
648 
649 ! Calcul du transport de V tenant compte d'echange par gradient
650 ! de pression horizontal avec l'environnement
651 
652  call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
653 ! & ,fraca*dvdq,zmax &
654  & ,fraca,zmax &
655  & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
656 
657  else
658 
659 ! calcul purement conservatif pour le transport de V
660  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
661  & ,zu,pduadj,zua,lev_out)
662  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
663  & ,zv,pdvadj,zva,lev_out)
664 
665  endif
666 
667 ! print*,'13 OK convect8'
668  do l=1,nlay
669  do ig=1,ngrid
670  pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
671  enddo
672  enddo
673 
674  if (prt_level.ge.1) print*,'14 OK convect8'
675 !------------------------------------------------------------------
676 ! Calculs de diagnostiques pour les sorties
677 !------------------------------------------------------------------
678 !calcul de fraca pour les sorties
679 
680  if (sorties) then
681  if (prt_level.ge.1) print*,'14a OK convect8'
682 ! calcul du niveau de condensation
683 ! initialisation
684  do ig=1,ngrid
685  nivcon(ig)=0
686  zcon(ig)=0.
687  enddo
688 !nouveau calcul
689  do ig=1,ngrid
690  chi=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
691  pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**chi
692  enddo
693 !IM do k=1,nlay
694  do k=1,nlay-1
695  do ig=1,ngrid
696  if ((pcon(ig).le.pplay(ig,k)) &
697  & .and.(pcon(ig).gt.pplay(ig,k+1))) then
698  zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(rg*rho(ig,k))/100.
699  endif
700  enddo
701  enddo
702 !IM
703  ierr=0
704  do ig=1,ngrid
705  if (pcon(ig).le.pplay(ig,nlay)) then
706  zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(rg*rho(ig,nlay))/100.
707  ierr=1
708  endif
709  enddo
710  if (ierr==1) then
711  abort_message = 'thermcellV0_main: les thermiques vont trop haut '
712  CALL abort_gcm(modname,abort_message,1)
713  endif
714 
715  if (prt_level.ge.1) print*,'14b OK convect8'
716  do k=nlay,1,-1
717  do ig=1,ngrid
718  if (zqla(ig,k).gt.1e-10) then
719  nivcon(ig)=k
720  zcon(ig)=zlev(ig,k)
721  endif
722  enddo
723  enddo
724  if (prt_level.ge.1) print*,'14c OK convect8'
725 !calcul des moments
726 !initialisation
727  do l=1,nlay
728  do ig=1,ngrid
729  q2(ig,l)=0.
730  wth2(ig,l)=0.
731  wth3(ig,l)=0.
732  ratqscth(ig,l)=0.
733  ratqsdiff(ig,l)=0.
734  enddo
735  enddo
736  if (prt_level.ge.1) print*,'14d OK convect8'
737  if (prt_level.ge.10)write(lunout,*) &
738  & 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
739  do l=1,nlay
740  do ig=1,ngrid
741  zf=fraca(ig,l)
742  zf2=zf/(1.-zf)
743 !
744  thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
745  if(zw2(ig,l).gt.1.e-10) then
746  wth2(ig,l)=zf2*(zw2(ig,l))**2
747  else
748  wth2(ig,l)=0.
749  endif
750  wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) &
751  & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
752  q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
753 !test: on calcul q2/po=ratqsc
754  ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
755  enddo
756  enddo
757 !calcul des flux: q, thetal et thetav
758  do l=1,nlay
759  do ig=1,ngrid
760  wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
761  wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
762  wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
763  enddo
764  enddo
765 !
766 
767 !!! nrlmd le 10/04/2012
768 
769 !------------Test sur le LCL des thermiques
770  do ig=1,ngrid
771  ok_lcl(ig)=.false.
772  if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
773  enddo
774 
775 !------------Localisation des niveaux entourant le LCL et du coef d'interpolation
776  do l=1,nlay-1
777  do ig=1,ngrid
778  if (ok_lcl(ig)) then
779  if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
780  klcl(ig)=l
781  interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
782  endif
783  endif
784  enddo
785  enddo
786 
787 !------------Hauteur des thermiques
788 !!jyg le 27/04/2012
789 !! do ig =1,ngrid
790 !! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
791 !! & -rhobarz(ig,klcl(ig)))*interp(ig)
792 !! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
793 !! zmax(ig)=pphi(ig,lmax(ig))/rg
794 !! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
795 !! enddo
796  do ig =1,ngrid
797  zmax(ig)=pphi(ig,lmax(ig))/rg
798  if (ok_lcl(ig)) then
799  rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
800  & -rhobarz(ig,klcl(ig)))*interp(ig)
801  zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*rg)
802  zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
803  else
804  rhobarz0(ig)=0.
805  zlcl(ig)=zmax(ig)
806  endif
807  enddo
808 !!jyg fin
809 
810 !------------Calcul des propriétés du thermique au LCL
811  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
812 
813  !-----Initialisation de la TKE moyenne
814  do l=1,nlay
815  do ig=1,ngrid
816  pbl_tke_max(ig,l)=0.
817  enddo
818  enddo
819 
820 !-----Calcul de la TKE moyenne
821  do nsrf=1,nbsrf
822  do l=1,nlay
823  do ig=1,ngrid
824  pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
825  enddo
826  enddo
827  enddo
828 
829 !-----Initialisations des TKE dans et hors des thermiques
830  do l=1,nlay
831  do ig=1,ngrid
832  therm_tke_max(ig,l)=pbl_tke_max(ig,l)
833  env_tke_max(ig,l)=pbl_tke_max(ig,l)
834  enddo
835  enddo
836 
837 !-----Calcul de la TKE transportée par les thermiques : therm_tke_max
838  call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
839  & rg,pplev,therm_tke_max)
840 ! print *,' thermcell_tke_transport -> ' !!jyg
841 
842 !-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
843  do l=1,nlay
844  do ig=1,ngrid
845  pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l) ! Recalcul de TKE moyenne aprés transport de TKE_TH
846  env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l)) ! Recalcul de TKE dans l'environnement aprés transport de TKE_TH
847  w_ls(ig,l)=-1.*omega(ig,l)/(rg*rhobarz(ig,l)) ! Vitesse verticale de grande échelle
848  enddo
849  enddo
850 ! print *,' apres w_ls = ' !!jyg
851 
852  do ig=1,ngrid
853  if (ok_lcl(ig)) then
854  fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
855  & -fraca(ig,klcl(ig)))*interp(ig)
856  w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
857  & -zw2(ig,klcl(ig)))*interp(ig)
858  w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
859  & -w_ls(ig,klcl(ig)))*interp(ig)
860  therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
861  & +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
862  env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
863  & -env_tke_max(ig,klcl(ig)))*interp(ig)
864  pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
865  & -pbl_tke_max(ig,klcl(ig)))*interp(ig)
866  if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
867  if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
868  if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
869  else
870  fraca0(ig)=0.
871  w0(ig)=0.
872 !!jyg le 27/04/2012
873 !! zlcl(ig)=0.
874 !!
875  endif
876  enddo
877 
878  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
879 ! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
880 
881 !------------Triggering------------------
882  IF (iflag_trig_bl.ge.1) THEN
883 
884 !-----Initialisations
885  depth(:)=0.
886  n2(:)=0.
887  s2(:)=0.
888  s_max(:)=0.
889 
890 !-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
891  do ig=1,ngrid
892  zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
893  depth(ig)=zmax_moy(ig)-zlcl(ig)
894  hmin(ig)=hmincoef*zlcl(ig)
895  if (depth(ig).ge.10.) then
896  s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
897  n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
898 !!
899 !!jyg le 27/04/2012
900 !! s_max(ig)=s2(ig)*log(n2(ig))
901 !! if (n2(ig) .lt. 1) s_max(ig)=0.
902  s_max(ig)=s2(ig)*log(max(n2(ig),1.))
903 !!fin jyg
904  else
905  s2(ig)=0.
906  n2(ig)=0.
907  s_max(ig)=0.
908  endif
909  enddo
910 ! print *,'avant Calcul de Wmax ' !!jyg
911 
912 !-----Calcul de Wmax et ALE_BL_STAT associée
913 !!jyg le 30/04/2012
914 !! do ig=1,ngrid
915 !! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
916 !! w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14))))
917 !! ale_bl_stat(ig)=0.5*w_max(ig)**2
918 !! else
919 !! w_max(ig)=0.
920 !! ale_bl_stat(ig)=0.
921 !! endif
922 !! enddo
923  susqr2pi=su*sqrt(2.*rpi)
924  reuler=exp(1.)
925  do ig=1,ngrid
926  if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
927  w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
928  ale_bl_stat(ig)=0.5*w_max(ig)**2
929  else
930  w_max(ig)=0.
931  ale_bl_stat(ig)=0.
932  endif
933  enddo
934 
935  ENDIF ! iflag_trig_bl
936 ! print *,'ENDIF iflag_trig_bl' !!jyg
937 
938 !------------Closure------------------
939 
940  IF (iflag_clos_bl.ge.1) THEN
941 
942 !-----Calcul de ALP_BL_STAT
943  do ig=1,ngrid
944  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
945  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
946  & (w0(ig)**2)
947  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
948  & +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
949  if (iflag_clos_bl.ge.2) then
950  alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
951  & (w0(ig)**2)
952  else
953  alp_bl_conv(ig)=0.
954  endif
955  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
956  enddo
957 
958 !-----Sécurité ALP infinie
959  do ig=1,ngrid
960  if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
961  enddo
962 
963  ENDIF ! (iflag_clos_bl.ge.1)
964 
965 !!! fin nrlmd le 10/04/2012
966 
967  if (prt_level.ge.10) then
968  ig=igout
969  do l=1,nlay
970  print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
971  print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
972  enddo
973  endif
974 
975 ! print*,'avant calcul ale et alp'
976 !calcul de ALE et ALP pour la convection
977  alp_bl(:)=0.
978  ale_bl(:)=0.
979 ! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
980  do l=1,nlay
981  do ig=1,ngrid
982  alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
983  ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
984 ! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
985  enddo
986  enddo
987 
988 !test:calcul de la ponderation des couches pour KE
989 !initialisations
990 
991  fm_tot(:)=0.
992  wght_th(:,:)=1.
993  lalim_conv(:)=lalim(:)
994 
995  do k=1,klev
996  do ig=1,ngrid
997  if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
998  enddo
999  enddo
1000 
1001 ! assez bizarre car, si on est dans la couche d'alim et que alim_star et
1002 ! plus petit que 1.e-10, on prend wght_th=1.
1003  do k=1,klev
1004  do ig=1,ngrid
1005  if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1006  wght_th(ig,k)=alim_star(ig,k)
1007  endif
1008  enddo
1009  enddo
1010 
1011 ! print*,'apres wght_th'
1012 !test pour prolonger la convection
1013  do ig=1,ngrid
1014 !v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1015  if ((alim_star(ig,1).lt.1.e-10)) then
1016  lalim_conv(ig)=1
1017  wght_th(ig,1)=1.
1018 ! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1019  endif
1020  enddo
1021 
1022 !------------------------------------------------------------------------
1023 ! Modif CR/FH 20110310 : Alp integree sur la verticale.
1024 ! Integrale verticale de ALP.
1025 ! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1026 ! couches
1027 !------------------------------------------------------------------------
1028 
1029  alp_int(:)=0.
1030  dp_int(:)=0.
1031  do l=2,nlay
1032  do ig=1,ngrid
1033  if(l.LE.lmax(ig)) THEN
1034  zdp=pplay(ig,l-1)-pplay(ig,l)
1035  alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
1036  dp_int(ig)=dp_int(ig)+zdp
1037  endif
1038  enddo
1039  enddo
1040 
1041  if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1042  do ig=1,ngrid
1043 !valeur integree de alp_bl * 0.5:
1044  if (dp_int(ig)>0.) then
1045  alp_bl(ig)=alp_int(ig)/dp_int(ig)
1046  endif
1047  enddo!
1048  endif
1049 
1050 
1051 ! Facteur multiplicatif sur Alp_bl
1052  alp_bl(:)=alp_bl_k*alp_bl(:)
1053 
1054 !------------------------------------------------------------------------
1055 
1056 
1057 !calcul du ratqscdiff
1058  if (prt_level.ge.1) print*,'14e OK convect8'
1059  var=0.
1060  vardiff=0.
1061  ratqsdiff(:,:)=0.
1062 
1063  do l=1,klev
1064  do ig=1,ngrid
1065  if (l<=lalim(ig)) then
1066  var=var+alim_star(ig,l)*zqta(ig,l)*1000.
1067  endif
1068  enddo
1069  enddo
1070 
1071  if (prt_level.ge.1) print*,'14f OK convect8'
1072 
1073  do l=1,klev
1074  do ig=1,ngrid
1075  if (l<=lalim(ig)) then
1076  zf=fraca(ig,l)
1077  zf2=zf/(1.-zf)
1078  vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
1079  endif
1080  enddo
1081  enddo
1082 
1083  if (prt_level.ge.1) print*,'14g OK convect8'
1084  do l=1,nlay
1085  do ig=1,ngrid
1086  ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1087 ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1088  enddo
1089  enddo
1090 !--------------------------------------------------------------------
1091 !
1092 !ecriture des fichiers sortie
1093 ! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1094 
1095 #ifdef wrgrads_thermcell
1096  if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1097 #include "thermcell_out3d.h"
1098 #endif
1099 
1100  endif
1101 
1102  if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1103 
1104  return
1105  end
1106 
1107 !-----------------------------------------------------------------------------
1108 
1109  subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1110  IMPLICIT NONE
1111 #include "iniprint.h"
1112 
1113  integer i, k, klon,klev
1114  real pplev(klon,klev+1),pplay(klon,klev)
1115  real ztv(klon,klev)
1116  real po(klon,klev)
1117  real ztva(klon,klev)
1118  real zqla(klon,klev)
1119  real f_star(klon,klev)
1120  real zw2(klon,klev)
1121  integer long(klon)
1122  real seuil
1123  character*21 comment
1124 
1125  if (prt_level.ge.1) THEN
1126  print*,'WARNING !!! TEST ',comment
1127  endif
1128  return
1129 
1130 ! test sur la hauteur des thermiques ...
1131  do i=1,klon
1132 !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1133  if (prt_level.ge.10) then
1134  print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1135  print*,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2'
1136  do k=1,klev
1137  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
1138  enddo
1139  endif
1140  enddo
1141 
1142 
1143  return
1144  end
1145 
1146 !!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1147 ! On transporte pbl_tke pour donner therm_tke
1148 ! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1149  subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1150  & rg,pplev,therm_tke_max)
1151  implicit none
1152 
1153 #include "iniprint.h"
1154 !=======================================================================
1155 !
1156 ! Calcul du transport verticale dans la couche limite en presence
1157 ! de "thermiques" explicitement representes
1158 ! calcul du dq/dt une fois qu'on connait les ascendances
1159 !
1160 !=======================================================================
1161 
1162  integer ngrid,nlay,nsrf
1163 
1164  real ptimestep
1165  real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1166  real entr0(ngrid,nlay),rg
1167  real therm_tke_max(ngrid,nlay)
1168  real detr0(ngrid,nlay)
1169 
1170 
1171  real masse(ngrid,nlay),fm(ngrid,nlay+1)
1172  real entr(ngrid,nlay)
1173  real q(ngrid,nlay)
1174  integer lev_out ! niveau pour les print
1175 
1176  real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1177 
1178  real zzm
1179 
1180  integer ig,k
1181  integer isrf
1182 
1183 
1184  lev_out=0
1185 
1186 
1187  if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1188 
1189 ! calcul du detrainement
1190  do k=1,nlay
1191  detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1192  masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/rg
1193  enddo
1194 
1195 
1196 ! Decalage vertical des entrainements et detrainements.
1197  masse(:,1)=0.5*masse0(:,1)
1198  entr(:,1)=0.5*entr0(:,1)
1199  detr(:,1)=0.5*detr0(:,1)
1200  fm(:,1)=0.
1201  do k=1,nlay-1
1202  masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1203  entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1204  detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1205  fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1206  enddo
1207  fm(:,nlay+1)=0.
1208 
1209 !!! nrlmd le 16/09/2010
1210 ! calcul de la valeur dans les ascendances
1211 ! do ig=1,ngrid
1212 ! qa(ig,1)=q(ig,1)
1213 ! enddo
1214 !!!
1215 
1216 !do isrf=1,nsrf
1217 
1218 ! q(:,:)=therm_tke(:,:,isrf)
1219  q(:,:)=therm_tke_max(:,:)
1220 !!! nrlmd le 16/09/2010
1221  do ig=1,ngrid
1222  qa(ig,1)=q(ig,1)
1223  enddo
1224 !!!
1225 
1226  if (1==1) then
1227  do k=2,nlay
1228  do ig=1,ngrid
1229  if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1230  & 1.e-5*masse(ig,k)) then
1231  qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) &
1232  & /(fm(ig,k+1)+detr(ig,k))
1233  else
1234  qa(ig,k)=q(ig,k)
1235  endif
1236  if (qa(ig,k).lt.0.) then
1237 ! print*,'qa<0!!!'
1238  endif
1239  if (q(ig,k).lt.0.) then
1240 ! print*,'q<0!!!'
1241  endif
1242  enddo
1243  enddo
1244 
1245 ! Calcul du flux subsident
1246 
1247  do k=2,nlay
1248  do ig=1,ngrid
1249  wqd(ig,k)=fm(ig,k)*q(ig,k)
1250  if (wqd(ig,k).lt.0.) then
1251 ! print*,'wqd<0!!!'
1252  endif
1253  enddo
1254  enddo
1255  do ig=1,ngrid
1256  wqd(ig,1)=0.
1257  wqd(ig,nlay+1)=0.
1258  enddo
1259 
1260 ! Calcul des tendances
1261  do k=1,nlay
1262  do ig=1,ngrid
1263  q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) &
1264  & -wqd(ig,k)+wqd(ig,k+1)) &
1265  & *ptimestep/masse(ig,k)
1266  enddo
1267  enddo
1268 
1269  endif
1270 
1271  therm_tke_max(:,:)=q(:,:)
1272 
1273  return
1274 !!! fin nrlmd le 10/04/2012
1275  end
1276