subroutine ll2xy (lon,lat,coord,nx,ny,gridx,gridy,rlngd,rlatd,
     &                  delx,dely,thetad,x,y)
c
c=======================================================================
c                                                                    ===
c  This routine computes the grid coordinates from the given         ===
c  longitude and latitude.                                           ===
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  LON, LAT  Longitude and latitude of the point.      (real; 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                                                                    ===
c  -------                                                           ===
c  Output:                                                           ===
c  -------                                                           ===
c                                                                    ===
c     X, Y      the corresponding grid coordinates (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 new
      FLOAT
     &      a,c2pi,cos1,cosl,cosp,cterm,deg2cm,delx,delxo,dely,delyo,
     &      dlat,dlon,f1x,f1y,f2x,f2y,gdxo,gdyo,gridx,gridy,lat,lon,plat
     &     ,plon,rlat,rlatd,rlato,rlon,rlngd,rlngo,rnlon,sin1,sinl,sinp,
     &      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,f1x,f1y,f2x,f2y,gdxo,gdyo,nxo,
     &     nyo,plat,plon,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
          f1x   = deg2cm/gridx
          f1y   = deg2cm/gridy
          f2x   = f1x*cos(rlatd*deg2rad)
          f2y   = f1y*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
        dlat = lat-rlatd
        dlon = lon-rlngd
        x    = xc+f1x*dlat*sterm+f2x*dlon*cterm
        y    = yc+f1y*dlat*cterm-f2y*dlon*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
        x    = xc+lon/gridx
        y    = yc+lat/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
        cosl = cos(lat*deg2rad)
        sinl = sin(lat*deg2rad)
        rlat = asin(sinl*sinp+cosl*cosp*cos(lon*deg2rad-plon))
        sin1 = cosl*sin(lon*deg2rad-plon)
        cos1 = sinl*cosp-sinp*cosl*cos(lon*deg2rad-plon)
        rlon = rnlon-atan2(sin1,cos1)
        rlon = lon_chk(rlon)
        rlon = rlon*rad2deg
        rlat = rlat*rad2deg
        x    = xc+rlon/gridx
        y    = yc+rlat/gridy
      endif
c
      return
      end