!! -------------------------------------------------------
!! This is a set of routines that the user can supply to
!! determine the flow inside their own object
!! An simple example is already set up here to illustrate
!! how the routines work.  The geometry for this problem
!! is solid body rotation inside of a circle.  Flow is
!! initialized to a strip in circle.
!!
!! Details of each routine are provided in the following
!! comments
!! -------------------------------------------------------


!! -------------------------------------------------------
!! user_initialize : Initialize geometry
!!
!! Input parameters :
!!             mx, my  : number of grid cells in x and y
!!                       directions
!!             dx, dy  : mesh sizes in x and y dir.
!!              x,y    : vector of x, y locations of center
!!                       of mesh cells.
!!        ax,ay, bx,by : Lower left corner (ax,ay) and upper
!!                       right corner (bx,by) of computational
!!                       domain.
!! -------------------------------------------------------
SUBROUTINE user_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 /user_data/ r_user

  r_user = 1.5  !! Radius of circle for solid body  rotation
                !! Center is assumed to be at (0,0)

END SUBROUTINE user_initialize


!! -----------------------------------------------------
!! user_fbody : Define level set used to determine
!!              irregular physical domain.
!!
!! Input parameters :
!!            x,y :  point (x,y) at which level set is to
!!                   be evaluated.
!!
!! Values returned : 1  (x,y) is in flow-domain
!!                  -1  (x,y) is NOT in flow-domain
!! -----------------------------------------------------
REAL FUNCTION user_fbody(x,y)
  IMPLICIT REAL(a-h,o-z)
  COMMON /user_data/ r_user

  IF (SQRT(x**2 + y**2) < r_user) THEN
     user_fbody = 1
  ELSE
     user_fbody = -1
  END IF
END FUNCTION user_fbody


!! -------------------------------------------------------
!! user_ic : This defines initial function q(x,y,0).
!! parameters are :
!!             mx, my  : number of grid cells in x and y
!!                       directions
!!             dx, dy  : mesh sizes in x and y dir.
!!             centers : (x,y) location of Centers of mass
!!                       of cell (i,j).
!!
!!                        centers(i,j,1) : x location
!!                        centers(i,j,2) : y location
!!             kappa   : capacity function.  kappa(i,j) is
!!                       capacity of cell (i,j).  Note that
!!                       this is NOT the modified kappa, so
!!                       it is set to 0 for cells outside
!!                       of the body.
!!  Output parameters :
!!             q       : (array) advected quantity.
!!
!! Function calls : user_ic2 : This specifies value for cell
!!                  (i,j).  In most cases, this is the function
!!                  that the user should modify.
!! NOTE : q is a three-dimensional array.
!! -------------------------------------------------------
SUBROUTINE user_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) = user_ic2(xc,yc,kappa(i,j))
     END DO
  END DO

END SUBROUTINE user_ic

!! --------------------------------------------------------
!! user_ic2 : Definition of q(x,y,0).
!!
!! Input parameters :
!!           xci,ycj : location at which function q(xc,yc,0)
!!                   is to be evaluated.  This corresponds
!!                   to the center of mass of cell (i,j)
!!            kij  : capacity of cell (i,j).
!!
!! output parameters :
!!            user_ic2 : value of q(xci,ycj,0)
!! --------------------------------------------------------
REAL FUNCTION user_ic2(xci,ycj,kij)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"
  COMMON /user_data/ r_user

  REAL kij

  IF (ABS(ycj) <= r_user/2.0 .AND. ABS(kij) > 1e-13) THEN
     user_ic2 = 1.0
  ELSE
     user_ic2 = 0.0
  END IF

END FUNCTION user_ic2

!! -------------------------------------------------------
!! user_stream : Set up stream function.
!! Input values :
!!          mx, my : Number of grid cells in x and y
!!                   directions
!!          dx, dy : mesh size in x,y directions
!!          x,y    : vector of x, y locations of center
!!                   of cells
!! Output values :
!!           psi   : array of stream function values.
!! -------------------------------------------------------
SUBROUTINE user_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) = user_stream2(xi,yj)
     END DO
  END DO
END SUBROUTINE user_stream

!! -------------------------------------------------------
!! user_stream2 : value of stream function psi(x,y)
!! Input values :
!!
!!          x,y : scalar value at which psi is to be
!!                    evaluated.
!! output values :
!!      psi(x,y) : value of stream function at (x,y)
!!
!! -------------------------------------------------------
REAL FUNCTION user_stream2(x,y)
  COMMON /user_data/ r_user
  COMMON /compi/ pi
  r2 = x**2 + y**2
  user_stream2 = 0.2*pi*MIN(r2,r_user**2)
END FUNCTION user_stream2


!! -----------------------------------------------------
!! bc2_user :  Set up boundary conditions on
!!             computatinal domain for advection term.
!! Input parameters :
!!             mx, my  : number of grid cells in x and y
!!                       directions
!!             dx, dy  : mesh sizes in x and y dir.
!!              x,y    : vector of x, y locations of center
!!                       of mesh cells.
!!                  q  : (array) current value of advected
!!                       quantity.
!!                aux  : Three dimensional array containing
!!                       the horizontal and vertical velocities
!!                       and the modified capacity function values
!!                          aux(i,j,1) : horizontal velocity
!!                                       in cell (i,j)
!!                          aux(i,j,2) : vertical velocity
!!                                       in cell (i,j)
!!                          aux(i,j,3) : capacity function value
!!                                       in cell (i,j)
!!                   t  : Current time.
!!
!! Output :
!!                 q  : Modified so that ghost cells have been
!!                      filled in correctly.  Unless the physical
!!                      irregular boundary intersects the
!!                      computational domain, ghost cell values
!!                      can all be set to 0 (as shown below).
!!                      In some cases (channel flow, for example)
!!                      outflow boundary conditions must be specified
!!                      See CLAWPACK documentation for the best way
!!                      to do this.
!! -------------------------------------------------------------
SUBROUTINE bc2_user(mx,my,dx,dy,x,y,q,aux,t)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:mx+mbc,1-mbc:my+mbc,*) :: q, aux
  REAL, DIMENSION(-1:maxgrid) :: x,y


  !! Nothing needs to happen here, since flow domain is contained
  !! completely within computational domain.
  DO i = 1-mbc, mx+mbc
     DO jbc=1,mbc
        q(i,1-jbc,1) = 0.0  !! bottom edge
        q(i,my+jbc,1) = 0.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_user
