LMDZ
thermcell_main.F90
Go to the documentation of this file.
1 !
2 ! $Id: thermcell_main.F90 2387 2015-11-07 09:27:40Z fhourdin $
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  & ,ztva )
22 
23  USE dimphy
24  USE ioipsl
25  USE indice_sol_mod
26  USE print_control_mod, ONLY: lunout,prt_level
27  IMPLICIT NONE
28 
29 !=======================================================================
30 ! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
31 ! Version du 09.02.07
32 ! Calcul du transport vertical dans la couche limite en presence
33 ! de "thermiques" explicitement representes avec processus nuageux
34 !
35 ! Reecriture a partir d'un listing papier a Habas, le 14/02/00
36 !
37 ! le thermique est suppose homogene et dissipe par melange avec
38 ! son environnement. la longueur l_mix controle l'efficacite du
39 ! melange
40 !
41 ! Le calcul du transport des differentes especes se fait en prenant
42 ! en compte:
43 ! 1. un flux de masse montant
44 ! 2. un flux de masse descendant
45 ! 3. un entrainement
46 ! 4. un detrainement
47 !
48 ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
49 ! Introduction of an implicit computation of vertical advection in
50 ! the environment of thermal plumes in thermcell_dq
51 ! impl = 0 : explicit, 1 : implicit, -1 : old version
52 ! controled by iflag_thermals =
53 ! 15, 16 run with impl=-1 : numerical convergence with NPv3
54 ! 17, 18 run with impl=1 : more stable
55 ! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
56 !
57 !=======================================================================
58 
59 
60 !-----------------------------------------------------------------------
61 ! declarations:
62 ! -------------
63 
64 #include "YOMCST.h"
65 #include "YOETHF.h"
66 #include "FCTTRE.h"
67 #include "thermcell.h"
68 
69 ! arguments:
70 ! ----------
71 
72 !IM 140508
73  INTEGER itap
74 
75  INTEGER ngrid,nlay
76  real ptimestep
77  REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
78  REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
79  REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
80  REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
81  REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
82  real pphi(ngrid,nlay)
83  LOGICAL debut
84 
85 ! local:
86 ! ------
87 
88  integer icount
89 
90  integer, save :: dvdq=1,dqimpl=-1
91 !$OMP THREADPRIVATE(dvdq,dqimpl)
92  data icount/0/
93  save icount
94 !$OMP THREADPRIVATE(icount)
95 
96  integer,save :: igout=1
97 !$OMP THREADPRIVATE(igout)
98  integer,save :: lunout1=6
99 !$OMP THREADPRIVATE(lunout1)
100  integer,save :: lev_out=10
101 !$OMP THREADPRIVATE(lev_out)
102 
103  REAL susqr2pi, reuler
104 
105  INTEGER ig,k,l,ll,ierr
106  real zsortie1d(klon)
107  INTEGER lmax(klon),lmin(klon),lalim(klon)
108  INTEGER lmix(klon)
109  INTEGER lmix_bis(klon)
110  real linter(klon)
111  real zmix(klon)
112  real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
113 ! real fraca(klon,klev)
114 
115  real zmax_sec(klon)
116 !on garde le zmax du pas de temps precedent
117  real zmax0(klon)
118 !FH/IM save zmax0
119 
120  real lambda
121 
122  real zlev(klon,klev+1),zlay(klon,klev)
123  real deltaz(klon,klev)
124  REAL zh(klon,klev)
125  real zthl(klon,klev),zdthladj(klon,klev)
126  REAL ztv(klon,klev)
127  real zu(klon,klev),zv(klon,klev),zo(klon,klev)
128  real zl(klon,klev)
129  real zsortie(klon,klev)
130  real zva(klon,klev)
131  real zua(klon,klev)
132  real zoa(klon,klev)
133 
134  real zta(klon,klev)
135  real zha(klon,klev)
136  real fraca(klon,klev+1)
137  real zf,zf2
138  real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
139  real q2(klon,klev)
140 ! FH probleme de dimensionnement avec l'allocation dynamique
141 ! common/comtherm/thetath2,wth2
142  real wq(klon,klev)
143  real wthl(klon,klev)
144  real wthv(klon,klev)
145 
146  real ratqscth(klon,klev)
147  real var
148  real vardiff
149  real ratqsdiff(klon,klev)
150 
151  logical sorties
152  real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
153  real zpspsk(klon,klev)
154 
155  real wmax(klon)
156  real wmax_tmp(klon)
157  real wmax_sec(klon)
158  real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
159  real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
160 
161  real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
162 !niveau de condensation
163  integer nivcon(klon)
164  real zcon(klon)
165  REAL chi
166  real zcon2(klon)
167  real pcon(klon)
168  real zqsat(klon,klev)
169  real zqsatth(klon,klev)
170 
171  real f_star(klon,klev+1),entr_star(klon,klev)
172  real detr_star(klon,klev)
173  real alim_star_tot(klon)
174  real alim_star(klon,klev)
175  real alim_star_clos(klon,klev)
176  real f(klon), f0(klon)
177 !FH/IM save f0
178  real zlevinter(klon)
179  real seuil
180  real csc(klon,klev)
181 
182 !!! nrlmd le 10/04/2012
183 
184 !------Entrées
185  real pbl_tke(klon,klev+1,nbsrf)
186  real pctsrf(klon,nbsrf)
187  real omega(klon,klev)
188  real airephy(klon)
189 !------Sorties
190  real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
191  real therm_tke_max0(klon),env_tke_max0(klon)
192  real n2(klon),s2(klon)
193  real ale_bl_stat(klon)
194  real therm_tke_max(klon,klev),env_tke_max(klon,klev)
195  real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
196 !------Local
197  integer nsrf
198  real rhobarz0(klon) ! Densité au LCL
199  logical ok_lcl(klon) ! Existence du LCL des thermiques
200  integer klcl(klon) ! Niveau du LCL
201  real interp(klon) ! Coef d'interpolation pour le LCL
202 !--Triggering
203  real su ! Surface unité: celle d'un updraft élémentaire
204  parameter(su=4e4)
205  real hcoef ! Coefficient directeur pour le calcul de s2
206  parameter(hcoef=1)
207  real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
208  parameter(hmincoef=0.3)
209  real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
210  parameter(eps1=0.3)
211  real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2
212  real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
213  real zmax_moy_coef
214  parameter(zmax_moy_coef=0.33)
215  real depth(klon) ! Epaisseur moyenne du cumulus
216  real w_max(klon) ! Vitesse max statistique
217  real s_max(klon)
218 !--Closure
219  real pbl_tke_max(klon,klev) ! Profil de TKE moyenne
220  real pbl_tke_max0(klon) ! TKE moyenne au LCL
221  real w_ls(klon,klev) ! Vitesse verticale grande échelle (m/s)
222  real coef_m ! On considère un rendement pour alp_bl_fluct_m
223  parameter(coef_m=1.)
224  real coef_tke ! On considère un rendement pour alp_bl_fluct_tke
225  parameter(coef_tke=1.)
226 
227 !!! fin nrlmd le 10/04/2012
228 
229 !
230  !nouvelles variables pour la convection
231  real ale_bl(klon)
232  real alp_bl(klon)
233  real alp_int(klon),dp_int(klon),zdp
234  real ale_int(klon)
235  integer n_int(klon)
236  real fm_tot(klon)
237  real wght_th(klon,klev)
238  integer lalim_conv(klon)
239 !v1d logical therm
240 !v1d save therm
241 
242  character*2 str2
243  character*10 str10
244 
245  character (len=20) :: modname='thermcell_main'
246  character (len=80) :: abort_message
247 
248  EXTERNAL scopy
249 !
250 
251 !-----------------------------------------------------------------------
252 ! initialisation:
253 ! ---------------
254 !
255 
256  seuil=0.25
257 
258  if (debut) then
259  if (iflag_thermals==15.or.iflag_thermals==16) then
260  dvdq=0
261  dqimpl=-1
262  else
263  dvdq=1
264  dqimpl=1
265  endif
266 
267  fm0=0.
268  entr0=0.
269  detr0=0.
270  endif
271  fm=0. ; entr=0. ; detr=0.
272  icount=icount+1
273 
274 !IM 090508 beg
275 !print*,'====================================================================='
276 !print*,'====================================================================='
277 !print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
278 !print*,'====================================================================='
279 !print*,'====================================================================='
280 !IM 090508 end
281 
282  if (prt_level.ge.1) print*,'thermcell_main V4'
283 
284  sorties=.true.
285  IF(ngrid.NE.klon) THEN
286  print*
287  print*,'STOP dans convadj'
288  print*,'ngrid =',ngrid
289  print*,'klon =',klon
290  ENDIF
291 !
292 ! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
293  do ig=1,klon
294  f0(ig)=max(f0(ig),1.e-2)
295  zmax0(ig)=max(zmax0(ig),40.)
296 !IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
297  enddo
298 
299  if (prt_level.ge.20) then
300  do ig=1,ngrid
301  print*,'th_main ig f0',ig,f0(ig)
302  enddo
303  endif
304 !-----------------------------------------------------------------------
305 ! Calcul de T,q,ql a partir de Tl et qT dans l environnement
306 ! --------------------------------------------------------------------
307 !
308  CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
309  & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
310 
311  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
312 
313 !------------------------------------------------------------------------
314 ! --------------------
315 !
316 !
317 ! + + + + + + + + + + +
318 !
319 !
320 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
321 ! wh,wt,wo ...
322 !
323 ! + + + + + + + + + + + zh,zu,zv,zo,rho
324 !
325 !
326 ! -------------------- zlev(1)
327 ! \\\\\\\\\\\\\\\\\\\\
328 !
329 !
330 
331 !-----------------------------------------------------------------------
332 ! Calcul des altitudes des couches
333 !-----------------------------------------------------------------------
334 
335  do l=2,nlay
336  zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/rg
337  enddo
338  zlev(:,1)=0.
339  zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/rg
340  do l=1,nlay
341  zlay(:,l)=pphi(:,l)/rg
342  enddo
343 !calcul de l epaisseur des couches
344  do l=1,nlay
345  deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
346  enddo
347 
348 ! print*,'2 OK convect8'
349 !-----------------------------------------------------------------------
350 ! Calcul des densites
351 !-----------------------------------------------------------------------
352 
353  rho(:,:)=pplay(:,:)/(zpspsk(:,:)*rd*ztv(:,:))
354 
355  if (prt_level.ge.10)write(lunout,*) &
356  & 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
357  rhobarz(:,1)=rho(:,1)
358 
359  do l=2,nlay
360  rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
361  enddo
362 
363 !calcul de la masse
364  do l=1,nlay
365  masse(:,l)=(pplev(:,l)-pplev(:,l+1))/rg
366  enddo
367 
368  if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
369 
370 !------------------------------------------------------------------
371 !
372 ! /|\
373 ! -------- | F_k+1 -------
374 ! ----> D_k
375 ! /|\ <---- E_k , A_k
376 ! -------- | F_k ---------
377 ! ----> D_k-1
378 ! <---- E_k-1 , A_k-1
379 !
380 !
381 !
382 !
383 !
384 ! ---------------------------
385 !
386 ! ----- F_lmax+1=0 ---------- \
387 ! lmax (zmax) |
388 ! --------------------------- |
389 ! |
390 ! --------------------------- |
391 ! |
392 ! --------------------------- |
393 ! |
394 ! --------------------------- |
395 ! |
396 ! --------------------------- |
397 ! | E
398 ! --------------------------- | D
399 ! |
400 ! --------------------------- |
401 ! |
402 ! --------------------------- \ |
403 ! lalim | |
404 ! --------------------------- | |
405 ! | |
406 ! --------------------------- | |
407 ! | A |
408 ! --------------------------- | |
409 ! | |
410 ! --------------------------- | |
411 ! lmin (=1 pour le moment) | |
412 ! ----- F_lmin=0 ------------ / /
413 !
414 ! ---------------------------
415 ! //////////////////////////
416 !
417 !
418 !=============================================================================
419 ! Calculs initiaux ne faisant pas intervenir les changements de phase
420 !=============================================================================
421 
422 !------------------------------------------------------------------
423 ! 1. alim_star est le profil vertical de l'alimentation a la base du
424 ! panache thermique, calcule a partir de la flotabilite de l'air sec
425 ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
426 !------------------------------------------------------------------
427 !
428  entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
429  lmin=1
430 
431 !-----------------------------------------------------------------------------
432 ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
433 ! panache sec conservatif (e=d=0) alimente selon alim_star
434 ! Il s'agit d'un calcul de type CAPE
435 ! zmax_sec est utilise pour determiner la geometrie du thermique.
436 !------------------------------------------------------------------------------
437 !---------------------------------------------------------------------------------
438 !calcul du melange et des variables dans le thermique
439 !--------------------------------------------------------------------------------
440 !
441  if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
442 !IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
443 
444 ! Gestion temporaire de plusieurs appels à thermcell_plume au travers
445 ! de la variable iflag_thermals
446 
447 ! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
448  if (iflag_thermals_ed<=9) then
449 ! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
450  CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
451  & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
452  & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
453  & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
454  & ,lev_out,lunout1,igout)
455 
456  elseif (iflag_thermals_ed>9) then
457 ! print*,'THERM RIO et al 2010, version d Arnaud'
458  CALL thermcellv1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
459  & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
460  & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
461  & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
462  & ,lev_out,lunout1,igout)
463 
464  endif
465 
466  if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
467 
468  call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
469  call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
470 
471  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
472  if (prt_level.ge.10) then
473  write(lunout1,*) 'Dans thermcell_main 2'
474  write(lunout1,*) 'lmin ',lmin(igout)
475  write(lunout1,*) 'lalim ',lalim(igout)
476  write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
477  write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
478  & ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
479  endif
480 
481 !-------------------------------------------------------------------------------
482 ! Calcul des caracteristiques du thermique:zmax,zmix,wmax
483 !-------------------------------------------------------------------------------
484 !
485  CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
486  & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
487 ! Attention, w2 est transforme en sa racine carree dans cette routine
488 ! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
489  wmax_tmp=0.
490  do l=1,nlay
491  wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
492  enddo
493 ! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
494 
495 
496 
497  call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
498  call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
499  call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
500  call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
501 
502  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
503 
504 !-------------------------------------------------------------------------------
505 ! Fermeture,determination de f
506 !-------------------------------------------------------------------------------
507 !
508 !
509 !! write(lunout,*)'THERM NOUVEAU XXXXX'
510  CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
511  & lalim,lmin,zmax_sec,wmax_sec,lev_out)
512 
513 
514 call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
515 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
516 
517  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
518  if (prt_level.ge.10) then
519  write(lunout1,*) 'Dans thermcell_main 1b'
520  write(lunout1,*) 'lmin ',lmin(igout)
521  write(lunout1,*) 'lalim ',lalim(igout)
522  write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
523  write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
524  & ,l=1,lalim(igout)+4)
525  endif
526 
527 
528 
529 
530 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
531 ! Apparemment sans importance
532  alim_star_clos(:,:)=alim_star(:,:)
533  alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
534 !
535 !CR Appel de la fermeture seche
536  if (iflag_thermals_closure.eq.1) then
537 
538  CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
539  & zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
540 
541 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
542 ! Appel avec les zmax et wmax tenant compte de la condensation
543 ! Semble moins bien marcher
544  else if (iflag_thermals_closure.eq.2) then
545 
546  CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
547  & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
548 
549  endif
550 
551 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
552 
553  if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
554 
555  if (tau_thermals>1.) then
556  lambda=exp(-ptimestep/tau_thermals)
557  f0=(1.-lambda)*f+lambda*f0
558  else
559  f0=f
560  endif
561 
562 ! Test valable seulement en 1D mais pas genant
563  if (.not. (f0(1).ge.0.) ) then
564  abort_message = .not..ge.' (f0(1)0.)'
565  CALL abort_physic (modname,abort_message,1)
566  endif
567 
568 !-------------------------------------------------------------------------------
569 !deduction des flux
570 !-------------------------------------------------------------------------------
571 
572  CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
573  & lalim,lmax,alim_star, &
574  & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
575  & detr,zqla,lev_out,lunout1,igout)
576 !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
577 
578  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
579  call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
580  call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
581 
582 !------------------------------------------------------------------
583 ! On ne prend pas directement les profils issus des calculs precedents
584 ! mais on s'autorise genereusement une relaxation vers ceci avec
585 ! une constante de temps tau_thermals (typiquement 1800s).
586 !------------------------------------------------------------------
587 
588  if (tau_thermals>1.) then
589  lambda=exp(-ptimestep/tau_thermals)
590  fm0=(1.-lambda)*fm+lambda*fm0
591  entr0=(1.-lambda)*entr+lambda*entr0
592  detr0=(1.-lambda)*detr+lambda*detr0
593  else
594  fm0=fm
595  entr0=entr
596  detr0=detr
597  endif
598 
599 !c------------------------------------------------------------------
600 ! calcul du transport vertical
601 !------------------------------------------------------------------
602 
603  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
604  & zthl,zdthladj,zta,lev_out)
605  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
606  & po,pdoadj,zoa,lev_out)
607 
608 !------------------------------------------------------------------
609 ! Calcul de la fraction de l'ascendance
610 !------------------------------------------------------------------
611  do ig=1,klon
612  fraca(ig,1)=0.
613  fraca(ig,nlay+1)=0.
614  enddo
615  do l=2,nlay
616  do ig=1,klon
617  if (zw2(ig,l).gt.1.e-10) then
618  fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
619  else
620  fraca(ig,l)=0.
621  endif
622  enddo
623  enddo
624 
625 !------------------------------------------------------------------
626 ! calcul du transport vertical du moment horizontal
627 !------------------------------------------------------------------
628 
629 !IM 090508
630  if (dvdq == 0 ) then
631 
632 ! Calcul du transport de V tenant compte d'echange par gradient
633 ! de pression horizontal avec l'environnement
634 
635  call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
636 ! & ,fraca*dvdq,zmax &
637  & ,fraca,zmax &
638  & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
639 
640  else
641 
642 ! calcul purement conservatif pour le transport de V
643  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
644  & ,zu,pduadj,zua,lev_out)
645  call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
646  & ,zv,pdvadj,zva,lev_out)
647 
648  endif
649 
650 ! print*,'13 OK convect8'
651  do l=1,nlay
652  do ig=1,ngrid
653  pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
654  enddo
655  enddo
656 
657  if (prt_level.ge.1) print*,'14 OK convect8'
658 !------------------------------------------------------------------
659 ! Calculs de diagnostiques pour les sorties
660 !------------------------------------------------------------------
661 !calcul de fraca pour les sorties
662 
663  if (sorties) then
664  if (prt_level.ge.1) print*,'14a OK convect8'
665 ! calcul du niveau de condensation
666 ! initialisation
667  do ig=1,ngrid
668  nivcon(ig)=0
669  zcon(ig)=0.
670  enddo
671 !nouveau calcul
672  do ig=1,ngrid
673  chi=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
674  pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**chi
675  enddo
676 !IM do k=1,nlay
677  do k=1,nlay-1
678  do ig=1,ngrid
679  if ((pcon(ig).le.pplay(ig,k)) &
680  & .and.(pcon(ig).gt.pplay(ig,k+1))) then
681  zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(rg*rho(ig,k))/100.
682  endif
683  enddo
684  enddo
685 !IM
686  ierr=0
687  do ig=1,ngrid
688  if (pcon(ig).le.pplay(ig,nlay)) then
689  zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(rg*rho(ig,nlay))/100.
690  ierr=1
691  endif
692  enddo
693  if (ierr==1) then
694  abort_message = 'thermcellV0_main: les thermiques vont trop haut '
695  CALL abort_physic (modname,abort_message,1)
696  endif
697 
698  if (prt_level.ge.1) print*,'14b OK convect8'
699  do k=nlay,1,-1
700  do ig=1,ngrid
701  if (zqla(ig,k).gt.1e-10) then
702  nivcon(ig)=k
703  zcon(ig)=zlev(ig,k)
704  endif
705  enddo
706  enddo
707  if (prt_level.ge.1) print*,'14c OK convect8'
708 !calcul des moments
709 !initialisation
710  do l=1,nlay
711  do ig=1,ngrid
712  q2(ig,l)=0.
713  wth2(ig,l)=0.
714  wth3(ig,l)=0.
715  ratqscth(ig,l)=0.
716  ratqsdiff(ig,l)=0.
717  enddo
718  enddo
719  if (prt_level.ge.1) print*,'14d OK convect8'
720  if (prt_level.ge.10)write(lunout,*) &
721  & 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
722  do l=1,nlay
723  do ig=1,ngrid
724  zf=fraca(ig,l)
725  zf2=zf/(1.-zf)
726 !
727  thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
728  if(zw2(ig,l).gt.1.e-10) then
729  wth2(ig,l)=zf2*(zw2(ig,l))**2
730  else
731  wth2(ig,l)=0.
732  endif
733  wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) &
734  & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
735  q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
736 !test: on calcul q2/po=ratqsc
737  ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
738  enddo
739  enddo
740 !calcul des flux: q, thetal et thetav
741  do l=1,nlay
742  do ig=1,ngrid
743  wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
744  wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
745  wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
746  enddo
747  enddo
748 !
749 ! $Id: thermcell_main.F90 2387 2015-11-07 09:27:40Z fhourdin $
750 !
751  CALL thermcell_alp(ngrid,nlay,ptimestep &
752  & ,pplay,pplev &
753  & ,fm0,entr0,lmax &
754  & ,ale_bl,alp_bl,lalim_conv,wght_th &
755  & ,zw2,fraca &
756 !!! necessire en plus
757  & ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
758 !!! nrlmd le 10/04/2012
759  & ,pbl_tke,pctsrf,omega,airephy &
760  & ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
761  & ,n2,s2,ale_bl_stat &
762  & ,therm_tke_max,env_tke_max &
763  & ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
764  & ,alp_bl_conv,alp_bl_stat &
765 !!! fin nrlmd le 10/04/2012
766  & )
767 
768 
769 
770 !calcul du ratqscdiff
771  if (prt_level.ge.1) print*,'14e OK convect8'
772  var=0.
773  vardiff=0.
774  ratqsdiff(:,:)=0.
775 
776  do l=1,klev
777  do ig=1,ngrid
778  if (l<=lalim(ig)) then
779  var=var+alim_star(ig,l)*zqta(ig,l)*1000.
780  endif
781  enddo
782  enddo
783 
784  if (prt_level.ge.1) print*,'14f OK convect8'
785 
786  do l=1,klev
787  do ig=1,ngrid
788  if (l<=lalim(ig)) then
789  zf=fraca(ig,l)
790  zf2=zf/(1.-zf)
791  vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
792  endif
793  enddo
794  enddo
795 
796  if (prt_level.ge.1) print*,'14g OK convect8'
797  do l=1,nlay
798  do ig=1,ngrid
799  ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
800 ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
801  enddo
802  enddo
803 !--------------------------------------------------------------------
804 !
805 !ecriture des fichiers sortie
806 ! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
807 
808  endif
809 
810  if (prt_level.ge.1) print*,'thermcell_main FIN OK'
811 
812  return
813  end
814 
815 !-----------------------------------------------------------------------------
816 
817  subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
818  USE print_control_mod, ONLY: prt_level
819  IMPLICIT NONE
820 
821  integer i, k, klon,klev
822  real pplev(klon,klev+1),pplay(klon,klev)
823  real ztv(klon,klev)
824  real po(klon,klev)
825  real ztva(klon,klev)
826  real zqla(klon,klev)
827  real f_star(klon,klev)
828  real zw2(klon,klev)
829  integer long(klon)
830  real seuil
831  character*21 comment
832 
833  if (prt_level.ge.1) THEN
834  print*,'WARNING !!! TEST ',comment
835  endif
836  return
837 
838 ! test sur la hauteur des thermiques ...
839  do i=1,klon
840 !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
841  if (prt_level.ge.10) then
842  print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
843  print*,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2'
844  do k=1,klev
845  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)
846  enddo
847  endif
848  enddo
849 
850 
851  return
852  end
853 
854 !!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
855 ! On transporte pbl_tke pour donner therm_tke
856 ! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
857  subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
858  & rg,pplev,therm_tke_max)
859  USE print_control_mod, ONLY: prt_level
860  implicit none
861 
862 !=======================================================================
863 !
864 ! Calcul du transport verticale dans la couche limite en presence
865 ! de "thermiques" explicitement representes
866 ! calcul du dq/dt une fois qu'on connait les ascendances
867 !
868 !=======================================================================
869 
870  integer ngrid,nlay,nsrf
871 
872  real ptimestep
873  real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
874  real entr0(ngrid,nlay),rg
875  real therm_tke_max(ngrid,nlay)
876  real detr0(ngrid,nlay)
877 
878 
879  real masse(ngrid,nlay),fm(ngrid,nlay+1)
880  real entr(ngrid,nlay)
881  real q(ngrid,nlay)
882  integer lev_out ! niveau pour les print
883 
884  real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
885 
886  real zzm
887 
888  integer ig,k
889  integer isrf
890 
891 
892  lev_out=0
893 
894 
895  if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
896 
897 ! calcul du detrainement
898  do k=1,nlay
899  detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
900  masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/rg
901  enddo
902 
903 
904 ! Decalage vertical des entrainements et detrainements.
905  masse(:,1)=0.5*masse0(:,1)
906  entr(:,1)=0.5*entr0(:,1)
907  detr(:,1)=0.5*detr0(:,1)
908  fm(:,1)=0.
909  do k=1,nlay-1
910  masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
911  entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
912  detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
913  fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
914  enddo
915  fm(:,nlay+1)=0.
916 
917 !!! nrlmd le 16/09/2010
918 ! calcul de la valeur dans les ascendances
919 ! do ig=1,ngrid
920 ! qa(ig,1)=q(ig,1)
921 ! enddo
922 !!!
923 
924 !do isrf=1,nsrf
925 
926 ! q(:,:)=therm_tke(:,:,isrf)
927  q(:,:)=therm_tke_max(:,:)
928 !!! nrlmd le 16/09/2010
929  do ig=1,ngrid
930  qa(ig,1)=q(ig,1)
931  enddo
932 !!!
933 
934  if (1==1) then
935  do k=2,nlay
936  do ig=1,ngrid
937  if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
938  & 1.e-5*masse(ig,k)) then
939  qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) &
940  & /(fm(ig,k+1)+detr(ig,k))
941  else
942  qa(ig,k)=q(ig,k)
943  endif
944  if (qa(ig,k).lt.0.) then
945 ! print*,'qa<0!!!'
946  endif
947  if (q(ig,k).lt.0.) then
948 ! print*,'q<0!!!'
949  endif
950  enddo
951  enddo
952 
953 ! Calcul du flux subsident
954 
955  do k=2,nlay
956  do ig=1,ngrid
957  wqd(ig,k)=fm(ig,k)*q(ig,k)
958  if (wqd(ig,k).lt.0.) then
959 ! print*,'wqd<0!!!'
960  endif
961  enddo
962  enddo
963  do ig=1,ngrid
964  wqd(ig,1)=0.
965  wqd(ig,nlay+1)=0.
966  enddo
967 
968 ! Calcul des tendances
969  do k=1,nlay
970  do ig=1,ngrid
971  q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) &
972  & -wqd(ig,k)+wqd(ig,k+1)) &
973  & *ptimestep/masse(ig,k)
974  enddo
975  enddo
976 
977  endif
978 
979  therm_tke_max(:,:)=q(:,:)
980 
981  return
982 !!! fin nrlmd le 10/04/2012
983  end
984 
thermcell_env
subroutine thermcell_env(ngrid, nlay, po, pt, pu, pv, pplay, pplev, zo, zh, zl, ztv, zthl, zu, zv, zpspsk, pqsat, lev_out)
Definition: thermcell_env.F90:3
thermcell_height
subroutine thermcell_height(ngrid, nlay, lalim, lmin, linter, lmix, zw2, zlev, lmax, zmax, zmax0, zmix, wmax, lev_out)
Definition: thermcell_height.F90:3
itap
!$Id klon initialisation mois suivants day_rain itap
Definition: calcul_divers.h:18
l
c c $Id c c calculs statistiques distribution nuage ftion du regime dynamique c c Ce calcul doit etre fait a partir de valeurs mensuelles CALL nbregdyn DO kmaxm1 DO l
Definition: calcul_REGDYN.h:13
dimphy::klon
integer, save klon
Definition: dimphy.F90:3
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
airephy
!$Id!c c c Common de passage de la geometrie de la dynamique a la physique real airephy(klon)
scopy
subroutine scopy(n, sx, incx, sy, incy)
Definition: cray.F:9
thermcell_dry
subroutine thermcell_dry(ngrid, nlay, zlev, pphi, ztv, alim_star, lalim, lmin, zmax, wmax, lev_out)
Definition: thermcell_dry.F90:6
var
!$Header!integer nvarmx s s s var
Definition: gradsdef.h:20
thermcell_closure
subroutine thermcell_closure(ngrid, nlay, r_aspect, ptimestep, rho, zlev, lalim, alim_star, f_star, zmax, wmax, f, lev_out)
Definition: thermcell_closure.F90:6
pplay
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL pplay
Definition: calcul_STDlev.h:26
prt_level
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
Definition: write_paramLMDZ_dyn.h:9
parameter
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
nsrf
c c $Id c nsrf
Definition: ini_coord_REGDYN.h:4
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
ll
!$Id Turb_fcg!implicit none!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc!cette routine permet d obtenir hq et ainsi de!pouvoir calculer la convergence et le cisaillement dans la physiq!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER klev REAL ll
Definition: 1Dconv.h:27
thermcell_flux2
subroutine thermcell_flux2(ngrid, klev, ptimestep, masse, lalim, lmax, alim_star, entr_star, detr_star, f, rhobarz, zlev, zw2, fm, entr, detr, zqla, lev_out, lunout1, igout)
Definition: thermcell_flux2.F90:8
iflag_thermals
nsplit_thermals!nrlmd le iflag_clos_bl tau_trig_deep real::s_trig!fin nrlmd le fact_thermals_ed_dz iflag_wake iflag_thermals_closure common ctherm1 iflag_thermals
Definition: thermcell.h:12
indice_sol_mod::nbsrf
integer, parameter nbsrf
Definition: indice_sol_mod.F90:2
thermcell_dq
subroutine thermcell_dq(ngrid, nlay, impl, ptimestep, fm, entr, masse, q, dq, qa, lev_out)
Definition: thermcell_dq.F90:3
print_control_mod
Definition: print_control_mod.F90:2
thermcell_alp
subroutine thermcell_alp(ngrid, nlay, ptimestep, pplay, pplev, fm0, entr0, lmax, ale_bl, alp_bl, lalim_conv, wght_th, zw2, fraca, pcon, rhobarz, wth3, wmax_sec, lalim, fm, alim_star, zmax, pbl_tke, pctsrf, omega, airephy, zlcl, fraca0, w0, w_conv, therm_tke_max0, env_tke_max0, n2, s2, ale_bl_stat, therm_tke_max, env_tke_max, alp_bl_det, alp_bl_fluct_m, alp_bl_fluct_tke, alp_bl_conv, alp_bl_stat)
Definition: thermcell_alp.F90:19
tau_thermals
nsplit_thermals!nrlmd le iflag_clos_bl tau_trig_deep real::s_trig!fin nrlmd le fact_thermals_ed_dz iflag_wake iflag_thermals_closure common ctherm1 iflag_thermals_closure common ctherm2 tau_thermals
Definition: thermcell.h:12
iflag_thermals_ed
nsplit_thermals!nrlmd le iflag_clos_bl tau_trig_deep real::s_trig!fin nrlmd le fact_thermals_ed_dz iflag_wake iflag_thermals_closure common ctherm1 iflag_thermals_closure common ctherm2 fact_thermals_ed_dz common ctherm4 iflag_wake common ctherm5 iflag_thermals_ed
Definition: thermcell.h:12
abort_physic
subroutine abort_physic(modname, message, ierr)
Definition: abort_physic.F90:3
pt
INTERFACE SUBROUTINE RRTM_ECRT_140GP pt
Definition: rrtm_ecrt_140gp.intfb.h:3
thermcell_dv2
subroutine thermcell_dv2(ngrid, nlay, ptimestep, fm, entr, masse, fraca, larga, u, v, du, dv, ua, va, lev_out)
Definition: thermcell_dv2.F90:4
dimphy
Definition: dimphy.F90:1
indice_sol_mod
Definition: indice_sol_mod.F90:1
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7
rg
real rg
Definition: comcstphy.h:1
lmax
lmax
Definition: iotd.h:1