c********************************************************
c Name:  genhrap.f
c
c Purpose:  Write the HRAP coordinates for a selected region of 
c cells to a file and create a subsequent file in geographic 
c coordinates in a suitable format to serve as input to the 
c GENERATE (polygon) command in ARC/INFO.  This program is designed 
c to be followed by genhrap.aml.
c
c Two options are available for defining the region of cells to be 
c created --(1) Define the latitude and longitude extent of the region 
c to be mapped, or (2) Specify the SW corner of the grid to be created 
c and the number of colums and rows of cells to be created.  With 
c either option, the program computes the HRAP coordinates of the 
c SW corner (if necessary) and generates grid cells starting with 
c the bottom row, moving from left to right, and then 
c moving to the next row up and repeating.
c
c Comments:  Only the output files hrap.COD.dat and inputgc.COD are
c required as input to genhrap.aml.  Intermediate files and optional
c files that were created in an earlier version are also listed below.
c Code can be compiled on a UNIX F77 compiler.  I do not know whether 
c the syntax is compatible with other compilers.
c
c Calls subroutines:  wll, topoly, crdat(numx,numy,xstart,ystart)
c 
c Inputs: none
c Output: "COD" is a user defined suffix
c     hrap.COD = file of hrap coordinates  /*temporary
c     geoc.COD = file of geocentric coordinates  /*temporary 
c     hrap.COD.dat = file containing HRAP coordinates in a format that
c                    can be attached to the polygon attribute table
c     *pster.COD = file of polar stereographic coordinates
c     inputgc.COD = input file of geoc. coordinates to make a polgon 
c        coverage 
c     *inpster.COD = input file of polar stereographic coordinates to
c                    make a polygon coverage
c     *inhrap.COD = input file of HRAP coordinates to make a polygon 
c                   coverage
c
c A * denotes optional files -- lines corresponding to the creation 
c of optional files have been commented out 
c in this version of the code.
c**********************************************************************

      program genhrap

c     <<< Variable Declaration >>>
      parameter (maxcol = 336, maxrow = 160)
c    *** maxcol and maxrow are limited to the extent of HRAP cells for
c       which data is available in the Arkans.-Red River Basin

      integer xstart,ystart,numx,numy,numpts,numx1,numy1
      double precision xhrap(maxcol), yhrap(maxrow)
      integer count,bool,rfunit,wfunit
c    *** rfunit and wfunit store the readfile unit number and the 
c    *** writefile unit number to be passed to the subroutine topoly.
      character suff*3,file1*8,file2*8,file3*12,file5*11
c    ***
c    <<< End of variable declaration >>>

c    *** Allow two options for defining the study region.
      write(*,*) 'Enter 1 to specify your region by latitudes'
      write(*,*) 'and longitudes of its corners.'
      write(*,*) ''
      write(*,*) '-OR-'
      write(*,*) ''
      write(*,*) 'Enter 2 to specify your region by hrap coordinates'
      write(*,*) 'and number of columns and rows.'
     
      read(*,*) bool    
      if (bool.eq.1) then
         call llinput(xstart,ystart,numx1,numy1)
        else

      write(*,*) 'Enter the hrap(x,y) for the lower left hand'
      write(*,*) 'corner of the region of interest:' 
      
      read(*,*) xstart,ystart
      write(*,*) 'Enter the number of grid columns and rows to be' 
      write(*,*) 'created:'
      read(*,*) numx1,numy1
      endif

c    *** Number of points to write is one greater than the number of 
c    *** columns or rows.  The name numx1 can be thought of as number of
c    *** x coordinates - 1.
    
      numx = numx1 + 1
      numy = numy1 + 1
      write(*,*) 'Enter a 3 character suffix to uniquely identify'
      write(*,*) 'your grid:'
      read(*,*) suff

c    ***Create names for all of the output files.
c    *** file1 = file of hrap coordinates
c    *** file2 = file of geocentric coordinates
c    *** file3 = file containing HRAP coordinates in a format that can be 
c                attached to the polygon attribute table
c    *** file4 = file of polar stereographic coordinates
c    *** file5 = input file of geoc. coordinates to make a polgon
c                coverage
c    *** file6 = input file of p. stereographic coordinates to make a 
c                polgon coverage
c    *** file7 = input file of hrap coordinates to make a polgon coverage

      file1 = 'hrap.'//suff
      file2 = 'geoc.'//suff
      file3 = 'hrap.'//suff//'.dat'
c     file4 = 'pster.'//suff
      file5 = 'inputgc.'//suff
c     file6 = 'inpster.'//suff
c     file7 = 'inhrap.'//suff
    
 
      open(unit = 10, file = file1, status = 'unknown')
      open (unit = 20, file = file2, status = 'unknown')
      open (unit = 30, file = file3, status = 'unknown')
c      open (unit = 40, file = file4, status = 'unknown') 
      open (unit = 50, file = file5, status = 'unknown')
c      open (unit = 60, file = file6, status = 'unknown')
c      open (unit = 70, file = file7, status = 'unknown')
    
c    *** Compute the total number of cell corners
      numpts = numx*numy

      xnew = xstart

      do 100 i=1,numx
         xhrap(i) = xnew
         xtemp = xnew + 1.0
         xnew = xtemp
 100  continue
    
      ynew = ystart
      do 200 j=1,numy
         yhrap(j) = ynew
         ytemp = ynew + 1.0
         ynew = ytemp
 200  continue

      count = 1
      do 300 j=1,numy
         do 400 i=1,numx
c          ** with "count" in the file, this file can be used as input
c          ** for creating an ARC/INFO point coverage
            write(10,*) count,xhrap(i),yhrap(j)
            count = count + 1
 400     continue
 300     continue
c    write(*,*) ystart

      call wll(numpts)
      rfunit = 20
      wfunit = 50
      call topoly(numx,numy,rfunit,wfunit)
c     rfunit = 40
c     wfunit = 60
c     call topoly(numx,numy,rfunit,wfunit)
c     rfunit = 10 
c     wfunit = 70
c     call topoly(numx,numy,rfunit,wfunit)
      call crdat(numx,numy,xstart,ystart)
      close(10)
      close(20)
      close(30)
      close(50)
      stop
      end  

c<<<<<<<<<<<<<<<<<<<<<<<SUBROUTINES>>>>>>>>>>>>>>>>>>>>>>>>>
c *********************************************************************
c
c Purpose: Convert HRAP coordinates contained in file "hrap.COD" to 
c lat-long coordinates based on a spherical earth and write them to an 
c output file that can be used to generate a point coverage.  The 
c parameter "numpts" stores the number of entries that will be expected 
c from "hrap.COD"
c
c Inputs: file hrap.COD
c Outputs: file geoc.COD
c
c*********************************************************************
      subroutine wll(numpts)
      double precision xhrap, yhrap, x, y
      double precision bigr, arg, latd, lond, ang
      double precision stlatd,earthr,mesh,stlond 
      integer rec,numpts
c
c***    Define constants    
      stlond = -105.0
      stlatd = 60.0
c*** earthr,mesh, x, and y are in meters.
      earthr = 6371200.0
      mesh = 4762.5

      rewind(unit=10)
    
      do 100 i=1,numpts
       
        read(10,*) rec,xhrap, yhrap
        x = (xhrap - 401.0)*mesh
        y = (yhrap - 1601.0)*mesh
    
        bigr = (x*x + y*y)**0.5
        arg = bigr/(earthr*(1 + dsind(stlatd)))
        latd = 90.0 - 2*datand(arg)

        ang = datan2d(y,x)

        if (y.gt.0) then 
          ang = 270.0-stlond-ang
        else
          ang = -90.0-stlond-ang
        endif
        if (ang.lt.180) then
          lond = -1 * ang
        else
          lond = 360.0 - ang
        endif
c***    Write polar stereographic coordinates and geocentric 
c***    coordinates to a file.
c        write(40,*) i,x,y
        write(20,*) i,lond, latd
 100  continue
      return
      end

c123456789012345678901234567890123456789012345678901234567890123456789012
c**********************************************************************
c Purpose: Given a list of corner points for a grid (can be (ID,x,y) or
c          (ID, lon,lat) in which the coordinates for the bottom row are 
c          listed one per line followed by the coordinates for the next
c          row up, create a file that can be used to generate a polygon
c          coverage of the grid cells.
c
c Input: File of corner points (ID,x,y),
c Ouput: File with lines: "poly-id, ll,lr,ur,ul,ll,end" -- repeated for
c        each polygon.  ll = lower left, lr = lower right, ur = upper 
c        right,ul = upper left
c           
c**********************************************************************
      subroutine topoly(numx,numy,rfunit,wfunit)

c      <<< Variable Declaration >>>
c      parameter (numx = 20, numy = 20)
c***    The old number of x-coordinates was 336.
c***    The old number of y-coordinates was 160.

      double precision xrowa(336),yrowa(336),xrowb(336),yrowb(336)
c       ** xrowa, yrowa are x and y coordinates of points in row a
      character*3 end
      integer i,l,rcount,r,polynum,numx,numy
      integer rfunit,wfunit
c    <<< End of Variable Declaration >>>

      end = 'end'

      rewind(unit=rfunit)
      rcount = 1
      polynum = 1
    
      do 200 i=1,numx
         read(rfunit,*) rec,xrowa(i),yrowa(i)
 200  continue

 100  if (rcount.lt.numy) then

      do 250 i=1,numx
         read(rfunit,*) rec,xrowb(i),yrowb(i)
 250  continue

      l = 1
 300  if (l.lt.numx) then
         r = l + 1
         write(wfunit,*) polynum, xrowa(l), yrowa(l)
         write(wfunit,*) xrowa(l),yrowa(l)
         write(wfunit,*) xrowa(r),yrowa(r)
         write(wfunit,*) xrowb(r),yrowb(r)
         write(wfunit,*) xrowb(l),yrowb(l)
         write(wfunit,*) xrowa(l),yrowa(l)
         write(wfunit,*) end
         l = l + 1
         polynum = polynum + 1
         goto 300
         endif

         rcount = rcount + 1
         do 350 i=1,numx
           xrowa(i) = xrowb(i)
           yrowa(i) = yrowb(i)
 350     continue
       goto 100
      endif
      write(wfunit,*) end

      return
      end
   
c **********************************************************************
c Purpose: This subprogram will create a data file that can be joined to
c the projected "hrap" polygon coverage so that "hrap" coordinates of the
c lower left hand corner of each polygon will be added to the appropriate
c line in the PAT.
c
c Note: The only difference between "hrap.COD.dat" produced by this 
c subroutine and "hrap.COD" produced by the main program is that 
c hrap.COD.dat does not contain entries for the last column and last
c row of points.
c **********************************************************************
    
      subroutine crdat(numx,numy,xstart,ystart)
c***    Old value of numx was 336
c***    Old value of numy was 160
      double precision xhrap(336), yhrap(160)
      integer count,numx,numy,xstart,ystart,numx1,numy1
    
      numx1 = numx - 1
      numy1 = numy - 1
      xnew = xstart

      do 100 i=1,numx1
         xhrap(i) = xnew
         xtemp = xnew + 1.0
         xnew = xtemp
 100  continue
    
      ynew = ystart
      do 200 j=1,numy1
         yhrap(j) = ynew
         ytemp = ynew + 1.0
         ynew = ytemp
 200  continue

      count = 1
      do 300 j=1,numy1
         do 400 i=1,numx1
            write(30,*) count,xhrap(i),yhrap(j)
            count = count + 1
 400     continue
 300  continue
      return
      end  

c**********************************************************************
c    At user's request, allow the user to input the latitude and 
c    longitude of the four corners that are of interest in the 
c    study.
c  
c    Note: The user should input geodetic coordinates.  These 
c    geodetic coordinates will be interpreted as geocentric coordinates
c    to be consistent with methodology used by the 
c    National Weather Service.
c**********************************************************************
      subroutine llinput(xstart,ystart,numx1,numy1)
    
c    <<< Variable Declaration >>>
      parameter (stlat = 60.0)
c*** clon is a constant used to account for the standard longitude
c*** see eqn. in "Geographic Positioning of the HRAP"
      parameter (clon = 15.0)
      parameter (rad = 6371.2)

      integer xstart,ystart,numx1,numy1
      real lon(4), lat(4)
      real sfactor,R,x,y,hrapx(4),hrapy(4)
c*** Declare variables llhrapx and llhrapy to pick the hrap coordinates of 
c*** the lower left hand coordinates desired.
      real minhx,minhy,maxhx,maxhy
c    <<< End Variable Declaration >>>

      write(*,*) 'Enter the latitudes and longitudes of four' 
      write(*,*) 'corners of a rectangle that encloses the' 
      write(*,*) 'study region (in decimal degrees).' 
      write(*,*) ''
      write(*,*) 'Enter a longitude value and then a space' 
      write(*,*) 'and then a latitude value. Hit return after'
      write(*,*) 'each coordinate.  Remember to input West'
      write(*,*) 'longitude values as negative numbers.'

      do 100 i = 1,4

        read(*,*) lon(i),lat(i)  
        sfactor = (1+sind(stlat))/(1+sind(lat(i)))
c** x and y are in km
        R = rad*cosd(lat(i))*sfactor
        x = R*cosd(lon(i)+clon)
        y = R*sind(lon(i)+clon)
        hrapx(i) = x/4.7625 + 401
        hrapy(i) = y/4.7625 + 1601
        write(*,*) 'hrapx, hrapy:', hrapx(i), hrapy(i)
100   continue
        minhx = hrapx(1)
        minhy = hrapy(1)
        maxhx = hrapx(1)
        maxhy = hrapy(1)

      do 200 j = 2,4
       
         if (hrapx(j).lt.minhx) then
            minhx = hrapx(j)
         endif
         if (hrapy(j).lt.minhy) then
            minhy = hrapy(j)
         endif
         if (hrapx(j).gt.maxhx) then
            maxhx = hrapx(j)
         endif
         if (hrapy(j).gt.maxhy) then
            maxhy = hrapy(j)
         endif

 200  continue 
      xstart = minhx
      ystart = minhy

      numx1 = maxhx - minhx
      numy1 = maxhy - minhy
      write(*,*) 'Lower left, num rows, num columns'
      write(*,*) xstart,ystart,numx1,numy1
      return
      end 
       
c*********************************************************************