!> @dir OrbitingParticles
!! @brief Contains files necessary for the Orbiting Particles test problem

!> @file problem.f90
!! @brief Main file for module Problem

!> @defgroup OrbitingParticles Orbiting Particle Module
!! @brief Module for setting up orbiting particles
!! @ingroup Modules

!> Orbiting Particle Module
!! @ingroup OrbitingParticles
MODULE Problem
  USE CommonFunctions
  USE Outflows
  USE DataDeclarations
  USE CoolingSrc
!  USE AnalyticCooling2
  USE DMCooling
  USE ZCooling
  USE Chemistry
  USE Shapes
  USE Fields
  USE Ambients
  USE Refinements
  USE ProcessingDeclarations
  USE Totals
  USE TemplateControl 
  USE Winds
  USE Clumps
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  INTEGER :: itracer2, itracer3, itracer4, itracer5
  REAL(KIND=qPREC) :: namb, nwind, Ramb_inner, tamb, twind, vwind, beta, mach
  REAL(KIND=qPREC) :: nclump,tclump, Rclump
  TYPE(ClumpDef), POINTER :: clump1
  TYPE(ClumpDef), POINTER :: clump2
  TYPE(WindDef), POINTER :: Wind
  REAL(KIND=qPREC), DIMENSION(3) :: pos_clump1=(/0,0,1/) 
  REAL(KIND=qPREC), DIMENSION(3) :: pos_clump2=(/0,0,1/) 

CONTAINS

  !> Initializes module variables
   SUBROUTINE ProblemModuleInit()      
      TYPE(AmbientDef), POINTER :: Ambient
      INTEGER :: i, iErr

      NAMELIST /ProblemData/ namb, nwind, tamb, & 
          twind, vwind, beta, mach,&
          nclump, tclump, Rclump,&
          pos_clump1, pos_clump2

      OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
      READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
      CLOSE(PROBLEM_DATA_HANDLE)      


      !Convert input parameters to computational units
      namb    = namb/nScale
      tamb=tamb/TempScale 
      nwind    = nwind/nScale
      twind=twind/TempScale

      nclump = nclump/nScale
      Rclump = Rclump/lScale
      tclump = tclump/TempScale
      ! Create ambient and setup tracers/diagnostics vars.
      CALL CreateAmbient(Ambient)
      !CALL AddTracer(Ambient%iTracer,'Ambient')
      Ambient%density=namb
      Ambient%pressure=namb*tamb
      !CALL AddDiagnosticVar(CoolingStrength_Field)
      !CALL AddTracer(itracer4, 'Jet tracer')
      !CALL AddTracer(itracer5, 'DensityFlag tracer')
     
    CALL CreateWind(Wind)
    wind%density=nwind
    Wind%velocity=mach*sqrt(gamma*nwind*twind)
    Wind%temperature=twind
    Wind%B=sqrt(2d0/beta)*(/0d0,1d0,0d0/)


    CALL CreateClump(clump1)
    clump1%position=pos_clump1 
    clump1%density=nclump
    clump1%radius=Rclump
    clump1%temperature=tclump
    clump1%density_profile=1 ! PN_Profile, soft clump density
    clump1%BG_density=namb
    CALL UpdateClump(clump1)
      
    CALL CreateClump(clump2)
    clump2%position=pos_clump2 
    clump2%density=nclump
      clump2%radius=Rclump
      clump2%temperature=tclump
      clump2%density_profile=1 ! PN_Profile, soft clump density
      clump2%BG_density=namb
      CALL UpdateClump(clump2)
      

   END SUBROUTINE ProblemModuleInit



   SUBROUTINE SetupDiagnostics

   end SUBROUTINE SetupDiagnostics

  !> Applies initial conditions
  !! @param Info Info object
   SUBROUTINE ProblemGridInit(Info)
      TYPE(InfoDef) :: Info
      INTEGER :: i,j,k,l,m,ii,jj,kk
      INTEGER :: rmbc,zrmbc,level
      INTEGER :: mx, my, mz
      INTEGER :: iErr

      return
   END SUBROUTINE ProblemGridInit

!> Applies Boundary conditions
!! @param Info Info object
SUBROUTINE ProblemBeforeStep(Info)
      TYPE(InfoDef) :: Info
      INTEGER :: i,j,k,l,m,ii,jj,kk
      INTEGER :: rmbc,zrmbc,level
      INTEGER :: mx, my, mz
      INTEGER :: iErr
      LOGICAL, SAVE :: once=.true.
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:)     :: q, aux
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: density
      REAL (KIND=qPrec), SAVE :: dxFIXED
      REAL(KIND=xprec) :: x,y,z,x2,y2,z2,xl,yl,zl,dx,dy,dz,r,r2,dr,dt,r3,torus2,cos2theta,ex1,ex2,Ftheta,f,dec7, bmax
      LOGICAL, DIMENSION(:,:,:),ALLOCATABLE::  density_flag

      return
      if (lRestart) return

      level=Info%level
      q=>Info%q
      aux=>Info%aux

      ! Calculating the number of ghost cells on each side of the grid.
      rmbc=levels(level)%gmbc(levels(level)%step)
      mx = Info%mX(1)
      my = Info%mX(2)
      mz = Info%mX(3)
      dx=levels(level)%dX
      dy=dx
      dz=dx
      xl=Info%xBounds(1,1)
      yl=Info%xBounds(2,1)
      zl=Info%xBounds(3,1)
      dt=levels(level)%dt
      SELECT CASE(nDim)
      CASE(2)
         zrmbc=0;mz=1;dz=0d0;zl=0d0
      CASE(3)
         zrmbc=rmbc
      END SELECT

      DO i=1-rmbc, mx+rmbc         ; x = (xl+(REAL(i,xPrec)-half)*dx)
         DO j=1-rmbc, my+rmbc      ; y = (yl+(REAL(j,xPrec)-half)*dy)
            DO k=1-zrmbc, mz+zrmbc ; z = (zl+(REAL(k,xPrec)-half)*dz)
              if(y<=GxBounds(2,1)+2d0*dy) then
                Info%q(i,j,k,1) = nwind;
                Info%q(i,j,k,iVx) = nwind*vwind;  !!vy
                Info%q(i,j,k,iVy) = 0d0;  !!vx
                Info%q(i,j,k,iE) = nwind * vwind *gamma7
                Info%q(i,j,k,iBy) =  sqrt(2d0/beta)
              end if
          END DO
        END DO
      END DO

END SUBROUTINE ProblemBeforeStep

!> Could be used to update grids pre-output
!! @param Info Info Object
SUBROUTINE ProblemAfterStep(Info)
TYPE(InfoDef) :: Info
END SUBROUTINE ProblemAfterStep

!> Could be used to set force refinement
!! @param Info Info object
SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:)     :: q
    INTEGER :: i,j,k,mx,my,mz,rmbc,zrmbc,level
    REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,lowerZ,upperZ,centerZ

    return

    level=Info%level
    rmbc=levels(level)%gmbc(levels(level)%step)
    zrmbc=rmbc
    mx = Info%mX(1)
    my = Info%mX(2)
    mz = Info%mX(3)
    dx=levels(level)%dX
    dy=dx
    dz=dx
    !xl=Info%xBounds(1,1)
    !yl=Info%xBounds(2,1)
    zl=Info%xBounds(3,1) !! local
    lowerZ=GxBounds(3,1)  !! global
    upperZ=GxBounds(3,2)  !! global
    centerZ=(upperZ-lowerZ)/2d0


END SUBROUTINE ProblemSetErrFlag

SUBROUTINE ProblemBeforeGlobalStep(n)
  INTEGER :: n
  INTEGER :: ierr

END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem

