      PROGRAM contour

      INTEGER PGBEG
!
      IF (PGBEG(0,'?',1,1) .NE. 1) STOP

      WRITE (*,'(A)') ' Routine PGCONT with PGCONL labels'
      CALL PGEX36

      CALL PGEND
      END program contour




SUBROUTINE PGEX36

      integer,parameter :: points1=251, points2=26, levels=20
      integer node, plot
real,parameter :: x1start=0.0, x2start=-13.5
real x1, x2, u1(points1,points2), u2(points1,points2), &
 u1min, u1max, u2min, u2max, x1end

INTEGER I,J
      REAL F(points1,points2),FMIN,FMAX,ALEV(1),TR(6)
      CHARACTER*32 LABEL
      DATA TR /x1start, 0.5, 0.0, x2start, 0.0, 1.0/

! use my data instead
      open(unit=101,file="data",form="formatted",status="old")
      do j=1,points2
      do i=1,points1
       read (101,*) node, x1, x2, u1(i,j), u2(i,j)
      end do
      end do
      close(101)

! convert to micron

u1 = u1*1.0e3
u2 = u2*1.0e3

! get the min/max

u1min = minval(u1)
u1max = maxval(u1)
u2min = minval(u2)
u2max = maxval(u2)

! do 2 plots
plt: do plot=1,2

 ! Clear the screen. Set up window and viewport.
 CALL PGPAGE
 CALL PGSVP( 0.1, 0.95, 0.1, 0.9 )
! CALL PGMTXT('t',1.0,0.0,0.0, 'u1, displacement along x1, \(0638)m')

if (plot .eq. 1) then
 x1end = 20.0
 CALL PGSWIN(1.0, x1end, x2start+1.0, -x2start-1.0)
 CALL PGBOX('bcnst',0.0,0,'bcnst',0.0,0)
 call pglab("x1,mm","x2,mm","u1, displacement along x1, \(0638)m")
 f = u1
 fmin = u1min
 fmax = u1max
end if

if (plot .eq. 2) then
 x1end = 20.0
 CALL PGSWIN(1.0, x1end, x2start+1.0, -x2start-1.0)
 CALL PGBOX('bcnst',0.0,0,'bcnst',0.0,0)
 call pglab("x1,mm","x2,mm","u2, displacement along x2, \(0638)m")
 f = u2
 fmin = u2min
 fmax = u2max
end if

!
! Draw the map.  PGCONT is called once for each contour, using
! different line attributes to distinguish contour levels.
!
      CALL PGBBUF
      DO 40 I=1,levels
          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/(levels-1)
          IF (MOD(I,5).EQ.0) THEN
              CALL PGSLW(3)
          ELSE
              CALL PGSLW(1)
          END IF
          IF ( alev(1) .lt. 0.0) THEN
              CALL PGSCI(2)
              CALL PGSLS(2)
          ELSE
              CALL PGSCI(3)
              CALL PGSLS(1)
          END IF
          CALL PGCONT(F,points1,points2,1,int(x1end/tr(2)),1,points2,ALEV,-1,TR)
   40 CONTINUE
      CALL PGSLW(1)
      CALL PGSLS(1)
      CALL PGEBUF
!
! Label the contours with PGCONL. Only even-numbered contours
! are labelled.
!
      CALL PGBBUF
      DO 50 I=1,levels
          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/(levels-1)
!          WRITE (LABEL,'(I2)') I
         WRITE (LABEL,'(F8.2)') ALEV
          IF ( alev(1) .lt. 0.0 ) THEN
              CALL PGSCI(2)
          ELSE
              CALL PGSCI(3)
          END IF
       CALL PGCONL(F,points1,points2,1,int(x1end/tr(2)),1,points2,ALEV,TR,LABEL,16,4)
 50   CONTINUE
      CALL PGSCI(1)
      CALL PGEBUF

end do plt

END subroutine pgex36
