!*********************************************************************
!*             Solving the Advection-Diffusion                       *
!*              Equation in Irregular Regions                        *
!*                                                                   *
!*********************************************************************
!*                                                                   *
!*                                                                   *
!*  This code reproduces the numerical results from the paper        *
!*                                                                   *
!*   Calhoun, D., LeVeque, R.J.,'A Cartesian Grid Finite-Volume      *
!*	Method for the Advection-Diffusion Equation in Irregular     *
!*	Geometries', submitted to Journal of Computational Physics   *
!*	1998.                                                        *
!*                                                                   *
!*  In addition, the code can be used for user-supplied geometries   *
!*  for which a stream function is known.                            *
!*                                                                   *
!*  Please direct all questions and comments to either one of the    *
!*  authors :                                                        *
!*                                                                   *
!*  calhoun@amath.washington.edu                                     *
!*  rjl@amath.washington.edu                                         *
!*                                                                   *
!*  This code can be downloaded from                                 *
!*                                                                   *
!*  http://www.amath.washington.edu/~rjl/clawpack/cartesian/advdiff/ *
!*                                                                   *
!*  See README file for more details on how to use this code         *
!*                                                                   *
!*  September 14, 1998                                               *
!*  Department of Applied Mathematics                                *
!*  University of Washington                                         *
!*  Seattle, WA  98195                                               *
!*                                                                   *
!*********************************************************************


PROGRAM advdiff
  IMPLICIT REAL (a-h,o-z)
  INCLUDE "parms.h"


  EXTERNAL bc2,rpn2ad1,src2_irr,rpt2ad1

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

  !! variables related to Claw2
  INTEGER :: mthlim(mwaves),method(7), nv(2)
  INTEGER :: adv_term, diff_term
  REAL cflv(4)
  REAL work(mwork)

  !! Grid variables
  REAL, DIMENSION(-1:maxgrid) :: x, y

  !! Time variables
  REAL dtv(5)

  !! Geometry
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: kappa, bdry_seg
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,2) :: centers, frside
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,2,2) :: midpts
  INTEGER, DIMENSION(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,2,4,2) :: neighbors
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,2,4) :: coeffs
  CHARACTER(len=LS) fname_output, fname


  COMMON /grid/ x,y
  COMMON /compi/ pi,pi2

  !! For boundary
  COMMON /problem_type/ iprob
  COMMON /interp_data_real/ coeffs
  COMMON /interp_data_int/ neighbors

  !! Solution method for diffusion step.
  COMMON /method_choice/ method_src, matvec_choice
  COMMON /parms/ D
  COMMON /grid_data_real/ dT, dx, dy, ax,ay,bx, by, beta, aux
  COMMON /grid_data_int/ mx, my

  OPEN(20,file='fort.info')
  OPEN(21,file='fort.x')
  OPEN(22,file='fort.y')
  OPEN(99,file='fort.dump')

  pi = 4.d0*ATAN(1.d0)
  pi2 = 2.d0*pi

  !! Input parameters for this run.
  WRITE(6,*) 'Input dx and dy.  (grid spacing) : '
  READ(5,*) dx,dy
  WRITE(6,*) 'Input dt and tend.  (time step and final time) : '
  READ(5,*) dt, T_final, nout
  WRITE(6,*) 'Input D : '
  READ(5,*) D
  WRITE(6,*) 'Input ax, ay, bx, by : '
  READ(5,*) ax, ay, bx, by
  WRITE(6,*) 'Input problem geometry :  '
  READ(5,*) iprob
  WRITE(6,*) 'Include diffusion term? (0/1) : '
  READ(5,*) diff_term
  WRITE(6,*) 'Include advection term? (0/1) : '
  READ(5,*) adv_term
  WRITE(6,*) 'Input location of large output files : '
  READ(5,*) fname_output
  WRITE(6,*) 'Input time step choice : '
  READ(5,*) itstep

  !! Define method to use for implicit diffusion step.
  method_src = 2   !! 1=backward Euler; 2 = trap; 3 =TRBDF

  len_output = istrlen(fname_output)

  !! -----------------------------------------
  !!  Set various parameters for clawpack
  !! -----------------------------------------
  cflv(1) = 2.5d0  ! maximum Courant number allowed.
  cflv(2) = 0.90  !  Unused, since method(1) = 0
                   ! clv(3:4) : output variables

  dtv(1) = dt  ! This will be in each pass through inner loop
  dtv(2) = 1e6   ! Unused, since method(1) = 0
                ! dtv(3:4) : output variables

  nv(1) = 1000   ! Unused, since method(1) = 0
              ! nv(2) : output variables

!!        mthlim(mw) = 0 for no limiter
!!                   = 1 for minmod
!!                   = 2 for superbee
!!                   = 3 for van Leer
!!                   = 4 for monotonized centered
  mthlim(1) = 4  ! 0 = Don't use a limiter

  method(1) = 0  ! 0 = Take fixed time step
  method(2) = 2  ! Use first order increment waves
  method(3) = 2  ! tranverse correction waves
  method(4) = 0  ! Kill verbose mode
  method(5) = diff_term  ! no source term
  method(6) = 3  ! no capacity function
  method(7) = 3  ! three auxilliary variables

  !! ---------------------------------------
  !! Set up grid
  !! ---------------------------------------
  WRITE(6,*) 'Setting up grid...'
  mx = (bx-ax + 1d-10)/dx
  my = (by-ay + 1d-10)/dy

  IF (mx > maxmx .OR. my > maxmy) THEN
     WRITE(6,*) 'Advection : mx or my too large; ', mx, my
     STOP
  END IF

  DO i=1-mbc,  mx + mbc
     x(i) = ax + (i-0.5d0) * dx
     IF (i >= 1 .AND. i <= mx)  THEN
        WRITE(21,2001) x(i)
     END IF
  END DO

  DO j = 1-mbc,my+mbc
     y(j) = ay + (j-0.5d0) * dy
     IF (j >= 1 .AND. j <= my)  THEN
        WRITE(22,2001) y(j)
     END IF
  END DO
2001 FORMAT(e16.8)
  WRITE(6,*) 'done'

  !! ------------------------------------------------
  !! Write out data for use by matlab scripts
  !! ------------------------------------------------
  OPEN(10,file='fort.matlab')
  WRITE(10,'(2I8)') mx, my
  WRITE(10,'(2F16.8)') dx, dy
  WRITE(10,'(4F8.4)') ax, ay, bx, by
  WRITE(10,'(I5)') iprob
  WRITE(10,*) fname_output
  CLOSE(10)

  !! ---------------------------------------------------
  !! Set up irregular boundaries and other geometry things
  !! ---------------------------------------------------
  WRITE(6,*) 'Setting up code...'
  CALL initialize_bdry(mx,my,dx,dy,x,y,ax,ay,bx,by,iprob)

  CALL get_geometry(mx,my,dx,dy,x,y,frside,kappa, midpts, &
     & centers, bdry_seg)

  IF (diff_term == 1) THEN
     CALL set_diffusion(mx,my,dx,dy,x,y,frside, kappa, &
          & centers, midpts, neighbors, coeffs)
  END IF

  !! --------------------------------------------------
  !! Initialize stream function and auxiliary variables
  !! --------------------------------------------------
  IF (adv_term == 1) THEN
     CALL set_psi(mx,my,dx,dy,x,y,psi,iprob)
  ELSE
     psi = 0
  END IF

  CALL set_aux(mx,my,dx,dy,x,y,frside,kappa,psi,D,aux,beta)

  fname = fname_output(1:len_output)//'fort.kappa'
  CALL write_q(maxmx,maxmy,mbc,LS,mx,my,kappa,0.0,fname,0)

  !! ---------------------------------------------------
  !! Set time information.
  !! ---------------------------------------------------

  T_start = 0.0
  dT_outer = T_final/nout
  IF (adv_term == 1) THEN
     CALL compute_tstep(mx,my,dx,dy,x,y,aux,frside,itstep,dt_min)
  ELSE
     dt_min = dT_outer
  END IF
220 FORMAT(A,1PE12.4)

  OPEN(11,file='fort.nplot')
  WRITE(11,'(i5)') nout
  CLOSE(11)

  !! --------------------------------------------------
  !! Summarize setup
  !! --------------------------------------------------
  WRITE(6,*) ' '
  WRITE(6,120) 'mx =      ',mx
  WRITE(6,120) 'my =      ',my
  WRITE(6,121) 'dx =      ',dx
  WRITE(6,121) 'dy =      ',dy
  WRITE(6,123) 'dt =      ',dt_min
  WRITE(6,121) '(ax,bx) = ',ax, bx
  WRITE(6,121) '(ay,by) = ',ay, by
  SELECT CASE(iprob)
  CASE(channel)
     WRITE(6,122) 'Problem geometry : Channel'
  CASE(annulus)
     WRITE(6,122) 'Problem geometry : Annulus'
  CASE(irr_objs)
     WRITE(6,122) 'Problem geometry : Irregular objects '
  CASE(user_domain)
     WRITE(6,122) 'Problem geometry : User domain '
  END SELECT
  WRITE(6,120) 'Advection term? ', adv_term
  WRITE(6,120) 'Diffusion term? ', diff_term
  WRITE(6,122) 'Output location for files : ',fname_output
  WRITE(6,122) ' '

120 FORMAT(A,I12)
121 FORMAT(A,2F12.4)
122 FORMAT(A,A)
123 FORMAT(A,1PE12.4)


  !! ---------------------------------------------------
  !! Initialize q
  !! ---------------------------------------------------

  WRITE(6,*) 'Initializing q...'
  CALL ic(mx,my,dx,dy,x,y,iprob,centers,kappa,q)

  fname = fname_output(1:len_output)//'fort.'
  CALL write_q(maxmx,maxmy,mbc,LS,mx,my,q,0.0,fname,100)

  !! ---------------------------------------------------
  !! Begin time steps.
  !! ---------------------------------------------------

  WRITE(6,*) 'Beginning time steps..'
  outer : DO n = 1,nout  ! Begin outer loop

     t_nm1 = T_start + (n-1) * dT_outer
     t_n = t_nm1 + dT_outer
     icount = 0
     time = t_n

     totalTime = 0
     inner : DO  ! if totalTime < dT

        IF (dT_outer - totalTime <= dt_min) THEN
           dtv(1) = dT_outer - totalTime
        ELSE
           dtv(1) = dt_min
        END IF
        tstep = dtv(1)
        t_inner_0 = t_nm1 + totalTime
        t_inner_1 = t_inner_0 + dtv(1)
        CALL bc2(maxmx,maxmy,meqn,mbc,mx,my,q,aux,t_n)
        IF (adv_term == 1) THEN
           CALL claw2(maxmx,maxmy,meqn,mwaves,mbc,mx,my, &
                &  q,aux,dx,dy,t_inner_0,t_inner_1,dtv, &
                &  cflv,nv,method,mthlim,&
                &  work,mwork,info,bc2,rpn2ad1,rpt2ad1,src2_irr)
        ELSE
           CALL src2_irr(maxmx,maxmy,meqn,mbc,mx,my,q,aux,t_n,dtv(1))
        END IF

        IF (info /= 0) THEN
           WRITE(6,*) 'Info = ',info
           STOP
        END IF

        CALL bc2(maxmx,maxmy,meqn,mbc,mx,my,q,aux,t_n)

        icount = icount + 1
        totalTime = totalTime + dtv(1)
        IF (totalTime >= dT_outer) EXIT
     END DO inner ! End inner loop

     CALL write_q(maxmx,maxmy,mbc,LS,mx,my,q,t_n,fname,100+n)
     WRITE(6,601) icount,t_n,info
601  FORMAT('CLAW2 returning after ',i5,' steps,  t =',e10.5, &
          '  info =',i3,/)

200  FORMAT(1P3E16.8)

     WRITE(20,1010) t_n,info,dtv(3),dtv(4),dtv(5), &
          cflv(3),cflv(4),nv(2)

1010 FORMAT( 't_n =         ',F12.4,/, &
          &  'info =        ',I12,/,&
          &  'smallest dt = ',1PE12.4,/,&
          &  'largest dt  = ',1PE12.4,/,&
          &  'last dt =     ',1PE12.4,/,&
          &  'largest cfl = ',F12.2,/,&
          &  'last cfl =    ',F12.2,/,&
          &  'steps taken = ',I12,/)

  END DO outer ! end outer loop
END PROGRAM advdiff
