c  program prjpt
c  one shot return of projection coordinates
	double precision phi1,phi2
	character ians*4
	sca=.001
100   print 1
1     format(' Enter projection number (0 for list, 99 to exit) :'$)
	read(5,*,err=100) iproj
      if(iproj.lt.1) then
      print 2
2       format(' 1=polyconic, 2=utm, 3=mercator, 4=lambert',/,
     1 ' 5,6,7= albers-us alaska hawaii, 8=polar stereographic,',/,
     1 ' 9= transverse mercator')
	go to 100
      endif
      if(iproj.gt.10) stop
	if(iproj.ge.5 .and. iproj.le.7) call setalb(iproj)
	if(iproj.eq.4) then
	phi1=33.d0
	phi2=45.d0
      print 4
4       format(' do you want to change the default standard',/,
     1 ' parallels of 33 and 45 degrees ?'$)
	read(5,6) ians
6       format(a4)
	if(ians(1:1).eq.'y' .or. ians(1:1).eq.'Y') then
      print 7
7       format(' enter new standard parallels (decimal deg.) :'$)
	read(5,*) phi1,phi2
	endif
	phi1=phi1*1.74532925d-2
	phi2=phi2*1.74532925d-2
	call setlam(phi1,phi2)
	endif
75    print 76
76    format
     1 (' Select degrees and minutes (1), or decimal degrees (2):'$)
      read(5,*) ideg
      if(ideg.ne.1.and.ideg.ne.2) go to 75
      if(ideg.eq.1) then
        print 8
8       format
     1 (' Enter base latitude (d,m), and central meridian (d,m):'$)
	read(5,*,err=100) yd,ym,xd,xm
        bl=yd+sign(ym,yd)/60.
        cm=xd+sign(xm,xd)/60.
      else
        print 85
85      format
     1 (' Enter base latitude (deg), and central meridian (deg):'$)
        read(5,*,err=100) bl,cm
      endif
	call prjctl(bl,cm,xdum,yb,cm,sca,iproj)
      print*,'Note: enter Ctrl-C to exit.'
10    if(ideg.eq.1) then
      print 15
15      format(' enter test lat (d,m), lon (d,m) :'$)
	read(5,*,err=100) yd,ym,xd,xm
c        if(yd.gt.90. .or. yd.lt.-90.) stop
	yd=yd+sign(ym,yd)/60.
	xd=xd+sign(xm,xd)/60.
      else
      print 155
155     format(' enter test lat (deg), lon (deg) :'$)
        read(5,*,err=100) yd,xd
      endif
        if(yd.gt.90. .or. yd.lt.-90.) stop
        call prjctl(yd,xd,x,y,cm,sca,iproj)
	y=y-yb
      print 20,x,y
20      format(' x=',f10.3,' km,  y=',f10.3,' km')
c      print*,' Continue ? [y/n]:'
c      read(5,6) ians
c      if(ians(1:1).eq.'n'.or.ians(1:1).eq.'N') stop
	go to 10
	end
c*******************************************************************************
	subroutine prjctl(ylat,xlon,x,y,cm,sfac,iprojt)
c  call projection routines
	double precision yl,xl,xp,yp,degr
	data degr/.0174532925199433d0/
	yl=dble(ylat)*degr
	xl=dble(xlon-cm)*degr
	go to (1,2,3,4,5,5,5,8,9),iprojt
1       call poly(yl,xl,xp,yp)
	go to 20
2     call utmfwd(yl,xl,xp,yp)
	go to 20
3       call merctr(yl,xl,xp,yp)
	go to 20
4       call lambert(yl,xl,xp,yp)
	go to 20
5       call albers(yl,xl,xp,yp)
	go to 20
8       call polars(yl,xl,xp,yp)
	go to 20
9       call transm(yl,xl,xp,yp)
20      x=sngl(xp)*sfac
	y=sngl(yp)*sfac
	return
	end
	subroutine albers(ylat,xlon,x,y)
	double precision ylat,xlon,x,y,n,rho1sq,sinbt1,twoc2n,
     1 nus,nals,nhaw,rho295,rho55,rho8,tcnus,tcnals,tcnhaw,
     1 a1,b,c1,d,e1,f1,g,h,theta,rho,sinbet,sinphi,s2
	data a1/9.954804334645587d-1/,b/4.492024607745888d-3/,
     1 c1/2.736435989866449d-5/,d/1.763992166249299d-7/,
     1 e1/1.160814577272288d-9/,f1/7.714265487727804d-12/,
     1 g/5.154557173568170d-14/,h/3.455700205911349d-16/,
     1 nus/6.02903493787094d-1/,nals/8.627447947235633d-1/,
     1 nhaw/2.241096394314637d-1/,rho295/8.49196923967458d13/,
     1 rho55/1.806308673895081d13/,rho8/7.943986660586285d14/,
     1 tcnus/1.346470921892769d14/,tcnals/9.409410848453636d13/,
     1 tcnhaw/3.622298555079059d14/,sbt295/4.907351753179611d-1/,
     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*
     1 (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
      subroutine poly(phi,dlamb,x,y)
      double precision phi,dlamb,xc,x,y,a,b,t,q,phi2,c1,c2
      a(xc)=
     1   6.378206402718907d 06 +xc*(
     1  -3.167517353503576d 06 +xc*(
     1   2.478805037574243d 05 +xc*(
     1  -3.530710396439220d 03 +xc*(
     1  -6.565371848240127d 02 +xc*(
     1   6.822539551727124d 01 +xc*(
     1  -2.888860980506611d 00 ))))))
      b(xc)=
     1   3.189103200618349d 06 +xc*(
     1  -2.115275857345996d 06 +xc*(
     1   4.144758431728325d 05 +xc*(
     1  -3.625295427368928d 04 +xc*(
     1   1.322429746943889d 03 +xc*(
     1   4.427183005616853d 01 +xc*(
     1  -8.630738701658902d 00 +xc*(
     1   4.279183401413811d-01 )))))))
      q(xc)=
     1   6.335034386662446d 06 +xc*(
     1   2.144094496614083d 04 +xc*(
     1  -4.182226973512467d 03 +xc*(
     1   3.609995316780635d 02 +xc*(
     1  -1.346978283534684d 01 ))))
      t(xc)=
     1   9.999999957157490d-01+xc*(
     1  -1.666665796975878d-01+xc*(
     1   8.333050613721043d-03+xc*(
     1  -1.980904608528695d-04+xc*
     1   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
	subroutine utmfwd(phi,dlam,x,y)
c  developed by g.i. evenden, usgs
	double precision phi,dlam,x,y,dl2,p
	x=1.d30
	y=1.d30
	p=phi*phi
	dl2=dlam*dlam
	if (p .le.1.94965360d0.and.dl2.le.3.7319881d-3)
     1    go to 10
	ier=-1
	return
10      ier=0
	y=phi*(
     1     6332500.47d0+p*(21431.67d0+p*(-4179.269d0+p*(
     1     359.981d0-p*13.267d0)))
     1 -dl2*(-3187827d0+p*(2114440d0+p*(-414363.5d0+p*(
     1     36344.6d0-p*1420.3d0)))
     1  -dl2*(1334935d0+p*(-2356027d0+p*(1371758d0-p*
     1     267852d0)))))
	x=dlam*(6375655.2d0+p*(-3166253.6d0+p*(247800.26d0+p*(
     1     -3569.65d0+p*(-617.35d0+p*50.89d0))))
     1  -dl2*(-1069479d0+p*(2661562d0+p*(-1785956d0+p*(
     1    485045d0-p*52022d0)))))
	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/,
     1 b2da2/9.932315290818186d-1/
	x=a*long
c       compute z/2
	z=0.5d0*(halfpi-datan(b2da2*dsin(lat)/dcos(lat)))
	y=a*dlog(dcos(z)/dsin(z))
	return
	end
      subroutine lambert(ylat,xlon,x,y)
      implicit double precision(a-h,o-z)
      save fn,fk
c
c   Lambert Conformal Conic forward projection program using
c    Clark 1866 ellipsoidal earth.  (See Map Projections Used by
c    the U. S. Geological Survey, GS Bulletin 1532, pp. 107-108.)
c
c      Input:  ylat - Latitude in radians.
c              xlon - Longitude in radians, east of central meridian.
c      Output: x    - Distance in meters east of central meridan.
c              y    - Distance in meters (negative) south of 90 degrees lat.
c
c    Note that if setlam entry is not called prior to calling
c     subroutine lambert, the constants fn & fk are set to calculate
c     x & y for standard parallels of 33. & 45. degrees.
c
      data a/6378206.4d0/,e2/.006768657997291099d0/
      data pi4/.785398163397448310d0/,e/.82271854223003258d-1/
      data b2da2/.993231342002708901d0/
      data fn/.63049989185603457d0/,fk/.124526547337527528d8/
      z=pi4-.5d0*datan(b2da2*(dtan(ylat)))
      r=fk*dtan(z)**fn
      theta=fn*xlon
      x=r*dsin(theta)
      y=-r*dcos(theta)
      return
c
      entry setlam(ylat1,ylat2)
c
c   This entry sets up the Lambert Conformal Conic constants for
c     standard parallels ylat1 & ylat2, entered in radians, using the
c     Clark 1866 ellipsoidal earth parameters of a = 6378206.4 meters
c     & b = 6356583.8 meters.  Note that e**2 = (a*a - b*b) / (a*a).
c
	if(ylat1.ne.ylat2 .and. ylat1*ylat2.gt.0.d0) go to 2
      print  1,ylat1,ylat2
1     format(2f10.3,' incorrect standard parallels,
     1 33 & 45 will be used')
	ylat1=.575958656d0
	ylat2=.785398167d0
2     cos1=dcos(ylat1)
      cos2=dcos(ylat2)
      sin1=dsin(ylat1)
      sin2=dsin(ylat2)
      esin1=e*sin1
      esin2=e*sin2
      fm1=cos1/dsqrt(1.d0-e2*sin1*sin1)
      fm2=cos2/dsqrt(1.d0-e2*sin2*sin2)
      t1=dtan(pi4-.5d0*datan(b2da2*(sin1/cos1)))
      t2=dtan(pi4-.5d0*datan(b2da2*(sin2/cos2)))
      fn=(dlog(fm1)-dlog(fm2))/(dlog(t1)-dlog(t2))
      fk=a*fm1/(fn*t1**fn)
      return
      end
	subroutine polars(ylat,xlon,x,y)
c   polar stereographic projection
c  partially checked by mike webring
c  there is about 1% error with polarr calculated this way,
c  but the form should be (1+e)/(1-e) (?) which makes bigger errors.
c  basic routine is due to kay edwards, flagstaff
	implicit double precision (a-h,o-z)
c        data e/8.1820567882165d-2/,r/6378160d0/,rad45/.785398163397448d0/,
        data e/8.1820567882165d-2/,rad45/.785398163397448d0/,
     1 ee/4.09102839410826d-2/,polarr/12713644.51d0/
	dlat=dabs(ylat)
c       polarr=(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=dcos(tane2)/dsin(tane2)
	tanz=tane1*(sine**ee)
	p=tanz*polarr
	y=-(p*dcos(xlon))
	x=p*dsin(xlon)
	if(ylat.lt.0.d0) y=-y
c       scale=p*dsqrt(1.d0-e*e*dsin(dlat)**2)/dcos(dlat)/r
	return
	end
      subroutine transm(ylat,xlon,x,y)
c
c   Transverse Mercator forward projection program using spherical
c     earth of radius 6371204. meters.  This is the radius of a
c     spherical earth of equivalent surface area to the Clark 1866
c     ellipsoid.  Note that the constant rsf is the product of this
c     radius and a scale factor of 0.926.  The scale factor is present
c     to balance the scale errors over North America.
c
c     input:  ylat - Latitude in radians.
c             xlon - Longitude in radians, east of central meridian.
c     output: x    - Distance in meters east of central meridian.
c             y    - Distance in meters north of equator.
c
      double precision ylat,xlon,x,y,rsf,b
      data rsf/5899734.904d0/
      b=dcos(ylat)*dsin(xlon)
      x=0.5d0*rsf*dlog((1.d0+b)/(1.d0-b))
      y=rsf*datan(dtan(ylat)/dcos(xlon))
      return
      end
