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

  COMPLEX cu
  COMPLEX, DIMENSION(0:maxn-1,maxi) :: xk,yk, dsdtk
  REAL, DIMENSION(maxi) :: arclength, radius_irr
  REAL, DIMENSION(maxi,2) :: centers_irr
  REAL, DIMENSION(1:2*maxi) :: xkdat, ykdat

  COMMON /interface_int/  n, m_bodies
  COMMON /interface_cplx/ xk, yk, dsdtk
  COMMON /interface_real/ arclength,centers_irr, radius_irr

  cu = (0,1)

  OPEN(10,file='irr/xk')
  OPEN(11,file='irr/yk')
  OPEN(12,file='irr/arclength')
  OPEN(13,file='irr/dsdt')
  OPEN(14,file='irr/circles_irr.out')

  READ(12,*) n,m_bodies

  IF (m_bodies > maxi) THEN
     WRITE(6,*) 'initialize_interface : m_bodies > maxi; ', m_bodies
     STOP
  END IF
  IF (n > maxn) THEN
     WRITE(6,*) 'Initialize interface : n > maxn', n
     STOP
  END IF

  DO i = 0,n-1
     READ(10,*) (xkdat(j),j=1,2*m_bodies)
     READ(11,*) (ykdat(j),j=1,2*m_bodies)
     DO j = 1,m_bodies
        xk(i,j) = xkdat(2*j - 1) + cu*xkdat(2*j)
        yk(i,j) = ykdat(2*j - 1) + cu*ykdat(2*j)
     END DO
  END DO

  DO j = 1,m_bodies
     READ(12,*) arclength(j)
  END DO

  DO j = 1,m_bodies
     DO i = 0,n-1
        READ(13,*) dr, di
        dsdtk(i,j) = dr + cu*di
     END DO
  END DO

  DO j = 1,m_bodies
     READ(14,*) r1,c1,c2
     radius_irr(j) = r1
     centers_irr(j,1) = c1
     centers_irr(j,2) = c2
  END DO


  CLOSE(10)
  CLOSE(11)
  CLOSE(12)
  CLOSE(13)
  CLOSE(14)

  CALL generate_fine_mesh(dx,dy)

END SUBROUTINE irr_initialize

!! ----------------------------------------------------
!! This is needed to get an accurate winding number
!! calculation to determine whether we are inside of a
!! flow domain or not.
!! ----------------------------------------------------
SUBROUTINE generate_fine_mesh(dx,dy)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  REAL, DIMENSION(0:maxn-1) :: ki
  COMPLEX, DIMENSION(0:maxn-1,maxi) :: xk,yk, dsdtk
  REAL, DIMENSION(maxi) :: arclength
  COMPLEX, DIMENSION(0:maxnz,maxi) :: z, dz
  COMPLEX cu

  COMMON /interface_int/  n, m_bodies
  COMMON /interface_cplx/ xk, yk, dsdtk
  COMMON /interface_real/ arclength
  COMMON /interface_fine/ z, dz
  COMMON /interface_fine_int/ nz(maxi)
  COMMON /compi/ pi

  cu = (0,1)
  sca = 2*pi/n

  DO i = 0,n/2
     ki(i) = float(i)*sca
  END DO

  DO i = n-1,n/2+1,-1
     ki(i) = -float(n-i)*sca
  END DO

  DO j = 1,m_bodies
     nz(j) = CEILING(arclength(j)*10/MIN(dx,dy))
     IF (nz(j) > maxnz) THEN
        WRITE(6,*) 'generate_fine_mesh : nz(j) > maxnz;', &
             & nz(j), maxnz
        STOP
     END IF
     DO i = 0,nz(j)-1
        ti = float(i*n)/float(nz(j))
        sumx = 0
        sumy = 0
        sumxt = 0
        sumyt = 0
        DO k = 0,n-1
           sumx = sumx +               REAL(xk(k,j)*EXP(cu*ki(k)*ti))
           sumy = sumy +               REAL(yk(k,j)*EXP(cu*ki(k)*ti))
           sumxt = sumxt +   REAL( cu*ki(k)*xk(k,j)*EXP(cu*ki(k)*ti))
           sumyt = sumyt +   REAL( cu*ki(k)*yk(k,j)*EXP(cu*ki(k)*ti))
        END DO
        z(i,j) = sumx + sumy*cu
        dz(i,j) = sumxt + sumyt*cu
     END DO
  END DO
END SUBROUTINE generate_fine_mesh


!! -----------------------------------------------------
!! Define level set to describe body;
!! Values returned : 1  (x,y) is in flow-domain
!!                  -1  (x,y) is NOT in flow-domain
!! This routine relies on other routines, defined below.
!! -----------------------------------------------------
REAL FUNCTION irr_fbody(x,y)
  IMPLICIT REAL(a-h,o-z)
  IF (fdisc_irr(x,y) == 1) THEN
     irr_fbody = 1
  ELSE
     irr_fbody = -1
  END IF
END

!! Is (x,y) in body j?
REAL FUNCTION f_inbody_j(x,y,j)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMPLEX zp, cu
  COMPLEX, DIMENSION(0:maxnz,maxi) :: z, dz
  REAL, DIMENSION(maxi,2) :: centers_irr
  REAL, DIMENSION(maxi) :: arclength, radius_irr

  COMMON /interface_real/ arclength, centers_irr, radius_irr
  COMMON /interface_int/  n, m_bodies
  COMMON /interface_fine/ z, dz
  COMMON /interface_fine_int/ nz(maxi)

  cu = (0,1)
  zp = x + y*cu

  !! This function should return 1 if (x,y) is in jth body,
  !! and zero otherwise.
  d2 = (x - centers_irr(j,1))**2 + (y - centers_irr(j,2))**2
  IF (d2 > radius_irr(j)**2) THEN
     f_inbody_j = 0
     RETURN
  END IF

  !! Winding number calculation
  sum = 0
  DO i = 0,nz(j)-1
     IF (z(i,j) /= zp) THEN
        sum = sum + AIMAG(dz(i,j)/(z(i,j) - zp))
     END IF
  END DO
  f_inbody_j = MAX(MIN(sum,1.0),0.0)
END FUNCTION f_inbody_j


!! Returns -1 if (x,y) is in some inclusion, and 1 if
!! (x,y) is not in any inclusion
REAL FUNCTION fdisc_irr(x,y)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMMON /interface_int/  n, m_bodies

  DO j = 1,m_bodies
     IF (ABS(f_inbody_j(x,y,j)) > 0.35) THEN
        fdisc_irr = -1
        RETURN
     END IF
  END DO
  fdisc_irr = 1
  RETURN
END FUNCTION fdisc_irr

!! Returns which body (x,y) is in.
INTEGER FUNCTION jfdisc(x,y)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"

  COMMON /interface_int/  n, m_bodies
  DO j = 1,m_bodies
     IF (abs(f_inbody_j(x,y,j))  > 0.35) THEN
        jfdisc = j
        RETURN
     END IF
  END DO
  jfdisc = 0
  RETURN
END FUNCTION jfdisc

!! -----------------------------------------------------
!! Initialize q.  This calls irr_ic2, below
!! -----------------------------------------------------
SUBROUTINE irr_ic(mx,my,dx,dy,x,y,centers,kappa,q)
  IMPLICIT REAL(a-h,o-z)

  INCLUDE "parms.h"
  REAL irr_ic2
  REAL, DIMENSION(1-mbc:maxgrid) :: x,y
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: kappa
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,*) :: q, centers

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

END SUBROUTINE irr_ic

REAL FUNCTION irr_ic2(xc,yc,kij)
  IMPLICIT REAL(a-h,o-z)
  INCLUDE "parms.h"
  REAL kij
  irr_ic2 = 0
END FUNCTION irr_ic2


!! -----------------------------------------------------
!! Set up stream function.  This is read from a file.
!! Note that values in ghost cells must also be set.
!! -----------------------------------------------------
SUBROUTINE irr_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

  OPEN(10,file='irr/psi.dat')
  READ(10,*) mxp, myp, nzp
  rf = LOG(REAL((mxp-1)/mx))/LOG(2.0)
  IF (ABS(INT(rf) - rf) > 1e-12) THEN
     WRITE(6,'(A,2I5, F5.2)') 'lrefine is not an integer;', mxp, mx, rf
     STOP
  END IF
  CLOSE(10)

  WRITE(6,'(A,F4.0)') 'Level of refinement : ',rf
  OPEN(10,file='irr/psi')
  DO n = 1,nzp
     READ(10,*) ir,jr,pr
     csi = (ir - 1)/(2**rf) + 1
     csj = (jr - 1)/(2**rf) + 1
     IF (ABS(INT(csi) - csi) < 1e-12 .AND. ABS(INT(csj) - csj) < 1e-12) THEN
        ics = INT(csi)
        jcs = INT(csj)
        psi(ics,jcs) = pr
     END IF
  END DO
  CLOSE(10)
  !! Now set values in ghost cells to 0.
  DO i = 1-mbc,mx+mbc
     DO jbc = 1,mbc
        psi(i,1-jbc) = 0
        IF (jbc > 1) THEN
           psi(i,my+jbc) = 0
        END IF
     END DO
  END DO
  DO ibc = 1,mbc
     DO j = 1-mbc,my+mbc
        psi(1-ibc,j) = 0.0
        IF (ibc > 1) THEN
           psi(mx+ibc,j) = 0.0
        END IF
     END DO
  END DO
END SUBROUTINE irr_stream


!! -----------------------------------------------------
!! Set up boundary conditions for advection term.
!! -----------------------------------------------------
SUBROUTINE bc2_irr(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
        IF (t <= 20) THEN  !! Inject one unit per time step for 20 time units
           q(1-ibc,j,1) = 1.0
        ELSE
           q(1-ibc,j,1) = 0.0    !! left edge
        ENDIF
        q(mx+ibc,j,1) = q(mx,j,1)   !! outflow conditions
     END DO
  END DO
END SUBROUTINE bc2_irr
