!#########################################################################
!		
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    resistive.f90 is part of AstroBEAR.
!
!    AstroBEAR is free software: you can redistribute it and/or modify	  
!    it under the terms of the GNU General Public License as published by 
!    the Free Software Foundation, either version 3 of the License, or    
!    (at your option) any later version.
!
!    AstroBEAR is distributed in the hope that it will be useful, 
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with AstroBEAR.  If not, see <http://www.gnu.org/licenses/>.
!
!#########################################################################
MODULE ResistiveSrc
  USE GlobalDeclarations
  USE TreeDeclarations
  USE DataDeclarations
  USE ExplicitDeclarations
  USE PhysicsDeclarations
  USE EOS

CONTAINS

  SUBROUTINE ResistiveGridAdvance(n,stepping)
    INTEGER :: n,stepping
    TYPE(NodeDef), POINTER :: node
    TYPE(NodeDefList), POINTER :: nodelist        
    !alternative solution: REAL ALLOCATABLE :: q,aux
    IF(MPI_ID==0)   PRINT*, "Advancing Resistive Solver on Grid..."
    nodelist=>nodes(n)%p
    DO WHILE (associated(nodelist))
       node=>nodelist%self
       CALL ResistiveSource(n,node%info,stepping)
       !alternative solution: store q and aux
       nodelist=>nodelist%next
    END DO

  END SUBROUTINE ResistiveGridAdvance

  SUBROUTINE ResistiveGridGetMaxSpeed(n)
    INTEGER :: n
    TYPE(NodeDef), POINTER :: node
    TYPE(NodeDefList), POINTER :: nodelist
    
    nodelist=>Nodes(n)%p 
    DO WHILE (associated(nodelist))
       node=>nodelist%self 
       CALL ResistiveGetMaxSpeed(node%info)
       nodelist=>nodelist%next 
    END DO

  END SUBROUTINE ResistiveGridGetMaxSpeed

  SUBROUTINE ResistiveSource(n,Info,stepping)
    TYPE(InfoDef) :: Info
    LOGICAL :: consist_en=.true.
    INTEGER :: n
    INTEGER :: i,j,k,z_off,stepping,mx,my,mz,mbc,zmbc
    INTEGER :: msx(3,2),msy(3,2),msz(3,2), mjx(3,2), mjy(3,2), mjz(3,2)
    INTEGER, ALLOCATABLE, DIMENSION(:) :: iflux,bflux
    REAL(KIND=qPREC) :: dx,dt,check_en,off_en
    REAL(KIND=qPREC), ALLOCATABLE, DIMENSION(:) :: new_b
    REAL(KIND=qPREC), ALLOCATABLE, DIMENSION(:,:,:) :: jx, jy, jz
    REAL(KIND=qPREC), ALLOCATABLE, DIMENSION(:,:,:,:) :: xflux, yflux, zflux
    REAL(KIND=qPREC), ALLOCATABLE, DIMENSION(:,:,:,:) :: fjx, fjy

    ! bounds for edge centered current jx, jy, jz
    mjx=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc,my+mbc+1,mz+zmbc+z_off/),(/3,2/))
    mjy=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc+1,my+mbc,mz+zmbc+z_off/),(/3,2/))
    mjz=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc+1,my+mbc+1,mz+zmbc/),(/3,2/))

    ! bounds for face centered energy fluxes
    msx=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc+1,my+mbc,mz+zmbc/),(/3,2/))
    msy=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc,my+mbc+1,mz+zmbc/),(/3,2/))
    msz=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc,my+mbc,mz+zmbc+z_off/),(/3,2/))

    ALLOCATE(iflux(1))
    iflux(:)=iE

    !allocate memory for cpuJx, cpuJy, cpuJz, cpuxflux, cpuyflux, cpuzflux, cpuxaux, cpuyaux, cpuzaux
    !assign Info%mx array to mx array
    !assign Info%xBounds array to xBounds array
    CALL cudaCallKernel(cpuJx, cpuJy, cpuJz, cpuxflux, cpuyflux, cpuzflux, cpuxaux, cpuyaux, cpuzaux, mx, xBounds, dx, resistivity)
    !copy cpuxaux, cpuyaux, cpuzaux into Info%aux(:,:,:,1:3)

    !store emf fluxes
    CALL StoreEmfs(Info,mjx,1,cpuJx*dt/dx)
    CALL StoreEmfs(Info,mjy,2,cpuJy*dt/dx)
    CALL StoreEmfs(Info,mjz,2,cpuJz*dt/dx)

    ! store energy flux for synchronizing AMR levels
    CALL storefixupfluxes(Info,msx,1,cpuxflux,iflux)
    CALL storefixupfluxes(Info,msy,2,cpuyflux,iflux)
    CALL storefixupfluxes(Info,msz,3,cpuzflux,iflux)

    DEALLOCATE(iflux)!; DEALLOCATE(new_b)

  END SUBROUTINE ResistiveSource

 SUBROUTINE ResistiveGetMaxSpeed(Info)
    TYPE(InfoDef) :: Info
    INTEGER :: level, m
    INTEGER :: mjx(3,2), mjy(3,2), mjz(3,2)
    INTEGER :: z_off,mx,my,mz,mbc,zmbc
    REAL(KIND=qPREC) :: dx, maxres, tmp
   
    maxres = 0d0
    level=Info%level
    dx=levels(level)%dx
    mx=Info%mX(1);my=Info%mX(2);mz=Info%mX(3)
    mbc=levels(n)%ambc(levels(n)%step)-1  
    SELECT CASE(nDim)
    CASE(2)
       zmbc=0;z_off=0
    CASE(3)
       zmbc=mbc;z_off=1
    END SELECT
   

    mjx=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc,my+mbc+1,mz+zmbc+z_off/),(/3,2/))
    mjy=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc+1,my+mbc,mz+zmbc+z_off/),(/3,2/))
    mjz=reshape((/1-mbc,1-mbc,1-zmbc,mx+mbc+1,my+mbc+1,mz+zmbc/),(/3,2/))
    IF(ResType==2)THEN
       DO i=mjx(1,1),mjx(1,2); DO j=mjx(2,1),mjx(2,2); DO k=mjx(3,1),mjx(3,2)
          CALL RealResistivity(Info,i,j,k,tmp,dx,1)
          maxres = max(maxres, tmp);
       END DO; END DO; END DO
       DO i=mjy(1,1),mjy(1,2); DO j=mjy(2,1),mjy(2,2); DO k=mjy(3,1),mjy(3,2)
          CALL RealResistivity(Info,i,j,k,tmp,dx,2)
          maxres = max(maxres, tmp);
       END DO; END DO; END DO
       DO i=mjz(1,1),mjz(1,2); DO j=mjz(2,1),mjz(2,2); DO k=mjz(3,1),mjz(3,2)
          CALL RealResistivity(Info,i,j,k,tmp,dx,3)
          maxres = max(maxres, tmp);
       END DO; END DO; END DO
    END IF
    
    explicit_maxspeed(level) = max(explicit_maxspeed(level),2d0*maxres/dx)

  END SUBROUTINE ResistiveGetMaxSpeed


  SUBROUTINE getT(Info,i,j,k,temp)
     TYPE(InfoDef) :: Info
     INTEGER :: i,j,k
     REAL(KIND=qPREC) :: temp
  
     temp=(gamma-1d0)*(Info%q(i,j,k,iE)-0.5*SUM(Info%q(i,j,k,2:4)**2)/Info%q(i,j,k,1)&
      -0.5*SUM(Info%q(i,j,k,iBx:iBz)**2))/Info%q(i,j,k,1)
  
  END SUBROUTINE getT
  
  
END MODULE ResistiveSrc



