c----------------------------------------------------------------------
      subroutine getproj(iprj,cm,bl)
c     Requests projection information including cental meridian and
c      base latitude.
      print '(/,(1x,a))',
     &  ' Available Projections...',
     &  ' ',
     &  ' 0=none, 1=poly, 2=transverse mercator="UTM"',
     &  '  3=merc, 4=lambert US, 14=general lambert',
     &  '   5=albers US, 6=albers Alaska, 7=albers Hawaii, 12=DNAG',
     &  '    99=proj not known'
      call ask(' Projection (0 to 14):  ')
      read *, iprj
      if(iprj.gt.0.and.iprj.lt.99) then
        call ask(' CM, BL (d m, d m) =       ')
        read *, cmd,cmm,bld,blm
        cm=-abs(decdeg(cmd,cmm))
c        this forces west longitude
        bl=decdeg(bld,blm)
      endif
      return
      end
c-----------------------------------------------------------------------------
      function int_ask(request)
      character*(*) request
    2 write(6,100)request
  100 format(a,' ',$)
      read(5,*,err=1)int_ask
      return
    1 write(6,101)
  101 format(/,' **ERROR - try again')
      call bell
      go to 2
      end
c--------------------------------------------------------------------------
      subroutine prheader(unit,id,pgm,ncol,nrow,nz,xo,dx,yo,dy,
     &    iproj,cm,bl)
      character id*56, pgm*8
      integer unit
      write(unit,1001) id,pgm
 1001 format(/,'    ID=',a56,'   pgm=',a8)
      write(unit,1002) 'nz  =',nz
      write(unit,1002) 'ncol=',ncol,'xo=',xo,'dx=',dx
      write(unit,1002) 'nrow=',nrow,'yo=',yo,'dy=',dy
 1002 format(1x,a8,i8,a8,g14.8,a8,g14.8)
      write(unit,'(1h0,a,i2,2(a,f9.4))')
     & '   Projection=',iproj,'     cm=',cm,'          bl=',bl
      write(unit,1003)
 1003 format(1h )
      return
      end
c------------------------------------------------------------------------
      subroutine putrow(iounit,ncol,row)
      real row(ncol)
      parameter (dummy = 0.0)
      write(iounit) dummy,row
      return
      end
c-----------------------------------------------------------------------
      subroutine wrheader(unit,id,pgm,ncol,nrow,nz,
     & xo,dx,yo,dy,iproj,cm,bl)
      character id*56, pgm*8
      integer unit
      write(unit) id,pgm,ncol,nrow,nz,xo,dx,yo,dy,iproj,cm,bl
      return
      end
c-----------------------------------------------------------------------
      subroutine proj
c      This subroutine is called thru its 3 entry points:
c           1. projset(iproj,cmdeg,bldeg) - sets up both fwd and inverse prjts.
c                iproj = projection (see list below)
c                cm = central meridian in degrees, negative for west lon.
c                bl = base latitude in degrees, positive for north lat.
c           2. projfwd(xdeg,ydeg,xkm,ykm) - does forward projection.
c                xdeg,ydeg = lon,lat in decimal degrees (- for west lon).
c                xkm,ykm = projected coordinates in km.
c           3. projinv(xkm,ykm,xdeg,ydeg) - does inverse projection.
c
c      Projections available are:
c            1 = polyconic
c            2 = ellipsoidal transverse mercator
c            3 = mercator
c            4 = Lambert (standard U.S. parallels at 33N and 45N.
c            5 = Albers US (parallels at 29.5 and 45.5)
c            6 = Albers Alaska (         55.0 and 65.0)
c            7 = Albers Hawaii (          8.0 and 18.0)
c           12 = DNAG = spherical transverse mercator
c           14 = Lambert with standard parallels requested.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      common/projconsts/ ipr,cm,bl,xcm,ybl
      data sca/0.001/
      decdeg(deg,xmin)= deg + sign(xmin,deg)/60.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry projset(iproj,cmdeg,bldeg)
      ipr=iproj
      cm=cmdeg
      bl=bldeg
c      if(ipr.eq.0) call getproj(ipr,cm,bl)
      if(ipr.eq.0) return
      goto (10,20,90,40,50,60,70,90,5,5,5,22,5,42), ipr
  5   print *, 'PROJECTION NUMBER GIVEN NOT SUPPORTED...'
      stop
c.....Polyconic
 10   call polyset(cm,bl)
       return
c.....Ellipsoidal transverse mercator
 20   call etmset(cm,bl)
       return
c.....DNAG = spherical transverse mercator
 22   call stmset(cm,bl,0.926)
c      scale factor set for DNAG base.
       return
c.....Lambert conformal conic (ellipsoidal) with standard US parallels
 40   call lambset(cm,bl,33.,45.)
c      parallels set to standard for US
       return
c.....Lambert conformal conic (ellipsoidal) with parallels requested
 42   print '(/,1h$,a)', '     *Give southern standard parallel (d m): '
      read *, p1d,p1m
      p1=decdeg(p1d,p1m)
      print '(/,1h$,a)', '     *Give northern standard parallel (d m): '
      read *, p2d,p2m
      p2=decdeg(p2d,p2m)
      call lambset(cm,bl,p1,p2)
       return
c.....Albers conical equal-area (ellipsoidal) - US standard parallels
 50   call albset(cm,bl,29.5,45.5)
       return
c.....Albers conical equal-area (ellipsoidal) - Alaska standard parallels
 60   call albset(cm,bl,55.0,65.0)
       return
c.....Albers conical equal-area (ellipsoidal) - Hawaii standard parallels
 70   call albset(cm,bl,8.0,18.0)
       return
c.....All others
 90   call prjctl(bl,cm,xcm,ybl,cm,sca,ipr)
       return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry projfwd(xdeg,ydeg,xkm,ykm)
      if(ipr.eq.0) then
        xkm=xdeg
        ykm=ydeg
        return
      endif
      goto (110,120,190,140,150,150,150,190,105,105,105,122,105,140),
     &   ipr
 105  print *, 'PROJECTION NUMBER GIVEN NOT SUPPORTED...'
      stop
 110  call polyfwd(xdeg,ydeg,xkm,ykm)
       return
 120  call etmfwd(xdeg,ydeg,xkm,ykm)
       return
 122  call stmfwd(xdeg,ydeg,xkm,ykm)
       return
 140  call lambfwd(xdeg,ydeg,xkm,ykm)
       return
 150  call albfwd(xdeg,ydeg,xkm,ykm)
       return
 190  call prjctl(ydeg,xdeg,xkm,y,cm,sca,ipr)
      ykm=y-ybl
       return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry projinv(xkm,ykm,xdeg,ydeg)
      if(ipr.eq.0) then
        xdeg=xkm
        ydeg=ykm
        return
      endif
      goto (210,220,205,240,250,250,250,205,205,205,205,222,205,240),
     &    ipr
 205  print *, 'THIS INVERSE PROJECTION NOT SUPPORTED'
      stop
 210  call polyinv(xkm,ykm,xdeg,ydeg)
       return
 220  call etminv(xkm,ykm,xdeg,ydeg)
       return
 222  call stminv(xkm,ykm,xdeg,ydeg)
       return
 240  call lambinv(xkm,ykm,xdeg,ydeg)
       return
 250  call albinv(xkm,ykm,xdeg,ydeg)
       return
      end
c--------------------------------------------------------------------------
      subroutine polyconic
c
c      Ellipsoidal form of the polyconic projection,
c       based on the Clarke 1866 ellipsoid.
c      Single precision used ....
c      See J.P. Snyder, 1982, Map projections used by the USGS:
c       USGS Bulletin no. 1532, p. 129-131, p. 256-258.
c      Programmed by R.Simpson 5/16/84
c
c      Used by calling 1st entry point to set constants, then 2nd or 3rd
c       entry points to do forward or inverse projection respectively.
c      For example:
c                         .
c                         .
c               call polyset(cm,bl)
c                         .
c                         .
c               call polyfwd(xdeg,ydeg,xkm,ykm)
c                         .
c                         .
c         (or)  call polyinv(xkm,ykm,xdeg,ydeg)
c                         .
c                         .
c
c      where: cm = central meridian in decimal degrees (-for west)
c             bl = base latidude in decimal degrees (+ for north)
c             xdeg = longitude in decimal degrees
c             ydeg = latitude in decimal degrees
c             xkm,ykm = projected coordinates in km.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real n,mfcn,mpfcn,m,m0,mn,mpn,ma,lam,lam0
      parameter (maxiter=10, tol=1.0E-6)
      parameter (a=6378.2064)
      parameter (e=0.0822719, esq=0.676866E-2)
      parameter (deg2rad=0.17453293E-1, rad2deg=57.295778)
      common/pconsts/ m0,c0,c2,c4,c6,lam0,phi0
c     functions...
c      (eqn 3-21, p. 20)....
      mfcn(phi) =
     &    a*(c0*phi-c2*sin(2.*phi)+c4*sin(4.*phi)-c6*sin(6.*phi))
c      (eqn 15-17, p. 130)....
      mpfcn(phi) =
     &    c0-2.0*c2*cos(2.0*phi)+4.0*c4*cos(4.0*phi)-6.0*c6*cos(6.0*phi)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry polyset(cm,bl)
c      cm= central meridian in degrees (- for w)
c      bl= base latitude in degrees
      lam0=cm*deg2rad
      phi0=bl*deg2rad
      c0=1.-e**2/4.-3.*e**4/64.-5.*e**6/256.
      c2=3.*e**2/8.+3.*e**4/32.+45.*e**6/1024.
      c4=15.*e**4/256.+45.*e**6/1024.
      c6=35.*e**6/3072.
      m0=mfcn(phi0)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry polyfwd(xdeg,ydeg,xkm,ykm)
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
c      xkm,ykm = projected coordinates in km.
      lam=xdeg*deg2rad
      phi=ydeg*deg2rad
      n=a/sqrt(1.-esq*sin(phi)**2)
      cotphi=1.0/tan(phi)
      m=mfcn(phi)
      if(phi.ne.0) then
        bige=(lam-lam0)*sin(phi)
        xkm=n*sin(bige)*cotphi
        ykm=m-m0+n*(1.0-cos(bige))*cotphi
      else
        xkm=a*(lam-lam0)
        ykm=-m0
      endif
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry polyinv(xkm,ykm,xdeg,ydeg)
c      xkm,ykm = projected coordinates in km.
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
      if(ykm.ne.-m0) then
        biga=(m0+ykm)/a
        bigb=(xkm/a)**2 + biga**2
        phin=biga
        do 40 i=1,maxiter
         bigc=sqrt(1.-esq*sin(phin)**2)*tan(phin)
         mn=mfcn(phin)
         mpn=mpfcn(phin)
         ma=mn/a
         phinp1=phin - (biga*(bigc*ma+1.)-ma-0.5*(ma**2+bigb)*bigc)/
     &   (esq*sin(2.*phin)*(ma**2+bigb-2.*biga*ma)/(4.*bigc)
     &    +(biga-ma)*(bigc*mpn-2./sin(2.*phin))-mpn)
         if(abs(phinp1-phin).le.tol) go to 60
         phin=phinp1
 40     continue
c.....  Upon normal exit from DO loop print...
        print *, 'WARNING...no convergence to tolerance'
        print *, '  after maxiterations =', maxiter
        print *, '    for xkm,ykm =', xkm,ykm
 60     phi=phinp1
        lam=asin(xkm*bigc/a)/sin(phi) + lam0
      else
        phi=0.0
        lam=xkm/a + lam0
      endif
      xdeg=lam*rad2deg
      ydeg=phi*rad2deg
      return
      end
c--------------------------------------------------------------------------
      subroutine etm
c      Ellipsoidal form of transverse mercator,
c       based on the Clarke 1866 ellipsoid.
c      Single precision used ....
c      See J.P. Snyder, 1982, Map projections used by the USGS:
c       USGS Bulletin no. 1532, p. 67-69, p.232-235.
c      Formulas not accurate for longitudes greater than 4 degrees from
c       the chosen central meridian (Snyder, p. 67).
c      Programmed by R.Simpson 12/21/83.
c      Used by calling 1st entry point to set constants, then 2nd or 3rd
c       entry points to do forward or inverse projection respectively.
c      For example:
c                         .
c                         .
c               call etmset(cm,bl)
c                         .
c                         .
c               call etmfwd(xdeg,ydeg,xkm,ykm)
c                         .
c                         .
c         (or)  call etminv(xkm,ykm,xdeg,ydeg)
c                         .
c                         .
c
c      where: cm = central meridian in decimal degrees (-for west)
c             bl = base latidude in decimal degrees (+ for north)
c             xdeg = longitude in decimal degrees
c             ydeg = latitude in decimal degrees
c             xkm,ykm = projected coordinates in km.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real k0,n,n1,mu,mfcn,m,m0,lam,lam0
      parameter (a=6378.2064)
      parameter (e=0.0822719, esq=0.676866E-2)
      parameter (k0=0.9996)
      parameter (deg2rad=0.17453293E-1, rad2deg=57.295778)
      common/etmconsts/ iwarn,epsq,m0,c0,c2,c4,c6,lam0,phi0,cmdeg
c     functions....
      mfcn(phi) =
     &    a*(c0*phi-c2*sin(2.*phi)+c4*sin(4.*phi)-c6*sin(6.*phi))
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry etmset(cm,bl)
c      cm= central meridian in degrees (- for w)
c      bl= base latitude in degrees
      iwarn=0
      cmdeg=cm
      lam0=cm*deg2rad
      phi0=bl*deg2rad
      c0=1.-e**2/4.-3.*e**4/64.-5.*e**6/256.
      c2=3.*e**2/8.+3.*e**4/32.+45.*e**6/1024.
      c4=15.*e**4/256.+45.*e**6/1024.
      c6=35.*e**6/3072.
      m0=mfcn(phi0)
      epsq=esq/(1.0-esq)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry etmfwd(xdeg,ydeg,xkm,ykm)
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
c      xkm,ykm = projected coordinates in km.
      if(abs(xdeg-cmdeg).gt.4.0.and.iwarn.eq.0) then
       print *, 'WARNING...accuracy drops for lons greater than'
       print *, '    4 degrees from the central meridian...'
       iwarn=iwarn+1
      endif
      lam=xdeg*deg2rad
      phi=ydeg*deg2rad
      n=a/sqrt(1.-esq*sin(phi)**2)
      tnphi=tan(phi)
      csphi=cos(phi)
      t=tnphi**2
      c=epsq*csphi**2
      biga=cos(phi)*(lam-lam0)
      m=mfcn(phi)
      xkm=k0*n*(biga+(1.-t+c)*biga**3/6. +
     &     (5.-18.*t+t**2+72.*c-58.*epsq)*biga**5/120.)
      ykm=k0*(m-m0+n*tnphi*(.5*biga**2 +
     &     (5.-t+9.*c+4.*c**2)*biga**4/24. +
     &     (61.-58.*t+t**2+600.*c-330.*epsq)*biga**6/720.))
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry etminv(xkm,ykm,xdeg,ydeg)
c      xkm,ykm = projected coordinates in km.
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
      m=m0 + ykm/k0
      e1=(1.-sqrt(1.-esq))/(1.+sqrt(1.-esq))
      mu=m/(a*(1.-esq/4.-3.*esq**2/64.-5.*esq**3/256.))
      phi1=mu+(3.*e1/2.-27.*e1**3/32.)*sin(2.*mu) +
     &      (21.*e1**2/16.-55.*e1**4/32.)*sin(4.*mu) +
     &      (151.*e1**3/96.)*sin(6.*mu)
      c1=epsq*cos(phi1)**2
      t1=tan(phi1)**2
      n1=a/sqrt(1.-esq*sin(phi1)**2)
      r1=a*(1-esq)/((1-esq*sin(phi1)**2)*sqrt(1-esq*sin(phi1)**2))
      d=xkm/(n1*k0)
      phi=phi1-(n1*tan(phi1)/r1)*
     &     (d**2/2.-(5.+3.*t1+10.*c1-4.*c1**2-9.*epsq)*d**4/24.
     &     + (61.+90*t1+298*c1+45*t1**2-252.*epsq-3.*c1**2)*d**6/720.)
      lam=lam0+(d-(1.+2.*t1+c1)*d**3/6.
     &     +(5.-2.*c1+28.*t1-3.*c1**2+8.*epsq+24.*t1**2)*d**5/120.)
     &     /cos(phi1)
      xdeg=lam*rad2deg
      ydeg=phi*rad2deg
      return
      end
c--------------------------------------------------------------------------
      subroutine stm
c      Spherical form of transverse mercator.
c      Single precision used ....
c      See J.P. Snyder, 1982, Map projections used by the USGS:
c       USGS Bulletin no. 1532, p. 64-67, p.231-232.
c      Programmed by R.Simpson 1/10/84.
c      Used by calling 1st entry point to set constants, then 2nd or 3rd
c       entry points to do forward or inverse projection respectively.
c      For example:
c                         .
c                         .
c               call stmset(cm,bl,scalefac)
c                         .
c                         .
c               call stmfwd(xdeg,ydeg,xkm,ykm)
c                         .
c                         .
c         (or)  call stminv(xkm,ykm,xdeg,ydeg)
c                         .
c                         .
c
c      where: cm = central meridian in decimal degrees (-for west)
c             bl = base latidude in decimal degrees (+ for north)
c             scalefac = central scale factor
c             xdeg = longitude in decimal degrees
c             ydeg = latitude in decimal degrees
c             xkm,ykm = projected coordinates in km.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real k0,phi,phi0,lam,lam0
      double precision d,e,dphi,dlam
      parameter (r=6371.204)
c      parameter (r=1.)
      parameter (deg2rad=0.17453293E-1, rad2deg=57.295778)
      common/stmconsts/ cmdeg,lam0,phi0,k0
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry stmset(cm,bl,scalefac)
c      scalefac = central scale factor
c      cm= central meridian in degrees (- for w)
c      bl= base latitude in degrees
      cmdeg=cm
      k0=scalefac
      lam0=cm*deg2rad
      phi0=bl*deg2rad
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry stmfwd(xdeg,ydeg,xkm,ykm)
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
c      xkm,ykm = projected coordinates in km.
      lam=xdeg*deg2rad
      phi=ydeg*deg2rad
      b = cos(phi)*sin(lam-lam0)
      xkm=0.5*r*k0*log((1.+b)/(1.-b))
      ykm=r*k0*(atan(tan(phi)/cos(lam-lam0)) - phi0)
c      k=k0/sqrt(1.-b**2)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry stminv(xkm,ykm,xdeg,ydeg)
c      xkm,ykm = projected coordinates in km.
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
      d=dble(ykm)/dble(r*k0) + dble(phi0)
      e=dble(xkm)/dble(r*k0)
      dphi=asin(sin(d)/cosh(e))
      dlam=dble(lam0) + atan(sinh(e)/cos(d))
      xdeg=dlam*rad2deg
      ydeg=dphi*rad2deg
      return
      end
c--------------------------------------------------------------------------
      subroutine lambert
c
c      Lambert conformal conic projection
c            based on the Clarke 1866 ellipsoid.
c      Single precision used .... max error probably < 5 meters.
c      See J.P. Snyder, 1982, Map projections used by the USGS:
c            USGS Bulletin no. 1532, p. 101-109, p.250-252.
c      Normal standard parallels for US are p1=33N and p2=45N.
c       See Snyder Table 8, p. 60-62 for state plane coordinate systems.
c      Programmed by R.Simpson 8/4/1983.
c
c      Used by calling 1st entry point to set constants, then 2nd or 3rd
c       entry points to do forward or inverse projection respectively.
c      For example:
c                         .
c                         .
c               call lambset(-96., 30.,p1,p2)
c                         .
c                         .
c               call lambfwd(-123., 45.,xkm,ykm)
c                         .
c                         .
c         (or)  call lambinv(xkm,ykm,dlon,dlat)
c                         .
c                         .
c      This will give x,y in km... (or) call to inverse will give back
c       dlon=-123., dlat=45.  p1 and p2 are the chosen parallels.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      parameter (a=6378.2064)
      parameter (e=0.0822719, esq=0.676866E-2)
      parameter (deg2rad=0.17453293E-1, rad2deg=57.295778, pi=3.1415927)
      real m1,m2,mfcn,n,lam,lam0
      common/lambconsts/ n,bigf,lam0,phi0,rho0
c     functions....
c.....(eqn 12-15, p. 97)...
      mfcn(cphi,sphi) = cphi/sqrt(1.0-esq*sphi**2)
c.....(eqn 13-9, p. 107)...
      tfcn(phi)=
     &   tan(0.25*pi-0.5*phi)/((1.-e*sin(phi))/(1.+e*sin(phi)))**(0.5*e)
c.....(eqn 13-7, p 107)...
      rhofcn(t) = a*bigf*(t)**n
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry lambset(cm,bl,p1,p2)
c      cm= central meridian in degrees (- for w)
c      bl= base latitude in degrees
c      p1,p2= parallels chosen (33 and 45 for U.S. State Maps)
c              (p1=south parallel, p2=north parallel)
      lam0=cm*deg2rad
      phi0=bl*deg2rad
      phi1=p1*deg2rad
      phi2=p2*deg2rad
      csphi1=cos(phi1)
      snphi1=sin(phi1)
      csphi2=cos(phi2)
      snphi2=sin(phi2)
      m1=mfcn(csphi1,snphi1)
      m2=mfcn(csphi2,snphi2)
      t0=tfcn(phi0)
      t1=tfcn(phi1)
      t2=tfcn(phi2)
      n=(log(m1)-log(m2))/(log(t1)-log(t2))
      bigf=m1/(n*(t1)**n)
      rho0=rhofcn(t0)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry lambfwd(xdeg,ydeg,xkm,ykm)
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
c      xkm,ykm = projected coordinates in km.
      lam=xdeg*deg2rad
      phi=ydeg*deg2rad
      t=tfcn(phi)
      rho=rhofcn(t)
      th=n*(lam-lam0)
      xkm=rho*sin(th)
      ykm=rho0-rho*cos(th)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry lambinv(xkm,ykm,xdeg,ydeg)
c      xkm,ykm = projected coordinates in km.
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
      rho=sign(sqrt(xkm**2+(rho0-ykm)**2),n)
      th=atan2(xkm*sign(1.,n),(rho0-ykm)*sign(1.,n))
      t=(rho/(a*bigf))**(1./n)
      chi=0.5*pi - 2.0*atan(t)
      ydeg=(chi+(esq/2+5*esq**2/24+esq**3/12)*sin(2*chi)
     &    +(7*esq**2/48+29*esq**3/240)*sin(4*chi)
     &    +(7*esq**3/160)*sin(6*chi))*rad2deg
      xdeg=(th/n+lam0)*rad2deg
      return
      end
c--------------------------------------------------------------------------
      subroutine alber
c      Albers equal-area conic projection
c       based on the Clarke 1866 ellipsoid.
c      Single precision used .... max error probably < 5 meters.
c      See J.P. Snyder, 1982, Map projections used by the USGS:
c       USGS Bulletin no. 1532, p. 93-99, p.245-248.
c      Set up to take standard parallels for US, Alaska, or Hawaii.
c      Programmed by R.Simpson 8/4/1983.
c      Used by calling 1st entry point to set constants, then 2nd or 3rd
c       entry points to do forward or inverse projection respectively.
c      For example:
c                         .
c                         .
c               call albset(-96., 30., p1,p2)
c                         .
c                         .
c               call albfwd(-123., 45., xkm,ykm)
c                         .
c                         .
c         (or)  call albinv(xkm,ykm,dlon,dlat)
c                         .
c                         .
c      This will give x,y in km... call to inverse will give back
c       dlon=-123., dlat=45.  p1 and p2 are the chosen standard parallels.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      parameter (a=6378.2064)
      parameter (e=0.0822719, esq=0.676866E-2)
      parameter (deg2rad=0.17453293E-1, rad2deg=57.295778)
      real m1,m2,mfcn,n,lam,lam0
      common/albconsts/ c,n,lam0,phi0,rho0
c     functions....
      mfcn(cphi,sphi) = cphi/sqrt(1.0-esq*sphi**2)
      qfcn(sphi) = (1.-esq)*(sphi/(1.-esq*sphi**2) - (0.5/e)*
     &                 log((1.-e*sphi)/(1.+e*sphi)))
      rhofcn(q) = a*sqrt(c-n*q)/n
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry albset(cm,bl,p1,p2)
c      cm= central meridian in degrees (- for w)
c      bl= base latitude in degrees
c      p1= southern standard parallel
c      p2= northern standard parallel
c               p1,p1= 29.5, 45.5 for U.S.
c                    = 55.0, 65.0 for Alaska
c                    =  8.0, 18.0 for Hawaii
      lam0=cm*deg2rad
      phi0=bl*deg2rad
      phi1=p1*deg2rad
      phi2=p2*deg2rad
      snphi0=sin(phi0)
      csphi1=cos(phi1)
      snphi1=sin(phi1)
      csphi2=cos(phi2)
      snphi2=sin(phi2)
      m1=mfcn(csphi1,snphi1)
      m2=mfcn(csphi2,snphi2)
      q1=qfcn(snphi1)
      q2=qfcn(snphi2)
      q0=qfcn(snphi0)
      n=(m1**2-m2**2)/(q2-q1)
      c=m1**2+n*q1
      rho0=rhofcn(q0)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry albfwd(xdeg,ydeg,xkm,ykm)
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
c      xkm,ykm = projected coordinates in km.
      lam=xdeg*deg2rad
      phi=ydeg*deg2rad
      snphi=sin(phi)
      csphi=cos(phi)
      q=qfcn(snphi)
      rho=rhofcn(q)
      th=n*(lam-lam0)
      xkm=rho*sin(th)
      ykm=rho0-rho*cos(th)
      return
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry albinv(xkm,ykm,xdeg,ydeg)
c      xkm,ykm = projected coordinates in km.
c      xdeg,ydeg = lon (- for w) and lat in decimal degrees.
      rho=sqrt(xkm**2+(rho0-ykm)**2)
      th=atan2(xkm*sign(1.,n),(rho0-ykm)*sign(1.,n))
      q=(c-(rho*n/a)**2)/n
      beta=asin(q/(1-((1-esq)/(2*e))*log((1-e)/(1+e))))
      ydeg=(beta+(esq/3+31*esq**2/180+517*esq**3/5040)*sin(2*beta)
     &    +(23*esq**2/360+251*esq**3/3780)*sin(4*beta)
     &    +(761*esq**3/45360)*sin(6*beta))*rad2deg
      xdeg=(lam0+th/n)*rad2deg
      return
      end
c--------------------------------------------------------------------------
      subroutine prjctl(ylat,xlon,x,y,cm,sfac,iprojt)
c
c     Programmed by Dick Godson, USGS Denver...
c
      double precision yl,xl,x1,y1
      yl=dble(ylat)*.0174532925199433d0
      xl=dble(abs(cm)-abs(xlon))*.0174532925199433d0
      go to (10,20,30,40,50,50,50,60),iprojt
   10 call poly(yl,xl,x1,y1)
      go to 90
   20 call utmfwd(yl,xl,x1,y1)
      go to 90
   30 call merctr(yl,xl,x1,y1)
      go to 90
   40 call lambrt(yl,xl,x1,y1)
      go to 90
   50 call albers(yl,xl,x1,y1)
      go to 90
   60 call pstereo(yl,xl,x1,y1)
      go to 90
   90 x=sngl(x1)*sfac
      y=sngl(y1)*sfac
      return
      end
c--------------------------------------------------------------------------
      subroutine poly(phi,dlamb,x,y)
      double precision phi,dlamb,xc,x,y,a,b,t,q,phi2,c1,c2
      a(xc)=
     &   6.378206402718907d 06 +xc*(
     &  -3.167517353503576d 06 +xc*(
     &   2.478805037574243d 05 +xc*(
     &  -3.530710396439220d 03 +xc*(
     &  -6.565371848240127d 02 +xc*(
     &   6.822539551727124d 01 +xc*(
     &  -2.888860980506611d 00 ))))))
      b(xc)=
     &   3.189103200618349d 06 +xc*(
     &  -2.115275857345996d 06 +xc*(
     &   4.144758431728325d 05 +xc*(
     &  -3.625295427368928d 04 +xc*(
     &   1.322429746943889d 03 +xc*(
     &   4.427183005616853d 01 +xc*(
     &  -8.630738701658902d 00 +xc*(
     &   4.279183401413811d-01 )))))))
      q(xc)=
     &   6.335034386662446d 06 +xc*(
     &   2.144094496614083d 04 +xc*(
     &  -4.182226973512467d 03 +xc*(
     &   3.609995316780635d 02 +xc*(
     &  -1.346978283534684d 01 ))))
      t(xc)=
     &   9.999999957157490d-01+xc*(
     &  -1.666665796975878d-01+xc*(
     &   8.333050613721043d-03+xc*(
     &  -1.980904608528695d-04+xc*
     &   2.605165638554101d-06 )))
      phi2=phi**2
      c1=dlamb*phi
      c2=(c1*t(phi2))**2
      x=dlamb*a(phi2)*t(c2)
      y=dlamb*c1   *b(phi2)*t(.25d0*c2)**2 + phi*q(phi2)
      return
      end
c--------------------------------------------------------------------------
      subroutine utmfwd(phi,dlam,x,y)
c
c  SUBROUTINE UTMFWD
c
c  PURPOSE..
c    TO DETERMINE BASIC UTM COORDIATES FROM LATITUDE
c    AND LONGITUDE INPUT.  CLARK 1866 SPHEROID CONSTANTS.
c
c  USAGE..
c    CALL UTMFWD(PHI,DLAM,X,Y,IER)
c
c  PARAMETER DESCRIPTION..
c    PHI    - INPUT LATITUDE IN DOUBLE PRECISION RADIANS.
c             DABS(PHI).LT.1.3963 (80 DEGREES).
c    DLAM   - INPUT LONGITUDE IN DOUBLE PRECISION RADIANS.
c             DABS(DLAM).LT.0.06109 (3.5 DEGREES).
c             NOTE.. DLAM MEASURED FROM CENTRAL MERIDIAN.
c    X      - PROJECTED COORDINATE IN METERS FROM THE
c             CENTRAL MERIDIAN (DOUBLE PRECISION).
c    Y      - PROJECTED COORDINATE IN METERS FROM THE
c             EQUATOR (DOUBLE PRECISION).
c    IER    - ERROR RETURN CODE.
c             IER=0 SUCCESSFUL COMPLETION.
c             IER.NE.0 ERROR IN PHI OR DLAM INPUT.
c
c  REMARKS..
c    THIS ROUTINE IS ACCURATE TO 1 METER.  ABSOLUTE
c    ERROR LESS THAN .5 METERS.
c    NEAREST CENTRAL MERIDIAN MAY BE DETERMINED EXTERNALLY
c    BY..
c    CN=DSIGN(10800D0,LON)+AINT(SNGL(LON/21600D0))*21600D0
c      WHERE LON IS LONGITUDE IN DOUBLE PRECISION SECONDS.
c    FALSE EASTING (500000 METERS) AND NORTHING (0 FOR
c    NORTHERN LATITUDES, 10000000 METERS FOR SOUTHERN
c    LATITUDES) FOR TRUE UTM COORDINATES MUST BE COMPUTED
c    EXTERNALLY.
c
c  SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED..
c    NONE
c
c  DEVELOPED AND CODED BY..
c    GERALD IAN EVENDEN
c    U. S. GEOLOGICAL SURVEY
c
c
          double precision phi,dlam,x,y,dl2,p
c
c  CHECK INPUT PRAMETERS..
          p=phi*phi
          dl2=dlam*dlam
c         if (p .le.1.94965360d0.and.dl2.le.3.7319881d-3)
c             go to 10
c         ier=-1
c         return
c
10        ier=0
c
c  DETERMINE PROJECTED COORDINATES..
          y=phi*(
     &          6332500.47d0+p*(21431.67d0+p*(-4179.269d0+p*(
     &          359.981d0-p*13.267d0)))
     &      -dl2*(-3187827d0+p*(2114440d0+p*(-414363.5d0+p*(
     &          36344.6d0-p*1420.3d0)))
     &       -dl2*(1334935d0+p*(-2356027d0+p*(1371758d0-p*
     &          267852d0)))))
          x=dlam*(6375655.2d0+p*(-3166253.6d0+p*(247800.26d0+p*(
     &          -3569.65d0+p*(-617.35d0+p*50.89d0))))
     &       -dl2*(-1069479d0+p*(2661562d0+p*(-1785956d0+p*(
     &         485045d0-p*52022d0)))))
c
          return
c
          end
c-------------------------------------------------------------------------
          subroutine merctr(lat,long,x,y)
          double precision lat,long,x,y,a,halfpi,b2da2,z
          data a/6378206.d0/,halfpi/1.57079632679489d0/,
     &      b2da2/9.932315290818186d-1/
          x=a*long
c
c         compute z/2
c
          z=0.5d0*(halfpi-datan(b2da2*dsin(lat)/dcos(lat)))
          y=a*dlog(dcos(z)/dsin(z))
          return
          end
c--------------------------------------------------------------------------
          subroutine lambrt(lat,long,x,y)
          double precision lat,long,x,y,z,k,l,r,b2da2,halfpi,theta
          data k/1.245265545827296d7/,l/6.304998344052458d-1/,
     &      b2da2/9.932315290818186d-1/,halfpi/1.57079632679489662d0/
c
c         compute z/2
c
          z=0.5d0*(halfpi-datan(b2da2*dsin(lat)/dcos(lat)))
c
c         compute r
c
          r=k*(dsin(z)/dcos(z))**l
c
c         compute coordinates
c
          theta=l*long
          x=r*dsin(theta)
          y=-r*dcos(theta)
          return
          end
c--------------------------------------------------------------------------
          subroutine albers(ylat,xlon,x,y)
          common/albcom/n,rho1sq,sinbt1,twoc2n
          double precision ylat,xlon,x,y,n,rho1sq,sinbt1,twoc2n,
     &      nus,nals,nhaw,rho295,rho55,rho8,tcnus,tcnals,tcnhaw,
     &      a1,b,c1,d,e1,f1,g,h,theta,rho,sinbet,sinphi,s2,
     &      sbt295,sbt55,sbt8
          data a1/9.954804334645587d-1/,b/4.492024607745888d-3/,
     &      c1/2.736435989866449d-5/,d/1.763992166249299d-7/,
     &      e1/1.160814577272288d-9/,f1/7.714265487727804d-12/,
     &      g/5.154557173568170d-14/,h/3.455700205911349d-16/,
     &      nus/6.02903493787094d-1/,nals/8.627447947235633d-1/,
     &      nhaw/2.241096394314637d-1/,rho295/8.49196923967458d13/,
     &      rho55/1.806308673895081d13/,rho8/7.943986660586285d14/,
     &      tcnus/1.346470921892769d14/,tcnals/9.409410848453636d13/,
     &      tcnhaw/3.622298555079059d14/,sbt295/4.907351753179611d-1/,
     &      sbt55/8.1792905450587868d-1/,sbt8/1.385562096187223d-1/
          sinphi=dsin(ylat)
          s2=sinphi*sinphi
          sinbet=sinphi*(a1+s2*(b+s2*(c1+s2*(d+s2*
     &      (e1+s2*(f1+s2*(g+s2*h)))))))
          rho=dsqrt(rho1sq+twoc2n*(sinbt1-sinbet))
          theta=n*xlon
          x=rho*dsin(theta)
          y=-rho*dcos(theta)
          return
c
c set up constants
          entry setalb(iproj)
          if(iproj-6)500,510,520
500       n=nus
          rho1sq=rho295
          sinbt1=sbt295
          twoc2n=tcnus
          return
510       n=nals
          rho1sq=rho55
          sinbt1=sbt55
          twoc2n=tcnals
          return
520       n=nhaw
          rho1sq=rho8
          sinbt1=sbt8
          twoc2n=tcnhaw
          return
          end
c
c----------------------------------------------------------
c
c  SUBROUTINE UTMINV
c
c  PURPOSE..
c    TO DETERMINE LATITUDE AND LONGITUDE FROM UTM
c    COORDINATES.  CLARK 1866 SPHEROID CONSTANTS.
c
c  USAGE..
c    CALL UTMINV(PHI,DLAM,X,Y,IER)
c
c  DESCRIPTION OF PARAMETERS..
c    PHI    - OUTPUT LATITUDE IN DOUBLE PRECISION RADIANS.
c    DLAM   - OUTPUT LONGITUDE IN DOUBLE PRECISION RADIANS
c             FROM CENTRAL MERIDIAN.
c    X      - DOUBLE PRECISION UTM METERS FROM THE CENTRAL
c             MERIDIAN.  FALSE EASTING MUST BE
c             REMOVE EXTERNALLY.
c             DABS(X).LE.390000 METERS.
c    Y      - DOUBLE PRECISION UTM METERS FROM THE EQUATOR.
c             FALSE NORTHING (SOUTHERN LATITUDES) MUST BE
c             REMOVED EXTERNALLY.
c             DABS(Y).LE.8900000 METERS.
c    IER    - ERROR RETURN CODE.
c             IER=0 SUCCESSFUL COMPLETION.
c             IER.NE.0 ERROR IN X OR Y INPUT.
c
c  REMARKS..
c    ABSOLUTE ERROR OF PHI AND DLAM LESS THAN 10**-8
c    RADIANS.
c
c  SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED..
c    NONE
c
c  DEVELOPED AND CODED BY..
c    GERALD IAN EVENDEN
c    U. S. GEOLOGICAL SURVEY
c
c
c--------------------------------------------------------------------------
          subroutine utminv(phi,dlam,x,y)
c
          double precision dlam,phi,x,y,x1,x2,y1,y2,
     &       ph1,ph2,dt
c
          equivalence (dt,x1,y1),(x2,y2)
c
          y1=y*1d-7
          y2=y1*y1
c         if (y2.le..7921d0.and.dabs(x).le.390000d0) go to 10
c         ier=-1
c         return
c
c  COMPUTE INVERSE PROJECTION
c 10       ier=0
          ph1=y1*(1.579155027d0+y2*(-.0133273056d0+y2*(
     &         .006813653d0+y2*(-.0018179637d0+y2*(
     &         .000352375d0-y2*.0000494195d0)))))
          ph2=ph1*ph1
          x1=x*1d-7/(.637565512d0+ph2*(-.318782754d0+ph2*(
     &         .0265652212d0+ph2*(-.0008854896d0+ph2*(
     &         .00001579516d0-ph2*.00000016787d0)))))
          x2=x1*x1
          dlam=x1*(.99999998d0+ph2*(-.00338386d0+ph2*(
     &         .001120504d0+ph2*(-.000143997d0+ph2*
     &         .0000081386d0)))
     &      -x2*(.167804d0+ph2*(.162647d0+ph2*(
     &         -.054578d0+ph2*(.0073115d0-ph2*.00045341d0)))
     &      -x2*(.04297d0+ph2*(.1402d0-ph2*.03291d0))))
          dt=-ph1*x2*(.4999998d0+ph2*(-.08671332d0+ph2*(
     &         .00584112d0+ph2*(-.000440204d0+ph2*
     &         .000031613d0)))
     &      -x2*(.20928d0+ph2*(-.0028705d0-ph2*.0081252d0)))
          phi=ph1+dt*(1.0068147d0+ph2*(-.51022062d0+ph2*(
     &         .047620417d0+ph2*(-.003100987d0+ph2*(
     &         .000282232d0-ph2*.0000185462d0))))
     &      +dt*(ph1*(-.51358d0+ph2*(.354256d0+ph2*(
     &         -.0770704d0+ph2*.0073699d0)))))
c
          return
c
          end
c--------------------------------------------------------------------------
          subroutine pstereo(ylat,xlon,x,y)
c
c         polar stereographic projection
c
          double precision ylat,xlon,x,y,sinc,sine,tane1,rad45,e,ee,
     &      tanz,polarc,a,dlat2,tane2,dlat
c          data e/.08182d0/,r/6378160d0/,rad45/.7854d0/,
          data e/.08182d0/,rad45/.7854d0/,
     &      ee/.04091d0/,polarc/12713645.1d0/
          dlat=dabs(ylat)
c         polarc=2.d0*r/dsqrt(1.d0-e*e)*((1.d0-e)/(1.d0+e))**ee
          sinc=e*dsin(dlat)
          sine=(1.d0+sinc)/(1.d0-sinc)
          dlat2=dlat*.5d0
          tane2=rad45+dlat2
          tane1=cos(tane2)/sin(tane2)
          tanz=tane1*(sine**ee)
          a=tanz*polarc
          y=-(a*dcos(xlon))
          x=a*dsin(xlon)
          if(ylat.lt.0.d0) y=-y
c         scale=a*dsqrt(1.d0-e*e*dsin(dlat)**2)/dcos(dlat)/r
          return
          end
c--------------------------------------------------------------------------
      function   decdeg(deg,min)
c     Converts degrees and minutes to decimal degrees...
      real decdeg, deg, min
      decdeg = deg + sign(min,deg)/60.
      return
      end     
tes to decimal degrees...

