C  Conversions to and from the geo (geographic) and the map
C  (Map-oriented)  coordinate systems.
C
C  The geo system is a 3-D Cartesian coordinate system, with origin
C  at the Center of the Earth, with the x_3 axis pointing to the North
C  Pole, the x_1 axis pointing to the Equator at the Greenwich Meridian,
C  and the x_2 axis by the right hand rule pointing to the Equator
C  at 90 E.
C
C  In the map system, the axis of the map projection passes through the
C  center of the Earth, and through a point on the surface that acts as the
c  "Pole" of the map (Which will coincide with the North or South Pole,
C  unless the projection is Oblique or Transverse).  The x_3 axis of the
C  map coordinate system is aligned with this map pole.  In Lambert and
C  Mercator projections, a "Cut" extends from the map pole along a great
C  circle to its antipode; the unrolled map is discontinuous there.  The
C  x_1 axis of the map coordinate system is diametrically opposite this
C  cut, and 90 degrees from the pole; the x_2 axis is selected to complete
C  a right hand rule.
C
C  The coefficients of the map coordinate system relative to the geo system
C  are given in elements 2 through 10 of the stcprm array.
C
      subroutine basg2m(stcprm, xin, xout)
C  receives the vector xin, giving the components of a point in the geo
C  system; returns xout, the components of the same point in the map system.
      dimension stcprm(15),xin(3),xout(3)
      ifind(l,k) = -2 + 3*l + k
        do k=1,3
          xout(k) = 0.
        enddo
        do l=1,3
          do k=1,3
            xout(l) = xout(l) + xin(k) * stcprm(ifind(l,k))
          enddo
        enddo
        return
      end

      subroutine basm2g(stcprm, xin, xout)
C  receives the vector xin, giving the components of a point in the map
C  system; returns xout, the components of the same point in the geo system.
      dimension stcprm(15),xin(3),xout(3)
      ifind(l,k) = -2 + 3*l + k
        do k=1,3
          xout(k) = 0.
        enddo
        do l=1,3
          do k=1,3
            xout(l) = xout(l) + xin(k) * stcprm(ifind(k,l))
          enddo
        enddo
        return
      end
      subroutine cc2gll(stcprm, alat, along, ue,vn, ug,vg)
      parameter (sin1dg=.017452406437284)
      real stcprm(15)
      real norm
      call cpolll(stcprm, alat,along, enx,eny,enz)
      norm = sqrt(enx*enx + eny*eny)
      if (norm .le. SIN1DG) then
        call cgrnll(stcprm, alat,along, enx,eny,enz)
        norm = sqrt(enx*enx + eny*eny)
      endif
      enx = enx / norm
      eny = eny / norm
C* CHANGE MADE 3/9/99 TO ALLOW UG,VG AND UE,VN TO USE SAME STORAGE
      temp = enx * vn + eny * ue
      vg = eny * vn - enx * ue
      ug = temp
C* WITHOUT INTERFERING WITH OTHER
      return
      end
      subroutine cc2gxy(stcprm, x,y, ue,vn, ug,vg)
      parameter (sin1dg=.017452406437284)
      real stcprm(15)
      real norm
      call cpolxy(stcprm, x,y, enx,eny,enz)
      norm = sqrt(enx*enx + eny*eny)
      if (norm .le. sin1dg) then
        call cgrnxy(stcprm, x,y, enx,eny,enz)
        norm = sqrt(enx*enx + eny*eny)
      endif
      enx = enx/norm
      eny = eny/norm
C* CHANGE MADE 3/9/99 TO ALLOW UG,VG AND UE,VN TO USE SAME STORAGE
      temp = enx * vn + eny * ue
      vg = eny * vn - enx * ue
      ug = temp
C* WITHOUT INTERFERING WITH OTHER
      return
      end
      subroutine ccrvll(stcprm, alat, along, gx, gy)
      parameter (REARTH = 6371.2)
      real stcprm(15)
      real map(3),geog(3)
      call ll_geo(alat,along, geog)
      call basg2m(stcprm, geog, map)
      if ( abs(map(3)) .ge. 1.) then
        if (stcprm(1) .eq. map(3)) then
          gx = 0.
          gy = 0.
          return
        else
          xi = 0.
          eta = sign(.5e21,map(3)) / REARTH
        endif
      else
        fact = (stcprm(1) - map(3)) /
     C    (1. - map(3)) / (1. + map(3)) / REARTH
        xi = - map(2) * fact
        eta = map(1) * fact
        if (abs(stcprm(1)) .lt. 1.) then
          glambda = (stcprm(1) - 1.) * atan2(map(2), map(1))
          clambda = cos(glambda)
          slambda = sin(glambda)
          fact = xi * clambda - eta * slambda
          eta = xi * slambda + eta * clambda
          xi = fact
        endif
      endif
      gx = xi * stcprm(13) + eta * stcprm(14)
      gy = eta * stcprm(13) - xi * stcprm(14)
      return
      end
      subroutine ccrvxy(stcprm, x, y, gx, gy)
      parameter (REARTH = 6371.2)
      real stcprm(15)
      real map(3)
      call xy_map(stcprm, x,y, map)
      if (abs(map(3)) .ge. 1.) then
        if (stcprm(1) .eq. map(3)) then
          gx = 0.
          gy = 0.
        else
          xi = 0.
          eta = sign(.5e21,map(3)) / REARTH
        endif
      else
        xi0 = (x - stcprm(11)) * stcprm(15) / REARTH
        eta0 = (y - stcprm(12)) * stcprm(15) / REARTH
        xi = xi0 * stcprm(13) - eta0 * stcprm(14)
        eta = xi0 * stcprm(14) + eta0 * stcprm(13)
        xi = - stcprm(1) * xi
        eta = 1. - stcprm(1) * eta
        fact = sqrt((xi*xi + eta*eta) / (1. - map(3)) / (1. + map(3)))
        xi = xi/fact
        eta = eta/fact
        fact = (stcprm(1) - map(3)) / (1. - map(3)) / (1. + map(3)) 
     C     / REARTH
        xi = xi * fact
        eta = eta * fact
      endif
      gx = xi * stcprm(13) + eta * stcprm(14)
      gy = eta * stcprm(13) - xi * stcprm(14)
      return
      end
      subroutine cg2cll(stcprm, alat,along, ug,vg, ue,vn)
      parameter (sin1dg=.017452406437284)
      real stcprm(15)
      real norm
      call cpolll(stcprm, alat,along ,enx,eny,enz)
      if (norm .le. SIN1DG) then
        call cgrnll(stcprm, alat,along, enx,eny,enz)
        norm = sqrt(enx*enx + eny*eny)
      endif
      enx = enx / norm
      eny = eny / norm
C* CHANGE MADE 3/9/99 TO ALLOW UG,VG AND UE,VN TO USE SAME STORAGE
      temp = - enx * vg + eny * ug
      vn = eny * vg + enx * ug
      ue = temp
C* WITHOUT INTERFERING WITH EACH OTHER
      return
      end
      subroutine cg2cxy(stcprm, x,y, ug,vg, ue,vn)
      parameter (sin1dg=.017452406437284)
      real stcprm(15)
      real norm
      call cpolxy(stcprm, x,y, enx,eny,enz)
      norm = sqrt(enx*enx + eny*eny)
      if (norm .le. sin1dg) then
        call cgrnxy(stcprm, x,y, enx,eny,enz)
        norm = sqrt(enx*enx + eny*eny)
      endif
      enx = enx / norm
      eny = eny / norm
C* CHANGE MADE 3/9/99 TO ALLOW UG,VG AND UE,VN TO USE SAME STORAGE
      temp = - enx * vg + eny * ug
      vn = eny * vg + enx * ug
      ue = temp
C* WITHOUT INTERFERING WITH EACH OTHER
      return
      end
      subroutine cgrnll(stcprm, alat,along, enx,eny,enz)
      real stcprm(15)
      real map(3),pole(3),geog(3)
      call ll_geo(alat,along, geog)
      call basg2m(stcprm, geog, map)
      do k=1,3
        pole(k) = stcprm(3*k - 1)
      enddo
      call proj_3d(stcprm, map, pole, enx,eny,enz)
      return
      end
      subroutine cgrnxy(stcprm, x,y, enx,eny,enz)
      real stcprm(15)
      real map(3),pole(3)
      call xy_map(stcprm, x,y, map)
      do k=1,3
        pole(k) = stcprm(3*k - 1)
      enddo
      call proj_3d(stcprm, map, pole, enx,eny,enz)
      return
      end
      real function cgszll(stcprm, alat, along)
      real stcprm(15)
      real map(3)
      real geog(3)
      call ll_geo(alat, along, geog)
      call basg2m(stcprm, geog, map)
      if (map(3) .ge. 1.) then
        if (stcprm(1) .ge. 1.) then
          cgszll = 2. * stcprm(15)
        else
          cgszll = 0.
        endif
      else if (map(3) .le. -1.) then
        if (stcprm(1) .le. -1.) then
          cgszll = 2. * stcprm(15)
        else
          cgszll = 0.
        endif
      else if (abs(stcprm(1)) .ge. 1.) then
        cgszll = stcprm(15) * (1. + sign(map(3),stcprm(1)))
      else
        ymerc = -.5 * log( (1. - map(3))/(1. + map(3)) )
        cgszll = stcprm(15) * exp( - (1.-stcprm(1)) * ymerc) * 
     C   (1. + map(3))
      endif
      return
      end
      function cgszxy(stcprm, x, y)
      real stcprm(15)
      real map(3)
      call xy_map(stcprm, x,y,map)
      if (map(3) .ge. 1.)  then
        if (stcprm(1) .ge. 1.) then
          cgszxy = 2.*stcprm(15)
        else
          cgszxy = 0.
        endif
      else if (map(3) .le. -1.) then
        if (stcprm(1) .le. -1.) then
          cgszxy = 2.*stcprm(15)
        else
          cgszxy = 0.
        endif
      else if (abs(stcprm(1)) .ge. 1.) then
        cgszxy = stcprm(15) * (1. + sign(map(3),stcprm(1)))
      else
        ymerc = -.5 * log( (1. - map(3))/(1. + map(3)) )
        cgszxy = stcprm(15) * exp( - (1. - stcprm(1)) * ymerc) *
     C      (1. + map(3))
      endif
      return
      end
      subroutine cll2xy(stcprm, alat, along, x, y )
      real stcprm(15)
      real geog(3),temp(3)
      call ll_geo(alat,along, geog)
      call basg2m(stcprm, geog, temp)
      call map_xy(stcprm, temp, x,y)
      return
      end
      integer function cml2xy(stcprm, alat, along, lsiz, x, y )
      real stcprm(15)
      real alat(*),along(*),x(*),y(*)
      real geog(3), mapold(3), mapnew(3), mapmid(3)
      logical cross
      cross(a,b) = ((a .le. 0.) .and. (b .gt. 0.))
     c .or. ((b.le.0.) .and. (a .gt. 0.))
C* Given a string of latitudes and longitudes, returns a string
C* of x-y values.  In the event that the *first* segment crosses the
C* 'cut' of the map, from side a to side b, the first x-y returned
C* will be the side b version of the cut-crossing point.  (If size_l<0,
C* then only the first part of the segment, on side a, will be returned.)
C*  If *any other* segment crosses the cut, this ends the transfer, the
C* last x-y value will be the near side of cut, and the returned value
C* will be negative.  Returns an int, the absolute value of which is
C* the number of points transformed, and the sign of which indicates
C* (if negative) that more points remain to be transformed.*/
      if (lsiz .eq. 0) then
C  No points provided; none returned.
        cml2xy = 0
        return
      endif
      kount=1
      npoints = iabs(lsiz)
      call ll_geo(alat(1),along(1), geog)
      call basg2m(stcprm, geog, mapold)
      call map_xy(stcprm, mapold, x(1),y(1))
      if (npoints .eq. 1) then
C* Only one point provided; one point returned.
        cml2xy = 1
        return
      endif
      if (lsiz .gt. 0) then
C* If first segment crossed, truncate first section.
        call ll_geo(alat(2),along(2), geog)
        call basg2m(stcprm, geog, mapnew)
        if (cross(mapold(2),mapnew(2))) then
          fact = mapnew(2) / (mapnew(2) - mapold(2))
          mapmid(2)=0.
          mapmid(1) = mapnew(1) + fact * (mapold(1) - mapnew(1))
          mapmid(3) = mapnew(3) + fact * (mapold(3) - mapnew(3))
          if (mapmid(1) .le. 0. ) then
C* only in this case is there a crossing of the cut.
            sizmid = sqrt(mapmid(1)**2 + mapmid(3)**2)
            if (sizmid .gt. 0.) then
              mapmid(1) = mapmid(1) / sizmid
              mapmid(3) = mapmid(3) / sizmid
            else
              mapmid(1) = -1.
            endif
            if (mapnew(2) .gt. 0.) then
              call map_xe(stcprm, mapmid, xi,eta, 'e')
            else 
              call map_xe(stcprm, mapmid, xi,eta, 'w')
            endif
            call xe_xy(stcprm, xi,eta, x(1),y(1))
          endif
        endif
        call map_xy(stcprm, mapnew, x(2), y(2))
        do k=1,3
          mapold(k) = mapnew(k)
        enddo
        kount=kount+1
      endif
C* geog values indicate v[0] "Toward lat0,lon0", v[1] East @0,0,
C* v[2] North @0,0.  Thus cut half-plane is v[1]=0, v[0]<=0.
C* Interpolate.*/
      do kount=kount+1,npoints
        call ll_geo(alat(kount),along(kount),geog)
        call basg2m(stcprm, geog, mapnew)
        if ( cross(mapold(2),mapnew(2)) ) then
C* Place interpolated value in location count and exit*/
C*If cut crossed, replace xy-point in location 0 with interpolated value*/
          fact = mapnew(2) / (mapnew(2) - mapold(2))
          mapmid(2)=0.
          mapmid(1) = mapnew(1) + fact * (mapold(1) - mapnew(1))
          mapmid(3) = mapnew(3) + fact * (mapold(3) - mapnew(3))
          if (mapmid(1) .le. 0. ) then
C* only in this case is there a crossing of the cut.
            sizmid = sqrt(mapmid(1)**2 + mapmid(3)**2)
            if (sizmid .gt. 0.) then
              mapmid(1) = mapmid(1) / sizmid
              mapmid(3) = mapmid(3) / sizmid
            else
              mapmid(1) = -1.
            endif
          endif
          if (mapold(2) .gt. 0.) then
            call map_xe(stcprm, mapmid, xi,eta, 'e')
          else
            call map_xe(stcprm, mapmid, xi,eta, 'w')
          endif
          call xe_xy(stcprm, xi,eta, x(kount),y(kount))
          ml2xy = -kount
          return
        endif
        call map_xy(stcprm, mapnew, x(kount),y(kount))
        do k=1,3
          mapold(k) = mapnew(k)
        enddo
      enddo
      ml2xy = npoints
      return
      end
      subroutine cpolll(stcprm, alat,along, enx,eny,enz)
      real stcprm(15)
      real map(3),pole(3),geog(3)
      call ll_geo(alat,along, geog)
      call basg2m(stcprm, geog, map)
      do k=1,3
        pole(k) = stcprm(3*k + 1)
      enddo
      call proj_3d(stcprm, map, pole, enx,eny,enz)
      return
      end
      subroutine cpolxy(stcprm, x,y, enx,eny,enz)
      real stcprm(15)
      real map(3),pole(3)
      call xy_map(stcprm, x,y, map)
      do k=1,3
        pole(k) = stcprm(3*k + 1)
      enddo
      call proj_3d(stcprm, map, pole, enx,eny,enz)
      return
      end
      subroutine cxy2ll(stcprm , x, y, alat, along)
      real stcprm(15)
      real map(3),geog(3)
      call xy_map(stcprm, x,y, map)
      call basm2g(stcprm, map, geog)
      call geo_ll(geog, alat,along)
      return
      end
      subroutine prmprt(stcprm)
      real stcprm(15)
      ifind(l,k) = -2 + 3*l +k
      write(*,*)' Gamma = ',stcprm(1)
      do k=1,3
        write(*,'('' line '',i1,'':  '',$)')k-1
        write(*,'(3f10.6)') (stcprm(ifind(l,k)),l=1,3)
      enddo
      write(*,*)'x0 = ',stcprm(11),', y0 = ',stcprm(12)
      write(*,*)'cos(th) = ',stcprm(13),', sin(th) = ',stcprm(14)
      write(*,*)'gridsize = ',stcprm(15)
      write(*,*)
      return
      end

      integer function ans_pt(prompt,values)
      integer scan
      character * (*) prompt,values
      character*80 instr
        write(*,'(1x,a,$)') prompt
        read(*,'(a)') instr
        ans_pt = (scan(values,instr)+1)/2
      return
      end

C  void main() {
      parameter (REARTH=6371.2)
      real stcprm(15)
      integer ano_pt,ans_pt,choice,ano_pj
      alat = 0.
      ano_pj = 1
      do  while (ano_pj .gt. 0)
        write(*,*) 'O - Oblique Stereographic'
        write(*,*) 'L - Lambert Polar'
        write(*,*) 'T - Transverse Mercator'
        write(*,*) 'U - Oblique Mercator'
        write(*,*) 'C - Oblique Lambert'
        choice = ans_pt('Enter Choice: ','OoLlTtUuCc')
        if (choice .eq. 1) then
          write(*,'('' Enter Latitude and Longitude: '',$)')
          read(*,*)alat,along
          call sobstr(stcprm,alat,along)
        else if (choice .eq. 2) then
          write(*,'('' Enter Two Ref. Latitudes and One Longitude: '',
     C            $)')
          read(*,*)alat1,alat2,along
          write(*,*)'Eqvlat = ',eqvlat(alat1,alat2)
          call stlmbr(stcprm,eqvlat(alat1,alat2),along)
        else if (choice .eq. 3) then
          write(*,'('' Enter Latitude and Longitude: '',$)')
          read(*,*)alat,along
          call stvmrc(stcprm,alat,along)
        else if (choice .eq. 4) then
          write(*,'(a,$)') 'Enter Central latitude and longitude: '
          read(*,*) alat1,along1
          write(*,'(a,$)') 'Enter Secondary latitude and longitude: '
          read(*,*) alat2,along2
          call sobmrc(stcprm, alat1, along1, alat2, along2)
        else if (choice .eq. 5) then
          write(*,'(a,$)') 'Enter Central latitude and longitude: '
          read(*,*) alat,along
          write(*,'(a,$)') 'Enter lat & lon of second point on circle: '
          read(*,*) alat1,along1
          write(*,'(a,$)') 'Enter lat & lon of third point on circle: '
          read(*,*) alat2,along2
          call soblmb(stcprm, alat1,along1, alat,along, alat2,along2)
        endif
        choice = ans_pt('1-point or 2-point scaling? ','1o2t')
        if (choice .eq. 1) then
          write(*,'(a,$)') 'Enter x,y of anchor point: '
          read(*,*) x,y
          write(*,'(a,$)') 'Enter lat,long of anchor point: '
          read(*,*) alat,along
          write(*,'(a,$)') 'Enter lat,long of reference point: '
          read (*,*) rlat,rlong
          write(*,'(a,$)') 'Enter grid size at reference point: '
          read(*,*) grdsiz
          write(*,'(a,$)')
     C       'Enter y-axis orientation at reference point: '
          read(*,*) orient
          call stcm1p(stcprm, x,y, alat,along,
     C        rlat,rlong, grdsiz,orient)
        else
          write(*,'(a,$)') 'Enter x,y of first anchor point: '
          read(*,*) x1,y1
          write(*,'(a,$)') 'Enter lat,long of first anchor point: '
          read(*,*) alat1,along1
          write(*,'(a,$)') 'Enter x,y of second anchor point: '
          read(*,*) x2,y2
          write(*,'(a,$)') 'Enter lat,long of second anchor point: '
          read(*,*) alat2,along2
          call stcm2p(stcprm,x1,y1,alat1,along1, x2,y2,alat2,along2)
        endif
        call prmprt(stcprm)
        ano_pt = 1
        do while (ano_pt .ne. 0)
          ano_pt = ans_pt('Translate x,y point? (y/n) ','yY')
          if (ano_pt .ne. 0) then
            write(*,'(1x,a,$)') 'Enter x,y: '
            read(*,*) x,y
            call cxy2ll(stcprm, x,y, alat,along)
            call cll2xy(stcprm, alat,along, oldx,oldy)
            write(*,*) 'x,y = (',oldx,',',oldy,'), lat,long = (',alat,
     C         ',',along,').'
            write(*,*) 'gridsize(x,y) = ', cgszxy(stcprm,x,y),
     C         'gridsize(l,l) = ', cgszll(stcprm,alat,along)
            call cpolxy(stcprm, x,y, enx,eny,enz)
            write(*,*) 'Polar axis from x,y = (',
     C         enx,',',eny,',',enz,')'
            call cpolll(stcprm, alat,along, enx,eny,enz)
            write(*,*) 'Polar axis from lat,long = (',
     C         enx,',',eny,',',enz,')'
            call cgrnxy(stcprm, x,y, enx,eny,enz)
            write(*,*) 'Greenwich axis from x,y = (',
     C         enx,',',eny,',',enz,')'
            call cgrnll(stcprm, alat,along, enx,eny,enz)
            write(*,*) 'Greenwich axis from lat,long = (',
     C            enx,',',eny,',',enz,')'
            call cc2gxy(stcprm, x,y, 0.,10., ug,vg)
            call cg2cxy(stcprm, x,y, ug,vg, ue,vn)
            write(*,*) 'x,y winds from (E,N) to (Ug,Vg):(',ue,',',vn,
     C         ') to (',ug,',',vg,')'
            call cc2gll(stcprm, alat,along, 0.,10., ug,vg)
            call cg2cll(stcprm, alat,along, ug,vg, ue,vn)
            write(*,*) 'l,l winds from (E,N) to (Ug,Vg):(',ue,',',vn,
     C            ') to (',ug,',',vg,')'
            call ccrvxy(stcprm, x,y, gx,gy)
            write(*,*) 'x,y curvature vector (gx,gy):(',gx,',',gy,')'
            call ccrvll(stcprm, alat,along, gx,gy)
            write(*,*) 'lat,long curvature vector (gx,gy):('
     C         ,gx,',',gy,')'
           endif
         enddo
      ano_pj = ans_pt('Another Projection? (y/n) ','yY')
      enddo
      stop
      end
      function eqvlat(lat1,lat2)
C*  Written 12/21/94 by Dr. Albion Taylor
C*
C*    This function is provided to assist in finding the tangent latitude
C*    equivalent to the 2-reference latitude specification in the legend
C*    of most lambert conformal maps.  If the map specifies "scale
C*    1:xxxxx true at 40N and 60N", then eqvlat(40.,60.) will return the
C*    equivalent tangent latitude.
C*  INPUTS:
C*    lat1, lat2:  The two latitudes specified in the map legend
C*  RETURNS:
C*    the equivalent tangent latitude
C*  EXAMPLE:  stcmap(& strcmp, eqvlat(40.,60.), 90.)
C*/
      real lat1,lat2
      parameter (pi=3.14159265358979323846,radpdg=pi/180.)
      parameter (dgprad=180./pi)
      parameter (fsm = 1.e-3)
        slat1 = sin(radpdg * lat1)
        slat2 = sin(radpdg * lat2)
        if (slat1 + slat2 .eq. 0.) then
          eqvlat = 0.
          return
        endif
        if (abs(slat1) .ge. 1.) then
          eqvlat = sign(90.,slat1)
          return
        endif
        if (abs(slat2) .ge. 1.) then
          eqvlat = sign(90.,slat2)
        endif
        if ( abs(slat1 - slat2) .gt. fsm) then
          al1 = log( (1. - slat1) / (1. - slat2) )
          al2 = log( (1. + slat1) / (1. + slat2) )
        else
          tau = - (slat1 - slat2)/(2. - slat1 - slat2)
          tau = tau * tau
          al1 = 2./(2. - slat1 - slat2) * (1.    + tau *
     a                                    (1./3. + tau *
     b                                    (1./5. + tau *
     c                                    (1./7. + tau))))
          tau = (slat1 - slat2)/(2. + slat1 + slat2)
          tau = tau * tau
          al2 = -2./(2. + slat1 + slat2) * (1.    + tau *
     a                                     (1./3. + tau *
     b                                     (1./5. + tau *
     c                                     (1./7. + tau))))
          endif
          eqvlat = asin ((al1 + al2) / (al1 - al2)) * DGPRAD
          return
        end

      subroutine geo_ll(geog,plat,plon)
C  given a vector geog in the geo (geographic) coordinate system,
C  returns as plat and plon the latitude and longitude of the 
C  corresponding point.
C
C  The geo system is a 3-D Cartesian coordinate system, with origin
C  at the Center of the Earth, with the x_3 axis pointing to the North
C  Pole, the x_1 axis pointing to the Equator at the Greenwich Meridian,
C  and the x_2 axis by the right hand rule pointing to the Equator
C  at 90 E.
      parameter (pi=3.14159265358979323846,radpdg=pi/180.)
      parameter (dgprad=180./pi)
      dimension geog(3)
      fact = geog(1)*geog(1) + geog(2)*geog(2)
      if (fact .le. 0.) then
        plon = 0.
        plat = sign(90.,geog(3))
      else
        fact = sqrt(fact)
        plat = DGPRAD * atan2(geog(3),fact)
        plon = DGPRAD * atan2(geog(2),geog(1))
      endif
      return
      end

      function csabva(a, b)
C* returns (1. - cos(a * b) ) / a / a, or its limit as a -> 0
        term = snabva(.5 * a, b)
        csabva = .5 * term * term
        return
      end

      real function lnabva(a, b)
C* returns ln (1. + a * b) / a, or its limit as a -> 0.
        c = a * b
        if (abs(c) .gt. 0.01) then
          lnabva = log( 1. + c ) / a
        else
C* using (1. + t) / (1. - t) = 1. + c when t = c / (2. + c)
          t = c/(2. + c)
          t = t * t
          lnabva =  b / (2. + c) * (2.    + t *
     a                             (2./3. + t *
     b                             (2./5. + t *
     c                             (2./7. ))))
        endif
        return
      end

      function snabva(a, b)
C* returns sin(a * b) / a, or its limit as a -> 0.
        c = a * b
        csq = c * c
        if (csq .gt. .001) then
          snabva = sin( c )  / a
        else
          snabva = b * ( 1. - csq / 6. *
     a                 ( 1. - csq / 20. *
     b                 ( 1. - csq / 42. )))
        endif
        return
      end

C  A series of functions which may have problems in certain ranges
C  of one of their parameters.  These functions, if written normally,
C  may suffer round-off problems.  They are used in particular for
C  Lambert Conformal projections whose parameters make them almost
C  the same as a Mercator Projection.

      function xpabva(a, b)
C* returns (exp (a * b) - 1. ) / a, or its limit as a -> 0.
        c = a * b
        csq = .25 * c * c
        if ( csq  .gt. .001) then
          xpabva = (exp( c ) - 1.) / a
        else
          xpabva = exp(.5 * c) * b * (1. + csq / 6. *
     a                               (1. + csq / 20. *
     b                               (1. + csq / 42. )))
        endif
        return
      end

      function atnabv(a, b, c)
C* returns atan2(a*b,1-a*c)/a, or its limit as a -> 0
        xi = a * b
        eta = a * c
        vsq = xi * xi + eta * eta
        if ( vsq - 2. * eta + 1. .le. 0.) then
          atnabv = 0.
          return
        endif
        if (abs(xi) .gt. .01 * (1. - eta) ) then
          atnabv = atan2(xi, 1. - eta) / a
        else
          t = xi/(1. - eta)
          t = t * t
          atnabv = b / (1. - eta) * (1.    - t *
     a                              (1./3. - t *
     b                              (1./5. - t *
     c                              (1./7. ))))
        endif
        return
      end
      subroutine ll_geo(xlat, xlong, vector)
C  Given a latitude xlat and longitude xlong, returns the unit vector
C  directed to the given point in the geo (geographic) coordinate system.
C
C  The geo system is a 3-D Cartesian coordinate system, with origin
C  at the Center of the Earth, with the x_3 axis pointing to the North
C  Pole, the x_1 axis pointing to the Equator at the Greenwich Meridian,
C  and the x_2 axis by the right hand rule pointing to the Equator
C  at 90 E.
      parameter (pi=3.14159265358979323846,radpdg=pi/180.)
      parameter (dgprad=180./pi)
      dimension vector(3)
        vector(3) = sin(RADPDG * xlat)
        clat = cos(RADPDG * xlat)
        vector(1) = clat * cos(RADPDG * xlong)
        vector(2) = clat * sin(RADPDG * xlong)
      return
      end
      subroutine map_xe(stcprm, x_map, xi, eta,flag)
      parameter (pi=3.14159265358979323846,pi_2=pi/2.)
      character * 1 flag
      dimension stcprm(15),x_map(3)
        if (abs(x_map(3)) .ge. 1.) then
C* Projection pole or its antipodes
            xi = 0.
          if (stcprm(1) * x_map(3) .gt. 0.) then
            eta = 1./stcprm(1)
          else
            if (x_map(3) .gt. 0) then
              eta = 1.0e10
            else
              eta = -1.0e10
            endif
          endif
          return
        else
          if (abs(stcprm(1)) .eq. 1.) then
C* Stereographic Case, away from pole
            fact = 1. / (stcprm(1) + x_map(3) )
            xi = x_map(2)*fact
            eta = 1./stcprm(1) - x_map(1) * fact
          else
            ymerc = .5 * log( ( 1. + x_map(3) ) / (1. - x_map(3) ) )
C* This Projection has a cut.  Select theta according to the rule
C* If cutflag = 'e', -PI/2 <= theta < 3PI/2, if cutflag = 'E',
C* 0 <= theta < 2PI, if cutflag = 'w', -3Pi/2 <= theta < PI/2,
C* cutflag = 'W', -2PI <= theta < 0., else -PI<=theta<PI. 
          if (flag .eq. 'E') then
            theta = atan2(- x_map(2), - x_map(1) ) + pi
          else if (flag .eq. 'W') then
            theta = atan2(- x_map(2), - x_map(1) ) - pi
          else if (flag .eq. 'e') then
            theta = atan2(- x_map(1), x_map(2) )  + pi_2
          else if (flag .eq. 'w') then
            theta = atan2(x_map(1), - x_map(2) )  - pi_2
          else
            theta = atan2(x_map(2), x_map(1) )
          endif
            rhog = xpabva( stcprm(1), - ymerc)
C* rhog = ( exp( - gamma * ymerc ) - 1. ) / gamma */
            xi = (1. + stcprm(1) * rhog) * snabva(stcprm(1), theta)
            eta = stcprm(1) * (1. + stcprm(1) * rhog) *
     A                        csabva(stcprm(1), theta) - rhog
          end if
        end if
        return
      end
      subroutine map_xy(stcprm, x_map, x, y)
      dimension stcprm(15),x_map(3)
        call map_xe(stcprm, x_map, xi, eta, 'c')
        call xe_xy(stcprm, xi, eta, x, y)
        return
      end
      subroutine mpstrt(stcprm, conang, p_lat, p_long, r_lat, r_long)
      parameter (pi=3.14159265358979323846,radpdg=pi/180.)
      parameter (dgprad=180./pi)
      parameter (REARTH=6371.2)
      dimension stcprm(15)
Ctypedef struct {
C 1:     gamma,
C 2-10:  rotate(3,3),
C 11-12: x0,y0,
C 13-14: crotate,srotate,
C  15:   grdszq
C

C*
C*  General Purpose routine to Set Map Parameters.  Called by the
C*  special purpose routines for oblique stereographic, oblique and
C*  transverse Mercator, oblique Lambert Conformal, etc. Projections
C*  Inputs: p_lat,p_long - Latitude and Longitude of the Pole
C*            Point for the Projection
C*          r_lat,r_long - Latitude and Longitude of a
C*            reference point - 180degrees around the Pole Point
C*            from the cut
C*          conang - angle between the Projection Pole axis and
C*            the generators (sides) of the cone on which the Earth
C*            is projected. + or - 90 degrees indicates a plane,
C*            hence Stereographic; 0 degrees indicates a cylinder
C*            and hence Mercator.
C*  Outputs: stcprm - map parameters
C*/
      dimension temp(3)
      ifind(l,k) = -2 + 3*l + k
        call ll_geo(p_lat,p_long,temp)
        do k=1,3
          stcprm(ifind(3,k)) = temp(k)
        enddo
C  for (k=0;k<3;k++) stcprm->rotate[2][k] = temp.v[k];
        call ll_geo(r_lat,r_long,temp)
        do k=1,3
          stcprm(ifind(1,k)) = temp(k)
        enddo
C  for (k=0;k<3;k++) stcprm->rotate[0][k] = temp.v[k];
        tnorm = x_prod(stcprm(ifind(3,1)),stcprm(ifind(1,1)),
     C     stcprm(ifind(2,1)))
C  norm =
C  x_product (stcprm->rotate[2],stcprm->rotate[0],stcprm->rotate[1]);
        do k=1,3
          stcprm(ifind(2,k)) = stcprm(ifind(2,k)) /tnorm
        enddo
C  for (k=0;k<3;k++) stcprm->rotate[1][k] /= norm
        tnorm = x_prod(stcprm(ifind(2,1)),stcprm(ifind(3,1)),
     C     stcprm(ifind(1,1)))
C  x_product (stcprm->rotate[1],stcprm->rotate[2],stcprm->rotate[0]);
C  stcprm->x0 = stcprm->y0 = stcprm->srotate = 0;
        stcprm(11)=0.
        stcprm(12)=0.
        stcprm(14)=0.
        stcprm(13)=1.
        stcprm(15) = REARTH
        stcprm(1) = sin(RADPDG * conang )
C stcprm->crotate = 1.;
C stcprm->gridszeq = REARTH;
C stcprm->gamma = sin(RADPDEG * cone_ang);
C*geographic triple : i = equator @ std merid,
C       j = equator @ 90E,
C       k = North Pole */
C*map base triple : i' = M-prime meridian @ M-equator,
C     j' = M-East at location i',
C     k' = M-Pole */
C
      return
      end

      subroutine proj_3d(stcprm, point, vect, enx,eny,enz)
C/*
C *  At a given point, resolves the components of a vector vect in the
C *  local coordinate system of the map projection.  It is assumed that
C *  vect and point is given in the 3-dimensional geocentric _map_
C *  coordinate system, rather than the North centered _geo_ system.
C *  returns these components as enx, eny, and enz.  Replaces vect with
C *  its projection on the tangent plane at point.
C */
      real stcprm(15)
      real point(3),vect(3)
      dot_pr = 0.
      do k=1,3
        dot_pr = dot_pr + point(k)*vect(k)
      enddo
C/*
C *  dot_prod is the local vertical component of vect.  Note, point is
C *  assumed to be a unit vector.
C */
      enz = dot_pr
      do k=1,3
        vect(k) = vect(k) - dot_pr * point(k)
      enddo
C/*  vect has now been projected to a plane tangent to point.
C */
      fact = 1. + point(3)
      xi = vect(1) - vect(3) * point(1) / fact
      eta = vect(2) - vect(3) * point(2) / fact
C/*
C *  xi, eta represent components of vector on the projection plane,
C *  i.e. the 2-dimensional map system.
C */
      fact = stcprm(1) -1.
      if (fact .lt. 0.) then
C/*
C *  For Lambert Conformal and Mercator projections (gamma < 1.0) ,
C *  a rotation of the vector components is needed.
C */
        if ( abs(point(3)) .lt. 1.) then
          theta = fact * atan2(point(2),point(1))
          cgthta = cos(theta)
          sgthta = sin(theta)
          fact = xi * cgthta - eta * sgthta
          eta = eta * cgthta + xi * sgthta
          xi = fact
        endif
      endif
C/*
C *  Now rotate xi, eta to the final map projection direction.
C */
      enx = eta * stcprm(13) - xi * stcprm(14)
      eny = - xi * stcprm(13) - eta * stcprm(14)
      return
      end
      integer function scan(char,charset)
      character * (*) char,charset
      ilen1 =len(char)
      ilen2 = len(charset)
      do i = 1,ilen1
        if (index(charset,char(i:i)) .ne. 0) then
          scan = i
          return
        endif
      enddo
      scan = 0
      return
      end
      subroutine soblmb(stcprm, r_lat, r_lng, a_lat1, a_lng1,
     C    a_lat2, a_lng2)
      parameter (pi=3.14159265358979323846,radpdg=pi/180.)
      parameter (dgprad=180./pi)
      real a_lat1,a_lng1,r_lat,r_lng,a_lat2,a_lng2
      real stcprm(15)
C/*
C *  Set Map Parameters for an Oblique LaMBeRt Conic Conformal
C *          Projection.
C *  Inputs: r_lat,r_lng,
C *          a_lat1,a_lng1,
C *          a_lat2,a_lng2  - latitudes and longitudes of three
C *              points on the circle (great or small) which is
C *              tangent to the projection cone.  Point r_lat,r_lng is
C *              the reference point; 180degrees away from the cut.
C *  Outputs: stcprm - map parameters
C */
      real pnt_a1(3),pnt_rf(3),pnt_a2(3),pol_pt(3)
      real b_a(3),c_b(3),norm,sin_cn,cos_cn,temp,lat_pl,lng_pl
      call ll_geo(a_lat1,a_lng1,pnt_a1)
      call ll_geo(r_lat,r_lng,pnt_rf)
      call ll_geo(a_lat2,a_lng2,pnt_a2)
      do k=1,3
        b_a(k) = pnt_rf(k) - pnt_a1(k)
        c_b(k) = pnt_a2(k) - pnt_rf(k)
      enddo
      norm = x_prod(b_a,c_b,pol_pt)
      if (norm .ne. 0.) then
        call geo_ll(pol_pt,lat_pl,lng_pl)
        sin_cn = 0.
        do k=1,3
          pol_pt(k) = pol_pt(k)/norm
          sin_cn = sin_cn + pnt_rf(k) * pol_pt(k)
        enddo
        cos_cn = 0.
        do k=1,3
          temp = sin_cn *pol_pt(k) - pnt_rf(k)
          cos_cn = cos_cn + temp*temp
        enddo
        call mpstrt(stcprm,dgprad*atan2(sin_cn,cos_cn),lat_pl,lng_pl,
     C        r_lat,r_lng)
      else
        temp = 0.
        do k=1,3
          temp = temp + b_a(k) * b_a(k)
        enddo
        if (temp .ne. 0.) then
          call sobmrc(stcprm,r_lat,r_lng,a_lat1,a_lng1)
        else
          do k=1,3
            temp = temp + c_b(k)*c_b(k)
          enddo
          if (temp .ne. 0.) then
            call sobmrc(stcprm,r_lat,r_lng, a_lat1,a_lng1)
          else
            call sobstr(stcprm, r_lat,r_lng)
          endif
        endif
      endif
      return
      end
      subroutine sobmrc(stcprm, r_lat, r_long, an_lat, an_lng)
C/*
C *  Set Map Parameters for an Oblique MeRCator Projection
C *  Inputs: r_lat,r_long - latitude and longitude of reference
C *            point on the great circle tangent to the cylinder,
C *            and 180 degrees from the cut.
C *       an_lat,an_lng - latitude and longitude of
C *            an anchor point elsewhere on the same great circle.
C *  Outputs: stcprm - map parameters
C */
      real refpnt(3),anchpt(3),prjpol(3)
      real stcprm(15)
      real*8  p_lat,p_long,norm
      call ll_geo(r_lat,r_long,refpnt)
      call ll_geo(an_lat,an_lng,anchpt)
      norm = x_prod(anchpt,refpnt,prjpol)
      if (norm .ne. 0.) then
        call geo_ll(prjpol,p_lat,p_long)
        call mpstrt(stcprm, 0., p_lat,p_long, r_lat,r_long)
      else
        call stvmrc(stcprm, r_lat, r_long)
      endif
      return
      end

      subroutine sobstr(stcprm, p_lat, p_lon)
C*
C*  Set Map Parameters for an OBlique STeReographic Projection
C*  Inputs: p_lat,P_lon - latitude and longitude of the
C*          projection pole (tangent point)
C*  Outputs: stcprm - map parameters
C*/
      dimension stcprm(15)
      call mpstrt(stcprm, 90., p_lat,p_lon, p_lat-90.,p_lon)
      return
      end
      subroutine stcm1p(stcprm, x1, y1, xlat1, xlong1,
     C xlatrf, xlonrf, gridsz, orient)
      parameter (pi=3.14159265358979323846,radpdg=pi/180.)
      parameter (dgprad=180./pi)
      real stcprm(15)
      real enx,eny,enz,norm,x1a,y1a
      c_or = cos(RADPDG * orient)
      s_or = - sin(RADPDG * orient)
C  stcprm->x0 = stcprm->y0 = stcprm->srotate = 0;
      stcprm(11) = 0.
      stcprm(12) = 0.
      stcprm(13) = 1.
      stcprm(14) = 0.
      stcprm(15) = 1.
C  stcprm->crotate = stcprm -> gridszeq = 1.0;
      call cpolll(stcprm, xlatrf, xlonrf, enx, eny, enz)
      norm = sqrt(enx*enx + eny*eny)
      if (norm .eq. 0.) then
        call cgrnll(stcprm,xlatrf,xlonrf,enx,eny,enz)
        norm = sqrt (enx* enx + eny*eny)
      endif
      enx = enx/norm
      eny = eny/norm
C  stcprm->gridszeq *= gridsz / cgszll(stcprm, xlatrf,xlonrf);
      stcprm(15) = stcprm(15) * gridsz / cgszll(stcprm,xlatrf,xlonrf)
C  stcprm -> crotate = eny * c_or - enx * s_or;
      stcprm(13) = eny * c_or - enx * s_or
C  stcprm -> srotate = -eny * s_or - enx * c_or;
      stcprm(14) = -eny * s_or - enx * c_or
      call cll2xy(stcprm, xlat1,xlong1, x1a,y1a)
C  stcprm->x0 += x1 - x1a;
C  stcprm->y0 += y1 - y1a;
      stcprm(11) = stcprm(11) + x1 - x1a
      stcprm(12) = stcprm(12) + y1 - y1a
      return
      end
      subroutine stcm2p(stcprm, x1, y1, xlat1, xlong1,
     C  x2, y2, xlat2, xlong2)
      real stcprm(15)
      real x1a,y1a, x2a,y2a, den,dena
C  stcprm->x0 = stcprm->y0 = stcprm->srotate = 0;
      stcprm(11) = 0.
      stcprm(12) = 0.
      stcprm(14) = 0.
C  stcprm->crotate = stcprm->gridszeq = 1.;
      stcprm(13) = 1.
      stcprm(15) = 1.
      call cll2xy(stcprm, xlat1,xlong1, x1a,y1a)
      call cll2xy(stcprm, xlat2,xlong2, x2a,y2a)
      den = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2))
      dena = sqrt((x1a - x2a)*(x1a - x2a) + (y1a - y2a)*(y1a - y2a))
      stcprm(13) = ((x1a - x2a)*(x1 - x2) + (y1a - y2a)*(y1 - y2) )
     C     /den /dena
      stcprm(14) = ((y1a - y2a)*(x1 - x2) - (x1a - x2a)*(y1 - y2) )
     C     /den /dena
C  stcprm->gridszeq *= dena /den;
      stcprm(15) = stcprm(15) * dena / den
      call cll2xy(stcprm, xlat1,xlong1, x1a,y1a)
C  stcprm->x0 += x1 - x1a;
C  stcprm->y0 += y1 - y1a;
      stcprm(11) = stcprm(11) + x1 - x1a
      stcprm(12) = stcprm(12) + y1 - y1a
      return
      end
      subroutine stlmbr(stcprm, reflat, reflon)
C*
C*  Set Map Parameters for a North Polar LaMBeRt Conic Conformal
C*    Projection
C*  Inputs: reflat - tangent latitude of the tangent latitude of
C*                   cone.
C*          reflon - midrange longitude (180 degrees from cut)
C*  Outputs: stcprm - map parameters
C*/
      real stcprm(15)
      CALL mpstrt(stcprm,reflat,90.,reflon,0.,reflon)
      return
      end

      subroutine stcmap(stcprm, tnglat, reflon)
C*
C*  Set Map Parameters for a North Polar LaMBeRt Conic Conformal
C*    Projection
C*  included for compatibliity with previous version
C*  Inputs: tnglat - tangent latitude of the tangent latitude of
C*                   cone.
C*          reflon - midrange longitude (180 degrees from cut)
C*  Outputs: stcprm - map parameters
C*/
      real stcprm(15)
      CALL mpstrt(stcprm,tnglat,90.,reflon,0.,reflon)
      return
      end
      subroutine stvmrc(stcprm, r_lat, r_long)
      real stcprm(15)
C*
C*  Set Map Parameters for a TransVerse MeRCator Projection
C*  Inputs: r_lat,r_longit - latitude and longitude of the reference
C*            point on the meridian tangent to the cylinder,
C*            and 180 degrees from the cut.
C*  Outputs: stcprm - map parameters
C*/
      real refpnt(3),anchrp(3),projpl(3)
      call ll_geo(r_lat, r_long, refpnt)
      call ll_geo(r_lat - 90.,r_long,anchrp)
      anorm = x_prod(anchrp, refpnt, projpl)
      call geo_ll(projpl, p_lat, p_long)
      call mpstrt(stcprm, 0., p_lat,p_long, r_lat, r_long)
      return
      end

      function x_prod(A,B,C)
C  Returns as C the vector cross product of A and B:
C  C = A x B.  Also returns as function value the absolute value of C.
      dimension a(3),b(3),C(3),indx(4)
      data indx/2,3,1,2/
        D = 0.
        do k=1,3
          C(k) = A(indx(k))*B(indx(k+1)) - A(indx(k+1))*B(indx(k))
          D = D + C(k) * C(k)
        enddo
        x_prod = SQRT(D)
        return
      end

      subroutine xe_xy(stcprm, xi, eta, x, y)
      parameter (REARTH=6371.2)
      dimension stcprm(15)
        x = stcprm(11) + REARTH / stcprm(15) *
     A  (stcprm(13) * xi + stcprm(14) * eta)
        y = stcprm(12) + REARTH / stcprm(15) *
     B  (stcprm(13) * eta - stcprm(14) * xi)
        return
      end
      subroutine xy_map(stcprm, x,y, map)
      real stcprm(15)
      real map(3)
      real lnabva
      parameter (REARTH=6371.2)
        xi0 = (x - stcprm(11)) * stcprm(15) / REARTH
        eta0 = (y - stcprm(12)) * stcprm(15) / REARTH
        xi = xi0 * stcprm(13) - eta0 * stcprm(14)
        eta = xi0 * stcprm(14) + eta0 * stcprm(13)
        vsq = stcprm(1) * (xi*xi + eta*eta) - 2. * eta
        if (stcprm(1) * vsq + 1. .le. 3.e-8) then
C/* Case of point at Projection Pole */
          map(1) = 0.
          map(2) = 0.
          map(3) = sign(1.,stcprm(1))
        else
          if ( stcprm(1) .ge. 1.) then
C/* Stereographic Case, std */
            fact = 2. / (2. + vsq)
            map(2) = xi * fact
            map(1) = (1. - eta) * fact
            map(3) = fact - 1.
          else if (stcprm(1) .le. -1.) then
C/* Stereographic Case, centered at antipodes */
            fact = -2. / (2. - vsq)
            map(2) = xi * fact
            map(1) = (- 1. - eta) * fact
            map(3) = fact + 1.
          else
C/* Mercator or Lambert Case */
            ymerc = -.5 * lnabva(stcprm(1), vsq)
            if ( ymerc .gt. 0.) then
              fact = exp(-ymerc)
              cosphi = fact / (1. + fact * fact)
              cosphi = cosphi + cosphi
              map(3) = 1. - fact * cosphi
            else
              fact = exp(ymerc)
              cosphi = fact / (1. + fact * fact)
              cosphi = cosphi + cosphi
              map(3) = fact * cosphi - 1.
            endif
C/*      map.v[2] = 1. - fact * cosphi; */
            if (abs(map(3)) .lt. 1.) then
              theta = atnabv(stcprm(1), xi, eta)
              map(1) = cosphi * cos(theta)
              map(2) = cosphi * sin(theta)
            else
              map(1) = 0.
              map(2) = 0.
            endif
          endif
        endif
      return
      end
      subroutine xy_xe(stcprm,  x,y, xi,eta)
      real stcprm(15)
      parameter (REARTH=6371.2)
        xtmp = x - stcprm(11)
        ytmp = y - stcprm(12)
        xi = stcprm(15) / REARTH *
     C    (stcprm(13) * xtmp - stcprm(14) * ytmp)
        eta = stcprm(15) / REARTH *
     C    (stcprm(13) * ytmp + stcprm(14) * xtmp)
        return
      end