!#########################################################################
!
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    problem.f90 of module Template 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/>.
!
!#########################################################################
!> @dir Template
!! @brief Contains files necessary for the Template Calculation

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

!> @defgroup Template Template Module
!! @brief Module for calculating collapse of a uniform cloud
!! @ingroup Modules

!> Template Module
!! @ingroup Template
MODULE Problem
  USE GlobalDeclarations
  USE DataDeclarations
  USE Clumps
  USE Ambients
  USE Fields
  USE Refinements
  USE Shapes
  USE Clumps
  USE ProcessingDeclarations
  USE Winds
  USE Tests
  USE TemplateControl
  USE Outflows
  USE InternalBoundary
  USE PhysicsDeclarations
  USE SweepDeclarations
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
       ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep

  PRIVATE

  ! distRatio: core_distance/wire radius
  REAL(KIND=qPREC) :: pos(3), rWire, coreDistRatio = 5d0, rhoBoundary, tBoundary, rhoWire, tWire, rhoAmb, tAmb, rhoWind, tWind, betaWind, Bdir(3), bwind, vwind
  REAL(KIND=qPREC) :: bodyRes=1d6, ambRes=1d-3, shellFrac=0.1, shellResFrac=0.1, coreFrac=0d0, coreResFrac=1d-8, criticalRes=-1.0, rExcl=0d0
  INTEGER :: refFlag = 1, wAxis = 1, rwFlag = 0, RESRAMP = 3, excludeBField = 0
  REAL(KIND=qPREC) :: wDir(3) = (/1d0,0d0,0d0/)
  REAL(KIND=qPREC) :: tRamp = 1.224d0, vAmb, delta_vel=1d10, delta_T, delta_rho, bAmb(3), shellRes, OKres, diff_alpha2_0, rampstart, ramptime
  TYPE(InternalBoundaryDef), POINTER :: rInternalBoundary, lInternalBoundary
  TYPE(WindDef), POINTER :: Wind

CONTAINS

  SUBROUTINE ProblemModuleInit() BIND(C)
    TYPE(AmbientDef), POINTER :: Ambient
    TYPE(RefinementDef), POINTER :: Refinement			! For testing AMR (1 Aug'19)
    INTEGER :: itrace
    REAL(KIND=qPREC) ::cs, vxnorm=-1, vxnormR=1
    real, pointer :: a
    NAMELIST/ProblemData/ pos, coreDistRatio, rWire, rhoBoundary, tBoundary, rhoAmb, tAmb, rhoWind, tWind, bWind, Bdir, vWind, vxnorm, vxnormR, tRamp, &
	  & bodyRes, ambRes, shellFrac, shellResFrac, coreFrac, criticalRes, refFlag, waxis, bAmb, rwFlag, rhoWire, tWire, vAmb

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


    rhoWind=rhoWind/nScale
    tWind=tWind/TempScale
    !bWind=bWind*1e-9*tesla/BScale
    bWind=bWind/BScale
    vWind=vWind*1e5/velScale
    !bWind=sqrt(2d0*rhoWind*tWind/betaWind)
    !vWind=machWind*sqrt(gamma*tWind)
    
    vAmb=vAmb*1e5/velScale
    bAmb=bAmb/BScale

    rhoAmb=rhoAmb/nScale
    tAmb=tAmb/TempScale
    rhoBoundary=rhoBoundary/nScale
    tBoundary=tBoundary/TempScale

    rhoWire=rhoWire/nScale  !rhoWire=8.73 g/cm^3
    tWire=tWire/TempScale
    rWire=rWire/lScale

    wdir=0d0
    wdir(waxis)=1


    CALL CreateAmbient(Ambient, rhoAmb, rhoAmb*tAmb)
    Ambient%velocity=vAmb*wDir
    Ambient%B = bAmb
    ! Make 0 Ambient for case with velocity Ramping up or Excluded Clump
    CALL UpdateAmbient(Ambient)

    CALL CreateWind(Wind, rhoWind, tWind, vWind)
    Wind%dir=wAxis

! Change to 1, 2, 3 (x, y, z)
    Wind%edge=1				! Left or Right edge
    Wind%B=bWind*Bdir
    CALL UpdateWind(Wind)

    CALL CreateInternalBoundary(lInternalBoundary)
    lInternalBoundary%Shape%type=CYLINDER
    lInternalBoundary%Shape%position=pos
    lInternalBoundary%Shape%size_param=rWire*(/1d0,1d0,1d30/)
    CALL SetShapeOrientation(lInternalBoundary%Shape, 0d0, 0d0, 0d0) !creates rotation matrix
    lInternalBoundary%density=rhoWire
    lInternalBoundary%pressure=rhoWire*tWire
    lInternalBoundary%type=REFLECTING
    lInternalBoundary%vnormfactor=-1
    lInternalBoundary%MaxLevel=MaxLevel
    CALL UpdateInternalBoundary(lInternalBoundary)


    CALL CreateInternalBoundary(lInternalBoundary)
    lInternalBoundary%Shape%type=CYLINDER
    lInternalBoundary%Shape%position=pos+(/0d0,coreDistRatio,0d0/)*rWire
    lInternalBoundary%Shape%size_param=rWire*(/1d0,1d0,1d30/)
    CALL SetShapeOrientation(lInternalBoundary%Shape, 0d0, 0d0, 0d0) !creates rotation matrix
    lInternalBoundary%density=rhoWire
    lInternalBoundary%pressure=rhoWire*tWire
    lInternalBoundary%type=REFLECTING
    lInternalBoundary%vnormfactor=-1
    lInternalBoundary%MaxLevel=MaxLevel
    CALL UpdateInternalBoundary(lInternalBoundary)
  END SUBROUTINE ProblemModuleInit

  SUBROUTINE ProblemGridInit(Info) BIND(C)
    TYPE(InfoDef) :: Info

  END SUBROUTINE ProblemGridInit


  SUBROUTINE ProblemBeforeStep(Info) BIND(C)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemBeforeStep

  SUBROUTINE ProblemAfterStep(Info) BIND(C)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  SUBROUTINE ProblemSetErrFlag(Info) BIND(C)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemSetErrFlag


  SUBROUTINE ProblemBeforeGlobalStep(n)		 BIND(C)! Modify time variable parameter here
    INTEGER :: n
  END SUBROUTINE ProblemBeforeGlobalStep


END MODULE Problem

