GEO_ZENITH_ANGLE.f

    SUBROUTINE geo_zenith_angle(i,j,RLAT,RLON,SLAT,SLON,ZA)
!
!************************************************************************
!*                                                    *
!*  Module Name:    GOES_ZA                                   *
!*                                                    *
!*  Language:   Fortran-90     Library:                   *
!*  Version.Rev:    1.0 16 JUL 87   Programmer: KLEESPIES         *
!*          1.1 16 Mar 89           Kleespies             *
!*              Fix bug in ZA computation.                *
!*              Sep. 2010       Chuang: convert to Fortran 90     *
!                                       and incoporate to unified post    *
!*  Calling Seq:    CALL GOES_ZA(RLAT,RLON,SLAT,SLON,ZA)          *
!*  Description:    Computes local zenith angle from a point on the   *
!*      earth's surface to geosynchronous altitude.  This uses    *
!*      spherical geometry, so it isn't fancy.  A ZA gt 90 deg    *
!*      means the earth station is over the horizon from the      *
!*      satellite.                                            *
!*  Input Args: R*4  RLAT  latitude of earth point +- 90 deg      *
!*          R*4  RLON  longitude of earth point deg + E       *
!*          R*4  SLAT  latitude of subsatellite point     *
!*          R*4  SLON  longitude of subsatellite point    *
!*                                    *
!*  Output Args:    R*4  ZA    zenith angle to geostationary      *
!*                 altitude.                  *
!*                                    *
!*  Common Blks:    none                          *
!*                                    *
!*  Include:    none                          *
!*                                    *
!*  Externals:  none                          *
!*                                    *
!*  Data Files: none                          *
!*                                    *
!*  Restrictions:   Uses spherical geometry, so is not totally    *
!*      accurate.  Makes no check for quadrant changes, so may    *
!*      have problems around the prime meridian and the date      *
!*      line.  Also the poles.                    *
!*                                    *
!*  Error Codes:    none                          *
!*                                    *
!*  Error Messages: none                          *
!*                                    *
!************************************************************************

!        use params_mod, only: DTR,ERAD,PI,RTD
    IMPLICIT NONE
    REAL*8 dtr/0.017453293 /   ! degrees to radians
    REAL*8 re / 6370.0/        ! radius of earth
    REAL*8 pi /3.1415926/      ! pi
    REAL*8 rtd /57.29577951 /  ! radians to degrees
!   REAL*8 RZS /1786245696.0/  ! (geosynchronous height)**2
    REAL*8 res /40589641.0  /  ! (earth radius)**2
    REAL*8 c1  /1826835337.0/  ! RES + RZS
    REAL*8 c2  /538527888.0 /  ! 2*RE*RES
    
    REAL rzs /1786245696.0/  ! (geosynchronous height)**2 for GOES

!   REAL*8 A,B,C,COSD,COSE,E,P,PP
        INTEGER, INTENT(IN):: i,j
        REAL,INTENT(IN):: slon,rlat,rlon,slat
    REAL,INTENT(OUT):: za
!        REAL C1,C2,RES,RE
    REAL a,b,c,cosd,cose,e,p,pp
    
!   RE=ERAD/1.0E3
!   RES=RE**2.0
!   C1=RES+RZS
!   C2=2.0*RE*RES/1000.

    a = (90.0 - rlat) * dtr  !  colatitude of point
    b = (90.0 - slat) * dtr  !  colatitude of satellite
    
    IF(rlon>180.)THEN
     c = abs(rlon- 360. - slon) * dtr
    ELSE 
     c = abs(rlon - slon) * dtr 
    END IF 
    
    cosd = cos(a)*cos(b) + sin(a)*sin(b)*cos(c) ! great circle arc

    pp = c1 - c2*cosd
    p = sqrt(pp)

    cose = (pp + res - rzs) / (2.*re*p)
        cose=max(min(cose,1.),-1.)
    e = acos(cose)

    za = pi - e
    za = za * rtd
        za=max(za,0.)
    if(abs(rlon-360.-slon)<1. .and. abs(rlat-30.)<1.)print*,'Debug GEO_ZENITH',  &
    rlon,rlat,res,c1,c2,a,b,c,cosd,pp,p,cose,e,za

    RETURN
    END