subroutine xy2ll (x,y,coord,nx,ny,gridx,gridy,rlngd,rlatd,delx,
     &                  dely,thetad,lon,lat)
c
c=======================================================================
c                                                                    ===
c  This routine computes the longitude and latitude of a given       ===
c  point from the supplied grid coordinates.                         ===
c                                                                    ===
c  ------                                                            ===
c  Input:                                                            ===
c  ------                                                            ===
c                                                                    ===
c  COORD     flag for coordinate type.                 (integer)     ===
c               [0]  Cartesian (tangent plane) grid.                 ===
c               [1]  unrotated spherical grid.                       ===
c               [2]  rotated spherical grid.                         ===
c  DELX,                                                             ===
c  DELY      Offset from transformation point to                     ===
c            grid center.                              (real; cm|deg)===
c  GRIDX     grid spacing in x-direction               (real; cm|deg)===
c  GRIDY     grid spacing in y-direction               (real; cm|deg)===
c  NX, NY    number of grid points in x & y dir.s.     (integer)     ===
c  RLNGD     longitude at center of grid               (real; deg)   ===
c  RLATD     latitude at center of grid                (real; deg)   ===
c  THETAD    angle grid is rotated WRT east.           (real; deg)   ===
c     X, Y      the grid coordinates (real)                          ===
c                                                                    ===
c  -------                                                           ===
c  Output:                                                           ===
c  -------                                                           ===
c                                                                    ===
c  LON, LAT  Longitude and latitude of the point.      (real; deg)   ===
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 new
      FLOAT
     &      a,c2pi,cos1,cosp,cosr,cterm,deg2cm,delx,delxo,dely,delyo,dx,
     &      dy,gdxo,gdyo,gridx,gridy,lat,lon,plat,plon,rf1x,rf1y,rf2x,
     &      rf2y,rlat,rlatd,rlato,rlon,rlngd,rlngo,rnlon,sin1,sinp,sinr,
     &      sterm,thetad,thetao,x,xc,y,yc
      FLOAT
     &        lon_chk
      integer coord,coordo,nx,ny,nxo,nyo
      parameter (c2pi=pi+pi,deg2cm=re*deg2rad)
c
      save coordo,cosp,cterm,delxo,delyo,gdxo,gdyo,nxo,nyo,plat,plon,
     &     rf1x,rf1y,rf2x,rf2y,rlato,rlngo,rnlon,sinp,sterm,thetao,xc,yc
c
      data coordo,nxo,nyo /3*0/
      data delxo,delyo,gdxo,gdyo,rlato,rlngo,thetao /7*c0/
c
c-----------------------------------------------------------------------
c  Longitude range check function.
c-----------------------------------------------------------------------
c
      lon_chk(a)=sign(c1,a)*
     &           (mod(abs(a),c2pi)-c2pi*int(mod(abs(a),c2pi)/pi))
c
c=======================================================================
c Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c Determine if working with a new coordinate system.
c-----------------------------------------------------------------------
c
      new = (nx.ne.nxo).or.(ny.ne.nyo).or.(thetad.ne.thetao).or.
     &      (rlatd.ne.rlato).or.(rlngd.ne.rlngo).or.(gridx.ne.gdxo).or.
     &      (gridy.ne.gdyo).or.(coord.ne.coordo).or.(delx.ne.delxo).or.
     &         (dely.ne.delyo)
      if (new) then
        coordo = coord
        gdxo   = gridx
        gdyo   = gridy
        delxo  = delx
        delyo  = dely
        nxo    = nx
        nyo    = ny
        rlato  = rlatd
        rlngo  = rlngd
        thetao = thetad
      endif
c
c-----------------------------------------------------------------------
c  Cartesian (tangent plane) system.
c-----------------------------------------------------------------------
c
      if (coord.eq.0) then
        if (new) then
          rf1x  = gridx/deg2cm
          rf1y  = gridy/deg2cm
          rf2x  = rf1x/cos(rlatd*deg2rad)
          rf2y  = rf1y/cos(rlatd*deg2rad)
          sterm = sin(thetad*deg2rad)
          cterm = cos(thetad*deg2rad)
          xc    = p5*FLoaT(nx+1)-delx/gridx
          yc    = p5*FLoaT(ny+1)-dely/gridy
        endif
        dx  = x-xc
        dy  = y-yc
        lon = rlngd+rf2x*dx*cterm-rf2y*dy*sterm
        lat = rlatd+rf1y*dy*cterm+rf1x*dx*sterm
c
c-----------------------------------------------------------------------
c  Unrotated spherical system.
c-----------------------------------------------------------------------
c
      elseif (coord.eq.1) then
        if (new) then
          xc = p5*FLoaT(nx+1)-delx/gridx
          yc = p5*FLoaT(ny+1)-dely/gridy
        endif
        lon = (x-xc)*gridx
        lat = (y-yc)*gridy
c
c-----------------------------------------------------------------------
c  Rotated spherical system.
c-----------------------------------------------------------------------
c
      elseif (coord.eq.2) then
        if (new) then
          call rotparm (rlngd*deg2rad,rlatd*deg2rad,thetad*deg2rad,plat,
     &                  plon,rnlon)
          xc   = p5*FLoaT(nx+1)-delx/gridx
          yc   = p5*FLoaT(ny+1)-dely/gridy
          cosp = cos(plat)
          sinp = sin(plat)
        endif
        rlon = (x-xc)*gridx
        rlat = (y-yc)*gridy
        cosr = cos(rlat*deg2rad)
        sinr = sin(rlat*deg2rad)
        lat  = asin(sinr*sinp+cosr*cosp*cos(rlon*deg2rad-rnlon))
        sin1 = cosr*sin(rlon*deg2rad-rnlon)
        cos1 = sinr*cosp-sinp*cosr*cos(rlon*deg2rad-rnlon)
        lon  = plon-atan2(sin1,cos1)
        lon  = lon_chk(lon)
        lat  = lat*rad2deg
        lon  = lon*rad2deg
      endif
c
      return
      end