!#########################################################################
!		
!    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
    
    PRIVATE 
    type(ambientdef), pointer :: ambient
    TYPE(ParticleDef), POINTER :: Particle
    TYPE(WindDef), POINTER :: Wind
    TYPE(ParticleListDef), POINTER :: particlelist
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
    open(unit=problem_data_handle,file='problem.data',status='old')
    read(problem_data_handle,nml=AmbientData)
!    DO i=1,nDim
!        DO edge=1,2                   
!            IF (MPI_ID == 0 .AND. Gmthbc(i,edge) == 1) THEN
!            If (ASSOCIATED(Wind)) Then
!                Print *, "Yes Wind"
!            Else
!                Print *,"No Wind"
!            End If 
!                NULLIFY(Wind)
!                CALL CreateWind(Wind)
!                Wind%density=rhoambient/rscale
!                Wind%temperature=pambient/(rhoambient/amu/xmu)/boltzmann/tempscale
!                Wind%dir=i
!                Wind%edge=edge
!                Wind%type=OUTFLOW_ONLY
!            END IF
!        END DO
!    END DO
    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)
        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=5e-10/rscale
        Particle%OutflowObj%temperature=30000/tempscale
        Particle%OutflowObj%speed=1e8/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),a(5,128),b(5,192)
    INTEGER :: i,j,k,nx,ny,nz,rmbc,zrmbc,level,ii,jj,kk
    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')
    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')
    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')
    END SELECT
    do i=1,19
       read(195,*)
       read(196,*)
       read(197,*)
       read(198,*)
       read(199,*)
    end do
    select case(info%level)
    case(0)
        do k=1,128
           do j=1,128
              read(195,*),a(1,:)
              read(196,*),a(2,:)
              read(197,*),a(3,:)
              read(198,*),a(4,:)
              read(199,*),a(5,:)
              if (j.ge.57.and.j.le.72) then
                  do i=1,16
                      Info%q(i,j-56,k,1)=a(1,56+i)/rscale
                      Info%q(i,j-56,k,2)=a(1,56+i)*a(2,56+i)/rscale/velscale
                      Info%q(i,j-56,k,3)=a(1,56+i)*a(3,56+i)/rscale/velscale
                      Info%q(i,j-56,k,4)=a(1,56+i)*a(4,56+i)/rscale/velscale
                      Info%q(i,j-56,k,5)=a(1,56+i)*(a(5,56+i)+0.5*(a(2,56+i)&
                      *a(2,56+i)+a(3,56+i)*a(3,56+i)+a(4,56+i)*a(4,56+i)))/pscale
                  end do
              end if
           end do
        end do
    case(4)
        do k=1,128
           do j=1,128
              read(195,*),a(1,:)
              read(196,*),a(2,:)
              read(197,*),a(3,:)
              read(198,*),a(4,:)
              read(199,*),a(5,:)
              Info%q(1:128,j,k,1)=a(1,:)/rscale
              Info%q(1:128,j,k,2)=a(1,1:128)*a(2,1:128)/rscale/velscale
              Info%q(1:128,j,k,3)=a(1,1:128)*a(3,1:128)/rscale/velscale
              Info%q(1:128,j,k,4)=a(1,1:128)*a(4,1:128)/rscale/velscale
              Info%q(1:128,j,k,5)=a(1,1:128)*(a(5,1:128)+0.5*(a(2,1:128)*a(2,1:128)+a(3,1:128)*a(3,1:128)+a(4,1:128)*a(4,1:128)))/pscale
           end do
        end do
    case(7)
        do k=1,192
           do j=1,192
              read(195,*),b(1,:)
              read(196,*),b(2,:)
              read(197,*),b(3,:)
              read(198,*),b(4,:)
              read(199,*),b(5,:)
              Info%q(1:192,j,k,1)=b(1,:)/rscale
              Info%q(1:192,j,k,2)=b(1,1:192)*b(2,1:192)/rscale/velscale
              Info%q(1:192,j,k,3)=b(1,1:192)*b(3,1:192)/rscale/velscale
              Info%q(1:192,j,k,4)=b(1,1:192)*b(4,1:192)/rscale/velscale
              Info%q(1:192,j,k,5)=b(1,1:192)*(b(5,1:192)+0.5*(b(2,1:192)*b(2,1:192)+b(3,1:192)*b(3,1:192)+b(4,1:192)*b(4,1:192)))/pscale
           end do
        end do
    end select
    close(195)
    close(196)
    close(197)
    close(198)
    close(199)
END SUBROUTINE ProblemGridInit

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

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
    REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,pos(3)
    if (lrestart) then
        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)
         
        DO k=1, mz
            z = (zl+(REAL(k,xPrec)-half)*dz)
            DO i=1, mx
                x = (xl+(REAL(i,xPrec)-half)*dx)
                DO j=1, my
                    y = (yl+(REAL(j,xPrec)-half)*dy)
                    if (Info%level<4) then
                        if (x>-4000 .and. x<4000 .and. y>-4000 .and. y<4000 .and. z>-4000 .and. z<4000) then !used for 6000d quiescent phase
                            !.or. (x>-1000 .and. x<1000 .and. y>-1000 .and. y<1000 .and. &
                            !z>-5000 .and. z<5000)) then
                            Info%ErrFlag(i,j,k)=0   !1
                        else
                            Info%ErrFlag(i,j,k)=0
                        end if
                    else if (Info%level<8) then
                        if (x>-750 .and. x<750 .and. y>-750 .and. y<750 .and. z>-750 .and. z<750) then ! used for 6000d quiescent phase  
                            !.or. (x>-600 .and. x<600 .and. y>-600 .and. y<600 .and. &
                            !z>-4000 .and. z<4000)) then
                            Info%ErrFlag(i,j,k)=0  !1
                        else
                            Info%ErrFlag(i,j,k)=0
                        end if




                        if (x>4000 .or. x<-4000 .or. y>4000 .or. y<-4000 .or. &
                            z>4000 .or. z<-4000) then
                            Info%ErrFlag(i,j,k)=0
                        end if
                    end if
                END DO
            END DO
        END DO
    end if
    IF (lRestart) RETURN
    SELECT CASE(Info%level)
    CASE(0)
       Info%ErrFlag=0
        Info%ErrFlag(5:12,5:12,61:68)=1
        !Info%ErrFlag(61:68,61:68,61:68)=1
        !Info%ErrFlag(64:65,64:65,58:71)=1
    CASE(1:3)
        Info%ErrFlag=1
    CASE(4)
       Info%ErrFlag=0
        Info%ErrFlag(53:76,53:76,53:76)=1
        !Info%ErrFlag(64:65,64:65,50:79)=1
    case(5:6)
        Info%ErrFlag=1
    case(7)
    END SELECT
END SUBROUTINE ProblemSetErrFlag

SUBROUTINE ProblemBeforeGlobalStep(n)
    INTEGER :: n,i
    real :: starttime,endtime
    starttime=6000.0
    endtime=20000.0
    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<starttime .or. time>endtime) then
        particle%outflowobj%speed=0/velscale
        Particle%OutflowObj%density=5e-13/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=1e-11/rscale              
        ! 1e-11 for high momentum, 5e-13 for low momentum 
        Particle%OutflowObj%temperature=30000/tempscale
        Particle%OutflowObj%open_angle=Pi/2
        Particle%OutflowObj%fade=0d0
        Particle%OutflowObj%duration=3e9/timescale
    end if
END SUBROUTINE ProblemBeforeGlobalStep

! Routines that set initial and boundary conditions for stellar wind

END MODULE Problem
