PROGRAM GRIDL
C      PARAMETER (M=41,N=38)
C     PARAMETER (M=513,N=513)
C     PARAMETER (M=63,N=63)
C     PARAMETER (M=65,N=65)
      PARAMETER (M=47,N=51)
      DIMENSION ALA(M,N),ALO(M,N)
      XP=10.
      YP=13.
      CALL GRDLOC (XP,YP,XLO,XLA)
      PRINT 1004,XP,YP,XLO,XLA
 1004 FORMAT(' AT ',2F8.2,' LON,LAT IS ',2F16.6)
      DO 10 I=1,M
      XP=I
      DO 10 J=1,N
      YP=J
      CALL GRDLOC (XP,YP,XLO,XLA)
      ALA(I,J)=XLA
      ALO(I,J)=XLO
   10 CONTINUE
      PRINT 1000,M,N
 1000 FORMAT(' GRID POINT LOCATIONS - ',I4,'  X ',I4)
      DO 20 IS=1,M,15
      IE=IS+14
      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 ',15F8.2)
      PRINT 1002,  (ALO(I,J),I=IS,IE)
 1002 FORMAT(4X,' LON ',15F8.2)
C
   18 CONTINUE
   20 CONTINUE
      END
      SUBROUTINE GRDLOC(XP,YP,XLO,XLA)
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       AUSTRALIAN 47X47     23.403269292      24,24     0     47,47     -1
C     PARAMETER(XPP=32.,YPP=32.,XAX=10.,ER=31.204359052)
C     PARAMETER(XPP=33.,YPP=33.,XAX=10.,ER=31.204359052)
C     PARAMETER(XPP=19.0,YPP=42.0,XAX=-15.,ER=2.*31.204359052)
      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
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       LON,LAT LONGITUDE AND LATITUDE OF SPECIFIED GRID 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,A4=ER*ER)
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
      XLA=RADDEG*ASIN((A4-R2)/(A4+R2))
      RETURN
      END