SUBROUTINE VROUT (LUN,N,NT,LISTC,LPTR,LEND,XLON,XLAT,IER)
      INTEGER LUN, N, NT, LISTC(*), LPTR(*), LEND(N), IER
      REAL    XLON(NT), XLAT(NT)

!***********************************************************
!
!      Modified by Fred Fox, WERU 04/02/2010 based on VRPLOT
!                                              From STRIPACK
!                                            Robert J. Renka
!                                  Dept. of Computer Science
!                                       Univ. of North Texas
!                                           renka@cs.unt.edu
!                                                   07/16/98
!
!   This subroutine writes out a list of point pairs containing all
! the voronoi region boundary lines. Each pair of lines is preceded
! by the GRASS standard ASCII line pair specifier 'L 2', which can
! be easily removed by sed to create a gnuplot compatable input.
!
!   The parameters defining the Voronoi diagram may be com-
! puted by Subroutine CRLIST.
!
! On input:
!
!       LUN = Logical unit number in the range 0 to 99.
!             The unit should be opened with an appropriate
!             file name before the call to this routine.
!
!       N = Number of nodes (Voronoi centers) and Voronoi
!           regions.  N .GE. 3.
!
!       NT = Number of Voronoi region vertices (triangles,
!            including those in the extended triangulation
!            if the number of boundary nodes NB is nonzero):
!            NT = 2*N-4.
!
!       LISTC = Array of length 3*NT containing triangle
!               indexes (indexes to XLON and XLAT) stored
!               in 1-1 correspondence with LIST/LPTR entries
!               (or entries that would be stored in LIST for
!               the extended triangulation):  the index of
!               triangle (N1,N2,N3) is stored in LISTC(K),
!               LISTC(L), and LISTC(M), where LIST(K),
!               LIST(L), and LIST(M) are the indexes of N2
!               as a neighbor of N1, N3 as a neighbor of N2,
!               and N1 as a neighbor of N3.  The Voronoi
!               region associated with a node is defined by
!               the CCW-ordered sequence of circumcenters in
!               one-to-one correspondence with its adjacency
!               list (in the extended triangulation).
!
!       LPTR = Array of length 3*NT = 6*N-12 containing a
!              set of pointers (LISTC indexes) in one-to-one
!              correspondence with the elements of LISTC.
!              LISTC(LPTR(I)) indexes the triangle which
!              follows LISTC(I) in cyclical counterclockwise
!              order (the first neighbor follows the last
!              neighbor).
!
!       LEND = Array of length N containing a set of
!              pointers to triangle lists.  LP = LEND(K)
!              points to a triangle (indexed by LISTC(LP))
!              containing node K for K = 1 to N.
!
!       XLON,XLAT = Arrays of length NT containing the Longitude, 
!                 Latitude coordinates of the triangle
!                  circumcenters (Voronoi vertices).
!
! Input parameters are not altered by this routine.
!
! On output:
!
!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 if LUN, , N, or NT is outside
!                     its valid range.
!             IER = 3 if an error was encountered in writing
!                     to unit LUN.
!
! Modules required by VRPLOT:  None
!
! Intrinsic functions called by VRPLOT:  ABS, ATAN, COS,
!                                          NINT, REAL, SIN,
!                                          SQRT
!
!***********************************************************
!
      INTEGER KV1, KV2, LP, LPL, N0

! Local parameters:

! KV1,KV2 =   Endpoint indexes of a Voronoi edge
! LP =        LIST index (pointer)
! LPL =       Pointer to the last neighbor of N0
! N0 =        Index of a node

! Test for invalid parameters.

      IF (LUN .LT. 0  .OR.  LUN .GT. 99  .OR.                           &
     &    N .LT. 3  .OR.  NT .NE. 2*N-4)                                &
     &  GO TO 11

! Loop on nodes (Voronoi centers) N0.
!   LPL indexes the last neighbor of N0.

      DO 3 N0 = 1,N
        LPL = LEND(N0)

! Set KV2 to the first (and last) vertex index

        KV2 = LISTC(LPL)

! Loop on neighbors N1 of N0.  For each triangulation edge
!   N0-N1, KV1-KV2 is the corresponding Voronoi edge.

        LP = LPL
    1   LP = LPTR(LP)
          KV1 = KV2
          KV2 = LISTC(LP)

! Add edge KV1-KV2 to the path iff KV2 > KV1
!   (so that the edge is drawn only once).

          IF (KV2 .GT. KV1) THEN
              WRITE (LUN,1010,ERR=13) 'L 2'
              WRITE (LUN,1020,ERR=13) XLON(KV1), XLAT(KV1)
              WRITE (LUN,1020,ERR=13) XLON(KV2), XLAT(KV2)
          END IF

! Bottom of loops.

    2     IF (LP .NE. LPL) GO TO 1
    3   CONTINUE

!add format statements to eliminate space at beginning of line
 1010 format(a3)
 1020 format(f20.15,1x,f20.15)

! No error encountered.

      IER = 0
      RETURN

! Invalid input parameter LUN, , N, or NT.

   11 IER = 1
      write(*,*) 'Bad file handle or bad number'
      RETURN

! Error writing to unit LUN.

   13 IER = 3
      RETURN
      END