| 1 | !#########################################################################
|
|---|
| 2 | !
|
|---|
| 3 | ! Copyright (C) 2003-2012 Department of Physics and Astronomy,
|
|---|
| 4 | ! University of Rochester,
|
|---|
| 5 | ! Rochester, NY
|
|---|
| 6 | !
|
|---|
| 7 | ! problem.f90 of module PlanetaryAtmosphere is part of AstroBEAR.
|
|---|
| 8 | !
|
|---|
| 9 | ! AstroBEAR is free software: you can redistribute it and/or modify
|
|---|
| 10 | ! it under the terms of the GNU General Public License as published by
|
|---|
| 11 | ! the Free Software Foundation, either version 3 of the License, or
|
|---|
| 12 | ! (at your option) any later version.
|
|---|
| 13 | !
|
|---|
| 14 | ! AstroBEAR is distributed in the hope that it will be useful,
|
|---|
| 15 | ! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 16 | ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 17 | ! GNU General Public License for more details.
|
|---|
| 18 | !
|
|---|
| 19 | ! You should have received a copy of the GNU General Public License
|
|---|
| 20 | ! along with AstroBEAR. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 21 | !
|
|---|
| 22 | !#########################################################################
|
|---|
| 23 | !> @dir PlanetaryAtmosphere
|
|---|
| 24 | !! @brief Contains files necessary for the PlanetaryAtmosphere Calculation
|
|---|
| 25 |
|
|---|
| 26 | !> @file problem.f90
|
|---|
| 27 | !! @brief Main file for module Problem
|
|---|
| 28 |
|
|---|
| 29 | !> @defgroup PlanetaryAtmosphere PlanetaryAtmosphere Module
|
|---|
| 30 | !! @brief Module for calculating collapse of a uniform cloud
|
|---|
| 31 | !! @ingroup Modules
|
|---|
| 32 |
|
|---|
| 33 | !> PlanetaryAtmosphere Module
|
|---|
| 34 | !! @ingroup PlanetaryAtmosphere
|
|---|
| 35 | MODULE Problem
|
|---|
| 36 | USE GlobalDeclarations
|
|---|
| 37 | USE DataDeclarations
|
|---|
| 38 | USE Clumps
|
|---|
| 39 | USE Ambients
|
|---|
| 40 | USE Refinements
|
|---|
| 41 | USE CommonFunctions
|
|---|
| 42 | IMPLICIT NONE
|
|---|
| 43 | SAVE
|
|---|
| 44 | PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
|
|---|
| 45 | ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
|
|---|
| 46 | PRIVATE
|
|---|
| 47 |
|
|---|
| 48 | CONTAINS
|
|---|
| 49 |
|
|---|
| 50 | SUBROUTINE ProblemModuleInit()
|
|---|
| 51 | TYPE(AmbientDef), POINTER :: Ambient
|
|---|
| 52 | TYPE(ClumpDef), POINTER :: Clump
|
|---|
| 53 | TYPE(ParticleDef), POINTER :: StellarParticle
|
|---|
| 54 | TYPE(RefinementDef), POINTER :: Refinement
|
|---|
| 55 | LOGICAL :: lFixed, lPlanet
|
|---|
| 56 | REAL(KIND=qPREC) :: Mp, Rp, a, OrbitalPeriod, RotationPeriod, Ms, Rs, Ts, rho_amb, Omega_amb, T_amb, ecc, soft_frac, P0
|
|---|
| 57 | REAL(KIND=qPREC) :: Vp, Vs, Xp, Xs
|
|---|
| 58 | REAL(KIND=qPREC) :: M, r, P, T, rho
|
|---|
| 59 | INTEGER :: i, nEntries
|
|---|
| 60 |
|
|---|
| 61 |
|
|---|
| 62 | NAMELIST/ProblemData/ Mp, Rp, a, RotationPeriod, Ms, Rs, Ts, rho_amb, T_amb, Omega_amb, lFixed, ecc, lPlanet, soft_frac, P0
|
|---|
| 63 | OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
|
|---|
| 64 | READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
|
|---|
| 65 | CLOSE(PROBLEM_DATA_HANDLE)
|
|---|
| 66 |
|
|---|
| 67 | !Rescale all parameters to computational units
|
|---|
| 68 | Mp=Mp*MJupiter/mScale
|
|---|
| 69 | Rp=Rp*RJupiter/lScale
|
|---|
| 70 | Ms=Ms*MSolar/mScale
|
|---|
| 71 | Rs=Rs*RSolar/lScale
|
|---|
| 72 | ! OrbitalPeriod=OrbitalPeriod*day/TimeScale
|
|---|
| 73 | RotationPeriod=RotationPeriod*day/TimeScale
|
|---|
| 74 | a=a*AU/lScale
|
|---|
| 75 | rho_amb=rho_amb/rScale
|
|---|
| 76 | T_amb=T_amb/TempScale
|
|---|
| 77 | Omega_amb=Omega_amb*TimeScale/day
|
|---|
| 78 | P0=P0/pscale
|
|---|
| 79 |
|
|---|
| 80 | ! Calculate Star and planet's positions based on orbital parameters.
|
|---|
| 81 | ! Rotation about z axis
|
|---|
| 82 | ! Assume start from planet and star along x axis at perihelion
|
|---|
| 83 | ! see http://astrobear.pas.rochester.edu/trac/astrobear/blog/johannjc09232013 for the derivation
|
|---|
| 84 | Vp=sqrt(ScaleGrav*(Mp+Ms)*(1d0-ecc**2)/a)/(1d0+Mp/Ms)
|
|---|
| 85 | Vs=sqrt(ScaleGrav*(Mp+Ms)*(1d0-ecc**2)/a)/(1d0+Ms/Mp)
|
|---|
| 86 | Xp=a*Ms/(Mp+Ms)
|
|---|
| 87 | Xs=a*Mp/(Mp+Ms)
|
|---|
| 88 |
|
|---|
| 89 |
|
|---|
| 90 | IF (MPI_ID == 0) THEN
|
|---|
| 91 | WRITE(*,'(4A25)') 'parameter', 'value in cu', 'physical value', 'unit'
|
|---|
| 92 | WRITE(*,'(100A1)') (/('-',i=1,100)/)
|
|---|
| 93 | 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"
|
|---|
| 94 | 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"
|
|---|
| 95 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Planet mass = ", Mp, Mp*mScale/MJupiter, " Jupiter masses"
|
|---|
| 96 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Planet radius = ", Rp, Rp*lScale/RJupiter, " Jupiter radii"
|
|---|
| 97 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Solar mass = ", Ms, Ms*mScale/MSolar, " Solar masses"
|
|---|
| 98 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Solar radius = ", Rs, Rs*lScale/RSolar, " Solar radii"
|
|---|
| 99 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Rotation Period = ", RotationPeriod, RotationPeriod*TimeScale/day, " days"
|
|---|
| 100 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Ambient density = ", rho_amb, rho_amb*rScale, " g/cc"
|
|---|
| 101 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Ambient Temperature = ", T_amb, T_amb*TempScale, " K"
|
|---|
| 102 | WRITE(*,'(A25,E25.15,F25.15,A25)') "Semi Major Axis = ", a, a*lScale/AU, " AU"
|
|---|
| 103 | WRITE(*,'(100A1)') (/('-',i=1,100)/)
|
|---|
| 104 | 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,")"
|
|---|
| 105 | 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,")"
|
|---|
| 106 |
|
|---|
| 107 |
|
|---|
| 108 | END IF
|
|---|
| 109 |
|
|---|
| 110 |
|
|---|
| 111 | ! Create Particle for star
|
|---|
| 112 | IF (.NOT. lRestart) THEN
|
|---|
| 113 | CALL CreateParticle(StellarParticle)
|
|---|
| 114 | StellarParticle%q(1)=Ms
|
|---|
| 115 | StellarParticle%xloc=(/-Xs,0d0,0d0/)
|
|---|
| 116 | StellarParticle%lFixed=lFixed
|
|---|
| 117 | StellarParticle%buffer=-1
|
|---|
| 118 | IF (.NOT. lFixed) StellarParticle%q(imom(1:3))=(/0d0,-Vs,0d0/) - Cross3D((/0d0,0d0,OmegaRot/), StellarParticle%xloc) !Adjust local velocity for rotation.
|
|---|
| 119 |
|
|---|
| 120 | CALL UpdateParticle(StellarParticle) !creates point gravity object
|
|---|
| 121 | StellarParticle%PointGravityObj%soft_function = SPLINESOFT
|
|---|
| 122 | StellarParticle%PointGravityObj%soft_length = soft_frac*(Xp+Xs)
|
|---|
| 123 | CALL AddSinkParticle(StellarParticle)
|
|---|
| 124 | ELSE
|
|---|
| 125 | StellarParticle=>SinkParticles%self
|
|---|
| 126 | END IF
|
|---|
| 127 |
|
|---|
| 128 | CALL CreateAmbient(Ambient)
|
|---|
| 129 | Ambient%density=rho_amb
|
|---|
| 130 | Ambient%pressure=rho_amb*T_amb
|
|---|
| 131 | ! Ambient%PointGravity=>StellarParticle%PointGravityObj !this will automatically create hydrostatic ambient
|
|---|
| 132 | Ambient%Omega=(/0d0,0d0,Omega_amb-OmegaRot/) !Shift ambient omega by frame rotation
|
|---|
| 133 | Ambient%PersistInBoundaries=.true.
|
|---|
| 134 | CALL UpdateAmbient(Ambient)
|
|---|
| 135 |
|
|---|
| 136 |
|
|---|
| 137 | IF (lPlanet) THEN
|
|---|
| 138 | CALL CreateClump(Clump)
|
|---|
| 139 | Clump%position=(/Xp,0d0,0d0/)
|
|---|
| 140 | Clump%omega=2d0*Pi/RotationPeriod-OmegaRot !Adjust planet rotation by frame rotation
|
|---|
| 141 | IF (.NOT. lFixed) Clump%velocity=(/0d0,Vp,0d0/) - Cross3D((/0d0,0d0,OmegaRot/), Clump%position) !Adjust local velocity for rotation.
|
|---|
| 142 |
|
|---|
| 143 | OPEN (UNIT=PROBLEM_DATA_HANDLE, FILE='planet_profile.dat', STATUS="OLD")
|
|---|
| 144 | READ (PROBLEM_DATA_HANDLE, *) nEntries
|
|---|
| 145 | READ (PROBLEM_DATA_HANDLE, *) !Skip the header column
|
|---|
| 146 | CALL CreateProfile(Clump%Profile, nEntries, (/Mass_Field, P_Field/), RADIAL) !Mass field - is gas density
|
|---|
| 147 | DO i=1,nEntries
|
|---|
| 148 | READ(PROBLEM_DATA_HANDLE, *) M, r, P, T, rho
|
|---|
| 149 | Clump%Profile%data(i,:)=(/ r/lScale, rho/rScale, P/pScale /)
|
|---|
| 150 | END DO
|
|---|
| 151 | ! CALL Profile_SelfGravityHSE(Clump%Profile, P0) !Use internal pressure from profile data file to recalculate pressure profile in HSE
|
|---|
| 152 | CALL UpdateProfile(Clump%Profile)
|
|---|
| 153 | ! write(*,'(3E25.15)') transpose(clump%profile%data)
|
|---|
| 154 | Clump%radius=Clump%Profile%data(nEntries,1)
|
|---|
| 155 |
|
|---|
| 156 | Clump%thickness=0d0
|
|---|
| 157 | Clump%subsample=1
|
|---|
| 158 | CALL UpdateClump(Clump)
|
|---|
| 159 |
|
|---|
| 160 |
|
|---|
| 161 | CALL ClearAllRefinements()
|
|---|
| 162 | CALL CreateRefinement(Refinement)
|
|---|
| 163 | Refinement%Type=REFINE_INSIDE
|
|---|
| 164 | CALL CreateShape(Refinement%Shape,SPHERE,1.1*(/Rp,Rp,Rp/))
|
|---|
| 165 | Refinement%Shape%position=Clump%position
|
|---|
| 166 | CALL UpdateShape(Refinement%Shape)
|
|---|
| 167 | CALL UpdateRefinement(Refinement)
|
|---|
| 168 |
|
|---|
| 169 |
|
|---|
| 170 | CALL CreateRefinement(Refinement)
|
|---|
| 171 | Refinement%Type=DEREFINE_INSIDE
|
|---|
| 172 | CALL CreateShape(Refinement%Shape,SPHERE,.6*(/Rp,Rp,Rp/))
|
|---|
| 173 | Refinement%Shape%position=Clump%position
|
|---|
| 174 | CALL UpdateShape(Refinement%Shape)
|
|---|
| 175 | CALL UpdateRefinement(Refinement)
|
|---|
| 176 |
|
|---|
| 177 |
|
|---|
| 178 | END IF
|
|---|
| 179 | END SUBROUTINE ProblemModuleInit
|
|---|
| 180 |
|
|---|
| 181 | SUBROUTINE ProblemGridInit(Info)
|
|---|
| 182 | TYPE(InfoDef) :: Info
|
|---|
| 183 | END SUBROUTINE ProblemGridInit
|
|---|
| 184 |
|
|---|
| 185 | SUBROUTINE ProblemBeforeStep(Info)
|
|---|
| 186 | TYPE(InfoDef) :: Info
|
|---|
| 187 | END SUBROUTINE ProblemBeforeStep
|
|---|
| 188 |
|
|---|
| 189 | SUBROUTINE ProblemAfterStep(Info)
|
|---|
| 190 | TYPE(InfoDef) :: Info
|
|---|
| 191 | END SUBROUTINE ProblemAfterStep
|
|---|
| 192 |
|
|---|
| 193 | SUBROUTINE ProblemSetErrFlag(Info)
|
|---|
| 194 | TYPE(InfoDef) :: Info
|
|---|
| 195 | END SUBROUTINE ProblemSetErrFlag
|
|---|
| 196 |
|
|---|
| 197 | SUBROUTINE ProblemBeforeGlobalStep(n)
|
|---|
| 198 | INTEGER :: n
|
|---|
| 199 | END SUBROUTINE ProblemBeforeGlobalStep
|
|---|
| 200 |
|
|---|
| 201 | END MODULE
|
|---|
