Blog: Meeting Update June-17th Erini: problem.f90

!#########################################################################
!               
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    problem.f90 of module BonnorEbertSphere 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/>.
!
!#########################################################################
!============================================================================
! This problem places a bonner ebert sphere in the center of the grid. 
!============================================================================


MODULE PROBLEM
  USE DataDeclarations
  USE GlobalDeclarations
  USE PhysicsDeclarations
  USE LE_MODULE
  USE Profiles
  USE WINDS
  USE CLUMPS
  USE Ambients
  USE Totals

  IMPLICIT NONE    ! It's safer to require explicit declarations
  SAVE             ! Save module information
  PRIVATE

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

  TYPE(WindDef), POINTER :: Wind
  TYPE(ClumpDef), POINTER :: MyClump
  TYPE(AmbientDef), POINTER :: Ambient
  TYPE(TotalDef), POINTER :: Total
  TYPE(InfoDef), POINTER :: Info
  TYPE(ProfileDef), Pointer :: Profile 
 INTEGER :: nWinds
  REAL(KIND=qPrec), DIMENSION(3) :: xloc=(/0,0,0/) 
REAL(KIND=qPrec) :: R_outer,h,dtheta,a_scale,p_cent,temp_cent
REAL :: xi_last,xi, at,poly_index,M_star,rho_cent
REAL(KIND=qPrec) :: rhoOut, pOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut, crit_rad, p_crit
REAL(KIND=qPrec) :: Omega, tff
INTEGER :: npoints
NAMELIST /ProblemData/ poly_index,rho_cent,M_star
NAMELIST /AmbientData/ rhoOut, pOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut


CONTAINS

SUBROUTINE ProblemModuleInit() 
  
  INTEGER :: i, edge
  
  
    OPEN(UNIT = PROBLEM_DATA_HANDLE, FILE = 'problem.data')
    READ(PROBLEM_DATA_HANDLE, NML=ProblemData)
    READ(PROBLEM_DATA_HANDLE, NML=AmbientData)
    CLOSE(PROBLEM_DATA_HANDLE)
    rho_cent = rho_cent*hmass*Xmu ! now is g/cc
    
    CALL CreateAmbient(Ambient)
    Ambient%density=1d0
    Ambient%pressure=1d0
    Ambient%B(:)=(/BxOut,ByOut,BzOut/)
    Ambient%velocity(:)=(/vxOut, vyOut, vzOut/)
       
    CALL CreateProfile(Profile, nPoints, (/Mass_Field, P_Field/), RADIAL)
    CALL ConstructPolytropeProfile(Profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent,temp_cent)
      DO i=1,nPoints
         Profile%data(i,:)=Profile%data(i,:) / (/lScale, rScale, pScale, TempScale/) !rescale to cu
      END DO

CLOSE(PROBLEM_DATA_HANDLE)


  END SUBROUTINE ProblemModuleInit

SUBROUTINE ProblemGridInit(Info) !Specific for each info structure
    Type(InfoDef) :: Info
    INTEGER        :: i, j, k, mx, my, mz
    INTEGER        :: mbc, rmbc, zrmbc
    INTEGER        :: nx, ny, nz
    INTEGER        :: ndim
    REAL(KIND=xPrec):: x, y, z, xlower, ylower, zlower,dx, dy, dz, r,rho_profile
RETURN

! OPEN(UNIT = PROBLEM_DATA_HANDLE, FILE = 'problem.data')
!    READ(PROBLEM_DATA_HANDLE, NML=ProblemData)
! CALL ConstructPolytropeProfile(Profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent,temp_cent)
!Profile%data(:,2)=rho_profile(r)

!rmbc=levels(Info%level)%gmbc(levels(Info%level)%step)
!    mx=Info%mx(1);  xlower=Info%Xbounds(1,1) 
!    my=Info%mx(2);  ylower=Info%Xbounds(2,1)
!    mz=Info%mx(3);  zlower=Info%Xbounds(3,1)
!dx=levels(Info%level)%dx
!SELECT CASE(nDim)
   ! CASE(2)
      ! zrmbc=0
      ! mz=1
     !  zlower=0
    !   dz=0  
   ! CASE(3)
  !     zrmbc=rmbc
 !   END SELECT
!Do k=1-zrmbc, mz+zrmbc
      ! Do j=1-rmbc, my+rmbc
         ! Do i=1-rmbc, mx+rmbc
           !  x=(xlower + (REAL(i, xPrec)-half)*dx)
          !   y=(ylower + (REAL(j, xPrec)-half)*dx)
         !    z=(zlower + (REAL(k, xPrec)-half)*dx) 
        !     r = Sqrt(x**2+y**2+z**2) 
       !      r = (r*lscale)*a_scale                 ! takes r to cgs to xi 
      !        at = (2*R_outer*SQRT(PI*G*rho_cent))/xi_last            
     !        Info%q(i,j,k,1)=rho_profile(r)
    !       Info%q(i,j,k,iE) =(((Info%q(i,j,k,1)*rScale*(aT**2))/gamma)*gamma7)/pScale             
                                                             
             
!Info%q(i,j,k,iE) = (((Info%q(i,j,k,1)*rscale*(aT**2))/gamma)*gamma7)/pScaleA
   !       END DO
  !     END DO
 !   END DO
!CLOSE(PROBLEM_DATA_HANDLE)

  END SUBROUTINE ProblemGridInit
SUBROUTINE ProblemBeforeStep(INFO)

    TYPE(INFODEF) :: INFO
    INTEGER :: i

  END SUBROUTINE ProblemBeforeStep

SUBROUTINE ProblemAfterStep(INFO)
    TYPE(INFODEF) :: INFO
    ! No special after step instructions needed; this is a stub.

  END SUBROUTINE ProblemAfterStep
SUBROUTINE ProblemSetErrFlag(INFO)
    TYPE(INFODEF) :: INFO
    ! No special instructions needed; this is a stub.

  END SUBROUTINE ProblemSetErrFlag

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


END MODULE PROBLEM