function rotangle(lon,lat,coord,rlngd,rlatd,thetad)
c
c=======================================================================
c                                                                    ===
c  This function computes the local angle that the x-coordinate line ===
c  makes with respest to East.                                       ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     LON, LAT  the longitude and latitude of the point (degrees;    ===
c               real)                                                ===
c     COORD     flag for coordinate type (integer)                   ===
c                  COORD = 0  ->  Cartesian (tangent plane) grid.    ===
c                  COORD = 1  ->  unrotated spherical grid.          ===
c                  COORD = 2  ->  rotated spherical grid.            ===
c     RLNGD     longitude at center of grid (degrees; real)          ===
c     RLATD     latitude at center of grid (degrees; real)           ===
c     THETAD    angle grid is rotates with respect to east           ===
c               (degrees; real)                                      ===
c                                                                    ===
c  On Output:                                                        ===
c                                                                    ===
c     ROTANGLE  the local rotation angle  (radians; real)            ===
c                                                                    ===
c     Calls:  ROTPARM                                                ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      logical first,old
      FLOAT
     *      cos1,cosp0,cosp1,lat,lon,londif,plat,plon0,plon1,rlat,rlatd,
     *      rlat0,rlat1,rlngd,rlng0,rlng1,rnlon,rotangle,sin1,sinp0,
     *      sinp1,thetad,theta0,theta1
      integer coord
c
      save cosp0,cosp1,first,plon0,plon1,rlat0,rlat1,rlng0,rlng1,sinp0,
     *     sinp1,theta0,theta1
c
      data first /.true./
      data rlat0,rlng0,theta0,rlat1,rlng1,theta1 /6*c0/
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c  Cartesian (tangent plane) system.
c-----------------------------------------------------------------------
c
      if(coord.eq.0) then
        rotangle=thetad*deg2rad
c
c-----------------------------------------------------------------------
c  Unrotated Spherical system.
c-----------------------------------------------------------------------
c
      elseif(coord.eq.1) then
        rotangle=c0
c
c-----------------------------------------------------------------------
c  Rotated spherical system.
c-----------------------------------------------------------------------
c
      elseif(coord.eq.2) then
c
c  Determine if working with the base rotated spherical system.
c
        old=(thetad.eq.theta0).and.
     *      (rlatd.eq.rlat0).and.(rlngd.eq.rlng0)
        if(old.or.first) then
c
c  Working with base system.  Initialize on first call.
c
          if(first) then
            rlat0=rlatd
            rlng0=rlngd
            theta0=thetad
            first=.false.
            call rotparm(rlngd*deg2rad,rlatd*deg2rad,thetad*deg2rad,
     *                   plat,plon0,rnlon)
            cosp0=cos(plat)
            sinp0=sin(plat)
          endif
          londif=lon*deg2rad-plon0
          rlat=lat*deg2rad
          sin1=cosp0*sin(londif)
          cos1=sinp0*cos(rlat)-sin(rlat)*cosp0*cos(londif)
          rotangle=atan2(sin1,cos1)
        else
c
c  Working with secondary rotated spherical system.  See if its already
c  in storage.
c
          old=(thetad.eq.theta1).and.
     *        (rlatd.eq.rlat1).and.(rlngd.eq.rlng1)
          if(.not.old) then
            rlat1=rlatd
            rlng1=rlngd
            theta1=thetad
            call rotparm(rlngd*deg2rad,rlatd*deg2rad,thetad*deg2rad,
     *                   plat,plon1,rnlon)
            cosp1=cos(plat)
            sinp1=sin(plat)
          endif
          londif=lon*deg2rad-plon1
          rlat=lat*deg2rad
          sin1=cosp1*sin(londif)
          cos1=sinp1*cos(rlat)-sin(rlat)*cosp1*cos(londif)
          rotangle=atan2(sin1,cos1)
        endif
      endif
      return
      end