program poles_frc_therm_grid
c  Functionality: 
C     PROGRAM which demononstrates how to  get (x,y) and (lat.,lon.)
c     for the 43x42 grid of the: Thermal Frocing Fieldsc   
c
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
c Parameters for the Thermal Grid  specic grid
c     DX, DY: Grid interval
c     NX, NY: Number of grid cells in x and y direction. 

      PARAMETER (NX=44,NY=43,DX=160.,DY=160.)
      DIMENSION DLON(NX,NY),DLAT(NX,NY),X(NX),Y(NY)

      radius=6370.
      iday=0
      do 2 i=1,nx
         do 3 j=1,ny
            x(i)=((i-1)*1.0)*dx-4080.0
            y(j)=((j-1)*1.0)*dy-(3280.0+160.)

            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)))
c            print*,'i,j,x,y,lat,lon'
            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