scalar Module Reference

Functions/Subroutines
subroutine	adv_scal (f, f0, f2, fn, d_fdis, fdis, d_fflux, fflux_obc, deltat, source)
subroutine	vdif_scal (f, deltat)
subroutine	bcond_scal_ptsource (f, fn, fdis)
subroutine	bcond_scal_obc (f, fn, fflux_obc, f_obc, deltat, alpha_nudge)
subroutine	fct_sed (f, fn)
subroutine	calc_vflux (n, c, w, zk, flux)
real(sp) function	lim (a, b)
subroutine	limit_hor_flux_scal (fka, fkb, ia, ib, Da, Db, ele_num, delt, exflux)

Variables
logical, parameter	debug = .true.

Function/Subroutine Documentation

adv_scal()

subroutine scalar::adv_scal	(	real(sp), dimension(0:mt,kb), intent(in)	f,
		real(sp), dimension(0:mt,kb), intent(in)	f0,
		real(sp), dimension(0:mt,kb), intent(in)	f2,
		real(sp), dimension(0:mt,kb), intent(out)	fn,
		integer, intent(in)	d_fdis,
		real(sp), dimension(d_fdis), intent(in)	fdis,
		integer, intent(in)	d_fflux,
		real(sp), dimension(d_fflux,kbm1), intent(out)	fflux_obc,
		real(sp), intent(in)	deltat,
		logical, intent(in)	source
	)

Definition at line 60 of file mod_scal.f90.

!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

Here is the call graph for this function:

bcond_scal_obc()

subroutine scalar::bcond_scal_obc	(	real(sp), dimension(0:mt,kb), intent(in)	f,
		real(sp), dimension(0:mt,kb), intent(inout)	fn,
		real(sp), dimension(iobcn+1,kbm1), intent(in)	fflux_obc,
		real(sp), dimension(iobcn ), intent(in)	f_obc,
		real(sp), intent(in)	deltat,
		real(sp), intent(in)	alpha_nudge
	)

Definition at line 650 of file mod_scal.f90.

!------------------------------------------------------------------------------|
  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

bcond_scal_ptsource()

subroutine scalar::bcond_scal_ptsource	(	real(sp), dimension(0:mt,kb), intent(in)	f,
		real(sp), dimension(0:mt,kb), intent(out)	fn,
		real(sp), dimension(numqbc ), intent(in)	fdis
	)

Definition at line 601 of file mod_scal.f90.

!------------------------------------------------------------------------------|
  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

calc_vflux()

subroutine scalar::calc_vflux	(	integer, intent(in)	n,
		real(sp), dimension(n), intent(in)	c,
		real(sp), dimension(n+1), intent(in)	w,
		real(sp), dimension(n), intent(in)	zk,
		real(sp), dimension(n+1), intent(out)	flux
	)

Definition at line 790 of file mod_scal.f90.

  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

Here is the call graph for this function:

Here is the caller graph for this function:

fct_sed()

subroutine scalar::fct_sed	(	real(sp), dimension(0:mt,kb), intent(in)	f,
		real(sp), dimension(0:mt,kb), intent(inout)	fn
	)

Definition at line 703 of file mod_scal.f90.

  !==============================================================================|
  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"

Here is the call graph for this function:

lim()

real(sp) function scalar::lim	(	real(sp)	a,
		real(sp)	b
	)

Definition at line 842 of file mod_scal.f90.

  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)

limit_hor_flux_scal()

subroutine scalar::limit_hor_flux_scal	(	real(sp), intent(in)	fka,
		real(sp), intent(in)	fkb,
		integer, intent(in)	ia,
		integer, intent(in)	ib,
		real(sp), intent(in)	Da,
		real(sp), intent(in)	Db,
		integer, intent(in)	ele_num,
		real(sp), intent(in)	delt,
		real(sp), intent(inout)	exflux
	)

Definition at line 871 of file mod_scal.f90.

!
! 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

Here is the caller graph for this function:

vdif_scal()

subroutine scalar::vdif_scal	(	real(sp), dimension(0:mt,kb), intent(inout)	f,
		real(sp), intent(in)	deltat
	)

Definition at line 528 of file mod_scal.f90.

  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

Here is the call graph for this function:

Variable Documentation

debug

logical, parameter scalar::debug = .true.

Definition at line 51 of file mod_scal.f90.

  logical, parameter :: debug = .true.