source: trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F @ 2176

Last change on this file since 2176 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

  • Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
  • Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
  • Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

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