program poles_frc_wind_grid
c  Functionality: 
C     PROGRAM which demononstrates how to  get (x,y) and (lat.,lon.)
c     for the 43x42 grid of the: 
c     Geostrophic wind fields. 
c   Output: (users have to modify to get output) 
c     X, Y      -  Map projection coordinates in kilometers
c                  for each grid cell center. 
c     Dlat,Dlon -  Latitude, Longitude coordinates in degrees
c                  for each grid cell center. 
c Parameters for this specic grid
c     DX, DY: Grid interval
c     NX, NY: Number of grid cells in x and y direction. 
      PARAMETER (NX=43,NY=42,DX=160.,DY=160.)
      DIMENSION DLON(NX,NY),DLAT(NX,NY),X(NX),Y(NY)

C   Use Lambert Equival. Projec.
C   Generate Grid Coordinates (x-y) and lon./lat. (x,y) are in km

c         
      radius=6370.
      iday=0
      do 2 i=1,nx
         do 3 j=1,ny
            x(i)=((i-1)*1.0)*dx-4080.0+80.0
            y(j)=((j-1)*1.0)*dy-(3280.0+160.)+80.0
            rho=sqrt(x(i)*x(i)+y(j)*y(j))
            if (x(i).ge.0.) then
               dlon(i,j)=35.+(asind(y(j)/rho))
            else if (x(i).lt.0..and.y(j).ge.0.) then
               xita=asind(y(j)/rho)
               dlon(i,j)=35.+(180.-xita)
               if ( dlon(i,j).gt.180. ) then
                  dlon(i,j)=35.+(-180.-xita)
               end if
            else if (x(i).lt.0..and.y(j).lt.0.) then
               dlon(i,j)=35.+(-180.-asind(y(j)/rho))
            end if
            dlat(i,j)=90.-2*(asind(rho/(2*radius)))
            print 99,i,j,x(i), y(i),dlat(i,j),dlon(i,j)
 99         format(2(i2,1x),2(f8.2,1x),2(f7.2,1x))
 3       continue
 2    continue
      end