planet2d.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: planet2d.f90 #
5!# #
6!# Copyright (C) 2006-2019 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
9!# Jubin Lirawi <jlirawi@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
27
28!----------------------------------------------------------------------------!
38!----------------------------------------------------------------------------!
39PROGRAM planet2d
40 USE fosite_mod
43 USE constants_si_mod, ONLY : c, gn, kb, na, sb, ke, au, &
45 USE common_dict
46#include "tap.h"
47 IMPLICIT NONE
48 !--------------------------------------------------------------------------!
49 ! general constants ! !
50! REAL, PARAMETER :: GN = 6.6742D-11 ! Newtons grav. constant !
51! REAL, PARAMETER :: CC = 2.99792458D+8 ! speed of light [m/s] !
52! REAL, PARAMETER :: RG = 8.31447 ! molar gas constant !
53! REAL, PARAMETER :: SBconst = 5.670367D-8 ! Stefan-Boltzmann const !
54! REAL, PARAMETER :: AU = 1.49597870691E+11 ! astronomical unit [m] !
55! REAL, PARAMETER :: YEAR = 3.15576E+7 ! Julian year [s] !
56! REAL, PARAMETER :: DAY = 8.6400E+4 ! Solar Day [s] !
57! REAL, PARAMETER :: REARTH = 6.371E+6 ! radius of earth [m] !
58! REAL, PARAMETER :: RJUP = 71.492E+6 ! radius of jupiter [m] !
59! REAL, PARAMETER :: RSUN = 6.957E+8 ! radius of the sun [m] !
60! REAL, PARAMETER :: MEARTH = 5.9723D+24 ! mass of the earth [kg] !
61 ! planetary parameters ! !
62 REAL, PARAMETER :: tyear = 4.05 * day ! tropical year [s] !
63 REAL, PARAMETER :: rplanet = 0.788*rearth ! planetary radius [m] !
64 REAL, PARAMETER :: theta0 = 0.0 ! axis-plane-angle [rad] !
65 REAL, PARAMETER :: phi0 = 0.0 ! [rad] !
66 REAL, PARAMETER :: omega = 2*pi/tyear/8 ! ang. rotation [rad/s] !
67 REAL, PARAMETER :: fsun = omega/(2*pi) - & ! freq.(!) day-night[/s] !
68 1.0/tyear ! !
69 REAL, PARAMETER :: mass = 0.297*mearth ! mass of the planet [kg]!
70 REAL, PARAMETER :: gacc = gn*mass/rplanet**2! grav. accel. [m/s**2] !
71 ! orbital parameters ! !
72 REAL, PARAMETER :: sm_axis = 0.02219*au ! semi major axis [m] !
73 REAL, PARAMETER :: excent = 0.0 ! numerical eccentricity !
74 ! stellar parameters ! !
75 REAL, PARAMETER :: rstar = 0.121*rsun ! radius of the star [m] !
76 REAL, PARAMETER :: tstar = 2511 ! eff. temperature [K] !
77 REAL, PARAMETER :: fstar = (rstar/au)**2 & ! rad flux density @ 1 AU!
78 * sb*tstar**4 ! [W/m**2] !
79 ! simulation parameters ! !
80 REAL, PARAMETER :: tsim = 10 * tyear ! simulation stop time !
81 INTEGER, PARAMETER :: xres = 1 ! radius-resolution [px] !
82 INTEGER, PARAMETER :: yres = 16 ! theta-resolution [px] !
83 INTEGER, PARAMETER :: zres = 32 ! phi-resolution [px] !
84 ! gas properties of the atmosphere
85 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats!
86 REAL, PARAMETER :: p0 = 1.0e7 ! surf. press. [Pa] !
87 REAL, PARAMETER :: mu = 2.8586e-2 ! molar mass [kg/mol] !
88 REAL, PARAMETER :: t0 = 258.1 ! mean equil. temp. [K] !
89 ! optical properties of the atmosphere
90 REAL, PARAMETER :: albedo = 0.306 ! albedo of the planet !
91 ! output parameters ! !
92 INTEGER, PARAMETER :: onum = 10 ! num. output data sets !
93 CHARACTER(LEN=256), PARAMETER & ! output data dir !
94 :: odir = './' ! !
95 CHARACTER(LEN=256), PARAMETER & ! output data file name !
96 :: ofname = 'planet2d' ! !
97 !--------------------------------------------------------------------------!
98 CLASS(fosite), ALLOCATABLE :: sim
99 CLASS(sources_base), POINTER :: sp => null()
100 CLASS(sources_planetcooling), POINTER :: spcool => null()
101 REAL :: tmean
102 LOGICAL :: ok
103 !--------------------------------------------------------------------------!
104
105 tap_plan(1)
106
107 ALLOCATE(sim)
108 CALL sim%InitFosite()
109 CALL makeconfig(sim, sim%config)
110 CALL sim%Setup()
111 ! get pointer on planetary cooling source term
112 IF (ALLOCATED(sim%Sources)) THEN
113 sp => sim%Sources%GetSourcesPointer(planet_cooling)
114 IF (ASSOCIATED(sp)) THEN
115 SELECT TYPE(sp)
116 CLASS IS(sources_planetcooling)
117 spcool => sp
118 CLASS DEFAULT
119 NULLIFY(spcool)
120 END SELECT
121 END IF
122 END IF
123 IF (.NOT.ASSOCIATED(spcool)) &
124 CALL sim%Error("planet2d","planetary cooling sources not initialized")
125 ! set initial condition
126 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc%pvar, sim%Timedisc%cvar)
127 ! run simulation
128 CALL sim%Run()
129 ! compute final global mean surface temperature using local surface temperatures
130 ! computed by the planetary cooling module
131 tmean = sum(spcool%T_s%data1d(:)*sim%Mesh%volume%data1d(:),mask=sim%Mesh%without_ghost_zones%mask1d(:)) &
132 / sum(sim%Mesh%volume%data1d(:),mask=sim%Mesh%without_ghost_zones%mask1d(:))
133 CALL sim%Finalize()
134 DEALLOCATE(sim)
135
136 tap_check_close(t0,tmean,t0*0.01,"deviation of final mean surface temperature from initial T0 < 1%")
137 tap_done
138
139CONTAINS
140
141 SUBROUTINE makeconfig(Sim, config)
142 IMPLICIT NONE
143 !------------------------------------------------------------------------!
144 CLASS(fosite) :: Sim
145 TYPE(dict_typ), POINTER :: config
146 !------------------------------------------------------------------------!
147 ! Local variable declaration
148 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
149 timedisc, fluxes, heating, sources, &
150 grav, pmass, cooling, rotframe
151 !------------------------------------------------------------------------!
152 !mesh settings
153 mesh => dict( &
154 "meshtype" / midpoint, &
155 "geometry" / spherical_planet, &
156 "decomposition" / (/ 1, -1, 1/), &
157 "omega" / omega, &
158 "inum" / yres, &
159 "jnum" / zres, &
160 "knum" / xres, &
161 "xmin" / (0.05), &
162 "xmax" / (pi-0.05), &
163 "ymin" / (0.0), &
164 "ymax" / (2.0*pi), &
165 "zmin" / rplanet, &
166 "zmax" / rplanet, &
167 "gparam" / rplanet, &
168 "output/rotation" / 0, &
169 "output/volume" / 1, &
170 "output/dAz" / 1)
171
172 ! boundary conditions
173 boundary => dict( &
174 "western" / reflecting, &
175 "eastern" / reflecting, &
176 "southern" / periodic, &
177 "northern" / periodic, &
178 "bottomer" / reflecting, &
179 "topper" / reflecting)
180
181 ! physics settings
182 physics => dict( &
183 "problem" / euler, &
184 "mu" / mu, &
185 "gamma" / gamma)
186
187 ! flux calculation and reconstruction method
188 fluxes => dict( &
189 "fluxtype" / kt, &
190 "order" / linear, &
191! "variables" / CONSERVATIVE, &
192 "variables" / primitive, &
193 "limiter" / vanleer)
194
195 ! rotating frame with disabled centrifugal forces
196 ! Rotating planets are usually in an equilibrium state in which
197 ! gravitational and centrifugal forces balance, so that the planet
198 ! becomes a spheroid instead of a sphere. On the spheroid the
199 ! tangential components of the centrifugal forces vanish and
200 ! one can account for the vertical component defining an
201 ! effective gravitational acceleration (in 3D simulations).
202 rotframe => dict( &
203 "stype" / rotating_frame, &
204 "disable_centaccel" / 1)
205
206 ! cooling in infrared
207 cooling => dict( &
208 "stype" / planet_cooling, &
209 "output/Qcool" / 1, &
210 "output/T_s" / 1, &
211 "output/RHO_s" / 1, &
212 "output/P_s" / 1, &
213 "radflux" / fstar, &
214 "albedo" / albedo, &
215 "mean_distance" / (sm_axis*(1. + 0.5*excent**2)), &
216 "T_0" / t0, &
217 "cvis" / 0.1, &
218 "gacc" / gacc)
219
220 ! heating by a star
221 heating => dict( &
222 "stype" / planet_heating, &
223 "output/Qstar" / 1, &
224 "year" / tyear, &
225 "theta0" / theta0, &
226 "phi0" / phi0, &
227 "omegasun" / fsun, &
228 "albedo" / albedo, &
229 "radflux" / fstar,&
230 "semimajoraxis" / sm_axis, &
231 "excentricity" / excent, &
232 "cvis" / 0.1)
233
234 sources => dict( &
235 "rotframe" / rotframe, &
236 "cooling" / cooling, &
237 "heating" / heating)
238
239 ! time discretization settings
240 timedisc => dict( &
241 "method" / ssprk, &
242 "cfl" / 0.4, &
243 "stoptime" / tsim, &
244 "tol_rel" / 0.0095, &
245 "dtlimit" / 1.0e-15, &
246 "maxiter" / 1000000000)
247
248 datafile => dict(&
249 "fileformat" / vtk, &
250 "filename" / (trim(odir) // trim(ofname)), &
251 "count" / onum)
252
253 config => dict( &
254 "mesh" / mesh, &
255 "physics" / physics, &
256 "boundary" / boundary, &
257 "fluxes" / fluxes, &
258 "timedisc" / timedisc, &
259 "sources" / sources, &
260 "datafile" / datafile)
261 END SUBROUTINE makeconfig
262
263
264 SUBROUTINE initdata(Mesh,Physics,pvar,cvar)
265 IMPLICIT NONE
266 !------------------------------------------------------------------------!
267 CLASS(physics_base), INTENT(IN) :: Physics
268 CLASS(mesh_base), INTENT(IN) :: Mesh
269 CLASS(marray_compound), POINTER, INTENT(INOUT) :: pvar,cvar
270 !------------------------------------------------------------------------!
271 ! Local variable declaration
272 !------------------------------------------------------------------------!
273 SELECT TYPE(p => pvar)
274 CLASS IS(statevector_euler)
275 ! assuming hydrostatic equillibrium with constant (effective) gravitational
276 ! acceleration yields constant column density, i.e. vertically integrated density,
277 ! given by the ratio of mean surface pressure P0 and (effective) grav. acceleration GACC
278 p%density%data1d(:) = p0/gacc
279
280 ! assuming hydrostic equillibrium with mean surface
281 ! temperature T0 yields the initial condition for the 2D vertically
282 ! integrated pressure
283 p%pressure%data1d(:) = gamma*physics%Constants%RG / mu * t0 * p%density%data1d(:)
284
285 ! vanishing velocities in comoving frame
286 p%velocity%data1d(:) = 0.0
287 CLASS DEFAULT
288 CALL mesh%Error("planet2d::InitData","state vector not supported")
289 END SELECT
290
291 CALL physics%Convert2Conservative(pvar,cvar)
292 CALL mesh%Info(" DATA-----> initial condition: 2D planetary atmosphere")
293
294 END SUBROUTINE initdata
295
296END PROGRAM planet2d
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
Dictionary for generic data types.
Definition: common_dict.f90:61
type(dict_typ) function, pointer, public dict(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16, n17, n18, n19, n20)
Construct a new dictionary from several key/value pairs. Together with the Assign subroutine and over...
module for SI units and physical constants
real, parameter, public mjupiter
jupiter mass [kg]
real, parameter, public mearth
earth mass [kg]
real, parameter, public gn
gravitational constant [m^3/kg/s^2]
real, parameter, public rsun
sun radius [m]
real, parameter, public ke
electron scattering opacity [m^2/kg]
real, parameter, public c
vacuum speed of light [m/s]
real, parameter, public day
length of a day [s]
real, parameter, public sb
Stefan-Boltzmann constant [W/m^2/K^4].
real, parameter, public au
astronomical unit [m]
real, parameter, public na
Avogadro constant [1/mol].
real, parameter, public kb
Boltzmann constant [J/K].
real, parameter, public pi
real, parameter, public rjupiter
jupiter radius [m]
real, parameter, public rearth
earth radius [m]
real, parameter, public msun
sun mass [kg]
generic source terms module providing functionaly common to all source terms
gray cooling of planetary atmospheres
program planet2d
Definition: planet2d.f90:39
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
main fosite class
Definition: fosite.f90:71