!! -----------------------------------------------------
!! Initialize channel;  Set up any common blocks needed
!! by other routines
!! -----------------------------------------------------
SUBROUTINE channel_initialize(mx,my,dx,dy,x,y,ax,ay,bx,by)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMPLEX sz
  REAL, DIMENSION(-1:maxgrid) :: x, y
  COMMON /channel/ alpha,width,sz

  OPEN(12,file='channel.dat')
  READ(12,*) alpha,width
  READ(12,*) sx,sy
  sz = CMPLX(sx,sy)
  CLOSE(12)

END SUBROUTINE channel_initialize


!! -----------------------------------------------------
!! Define level set to describe body;
!! Values returned : 1  (x,y) is in flow-domain
!!                  -1  (x,y) is NOT in flow-domain
!! -----------------------------------------------------
REAL FUNCTION channel_fbody(x,y)

  COMPLEX sz, z,zi,w
  COMMON /channel/ alpha,width,sz
  COMMON /compi/ pi

  w = CMPLX(x,y)
  zi = (0,1)
  z = (w-sz)*EXP(-zi*alpha)
  IF (ABS(imag(z)) <= width) THEN
     channel_fbody = 1
  ELSE
     channel_fbody = -1
  END IF

END FUNCTION channel_fbody



!! -----------------------------------------------------
!! Initialize q.  This calls channel_ic2, below
!! -----------------------------------------------------
SUBROUTINE channel_ic(mx,my,dx,dy,x,y,centers,kappa,q)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMPLEX sz, zi, w, z

  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: kappa
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,*) :: q, centers
  REAL, DIMENSION(1-mbc:maxgrid) :: x,y

  DO i = 1-mbc,mx+mbc
     DO j = 1-mbc,my+mbc
        xc = centers(i,j,1)
        yc = centers(i,j,2)
        q(i,j,1) = channel_ic2(xc,yc,kappa(i,j))
     END DO
  END DO

END SUBROUTINE channel_ic


REAL FUNCTION channel_ic2(xc,yc,kij)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMPLEX sz, zi, w,z
  REAL kij
  COMMON /channel/ alpha,width,sz
  COMMON /parms/ D
  COMMON /compi/ pi

  t0 = 1.0
  zi = (0,1)
  w = CMPLX(xc,yc)
  z = (w-sz)*EXP(-zi*alpha)
  sij = REAL(z)
  wij = AIMAG(z)
  IF (ABS(wij) <= width) THEN
     channel_ic2 = 0.5*(erf((0.75 - (sij+3.25))/SQRT(4*D*t0)) + &
          &          erf((0.75 + sij+3.25)/SQRT(4*D*t0)))
  ELSE
     channel_ic2 = 0.0
  END IF
END FUNCTION channel_ic2



!! -----------------------------------------------------
!! Set up stream function.  Note that values in ghost
!! cells must also be set.
!! -----------------------------------------------------
SUBROUTINE channel_stream(mx,my,dx,dy,x,y,psi)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: psi
  REAL, DIMENSION(-1:maxgrid) :: x,y

  DO i = 1-mbc,mx+mbc
     DO j = 1-mbc,my+mbc
        xi = x(i) - dx/2
        yj = y(j) - dy/2
        psi(i,j) = channel_stream2(xi,yj)
     END DO
  END DO
END SUBROUTINE channel_stream

REAL FUNCTION channel_stream2(x,y)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMPLEX w,zi,z,sz
  COMMON/channel/ alpha,width,sz

  w = CMPLX(x,y)
  zi = (0,1)
  z = (w-sz)*EXP(-zi*alpha)
  IF (ABS(imag(z)) <= width) THEN
     channel_stream2 = -imag(z);
  ELSEIF (imag(z) < width) THEN
     channel_stream2 = width
  ELSE
     channel_stream2 = -width
  END IF
  channel_stream2 = -channel_stream2

END FUNCTION channel_stream2


!! -----------------------------------------------------
!! Set up boundary conditions for advection term.
!! -----------------------------------------------------
SUBROUTINE bc2_channel(mx,my,dx,dy,x,y,q,aux,t)
  IMPLICIT REAL(a-h,o-z)

  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,*) :: q, aux

  DO i = 1-mbc, mx+mbc
     DO jbc=1,mbc
        q(i,1-jbc,1) = 0  !! bottom edge
        q(i,my+jbc,1) = 0 !! top edge
     END DO
  END DO

  DO j = 1-mbc, my+mbc
     DO ibc=1,mbc
        q(1-ibc,j,1) = 0.0    !! left edge
        q(mx+ibc,j,1) = q(mx,j,1)   !! outflow conditions
     END DO
  END DO
END SUBROUTINE bc2_channel
