source: LMDZ4/trunk/libf/phylmd/thermcell_dry.F90 @ 923

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

Bascule de la physique du LMD vers la physique avec thermiques
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.0 KB
Line 
1       SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
2     &                            lalim,lmin,zmax,wmax,lev_out)
3
4!--------------------------------------------------------------------------
5!thermcell_dry: calcul de zmax et wmax du thermique sec
6!--------------------------------------------------------------------------
7       IMPLICIT NONE
8#include "YOMCST.h"       
9       INTEGER l,ig
10
11       INTEGER ngrid,nlay
12       REAL zlev(ngrid,nlay+1)
13       REAL pphi(ngrid,nlay)
14       REAl ztv(ngrid,nlay)
15       REAL alim_star(ngrid,nlay)
16       INTEGER lalim(ngrid)
17      integer lev_out                           ! niveau pour les print
18
19       REAL zmax(ngrid)
20       REAL wmax(ngrid)
21
22!variables locales
23       REAL zw2(ngrid,nlay+1)
24       REAL f_star(ngrid,nlay+1)
25       REAL ztva(ngrid,nlay+1)
26       REAL wmaxa(ngrid)
27       REAL wa_moy(ngrid,nlay+1)
28       REAL linter(ngrid),zlevinter(ngrid)
29       INTEGER lmix(ngrid),lmax(ngrid),lmin(ngrid)
30
31!initialisations
32       do ig=1,ngrid
33          do l=1,nlay+1
34             zw2(ig,l)=0.
35             wa_moy(ig,l)=0.
36          enddo
37       enddo
38       do ig=1,ngrid
39          do l=1,nlay
40             ztva(ig,l)=ztv(ig,l)
41          enddo
42       enddo
43       do ig=1,ngrid
44          wmax(ig)=0.
45          wmaxa(ig)=0.
46       enddo
47!calcul de la vitesse a partir de la CAPE en melangeant thetav
48
49!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
50! A eliminer
51! Ce if complique etait fait pour reperer la premiere couche instable
52! Ici, c'est lmin.
53!
54!       do l=1,nlay-2
55!         do ig=1,ngrid
56!            if (ztv(ig,l).gt.ztv(ig,l+1)  &
57!     &         .and.alim_star(ig,l).gt.1.e-10  &
58!     &         .and.zw2(ig,l).lt.1e-10) then
59!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
60
61
62! Calcul des F^*, integrale verticale de E^*
63       f_star(:,1)=0.
64       do l=1,nlay
65          f_star(:,l+1)=f_star(:,l)+alim_star(:,l)
66       enddo
67
68! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
69       linter(:)=0.
70
71! couche la plus haute concernee par le thermique.
72       lmax(:)=1
73
74! Le niveau linter est une variable continue qui se trouve dans la couche
75! lmax
76
77       do l=1,nlay-2
78         do ig=1,ngrid
79            if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
80
81!------------------------------------------------------------------------
82!  Calcul de la vitesse en haut de la premiere couche instable.
83!  Premiere couche du panache thermique
84!------------------------------------------------------------------------
85               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
86     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
87     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
88
89!------------------------------------------------------------------------
90! Tant que la vitesse en bas de la couche et la somme du flux de masse
91! et de l'entrainement (c'est a dire le flux de masse en haut) sont
92! positifs, on calcul
93! 1. le flux de masse en haut  f_star(ig,l+1)
94! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
95! 3. la vitesse au carré en haut zw2(ig,l+1)
96!------------------------------------------------------------------------
97
98!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
99!  A eliminer : dans cette version, si zw2 est > 0 on a un therique.
100!  et donc, au dessus, f_star(ig,l+1) est forcement suffisamment
101!  grand puisque on n'a pas de detrainement.
102!  f_star est une fonction croissante.
103!  c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
104!           else if ((zw2(ig,l).ge.1e-10).and.  &
105!    &               (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
106!              f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
107!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
108
109            else if (zw2(ig,l).ge.1e-10) then
110
111               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l)  &
112     &                    *ztv(ig,l))/f_star(ig,l+1)
113               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+  &
114     &                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
115     &                     *(zlev(ig,l+1)-zlev(ig,l))
116            endif
117! determination de zmax continu par interpolation lineaire
118!------------------------------------------------------------------------
119
120            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
121!               stop'On tombe sur le cas particulier de thermcell_dry'
122                print*,'On tombe sur le cas particulier de thermcell_dry'
123                zw2(ig,l+1)=0.
124                linter(ig)=l+1
125                lmax(ig)=l
126            endif
127
128            if (zw2(ig,l+1).lt.0.) then
129               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
130     &           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
131               zw2(ig,l+1)=0.
132               lmax(ig)=l
133            endif
134
135               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
136
137            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
138!   lmix est le niveau de la couche ou w (wa_moy) est maximum
139               lmix(ig)=l+1
140               wmaxa(ig)=wa_moy(ig,l+1)
141            endif
142         enddo
143      enddo
144       if (lev_out.ge.1) print*,'fin calcul zw2'
145!
146!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
147! A eliminer :
148! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
149! Calcul de la couche correspondant a la hauteur du thermique
150!      do ig=1,ngrid
151!         lmax(ig)=lalim(ig)
152!      enddo
153!      do ig=1,ngrid
154!         do l=nlay,lalim(ig)+1,-1
155!            if (zw2(ig,l).le.1.e-10) then
156!               lmax(ig)=l-1
157!            endif
158!         enddo
159!      enddo
160!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
161
162!   
163! Determination de zw2 max
164      do ig=1,ngrid
165         wmax(ig)=0.
166      enddo
167
168      do l=1,nlay
169         do ig=1,ngrid
170            if (l.le.lmax(ig)) then
171                zw2(ig,l)=sqrt(zw2(ig,l))
172                wmax(ig)=max(wmax(ig),zw2(ig,l))
173            else
174                 zw2(ig,l)=0.
175            endif
176          enddo
177      enddo
178
179!   Longueur caracteristique correspondant a la hauteur des thermiques.
180      do  ig=1,ngrid
181         zmax(ig)=0.
182         zlevinter(ig)=zlev(ig,1)
183      enddo
184      do  ig=1,ngrid
185! calcul de zlevinter
186
187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188! FH A eliminer
189! Simplification
190!          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
191!     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
192!     &    -zlev(ig,lmax(ig)))
193!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194
195          zlevinter(ig)=zlev(ig,lmax(ig)) + &
196     &    (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
197           zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
198      enddo
199
200! Verification que lalim<=lmax
201      do ig=1,ngrid
202         if(lalim(ig)>lmax(ig)) then
203            print*,'WARNING thermcell_dry ig=',ig,'  lalim=',lalim(ig),'  lmax(ig)=',lmax(ig)
204            lmax(ig)=lalim(ig)
205         endif
206      enddo
207     
208      RETURN
209      END
Note: See TracBrowser for help on using the repository browser.