/*@@
  @file      ADM_IterativeCN.F
  @date      August 1998
  @desc 
     Iterative crank nicholson in 3D
  @enddesc 
@@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_DefineThorn.h"

#include "CactusEinstein/Einstein/src/Einstein.h"


      subroutine ADM_ICN_Init(CCTK_ARGUMENTS)

      implicit  none

      DECLARE_CCTK_ARGUMENTS

      icn_iteration = 1

      return
      end


/*@@
  @routine    ADM_IterativeCN
  @date       August 1998
  @desc 
     Iterative Crank Nicholson in 3D
  @enddesc 
  @calls     
  @calledby   
  @history 
  @endhistory  
@@*/

      subroutine ADM_IterativeCN(CCTK_ARGUMENTS)

      USE ADM_Scalars

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer :: i,j,k,nx,ny,nz
      integer :: icn_maxit
      integer :: ierror

      CCTK_REAL :: dt,dx,dy,dz,halfdt
      CCTK_REAL :: 
     &     adm_rhsK_metric_xx,adm_rhsK_metric_xy,adm_rhsK_metric_xz,
     &     adm_rhsK_metric_yy,adm_rhsK_metric_yz,adm_rhsK_metric_zz,
     &     adm_rhsK_matter_xx,adm_rhsK_matter_xy,adm_rhsK_matter_xz,
     &     adm_rhsK_matter_yy,adm_rhsK_matter_yz,adm_rhsK_matter_zz

      CCTK_REAL :: zero,half,one,pi

      CCTK_REAL, DIMENSION(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: mask

c     Declarations for macros from ADM Evolution

#include "macro/ALPSOURCES_declare.h"
#include "macro/GSOURCES_declare.h"
#include "macro/KSOURCES_declare.h"


c     ==================================================================
c
c     INITIAL STUFF
c
c     ==================================================================

c     Numbers.

      zero = 0.0D0
      half = 0.5D0
      one  = 1.0D0

      pi = acos(-1.0D0)

c     Grid parameters.

      dt = cctk_delta_time

      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)

      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)


c     ==================================================================
c
c     SAVE OLD TIME STEP
c
c     ==================================================================

      if (icn_iteration == 1) then

         if (ADM_evolveLapse) then
            ADM_alp_p = alp
         end if

         ADM_gxx_p = gxx
         ADM_gxy_p = gxy
         ADM_gxz_p = gxz
         ADM_gyy_p = gyy
         ADM_gyz_p = gyz
         ADM_gzz_p = gzz

         ADM_kxx_p = kxx
         ADM_kxy_p = kxy
         ADM_kxz_p = kxz
         ADM_kyy_p = kyy
         ADM_kyz_p = kyz
         ADM_kzz_p = kzz

      end if


c     ==================================================================
c
c     THE ITERATIVE LOOP
c
c     ==================================================================

c     ICN parameters.

      if (ICN_itnum < 0) then
         icn_maxit = 50
      else
         icn_maxit = ICN_itnum
      end if

c     Do a new iteration if required.

      if (icn_iteration <= icn_maxit) then

c        Make the sources.

        CALL ADM_Lapse(CCTK_ARGUMENTS)

#include "adm_icn_inloop.inc"

         do k=2,nz-1
            do j=2,ny-1
               do i=2,nx-1

#include "macro/KSOURCES_guts.h"

                  if (ADM_evolveLapse) then
#include "macro/ALPSOURCES_guts.h"
                  end if

#include "macro/KSOURCES_undefine.h"
#include "macro/ALPSOURCES_undefine.h"
                  
               end do
            end do
         end do

         if (evolveg_boundary.eq.1) then

            if ((shift_state.ne.SHIFT_INACTIVE)) then
               call CCTK_WARN(0,"evolveg_boundary = yes  is incompatible with a non-zero shift")
            end if

            do k=1,nz
               do j=1,ny
                  do i=1,nx

#include "macro/GSOURCES_guts.h"
#include "macro/GSOURCES_undefine.h"

                  end do
               end do
            end do

         else

            do k=2,nz-1
               do j=2,ny-1
                  do i=2,nx-1

#include "macro/GSOURCES_guts.h"
#include "macro/GSOURCES_undefine.h"

                  end do
               end do
            end do

         end if

#include "adm_preloop_ADM.inc"

c        In ICN, all iterations except the last one
c        jump only half a time step.

         if (icn_iteration < icn_maxit) then
            halfdt = half*dt
         else
            halfdt = dt
         end if

c        Update variables.

         if (ADM_evolveLapse) then

            do k=2,nz-1
               do j=2,ny-1
                  do i=2,nx-1
                     alp(i,j,k) = ADM_alp_p(i,j,k) + halfdt*adms_alp(i,j,k)
                  end do
               end do
            end do

         end if

         if (evolveg_boundary.eq.1) then

            do k=1,nz
               do j=1,ny
                  do i=1,nx

                     gxx(i,j,k) = ADM_gxx_p(i,j,k) + halfdt*adms_gxx(i,j,k)
                     gxy(i,j,k) = ADM_gxy_p(i,j,k) + halfdt*adms_gxy(i,j,k)
                     gxz(i,j,k) = ADM_gxz_p(i,j,k) + halfdt*adms_gxz(i,j,k)
                     gyy(i,j,k) = ADM_gyy_p(i,j,k) + halfdt*adms_gyy(i,j,k)
                     gyz(i,j,k) = ADM_gyz_p(i,j,k) + halfdt*adms_gyz(i,j,k)
                     gzz(i,j,k) = ADM_gzz_p(i,j,k) + halfdt*adms_gzz(i,j,k)

                  end do
               end do
            end do

         else

            do k=2,nz-1
               do j=2,ny-1
                  do i=2,nx-1

                     gxx(i,j,k) = ADM_gxx_p(i,j,k) + halfdt*adms_gxx(i,j,k)
                     gxy(i,j,k) = ADM_gxy_p(i,j,k) + halfdt*adms_gxy(i,j,k)
                     gxz(i,j,k) = ADM_gxz_p(i,j,k) + halfdt*adms_gxz(i,j,k)
                     gyy(i,j,k) = ADM_gyy_p(i,j,k) + halfdt*adms_gyy(i,j,k)
                     gyz(i,j,k) = ADM_gyz_p(i,j,k) + halfdt*adms_gyz(i,j,k)
                     gzz(i,j,k) = ADM_gzz_p(i,j,k) + halfdt*adms_gzz(i,j,k)

                  end do
               end do
            end do

         end if

         do k=2,nz-1
            do j=2,ny-1
               do i=2,nx-1

                  kxx(i,j,k) = ADM_kxx_p(i,j,k) + halfdt*adms_kxx(i,j,k)
                  kxy(i,j,k) = ADM_kxy_p(i,j,k) + halfdt*adms_kxy(i,j,k)
                  kxz(i,j,k) = ADM_kxz_p(i,j,k) + halfdt*adms_kxz(i,j,k)
                  kyy(i,j,k) = ADM_kyy_p(i,j,k) + halfdt*adms_kyy(i,j,k)
                  kyz(i,j,k) = ADM_kyz_p(i,j,k) + halfdt*adms_kyz(i,j,k)
                  kzz(i,j,k) = ADM_kzz_p(i,j,k) + halfdt*adms_kzz(i,j,k)

               end do
            end do
         end do

c        Excision.

         if (excise.eq.1) then

#ifdef DEVELOPMENT_SIMPLEEXCISION

            if (ADM_evolveLapse) then
               call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .            alp,ADM_alp_p,one,one,'dt_extrapolation')
            end if

            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         gxx,ADM_gxx_p,one,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         gyy,ADM_gyy_p,one,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         gzz,ADM_gzz_p,one,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         gxy,ADM_gxy_p,zero,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         gxz,ADM_gxz_p,zero,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         gyz,ADM_gyz_p,zero,one,'dt_extrapolation')

            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         kxx,ADM_kxx_p,one,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         kyy,ADM_kyy_p,one,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         kzz,ADM_kzz_p,one,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         kxy,ADM_kxy_p,zero,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         kxz,ADM_kxz_p,zero,one,'dt_extrapolation')
            call SimpleExcision(_CCTK_ARGUMENTS,x,y,z,r,mask,
     .         kyz,ADM_kyz_p,zero,one,'dt_extrapolation')

#else
            call CCTK_WARN(0,"You have not compiled with SimpleExcision")
#endif

         end if

c        Staticrad boundary.  This is a hack, it really
c        shouldnot be done this way.  I need to talk to
c        Gab or Gerd later about it.

         if (CCTK_EQUALS(bound,'staticrad')) then

            if (ADM_evolveLapse) then
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         alp,ADM_alp_p,ADM_alp0,1.0D0)
            end if

            if (evolveg_boundary.eq.0) then
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         gxx,ADM_gxx_p,ADM_gxx0,1.0D0)
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         gyy,ADM_gyy_p,ADM_gyy0,1.0D0)
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         gzz,ADM_gzz_p,ADM_gzz0,1.0D0)
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         gxy,ADM_gxy_p,ADM_gxy0,1.0D0)
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         gxz,ADM_gxz_p,ADM_gxz0,1.0D0)
               call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .         gyz,ADM_gyz_p,ADM_gyz0,1.0D0)
            end if

            call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .      kxx,ADM_kxx_p,ADM_kxx0,1.0D0)
            call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .      kyy,ADM_kyy_p,ADM_kyy0,1.0D0)
            call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .      kzz,ADM_kzz_p,ADM_kzz0,1.0D0)
            call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .      kxy,ADM_kxy_p,ADM_kxy0,1.0D0)
            call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .      kxz,ADM_kxz_p,ADM_kxz0,1.0D0)
            call staticrad(cctk_bbox,nx,ny,nz,dt,dx,dy,dz,x,y,z,r,
     .      kyz,ADM_kyz_p,ADM_kyz0,1.0D0)

         end if

c        Boundaries.

         if (ADM_evolveLapse) then
            call ADM_Boundary(cctkGH,"einstein::lapse")
         end if

         if (evolveg_boundary.eq.0) then
            call ADM_Boundary(cctkGH,"einstein::metric")
         end if

         call ADM_Boundary(cctkGH,"einstein::curv")

#include "adm_icn_bound1.inc"

c        Synchronize.

         if (ADM_evolveLapse) then
            call CCTK_SyncGroup(ierror,cctkGH,"einstein::lapse")
         end if

         call CCTK_SyncGroup(ierror,cctkGH,"einstein::metric")
         call CCTK_SyncGroup(ierror,cctkGH,"einstein::curv")

c        Increment icn_iteration.

         icn_iteration = icn_iteration + 1

c     If we have done enough iterations, finish.

      else

c        Return the value of icn_iteration to 0.

         icn_iteration = 0

c        Find change in lapse.

         if (ADM_evolveLapse) then
            ADM_dtalp = alp - ADM_alp_p
         end if

      end if


c     ==================================================================
c
c     END
c
c     ==================================================================

      end subroutine ADM_IterativeCN