source: lmdz_wrf/trunk/WRFV3/lmdz/thermcell_main.F90 @ 1531

Last change on this file since 1531 was 179, checked in by lfita, 10 years ago

No Ale, Alp, ... values because pphi is at the intermediate level!!!
Removing raw prints of checkig for Ale, Alp

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