!! -----------------------------------------------------
!! Initialize flow-domain; set up any common blocks
!! needed by other routines
!! -----------------------------------------------------

SUBROUTINE annulus_initialize(mx,my,dx,dy,x,y,ax,ay,bx,by)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  REAL, DIMENSION(-1:maxgrid) :: x,y

  COMMON /annulus/ radius_annl(2), center_annl(2)

  OPEN(12,file='annulus.dat')
  READ(12,*) radius_annl(1), radius_annl(2)
  READ(12,*) center_annl(1), center_annl(2)
  CLOSE(12)

END SUBROUTINE annulus_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 annulus_fbody(x,y)
  IMPLICIT REAL(a-h,o-z)

  COMMON /annulus/ radius_annl(2), center_annl(2)
  COMMON /compi/ pi

  d2 = (x - center_annl(1))**2 + (y - center_annl(2))**2
  IF (d2 <= radius_annl(2)**2 .AND. d2 >= radius_annl(1)**2) THEN
     annulus_fbody = 1
  ELSE
     annulus_fbody = -1
  END IF

END FUNCTION annulus_fbody


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

  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) = annulus_ic2(xc,yc,kappa(i,j))
     END DO
  END DO
END SUBROUTINE annulus_ic

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

  REAL kij
  COMMON /annulus/ radius_annl(2), center_annl(2)
  COMMON /compi/ pi

  xc = xc - center_annl(1)
  yc = yc - center_annl(2)
  rij = SQRT(xc**2 + yc**2)
  thij = ATAN2(yc,xc)
  IF (ABS(thij-pi/2) <= pi/6 .AND. rij >= radius_annl(1) &
       & .AND. rij <= radius_annl(2)) THEN
     annulus_ic2 = 1.0
  ELSE
     annulus_ic2 = 0.0
  END IF
END FUNCTION annulus_ic2


!! -----------------------------------------------------
!! Set up stream function.  Note that values in ghost
!! cells must also be set.
!! -----------------------------------------------------
SUBROUTINE annulus_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) = annulus_stream2(xi,yj)
     END DO
  END DO

END SUBROUTINE annulus_stream

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

  COMMON /annulus/ radius_annl(2), center_annl(2)
  COMMON /compi/ pi

  u0 = 0.2*pi
  r1sq = radius_annl(1)**2
  r2sq = radius_annl(2)**2
  rsq = (x-center_annl(1))**2 + (y-center_annl(2))**2
  IF (rsq < r1sq) THEN
     annulus_stream2 = u0*(r2sq-r1sq)
  ELSEIF (rsq <= r2sq .AND. rsq >= r1sq) THEN
     annulus_stream2 = u0*(r2sq-rsq)
  ELSEIF (rsq > r2sq) THEN
     annulus_stream2 = 0.0
  END IF

  RETURN
END FUNCTION annulus_stream2


!! -----------------------------------------------------
!! Set up boundary conditions for advection term.
!! -----------------------------------------------------
SUBROUTINE bc2_annulus(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) = 0.0   !! right edge
     END DO
  END DO
END SUBROUTINE bc2_annulus
