!! -------------------------------------------------------------------
!! src2_irr.f90
!!
!! This solves the parabolic equation in the advection-diffusion equation
!!          q_t = \nabla (D \nabla q)
!! over the times step dt_claw.  One of three implicit algorithms may
!! be chosen : Backward Euler; Crank-Nicolson; or Trapazoidal-BDF2.
!! To solve the resulting linear system, the stabilized biconjugate
!! gradient algorithm (BICG-stab) of Van der Vorst is used.
!!
!! Boundary conditions on the computational domain are handled
!! through the call to 'define_bdry', also in this file.  The current
!! implementation only handles flux (Neumann) boundary conditions,
!! (although it is essentially set up to handle Dirichlet as well).
!! Furthermore, only static (not time dependent) boundary conditions
!! are implemented in the present version.
!!
!! At the boundary of the physical (flow) domain, the code is only
!! designed to handle flux conditions at this time.  At this time,
!! only homogeneous flux (insulated wall) are properly handled.
!!
!! This routine is called from CLAWPACK in the case that the user
!! also specifies an advection term, or alone, if the user only wishes
!! to solve the diffusion equation in complex geometry. Fluxes at full
!! and partial cell edges are approximated using coefficients determined
!! at setup time.  These fluxes are computed in the file 'compute_fluxes'
!! -------------------------------------------------------------------

SUBROUTINE src2_irr(maxmx_arg,maxmy_arg,meqn_arg,mbc_arg,mx,my,qn, &
     & aux_arg,t,dt_claw)
  IMPLICIT REAL (a-h,o-z)

  INCLUDE "parms.h"

  EXTERNAL apk_flux

  !! input arguments
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: C, H, HI, bdry_seg, WI
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,meqn) :: qn, qnp1, Wij
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,3) :: beta
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux) :: aux, aux_apk
  REAL, DIMENSION(1:maxmx+1,1:maxmy+1) :: F,G, FI, GI

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

  COMMON /grid/ x,y
  COMMON /grid_data_real/ dT, dx, dy, ax,ay,bx, by, beta, aux
  COMMON /grid_data_int/ mx_apk,my_apk
  COMMON /method_choice/ method, matvec_choice


  eps = 1e-12
  dt_outer = dt_claw
  mx_apk = mx
  my_apk = my

  !! Get inhomogenous boundary conditions
  !! NOTE : We are assuming that the boundary conditions are not time
  !! dependent.

  !! Define normal derivatives at boundaries of computational domain
  !! and on irregular objects.
  CALL define_bdry(mx,my,dx,dy,C,H)
  !!
  FI(1:mx,1:my) = 0.0
  GI(1:mx,1:my) = 0.0
  DO j = 1,my
     FI(1,j) = -beta(1,j,1)*C(1,j)
     FI(mx+1,j) = -beta(mx+1,j,1)*C(mx+1,j)
  END DO

  DO i = 1,mx
     GI(i,1) =    -beta(i,1,2)*C(i,1)
     GI(i,my+1) = -beta(i,my+1,2)*C(i,my+1)
  END DO
  GI(1,1) = 0.0  !! Kludge needed for now!
  GI(1,my) = 0.0

  !! Only homogenous (H = 0) flux conditions are handled at this time.
  DO i = 1,mx
     DO j = 1,my
        HI(i,j) = -beta(i,j,3)*H(i,j)
     END DO
  END DO

  WI(1:mx,1:my) = 0.0
  DO i = 1,mx
     DO j = 1,my
        WI(i,j) = -(dy*(FI(i+1,j) - FI(i,j)) + dx*(GI(i,j+1) - GI(i,j)) + &
             & bdry_seg(i,j)*HI(i,j))
     END DO
  END DO

  SELECT CASE(method)
  CASE(1)  !! Backward Euler
     a1 = 1.0
     a2 = 0.0
  CASE(2)  !! Crank-Nicolson
     a1 = 0.5
     a2 = 0.5
  CASE(3)  !! Trapazoidal-BDF
     a1 = 0.25
     a2 = 0.25
  END SELECT

  totalTime = 0.0
  inner : DO
     dt = MIN(dt_claw,dx)
     IF (dt_outer - totaltime <= dt) THEN
        dt = dt_outer - totalTime
     END IF

     CALL compute_fluxes(mx,my,qn,F,G)

     DO i = 1,mx
        DO j = 1,my
           dtkpa = dt/(aux(i,j,3)*dx*dy)
           Wij(i,j,1) = -(dy*(F(i+1,j) - F(i,j)) + dx*(G(i,j+1) - G(i,j)))
           qn(i,j,1) = qn(i,j,1) + dtkpa*(a2*Wij(i,j,1) + (a1 + a2)*WI(i,j))
        END DO
     END DO

     resid = eps
     kmax = mx*my
     iter = kmax
     matvec_choice = 1
     CALL bicg(mx,my,qn,qnp1,apk_flux,iter,resid,info)
     IF (info /= 0) THEN
        WRITE(6,'(A,I5)') 'src2_irr : Info /= 0; info = ',info
        STOP
     END IF

     IF (method == 3) THEN
        DO i = 1,mx
           DO j = 1,my
              dtkpa = dt/(aux(i,j,3)*dx*dy)
              qn(i,j,1) = -qn(i,j,1)/3 + 4/3*qnp1(i,j,1) + dtkpa/3*WI(i,j)
           END DO
        END DO

        resid = eps
        kmax  = mx*my
        iter = kmax
        matvec_choice = 2
        CALL bicg(mx,my,qn,qnp1,apk_flux,iter,resid,info)
        IF (info /= 0) THEN
           WRITE(6,'(A,I5)') 'src2_irr : Info /= 0; info = ',info
           STOP
        END IF
     END IF
     qn(1:mx,1:my,1) = qnp1(1:mx,1:my,1)

     totalTime = totalTime + dt
     IF (ABS(totalTime - dT_outer) < 1e-14) EXIT
  END DO inner

  RETURN
END SUBROUTINE src2_irr


!! ---------------------------------------------------------
!! Compute fluxes F and G at full cell edges.  Note that
!! length of partial edges is taken into account in the
!! definition of beta(i,j).
!! ---------------------------------------------------------
SUBROUTINE compute_fluxes(mx,my, q, F, G)
  IMPLICIT REAL (A-H, O-Z)
  INCLUDE "parms.h"

  REAL, DIMENSION(1:maxmx+1, 1:maxmy+1) :: F,G
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux) :: aux
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,meqn) :: q
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,2,4) :: coeffs
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,3) :: beta
  INTEGER, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+Mbc,2,4,2) :: neighbors

  COMMON /interp_data_real/ coeffs
  COMMON /interp_data_int/ neighbors
  COMMON /grid_data_real/ dt, dx, dy, ax,ay,bx,by,beta,aux

  !! Fill in ghost cells;  Since we only have flux conditions, there is no
  !! Coefficient in front of the q.
  DO i = 1-mbc, mx + mbc
     q(i,0,1) = q(i,1,1)   !! bottom edge
     q(i,my+1,1) = q(i,my,1) !! top edge
  END DO

  DO j = 1-mbc, my + mbc
     q(0,j,1) = q(1,j,1)   !! left edge
     q(mx+1,j,1) = q(mx,j,1) !! right edge
  END DO


  !! Do two separate loops to avoid any problems computing
  !! vertical fluxes in bottom and top borders, and horizontal
  !! fluxes in left and right borders.
  iside = 1
  DO i = 1,mx+1
     DO j = 1,my
        px = 0
        DO k = 1,4
           ni = neighbors(i,j,iside,k,1)
           nj = neighbors(i,j,iside,k,2)
           px = px + q(ni,nj,1)*coeffs(i,j,iside,k)
        END DO
        F(i,j) = -beta(i,j,1)*px
     END DO
  END DO

  iside = 2
  DO i = 1,mx
     DO j = 1,my+1
        py = 0
        DO k = 1,4
           ni = neighbors(i,j,iside,k,1)
           nj = neighbors(i,j,iside,k,2)
           py = py + q(ni,nj,1)*coeffs(i,j,iside,k)
        END DO
        G(i,j) = -beta(i,j,2)*py
     END DO
  END DO
END SUBROUTINE compute_fluxes


!! ---------------------------------------------------------
!! Matrix vector multiply.
!! ---------------------------------------------------------
SUBROUTINE apk_flux(q,mx,my,apk)

  IMPLICIT REAL (A-H, O-Z)
  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,meqn) :: apk, q
  REAL, DIMENSION(1:maxmx+1, 1:maxmy+1) :: F, G
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,maux) :: aux
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,3) :: beta

  COMMON /grid_data_real/ dt, dx, dy,ax,ay,bx,by,beta,aux
  COMMON /method_choice/ method, matvec_choice

  CALL compute_fluxes(mx,my,q,F,G)

  !! Coefficients for multistep methods.
  SELECT CASE(method)
  CASE(1)
     a1 = 1.0  !! Backward Euler
  CASE(2)
     a1 = 0.5  !! Crank-Nicolson
  CASE(3)
     SELECT CASE(matvec_choice)
     CASE(1)
        a1 = 0.25  !! First pass of TR-BDF2
     CASE(2)
        a1 = 1.0/3.0  !! Second pass
     END SELECT
  END SELECT

  DO j = 1,my
     DO i = 1,mx
        dtkpa = dt/(aux(i,j,3)*dx*dy)
        Wij = -(dy*(F(i+1,j) - F(i,j)) + dx*(G(i,j+1) - G(i,j)))
        apk(i,j,1) = q(i,j,1) - dtkpa*a1*Wij
     END DO
  END DO

  RETURN
END SUBROUTINE apk_flux


!! ----------------------------------------------------
!! This is not set up correctly to handle inhomogenous
!! boundary flux conditions on physical (flow) domain.
!! ----------------------------------------------------
SUBROUTINE define_bdry(mx,my,dx,dy,C,H)
  IMPLICIT REAL(a-h,o-z)

  INCLUDE "parms.h"

  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: C, H

  !! Specify flux conditions on computational domain
  !! This is the inhomogenous part and may be set to
  !! a non-zero value.
  C(1:mx+1,1:my+1) = 0.0
  DO i = 1,mx
     C(i,1) = 0.0
     C(i,my+1) = 0.0
  END DO

  DO j = 1,my
     C(1,j) = 0.0
     C(mx+1,j) = 0.0
  END DO

  !! Now specify conditions on boudary segments.  Here, we
  !! need to know which cell Cij has a boundary segment,
  !! and what the midpoint of this segment is.  Currently,
  !! this information is not provided, so only homogenous
  !! (zero) flux conditions will work.
  !! To do this properly, the <geo>_setup.f90 files should
  !! probably handle the details of this call.
  DO i = 1-mbc,mx+mbc
     DO j = 1-mbc,my+mbc
        H(i,j) = 0.0
     END DO
  END DO

END SUBROUTINE define_bdry
