u/erica/CoreCollapseBlog: problem.f90_mar15

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

!> @file problem.f90
!! @brief Main file for module Problem
MODULE Problem
  USE DataDeclarations
  USE Ambients
  !USE PointGravitySrc
  !USE ParticleDeclarations !n
  USE Clumps
  USE Profiles
  USE ProfileFunctions
  USE Winds
  IMPLICIT NONE
  SAVE

  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep

CONTAINS

  !> Initializes module variables
  SUBROUTINE ProblemModuleInit()

    TYPE(AmbientDef), POINTER :: Ambient       
    TYPE(WindDef), POINTER :: Wind
    TYPE(ClumpDef), POINTER :: Clump
    INTEGER :: edge
    INTEGER :: i, nEntries !soft_function
    REAL(KIND=qPREC), DIMENSION(3) :: xloc=(/0,0,0/), velocity= (/0,0,0/)  
    REAL(KIND=qPREC) :: rhoOut, pOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut, r_cutoff


    NAMELIST /ClumpData/ xloc, velocity!, soft_radius, soft_function!, mass
    NAMELIST /AmbientData/ rhoOut, pOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut

    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
    READ(PROBLEM_DATA_HANDLE,NML=AmbientData)
    CALL CreateAmbient(Ambient)
    Ambient%density=rhoOut/rscale
    !Ambient%density=.1
    !Ambient%pressure=10
    Ambient%pressure=pOut/pscale
    Ambient%PersistInBoundaries=.true.
    CALL UpdateAmbient(Ambient)

    READ(PROBLEM_DATA_HANDLE,NML=ClumpData)
    CALL CreateClump(Clump)
    Clump%position=xloc
    Clump%velocity=velocity
    !Clump%density=1 these things are set using the profile object
    !Clump%temperature=1 ""
    Clump%thickness=0d0
    Clump%subsample=1

    CLOSE(PROBLEM_DATA_HANDLE)

    OPEN (UNIT=PROBLEM_DATA_HANDLE, FILE='progenitor_profiles.dat', STATUS="OLD", POSITION='REWIND')
!    PRINT *, 'OPENED FILE OKAY'
    READ (PROBLEM_DATA_HANDLE, *) nEntries
 !   PRINT *, 'nEntries=', nEntries
    READ (PROBLEM_DATA_HANDLE, *)  !Skip the header column      
    CALL CreateProfile(Clump%Profile, nEntries, (/Mass_Field, P_Field/), RADIAL) !Mass field - is gas density
    DO i=1,nEntries
       READ(PROBLEM_DATA_HANDLE, *) Clump%profile%data(i,:) ! position, density, pressure in cgs
  !         PRINT *, 'CLUMP%PROFILE%DATA(I,1)=radius (cm)', 'CLUMP%PROFILE%DATA(I,2)=density (g/cc)', 'CLUMP%PROFILE%DATA(I,3)=press(dyn/cc)'
  !         PRINT *, CLUMP%PROFILE%DATA(I,1), CLUMP%PROFILE%DATA(I,2), CLUMP%PROFILE%DATA(I,3)
       Clump%profile%data(i,:)=Clump%profile%data(i,:) / (/lScale, rScale, pScale/) !rescale to cu
    !       PRINT *, 'CLUMP%PROFILE%DATA(I,1)=radius (cu)', 'CLUMP%PROFILE%DATA(I,2)=density (cu)', 'CLUMP%PROFILE%DATA(I,3)=press(cu)'
    !       PRINT *, CLUMP%PROFILE%DATA(I,1), CLUMP%PROFILE%DATA(I,2), CLUMP%PROFILE%DATA(I,3)
    END DO
    CLOSE(PROBLEM_DATA_HANDLE)


    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='profile_init.out')
    WRITE(PROBLEM_DATA_HANDLE,'(3E25.15)') (/(clump%profile%data(i,:)*(/lscale, rScale, pScale/), i=1,nEntries)/) 
    CLOSE(PROBLEM_DATA_HANDLE)

    pOut=pOut/pscale

    CALL Profile_SelfGravityHSE_Outside(Clump%Profile, pOut) !Use scaled ambient pressure from profile data file to recalculate pressure profile in HSE
    CALL UpdateProfile(Clump%Profile) !don't think this is actually needed -- subroutine is a stub in the code and data values are updated in prev func
  
    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='profile.out')
    WRITE(PROBLEM_DATA_HANDLE,'(3E25.15)') (/(clump%profile%data(i,:)*(/lscale, rScale, pScale/), i=1,nEntries)/) 
    CLOSE(PROBLEM_DATA_HANDLE)

    r_cutoff=Clump%profile%data(nEntries,1)
    Clump%radius=r_cutoff
    print *, 'clump%radius=', clump%radius, 'clump%radius (cm)=', clump%radius*lscale
    CALL UpdateClump(Clump)

    DO i=1,nDim
       DO edge=1,2
          IF (Gmthbc(i,edge) == 1) THEN 
             CALL CreateWind(Wind)
             Wind%density=rhoOut
             wind%temperature=tempOut
             wind%type=OUTFLOW_ONLY
             Wind%dir=i
             Wind%edge=edge
          END IF
       END DO
    END DO

  END SUBROUTINE ProblemModuleInit

  !> Applies initial conditions
  !! @param Info Info object
  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
  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
      INTEGER :: i,j,k
      Info%ErrFlag(:,:,:)=1  ! refines over entire grid
  END SUBROUTINE ProblemSetErrFlag

  SUBROUTINE ProblemBeforeGlobalStep(n)
     INTEGER :: n
  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem