source: LMDZ4/branches/LMDZ4V5.0-LF/libf/phylmd/thermcell_main.F90 @ 5360

Last change on this file since 5360 was 1311, checked in by Laurent Fairhead, 15 years ago

Modifications to thermals for TKE transport


Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.3 KB
Line 
1!
2! $Id: thermcell_main.F90 1311 2010-02-18 13:14:02Z abarral $
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     &                  ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed &
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!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
25!
26!   le thermique est suppose homogene et dissipe par melange avec
27!   son environnement. la longueur l_mix controle l'efficacite du
28!   melange
29!
30!   Le calcul du transport des differentes especes 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      integer iflag_thermals_ed
58      real ptimestep,l_mix,r_aspect
59      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
60      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
61      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
62      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
63      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
64      real pphi(ngrid,nlay)
65
66!   local:
67!   ------
68
69      integer icount
70      data icount/0/
71      save icount
72!$OMP THREADPRIVATE(icount)
73
74      integer,save :: igout=1
75!$OMP THREADPRIVATE(igout)
76      integer,save :: lunout1=6
77!$OMP THREADPRIVATE(lunout1)
78      integer,save :: lev_out=10
79!$OMP THREADPRIVATE(lev_out)
80
81      INTEGER ig,k,l,ll
82      real zsortie1d(klon)
83      INTEGER lmax(klon),lmin(klon),lalim(klon)
84      INTEGER lmix(klon)
85      INTEGER lmix_bis(klon)
86      real linter(klon)
87      real zmix(klon)
88      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
89!      real fraca(klon,klev)
90
91      real zmax_sec(klon)
92!on garde le zmax du pas de temps precedent
93      real zmax0(klon)
94!FH/IM     save zmax0
95
96      real lambda
97
98      real zlev(klon,klev+1),zlay(klon,klev)
99      real deltaz(klon,klev)
100      REAL zh(klon,klev)
101      real zthl(klon,klev),zdthladj(klon,klev)
102      REAL ztv(klon,klev)
103      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
104      real zl(klon,klev)
105      real zsortie(klon,klev)
106      real zva(klon,klev)
107      real zua(klon,klev)
108      real zoa(klon,klev)
109
110      real zta(klon,klev)
111      real zha(klon,klev)
112      real fraca(klon,klev+1)
113      real zf,zf2
114      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
115      real q2(klon,klev)
116! FH probleme de dimensionnement avec l'allocation dynamique
117!     common/comtherm/thetath2,wth2
118      real wq(klon,klev)
119      real wthl(klon,klev)
120      real wthv(klon,klev)
121   
122      real ratqscth(klon,klev)
123      real var
124      real vardiff
125      real ratqsdiff(klon,klev)
126
127      logical sorties
128      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
129      real zpspsk(klon,klev)
130
131      real wmax(klon)
132      real wmax_sec(klon)
133      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
134      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
135
136      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
137!niveau de condensation
138      integer nivcon(klon)
139      real zcon(klon)
140      REAL CHI
141      real zcon2(klon)
142      real pcon(klon)
143      real zqsat(klon,klev)
144      real zqsatth(klon,klev)
145
146      real f_star(klon,klev+1),entr_star(klon,klev)
147      real detr_star(klon,klev)
148      real alim_star_tot(klon)
149      real alim_star(klon,klev)
150      real alim_star_clos(klon,klev)
151      real f(klon), f0(klon)
152!FH/IM     save f0
153      real zlevinter(klon)
154      logical debut
155       real seuil
156      real csc(klon,klev)
157
158!
159      !nouvelles variables pour la convection
160      real Ale_bl(klon)
161      real Alp_bl(klon)
162      real alp_int(klon)
163      real ale_int(klon)
164      integer n_int(klon)
165      real fm_tot(klon)
166      real wght_th(klon,klev)
167      integer lalim_conv(klon)
168!v1d     logical therm
169!v1d     save therm
170
171      character*2 str2
172      character*10 str10
173
174      character (len=20) :: modname='thermcell_main'
175      character (len=80) :: abort_message
176
177      EXTERNAL SCOPY
178!
179
180!-----------------------------------------------------------------------
181!   initialisation:
182!   ---------------
183!
184
185      seuil=0.25
186
187      if (debut)  then
188         fm0=0.
189         entr0=0.
190         detr0=0.
191
192
193! #define wrgrads_thermcell
194#undef wrgrads_thermcell
195#ifdef wrgrads_thermcell
196! Initialisation des sorties grads pour les thermiques.
197! Pour l'instant en 1D sur le point igout.
198! Utilise par thermcell_out3d.h
199         str10='therm'
200         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
201     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
202     &   ptimestep,str10,'therm ')
203#endif
204
205
206
207      endif
208
209      fm=0. ; entr=0. ; detr=0.
210
211      icount=icount+1
212
213!IM 090508 beg
214!print*,'====================================================================='
215!print*,'====================================================================='
216!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
217!print*,'====================================================================='
218!print*,'====================================================================='
219!IM 090508 end
220
221      if (prt_level.ge.1) print*,'thermcell_main V4'
222
223       sorties=.true.
224      IF(ngrid.NE.klon) THEN
225         PRINT*
226         PRINT*,'STOP dans convadj'
227         PRINT*,'ngrid    =',ngrid
228         PRINT*,'klon  =',klon
229      ENDIF
230!
231!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
232     do ig=1,klon
233      if (prt_level.ge.20) then
234       print*,'th_main ig f0',ig,f0(ig)
235      endif
236         f0(ig)=max(f0(ig),1.e-2)
237         zmax0(ig)=max(zmax0(ig),40.)
238!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
239     enddo
240
241!-----------------------------------------------------------------------
242! Calcul de T,q,ql a partir de Tl et qT dans l environnement
243!   --------------------------------------------------------------------
244!
245      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
246     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
247       
248      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
249
250!------------------------------------------------------------------------
251!                       --------------------
252!
253!
254!                       + + + + + + + + + + +
255!
256!
257!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
258!  wh,wt,wo ...
259!
260!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
261!
262!
263!                       --------------------   zlev(1)
264!                       \\\\\\\\\\\\\\\\\\\\
265!
266!
267
268!-----------------------------------------------------------------------
269!   Calcul des altitudes des couches
270!-----------------------------------------------------------------------
271
272      do l=2,nlay
273         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
274      enddo
275         zlev(:,1)=0.
276         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
277      do l=1,nlay
278         zlay(:,l)=pphi(:,l)/RG
279      enddo
280!calcul de l epaisseur des couches
281      do l=1,nlay
282         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
283      enddo
284
285!     print*,'2 OK convect8'
286!-----------------------------------------------------------------------
287!   Calcul des densites
288!-----------------------------------------------------------------------
289
290      do l=1,nlay
291         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
292      enddo
293
294!IM
295     if (prt_level.ge.10)write(lunout,*)                                &
296    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
297      rhobarz(:,1)=rho(:,1)
298
299      do l=2,nlay
300         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
301      enddo
302
303!calcul de la masse
304      do l=1,nlay
305         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
306      enddo
307
308      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
309
310!------------------------------------------------------------------
311!
312!             /|\
313!    --------  |  F_k+1 -------   
314!                              ----> D_k
315!             /|\              <---- E_k , A_k
316!    --------  |  F_k ---------
317!                              ----> D_k-1
318!                              <---- E_k-1 , A_k-1
319!
320!
321!
322!
323!
324!    ---------------------------
325!
326!    ----- F_lmax+1=0 ----------         \
327!            lmax     (zmax)              |
328!    ---------------------------          |
329!                                         |
330!    ---------------------------          |
331!                                         |
332!    ---------------------------          |
333!                                         |
334!    ---------------------------          |
335!                                         |
336!    ---------------------------          |
337!                                         |  E
338!    ---------------------------          |  D
339!                                         |
340!    ---------------------------          |
341!                                         |
342!    ---------------------------  \       |
343!            lalim                 |      |
344!    ---------------------------   |      |
345!                                  |      |
346!    ---------------------------   |      |
347!                                  | A    |
348!    ---------------------------   |      |
349!                                  |      |
350!    ---------------------------   |      |
351!    lmin  (=1 pour le moment)     |      |
352!    ----- F_lmin=0 ------------  /      /
353!
354!    ---------------------------
355!    //////////////////////////
356!
357!
358!=============================================================================
359!  Calculs initiaux ne faisant pas intervenir les changements de phase
360!=============================================================================
361
362!------------------------------------------------------------------
363!  1. alim_star est le profil vertical de l'alimentation a la base du
364!     panache thermique, calcule a partir de la flotabilite de l'air sec
365!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
366!------------------------------------------------------------------
367!
368      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
369      lmin=1
370
371!-----------------------------------------------------------------------------
372!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
373!     panache sec conservatif (e=d=0) alimente selon alim_star
374!     Il s'agit d'un calcul de type CAPE
375!     zmax_sec est utilise pour determiner la geometrie du thermique.
376!------------------------------------------------------------------------------
377!---------------------------------------------------------------------------------
378!calcul du melange et des variables dans le thermique
379!--------------------------------------------------------------------------------
380!
381      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
382!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
383
384! Gestion temporaire de plusieurs appels à thermcell_plume au travers
385! de la variable iflag_thermals
386
387!      print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
388      if (iflag_thermals_ed<=9) then
389!         print*,'THERM NOVUELLE/NOUVELLE/ANCIENNE'
390         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
391     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
392     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
393     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
394     &    ,lev_out,lunout1,igout)
395
396      elseif (iflag_thermals_ed<=19) then
397! Version d'Arnaud Jam
398!         print*,'THERM RIO et al 2010'
399         CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
400     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
401     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
402     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
403     &    ,lev_out,lunout1,igout)
404
405      endif
406
407      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
408
409      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
410      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
411
412      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
413      if (prt_level.ge.10) then
414         write(lunout1,*) 'Dans thermcell_main 2'
415         write(lunout1,*) 'lmin ',lmin(igout)
416         write(lunout1,*) 'lalim ',lalim(igout)
417         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
418         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
419     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
420      endif
421
422!-------------------------------------------------------------------------------
423! Calcul des caracteristiques du thermique:zmax,zmix,wmax
424!-------------------------------------------------------------------------------
425!
426      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
427     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
428
429
430
431      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
432      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
433      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
434      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
435
436      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
437
438!-------------------------------------------------------------------------------
439! Fermeture,determination de f
440!-------------------------------------------------------------------------------
441!
442!
443!!      write(lunout,*)'THERM NOUVEAU XXXXX'
444      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
445    &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
446
447call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
448call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
449
450      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
451      if (prt_level.ge.10) then
452         write(lunout1,*) 'Dans thermcell_main 1b'
453         write(lunout1,*) 'lmin ',lmin(igout)
454         write(lunout1,*) 'lalim ',lalim(igout)
455         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
456         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
457     &    ,l=1,lalim(igout)+4)
458      endif
459
460
461
462!print*,'THERM 26JJJ'
463
464! Choix de la fonction d'alimentation utilisee pour la fermeture.
465! Apparemment sans importance
466      alim_star_clos(:,:)=alim_star(:,:)
467      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
468
469! Appel avec la version seche
470      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
471     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
472
473!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
474! Appel avec les zmax et wmax tenant compte de la condensation
475! Semble moins bien marcher
476!     CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
477!    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
478!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
479
480      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
481
482      if (tau_thermals>1.) then
483         lambda=exp(-ptimestep/tau_thermals)
484         f0=(1.-lambda)*f+lambda*f0
485      else
486         f0=f
487      endif
488
489! Test valable seulement en 1D mais pas genant
490      if (.not. (f0(1).ge.0.) ) then
491              abort_message = '.not. (f0(1).ge.0.)'
492              CALL abort_gcm (modname,abort_message,1)
493      endif
494
495!-------------------------------------------------------------------------------
496!deduction des flux
497!-------------------------------------------------------------------------------
498
499      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
500     &       lalim,lmax,alim_star,  &
501     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
502     &       detr,zqla,lev_out,lunout1,igout)
503!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
504
505      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
506      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
507      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
508
509!------------------------------------------------------------------
510!   On ne prend pas directement les profils issus des calculs precedents
511!   mais on s'autorise genereusement une relaxation vers ceci avec
512!   une constante de temps tau_thermals (typiquement 1800s).
513!------------------------------------------------------------------
514
515      if (tau_thermals>1.) then
516         lambda=exp(-ptimestep/tau_thermals)
517         fm0=(1.-lambda)*fm+lambda*fm0
518         entr0=(1.-lambda)*entr+lambda*entr0
519         detr0=(1.-lambda)*detr+lambda*detr0
520      else
521         fm0=fm
522         entr0=entr
523         detr0=detr
524      endif
525
526!c------------------------------------------------------------------
527!   calcul du transport vertical
528!------------------------------------------------------------------
529
530      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
531     &                    zthl,zdthladj,zta,lev_out)
532      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
533     &                   po,pdoadj,zoa,lev_out)
534
535!------------------------------------------------------------------
536! Calcul de la fraction de l'ascendance
537!------------------------------------------------------------------
538      do ig=1,klon
539         fraca(ig,1)=0.
540         fraca(ig,nlay+1)=0.
541      enddo
542      do l=2,nlay
543         do ig=1,klon
544            if (zw2(ig,l).gt.1.e-10) then
545            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
546            else
547            fraca(ig,l)=0.
548            endif
549         enddo
550      enddo
551     
552!------------------------------------------------------------------
553!  calcul du transport vertical du moment horizontal
554!------------------------------------------------------------------
555
556!IM 090508 
557      if (1.eq.1) then
558!IM 070508 vers. _dq       
559!     if (1.eq.0) then
560
561
562! Calcul du transport de V tenant compte d'echange par gradient
563! de pression horizontal avec l'environnement
564
565         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
566     &    ,fraca,zmax  &
567     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
568
569      else
570
571! calcul purement conservatif pour le transport de V
572         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
573     &    ,zu,pduadj,zua,lev_out)
574         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
575     &    ,zv,pdvadj,zva,lev_out)
576      endif
577
578!     print*,'13 OK convect8'
579      do l=1,nlay
580         do ig=1,ngrid
581           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
582         enddo
583      enddo
584
585      if (prt_level.ge.1) print*,'14 OK convect8'
586!------------------------------------------------------------------
587!   Calculs de diagnostiques pour les sorties
588!------------------------------------------------------------------
589!calcul de fraca pour les sorties
590     
591      if (sorties) then
592      if (prt_level.ge.1) print*,'14a OK convect8'
593! calcul du niveau de condensation
594! initialisation
595      do ig=1,ngrid
596         nivcon(ig)=0
597         zcon(ig)=0.
598      enddo
599!nouveau calcul
600      do ig=1,ngrid
601      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
602      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
603      enddo
604      do k=1,nlay
605         do ig=1,ngrid
606         if ((pcon(ig).le.pplay(ig,k))  &
607     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
608            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
609         endif
610         enddo
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      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
642!
643      if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
644            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
645            if(zw2(ig,l).gt.1.e-10) then
646             wth2(ig,l)=zf2*(zw2(ig,l))**2
647            else
648             wth2(ig,l)=0.
649            endif
650!           print*,'wth2=',wth2(ig,l)
651            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
652     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
653      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
654            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
655!test: on calcul q2/po=ratqsc
656            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
657         enddo
658      enddo
659!calcul des flux: q, thetal et thetav
660      do l=1,nlay
661         do ig=1,ngrid
662      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
663      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
664      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
665         enddo
666      enddo
667!calcul de ale_bl et alp_bl
668!pour le calcul d'une valeur integree entre la surface et lmax
669      do ig=1,ngrid
670      alp_int(ig)=0.
671      ale_int(ig)=0.
672      n_int(ig)=0
673      enddo
674!
675      do l=1,nlay
676      do ig=1,ngrid
677       if(l.LE.lmax(ig)) THEN
678        alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
679        ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
680        n_int(ig)=n_int(ig)+1
681       endif
682      enddo
683      enddo
684!      print*,'avant calcul ale et alp'
685!calcul de ALE et ALP pour la convection
686      do ig=1,ngrid
687!      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
688!          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
689!      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig))
690!     &           *0.1
691!valeur integree de alp_bl * 0.5:
692       if (n_int(ig).gt.0) then
693       Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
694!       if (Alp_bl(ig).lt.0.) then
695!       Alp_bl(ig)=0.
696       endif
697!       endif
698!         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
699!     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
700!            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
701!      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
702!      if (nivcon(ig).eq.1) then
703!       Ale_bl(ig)=0.
704!       else
705!valeur max de ale_bl:
706       Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2
707!     & /2.
708!     & *0.1
709!        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
710!       if (n_int(ig).gt.0) then
711!       Ale_bl(ig)=ale_int(ig)/n_int(ig)
712!        Ale_bl(ig)=4.
713!       endif
714!       endif
715!            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
716!          Ale_bl(ig)=wth2(ig,nivcon(ig))
717!          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
718      enddo
719!test:calcul de la ponderation des couches pour KE
720!initialisations
721!      print*,'ponderation'
722      do ig=1,ngrid
723           fm_tot(ig)=0.
724      enddo
725       do ig=1,ngrid
726        do k=1,klev
727           wght_th(ig,k)=1.
728        enddo
729       enddo
730       do ig=1,ngrid
731!         lalim_conv(ig)=lmix_bis(ig)
732!la hauteur de la couche alim_conv = hauteur couche alim_therm
733         lalim_conv(ig)=lalim(ig)
734!         zentr(ig)=zlev(ig,lalim(ig))
735      enddo
736      do ig=1,ngrid
737        do k=1,lalim_conv(ig)
738           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
739        enddo
740      enddo
741      do ig=1,ngrid
742        do k=1,lalim_conv(ig)
743           if (fm_tot(ig).gt.1.e-10) then
744!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
745           endif
746!on pondere chaque couche par a*
747             if (alim_star(ig,k).gt.1.e-10) then
748             wght_th(ig,k)=alim_star(ig,k)
749             else
750             wght_th(ig,k)=1.
751             endif
752        enddo
753      enddo
754!      print*,'apres wght_th'
755!test pour prolonger la convection
756      do ig=1,ngrid
757!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
758      if ((alim_star(ig,1).lt.1.e-10)) then
759      lalim_conv(ig)=1
760      wght_th(ig,1)=1.
761!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
762      endif
763      enddo
764
765!calcul du ratqscdiff
766      if (prt_level.ge.1) print*,'14e OK convect8'
767      var=0.
768      vardiff=0.
769      ratqsdiff(:,:)=0.
770      do ig=1,ngrid
771         do l=1,lalim(ig)
772            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
773         enddo
774      enddo
775      if (prt_level.ge.1) print*,'14f OK convect8'
776      do ig=1,ngrid
777          do l=1,lalim(ig)
778          zf=fraca(ig,l)
779          zf2=zf/(1.-zf)
780          vardiff=vardiff+alim_star(ig,l)  &
781     &           *(zqta(ig,l)*1000.-var)**2
782!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
783!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
784          enddo
785      enddo
786      if (prt_level.ge.1) print*,'14g OK convect8'
787      do l=1,nlay
788         do ig=1,ngrid
789            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
790!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
791         enddo
792      enddo
793!--------------------------------------------------------------------   
794!
795!ecriture des fichiers sortie
796!     print*,'15 OK convect8'
797
798#ifdef wrgrads_thermcell
799      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
800#include "thermcell_out3d.h"
801#endif
802
803      endif
804
805      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
806
807      return
808      end
809
810!-----------------------------------------------------------------------------
811
812      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
813      IMPLICIT NONE
814#include "iniprint.h"
815
816      integer i, k, klon,klev
817      real pplev(klon,klev+1),pplay(klon,klev)
818      real ztv(klon,klev)
819      real po(klon,klev)
820      real ztva(klon,klev)
821      real zqla(klon,klev)
822      real f_star(klon,klev)
823      real zw2(klon,klev)
824      integer long(klon)
825      real seuil
826      character*21 comment
827
828      if (prt_level.ge.1) THEN
829       print*,'WARNING !!! TEST ',comment
830      endif
831      return
832
833!  test sur la hauteur des thermiques ...
834         do i=1,klon
835!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
836           if (prt_level.ge.10) then
837               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
838               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
839               do k=1,klev
840                  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)
841               enddo
842           endif
843         enddo
844
845
846      return
847      end
848
Note: See TracBrowser for help on using the repository browser.