source: trunk/LMDZ.MARS/libf/phymars/vlz_fi.F @ 3567

Last change on this file since 3567 was 2448, checked in by cmathe, 4 years ago

Delete test_vmr_co2.F, old and unused anymore

File size: 5.8 KB
Line 
1      SUBROUTINE vlz_fi(ngrid,nlay,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
19c
20c
21c   Arguments:
22c   ----------
23      integer,intent(in) :: ngrid ! number of atmospheric columns
24      integer,intent(in) :: nlay ! number of atmospheric layers
25      real masse(ngrid,nlay),pente_max
26      REAL q(ngrid,nlay)
27      REAL w(ngrid,nlay)
28      REAL wq(ngrid,nlay+1)
29c
30c      Local
31c   ---------
32c
33      INTEGER i,ij,l,j,ii
34c
35
36      real dzq(ngrid,nlay),dzqw(ngrid,nlay),adzqw(ngrid,nlay),dzqmax
37      real newmasse
38      real sigw, Mtot, MQtot
39      integer m
40
41      REAL      SSUM,CVMGP,CVMGT
42      integer ismax,ismin
43
44
45c    On oriente tout dans le sens de la pression c'est a dire dans le
46c    sens de W
47
48      do l=2,nlay
49         do ij=1,ngrid
50            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
51            adzqw(ij,l)=abs(dzqw(ij,l))
52         enddo
53      enddo
54
55      do l=2,nlay-1
56         do ij=1,ngrid
57#ifdef CRAY
58            dzq(ij,l)=0.5*
59     ,      cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
60#else
61            if(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) then
62                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
63            else
64                dzq(ij,l)=0.
65            endif
66#endif
67            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
68            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
69         enddo
70      enddo
71
72      do ij=1,ngrid
73         dzq(ij,1)=0.
74         dzq(ij,nlay)=0.
75      enddo
76c ---------------------------------------------------------------
77c   .... calcul des termes d'advection verticale  .......
78c ---------------------------------------------------------------
79
80c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
81c
82c      No flux at the model top:
83       do ij=1,ngrid
84          wq(ij,nlay+1)=0.
85       enddo
86
87c      1) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)     
88c      ===============================
89
90       do l = 1,nlay          ! loop different than when w<0
91        do  ij = 1,ngrid
92         if(w(ij,l).gt.0.)then
93
94c         Regular scheme (transfered mass < 1 layer)
95c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
96          if(w(ij,l).le.masse(ij,l))then
97            sigw=w(ij,l)/masse(ij,l)
98            wq(ij,l)=w(ij,l)*(q(ij,l)+0.5*(1.-sigw)*dzq(ij,l))
99           
100
101c         Extended scheme (transfered mass > 1 layer)
102c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
103          else
104            m=l
105            Mtot = masse(ij,m)
106            MQtot = masse(ij,m)*q(ij,m)
107            if(m.ge.nlay)goto 88
108            do while(w(ij,l).gt.(Mtot+masse(ij,m+1)))
109                m=m+1
110                Mtot = Mtot + masse(ij,m)
111                MQtot = MQtot + masse(ij,m)*q(ij,m)
112                if(m.ge.nlay)goto 88
113            end do
114 88         continue
115            if (m.lt.nlay) then
116                sigw=(w(ij,l)-Mtot)/masse(ij,m+1)
117                wq(ij,l)=(MQtot + (w(ij,l)-Mtot)*
118     &          (q(ij,m+1)+0.5*(1.-sigw)*dzq(ij,m+1)) )
119            else
120                w(ij,l) = Mtot
121                wq(ij,l) = Mqtot
122            end if
123
124          end if
125         end if
126        enddo
127       enddo
128
129c      2) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)     
130c      ===============================
131       goto 99 ! SKIPPING THIS PART FOR SEDIMENTATION
132
133c      Surface flux up:
134       do  ij = 1,ngrid
135         if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid
136       end do
137
138       do l = 1,nlay-1  ! loop different than when w>0
139        do  ij = 1,ngrid
140         if(w(ij,l+1).le.0)then
141
142c         Regular scheme (transfered mass < 1 layer)
143c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144          if(-w(ij,l+1).le.masse(ij,l))then
145            sigw=w(ij,l+1)/masse(ij,l)
146            wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
147c         Extended scheme (transfered mass > 1 layer)
148c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
149          else
150             m = l-1
151             Mtot = masse(ij,m+1)
152             MQtot = masse(ij,m+1)*q(ij,m+1)
153             if (m.le.0)goto 77
154             do while(-w(ij,l+1).gt.(Mtot+masse(ij,m)))
155                m=m-1
156                Mtot = Mtot + masse(ij,m+1)
157                MQtot = MQtot + masse(ij,m+1)*q(ij,m+1)
158                if (m.le.0)goto 77
159             end do
160 77          continue
161
162             if (m.gt.0) then
163                sigw=(w(ij,l+1)+Mtot)/masse(ij,m)
164                wq(ij,l+1)= - (MQtot + (-w(ij,l+1)-Mtot)*
165     &          (q(ij,m)-0.5*(1.+sigw)*dzq(ij,m))  )
166             else
167c               wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*qm(ij,1))
168                write(*,*) 'a rather weird situation in vlz_fi !'
169                call abort_physic("vlz_fi","weird situation",1)
170             end if
171
172          endif
173         endif
174        enddo
175       enddo
176 99    continue
177
178      do l=1,nlay
179         do ij=1,ngrid
180
181cccccccc lines below not used for sedimentation (No real flux)
182ccccc       newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
183ccccc       q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
184ccccc&         /newmasse
185ccccc       masse(ij,l)=newmasse
186
187c            it cannot entrain more than available mass !
188            if ((wq(ij,l+1)-wq(ij,l)) .lt. -(masse(ij,l)*q(ij,l))) then
189              wq(ij,l+1) = wq(ij,l)-masse(ij,l)*q(ij,l)
190            end if
191
192            q(ij,l)=q(ij,l) +  (wq(ij,l+1)-wq(ij,l))/masse(ij,l) 
193         enddo
194      enddo
195
196      return
197      end
Note: See TracBrowser for help on using the repository browser.