##==============================================================================
## Transport in a one-dimensional finite difference grid
##==============================================================================

tran.1D <- function(C, C.up=C[1], C.down=C[length(C)],
     flux.up=NULL, flux.down=NULL, a.bl.up=NULL, C.bl.up=NULL,
		 a.bl.down=NULL, C.bl.down=NULL,
     D=0, v=0, AFDW=1, VF=1, A=1,
     dx, full.check = FALSE, full.output = FALSE) {


### INPUT CHECKS
  if (is.null(dx))
    stop("error: dx should be inputted ")

  N <- length(C)
  if (N == 0)
    stop("C should be a vector with numeric values")

  if (is.null(C.up))   C.up <- C[1]
  if (is.null(C.down)) C.down <- C[N]

  ## default infilling of grid parameters

  if (!is.list(AFDW)) AFDW <- list(int=AFDW)
  if (!is.list(D)) D <- list(int=D)
  if (!is.list(v)) v <- list(int=v)

  if (!is.list(VF))
    VF <- list(int=rep(VF,length.out=(N+1)),
                mid=0.5*(rep(VF,length.out=(N+1))[1:N]+
                rep(VF,length.out=(N+1))[2:(N+1)]))
  if (!is.list(A))
    A <- list(int=rep(A,length.out=(N+1)),
              mid=0.5*(rep(A,length.out=(length(C)+1))[1:N]+
              rep(A,length.out=(N+1))[2:(N+1)]))

  if (is.list(dx)) grid <- dx
  if (!is.list(dx))
    grid <- list(dx=rep(dx,length.out=N),
                 dx.aux=0.5*(c(0,rep(dx,length.out=N))+
                 c(rep(dx,length.out=N),0)))

  ## check dimensions of input arguments

  if (full.check) {

    ## check input of AFDW
    gn <- names(AFDW)
    if (! "int" %in% gn)
      stop("error: AFDW should be a list that contains 'int', the AFDW values at the interface of the grid cell ")
    if (is.null(AFDW$int))
      stop("error: AFDW is NULL, should contain (numeric) values")
    if (!is.null(AFDW$int)) {
      if (!((length(AFDW$int)==1) || (length(AFDW$int)==(N+1))))
        stop("error: AFDW should be a vector of length 1 or N+1")
    }
    if (any(AFDW$int < 0)||any(AFDW$int > 1))
	    stop("error: the AFDW should always range between 0 and 1")

    ## check input of D
    gn <- names(D)
    if (! "int" %in% gn)
      stop("error: D should be a list that contains 'int', the D values at the interface of the grid cell ")
    if (is.null(D$int))
      stop("error: D is NULL, should contain (numeric) values")
    if (!is.null(D$int)) {
      if (!((length(D$int)==1) || (length(D$int)==(N+1))))
        stop("error: D should be a vector of length 1 or N+1")
    }
    if (any(D$int < 0))
	    stop("error: the diffusion coefficient should always be positive")

    ## check input of v
    gn <- names(v)
    if (! "int" %in% gn)
      stop("error: v should be a list that contains 'int', the v values at the interface of the grid cell ")
    if (is.null(v$int))
      stop("error: the advective velocity v is NULL, should contain (numeric) values")
    if (!is.null(v$int)) {
      if (!((length(v$int)==1) || (length(v$int)==(N+1))))
        stop("error: v should be a vector of length 1 or N+1")
    }
## This constraint no longer holds I presume (filip)
##    if (any (v$int < 0) & any (v$int > 0))
##	    stop("error: the advective velocity cannot be both positive and negative within the same domain")

    ## check input of VF
    gn <- names(VF)
    if (! "int" %in% gn)
      stop("error: VF should be a list that contains 'int', the area values at the interface of the grid cell ")
    if (! "mid" %in% gn)
      stop("error: VF should be a list that contains 'mid', the area at the middle of the grid cells")
    if (is.null(VF$int) || is.null(VF$mid))
      stop("error: the volume fraction VF should contain (numeric) values")
    if (!is.null(VF$int)) {
      if (!((length(VF$int)==1) || (length(VF$int)==(N+1))))
        stop("error: VF$int should be a vector of length 1 or N+1")
    }
    if (!is.null(VF$mid)) {
      if (!((length(VF$mid)==1) || (length(VF$mid)==(N))))
        stop("error: VF$mid should be a vector of length 1 or N")
    }
    if (any(VF$int < 0) || any(VF$mid < 0) || any(VF$int > 1) ||any(VF$mid > 1))
      stop("error: the volume fraction should range between 0 and 1")

    ## check input of A
    gn <- names(A)
    if (! "int" %in% gn)
      stop("error: A should be a list that contains 'int', the area values at the interface of the grid cell ")
    if (! "mid" %in% gn)
      stop("error: A should be a list that contains 'mid', the area at the middle of the grid cells")
    if (is.null(A$int) || is.null(A$mid))
      stop("error: the surface area A is NULL, should contain (numeric) values")
    if (!is.null(A$int)) {
      if (!((length(A$int)==1) || (length(A$int)==(N+1))))
        stop("error: A$int should be a vector of length 1 or N+1")
    }
    if (!is.null(A$mid)) {
      if (!((length(A$mid)==1) || (length(A$mid)==(N))))
        stop("error: A$mid should be a vector of length 1 or N")
    }
    if (any (A$int < 0) || any (A$mid < 0))
      stop("error: the area A should always be positive")

    ## check input of grid
    gn <- names(grid)
    if (! "dx" %in% gn)
      stop("error: grid should be a list that contains 'dx' ")
    if (! "dx.aux" %in% gn)
	    stop("error: grid should be a list that contains 'dx.aux' ")
    if (is.null(grid$dx) || is.null(grid$dx.aux))
    	stop("error: the grid should be a list with (numeric) values for 'dx' and 'dx.aux' ")
    if (!is.null(grid$dx)) {
      if (!((length(grid$dx)==1) || (length(grid$dx)==N)))
        stop("error: dx should be a vector of length 1 or N")
    }
    if (!is.null(grid$dx.aux)) {
      if (!((length(grid$dx.aux)==1) || (length(grid$dx.aux)==(N+1))))
        stop("error: dx.aux should be a vector of length 1 or N+1")
    }
    if (any(grid$dx <= 0) || any(grid$dx.aux <= 0) )
    	stop("error: the grid distances dx and dx.aux should always be positive")

    ## check input of boundary layer parameters
    if (!is.null(a.bl.up) & !is.null(C.bl.up)) {
    	if (a.bl.up < 0)
        stop("error: the boundary layer transfer coefficient should be positive")
    }

    if (!is.null(a.bl.down) & !is.null(C.bl.down)) {
	    if (a.bl.down < 0)
        stop("error: the boundary layer transfer coefficient should be positive")
    }
  } # end full.check


### FUNCTION BODY: CALCULATIONS


  if (full.output) {
    ## Impose boundary flux at upper boundary when needed
    ## Default boundary condition is zero gradient
    if (! is.null (flux.up)) {
      if (v$int[1] >= 0) {
      ## advection directed downwards
        nom <- flux.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
               (1-AFDW$int[1])*v$int[1])*C[1]
        denom <- VF$int[1]*(D$int[1]/grid$dx.aux[1]+
                 AFDW$int[1]*v$int[1])
	    } else	{
      ## advection directed upwards
        nom <- flux.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
               AFDW$int[1]*v$int[1])*C[1]
        denom <- VF$int[1]*(D$int[1]/grid$dx.aux[1]+
                (1-AFDW$int[1])*v$int[1])
	    }
      C.up <- nom/denom
    }
    ## Impose boundary flux at lower boundary when needed
    ## Default boundary condition is no gradient
    if (! is.null (flux.down)) {
      if (v$int[N+1] >= 0) {
        ## advection directed downwards
	      nom <- flux.down - VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1] +
               AFDW$int[N+1]*v$int[N+1])*C[N]
        denom <- -VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1]+
                (1-AFDW$int[N+1])*v$int[N+1])
      } else {
        ## advection directed downwards
	      nom <- flux.down - VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1] +
               (1-AFDW$int[N+1])*v$int[N+1])*C[N]
        denom <- -VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1]+
                 AFDW$int[N+1]*v$int[N+1])
      }
      C.down <- nom/denom
    }

    ## when upstream boundary layer is present, calculate new C.up
    if (!is.null(a.bl.up) & !is.null(C.bl.up)) {
      if (v$int[1] >= 0)	{
        ## advection directed downwards
        nom <- a.bl.up*C.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
              (1-AFDW$int[1])*v$int[1])*C[1]
        denom <- a.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
                 AFDW$int[1]*v$int[1])
	    } else	{
        ## advection directed upwards
        nom <- a.bl.up*C.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
               AFDW$int[1]*v$int[1])*C[1]
        denom <- a.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
                (1-AFDW$int[1])*v$int[1])
	    }
      C.up <- nom/denom
    }

    ## when downstream boundary layer is present, calculate new C.up
    if (!is.null(a.bl.down) & !is.null(C.bl.down)) {
      if (v$int[N+1] >= 0)	{
        ## advection directed downwards
       nom <- a.bl.down*C.bl.down + VF$int[N+1]*(D$int[N+1]/
              grid$dx.aux[N+1] + (1-AFDW$int[N+1])*v$int[N+1])*C[N]
       denom <- a.bl.down + VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1] +
                AFDW$int[N+1]*v$int[N+1])
     	} else	{
        ## advection directed upwards
        nom <- a.bl.down*C.bl.down + VF$int[N+1]*(D$int[N+1]/
               grid$dx.aux[N+1] + AFDW$int[N+1]*v$int[N+1])*C[N]
        denom <- a.bl.down + VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1] +
                 (1-AFDW$int[N+1])*v$int[N+1])
	    }
      C.down <- nom/denom
    }

    ## Calculate diffusive part of the flux

    dif.flux <- as.vector(-VF$int*D$int*diff(c(C.up,C,C.down))/
                          grid$dx.aux)
    adv.flux <- rep(0,length.out=length(dif.flux))

    ## Add advective part of the flux if needed
      if (any (v$int > 0)) {   # advection directed downwards
       vv <- v$int
       vv[v$int<0] <- 0
	     conc <- AFDW$int*c(C.up,C)
	     if (any (AFDW$int < 1))
         conc <- conc +(1-AFDW$int)*c(C,C.down)
	     adv.flux <- adv.flux + as.vector(VF$int*vv*conc)
    }
    if (any (v$int < 0)) {   # advection directed upwards
       vv <- v$int
       vv[v$int>0] <- 0
       conc <- AFDW$int*c(C,C.down)
	     if (any (AFDW$int < 1))
         conc <- conc +(1-AFDW$int)*c(C.up,C)
	     adv.flux <- adv.flux + as.vector(VF$int*vv*conc)
    }

    flux <- dif.flux + adv.flux

  } else { # not full.output

    ## when upstream boundary layer is present, calculate new C.up
    if (!is.null(a.bl.up) & !is.null(C.bl.up))  {
      if (v$int[1] >= 0) { # advection directed downwards
        nom <- a.bl.up*C.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
               (1-AFDW$int[1])*v$int[1])*C[1]
        denom <- a.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
               AFDW$int[1]*v$int[1])
	    } else	{  # advection directed upwards
        nom <- a.bl.up*C.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
               AFDW$int[1]*v$int[1])*C[1]
        denom <- a.bl.up + VF$int[1]*(D$int[1]/grid$dx.aux[1] +
                 (1-AFDW$int[1])*v$int[1])
	    }
      C.up <- nom/denom
    }

    ## when upstream boundary layer is present, calculate new C.up
    if (!is.null(a.bl.down) & !is.null(C.bl.down)) {
      if (v$int[N+1] >= 0)	{   # advection directed downwards
        nom <- a.bl.down*C.bl.down + VF$int[N+1]*(D$int[N+1]/
               grid$dx.aux[N+1] + (1-AFDW$int[N+1])*v$int[N+1])*C[N]
        denom <- a.bl.down + VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1] +
                 AFDW$int[N+1]*v$int[N+1])
	    } else	{  # advection directed upwards
        nom <- a.bl.down*C.bl.down + VF$int[N+1]*(D$int[N+1]/
                 grid$dx.aux[N+1] + AFDW$int[N+1]*v$int[N+1])*C[N]
        denom <- a.bl.down + VF$int[N+1]*(D$int[N+1]/grid$dx.aux[N+1] +
                (1-AFDW$int[N+1])*v$int[N+1])
	    }
      C.down <- nom/denom
    }

    ## Calculate diffusive part of the flux
	  flux <- as.vector(-(VF$int)*D$int*
                      diff(c(C.up,C,C.down))/grid$dx.aux)
    ## Add advective part of the flux if needed
    if (any (v$int > 0)) {   # advection directed downwards
       vv <- v$int
       vv[v$int<0] <- 0
	     conc <- AFDW$int*c(C.up,C)
	     if (any (AFDW$int < 1))
         conc <- conc +(1-AFDW$int)*c(C,C.down)
	     flux <- flux + as.vector(VF$int*vv*conc)
    }
    if (any (v$int < 0)) {   # advection directed upwards
       vv <- v$int
       vv[v$int>0] <- 0
       conc <- AFDW$int*c(C,C.down)
	     if (any (AFDW$int < 1))
         conc <- conc +(1-AFDW$int)*c(C.up,C)
	     flux <- flux + as.vector(VF$int*vv*conc)
    }

  }

  if (! is.null (flux.up)) flux[1] <- flux.up
  if (! is.null (flux.down)) flux[N+1] <- flux.down

    
## Calculate rate of change = Flux gradient       
  dC <- -diff(A$int*flux)/A$mid/VF$mid/grid$dx

  if (!full.output){
    return (list (dC = dC,                   # Rate of change due to advective-diffusive transport in each grid cell
  					flux.up = flux[1],               # Flux across lower boundary interface; positive = IN
	  				flux.down = flux[length(flux)])) # Flux across lower boundary interface; positive = OUT

  } else {
  	return (list (dC = dC,                   # Rate of change due to advective-diffusive transport in each grid cell
	  				C.up = C.up,                     # Concentration at upper interface
		  			C.down = C.down,                 # Concentration at lower interface
			  		dif.flux = dif.flux,             # Diffusive flux across at the interface of each grid cell
				  	adv.flux = adv.flux,             # Advective flux across at the interface of each grid cell
					  flux = flux,                     # Flux across at the interface of each grid cell
  					flux.up = flux[1],               # Flux across lower boundary interface; positive = IN
	  				flux.down = flux[length(flux)])) # Flux across lower boundary interface; positive = OUT
  }

} # end tran.1D