source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 @ 1299

Last change on this file since 1299 was 1299, checked in by Laurent Fairhead, 14 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.6 KB
Line 
1!
2! $Id: thermcell_plume.F90 1299 2010-01-20 14:27:21Z fairhead $
3!
4      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
5     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
6     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
7     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8     &           ,lev_out,lunout1,igout)
9
10!--------------------------------------------------------------------------
11!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12!--------------------------------------------------------------------------
13
14      IMPLICIT NONE
15
16#include "YOMCST.h"
17#include "YOETHF.h"
18#include "FCTTRE.h"
19#include "iniprint.h"
20#include "thermcell.h"
21
22      INTEGER itap
23      INTEGER lunout1,igout
24      INTEGER ngrid,klev
25      REAL ptimestep
26      REAL ztv(ngrid,klev)
27      REAL zthl(ngrid,klev)
28      REAL po(ngrid,klev)
29      REAL zl(ngrid,klev)
30      REAL rhobarz(ngrid,klev)
31      REAL zlev(ngrid,klev+1)
32      REAL pplev(ngrid,klev+1)
33      REAL pphi(ngrid,klev)
34      REAL zpspsk(ngrid,klev)
35      REAL alim_star(ngrid,klev)
36      REAL f0(ngrid)
37      INTEGER lalim(ngrid)
38      integer lev_out                           ! niveau pour les print
39      real zcon2(ngrid)
40   
41      real alim_star_tot(ngrid)
42
43      REAL ztva(ngrid,klev)
44      REAL ztla(ngrid,klev)
45      REAL zqla(ngrid,klev)
46      REAL zqta(ngrid,klev)
47      REAL zha(ngrid,klev)
48
49      REAL detr_star(ngrid,klev)
50      REAL coefc
51      REAL entr_star(ngrid,klev)
52      REAL detr(ngrid,klev)
53      REAL entr(ngrid,klev)
54
55      REAL csc(ngrid,klev)
56
57      REAL zw2(ngrid,klev+1)
58      REAL w_est(ngrid,klev+1)
59      REAL f_star(ngrid,klev+1)
60      REAL wa_moy(ngrid,klev+1)
61
62      REAL ztva_est(ngrid,klev)
63      REAL zqla_est(ngrid,klev)
64      REAL zqsatth(ngrid,klev)
65      REAL zta_est(ngrid,klev)
66      REAL zdw2
67      REAL zw2modif
68      REAL zeps
69
70      REAL linter(ngrid)
71      INTEGER lmix(ngrid)
72      INTEGER lmix_bis(ngrid)
73      REAL    wmaxa(ngrid)
74
75      INTEGER ig,l,k
76
77      real zdz,zfact,zbuoy,zalpha
78      real zcor,zdelta,zcvm5,qlbef
79      real Tbef,qsatbef
80      real dqsat_dT,DT,num,denom
81      REAL REPS,RLvCp,DDT0
82      PARAMETER (DDT0=.01)
83      logical Zsat
84      LOGICAL active(ngrid),activetmp(ngrid)
85      REAL fact_gamma,fact_epsilon,fact_gamma2
86      REAL c2(ngrid,klev)
87
88      REAL zw2fact,expa
89      Zsat=.false.
90! Initialisation
91      RLvCp = RLVTT/RCPD
92     
93      if (iflag_thermals_ed==0) then
94         fact_gamma=1.
95         fact_epsilon=1.
96      else if (iflag_thermals_ed==1)  then
97! Valeurs au moment de la premiere soumission des papiers
98         fact_gamma=1.
99         fact_epsilon=0.002
100         fact_gamma2=0.6
101! Suggestions des Fleurs, Septembre 2009
102         fact_epsilon=0.015
103!test cr
104!         fact_epsilon=0.002
105         fact_gamma=0.9
106         fact_gamma2=0.7
107
108      else if (iflag_thermals_ed==2)  then
109         fact_gamma=1.
110         fact_epsilon=2.
111      else if (iflag_thermals_ed==3)  then
112         fact_gamma=3./4.
113         fact_epsilon=3.
114      endif
115
116      write(lunout,*)'THERM 31H '
117
118! Initialisations des variables reeles
119if (1==0) then
120      ztva(:,:)=ztv(:,:)
121      ztva_est(:,:)=ztva(:,:)
122      ztla(:,:)=zthl(:,:)
123      zqta(:,:)=po(:,:)
124      zha(:,:) = ztva(:,:)
125else
126      ztva(:,:)=0.
127      ztva_est(:,:)=0.
128      ztla(:,:)=0.
129      zqta(:,:)=0.
130      zha(:,:) =0.
131endif
132
133      zqla_est(:,:)=0.
134      zqsatth(:,:)=0.
135      zqla(:,:)=0.
136      detr_star(:,:)=0.
137      entr_star(:,:)=0.
138      alim_star(:,:)=0.
139      alim_star_tot(:)=0.
140      csc(:,:)=0.
141      detr(:,:)=0.
142      entr(:,:)=0.
143      zw2(:,:)=0.
144      w_est(:,:)=0.
145      f_star(:,:)=0.
146      wa_moy(:,:)=0.
147      linter(:)=1.
148      linter(:)=1.
149
150! Initialisation des variables entieres
151      lmix(:)=1
152      lmix_bis(:)=2
153      wmaxa(:)=0.
154      lalim(:)=1
155
156!-------------------------------------------------------------------------
157! On ne considere comme actif que les colonnes dont les deux premieres
158! couches sont instables.
159!-------------------------------------------------------------------------
160      active(:)=ztv(:,1)>ztv(:,2)
161
162!-------------------------------------------------------------------------
163! Definition de l'alimentation a l'origine dans thermcell_init
164!-------------------------------------------------------------------------
165      do l=1,klev-1
166         do ig=1,ngrid
167            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
168               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
169     &                       *sqrt(zlev(ig,l+1))
170               lalim(:)=l+1
171               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
172            endif
173         enddo
174      enddo
175      do l=1,klev
176         do ig=1,ngrid
177            if (alim_star_tot(ig) > 1.e-10 ) then
178               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
179            endif
180         enddo
181      enddo
182      alim_star_tot(:)=1.
183
184
185!------------------------------------------------------------------------------
186! Calcul dans la premiere couche
187! On decide dans cette version que le thermique n'est actif que si la premiere
188! couche est instable.
189! Pourrait etre change si on veut que le thermiques puisse se déclencher
190! dans une couche l>1
191!------------------------------------------------------------------------------
192do ig=1,ngrid
193! Le panache va prendre au debut les caracteristiques de l'air contenu
194! dans cette couche.
195    if (active(ig)) then
196    ztla(ig,1)=zthl(ig,1)
197    zqta(ig,1)=po(ig,1)
198    zqla(ig,1)=zl(ig,1)
199!cr: attention, prise en compte de f*(1)=1
200    f_star(ig,2)=alim_star(ig,1)
201    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
202&                     *(zlev(ig,2)-zlev(ig,1))  &
203&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
204    w_est(ig,2)=zw2(ig,2)
205    endif
206enddo
207!
208
209!==============================================================================
210!boucle de calcul de la vitesse verticale dans le thermique
211!==============================================================================
212do l=2,klev-1
213!==============================================================================
214
215
216! On decide si le thermique est encore actif ou non
217! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
218    do ig=1,ngrid
219       active(ig)=active(ig) &
220&                 .and. zw2(ig,l)>1.e-10 &
221&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
222    enddo
223
224
225
226! Premier calcul de la vitesse verticale a partir de la temperature
227! potentielle virtuelle
228!     if (1.eq.1) then
229!         w_est(ig,3)=zw2(ig,2)* &
230!    &      ((f_star(ig,2))**2) &
231!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
232!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
233!    &      *(zlev(ig,3)-zlev(ig,2))
234!      endif
235
236
237!---------------------------------------------------------------------------
238! calcul des proprietes thermodynamiques et de la vitesse de la couche l
239! sans tenir compte du detrainement et de l'entrainement dans cette
240! couche
241! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
242! avant) a l'alimentation pour avoir un calcul plus propre
243!---------------------------------------------------------------------------
244
245     call thermcell_condens(ngrid,active, &
246&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
247
248    do ig=1,ngrid
249        if(active(ig)) then
250        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
251        zta_est(ig,l)=ztva_est(ig,l)
252        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
253        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
254     &      -zqla_est(ig,l))-zqla_est(ig,l))
255
256             w_est(ig,l+1)=zw2(ig,l)*  &
257     &                   ((f_star(ig,l))**2)  &
258     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
259     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
260     &                   *(zlev(ig,l+1)-zlev(ig,l))
261             if (w_est(ig,l+1).lt.0.) then
262                w_est(ig,l+1)=zw2(ig,l)
263             endif
264       endif
265    enddo
266
267!-------------------------------------------------
268!calcul des taux d'entrainement et de detrainement
269!-------------------------------------------------
270
271     do ig=1,ngrid
272        if (active(ig)) then
273          zdz=zlev(ig,l+1)-zlev(ig,l)
274          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
275          zfact=fact_gamma/(1.+fact_gamma)
276 
277! estimation de la fraction couverte par les thermiques
278          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
279
280!calcul de la soumission papier
281          if (1.eq.1) then
282             fact_epsilon=0.0007
283!             zfact=0.9/(1.+0.9)
284             zfact=0.3
285             fact_gamma=0.7
286             fact_gamma2=0.6
287             expa=0.25
288! Calcul  du taux d'entrainement entr_star (epsilon)
289           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
290     &     zbuoy/w_est(ig,l+1) )&
291!- fact_epsilon/zalpha**0.25  ) &
292     &     +0.000 )
293
294!           entr_star(ig,l)=f_star(ig,l)*zdz * (  1./3 * MAX(0.,  &     
295!     &     zbuoy/w_est(ig,l+1) - 1./zalpha**0.25  ) &
296!     &     +0.000 )
297! Calcul du taux de detrainement detr_star (delta)
298!           if (zqla_est(ig,l).lt.1.e-10) then
299!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
300!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
301!     &     +0.0006 )
302!           else
303!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
304!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
305!     &     +0.002 )
306!           endif
307           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
308     &     fact_gamma2 * MAX(0., &
309!fact_epsilon/zalpha**0.25
310     &      -zbuoy/w_est(ig,l+1) )       &
311!     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
312!test
313     &     +  0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &   
314     &     +0.0000 )
315          else
316
317! Calcul  du taux d'entrainement entr_star (epsilon)
318           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
319     &     zbuoy/w_est(ig,l+1) - fact_epsilon  ) &
320     &     +0.0000 )
321
322! Calcul du taux de detrainement detr_star (delta)
323           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
324     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
325     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
326     &     +0.0000 )
327
328          endif
329!endif choix du calcul de E* et D*
330
331!cr test
332!           entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)     
333
334! Prise en compte de la fraction
335!      detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
336
337! En dessous de lalim, on prend le max de alim_star et entr_star pour
338! alim_star et 0 sinon
339        if (l.lt.lalim(ig)) then
340          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
341          entr_star(ig,l)=0.
342        endif
343
344! Calcul du flux montant normalise
345      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
346     &              -detr_star(ig,l)
347
348      endif
349   enddo
350
351!----------------------------------------------------------------------------
352!calcul de la vitesse verticale en melangeant Tl et qt du thermique
353!---------------------------------------------------------------------------
354   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
355   do ig=1,ngrid
356       if (activetmp(ig)) then
357           Zsat=.false.
358           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
359     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
360     &            /(f_star(ig,l+1)+detr_star(ig,l))
361           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
362     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
363     &            /(f_star(ig,l+1)+detr_star(ig,l))
364
365        endif
366    enddo
367
368   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
369
370
371   do ig=1,ngrid
372      if (activetmp(ig)) then
373        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
374! on ecrit de maniere conservative (sat ou non)
375!          T = Tl +Lv/Cp ql
376           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
377           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
378!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
379           zha(ig,l) = ztva(ig,l)
380           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
381     &              -zqla(ig,l))-zqla(ig,l))
382
383!on ecrit zqsat
384           zqsatth(ig,l)=qsatbef 
385
386!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
387!          zw2(ig,l+1)=&
388!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
389!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
390!     &                 *1./(1.+fact_gamma) &
391!     &                 *(zlev(ig,l+1)-zlev(ig,l))
392!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
393! La meme en plus modulaire :
394           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
395           zdz=zlev(ig,l+1)-zlev(ig,l)
396
397
398           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
399
400if (1==0) then
401           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
402           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
403           zw2(ig,l+1)=zw2modif+zdw2
404else
405! Tentative de reecriture de l'equation de w2. A reprendre ...
406!           zdw2=2*zdz*zbuoy
407!           zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon))
408!!!!!      zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
409!formulation Arnaud
410!           zw2fact=zbuoy*zalpha**expa/fact_epsilon
411!           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) &
412!      &    +zw2fact
413           if (zbuoy.gt.1.e-10) then
414           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact)
415           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
416      &    +zw2fact
417           else
418           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma)
419           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
420      &    +zw2fact
421
422           endif
423
424endif
425!           zw2(ig,l+1)=zw2modif+zdw2
426      endif
427   enddo
428
429        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
430!
431!---------------------------------------------------------------------------
432!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
433!---------------------------------------------------------------------------
434
435   do ig=1,ngrid
436            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
437!               stop'On tombe sur le cas particulier de thermcell_dry'
438                write(lunout,*)                                         &
439     &          'On tombe sur le cas particulier de thermcell_plume'
440                zw2(ig,l+1)=0.
441                linter(ig)=l+1
442            endif
443
444        if (zw2(ig,l+1).lt.0.) then
445           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
446     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
447           zw2(ig,l+1)=0.
448        endif
449
450           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
451
452        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
453!   lmix est le niveau de la couche ou w (wa_moy) est maximum
454!on rajoute le calcul de lmix_bis
455            if (zqla(ig,l).lt.1.e-10) then
456               lmix_bis(ig)=l+1
457            endif
458            lmix(ig)=l+1
459            wmaxa(ig)=wa_moy(ig,l+1)
460        endif
461   enddo
462
463!=========================================================================
464! FIN DE LA BOUCLE VERTICALE
465      enddo
466!=========================================================================
467
468!on recalcule alim_star_tot
469       do ig=1,ngrid
470          alim_star_tot(ig)=0.
471       enddo
472       do ig=1,ngrid
473          do l=1,lalim(ig)-1
474          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
475          enddo
476       enddo
477       
478
479        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
480
481!     print*,'thermcell_plume OK'
482
483      return
484      end
Note: See TracBrowser for help on using the repository browser.