!#########################################################################
!		
!    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 PDFs
    USE Projections
    USE Refinements
    USE ProcessingDeclarations
    USE Winds
    USE TemplateControl
    IMPLICIT NONE
    SAVE
    PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
       ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
    REAL :: time, efact
    
    PRIVATE 
    type(ambientdef), pointer :: ambient
    TYPE(ParticleDef), POINTER :: Particle
    TYPE(WindDef), POINTER :: Wind
    TYPE(ParticleListDef), POINTER :: particlelist
    type(refinementdef), pointer :: refinejet
    integer, dimension(3,2) :: data_offsets_0, data_offset_4, data_offset_7
    integer, dimension(3,2) :: ibox_0, ibox_4, ibox_7, ref_0, ref_4
    real, dimension(3,2):: box_0, box_4, box_7
    real :: outflowstarttime, outflowendtime, outflowdensity 
    real(kind=qprec) :: outflow_extent=0
    real :: mintracer=.1
CONTAINS

SUBROUTINE ProblemModuleInit()
    INTEGER :: n,i,j,k,edge
    real(KIND=qPREC) :: rhoambient,pambient,mass,xloc(3),vel(3),alpha,radiusw
    integer :: grav_soft_rad
    TYPE(ProjectionDef), POINTER :: Projection
    NAMELIST/AmbientData/ rhoambient,pambient
    NAMELIST /ParticleData/ mass,xloc,vel,alpha,radiusw,grav_soft_rad
    NAMELIST/ProblemData/ outflowstarttime, outflowendtime, outflowdensity, mintracer 
    open(unit=problem_data_handle,file='problem.data',status='old')
    read(problem_data_handle,nml=ProblemData)

    ! Make sure interpopts starts at zero for initial grid setup
    if (.not. lRestart) then 
        InterpOpts=0
        MinDensity=4d-6  ! << set mindensity higher to overwrite the odd spherical region, change to 1d-10 later in ProblemBeforeGlobalStep 
        MinTemp=20       ! << set mintemp higher for frame 0, change to 10 later in ProblemBeforeGlobalStep 
    end if 

    if (lRestart) then 
        particle=>SinkParticles%self                ! handle to the particle 
        Call addtracer(particle%outflowobj%iTracer,'Outflow_Tracer')  ! <<< tracers need to be added in the same order when the code is restarted, but don't re-creat the particle - or the outflow object   
    end if 

    ! Physical bounds of boxes
    box_0=64000*reshape((/-1,-1,-1,1,1,1/),(/3,2/))
    box_4=4000*reshape((/-1,-1,-1,1,1,1/),(/3,2/))
    box_7=750*reshape((/-1,-1,-1,1,1,1/),(/3,2/))

    ! index of boxes in each level's index space
    ibox_0=nint(box_0-spread(GxBounds(:,1),2,2))/levels(0)%dx + spread((/1,0/),1,3)
    if (maxlevel >= 4) then
       ibox_4=nint(box_4-spread(GxBounds(:,1),2,2))/levels(4)%dx + spread((/1,0/),1,3)
       if (maxlevel >= 7) then
          ibox_7=nint(box_7-spread(GxBounds(:,1),2,2))/levels(7)%dx + spread((/1,0/),1,3)
       end if
    end if
    
    CALL CreateRefinement(RefineJet)
    CALL CreateShape(RefineJet%shape, RECTANGULAR_PRISM, (/750d0,750d0,0d0/)) !position defaults to 0,0,0
    RefineJet%type=REFINE_INSIDE

    !> internal energy of ideal gas goes as T/(gamma-1)
    !> So to convert internal energy from gamma=5/3 to the current simulation gamma,  
    !> we need to multiply by (5/3-1) and divide by (gamma-1)
    efact = (5.0/3.0-1.0)/(gamma-1.0) !< factor for scaling internal energy from input data

    !    call createambient(ambient)
    !    ambient%density=rhoambient
    !    ambient%pressure=pambient
    READ(problem_data_handle,nml=ParticleData)
    if (.not. lrestart) then
       NULLIFY(Particle)
       CALL CreateParticle(Particle)
       Particle%Q(1)=mass/mscale
       Particle%xloc=xloc/lscale
       Particle%Q(imom(1:nDim))=vel(1:nDim)/velscale
       CALL AddSinkParticle(Particle)
       Particle%iAccrete=NOACCRETION
       Particle%lfixed=.true.
       CALL CreatePointGravityObject(Particle%PointGravityObj)
       Particle%PointGravityObj%alpha=1
       Particle%PointGravityObj%soft_length=REAL(grav_soft_rad,qPREC)*sink_dx
       Particle%PointGravityObj%soft_function=SPLINESOFT
       Particle%PointGravityObj%Mass=Particle%Q(1)
       Particle%PointGravityObj%v0(1:nDim)=Particle%Q(imom(1:nDim))
       Particle%PointGravityObj%x0=Particle%xloc
       !        Particle%PointGravityObj%luminosity=3300*lsun/(pscale/timescale)
       CALL CreateOutflow(Particle%OutflowObj)
       call addtracer(particle%outflowobj%iTracer,'Outflow_Tracer')
       Particle%OutflowObj%duration=3e9/timescale
       Particle%OutflowObj%radius=REAL(radiusw,qPREC)*sink_dx
       Particle%OutflowObj%thickness=REAL(radiusw,qPREC)*sink_dx
       Particle%OutflowObj%open_angle=Pi/2
       Particle%OutflowObj%theta=0
       Particle%OutflowObj%phi=0
       Particle%OutflowObj%density=outflowdensity/rscale
       Particle%OutflowObj%temperature=30000/tempscale
       Particle%OutflowObj%speed=0/velscale
       Particle%OutflowObj%fade=0d0
       Particle%OutflowObj%velocity(1:nDim)=Particle%Q(imom(1:nDim))
       Particle%OutflowObj%position= Particle%xloc
       Particle%OutflowObj%begin=-1d0
       CALL UpdateOutflow(Particle%OutflowObj)
    end if
    close(problem_data_handle)
    !CALL PostProcessInit()
END SUBROUTINE ProblemModuleInit

SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,pos(3)
    INTEGER :: i,j,k,nx,ny,nz,rmbc,zrmbc,level,ii,jj,kk, ibox(3,2), mb(3,2), mc(3,2)
    real(kind=qprec), dimension(:,:,:,:), allocatable :: a

    nx=Info%mx(1)
    ny=Info%mx(2)
    nz=Info%mx(3)
    IF (lRestart) RETURN
    IF (.NOT. ANY(Info%level == (/0,4,7/))) RETURN
    SELECT CASE(Info%level)
    CASE(0)
       open(unit=195,file='NewLargeGrid_density_g_cm3_grid.dat',status='old',form='formatted',action='read')
       open(unit=196,file='NewLargeGrid_v_x_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=197,file='NewLargeGrid_v_y_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=198,file='NewLargeGrid_v_z_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=199,file='NewLargeGrid_u_erg_g_grid.dat',status='old',form='formatted',action='read')
       ibox=ibox_0
    CASE(4)
       open(unit=195,file='NewMedGrid_density_g_cm3_grid.dat',status='old',form='formatted',action='read')
       open(unit=196,file='NewMedGrid_v_x_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=197,file='NewMedGrid_v_y_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=198,file='NewMedGrid_v_z_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=199,file='NewMedGrid_u_erg_g_grid.dat',status='old',form='formatted',action='read')
       ibox=ibox_4
    CASE(7)
       open(unit=195,file='NewSmallGrid_density_g_cm3_grid.dat',status='old',form='formatted',action='read')
       open(unit=196,file='NewSmallGrid_v_x_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=197,file='NewSmallGrid_v_y_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=198,file='NewSmallGrid_v_z_cm_s_grid.dat',status='old',form='formatted',action='read')
       open(unit=199,file='NewSmallGrid_u_erg_g_grid.dat',status='old',form='formatted',action='read')
       ibox=ibox_7
    END SELECT    
 
   do i=1,19  !< Reads the first 19 lines from each file - presumably containing meta data
       read(195,*)
       read(196,*)
       read(197,*)
       read(198,*)
       read(199,*)
    end do

    ! Allocate space to store local box
    allocate(a(ibox(1,1):ibox(1,2), ibox(2,1):ibox(2,2), ibox(3,1):ibox(3,2), 5))

    ! Read in data
    read(195,*),a(:,:,:,1) !< must be density in g
    read(196,*),a(:,:,:,2) !< x-velocity in cm/s
    read(197,*),a(:,:,:,3) !< y-velocity in cm/s
    read(198,*),a(:,:,:,4) !< z-velocity in cm/s
    read(199,*),a(:,:,:,5) !< internal energy in erg/g
    close(195)
    close(196)
    close(197)
    close(198)
    close(199)

    a(:,:,:,5)=a(:,:,:,5)*a(:,:,:,1)*efact / pscale !< rescale internal energy for simulation gamma
    a(:,:,:,1)=a(:,:,:,1)/rscale
    a(:,:,:,2:4)=a(:,:,:,2:4)/velscale * spread(a(:,:,:,1),4,3)
    a(:,:,:,5)=a(:,:,:,5)+.5*sum(a(:,:,:,2:4)**2,4)/a(:,:,:,1)

    ! Calculate overlap
    mb(:,2)=min(ibox(:,2), info%mglobal(:,2))
    mb(:,1)=max(ibox(:,1), info%mglobal(:,1))
    mc=mb-spread(info%mglobal(:,1),2,2)+1
    
    ! Update Info%q
    Info%q(mc(1,1):mc(1,2), mc(2,1):mc(2,2), mc(3,1):mc(3,2),1:5)=a(mb(1,1):mb(1,2), mb(2,1):mb(2,2), mb(3,1):mb(3,2), :)
    deallocate(a)
END SUBROUTINE ProblemGridInit

SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
    integer :: k
    real(KIND=qPREC) :: y
    if (outflow_extent < info%xbounds(3,2) .or. -outflow_extent > info%xbounds(3,1)) then
       do k=1,Info%mx(3)
          y=Info%xBounds(3,1)+(real(k)-.5)*levels(info%level)%dx
          if (outflow_extent<abs(y)) then
             if (maxval(Info%q(1:Info%mX(1),1:Info%mX(2),k,particle%outflowobj%itracer)) > mintracer) then
                outflow_extent=abs(y)
             end if
          end if
       end do
    end if
END SUBROUTINE ProblemBeforeStep

!> on Jan 31, 2019 
!> Added the following subroutine to ensure that 
!> the outflow is not cooled to minTemp.      
SUBROUTINE ProtectOutflowTemp(pos,t,dt,q,Ghost)
  USE GlobalDeclarations
  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
  if (q(particle%outflowobj%itracer) > 1e-4*particle%outflowobj%density) then
     q(iE)=max(q(1)*FloorTemp/TempScale, q(iE))
  end if
END SUBROUTINE ProtectOutflowTemp

SUBROUTINE PostProcessInit()
    TYPE(ProjectionDef), POINTER :: Projection
    CALL CreateProjection(Projection)
    Projection%field%id=MASS_FIELD
    Projection%dim=1

    CALL CreateProjection(Projection)
    Projection%field%id=MASS_FIELD
    Projection%dim=2

    CALL CreateProjection(Projection)
    Projection%field%id=MASS_FIELD
    Projection%dim=3
END SUBROUTINE PostProcessInit

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

SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
    REAL(KIND=xprec) :: vrot(3),xloc_primary(3),xloc_secondary(3),iparticle
    INTEGER :: n,i,j,k,mx,my,mz,rmbc,zrmbc,level,ii,jj,kk, mb(3,2), mc(3,2)
    REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,pos(3)
    mb(:,1)=1
    mb(:,2)=0
    select case(info%level)
    case(0:3)
       if (maxlevel >= 4) mb=leveldown(ibox_4,4,info%level)
    case(4:6)
       if (maxlevel >= 7) mb=leveldown(ibox_7,7,info%level)
    end select
    mb(:,1)=max(mb(:,1), info%mglobal(:,1))
    mb(:,2)=min(mb(:,2), info%mglobal(:,2))
    mc=mb-spread(info%mglobal(:,1),2,2)+1
    info%errflag(mc(1,1):mc(1,2), mc(2,1):mc(2,2), mc(3,1):mc(3,2))=1
END SUBROUTINE ProblemSetErrFlag

SUBROUTINE ProblemBeforeGlobalStep(n)      
    INTEGER, INTENT(IN) :: n
    INTEGER :: i
    INTEGER :: ierr

    ! recalculate outflow_extent before updating the refinement shape's extends 
    ! prior to regridding - and should preserve the refined data 
    TYPE(nodedeflist), pointer :: nodelist 
    TYPE(nodedef), pointer :: node 
    if (lRestart .and. outflow_extent==0) then
        DO i=0, maxlevel 
            nodelist=>Nodes(i)%p 
            DO WHILE (associated(nodelist))
                node=>nodelist%self
                CALL ProblemBeforeStep(node%info)
                nodelist=>nodelist%next
            end DO
        end DO 
    end if 

    if (levels(n)%dt>0) then  ! make sure there is a time step associated with the step   
    ! << set mindensity and mintemp lower for the rest of simulation 
        MinDensity=1d-10
        MinTemp=10 
    !Set interpopts back to 1 for rest of simulation
        InterpOpts=1
    end if 

    time=levels(n)%tnow
    i=1
    particlelist=>sinkparticles
    DO WHILE (ASSOCIATED(particlelist))
        particle=>particlelist%self
        particlelist=>particlelist%next
        i=i+1
    END DO
    if (time<outflowstarttime .or. time>outflowendtime) then
        particle%outflowobj%speed=0/velscale
        Particle%OutflowObj%density=outflowdensity/rscale
        Particle%OutflowObj%temperature=30000/tempscale
        Particle%OutflowObj%open_angle=Pi/2
        Particle%OutflowObj%fade=0d0
        Particle%OutflowObj%duration=3e9/timescale
    else
        particle%outflowobj%speed=3e7/velscale
        Particle%OutflowObj%density=outflowdensity/rscale  
!1e-11 for hign momentum flux, and 5e-13 for low momentum flux  
        Particle%OutflowObj%temperature=30000/tempscale
        Particle%OutflowObj%open_angle=Pi/2
        Particle%OutflowObj%fade=0d0
        Particle%OutflowObj%duration=3e9/timescale
    end if
#if defined _MPI
    call mpi_allreduce(MPI_IN_PLACE, outflow_extent, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)
#endif

    RefineJet%shape%size_param(3)=outflow_extent
    CALL UpdateShape(RefineJet%shape)
    
END SUBROUTINE ProblemBeforeGlobalStep

! Routines that set initial and boundary conditions for stellar wind

END MODULE Problem