source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 @ 1361

Last change on this file since 1361 was 1338, checked in by idelkadi, 15 years ago

Nouvelle version des thermiques

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