C
C________________________________________________________________
C
C     SUBROUTINE  L L P O S T
C________________________________________________________________
C
        subroutine llpost
        common /llp/ latm(3),latx(3),longm(3),longx(3),cm(3),
     1               baslat(3),
     1    iproj,xxx(2),yyy(2),sizep,unit,iplotr,neat,tint,itpost,
     1    ib,phi1,phi2
        character*1 ctype
        common/captur/icapt,iwhite,ctype
        dimension scaf(3),range(2),xx(2),yy(2)
        real latm,latx,longm,longx
        character fmt*16
        logical nos,nom
        data d2r/1.7453292e-2/, conm/1.666667e-2/,
     1       cons/2.777778e-4/,
     1    scaf/39.370079,1.,.001/, icross/6/,
     1    nmax/10/
        dms2d(a,b,c)=a+sign(b,a)*conm+sign(c,a)*cons
c
        if(ctype.eq.'l') then
        icol=1
        else
        icol=-1
        endif

        if(iproj.lt.1 .or.iproj.gt.10) return
        sca=scaf(int(unit)+1)
        do 1 i=1,2
        xx(i)=xxx(i)/(1000.*sca)
1       yy(i)=yyy(i)/(1000.*sca)
c       changes made by r.godson 3/90 to substitute xxx & yyy
c       for xx & yy in several places below to allow proper
c       labelling when units are not in kilometers.
        hfx=.5*(xxx(1)+xxx(2))
        hfy=.5*(yyy(1)+yyy(2))
        bld=dms2d(baslat(1),baslat(2),baslat(3))
        cmd=dms2d(cm(1),cm(2),cm(3))
        if(abs(bld).gt.90. .or. abs(cmd).gt.360.) return
        call prjctl(bld,cmd,xb,yb,cmd,sca,iproj)
        ydm=dms2d(latm(1),latm(2),latm(3))
        ydx=dms2d(latx(1),latx(2),latx(3))
        xdx=dms2d(longm(1),longm(2),longm(3))
        xdm=dms2d(longx(1),longx(2),longx(3))
c  adjust cem for conic projection distance to equator
        cem=4.0e7
        if(iproj.eq.4 .and. abs(bld).lt.15.) cem=4.21e7
        if(iproj.gt.4 .and. abs(bld).lt.15.) cem=3.84e7
        ce=cem*sca
        dt=(yy(2)-yy(1))*360./(ce*5.)
c  min-max latitude of area
        if(ydm.ge.ydx) go to 9
        if(ydm.ne.0.0 .and. ydx.ne.0.0) go to 10
9       ydm=bld+(yy(1)*360.)/ce-dt
        ydx=bld+(yy(2)*360.)/ce+dt
        if(ydm.lt.-90.0) ydm=-90.0
        if(ydx.gt.90.0)  ydx=90.0
        latx(1)=ydx
        latm(1)=ydm
c  min-max longitude of area
10      if(xdm.ne.0.0 .and. xdx.ne.0.0) go to 20
        cosp=cos(d2r*amax1(abs(ydx),abs(ydm)))
        if(cosp.lt..17364) cosp=.17364
        xdm=cmd+(xx(1)*360.)/(cosp*ce)-dt
        xdx=cmd+(xx(2)*360.)/(cosp*ce)+dt
        longm(1)=xdx
        longx(1)=xdm
c  default tick interval for less than nmax ticks per side
c  result is multiple of 1, 2 or 5 deg, min, or, sec
20      if(tint.gt.0.0) go to 30
        range(1)=0.
        range(2)=amax1((xdx-xdm),(ydx-ydm))
        tint=0.
        dr=range(2)*60./float(nmax)
        if(dr.ge.60.) then
         call setax(range,tint,nmax,ndum,fmt)
         tint=tint*60.
        else if(dr.ge.1. .and. dr.lt.60.) then
         range(2)=range(2)*60.
         call setax(range,tint,nmax,ndum,fmt)
        else
         range(2)=range(2)*3600.
         call setax(range,tint,nmax,ndum,fmt)
         tint=tint/60.
        endif
c  final initialization
30      if(itpost.le.0) itpost=2
        dt=(aint(tint*60.+.0001))/3600.
        dt2=dt*float(itpost)
        dt3=dt2*60.0
        nos=.false.
        if(dt3-aint(dt3).eq.0.0) nos=.true.
        nom=.false.
        if(amod(dt3,60.0).eq.0.0) nom=.true.
        if(yy(2)-yy(1).lt.20.)  ydm=aint(60.*ydm)/60.
        if(yy(2)-yy(1).ge.20.)  ydm=aint(ydm+90.)-90.
        if(yy(2)-yy(1).gt.500.) ydm=5.*aint((ydm+90.)/5.) - 90.
        if(xx(2)-xx(1).lt.20.)  xdm=aint(60.*xdm)/60.
        if(xx(2)-xx(1).ge.20.)  xdm=aint(xdm+360.)-360.
c     utm has a limit of 3 degrees
        if(xx(2)-xx(1).gt.500. .and. iproj.ne.2)
     1    xdm=5.*aint((xdm+360.)/5.) - 360.
c  plot longitude labels
        if(sizep.le.0.0) return
        xd=xdm
100     if(xd.gt.xdx) go to 190
        call dms(xd,id,im,is)
        id=id
        yd=ydm
105     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        if(x.ge.xxx(1).and.x.le.xxx(2) .and.
     1    y.ge.yyy(1).and.y.lt.hfy) go to 109
        yd=yd+dt
        if(yd-ydx) 105,105,110
109     call latlab(-1.,0,sizep,x,yyy(1),id,im,is,nos,nom)
110     if(iplotr.eq.1 .or. iplotr.eq.4.or.iplotr.gt.7) go to 120
        yd=ydx
115     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        if(x.ge.xxx(1).and.x.le.xxx(2) .and.
     1    y.gt.hfy.and.y.le.yyy(2)) go to 119
        yd=yd-dt
        if(yd-ydm) 120,115,115
119     call latlab(1.,0,sizep,x,yyy(2),id,im,is,nos,nom)
120     xd=xd+dt2
        go to 100
c  plot latitude labels
190     yd=ydm
200     if(yd.gt.ydx) go to 300
        call dms(yd,id,im,is)
        xd=xdm
205     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        if(y.ge.yyy(1).and.y.le.yyy(2) .and.
     1    x.ge.xxx(1).and.x.lt.hfx) go to 209
        xd=xd+dt
        if(xd-xdx) 205,205,210
209     call latlab(-1.,1,sizep,xxx(1),y,id,im,is,nos,nom)
210     if(iplotr.eq.1 .or. iplotr.eq.4.or.iplotr.gt.7) go to 220
        xd=xdx
215     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        if(y.ge.yyy(1).and.y.le.yyy(2) .and.
     1    x.gt.hfx.and.x.le.xxx(2)) go to 219
        xd=xd-dt
        if(xd-xdm) 220,215,215
219     call latlab(1.,1,sizep,xxx(2),y,id,im,is,nos,nom)
220     yd=yd+dt2
        go to 200
c  plot interior lat-lon ticks
300     yd=ydm
        do 360 j=1,int((ydx-yd)/dt)+1
        xd=xdm
        do 350 i=1,int((xdx-xd)/dt)+1
        call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        if(x.lt.xxx(1) .or. x.gt.xxx(2)) go to 350
        if(y.lt.yyy(1) .or. y.gt.yyy(2)) go to 350
        if(sizep.lt..00025) then
          call vchar(x,y,icross,1,2*icol,2000.0*sizep,0.,0.,0.)
        else
          call vchar(x,y,icross,1,2*icol,2.0*sizep,0.,0.,0.)
        end if
350     xd=xd+dt
360     yd=yd+dt
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  D M S
C________________________________________________________________
C
        subroutine dms(d,id,im,is)
        is=int(abs(d)*3600.+.5)
        id=is/3600
        im=(is-id*3600)/60
        is=is-(id*3600+im*60)
        id=isign(id,int(d))
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  L A T L A B
C________________________________________________________________
C
        subroutine latlab(side,lat,sz,x,y,jd,jm,js,nosec,nomin)
        character*1 ctype
        common/captur/icapt,iwhite,ctype
        dimension labint(3)
        character label*12,jd2*1,jm2*1,js2*1
        equivalence (labint,label)
        logical nosec,nomin
        jd2=char(24)
        jm2=char(39)
        js2=char(34)
        write(label,50) iabs(jd),jd2,jm,jm2,js,js2
50      format(i4,a1,2(i2,a1))
        if(jm.ne.0 .or. js.ne.0) nomin=.false.
        if(js.ne.0) nosec=.false.
        nch=11
        if(nosec) nch=8
        if(nomin) nch=5
        xl=nch*sz
        yof=0.
        if(lat.eq.1) then
        if(side.lt.0.) xof=-(2.*sz+xl)
        if(side.gt.0.) xof=2.*sz
        else
        xof=-xl/2.
        yof=2.*sz*side
        endif
        call vchar(x,y,labint,nch,(iwhite*100)+2,sz,0.,xof,yof)
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  P R J C T L
C________________________________________________________________
C
        subroutine prjctl(ylat,xlon,x,y,cmer,sfac,iprojt)
        common /llp/ latm(3),latx(3),longm(3),longx(3),cm(3),
     1               baslat(3),
     1    iproj,xxx(2),yyy(2),sizep,unit,iplotr,neat,tint,itpost,
     1    ib,phi1,phi2
c  call projection routines
        double precision yl,xl,x1,y1,degr,p1,p2
        data degr/.0174532925199433d0/
        yl=dble(ylat)*degr
        xl=dble(xlon-cmer)*degr
        go to (1,2,3,4,5,5,5,8,9),iprojt
1       call poly(yl,xl,x1,y1)
        go to 20
2       call utmfwd(yl,xl,x1,y1)
        go to 20
3       call merctr(yl,xl,x1,y1)
        go to 20
4       p1=dble(phi1)*degr
        p2=dble(phi2)*degr
        call setlam(p1,p2)
        call lamber(yl,xl,x1,y1)
        go to 20
5       call setalb(iprojt)
        call albers(yl,xl,x1,y1)
        go to 20
8       call polars(yl,xl,x1,y1)
        go to 20
9       call transm(yl,xl,x1,y1)
20      x=sngl(x1)*sfac
        y=sngl(y1)*sfac
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  P O L Y
C________________________________________________________________
C
      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
C
C________________________________________________________________
C
C     SUBROUTINE  U T M F W D
C________________________________________________________________
C
        subroutine utmfwd(phi,dlam,x,y)
c  developed by g.i. evenden, usgs
        double precision phi,dlam,x,y,dl2,p
        p=phi*phi
        dl2=dlam*dlam
        if (p .le.1.94965360d0.and.dl2.le.3.7319881d-3)
     1       go to 10
        x=1.d30
        y=1.d30
        return
10      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
C
C________________________________________________________________
C
C     SUBROUTINE  M E R C T R
C________________________________________________________________
C
        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
C
C________________________________________________________________
C
C     SUBROUTINE  L A M B E R
C________________________________________________________________
C
      subroutine lamber(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
c                     central meridian.
c      Output: x    - Distance in meters east of central meridan.
c              y    - Distance in meters (negative) south of
c                     90 degrees lat.
c
c    Note that if setlam entry is not called prior to calling
c     subroutine lamber, the constants fn & fk are set
c     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
C________________________________________________________________
C
C     ENTRY  S E T L A M
C________________________________________________________________
C
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
c     the
c     Clark 1866 ellipsoidal earth parameters of a = 6378206.4
c     meters
c     & b = 6356583.8 meters.
c     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, 33 & 45',
     1                ' 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
C
C________________________________________________________________
C
C     SUBROUTINE  A L B E R S
C________________________________________________________________
C
        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/,
     1    tcnals/9.409410848453636d13/,
     1    tcnhaw/3.622298555079059d14/,
     1    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________________________________________________________________
C
C     ENTRY  S E T A L B
C________________________________________________________________
C
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  P O L A R S
C________________________________________________________________
C
        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
c  errors.
c  basic routine is due to kay edwards, flagstaff
        implicit double precision (a-h,o-z)
        data e/8.1820567882165d-2/,
     1       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
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  T R A N S M
C________________________________________________________________
C
      subroutine transm(ylat,xlon,x,y)
c
c   Transverse Mercator forward projection program using
c   spherical
c     earth of radius 6371204. meters.  This is the radius of a
c     spherical earth of equivalent surface area to the Clark
c     1866
c     ellipsoid.  Note that the constant rsf is the product of
c     this
c     radius and a scale factor of 0.926.  The scale factor is
c     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
c                    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
