!######################################################################### ! ! Copyright (C) 2003-2012 Department of Physics and Astronomy, ! University of Rochester, ! Rochester, NY ! ! problem.f90 of module PlanetaryAtmosphere 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 . ! !######################################################################### !> @dir PlanetaryAtmosphere !! @brief Contains files necessary for the PlanetaryAtmosphere Calculation !> @file problem.f90 !! @brief Main file for module Problem !> @defgroup PlanetaryAtmosphere PlanetaryAtmosphere Module !! @brief Module for calculating collapse of a uniform cloud !! @ingroup Modules !> PlanetaryAtmosphere Module !! @ingroup PlanetaryAtmosphere MODULE Problem USE GlobalDeclarations USE DataDeclarations USE Clumps USE Ambients USE Refinements USE CommonFunctions IMPLICIT NONE SAVE PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, & ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep PRIVATE CONTAINS SUBROUTINE ProblemModuleInit() TYPE(AmbientDef), POINTER :: Ambient TYPE(ClumpDef), POINTER :: Clump TYPE(ParticleDef), POINTER :: StellarParticle TYPE(RefinementDef), POINTER :: Refinement LOGICAL :: lFixed, lPlanet REAL(KIND=qPREC) :: Mp, Rp, a, OrbitalPeriod, RotationPeriod, Ms, Rs, Ts, rho_amb, Omega_amb, T_amb, ecc, soft_frac, P0 REAL(KIND=qPREC) :: Vp, Vs, Xp, Xs REAL(KIND=qPREC) :: M, r, P, T, rho INTEGER :: i, nEntries NAMELIST/ProblemData/ Mp, Rp, a, RotationPeriod, Ms, Rs, Ts, rho_amb, T_amb, Omega_amb, lFixed, ecc, lPlanet, soft_frac, P0 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD") READ(PROBLEM_DATA_HANDLE,NML=ProblemData) CLOSE(PROBLEM_DATA_HANDLE) !Rescale all parameters to computational units Mp=Mp*MJupiter/mScale Rp=Rp*RJupiter/lScale Ms=Ms*MSolar/mScale Rs=Rs*RSolar/lScale ! OrbitalPeriod=OrbitalPeriod*day/TimeScale RotationPeriod=RotationPeriod*day/TimeScale a=a*AU/lScale rho_amb=rho_amb/rScale T_amb=T_amb/TempScale Omega_amb=Omega_amb*TimeScale/day P0=P0/pscale ! Calculate Star and planet's positions based on orbital parameters. ! Rotation about z axis ! Assume start from planet and star along x axis at perihelion ! see http://astrobear.pas.rochester.edu/trac/astrobear/blog/johannjc09232013 for the derivation Vp=sqrt(ScaleGrav*(Mp+Ms)*(1d0-ecc**2)/a)/(1d0+Mp/Ms) Vs=sqrt(ScaleGrav*(Mp+Ms)*(1d0-ecc**2)/a)/(1d0+Ms/Mp) Xp=a*Ms/(Mp+Ms) Xs=a*Mp/(Mp+Ms) IF (MPI_ID == 0) THEN WRITE(*,'(4A25)') 'parameter', 'value in cu', 'physical value', 'unit' WRITE(*,'(100A1)') (/('-',i=1,100)/) WRITE(*,'(A25,E25.15,F25.15,A25)') "Orbit Period = ", 2*Pi*sqrt(a**3/ScaleGrav/(Ms+Mp)), 2*Pi*sqrt(a**3/ScaleGrav/(Ms+Mp)) * TimeScale/day, " days" WRITE(*,'(A25,E25.15,F25.15,A25)') "Orbital Omega = ",sqrt(ScaleGrav*(Ms+Mp)/a**3), sqrt(ScaleGrav*(Ms+Mp)/a**3)/TimeScale*day, " radians/day" WRITE(*,'(A25,E25.15,F25.15,A25)') "Planet mass = ", Mp, Mp*mScale/MJupiter, " Jupiter masses" WRITE(*,'(A25,E25.15,F25.15,A25)') "Planet radius = ", Rp, Rp*lScale/RJupiter, " Jupiter radii" WRITE(*,'(A25,E25.15,F25.15,A25)') "Solar mass = ", Ms, Ms*mScale/MSolar, " Solar masses" WRITE(*,'(A25,E25.15,F25.15,A25)') "Solar radius = ", Rs, Rs*lScale/RSolar, " Solar radii" WRITE(*,'(A25,E25.15,F25.15,A25)') "Rotation Period = ", RotationPeriod, RotationPeriod*TimeScale/day, " days" WRITE(*,'(A25,E25.15,F25.15,A25)') "Ambient density = ", rho_amb, rho_amb*rScale, " g/cc" WRITE(*,'(A25,E25.15,F25.15,A25)') "Ambient Temperature = ", T_amb, T_amb*TempScale, " K" WRITE(*,'(A25,E25.15,F25.15,A25)') "Semi Major Axis = ", a, a*lScale/AU, " AU" WRITE(*,'(100A1)') (/('-',i=1,100)/) WRITE(*,'(A25,A1,F5.2,A1,F5.2,A1,F5.2,A1,A25,A1,F5.2,A1,F5.2,A1,F5.2,A1)') "Planet located at ", "(", Xp,",",0d0,",",0d0,")", "with velocity ", "(", 0d0,",",Vp,",",0d0,")" WRITE(*,'(A25,A1,F5.2,A1,F5.2,A1,F5.2,A1,A25,A1,F5.2,A1,F5.2,A1,F5.2,A1)') "Star located at ", "(", -Xs,",",0d0,",",0d0,")", "with velocity ", "(", 0d0,",",-Vs,",",0d0,")" END IF ! Create Particle for star IF (.NOT. lRestart) THEN CALL CreateParticle(StellarParticle) StellarParticle%q(1)=Ms StellarParticle%xloc=(/-Xs,0d0,0d0/) StellarParticle%lFixed=lFixed StellarParticle%buffer=-1 IF (.NOT. lFixed) StellarParticle%q(imom(1:3))=(/0d0,-Vs,0d0/) - Cross3D((/0d0,0d0,OmegaRot/), StellarParticle%xloc) !Adjust local velocity for rotation. CALL UpdateParticle(StellarParticle) !creates point gravity object StellarParticle%PointGravityObj%soft_function = SPLINESOFT StellarParticle%PointGravityObj%soft_length = soft_frac*(Xp+Xs) CALL AddSinkParticle(StellarParticle) ELSE StellarParticle=>SinkParticles%self END IF CALL CreateAmbient(Ambient) Ambient%density=rho_amb Ambient%pressure=rho_amb*T_amb ! Ambient%PointGravity=>StellarParticle%PointGravityObj !this will automatically create hydrostatic ambient Ambient%Omega=(/0d0,0d0,Omega_amb-OmegaRot/) !Shift ambient omega by frame rotation Ambient%PersistInBoundaries=.true. CALL UpdateAmbient(Ambient) IF (lPlanet) THEN CALL CreateClump(Clump) Clump%position=(/Xp,0d0,0d0/) Clump%omega=2d0*Pi/RotationPeriod-OmegaRot !Adjust planet rotation by frame rotation IF (.NOT. lFixed) Clump%velocity=(/0d0,Vp,0d0/) - Cross3D((/0d0,0d0,OmegaRot/), Clump%position) !Adjust local velocity for rotation. OPEN (UNIT=PROBLEM_DATA_HANDLE, FILE='planet_profile.dat', STATUS="OLD") READ (PROBLEM_DATA_HANDLE, *) 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, *) M, r, P, T, rho Clump%Profile%data(i,:)=(/ r/lScale, rho/rScale, P/pScale /) END DO ! CALL Profile_SelfGravityHSE(Clump%Profile, P0) !Use internal pressure from profile data file to recalculate pressure profile in HSE CALL UpdateProfile(Clump%Profile) ! write(*,'(3E25.15)') transpose(clump%profile%data) Clump%radius=Clump%Profile%data(nEntries,1) Clump%thickness=0d0 Clump%subsample=1 CALL UpdateClump(Clump) CALL ClearAllRefinements() CALL CreateRefinement(Refinement) Refinement%Type=REFINE_INSIDE CALL CreateShape(Refinement%Shape,SPHERE,1.1*(/Rp,Rp,Rp/)) Refinement%Shape%position=Clump%position CALL UpdateShape(Refinement%Shape) CALL UpdateRefinement(Refinement) CALL CreateRefinement(Refinement) Refinement%Type=DEREFINE_INSIDE CALL CreateShape(Refinement%Shape,SPHERE,.6*(/Rp,Rp,Rp/)) Refinement%Shape%position=Clump%position CALL UpdateShape(Refinement%Shape) CALL UpdateRefinement(Refinement) END IF END SUBROUTINE ProblemModuleInit SUBROUTINE ProblemGridInit(Info) TYPE(InfoDef) :: Info END SUBROUTINE ProblemGridInit SUBROUTINE ProblemBeforeStep(Info) TYPE(InfoDef) :: Info END SUBROUTINE ProblemBeforeStep SUBROUTINE ProblemAfterStep(Info) TYPE(InfoDef) :: Info END SUBROUTINE ProblemAfterStep SUBROUTINE ProblemSetErrFlag(Info) TYPE(InfoDef) :: Info END SUBROUTINE ProblemSetErrFlag SUBROUTINE ProblemBeforeGlobalStep(n) INTEGER :: n END SUBROUTINE ProblemBeforeGlobalStep END MODULE