source: LMDZ4/trunk/libf/phytherm/thermcell_main.F90 @ 834

Last change on this file since 834 was 814, checked in by Laurent Fairhead, 17 years ago

Rajout de la physique utilisant les thermiques FH
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.8 KB
Line 
1      SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
2     &                  ,pplay,pplev,pphi,debut  &
3     &                  ,pu,pv,pt,po  &
4     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
5     &                  ,fm0,entr0,zqla,lmax  &
6     &                  ,ratqscth,ratqsdiff,zqsatth  &
7     &                  ,r_aspect,l_mix,w2di,tho)
8
9      IMPLICIT NONE
10
11!=======================================================================
12!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
13!   Version du 09.02.07
14!   Calcul du transport vertical dans la couche limite en presence
15!   de "thermiques" explicitement representes avec processus nuageux
16!
17!   Réécriture à partir d'un listing papier à Habas, le 14/02/00
18!
19!   le thermique est supposé homogène et dissipé par mélange avec
20!   son environnement. la longueur l_mix contrôle l'efficacité du
21!   mélange
22!
23!   Le calcul du transport des différentes espèces se fait en prenant
24!   en compte:
25!     1. un flux de masse montant
26!     2. un flux de masse descendant
27!     3. un entrainement
28!     4. un detrainement
29!
30!=======================================================================
31
32!-----------------------------------------------------------------------
33!   declarations:
34!   -------------
35
36#include "dimensions.h"
37#include "dimphy.h"
38#include "YOMCST.h"
39#include "YOETHF.h"
40#include "FCTTRE.h"
41
42!   arguments:
43!   ----------
44
45      INTEGER ngrid,nlay,w2di,tho
46      real ptimestep,l_mix,r_aspect
47      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
48      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
49      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
50      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
51      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
52      real pphi(ngrid,nlay)
53
54!   local:
55!   ------
56
57      integer,save :: igout=871
58      integer,save :: lunout=6
59      integer,save :: lev_out=0
60
61      INTEGER ig,k,l,ll
62      real zsortie1d(klon)
63      INTEGER lmax(klon),lmin(klon),lalim(klon)
64      INTEGER lmix(klon)
65      real linter(klon)
66      real zmix(klon)
67      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev)
68      real zmax_sec(klon)
69      real w_est(klon,klev+1)
70!on garde le zmax du pas de temps precedent
71      real zmax0(klon)
72      save zmax0
73
74      real zlev(klon,klev+1),zlay(klon,klev)
75      real deltaz(klon,klev)
76      REAL zh(klon,klev),zdhadj(klon,klev)
77      real zthl(klon,klev),zdthladj(klon,klev)
78      REAL ztv(klon,klev)
79      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
80      real zl(klon,klev)
81      real zsortie(klon,klev)
82      real zva(klon,klev)
83      real zua(klon,klev)
84      real zoa(klon,klev)
85
86      real zta(klon,klev)
87      real zha(klon,klev)
88      real fraca(klon,klev+1)
89      real zf,zf2
90      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
91      real q2(klon,klev)
92      common/comtherm/thetath2,wth2
93   
94      real ratqscth(klon,klev)
95      real var
96      real vardiff
97      real ratqsdiff(klon,klev)
98      integer isplit,nsplit
99      parameter (nsplit=10)
100      data isplit/0/
101      save isplit
102
103      logical sorties
104      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
105      real zpspsk(klon,klev)
106
107      real wmax(klon)
108      real wmax_sec(klon)
109      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
110      real detr0(klon,klev)
111      real fm(klon,klev+1),entr(klon,klev)
112
113      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
114!niveau de condensation
115      real nivcon(klon)
116      real zcon(klon)
117      REAL CHI
118      real zcon2(klon)
119      real pcon(klon)
120      real zqsat(klon,klev)
121      real zqsatth(klon,klev)
122
123      real f_star(klon,klev+1),entr_star(klon,klev)
124      real detr_star(klon,klev)
125      real alim_star_tot(klon),alim_star2(klon)
126      real alim_star(klon,klev)
127      real f(klon), f0(klon)
128      save f0
129      real zlevinter(klon)
130      logical debut
131
132!
133
134      character*2 str2
135      character*10 str10
136
137      EXTERNAL SCOPY
138!
139
140!-----------------------------------------------------------------------
141!   initialisation:
142!   ---------------
143!
144      if (lev_out.ge.1) print*,'thermcell_main V4'
145
146       sorties=.true.
147      IF(ngrid.NE.klon) THEN
148         PRINT*
149         PRINT*,'STOP dans convadj'
150         PRINT*,'ngrid    =',ngrid
151         PRINT*,'klon  =',klon
152      ENDIF
153!
154!Initialisation
155!
156      do ig=1,klon     
157      if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
158            f0(ig)=1.e-5
159            zmax0(ig)=40.
160      endif
161      enddo
162
163
164
165!-----------------------------------------------------------------------
166! Calcul de T,q,ql a partir de Tl et qT dans l environnement
167!   --------------------------------------------------------------------
168!
169      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
170     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
171       
172      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_env'
173
174!------------------------------------------------------------------------
175!                       --------------------
176!
177!
178!                       + + + + + + + + + + +
179!
180!
181!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
182!  wh,wt,wo ...
183!
184!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
185!
186!
187!                       --------------------   zlev(1)
188!                       \\\\\\\\\\\\\\\\\\\\
189!
190!
191
192!-----------------------------------------------------------------------
193!   Calcul des altitudes des couches
194!-----------------------------------------------------------------------
195
196      do l=2,nlay
197         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
198      enddo
199         zlev(:,1)=0.
200         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
201      do l=1,nlay
202         zlay(:,l)=pphi(:,l)/RG
203      enddo
204!calcul de l epaisseur des couches
205      do l=1,nlay
206         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
207      enddo
208
209!     print*,'2 OK convect8'
210!-----------------------------------------------------------------------
211!   Calcul des densites
212!-----------------------------------------------------------------------
213
214      do l=1,nlay
215         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
216      enddo
217
218      do l=2,nlay
219         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
220      enddo
221
222!calcul de la masse
223      do l=1,nlay
224         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
225      enddo
226
227      if (lev_out.ge.1) print*,'thermcell_main apres initialisation'
228
229!------------------------------------------------------------------
230!   Calcul de w2, carre de w a partir de la cape
231!
232!   Indicages:
233!   l'ascendance provenant du niveau k traverse l'interface l avec
234!   une vitesse wa(k,l).
235!
236!                       --------------------
237!
238!                       + + + + + + + + + +
239!
240!  wa(k,l)   ----       --------------------    l
241!             /\
242!            /||\       + + + + + + + + + +
243!             ||
244!             ||        --------------------
245!             ||
246!             ||        + + + + + + + + + +
247!             ||
248!             ||        --------------------
249!             ||__
250!             |___      + + + + + + + + + +     k
251!
252!                       --------------------
253!
254!
255!
256!------------------------------------------------------------------
257!definition du profil d alimentation a partir de la flottabilite:
258!alim_star, alim_star_tot, lalim et lmin
259!------------------------------------------------------------------
260!
261      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
262      CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
263     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
264
265      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init'
266      if (lev_out.ge.10) then
267         write(lunout,*) 'Dans thermcell_main 1'
268         write(lunout,*) 'lmin ',lmin(igout)
269         write(lunout,*) 'lalim ',lalim(igout)
270      endif
271
272!-----------------------------------------------------------------------------------
273!calcul des caracteristiques du thermique sec pour le calcul de detr et la fermeture
274!-----------------------------------------------------------------------------------
275!
276      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
277     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
278
279      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_dry'
280
281
282
283!---------------------------------------------------------------------------------
284!calcul du melange et des variables dans le thermique
285!--------------------------------------------------------------------------------
286!
287      CALL thermcell_plume(ngrid,nlay,ztv,zthl,po,zl,rhobarz,  &
288     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
289     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
290     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
291
292      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume'
293      if (lev_out.ge.10) then
294         write(lunout,*) 'Dans thermcell_main 2'
295         write(lunout,*) 'lmin ',lmin(igout)
296         write(lunout,*) 'lalim ',lalim(igout)
297         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
298         write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
299     &    ,f_star(igout,l+1),l=1,lmax(igout))
300      endif
301
302!-------------------------------------------------------------------------------
303! Calcul des caracteristiques du thermique:zmax,zmix,wmax
304!-------------------------------------------------------------------------------
305!
306      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
307     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
308
309      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height'
310
311!-------------------------------------------------------------------------------
312! Fermeture,determination de f
313!-------------------------------------------------------------------------------
314
315      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
316     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
317
318      if(lev_out.ge.1)print*,'thermcell_closure apres thermcell_closure'
319
320!-------------------------------------------------------------------------------
321!deduction des flux
322!-------------------------------------------------------------------------------
323
324      CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
325     &       lalim,lmax,alim_star,  &
326     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
327     &       detr,zqla,zmax,lev_out,lunout,igout)
328
329      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_flux'
330
331!c------------------------------------------------------------------
332!   calcul du transport vertical
333!------------------------------------------------------------------
334
335      if (w2di.eq.1) then
336         fm0=fm0+ptimestep*(fm-fm0)/float(tho)
337         entr0=entr0+ptimestep*(entr-entr0)/float(tho)
338      else
339         fm0=fm
340         entr0=entr
341         detr0=detr
342      endif
343
344      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
345     &                    zthl,zdthladj,zta,lev_out)
346      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
347     &                   po,pdoadj,zoa,lev_out)
348
349      if (1.eq.0) then
350
351! Calcul du transport de V tenant compte d'echange par gradient
352! de pression horizontal avec l'environnement
353
354         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
355     &    ,fraca,zmax  &
356     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
357      else
358
359! calcul purement conservatif pour le transport de V
360         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
361     &    ,zu,pduadj,zua,lev_out)
362         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
363     &    ,zv,pdvadj,zva,lev_out)
364      endif
365
366!     print*,'13 OK convect8'
367      do l=1,nlay
368         do ig=1,ngrid
369           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
370         enddo
371      enddo
372
373      print*,'14 OK convect8'
374!------------------------------------------------------------------
375!   Calculs de diagnostiques pour les sorties
376!------------------------------------------------------------------
377!calcul de fraca pour les sorties
378     
379      if (sorties) then
380      do ig=1,klon
381         fraca(ig,1)=0.
382      enddo
383      do l=2,nlay
384         do ig=1,klon
385            if (zw2(ig,l).gt.1.e-10) then
386            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
387            else
388            fraca(ig,l)=0.
389            endif
390         enddo
391      enddo
392     
393! calcul du niveau de condensation
394! initialisation
395      do ig=1,ngrid
396         nivcon(ig)=0.
397         zcon(ig)=0.
398      enddo
399!nouveau calcul
400      do ig=1,ngrid
401      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
402      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
403      enddo
404      do k=1,nlay
405         do ig=1,ngrid
406         if ((pcon(ig).le.pplay(ig,k))  &
407     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
408            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
409         endif
410         enddo
411      enddo
412      do k=nlay,1,-1
413         do ig=1,ngrid
414            if (zqla(ig,k).gt.1e-10) then
415               nivcon(ig)=k
416               zcon(ig)=zlev(ig,k)
417            endif
418         enddo
419      enddo
420!calcul des moments
421!initialisation
422      do l=1,nlay
423         do ig=1,ngrid
424            q2(ig,l)=0.
425            wth2(ig,l)=0.
426            wth3(ig,l)=0.
427            ratqscth(ig,l)=0.
428            ratqsdiff(ig,l)=0.
429         enddo
430      enddo     
431      do l=1,nlay
432         do ig=1,ngrid
433            zf=fraca(ig,l)
434            zf2=zf/(1.-zf)
435            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
436            wth2(ig,l)=zf2*(zw2(ig,l))**2
437!           print*,'wth2=',wth2(ig,l)
438            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
439     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
440            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
441!test: on calcul q2/po=ratqsc
442            ratqscth(ig,l)=sqrt(q2(ig,l))/(po(ig,l)*1000.)
443         enddo
444      enddo
445!calcul du ratqscdiff
446      var=0.
447      vardiff=0.
448      ratqsdiff(:,:)=0.
449      do ig=1,ngrid
450         do l=1,lalim(ig)
451            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
452         enddo
453      enddo
454      do ig=1,ngrid
455          do l=1,lalim(ig)
456          zf=fraca(ig,l)
457          zf2=zf/(1.-zf)
458          vardiff=vardiff+alim_star(ig,l)  &
459     &           *(zqta(ig,l)*1000.-var)**2
460!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
461!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
462          enddo
463      enddo
464      do l=1,nlay
465         do ig=1,ngrid
466            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
467!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
468         enddo
469      enddo
470!--------------------------------------------------------------------   
471!
472!ecriture des fichiers sortie
473!     print*,'15 OK convect8'
474
475      isplit=isplit+1       
476
477
478#ifdef und
479      if (lev_out.ge.1) print*,'thermcell_main sorties 1D'
480#include "thermcell_out1d.h"
481#endif
482
483
484! #define troisD
485      if (lev_out.ge.1) print*,'thermcell_main sorties 3D'
486#ifdef troisD
487#include "thermcell_out3d.h"
488#endif
489
490      endif
491
492      if (lev_out.ge.1) print*,'thermcell_main FIN  OK'
493
494      return
495      end
496
497!-----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.