c program prjrec
c
      dimension boug(6)
	double precision d2r,dp1,dp2
	character ifile*56,ofile*56,ifmt*56,ofmt*56,blank*56,
     1 pgm*8,dcon*3,bid*8,cf*50
      common/cparm/ifile,ofile,ifmt,ofmt,dcon
      common/rparm/baslat(3),cm(3),iproj,phi1,phi2,sca,s1,s2,
     & ianom,janom,ieast,irec,jrec
	data ind/5/,blank/' '/,pgm/'prjrec'/
	d2r=3.14159265d0/180.d0
        ifile=' '
        ofile=' '
        ifmt= ' '
        ofmt= ' '
        irec=2
        jrec=2
        dcon='deg'
        ieast=0
        ianom=1
        janom=0
        sca=.001
        s1=0.
        s2=0.
        iproj=0
        baslat(1)=0.
        baslat(2)=0.
        baslat(3)=0.
        cm(1)=0.
        cm(2)=0.
        cm(3)=0.
        phi1=33.
        phi2=45.
        do 121 i=1,6
        boug(i)=0.
121     continue
c
	print  1
1	format(' enter command filename or ''help'':')
	read(ind,2)cf
2	format(a)
	if(cf.eq.'help') print  130
130	format(' input keys (irec)',/,' 1 bin lon,lat,z',/,' 2 bin post',/,
     1 ' 3 formatted lon,lat,z',/,' 4 formatted ascii post',/,
     1 ' 5 formatted deg&min',/,' output keys (jrec) same as input',/,
     1 ' both default to 2')

	if (cf.eq.'help') stop
	open(unit=9,file=cf,form='formatted',status='old')
      call namepr(9)
	close(9)
	if(ifile.ne.' ') go to 40
	print  41
41	format(' input file :'$)
	read(ind,2) ifile
40	if(ofile.ne.' ') go to 42
	print  43
43	format(' output file :'$)
	read(ind,2) ofile
42	continue
	a=-1.
	if(ieast.ne.0) a=1.
	assign 51 to icon
	if(irec.eq.5) go to 30
	cdeg=1.
	if(dcon.eq.'min') cdeg=1./60.
	if(dcon.eq.'sec') cdeg=1./3600.
	if(cdeg.ne.1.0) assign 50 to icon
30	continue
c
	if(irec.lt.1 .or. irec.gt.5) stop ' 1>irec>5'
	if(jrec.lt.1 .or. jrec.gt.4) stop ' 1>jrec>4'
	if(janom.ne.0 .and. (jrec.eq.2 .or. jrec.eq.4))
     1 stop ' janom and jrec are inconsistent'
	if(irec.ge.3 .and. ifmt.eq.blank) stop ' need input format'
	if(jrec.ge.3 .and. ofmt.eq.blank) stop ' need output format'
	if(janom.le.0) janom=2
c
	if(irec.lt.3) open(10,file=ifile,form='unformatted',status='old')
	if(irec.ge.3) open(10,file=ifile,form='formatted',status='old')
	if(jrec.lt.3) open(unit=11,file=ofile,form='unformatted',
     1          status='unknown')
	if(jrec.ge.3) open(unit=11,file=ofile,form='formatted',
     1          status='unknown')
	if(iproj.ge.5 .and. iproj.le.7) call setalb(iproj)
	dp1=dble(phi1)*d2r
	dp2=dble(phi2)*d2r
	if(iproj.eq.4) call setlam(dp1,dp2)
	yb=baslat(1)+baslat(2)*1.6666667e-2+baslat(3)*2.77777e-4
	cmd=abs(cm(1))+cm(2)*1.6666667e-2+cm(3)*2.77777e-4
	cmd=sign(cmd,a)
	call prjctl(yb,cmd,xc,yc,cmd,sca,iproj)
	it=0
	ic=0
c
9	go to (21,22,23,24,25),irec
21	read(10,end=10)xd,yd,z
	boug(1)=z
	go to 13
22	read(10,end=10) bid,xd,yd,boug
	go to 12
23	read(10,ifmt,end=10) xd,yd,z
	boug(1)=z
	go to 13
24	read(10,ifmt,end=10)bid,xd,yd,(boug(i),i=1,ianom)
	go to 12
25	read(10,ifmt,end=10)bid,yd,ym,xd,xm,(boug(i),i=1,ianom)
	yd=yd+ym*1.666667e-2
	xd=abs(xd)+xm*1.666667e-2
12	z=boug(janom)
13	it=it+1
	xd=sign(xd,a)
	go to icon,(50,51)
50	xd=xd*cdeg
	yd=yd*cdeg
51	if(abs(yd).gt.90.0) go to 9
	call prjctl(yd,xd,x,y,cmd,sca,iproj)
	if(x.lt.-1.e15) go to 9
	y=y-yc
	go to (3,4,5,6),jrec
3	write(11)x,y,z
	go to 7
4	write(11)bid,x,y,boug
	go to 7
5	write(11,ofmt)x,y,z
	go to 7
6	write(11,ofmt) bid,x,y,(boug(i),i=1,ianom)
	go to 7
7	ic=ic+1
	go to 9
10	print  11,ic,it
11	format(i7,' stations output from',i7,' input')
	close(10)
	close(11)
	stop
	end
	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
	return
10	continue
	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,phi
	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 
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.
	implicit double precision (a-h,o-z)
	data e/8.1820567882165d-2/,r/6378160d0/,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
      subroutine namepr(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
      if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if((var(l).eq.blank.or.var(l).eq.comma).and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank.or.var(j).eq.comma) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 15
  120 m=l
      go to 15
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
c*********
      subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c     assigns values to proper variable
c     variables are passed to program prjrec through common blocks
c     numr=position in the array var where real variables start
c     numa=position in the array var where arrays start
c     nnvar=number of variables in program prjrec
c
      parameter(nnvar=18,numr=7)
      character*6 pvar,var(nnvar)
      character*56 tvar,kvar,cfmt
      logical chv
      character ifile*56,ofile*56,ifmt*56,ofmt*56,dcon*3
      common/cparm/ifile,ofile,ifmt,ofmt,dcon
      common/rparm/baslat(3),cm(3),iproj,phi1,phi2,sca,
     & s1,s2,ianom,janom,ieast,irec,jrec
      data var/'iproj','ianom','janom','ieast','irec','jrec',
     & 'phi1','phi2','sca','s1','s2','ifile','ofile','ifmt',
     & 'ofmt','dcon','baslat','cm'/
      numa=17
      do 190 i=1,nnvar
      if(pvar.ne.var(i)) go to 190
      if(.not.chv) then
c
c     noncharacter value
c     right justify then number in variable kvar
c
      m=57-nn
      im=m-1
      kvar(m:56)=tvar(1:nn)
      if(i.lt.numr) then
c
c     integer value
c
      write(cfmt,50) im,nn
   50 format('(',i2,'x,i',i2,')')
      read(kvar,cfmt) jvar
      else
c
c     real value
c
      write(cfmt,60) im,nn
   60 format('(',i2,'x,g',i2,'.0)')
      read(kvar,cfmt) xvar
      endif
      endif
      go to (101,102,103,104,105,106,107,108,109,110,111,112,
     & 113,114,115,116,117,118),i
  101 iproj=jvar
      go to 200
  102 ianom=jvar
      go to 200
  103 janom=jvar
      go to 200
  104 ieast=jvar
      go to 200
  105 irec=jvar
      go to 200
  106 jrec=jvar
      go to 200
  107 phi1=xvar
      go to 200
  108 phi2=xvar
      go to 200
  109 sca=xvar
      go to 200
  110 s1=xvar
      go to 200
  111 s2=xvar
      go to 200
  112 ifile=tvar(1:nn)
      go to 200
  113 ofile=tvar(1:nn)
      go to 200
  114 ifmt=tvar(1:nn)
      go to 200
  115 ofmt=tvar(1:nn)
      go to 200
  116 dcon=tvar(1:nn)
      go to 200
  117 baslat(inum)=xvar
      go to 200
  118 cm(inum)=xvar
      go to 200
  190 continue
      write(*,*) ' error in namelist - ',pvar,' variable not included'
      stop
  200 nvar=i
      return
      end
