source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/vlz_fi.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 6.4 KB
Line 
1      SUBROUTINE vlz_fi(ngrid,q,pente_max,masse,w,wq)
2c
3c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
4c
5c    ********************************************************************
6c     Shema  d'advection " pseudo amont " dans la verticale
7c    pour appel dans la physique (sedimentation)
8c    ********************************************************************
9c    q rapport de melange (kg/kg)...
10c    masse : masse de la couche Dp/g
11c    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
12c    pente_max = 2 conseillee
13c
14c
15c   --------------------------------------------------------------------
16      IMPLICIT NONE
17c
18!-----------------------------------------------------------------------
19!   INCLUDE 'dimensions.h'
20!
21!   dimensions.h contient les dimensions du modele
22!   ndm est tel que iim=2**ndm
23!-----------------------------------------------------------------------
24
25      INTEGER iim,jjm,llm,ndm
26
27      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
28
29!-----------------------------------------------------------------------
30!-----------------------------------------------------------------------
31!   INCLUDE 'dimphys.h'
32
33! ngridmx : number of horizontal grid points
34! note: the -1/jjm term will be 0; unless jj=1
35      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)   
36! nlayermx : number of atmospheric layers
37      integer, parameter :: nlayermx = llm 
38! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
39      !integer, parameter :: nsoilmx = 4 ! for a test
40      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
41      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
42!-----------------------------------------------------------------------
43
44
45c
46c
47c   Arguments:
48c   ----------
49      integer ngrid
50      real masse(ngrid,llm),pente_max
51      REAL q(ngrid,llm)
52      REAL w(ngrid,llm)
53      REAL wq(ngrid,llm+1)
54c
55c      Local 
56c   ---------
57c
58      INTEGER i,ij,l,j,ii
59c
60
61      real dzq(ngrid,llm),dzqw(ngrid,llm),adzqw(ngrid,llm),dzqmax
62      real newmasse
63      real sigw, Mtot, MQtot
64      integer m
65
66      REAL      SSUM,CVMGP,CVMGT
67      integer ismax,ismin
68
69
70c    On oriente tout dans le sens de la pression c'est a dire dans le
71c    sens de W
72
73      do l=2,llm
74         do ij=1,ngrid
75            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
76            adzqw(ij,l)=abs(dzqw(ij,l))
77         enddo
78      enddo
79
80      do l=2,llm-1
81         do ij=1,ngrid
82            if(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) then
83                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
84            else
85                dzq(ij,l)=0.
86            endif
87            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
88            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
89         enddo
90      enddo
91
92      do ij=1,ngrid
93         dzq(ij,1)=0.
94         dzq(ij,llm)=0.
95      enddo
96c ---------------------------------------------------------------
97c   .... calcul des termes d'advection verticale  .......
98c ---------------------------------------------------------------
99
100c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
101c
102c      No flux at the model top:
103       do ij=1,ngrid
104          wq(ij,llm+1)=0.
105       enddo
106
107c      1) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)     
108c      ===============================
109
110       do l = 1,llm          ! loop different than when w<0
111        do ij=1,ngrid
112
113         if(w(ij,l).gt.0.)then
114
115c         Regular scheme (transfered mass < 1 layer)
116c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
117          if(w(ij,l).le.masse(ij,l))then
118            sigw=w(ij,l)/masse(ij,l)
119            wq(ij,l)=w(ij,l)*(q(ij,l)+0.5*(1.-sigw)*dzq(ij,l))
120           
121
122c         Extended scheme (transfered mass > 1 layer)
123c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
124          else
125            m=l
126            Mtot = masse(ij,m)
127            MQtot = masse(ij,m)*q(ij,m)
128            if(m.ge.llm)goto 88
129            do while(w(ij,l).gt.(Mtot+masse(ij,m+1)))
130                m=m+1
131                Mtot = Mtot + masse(ij,m)
132                MQtot = MQtot + masse(ij,m)*q(ij,m)
133                if(m.ge.llm)goto 88
134            end do
135 88         continue
136            if (m.lt.llm) then
137                sigw=(w(ij,l)-Mtot)/masse(ij,m+1)
138                wq(ij,l)=(MQtot + (w(ij,l)-Mtot)*
139     &          (q(ij,m+1)+0.5*(1.-sigw)*dzq(ij,m+1)) )
140            else
141                w(ij,l) = Mtot
142                wq(ij,l) = Mqtot
143            end if
144          end if
145         end if
146        enddo
147       enddo
148
149c      2) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)     
150c      ===============================
151       goto 99 ! SKIPPING THIS PART FOR SEDIMENTATION
152
153c      Surface flux up:
154       do ij=1,ngrid
155         if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid
156       end do
157
158       do l = 1,llm-1  ! loop different than when w>0
159        do ij=1,ngrid
160         if(w(ij,l+1).le.0)then
161
162c         Regular scheme (transfered mass < 1 layer)
163c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
164          if(-w(ij,l+1).le.masse(ij,l))then
165            sigw=w(ij,l+1)/masse(ij,l)
166            wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
167c         Extended scheme (transfered mass > 1 layer)
168c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
169          else
170             m = l-1
171             Mtot = masse(ij,m+1)
172             MQtot = masse(ij,m+1)*q(ij,m+1)
173             if (m.le.0)goto 77
174             do while(-w(ij,l+1).gt.(Mtot+masse(ij,m)))
175                m=m-1
176                Mtot = Mtot + masse(ij,m+1)
177                MQtot = MQtot + masse(ij,m+1)*q(ij,m+1)
178                if (m.le.0)goto 77
179             end do
180 77          continue
181
182             if (m.gt.0) then
183                sigw=(w(ij,l+1)+Mtot)/masse(ij,m)
184                wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*
185     &          (q(ij,m)-0.5*(1.+sigw)*dzq(ij,m))  )
186             else
187c               wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*qm(ij,1))
188                write(*,*) 'a rather weird situation in vlz_fi !'
189                stop
190             end if
191          endif
192         endif
193        enddo
194       enddo
195 99    continue
196
197      do l=1,llm
198         do ij=1,ngrid
199
200cccccccc lines below not used for sedimentation (No real flux)
201ccccc       newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l) 
202ccccc       q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
203ccccc&         /newmasse
204ccccc       masse(ij,l)=newmasse
205
206            q(ij,l)=q(ij,l) +  (wq(ij,l+1)-wq(ij,l))/masse(ij,l)
207
208         enddo
209      enddo
210
211
212
213      return
214      end
Note: See TracBrowser for help on using the repository browser.