!> @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 Fields
  USE Ambients
  USE Refinements
  USE ProcessingDeclarations
  USE Clumps
  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
  REAL(KIND=qPREC) :: Bjet, Bclump
  real(kind=qprec) :: vwind=10 !km/s
  REAL(KIND=qPREC) :: nPower_amb !!ambient density distribution r^{n}
  REAL(KIND=qPREC) :: nClump, RClump, tClump, vClump, rho_clump, sigma_Clump, sigmaFACTOR_Clump, thicknessClump=.1 
  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_clump=(/0,0,1/) !Clump starting location 
  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

 INTEGER :: buff(0:MaxDepth) !number of extra cells of buffer around sink region to set err flags

  LOGICAL :: rings, stratified, torus,lform
  TYPE(ParticleDef), POINTER :: Particle
  INTEGER, PARAMETER :: collimated=1, conical=2, clump=3
  INTEGER :: outflowType=2
  INTEGER :: oldframe = -1 
  REAl(KIND=qPrec) :: AGB_speed
  REAL(KIND=qPrec) :: Clump_StartTime, Jet_OffTime
  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(ClumpDef), POINTER :: myClump
  TYPE(TotalDef), POINTER :: Total
  LOGICAL :: lClumpAdded = .false., lPlaceClump=.false., lRefineClump=.false., lPlaceClumpInBoundaries=.false., lprintTotal=.false.
  TYPE(RefinementDef), POINTER :: Refinement, ClumpRefinement
  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

      NAMELIST /ProblemData/ namb, njet, Bamb, tamb, Ramb_inner,Rjet, pos_amb,&
           tjet, vjet, nClump, RClump, tClump, vClump, Bstar, magnetic_exponent, &
            rings, stratified, nPower_amb, AGB_speed, &
           torus,torusalpha,torusbeta, tt, sigmaB, sigma_Clump, MDOT_jet, &
           vwind, buff, jet_width, open_angle, tf, pos_jet, pos_clump, &
           Clump_StartTime, Jet_OffTime, Jet_DimTime,lprintTotal

      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
      vwind      = vwind*1e5/VelScale !from km/s to cu
      vjet       = vjet*1e5/velscale
      vclump     = vclump*1e5/velscale
      AGB_speed = AGB_speed*1e5/VelScale ! convert AGB_speed to comp unit
      rho_amb    = namb/nScale
      tamb=tamb/TempScale 
      tjet=tjet/TempScale
      Bamb=Bamb/BScale
      Bstar=Bstar/BScale
      Rjet=Rjet*au/lScale

!      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_clump  = nclump/nScale
      !rho_wind   = nwind/nScale
      !MDOT_wind = 4d0*Pi*Rjet**2*rho_wind*vwind
      MDOT_jet =  4d0*Pi*Rjet**2*rho_jet*vjet

      IF (MPI_ID == 0) THEN
      !   write(*,*) 'MDOT_wind = ', MDOT_wind*(mScale/TimeScale)/(mSolar/yr) !from cu to Msun/yr 
         write(*,*) 'MDOT_jet = ', MDOT_jet * (mScale/TimeScale)/(mSolar/yr) !from cu to Msun/yr
      end IF


      IF (Bstar /= 0d0) THEN
         Bjet = Bstar/(Rjet**magnetic_exponent)
      ELSE
         Bjet=0d0
      END IF

     ! Bclump=sqrt(sigma_Clump*rho_clump*vClump**2d0)
      IF (MPI_ID == 0) THEN
         write(*,*) 'BJet = ', Bjet
      END IF

      
      ! Create ambient and setup tracers/diagnostics vars.
      CALL CreateAmbient(Ambient)
      CALL AddDiagnosticVar(CoolingStrength_Field)
      CALL AddTracer(itracer3, 'Clump tracer')
      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/))
      Refinement%BufferCells=buff(0)

      ! Create clump for model with clump
      CALL CreateClump(myClump)
      myClump%position=pos_Clump!(/0d0,-1.0*(vClump*Clump_Starttime/velscale),0d0/)
      myClump%density=nClump/nscale
      myClump%radius=RClump
      myClump%temperature=tClump/tempScale
      myClump%density_profile=2 ! PN_Profile, soft clump density
      myClump%BG_density=namb/nScale
      CALL UpdateClump(myClump)

      CALL SetupDiagnostics

   END SUBROUTINE ProblemModuleInit



   SUBROUTINE SetupDiagnostics

     !!=== print out the total Mass and Momentum
     if(lprintTotal) then
        lSolarMassTotal=.true.
        CALL CreateTotal(Total)
        Total%Field%id=Mass_Field
        Total%Field%component=GASCOMP
        Total%Field%name='Mass'

        CALL CreateTotal(Total)
        Total%Field%id=Px_Field
        Total%Field%component=GASCOMP
        Total%Field%name='Px'

        CALL CreateTotal(Total)
        Total%Field%id=Py_Field
        Total%Field%component=GASCOMP
        Total%Field%name='Py'

        CALL CreateTotal(Total)
        Total%Field%id=Pz_Field
        Total%Field%component=GASCOMP
        Total%Field%name='Pz'

     end if
   end SUBROUTINE SetupDiagnostics




  !> Applies initial conditions
  !! @param Info Info object
   SUBROUTINE ProblemGridInit(Info)
      TYPE(InfoDef) :: Info
      if (lrestart) return
      CALL ApplyFunction(SetupAmbient, Info)
      CALL ApplyFunction(SetupJet, Info)        
      IF (nDim == 3 .AND. lMHD) CALL ApplyEMF(JetInitEmf, Info)
!      IF (nDim == 3 .AND. lMHD) CALL ApplyBField(JetInitEmf, Info)
!      write(*,*) 'gridinit applied', Info%q(1,9,1,1)
    END SUBROUTINE ProblemGridInit

    SUBROUTINE JetInitEmf(pos, t, dt, q)!,Ghost)
      REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
      REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
      REAL(KIND=qPREC), INTENT(IN) :: t, dt
      !LOGICAL, INTENT(IN) :: Ghost
      q=0
      q(2)=JetVectorPotentialZ(pos-pos_jet)
    END SUBROUTINE JetInitEmf


    FUNCTION JetVectorPotentialZ(pos)
      REAL(KIND=qPREC) :: JetVectorPotentialZ, pos(3), r, s, ss, a, c, z, const
      s = sqrt(pos(1)**2+pos(3)**2)
      z = pos(2)

      !write(*,*) 'need to calculate vector potential for b field'
      !stop     

      a=1.0/tan(open_angle)
      C=Bjet/Rjet
      
      

!      a = Bjet/Rjet/6
!      c = tan(open_angle) 
!      const=0d0

      r=sqrt(s**2+z**2)
      if (r < Rjet) then
         if (s .le.  a*z) then
!            JetVectorPotentialZ=-a*(r*(s*(3d0-2d0*c*s/z)-2d0*z*s*c)+3d0*z**2*log(r+s))
            JetVectorPotentialZ=-C*(r*(3*a*s*z-2*s**2-2*z**2)+3*a*z**3*log(r+s))/(6*a*z)

            s=a*z
            r=sqrt(s**2+z**2)
            JetVectorPotentialZ=JetVectorPotentialZ + C*(r*(3*a*s*z-2*s**2-2*z**2)+3*a*z**3*log(r+s))/(6*a*z)
         else 
            JetVectorPotentialZ=0
         end if
      else
         JetVectorPotentialZ=0         
      end if
    END FUNCTION  JetVectorPotentialZ

    FUNCTION JetBPhi(pos)
      REAL(KIND=qPREC) :: s, ss, JetBPhi, pos(3), r, snorm
      r=sqrt(sum(pos**2))
      IF (pos(2) > 0d0 .and. r .le. Rjet) THEN         
         s=sqrt(pos(1)**2+pos(3)**2)
         snorm=min(Rjet**2-pos(2)**2,pos(2)/tan(open_angle))
         snorm=pos(2)/tan(open_angle)
         ss=s/snorm
         if (ss > 1) then
            JetBphi=Bamb
         else
            JetBphi=r/Rjet*Bjet*(1d0-ss)+Bamb
         end if
        ! JetBphi=Bjet*pos(2)/r
      ELSE
         JetBPhi=Bamb
      END IF
    END FUNCTION JetBPhi


    SUBROUTINE JetSourceEmf(pos, t, dt, q)
      REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
      REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
      REAL(KIND=qPREC), INTENT(IN) :: t, dt
      REAL(KIND=qPREC), DIMENSION(MaxVars) :: qjet
      REAL(KIND=qPREC) :: BField, pos_rel(3), r, s, ss, Bphi, Btarget(3)
      qjet=0d0
      pos_rel=pos-pos_jet
      r=sqrt(sum(pos_rel**2))
      s=sqrt(pos_rel(1)**2+pos_rel(3)**2)
      if (s < tiny(s)) then
         q=0
      else
         IF (pos_rel(2) > 0d0 .and. r .le. Rjet) THEN         
            Btarget=JetBPhi(pos_rel)*(/-pos_rel(3),0d0,pos_rel(1)/)/s ! phi unit vector        
            CALL SetupJet(pos, t, dt, qjet, .false.)
            q=Cross3D(qjet(imom(1:3)), Btarget)*dt !estimated emf to counter induction equation         
            !         write(*,*) pos
            !         write(*,*) Btarget
            !         write(*,*) qjet(imom(1:3))
            !         write(*,*) t, dt, s, r
         ELSE
            q=0
         END IF
      end if
    END SUBROUTINE JetSourceEmf


    SUBROUTINE SetupAmbient(pos, t, dt, q, Ghost)
      REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
      REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
      REAL(KIND=qPREC), INTENT(IN) :: t, dt
      LOGICAL, INTENT(IN) :: Ghost
      REAL(KIND=xprec) :: r,torus2,cos2theta,ex1,ex2,Ftheta, pos_rel(3), theta
      
      pos_rel=pos-pos_amb
      r=sqrt(sum(pos_rel**2)) ! distance from ambient center
      
      
      !A M B I E N T 
      q(1) = rho_amb
      q(imom(1:nDim)) = 0d0
      if (itracer3 /= 0) q(itracer3)=0d0
      if (itracer5 /= 0) q(itracer5)=0d0

      if (stratified) then !27jan'12
         q(1) = q(1) * ((Ramb_inner/r)**nPower_amb) ! power law profile
         q(imom(1:nDim)) = AGB_speed*pos_rel(1:nDim)/r !radial outflow
      end if

       !!! TORUS
      !From Frnak & Mellema, 1994ApJ...430..800F: produced a pole-to-equator density contrast
      torus2=1d0
      if ( torus ) then
         theta=acos(pos_rel(2)/r)
         cos2theta=cos(2d0*theta) !cos 2* declination angle
         ex1=-2d0*torusbeta
         ex2=torusbeta*cos2theta-torusbeta
         Ftheta=1d0-(torusalpha*(exp(ex2)-1d0)/(exp(ex1)-1d0)  )
         torus2=1d0/Ftheta

         if (r.ge.Ramb_inner) then
            q(1) = q(1) * ((Ramb_inner/r)**nPower_amb) ! power law profile
            q(imom(1:nDim)) = AGB_speed*pos_rel(1:nDim)/r !radial outflow
         end if

         q(1)=q(1)*torus2
                
      end if !if torus
      q(iBz) = Bamb
      q(iE) = q(1)*tamb
    END SUBROUTINE SetupAmbient


    SUBROUTINE SetupJet(pos, time, dt, q, Ghost)
      REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
      REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
      REAL(KIND=qPREC), INTENT(IN) :: time, dt
      LOGICAL, INTENT(IN) :: Ghost
      REAL(KIND=qPREC) :: fact, fact2, pos_rel(3), timef, v(3), theta, r, s, ss, Bphi
      !INITIAL ACCELERATION OF JET/WIND
      IF (time < tt) THEN
         timef = sin(Pi/2d0*time/tt) !rising to 1 following a sin trajectory over tt
      ELSE
         timef=1d0
      END IF

      pos_rel=pos-pos_jet

      r=sqrt(sum(pos_rel**2))
!      q=0
      IF (pos_rel(2)>0d0 .and. r.le.Rjet) THEN !20may2014 Martin version, before 1.2.2015
        ! !======== J E T=========:
        ! IF (outflowType == collimated) then
        !    q(1)=rho_jet
        !    v=vjet*(/0d0,1d0,0d0/)*timef !7 Nov 2014 get rid the Rjet for vjet
        !    q(imom(1:nDim))=v(1:nDim)
        !    q(iE)=q(1)*tjet

        !    !========= W I N D ===========:							
        ! ELSEIF (outflowType == conical) THEN 
        theta=acos(pos_rel(2)/r)
        if (abs(theta) > open_angle) then
           fact=1.00d0*(Pi*half-abs(theta))/ (Pi*half-open_angle) ! ramp down density outside of cone
        else
           fact=1.0
        end if
        q(1)=rho_jet*fact!+(1d0-fact)*q(1)
        v=vjet*pos_rel/Rjet*timef ! ramped up radial velocity (v ~ r)

        !tempering, a la Lee & Sahai, 2003
        if (pos_rel(2) > 0) then
           if (theta < open_angle) then
              fact2=1d0
           else
              fact2=0d0
           end if
           if (tf.ne.0d0) then
              !q(1) = rho_jet*max(0.001d0,exp(-(theta/tf)**2d0))*fact+(1d0-fact)*rho_amb
              q(1) = q(1)*max(0.001d0,exp(-(theta/tf)**2d0))
              fact2=fact2*max(0.001d0,dexp(-(theta/tf)**2d0))
           end if
        else
           fact2=0d0           
        end if
        
        q(imom(1:nDim))=fact2*v

        q(iE)=q(1)*tjet
        ! END IF
        
        
        IF (lMHD) THEN            
           IF (nDim == 2) then ! B_phi in cylindrical geometry is a scalar (not aux field)
              q(iBz)=JetBPhi(pos_rel)
           END IF
        END IF
        if (itracer4 /= 0) q(itracer4)=q(1)

      !!!===== Commented out by B.Liu on 06/08/2018 
       !! necessary to keep Clump on the Jet launch region
       !CALL SetupClump(pos, time, dt, q, ghost)
      ! CALL SetupClumpInJet(pos, time, dt, q, ghost)
            
     !   if (time > clump_starttime) then
     !      if (sqrt(sum((pos-pos_clump+vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2)) < RClump) then
     !         !           q(imom(1:nDim))=q(imom(1:nDim))*q(1)
     !         !           q(1)=q(1)+rho_clump
     !         !           q(imom(2))=q(imom(2))+rho_clump*vClump
     !         !           q(imom(1:nDim))=q(imom(1:nDim))/q(1)
     !         
     !         q(1)=rho_clump
     !         q(iE)=rho_clump*tClump
     !         q(imom(1:nDim))=0d0
     !         q(imom(2))=vClump
     !      end if
     !   end if
     END IF

    END SUBROUTINE SetupJet

    SUBROUTINE SetupClump(pos, time, dt, q, Ghost)
      REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
      REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
      REAL(KIND=qPREC), INTENT(IN) :: time, dt
      LOGICAL, INTENT(IN) :: Ghost

      !if (time >= tplace) then
      if (sqrt(sum((pos-pos_clump-vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2)) <= RClump) then
            !if (sqrt(sum((pos-pos_clump)**2)) <= RClump) then
            !   print*,"vClump=", vClump, "time=", time
            !   print*, "sqrt(sum((pos-pos_clump+vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2))=",sqrt(sum((pos-pos_clump+vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2)), "RClump=", RClump
            !q(imom(1:nDim))=q(imom(1:nDim))*q(1)
            !q(1)=q(1)+rho_clump
            !q(imom(2))=q(imom(2))+rho_clump*vClump
            !q(imom(1:nDim))=q(imom(1:nDim))/q(1)
            q(1)=rho_clump
            q(iE)=rho_clump*tClump
            q(imom(1:nDim))=0d0
            q(imom(2))=vClump
       end if
    END SUBROUTINE SetupClump



    SUBROUTINE SetupClumpInJet(pos, time, dt, q, Ghost)
      REAL(KIND=qPREC), DIMENSION(:), INTENT(IN) :: pos
      REAL(KIND=qPREC), DIMENSION(:), INTENT(INOUT) :: q
      REAL(KIND=qPREC), INTENT(IN) :: time, dt
      LOGICAL, INTENT(IN) :: Ghost
      real(kind=qprec) :: rho_smooth, f, r
!      if (time >= tplace) then
         r=sqrt(sum((pos-pos_clump-vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2))
         if (r <= RClump) then
            !if (sqrt(sum((pos-pos_clump)**2)) <= RClump) then
            !   print*,"vClump=", vClump, "time=", time
            !   print*, "sqrt(sum((pos-pos_clump+vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2))=",sqrt(sum((pos-pos_clump+vClump*(/0d0,1d0,0d0/)*(time-clump_starttime))**2)), "RClump=", RClump
            f=smooth_tanh(r/RClump, thicknessClump)
            rho_smooth=rho_clump*f

            q(imom(1:nDim))=q(imom(1:nDim))*q(1)
            q(1)=q(1)+rho_smooth
            q(imom(2))=q(imom(2))+rho_smooth*vClump
            q(imom(1:nDim))=q(imom(1:nDim))/q(1)
            q(iE)=q(iE)+rho_smooth*tClump
!            q(1)=rho_clump
!            q(iE)=rho_clump*tClump
!            q(imom(1:nDim))=0d0
!            q(imom(2))=vClump
         end if
!      end if
    END SUBROUTINE SetupClumpInJet


!> Applies Boundary conditions
!! @param Info Info object
SUBROUTINE ProblemBeforeStep(Info)
      TYPE(InfoDef) :: Info
      INTEGER :: i,j,k,ii,jj,kk,mx,my,mz,rmbc,zrmbc,level,iErr
      REAL(KIND=xprec) :: x,y,z,x2,y2,z2,xl,yl,zl,zzl,dx,dy,dz,time, r,r2, dt, dr,velFactor,torus2,cos2theta,ex1,ex2,Ftheta, drho
      REAL(KIND=qPREC) :: qnew(NrHydroVars), qwind(NrHydroVars), qjet(NrHydroVars), v(3), rhoFact, fact, fact2, f, timef, bmax, theta
      REAL(KIND=qPREC) :: pressure
      LOGICAL, SAVE :: once=.true.
      REAL(KIND=xprec), SAVE :: dxFIXED=.5d0
      INTEGER, POINTER, DIMENSION(:,:,:) :: mSs
      REAL(KIND=qPREC), POINTER, DIMENSION(:,:) :: offsets
      INTEGER :: nOverlaps

      level=Info%level
      dt=levels(level)%dt
      time=levels(level)%tnow

!      IF (time == tplace .and. level==MaxLevel) then
!         !write (*,*) 'placing clump now'
!         CALL ApplyFunction(SetupClump, Info)
!      end IF
      IF (lPlaceClump) THEN
         myClump%position=pos_Clump
         myClump%velocity=(/0d0,vClump,0d0/)
         CALL UpdateClump(myClump)
         CALL CalcPhysicalOverlaps(Info, myClump%Shape%xBounds, mSs, nOverlaps, offsets, IEVERYWHERE, lHydroPeriodic)
         IF (nOverlaps > 0) THEN
            CALL PlaceClump(Info, myClump, nOverlaps, mSs, offsets)
            
         !   CALL ApplyFunction(SetupClump, Info)
            DEALLOCATE(mSs, offsets)
         END IF
      ELSEIF (lPlaceClumpInBoundaries) THEN
         CALL CalcPhysicalOverlaps(Info, myClump%Shape%xBounds, mSs, nOverlaps, offsets, IBOUNDARY(2,1), lHydroPeriodic)
         IF (nOverlaps > 0) THEN
            CALL PlaceClump(Info, myClump, nOverlaps, mSs, offsets)
            !CALL ApplyFunction(SetupClump, Info)
            DEALLOCATE(mSs, offsets)
         END IF
      END IF
      
     !! Jet_OffTime tested and works
      !IF (time + dt .le. Jet_OffTime) THEN
      IF (time  .le. Jet_OffTime) THEN
 !! commented CALL SetupJet line  out by Baowei on June 8th to avoid the stepping of Clump zone.
 !! 
!         CALL ApplyFunction(SetupJet, Info)        
         IF (nDim == 3 .AND. lMHD) CALL ApplyEMF(JetSourceEmf, Info)
         !write(*,*) 'beforestep applied', Info%q(1,9,1,1)
         !IF (nDim == 3 .AND. lMHD) CALL ApplyEMF(JetInitEmf, Info)
         !if (time > clump_starttime) then
         !   CALL ApplyFunction(SetupClump, Info)
         !   print*, "Clump setup"
         !   stop
         !end if
      END IF

            
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
  INTEGER :: massPrintTime=0 !! only need to print the initial mass once.

  !  IF (n == 0 .AND. abs(levels(n)%tnow - Clump_StartTime) &
  !   < tiny(Clump_StartTime) .AND. .NOT. lClumpAdded ) THEN
  !    CALL AddClump()
  !    lClumpAdded=.true.
  ! END IF

  IF (n == 0) THEN
     IF (levels(n)%tnow .le. Clump_StartTime .AND. levels(n)%tnow + levels(n)%dt > Clump_StartTime) THEN
        !write(*,*) 'turning on place clump'
        CALL CreateRefinement(ClumpRefinement)
        CALL CreateShape(ClumpRefinement%Shape, RECTANGULAR_PRISM, (/RClump,RClump,RClump/))
        ClumpRefinement%Shape%position=pos_clump
        CALL UpdateShape(ClumpRefinement%Shape)
        ClumpRefinement%type=REFINE_INSIDE
        CALL UpdateRefinement(ClumpRefinement)
        tplace=levels(n)%tnow+levels(n)%dt
        lPlaceClump=.true.
     ELSE IF (abs(levels(n)%tnow-tplace) < 1e-10) THEN
        !write(*,*) 'PlaceClump is true'
        lPlaceClump=.true.
        lRefineClump=.false.
     ELSE IF (lPlaceClump) THEN
        !write(*,*) 'destroying refinement'
        lPlaceClumpInBoundaries=.true.
        lPlaceClump=.false.
        IF (ASSOCIATED(ClumpRefinement)) CALL DestroyRefinement(ClumpRefinement)
        
     END IF
  END IF
END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem

