source: LMDZ4/trunk/libf/phylmd/thermcell_main.F90 @ 1080

Last change on this file since 1080 was 1026, checked in by lmdzadmin, 16 years ago

Modifs thermiques
FH/IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.7 KB
Line 
1!
2! $Header$
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 &
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!
153      !nouvelles variables pour la convection
154      real Ale_bl(klon)
155      real Alp_bl(klon)
156      real alp_int(klon)
157      real ale_int(klon)
158      integer n_int(klon)
159      real fm_tot(klon)
160      real wght_th(klon,klev)
161      integer lalim_conv(klon)
162!v1d     logical therm
163!v1d     save therm
164
165      character*2 str2
166      character*10 str10
167
168      EXTERNAL SCOPY
169!
170
171!-----------------------------------------------------------------------
172!   initialisation:
173!   ---------------
174!
175
176      seuil=0.25
177
178      if (debut)  then
179         fm0=0.
180         entr0=0.
181         detr0=0.
182
183
184! #define wrgrads_thermcell
185#undef wrgrads_thermcell
186#ifdef wrgrads_thermcell
187! Initialisation des sorties grads pour les thermiques.
188! Pour l'instant en 1D sur le point igout.
189! Utilise par thermcell_out3d.h
190         str10='therm'
191         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
192     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
193     &   ptimestep,str10,'therm ')
194#endif
195
196
197
198      endif
199
200      fm=0. ; entr=0. ; detr=0.
201
202      icount=icount+1
203
204!IM 090508 beg
205!print*,'====================================================================='
206!print*,'====================================================================='
207!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
208!print*,'====================================================================='
209!print*,'====================================================================='
210!IM 090508 end
211
212      if (prt_level.ge.1) print*,'thermcell_main V4'
213
214       sorties=.true.
215      IF(ngrid.NE.klon) THEN
216         PRINT*
217         PRINT*,'STOP dans convadj'
218         PRINT*,'ngrid    =',ngrid
219         PRINT*,'klon  =',klon
220      ENDIF
221!
222!Initialisation
223!
224!    IF (1.eq.0) THEN
225!     do ig=1,klon     
226!FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
227!     if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
228!           f0(ig)=1.e-5
229!           zmax0(ig)=40.
230!v1d        therm=.false.
231!     endif
232!     enddo
233!    ENDIF !(1.eq.0) THEN
234     print*,'WARNING thermcell_main f0=max(f0,1.e-2)'
235     do ig=1,klon
236      if (prt_level.ge.20) then
237       print*,'th_main ig f0',ig,f0(ig)
238      endif
239         f0(ig)=max(f0(ig),1.e-2)
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      print*,'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
298      rhobarz(:,1)=rho(:,1)
299
300      do l=2,nlay
301         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
302      enddo
303
304!calcul de la masse
305      do l=1,nlay
306         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
307      enddo
308
309      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
310
311!------------------------------------------------------------------
312!
313!             /|\
314!    --------  |  F_k+1 -------   
315!                              ----> D_k
316!             /|\              <---- E_k , A_k
317!    --------  |  F_k ---------
318!                              ----> D_k-1
319!                              <---- E_k-1 , A_k-1
320!
321!
322!
323!
324!
325!    ---------------------------
326!
327!    ----- F_lmax+1=0 ----------         \
328!            lmax     (zmax)              |
329!    ---------------------------          |
330!                                         |
331!    ---------------------------          |
332!                                         |
333!    ---------------------------          |
334!                                         |
335!    ---------------------------          |
336!                                         |
337!    ---------------------------          |
338!                                         |  E
339!    ---------------------------          |  D
340!                                         |
341!    ---------------------------          |
342!                                         |
343!    ---------------------------  \       |
344!            lalim                 |      |
345!    ---------------------------   |      |
346!                                  |      |
347!    ---------------------------   |      |
348!                                  | A    |
349!    ---------------------------   |      |
350!                                  |      |
351!    ---------------------------   |      |
352!    lmin  (=1 pour le moment)     |      |
353!    ----- F_lmin=0 ------------  /      /
354!
355!    ---------------------------
356!    //////////////////////////
357!
358!
359!=============================================================================
360!  Calculs initiaux ne faisant pas intervenir les changements de phase
361!=============================================================================
362
363!------------------------------------------------------------------
364!  1. alim_star est le profil vertical de l'alimentation à la base du
365!     panache thermique, calculé à partir de la flotabilité de l'air sec
366!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
367!------------------------------------------------------------------
368!
369      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
370      CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
371     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
372
373call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
374call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
375
376
377      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
378      if (prt_level.ge.10) then
379         write(lunout1,*) 'Dans thermcell_main 1'
380         write(lunout1,*) 'lmin ',lmin(igout)
381         write(lunout1,*) 'lalim ',lalim(igout)
382         write(lunout1,*) ' ig l alim_star thetav'
383         write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
384     &   ,ztv(igout,l),l=1,lalim(igout)+4)
385      endif
386
387!v1d      do ig=1,klon
388!v1d     if (alim_star(ig,1).gt.1.e-10) then
389!v1d     therm=.true.
390!v1d     endif
391!v1d      enddo
392!-----------------------------------------------------------------------------
393!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
394!     panache sec conservatif (e=d=0) alimente selon alim_star
395!     Il s'agit d'un calcul de type CAPE
396!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
397!------------------------------------------------------------------------------
398!
399      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
400     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
401
402call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
403call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
404
405      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
406      if (prt_level.ge.10) then
407         write(lunout1,*) 'Dans thermcell_main 1b'
408         write(lunout1,*) 'lmin ',lmin(igout)
409         write(lunout1,*) 'lalim ',lalim(igout)
410         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
411         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
412     &    ,l=1,lalim(igout)+4)
413      endif
414
415
416
417!---------------------------------------------------------------------------------
418!calcul du melange et des variables dans le thermique
419!--------------------------------------------------------------------------------
420!
421      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
422!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
423      CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
424     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
425     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
426     &           ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
427     &            ,lev_out,lunout1,igout)
428      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
429
430      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
431      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
432
433      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
434      if (prt_level.ge.10) then
435         write(lunout1,*) 'Dans thermcell_main 2'
436         write(lunout1,*) 'lmin ',lmin(igout)
437         write(lunout1,*) 'lalim ',lalim(igout)
438         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
439         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
440     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
441      endif
442
443!-------------------------------------------------------------------------------
444! Calcul des caracteristiques du thermique:zmax,zmix,wmax
445!-------------------------------------------------------------------------------
446!
447      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
448     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
449
450
451      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
452      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
453      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
454      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
455
456      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
457
458!-------------------------------------------------------------------------------
459! Fermeture,determination de f
460!-------------------------------------------------------------------------------
461!
462!avant closure: on redéfinit lalim, alim_star_tot et alim_star
463!       do ig=1,klon
464!       do l=2,lalim(ig)
465!       alim_star(ig,l)=entr_star(ig,l)
466!       entr_star(ig,l)=0.
467!       enddo
468!       enddo
469
470      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
471     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
472
473      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
474
475      if (tau_thermals>1.) then
476         lambda=exp(-ptimestep/tau_thermals)
477         f0=(1.-lambda)*f+lambda*f0
478      else
479         f0=f
480      endif
481
482! Test valable seulement en 1D mais pas genant
483      if (.not. (f0(1).ge.0.) ) then
484           stop'Dans thermcell_main'
485      endif
486
487!-------------------------------------------------------------------------------
488!deduction des flux
489!-------------------------------------------------------------------------------
490
491      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
492     &       lalim,lmax,alim_star,  &
493     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
494     &       detr,zqla,lev_out,lunout1,igout)
495!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
496
497      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
498      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
499      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
500
501!------------------------------------------------------------------
502!   On ne prend pas directement les profils issus des calculs precedents
503!   mais on s'autorise genereusement une relaxation vers ceci avec
504!   une constante de temps tau_thermals (typiquement 1800s).
505!------------------------------------------------------------------
506
507      if (tau_thermals>1.) then
508         lambda=exp(-ptimestep/tau_thermals)
509         fm0=(1.-lambda)*fm+lambda*fm0
510         entr0=(1.-lambda)*entr+lambda*entr0
511!        detr0=(1.-lambda)*detr+lambda*detr0
512      else
513         fm0=fm
514         entr0=entr
515         detr0=detr
516      endif
517
518!c------------------------------------------------------------------
519!   calcul du transport vertical
520!------------------------------------------------------------------
521
522      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
523     &                    zthl,zdthladj,zta,lev_out)
524      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
525     &                   po,pdoadj,zoa,lev_out)
526
527!------------------------------------------------------------------
528! Calcul de la fraction de l'ascendance
529!------------------------------------------------------------------
530      do ig=1,klon
531         fraca(ig,1)=0.
532         fraca(ig,nlay+1)=0.
533      enddo
534      do l=2,nlay
535         do ig=1,klon
536            if (zw2(ig,l).gt.1.e-10) then
537            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
538            else
539            fraca(ig,l)=0.
540            endif
541         enddo
542      enddo
543     
544!------------------------------------------------------------------
545!  calcul du transport vertical du moment horizontal
546!------------------------------------------------------------------
547
548!IM 090508 
549      if (1.eq.1) then
550!IM 070508 vers. _dq       
551!     if (1.eq.0) then
552
553
554! Calcul du transport de V tenant compte d'echange par gradient
555! de pression horizontal avec l'environnement
556
557         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
558     &    ,fraca,zmax  &
559     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
560!IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
561      else
562
563! calcul purement conservatif pour le transport de V
564         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
565     &    ,zu,pduadj,zua,lev_out)
566         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
567     &    ,zv,pdvadj,zva,lev_out)
568      endif
569
570!     print*,'13 OK convect8'
571      do l=1,nlay
572         do ig=1,ngrid
573           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
574         enddo
575      enddo
576
577      if (prt_level.ge.1) print*,'14 OK convect8'
578!------------------------------------------------------------------
579!   Calculs de diagnostiques pour les sorties
580!------------------------------------------------------------------
581!calcul de fraca pour les sorties
582     
583      if (sorties) then
584      if (prt_level.ge.1) print*,'14a OK convect8'
585! calcul du niveau de condensation
586! initialisation
587      do ig=1,ngrid
588         nivcon(ig)=0
589         zcon(ig)=0.
590      enddo
591!nouveau calcul
592      do ig=1,ngrid
593      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
594      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
595      enddo
596      do k=1,nlay
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      if (prt_level.ge.1) print*,'14b OK convect8'
605      do k=nlay,1,-1
606         do ig=1,ngrid
607            if (zqla(ig,k).gt.1e-10) then
608               nivcon(ig)=k
609               zcon(ig)=zlev(ig,k)
610            endif
611         enddo
612      enddo
613      if (prt_level.ge.1) print*,'14c OK convect8'
614!calcul des moments
615!initialisation
616      do l=1,nlay
617         do ig=1,ngrid
618            q2(ig,l)=0.
619            wth2(ig,l)=0.
620            wth3(ig,l)=0.
621            ratqscth(ig,l)=0.
622            ratqsdiff(ig,l)=0.
623         enddo
624      enddo     
625      if (prt_level.ge.1) print*,'14d OK convect8'
626      print*,'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
627      do l=1,nlay
628         do ig=1,ngrid
629            zf=fraca(ig,l)
630            zf2=zf/(1.-zf)
631!
632      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
633!
634      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)
635            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
636            if(zw2(ig,l).gt.1.e-10) then
637             wth2(ig,l)=zf2*(zw2(ig,l))**2
638            else
639             wth2(ig,l)=0.
640            endif
641!           print*,'wth2=',wth2(ig,l)
642            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
643     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
644      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
645            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
646!test: on calcul q2/po=ratqsc
647            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
648         enddo
649      enddo
650!calcul de ale_bl et alp_bl
651!pour le calcul d'une valeur intégrée entre la surface et lmax
652      do ig=1,ngrid
653      alp_int(ig)=0.
654      ale_int(ig)=0.
655      n_int(ig)=0
656      enddo
657!
658      do l=1,nlay
659      do ig=1,ngrid
660       if(l.LE.lmax(ig)) THEN
661        alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
662        ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
663        n_int(ig)=n_int(ig)+1
664       endif
665      enddo
666      enddo
667!      print*,'avant calcul ale et alp'
668!calcul de ALE et ALP pour la convection
669      do ig=1,ngrid
670!      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
671!          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
672!      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig))
673!     &           *0.1
674!valeur integree de alp_bl * 0.5:
675       if (n_int(ig).gt.0) then
676       Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
677!       if (Alp_bl(ig).lt.0.) then
678!       Alp_bl(ig)=0.
679       endif
680!       endif
681!         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
682!     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
683!            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
684!      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
685!      if (nivcon(ig).eq.1) then
686!       Ale_bl(ig)=0.
687!       else
688!valeur max de ale_bl:
689       Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2
690!     & /2.
691!     & *0.1
692!        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
693!       if (n_int(ig).gt.0) then
694!       Ale_bl(ig)=ale_int(ig)/n_int(ig)
695!        Ale_bl(ig)=4.
696!       endif
697!       endif
698!            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
699!          Ale_bl(ig)=wth2(ig,nivcon(ig))
700!          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
701      enddo
702!test:calcul de la ponderation des couches pour KE
703!initialisations
704!      print*,'ponderation'
705      do ig=1,ngrid
706           fm_tot(ig)=0.
707      enddo
708       do ig=1,ngrid
709        do k=1,klev
710           wght_th(ig,k)=1.
711        enddo
712       enddo
713       do ig=1,ngrid
714!         lalim_conv(ig)=lmix_bis(ig)
715!la hauteur de la couche alim_conv = hauteur couche alim_therm
716         lalim_conv(ig)=lalim(ig)
717!         zentr(ig)=zlev(ig,lalim(ig))
718      enddo
719      do ig=1,ngrid
720        do k=1,lalim_conv(ig)
721           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
722        enddo
723      enddo
724      do ig=1,ngrid
725        do k=1,lalim_conv(ig)
726           if (fm_tot(ig).gt.1.e-10) then
727!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
728           endif
729!on pondere chaque couche par a*
730             if (alim_star(ig,k).gt.1.e-10) then
731             wght_th(ig,k)=alim_star(ig,k)
732             else
733             wght_th(ig,k)=1.
734             endif
735        enddo
736      enddo
737!      print*,'apres wght_th'
738!test pour prolonger la convection
739      do ig=1,ngrid
740!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
741      if ((alim_star(ig,1).lt.1.e-10)) then
742      lalim_conv(ig)=1
743      wght_th(ig,1)=1.
744!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
745      endif
746      enddo
747
748!calcul du ratqscdiff
749      if (prt_level.ge.1) print*,'14e OK convect8'
750      var=0.
751      vardiff=0.
752      ratqsdiff(:,:)=0.
753      do ig=1,ngrid
754         do l=1,lalim(ig)
755            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
756         enddo
757      enddo
758      if (prt_level.ge.1) print*,'14f OK convect8'
759      do ig=1,ngrid
760          do l=1,lalim(ig)
761          zf=fraca(ig,l)
762          zf2=zf/(1.-zf)
763          vardiff=vardiff+alim_star(ig,l)  &
764     &           *(zqta(ig,l)*1000.-var)**2
765!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
766!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
767          enddo
768      enddo
769      if (prt_level.ge.1) print*,'14g OK convect8'
770      do l=1,nlay
771         do ig=1,ngrid
772            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
773!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
774         enddo
775      enddo
776!--------------------------------------------------------------------   
777!
778!ecriture des fichiers sortie
779!     print*,'15 OK convect8'
780
781      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
782#ifdef wrgrads_thermcell
783#include "thermcell_out3d.h"
784#endif
785
786      endif
787
788      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
789
790!     if(icount.eq.501) stop'au pas 301 dans thermcell_main'
791      return
792      end
793
794!-----------------------------------------------------------------------------
795
796      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
797      IMPLICIT NONE
798#include "iniprint.h"
799
800      integer i, k, klon,klev
801      real pplev(klon,klev+1),pplay(klon,klev)
802      real ztv(klon,klev)
803      real po(klon,klev)
804      real ztva(klon,klev)
805      real zqla(klon,klev)
806      real f_star(klon,klev)
807      real zw2(klon,klev)
808      integer long(klon)
809      real seuil
810      character*21 comment
811
812      if (prt_level.ge.1) THEN
813       print*,'WARNING !!! TEST ',comment
814      endif
815      return
816
817!  test sur la hauteur des thermiques ...
818         do i=1,klon
819!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
820           if (prt_level.ge.10) then
821               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
822               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
823               do k=1,klev
824                  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)
825               enddo
826!              stop
827           endif
828         enddo
829
830
831      return
832      end
833
Note: See TracBrowser for help on using the repository browser.