gres02 - Compute Greens functions from load Love numbers
Main program is located in ~/Green/p/m/gres02.f and
subroutines in ~/Green/p/*.f
Program needs an instruction file (INS, e.g. gres02.ins)
and a file with load Love numbers (LLN).
The LLN's can be given sparsely. Example ~/Green/loadf.dat
Program prints on unit 6. A second output unit can be defined for
more extensive printing.
This is a very early program from the Uppsala (Hällby) period of the
author.
Some special pieces of code are sidetracks and will not be explained
here ("Newton numbers",
more than one LLN file). In essence it follows the maths of W.E.
Farrell (1972).
Not recommended: Use of Euler-Maclaurin for summation.
Structure of the INS-file:
Open-file block
Namelist &PARA
Distances θ if QSPECA .or. QSPECH
Open-file block
21 R LLN input
31 B greens.dat output
#n B protocol
#n B "Newton numbers" if #n
= IUNN .gt. 0
Q
Namelist PARA
Variable - Type -
Explanation
[default]
-----------------------------------------------------------------------------------------------------
LOVE
NUMBERS INPUT
ITAG
- char*8 - Tag for the LLN's, e.g. 'L00EGBC', to
locate the data section
LVNZB,LVNZE
DISTANCES WHERE GREENS IS COMPUTED
QSPECA -
logical - Read all distances after namelist &end
[.false.]
False means that distances are generated autmomatically
QSPECH - logical
- Read distances of the last decade after
&PARA .. &END
[.false.]
NRES
- integer - If .NOT. QSPECA, divide each decade into NRES
intervals
[10]
else read NRES distances from the section after &PARA ..
&END
THMIN,THMAX - real*8 -
Distance range in
degrees
[0.1, 179.9]
"NEWTON NUMBERS" (expansion of a point load that
is not located at
exactly earth's radius)
MTP
IUNN
TPD,TPMIN
NNSTYP
GREENS FUNCTION, output
OTAG -
char*8 - The first character will always be changed to
'G'
[' n.n. ']
The string is used to make the output section searchable
for the Greens function utility ~/Oload/afor/p/greensf.f
NEWTAG -
logical - Use OTAG, else derive from
ITAG
[.false.]
QHEX -
logical - .true.: Write one output line per
distance
[.true.]
.false.: Write two lines if more than 3 components
GREENS FUNCTION, select components:
QREAL
- logcial - compute only the real part even if LLN's are
complex.
[.false.]
IOUTM
- integer - number of components to compute (max
8.)
[4]
IOUT(8) - integer -
component
codes:
[1,2,3,4,5,6,7,0]
1: Potential (1+k)
2: Radial (h)
3: Tangential (l)
4: Gravity (1 + 2h -(n+1)k)
5: Tilt
6: Astronomical deflection
7: Ocean tide-effective potential (1+k-h)
8: Truncated Newton numbers (???)
9: Newton (???)
10: Covariance function itself (???)
11: Areal Strain
12: Strain along-distance
13: Strain cross-distance
14: Like 1, with Newtonian added ("NPOT")
15: Like 4, with Newtonian added ("NGRA")
16: Displacement potential of 3
17: Reserved for special tests
???: deserves further explanation.
QANEWT -
logical - add the Newtonian effect if applicable (???)
[.false.]
The mass is located at Earth's radius
COMPUTATION OF THE INFINITE SUMS: Tilt
TILTD -
real*4 - Below this distance, compute tilt using the
direct formula,
beyond by derivative of tide-effective
potential
[181.0]
COMPUTATION OF THE INFINITE SUMS: Kummer transform
QKUMMER - logical
- Do Kummer transform with the asymptotic
LLN's
[.true.]
COMPUTATION OF THE INFINITE SUMS: van Wijngaarden algorithm
MSEUWI -
integer - Use Euler-Maclaurin if the largest degree is less
than
[11]
MADDEU -
integer - ... MSEUWI+MADDEU
[0]
else van Wijngaarden
NDIRCT
- integer - Direct summation of the first NDIRCT
terms
[100]
EPSR,EPSI -
real*8 - Precision of the algorithm, real and
imaginary
parts
[1.d05, . ]
EPSI will default to 10*EPSR
ILIM -
integer - The maximum number of Euler transformations in
the
[100]
van Wijngaarden algorithm. Actually used is
MIN(ILIM,9995,M+5) where M is the maximum degree of the
LLN's
IPRSST -
integer - How many times the error has to stay below EPSR
and EPSI
[5]
QSUMD -
integer - Use direct summation as the exit strategy if van
Wijngaarden [.false.]
fails.
COMPUTATION OF THE INFINITE SUMS: Optional disk
factor
DSCALL -
logical - If .true., use disk factor on all
components.
[.false.]
If .false., use it only where sums converge badly (like tilt)
QLIMDISC - logical -
.FALSE.: Use fixed disk radius RDISC
else
[.true.]
RDISC
- real*8 - MIN(distance/FRDISC,
RDISC)
[1.0d0]
FRDISC -
real*8
-
[30.0d0]
ASYMPTOTIC LOVE NUMBERS n -> inf.
QINFLN -
logical - request from subroutine RELOVE to return the
asymptotic LLN's [.true.]
TINFH, -K, -L - real*8 - factor to multiply the
asymptotic
LLN's
[3*1.0d0]
This may be useful if the asymptotic values are not in the
file. The max-degree values are substituted, and with these
factors, they can be tuned according to experience.
NORMALISATION:
QFARRL -
logical - only used in printout: use Farrell's normalisation
(against
distance), else normalise with the trigonometric factor
[.false.]
PRINT
AND TRACE:
IUPTC, QTRACE, QTRACP
LTST1,LTST2 - integer - Monitor LLN
interpolation in this degree
interval
[-1,-1]
OTHER:
QAUTOD
(not needed, QSPECA will set it)
IULOVE(2),PTAG,LNPRC
(a second LLN file to modify the first set by percentage of
change)
TVAL
(Legendre function evaluation/iteration test)
TFIX,QCWT,QCWA
(convolution with one of two types of covariances)
.bye para
Example of an instruction file, Farrell's LLN's with disk
factors of fixed diameter
components RADIAL, GRAVITY, OTEPOTNL, NOTP
8 B GREPTC.TXT
31 B
MD00EGBC-0.125.DAT
Trunc, disk, newton cmps
21 R loadf.dat
Q
&PARA
ITAG='L00EGBC ', PTAG='P00I% ', OTAG='
00EGBV ', NEWTAG=.TRUE.
RDISC=0.125, DSCALL=.TRUE., QLIMDISC=.FALSE.
IUPTC=8
THMIN=0.01, THMAX=179.9, NRES=114
EPSR=1.D-6, EPSI=2.D-6, NDIRCT=0, IPRSST=9,
QSUMD=.TRUE.
QTRACE=.FALSE.
ILIM=10000
QREAL=.TRUE.
MSEUWI=13, MADDEU=0
IOUTM=4, IOUT=2,4,7,14,3*0
LVNZB=0, LVNZE=0
LTST1=25,LTST2=40
QTRACP=.FALSE.
QSPECA=.TRUE.
QSPECH=.FALSE.
QFARRL=.FALSE.
QINFLN=.TRUE.
QHEX=.TRUE.
LNPRC=.FALSE.
&END
1.e-2, 1.166e-2, 1.359e-2, 1.585e-2, 1.848e-2, 2.154e-2,
2.512e-2, 2.929e-2,
3.415e-2, 3.981e-2, 4.642e-2, 5.412e-2, 6.31e-2,
7.356e-2, 8.577e-2,
1.e-1, 1.166e-1, 1.359e-1, 1.585e-1, 1.848e-1, 2.154e-1,
2.512e-1, 2.929e-1,
3.415e-1, 3.981e-1, 4.642e-1, 5.412e-1, 6.31e-1,
7.356e-1, 8.577e-1,
1., 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,
3.2, 3.4, 3.6, 3.8,
4., 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0,
6.2, 6.4, 6.6, 6.8,
7., 7.2, 7.4, 7.6, 7.8, 8.0, 8.2, 8.4, 8.6, 8.8, 9.0,
9.2, 9.4, 9.6, 9.8,
10., 10.25, 10.5, 10.75, 11., 11.25, 11.5, 11.75, 12.,
12.25, 12.5, 12.75,
13., 13.25, 13.5,
13.75, 14., 14.25, 14.5, 14.75, 15., 15.25, 15.5, 15.75,
16., 16.25, 16.5,
16.75, 17., 17.25,
18., 18.25, 18.5, 18.75, 19., 19.25, 19.5, 19.75, 20.
.bye