#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: Advection schemes --- grid centers\label{sec:advectionMean}
!
! !INTERFACE:
   subroutine adv_center(N,dt,h,ho,ww,Bcup,Bcdw,Yup,Ydw,method,mode,Y)
!
! !DESCRIPTION:
!
! This subroutine solves a one-dimensional advection equation. There are two
! options, depending whether the advection should be conservative or not.
! Conservative advection has to be applied when settling of sediment or
! rising of phytoplankton is considered. In this case the advection is of
! the form
!  \begin{equation}
!   \label{Yadvection_cons}
!    \partder{Y}{t} = - \partder{F}{z}
!    \comma
!  \end{equation}
! where $F=wY$ is the flux caused by the advective velocity, $w$.
!
! Non-conservative advective transport has to be applied, when the water
! has a non-zero vertical velocity. In three-dimensional applications,
! this transport would be conservative, since vertical divergence would be
! compensated by horizontal convergence and vice versa. However, the
! key assumption of one-dimensional modelling is horizontal homogeneity,
! such that we indeed have to apply a vertically non-conservative method,
! which is of the form
!  \begin{equation}
!   \label{Yadvection_noncons}
!    \partder{Y}{t} = - w\partder{Y}{z}
!                   = - \left(\partder{F}{z} - Y\partder{w}{z} \right).
!  \end{equation}
!
! The discretized form of \eq{Yadvection_cons} is
!  \begin{equation}
!   \label{advDiscretized_cons}
!   Y_i^{n+1} = Y_i^n
!   - \dfrac{\Delta t}{h_i}
!    \left( F^n_{i} - F^n_{i-1} \right)
!   \comma
!  \end{equation}
! where the integers $n$ and $i$ correspond to the present time and space
! level, respectively.
!
! For the non-conservative form \eq{Yadvection_noncons},
! an extra term needs to be included:
!  \begin{equation}
!   \label{advDiscretized_noncons}
!   Y_i^{n+1} = Y_i^n
!   - \dfrac{\Delta t}{h_i}
!    \left( F^n_{i} - F^n_{i-1} -Y_i^n \left(w_k-w_{k-1}  \right)\right).
!  \end{equation}
!
!  Which advection method is applied is decided by the flag {\tt mode},
!  which gives conservative advection \eq{advDiscretized_cons}
!  for {\tt mode=1} and
!  non-conservative advection \eq{advDiscretized_noncons} for {\tt mode=0}.
!
! Fluxes are defined at the grid faces, the variable $Y_i$ is defined at the
!  grid centers. The fluxes are computed in an upstream-biased way,
!  \begin{equation}
!   \label{upstream}
!   F^n_{i} = \dfrac{1}{\Delta t}
!   \int_{z^\text{Face}_{i} - w \Delta t}^{z^\text{Face}_{i}} Y(z') dz'
!   \point
!  \end{equation}
!  For a third-order polynomial approximation of $Y$ (see \cite{Pietrzak98}),
!  these fluxes can be written the in so-called Lax-Wendroff form as
!  \begin{equation}
!   \label{fluxDiscretized}
!    \begin{array}{rcll}
!      F_{i} &=& w_{i} \left(Y_i +  \dfrac{1}{2} \Phi^+_{i}
!      \left(1-\magn{c_{i}} \right) \left( Y_{i+1} - Y_i \right) \right)
!      \quad & \text{for} \quad w_{i} > 0
!      \comma  \\[5mm]
!      F_{i} &=& w_{i} \left(Y_{i+1} +  \dfrac{1}{2} \Phi^-_{i}
!      \left(1-\magn{c_{i}} \right) \left( Y_i - Y_{i+1} \right) \right)
!      \quad & \text{for} \quad w_{i} < 0
!      \comma
!    \end{array}
!  \end{equation}
!  where $c_{i} = 2 w_{i} \Delta t / (h_i+h_{i+1})$ is the Courant number.
!  The factors appearing in \eq{fluxDiscretized} are defined as
!  \begin{equation}
!   \label{phiDiscretized}
!  \Phi^+_{i} =  \alpha_{i} +  \beta_{i}  r^+_{i}
!  \comma
!  \Phi^-_{i} =  \alpha_{i} +  \beta_{i}  r^-_{i}
!  \comma
!  \end{equation}
!  where
!  \begin{equation}
!   \label{alphaDiscretized}
!   \alpha_{i} = \dfrac{1}{2}
!    + \dfrac{1}{6} \left( 1- 2 \magn{c_{i}} \right) \comma
!   \beta_{i} = \dfrac{1}{2}
!    - \dfrac{1}{6} \left( 1- 2 \magn{c_{i}} \right)
!  \point
!  \end{equation}
!  The upstream and downstream slope parameters are
!  \begin{equation}
!   \label{slopeDiscretized}
!   r^+_{i} = \dfrac{Y_i - Y_{i-1}}{Y_{i+1}-Y_{i}}  \comma
!   r^-_{i} = \dfrac{Y_{i+2} - Y_{i+1}}{Y_{i+1}-Y_{i}}
!  \point
!  \end{equation}
!
!  To obtain monotonic and positive schemes also in the presence of strong
!  gradients, so-called slope limiters are aplied for the factors $\Phi^+_{i}$
!  and $\Phi^-_{i}$. The two most obvious cases are
!  the first-order upstream discretisation with $\Phi^+_{i}=\Phi^-_{i}=0$
!  and the Lax-Wendroff scheme with  $\Phi^+_{i}=\Phi^-_{i}=1$.
!  The subroutine {\tt adv\_center.F90} provides six different slope-limiters,
!  all discussed in detail by \cite{Pietrzak98}:
!
! \begin{itemize}
!  \item first-order upstream ({\tt method=UPSTREAM})
!  \item second-order upstream-biased polynomial scheme ({\tt method=P1},
!        not yet implemented)
!  \item third-order upstream-biased polynomial scheme ({\tt method=P2})
!  \item third-order scheme (TVD) with Superbee limiter ({\tt method=Superbee})
!  \item third-order scheme (TVD) with MUSCL limiter ({\tt method=MUSCL})
!  \item third-order scheme (TVD) with ULTIMATE QUICKEST limiter
!        ({\tt method=P2\_PDM})
! \end{itemize}
!
!
! If during a certain time step the maximum Courant number is larger
! than one, a split iteration will be carried out which guarantees that the
! split step Courant numbers are just smaller than 1.
!
! Several kinds of boundary conditions are implemented for the upper
! and lower boundaries. They are set by the integer values {\tt Bcup}
! and {\tt Bcdw}, that have to correspond to the parameters defined
! in the module {\tt util}, see \sect{sec:utils}. The
! following choices exist at the moment:
!
! For the value {\tt flux}, the boundary values {\tt Yup} and {\tt Ydw} are
! interpreted as specified fluxes at the uppermost and lowest interface.
! Fluxes into the boundary cells are counted positive by convention.
! For the value {\tt value}, {\tt Yup} and {\tt Ydw} specify the value
! of $Y$ at the interfaces, and the flux is computed by multiplying with
! the (known) speed  at the interface. For the value {\tt oneSided},
! {\tt Yup} and {\tt Ydw} are ignored and the flux is computed
! from a one-sided first-order upstream discretisation using the speed
! at the interface and the value of $Y$ at the center of the boundary cell.
! For the value {\tt zeroDivergence}, the fluxes into and out of the
! respective boundary cell are set equal.
! This corresponds to a zero-gradient formulation, or to zero
! flux divergence in the boundary cells.
!
! Be careful that your boundary conditions are mathematically well defined.
! For example, specifying an inflow into the boundary cell with the
! speed at the boundary being directed outward does not make sense.
!
!
! !USES:
   use util
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  number of vertical layers
   integer,  intent(in)                :: N

!  time step (s)
   REALTYPE, intent(in)                :: dt

!  layer thickness (m)
   REALTYPE, intent(in)                :: h(0:N)

!  old layer thickness (m)
   REALTYPE, intent(in)                :: ho(0:N)

!  vertical advection speed
   REALTYPE, intent(in)                :: ww(0:N)

!  type of upper BC
   integer,  intent(in)                :: Bcup

!  type of lower BC
   integer,  intent(in)                :: Bcdw

!  value of upper BC
   REALTYPE, intent(in)                :: Yup

!  value of lower BC
   REALTYPE, intent(in)                :: Ydw

!  type of advection scheme
   integer,  intent(in)                :: method

!  advection mode (0: non-conservative, 1: conservative)
   integer,  intent(in)                :: mode
!
! !INPUT/OUTPUT PARAMETERS:
   REALTYPE                            :: Y(0:N)
!
! !DEFINED PARAMETERS:
   REALTYPE,     parameter             :: one6th=1.0d0/6.0d0
   integer,      parameter             :: itmax=100
!
! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!  $Log: adv_center.F90,v $
!  Revision 20.0  2013/12/14 00:13:57  fms
!  Merged revision 1.1.2.1 onto trunk
!
!  Revision 1.1.2.1  2012/05/15 16:01:18  smg
!  initial cvs ci to mom5
!  AUTHOR:Griffies
!  REVIEWERS:
!  TEST STATUS:
!  CHANGES PUBLIC INTERFACES?
!  CHANGES ANSWERS?
!
!  Revision 1.1.2.1.390.1  2012/04/23 20:33:41  smg
!  updated to gotm-2012.03.09 CVS tagged code.
!  AUTHOR:Martin Schmidt
!  REVIEWERS:Griffies
!  TEST STATUS:
!  CHANGES PUBLIC INTERFACES?
!  CHANGES ANSWERS?
!
!  Revision 1.5  2010-09-17 12:53:52  jorn
!  extensive code clean-up to ensure proper initialization and clean-up of all variables
!
!  Revision 1.4  2006-11-06 13:36:46  hb
!  Option for conservative vertical advection added to adv_center
!
!  Revision 1.3  2006-03-20 09:06:38  kbk
!  removed explicit double precission dependency
!
!  Revision 1.2  2005/11/18 10:59:34  kbk
!  removed unused variables - some left in parameter lists
!
!  Revision 1.1  2005/06/27 10:54:33  kbk
!  new files needed
!
!
!
!EOP
!
! !LOCAL VARIABLES:
   integer                              :: i,k,it
   REALTYPE                             :: x,r,Phi,limit
   REALTYPE                             :: Yu,Yc,Yd
   REALTYPE                             :: c,cmax
   REALTYPE                             :: cu(0:N)
!
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(*,*) 'adv_center # ',Ncall
#endif

!  initialize interface fluxes with zero
   cu   = _ZERO_

!  initialize maximum Courant number
   cmax = _ZERO_

!  compute maximum Courant number
   do k=1,N-1
      c=abs(ww(k))*dt/(0.5*(h(k)+h(k+1)))
      if (c.gt.cmax) cmax=c
   enddo

   it=min(itmax,int(cmax)+1)

!#ifdef DEBUG
   if (it .gt. 1) then
      STDERR 'In adv_center():'
      STDERR 'Maximum Courant number is ',cmax
      STDERR it,' iterations used for vertical advection'
   endif
!#endif

!  splitting loop
   do i=1,it

!     vertical loop
      do k=1,N-1

!        compute the slope ration
         if (ww(k) .gt. _ZERO_) then

!           compute Courant number
            c=ww(k)/float(it)*dt/(0.5*(h(k)+h(k+1)))

            if (k .gt. 1) then
               Yu=Y(k-1)                              ! upstream value
            else
               Yu=Y(k)
            end if
            Yc=Y(k  )                                 ! central value
            Yd=Y(k+1)                                 ! downstream value

!           compute slope ration
            if (abs(Yd-Yc) .gt. 1e-10) then
               r=(Yc-Yu)/(Yd-Yc)
            else
               r=(Yc-Yu)*1.e10
            end if

!        negative speed
         else

!           compute Courant number
            c=-ww(k)/float(it)*dt/(0.5*(h(k)+h(k+1)))

            if (k .lt. N-1) then
               Yu=Y(k+2)                              ! upstream value
            else
               Yu=Y(k+1)
            end if
            Yc=Y(k+1)                                 ! central value
            Yd=Y(k  )                                 ! downstream value


!           compute slope ratio
            if (abs(Yc-Yd) .gt. 1e-10) then
               r=(Yu-Yc)/(Yc-Yd)
            else
               r=(Yu-Yc)*1.e10
            end if

         end if

!        compute the flux-factor phi
         x    =  one6th*(1.-2.0*c)
         Phi  =  (0.5+x)+(0.5-x)*r

!        limit the flux according to different suggestions
         select case (method)
            case (UPSTREAM)
               limit=_ZERO_
            case (P1)
               FATAL "P1 advection method not yet implemented, choose other method"
               stop  "adv_center.F90"
            case ((P2),(P2_PDM))
               if (method.eq.P2) then
                  limit=Phi
               else
                  limit=max(_ZERO_,min(Phi,2./(1.-c),2.*r/(c+1.e-10)))
               end if
            case (Superbee)
               limit=max(_ZERO_, min(_ONE_, 2.0*r), min(r,2.*_ONE_) )
            case (MUSCL)
               limit=max(_ZERO_,min(2.*_ONE_,2.0*r,0.5*(1.0+r)))
            case default
               LEVEL3 method
               FATAL 'unkown advection method in adv_center()'
               stop
          end select

!        compute the limited flux
         cu(k)=ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc))

      end do

!     do the upper boundary conditions
      select case (Bcup)
      case (flux)
         cu(N) = - Yup              ! flux into the domain is positive
      case (value)
         cu(N) =  ww(N)*Yup
      case (oneSided)
         if (ww(N).ge._ZERO_) then
            cu(N) =  ww(N)*Y(N)
         else
            cu(N) = _ZERO_
         end if
      case (zeroDivergence)
         cu(N) = cu(N-1)
      case default
         FATAL 'unkown upper boundary condition type in adv_center()'
         stop
      end select


!     do the lower boundary conditions
      select case (Bcdw)
      case (flux)
         cu(0) =   Ydw               ! flux into the domain is positive
      case (value)
         cu(0) =  ww(0)*Ydw
      case (oneSided)
         if(ww(0).le._ZERO_) then
            cu(0) =  ww(0)*Y(1)
         else
            cu(0) = _ZERO_
         end if
      case (zeroDivergence)
         cu(0) = cu(1)
      case default
         FATAL 'unkown lower boundary condition type in adv_center()'
         stop
      end select

!     do the vertical advection step which will be used for prescribed
!     vertical flow velocity and for settling of suspended matter.

      if (mode.eq.0) then ! non-conservative
         do k=1,N
            Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/        &
                 h(k)-Y(k)*(ww(k)-ww(k-1))/h(k))
         enddo
      else                ! conservative
         do k=1,N
            Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
         enddo
      end if

   end do ! end of the iteration loop


#ifdef DEBUG
   write(*,*) 'Leaving adv_center()'
   write(*,*)
#endif

   return
   end subroutine adv_center
!EOC

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------
