!> @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 ParticleDeclarations
  USE CoolingSrc
  USE DMCooling
  USE ZCooling
  USE Chemistry
  USE Shapes
  USE Fields
  USE Ambients
  USE Refinements
  USE ProcessingDeclarations
  USE Totals
  USE TemplateControl 
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  INTEGER :: itracer2, itracer3, itracer4, itracer5
  REAL(KIND=qPREC) :: namb, njet, Bamb, Ramb_inner, tamb, Rjet, tjet, vjet,torusalpha,torusbeta, tt, sigmaB, MDOT_jet, rho_amb, rho_jet, jet_width, open_angle, sigmaFACTOR, tf, Bstar, Rstar, njet2, Rjet2, tjet2, vjet2, sigmaB2, rho_jet2
  REAL(KIND=qPREC) :: Bjet, Bjet2
  real(kind=qprec) :: vwind=10 !km/s
  REAL(KIND=qPREC) :: nPower_amb !!ambient density distribution r^{n}
  REAL(KIND=qPREC) :: Jet_DimTime
  REAL(KIND=qPREC), SAVE :: delta_mass, delta_mass_3D, local_mass, local_mass_3D
  REAL(KIND=qPREC), DIMENSION(3) :: pos_amb=(/0,0,1/) !Ambient center when stratified
  REAL(KIND=qPREC), DIMENSION(3) :: pos_jet=(/0,0,1/) !Jet starting location
  REAL(KIND=qPREC), DIMENSION(3) :: pos_jet2=(/0,0,1/) !Jet starting location

  REAL(KIND=qPREC) :: WedgeAngle=0d0 ! WedgeAngle for Wedge tip
 INTEGER :: buff(0:MaxDepth) !number of extra cells of buffer around sink region to set err flags

  LOGICAL :: rings, stratified, torus,lform
  LOGICAL :: lSolarMassTotal 
  INTEGER :: outflowType=2
  INTEGER :: oldframe = -1 
  REAl(KIND=qPrec) :: AGB_speed=0.0
  REAL(KIND=qPrec) :: ne=1.0, Temp ! for DMcooling, ne is not used if lTrackHydrgoen and lTrackHelium are not turned on in physics.data
  TYPE(TotalDef), POINTER :: Total
  TYPE(RefinementDef), POINTER :: Refinement
  TYPE(OutflowDef), POINTER :: Jet1, Jet2
  REAL(KIND=qPREC) :: tplace=-1d30
  LOGICAL :: lTest3DMag = .true.
  REAL :: RSurface=1, BSurface, magnetic_exponent

CONTAINS

  !> Initializes module variables
   SUBROUTINE ProblemModuleInit()      
      REAL(KIND=qPREC) :: position(3)
      REAL(KIND=qPREC) :: vel(3)
      TYPE(AmbientDef), POINTER :: Ambient
      INTEGER :: i, iErr

      NAMELIST /ProblemData/ namb, njet, Bamb, tamb, Ramb_inner,Rjet, pos_amb,&
           tjet, vjet, Bjet, &
            rings, stratified, nPower_amb, AGB_speed, &
           torus,torusalpha,torusbeta, tt,  &
           vwind, buff, jet_width, open_angle, tf, WedgeAngle, pos_jet, &
           Jet_DimTime,tjet2, vjet2, njet2, Rjet2, pos_jet2, Bjet2

      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

      tf         = tf*Pi/180d0 !convert to radians
      open_angle = open_angle*Pi/180d0 ! convert to radians
      WedgeAngle=WedgeAngle*pi/180d0

      vjet       = vjet*1e5/velscale
      rho_amb    = namb/nScale
      tamb=tamb/TempScale 
      tjet=tjet/TempScale
      Bamb=Bamb/BScale
      Bjet=Bjet/BScale
      Rjet=Rjet/lScale


      tjet2=tjet2/TempScale
      vjet2 = vjet2*1e5/velscale;
      Rjet2=Rjet2/lScale
      Bjet2=Bjet2/BScale

!      To calculate densities from mdot
!      MDOT_wind  = MDOT_wind*(mSolar/yr)/ (mScale/TimeScale) !from Msun/yr to cu
!      MDOT_jet   = MDOT_jet*(mSolar/yr)/ (mScale/TimeScale) !from Msun/yr to cu
!      rho_wind= MDOT_wind/ (4d0*Pi*Rjet**2*vwind)
!      rho_jet= MDOT_jet/ (4d0*Pi*Rjet**2*vjet)
      
!     Calculate mdot from densities
      rho_jet    = njet/nScale
      rho_jet2    = njet2/nScale


      !IF (sigmaB /=0d0) THEN
      !     IF(MPI_ID==0) print*, "sigmaB !=0, Use sigmaB to set up the jet launching zone magnetic field"
      !     Bjet=sqrt(sigmaB*rho_jet*(vjet)**2d0)

      !     !IF (MPI_ID == 0) THEN
      !     !   !write(*,*) 'BJet = ', Bjet
      !     !   write(135,*) 'sigmaB  = ', sigmaB 
      !     !   write(135,*) 'Calculated Magnetic field at launch zone surface: Bjet = ', Bjet*Bscale, " Gauss" 
      !     !END IF
      !END IF
      

      ! Create ambient and setup tracers/diagnostics vars.
      CALL CreateAmbient(Ambient)
      Ambient%density=rho_amb
      !CALL AddDiagnosticVar(CoolingStrength_Field)
      !CALL AddTracer(itracer4, 'Jet tracer')
      !CALL AddTracer(itracer5, 'DensityFlag tracer')
     
     !  Create refinement around base of jet
      !CALL CreateRefinement(Refinement)
      !CALL CreateShape(Refinement%Shape, RECTANGULAR_PRISM, (/Rjet,Rjet,Rjet/))
      !CALL CreateShape(Refinement%Shape, RECTANGULAR_PRISM, (/Rjet2,Rjet2,Rjet2/))
      !Refinement%BufferCells=buff(0)

      
      CALL CreateOutflow(Jet1)
      !Jet1%duration=20
      !Jet1%base=0.5*Rjet !cylindrical radius of jet
      !Jet1%position=(/0,0,0/) !pos_jet
      Jet1%density=rho_jet
      Jet1%radius=Rjet!levels(MaxLevel)%dx*36
      Jet1%thickness=0.5*Rjet!levels(MaxLevel)%dx*32
      Jet1%temperature=tjet
      !Jet1%decay=10.0
      !Jet1%open_angle=open_angle
      Jet1%theta=open_angle
      Jet1%phi=0.5*pi !d0
      Jet1%velocity=(/0d0,0d0,vjet/)
      !Jet1%speed=vjet
      Jet1%B=Bjet
      Jet1%begin=0d0
      CALL UpdateOutflow(Jet1)

      !nullify(Jet1);
      !Jet1%position=pos_jet2
      !Jet1%open_angle=half*Pi
      !Jet1%radius=Rjet2
      !Jet1%thickness=Rjet2
      !Jet1%begin=-1d0
      !Jet1%velocity=vjet2
      !Jet1%duration=200000
      !!Jet1%decay=outflow_decay
      !Jet1%density=rho_jet2
      !Jet1%temperature=tjet2
      !CALL UpdateOutflow(Jet1)
      CALL CreateOutflow(Jet2)
      !Jet2%duration=20000
      Jet2%position=pos_jet2
      Jet2%density=rho_jet2
      Jet2%radius=Rjet2
      Jet2%thickness=0.5*Rjet2
      !Jet2%base=Rjet2 !cylindrical radius of jet
      Jet2%temperature=tjet2
      !Jet2%open_angle=open_angle
      Jet1%theta=pi-open_angle  ! open angle from z-axis
      Jet1%phi=0.5*pi  !open angle around z axis from x
      Jet2%velocity=(/0d0,0d0,vjet2/)
      Jet2%B=Bjet2
      !Jet2%speed=vjet2
      Jet2%begin=0d0
      CALL UpdateOutflow(Jet2)

  !    CALL SetupDiagnostics

   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, rmbc, zrmbc, mx, my, mz, level
      REAL(KIND=qPREC), POINTER, DIMENSION(:,:,:,:) :: q
      REAL(KIND=qPREC) :: nH, nHe, SIIFrac

      if (lrestart) return

      IF (iCooling == 4) THEN
         level = Info%level
         q=>Info%q
         rmbc = levels(level)%gmbc(levels(level)%step)
         zrmbc = 0
         mx = Info%mX(1)
         my = Info%mX(2)
         mz = Info%mX(3)
         IF (nDim == 3) zrmbc = rmbc
         SIIFrac = 10d0**(7.2d0 - 12d0)

         DO i=1-rmbc, mx+rmbc
            DO j=1-rmbc, my+rmbc
               DO k=1-zrmbc, mz+zrmbc
                  IF (lTrackHydrogen) THEN
                     nH = q(i,j,k,1)/(1d0+xH+HeFrac*(1d0+xHe+2d0*xHeII)+ZFrac*(1d0+xZ))
                     q(i,j,k,iH) = (1d0-xH)*nH
                     q(i,j,k,iHII) = xH*nH
                     IF (lTrackHelium) THEN
                        nHe = HeFrac*nH
                        q(i,j,k,iHe) = (1d0-xHe-xHeII)*nHe
                        q(i,j,k,iHeII) = xHe*nHe
                        q(i,j,k,iHeIII) = xHeII*nHe
                     END IF
                     IF (lTrackSulfur) THEN
                        q(i,j,k,iSII) = SIIFrac*nH
                        q(i,j,k,iSIII) = 0d0
                        q(i,j,k,iSIV) = 0d0
                     END IF
                  END IF
               END DO
            END DO
         END DO
      END IF

     ! CALL ApplyFunction(SetupAmbient, Info)
    !  CALL ApplyFunction(SetupJet, Info)        
    !  CALL ApplyFunction(SetupJet2, Info)        
    !  IF (nDim == 3 .AND. lMHD) THEN
    !     CALL ApplyEMF(JetInitEmf, Info)
    !     CALL ApplyEMF(JetInitEmf2, Info)
    !  END IF
!      IF (nDim == 3 .AND. lMHD) CALL ApplyBField(JetInitEmf, Info)
!      write(*,*) 'gridinit applied', Info%q(1,9,1,1)
    END SUBROUTINE ProblemGridInit

!> Applies Boundary conditions
!! @param Info Info object
SUBROUTINE ProblemBeforeStep(Info)
      TYPE(InfoDef) :: Info
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
return
END SUBROUTINE ProblemSetErrFlag

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

END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem

