PROGRAM GRIDL
C
C Compute and print the lat, lon, and grid spacing for grid points
C in various polar stereographic grids.  To run for a specific
C grid the proper array size must be selected below and the proper
C grid definition parameters must be selected in the GRDLOC
C subroutine.
C
C     PARAMETER (M=41,N=38)
C     PARAMETER (M=513,N=513)
      PARAMETER (M=65,N=65)
C     PARAMETER (M=63,N=63)
C     PARAMETER (M=47,N=51)
C     PARAMETER (M=147,N=110)
      DIMENSION ALA(M,N),ALO(M,N),AGSP(M,N)
      XP=21.
      YP=18.
C GRID LOCATION OF ARBITRARY POINT
      CALL GRDLOC (XP,YP,XLO,XLA,GSP)
      PRINT 1004,XP,YP,XLO,XLA,GSP
 1004 FORMAT(' AT ',2F8.2,' LON,LAT IS ',3F12.4)
C GRID SPACING AT VARIOUS LATITUDES
      DO 8 LAT=0,90,10
      XLAT=LAT
      CALL LATGSP(XLAT,GSP)
      PRINT 1006,XLAT,GSP
 1006 FORMAT(' LAT, SPACING ',2F12.3)
    8 CONTINUE
C GRID LOCATION AND SPACING FOR POINTS ON THE GRID
      DO 10 I=1,M
      XP=I
      DO 10 J=1,N
      YP=J
      CALL GRDLOC (XP,YP,XLO,XLA,GSP)
      ALA(I,J)=XLA
      ALO(I,J)=XLO
      AGSP(I,J)=GSP
   10 CONTINUE
      PRINT 1000,M,N
 1000 FORMAT(' GRID POINT LOCATIONS - ',I4,'  X ',I4)
      DO 20 IS=1,M,14
      IE=IS+13
      IF(IE.GT.M) IE=M
      PRINT 1003,(I,I=IS,IE)
 1003 FORMAT(9X,16I8)
      DO 18 JJ=1,N
      J=N-JJ+1
      PRINT 1001,J,(ALA(I,J),I=IS,IE)
 1001 FORMAT(1X,I3,' LAT ',14F8.2)
      PRINT 1002,  (ALO(I,J),I=IS,IE)
 1002 FORMAT(4X,' LON ',14F8.2)
      PRINT 1005,  (AGSP(I,J),I=IS,IE)
 1005 FORMAT(4X,' SPC ',14F8.2)
C
   18 CONTINUE
   20 CONTINUE
      END
      SUBROUTINE GRDLOC(XP,YP,XLO,XLA,GSP)
C FIND LON,LAT OF GIVEN GRID POINTS (XP,YP)
C
C GRID DEPENDENT PARAMETERS FOLLOW -
C
C    THE POLE POINT OF THE GRID IS XPP,YPP
C    ER IS EARTH'S RADIUS IN GRID UNITS
C    XAX IS THE LONGITUDE (W IS NEGATIVE) OF THE +X AXIS OF THE GRID
C
C NAVY 63x63 grid parameters
C     PARAMETER(XPP=32.,YPP=32.,XAX=10.,ER=31.204359052)
C NMC 65x65 grid parameters
      PARAMETER(XPP=33.,YPP=33.,XAX=10.,ER=31.204359052)
C 190km lfm and ngm paramters
C     PARAMETER(XPP=19.0,YPP=42.0,XAX=-15.,ER=2.*31.204359052)
C 47x51 grid parameters
C     PARAMETER(XPP=24.,YPP=26.,XAX=10.,ER=31.204359052)
C NEPH GRID - GIVES LAT LON OF CORNER OF GRID BOX
C     PARAMETER(XPP=256.,YPP=256.,XAX= 10.,ER=8.*31.204359052)
C     PARAMETER(XPP=24.0,YPP=24.0,XAX=000.,ER=23.403269292)
C NCEP GRID 104
C     PARAMETER(XPP=75.5,YPP=109.5,XAX=-15.,ER=131.0)
C
      PARAMETER(ER2=ER*ER,EAR=6371.2213)
C INPUT
C
C       XP,YP GRID POINT FOR WHICH LON,LAT IS DESIRED
C
C NOTE --
C  IF PROJECTION IS LOOKING DOWN ON NORTH POLE
C       (COUNTERCLOCKWISE AROUND POLE IS EAST)
C       EAST LONGITUDES ARE CONSIDERED POSITIVE (EASTWARD INCREMENTS POSITIVE)
C  IF PROJECTION IS LOOKING DOWN ON SOUTH POLE
C       (COUNTERCLOCKWISE AROUND POLE IS WEST)
C       WEST LONGITUDES ARE CONSIDERED POSITIVE (WESTWARD INCREMENTS POSITIVE)
C
C
C OUTPUT
C       XLON,XLAT LONGITUDE AND LATITUDE OF SPECIFIED GRID POINT.
C       GSP IS GRID SPACING AT THIS POINT
C
C INTERNAL
C
C  R IS RADIUS IN GRID UNITS FROM POLE POINT OF DESIRED POINT
C  PHI IS LATITUDE IN RADIANS
C  ELONG IS LONG EAST OF XAX (+X AXIS ON GRID) -- IN RADIANS
C
      PARAMETER (PI=3.14159265,PI4=PI/4,RADDEG=180./PI)
C
      YY=YP-YPP
      XX=XP-XPP
      XLO=0.
      IF(YY.EQ.0. .AND. XX.EQ.0.) GO TO 10
      ELONG=RADDEG*ATAN2(YY,XX)
      XLO=ELONG+XAX
      IF(XLO.GT.180.) XLO=XLO-360.
      IF(XLO.LT.-180.) XLO=XLO+360.
   10 CONTINUE
      R2=XX*XX+YY*YY
      SS=(ER2-R2)/(ER2+R2)
      XLA=RADDEG*ASIN(SS)
      GSP=EAR*(1.+SS)/ER
      RETURN
      ENTRY LATGSP(XLAT,GSP)
C RETURN GRID SPACING IN KM (GSP) AT A GIVEN LATITUDE (XLAT)
      SS=SIN(XLAT/RADDEG)
      GSP=EAR*(1.+SS)/ER
      RETURN
      END