agndisk.f90
Program and data initialization for 3D rotating flow.
Program and data initialization for 3D rotating flow initial condition:
hydrostatic equilibrium for slightly non-keplerian flow deviation from Keplerian:
\begin{eqnarray*} v_{\varphi}(s,z) &=& \frac{v_{\mathrm{kep}}(s,z)}{\lambda (z)} \\ \lambda(z) &=& 1 + \frac{A_0}{1+(z/Z_0)^2} \\ P_\mathrm{c}(s) &=& P_{\mathrm{0}}\exp{(-s/S_0)}, \end{eqnarray*}
with \( P_\mathrm{c}(s) \) the pressure in equatorial plane and constant \( A_0 = 0.1 \), \( Z_0 = 300 \, L_{\mathrm{scale}} \), \( S_0 = 150 L_{\mathrm{scale}} \).
1!#############################################################################
2!# #
3!# fosite - 2D hydrodynamical simulation program #
4!# module: agndisk.f90 #
5!# #
6!# Copyright (C) 2006-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program 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 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
25
26!----------------------------------------------------------------------------!
43!----------------------------------------------------------------------------!
46#ifdef PARALLEL
47#ifdef HAVE_MPI_MOD
48 USE mpi
49#endif
50#endif
51 IMPLICIT NONE
52#ifdef PARALLEL
53#ifdef HAVE_MPIF_H
54 include 'mpif.h'
55#endif
56#endif
57 !--------------------------------------------------------------------------!
58 ! general constants
59 REAL, PARAMETER :: GN = 6.6742d-11 ! Newtons grav. constant !
60 REAL, PARAMETER :: CC = 2.99792458d+8 ! speed of light !
61 REAL, PARAMETER :: MSUN = 1.989d+30 ! solar mass [kg] !
62 REAL, PARAMETER :: AU = 1.49597870691e+11 ! astronomical unit [m] !
63 REAL, PARAMETER :: PARSEC = 0.5/pi*1.296e+6 * au ! parsec [m] !
64 REAL, PARAMETER :: YEAR = 3.15576e+7 ! Julian year [sec] !
65 REAL, PARAMETER :: ETA_ACC = 0.1 ! accretion efficiency !
66 ! simulation parameters
67 REAL, PARAMETER :: TSIM = 1.0e+4*year ! simulation time !
68 REAL, PARAMETER :: M_BH = 1.0e+6 * msun ! central point mass !
69 REAL, PARAMETER :: M_DISK = 1.0e+0 * m_bh ! accretion disk mass !
70 REAL, PARAMETER :: RS = 2 * gn * m_bh / cc**2 ! Schwarzschild radius !
71 REAL, PARAMETER :: LSCALE = 1.0e+3 * rs ! typical length scale !
72 REAL, PARAMETER :: ETASGS = 1.0e-6
73 REAL, PARAMETER :: gamma = 1.4
74 ! computational domain
75 REAL, PARAMETER :: RMIN = 5.0d+0 * lscale
76 REAL, PARAMETER :: RMAX = 3.0d+2 * lscale
77 REAL, PARAMETER :: ZMIN = -1.5d+2 * lscale
78 REAL, PARAMETER :: ZMAX = 1.5d+2 * lscale
79 REAL, PARAMETER :: ZTAN = 1.0e+2 * lscale ! scale for tancyl. coord.!
80 ! initial condition:
81 !! hydrostatic equilibrium for slightly non-keplerian flow
82 !! deviation from Keplerian: v_phi(s,z) = v_kep(s,z) / lambda(z)
83 !! lambda(z) = 1 + A0/(1+(z/Z0)**2)
84 !! pressure in equatorial plane: Pc(s) = P0 * EXP(-s/S0)
85 REAL, PARAMETER :: Z0 = 3.0e+02 * lscale ! vertical scale !
86 REAL, PARAMETER :: A0 = 1.0e-01 ! velocity deviation !
87 REAL, PARAMETER :: S0 = 1.5e+02 * lscale ! radial scale !
88 REAL, PARAMETER :: B0 = -1.0 ! <1 ! ! slope den. power law !
89 REAL, PARAMETER :: RHO_INF = 1.0e-15 ! den. at infinity !
90 REAL, PARAMETER :: P_INF = 1.4e-10 ! pressure at infinity !
91 INTEGER, PARAMETER :: PHYS = euler ! transport model !
92 ! resolution
93 INTEGER, PARAMETER :: MGEO = cylindrical
94 INTEGER, PARAMETER :: RRES = 100 ! res. in r-direction !
95 INTEGER, PARAMETER :: PHIRES = 1 ! res. in phi-direction !
96 INTEGER, PARAMETER :: ZRES = 101 ! res. in z-direction !
97 ! output file parameter
98 INTEGER, PARAMETER :: ONUM = 1000 ! number of outputs !
99 CHARACTER(LEN=256), PARAMETER :: ODIR &! output directory
100 = "./"
101 CHARACTER(LEN=256), PARAMETER & ! output data file name !
102 :: OFNAME = 'agndisk'
103 !--------------------------------------------------------------------------!
104 CLASS(fosite), ALLOCATABLE :: Sim
105 !--------------------------------------------------------------------------!
106
107 ALLOCATE(sim)
108
109 CALL sim%InitFosite()
110
112
113 CALL sim%Setup()
114
115 ! set initial condition
117
118 CALL sim%Run()
119
120 CALL sim%Finalize()
121
122CONTAINS
123
125 IMPLICIT NONE
126 !------------------------------------------------------------------------!
127 TYPE(Dict_TYP),POINTER :: config
128 !------------------------------------------------------------------------!
129 ! Local variable declaration
130 TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, logfile, sources, &
131 timedisc,fluxes,stemp, grav
132 !------------------------------------------------------------------------!
133 ! mesh settings
134 SELECT CASE(mgeo)
135 CASE(cylindrical)
136 mesh => dict( &
137 "meshtype" / midpoint, &
138 "geometry" / mgeo, &
139 "inum" / rres, &
140 "jnum" / phires, &
141 "knum" / zres, &
142 "xmin" / rmin, &
143 "xmax" / rmax, &
144 "ymin" / (0.0), &
145 "ymax" / (2.0*pi), &
146 "zmin" / zmin, &
147 "zmax" / zmax)
148
149 ! boundary conditions
150 boundary => dict( &
151! "western" / REFLECTING, &
152! "eastern" / REFLECTING, &
153 "western" / custom, &
154! "eastern" / CUSTOM, &
155! "western" / NO_GRADIENTS, &
156! "eastern" / NO_GRADIENTS, &
157! "western" / FARFIELD, &
158 "eastern" / farfield, &
159 "bottomer" / farfield, &
160 "topper" / farfield, &
161! "southern" / AXIS, &
162 "southern" / periodic, &
163! "southern" / REFLECTING, &
164 "northern" / periodic)!, &
165! "northern" / REFLECTING)
166! "bottomer" / NO_GRADIENTS, &
167! "northern" / CUSTOM)
168! "southern" / NO_GRADIENTS, &
169! "topper" / NO_GRADIENTS)
170! "southern" / FARFIELD, &
171! "northern" / FARFIELD)
172! "northern" / REFLECTING)
173
174 CASE(spherical)
175 mesh => dict( &
176 "meshtype" / midpoint, &
177 "geometry" / mgeo, &
178 "inum" / rres, &
179 "jnum" / zres, &
180 "knum" / phires, &
181 "xmin" / rmin, &
182 "xmax" / rmax, &
183 "ymin" / (0.0), &
184 "ymax" / (pi), &
185 "zmin" / (0.0), &
186 "zmax" / (2.0*pi))
187
188 ! boundary conditions
189 boundary => dict( &
190! "western" / REFLECTING, &
191! "eastern" / REFLECTING, &
192 "western" / custom, &
193! "eastern" / CUSTOM, &
194! "western" / NO_GRADIENTS, &
195! "eastern" / NO_GRADIENTS, &
196! "western" / FARFIELD, &
197 "eastern" / farfield, &
198 "bottomer" / periodic, &
199 "topper" / periodic, &
200! "southern" / AXIS, &
201 "southern" / axis, &
202! "southern" / REFLECTING, &
203 "northern" / axis)!, &
204! "northern" / REFLECTING)
205! "bottomer" / NO_GRADIENTS, &
206! "northern" / CUSTOM)
207! "southern" / NO_GRADIENTS, &
208! "topper" / NO_GRADIENTS)
209! "southern" / FARFIELD, &
210! "northern" / FARFIELD)
211! "northern" / REFLECTING)
212
213
214 END SELECT
215
216
217 ! physics settings
218 physics => dict( &
219 "problem" / phys, &
220 "gamma" / gamma, & ! ratio of specific heats !
221 "mu" / 6.02e-4, & ! mean molecular weight !
222 "dpmax" / 1.0e+3) ! for advanced time step control !
223
224 ! numerical scheme for flux calculation
225 fluxes => dict(&
226 "fluxtype" / kt, &
227 "order" / linear, &
228 "variables" / primitive, & ! vars. to use for reconstruction!
229 "limiter" / minmod, & ! one of: minmod, monocent,... !
230 "theta" / 1.2) ! optional parameter for limiter !
231
232 ! source term due to a point mass
233 grav => dict("stype" / gravity, &
234 "output/accel" / 1, &
235 "pmass/gtype" / pointmass, & ! grav. accel. of a point mass !
236 "pmass/potential" / newton, &
237 "pmass/outbound" / 1, &
238! "pmass/acclimit" / (2.2E-9/ETA_ACC/YEAR), &
239 "pmass/mass" / m_bh) ! mass of the accreting object[kg] !
240
241
242 ! account for energy losses due to radiative cooling
243!!$ stemp => Dict(, &
244!!$ "stype" / COOLING,, &
245!!$ "cvis" / 0.9) ! CFL number for cooling time !
246!!$ CALL SetAttr(sources, "cooling", stemp)
247
248
249! ! viscosity source term
250! stemp => Dict("stype" / VISCOSITY, &
251! "cvis" / 0.6, &
252!! "vismodel" / ALPHA, &
253!! "dynconst" / 0.1, &
254! "vismodel" / BETA, &
255! "dynconst" / 0.001, &
256! "output/stress" / 1, &
257! "output/kinvis" / 1, &
258! "output/dynvis" / 1)
259!!$ CALL SetAttr(sources, "viscosity", stemp)
260
261
262 ! time discretization settings
263 timedisc => dict(&
264 "method" / ssprk, &
265 "order" / 5, &
266 "cfl" / 0.4, &
267 "stoptime" / tsim, &
268 "dtlimit" / 1.0e-3, &
269 "maxiter" / 1000000000, &
270 "output/density" / 1, &
271 "output/xvelocity" / 1, &
272 "output/yvelocity" / 1, &
273 "output/zvelocity" / 1, &
274 "output/pressure" / 1, &
275 "output/fluxes" / 1, &
276 "output/rhs" / 1, &
277 "output/sgspressure" / 1, &
278 "output/geometrical_sources" / 1, &
279 "output/external_sources" / 1)
280
281 sources => dict( &
282 "grav" / grav)
283
284! ! initialize log input/output
285! logfile => Dict(, &
286! "fileformat" / BINARY, &
287! +"filename" / (TRIM(ODIR) // TRIM(OFNAME) // "log"), &
288! +"dtwall" / 3600, &
289! +"filecycles" / 1)
290
291 ! initialize data input/output
292 datafile => dict( &
293 "fileformat" / vtk, &
294 ! "fileformat" / XDMF, &
295 "unit" / 5555, &
296 "filename" / (trim(odir) // trim(ofname)), &
297 "stoptime" / tsim, &
298 "count" / onum, &
299 "filecycles" / (onum+1))
300
301 config => dict("mesh" / mesh, &
302 "physics" / physics, &
303 "boundary" / boundary, &
304 "fluxes" / fluxes, &
305 "sources" / sources, &
306 "timedisc" / timedisc, &
307! "logfile" / logfile, &
308 "datafile" / datafile)
309
311
315 IMPLICIT NONE
316 !------------------------------------------------------------------------!
317 CLASS(Mesh_base), INTENT(IN) :: Mesh
318 CLASS(Physics_base), INTENT(INOUT) :: Physics
319 CLASS(Fluxes_base), INTENT(INOUT) :: Fluxes
320 CLASS(Timedisc_base), INTENT(INOUT) :: Timedisc
321 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
322 !------------------------------------------------------------------------!
323 ! Local variable declaration
324 CLASS(sources_base), POINTER :: sp => null()
325 CLASS(sources_gravity), POINTER :: gp => null()
326 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,2) :: dv
327 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX) :: velo_periodic
328 REAL :: s,z,tau,lambda,tau1,tau2,dPs,dPz,s1,s2,z1,z2,P0,rho0,cs_inf
329 REAL :: mdisk,am
330#ifdef PARALLEL
331 INTEGER :: ierror
332 REAL :: mdisk_local,am_local
333#endif
334 INTEGER :: i,j,k,n
335 CHARACTER(LEN=20) :: mdisk_str,am_str
336 !------------------------------------------------------------------------!
337 ! get gravitational acceleration
338 IF (ALLOCATED(sources)) &
339 sp => sources%GetSourcesPointer(gravity)
340 IF (.NOT.ASSOCIATED(sp)) &
341 CALL physics%Error("agndisk::InitData","no gravity term initialized")
342
343 ! initial condition (computational domain)
344 ! the total disk mass is proportional to the central pressure;
345 ! to determine the necessary central pressure for a given
346 ! disk mass, the initial condition is computed twice
347 timedisc%pvar%data2d(:,physics%XVELOCITY) = 0.0
348 timedisc%pvar%data2d(:,physics%ZVELOCITY) = 0.0
349 p0 = 1.0e+10
350 rho0 = 1.0e-09
351 cs_inf = sqrt(gamma*p_inf/rho_inf)
352 DO n=1,2
353 DO k=mesh%KGMIN,mesh%KGMAX
354 DO j=mesh%JGMIN,mesh%JGMAX
355 DO i=mesh%IGMIN,mesh%IGMAX
356 s = abs(mesh%bccart(i,j,k,1))
357 z = mesh%bccart(i,j,k,3)
358!if (s .eq. 0.0 .or. z .eq. 0.0) print *,s,z,i,j
359
360 tau = func_tau(s,z)
361 lambda = func_lambda(z)
362 timedisc%pvar%data4d(i,j,k,physics%PRESSURE) = pc(tau,rho0,cs_inf,gamma)
363 z1 = mesh%cart%faces(i,j,k,5,3) ! western cell faces z coordinate
364 z2 = mesh%cart%faces(i,j,k,6,3) ! eastern cell faces z coordinate
365 tau1 = func_tau(s,z1)
366 tau2 = func_tau(s,z2)
368 IF (abs(z).GT.0.0) THEN
369!!$ Timedisc%pvar(i,j,Physics%DENSITY) = rho0 * rho(s,z,tau,lambda)
370 timedisc%pvar%data4d(i,j,k,physics%DENSITY) = dpz / gz(s,z)
371 ELSE
372 timedisc%pvar%data4d(i,j,k,physics%DENSITY) = rho0 * (s/s0)**b0
373 END IF
374 s1 = abs(mesh%cart%faces(i,j,k,1,1)) ! southern cell faces s coordinate
375 s2 = mesh%cart%faces(i,j,k,2,1) ! northern cell faces s coordinate
376 tau1 = func_tau(s1,z)
377 tau2 = func_tau(s2,z)
378if (tau1 .eq. 0.0 .or. tau2 .eq. 0.0) print *,tau1,tau2,i,j
380! Timedisc%pvar%data4d(i,j,k,Physics%YVELOCITY) = SQRT(MAX(0.0, &
381! s*(dPs/Timedisc%pvar%data4d(i,j,k,Physics%DENSITY) - gs(s,z))))
382 timedisc%pvar%data4d(i,j,k,physics%YVELOCITY) = sqrt(max(0.0,-s*gs(s,z) / lambda))
383 END DO
384 END DO
385 END DO
386 ! 2. initialize gravitational acceleration
387! CALL gp%UpdateGravity(Mesh,Physics,Fluxes,Timedisc%pvar,0.0,0.0)
388! Timedisc%pvar%data4d(:,:,:,2:2+Physics%VDIM) = &
389! Timedisc%GetCentrifugalVelocity(Mesh,Physics,Fluxes,Sources,(/0.,0.,1./),gp%accel%data4d)
390 ! compute disk mass
391 mdisk = sum(timedisc%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%DENSITY) &
392 * mesh%volume%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
393 ! compute angular momentum
394 am = sum(mesh%hy%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) &
395 * timedisc%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%YVELOCITY) &
396 * timedisc%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%DENSITY) &
397 * mesh%volume%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
398
399#ifdef PARALLEL
400 mdisk_local = mdisk
401 CALL mpi_allreduce(mdisk_local,mdisk,1,default_mpi_real,mpi_sum, &
402 mesh%comm_cart,ierror)
403 am_local = am
404 CALL mpi_allreduce(am_local,am,1,default_mpi_real,mpi_sum, &
405 mesh%comm_cart,ierror)
406#endif
407!!$ ! adjust central pressure
408!!$ P0 = M_DISK / mdisk * P0
409 ! adjust density at (S0,Z=0)
410 rho0 = m_disk / mdisk * rho0
411 END DO
412
413 ! add velocity perturbations
414!!$ CALL RANDOM_SEED
415!!$ CALL RANDOM_NUMBER(dv)
416!!$ Timedisc%pvar(:,:,Physics%XVELOCITY) = Timedisc%pvar(:,:,Physics%XVELOCITY) &
417!!$ + (dv(:,:,1)-0.5)*2.0E+2
418!!$ Timedisc%pvar(:,:,Physics%YVELOCITY) = Timedisc%pvar(:,:,Physics%YVELOCITY) &
419!!$ + (dv(:,:,2)-0.5)*2.0E+2
420
421
422! ! eastern
423 SELECT TYPE(beast => timedisc%Boundary%Boundary(east)%p)
424 CLASS IS (boundary_farfield)
425 DO k=mesh%KMIN,mesh%KMAX
426 DO j=mesh%JMIN,mesh%JMAX
427 DO i=1,mesh%GINUM
428 beast%data(i,j,k,:) = timedisc%pvar%data4d(mesh%IMAX+i,j,k,:)
429 END DO
430 END DO
431 END DO
432 END SELECT
433
434 ! southern
435 SELECT TYPE(bwest => timedisc%Boundary%Boundary(west)%p)
436 CLASS IS(boundary_custom)
437 CALL bwest%SetCustomBoundaries(mesh,physics, &
438 (/custom_nograd,custom_outflow,custom_kepler,custom_nograd,custom_nograd/))
439 END SELECT
440
441 ! bottom
442 SELECT TYPE (bbottom => timedisc%Boundary%Boundary(bottom)%p)
443 CLASS IS(boundary_farfield)
444 DO k=1,mesh%GKNUM
445 DO j=mesh%JMIN,mesh%JMAX
446 DO i=mesh%IMIN,mesh%IMAX
447 bbottom%data(i,j,k,:) = timedisc%pvar%data4d(i,j,mesh%KMIN-k,:)
448 END DO
449 END DO
450 END DO
451 END SELECT
452
453 ! Top
454 SELECT TYPE (btop => timedisc%Boundary%Boundary(top)%p)
455 CLASS IS(boundary_farfield)
456 DO k=1,mesh%GKNUM
457 DO j=mesh%JMIN,mesh%JMAX
458 DO i=mesh%IMIN,mesh%IMAX
459 btop%data(i,j,k,:) = timedisc%pvar%data4d(i,j,mesh%KMAX+k,:)
460 END DO
461 END DO
462 END DO
463 END SELECT
464
465 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
466
467 CALL timedisc%Boundary%CenterBoundary(mesh,physics,0.0,timedisc%pvar, &
468 timedisc%cvar)
469
470 ! print some information
471 WRITE (mdisk_str, '(ES10.2)') mdisk/msun
472 WRITE (am_str, '(ES10.2)') am
473 CALL mesh%Info(" DATA-----> initial condition: non Keplerian flow")
474 CALL mesh%Info(" disk mass: " // trim(mdisk_str) // " M_sun")
475 CALL mesh%Info(" angular momentum: " // trim(am_str) // " kg/m^2/s")
476
478
480 IMPLICIT NONE
481 REAL,INTENT(IN) :: sc,zc
482 REAL :: res
483 res = sqrt(gn*m_bh/sc) * (1.0d0+(zc/sc)**2)**(-0.75)
485
487 IMPLICIT NONE
488 REAL,INTENT(IN) :: sc,zc
489 REAL :: res
490 res = -(gn*m_bh/sc**2) / sqrt(1.0d0+(zc/sc)**2)**3
492
494 IMPLICIT NONE
495 REAL,INTENT(IN) :: sc,zc
496 REAL :: res
497 res = -(gn*m_bh/sc**2) * (zc/sc) / sqrt(1.0d0+(zc/sc)**2)**3
499
501 IMPLICIT NONE
502 REAL,INTENT(IN) :: zc
503 REAL :: res
504 res = 1.0d0 + a0/(1.0d0+(zc/z0)**2)
506
508 IMPLICIT NONE
509 REAL,INTENT(IN) :: sc,zc
510 REAL :: res
511 res = sc*sqrt(1.0d0 + (zc/sc)**2*(1.0d0 + (1.0d0+0.5d0*(zc/z0)**2)/a0))
513
515 IMPLICIT NONE
516 REAL,INTENT(IN) :: tau,rho_s0,cs_inf,gamma
517 REAL :: res
518! res = P0 * EXP(-tau/S0)
519 res = p_inf * (1.0 + a0/(1.0+a0) * 1.0/(1.0-b0) * 0.5 * gamma * &
520 (cc/cs_inf)**2 * rho_s0/rho_inf * (rs/s0) * (tau/s0)**(b0-1.0))
522
523!!$ ELEMENTAL FUNCTION dPc(tau) RESULT(res)
524!!$ IMPLICIT NONE
525!!$ REAL,INTENT(IN) :: tau
526!!$ REAL :: res
527!!$ res = -1.0D0/S0 * Pc(tau)
528!!$ END FUNCTION dPc
529
530! ELEMENTAL FUNCTION rho(sc,zc,tau,lambda) RESULT(res)
531! IMPLICIT NONE
532! REAL,INTENT(IN) :: sc,zc,tau,lambda
533! REAL :: res
534! res = A0/(1.0+A0) * lambda/(lambda-1.0) * sqrt(1.0+(zc/sc)**2)**3 &
535! * (tau/S0)**B0 * (s/tau)**3
536! END FUNCTION rho
537
538
Definition: fosite.f90:41
generic source terms module providing functionaly common to all source terms
Definition: sources_base.f90:42
generic gravity terms module providing functionaly common to all gravity terms
Definition: sources_gravity.f90:41
