MODULE MTI

  USE NodeInfoDef

  USE GlobalDeclarations
  IMPLICIT NONE    
  SAVE
  PRIVATE
 
  PUBLIC :: scaleMTI,qinitMTI,b4stepMTI
  INTEGER :: MTItype, nmodex, nmodey, bc
  REAL(KIND=qPrec) :: rho0, p0, y0, v0, B0, B_static, lx
  NAMELIST/MTIdata/MTItype, nmodex, nmodey, rho0, p0, y0, v0, B0, B_static, lx, bc

CONTAINS

  SUBROUTINE scaleMTI
    ! Reads in the MTI data file
    INTEGER :: iErr
    OPEN(UNIT=10,IOSTAT=iErr,ERR=50,FILE='MTI.data',STATUS='old',      &
         FORM='formatted')
50  IF (iErr /= 0) THEN
       PRINT *,'RUNNING MTI: ERROR opening MTI.data!'
       STOP
    END IF
    READ(10,NML=MTIdata,IOSTAT=iErr,ERR=51)
    CLOSE(10)
51  IF (iErr /= 0) THEN
       PRINT *,'RUNNING MTI: ERROR reading from MTI.data! Please check the data file'
       STOP
    END IF
  END SUBROUTINE scaleMTI

  SUBROUTINE qinitMTI(Info)
    ! Interface declarations
    TYPE (NodeInfo) :: Info  ! Data associated with this grid
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:,:)     :: q, aux
    INTEGER :: mx,my,mz,mbc,rmbc,zrmbc
    INTEGER :: i, j, k
    REAL(KIND=xPrec) :: x, x0,dx, dy, dz, xpos, ypos, zpos, &
         xl, yl, zl, r, r0, T, t0
    REAL(KIND=qPrec) :: rho, vy, p, B, bbx, bby, phi, phix, phiy
    REAL(KIND=qPrec), DIMENSION(:,:), ALLOCATABLE :: VectorPotential
    CHARACTER :: MTIname

    q=>Info%q
    aux=>Info%aux
    mbc=Info%mbc
    IF(Info%ErrorEstimation) THEN
       rmbc=(Info%AMRSteps-Info%AMRStep+1)*mbc
    ELSE
       rmbc=Info%r*mbc
    END IF


    mx=Info%mX(1); my=Info%mX(2); mz=Info%mX(3)
    dx=Info%dX(1); dy=Info%dX(2); dz=Info%dX(3)
    xl=Info%Xlower(1); yl=Info%Xlower(2); zl=Info%Xlower(3)
    r0 = lx/1.0
    x0 = 0.05
    t0 = 1.0e+2
    
    SELECT CASE(Info%nDim)
    CASE(2)
       zrmbc=0;mz=1;zl=0;dz=0
    CASE(3)
       zrmbc=rmbc
    END SELECT
    ALLOCATE(VectorPotential(1-rmbc:mx+rmbc,1-rmbc:my+rmbc))
    q(:,:,:,:,:) = 0d0; aux(:,:,:,:,:) = 0d0
    VectorPotential = 0d0
    
    SELECT CASE(MTItype)
    CASE(1) ! Simple 2D Vy Perturbation
       MTIname = 'Simple 2D Vy'
       DO j = 1-rmbc, my+rmbc
          DO i = 1-rmbc, mx+rmbc
             xpos = xl+(real(i,xPrec)-0.5)*dx
             ypos = yl+(real(j,xPrec)-0.5)*dy
             rho = rho0*(1d0-ypos/y0)**2
             p = p0*(1d0-ypos/y0)**3
             vy = v0*sin(real(nmodex)*pi*xpos/lx)
             q(i,j,1,1,1) = rho
             q(i,j,1,1,3) = rho*vy
             q(i,j,1,1,iBx) = B0
             aux(i,j,1,1,1) = B0
             q(i,j,1,1,iE) = p/(gamma-1d0)+0.5*rho*vy**2+0.5*B0**2
          END DO
       END DO
       PRINT*, 'p0 = ', p0
       PRINT*, 'rho0 = ', rho0
       PRINT*, 'temp = ', p0/rho0
       PRINT*, 'mCentral = ', mCentral
    CASE(2) ! Empty Grid MHD Conduction for test.
       MTIname = 'Empty Grid Test'
       DO k = 1-zrmbc, mz+zrmbc
       DO j = 1-rmbc, my+rmbc
       DO i = 1-rmbc, mx+rmbc
          xpos = xl+(real(i,xPrec)-0.5)*dx
          ypos = yl+(real(j,xPrec)-0.5)*dy
          r = xpos**2+ypos**2
          IF(nDim==3)THEN
             zpos = zl+(real(k,xPrec)-0.5)*dz
             r = r+zpos**2
          END IF
          r = SQRT(r)
          rho = exp(r**2/r0**2)
          p = 1.0
          q(i,j,k,1,1) = rho
          q(i,j,k,1,iBx) = B0
          aux(i,j,k,1,1) = B0
          q(i,j,k,1,iE) = p/(gamma-1d0)+0.5*B0**2
       END DO
       END DO
       END DO
    CASE(3) ! Sine B Field Diffusion.
       MTIname = 'Sine Field Test'
       PRINT*, "B0*lx = ", B0*lx
       DO j = 1-rmbc, my+rmbc
          DO i = 1-rmbc, mx+rmbc
             xpos = xl+(real(i,xPrec)-0.5)*dx
             ypos = yl+(real(j,xPrec)-0.5)*dy
             IF((xpos-x0).le.-5e-5)THEN
                T = t0
                IF(xpos.gt.0d0) T = t0*(1-xpos**2/x0**2)**(1/2.5)
             ELSE
                T = 8.3240004
             END IF
             p = p0
             rho = p/T
             q(i,j,1,1,1) = rho
             VectorPotential(i,j)=0d0
             IF(xpos.gt.-0.05.and.xpos.lt.0.15)THEN
                phi = pi*(xpos-0.5*lx)/lx
                phix= 8.0*pi*xpos/lx; phiy = 8.0*pi*ypos/lx
                VectorPotential(i,j)=(B0*lx/(16d0*pi))*(cos(phix)-cos(phiy))*(1d0+cos(phi))
             END IF
          END DO
       END DO
       DO j = 0,my+1
          DO i = 0,mx+1
             bbx = (VectorPotential(i,j+1)-VectorPotential(i,j))/dy
             bby = -(VectorPotential(i+1,j)-VectorPotential(i,j))/dx
             aux(i,j,1,1,1) = bbx; aux(i,j,1,1,2) = bby
             Info%q(i,j,1,1,iBx) = bbx; Info%q(i,j,1,1,iBy) = bby
             q(i,j,1,1,iE) = p/(gamma-1d0)+0.5*(SUM(Info%q(i,j,1,1,iBx:iBy)**2))
          END DO
       END DO
       
    CASE DEFAULT
       MTIname = 'UNKNOWN'
       PRINT*, 'RUNNING MTI: UNKNOWN MTI type. Please check the data file.'
       STOP
    END SELECT
    DEALLOCATE(VectorPotential)
    
    IF((Info%MPI_id == 0).AND.(Info%level == 0))THEN
       PRINT*
       PRINT*, '=================================================='
       PRINT*, 'MTI PROBLEM INITIALIZED:'
       PRINT*, 'PERTURBATION TYPE:', MTIname, 'Perturbation'
       PRINT*, 'MODE NUMBER: ', nmodex, nmodey
       PRINT*, 'BOTTOM DENSITY: ', rho0
       PRINT*, 'BOTTOM PRESSURE: ', p0
       PRINT*, 'FIELD STRENGTH: ', B0**2
       PRINT*, 'SCALE HIGHT: ', y0
       PRINT*, 'ACCELERATION: ', mCentral
       PRINT*, '=================================================='
       PRINT*
    END IF
  END SUBROUTINE qinitMTI
  
  SUBROUTINE b4stepMTI(Info)
    ! Interface declarations
    TYPE (NodeInfo) :: Info  ! Data associated with this grid
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:,:)     :: q,aux
    INTEGER :: i,j,k,mx,my,mz,mbc,rmbc,zrmbc, ibc
    REAL (KIND=xprec) :: dx, dy
    REAL (KIND=qPrec), ALLOCATABLE, DIMENSION(:) :: ptemp

    q=>Info%q

    aux=>Info%aux
    mbc=Info%mbc
    IF(Info%ErrorEstimation) THEN
       rmbc=(Info%AMRSteps-Info%AMRStep+1)*mbc
    ELSE
       rmbc=Info%r*mbc
    END IF

    mx=Info%mX(1); my=Info%mX(2); mz=Info%mX(3)
    dx=Info%dX(1); dy=Info%dX(2)

    SELECT CASE(bc)
    CASE(0)
       CONTINUE
    CASE(1) ! fixed boundary.
       IF(Info%mthbc(3)==3)THEN
          ALLOCATE(ptemp(1-rmbc:mx+rmbc))          
          DO ibc=1,rmbc
             ptemp(:) = p0*(1d0-(real(1-ibc)-0.5)*dy/y0)**3
             q(:,1-ibc,1,1,iE) = ptemp(:)/(gamma-1d0)+0.5*(q(:,1-ibc,1,1,2)**2+&
                  q(:,1-ibc,1,1,3)**2)/q(:,1-ibc,1,1,1)
          END DO
          DEALLOCATE(ptemp)
       END IF

       IF(Info%mthbc(4)==3)THEN
          ALLOCATE(ptemp(1-rmbc:mx+rmbc))
          DO ibc=1,rmbc
             ptemp(:) = p0*(1d0-(real(my+ibc)-0.5)*dy/y0)**3
             q(:,my+ibc,1,1,iE) = ptemp(:)/(gamma-1d0)+0.5*(q(:,my+ibc,1,1,2)**2+&
                  q(:,my+ibc,1,1,3)**2)/q(:,my+ibc,1,1,1)
          END DO
          DEALLOCATE(ptemp)
       END IF
       
    CASE(2) ! extrapolation
       IF(Info%mthbc(3)==3)THEN
          ALLOCATE(ptemp(1-rmbc:mx+rmbc))
          ptemp(:) = (gamma-1d0)*(q(:,1,1,1,iE)-0.5*(SUM(q(:,1,1,1,2:3)**2))/q(:,1,1,1,1)-0.5*(SUM(q(:,1,1,1,iBx:iBy))**2)) ! obtain the bottom line pressure.
          DO ibc=1,rmbc
             ptemp(:) = ptemp(:)+dy*0.5*(q(:,2-ibc,1,1,1)+q(:,1-ibc,1,1,1))*mCentral ! increase by the amount of rho*g*dy
             q(:,1-ibc,1,1,iE) = ptemp(:)/(gamma-1d0)+0.5*(SUM(q(:,1-ibc,1,1,2:3)**2))/q(:,1-ibc,1,1,1)+0.5*(SUM(q(:,1,1,1,iBx:iBy))**2)
          END DO
          DEALLOCATE(ptemp)
       END IF

       IF(Info%mthbc(4)==3)THEN
          ALLOCATE(ptemp(1-rmbc:mx+rmbc))
          ptemp(:) = (gamma-1d0)*(q(:,my,1,1,iE)-0.5*(SUM(q(:,my,1,1,2:3)**2))/q(:,my,1,1,1)-0.5*(SUM(q(:,my,1,1,iBx:iBy))**2)) ! obtain the top line pressure.
          DO ibc=1,rmbc
             ptemp(:) = ptemp(:)-dy*0.5*(q(:,my+ibc-1,1,1,1)+q(:,my+ibc,1,1,1))*mCentral ! decrease by the amount of rho*g*dy
             q(:,my+ibc,1,1,iE) = ptemp(:)/(gamma-1d0)+0.5*(SUM(q(:,my+ibc,1,1,2:3)**2))/q(:,my+ibc,1,1,1)+0.5*(SUM(q(:,my,1,1,iBx:iBy))**2)
          END DO
          DEALLOCATE(ptemp)
       END IF
    CASE(3) ! Reservoire on Two Ends
       IF(Info%mthbc(1)==1)THEN
          q(1-rmbc:0,:,:,1,1) = p0/100.000000
          q(1-rmbc:0,:,:,1,2:4) = 0.0
          q(1-rmbc:0,:,:,1,iE) = p0/(gamma-1.0)+0.5*B0**2
          !PRINT*
          !PRINT*, "ghost and edge ie:", q(0,8,1,1,iE), q(1,8,1,1,iE)
          !PRINT*, "ghost and edge px:", q(0,8,1,1,2), q(1,8,1,1,2)
          !PRINT*
       END IF
       IF(Info%mthbc(2)==1)THEN
          q(mx+1:mx+rmbc,:,:,1,1) = p0/8.3240004
          q(mx+1:mx+rmbc,:,:,1,2:4) = 0.0
          q(mx+1:mx+rmbc,:,:,1,iE) = p0/(gamma-1.0)+0.5*B0**2
       END IF
    CASE DEFAULT
       PRINT*, 'UNKNOWN BOUNDARY CONDITION IN B4STEP. HALTING... '
       STOP
    END SELECT

  END SUBROUTINE b4stepMTI
  
END MODULE MTI
