c Last Update to Flux Routines: 28 Apr 1994
c
      SUBROUTINE albedo(iy,m,id,ihr,min,poslat,poslong,solar_in,
     *                  albedo_payne,albedo_fresnel,ierr)
                        
c Peter A. Coppin
c CSIRO Centre for Environmental Mechanics
c email: coppin@enmech.csiro.au
c version 1.0, 28-Apr-94

c this subroutine calculates the albedo of sea water using two methods;
c
c the first is via the model of PAYNE (1972)
c       Payne, R.E. (1972) "Albedo of the sea Surface"
c       J. Atm. Sci. 29: 959-970
c which is essentially a fit to data obtained over open water.
c It uses only date, time, position, the sun's elevation and
c the total short wave solar radiation as inputs. The latter is compared 
c to the maximum solar input at that sun angle with no atmosphere. This is
c a simple parameterisation for cloudiness. The values from the model fit
c well with observations in a statistical sense, but due to the simplicity
c of the model individual values may be incorrect.
c
c the second method simply calculates the theoretical albedo for sea water
c with properties encountered in TOGA-COARE using Fresnel theory
c i.e. point radiation source, perfectly flat surface and
c perfectly clear water. This provides an interesting theoretical
c limit to albedo values. Real values will be higher in the middle of the day
c and lower at the extremes of the day due to the presence of diffuse 
c radiation and waves on the sea surface.

c NOTE:
c the time used is important. It should be the centre time of any averaging
c period. The routine gives an error if the sun is below the horizon, so
c adjustments to the nominated time may be necessary. 
      
c INPUTS
c iy - (integer) year between 50 and 99 (1950 and 1999)
c m - (integer) month between 1 and 12
c id - (integer) day between 1 and end of month (<31!)
c ihr - (integer) hour between 1 and 24 (this must be UTC!!)
c min - (integer) minute between 1 and 60 (also UTC)
c poslat - (real) latitude in decimal degrees, S is negative
c poslong - (real) longitude in decimal degrees, E is positive
c solar_in - (real) incoming solar short wave radiation in W/m^2
c
c OUTPUTS
c albedo_payne - (real) albedo from PAYNE model ( value between 0.0 and 1.0)
c             sometimes returned as nonsense value when error
c albedo_fresnel - (real) albedo from Fresnel theory ( value between 0.0 and 1.0)
c             sometimes returned as nonsense value when error
c ierr - error flag, see following comments
c  
c users MUST check that the error code is clear (ierr=0) as there is a multitude of
c ways for the routine to fail
c
c ERROR CODES
c ierr=1 (from SUNVEC) year not between 50 and 99
c ierr=2 (from SUNVEC) month > 12 or day >31
c ierr=3 (from SUNVEC) month or day negative
c ierr=4  solar zenith distance above 90 deg i.e. NIGHT TIME!!
c ierr=5  atmospheric transmission negative of exceeds 1, check solar radiation
c          value or time of day used is not appropriate.                      
c ierr=6  error in PAYNE model lookup table (routine should not get this far)

c there is some double precision arithmetic in the SUNVEC subroutine which
c               should be OK on most FORTRAN77 compilers
c
      pi = 4. * ATAN(1.)
      drad = pi / 180.
      ierr=0
      albedo_payne=999.99
      albedo_fresnel=999.99
      
c call solar angle calculator
         
      call SUNVEC(poslat,poslong,iy,m,id,ihr,min,SL,SM,SN,ierrsv)
          
      IF (ierrsv.ne.0) THEN
        ierr=ierrsv
        return
      END IF
        
        azimuth = ATAN2(-sm , sl) / drad
       
      IF( azimuth .lt. 0) THEN
        azimuth = 360 + azimuth
      END IF
      
      zenangle = ATAN2(SQRT(1 - sn ** 2) , sn) / drad
        
      if (zenangle .gt. 90.0)then
        ierr=4
        return
      else 
        elevation=90.0-zenangle
c calculate fresnel limit albedo

        call fresnel(zenangle,albedo_fresnel)
        
c  calculate atmospheric transmissivity

        trans=solar_in/(1353.0*cos(zenangle*drad))

        if(trans.ge.0.0.and.trans.le.1.0)then      
c call PAYNE model
          call PAYNE(trans,elevation,albedo_payne,ierrp)      
          if(ierrp.ne.0)then
            ierr=6
            return
          end if
        else
          ierr=5
          return
        end if
      
      end if
      
      RETURN
      END



c-----------------------------------------------------------------------
      SUBROUTINE PAYNE(trans,elevation,albedo,ierr)
      
      DIMENSION payneobs(21,46)
      
      
      DATA (payneobs(i,1),i=1,21)/0.061,0.062,0.072,0.087,0.115,0.163,
     *0.235,0.318,0.395,0.472,0.542,0.604,0.655,0.693,0.719,0.732,
     *0.73, 0.681,0.581,0.453,0.425/
      DATA (payneobs(i,2),i=1,21)/0.061,0.062,0.07, 0.083,0.108,0.145,
     *0.198,0.263,0.336,0.415,0.487,0.547,0.595,0.631,0.656,0.67,
     *0.652,0.602,0.494,0.398,0.37/
      DATA (payneobs(i,3),i=1,21)/0.061,0.061,0.068,0.079,0.098,0.13,
     *0.174,0.228,0.29, 0.357,0.424,0.498,0.556,0.588,0.603,0.592,
     *0.556,0.488,0.393,0.342,0.325/
      DATA (payneobs(i,4),i=1,21)/0.061,0.061,0.065,0.073,0.086,0.11,
     * 0.15, 0.192,0.248,0.306,0.36, 0.407,0.444,0.469,0.48, 0.474,
     *0.444,0.386,0.333,0.301,0.29/
      DATA (payneobs(i,5),i=1,21)/0.061,0.061,0.065,0.07, 0.082,0.101,
     *0.131,0.168,0.208,0.252,0.295,0.331,0.358,0.375,0.385,0.377,
     *0.356,0.32, 0.288,0.266,0.255/
      DATA (payneobs(i,6),i=1,21)/0.061,0.061,0.063,0.068,0.077,0.092,
     *0.114,0.143,0.176,0.21, 0.242,0.272,0.288,0.296,0.3,  0.291,
     *0.273,0.252,0.237,0.266,0.22/
      DATA (payneobs(i,7),i=1,21)/0.061,0.061,0.062,0.066,0.072,0.084,
     *0.103,0.127,0.151,0.176,0.198,0.219,0.236,0.245,0.25, 0.246,
     *0.235,0.222,0.211,0.205,0.2/
      DATA (payneobs(i,8),i=1,21)/0.061,0.061,0.061,0.065,0.071,0.079,
     *0.094,0.113,0.134,0.154,0.173,0.185,0.19, 0.193,0.193,0.19,
     *0.188,0.185,0.182,0.18, 0.178/
      DATA (payneobs(i,9),i=1,21)/0.061,0.061,0.061,0.064,0.067,0.072,
     *0.083,0.099,0.117,0.135,0.15, 0.16, 0.164,0.165,0.164,0.162,
     *0.16, 0.159,0.158,0.157,0.157/
      DATA (payneobs(i,10),i=1,21)/0.061,0.061,0.061,0.063,0.067,0.072,
     *0.08, 0.092,0.107,0.125,0.136,0.141,0.145,0.145,0.145,0.144,
     *0.143,0.142,0.141,0.14, 0.14/
      DATA (payneobs(i,11),i=1,21)/0.061,0.061,0.061,0.062,0.065,0.068,
     *0.074,0.084,0.097,0.111,0.121,0.127,0.13, 0.131,0.131,0.13,
     *0.129,0.127,0.126,0.125,0.122/
      DATA (payneobs(i,12),i=1,21)/0.061,0.061,0.061,0.061,0.063,0.067,
     *0.074,0.082,0.091,0.102,0.11, 0.116,0.119,0.118,0.116,0.114,
     *0.113,0.111,0.11, 0.109,0.108/
      DATA (payneobs(i,13),i=1,21)/0.061,0.061,0.061,0.061,0.062,0.064,
     *0.07, 0.076,0.085,0.094,0.101,0.105,0.107,0.106,0.103,0.1,  
     *0.097,0.096,0.095,0.095,0.095/
      DATA (payneobs(i,14),i=1,21)/0.061,0.061,0.061,0.06, 0.061,0.063,
     *0.067,0.072,0.079,0.086,0.093,0.097,0.098,0.097,0.092,0.088,
     *0.086,0.084,0.083,0.083,0.083/
      DATA (payneobs(i,15),i=1,21)/0.061,0.061,0.061,0.06, 0.061,0.062,
     *0.065,0.07, 0.075,0.081,0.086,0.089,0.09, 0.088,0.084,0.08, 
     *0.077,0.075,0.074,0.074,0.074/
      DATA (payneobs(i,16),i=1,21)/0.061,0.061,0.061,0.06, 0.06, 0.061,
     *0.064,0.067,0.071,0.076,0.081,0.083,0.084,0.081,0.076,0.072,
     *0.069,0.067,0.066,0.065,0.065/
      DATA (payneobs(i,17),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.061,
     *0.063,0.065,0.068,0.072,0.076,0.077,0.076,0.074,0.071,0.067,
     *0.064,0.062,0.061,0.061,0.061/
      DATA (payneobs(i,18),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.061,
     *0.062,0.064,0.067,0.071,0.073,0.074,0.073,0.069,0.065,0.062,
     *0.06, 0.058,0.057,0.057,0.056/
      DATA (payneobs(i,19),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, 
     *0.061,0.062,0.065,0.068,0.069,0.069,0.068,0.065,0.061,0.058,
     *0.055,0.054,0.053,0.052,0.052/
      DATA (payneobs(i,20),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, 
     *0.06, 0.062,0.063,0.066,0.067,0.066,0.064,0.061,0.057,0.054,
     *0.051,0.05, 0.049,0.048,0.048/
      DATA (payneobs(i,21),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, 
     *0.06, 0.06, 0.062,0.065,0.065,0.063,0.06, 0.057,0.054,0.05, 
     *0.047,0.046,0.045,0.044,0.044/
      DATA (payneobs(i,22),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, 
     *0.06, 0.06, 0.061,0.063,0.064,0.061,0.058,0.055,0.051,0.047,
     *0.044,0.042,0.041,0.04, 0.04/
      DATA (payneobs(i,23),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, 
     *0.059,0.06, 0.06, 0.062,0.062,0.059,0.056,0.052,0.049,0.045,
     *0.042,0.04, 0.039,0.038,0.038/
      DATA (payneobs(i,24),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.059,0.06, 0.061,0.06, 0.057,0.054,0.05, 0.047,0.043,
     *0.039,0.036,0.034,0.033,0.033/
      DATA (payneobs(i,25),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059,
     *0.059,0.059,0.06, 0.06, 0.059,0.056,0.053,0.049,0.045,0.041,
     *0.037,0.035,0.033,0.032,0.032/
      DATA (payneobs(i,26),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.059,0.059,0.059,0.058,0.055,0.051,0.047,0.043,0.039,
     *0.035,0.033,0.032,0.031,0.031/
      DATA (payneobs(i,27),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.059,0.059,0.058,0.057,0.054,0.05, 0.046,0.043,0.039,
     *0.035,0.032,0.031,0.03, 0.03/
      DATA (payneobs(i,28),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.059,0.058,0.057,0.056,0.053,0.049,0.046,0.042,0.038,
     *0.035,0.032,0.03, 0.029,0.029/
      DATA (payneobs(i,29),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.059,0.058,0.057,0.055,0.053,0.048,0.044,0.041,0.037,
     *0.034,0.031,0.029,0.028,0.028/
      DATA (payneobs(i,30),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.059,0.058,0.057,0.055,0.052,0.048,0.044,0.04, 0.036,
     *0.033,0.03, 0.028,0.027,0.027/
      DATA (payneobs(i,31),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.056,0.054,0.051,0.047,0.043,0.039,0.036,
     *0.033,0.03, 0.028,0.027,0.026/
      DATA (payneobs(i,32),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.055,0.053,0.05, 0.046,0.042,0.039,0.035,
     *0.032,0.03, 0.028,0.026,0.026/
      DATA (payneobs(i,33),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.055,0.053,0.05, 0.046,0.042,0.038,0.035,
     *0.032,0.03, 0.028,0.026,0.026/
      DATA (payneobs(i,34),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.055,0.052,0.049,0.045,0.041,0.038,0.034,
     *0.032,0.029,0.027,0.026,0.026/
      DATA (payneobs(i,35),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.055,0.052,0.049,0.045,0.041,0.037,0.033,
     *0.029,0.027,0.026,0.025,0.025/
      DATA (payneobs(i,36),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.055,0.052,0.049,0.045,0.04, 0.036,0.032,
     *0.029,0.027,0.026,0.025,0.025/
      DATA (payneobs(i,37),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.057,0.055,0.051,0.048,0.044,0.04, 0.036,0.032,
     *0.029,0.027,0.026,0.025,0.025/
      DATA (payneobs(i,38),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.056,0.054,0.051,0.047,0.043,0.039,0.035,0.032,
     *0.029,0.027,0.025,0.025,0.025/
      DATA (payneobs(i,39),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.057,0.056,0.053,0.05, 0.047,0.043,0.039,0.035,0.031,
     *0.028,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,40),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.056,0.054,0.05, 0.047,0.043,0.039,0.034,0.031,
     *0.028,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,41),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.056,0.053,0.05, 0.046,0.042,0.038,0.034,0.031,
     *0.028,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,42),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.056,0.053,0.05, 0.046,0.042,0.038,0.034,0.03, 
     *0.027,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,43),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059,
     *0.059,0.058,0.056,0.054,0.051,0.047,0.043,0.038,0.034,0.03, 
     *0.027,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,44),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059,
     *0.058,0.057,0.056,0.054,0.05, 0.047,0.042,0.038,0.034,0.03, 
     *0.027,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,45),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059,
     *0.058,0.057,0.056,0.053,0.05, 0.046,0.042,0.038,0.034,0.03, 
     *0.028,0.026,0.025,0.025,0.025/
      DATA (payneobs(i,46),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059,
     *0.058,0.057,0.055,0.053,0.05, 0.046,0.042,0.038,0.034,0.03, 
     *0.028,0.026,0.025,0.025,0.025/

      ierr=0
      
      if(elevation.ge.0.0.and.elevation.le.90.0)then
        ielev=nint(elevation/2.)+1
      else
        ierr=1
      end if
      
      if(trans.ge.0.0.and.trans.le.1.0)then
        itrans=nint(trans*20.)+1
      else
        ierr=1
      end if

      if(ierr.eq.0)then
        albedo=payneobs(itrans,ielev)
      else
        albedo=999.99
      end if
      
      RETURN
      END
      
c---------------------------------------------------------------------
      subroutine fresnel(zenangle,albedo_fresnel)
      
c this routine calculates the Fresnel limit albedo for sea water
c as a funstion of solar zenith distance in degrees
c for water at 35 ppt salinity and 30 deg C i.e.ref index of 1.339675

      pi = 4. * ATAN(1.)
      drad = pi / 180.
      refindex=1.339675
      alpha=zenangle*drad
      beta=asin(sin(alpha)/refindex)
      albedo_fresnel=0.5*((tan(alpha-beta)**2/tan(alpha+beta)**2)+
     * (sin(alpha-beta)**2/sin(alpha+beta)**2))
     
      return
      end


c---------------------------------------------------------------------------------        

        SUBROUTINE SUNVEC(poslat,poslong,iy,m,id,ihr,min,SL,SM,SN,ierr)

C       SHELL 14.8.73   FOLLOWS GOODSPEED.
C       SUBROUTINE VERSION FOR CALCULATING VECTORS AS REQUIRED.
C       CALCULATES SUN'S VECTOR DIRECTION COSINES(SL,SM,SN)
C  TAKES IN TIME,DAY,MON,YEAR.  CONV. TIME TO HR AND MIN
C
C       RESULTS HAVE BEEN CHECKED AGAINST GOODSPEED'S PROGRAM
C       ON THE 3600, AND THE ANGLES GIVEN BY THIS PROGRAM
C       ON THE PDP-9 ARE CORRECT TO WITHIN 0.01 DEGREES.
C  MOD  DE  14.8.73
C
C
        DOUBLE PRECISION DAY,DAYT,DELTDY,TWOPI,T0,DAY2,A,ARG
        REAL LONG,LAT,NORTH
C
        DIMENSION B(3),OM(3),EAST(3),NORTH(3),GL(3),GN(3),OVERH(3)
        DIMENSION SINTOM(3),RADVEC(3),MON(12)
        DATA GL(1),GL(2),GL(3)/-0.226715,-0.973961,0.0/
        DATA GN(1),GN(2),GN(3)/0.89356453,-0.2080006,-0.3978419/
        DATA OM(1),OM(2),OM(3)/0.3874825,-0.090196746,0.917454/
        TWOPI=.6283185308D01
        T0=60.0 * 24.0 * 1.90123508D - 06
        PI=3.141593

        ierr=0
        LONG=poslong*TWOPI/360.0
        LAT=poslat*TWOPI/360.0
        SINLAT=SIN(LAT)
        COSLAT=COS(LAT)
C
C  OLD SUBR GSDATE
C       SHELL 21.6.73  CALCULATES NUMBER OF DAYS FROM 4.1.50 FOR GSSUN 
C
        DO 3 I=1,12
3       MON(I)=31
        HR=IHR
        RMIN=MIN
C
        MON(4)=30
        MON(6)=30
        MON(9)=30
        MON(11)=30
        MON(2)=28

        IF(IY.LT.50.OR.IY.GT.100) then
          ierr=1
          return
        END IF
        
        Y=IY
        
        IF(ID.GT.31.OR.M.GT.12) then
          ierr=2
          return
        END IF
        
        IF(ID.LE.0.OR.M.LE.0) then
           ierr=3
           return
        END IF
         
C       NOW CALCULATE NUMBER OF DAYS SINCE 1.1.IY
C
        IDAY=0
        DO 10 K=1,12
        IF(K.EQ.M) GO TO 20
10      IDAY=IDAY+MON(K)
20      IDAY=IDAY+ID
C       NOW ADD NUMBER OF DAYS FROM 4.1.50 TO 1.1.IY
        IIY=IY-50
        IDAY=IDAY+365*IIY
C       ADD 1 DAY FOR EACH LEAP YEAR BETWEEN 1950 AND 19IY
        IF(IY.LT.52) GO TO 31
        DO 35 I=52,IY,4
        IF(I.EQ.IY.AND.M.LT.3) GO TO 31
        IDAY=IDAY+1
35      CONTINUE
C
31      IDAY=IDAY-4
        DAY=IDAY
        HOUR=HR+RMIN/60.0
c change to zero hour zone offset to suit UTC
        DELTDY=(HOUR-(1.0+40.0/60.0))/24.0
c        DELTDY=(HOUR-(1.0+10.0+40.0/60.0))/24.0
C       ALLOWS FOR PERIHELION AT 0140HRS, 1950, AND HOURZONE.
        DAY=DAY+DELTDY
        DAYT=DAY
C
C  OLD SUBR GSNUSS
C       SHELL  6.7.73    FOLLOWS GOODSPEED
C       CALCULATES VECTOR FROM SUN TO EARTH AT GIVEN TIME.
C       EC IS ECCENTRICITY OF EARTH'S ORBIT AROUND SUN.
C
C
        EC=0.016751
        DAY2=DAY*T0
        DAY2=DMOD(DAY2,1.0D0)
C       THIS TAKES THE REMAINDER ONLY !!!
        A=TWOPI*DAY2
C       DAY2 IS EXCESS NUMBER OF DAYS IN CURRENT "YEAR".
C       NEWTON RAPHSON SOLUTION FOR E, THE ECCENTRIC ANOMOLY.
        N=0
C       FIRST APPROXIMATION
        E=A
        SINEN=SIN(E)
        HOLD=SINEN
        COSEN=COS(E)
51      N=N+1
        F=E-EC*SINEN-A
        F1=1.-EC*COSEN
        RF=F/F1
        SINEN=-RF*COSEN+SINEN
        COSEN=RF*HOLD+COSEN
        E=E-RF
        RE=ABS(RF/E)
        IF(RE.GT.1.E-05) GO TO 51
C       SINTH=SQRT(1.0-EC*EC)*SINEN
        SINTH=0.99985886*SINEN
C       THIS IS THE "TRUE VALUE"
C*      SINTH=.9997194*SINEN
C       THIS IS GOODSPEED'S VALUE--INCORRECT!
C  HOWEVER THE RESULTANT ERROR IN SUN'S ANGLES IS LESS THAN 0.01 DEG.
        COSTH=COSEN-EC
        FAC=1.0/SQRT(SINTH*SINTH+COSTH*COSTH)
        RADVEC(1)=FAC*COSTH
C       X-COMPONENT OF RADVEC,  E.T.C......
        RADVEC(2)=FAC*SINTH
        RADVEC(3)=0
C
C  OLD SUBR GSOLVC
C       SHELL   6.7.73  FOLLOWS GOODSPEED 5.4.72
CALCULATES SUN'S DIRECTIONAL COSINES IN LOCAL COORDINATE FRAME
C       CODING SLIGHTLY ALTERED TO GIVE CORRECT(!!!) RESULTS.
C       FOR GRIFFITH (NSW, AUSTRALIA)
c but good enough for TOGA-COARE water albedos (PAC)
C
        ARG=60.0*24.0*.4375269285D-02*DAYT+LONG+2.2312617D0
        ARG=ARG/TWOPI
        ARG=DMOD(ARG,1.0D0)
        ARG=ARG*TWOPI
        COSA=COS(ARG)
        SINA=SIN(ARG)
C
C
        DO 70 I=1,3
        OVERH(I)=COSA*GL(I)+SINA*GN(I)
        SINTOM(I)=SINLAT*OM(I)
70      OVERH(I)=OVERH(I)*COSLAT+SINTOM(I)
        CALL GSVEC(OM,OVERH,EAST,4)
        CALL GSVEC(EAST,B,EAST,8)
        CALL GSVEC(OVERH,EAST,NORTH,4)
        CALL GSVEC(RADVEC,NORTH,B,3)
        SL=-B(1)
        CALL GSVEC(RADVEC,EAST,B,3)
        SM=B(1)
        CALL GSVEC(RADVEC,OVERH,B,3)
        SN=-B(1)
C
C
        RETURN
        END
      FUNCTION ACOS(A)
C     COS -VE IN 2,3 QUAD.RETURN PI + ATAN(ANGLE)
        B=1.0-A*A
        IF (B) 5,6,6
6       IF (A) 1,2,3
1     ACOS= 3.14159 + ATAN(SQRT(1.0- A*A)/A)
      RETURN
2     ACOS = 3.14159/2.0
      RETURN
3     ACOS = ATAN(SQRT(1.0-A*A)/A)
      RETURN
5       ACOS=999.
C     THIS INTENDED AS DIAGNOSTIC
      RETURN
      END

C------------------------------------------------------

        SUBROUTINE GSVEC(A,B,C,IOP)

C       SUBROUTINE TO PERFORM VECTOR OPERATIONS
C       SHELL   18.5.73    MODIFIED  8.6.73  CHECKED 12.6.73
C       A,B,C ARE VECTORS, IOP IS VECTOR OPERATION
C       1=ADDITION
C       2=SUBTRACTION
C       3=DOT PRODUCT
C       4=CROSS PRODUCT
C       5=MAGNITUDE
C       6=MULTIPLICATION OF VECTOR BY A SCALAR
C       7=COS(ANGLE)BETWEEN 2 VECTORS, RETURNING C(1) AS ANGLE,
C       AND C(2) AS COS OF THE ANGLE.
C       8=NORMALIZATION OF VECTOR A, RETURNING UNDER SAME NAME
C
C
C
        DIMENSION A(3),B(3),C(3)
        GO TO (1,2,3,4,5,6,7,8), IOP
2       DO 12 I=1,3
        C(I)=A(I)-B(I)
12      CONTINUE
        RETURN
6       DO 16 I=1,3
C  B(1) IS THE SCALAR
16      C(I)=B(1)*A(I)
        RETURN
C
C
7       C(2)=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
C  C(2) IS DOT PRODUCT
C  RMA AND RMB ARE MAGNITUDES OF A AND B RESPECTIVELY
        RMA=SQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3))
        RMB=SQRT(B(1)*B(1)+B(2)*B(2)+B(3)*B(3))
        C(2)=C(2)/(RMA*RMB)
C  C(2) IS COS(ANGLE) BETWEEN VECTORS A AND B
        C(3)=1.0E+38
C       INDEX TO 38 FROM 70 FOR PDP 11/10
C  THIS IS A DIAGNOSTIC
        CTH=C(2)
        C(1)=ACOS(CTH)
C       C(1) LIES BETWEEN ZERO AND PI (3.141593)
        RETURN
C
C
1       DO 11 I=1,3
        C(I)=A(I)+B(I)
11      CONTINUE
        RETURN
3       C(1)=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
        C(3)=1.0E+38
        C(2)=C(3)
C       THIS IS A DIAGNOSTIC WHICH GIVES NONSENSE VALUES
C       IF THE RESULT OF THE DOT PRODUCT IS USED AS A VECTOR.
        RETURN
4       C(1)=A(2)*B(3)-A(3)*B(2)
        C(2)=A(3)*B(1)-A(1)*B(3)
        C(3)=A(1)*B(2)-A(2)*B(1)
        RETURN
5       C(1)=SQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3))
        C(3)=1.0E+38
        C(2)=C(3)
C       THIS IS A DIAGNOSTIC.
        RETURN
8       D=0.0
        DO 18 I=1,3
        D=D+A(I)*A(I)
18      CONTINUE
        IF(D.EQ.0.0) GO TO 80
        D=SQRT(D)
        DO 19 I=1,3
        A(I)=A(I)/D
        C(I)=A(I)
19      CONTINUE
80      RETURN
        END
        
c----------------------------------------------------------------------