source: LMDZ5/trunk/libf/phylmd/thermcell_main.F90 @ 2172

Last change on this file since 2172 was 2079, checked in by lguez, 10 years ago

Protect against division by 0 in computation of proba_notrig in procedure physiq: initialize s2 at some low value instead of 0.

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