path: root/src/slalib/refro.f

SUBROUTINE sla_REFRO (ZOBS, HM, TDK, PMB, RH, WL, PHI, TLR,
     :                      EPS, REF)
*+
*     - - - - - -
*      R E F R O
*     - - - - - -
*
*  Atmospheric refraction for radio and optical/IR wavelengths.
*
*  Given:
*    ZOBS    d  observed zenith distance of the source (radian)
*    HM      d  height of the observer above sea level (metre)
*    TDK     d  ambient temperature at the observer (deg K)
*    PMB     d  pressure at the observer (millibar)
*    RH      d  relative humidity at the observer (range 0-1)
*    WL      d  effective wavelength of the source (micrometre)
*    PHI     d  latitude of the observer (radian, astronomical)
*    TLR     d  temperature lapse rate in the troposphere (degK/metre)
*    EPS     d  precision required to terminate iteration (radian)
*
*  Returned:
*    REF     d  refraction: in vacuo ZD minus observed ZD (radian)
*
*  Notes:
*
*  1  A suggested value for the TLR argument is 0.0065D0.  The
*     refraction is significantly affected by TLR, and if studies
*     of the local atmosphere have been carried out a better TLR
*     value may be available.
*
*  2  A suggested value for the EPS argument is 1D-8.  The result is
*     usually at least two orders of magnitude more computationally
*     precise than the supplied EPS value.
*
*  3  The routine computes the refraction for zenith distances up
*     to and a little beyond 90 deg using the method of Hohenkerk
*     and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
*     in the Explanatory Supplement, 1992 edition - see section 3.281).
*
*  4  The code is a development of the optical/IR refraction subroutine
*     AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to
*     support the radio case.  Apart from merely cosmetic changes, the
*     following modifications to the original HMNAO optical/IR refraction
*     code have been made:
*
*     .  The angle arguments have been changed to radians.
*
*     .  Any value of ZOBS is allowed (see note 6, below).
*
*     .  Other argument values have been limited to safe values.
*
*     .  Murray's values for the gas constants have been used
*        (Vectorial Astrometry, Adam Hilger, 1983).
*
*     .  The numerical integration phase has been rearranged for
*        extra clarity.
*
*     .  A better model for Ps(T) has been adopted (taken from
*        Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
*
*     .  More accurate expressions for Pwo have been adopted
*        (again from Gill 1982).
*
*     .  Provision for radio wavelengths has been added using
*        expressions devised by A.T.Sinclair, RGO (private
*        communication 1989), based on the Essen & Froome
*        refractivity formula adopted in Resolution 1 of the
*        13th International Geodesy Association General Assembly
*        (Bulletin Geodesique 70 p390, 1963).
*
*     .  Various small changes have been made to gain speed.
*
*     None of the changes significantly affects the optical/IR results
*     with respect to the algorithm given in the 1992 Explanatory
*     Supplement.  For example, at 70 deg zenith distance the present
*     routine agrees with the ES algorithm to better than 0.05 arcsec
*     for any reasonable combination of parameters.  However, the
*     improved water-vapour expressions do make a significant difference
*     in the radio band, at 70 deg zenith distance reaching almost
*     4 arcsec for a hot, humid, low-altitude site during a period of
*     low pressure.
*
*  5  The radio refraction is chosen by specifying WL > 100 micrometres.
*     Because the algorithm takes no account of the ionosphere, the
*     accuracy deteriorates at low frequencies, below about 30 MHz.
*
*  6  Before use, the value of ZOBS is expressed in the range +/- pi.
*     If this ranged ZOBS is -ve, the result REF is computed from its
*     absolute value before being made -ve to match.  In addition, if
*     it has an absolute value greater than 93 deg, a fixed REF value
*     equal to the result for ZOBS = 93 deg is returned, appropriately
*     signed.
*
*  7  As in the original Hohenkerk and Sinclair algorithm, fixed values
*     of the water vapour polytrope exponent, the height of the
*     tropopause, and the height at which refraction is negligible are
*     used.
*
*  8  The radio refraction has been tested against work done by
*     Iain Coulson, JACH, (private communication 1995) for the
*     James Clerk Maxwell Telescope, Mauna Kea.  For typical conditions,
*     agreement at the 0.1 arcsec level is achieved for moderate ZD,
*     worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg.  At hot and
*     humid sea-level sites the accuracy will not be as good.
*
*  9  It should be noted that the relative humidity RH is formally
*     defined in terms of "mixing ratio" rather than pressures or
*     densities as is often stated.  It is the mass of water per unit
*     mass of dry air divided by that for saturated air at the same
*     temperature and pressure (see Gill 1982).
*
*  Called:  sla_DRANGE, sla__ATMT, sla__ATMS
*
*  P.T.Wallace   Starlink   3 June 1997
*
*  Copyright (C) 1997 Rutherford Appleton Laboratory
*-

      IMPLICIT NONE

      DOUBLE PRECISION ZOBS,HM,TDK,PMB,RH,WL,PHI,TLR,EPS,REF

*
*  Fixed parameters
*
      DOUBLE PRECISION D93,GCR,DMD,DMW,S,DELTA,HT,HS
*  93 degrees in radians
      PARAMETER (D93=1.623156204D0)
*  Universal gas constant
      PARAMETER (GCR=8314.32D0)
*  Molecular weight of dry air
      PARAMETER (DMD=28.9644D0)
*  Molecular weight of water vapour
      PARAMETER (DMW=18.0152D0)
*  Mean Earth radius (metre)
      PARAMETER (S=6378120D0)
*  Exponent of temperature dependence of water vapour pressure
      PARAMETER (DELTA=18.36D0)
*  Height of tropopause (metre)
      PARAMETER (HT=11000D0)
*  Upper limit for refractive effects (metre)
      PARAMETER (HS=80000D0)

      INTEGER IS,K,N,I,J
      LOGICAL OPTIC,LOOP
      DOUBLE PRECISION ZOBS1,ZOBS2,HMOK,TDKOK,PMBOK,RHOK,WLOK,ALPHA,
     :                 TOL,WLSQ,GB,A,GAMAL,GAMMA,GAMM2,DELM2,
     :                 TDC,PSAT,PWO,W,
     :                 C1,C2,C3,C4,C5,C6,R0,TEMPO,DN0,RDNDR0,SK0,F0,
     :                 RT,TT,DNT,RDNDRT,SINE,ZT,FT,DNTS,RDNDRP,ZTS,FTS,
     :                 RS,DNS,RDNDRS,ZS,FS,REFOLD,Z0,ZRANGE,FB,FF,FO,FE,
     :                 H,R,SZ,RG,DR,TG,DN,RDNDR,T,F,REFP,REFT

      DOUBLE PRECISION sla_DRANGE

*  The refraction integrand
      DOUBLE PRECISION REFI
      REFI(R,DN,RDNDR) = RDNDR/(DN+RDNDR)



*  Transform ZOBS into the normal range.
      ZOBS1 = sla_DRANGE(ZOBS)
      ZOBS2 = MIN(ABS(ZOBS1),D93)

*  Keep other arguments within safe bounds.
      HMOK = MIN(MAX(HM,-1D3),10D3)
      TDKOK = MIN(MAX(TDK,100D0),500D0)
      PMBOK = MIN(MAX(PMB,0D0),10000D0)
      RHOK = MIN(MAX(RH,0D0),1D0)
      WLOK = MAX(WL,0.1D0)
      ALPHA = MIN(MAX(ABS(TLR),0.001D0),0.01D0)

*  Tolerance for iteration.
      TOL = MIN(MAX(ABS(EPS),1D-12),0.1D0)/2D0

*  Decide whether optical/IR or radio case - switch at 100 microns.
      OPTIC = WLOK.LE.100D0

*  Set up model atmosphere parameters defined at the observer.
      WLSQ = WLOK*WLOK
      GB = 9.784D0*(1D0-0.0026D0*COS(PHI+PHI)-0.00000028D0*HMOK)
      IF (OPTIC) THEN
         A = (287.604D0+(1.6288D0+0.0136D0/WLSQ)/WLSQ)
     :                                              *273.15D-6/1013.25D0
      ELSE
         A = 77.624D-6
      END IF
      GAMAL = (GB*DMD)/GCR
      GAMMA = GAMAL/ALPHA
      GAMM2 = GAMMA-2D0
      DELM2 = DELTA-2D0
      TDC = TDKOK-273.15D0
      PSAT = 10D0**((0.7859D0+0.03477D0*TDC)/(1D0+0.00412D0*TDC))*
     :                                (1D0+PMBOK*(4.5D-6+6D-10*TDC*TDC))
      IF (PMBOK.GT.0D0) THEN
         PWO = RHOK*PSAT/(1D0-(1D0-RHOK)*PSAT/PMBOK)
      ELSE
         PWO = 0D0
      END IF
      W = PWO*(1D0-DMW/DMD)*GAMMA/(DELTA-GAMMA)
      C1 = A*(PMBOK+W)/TDKOK
      IF (OPTIC) THEN
         C2 = (A*W+11.2684D-6*PWO)/TDKOK
      ELSE
         C2 = (A*W+12.92D-6*PWO)/TDKOK
      END IF
      C3 = (GAMMA-1D0)*ALPHA*C1/TDKOK
      C4 = (DELTA-1D0)*ALPHA*C2/TDKOK
      IF (OPTIC) THEN
         C5 = 0D0
         C6 = 0D0
      ELSE
         C5 = 371897D-6*PWO/TDKOK
         C6 = C5*DELM2*ALPHA/(TDKOK*TDKOK)
      END IF

*  Conditions at the observer.
      R0 = S+HMOK
      CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,
     :                                              R0,TEMPO,DN0,RDNDR0)
      SK0 = DN0*R0*SIN(ZOBS2)
      F0 = REFI(R0,DN0,RDNDR0)

*  Conditions in the troposphere at the tropopause.
      RT = S+HT
      CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,
     :                                                 RT,TT,DNT,RDNDRT)
      SINE = SK0/(RT*DNT)
      ZT = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0)))
      FT = REFI(RT,DNT,RDNDRT)

*  Conditions in the stratosphere at the tropopause.
      CALL sla__ATMS(RT,TT,DNT,GAMAL,RT,DNTS,RDNDRP)
      SINE = SK0/(RT*DNTS)
      ZTS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0)))
      FTS = REFI(RT,DNTS,RDNDRP)

*  Conditions at the stratosphere limit.
      RS = S+HS
      CALL sla__ATMS(RT,TT,DNT,GAMAL,RS,DNS,RDNDRS)
      SINE = SK0/(RS*DNS)
      ZS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0)))
      FS = REFI(RS,DNS,RDNDRS)

*
*  Integrate the refraction integral in two parts;  first in the
*  troposphere (K=1), then in the stratosphere (K=2).
*

*  Initialize previous refraction to ensure at least two iterations.
      REFOLD = 1D6

*  Start off with 8 strips for the troposphere integration, and then
*  use the final troposphere value for the stratosphere integration,
*  which tends to need more strips.
      IS = 8

*  Troposphere then stratosphere.
      DO K = 1,2

*     Start Z, Z range, and start and end values.
         IF (K.EQ.1) THEN
            Z0 = ZOBS2
            ZRANGE = ZT-Z0
            FB = F0
            FF = FT
         ELSE
            Z0 = ZTS
            ZRANGE = ZS-Z0
            FB = FTS
            FF = FS
         END IF

*     Sums of odd and even values.
         FO = 0D0
         FE = 0D0

*     First time through the loop we have to do every point.
         N = 1

*     Start of iteration loop (terminates at specified precision).
         LOOP = .TRUE.
         DO WHILE (LOOP)

*        Strip width.
            H = ZRANGE/DBLE(IS)

*        Initialize distance from Earth centre for quadrature pass.
            IF (K.EQ.1) THEN
               R = R0
            ELSE
               R = RT
            END IF

*        One pass (no need to compute evens after first time).
            DO I=1,IS-1,N

*           Sine of observed zenith distance.
               SZ = SIN(Z0+H*DBLE(I))

*           Find R (to the nearest metre, maximum four iterations).
               IF (SZ.GT.1D-20) THEN
                  W = SK0/SZ
                  RG = R
                  DR = 1D6
                  J = 0
                  DO WHILE (ABS(DR).GT.1D0.AND.J.LT.4)
                     J=J+1
                     IF (K.EQ.1) THEN
                        CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,
     :                                 C1,C2,C3,C4,C5,C6,RG,TG,DN,RDNDR)
                     ELSE
                        CALL sla__ATMS(RT,TT,DNT,GAMAL,RG,DN,RDNDR)
                     END IF
                     DR = (RG*DN-W)/(DN+RDNDR)
                     RG = RG-DR
                  END DO
                  R = RG
               END IF

*           Find the refractive index and integrand at R.
               IF (K.EQ.1) THEN
                  CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,
     :                                   C1,C2,C3,C4,C5,C6,R,T,DN,RDNDR)
               ELSE
                  CALL sla__ATMS(RT,TT,DNT,GAMAL,R,DN,RDNDR)
               END IF
               F = REFI(R,DN,RDNDR)

*           Accumulate odd and (first time only) even values.
               IF (N.EQ.1.AND.MOD(I,2).EQ.0) THEN
                  FE = FE+F
               ELSE
                  FO = FO+F
               END IF
            END DO

*        Evaluate the integrand using Simpson's Rule.
            REFP = H*(FB+4D0*FO+2D0*FE+FF)/3D0

*        Has the required precision been achieved?
            IF (ABS(REFP-REFOLD).GT.TOL) THEN

*           No: prepare for next iteration.

*           Save current value for convergence test.
               REFOLD = REFP

*           Double the number of strips.
               IS = IS+IS

*           Sum of all current values = sum of next pass's even values.
               FE = FE+FO

*           Prepare for new odd values.
               FO = 0D0

*           Skip even values next time.
               N = 2
            ELSE

*           Yes: save troposphere component and terminate the loop.
               IF (K.EQ.1) REFT = REFP
               LOOP = .FALSE.
            END IF
         END DO
      END DO

*  Result.
      REF = REFT+REFP
      IF (ZOBS1.LT.0D0) REF = -REF

      END