c  program to map a grid incremented in
c  lat-lon onto a rectangular coordinate system.
c  restriction on projection types: meridians must be
c  straight lines, with a unit of latitude represented
c  by a constant amount of y displacement.
c  implemented projections are: utm, lambert, and albers.
c
c  the central meridian affects the rotation of the data
c  and should match published maps.  for instance
c  1:250,000 topo sheets might be 1 degree from the center
c  of a utm zone, but to produce a map squared on a plot
c  project from the center.  grids that will be merged
c  together should be projected with the same central meridian.
c
c    coded by Mike Webring
c    Geophysics branch, USGS
c
	common /prjgrd/ buf(30000),x(1500),z(1500),zint(1500)
	common /spl15/ wrk(7500)
	double precision xd,yd,ryd,d,ddr,ddr2,y2or,xx,xa,phi1,phi2
	double precision cord(2,1500)
c
      character    ians*16,ians1*1
	character    id*56, blank*56, id2*56, p*8, p2*8
        character*56 ifile, jfile
	logical ysplin
c
        data inp/5/, p2/'prjgrd'/, sca/.001/,
     1  maxcr/1500/, dval/0.170412e+39/, isw/12/,jsw/13/,ksw/14/
	blank = ' '
c
1     print  5
5	format(' enter projection number :'$)
	read *, ipr
      if(ipr.lt.1.or.ipr.gt.9) then
      print 10
10	format(' polyconic=1, utm=2, lambert conformal=4,',/
     1 ' albers equal area us=5, alaska=6, hawaii=7,',/
     1 ' polar stereographic=8, transverse mercator=9')
      go to 1
      else
      endif
	if(ipr.eq.3) stop ' mercator not allowed'
	if(ipr.ge.5 .and. ipr.le.7) call setalb(ipr)
	phi1 = 33.d0
	phi2 = 45.d0
	if(ipr.ne.4) go to 40
15    print 20
20	format(' do you want to change the lambert standard parallels ?'$)

      read(inp,37) ians
37    format(a)
	if(ians(1:1).eq.'n' .or. ians(1:1).eq.'N') go to 40
	if(.not.(ians(1:1).eq.'y' .or. ians(1:1).eq.'Y')) go to 15
      print  30
30	format(' enter phi1 and phi2 (deg) :'$)
	read *, phi1,phi2
40	phi1=phi1*3.1415927d0/180.
	phi2=phi2*3.1415927d0/180.
	call setlam(phi1,phi2)
c
      print  50
50	format(' input grid : '$)
	read(inp,37)ifile
	open (unit=isw,file=ifile,form='unformatted',status='old')
	read(isw) id,p, nc,nr,nz, xo,dx,yo,dy
      print 60,  id,p, nc,nr,    xo,yo,dx,dy
60	format(/,1x,a56,'pgm:',a8,/,' ncol=',i5,'  nrow=',i5,/,
     1 ' xo,yo= ',1p2g13.5,'  dx,dy= ',2g11.3/)
c
70    print  80
80    format(' input units (1=deg,2=min,3=sec) :'$)
	unit=0.
      read*, iunit
      if(iunit.eq.1)unit=1.
      if(iunit.eq.2)unit=1./60.
      if(iunit.eq.3)unit=1./3600.
	if(unit.eq.0.0) go to 70
	xo=xo*unit
	yo=yo*unit
	dx=dx*unit
	dy=dy*unit
c
      print *, ' projected output grid :'
	read(inp,37) jfile
	open (unit=jsw, file=jfile, form='unformatted', status='new' )
      print *,' title, <car ret> to leave the same  '
	read(inp,37) id2
	if ( id2 .eq. blank ) id2 = id
      print *, ' output dx,dy in km: '
	read *, dx2, dy2
c
      print *, ' central meridian (NEGATIVE is west longitude)'
      read *, cm
	if ( xo .gt. 0 ) then
        print *, ' grid origin is', xo
        print *, ' is this west longitude ?'
	  if ( noyes(ldum) .eq. 1 ) xo = -xo
	endif
      print *, ' base latitude :'
      read *, blat
	if ( abs( xo - cm ) .lt. 1.e-4) xo = cm
	dx = abs(dx)
	call prjctl( blat, cm, dum, yb, cm, sca, ipr )
c
c  get output dimensions
	xlon = xo
	rlt2 = yo + (nr-1)*dy
139   print  140
140	format(' change area covered ?')
      read(inp,37)ians1
      if(ians1.eq.'y' .or. ians1.eq.'Y') go to 150
      if(.not.(ians1.eq.'n' .or. ians1.eq.'N')) go to 139
	call prjctl(yo,cm,dum,yo2,cm,sca,ipr)
	yo3=yo2-yb
	call prjctl(yo,xo,xo2,dum,cm,sca,ipr)
	rln2=xo+(nc-1)*dx
	call prjctl(yo,rln2,xmax,dum,cm,sca,ipr)
	call prjctl(rlt2,xo,dum,ymax,cm,sca,ipr)
	call prjctl(rlt2,rln2,dum,ymax2,cm,sca,ipr)
	if(ymax2.gt.ymax) ymax=ymax2
	nr2=((ymax-yo2)/dy2)+1
	nc2=((xmax-xo2)/dx2)+1
      print 145,nc2,nr2
145	format(' output dimensions ',i5,' x',i5)
	go to 170
150	continue
      print  160
160	format(' enter new xo (km), yo (km), ncol, nrow :')
	read *, xo2,yo3,nc2,nr2
	yo2=yo3+yb
c
170	if(nc2.gt.maxcr .or. nr2.gt.maxcr) go to 890
	n2=int(float(nr2)/20.+.999)
	n=int(float(nr)/20.+.999)
	m=int(float(nc)/20.+.999)
	n10=n*20
	m10=m*20
	nmax=n
	if(n2.gt.n) nmax=n2
	nx10=nmax*20
      open (ksw,access='direct',form='unformatted',
     1 status='scratch',recl=1600)
c
c  chk utm width
	if(ipr.ne.2) go to 200
	if(abs(xo-cm).gt.3.5) go to 180
	chk=xo+dx*float(nc)
	if(abs(chk-cm).lt.3.5) go to 200
180   print  190
190	format(' error: grid extends 3 degrees beyond central meridian')
	go to 999
c
200	call rio(buf,m,nmax,1,1,m10,-1,isw,ksw)
c
c  spline in y
 	ysplin = .true.
	ic    = 1
	rnr   = 1.0 / float(nr-1)
	do 300 ibuf = 1, m
	loc = ibuf
	call rio(buf,nmax,0,loc,m,nx10,1,isw,ksw)
	iof = 0
	do 290 ir = 1, 20
c
	if (ysplin) then
	call prjctl(yo,xlon,px1,py1,cm,sca,ipr)
	call prjctl(rlt2,xlon,px2,py2,cm,sca,ipr)
	cord(1,ic)=dble(px1)
	xd=dble(px2)-cord(1,ic)
	yd=dble(py2)-dble(py1)
	if(yd.eq.0.d0) stop ' zero width input'
	ryd=1.d0/yd
	d=dsqrt(xd*xd+yd*yd)
	ddr=d*dble(rnr)
	dr=sngl(ddr)
	cord(2,ic)=dble(dy2)*xd*ryd
	ddr2=dble(dy2)*d*ryd
	dr2=sngl(ddr2)
	y2or=dble(yo2)-dble(py1)
	xx=y2or*xd*ryd
	cord(1,ic)=cord(1,ic)+xx
	xa=-y2or*d*ryd
	x(1)=sngl(xa)
	z(1)=buf(1+iof)
        do 90 i = 2, nr
	  z(i) = buf(i+iof)
	  xa = xa + ddr
	  x(i) = sngl(xa)
90      continue
	call spl1k(nr,x,z,dr,0.,nx10,buf(1+iof),dr2)
	if (ic .eq. nc) ysplin = .false.
c
	else
        do 91 i = 1 + iof, nx10 + iof
	  buf(i) = dval
91      continue
	endif
c
	ic   = ic + 1
	iof  = iof + nx10
290	xlon = xlon + dx
	loc  = ibuf
300	call rio(buf,n2,0,loc,m,nx10,0,isw,ksw)
c
c  spline in x
	ir = 1
	loc = 1
	write(jsw) id2,p2, nc2,nr2,nz, xo2,dx2,yo3,dy2, ipr,cm,blat
        do 94 ibuf = 1, n2
	  call rio(buf,m,0,loc,1,m10,1,isw,ksw)
          do 93 j = 1, 20
	    iof = (j-1) * m10
            do 92 i = 1, nc
	      x(i) = sngl(cord(1,i))
	      cord(1,i) = cord(1,i) + cord(2,i)
92          continue
	    call spl1k(nc,x,buf(1+iof),0.,xo2,nc2,zint,dx2)
	    call rowio(nc2,zint,0,isw,jsw,iend)
	    if (ir.eq.nr2) go to 999
	    ir = ir + 1
93        continue
	  loc=loc+m
94      continue
c
999	continue
	close(ksw)
	close(isw)
	close(jsw)
	stop
890   print  891,nc2,nr2
891	format(' exceeded output dimension with ncol=',i6,' nrow=',i6)
9999  print  9998,py2,py1,yo,rlt2,cord(1,ic),px2
9998	format(6e12.4)
	stop
	end
c*******************************************************************************
	subroutine rio(z,nblk,nbuf,loc2,idl,n10,iop,isw,ksw)
	dimension hdr(16),w(400),z(1)
	logical noread
      data dval/0.170412e+39/
	loc=loc2
	assign 999 to ir
	if(iop)1,2,3
1	assign 99 to ir
	noread=.false.
	itot=0
	rewind(isw)
	read(isw) hdr,nc,nr,nz,xo,dx,yo,dy
	do 100 ibuf=1,nbuf
	do 5 i=1,20*n10
5	z(i)=dval
	ii=1
	do 10 irow=1,20
	if(noread) go to 10
	call rowio(nc,z(ii),-1,isw,isw,ie)
	if(ie.eq.0) itot=itot+1
	if(ie.eq.1 .or. itot.ge.nr) noread=.true.
10	ii=ii+n10
	go to 2
99	continue
100	continue
      if(itot.ne.nr) print  101,itot,nr
101	format(' found only ',i4,' of ',i4,' rows')
	return
c
2	do 30 iblk=1,nblk
	lw=1
	ls=(iblk-1)*20+1
	do 21 j=1,20
	l2=ls
	l=lw
	do 20 i=1,20
	w(l)=z(l2)
	l=l+20
20	l2=l2+1
	lw=lw+1
21	ls=ls+n10
      write(ksw,rec=loc)w
30	loc=loc+idl
	go to ir,(99,999)
c
3	do 60 iblk=1,nblk
      read(ksw,rec=loc) w
	loc=loc+idl
	l=1
	ls=(iblk-1)*20+1
	do 55 j=1,20
	l2=ls
	do 50 i=1,20
	z(l2)=w(l)
	l=l+1
50	l2=l2+1
55	ls=ls+n10
60	continue
999	return
	end
c*******************************************************************************
	subroutine spl1k(na,x,z,dx,xo,na2,zint,dx2)
	dimension d(2),z(na),x(na),zint(na2)
	common /spl15/a(1500),b(1500),c(1500)
     1  ,w1(1500),w2(1500)
	logical if
      data d/0.,0./,dval/0.170412e+39/
	if=.false.
	do 10 i=1,na2
10	zint(i)=dval
	do 100 i=1,na
	if(.not.if) go to 1
	if(z(i).gt.1.e37) go to 2
3	ni=ni+1
	if(i.eq.na) go to 4
	go to 100
1	if(z(i).gt.1.e37) go to 100
	if=.true.
	ni=0
	ib=i
	go to 3
2	if=.false.
4	if(ni.lt.3) go to 100
	call spln1(ni,dx,x(ib),z(ib),a,b,c,0,d,w1,w2)
	ib3=ib-1
	m=ni-1+ib
	xd=abs(xo-x(ib))
	if(xo.ge.x(ib)) go to 22
	i2=ib
	ib2=int(xd/dx2+1.999999)
	xx=xo+(ib2-1)*dx2
24	if(xx.ge.x(i2)) go to 30
c     print  26
26	format(' index error')
	xx=xx+dx2
	go to 24
22	ib2=1
	xx=xo
	i2=1
30	if(xx.gt.x(m)) go to 100
31	if(xx.gt.x(i2+1)) i2=i2+1
	if(xx.gt.x(i2+1)) go to 31
	xd=xx-x(i2)
	j=i2-ib3
	if(ib2.gt.na2) return
	zint(ib2)=z(i2)+((c(j)*xd+b(j))*xd+a(j))*xd
	ib2=ib2+1
	xx=xx+dx2
	go to 30
100	continue
	return
	end
c*******************************************************************************
	subroutine rowio(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
c*******************************************************************************
	subroutine spln1(m,h,x,y,a,b,c,t,d,p,s)
c	coded by W. L. Anderson, U.S.G.S.
	dimension x(m),y(m),a(m),b(m),c(m),d(2),p(m),s(m)
	real mul
	integer t
	if(t.lt.0 .or. t.gt.1 .or. h.lt.0. .or. m.lt.3) go to 999
	n=m-1
	if(t.eq.0) go to 20
c  1ST DERIV BOUNDS GIVEN
	ne=n-1
	if(h)999,11,1
c EQUAL SPACING H>0  T=1
1	hh=3./h
	do 2 i=1,ne
	b(i)=4.
	c(i)=1.
	a(i)=1.
2	p(i)=hh*(y(i+2)-y(i))
	p(1)=p(1)-d(1)
	p(ne)=p(ne)-d(2)
c SOLN OF TRIDIAGONAL MATRIX ORDER NE
3	c(1)=c(1)/b(1)
	p(1)=p(1)/b(1)
	do 4 i=2,ne
	mul=1./(b(i)-a(i)*c(i-1))
	c(i)=mul*c(i)
4	p(i)=mul*(p(i)-a(i)*p(i-1))
c OBTAIN SPLINE COEF
	a(ne+t)=p(ne)
	i=ne-1
5	a(i+t)=p(i)-c(i)*a(i+t+1)
	i=i-1
	if(i.ge.1) go to 5
	if(t.eq.0) go to 6
	a(1)=d(1)
	a(m)=d(2)
6	if(h.eq.0.) go to 14
	hh=1./h
	do 7 i=1,n
	i1=i+1
	mul=hh*(y(i1)-y(i))
	b(i)=hh*(3.*mul-(a(i1)+2.*a(i)))
7	c(i)=hh*hh*(-2.*mul+a(i1)+a(i))
	return
c UNEQUAL SPACING H=0  T=1
11	do 12 i=1,n
	i1=i+1
12	s(i1)=x(i1)-x(i)
	do 13 i=1,ne
	i1=i+1
	i2=i+2
	b(i)=2.*(s(i1)+s(i2))
	c(i)=s(i1)
	a(i)=s(i2)
	sq1=s(i1)*s(i1)
	sq2=s(i2)*s(i2)
13	p(i)=3.*(sq1*(y(i2)-y(i1))+sq2*(y(i1)-y(i)))/(s(i1)*s(i2))
c       B(I)=2.*(S(I+1)+S(I+2))
c       C(I)=S(I+1)
c       A(I)=S(I+2)
c13     P(I)=3.*(S(I+1)**2*(Y(I+2)-Y(I+1))+S(I+2)**2*(Y(I+1)-Y(I)))/
c      1  (S(I+1)*S(I+2))
	p(1)=p(1)-s(3)*d(1)
	p(ne)=p(ne)-s(n)*d(2)
	go to 3
14	do 15 i=1,n
	i1=i+1
	hh=1./s(i1)
	hs=hh*hh
	mul=(y(i1)-y(i))*hs
	b(i)=3.*mul-(a(i1)+2.*a(i))*hh
15	c(i)=-2.*mul*hh+(a(i1)+a(i))*hs
	return
c  2ND DERIV BOUNDS GIVEN
20	ne=n+1
	if(h)999,31,21
c  EQUAL SPACING H>0  T=0
21	hh=3./h
	do 22 i=2,n
	b(i)=4.
	c(i)=1.
	a(i)=1.
22	p(i)=hh*(y(i+1)-y(i-1))
	b(1)=2.
	b(ne)=2.
	c(1)=1.
	c(ne)=1.
	a(ne)=1.
	p(1)=hh*(y(2)-y(1))-.5*h*d(1)
	p(ne)=hh*(y(m)-y(n))+.5*h*d(2)
	go to 3
c  UNEQUAL SPACING H=0  T=0
31	do 32 i=1,n
	i1=i+1
32	s(i1)=x(i1)-x(i)
	n1=n-1
	do 33 i=1,n1
	i1=i+1
	i2=i+2
	b(i1)=2.*(s(i1)+s(i2))
	c(i1)=s(i1)
	a(i1)=s(i2)
	sq1=s(i1)*s(i1)
	sq2=s(i2)*s(i2)
33	p(i1)=3.*(sq1*(y(i2)-y(i1))+sq2*(y(i1)-y(i)))/(s(i1)*s(i2))
c       B(I+1)=2.*(S(I+1)+S(I+2))
c       C(I+1)=S(I+1)
c       A(I+1)=S(I+2)
c33     P(I+1)=3.*(S(I+1)**2*(Y(I+2)-Y(I+1))+S(I+2)**2*(Y(I+1)-Y(I)))/
c      1 (S(I+1)*S(I+2))
	b(1)=2.
	b(ne)=2.
	c(1)=1.
	c(ne)=1.
	a(ne)=1.
	p(1)=3.*(y(2)-y(1))/s(2)-.5*s(2)*d(1)
	p(ne)=3.*(y(m)-y(n))/s(m)+.5*s(m)*d(2)
	go to 3
999	m=-abs(m)
	return
	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 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 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 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
c******************************************************************************
	function noyes()
c  returns -1 for indeterminate, 0 for no, 1 for yes
	character  inp*4, blank*4
	data       blank/ '    ' /

	icnt  = 0
	noyes = -1

1	inp   = ' '
	read ( 5, 2, end=5, err=5 ) inp
2	format( a )
	if ( inp .eq. blank ) return

	do 3 n = 1, 4
	  if ( inp(n:n) .ne. ' ' ) then
	    if ( inp(n:n).eq.'y' .or. inp(n:n).eq.'Y' ) noyes = 1
	    if ( inp(n:n).eq.'n' .or. inp(n:n).eq.'N' ) noyes = 0
	    if ( noyes .gt. -1 ) return
	    go to 5
	  endif
3	continue

5	if ( icnt .lt. 3 ) then
	  icnt = icnt + 1
	  write( 6, * ) ' <cr>, no, or yes :'
	  go to 1
	endif

	write( 6, * ) ' count exceeded, returning indeterminate state '
9	return
	end
