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