New Greens function ~/Oload/afor/p/greens-v2.f

The main idea is to sample a Greens function in a long array and interpolate quadratically for further precision when zgreen(distance) is called.
The Greens function resulting from astyanax and greensm (or gres02), often sparsely sampled, is interpolated along a logarithmic scale after removing the singular factor.
The source file includes a function gnrm(distance) for division of the result of zgreen. 

complex*16 function zgreen (dis)
dis        - real*8   - distance in radians

Returns the complex value of the Greens function.


real*8 function physf_zgreen ()
Returns the physical factor, either effect per sea volume (m^3) or
per meter sea level elevation in an equatorial grid cell.

real*8 function v_newton()
Returns the Newtonian effect, a value of 1.0d0 or 0.0d0 depending on the component.

real*8 function gnrm (dis)
Returns the value of the normalisation function. The physical loading effect zeff due to
one meter of sea level elevation is given by
      xnewton=v_newton()
      pf=physf_zgreen()
      zeff=pf*(xnewton+zgreen(dis2load))/grnm(dis2load)*dcos(ylat_load)

A sequence of initialisation and parameter settings must precede first use.

A typical sequence may look like this:
      call zgreen_init (file,comp,npoli,ntbn)
      call zgreen_check ()
      call zgreen_print (what)
      call set_gpars (-1.d0,-1.d0,-1.d0,'mm')
      call use_gridres (.true.,1,1)
      call zgreen_extrapolate (qxtr)
and numerical results are obtained with e.g.
      do i=1,nd-1
         if (dis(i).gt.0.d0) print '(f10.6,2x,1p,2e12.4,2x,e12.4)'
     &,      dis(i), zgreen(dis(i)*rad), 1.d0/gnrm(dis(i)*rad)
      enddo

Parameter setting subroutines

subroutine set_gpars (gsurface,rsurface,rhowater,option)
gsurface  - real*8    - little-g at the earth's surface (9.81 m/s²) ¹
rsurface  - real*8    - earth radius (6.371d6 m) ¹
rhowater  - real*8    - density of sea water  (1023 kg/m³) ¹
                        The constant of gravity is hard-wired (6.6743d-11 m³/kg/s²)
option    - character - to scale the results, any of the following strings
                        ' '  'reset'  'unity'  set the scale to 1
                        'nano'  'nrad'  'nm/s**2'  'nm/s^2'  to 10⁹
                        'micro' 'urad'  'um/s**2'  'um/s^2'  to 10⁶
                        'mm'  'milli'  'cm'  'centi'  to 10³ or 10², respectively.
                        As a final alternative a numerical value in string format
                        can be specified.
                        If the scale is not to be changed, specify 'keep'.
                        'reset' will yield unity scaling.
                        The strings are not case-sensitive
¹A value < 0 defaults to a previously substituted parameter or the internal standard value.


subroutine zgreen_init (greensfn,component,npoli,lentbl)
greensfn  - character - File name of the Greens function to be input.
component - char*8    - At least four characters, any of the following strings
                        'RADI'  'TANG'  'PTAN'  'NN-S'  'XX-S'  'ARES'  'AREA'
                        'GRAV'  'TILT'  'ASTR'  'POTE'  'OTEP'
npoli     - integer   - Order of the interpolation polynomial.
lentbl    - integer   - Length of the interpolated table, <= 100,000


entry zgreen_extrapolate (qq)
qq        - logical   - Allow (if .true.) extrapolation of the tabled function
                        to very small distances. Quadratic extrapolation.

entry use_gridconst (qq,dx,dy,opt)
entry use_gridres (qq,nx,ny)
qq        - logical   - Include (if .true.) the area of an equatorial quadrangle in the
                        physical factor physf_zgreen()
                        Then specify the cell width and height in
dx,dy     - real*8    - in radians (default or opt='rad')
                        in degrees (opt='deg')
                        or in number of subintervals per degree (opt='/deg')
nx,ny     - integer   - like under opt='/deg'
opt       - character - see dx,dy above.


Convenience subroutines:

subroutine zgreen_check ()
finds the value and the distance (table entry) of the maximum successive difference in the
interpolated table. Increasing the polynomial order can give lower such changes but only
to a limited extent. The numerical precision in the input file will set a limit.
Increasing the length of the table might give lower but betraying differences.
 
subroutine zgreen_print (what)
Prints the input values (what='I') or the interpolated table ('T') or both ('B').
what - char*1