source: trunk/LMDZ.EARTH/libf/phylmd/thermcell_main.F90 @ 357

Last change on this file since 357 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5

File size: 27.0 KB
Line 
1!
2! $Id: thermcell_main.F90 1403 2010-07-01 09:02:53Z fairhead $
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,ztv &
13     &                  ,zpspsk,ztla,zthl)
14
15      USE dimphy
16      USE comgeomphy , ONLY:rlond,rlatd
17      IMPLICIT NONE
18
19!=======================================================================
20!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
21!   Version du 09.02.07
22!   Calcul du transport vertical dans la couche limite en presence
23!   de "thermiques" explicitement representes avec processus nuageux
24!
25!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
26!
27!   le thermique est suppose homogene et dissipe par melange avec
28!   son environnement. la longueur l_mix controle l'efficacite du
29!   melange
30!
31!   Le calcul du transport des differentes especes se fait en prenant
32!   en compte:
33!     1. un flux de masse montant
34!     2. un flux de masse descendant
35!     3. un entrainement
36!     4. un detrainement
37!
38!=======================================================================
39
40!-----------------------------------------------------------------------
41!   declarations:
42!   -------------
43
44#include "dimensions.h"
45#include "YOMCST.h"
46#include "YOETHF.h"
47#include "FCTTRE.h"
48#include "iniprint.h"
49
50!   arguments:
51!   ----------
52
53!IM 140508
54      INTEGER itap
55
56      INTEGER ngrid,nlay,w2di
57      real tau_thermals
58      integer iflag_thermals_ed
59      real ptimestep,l_mix,r_aspect
60      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
61      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
62      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
63      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
64      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
65      real pphi(ngrid,nlay)
66
67!   local:
68!   ------
69
70      integer icount
71      data icount/0/
72      save icount
73!$OMP THREADPRIVATE(icount)
74
75      integer,save :: igout=1
76!$OMP THREADPRIVATE(igout)
77      integer,save :: lunout1=6
78!$OMP THREADPRIVATE(lunout1)
79      integer,save :: lev_out=10
80!$OMP THREADPRIVATE(lev_out)
81
82      INTEGER ig,k,l,ll
83      real zsortie1d(klon)
84      INTEGER lmax(klon),lmin(klon),lalim(klon)
85      INTEGER lmix(klon)
86      INTEGER lmix_bis(klon)
87      real linter(klon)
88      real zmix(klon)
89      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
90!      real fraca(klon,klev)
91
92      real zmax_sec(klon)
93!on garde le zmax du pas de temps precedent
94      real zmax0(klon)
95!FH/IM     save zmax0
96
97      real lambda
98
99      real zlev(klon,klev+1),zlay(klon,klev)
100      real deltaz(klon,klev)
101      REAL zh(klon,klev)
102      real zthl(klon,klev),zdthladj(klon,klev)
103      REAL ztv(klon,klev)
104      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
105      real zl(klon,klev)
106      real zsortie(klon,klev)
107      real zva(klon,klev)
108      real zua(klon,klev)
109      real zoa(klon,klev)
110
111      real zta(klon,klev)
112      real zha(klon,klev)
113      real fraca(klon,klev+1)
114      real zf,zf2
115      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
116      real q2(klon,klev)
117! FH probleme de dimensionnement avec l'allocation dynamique
118!     common/comtherm/thetath2,wth2
119      real wq(klon,klev)
120      real wthl(klon,klev)
121      real wthv(klon,klev)
122   
123      real ratqscth(klon,klev)
124      real var
125      real vardiff
126      real ratqsdiff(klon,klev)
127
128      logical sorties
129      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
130      real zpspsk(klon,klev)
131
132      real wmax(klon)
133      real wmax_tmp(klon)
134      real wmax_sec(klon)
135      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
136      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
137
138      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
139!niveau de condensation
140      integer nivcon(klon)
141      real zcon(klon)
142      REAL CHI
143      real zcon2(klon)
144      real pcon(klon)
145      real zqsat(klon,klev)
146      real zqsatth(klon,klev)
147
148      real f_star(klon,klev+1),entr_star(klon,klev)
149      real detr_star(klon,klev)
150      real alim_star_tot(klon)
151      real alim_star(klon,klev)
152      real alim_star_clos(klon,klev)
153      real f(klon), f0(klon)
154!FH/IM     save f0
155      real zlevinter(klon)
156      logical debut
157       real seuil
158      real csc(klon,klev)
159
160!
161      !nouvelles variables pour la convection
162      real Ale_bl(klon)
163      real Alp_bl(klon)
164      real alp_int(klon)
165      real ale_int(klon)
166      integer n_int(klon)
167      real fm_tot(klon)
168      real wght_th(klon,klev)
169      integer lalim_conv(klon)
170!v1d     logical therm
171!v1d     save therm
172
173      character*2 str2
174      character*10 str10
175
176      character (len=20) :: modname='thermcell_main'
177      character (len=80) :: abort_message
178
179      EXTERNAL SCOPY
180!
181
182!-----------------------------------------------------------------------
183!   initialisation:
184!   ---------------
185!
186
187      seuil=0.25
188
189      if (debut)  then
190         fm0=0.
191         entr0=0.
192         detr0=0.
193
194
195#undef wrgrads_thermcell
196#ifdef wrgrads_thermcell
197! Initialisation des sorties grads pour les thermiques.
198! Pour l'instant en 1D sur le point igout.
199! Utilise par thermcell_out3d.h
200         str10='therm'
201         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
202     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
203     &   ptimestep,str10,'therm ')
204#endif
205
206
207
208      endif
209
210      fm=0. ; entr=0. ; detr=0.
211
212
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 NOUVELLE/NOUVELLE Arnaud'
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>9) 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
471! Choix de la fonction d'alimentation utilisee pour la fermeture.
472! Apparemment sans importance
473      alim_star_clos(:,:)=alim_star(:,:)
474      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
475
476! Appel avec la version seche
477      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
478     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
479
480!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481! Appel avec les zmax et wmax tenant compte de la condensation
482! Semble moins bien marcher
483!     CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
484!    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
485!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
486
487      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
488
489      if (tau_thermals>1.) then
490         lambda=exp(-ptimestep/tau_thermals)
491         f0=(1.-lambda)*f+lambda*f0
492      else
493         f0=f
494      endif
495
496! Test valable seulement en 1D mais pas genant
497      if (.not. (f0(1).ge.0.) ) then
498              abort_message = '.not. (f0(1).ge.0.)'
499              CALL abort_gcm (modname,abort_message,1)
500      endif
501
502!-------------------------------------------------------------------------------
503!deduction des flux
504!-------------------------------------------------------------------------------
505
506      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
507     &       lalim,lmax,alim_star,  &
508     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
509     &       detr,zqla,lev_out,lunout1,igout)
510!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
511
512      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
513      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
514      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
515
516!------------------------------------------------------------------
517!   On ne prend pas directement les profils issus des calculs precedents
518!   mais on s'autorise genereusement une relaxation vers ceci avec
519!   une constante de temps tau_thermals (typiquement 1800s).
520!------------------------------------------------------------------
521
522      if (tau_thermals>1.) then
523         lambda=exp(-ptimestep/tau_thermals)
524         fm0=(1.-lambda)*fm+lambda*fm0
525         entr0=(1.-lambda)*entr+lambda*entr0
526         detr0=(1.-lambda)*detr+lambda*detr0
527      else
528         fm0=fm
529         entr0=entr
530         detr0=detr
531      endif
532
533!c------------------------------------------------------------------
534!   calcul du transport vertical
535!------------------------------------------------------------------
536
537      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
538     &                    zthl,zdthladj,zta,lev_out)
539      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
540     &                   po,pdoadj,zoa,lev_out)
541
542!------------------------------------------------------------------
543! Calcul de la fraction de l'ascendance
544!------------------------------------------------------------------
545      do ig=1,klon
546         fraca(ig,1)=0.
547         fraca(ig,nlay+1)=0.
548      enddo
549      do l=2,nlay
550         do ig=1,klon
551            if (zw2(ig,l).gt.1.e-10) then
552            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
553            else
554            fraca(ig,l)=0.
555            endif
556         enddo
557      enddo
558     
559!------------------------------------------------------------------
560!  calcul du transport vertical du moment horizontal
561!------------------------------------------------------------------
562
563!IM 090508 
564      if (1.eq.1) then
565!IM 070508 vers. _dq       
566!     if (1.eq.0) then
567
568
569! Calcul du transport de V tenant compte d'echange par gradient
570! de pression horizontal avec l'environnement
571
572         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
573     &    ,fraca,zmax  &
574     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
575
576      else
577
578! calcul purement conservatif pour le transport de V
579         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
580     &    ,zu,pduadj,zua,lev_out)
581         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
582     &    ,zv,pdvadj,zva,lev_out)
583      endif
584
585!     print*,'13 OK convect8'
586      do l=1,nlay
587         do ig=1,ngrid
588           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
589         enddo
590      enddo
591
592      if (prt_level.ge.1) print*,'14 OK convect8'
593!------------------------------------------------------------------
594!   Calculs de diagnostiques pour les sorties
595!------------------------------------------------------------------
596!calcul de fraca pour les sorties
597     
598      if (sorties) then
599      if (prt_level.ge.1) print*,'14a OK convect8'
600! calcul du niveau de condensation
601! initialisation
602      do ig=1,ngrid
603         nivcon(ig)=0
604         zcon(ig)=0.
605      enddo
606!nouveau calcul
607      do ig=1,ngrid
608      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
609      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
610      enddo
611!IM   do k=1,nlay
612      do k=1,nlay-1
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!IM
621      do ig=1,ngrid
622        if (pcon(ig).le.pplay(ig,nlay)) then
623           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
624           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
625           CALL abort_gcm (modname,abort_message,1)
626        endif
627      enddo
628      if (prt_level.ge.1) print*,'14b OK convect8'
629      do k=nlay,1,-1
630         do ig=1,ngrid
631            if (zqla(ig,k).gt.1e-10) then
632               nivcon(ig)=k
633               zcon(ig)=zlev(ig,k)
634            endif
635         enddo
636      enddo
637      if (prt_level.ge.1) print*,'14c OK convect8'
638!calcul des moments
639!initialisation
640      do l=1,nlay
641         do ig=1,ngrid
642            q2(ig,l)=0.
643            wth2(ig,l)=0.
644            wth3(ig,l)=0.
645            ratqscth(ig,l)=0.
646            ratqsdiff(ig,l)=0.
647         enddo
648      enddo     
649      if (prt_level.ge.1) print*,'14d OK convect8'
650      if (prt_level.ge.10)write(lunout,*)                                &
651    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
652      do l=1,nlay
653         do ig=1,ngrid
654            zf=fraca(ig,l)
655            zf2=zf/(1.-zf)
656!
657      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
658!
659      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)
660            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
661            if(zw2(ig,l).gt.1.e-10) then
662             wth2(ig,l)=zf2*(zw2(ig,l))**2
663            else
664             wth2(ig,l)=0.
665            endif
666!           print*,'wth2=',wth2(ig,l)
667            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
668     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
669      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
670            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
671!test: on calcul q2/po=ratqsc
672            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
673         enddo
674      enddo
675!calcul des flux: q, thetal et thetav
676      do l=1,nlay
677         do ig=1,ngrid
678      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
679      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
680      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
681         enddo
682      enddo
683!
684!      print*,'avant calcul ale et alp'
685!calcul de ALE et ALP pour la convection
686      Alp_bl(:)=0.
687      Ale_bl(:)=0.
688!          print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
689      do l=1,nlay
690      do ig=1,ngrid
691           Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
692           Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2)
693!          print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
694      enddo
695      enddo
696
697!     print*,'AAAAAAA ',Alp_bl,Ale_bl,lmix
698
699
700! TEST. IL FAUT REECRIRE LES ALE et ALP
701!     Ale_bl(:)=0.5*wmax(:)*wmax(:)
702!     Alp_bl(:)=0.1*wmax(:)*wmax(:)*wmax(:)
703
704!test:calcul de la ponderation des couches pour KE
705!initialisations
706!      print*,'ponderation'
707      do ig=1,ngrid
708           fm_tot(ig)=0.
709      enddo
710       do ig=1,ngrid
711        do k=1,klev
712           wght_th(ig,k)=1.
713        enddo
714       enddo
715       do ig=1,ngrid
716!         lalim_conv(ig)=lmix_bis(ig)
717!la hauteur de la couche alim_conv = hauteur couche alim_therm
718         lalim_conv(ig)=lalim(ig)
719!         zentr(ig)=zlev(ig,lalim(ig))
720      enddo
721      do ig=1,ngrid
722        do k=1,lalim_conv(ig)
723           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
724        enddo
725      enddo
726      do ig=1,ngrid
727        do k=1,lalim_conv(ig)
728           if (fm_tot(ig).gt.1.e-10) then
729!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
730           endif
731!on pondere chaque couche par a*
732             if (alim_star(ig,k).gt.1.e-10) then
733             wght_th(ig,k)=alim_star(ig,k)
734             else
735             wght_th(ig,k)=1.
736             endif
737        enddo
738      enddo
739!      print*,'apres wght_th'
740!test pour prolonger la convection
741      do ig=1,ngrid
742!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
743      if ((alim_star(ig,1).lt.1.e-10)) then
744      lalim_conv(ig)=1
745      wght_th(ig,1)=1.
746!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
747      endif
748      enddo
749
750!calcul du ratqscdiff
751      if (prt_level.ge.1) print*,'14e OK convect8'
752      var=0.
753      vardiff=0.
754      ratqsdiff(:,:)=0.
755      do ig=1,ngrid
756         do l=1,lalim(ig)
757            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
758         enddo
759      enddo
760      if (prt_level.ge.1) print*,'14f OK convect8'
761      do ig=1,ngrid
762          do l=1,lalim(ig)
763          zf=fraca(ig,l)
764          zf2=zf/(1.-zf)
765          vardiff=vardiff+alim_star(ig,l)  &
766     &           *(zqta(ig,l)*1000.-var)**2
767!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
768!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
769          enddo
770      enddo
771      if (prt_level.ge.1) print*,'14g OK convect8'
772      do l=1,nlay
773         do ig=1,ngrid
774            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
775!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
776         enddo
777      enddo
778!--------------------------------------------------------------------   
779!
780!ecriture des fichiers sortie
781!     print*,'15 OK convect8'
782
783#ifdef wrgrads_thermcell
784      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
785#include "thermcell_out3d.h"
786#endif
787
788      endif
789
790      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
791
792      return
793      end
794
795!-----------------------------------------------------------------------------
796
797      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
798      IMPLICIT NONE
799#include "iniprint.h"
800
801      integer i, k, klon,klev
802      real pplev(klon,klev+1),pplay(klon,klev)
803      real ztv(klon,klev)
804      real po(klon,klev)
805      real ztva(klon,klev)
806      real zqla(klon,klev)
807      real f_star(klon,klev)
808      real zw2(klon,klev)
809      integer long(klon)
810      real seuil
811      character*21 comment
812
813      if (prt_level.ge.1) THEN
814       print*,'WARNING !!! TEST ',comment
815      endif
816      return
817
818!  test sur la hauteur des thermiques ...
819         do i=1,klon
820!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
821           if (prt_level.ge.10) then
822               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
823               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
824               do k=1,klev
825                  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)
826               enddo
827           endif
828         enddo
829
830
831      return
832      end
833
Note: See TracBrowser for help on using the repository browser.