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

Last change on this file since 1638 was 1638, checked in by idelkadi, 12 years ago

Introduction du declenchement stochastique de la convection

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