source: LMDZ5/branches/LF-private/libf/phylmd/thermcell_main.F90 @ 3839

Last change on this file since 3839 was 1790, checked in by Ehouarn Millour, 11 years ago

Ré-implémentation des commissions écrasées par la rev 1785. (Notamment lmdz1d.F)
UG
...................................................

Reimplementation of the commitions undonned by rev 1785. (Especially lmdy1d.F)
UG

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