mod_scal.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!=======================================================================
! FVCOM Scalar Module  
!
!    contains methods:
!        Adv_Scal            => Advect a Scalar Quantity 
!        Vdif_Scal           => Vertical Diffusion of Scalar Quantity
!        Bcond_Scal_OBC      => Open Boundary Condition for Scalar
!        Bcond_Scal_PTsource => Point Sources of Scalar
!=======================================================================
Module scalar

  logical, parameter :: debug = .true. 

  contains
!==============================================================================|
! Calculate Horizontal Advection and Diffusion For Scalar (f)                  |
!==============================================================================|
!  Subroutine Adv_Scal(f,fn,d_fdis,fdis,d_fflux,fflux_obc,deltat,source)
!J. Ge for tracer advection
  Subroutine adv_scal(f,f0,f2,fn,d_fdis,fdis,d_fflux,fflux_obc,deltat,source)
!J. Ge for tracer advection

!------------------------------------------------------------------------------|

  use all_vars
  use lims, only: m,mt,n,nt,kbm1,kb
  use bcs
  use mod_obcs
  use mod_wd


  implicit none
  real(sp), intent(in ), dimension(0:mt,kb)      :: f
!J. Ge for tracer advection
  real(sp), intent(in ), dimension(0:mt,kb)      :: f0,f2
!J. Ge for tracer advection 
  real(sp), intent(out), dimension(0:mt,kb)      :: fn
  integer , intent(in )                          :: d_fdis
  real(sp), intent(in ), dimension(d_fdis)       :: fdis
  integer , intent(in )                          :: d_fflux
  real(sp), intent(out), dimension(d_fflux,kbm1) :: fflux_obc 
  real(sp), intent(in )                          :: deltat
  logical , intent(in )                          :: source

  !----------------local--------------------------------------
  real(sp), dimension(0:mt,kb)   :: xflux,xflux_adv
  real(sp), dimension(m)         :: pupx,pupy,pvpx,pvpy  
  real(sp), dimension(m)         :: pfpx,pfpy,pfpxd,pfpyd,viscoff
  real(sp), dimension(3*nt,kb)      :: dtij 
  real(sp), dimension(3*nt,kbm1) :: uvn
  real(sp), dimension(kb)        :: vflux
  real(sp) :: utmp,vtmp,sitai,ffd,ff1,x11,y11,x22,y22,x33,y33
  real(sp) :: tmp1,tmp2,xi,yi
  real(sp) :: dxa,dya,dxb,dyb,fij1,fij2,un
  real(sp) :: txx,tyy,fxx,fyy,viscof,exflux,temp,fpoint
  real(sp) :: fact,fm1,fmean
  integer  :: i,i1,i2,ia,ib,j,j1,j2,k,jtmp,jj
!J. Ge for tracer advectio
  real(sp),dimension(kb) :: ftmp
!J. Ge for tracer advection

  real(sp) :: da,db,ds



!------------------------------------------------------------------------------!

!-------------------------------------------------------
!Calculate Mean Values
!-------------------------------------------------------

  fmean = sum(f(1:m,1:kbm1))/float(m*kbm1)

!-------------------------------------------------------
!Initialize Multipliers to Control Horizontal Diff
!-------------------------------------------------------

  fact = 0.0_sp
  fm1  = 1.0_sp
  if(horizontal_mixing_type == 'closure') then
    fact = 1.0_sp
    fm1  = 0.0_sp
  end if
     
!-------------------------------------------------------
!Initialize Fluxes
!-------------------------------------------------------
  xflux     = 0.0_sp
  xflux_adv = 0.0_sp

!-------------------------------------------------------
!Calculate Normal Velocity on Control Volume Edges
!-------------------------------------------------------
!!# if !defined (1)
  do i=1,ncv
    i1=ntrg(i)
    !dtij(i)=dt1(i1)
    do k=1,kbm1
      dtij(i,k)=dt1(i1)*dz1(i1,k) 
      uvn(i,k) = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i)


    end do
  end do
!!# else
!!  do i=1,ncv
!!    i1=ntrg(i)
!!    dtij(i)=dt1(i1)
!!    do k=1,kbm1
!!      uvn(i,k) = vs(i1,k)*dltxe(i) - us(i1,k)*dltye(i)
!!    end do
!!  end do
!!# endif

!
!--Calculate the Advection and Horizontal Diffusion Terms----------------------!
!

   do k=1,kbm1
      pfpx  = 0.0_sp 
      pfpy  = 0.0_sp 
      pfpxd = 0.0_sp 
      pfpyd = 0.0_sp
     do i=1,m
       do j=1,ntsn(i)-1
         i1=nbsn(i,j)
         i2=nbsn(i,j+1)

!         IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 1)THEN
!          FFD=0.5_SP*(f(I,K)+f(I2,K))
!          FF1=0.5_SP*(f(I,K)+f(I2,K))
!    ELSE IF(ISWETN(I1) == 1 .AND. ISWETN(I2) == 0)THEN
!          FFD=0.5_SP*(f(I1,K)+f(I,K))
!          FF1=0.5_SP*(f(I1,K)+f(I,K))
!    ELSE IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 0)THEN
!          FFD=0.5_SP*(f(I,K)+f(I,K))
!          FF1=0.5_SP*(f(I,K)+f(I,K))
!    ELSE
!          FFD=0.5_SP*(f(I1,K)+f(I2,K))
!          FF1=0.5_SP*(f(I1,K)+f(I2,K))
!    END IF
 
!J. Ge for tracer advection
         IF(backward_advection/=.true.)THEN
           IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
            ffd=0.5_sp*(f(i,k)+f(i2,k))
            ff1=0.5_sp*(f(i,k)+f(i2,k))
       ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
            ffd=0.5_sp*(f(i1,k)+f(i,k))
            ff1=0.5_sp*(f(i1,k)+f(i,k))
       ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
            ffd=0.5_sp*(f(i,k)+f(i,k))
            ff1=0.5_sp*(f(i,k)+f(i,k))
       ELSE
            ffd=0.5_sp*(f(i1,k)+f(i2,k))
            ff1=0.5_sp*(f(i1,k)+f(i2,k))
       END IF
         ELSE
           IF(backward_step==1)THEN
             IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
               ffd=0.5_sp*((f0(i,k)+f(i,k))*0.5+(f0(i2,k)+f(i2,k))*0.5)
               ff1=0.5_sp*((f0(i,k)+f(i,k))*0.5+(f0(i2,k)+f(i2,k))*0.5)
             ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((f0(i1,k)+f(i1,k))*0.5+(f0(i,k)+f(i,k))*0.5)
               ff1=0.5_sp*((f0(i1,k)+f(i1,k))*0.5+(f0(i,k)+f(i,k))*0.5)
         ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((f0(i,k)+f(i,k))*0.5+(f0(i,k)+f(i,k))*0.5)
               ff1=0.5_sp*((f0(i,k)+f(i,k))*0.5+(f0(i,k)+f(i,k))*0.5)
         ELSE
               ffd=0.5_sp*((f0(i1,k)+f(i1,k))*0.5+(f0(i2,k)+f(i2,k))*0.5)
               ff1=0.5_sp*((f0(i1,k)+f(i1,k))*0.5+(f0(i2,k)+f(i2,k))*0.5)
             END IF
           ELSEIF(backward_step==2)THEN
             IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
               ffd=0.5_sp*((f2(i,k)+f0(i,k)+f(i,k))/3.0_sp+(f2(i2,k)+f0(i2,k)+f(i2,k))/3.0_sp)
               ff1=0.5_sp*((f2(i,k)+f0(i,k)+f(i,k))/3.0_sp+(f2(i2,k)+f0(i2,k)+f(i2,k))/3.0_sp)
             ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((f2(i1,k)+f0(i1,k)+f(i1,k))/3.0_sp+(f2(i,k)+f0(i,k)+f(i,k))/3.0_sp)
               ff1=0.5_sp*((f2(i1,k)+f0(i1,k)+f(i1,k))/3.0_sp+(f2(i,k)+f0(i,k)+f(i,k))/3.0_sp)
         ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
               ffd=0.5_sp*((f2(i,k)+f0(i,k)+f(i,k))/3.0_sp+(f2(i,k)+f0(i,k)+f(i,k))/3.0_sp)
               ff1=0.5_sp*((f2(i,k)+f0(i,k)+f(i,k))/3.0_sp+(f2(i,k)+f0(i,k)+f(i,k))/3.0_sp)
         ELSE
               ffd=0.5_sp*((f2(i1,k)+f0(i1,k)+f(i1,k))/3.0_sp+(f2(i2,k)+f0(i2,k)+f(i2,k))/3.0_sp)
               ff1=0.5_sp*((f2(i1,k)+f0(i1,k)+f(i1,k))/3.0_sp+(f2(i2,k)+f0(i2,k)+f(i2,k))/3.0_sp)
         END IF  
           ENDIF

         END IF
!J. Ge for tracer advection
     
         pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2))
         pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1))
         pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2))
         pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1))
       end do

! gather all neighboring control volumes connecting at dam node 

       pfpx(i)  =pfpx(i )/art2(i)
       pfpy(i)  =pfpy(i )/art2(i)
       pfpxd(i) =pfpxd(i)/art2(i)
       pfpyd(i) =pfpyd(i)/art2(i)

     end do
          
     if(k == kbm1)then
       do i=1,m
         pfpxb(i) = pfpx(i)
         pfpyb(i) = pfpy(i)
       end do
     end if

     do i=1,m
       pupx(i)=0.0_sp
       pupy(i)=0.0_sp
       pvpx(i)=0.0_sp
       pvpy(i)=0.0_sp
       j=1
       i1=nbve(i,j)
       jtmp=nbvt(i,j)
       j1=jtmp+1-(jtmp+1)/4*3
       j2=jtmp+2-(jtmp+2)/4*3
       x11=0.5_sp*(vx(i)+vx(nv(i1,j1)))
       y11=0.5_sp*(vy(i)+vy(nv(i1,j1)))
       x22=xc(i1)
       y22=yc(i1)
       x33=0.5_sp*(vx(i)+vx(nv(i1,j2)))
       y33=0.5_sp*(vy(i)+vy(nv(i1,j2)))

       pupx(i)=pupx(i)+u(i1,k)*(y11-y33)
       pupy(i)=pupy(i)+u(i1,k)*(x33-x11)
       pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33)
       pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11)

       if(isonb(i) /= 0) then
         pupx(i)=pupx(i)+u(i1,k)*(vy(i)-y11)
         pupy(i)=pupy(i)+u(i1,k)*(x11-vx(i))
         pvpx(i)=pvpx(i)+v(i1,k)*(vy(i)-y11)
         pvpy(i)=pvpy(i)+v(i1,k)*(x11-vx(i))
       end if

       do j=2,ntve(i)-1
         i1=nbve(i,j)
         jtmp=nbvt(i,j)
         j1=jtmp+1-(jtmp+1)/4*3
         j2=jtmp+2-(jtmp+2)/4*3
         x11=0.5_sp*(vx(i)+vx(nv(i1,j1)))
         y11=0.5_sp*(vy(i)+vy(nv(i1,j1)))
         x22=xc(i1)
         y22=yc(i1)
         x33=0.5_sp*(vx(i)+vx(nv(i1,j2)))
         y33=0.5_sp*(vy(i)+vy(nv(i1,j2)))

         pupx(i)=pupx(i)+u(i1,k)*(y11-y33)
         pupy(i)=pupy(i)+u(i1,k)*(x33-x11)
         pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33)
         pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11)
       end do
       j=ntve(i)
       i1=nbve(i,j)
       jtmp=nbvt(i,j)
       j1=jtmp+1-(jtmp+1)/4*3
       j2=jtmp+2-(jtmp+2)/4*3
       x11=0.5_sp*(vx(i)+vx(nv(i1,j1)))
       y11=0.5_sp*(vy(i)+vy(nv(i1,j1)))
       x22=xc(i1)
       y22=yc(i1)
       x33=0.5_sp*(vx(i)+vx(nv(i1,j2)))
       y33=0.5_sp*(vy(i)+vy(nv(i1,j2)))

       pupx(i)=pupx(i)+u(i1,k)*(y11-y33)
       pupy(i)=pupy(i)+u(i1,k)*(x33-x11)
       pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33)
       pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11)

       if(isonb(i) /= 0) then
         pupx(i)=pupx(i)+u(i1,k)*(y11-vy(i))
         pupy(i)=pupy(i)+u(i1,k)*(vx(i)-x11)
         pvpx(i)=pvpx(i)+v(i1,k)*(y11-vy(i))
         pvpy(i)=pvpy(i)+v(i1,k)*(vx(i)-x11)
       end if
       pupx(i)=pupx(i)/art1(i)
       pupy(i)=pupy(i)/art1(i)
       pvpx(i)=pvpx(i)/art1(i)
       pvpy(i)=pvpy(i)/art1(i)
       tmp1=pupx(i)**2+pvpy(i)**2
       tmp2=0.5_sp*(pupy(i)+pvpx(i))**2
       viscoff(i)=sqrt(tmp1+tmp2)*art1(i)
     end do
!     if(k == kbm1) then
!       ah_bottom(1:m) = horcon*(fact*viscoff(1:m) + fm1)
!     end if


     do i=1,ncv_i
       ia=niec(i,1)
       ib=niec(i,2)
       xi=0.5_sp*(xije(i,1)+xije(i,2))
       yi=0.5_sp*(yije(i,1)+yije(i,2))
       dxa=xi-vx(ia)
       dya=yi-vy(ia)
       dxb=xi-vx(ib)
       dyb=yi-vy(ib)

!       fij1=f(ia,k)+dxa*pfpx(ia)+dya*pfpy(ia)
!       fij2=f(ib,k)+dxb*pfpx(ib)+dyb*pfpy(ib)
!J. Ge for tracer advection
       IF(backward_advection/=.true.)THEN
         fij1=f(ia,k)+dxa*pfpx(ia)+dya*pfpy(ia)
         fij2=f(ib,k)+dxb*pfpx(ib)+dyb*pfpy(ib)
       ELSE
         IF(backward_step==1)THEN
           fij1=(f0(ia,k)+f(ia,k))*0.5+dxa*pfpx(ia)+dya*pfpy(ia)
           fij2=(f0(ib,k)+f(ib,k))*0.5+dxb*pfpx(ib)+dyb*pfpy(ib)
         ELSEIF(backward_step==2)THEN
           fij1=(f2(ia,k)+f0(ia,k)+f(ia,k))/3.0_sp+dxa*pfpx(ia)+dya*pfpy(ia)
           fij2=(f2(ib,k)+f0(ib,k)+f(ib,k))/3.0_sp+dxb*pfpx(ib)+dyb*pfpy(ib)
         ENDIF
       ENDIF
!J. Ge for tracer advection
       un=uvn(i,k)

!       viscof=horcon*(fact*(viscoff(ia)+viscoff(ib))*0.5_sp + fm1)
        viscof=(fact*0.5_sp*(viscoff(ia)*nn_hvc(ia)+viscoff(ib)*nn_hvc(ib)) + fm1*0.5_sp*(nn_hvc(ia)+nn_hvc(ib)))

       txx=0.5_sp*(pfpxd(ia)+pfpxd(ib))*viscof
       tyy=0.5_sp*(pfpyd(ia)+pfpyd(ib))*viscof

       fxx=-dtij(i,k)*txx*dltye(i)
       fyy= dtij(i,k)*tyy*dltxe(i)


       exflux=-un*dtij(i,k)* &
          ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy


       ! --------- new: 14-02-2012 ,K.Lettmann ------
       ! Limit the total flux according to the amount
       ! of mass within the respective neigbour control volumnes

       da = dz(ia,k)*dt(ia)
       db = dz(ib,k)*dt(ib)
       ds = 0.5_sp*( dz(ia,k) + dz(ib,k) )  ! Approximation of water depth for segment in sigma units

       exflux = exflux / ds ! trafo to old eflux unit (FVCOM 2.7)
       call limit_hor_flux_scal(f(ia,k),f(ib,k),ia,ib,da,db,ntrg(i),deltat,exflux)
       exflux = exflux * ds ! trafo back to new eflux unit (>FVCOM 3.1.4)
       ! --------------------------------------------


       xflux(ia,k)=xflux(ia,k)+exflux
       xflux(ib,k)=xflux(ib,k)-exflux

       xflux_adv(ia,k)=xflux_adv(ia,k)+(exflux-fxx-fyy)
       xflux_adv(ib,k)=xflux_adv(ib,k)-(exflux-fxx-fyy)

     end do
  end do !!sigma loop

!---------------------------------------------------------------------------------
! Accumulate Fluxes at Boundary Nodes
!---------------------------------------------------------------------------------
 

!---------------------------------------------------------------------------------
! Store Advective Fluxes at the Boundary
!---------------------------------------------------------------------------------
  do k=1,kbm1
     if(iobcn > 0) then
       do i=1,iobcn
         i1=i_obc_n(i)
         fflux_obc(i,k)=xflux_adv(i1,k)
       end do
     end if
  end do

!---------------------------------------------------------------------------------
! Calculate Vertical Advection Terms 
!---------------------------------------------------------------------------------

   do i=1,m 
     if(iswetn(i)*iswetnt(i) == 1) then
!J. Ge for tracer advection
       !w_tmp = wts(i,1:kb)
       !dz_tmp= dz(i,1:kbm1)
       IF(backward_advection/=.true.)THEN
         !ftmp = f(i,1:kbm1)
         !call calc_vflux(kbm1,ftmp,w_tmp,dz_tmp,vflux)
         call calc_vflux(kbm1,f(i,1:kbm1),wts(i,1:kb),dz(i,1:kbm1),vflux)
       ELSE
         IF(backward_step==1)THEN
           ftmp=(f0(i,1:kbm1)+f(i,1:kbm1))*0.5
         ELSEIF(backward_step==2)THEN
           ftmp=(f2(i,1:kbm1)+f0(i,1:kbm1)+f(i,1:kbm1))/3.0_sp
         ENDIF
         !call calc_vflux(kbm1,ftmp,w_tmp,dz_tmp,vflux)
         call calc_vflux(kbm1,ftmp,wts(i,1:kb),dz(i,1:kbm1),vflux)
       ENDIF
!J. Ge for tracer advection
!!!      vflux = 0.0_SP

     do k=1,kbm1
       if(isonb(i) == 2) then
         xflux(i,k)= vflux(k)*art1(i)   !/dz(i,k)
!JQI         xflux(i,k)= (vflux(k)-vflux(k+1))*art1(i)/dz(i,k)
       else
!JQI         xflux(i,k)=xflux(i,k)+ (vflux(k)-vflux(k+1))*art1(i)/dz(i,k)
         xflux(i,k)=xflux(i,k)+ vflux(k)*art1(i)    !/dz(i,k)
       end if
     end do
     end if
   end do

!-------------------------------------------------------
!Point Source                                      
!-------------------------------------------------------
  if(source)then  !!user specified

  if(river_ts_setting == 'calculated') then
    if(river_inflow_location == 'node') then
        do j=1,numqbc
          jj=inodeq(j)
          fpoint=fdis(j)
          do k=1,kbm1
            xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*fpoint !/dz(jj,k)
          end do
        end do
    else if(river_inflow_location == 'edge') then
      write(*,*)'scalar advection not setup for "edge" point source'
      stop
    end if
  end if

  else

  if(river_ts_setting == 'calculated')then   
    if(river_inflow_location == 'node') then
        do j=1,numqbc
          jj=inodeq(j)
          do k=1,kbm1
            fpoint = f(jj,k)
!J. Ge for tracer advection
            IF(backward_advection/=.true.)THEN
              fpoint = f(jj,k)
            ELSE
              IF(backward_step==1)THEN
                fpoint = (f0(jj,k)+f(jj,k))*0.5
              ELSEIF(backward_step==2)THEN
                fpoint = (f2(jj,k)+f0(jj,k)+f(jj,k))/3.0_sp
              ENDIF
            ENDIF
!J. Ge for tracer advection
            xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*fpoint! /dz(jj,k)
          end do
        end do
    else if(river_inflow_location == 'edge') then
      write(*,*)'scalar advection not setup for "edge" point source'
      stop
    end if
  end if

  endif
!------------------------------------------------------------------------
!Update Scalar Quantity
!------------------------------------------------------------------------

  do i=1,m
    if(iswetn(i)*iswetnt(i) == 1 )then
    do k=1,kbm1
      !fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i))
      fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/(dt(i)*dz(i,k))))*(dt(i)/dtfa(i)) 
    end do
    else
    do k=1,kbm1
      fn(i,k)=f(i,k)
    end do
    end if
  end do

  return
  End Subroutine adv_scal
!==============================================================================|

!==============================================================================|
! Vertical Diffusion of Scalar                                                 |
!==============================================================================|
  Subroutine vdif_scal(f,deltat)

  use mtridiagonal_scal
  use all_vars 

  Implicit None 
  Real(sp), intent(inout) :: f(0:mt,kb)
  Real(sp), intent(in   ) :: deltat
  !--local--------------------
  integer  :: i,k,ll
  real(sp) :: dsqrd,dfdz,visb
  real(sp) :: fsol(0:kb)


  call init_tridiagonal_scal(kb)

  Do i=1,m
     dsqrd = d(i)*d(i)

    !----------------------------------------------------------------
    !  Set up Diagonals of Matrix (lower=au,diag=bu,upper=cu)
    !----------------------------------------------------------------
    

    !Surface
    au(1) = 0.0
    cu(1)=      - deltat*(kh(i,2)+umol)/(dzz(i,1)*dz(i,1)*dsqrd)
    bu(1)=  1.0 - cu(1) 

    !Interior
    do k=2,kbm1-1
      au(k) =     - deltat*(kh(i,k  )+umol)/(dzz(i,k-1)*dz(i,k)*dsqrd)
      cu(k) =     - deltat*(kh(i,k+1)+umol)/(dzz(i,k  )*dz(i,k)*dsqrd)
      bu(k) = 1.0 - cu(k) - au(k) 
    end do

    !Bottom
     au(kbm1) =     - deltat*(kh(i,kbm1)+umol)/(dzz(i,kbm1-1)*dz(i,kbm1)*dsqrd)
     cu(kbm1) = 0.0
     bu(kbm1) = 1.0 - au(kbm1) 

    !----------------------------------------------------------------
    ! Set up RHS forcing vector and boundary conditions 
    !----------------------------------------------------------------
    do k=1,kbm1
      du(k) = f(i,k)
    end do

    !Free Surface: No flux

    !Bottom: No flux
      

    !----------------------------------------------------------------
    ! Solve 
    !----------------------------------------------------------------

     call tridiagonal_scal(kb,1,kbm1,fsol)
    
     !Transfer
     f(i,1:kbm1) = fsol(1:kbm1)

  End Do



  End Subroutine vdif_scal


!==============================================================================|
! Set Point Source Conditions for Scalar Function                              |
!==============================================================================|

  Subroutine bcond_scal_ptsource(f,fn,fdis)

!------------------------------------------------------------------------------|
  use all_vars
  use bcs
  use mod_obcs
  implicit none
  real(sp), intent(in ), dimension(0:mt,kb)      :: f 
  real(sp), intent(out), dimension(0:mt,kb)      :: fn
  real(sp), intent(in ), dimension(numqbc )      :: fdis
!--local-------------------------------------------
  integer  :: i,j,k,j1,j11,j22
!------------------------------------------------------------------------------|


!--------------------------------------------
! Set Source Terms
!--------------------------------------------
  if(river_ts_setting == 'specified') then
    if(numqbc > 0) then
      if(river_inflow_location == 'node') then
        do i=1,numqbc
          j11=inodeq(i)
          do k=1,kbm1
            fn(j11,k)=fdis(i)
          end do
        end do
      else if(river_inflow_location == 'edge') then
        do i=1,numqbc
          j11=n_icellq(i,1)
          j22=n_icellq(i,2)
          do k=1,kbm1
            fn(j11,k)=fdis(i)
            fn(j22,k)=fdis(i)
          end do
        end do
      end if
    end if
  end if

  return
  End Subroutine bcond_scal_ptsource 
!==============================================================================|
!==============================================================================|

!==============================================================================|
!   Set Boundary Conditions for Scalar Function on Open Boundary               |
!==============================================================================|

  Subroutine bcond_scal_obc(f,fn,fflux_obc,f_obc,deltat,alpha_nudge)

!------------------------------------------------------------------------------|
  use all_vars
  use bcs
  use mod_obcs
  implicit none
  real(sp), intent(in   ), dimension(0:mt,kb)      :: f 
  real(sp), intent(inout), dimension(0:mt,kb)      :: fn
  real(sp), intent(in   ), dimension(iobcn+1,kbm1) :: fflux_obc 
  real(sp), intent(in   ), dimension(iobcn       ) :: f_obc 
  real(sp), intent(in   )                          :: deltat
  real(sp), intent(in   )                          :: alpha_nudge 
!--local-------------------------------------------
  real(sp) :: f2d,f2d_next,f2d_obc,xflux2d,tmp
  integer  :: i,j,k,j1,j11,j22
!------------------------------------------------------------------------------|
       
!--------------------------------------------
! Set Scalar Value on Open Boundary
!--------------------------------------------
  if(iobcn > 0) then
    do i=1,iobcn
      j=i_obc_n(i)
      j1=next_obc(i)
      f2d=0.0_sp
      f2d_next=0.0_sp
      xflux2d=0.0_sp
      do k=1,kbm1
        f2d=f2d+f(j,k)*dz(j,k)
        f2d_next=f2d_next+fn(j1,k)*dz(j1,k)
        xflux2d=xflux2d+fflux_obc(i,k)*dz(j,k)
      end do
  
      if(uard_obcn(i) > 0.0_sp) then
        tmp=xflux2d+f2d*uard_obcn(i)
        f2d_obc=(f2d*dt(j)-tmp*deltat/art1(j))/d(j)
        do k=1,kbm1
          fn(j,k)=fn(j1,k) !f2d_obc+(fn(j1,k)-f2d_next)
        end do
      else
        do k=1,kbm1
          fn(j,k) = f(j,k)-alpha_nudge*(f(j,k)-f_obc(i))
        end do
      end if
    end do
  endif

  return
  End Subroutine bcond_scal_obc 
!==============================================================================|
!==============================================================================|

  Subroutine fct_sed(f,fn)
  !==============================================================================|
  USE all_vars
  USE mod_utils
  USE bcs
  USE mod_obcs
  IMPLICIT NONE
  real(sp), intent(inout), dimension(0:mt,kb)      :: fn
  real(sp), intent(in), dimension(0:mt,kb)      :: f
  REAL(SP):: SMAX,SMIN
  INTEGER :: I,J,K,K1
  !==============================================================================|
  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"Start: fct_sed"

  nodes: DO i=1,m

     ! SKIP OPEN BOUNDARY NODES
     IF(iobcn > 0)THEN
        DO j=1,iobcn
           IF(i == i_obc_n(j)) cycle nodes
        END DO
     END IF

     ! SKIP RIVER INFLOW POINTS
     IF(numqbc > 0)THEN
        DO j=1,numqbc
           IF(river_inflow_location == 'node')THEN
              IF(i == inodeq(j)) cycle nodes
           END IF
           IF(river_inflow_location == 'edge')THEN
              IF(i == n_icellq(j,1) .OR. i == n_icellq(j,2)) cycle nodes
           END IF
        END DO
     END IF

     ! SKIP GROUND WATER INFLOW POINTS
     IF(bfwdis(i) .GT. 0.0_sp .and. groundwater_salt_on) cycle nodes

     k1 = 1
     IF(precipitation_on) k1 = 2
!     DO K=1,KBM1
     DO k=k1,kbm1
        smax = maxval(f(nbsn(i,1:ntsn(i)),k))
        smin = minval(f(nbsn(i,1:ntsn(i)),k))

        IF(k == 1)THEN
           smax = max(smax,(f(i,k)*dz(i,k+1)+f(i,k+1)*dz(i,k))/  &
                (dz(i,k)+dz(i,k+1)))
           smin = min(smin,(f(i,k)*dz(i,k+1)+f(i,k+1)*dz(i,k))/  &
                (dz(i,k)+dz(i,k+1)))
        ELSE IF(k == kbm1)THEN
           smax = max(smax,(f(i,k)*dz(i,k-1)+f(i,k-1)*dz(i,k))/  &
                (dz(i,k)+dz(i,k-1)))
           smin = min(smin,(f(i,k)*dz(i,k-1)+f(i,k-1)*dz(i,k))/  &
                (dz(i,k)+dz(i,k-1)))
        ELSE
           smax = max(smax,(f(i,k)*dz(i,k-1)+f(i,k-1)*dz(i,k))/  &
                (dz(i,k)+dz(i,k-1)),                             &
                (f(i,k)*dz(i,k+1)+f(i,k+1)*dz(i,k))/           &
                (dz(i,k)+dz(i,k+1)))
           smin = min(smin,(f(i,k)*dz(i,k-1)+f(i,k-1)*dz(i,k))/  &
                (dz(i,k)+dz(i,k-1)),                             &
                (f(i,k)*dz(i,k+1)+f(i,k+1)*dz(i,k))/           &
                (dz(i,k)+dz(i,k+1)))
        END IF

        IF(smin-fn(i,k) > 0.0_sp)fn(i,k) = smin
        IF(fn(i,k)-smax > 0.0_sp)fn(i,k) = smax

     END DO
  END DO nodes

  WHERE(fn < 0.0_sp)fn=0.0_sp
  
  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"End: fct_sed"
  End Subroutine fct_sed

!==========================================================================
! Calculate Fluxes for Vertical Advection Equation                            
! n: number of cells
! c: scalar variable (1:n)
! w: velocity field at cell interfaces (1:n+1)
! note:  we dont use face normals to construct inflow/outflow
!        thus we add dissipation term instead of subtracting because 
!        positive velocity is up while computational coordinates increase
!        down towards bottom.  
!==========================================================================
  Subroutine calc_vflux(n,c,w,zk,flux) 
  use mod_prec
  use mod_utils, only : limled2
  implicit none
  integer , intent(in ) :: n
  real(sp), intent(in ) ::  c(n)
  real(sp), intent(in ) ::  w(n+1) 
  real(sp), intent(in ) ::  zk(n)
  real(sp), intent(out) ::  flux(n+1)
  real(sp) :: conv(n+1),diss(n+1)
  real(sp) :: cin(0:n+1)
  real(sp) :: sl_h(0:n+1)
  real(sp) :: dis4,sl_u,sl_f
  integer  :: i

  !transfer to working array
  cin(1:n) = c(1:n)
  sl_h(1:n) = zk(1:n)

  !surface bcs (no flux)
  cin(0)  =  -cin(1) 
!  cin(-1) =  -cin(2)
  sl_h(0) = zk(1)
  
  !bottom bcs (no flux)
  cin(n+1) = -cin(n) 
!  cin(n+2) = -cin(n-1)
  sl_h(n+1) = zk(n)

  !flux computation
  do i=2,n
    dis4    = .5*abs(w(i))
    conv(i) = w(i)*(cin(i)+cin(i-1))*0.5_sp
    sl_u = 2.0_sp*(cin(i)-cin(i+1))/(sl_h(i)+sl_h(i+1))
    sl_f = 2.0_sp*(cin(i-2)-cin(i-1))/(sl_h(i-2)+sl_h(i-1))
    diss(i) = dis4*(cin(i-1)-cin(i)-0.5_sp*limled2(sl_u,sl_f,2.0_sp)*(sl_h(i-1)+sl_h(i))) 
!!    flux(i) = conv(i)+diss(i)
  end do
  conv(1) = 0.0_sp
  diss(1) = 0.0_sp
  conv(n+1) = 0.0_sp
  diss(n+1) = 0.0_sp
  
  do i=1,n
    flux(i) = conv(i)-conv(i+1)+diss(i+1)-diss(i)
  end do

  End Subroutine calc_vflux
  
!==========================================================================
! Calculate LED Limiter L(u,v)  
!==========================================================================
  Function lim(a,b)
  use mod_prec
  real(sp) lim,a,b
  real(sp) q,r
  real(sp) eps
  eps = epsilon(eps)
  
  q = 0. !1st order
  q = 1. !minmod
  q = 2. !van leer

  r = abs(   (a-b)/(abs(a)+abs(b)+eps) )**q
  lim = .5*(1-r)*(a+b)

  End Function lim


! =============================================================================
! K.lettmann, Feb 2012, lettmann@icbm.de
!
! Methods for flux limiting for problems witht the diffusion step
! to omit negative concentrations
! ==============================================================================
!
!==============================================================================|
! Horizontal flux Limitation For Scalar (f)                  |
! 
! K.lettmann, Feb 2012, lettmann@icbm.de
!==============================================================================|
  Subroutine limit_hor_flux_scal(fka,fkb,ia,ib,Da,Db,ele_num,delt,exflux)
!
! This subroutine limits the total flux about a control volume segment
! in case there is not enough tracer in the source element, 
! the flux is coming from.
!------------------------------------------------------------------------------|

  use all_vars
  use lims, only: m,mt,n,nt,kbm1,kb
  use bcs
  use mod_obcs

  implicit none
  real(sp), intent(in )                  :: fka        ! scalar quantity on k-sigma level in cv a
  real(sp), intent(in )                  :: fkb        ! scalar quantity on k-sigma level in cv b
  integer , intent(in )                  :: ele_num    ! The element/triangle we are working
  integer , intent(in )                  :: ia,ib      ! Index of the left and right node to the SCV side S
  real(sp), intent(in )                  :: delt       ! time step for diffusion equation
  real(sp), intent(in )                  :: Da         ! water depth  of sigma layer [m] on side a of interface
  real(sp), intent(in )                  :: Db         ! water depth  of sigma layer [m] on side b of interface
  real(sp), intent(inout)                :: exflux     ! total flux over side S (diffusion + advection) []

  !----------------local--------------------------------------
  integer  :: i_source
  real(sp) :: x_ce,xa,xb,y_ce,ya,yb !coordinates of triangle nodes
  real(sp) :: f_source              !concentrations in source element
  real(sp) :: A_cv                  !area of part of control volume
  real(sp) :: delh                  !hight of sigma-layer
  real(sp) :: D_source              !water depth of sigma layer [m] at source node
  real(sp) :: amount_source         !amount of tracer in source node control volume
  real(sp) :: flux_max              !maximal possible outflux of tracer from source node control volume
!------------------------------------------------------------------------------!

  ! ----- the coordinates of the corner points of the double control volume ---
  x_ce = xc(ele_num) ;   y_ce = yc(ele_num)  ! center node
  xa = vx(ia) ;   ya = vy(ia)                ! left node
  xb = vx(ib) ;   yb = vy(ib)                ! right node
  !print*,'center node: ',x_ce,y_ce
  ! ---------------------------------------------------------------------------


  ! -------- calculate the area of the tracer control volume -------------------
  ! it is only the small part of the total control volume / not art1
   ! area for cartesian triangle on a plane
   a_cv = 0.5* abs( 0.5*( (xb-xa)*(y_ce-ya) - (x_ce-xa)*(yb-ya) ) )
  ! ----------------------------------------------------------------------------

  ! ------- find the source element, where the flux comes from -----------------
  if (exflux >= 0.0) then 
       i_source = ia
       f_source = fka
       d_source = da
  else
       i_source = ib
       f_source = fkb
       d_source = db
  end if
  !print*,'source node:',i_source
  ! -----------------------------------------------------------------------------

  ! ------- calculate the maximal possible outflow from source element ---------
  ! it cannot flow out more tracer than being in the element
  !print*,'f_source',f_source
  !print*,'D:',Dn(ia),Dn(ib)

  amount_source = f_source * a_cv
  flux_max = amount_source*d_source/delt

  !print*,'flux_max:',flux_max
  ! -----------------------------------------------------------------------------


  ! ------ limit the exflux if necessary ----------------------------------------
  if (exflux >= 0.0) then
    exflux = min(flux_max,exflux)
  else
    exflux = max(-flux_max,exflux)
  end if
  ! -----------------------------------------------------------------------------


  return
  End Subroutine limit_hor_flux_scal
!==============================================================================|

End Module scalar