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

      MODULE vlz_fi_mod
      
      IMPLICIT NONE
      
      CONTAINS
      
      SUBROUTINE vlz_fi(ngrid,nlay,q,pente_max,masse,w,wq)
c
c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
c
c    ********************************************************************
c    "pseudo upstream" Advection scheme along the vertical
c    to be used in the physics (sedimentation)
c    ********************************************************************

      IMPLICIT NONE
c

c
c
c   Arguments:
c   ----------
      integer,intent(in) :: ngrid ! number of atmospheric columns
      integer,intent(in) :: nlay ! number of atmospheric layers
      real,intent(in) :: masse(ngrid,nlay) ! mass of atmospheric layer delta(P)/g
      real,intent(in) :: pente_max ! maximum slope for the scheme (2 is recommended)
      real,intent(inout) :: q(ngrid,nlay) ! tracer mixing ratio (kg/kg)
      real,intent(inout) :: w(ngrid,nlay) ! mass of atmosphere "transfered" over the time step (kg.m-2)
      real,intent(out) :: wq(ngrid,nlay+1) ! trancer increment due to advection (kg)
c
c      Local 
c   ---------
c
      INTEGER i,ij,l,j,ii
c

      real dzq(ngrid,nlay),dzqw(ngrid,nlay),adzqw(ngrid,nlay),dzqmax
      real newmasse
      real sigw, Mtot, MQtot
      integer m


c    Orientation follows pressure, i.e. follows W

      do l=2,nlay
         do ij=1,ngrid
            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
            adzqw(ij,l)=abs(dzqw(ij,l))
         enddo
      enddo

      do l=2,nlay-1
         do ij=1,ngrid
            if(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) then
                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
            else
                dzq(ij,l)=0.
            endif

            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
         enddo
      enddo

      do ij=1,ngrid
         dzq(ij,1)=0.
         dzq(ij,nlay)=0.
      enddo
c ---------------------------------------------------------------
c   .... compute vertical advection terms   .......
c ---------------------------------------------------------------

c compute  - d( q   * w )/ d(sigma) , later added to dq to compute dq
c
c      No flux at the model top:
       do ij=1,ngrid
          wq(ij,nlay+1)=0.
       enddo

c      1) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)     
c      ===============================

       do l = 1,nlay          ! loop different than when w<0
        do  ij = 1,ngrid
         if(w(ij,l).gt.0.)then

c         Regular scheme (transfered mass < 1 layer)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if(w(ij,l).le.masse(ij,l))then
            sigw=w(ij,l)/masse(ij,l)
            wq(ij,l)=w(ij,l)*(q(ij,l)+0.5*(1.-sigw)*dzq(ij,l))
            

c         Extended scheme (transfered mass > 1 layer)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          else 
            m=l
            Mtot = masse(ij,m)
            MQtot = masse(ij,m)*q(ij,m)
            if(m.ge.nlay)goto 88
            do while(w(ij,l).gt.(Mtot+masse(ij,m+1)))
                m=m+1
                Mtot = Mtot + masse(ij,m)
                MQtot = MQtot + masse(ij,m)*q(ij,m)
                if(m.ge.nlay)goto 88
            end do
 88         continue
            if (m.lt.nlay) then
                sigw=(w(ij,l)-Mtot)/masse(ij,m+1)
                wq(ij,l)=(MQtot + (w(ij,l)-Mtot)*
     &          (q(ij,m+1)+0.5*(1.-sigw)*dzq(ij,m+1)) )
            else
                w(ij,l) = Mtot
                wq(ij,l) = Mqtot 
            end if

          end if
         end if
        enddo
       enddo

c      2) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)     
c      ===============================
       goto 99 ! SKIPPING THIS PART FOR SEDIMENTATION 

c      Surface flux up:
       do  ij = 1,ngrid
         if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid
       end do

       do l = 1,nlay-1  ! loop different than when w>0
        do  ij = 1,ngrid
         if(w(ij,l+1).le.0)then

c         Regular scheme (transfered mass < 1 layer)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if(-w(ij,l+1).le.masse(ij,l))then
            sigw=w(ij,l+1)/masse(ij,l)
            wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
c         Extended scheme (transfered mass > 1 layer)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          else 
             m = l-1
             Mtot = masse(ij,m+1)
             MQtot = masse(ij,m+1)*q(ij,m+1)
             if (m.le.0)goto 77
             do while(-w(ij,l+1).gt.(Mtot+masse(ij,m)))
                m=m-1
                Mtot = Mtot + masse(ij,m+1)
                MQtot = MQtot + masse(ij,m+1)*q(ij,m+1)
                if (m.le.0)goto 77
             end do
 77          continue

             if (m.gt.0) then
                sigw=(w(ij,l+1)+Mtot)/masse(ij,m)
                wq(ij,l+1)= - (MQtot + (-w(ij,l+1)-Mtot)*
     &          (q(ij,m)-0.5*(1.+sigw)*dzq(ij,m))  )
             else
c               wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*qm(ij,1))
                write(*,*) 'a rather weird situation in vlz_fi !'
                call abort_physic("vlz_fi","weird situation",1)
             end if

          endif
         endif
        enddo
       enddo
 99    continue

      do l=1,nlay
         do ij=1,ngrid

cccccccc lines below not used for sedimentation (No real flux)
ccccc       newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l) 
ccccc       q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
ccccc&         /newmasse
ccccc       masse(ij,l)=newmasse

c            it cannot entrain more than available mass !
            if ((wq(ij,l+1)-wq(ij,l)) .lt. -(masse(ij,l)*q(ij,l))) then
              wq(ij,l+1) = wq(ij,l)-masse(ij,l)*q(ij,l)
            end if

            q(ij,l)=q(ij,l) +  (wq(ij,l+1)-wq(ij,l))/masse(ij,l)  
         enddo
      enddo

      END SUBROUTINE vlz_fi

      END MODULE vlz_fi_mod