source: lmdz_wrf/trunk/WRFV3/lmdz/thermcellV0_main.F90 @ 1465

Last change on this file since 1465 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 68.4 KB
Line 
1!
2! $Id$
3!
4      SUBROUTINE thermcellV0_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     &                  ,r_aspect,l_mix,tau_thermals &
11     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
12     &                  ,zmax0, f0,zw2,fraca)
13
14      USE dimphy
15      USE comgeomphy , ONLY:rlond,rlatd
16      IMPLICIT NONE
17
18!=======================================================================
19!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
20!   Version du 09.02.07
21!   Calcul du transport vertical dans la couche limite en presence
22!   de "thermiques" explicitement representes avec processus nuageux
23!
24!   Réécriture à partir d'un listing papier à Habas, le 14/02/00
25!
26!   le thermique est supposé homogène et dissipé par mélange avec
27!   son environnement. la longueur l_mix contrôle l'efficacité du
28!   mélange
29!
30!   Le calcul du transport des différentes espèces se fait en prenant
31!   en compte:
32!     1. un flux de masse montant
33!     2. un flux de masse descendant
34!     3. un entrainement
35!     4. un detrainement
36!
37!=======================================================================
38
39!-----------------------------------------------------------------------
40!   declarations:
41!   -------------
42
43#include "dimensions.h"
44#include "YOMCST.h"
45#include "YOETHF.h"
46#include "FCTTRE.h"
47#include "iniprint.h"
48
49!   arguments:
50!   ----------
51
52!IM 140508
53      INTEGER itap
54
55      INTEGER ngrid,nlay,w2di
56      real tau_thermals
57      real ptimestep,l_mix,r_aspect
58      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
59      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
60      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
61      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
62      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
63      real pphi(ngrid,nlay)
64
65!   local:
66!   ------
67
68      integer icount
69      data icount/0/
70      save icount
71!$OMP THREADPRIVATE(icount)
72
73      integer,save :: igout=1
74!$OMP THREADPRIVATE(igout)
75      integer,save :: lunout1=6
76!$OMP THREADPRIVATE(lunout1)
77      integer,save :: lev_out=10
78!$OMP THREADPRIVATE(lev_out)
79
80      INTEGER ig,k,l,ll
81      real zsortie1d(klon)
82      INTEGER lmax(klon),lmin(klon),lalim(klon)
83      INTEGER lmix(klon)
84      INTEGER lmix_bis(klon)
85      real linter(klon)
86      real zmix(klon)
87      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1)
88!      real fraca(klon,klev)
89
90      real zmax_sec(klon)
91!on garde le zmax du pas de temps precedent
92      real zmax0(klon)
93!FH/IM     save zmax0
94
95      real lambda
96
97      real zlev(klon,klev+1),zlay(klon,klev)
98      real deltaz(klon,klev)
99      REAL zh(klon,klev)
100      real zthl(klon,klev),zdthladj(klon,klev)
101      REAL ztv(klon,klev)
102      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
103      real zl(klon,klev)
104      real zsortie(klon,klev)
105      real zva(klon,klev)
106      real zua(klon,klev)
107      real zoa(klon,klev)
108
109      real zta(klon,klev)
110      real zha(klon,klev)
111      real fraca(klon,klev+1)
112      real zf,zf2
113      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
114      real q2(klon,klev)
115! FH probleme de dimensionnement avec l'allocation dynamique
116!     common/comtherm/thetath2,wth2
117   
118      real ratqscth(klon,klev)
119      real var
120      real vardiff
121      real ratqsdiff(klon,klev)
122
123      logical sorties
124      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
125      real zpspsk(klon,klev)
126
127      real wmax(klon)
128      real wmax_sec(klon)
129      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
130      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
131
132      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
133!niveau de condensation
134      integer nivcon(klon)
135      real zcon(klon)
136      REAL CHI
137      real zcon2(klon)
138      real pcon(klon)
139      real zqsat(klon,klev)
140      real zqsatth(klon,klev)
141
142      real f_star(klon,klev+1),entr_star(klon,klev)
143      real detr_star(klon,klev)
144      real alim_star_tot(klon),alim_star2(klon)
145      real alim_star(klon,klev)
146      real f(klon), f0(klon)
147!FH/IM     save f0
148      real zlevinter(klon)
149      logical debut
150       real seuil
151
152! Declaration uniquement pour les sorties dans thermcell_out3d.
153! Inutilise en 3D
154      real wthl(klon,klev)
155      real wthv(klon,klev)
156      real wq(klon,klev)
157!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158
159
160!
161      !nouvelles variables pour la convection
162      real Ale_bl(klon)
163      real Alp_bl(klon)
164      real alp_int(klon)
165      real ale_int(klon)
166      integer n_int(klon)
167      real fm_tot(klon)
168      real wght_th(klon,klev)
169      integer lalim_conv(klon)
170!v1d     logical therm
171!v1d     save therm
172
173      character*2 str2
174      character*10 str10
175
176      character (len=20) :: modname='thermcellV0_main'
177      character (len=80) :: abort_message
178
179      EXTERNAL SCOPY
180!
181
182!-----------------------------------------------------------------------
183!   initialisation:
184!   ---------------
185!
186
187      seuil=0.25
188
189      if (debut)  then
190         fm0=0.
191         entr0=0.
192         detr0=0.
193
194
195#undef wrgrads_thermcell
196#ifdef wrgrads_thermcell
197! Initialisation des sorties grads pour les thermiques.
198! Pour l'instant en 1D sur le point igout.
199! Utilise par thermcell_out3d.h
200         str10='therm'
201         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
202     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
203     &   ptimestep,str10,'therm ')
204#endif
205
206
207
208      endif
209
210      fm=0. ; entr=0. ; detr=0.
211
212      icount=icount+1
213
214!IM 090508 beg
215!print*,'====================================================================='
216!print*,'====================================================================='
217!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
218!print*,'====================================================================='
219!print*,'====================================================================='
220!IM 090508 end
221
222      if (prt_level.ge.1) print*,'thermcell_main V4'
223
224       sorties=.true.
225      IF(ngrid.NE.klon) THEN
226         PRINT*
227         PRINT*,'STOP dans convadj'
228         PRINT*,'ngrid    =',ngrid
229         PRINT*,'klon  =',klon
230      ENDIF
231!
232!Initialisation
233!
234     if (prt_level.ge.10)write(lunout,*)                                &
235    &     'WARNING thermcell_main f0=max(f0,1.e-2)'
236     do ig=1,klon
237         f0(ig)=max(f0(ig),1.e-2)
238     enddo
239
240!-----------------------------------------------------------------------
241! Calcul de T,q,ql a partir de Tl et qT dans l environnement
242!   --------------------------------------------------------------------
243!
244      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
245     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
246       
247      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
248
249!------------------------------------------------------------------------
250!                       --------------------
251!
252!
253!                       + + + + + + + + + + +
254!
255!
256!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
257!  wh,wt,wo ...
258!
259!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
260!
261!
262!                       --------------------   zlev(1)
263!                       \\\\\\\\\\\\\\\\\\\\
264!
265!
266
267!-----------------------------------------------------------------------
268!   Calcul des altitudes des couches
269!-----------------------------------------------------------------------
270
271      do l=2,nlay
272         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
273      enddo
274         zlev(:,1)=0.
275         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
276      do l=1,nlay
277         zlay(:,l)=pphi(:,l)/RG
278      enddo
279!calcul de l epaisseur des couches
280      do l=1,nlay
281         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
282      enddo
283
284!     print*,'2 OK convect8'
285!-----------------------------------------------------------------------
286!   Calcul des densites
287!-----------------------------------------------------------------------
288
289      do l=1,nlay
290         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
291      enddo
292
293!IM
294     if (prt_level.ge.10)write(lunout,*)                                &
295    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
296      rhobarz(:,1)=rho(:,1)
297
298      do l=2,nlay
299         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
300      enddo
301
302!calcul de la masse
303      do l=1,nlay
304         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
305      enddo
306
307      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
308
309!------------------------------------------------------------------
310!
311!             /|\
312!    --------  |  F_k+1 -------   
313!                              ----> D_k
314!             /|\              <---- E_k , A_k
315!    --------  |  F_k ---------
316!                              ----> D_k-1
317!                              <---- E_k-1 , A_k-1
318!
319!
320!
321!
322!
323!    ---------------------------
324!
325!    ----- F_lmax+1=0 ----------         \
326!            lmax     (zmax)              |
327!    ---------------------------          |
328!                                         |
329!    ---------------------------          |
330!                                         |
331!    ---------------------------          |
332!                                         |
333!    ---------------------------          |
334!                                         |
335!    ---------------------------          |
336!                                         |  E
337!    ---------------------------          |  D
338!                                         |
339!    ---------------------------          |
340!                                         |
341!    ---------------------------  \       |
342!            lalim                 |      |
343!    ---------------------------   |      |
344!                                  |      |
345!    ---------------------------   |      |
346!                                  | A    |
347!    ---------------------------   |      |
348!                                  |      |
349!    ---------------------------   |      |
350!    lmin  (=1 pour le moment)     |      |
351!    ----- F_lmin=0 ------------  /      /
352!
353!    ---------------------------
354!    //////////////////////////
355!
356!
357!=============================================================================
358!  Calculs initiaux ne faisant pas intervenir les changements de phase
359!=============================================================================
360
361!------------------------------------------------------------------
362!  1. alim_star est le profil vertical de l'alimentation à la base du
363!     panache thermique, calculé à partir de la flotabilité de l'air sec
364!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
365!------------------------------------------------------------------
366!
367      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
368      CALL thermcellV0_init(ngrid,nlay,ztv,zlay,zlev,  &
369     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
370
371call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
372call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
373
374
375      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
376      if (prt_level.ge.10) then
377         write(lunout1,*) 'Dans thermcell_main 1'
378         write(lunout1,*) 'lmin ',lmin(igout)
379         write(lunout1,*) 'lalim ',lalim(igout)
380         write(lunout1,*) ' ig l alim_star thetav'
381         write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
382     &   ,ztv(igout,l),l=1,lalim(igout)+4)
383      endif
384
385!v1d      do ig=1,klon
386!v1d     if (alim_star(ig,1).gt.1.e-10) then
387!v1d     therm=.true.
388!v1d     endif
389!v1d      enddo
390!-----------------------------------------------------------------------------
391!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
392!     panache sec conservatif (e=d=0) alimente selon alim_star
393!     Il s'agit d'un calcul de type CAPE
394!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
395!------------------------------------------------------------------------------
396!
397      CALL thermcellV0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
398     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
399
400call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
401call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
402
403      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
404      if (prt_level.ge.10) then
405         write(lunout1,*) 'Dans thermcell_main 1b'
406         write(lunout1,*) 'lmin ',lmin(igout)
407         write(lunout1,*) 'lalim ',lalim(igout)
408         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
409         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
410     &    ,l=1,lalim(igout)+4)
411      endif
412
413
414
415!---------------------------------------------------------------------------------
416!calcul du melange et des variables dans le thermique
417!--------------------------------------------------------------------------------
418!
419      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
420!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
421      CALL thermcellV0_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
422     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
423     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
424     &           ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
425     &            ,lev_out,lunout1,igout)
426      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
427
428      call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
429      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
430
431      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
432      if (prt_level.ge.10) then
433         write(lunout1,*) 'Dans thermcell_main 2'
434         write(lunout1,*) 'lmin ',lmin(igout)
435         write(lunout1,*) 'lalim ',lalim(igout)
436         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
437         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
438     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
439      endif
440
441!-------------------------------------------------------------------------------
442! Calcul des caracteristiques du thermique:zmax,zmix,wmax
443!-------------------------------------------------------------------------------
444!
445      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
446     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
447
448
449      call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
450      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
451      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
452      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
453
454      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
455
456!-------------------------------------------------------------------------------
457! Fermeture,determination de f
458!-------------------------------------------------------------------------------
459!
460!avant closure: on redéfinit lalim, alim_star_tot et alim_star
461!       do ig=1,klon
462!       do l=2,lalim(ig)
463!       alim_star(ig,l)=entr_star(ig,l)
464!       entr_star(ig,l)=0.
465!       enddo
466!       enddo
467
468      CALL thermcellV0_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
469     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
470
471      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
472
473      if (tau_thermals>1.) then
474         lambda=exp(-ptimestep/tau_thermals)
475         f0=(1.-lambda)*f+lambda*f0
476      else
477         f0=f
478      endif
479
480! Test valable seulement en 1D mais pas genant
481      if (.not. (f0(1).ge.0.) ) then
482        abort_message = 'Dans thermcell_main f0(1).lt.0 '
483        CALL abort_gcm (modname,abort_message,1)
484      endif
485
486!-------------------------------------------------------------------------------
487!deduction des flux
488!-------------------------------------------------------------------------------
489
490      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
491     &       lalim,lmax,alim_star,  &
492     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
493     &       detr,zqla,lev_out,lunout1,igout)
494!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
495
496      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
497      call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
498      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
499
500!------------------------------------------------------------------
501!   On ne prend pas directement les profils issus des calculs precedents
502!   mais on s'autorise genereusement une relaxation vers ceci avec
503!   une constante de temps tau_thermals (typiquement 1800s).
504!------------------------------------------------------------------
505
506      if (tau_thermals>1.) then
507         lambda=exp(-ptimestep/tau_thermals)
508         fm0=(1.-lambda)*fm+lambda*fm0
509         entr0=(1.-lambda)*entr+lambda*entr0
510!        detr0=(1.-lambda)*detr+lambda*detr0
511      else
512         fm0=fm
513         entr0=entr
514         detr0=detr
515      endif
516
517!c------------------------------------------------------------------
518!   calcul du transport vertical
519!------------------------------------------------------------------
520
521      call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse,  &
522     &                    zthl,zdthladj,zta,lev_out)
523      call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse,  &
524     &                   po,pdoadj,zoa,lev_out)
525
526!------------------------------------------------------------------
527! Calcul de la fraction de l'ascendance
528!------------------------------------------------------------------
529      do ig=1,klon
530         fraca(ig,1)=0.
531         fraca(ig,nlay+1)=0.
532      enddo
533      do l=2,nlay
534         do ig=1,klon
535            if (zw2(ig,l).gt.1.e-10) then
536            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
537            else
538            fraca(ig,l)=0.
539            endif
540         enddo
541      enddo
542     
543!------------------------------------------------------------------
544!  calcul du transport vertical du moment horizontal
545!------------------------------------------------------------------
546
547!IM 090508 
548      if (1.eq.1) then
549!IM 070508 vers. _dq       
550!     if (1.eq.0) then
551
552
553! Calcul du transport de V tenant compte d'echange par gradient
554! de pression horizontal avec l'environnement
555
556         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
557     &    ,fraca,zmax  &
558     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
559!IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
560      else
561
562! calcul purement conservatif pour le transport de V
563         call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse  &
564     &    ,zu,pduadj,zua,lev_out)
565         call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse  &
566     &    ,zv,pdvadj,zva,lev_out)
567      endif
568
569!     print*,'13 OK convect8'
570      do l=1,nlay
571         do ig=1,ngrid
572           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
573         enddo
574      enddo
575
576      if (prt_level.ge.1) print*,'14 OK convect8'
577!------------------------------------------------------------------
578!   Calculs de diagnostiques pour les sorties
579!------------------------------------------------------------------
580!calcul de fraca pour les sorties
581     
582      if (sorties) then
583      if (prt_level.ge.1) print*,'14a OK convect8'
584! calcul du niveau de condensation
585! initialisation
586      do ig=1,ngrid
587         nivcon(ig)=0
588         zcon(ig)=0.
589      enddo
590!nouveau calcul
591      do ig=1,ngrid
592      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
593      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
594      enddo
595!IM   do k=1,nlay
596      do k=1,nlay-1
597         do ig=1,ngrid
598         if ((pcon(ig).le.pplay(ig,k))  &
599     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
600            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
601         endif
602         enddo
603      enddo
604!IM
605      do ig=1,ngrid
606        if (pcon(ig).le.pplay(ig,nlay)) then
607           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
608           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
609           CALL abort_gcm (modname,abort_message,1)
610        endif
611      enddo
612      if (prt_level.ge.1) print*,'14b OK convect8'
613      do k=nlay,1,-1
614         do ig=1,ngrid
615            if (zqla(ig,k).gt.1e-10) then
616               nivcon(ig)=k
617               zcon(ig)=zlev(ig,k)
618            endif
619         enddo
620      enddo
621      if (prt_level.ge.1) print*,'14c OK convect8'
622!calcul des moments
623!initialisation
624      do l=1,nlay
625         do ig=1,ngrid
626            q2(ig,l)=0.
627            wth2(ig,l)=0.
628            wth3(ig,l)=0.
629            ratqscth(ig,l)=0.
630            ratqsdiff(ig,l)=0.
631         enddo
632      enddo     
633      if (prt_level.ge.1) print*,'14d OK convect8'
634      if (prt_level.ge.10)write(lunout,*)                                &
635    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
636      do l=1,nlay
637         do ig=1,ngrid
638            zf=fraca(ig,l)
639            zf2=zf/(1.-zf)
640!
641            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
642            if(zw2(ig,l).gt.1.e-10) then
643             wth2(ig,l)=zf2*(zw2(ig,l))**2
644            else
645             wth2(ig,l)=0.
646            endif
647!           print*,'wth2=',wth2(ig,l)
648            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
649     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
650            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
651!test: on calcul q2/po=ratqsc
652            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
653         enddo
654      enddo
655
656      if (prt_level.ge.10) then
657          print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
658          ig=igout
659          do l=1,nlay
660             print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
661          enddo
662          do l=1,nlay
663             print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
664          enddo
665      endif
666
667      do ig=1,ngrid
668      alp_int(ig)=0.
669      ale_int(ig)=0.
670      n_int(ig)=0
671      enddo
672!
673      do l=1,nlay
674      do ig=1,ngrid
675       if(l.LE.lmax(ig)) THEN
676        alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
677        ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
678        n_int(ig)=n_int(ig)+1
679       endif
680      enddo
681      enddo
682!      print*,'avant calcul ale et alp'
683!calcul de ALE et ALP pour la convection
684      do ig=1,ngrid
685!      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
686!          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
687!      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig))
688!     &           *0.1
689!valeur integree de alp_bl * 0.5:
690       if (n_int(ig).gt.0) then
691       Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
692!       if (Alp_bl(ig).lt.0.) then
693!       Alp_bl(ig)=0.
694       endif
695!       endif
696!         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
697!     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
698!            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
699!      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
700!      if (nivcon(ig).eq.1) then
701!       Ale_bl(ig)=0.
702!       else
703!valeur max de ale_bl:
704       Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2
705!     & /2.
706!     & *0.1
707!        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
708!       if (n_int(ig).gt.0) then
709!       Ale_bl(ig)=ale_int(ig)/n_int(ig)
710!        Ale_bl(ig)=4.
711!       endif
712!       endif
713!            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
714!          Ale_bl(ig)=wth2(ig,nivcon(ig))
715!          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
716      enddo
717!test:calcul de la ponderation des couches pour KE
718!initialisations
719!      print*,'ponderation'
720      do ig=1,ngrid
721           fm_tot(ig)=0.
722      enddo
723       do ig=1,ngrid
724        do k=1,klev
725           wght_th(ig,k)=1.
726        enddo
727       enddo
728       do ig=1,ngrid
729!         lalim_conv(ig)=lmix_bis(ig)
730!la hauteur de la couche alim_conv = hauteur couche alim_therm
731         lalim_conv(ig)=lalim(ig)
732!         zentr(ig)=zlev(ig,lalim(ig))
733      enddo
734      do ig=1,ngrid
735        do k=1,lalim_conv(ig)
736           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
737        enddo
738      enddo
739      do ig=1,ngrid
740        do k=1,lalim_conv(ig)
741           if (fm_tot(ig).gt.1.e-10) then
742!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
743           endif
744!on pondere chaque couche par a*
745             if (alim_star(ig,k).gt.1.e-10) then
746             wght_th(ig,k)=alim_star(ig,k)
747             else
748             wght_th(ig,k)=1.
749             endif
750        enddo
751      enddo
752!      print*,'apres wght_th'
753!test pour prolonger la convection
754      do ig=1,ngrid
755!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
756      if ((alim_star(ig,1).lt.1.e-10)) then
757      lalim_conv(ig)=1
758      wght_th(ig,1)=1.
759!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
760      endif
761      enddo
762
763!calcul du ratqscdiff
764      if (prt_level.ge.1) print*,'14e OK convect8'
765      var=0.
766      vardiff=0.
767      ratqsdiff(:,:)=0.
768      do ig=1,ngrid
769         do l=1,lalim(ig)
770            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
771         enddo
772      enddo
773      if (prt_level.ge.1) print*,'14f OK convect8'
774      do ig=1,ngrid
775          do l=1,lalim(ig)
776          zf=fraca(ig,l)
777          zf2=zf/(1.-zf)
778          vardiff=vardiff+alim_star(ig,l)  &
779     &           *(zqta(ig,l)*1000.-var)**2
780!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
781!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
782          enddo
783      enddo
784      if (prt_level.ge.1) print*,'14g OK convect8'
785      do l=1,nlay
786         do ig=1,ngrid
787            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
788!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
789         enddo
790      enddo
791!--------------------------------------------------------------------   
792!
793!ecriture des fichiers sortie
794!     print*,'15 OK convect8'
795
796      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
797#ifdef wrgrads_thermcell
798#include "thermcell_out3d.h"
799#endif
800
801      endif
802
803      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
804
805!     if(icount.eq.501) stop'au pas 301 dans thermcell_main'
806      return
807      end
808
809!-----------------------------------------------------------------------------
810
811      subroutine testV0_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
812      IMPLICIT NONE
813#include "iniprint.h"
814
815      integer i, k, klon,klev
816      real pplev(klon,klev+1),pplay(klon,klev)
817      real ztv(klon,klev)
818      real po(klon,klev)
819      real ztva(klon,klev)
820      real zqla(klon,klev)
821      real f_star(klon,klev)
822      real zw2(klon,klev)
823      integer long(klon)
824      real seuil
825      character*21 comment
826
827      if (prt_level.ge.1) THEN
828       print*,'WARNING !!! TEST ',comment
829      endif
830      return
831
832!  test sur la hauteur des thermiques ...
833         do i=1,klon
834!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
835           if (prt_level.ge.10) then
836               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
837               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
838               do k=1,klev
839                  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)
840               enddo
841           endif
842         enddo
843
844
845      return
846      end
847
848!==============================================================================
849      SUBROUTINE thermcellV0_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
850     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
851
852!-------------------------------------------------------------------------
853!thermcell_closure: fermeture, determination de f
854!-------------------------------------------------------------------------
855      IMPLICIT NONE
856
857#include "iniprint.h"
858#include "thermcell.h"
859      INTEGER ngrid,nlay
860      INTEGER ig,k       
861      REAL r_aspect,ptimestep
862      integer lev_out                           ! niveau pour les print
863
864      INTEGER lalim(ngrid)
865      REAL alim_star(ngrid,nlay)
866      REAL alim_star_tot(ngrid)
867      REAL rho(ngrid,nlay)
868      REAL zlev(ngrid,nlay)
869      REAL zmax(ngrid),zmax_sec(ngrid)
870      REAL wmax(ngrid),wmax_sec(ngrid)
871      real zdenom
872
873      REAL alim_star2(ngrid)
874
875      REAL f(ngrid)
876
877      character (len=20) :: modname='thermcellV0_main'
878      character (len=80) :: abort_message
879
880      do ig=1,ngrid
881         alim_star2(ig)=0.
882      enddo
883      do ig=1,ngrid
884         if (alim_star(ig,1).LT.1.e-10) then
885            f(ig)=0.
886         else   
887             do k=1,lalim(ig)
888                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
889     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
890             enddo
891             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
892             if (zdenom<1.e-14) then
893                print*,'ig=',ig
894                print*,'alim_star2',alim_star2(ig)
895                print*,'zmax',zmax(ig)
896                print*,'r_aspect',r_aspect
897                print*,'zdenom',zdenom
898                print*,'alim_star',alim_star(ig,:)
899                print*,'zmax_sec',zmax_sec(ig)
900                print*,'wmax_sec',wmax_sec(ig)
901                abort_message = 'zdenom<1.e-14'
902                CALL abort_gcm (modname,abort_message,1)
903             endif
904             if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then
905             f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
906     &             *alim_star2(ig))
907!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
908!    &                     zmax_sec(ig))*wmax_sec(ig))
909             if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
910             else
911             f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
912!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
913!     &                     zmax(ig))*wmax(ig))
914             if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
915             endif
916         endif
917!         f0(ig)=f(ig)
918      enddo
919      if (prt_level.ge.1) print*,'apres fermeture'
920
921!
922      return
923      end
924!==============================================================================
925      SUBROUTINE thermcellV0_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
926     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
927     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
928     &           ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
929     &           ,lev_out,lunout1,igout)
930
931!--------------------------------------------------------------------------
932!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
933!--------------------------------------------------------------------------
934
935      IMPLICIT NONE
936
937#include "YOMCST.h"
938#include "YOETHF.h"
939#include "FCTTRE.h"
940#include "iniprint.h"
941#include "thermcell.h"
942
943      INTEGER itap
944      INTEGER lunout1,igout
945      INTEGER ngrid,klev
946      REAL ptimestep
947      REAL ztv(ngrid,klev)
948      REAL zthl(ngrid,klev)
949      REAL po(ngrid,klev)
950      REAL zl(ngrid,klev)
951      REAL rhobarz(ngrid,klev)
952      REAL zlev(ngrid,klev+1)
953      REAL pplev(ngrid,klev+1)
954      REAL pphi(ngrid,klev)
955      REAL zpspsk(ngrid,klev)
956      REAL alim_star(ngrid,klev)
957      REAL zmax_sec(ngrid)
958      REAL f0(ngrid)
959      REAL l_mix
960      REAL r_aspect
961      INTEGER lalim(ngrid)
962      integer lev_out                           ! niveau pour les print
963      real zcon2(ngrid)
964   
965      real alim_star_tot(ngrid)
966
967      REAL ztva(ngrid,klev)
968      REAL ztla(ngrid,klev)
969      REAL zqla(ngrid,klev)
970      REAL zqla0(ngrid,klev)
971      REAL zqta(ngrid,klev)
972      REAL zha(ngrid,klev)
973
974      REAL detr_star(ngrid,klev)
975      REAL coefc
976      REAL detr_stara(ngrid,klev)
977      REAL detr_starb(ngrid,klev)
978      REAL detr_starc(ngrid,klev)
979      REAL detr_star0(ngrid,klev)
980      REAL detr_star1(ngrid,klev)
981      REAL detr_star2(ngrid,klev)
982
983      REAL entr_star(ngrid,klev)
984      REAL entr_star1(ngrid,klev)
985      REAL entr_star2(ngrid,klev)
986      REAL detr(ngrid,klev)
987      REAL entr(ngrid,klev)
988
989      REAL zw2(ngrid,klev+1)
990      REAL w_est(ngrid,klev+1)
991      REAL f_star(ngrid,klev+1)
992      REAL wa_moy(ngrid,klev+1)
993
994      REAL ztva_est(ngrid,klev)
995      REAL zqla_est(ngrid,klev)
996      REAL zqsatth(ngrid,klev)
997      REAL zta_est(ngrid,klev)
998
999      REAL linter(ngrid)
1000      INTEGER lmix(ngrid)
1001      INTEGER lmix_bis(ngrid)
1002      REAL    wmaxa(ngrid)
1003
1004      INTEGER ig,l,k
1005
1006      real zcor,zdelta,zcvm5,qlbef
1007      real Tbef,qsatbef
1008      real dqsat_dT,DT,num,denom
1009      REAL REPS,RLvCp,DDT0
1010      PARAMETER (DDT0=.01)
1011      logical Zsat
1012      REAL fact_gamma,fact_epsilon
1013      REAL c2(ngrid,klev)
1014
1015      Zsat=.false.
1016! Initialisation
1017      RLvCp = RLVTT/RCPD
1018     
1019      if (iflag_thermals_ed==0) then
1020         fact_gamma=1.
1021         fact_epsilon=1.
1022      else if (iflag_thermals_ed==1)  then
1023         fact_gamma=1.
1024         fact_epsilon=1.
1025      else if (iflag_thermals_ed==2)  then
1026         fact_gamma=1.
1027         fact_epsilon=2.
1028      endif
1029
1030      do l=1,klev
1031         do ig=1,ngrid
1032            zqla_est(ig,l)=0.
1033            ztva_est(ig,l)=ztva(ig,l)
1034            zqsatth(ig,l)=0.
1035         enddo
1036      enddo
1037
1038!CR: attention test couche alim
1039!     do l=2,klev
1040!     do ig=1,ngrid
1041!        alim_star(ig,l)=0.
1042!     enddo
1043!     enddo
1044!AM:initialisations du thermique
1045      do k=1,klev
1046         do ig=1,ngrid
1047            ztva(ig,k)=ztv(ig,k)
1048            ztla(ig,k)=zthl(ig,k)
1049            zqla(ig,k)=0.
1050            zqta(ig,k)=po(ig,k)
1051!
1052            ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
1053            ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
1054            zha(ig,k) = ztva(ig,k)
1055!
1056         enddo
1057      enddo
1058      do k=1,klev
1059        do ig=1,ngrid
1060           detr_star(ig,k)=0.
1061           entr_star(ig,k)=0.
1062
1063           detr_stara(ig,k)=0.
1064           detr_starb(ig,k)=0.
1065           detr_starc(ig,k)=0.
1066           detr_star0(ig,k)=0.
1067           zqla0(ig,k)=0.
1068           detr_star1(ig,k)=0.
1069           detr_star2(ig,k)=0.
1070           entr_star1(ig,k)=0.
1071           entr_star2(ig,k)=0.
1072
1073           detr(ig,k)=0.
1074           entr(ig,k)=0.
1075        enddo
1076      enddo
1077      if (prt_level.ge.1) print*,'7 OK convect8'
1078      do k=1,klev+1
1079         do ig=1,ngrid
1080            zw2(ig,k)=0.
1081            w_est(ig,k)=0.
1082            f_star(ig,k)=0.
1083            wa_moy(ig,k)=0.
1084         enddo
1085      enddo
1086
1087      if (prt_level.ge.1) print*,'8 OK convect8'
1088      do ig=1,ngrid
1089         linter(ig)=1.
1090         lmix(ig)=1
1091         lmix_bis(ig)=2
1092         wmaxa(ig)=0.
1093      enddo
1094
1095!-----------------------------------------------------------------------------------
1096!boucle de calcul de la vitesse verticale dans le thermique
1097!-----------------------------------------------------------------------------------
1098      do l=1,klev-1
1099         do ig=1,ngrid
1100
1101
1102
1103! Calcul dans la premiere couche active du thermique (ce qu'on teste
1104! en disant que la couche est instable et que w2 en bas de la couche
1105! est nulle.
1106
1107            if (ztv(ig,l).gt.ztv(ig,l+1)  &
1108     &         .and.alim_star(ig,l).gt.1.e-10  &
1109     &         .and.zw2(ig,l).lt.1e-10) then
1110
1111
1112! Le panache va prendre au debut les caracteristiques de l'air contenu
1113! dans cette couche.
1114               ztla(ig,l)=zthl(ig,l)
1115               zqta(ig,l)=po(ig,l)
1116               zqla(ig,l)=zl(ig,l)
1117               f_star(ig,l+1)=alim_star(ig,l)
1118
1119               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
1120     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
1121     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1122               w_est(ig,l+1)=zw2(ig,l+1)
1123!
1124
1125
1126            else if ((zw2(ig,l).ge.1e-10).and.  &
1127     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
1128!estimation du detrainement a partir de la geometrie du pas precedent
1129!tests sur la definition du detr
1130!calcul de detr_star et entr_star
1131
1132
1133
1134!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1135! FH le test miraculeux de Catherine ? Le bout du tunel ?
1136!               w_est(ig,3)=zw2(ig,2)*  &
1137!    &                   ((f_star(ig,2))**2)  &
1138!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
1139!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
1140!    &                   *(zlev(ig,3)-zlev(ig,2))
1141!     if (l.gt.2) then
1142!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1143
1144
1145
1146! Premier calcul de la vitesse verticale a partir de la temperature
1147! potentielle virtuelle
1148
1149! FH CESTQUOI CA ????
1150#define int1d2
1151!#undef int1d2
1152#ifdef int1d2
1153      if (l.ge.2) then
1154#else
1155      if (l.gt.2) then
1156#endif
1157
1158      if (1.eq.1) then
1159          w_est(ig,3)=zw2(ig,2)* &
1160     &      ((f_star(ig,2))**2) &
1161     &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
1162     &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
1163!     &      *1./3. &
1164     &      *(zlev(ig,3)-zlev(ig,2))
1165       endif
1166
1167
1168!---------------------------------------------------------------------------
1169!calcul de l entrainement et du detrainement lateral
1170!---------------------------------------------------------------------------
1171!
1172!test:estimation de ztva_new_est sans entrainement
1173
1174               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
1175               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1176               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
1177               qsatbef=MIN(0.5,qsatbef)
1178               zcor=1./(1.-retv*qsatbef)
1179               qsatbef=qsatbef*zcor
1180               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
1181               if (Zsat) then
1182               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
1183               DT = 0.5*RLvCp*qlbef
1184               do while (abs(DT).gt.DDT0)
1185                 Tbef=Tbef+DT
1186                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1187                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
1188                 qsatbef=MIN(0.5,qsatbef)
1189                 zcor=1./(1.-retv*qsatbef)
1190                 qsatbef=qsatbef*zcor
1191                 qlbef=zqta(ig,l-1)-qsatbef
1192
1193                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1194                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
1195                 zcor=1./(1.-retv*qsatbef)
1196                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
1197                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
1198                 denom=1.+RLvCp*dqsat_dT
1199                 DT=num/denom
1200               enddo
1201                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
1202               endif
1203        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
1204        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
1205        zta_est(ig,l)=ztva_est(ig,l)
1206        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
1207     &      -zqla_est(ig,l))-zqla_est(ig,l))
1208
1209             w_est(ig,l+1)=zw2(ig,l)*  &
1210     &                   ((f_star(ig,l))**2)  &
1211     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
1212     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
1213!     &                   *1./3. &
1214     &                   *(zlev(ig,l+1)-zlev(ig,l))
1215             if (w_est(ig,l+1).lt.0.) then
1216                w_est(ig,l+1)=zw2(ig,l)
1217             endif
1218!
1219!calcul du detrainement
1220!=======================
1221
1222!CR:on vire les modifs
1223         if (iflag_thermals_ed==0) then
1224
1225! Modifications du calcul du detrainement.
1226! Dans la version de la these de Catherine, on passe brusquement
1227! de la version seche a la version nuageuse pour le detrainement
1228! ce qui peut occasioner des oscillations.
1229! dans la nouvelle version, on commence par calculer un detrainement sec.
1230! Puis un autre en cas de nuages.
1231! Puis on combine les deux lineairement en fonction de la quantite d'eau.
1232
1233#define int1d3
1234!#undef int1d3
1235#define RIO_TH
1236#ifdef RIO_TH
1237!1. Cas non nuageux
1238! 1.1 on est sous le zmax_sec et w croit
1239          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
1240     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
1241#ifdef int1d3
1242     &       (zqla_est(ig,l).lt.1.e-10)) then
1243#else
1244     &       (zqla(ig,l-1).lt.1.e-10)) then
1245#endif
1246             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
1247     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
1248     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
1249     &       /(r_aspect*zmax_sec(ig)))
1250             detr_stara(ig,l)=detr_star(ig,l)
1251
1252       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
1253
1254! 1.2 on est sous le zmax_sec et w decroit
1255          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
1256#ifdef int1d3
1257     &            (zqla_est(ig,l).lt.1.e-10)) then
1258#else
1259     &            (zqla(ig,l-1).lt.1.e-10)) then
1260#endif
1261             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
1262     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
1263     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
1264     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
1265     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
1266     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
1267     &       *((zmax_sec(ig)-zlev(ig,l))/  &
1268     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
1269             detr_starb(ig,l)=detr_star(ig,l)
1270
1271        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
1272
1273          else
1274
1275! 1.3 dans les autres cas
1276             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
1277     &                      *(zlev(ig,l+1)-zlev(ig,l))
1278             detr_starc(ig,l)=detr_star(ig,l)
1279
1280        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
1281             
1282          endif
1283
1284#else
1285
1286! 1.1 on est sous le zmax_sec et w croit
1287          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
1288     &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
1289             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
1290     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
1291     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
1292     &       /(r_aspect*zmax_sec(ig)))
1293
1294       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
1295
1296! 1.2 on est sous le zmax_sec et w decroit
1297          else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
1298             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
1299     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
1300     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
1301     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
1302     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
1303     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
1304     &       *((zmax_sec(ig)-zlev(ig,l))/  &
1305     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
1306       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
1307
1308          else
1309             detr_star=0.
1310          endif
1311
1312! 1.3 dans les autres cas
1313          detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
1314     &                      *(zlev(ig,l+1)-zlev(ig,l))
1315
1316          coefc=min(zqla(ig,l-1)/1.e-3,1.)
1317          if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
1318          coefc=1.
1319! il semble qu'il soit important de baser le calcul sur
1320! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
1321          detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
1322
1323        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
1324
1325#endif
1326
1327
1328        if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
1329!IM 730508 beg
1330!        if(itap.GE.7200) THEN
1331!         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
1332!        endif
1333!IM 730508 end
1334         
1335         zqla0(ig,l)=zqla_est(ig,l)
1336         detr_star0(ig,l)=detr_star(ig,l)
1337!IM 060508 beg
1338!         if(detr_star(ig,l).GT.1.) THEN
1339!          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
1340!   &      ,detr_starc(ig,l),coefc
1341!         endif
1342!IM 060508 end
1343!IM 160508 beg
1344!IM 160508       IF (f0(ig).NE.0.) THEN
1345           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
1346!IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
1347!IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
1348!IM 160508       ELSE
1349!IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
1350!IM 160508       ENDIF
1351!IM 160508 end
1352!IM 060508 beg
1353!        if(detr_star(ig,l).GT.1.) THEN
1354!         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
1355!   &     REAL(1)/f0(ig)
1356!        endif
1357!IM 060508 end
1358        if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
1359!
1360!calcul de entr_star
1361
1362! #undef test2
1363! #ifdef test2
1364! La version test2 destabilise beaucoup le modele.
1365! Il semble donc que ca aide d'avoir un entrainement important sous
1366! le nuage.
1367!         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
1368!          entr_star(ig,l)=0.4*detr_star(ig,l)
1369!         else
1370!          entr_star(ig,l)=0.
1371!         endif
1372! #else
1373!
1374! Deplacement du calcul de entr_star pour eviter d'avoir aussi
1375! entr_star > fstar.
1376! Redeplacer suite a la transformation du cas detr>f
1377! FH
1378
1379        if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
1380#define int1d
1381!FH 070508 #define int1d4
1382!#undef int1d4
1383! L'option int1d4 correspond au choix dans le cas ou le detrainement
1384! devient trop grand.
1385
1386#ifdef int1d
1387
1388#ifdef int1d4
1389#else
1390       detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
1391!FH 070508 plus
1392       detr_star(ig,l)=min(detr_star(ig,l),1.)
1393#endif
1394
1395       entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
1396
1397        if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
1398#ifdef int1d4
1399! Si le detrainement excede le flux en bas + l'entrainement, le thermique
1400! doit disparaitre.
1401       if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
1402          detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
1403          f_star(ig,l+1)=0.
1404          linter(ig)=l+1
1405          zw2(ig,l+1)=-1.e-10
1406       endif
1407#endif
1408
1409
1410#else
1411
1412        if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
1413        if(l.gt.lalim(ig)) then
1414         entr_star(ig,l)=0.4*detr_star(ig,l)
1415        else
1416
1417! FH :
1418! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
1419! en haut de la couche d'alimentation.
1420! A remettre en questoin a la premiere occasion mais ca peut aider a
1421! ecrire un code robuste.
1422! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
1423! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
1424! d'alimentation, ce qui n'est pas forcement heureux.
1425
1426        if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
1427#undef pre_int1c
1428#ifdef pre_int1c
1429         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
1430         detr_star(ig,l)=entr_star(ig,l)
1431#else
1432         entr_star(ig,l)=0.
1433#endif
1434
1435        endif
1436
1437#endif
1438
1439        if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
1440        entr_star1(ig,l)=entr_star(ig,l)
1441        detr_star1(ig,l)=detr_star(ig,l)
1442!
1443
1444#ifdef int1d
1445#else
1446        if (detr_star(ig,l).gt.f_star(ig,l)) then
1447
1448!  Ce test est là entre autres parce qu'on passe par des valeurs
1449!  delirantes de detr_star.
1450!  ca vaut sans doute le coup de verifier pourquoi.
1451
1452           detr_star(ig,l)=f_star(ig,l)
1453#ifdef pre_int1c
1454           if (l.gt.lalim(ig)+1) then
1455               entr_star(ig,l)=0.
1456               alim_star(ig,l)=0.
1457! FH ajout pour forcer a stoper le thermique juste sous le sommet
1458! de la couche (voir calcul de finter)
1459               zw2(ig,l+1)=-1.e-10
1460               linter(ig)=l+1
1461            else
1462               entr_star(ig,l)=0.4*detr_star(ig,l)
1463            endif
1464#else
1465           entr_star(ig,l)=0.4*detr_star(ig,l)
1466#endif
1467        endif
1468#endif
1469
1470      else !l > 2
1471         detr_star(ig,l)=0.
1472         entr_star(ig,l)=0.
1473      endif
1474
1475        entr_star2(ig,l)=entr_star(ig,l)
1476        detr_star2(ig,l)=detr_star(ig,l)
1477        if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
1478
1479       endif  ! iflag_thermals_ed==0
1480
1481!CR:nvlle def de entr_star et detr_star
1482      if (iflag_thermals_ed>=1) then
1483!      if (l.lt.lalim(ig)) then
1484!      if (l.lt.2) then
1485!        entr_star(ig,l)=0.
1486!        detr_star(ig,l)=0.
1487!      else
1488!      if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then
1489!         entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1490!      else
1491!         entr_star(ig,l)=  &
1492!     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
1493!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
1494!     &                *(zlev(ig,l+1)-zlev(ig,l))
1495
1496 
1497         entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &         
1498     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
1499     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)   &
1500     &                *(zlev(ig,l+1)-zlev(ig,l))) &
1501     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1502
1503        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
1504            alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
1505            lalim(ig)=lmix_bis(ig)
1506            if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
1507        endif
1508
1509        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
1510!        c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
1511         c2(ig,l)=0.001
1512         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
1513     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
1514     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
1515     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
1516     &                *(zlev(ig,l+1)-zlev(ig,l)))                    &
1517     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1518
1519       else
1520!         c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
1521          c2(ig,l)=0.003
1522
1523         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
1524     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
1525     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
1526     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
1527     &                *(zlev(ig,l+1)-zlev(ig,l))) &
1528     &                +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1529       endif
1530         
1531           
1532!        detr_star(ig,l)=detr_star(ig,l)*3.
1533!        if (l.lt.lalim(ig)) then
1534!          entr_star(ig,l)=0.
1535!        endif
1536!        if (l.lt.2) then
1537!          entr_star(ig,l)=0.
1538!          detr_star(ig,l)=0.
1539!        endif
1540
1541
1542!      endif
1543!      else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
1544!      entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1))        &
1545!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
1546!     &                *(zlev(ig,l+1)-zlev(ig,l))
1547!      detr_star(ig,l)=0.002*f_star(ig,l)                         &
1548!     &                *(zlev(ig,l+1)-zlev(ig,l))
1549!      else
1550!      entr_star(ig,l)=0.001*f_star(ig,l)                         &
1551!     &                *(zlev(ig,l+1)-zlev(ig,l))
1552!      detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1))       &
1553!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))       &
1554!     &                *(zlev(ig,l+1)-zlev(ig,l))                      &
1555!     &                +0.002*f_star(ig,l)                             &
1556!     &                *(zlev(ig,l+1)-zlev(ig,l))
1557!      endif
1558
1559      endif   ! iflag_thermals_ed==1
1560
1561!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1562! FH inutile si on conserve comme on l'a fait plus haut entr=detr
1563! dans la couche d'alimentation
1564!pas d entrainement dans la couche alim
1565!      if ((l.le.lalim(ig))) then
1566!           entr_star(ig,l)=0.
1567!      endif
1568!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1569!
1570!prise en compte du detrainement et de l entrainement dans le calcul du flux
1571
1572      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
1573     &              -detr_star(ig,l)
1574
1575!test sur le signe de f_star
1576        if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
1577       if (f_star(ig,l+1).gt.1.e-10) then
1578!----------------------------------------------------------------------------
1579!calcul de la vitesse verticale en melangeant Tl et qt du thermique
1580!---------------------------------------------------------------------------
1581!
1582       Zsat=.false.
1583       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
1584     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
1585     &            /(f_star(ig,l+1)+detr_star(ig,l))
1586!
1587       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
1588     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
1589     &            /(f_star(ig,l+1)+detr_star(ig,l))
1590
1591               Tbef=ztla(ig,l)*zpspsk(ig,l)
1592               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1593               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
1594               qsatbef=MIN(0.5,qsatbef)
1595               zcor=1./(1.-retv*qsatbef)
1596               qsatbef=qsatbef*zcor
1597               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
1598               if (Zsat) then
1599               qlbef=max(0.,zqta(ig,l)-qsatbef)
1600               DT = 0.5*RLvCp*qlbef
1601               do while (abs(DT).gt.DDT0)
1602                 Tbef=Tbef+DT
1603                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1604                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
1605                 qsatbef=MIN(0.5,qsatbef)
1606                 zcor=1./(1.-retv*qsatbef)
1607                 qsatbef=qsatbef*zcor
1608                 qlbef=zqta(ig,l)-qsatbef
1609
1610                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1611                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
1612                 zcor=1./(1.-retv*qsatbef)
1613                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
1614                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
1615                 denom=1.+RLvCp*dqsat_dT
1616                 DT=num/denom
1617              enddo
1618                 zqla(ig,l) = max(0.,qlbef)
1619              endif
1620!   
1621        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
1622! on ecrit de maniere conservative (sat ou non)
1623!          T = Tl +Lv/Cp ql
1624           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1625           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1626!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1627           zha(ig,l) = ztva(ig,l)
1628           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
1629     &              -zqla(ig,l))-zqla(ig,l))
1630
1631!on ecrit zqsat
1632           zqsatth(ig,l)=qsatbef 
1633!calcul de vitesse
1634           zw2(ig,l+1)=zw2(ig,l)*  &
1635     &                 ((f_star(ig,l))**2)  &
1636!  Tests de Catherine
1637!     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
1638     &      /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
1639     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
1640     &                 *fact_gamma &
1641     &                 *(zlev(ig,l+1)-zlev(ig,l))
1642!prise en compte des forces de pression que qd flottabilité<0
1643!              zw2(ig,l+1)=zw2(ig,l)*  &
1644!     &            1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &       
1645!     &                 (f_star(ig,l))**2 &
1646!     &                 /(f_star(ig,l)+entr_star(ig,l))**2+ &
1647!     &                 (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+  &       
1648!     &                 /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
1649!     &                 /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
1650!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
1651!     &                 *1./3. &
1652!     &                 *(zlev(ig,l+1)-zlev(ig,l))
1653         
1654!        write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
1655!     &              -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
1656!     &               2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1657
1658 
1659!             zw2(ig,l+1)=zw2(ig,l)*  &
1660!     &                 (2.-2.*entr_star(ig,l)/f_star(ig,l)) & 
1661!     &                 -zw2(ig,l-1)+  &       
1662!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
1663!     &                 *1./3. &
1664!     &                 *(zlev(ig,l+1)-zlev(ig,l))             
1665
1666            endif
1667        endif
1668        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1669!
1670!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1671
1672            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1673                print*,'On tombe sur le cas particulier de thermcell_plume'
1674                zw2(ig,l+1)=0.
1675                linter(ig)=l+1
1676            endif
1677
1678!        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
1679        if (zw2(ig,l+1).lt.0.) then
1680           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
1681     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1682           zw2(ig,l+1)=0.
1683        endif
1684
1685           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1686
1687        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1688!   lmix est le niveau de la couche ou w (wa_moy) est maximum
1689!on rajoute le calcul de lmix_bis
1690            if (zqla(ig,l).lt.1.e-10) then
1691               lmix_bis(ig)=l+1
1692            endif
1693            lmix(ig)=l+1
1694            wmaxa(ig)=wa_moy(ig,l+1)
1695        endif
1696        enddo
1697      enddo
1698
1699!on remplace a* par e* ds premiere couche
1700!      if (iflag_thermals_ed.ge.1) then
1701!       do ig=1,ngrid
1702!       do l=2,klev
1703!          if (l.lt.lalim(ig)) then
1704!             alim_star(ig,l)=entr_star(ig,l)
1705!          endif
1706!       enddo
1707!       enddo
1708!       do ig=1,ngrid
1709!          lalim(ig)=lmix_bis(ig)
1710!       enddo
1711!      endif
1712       if (iflag_thermals_ed.ge.1) then
1713          do ig=1,ngrid
1714             do l=2,lalim(ig)
1715                alim_star(ig,l)=entr_star(ig,l)
1716                entr_star(ig,l)=0.
1717             enddo
1718           enddo
1719       endif
1720        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1721
1722!     print*,'thermcell_plume OK'
1723
1724      return
1725      end
1726!==============================================================================
1727       SUBROUTINE thermcellV0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
1728     &                            lalim,lmin,zmax,wmax,lev_out)
1729
1730!--------------------------------------------------------------------------
1731!thermcell_dry: calcul de zmax et wmax du thermique sec
1732!--------------------------------------------------------------------------
1733       IMPLICIT NONE
1734#include "YOMCST.h"       
1735#include "iniprint.h"
1736       INTEGER l,ig
1737
1738       INTEGER ngrid,nlay
1739       REAL zlev(ngrid,nlay+1)
1740       REAL pphi(ngrid,nlay)
1741       REAl ztv(ngrid,nlay)
1742       REAL alim_star(ngrid,nlay)
1743       INTEGER lalim(ngrid)
1744      integer lev_out                           ! niveau pour les print
1745
1746       REAL zmax(ngrid)
1747       REAL wmax(ngrid)
1748
1749!variables locales
1750       REAL zw2(ngrid,nlay+1)
1751       REAL f_star(ngrid,nlay+1)
1752       REAL ztva(ngrid,nlay+1)
1753       REAL wmaxa(ngrid)
1754       REAL wa_moy(ngrid,nlay+1)
1755       REAL linter(ngrid),zlevinter(ngrid)
1756       INTEGER lmix(ngrid),lmax(ngrid),lmin(ngrid)
1757
1758!initialisations
1759       do ig=1,ngrid
1760          do l=1,nlay+1
1761             zw2(ig,l)=0.
1762             wa_moy(ig,l)=0.
1763          enddo
1764       enddo
1765       do ig=1,ngrid
1766          do l=1,nlay
1767             ztva(ig,l)=ztv(ig,l)
1768          enddo
1769       enddo
1770       do ig=1,ngrid
1771          wmax(ig)=0.
1772          wmaxa(ig)=0.
1773       enddo
1774!calcul de la vitesse a partir de la CAPE en melangeant thetav
1775
1776!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1777! A eliminer
1778! Ce if complique etait fait pour reperer la premiere couche instable
1779! Ici, c'est lmin.
1780!
1781!       do l=1,nlay-2
1782!         do ig=1,ngrid
1783!            if (ztv(ig,l).gt.ztv(ig,l+1)  &
1784!     &         .and.alim_star(ig,l).gt.1.e-10  &
1785!     &         .and.zw2(ig,l).lt.1e-10) then
1786!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1787
1788
1789! Calcul des F^*, integrale verticale de E^*
1790       f_star(:,1)=0.
1791       do l=1,nlay
1792          f_star(:,l+1)=f_star(:,l)+alim_star(:,l)
1793       enddo
1794
1795! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
1796       linter(:)=0.
1797
1798! couche la plus haute concernee par le thermique.
1799       lmax(:)=1
1800
1801! Le niveau linter est une variable continue qui se trouve dans la couche
1802! lmax
1803
1804       do l=1,nlay-2
1805         do ig=1,ngrid
1806            if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
1807
1808!------------------------------------------------------------------------
1809!  Calcul de la vitesse en haut de la premiere couche instable.
1810!  Premiere couche du panache thermique
1811!------------------------------------------------------------------------
1812               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
1813     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
1814     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1815
1816!------------------------------------------------------------------------
1817! Tant que la vitesse en bas de la couche et la somme du flux de masse
1818! et de l'entrainement (c'est a dire le flux de masse en haut) sont
1819! positifs, on calcul
1820! 1. le flux de masse en haut  f_star(ig,l+1)
1821! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
1822! 3. la vitesse au carré en haut zw2(ig,l+1)
1823!------------------------------------------------------------------------
1824
1825!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1826!  A eliminer : dans cette version, si zw2 est > 0 on a un therique.
1827!  et donc, au dessus, f_star(ig,l+1) est forcement suffisamment
1828!  grand puisque on n'a pas de detrainement.
1829!  f_star est une fonction croissante.
1830!  c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
1831!           else if ((zw2(ig,l).ge.1e-10).and.  &
1832!    &               (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
1833!              f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
1834!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1835
1836            else if (zw2(ig,l).ge.1e-10) then
1837
1838               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l)  &
1839     &                    *ztv(ig,l))/f_star(ig,l+1)
1840               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+  &
1841     &                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
1842     &                     *(zlev(ig,l+1)-zlev(ig,l))
1843            endif
1844! determination de zmax continu par interpolation lineaire
1845!------------------------------------------------------------------------
1846
1847            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1848!               print*,'On tombe sur le cas particulier de thermcell_dry'
1849                zw2(ig,l+1)=0.
1850                linter(ig)=l+1
1851                lmax(ig)=l
1852            endif
1853
1854            if (zw2(ig,l+1).lt.0.) then
1855               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
1856     &           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1857               zw2(ig,l+1)=0.
1858               lmax(ig)=l
1859            endif
1860
1861               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1862
1863            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1864!   lmix est le niveau de la couche ou w (wa_moy) est maximum
1865               lmix(ig)=l+1
1866               wmaxa(ig)=wa_moy(ig,l+1)
1867            endif
1868         enddo
1869      enddo
1870       if (prt_level.ge.1) print*,'fin calcul zw2'
1871!
1872!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1873! A eliminer :
1874! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
1875! Calcul de la couche correspondant a la hauteur du thermique
1876!      do ig=1,ngrid
1877!         lmax(ig)=lalim(ig)
1878!      enddo
1879!      do ig=1,ngrid
1880!         do l=nlay,lalim(ig)+1,-1
1881!            if (zw2(ig,l).le.1.e-10) then
1882!               lmax(ig)=l-1
1883!            endif
1884!         enddo
1885!      enddo
1886!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1887
1888!   
1889! Determination de zw2 max
1890      do ig=1,ngrid
1891         wmax(ig)=0.
1892      enddo
1893
1894      do l=1,nlay
1895         do ig=1,ngrid
1896            if (l.le.lmax(ig)) then
1897                zw2(ig,l)=sqrt(zw2(ig,l))
1898                wmax(ig)=max(wmax(ig),zw2(ig,l))
1899            else
1900                 zw2(ig,l)=0.
1901            endif
1902          enddo
1903      enddo
1904
1905!   Longueur caracteristique correspondant a la hauteur des thermiques.
1906      do  ig=1,ngrid
1907         zmax(ig)=0.
1908         zlevinter(ig)=zlev(ig,1)
1909      enddo
1910      do  ig=1,ngrid
1911! calcul de zlevinter
1912
1913!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1914! FH A eliminer
1915! Simplification
1916!          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
1917!     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
1918!     &    -zlev(ig,lmax(ig)))
1919!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1920
1921          zlevinter(ig)=zlev(ig,lmax(ig)) + &
1922     &    (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
1923           zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1924      enddo
1925
1926! Verification que lalim<=lmax
1927      do ig=1,ngrid
1928         if(lalim(ig)>lmax(ig)) then
1929           if ( prt_level > 1 ) THEN
1930            print*,'WARNING thermcell_dry ig=',ig,'  lalim=',lalim(ig),'  lmax(ig)=',lmax(ig)
1931           endif
1932           lmax(ig)=lalim(ig)
1933         endif
1934      enddo
1935     
1936      RETURN
1937      END
1938!==============================================================================
1939      SUBROUTINE thermcellV0_init(ngrid,nlay,ztv,zlay,zlev,  &
1940     &                  lalim,lmin,alim_star,alim_star_tot,lev_out)
1941
1942!----------------------------------------------------------------------
1943!thermcell_init: calcul du profil d alimentation du thermique
1944!----------------------------------------------------------------------
1945      IMPLICIT NONE
1946#include "iniprint.h"
1947#include "thermcell.h"
1948
1949      INTEGER l,ig
1950!arguments d entree
1951      INTEGER ngrid,nlay
1952      REAL ztv(ngrid,nlay)
1953      REAL zlay(ngrid,nlay)
1954      REAL zlev(ngrid,nlay+1)
1955!arguments de sortie
1956      INTEGER lalim(ngrid)
1957      INTEGER lmin(ngrid)
1958      REAL alim_star(ngrid,nlay)
1959      REAL alim_star_tot(ngrid)
1960      integer lev_out                           ! niveau pour les print
1961     
1962      REAL zzalim(ngrid)
1963!CR: ponderation entrainement des couches instables
1964!def des alim_star tels que alim=f*alim_star     
1965
1966      do l=1,nlay
1967         do ig=1,ngrid
1968            alim_star(ig,l)=0.
1969         enddo
1970      enddo
1971! determination de la longueur de la couche d entrainement
1972      do ig=1,ngrid
1973         lalim(ig)=1
1974      enddo
1975
1976      if (iflag_thermals_ed.ge.1) then
1977!si la premiÚre couche est instable, on declenche un thermique
1978         do ig=1,ngrid
1979            if (ztv(ig,1).gt.ztv(ig,2)) then
1980               lmin(ig)=1
1981               lalim(ig)=2
1982               alim_star(ig,1)=1.
1983               alim_star_tot(ig)=alim_star(ig,1)
1984               if(prt_level.GE.10) print*,'init',alim_star(ig,1),alim_star_tot(ig)
1985            else
1986                lmin(ig)=1
1987                lalim(ig)=1
1988                alim_star(ig,1)=0.
1989                alim_star_tot(ig)=0.
1990            endif
1991         enddo
1992 
1993         else
1994!else iflag_thermals_ed=0 ancienne def de l alim
1995
1996!on ne considere que les premieres couches instables
1997      do l=nlay-2,1,-1
1998         do ig=1,ngrid
1999            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
2000     &          ztv(ig,l+1).le.ztv(ig,l+2)) then
2001               lalim(ig)=l+1
2002            endif
2003          enddo
2004      enddo
2005
2006! determination du lmin: couche d ou provient le thermique
2007
2008      do ig=1,ngrid
2009! FH initialisation de lmin a nlay plutot que 1.
2010!        lmin(ig)=nlay
2011         lmin(ig)=1
2012      enddo
2013      do l=nlay,2,-1
2014         do ig=1,ngrid
2015            if (ztv(ig,l-1).gt.ztv(ig,l)) then
2016               lmin(ig)=l-1
2017            endif
2018         enddo
2019      enddo
2020!
2021      zzalim(:)=0.
2022      do l=1,nlay-1
2023         do ig=1,ngrid
2024             if (l<lalim(ig)) then
2025                zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
2026             endif
2027          enddo
2028      enddo
2029      do ig=1,ngrid
2030          if (lalim(ig)>1) then
2031             zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
2032          else
2033             zzalim(ig)=zlay(ig,1)
2034          endif
2035      enddo
2036
2037      if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))
2038
2039! definition de l'entrainement des couches
2040      if (1.eq.1) then
2041      do l=1,nlay-1
2042         do ig=1,ngrid
2043            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
2044     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
2045!def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
2046             alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
2047     &                       *sqrt(zlev(ig,l+1))
2048            endif
2049         enddo
2050      enddo
2051      else
2052      do l=1,nlay-1
2053         do ig=1,ngrid
2054            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
2055     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
2056             alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
2057     &        *(zlev(ig,l+1)-zlev(ig,l))
2058            endif
2059         enddo
2060      enddo
2061      endif
2062     
2063! pas de thermique si couche 1 stable
2064      do ig=1,ngrid
2065!CRnouveau test
2066        if (alim_star(ig,1).lt.1.e-10) then
2067            do l=1,nlay
2068                alim_star(ig,l)=0.
2069            enddo
2070            lmin(ig)=1
2071         endif
2072      enddo
2073! calcul de l alimentation totale
2074      do ig=1,ngrid
2075         alim_star_tot(ig)=0.
2076      enddo
2077      do l=1,nlay
2078         do ig=1,ngrid
2079            alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
2080         enddo
2081      enddo
2082!
2083! Calcul entrainement normalise
2084      do l=1,nlay
2085         do ig=1,ngrid
2086            if (alim_star_tot(ig).gt.1.e-10) then
2087               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
2088            endif
2089         enddo
2090      enddo
2091       
2092!on remet alim_star_tot a 1
2093      do ig=1,ngrid
2094         alim_star_tot(ig)=1.
2095      enddo
2096
2097      endif
2098!endif iflag_thermals_ed
2099      return
2100      end 
2101!==============================================================================
Note: See TracBrowser for help on using the repository browser.