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

Last change on this file since 2346 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM

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