c   GENPROJ.FOR
c     a general forward and inverse projection program
c     that works on binary xyz or post files
c     a xyz file can be created from a post file
c     important note**** input coordinates for forward
c     projection of west longtitudes must be negative
      character*1 ans
      character*1 type
      character*8 id
      CHARACTER*56 NAME,NAME1
      dimension z(24),cmbl(4)
      real*8 id8
c
        call askin
      write(*,*) 'GENERAL PROJECTION PROGRAM'
      write(*,*) 'west longtitudes must be negative'
C      WRITE(*,1000)
C 1000 FORMAT(' ENTER INPUT FILE NAME: '$)
C      READ(*,1100) NAME
C 1100 FORMAT(A)
C      OPEN(8,FILE=NAME,STATUS='OLD',form='unformatted')
        name=' '
400     call askc('Enter input file name',name,ierr)
        if(ierr.lt.0) stop
        call xopen(8,name,'old','read',itype,ierr)
        if(ierr.ne.0) go to 400
C      WRITE(*,1150)
C 1150 FORMAT(' ENTER OUTPUT FILE NAME: '$)
C      READ(*,1100) NAME1
C      OPEN(9,FILE=NAME1,STATUS='UNKNOWN',FORM='UNFORMATTED')
C      write(*,*) 'is input file xyz or pos? '
C      read(*,1100) type
        name1=' '
401     call askc('Enter output file name',name1,ierr)
        if(ierr.eq.-1) go to 401
        if(ierr.eq.-2) then
          call xclose(8,'delete')
          go to 400
        endif
        if(itype.ne.-1) then
          call xopen(9,name1,'new','write',itype,ierr)
        else
          open(9,file=name1,status='new',form='formatted',iostat=ierr)
          if(ierr.eq.6415) then
            ians=0
            call aski4l(' File exists -- OK to overwrite?',ians,ierr)
            if(ians.eq.1) then
              open(9,file=name1,status='unknown',form='formatted',
     1             iostat=ierr)
            endif
          endif
        endif
c        if(ierr.ne.0) go to 401
C      if(type.eq.'p'.or.type.eq.'P') then
C      ianom=0
C      write(*,*) 'convert post to xyz?(y/n) '
C      read(*,1100) ans
C      if(ans.eq.'y'.or.ans.eq.'Y') then
C      write(*,*) 'enter z number '
C      read(*,*) ianom
C      endif
C      endif
C      write(*,*) 'enter cm,baslat,proj no. & 1(forw) or 2(inv)'
C      read(*,*) cm,baslat,iproj,im
        call mphnam
        cmbl(1)=0.
        cmbl(2)=0.
        cmbl(3)=0.
        cmbl(4)=0.
402     call askf4a('Enter cm,baslat,proj no. & 1(forw) or 2(inv)',
     1              cmbl,4,ierr)
        if(ierr.eq.-1) go to 402
        if(ierr.eq.-2) then
          call xclose(9,'delete')
          go to 401
        endif
        cm=cmbl(1)
        baslat=cmbl(2)
        iproj=cmbl(3)
        im=cmbl(4)
        if(iproj.lt.1.or.iproj.gt.9) go to 402
        if(im.lt.1.or.im.gt.2) go to 402
      call projst(iproj,cm,baslat)
C      if(type.eq.'x'.or.type.eq.'X') then
        if(itype.eq.0) then
c
c     xyz file
c
      go to (10,20),im
C   10 READ(8,END=70) x,y,zz
   10 call xioxyz('r',8,z,3,ierr)
      if(ierr.ne.0) go to 90
      call projfw(z(1),z(2),xkm,ykm)
      z(1)=xkm
      z(2)=ykm
C      WRITE(9) xkm,ykm,zz
      call xioxyz('w',9,z,3,ierr)
      if(ierr.ne.0) stop 'error writing xyz file'
      go to 10
C   20 READ(8,END=70) x,y,zz
   20 call xioxyz('r',8,z,3,ierr)
      if(ierr.ne.0) go to 90
      x=z(1)
      y=z(2)
      call projin(x,y,xdeg,ydeg)
      z(1)=xdeg
      z(2)=ydeg
C      write(9) xdeg,ydeg,zz
      call xioxyz('w',9,z,3,ierr)
      if(ierr.ne.0) stop 'error writing xyz file'
      go to 20
      else if(itype.gt.0) then
c
c     post file
c
C      nword=lenrec(8)
C      if(ianom.eq.0) then
      go to (30,40),im
C   30 call io(z,nword,8,1,*70)
   30 call xiopst('r',8,id,z,itype,ierr)
      if(ierr.ne.0) go to 90
C      call projfw(z(3),z(4),xkm,ykm)
C      z(3)=xkm
C      z(4)=ykm
      call projfw(z(1),z(2),xkm,ykm)
      z(1)=xkm
      z(2)=ykm
C      call io(z,nword,9,2,*70)
      call xiopst('w',9,id,z,itype,ierr)
      if(ierr.ne.0) stop 'error writing post file'
      go to 30
C   40 call io(z,nword,8,1,*70)
   40 call xiopst('r',8,id,z,itype,ierr)
      if(ierr.ne.0) go to 90
C      call projin(z(3),z(4),xdeg,ydeg)
C      z(3)=xdeg
C      z(4)=ydeg
      call projin(z(1),z(2),xdeg,ydeg)
      z(1)=xdeg
      z(2)=ydeg
C      call io(z,nword,9,2,*70)
      call xiopst('w',9,id,z,itype,ierr)
      if(ierr.ne.0) stop 'error writing post file'
      go to 40
      else
c
c  ascii file
c
c      ians=0
c403   call aski4l('Do the ascii records include an id field?',ians,
c     1             ierr)
c      if(ierr.eq.-2) go to 402
      nchan=3
405   print*,'Enter the number of channels in the ascii file,'
      call aski4('  including x and y',nchan,ierr)
c   50 if(ians.eq.0) then
   50 read(8,*,end=90) (z(i),i=1,nchan)
c      else
c        read(5,100,end=90) (id,(z(i),i=1,nchan))
c 100  format(a,24g16.8)
c      endif
      go to (60,70) im
   60 call projfw(z(1),z(2),xkm,ykm)
      z(1)=xkm
      z(2)=ykm
      go to 80
   70 call projin(z(1),z(2),xdeg,ydeg)
      z(1)=xdeg
      z(2)=ydeg
c   80 if(ians.eq.0) then
   80 write(9,*) (z(i),i=1,nchan)
c      else
c        write(6,100) (id,(z(i),i=1,nchan))
c      endif
      go to 50
C      else
Cc
Cc     convert post to xyz file
Cc
C      ianom=ianom+4
C      go to (50,60),im
C   50 call io(z,nword,8,1,*70)
C      call projfw(z(3),z(4),xkm,ykm)
C      z(1)=xkm
C      z(2)=ykm
C      z(3)=z(ianom)
C      call io(z,3,9,2,*70)
C      go to 50
C   60 call io(z,nword,8,1,*70)
C      call projin(z(3),z(4),xdeg,ydeg)
C      z(1)=xdeg
C      z(2)=ydeg
C      z(3)=z(ianom)
C      call io(z,3,9,2,*70)
C      go to 60
C      endif
      endif
   90 CLOSE(8)
      STOP
      END
      subroutine proj
c      This subroutine is called thru its 3 entry points:
c           1. projst(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. projfw(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. projin(xkm,ykm,xdeg,ydeg) - does inverse projection.
c
c      Projections available are:
c            0 = none
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            8 = polar stereographic
c            9 = DNAG = spherical transverse mercator
c           14 = Lambert with standard parallels requested.
c           15 = Lambert with Antarctic parallels as used by the
c                German Antarctic North Victorialand Expedition
c                automatically shifts origin to x0=1444, y0=1123
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      common/prjcon/ ipr,cm,bl,xcm,ybl
      data sca/0.001/
      decdeg(deg,xmin)= deg + sign(xmin,deg)/60.


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry projst(iproj,cmdeg,bldeg)

      ipr=iproj
      cm=cmdeg
      bl=bldeg

      goto (10,20,90,40,50,60,70,90,22,5,5,5,5,42,44), ipr
      if (ipr.eq.0) return
  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.....Lambert conformal conic (ellipsoidal) with Antarctic parallels
 44   call lambset(cm,bl,72.66667,75.33333)
      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 projfw(xdeg,ydeg,xkm,ykm)

      goto (110,120,190,140,150,150,150,190,122,105,105,105,105,140,
     &   140),
     &   ipr
      if (ipr.eq.0) then
         xkm=xdeg
         ykm=ydeg
         return
      end if
 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)
      if (ipr.eq.15) then

c     shift origin for GANOVEX grids

      xkm=1444.-xkm
      ykm=1123.-ykm
      end if

       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 projin(xkm,ykm,xdeg,ydeg)

      goto (210,220,205,240,250,250,250,205,222,205,205,205,205,240,
     &    239),
     &    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

 239  continue

c     Shift origin of GANOVEX grids
      xkm=1444.-xkm
      ykm=1123.-ykm

 240  call lambinv(xkm,ykm,xdeg,ydeg)
       return

 250  call albinv(xkm,ykm,xdeg,ydeg)
       return

      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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

      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.2957795)
      real n,mfcn,mpfcn,m,m0,mn,mpn,ma,lam,lam0
      common/copoly/ 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


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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.2957795)
      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)
      t=tnphi**2
c     following line added 3/2/93
      csphi=cos(phi)
      c=epsq*csphi**2
      biga=csphi*(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



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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
      parameter (r=6371.204)
c      parameter (r=1.)
      parameter (deg2rad=0.17453293E-1, rad2deg=57.2957795)
      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=ykm/(r*k0) + phi0

      phi=asin(sin(d)/cosh(xkm/(r*k0)))
      lam=lam0 + atan(sinh(xkm/(r*k0))/cos(d))

      xdeg=lam*rad2deg
      ydeg=phi*rad2deg
      
      return
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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.2957795, 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


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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.2957795)
      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)
      double precision yl,xl,x1,y1
      yl=dble(ylat)*.0174532925199433d0
      xl=dble(abs(cm)-abs(xlon))*.0174532925199433d0
      go to (30,30,30,60,60,60,60,60),iprojt
   30 call merctr(yl,xl,x1,y1)
      go to 90
   60 call pstereo(yl,xl,x1,y1)
   90 x=sngl(x1)*sfac
      y=sngl(y1)*sfac
      return
      end
          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
          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
          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*********
c      subroutine io(z,n,idev,irw,*)
c      dimension z(n)
c      go to (10,20),irw
c   10 read(idev,end=30) z
c      return
c   20 write(idev) z
c      return
c   30 return 1
c      end
c*********
c      function lenrec(idev)
cc     find the number of words in a record
c      dimension iz(24)
c      do 20 lenrec=24,1,-1
c      read(idev,err=10) (iz(i),i=1,lenrec)
c      go to 30
c   10 rewind idev
c   20 continue
c   30 rewind idev
c      return
c      end

