Loading components


Gravity at height

Here is the code fragment used to compute the extra newtonian attraction in the vertical due to topographic height near a tidal coast.
Starting with height above ellipsoid of station IS

            hh = salt(is)/re
            h1 = 1.d0+hh
            h2 = hh*(2-hh/h1)
            hq = hh*hh
            sthq = 4*sthq
            dshq = h1+hq/sthq
            dsh = dsqrt(dshq)
            V_NEWTON = (1/h1+h2/sthq/dshq)/dsh
            gvalu = ( V_NEWTON + G(ITAB) )/sthu

where RE is earth radius, STHQ = STH2 at entry, STH = sin(θ/2), θ the distance angle between the field point on the sphere and the load point, and G() is the Greens function table


Spheroidal potential of tangential displacement

Rather than preparing two data set for horizontal (=tangential) displacement, a spheroidal potential suffices. The gridded PTAN, as it's named, is interpolated together with the gradient operator.

Since we interpolate scalar fields with a quadratic algorithm, the first derivatives are simple to incorporate. The associated function is zquint ( zquint.f ). There are parameter-setting entries. There is a real-valued equivalent in rquint.f with the same parameter-setting calls to quint_opt. Here is the description of zquint.

ZQUINT

      complex function zquint (z,m,n,x,y,zdx,zdy)
Interpolation * 2-D Quadratic polyn., full array.

RQUINT
      real function rquint (r,m,n,x,y,rdx,rdy)

And yet another interpolation routine. We describe the complex routine.
This time the complex array z(m,n) has data everywhere.
Returned value is interpolated quadratically at 1 < x < m, 1 < y < n.
Additionally returned: The gradient at (x,y), (zdx,zdy)

GEO mode:
zdx and zdy are in physical units; length unit=rad.
If this is not desired, call Quint_Opt ('+G-P')

GEO mode activation and parameter-setting calls:
Reading a global tide model should set common area /crsar/.
/crsar/ parameters can alternately be set using

    Call Set_Got_Ar (...)
c.f. oteu54.f If /crsar/ thus describes the geography of Z(M,N),

    Call Quint_Opt ('+G')
signifies GEO mode. Zquint will interpret X,Y as longitude,latitude.
Circular connection is determined by the GLOBAL parameter in common /crsar/.

N.B. in GEO mode Z(.,1) is the northernmost row of data.
N.B. there is no automatic /crsar/ validation; default values fit the Schwiderski grid.

In non-GEO mode,
    Call Quint_Opt ('+C')
will declare Z(M,N) as circularly connected in the first array subscript.

Default option is '-C-G'.

To process Catesian arrays with north row first
    Call Quint_Opt ('-G+R')

To process spherical grids with the south row first:
    Call Quint_Opt ('+G+R')

To enter Trace mode,
    Call Quint_Opt ('+T')

Parameter-setting subroutine     QUINT_OPT

 
      subroutine quint_opt (opt)
Options: The system wakes up with -G -S -P -C -R -T, ER=6.371e6 [m]

+C      - Assume array is circular connected.

+I      - Ignore flags, eventhough zquint_f is called

-G      - Assume plane grid, Cartesian relation between array row and y
          coordinate. That means the first array row is south of the rest.

+G      - Assume circular connection according to "/crsar/ ...,global".
          Assume northernmost row has number 1, and
          use area bounds according to "/crsar/..."
          Imply +P (qphys=.true.), i.e. gradients are per radian.
          To set common /crsar/ use  Set_Got_Ar (...)  or read a tide
          model with a global reading routine.

+S:Z or +S:M  Stereographic grid. Specify grid kind Z or M

+SGZ or +SGM  Stereographic grid, calling with coordinates on the grid
              (floating point, 0<=X<=M, 0<=Y<=N)
              You can code "=" or "g" for "G".

+R      - Reverse order north/south. Obs! Question: What's normal ?
          With -G, +R means that the first row of the array is the north
          of the others. With +G, +R means that the first row of the 
          array is south of the others.

-P      - Gradients are per unit step in the Cartesian grid.
+P      - Physical units of the gradients ON THE EARTH SURFACE.

+Punit  - The following units can be given:
           +P/Mm    - gradient per 1000 km
           +P/km    - gradient per km
           +P/m     - gradient per m
           +P/deg   - gradient per degree

+T      - Switch on trace option = extensive print.

ER=val, - Numeric value of the planet radius in m.

      Entry QUINT_ER (r)

r       - real - New earth radius if r>0
                 If called after setting option +T, quint_er will print the new (or old) radius in the protocol.

Examples:

Spherical grid  from Oload/p/mpp/grdol.f
         implicit complex (z)
         dimension zolm(3), zra(1440,720), zpt(1440,720)
         character gotar*4, opt_quint*8
...
C read two global loading maps with RADI=radial and PTAN components
         call getsiz (31,m,n)
         call getzm  (31,zra,m,n,typ,tag)
         call getzm  (32,zpt,m,n,typ,tag)
         call set_got_ar_for (gotar,n)
...
         call zquint_opt ('+C+G'//opt_quint)
...
C xi,yo = longitude,latitude in degrees
         zolm(1)=zquint(zra,m,n,xi,yi,zolm(2),zolm(3))
         zdummy =zquint(zpt,m,n,xi,yi,zolm(2),zolm(3))




Plane grid  from OTEQ/PROG/s/otes18.f
         dimension el(m,n,2)

C The array EL is regional and on the tangential plane.
C Stereographic projection is initialised.
         call zquint_opt ('+S:Z')
...
C xi,yo = longitude,latitude in degrees
         rz=rquint(el(1,1,ino),m,n,xd,yd,gx,gy)