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

Last change on this file since 1373 was 1372, checked in by musat, 15 years ago

Bug fix: overflow for pplay field
IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.9 KB
Line 
1!
2! $Id: thermcell_main.F90 1372 2010-04-27 10:56:22Z musat $
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!IM   do k=1,nlay
613      do k=1,nlay-1
614         do ig=1,ngrid
615         if ((pcon(ig).le.pplay(ig,k))  &
616     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
617            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
618         endif
619         enddo
620      enddo
621!IM
622      do ig=1,ngrid
623        if (pcon(ig).le.pplay(ig,nlay)) then
624           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
625           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
626           CALL abort_gcm (modname,abort_message,1)
627        endif
628      enddo
629      if (prt_level.ge.1) print*,'14b OK convect8'
630      do k=nlay,1,-1
631         do ig=1,ngrid
632            if (zqla(ig,k).gt.1e-10) then
633               nivcon(ig)=k
634               zcon(ig)=zlev(ig,k)
635            endif
636         enddo
637      enddo
638      if (prt_level.ge.1) print*,'14c OK convect8'
639!calcul des moments
640!initialisation
641      do l=1,nlay
642         do ig=1,ngrid
643            q2(ig,l)=0.
644            wth2(ig,l)=0.
645            wth3(ig,l)=0.
646            ratqscth(ig,l)=0.
647            ratqsdiff(ig,l)=0.
648         enddo
649      enddo     
650      if (prt_level.ge.1) print*,'14d OK convect8'
651      if (prt_level.ge.10)write(lunout,*)                                &
652    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
653      do l=1,nlay
654         do ig=1,ngrid
655            zf=fraca(ig,l)
656            zf2=zf/(1.-zf)
657!
658      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
659!
660      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)
661            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
662            if(zw2(ig,l).gt.1.e-10) then
663             wth2(ig,l)=zf2*(zw2(ig,l))**2
664            else
665             wth2(ig,l)=0.
666            endif
667!           print*,'wth2=',wth2(ig,l)
668            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
669     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
670      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
671            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
672!test: on calcul q2/po=ratqsc
673            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
674         enddo
675      enddo
676!calcul des flux: q, thetal et thetav
677      do l=1,nlay
678         do ig=1,ngrid
679      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
680      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
681      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
682         enddo
683      enddo
684!
685!      print*,'avant calcul ale et alp'
686!calcul de ALE et ALP pour la convection
687      Alp_bl(:)=0.
688      Ale_bl(:)=0.
689!          print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
690      do l=1,nlay
691      do ig=1,ngrid
692           Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
693           Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2)
694!          print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
695      enddo
696      enddo
697
698!     print*,'AAAAAAA ',Alp_bl,Ale_bl,lmix
699
700
701! TEST. IL FAUT REECRIRE LES ALE et ALP
702!     Ale_bl(:)=0.5*wmax(:)*wmax(:)
703!     Alp_bl(:)=0.1*wmax(:)*wmax(:)*wmax(:)
704
705!test:calcul de la ponderation des couches pour KE
706!initialisations
707!      print*,'ponderation'
708      do ig=1,ngrid
709           fm_tot(ig)=0.
710      enddo
711       do ig=1,ngrid
712        do k=1,klev
713           wght_th(ig,k)=1.
714        enddo
715       enddo
716       do ig=1,ngrid
717!         lalim_conv(ig)=lmix_bis(ig)
718!la hauteur de la couche alim_conv = hauteur couche alim_therm
719         lalim_conv(ig)=lalim(ig)
720!         zentr(ig)=zlev(ig,lalim(ig))
721      enddo
722      do ig=1,ngrid
723        do k=1,lalim_conv(ig)
724           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
725        enddo
726      enddo
727      do ig=1,ngrid
728        do k=1,lalim_conv(ig)
729           if (fm_tot(ig).gt.1.e-10) then
730!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
731           endif
732!on pondere chaque couche par a*
733             if (alim_star(ig,k).gt.1.e-10) then
734             wght_th(ig,k)=alim_star(ig,k)
735             else
736             wght_th(ig,k)=1.
737             endif
738        enddo
739      enddo
740!      print*,'apres wght_th'
741!test pour prolonger la convection
742      do ig=1,ngrid
743!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
744      if ((alim_star(ig,1).lt.1.e-10)) then
745      lalim_conv(ig)=1
746      wght_th(ig,1)=1.
747!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
748      endif
749      enddo
750
751!calcul du ratqscdiff
752      if (prt_level.ge.1) print*,'14e OK convect8'
753      var=0.
754      vardiff=0.
755      ratqsdiff(:,:)=0.
756      do ig=1,ngrid
757         do l=1,lalim(ig)
758            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
759         enddo
760      enddo
761      if (prt_level.ge.1) print*,'14f OK convect8'
762      do ig=1,ngrid
763          do l=1,lalim(ig)
764          zf=fraca(ig,l)
765          zf2=zf/(1.-zf)
766          vardiff=vardiff+alim_star(ig,l)  &
767     &           *(zqta(ig,l)*1000.-var)**2
768!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
769!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
770          enddo
771      enddo
772      if (prt_level.ge.1) print*,'14g OK convect8'
773      do l=1,nlay
774         do ig=1,ngrid
775            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
776!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
777         enddo
778      enddo
779!--------------------------------------------------------------------   
780!
781!ecriture des fichiers sortie
782!     print*,'15 OK convect8'
783
784#ifdef wrgrads_thermcell
785      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
786#include "thermcell_out3d.h"
787#endif
788
789      endif
790
791      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
792
793      return
794      end
795
796!-----------------------------------------------------------------------------
797
798      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
799      IMPLICIT NONE
800#include "iniprint.h"
801
802      integer i, k, klon,klev
803      real pplev(klon,klev+1),pplay(klon,klev)
804      real ztv(klon,klev)
805      real po(klon,klev)
806      real ztva(klon,klev)
807      real zqla(klon,klev)
808      real f_star(klon,klev)
809      real zw2(klon,klev)
810      integer long(klon)
811      real seuil
812      character*21 comment
813
814      if (prt_level.ge.1) THEN
815       print*,'WARNING !!! TEST ',comment
816      endif
817      return
818
819!  test sur la hauteur des thermiques ...
820         do i=1,klon
821!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
822           if (prt_level.ge.10) then
823               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
824               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
825               do k=1,klev
826                  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)
827               enddo
828           endif
829         enddo
830
831
832      return
833      end
834
Note: See TracBrowser for help on using the repository browser.