advection.1D <- function (C, C.up = C[1], C.down = C[length(C)],
                       flux.up = NULL, flux.down=NULL,
                       v, VF=1, A =1,  dx, dt.default = 1,
                       adv.method = c("muscl","super","quick","p3","up"),
                       full.check = FALSE) {

  # number of compartments
  n <- as.integer(length(C))

  adv.method <- match.arg(adv.method)
  advmet <- pmatch(adv.method,c("muscl","super","quick","p3","up"))
  if (is.na(advmet))
    stop ("'adv.method' not known: ",adv.method)
    
  # types of boundaries and inputs
  if (is.null(flux.up)) {
    upbnd  <- 2
    upval  <- C.up
  } else {
    upbnd  <- 1
    upval  <- flux.up
  }
  if (is.null(flux.down)) {
    dwnbnd  <- 2
    dwnval  <- C.down
  } else {
    dwnbnd  <- 1
    dwnval  <- flux.down
  }

  # timestep this works only for latest version of deSolve  !
  dt <- timestep(prev=FALSE)
#  print(dt)
  if (is.nan(dt)) dt <- dt.default
  if (dt < 1e-30) dt <- dt.default   # if dt is = 0 or ~0
  if (dt > 1e30)  dt <- dt.default

  # velocity, grid sizes, volume fractions, surface areas
  v     <- rep(v, length.out = n+1)
  
  if (is.list(dx)) {
    dx.aux <- dx$dx.aux
    dx <- dx$dx
  } else { 
    dx    <- rep(dx,length.out = n)
#    dx.aux<- 0.5*(c(0,rep(dx,length.out=n))+
#                  c(rep(dx,length.out=n),0))
    dx.aux <- rep(dx,length.out = n+1)

  }
  
  if (is.list(VF)) {
    VFint <- VF$int
    VFmid <- VF$mid
  } else { 
    VFint <- rep(VF,length.out=(n+1))
    VFmid <- 0.5*(rep(VF,length.out=(n+1))[1:n]+
                  rep(VF,length.out=(n+1))[2:(n+1)])
  }
  if (is.list(A)) {
    Aint <- A$int
    Amid <- A$mid
  } else { 
    Aint  <- rep(A,length.out=(n+1))
    Amid  <- 0.5*(rep(A,length.out=(length(C)+1))[1:n]+
                  rep(A,length.out=(n+1))[2:(n+1)])
  }
  storage.mode(dx) <- storage.mode(dx.aux) <- "double"
  storage.mode(VFint) <- storage.mode(VFmid) <- "double"
  storage.mode(Aint) <- storage.mode(Amid) <- "double"
  storage.mode(v) <- "double"

  # The advection routine ...
  O<-.Fortran("advection",n,as.double(C), as.double(dt),
   dx, dx.aux, v,
   as.integer(upbnd), as.integer(dwnbnd), as.double(upval), as.double (dwnval),
   VFint, VFmid, Aint, Amid,
   as.integer(advmet), as.integer(1), dy=as.double(rep(0.,n)),
   cu=as.double(rep(0.,n+1)), it=as.integer(0), package = "ReacTran")
   
  # return the "rate of change" and the fluxes
   list(dC=O$dy, adv.flux=O$cu, flux.up =O$cu[1], flux.down=O$cu[n+1], 
     it=O$it)
}

## =============================================================================
## volumetric advective transport
## =============================================================================

advection.volume.1D <- function (C, C.up = C[1], C.down = C[length(C)],
                       F.up = NULL, F.down=NULL,
                       flow, V, dt.default = 1,
                       adv.method = c("muscl","super","quick","p3","up"),
                       full.check = FALSE) {

  # number of compartments
  n <- as.integer(length(C))
  adv.method <- match.arg(adv.method)
  advmet <- pmatch(adv.method,c("muscl","super","quick","p3","up"))
  if (is.na(advmet))
    stop ("'adv.method' not known: ",adv.method)
    
  # types of boundaries and inputs
  if (is.null(F.up)) {
    upbnd  <- 2
    upval  <- C.up
  } else {
    upbnd  <- 1
    upval  <- F.up
  }
  if (is.null(F.down)) {
    dwnbnd  <- 2
    dwnval  <- C.down
  } else {
    dwnbnd  <- 1
    dwnval  <- F.down
  }

  # timestep - this works only for latest version of deSolve  !
  dt <- timestep(prev = FALSE)
  if (dt == 0)   dt <- dt.default
  if (dt > 1e30) dt <- dt.default
  
  # velocity, grid sizes, volume fractions, surface areas
  flow  <- rep(flow, length.out = n+1)
  V     <- rep(V,length.out = n)
  V.aux <- c(V[1],V)

  storage.mode(V) <- storage.mode(V.aux) <- "double"
  storage.mode(flow) <- "double"

  # The advection routine ...
  O<-.Fortran("advectvol",n,as.double(C), as.double(dt),
   V, V.aux, flow,
   as.integer(upbnd), as.integer(dwnbnd), 
   as.double(upval), as.double (dwnval),
   as.integer(advmet), as.integer(1), dy=as.double(rep(0.,n)),
   cu=as.double(rep(0.,n+1)), it=as.integer(0), package = "ReacTran")

  # return the "rate of cange" and the fluxes
   list(dC=O$dy, F=O$cu, F.up =O$cu[1], F.down=O$cu[n+1], it=O$it)
}