This procedure computes (ocean tide) surface loading effects.
 Observation components currently available:

 RADI (radial ~= vertical displacement)
 TANG (tangential = horizontal ~= transverse displacement)
 OTEP (Ocean Tide Effective Potential, converted to vertical displacement of level surface)
 GRAV (gravity)
 TILT (inclination of the vertical with respect to the solid earth surface)
 ASTR (astronomical inclination of the vertical)
 ASTN (=ASTR, different distance normalisation of Green's function, though.)

 Method: Point mass convolution, refinement of ocean grid in the near-zone
 of a site. The loading kernel (Green's) function is of complex type
 in order to allow for visco-elastic earth models.

 In one run, rsload processes a number of sites, one component, and one
 partial tide. TANG, TILT, ASTx (generally: effects along an azimuth) require
 only one run to obtain north and east components.

 Components that imply a Newtonian effect are always computed with that
 effect included (GRAV TILT AST* OTEP)

 GRAVity effects can be computed at topographic height above the sea surface.
 The option is set/removed via a subroutine entry. (The equivalent option
 for tilt is not implemented yet.)

 Rsload can except an ocean area, either a spherical quadrangle or a rect-
 angular tangent plane, stereographically projected
 (keyword: regional area clipping).

 Clipping modes: Either skip the tide inside a spherical quadrangle.
                 Or skip the tide outside a spherical quadrangle.

 If the quadrangle outlines an area for which a plane tide model has been
 computed, weights in the overlapping zones can be used to paste the
 areas together. The weight array can be computed using Oload/p/m/olstm.f

 If a site is near the coast, the near-zone loading effect can be excepted
 for post-processing by OLMPP (keywords: (local) exception processing, post-
 processing). The excepted area is 3x3 or 5x5 load cells wide. The support
 information passed to OLMPP is from a 5x5 or 7x7 area, respectively, always
 with a site at the centre.
   This program does not retrieve the TOPO data though ! (Deferred TOPO
 processing is achieved with the new version of OLMPP.)
 

 Prerequisits:
 Green's function must have been initiated. CALL INIGRE.

 Observing sites in  /CSITE/ - to be set by calling program.
                               Associated BLOCK DATA in ols05.f

 Near zone load exception, announced
                  in /CXCPT/ - array element QXCPT(isite) to be set true by
                               calling program.
                               Associated BLOCK DATA included below.

 Tide model grid information (global, spherical models)
                  in /CRSAR/ - to set, use GETFLA or GETGOT or GETCRT or ...
                               Associated BLOCK DATA in ols04.f

 Area clip limits in /CRSWI/ - to be set by calling program.
                               Associated BLOCK DATA in ols04.f
                               Default: Use all global data.

 Local area overlap weights (plane models)
           in /CGWA/ and  /CUGWA/ - to set, use GETGWA
                               to create, exec OLSTM.exp
                               Default: ignore.

 Refer to olsr01.fc-doc for documentation as to the named common areas.
 
 Calling parameters:
 -------------------

 Z(M,N)        - complex - tide array, organized north (.,1) to south (.,N).
 FLZ(M,N)      - integer - land/sea flags, tested on { .eq.0 | .ne.0 }.
 BSX(M),BSY(M) - real*8  - work arrays.
 CRTDIS        - real    - distance limit [deg] for integrated Green's fct.
 DISMAX        - real    - max.distance [deg] for Green's fct. table.
 QEXPG         - logical - .true.:  Make a new Green's fct. table.
                           .false.: Use the currently available one.
                           returned: .false.
 GDIM          - integer - length of Green's table. 2*(GDIM+1) four-byte
                           words must be allocated in Blank Common.
 CLO(mxs)      - complex - returned: loading effects, units = mm.
 CLX(mxs)      - complex - returned: loading effects in accompanying component
                           (cross-azimuth).
                           mxs  is a hard-wired parameter: number of sites.
 QAZR          - logical - returned: .true. if azimuthal component
 UNIT          - char*4  - returned: metric units
 QTR           - logical - Trace on.

 Blank Common:
 -------------
 To be allocated in Main:
       real*4 ctab,g
       common ctab, g(2*GDIM+2)
 GDIM must be 1000 at least.

 ctab        - real*4 - Look-up index distance scale; set by rsload.
 g           - real*4 - Green's table. Internal use: complex*8
 
     entry rsload_pic (d1,d2,id3,id4,qq)

 
     entry rsload_ortho (qq)
 
     entry rsload_hires (ires)
 
     entry rsload_amp (v)
 
     entry rsload_crit_alt (v)
 
     entry rsload_skip_xcpt
 
     entry rsload_init
 
     entry rsload_warning (qq)