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

      subroutine IVP(CCTK_FARGUMENTS)

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      INTEGER CCTK_Equals

c     local vars
      CCTK_REAL temp, max_del_psi
      integer i,pass
      integer, save:: didit=0

      integer int_ivp_poststep_every
      integer reduce_handle,ierr,ivp_del_psi_index


      if(CCTK_Equals(domain, "quadrant").eq.1) then
         if(ivp_ts.eq.0) then
            call CCTK_Warn(0,31,"IVP.F","IVP", "no support for quadrant");
         endif
      endif

#ifdef BMUTIL
      write(*,*)'BMutil user, IVP does not support your thorn yet.'
      STOP
c$$$      if(CONTAINS('ivp_set_bm_vars','yes')) then
c$$$         write(*,*) 'IVP: calling BM routines metricderiv and '
c$$$         write(*,*) '     vectorini.'
c$$$         call metricderiv(METRICDERIVVARS_ARGS)
c$$$         call vectorini(VECTORINIVARS_ARGS)
c$$$      endif
#endif

c FIXME: routine vectorini_point in BMA; is this the routine you want?


c-----------------------------------------------------------------------
c adding line to check errors.

c crazy t3e hack
c cctk_iteration is integer and parameters are integer*4
c t3e compiler complains about mod function
      int_ivp_poststep_every = int(ivp_poststep_every)

      if(ivp_poststep.ne.0) then
         write(*,*)'cctk_iteration = ',cctk_iteration
         if (mod(cctk_iteration,int_ivp_poststep_every).ne.0) then
            return
         else
            if (cctk_iteration.ne.0) then
               write(*,*)'Go crazy folks!  Go crazy!!!'
c FIXME: do I want to call an exit here???
            endif
         endif
      endif

c     need some special logic so the IVP does not get solved twice
c     on iteration 0
      if((ivp_poststep.ne.0).and.(cctk_iteration.ne.0))then
         if (didit.eq.1) return
      endif

c     get conformal metric
      if(CCTK_Equals(ivp_conformal_metric, "flat").ne.0) then
         write(*,*) 'IVP: setting conformal metric to flat'
         gxx = 1.0d0
         gyy = 1.0d0
         gzz = 1.0d0
         gxy = 0.0d0
         gxz = 0.0d0
         gyz = 0.0d0
      elseif(CCTK_Equals(ivp_conformal_metric, "input").ne.0) then
         write(*,*) 'IVP: using input metric as conformal metric'
      elseif(CCTK_Equals(ivp_conformal_metric,"test").ne.0) then
         write(*,*) 'IVP: using conformal gaussian as test'
         gxx = 1.0d0 + 2.0d-1 * exp(-r**2/5.0d0)
         gyy = gxx
         gzz = gxx
         gxy = 0.0d0
         gxz = 0.0d0
         gyz = 0.0d0
      else
         write(*,*) 'IVP:error,must set parameter ivp_conformal_metric'
         write(*,*) '           to something reasonable!'
         write(*,*) 'choose flat, input, test, etc.'
         STOP
      endif

c     get conformal trK and astar
      if(CCTK_Equals(ivp_trK_conformal_astar,"zero").ne.0) then
         write(*,*) 'IVP:setting trK and astar to zero'
         trK = 0.0d0
         ivp_astarxx = 0.0d0
         ivp_astarxy = 0.0d0
         ivp_astarxz = 0.0d0
         ivp_astaryy = 0.0d0
         ivp_astaryz = 0.0d0
         ivp_astarzz = 0.0d0
      elseif(CCTK_Equals(ivp_trK_conformal_astar,"input").ne.0) then
         write(*,*) 'IVP:ivp_trK_conformal_astar = traceless part of '
         write(*,*) ' input extrinsic curvature variables (kxx, etc.)'
          if(assume_astar_divfree.eq.1) then
            write(*,*) '----------------------------------------------'
            write(*,*) 'IVP:Warning!!!!!  Yorks decomposition assumes'
            write(*,*) ' the conformal astar is divergence free!'
            write(*,*) ' If your input astar is not divergence free,'
            write(*,*) ' you will NOT solve the momentum equations!'
            write(*,*) ' To allow non-divergence-free astar, set the '
            write(*,*) ' following in your parameter file:'
            write(*,*) '   '
            write(*,*) 'assume_astar_divfree = "no"'
            write(*,*) '----------------------------------------------'
c           set astar and trK below, after call to setupIVP
         endif
      elseif(CCTK_Equals(ivp_trK_conformal_astar,"test").ne.0) then
         trK = 4.0d-1 * exp(-r**2/5.0d0)
         ivp_astarxx = 0.0d0
         ivp_astarxy = 0.0d0
         ivp_astarxz = 0.0d0
         ivp_astaryy = 0.0d0
         ivp_astaryz = 0.0d0
         ivp_astarzz = 0.0d0
      else
         write(*,*) 'IVP:error,set parameter ivp_trK_conformal_astar'
         write(*,*) '           to something reasonable!'
         write(*,*) 'choose zero, input, test, etc.'
         STOP
      endif

      if(ivp_ts.ne.0) then
         if(CCTK_Equals(ivp_trK_conformal_astar, "zero").eq.0) then
            write(*,*) 'IVP: Hey, if you want time symmetry, then you '
            write(*,*) '     had better set "ivp_trK_conformal_astar" '
            write(*,*) '     to "zero"!!!'
            STOP
         endif
         write(*,*) 'IVP: assuming Time Symmetry!!!'
         ivp_wx = 0.0d0
         ivp_wy = 0.0d0
         ivp_wz = 0.0d0
      endif

c     calculate all the stuff that does not change during the solve,
c     like scalar curvature, inverse metric
      call setupIVP(CCTK_FARGUMENTS)

c     set astar and trK for input case
      if(CCTK_Equals(ivp_trk_conformal_astar,"input").eq.1) then
         call set_astar_trk(CCTK_FARGUMENTS)
      endif

c     last ditch chance to zero out the tr(K).
      if(ivp_force_trK_zero.eq.1) then
         trK =  0.0d0
      endif

c     initial guess for solution:
      ivp_psi = 1.0d0
      ivp_wx = 0.0d0
      ivp_wy = 0.0d0
      ivp_wz = 0.0d0

      ivp_del_psi = 1.0d10
      max_del_psi = 1.0d10

      i = 0

c     Do the solve in two passes.  The first pass uses a larger
c     tolerence for the elliptic solvers
      do pass = 1,2
         if (pass.eq.1) then
            write(*,*) 'IVP:starting pass 1 w/larger stopping criteria'
            write(*,*) '     for elliptic solvers.'
         endif
         if ((pass.eq.2).and.(ivp_2pass.ne.0)) then
            write(*,*) 'IVP:starting pass 2 w/smaller stopping criteria'
            write(*,*) '     for elliptic solvers.'
            ivp_del_psi = 1.0d10
            max_del_psi = 1.0d10
         endif
         do while(max_del_psi .gt. ivp_del_psi_max)
            if((pass.eq.1).or.((pass.eq.2).and.(ivp_2pass.ne.0))) then
               i=i+1
               write(*,*)'IVP:main iteration ',i,' with max_del_psi',
     .              max_del_psi
               ivp_del_psi = ivp_psi
               call solve_iter_IVP(CCTK_FARGUMENTS,pass)
               ivp_del_psi = abs(ivp_psi - ivp_del_psi)

c     reduce ivp_del_psi
         call CCTK_ReductionHandle(reduce_handle,"maximum")

         call CCTK_VarIndex (ivp_del_psi_index, 
     &        "ivp::ivp_del_psi")

         call CCTK_Reduce(ierr, cctkGH, -1, reduce_handle,
     $        1,CCTK_VARIABLE_REAL,max_del_psi, 
     $        1,ivp_del_psi_index) 

         if (ierr.ne.0) then
            call CCTK_WARN(1,"Reduction of ivp_del_psi");
            write(*,*)'ierr = ',ierr
            STOP
         endif

c     force us out of the loop if things are time symmetric
               if(ivp_ts.eq.1) then
                  max_del_psi = 0.0d0
               endif
            endif
         enddo
      enddo
      write(*,*) 'IVP: done with iterating, max change in psi: ',
     .   max_del_psi

c     removed IO_Write1D for ivp_psi, ivp_wi.  these routines now called
c     from Output1DGH and such, now called IOASCII_Write1D.

c     now, apply the solution!
      call ivp_conformal_trans(CCTK_FARGUMENTS)

      if (mahc_evolution.ne.0) then
         call mahc_ivp_conformal(CCTK_FARGUMENTS)
      endif

c     special symmetries to metric
      call CartSymGN(ierr,cctkGH,"einstein::metric")
c     special symmetries to extrinsic curvature
      call CartSymGN(ierr,cctkGH,"einstein::curv")
      if (mahc_evolution.ne.0) then
         call CartSymGN(ierr,cctkGH,"gen_rel_perfect_fluid::Mahc_evol_vars")
         call CartSymGN(ierr,cctkGH,"gen_rel_perfect_fluid::Mahc_prim_vars")
      endif

      didit = 1

      call CCTK_Info("IVP", "Done with IVP")

      return
      end