c
c     Copyright (c) 1999,2000, Compaq Computer Corporation
c
c     See the file "license.terms" for information on usage and redistribution
c     of this file, and for a DISCLAIMER OF ALL WARRANTIES.
c
c     plot multiple 1d cut lines
c
      subroutine compare_plot(lsave,lplotfirst,plot_mode)
      implicit none
c
c     Programmer: John Faricelli
c     Date: 10-Oct-85
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'compare.cmn'
      include 'units.cmn'
      include 'data.cmn'
      include 'defaults.cmn'
c      include 'annotations.cmn'
      include 'io.cmn'
      include 'user.cmn'
      include 'title.cmn'
      include 'callable.cmn'
c
      logical lsave, lplotfirst
      integer plot_mode
c
c     locals
c
      logical do_header, different,
     1        fake_ln03, data_selected, first_time, saw_log, saw_linear,
     1        x_char_data(2)
      character*40 first_quantity
      character*1  answer, first_axis, first_plane
      integer i, j, slot, ilen, ulen, string_len
      integer v_exponent(2), v_minor_tic_freq(2)
      real real_vec(10)
      real v_tic_major(2), delta
      real v_min(2), v_max(2)
c
CC BEGIN MENU
      real hmin_full(2), hmax_full(2), h_min(2), h_max(2)
      real prev_hmin_plot(2), prev_hmax_plot(2),
     1     prev_h_tic_major(2)
      integer prev_h_minor_tic_freq(2), prev_h_exponent(2)
c
      real vmin_full(2), vmax_full(2), vmin_plot(2), vmax_plot(2)
      real prev_vmin_plot(2), prev_vmax_plot(2),
     2     prev_v_tic_major(2)
      integer item_number, int_vec(2),
     1     prev_v_minor_tic_freq(2), iuser_hscale(2)
      integer iuser_vscale(2), prev_iuser_vscale(2),
     1        prev_v_exponent(2)
      integer ipage, itag, ierror, ixaxis, iyaxis
      character*80 string
      character*27 xaxis_name(2), yaxis_name(2)
      data xaxis_name /'BOTTOM (X) AXIS PARAMETERS:',
     1                 'TOP    (X) AXIS PARAMETERS:'/
      data yaxis_name /'LEFT   (Y) AXIS PARAMETERS:',
     1                 'RIGHT  (Y) AXIS PARAMETERS:'/
      character*31 axis_prompt(2)
      data axis_prompt /'21) Switch to top/right axes  ',
     1                  '21) Switch to left/bottom axes'/
      character*6 str_hlog, str_vlog
      character*3 str_habs, str_vabs, str_labels, str_aspect
      real h, h1(2), h2(2), h3(2), h4(2), x1(2), x2(2), y1(2), y2(2),
     1     lastx, lasty, delx, dely, temp_min, temp_max
      logical eof, do_vlog(2), do_vabs(2), do_hlog(2), do_habs(2),
     1        do_vuscale(2), do_huscale(2), prev_vabs(2), prev_habs(2),
     2        abort, is_graphics_terminal, lopen, delayed_message,
     3        do_zoom, do_unzoom, prev_v_saved(2), prev_h_saved(2),
     4        lannotate, do_probe, first_probe, lfirst, lreset,
     5        prev_h_used, prev_v_used, update_h_axis, update_v_axis,
     6        do_unzoom2, do_zoom2, do_draw_to_scale
CC END MENU
c
      real xndc0, yndc0, xndc1, yndc1, ndc_aspect_ratio, aspect_ratio
c
      real rdum, rdum2, rdum3
      integer idum, idum2, idum3
      character*1 cdum
c
      lopen = .FALSE.
      lannotate = .FALSE.
      do_header = .TRUE.
      do_draw_to_scale = .FALSE.
      delayed_message = .FALSE.
      ulen = string_len(pm_length_unit)
      lfirst = lplotfirst
c
c     Check for data selected
      data_selected = .FALSE.
      do i=1,n_cuts
         slot = slot_index(i)
         if( multi_1d(slot)%selected .eq. 1 ) data_selected = .TRUE.
      end do
      if( .NOT. data_selected ) then
         if( lbatch ) then
            call plt$error_message('compare_plot',
     1         'No data is selected')
         else
            call beep_terminal
            call get_char(answer,
     1      'Error: No data is selected! Hit return to continue...',0)
         endif
         return
      endif
c
c     Determine which, if any, datasets have X character data
      x_char_data(1) = .FALSE.
      x_char_data(2) = .FALSE.
      do i=1,n_cuts
         slot = slot_index(i)
         if( multi_1d(slot)%selected .eq. 1 ) then
            if( multi_1d(slot)%x_char_data ) then
               x_char_data(multi_1d(slot)%xaxis_tag) = .TRUE.
            endif
         endif
      end do
c
c     Initialize axis "page" number
c     ipage=1 lists left, lower axes
c     ipage=2 lists right, upper axes
      ipage = 1
c
c     Horizontal scale factor not implemented yet (if ever)
      do i=1,2
         do_huscale(i)   = .FALSE.
         iuser_hscale(i) = 0
      end do
c
c     Initialize data which gets saved over subroutine calls
c
      if( comp_lfirst ) then
         comp_lfirst = .FALSE.
         call pm_init_comp_sticky
c
         num_titles = 1
         title_dat(1) = ' '
c
         do i=1,2
            comp_h_used(i) = .FALSE.
            comp_v_used(i) = .FALSE.
            do j=1,n_cuts
               slot = slot_index(j)
               if( multi_1d(slot)%selected .eq. 1 .and.
     1             multi_1d(slot)%xaxis_tag .eq. i ) then
                   comp_h_used(i) = .TRUE.
               endif
               if( multi_1d(slot)%selected .eq. 1 .and.
     1             multi_1d(slot)%yaxis_tag .eq. i ) then
                  comp_v_used(i) = .TRUE.
               endif
            end do
         end do
c
c        Find min/max of coordinates, discard x bound information
         do i=1,2
            do_hlog(i)  = .FALSE.
            do_habs(i)  = .FALSE.
            call recompute_minmax_coords(comp_h_min(i), comp_h_max(i),
     1           rdum, rdum2, comp_h_tic_major(i),
     2           comp_h_minor_tic_freq(i),
     3           iuser_hscale(i), comp_h_exponent(i),
     4           do_habs(i), do_hlog(i), do_huscale(i), i,
     5           work_vec(1,1))
            if( comp_h_exponent(i) .ne. 0 ) then
               rdum = 10.0**comp_h_exponent(i)
               comp_h_min(i) = comp_h_min(i)/rdum
               comp_h_max(i) = comp_h_max(i)/rdum
            endif
            if( plot_mode .eq. 1 ) then
c              Check for funny limits. Avoid min = max.
               temp_min = comp_h_min(i)
               temp_max = comp_h_max(i)
               call pm_avoid_funny_limits(temp_min, temp_max,
     1                 comp_h_min(i), comp_h_max(i))
               endif
               if( x_char_data(i) ) then
                  comp_h_tic_major(i) = 1.0
                  comp_h_minor_tic_freq(i) = 1
               endif
            if( plot_mode .eq. 2 ) then  ! Barchart plot
               comp_h_min(i) = comp_h_min(i) - 1.0
               comp_h_max(i) = comp_h_max(i) + 1.0
               comp_h_tic_major(i) = 1.0
               comp_h_minor_tic_freq(i) = 1
            endif
         end do
c
c        Try to determine default horizontal axis labels,
c        depending on quantities involved.
c        Horizontal axis label
         different = .FALSE.
         first_time = .TRUE.
         do i=1,n_cuts
            slot = slot_index(i)
c           Only look at data which has been selected
            if( multi_1d(slot)%selected .eq. 1 ) then
               if( first_time ) then
                  first_axis = multi_1d(slot)%cut_line_axis
                  first_plane = multi_1d(slot)%cut_plane_axis
                  first_time = .FALSE.
              endif
               if( multi_1d(slot)%cut_line_axis .ne. first_axis ) 
     1            different = .TRUE.
               if( multi_1d(slot)%cut_plane_axis .ne. first_plane ) 
     1            different = .TRUE.
**               if( multi_1d(slot)%type .eq. COMP_USER_TYPE )
**     1            different = .TRUE.
           endif
         end do
         if( different ) then
            comp_hunits(1) = ' '
            comp_horiz_label(1) = ' '
         else
            comp_hunits(1) = pm_length_unit(1:ulen)
            if( binfiletype .eq. DT_TIF .and. tif_1d_data ) then
               comp_horiz_label(1) = tif_names(tif_ixcol)
               comp_hunits(1) = tif_units(tif_ixcol)
            else if( binfiletype .eq. DT_GRAPES) then
               comp_horiz_label(1) = grapes_x_quantity
               comp_hunits(1) = grapes_x_units
            else if( binfiletype .eq. DT_ASCII ) then
               comp_horiz_label(1) = user_1d_x_name
               comp_hunits(1) = user_1d_x_units
            else if( first_axis .eq. 'X' ) then
               if( first_plane .eq. 'Z' ) then
                  comp_horiz_label(1) = 'Depth'
               else if( first_plane .eq. 'Y' ) then
                  comp_horiz_label(1) = 'Width'
               else
                  comp_horiz_label(1) = 'Time'
                  comp_hunits(1) = 's'
               endif
            else if( first_axis .eq. 'Y' ) then
               if( first_plane .eq. 'Z' ) then
                  comp_horiz_label(1) = 'Length'
               else if( first_plane .eq. 'X') then
                  comp_horiz_label(1) = 'Width'
               else
                  comp_horiz_label(1) = 'Time'
                  comp_hunits(1) = 's'
               endif
            else if( first_axis .eq. 'Z' ) then
               if( first_plane .eq. 'Y' ) then
                  comp_horiz_label(1) = 'Length'
               else
                  comp_horiz_label(1) = 'Depth'
               endif
            else if( first_axis .eq. 'T' ) then
               if( first_plane .eq. 'Y' ) then
                  comp_horiz_label(1) = 'Length'
               else
                  comp_horiz_label(1) = 'Depth'
               endif
            endif
         endif
c
c        Vertical axis label
         different = .FALSE.
         first_time = .TRUE.
         do i=1,n_cuts
            slot = slot_index(i)
c           Only look at data which has been selected
            if( multi_1d(slot)%selected .eq. 1 ) then
               if( first_time ) then
                  first_quantity = multi_1d(slot)%quantity
                  first_time = .FALSE.
               endif
               if( multi_1d(slot)%quantity .ne. first_quantity ) 
     1            different = .TRUE.
**               if( multi_1d(slot)%type .eq. COMP_USER_TYPE )
**     1            different = .TRUE.
            endif
         end do
         if( different ) then
            comp_vunits(1) = ' '
            comp_vert_label(1) = ' '
         else
            comp_vert_label(1) = first_quantity
            slot = slot_index(1)
            if( multi_1d(slot)%type .eq. COMP_CUT_TYPE .or.
     1          multi_1d(slot)%type .eq. COMP_INTERFACE_TYPE .or.
     2          multi_1d(slot)%type .eq. COMP_USER_TYPE ) then
               comp_vunits(1) = multi_1d(slot)%units
            else
               comp_vunits(1) = multi_1d(slot)%int_units
            endif
         endif
c
c        Compute default vertical axis parameters
c
         do j=1,2
c           Vertical axis logarithmic?
            saw_log    = .FALSE.
            saw_linear = .FALSE.
            do i=1,n_cuts
               slot = slot_index(i)
c              Only look at data which has been selected
               if( multi_1d(slot)%selected .eq. 1 .and.
     1             multi_1d(slot)%yaxis_tag .eq. j ) then
                  if( multi_1d(slot)%log .eq. .TRUE. ) then
                     saw_log = .TRUE.
                  else
                     saw_linear = .TRUE.
                  endif
               endif
            end do
            if( saw_log .and. .NOT. saw_linear ) then
               do_vlog(j) = .TRUE.
               do_vabs(j) = .TRUE.
            else
               do_vlog(j) = .FALSE.
               do_vabs(j) = .FALSE.
            endif
c
            do_vuscale(j) = .FALSE.
            v_exponent(j) = 0
            iuser_vscale(j) = 0
            call recompute_minmax_comp(
     1                         vmin_full(j), vmax_full(j),
     2                         vmin_plot(j), vmax_plot(j),
     3                         v_tic_major(j), v_minor_tic_freq(j),
     4                         iuser_vscale(j), v_exponent(j),
     5                         do_vabs(j), do_vlog(j), do_vuscale(j),
     6                         j, work_vec(1,1), plot_mode)
         end do
      endif
c
      if( lsave ) then
c        Find min/max of coordinates
         do_hlog(1) = comp_h_log(1)
         do_habs(1) = comp_h_abs(1)
         do_hlog(2) = comp_h_log(2)
         do_habs(2) = comp_h_abs(2)
         do i=1,2
            if( comp_h_used(i) ) then
               idum2 = 0
               call recompute_minmax_coords(hmin_full(i), hmax_full(i),
     1              rdum, rdum2, rdum3, idum, idum2, idum3,
     2              do_habs(i), do_hlog(i), do_huscale(i), i,
     3              work_vec(1,1))
c
c              Set axis parameters if autoscale flag was set
               if( comp_h_autoscale_min(i) .or.
     1             comp_h_autoscale_max(i) ) then
                  if(comp_h_autoscale_min(i)) comp_h_min(i)=hmin_full(i)
                  if(comp_h_autoscale_max(i)) comp_h_max(i)=hmax_full(i)
c                 Sanity check the min/max
                  lreset = .FALSE.
                  if( comp_h_min(i) .ge. comp_h_max(i) ) then
                     comp_h_min(i) = hmin_full(i)
                     comp_h_max(i) = hmax_full(i)
                     lreset = .TRUE.
                     string = 'Invalid user specified min/max'
                     delayed_message = .TRUE.
                  endif
                  if( .NOT. do_hlog(i) ) then
                     comp_h_tic_major(i) = rdum3
                     comp_h_minor_tic_freq(i) = idum
                     comp_h_exponent(i) = idum3
                     if( comp_h_exponent(i) .ne. 0 ) then
                        rdum = 10.0**comp_h_exponent(i)
                        if(comp_h_autoscale_min(i) .or. lreset)
     1                     comp_h_min(i) = comp_h_min(i)/rdum
                        if(comp_h_autoscale_max(i) .or. lreset)
     1                     comp_h_max(i) = comp_h_max(i)/rdum
                     endif
                     if( plot_mode .eq. 2 .or. x_char_data(i) ) then
                        comp_h_tic_major(i) = 1.0
                        comp_h_minor_tic_freq(i) = 1
                     endif
                  endif
               endif
c              Compute defaults for "prev" quantities
               if( do_hlog(i) ) then
                  prev_h_saved(i) = .FALSE.
               else
                  prev_h_saved(i) = .TRUE.
                  prev_hmin_plot(i) = comp_h_min(i)
                  prev_hmax_plot(i) = comp_h_max(i)
                  prev_h_tic_major(i) = comp_h_tic_major(i)
                  prev_h_minor_tic_freq(i) = comp_h_minor_tic_freq(i)
                  prev_h_exponent(i) = comp_h_exponent(i)
               endif
            endif
         end do
c
c        Find min/max of data
         do_vlog(1) = comp_v_log(1)
         do_vabs(1) = comp_v_abs(1)
         do_vlog(2) = comp_v_log(2)
         do_vabs(2) = comp_v_abs(2)
         do_vuscale(1) = .FALSE.
         do_vuscale(2) = .FALSE.
         vmin_plot(1) = comp_v_min(1)
         vmax_plot(1) = comp_v_max(1)
         v_exponent(1) = comp_v_exponent(1)
         v_tic_major(1) = comp_v_tic_major(1)
         v_minor_tic_freq(1) = comp_v_minor_tic_freq(1)
         vmin_plot(2) = comp_v_min(2)
         vmax_plot(2) = comp_v_max(2)
         v_exponent(2) = comp_v_exponent(2)
         v_tic_major(2) = comp_v_tic_major(2)
         v_minor_tic_freq(2) = comp_v_minor_tic_freq(2)
         iuser_vscale(1) = 0
         iuser_vscale(2) = 0
         do i=1,2
            if( comp_v_used(i) ) then
               call recompute_minmax_comp(
     1                         vmin_full(i), vmax_full(i),
     2                         rdum, rdum2,
     3                         rdum3, idum,
     4                         iuser_vscale(i), idum2,
     5                         do_vabs(i), do_vlog(i), do_vuscale(i),
     6                         i, work_vec(1,1), plot_mode)
               if( comp_v_autoscale_min(i) .or.
     1             comp_v_autoscale_max(i)) then
c                 Set defaults min/max for vertical axis
                  if(comp_v_autoscale_min(i)) comp_v_min(i) = rdum
                  if(comp_v_autoscale_max(i)) comp_v_max(i) = rdum2
c                 Sanity check the min/max
                  if( comp_v_min(i) .ge. comp_v_max(i) ) then
                     comp_v_min(i) = rdum
                     comp_v_max(i) = rdum2
                     string = 'Invalid user specified min/max'
                     delayed_message = .TRUE.
                  endif
                  vmin_plot(i) = comp_v_min(i)
                  vmax_plot(i) = comp_v_max(i)
                  if( .NOT. do_vlog(i) ) then
                     comp_v_tic_major(i) = rdum3
                     comp_v_minor_tic_freq(i) = idum
                     comp_v_exponent(i) = idum2
                     v_tic_major(i) = comp_v_tic_major(i)
                     v_minor_tic_freq(i) = comp_v_minor_tic_freq(i)
                     v_exponent(i) = comp_v_exponent(i)
                  endif
               else
                  if( idum2 .ne. comp_v_exponent(i) ) then
                     do_vuscale(i) = .TRUE.
                     iuser_vscale(i) = comp_v_exponent(i)
                  endif
                  v_exponent(i) = comp_v_exponent(i)
               endif
c              Compute defaults for "prev" quantities
               if( do_vlog(i) ) then
                  prev_v_saved(i) = .FALSE.
               else
                  prev_v_saved(i) = .TRUE.
                  prev_vmin_plot(i) = comp_v_min(i)
                  prev_vmax_plot(i) = comp_v_max(i)
                  prev_v_tic_major(i) = comp_v_tic_major(i)
                  prev_v_minor_tic_freq(i) = comp_v_minor_tic_freq(i)
                  prev_iuser_vscale(i) = iuser_vscale(i)
                  prev_v_exponent(i) = comp_v_exponent(i)
               endif
            endif
         end do
      else
c        Update hmin_full, hmax_full, vmin_full, vmax_full,
c        in case new data has been added.
c        Also, check to see if we've added a new axis...
         do i=1,2
            prev_h_used = comp_h_used(i)
            prev_v_used = comp_h_used(i)
            update_h_axis = .FALSE.
            update_v_axis = .FALSE.
            comp_h_used(i) = .FALSE.
            comp_v_used(i) = .FALSE.
            do j=1,n_cuts
               slot = slot_index(j)
               if( multi_1d(slot)%selected .eq. 1 .and.
     1             multi_1d(slot)%xaxis_tag .eq. i ) then
                   comp_h_used(i) = .TRUE.
               endif
               if( multi_1d(slot)%selected .eq. 1 .and.
     1             multi_1d(slot)%yaxis_tag .eq. i ) then
                   comp_v_used(i) = .TRUE.
                endif
            end do
            if( comp_h_used(i) .and. .not. prev_h_used ) then
               update_h_axis = .TRUE.
            endif
            if( comp_v_used(i) .and. .not. prev_v_used ) then
               update_v_axis = .TRUE.
            endif
c
            call recompute_minmax_coords(hmin_full(i), hmax_full(i),
     1           rdum, rdum2, rdum3, idum, idum2, idum3,
     2           do_habs(i), do_hlog(i), do_huscale(i), i,
     3           work_vec(1,1))
            if( update_h_axis ) then
               comp_h_min(i) = hmin_full(i)
               comp_h_max(i) = hmax_full(i)
               comp_h_tic_major(i) = rdum3
               comp_h_minor_tic_freq(i) = idum
               comp_h_exponent(i)  = idum3
               if( comp_h_exponent(i) .ne. 0 ) then
                  rdum = 10.0**comp_h_exponent(i)
                  comp_h_min(i) = comp_h_min(i)/rdum
                  comp_h_max(i) = comp_h_max(i)/rdum
               endif
               if( x_char_data(i) .or. plot_mode .eq. 2 ) then
                  comp_h_tic_major(i) = 1.0
                  comp_h_minor_tic_freq(i) = 1
               endif
            endif
c
            call recompute_minmax_comp(
     1                         vmin_full(i), vmax_full(i),
     2                         rdum, rdum2,
     3                         rdum3, idum,
     4                         0, idum2,
     5                         do_vabs(i), do_vlog(i), .FALSE.,
     6                         i, work_vec(1,1), plot_mode)
            if( update_v_axis ) then
               vmin_plot(i) = rdum
               vmax_plot(i) = rdum2
               v_tic_major(i) = rdum3
               v_minor_tic_freq(i) = idum
               v_exponent(i) = idum2
            endif
c
            if( do_hlog(i) ) then
               prev_h_saved(i) = .FALSE.
            else
               prev_h_saved(i) = .TRUE.
               prev_hmin_plot(i) = comp_h_min(i)
               prev_hmax_plot(i) = comp_h_max(i)
               prev_h_tic_major(i) = comp_h_tic_major(i)
               prev_h_minor_tic_freq(i) = comp_h_minor_tic_freq(i)
               prev_h_exponent(i) = comp_h_exponent(i)
            endif
c
            if( do_vlog(i) ) then
               prev_v_saved(i) = .FALSE.
            else
               prev_v_saved(i) = .TRUE.
               prev_iuser_vscale(i) = 0
               prev_vmin_plot(i) = vmin_plot(i)
               prev_vmax_plot(i) = vmax_plot(i)
               prev_v_tic_major(i) = v_tic_major(i)
               prev_v_minor_tic_freq(i)  = v_minor_tic_freq(i)
               prev_v_exponent(i) = v_exponent(i)
            endif
         end do
      endif
c
      do i=1,2
         prev_habs(i)  = .FALSE.
         prev_vabs(i)  = .FALSE.
      end do
c
c     Take care of parameters set via callable interface
c     This routine may reset the X/Y axis limits
CC !!!NOTE WELL!!! This should work, but I should do something better
CC    about callable_title info
      if( callable_num_titles .eq. 1 ) then
         num_titles = 2
         title_dat(2) = callable_titles(1)
      endif
      call get_axis_callable(
     1      comp_h_min(1), comp_h_max(1),
     2      comp_h_exponent(1), comp_h_tic_major(1),
     3      comp_h_minor_tic_freq(1),
     4      comp_hunits(1),
     5      comp_horiz_label(1), do_hlog(1), do_habs(1),
     6      vmin_plot(1), vmax_plot(1),
     7      v_exponent(1), v_tic_major(1), v_minor_tic_freq(1),
     8      comp_vunits(1),
     9      comp_vert_label(1), do_vlog(1), do_vabs(1))
c
c     If batch processing, perform hardcopy and return
      if( lbatch .or. lrestore ) then
         previous_plot_device = plot_device
         if( lbatch ) then
            plot_device = hardcopy_device
            if( lplotfile ) then
               plot_file_name = given_plotfilename
            endif
         endif
c
c        Stuff axis information into save structures
         call copy_axis_to_save(PLOT_COMPARE_COMMAND,
     1      comp_h_min(1), comp_h_max(1),
     2      comp_h_exponent(1), comp_h_tic_major(1),
     3      comp_h_minor_tic_freq(1),
     4      comp_hunits(1), comp_horiz_label(1),
     5      do_hlog(1), do_habs(1), comp_h_used(1),
     6      vmin_plot(1), vmax_plot(1),
     7      v_exponent(1), v_tic_major(1), v_minor_tic_freq(1),
     8      comp_vunits(1), comp_vert_label(1),
     9      do_vlog(1), do_vabs(1), comp_v_used(1),
     A      rdum, rdum,
     B      idum, rdum, idum, cdum,
     C      cdum, .FALSE., .FALSE., .FALSE.)
         call copy_axis2_to_save(comp_h_min(2), comp_h_max(2),
     1      comp_h_exponent(2), comp_h_tic_major(2),
     2      comp_h_minor_tic_freq(2),
     3      comp_hunits(2), comp_horiz_label(2),
     4      do_hlog(2), do_habs(2), comp_h_used(2),
     5      vmin_plot(2), vmax_plot(2),
     6      v_exponent(2), v_tic_major(2),
     7      v_minor_tic_freq(2),
     8      comp_vunits(2), comp_vert_label(2),
     9      do_vlog(2), do_vabs(2), comp_v_used(2))
c
c        Stuff data information into save structures
         call copy_compare_to_save(plot_mode)
c
c        Copy data into temp store, take absolute value, log if necessary...
         do 115 i=1,n_cuts
            slot = slot_index(i)
            itag = multi_1d(slot)%yaxis_tag
            do 116 j=1,multi_1d(slot)%np
               if( do_vabs(itag) .or. do_vlog(itag) ) then
                  work_vec(j,slot) = ABS(multi_1d(slot)%data(j))
                  if( do_vlog(itag) ) then
c                    check for log(0.0), fudge
                     if( work_vec(j,slot) .eq. 0.0 )
     1                  work_vec(j,slot) = PM_MIN_FLOAT
                     work_vec(j,slot) = ALOG10(work_vec(j,slot))
                  endif
               else
                  work_vec(j,slot) = multi_1d(slot)%data(j)
               endif
116         continue
115      continue
c
         do i=1,2
            if( comp_v_used(i) ) then
               if( do_vlog(i) ) then
                  v_min(i) = ALOG10(vmin_plot(i))
                  v_max(i) = ALOG10(vmax_plot(i))
               else
                  v_min(i) = vmin_plot(i)
                  v_max(i) = vmax_plot(i)
               endif
            endif
         end do
c
         do i=1,2
            if( comp_h_used(i) ) then
               if( do_hlog(i) ) then
                  h_min(i) = ALOG10(comp_h_min(i))
                  h_max(i) = ALOG10(comp_h_max(i))
               else
                  h_min(i) = comp_h_min(i)
                  h_max(i) = comp_h_max(i)
               endif
            endif
         end do
c
         call do_multi_plot(work_vec, h_min, h_max,
     1      comp_h_exponent,
     2      comp_h_tic_major, comp_h_minor_tic_freq,
     3      v_min, v_max,
     4      v_exponent, v_tic_major, v_minor_tic_freq,
     5      comp_hunits, comp_vunits,
     6      comp_horiz_label, comp_vert_label, do_hlog, do_vlog,
     7      comp_h_used, comp_v_used, do_header, plot_mode,
     8      lsave, lopen, lannotate, ierror)
c
         plot_device = previous_plot_device
         return
      endif
c
c     Main return point for menu
c
4000  continue
c
      if( do_hlog(ipage) ) then
         str_hlog = 'log'
      else
         str_hlog = 'linear'
      endif
      if( do_habs(ipage) ) then
         str_habs = 'yes'
      else
         str_habs = 'no'
      endif
c
      if( do_vlog(ipage) ) then
         str_vlog = 'log'
      else
         str_vlog = 'linear'
      endif
      if( do_vabs(ipage) ) then
         str_vabs = 'yes'
      else
         str_vabs = 'no'
      endif
      if( comp_nolabels ) then
         str_labels = 'yes'
      else
         str_labels = 'no'
      endif
      if( do_draw_to_scale ) then
         str_aspect = 'yes'
      else
         str_aspect = ' no'
      endif
c
c     Write out menu form:
      call clear_screen
      call display_string(1,1 ,'1)  Title:',10, 0)
c
      call display_string(3 ,1 ,xaxis_name(ipage),27,0)
      if( .NOT. comp_h_used(ipage) )
     1   call display_string(3 ,30,
     2   '(inactive: no data associated with this axis)',45,1)
      call display_string(4 ,1 ,'Data range:',11,0)
      call display_string(4 ,24,' to ',4,0)
      call display_string(6 ,1 ,'2)  Plot range:',15,0)
      call display_string(6 ,24,' to ',4,0)
      call display_string(7 ,1 ,'3)  Major tics at:',18,0)
      call display_string(8 ,1 ,'4)  Minor tic freq:',19,0)
      call display_string(9 ,1 ,'5)  Scale factor:       10**',28,0)
      call display_string(4 ,43,'6)  Linear/log: ',16, 0)
      call display_string(5 ,43,'7)  Absolute value: ',20, 0)
      call display_string(6 ,43,'8)  Recompute default axis parameters',
     1                    37,0)
      call display_string(7 ,43,'9)  Label: ',11, 0)
      call display_string(8 ,43,'10) Units: ',11, 0)
c
      call display_string(11,1 ,yaxis_name(ipage),27,0)
      if( .NOT. comp_v_used(ipage) )
     1   call display_string(11,30,
     2   '(inactive: no data associated with this axis)',45,1)
      call display_string(12,1 ,'Data range:',11,0)
      call display_string(12,24,' to',3,0)
      call display_string(14,1 ,'11) Plot range:',15,0)
      call display_string(14,24,' to',3,0)
      call display_string(15,1 ,'12) Major tics at:',18,0)
      call display_string(16,1 ,'13) Minor tic freq:',19,0)
      call display_string(17,1 ,'14) Scale factor:       10**',28,0)
      call display_string(12,43,'15) Linear/log: ',16, 0)
      call display_string(13,43,'16) Absolute value: ',20, 0)
      call display_string(14,43,'17) Recompute default axis parameters',
     1                    37,0)
      call display_string(15,43,'18) Label: ',11, 0)
      call display_string(16,43,'19) Units: ',11, 0)
      call display_string(18,43,'20) Suppress labels: ',21, 0)
      call display_string(19,43,axis_prompt(ipage),31, 0)
      call display_string(20,43,'22) Draw to scale: ',19, 0)
c
      call display_string(23,1 ,'Messages:',9,0)
c
c     Title
      call display_string(1 , 12, title_dat(1), 69, 1)
c
c     X Range
      call display_real  (4 , 15, hmin_full(ipage), 1)
      call display_real  (4 , 28, hmax_full(ipage), 1)
      call display_real  (6 , 15, comp_h_min(ipage), 1)
      call display_real  (6 , 28, comp_h_max(ipage), 1)
      if( do_hlog(ipage) ) then
         call display_string(7 , 23, '       NA', 9, 1)
         call display_string(8 , 27, '   NA', 5, 1)
      else    
c        X major tics
         call display_real  (7 , 23, comp_h_tic_major(ipage), 1)
c        X minor tic freq
         call display_int   (8 , 27, comp_h_minor_tic_freq(ipage), 5)
      endif
c     X scale factor
      call display_int   (9 , 29, comp_h_exponent(ipage),3)
c     Linear/log
      call display_string(4 , 59, str_hlog, 6, 1)
c     Absolute value
      call display_string(5 , 63, str_habs, 3, 1)
c     Horizontal label
      call display_string(7 , 54, comp_horiz_label(ipage), 27, 1)
c     Horizontal Units
      call display_string(8 , 54, comp_hunits(ipage), 20, 1)
c
c     Y Range
      call display_real  (12, 15, vmin_full(ipage), 1)
      call display_real  (12, 28, vmax_full(ipage), 1)
      call display_real  (14, 15, vmin_plot(ipage), 1)
      call display_real  (14, 28, vmax_plot(ipage), 1)
      if( do_vlog(ipage) ) then
         call display_string(15, 23, '       NA', 9, 1)
         call display_string(16, 27, '   NA', 5, 1)
      else
c        Y major tics
         call display_real  (15, 23, v_tic_major(ipage), 1)
c        Y minor tic freq
         call display_int   (16, 27, v_minor_tic_freq(ipage), 5)
      endif
c     Y scale factor
      call display_int  (17, 29, v_exponent(ipage), 3)
c     Linear/log
      call display_string(12, 59, str_vlog, 6, 1)
c     Absolute value
      call display_string(13, 63, str_vabs, 3, 1)
c     Vertical label
      call display_string(15, 54, comp_vert_label(ipage), 27, 1)
c     Vertical Units
      call display_string(16, 54, comp_vunits(ipage), 20, 1)
c     Label switch
      call display_string(18, 64, str_labels, 3, 1)
c     Aspect switch
      call display_string(20, 64, str_aspect, 3, 1)
c
c     Display any messaged from before window refresh
      if( delayed_message ) then
         call display_message(string)
         delayed_message = .FALSE.
      endif
c
2000  call display_prompt(
     1 'Enter item num or '//
     2 'Annotate Plot Hardcopy Keepwin Repaint Zoom Unzoom Sample Exit')
      if( lfirst ) then
         lfirst = .FALSE.
         eof = .FALSE.
         string = 'P'
      else
         call read_from_input_area(string, eof)
      endif
      if( eof ) goto 2000
      call upcase(string)
      if( string(1:1) .eq. 'P' ) then
         call display_prompt('Making plot')
c
c        Stuff axis information into save structures
         call copy_axis_to_save(PLOT_COMPARE_COMMAND,
     1      comp_h_min(1), comp_h_max(1),
     2      comp_h_exponent(1), comp_h_tic_major(1),
     3      comp_h_minor_tic_freq(1),
     4      comp_hunits(1), comp_horiz_label(1),
     5      do_hlog(1), do_habs(1), comp_h_used(1),
     6      vmin_plot, vmax_plot,
     7      v_exponent(1), v_tic_major(1), v_minor_tic_freq(1),
     8      comp_vunits(1), comp_vert_label(1),
     9      do_vlog(1), do_vabs(1), comp_v_used(1),
     A      rdum, rdum,
     B      idum, rdum, idum, cdum,
     C      cdum, .FALSE., .FALSE., .FALSE.)
         call copy_axis2_to_save(comp_h_min(2), comp_h_max(2),
     1      comp_h_exponent(2), comp_h_tic_major(2),
     2      comp_h_minor_tic_freq(2),
     3      comp_hunits(2), comp_horiz_label(2),
     4      do_hlog(2), do_habs(2), comp_h_used(2),
     5      vmin_plot(2), vmax_plot(2),
     6      v_exponent(2), v_tic_major(2),
     7      v_minor_tic_freq(2),
     8      comp_vunits(2), comp_vert_label(2),
     9      do_vlog(2), do_vabs(2), comp_v_used(2))
c
c        Stuff data information into save structures
         call copy_compare_to_save(plot_mode)
c
c        Copy data into temp store, take absolute value, log if necessary...
         do 100 i=1,n_cuts 
            slot = slot_index(i)
            itag = multi_1d(slot)%yaxis_tag
            do 110 j=1,multi_1d(slot)%np
               if( do_vabs(itag) .or. do_vlog(itag) ) then
                  work_vec(j,slot) = ABS(multi_1d(slot)%data(j))
                  if( do_vlog(itag) ) then
c                    check for log(0.0), fudge
                     if( work_vec(j,slot) .eq. 0.0 )
     1                  work_vec(j,slot) = PM_MIN_FLOAT
                     work_vec(j,slot) = ALOG10(work_vec(j,slot))
                  endif
               else
                  work_vec(j,slot) = multi_1d(slot)%data(j)
               endif
110         continue
100      continue
c
         do i=1,2
            if( comp_v_used(i) ) then
               if( do_vlog(i) ) then
                  v_min(i) = ALOG10(vmin_plot(i))
                  v_max(i) = ALOG10(vmax_plot(i))
               else
                  v_min(i) = vmin_plot(i)
                  v_max(i) = vmax_plot(i)
               endif
            endif
         end do
c
         do i=1,2
            if( comp_h_used(i) ) then
               if( do_hlog(i) ) then
                  h_min(i) = ALOG10(comp_h_min(i))
                  h_max(i) = ALOG10(comp_h_max(i))
               else
                  h_min(i) = comp_h_min(i)
                  h_max(i) = comp_h_max(i)
               endif
            endif
         end do
c
         call do_multi_plot(work_vec, h_min, h_max,
     1      comp_h_exponent,
     2      comp_h_tic_major, comp_h_minor_tic_freq,
     3      v_min, v_max,
     4      v_exponent, v_tic_major, v_minor_tic_freq,
     5      comp_hunits, comp_vunits,
     6      comp_horiz_label, comp_vert_label, do_hlog, do_vlog,
     7      comp_h_used, comp_v_used, do_header, plot_mode,
     8      lsave, lopen, lannotate, ierror)
         if( ierror .eq. PLT$NORMAL ) then
            lopen = .TRUE.
         else
            call plt$get_error_message(ierror,string)
            delayed_message = .TRUE.
         endif
c
         goto 4000
c
      else if( string(1:1) .eq. 'Z' .or. string(1:1) .eq. 'U') then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: cannot zoom on this graphics device')
               goto 2000
            else
               do_zoom    = .FALSE.
               do_zoom2   = .FALSE.
               do_unzoom  = .FALSE.
               do_unzoom2 = .FALSE.
               if( string(1:1) .eq. 'Z' ) then
                  if( string(2:2) .eq. '2' ) then
                     do_zoom  = .TRUE.
                     do_zoom2 = .TRUE.
                  else
                     call display_prompt(
     1     'Enter lower left, upper right. ' // PM_CANCEL_BUTTON //
     1     ' button cancels input')
                  endif
                  do i=1,2
                     if( comp_h_used(i) ) then
                        if( do_hlog(i) ) then
                           h1(i) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                            comp_h_min(i))))
                           h2(i) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                            comp_h_max(i))))
                        else
c                          Unscale any remote exponent
                           if( comp_h_exponent(i) .ne. 0 ) then
                              h1(i) = comp_h_min(i)*
     1                                10.0**comp_h_exponent(i)
                              h2(i) = comp_h_max(i)*
     1                                10.0**comp_h_exponent(i)
                           else
                              h1(i) = comp_h_min(i)
                              h2(i) = comp_h_max(i)
                           endif
                        endif
                     endif
                     if( comp_v_used(i) ) then
                        if( do_vlog(i) ) then
                           h3(i) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                            vmin_plot(i))))
                           h4(i) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                            vmax_plot(i))))
                        else
c                          Unscale any remote exponent
                           if( v_exponent(i) .ne. 0 ) then
                              h3(i) = vmin_plot(i)*10.0**v_exponent(i)
                              h4(i) = vmax_plot(i)*10.0**v_exponent(i)
                           else
                              h3(i) = vmin_plot(i)
                              h4(i) = vmax_plot(i)
                           endif
                        endif
                     endif
                  end do
c                 Get plot area in NDC space
                  call pm_get_compare_plot_ndc_space(xndc0, yndc0,
     1                 xndc1, yndc1)
                  if( do_zoom2 ) then
c                    Do a 2X zoom
                     delta = 0.25*(h2(1) - h1(1))
                     x1(1) = h1(1) + delta
                     x2(1) = h2(1) - delta
                     delta = 0.25*(h4(1) - h3(1))
                     y1(1) = h3(1) + delta
                     y2(1) = h4(1) - delta
                     do_zoom = .TRUE.
                  else
                     call window_zoom(current_window_id,
     1                  h1(1), h2(1), h3(1), h4(1),
     2                  xndc0, xndc1, yndc0, yndc1, .FALSE.,
     3                  x1(1), x2(1), y1(1), y2(1), do_zoom)
                  endif
                  if( do_zoom .and. do_draw_to_scale ) then
                     ndc_aspect_ratio = (xndc1 - xndc0)/(yndc1 - yndc0)
                     if( ratio_y .eq. 0.0 ) then
c                       Don't have aspect ratio information yet,
c                       Use approximate ratio for PostScript device
                        aspect_ratio = ndc_aspect_ratio*1.30
                     else
                        aspect_ratio=ndc_aspect_ratio*(ratio_x/ratio_y)
                     endif
                     call pm_set_plot_bounds_to_scale(aspect_ratio,'E',
     1                  x1(1), x2(1), rdum,  idum,
     1                  y1(1), y2(1), rdum2, idum2)
                  endif
c                 Scale second axes proportionally to first axes
                  if( do_zoom ) then
                     if( comp_h_used(2) ) then
                        call scale_like_zoom(h1,h2,x1,x2)
                     endif
                     if( comp_v_used(2) ) then
                        call scale_like_zoom(h3,h4,y1,y2)
                     endif
                  endif
                  do i=1,2
                     if( do_zoom ) then
                        if( comp_h_used(i) .and. do_hlog(i) ) then
                           x1(i) = 10.0**(MIN(x1(i),38.0))
                           x2(i) = 10.0**(MIN(x2(i),38.0))
                        endif
                        if( comp_v_used(i) .and. do_vlog(i) ) then
                           y1(i) = 10.0**(MIN(y1(i),38.0))
                           y2(i) = 10.0**(MIN(y2(i),38.0))
                        endif
                     endif
                  enddo
               else
                  if( string(2:2) .eq. '2' ) then
c                    Do a 2x unzoom
                     do i=1,2
                        if( comp_h_used(i) ) then
                           delta = 0.5*(comp_h_max(i) - comp_h_min(i))
                           x1(i) = comp_h_min(i) - delta
                           x2(i) = comp_h_max(i) + delta
                           x1(i) = x1(i)*10.0**comp_h_exponent(i)
                           x2(i) = x2(i)*10.0**comp_h_exponent(i)
                        endif
                        if( comp_v_used(i) ) then
                           delta = 0.5*(vmax_plot(i) - vmin_plot(i))
                           y1(i) = vmin_plot(i) - delta
                           y2(i) = vmax_plot(i) + delta
                           y1(i) = y1(i)*10.0**v_exponent(i)
                           y2(i) = y2(i)*10.0**v_exponent(i)
                        endif
                     enddo
                     do_unzoom2 = .TRUE.
                  else
c                    Use full plot limits that we computed earlier
                     do i=1,2
                        if( comp_h_used(i) ) then
                           x1(i) = hmin_full(i)
                           x2(i) = hmax_full(i)
                           if( plot_mode .eq. 2 ) then
                              x1(i) = x1(i) - 1.0
                              x2(i) = x2(i) + 1.0
                           endif
                           comp_h_exponent(i) = 0
                        endif
                        if( comp_v_used(i) ) then
                           y1(i) = vmin_full(i)
                           y2(i) = vmax_full(i)
                        endif
                     enddo
                     if( do_draw_to_scale ) then
                        ndc_aspect_ratio=(xndc1 - xndc0)/(yndc1 - yndc0)
                        if( ratio_y .eq. 0.0 ) then
c                          Don't have aspect ratio information yet,
c                          Use approximate ratio for PostScript device
                           aspect_ratio = ndc_aspect_ratio*1.30
                        else
                           aspect_ratio = ndc_aspect_ratio*
     1                                    (ratio_x/ratio_y)
                        endif
                        call pm_set_plot_bounds_to_scale(
     1                            aspect_ratio,'E',
     1                            x1(1), x2(1), rdum,  idum,
     1                            y1(1), y2(1), rdum2, idum2)
                     endif
                  endif
                  do_unzoom = .TRUE.
               endif
               if( do_zoom .or. do_unzoom ) then
                  do i=1,2
                     if( comp_h_used(i) ) then
                        comp_h_min(i)  = x1(i)
                        comp_h_max(i)  = x2(i)
                     endif
                     if( comp_v_used(i) ) then
                        vmin_plot(i)  = y1(i)
                        vmax_plot(i)  = y2(i)
                     endif
                  end do
c
c                 Get tic mark parameters
                  do i=1,2
                     if( comp_h_used(i) ) then
                        if( .NOT. do_hlog(i) ) then
                           call get_default_limits(
     1                     comp_h_min(i), comp_h_max(i),
     2                     rdum, rdum2,
     3                     comp_h_tic_major(i),comp_h_minor_tic_freq(i),
     4                     comp_h_exponent(i), .FALSE.)
                           if( x_char_data(i) .or.
     1                         plot_mode .eq. 2 ) then
                              comp_h_tic_major(i)      = 1.0
                              comp_h_minor_tic_freq(i) = 1
                           endif
c                          Handle the case of unzoom
                           if( comp_h_exponent(i) .ne. 0 ) then
                              comp_h_min(i)=comp_h_min(i)/
     1                                      10.0**(comp_h_exponent(i))
                              comp_h_max(i)=comp_h_max(i)/
     1                                      10.0**(comp_h_exponent(i))
                           endif
                        endif
                     endif
                     if( comp_v_used(i) ) then
                        if( do_unzoom .and. .not. do_unzoom2 .and.
     1                      .not. do_draw_to_scale ) then
c                          Recompute all vertical axis parameters
                           do_vuscale(i) = .FALSE.
                           iuser_vscale(i) = 0
                           call recompute_minmax_comp(
     1                        rdum, rdum2,
     2                        vmin_plot(i), vmax_plot(i),
     3                        v_tic_major(i), v_minor_tic_freq(i),
     4                        iuser_vscale(i), v_exponent(i),
     5                        do_vabs(i), do_vlog(i), do_vuscale(i), i,
     6                        work_vec(1,1), plot_mode)
                        else if( .NOT. do_vlog(i) ) then
c                          Just get updated tic mark info
                           call get_default_limits(
     1                        vmin_plot(i), vmax_plot(i),
     2                        rdum, rdum2,
     3                        v_tic_major(i), v_minor_tic_freq(i),
     4                        v_exponent(i), .FALSE.)
                           if( v_exponent(i) .ne. 0 ) then
                           vmin_plot(i)=vmin_plot(i)/10.0**v_exponent(i)
                           vmax_plot(i)=vmax_plot(i)/10.0**v_exponent(i)
                           endif
                        endif
                     endif
                  end do
c
c                 Stuff axis information into save structures
                  call copy_axis_to_save(PLOT_COMPARE_COMMAND,
     1               comp_h_min(1), comp_h_max(1),
     2               comp_h_exponent(1), comp_h_tic_major(1),
     3               comp_h_minor_tic_freq(1),
     4               comp_hunits(1), comp_horiz_label(1),
     5               do_hlog(1), do_habs(1), comp_h_used(1),
     6               vmin_plot(1), vmax_plot(1),
     7               v_exponent(1), v_tic_major(1),
     8               v_minor_tic_freq(1),
     9               comp_vunits(1), comp_vert_label(1),
     A               do_vlog(1), do_vabs(1), comp_v_used(1),
     B               rdum, rdum,
     C               idum, rdum, idum, cdum,
     D               cdum, .FALSE., .FALSE., .FALSE.)
                  call copy_axis2_to_save(comp_h_min(2), comp_h_max(2),
     1               comp_h_exponent(2), comp_h_tic_major(2),
     2               comp_h_minor_tic_freq(2),
     3               comp_hunits(2), comp_horiz_label(2),
     4               do_hlog(2), do_habs(2), comp_h_used(2),
     5               vmin_plot(2), vmax_plot(2),
     6               v_exponent(2), v_tic_major(2),
     7               v_minor_tic_freq(2),
     8               comp_vunits(2), comp_vert_label(2),
     9               do_vlog(2), do_vabs(2), comp_v_used(2))
c
c                 Stuff data information into save structures
                  call copy_compare_to_save(plot_mode)
c
c                 Copy data into temp store, take absolute value,
c                 log if necessary...
                  do 140 i=1,n_cuts 
                     slot = slot_index(i)
                     itag = multi_1d(slot)%yaxis_tag
                     do 150 j=1,multi_1d(slot)%np
                        if( do_vabs(itag) .or. do_vlog(itag) ) then
                           work_vec(j,slot) =
     1                        ABS(multi_1d(slot)%data(j))
                           if( do_vlog(itag) ) then
c                             check for log(0.0), then fudge
                              if( work_vec(j,slot) .eq. 0.0 )
     1                           work_vec(j,slot) = PM_MIN_FLOAT
                              work_vec(j,slot) =
     1                           ALOG10(work_vec(j,slot))
                           endif
                        else
                           work_vec(j,slot) = multi_1d(slot)%data(j)
                        endif
150                  continue
140               continue
c
                  do i=1,2
                     if( do_vlog(i) ) then
                        v_min(i) = ALOG10(vmin_plot(i))
                        v_max(i) = ALOG10(vmax_plot(i))
                     else
                        v_min(i) = vmin_plot(i)
                        v_max(i) = vmax_plot(i)
                     endif
                  end do
c
                  do i=1,2
                     if( comp_h_used(i) ) then
                        if( do_hlog(i) ) then
                           h_min(i) = ALOG10(comp_h_min(i))
                           h_max(i) = ALOG10(comp_h_max(i))
                        else
                           h_min(i) = comp_h_min(i)
                           h_max(i) = comp_h_max(i)
                        endif
                     endif
                  end do
c
                  call do_multi_plot(work_vec,
     1               h_min, h_max, comp_h_exponent,
     2               comp_h_tic_major(1), comp_h_minor_tic_freq,
     3               v_min, v_max, v_exponent, 
     4               v_tic_major, v_minor_tic_freq,
     5               comp_hunits, comp_vunits,
     6               comp_horiz_label, comp_vert_label, do_hlog,do_vlog,
     7               comp_h_used, comp_v_used, do_header, plot_mode,
     8               lsave, lopen, lannotate, ierror)
                  if( ierror .eq. PLT$NORMAL ) then
                     lopen = .TRUE.
                  else
                     call plt$get_error_message(ierror,string)
                     delayed_message = .TRUE.
                  endif
               endif
            endif
            goto 4000
         else
            call display_message(
     1         'Plot must be displayed before you can zoom it')
         endif
         goto 2000
c
      else if( string(1:1) .eq. 'H' ) then
         fake_ln03 = .TRUE.
         call inquire_hardcopy_menu(fake_ln03, abort)
         if( abort ) goto 4000
         previous_plot_device = plot_device
         plot_device = hardcopy_device
         call display_prompt('Making hardcopy')
c
c        Stuff axis information into save structures
         call copy_axis_to_save(PLOT_COMPARE_COMMAND,
     1      comp_h_min(1), comp_h_max(1),
     2      comp_h_exponent(1), comp_h_tic_major(1),
     3      comp_h_minor_tic_freq(1),
     4      comp_hunits(1), comp_horiz_label(1),
     5      do_hlog(1), do_habs(1), comp_h_used(1),
     6      vmin_plot(1), vmax_plot(1),
     7      v_exponent(1), v_tic_major(1), v_minor_tic_freq(1),
     8      comp_vunits(1), comp_vert_label(1),
     9      do_vlog(1), do_vabs(1), comp_v_used(1),
     A      rdum, rdum,
     B      idum, rdum, idum, cdum,
     C      cdum, .FALSE., .FALSE., .FALSE.)
         call copy_axis2_to_save(comp_h_min(2), comp_h_max(2),
     1      comp_h_exponent(2), comp_h_tic_major(2),
     2      comp_h_minor_tic_freq(2),
     3      comp_hunits(2), comp_horiz_label(2),
     4      do_hlog(2), do_habs(2), comp_h_used(2),
     5      vmin_plot(2), vmax_plot(2),
     6      v_exponent(2), v_tic_major(2),
     7      v_minor_tic_freq(2),
     8      comp_vunits(2), comp_vert_label(2),
     9      do_vlog(2), do_vabs(2), comp_v_used(2))
c
c        Stuff data information into save structures
         call copy_compare_to_save(plot_mode)
c
c        Copy data into temp store, take absolute value, log if necessary...
         do 120 i=1,n_cuts
            slot = slot_index(i)
            itag = multi_1d(slot)%yaxis_tag
            do 130 j=1,multi_1d(slot)%np
               if( do_vabs(itag) .or. do_vlog(itag) ) then
                  work_vec(j,slot) = ABS(multi_1d(slot)%data(j))
                  if( do_vlog(itag) ) then
c                    check for log(0.0)
c                    A fudge
                     if( work_vec(j,slot) .eq. 0.0 )
     1                  work_vec(j,slot) = PM_MIN_FLOAT
                     work_vec(j,slot) = ALOG10(work_vec(j,slot))
                  endif
               else
                  work_vec(j,slot) = multi_1d(slot)%data(j)
               endif
130         continue
120      continue
c
         do i=1,2
            if( do_vlog(i) ) then
               v_min(i) = ALOG10(vmin_plot(i))
               v_max(i) = ALOG10(vmax_plot(i))
            else
               v_min(i) = vmin_plot(i)
               v_max(i) = vmax_plot(i)
            endif
         end do
c
         do i=1,2
            if( comp_h_used(i) ) then
               if( do_hlog(i) ) then
                  h_min(i) = ALOG10(comp_h_min(i))
                  h_max(i) = ALOG10(comp_h_max(i))
               else
                  h_min(i) = comp_h_min(i)
                  h_max(i) = comp_h_max(i)
               endif
            endif
         end do
c
         call do_multi_plot(work_vec, h_min, h_max,
     1      comp_h_exponent,
     2      comp_h_tic_major, comp_h_minor_tic_freq,
     3      v_min, v_max,
     4      v_exponent, v_tic_major, v_minor_tic_freq,
     5      comp_hunits, comp_vunits,
     6      comp_horiz_label, comp_vert_label, do_hlog, do_vlog,
     7      comp_h_used, comp_v_used, do_header, plot_mode,
     8      lsave, lopen, lannotate, ierror)
c
         plot_device = previous_plot_device
         if( ierror .eq. PLT$NORMAL ) then
            ilen = string_len(plot_file_name)
            write(string,200,err=201) plot_file_name(1:ilen)
200         format('Plot file ',a' created')
201         continue
         else
            call plt$get_error_message(ierror,string)
         endif
         delayed_message = .TRUE.
         goto 4000
c
      else if( string(1:1) .eq. 'A' ) then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: cannot annotate on this graphics device')
            else
               call anno_annotate_graph(current_window_id)
               lannotate = .TRUE.
               goto 4000
            endif
         else
            call display_message(
     1         'Plot must be displayed before you can annotate it')
         endif
         goto 2000
c     Special case to catch E or Q (which are "legal" floating point
c     numbers!)
      else if( string(1:1) .eq. 'E' .or. string(1:1) .eq. 'Q') then
         goto 5000
      else if( string(1:1) .eq. 'R') then ! Repaint screen
         goto 4000
      else if( string(1:1) .eq. 'K') then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: no window functions on this graphics device')
               goto 2000
            else
c              Set "save this window" flag
c              (only for workstations and when using retained window style)
               call retain_current_window
            endif
         else
            call display_message(
     1      'Plot must be displayed in a window before you can keep it')
         endif
         goto 2000
      else if( string(1:1) .eq. 'S') then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: can not sample on this graphics device')
               goto 2000
            else
               call display_prompt(
     1         'Enter sample point with left mouse button. ' //
     2         PM_CANCEL_BUTTON // ' button cancels input')
               first_probe = .TRUE.
c              Top of sample loop...
720            continue
c              Determine which axes to use, depending on current axes
c              displayed to user
               if( comp_h_used(ipage) ) then
                  ixaxis = ipage
               else
                  ixaxis = MOD(ipage,2)+1
               endif
               if( comp_v_used(ipage) ) then
                  iyaxis = ipage
               else
                  iyaxis = MOD(ipage,2)+1
               endif
               if( do_hlog(ixaxis) ) then
                  h1(1) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                   comp_h_min(ixaxis))))
                  h2(1) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                   comp_h_max(ixaxis))))
               else
                  h1(1) = comp_h_min(ixaxis)
                  h2(1) = comp_h_max(ixaxis)
               endif
               if( do_vlog(iyaxis) ) then
                  h3(1) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                   vmin_plot(iyaxis))))
                  h4(1) = ALOG10(ABS(MAX(PM_MIN_FLOAT,
     1                                   vmax_plot(iyaxis))))
               else
                  h3(1) = vmin_plot(iyaxis)
                  h4(1) = vmax_plot(iyaxis)
               endif
c              Set plot area in NDC space
               call pm_get_compare_plot_ndc_space(xndc0, yndc0,
     1              xndc1, yndc1)
               call get_probe_location(current_window_id,
     1            h1(1), h2(1), h3(1), h4(1),
     2            xndc0, xndc1, yndc0, yndc1, .FALSE.,
     3            x1(1), y1(1), do_probe)
               if( do_probe ) then
                  if( do_hlog(ixaxis) ) then
                     x1(1) = 10.0**(MIN(x1(1),38.0))
                  else if(comp_h_exponent(ixaxis) .ne. 0) then
                     x1(1) = x1(1)*10.0**(comp_h_exponent(ixaxis))
                  endif
                  if( do_vlog(iyaxis) ) then
                     y1(1) = 10.0**(MIN(y1(1),38.0))
                  else if(v_exponent(iyaxis) .ne. 0) then
                     y1(1) = y1(1)*10.0**(v_exponent(iyaxis))
                  endif
               endif
               if( .NOT. do_probe ) goto 710
               if( first_probe ) then
                  write(string,700,err=702) x1(1), y1(1)
700               format('Probe at (',1pe12.5,',',1pe12.5,')')
                  first_probe = .FALSE.
                  lastx = x1(1)
                  lasty = y1(1)
               else
                  delx = x1(1) - lastx
                  dely = y1(1) - lasty
                  write(string,701,err=702) x1(1), y1(1), delx, dely
701               format('Probe at (',1pe12.5,',',1pe12.5,') change: (',
     1                   1pe12.5,',',1pe12.5,')')
                  lastx = x1(1)
                  lasty = y1(1)
               endif
702            continue
               call display_message_quietly(string)
               goto 720
710            continue
            endif
         else
            call display_message(
     1      'Plot must be displayed in a window before you can sample')
         endif
         goto 2000
      else
         read(string,*,err=1000) item_number
         if( item_number .eq. 1 ) then
            call get_string_with_prompt(title_dat(1), 
     1              'Enter title (may use mathmode)')
            num_titles = MAX(1,num_titles)
            call display_string(1 , 12, title_dat(1), 69, 1)
c
         elseif ( item_number .eq. 2 ) then
c           X axis min/max
            real_vec(1) = comp_h_min(ipage)
            real_vec(2) = comp_h_max(ipage)
            call get_reals_with_prompt(real_vec, 2,
     1              'Enter min and max values for horizontal axis')
            if( real_vec(1) .ge. real_vec(2) ) then
               call display_message(
     1                 'Min value must be less than max')
            else if( do_hlog(ipage) .and.
     1         (real_vec(1) .le. 0.0 .or. real_vec(2) .le. 0.0) ) then
               call display_message(
     1                 'Log limits must be greater than zero')
            else
               comp_h_min(ipage) = real_vec(1)
               comp_h_max(ipage) = real_vec(2)
               if( do_draw_to_scale .and. .NOT. do_hlog(ipage) .and.
     1             .NOT. do_vlog(ipage) ) then
c                 Get plot aspect ratio
                  call pm_get_compare_plot_ndc_space(xndc0, yndc0,
     1                    xndc1, yndc1)
                  ndc_aspect_ratio = (xndc1 - xndc0)/(yndc1 - yndc0)
                  if( ratio_y .eq. 0.0 ) then
c                    Don't have aspect ratio information yet,
c                    Use approximate ratio for PostScript device
                      aspect_ratio = ndc_aspect_ratio*1.30
                  else
                     aspect_ratio = ndc_aspect_ratio*(ratio_x/ratio_y)
                  endif
c                 Reset Y plot bounds to maintain aspect ratio
                  call pm_set_plot_bounds_to_scale(aspect_ratio,'X',
     1                     comp_h_min(ipage),comp_h_max(ipage),
     2                     comp_h_tic_major(ipage),
     1                     comp_h_minor_tic_freq(ipage),
     3                     vmin_plot(ipage),vmax_plot(ipage),
     4                     v_tic_major(ipage),v_minor_tic_freq(ipage))
c                 Display Y axis parameters
                  call display_real(14, 15, vmin_plot(ipage), 1)
                  call display_real(14, 28, vmax_plot(ipage), 1)
                  call display_real(15, 23, v_tic_major(ipage), 1)
                  call display_int (16, 27, v_minor_tic_freq(ipage), 5)
                  call display_int (17, 29, v_exponent(ipage), 3)
               endif
               call display_real(6 , 15, comp_h_min(ipage), 1)
               call display_real(6 , 28, comp_h_max(ipage), 1)
               if( .NOT. do_hlog(ipage) ) then
c                 Generate compatible tic mark information
                  call get_default_limits(
     1                 comp_h_min(ipage), comp_h_max(ipage),
     2                 rdum, rdum2,
     3                 comp_h_tic_major(ipage),
     4                 comp_h_minor_tic_freq(ipage),
     5                 idum, .FALSE.)
                  comp_h_tic_major(ipage) =
     1               comp_h_tic_major(ipage)*10.0**idum
                  if( x_char_data(ipage) .or. plot_mode .eq. 2 ) then
                     comp_h_tic_major(ipage) = 1.0
                     comp_h_minor_tic_freq(ipage) = 1
                  endif
                  call display_real(7 , 23,comp_h_tic_major(ipage), 1)
                  call display_int (8 , 27,comp_h_minor_tic_freq(ipage),
     1                              5)
               endif
            endif
c
         elseif ( item_number .eq. 3 ) then
c           X axis major tic increment
            real_vec(1) = comp_h_tic_major(ipage)
            if( do_hlog(ipage) ) then
               call display_message(
     1            'Not used when plotting on log scale')
            else
               call get_reals_with_prompt(real_vec, 1,
     1                 'Enter major tic increment for horizontal axis')
               delta = comp_h_max(ipage) - comp_h_min(ipage)
               if( real_vec(1) .gt. delta .or.
     1            real_vec(1) .le. 0.0 ) then
                  call display_message(
     1                    'Increment < 0 or larger than plot range')
               else
                  comp_h_tic_major(ipage) = real_vec(1)
                  call display_real(7 , 23, comp_h_tic_major(ipage), 1)
               endif
           endif
c
         elseif ( item_number .eq. 4 ) then
c           Minor tic frequency
            int_vec(1) = comp_h_minor_tic_freq(ipage)
            if( do_hlog(ipage) ) then
               call display_message(
     1            'Not used when plotting on log scale')
            else
               call get_ints_with_prompt(int_vec, 1,
     1                 'Enter minor tic frequency for horizontal axis')
               if( int_vec(1) .lt. 1 ) then
                  call display_message('Illegal value')
               else
                  comp_h_minor_tic_freq(ipage) = int_vec(1)
                  call display_int (8 , 27, comp_h_minor_tic_freq(ipage)
     1                              ,5)
               endif
            endif
c
         elseif ( item_number .eq. 5 ) then
c           X axis scale factor
            call display_message('Sorry: not implemented yet')
c
         elseif ( item_number .eq. 6 ) then
c           X axis log
            do_hlog(ipage) = .NOT. do_hlog(ipage)
            if( do_hlog(ipage) ) then
               prev_habs(ipage) = do_habs(ipage)
               do_habs(ipage) = .TRUE.
               str_hlog = 'log'
               str_habs = 'yes'
c              Save previous values
               prev_h_exponent(ipage) = comp_h_exponent(ipage)
               prev_hmin_plot(ipage) = comp_h_min(ipage)
               prev_hmax_plot(ipage) = comp_h_max(ipage)
               prev_h_tic_major(ipage)  = comp_h_tic_major(ipage)
               prev_h_minor_tic_freq(ipage) =
     1            comp_h_minor_tic_freq(ipage)
               prev_h_saved(ipage) = .TRUE.
               call recompute_minmax_coords(
     1            hmin_full(ipage), hmax_full(ipage),
     2            comp_h_min(ipage), comp_h_max(ipage),
     3            comp_h_tic_major(ipage), comp_h_minor_tic_freq(ipage),
     4            iuser_hscale(ipage), comp_h_exponent(ipage),
     5            do_habs(ipage), do_hlog(ipage), do_huscale(ipage),
     6            ipage, work_vec(1,1))
               comp_h_exponent(ipage) = 0
c              Mark X axis tics as not applicable
               call display_string(7 , 23, '       NA', 9, 1)
               call display_string(8 , 27, '   NA', 5, 1)
            else
               do_habs(ipage) = prev_habs(ipage)
               str_hlog = 'linear'
               if( do_habs(ipage) ) then
                  str_habs = 'yes'
               else
                  str_habs = 'no'
               endif
c              Recompute hmin_full/hmax_full
               call recompute_minmax_coords(
     1            hmin_full(ipage), hmax_full(ipage),
     1            comp_h_min(ipage), comp_h_max(ipage),
     2            comp_h_tic_major(ipage), comp_h_minor_tic_freq(ipage),
     3            iuser_hscale(ipage), comp_h_exponent(ipage),
     4            do_habs(ipage), do_hlog(ipage), do_huscale(ipage),
     5            ipage, work_vec(1,1))
c              Restore previous axis values
               if( prev_h_saved(ipage) ) then
                  comp_h_exponent(ipage) = prev_h_exponent(ipage)
                  comp_h_min(ipage) = prev_hmin_plot(ipage)
                  comp_h_max(ipage) = prev_hmax_plot(ipage)
                  comp_h_tic_major(ipage) = prev_h_tic_major(ipage)
                  comp_h_minor_tic_freq(ipage) =
     1               prev_h_minor_tic_freq(ipage)
               endif
c              Re-display tic/scale data
               call display_real (7 , 23, comp_h_tic_major(ipage), 1)
               call display_int  (8 , 27, comp_h_minor_tic_freq(ipage),
     1                            5)
            endif
            call display_real  (4 , 15, hmin_full(ipage), 1)
            call display_real  (4 , 28, hmax_full(ipage), 1)
            call display_real  (6 , 15, comp_h_min(ipage), 1)
            call display_real  (6 , 28, comp_h_max(ipage), 1)
            call display_int   (9 , 29, comp_h_exponent(ipage), 3)
            call display_string(4 , 59, str_hlog, 6, 1)
            call display_string(5 , 63, str_habs, 3, 1)
c
         elseif ( item_number .eq. 7 ) then
c           X axis absolute value (not implented)
            call display_message('Sorry: not implemented yet')
c
         elseif ( item_number .eq. 8 ) then
c           X axis recompute defaults
            do_huscale(ipage) = .FALSE.
            iuser_hscale(ipage) = 0
            call recompute_minmax_coords(
     1           rdum, rdum2,
     2           comp_h_min(ipage), comp_h_max(ipage),
     3           comp_h_tic_major(ipage), comp_h_minor_tic_freq(ipage),
     4           iuser_hscale(ipage), comp_h_exponent(ipage),
     5           do_habs(ipage), do_hlog(ipage), do_huscale(ipage),
     6           ipage, work_vec(1,1))
            if( x_char_data(ipage) .or. plot_mode .eq. 2 ) then
               comp_h_min(ipage) = comp_h_min(ipage) - 1.0
               comp_h_max(ipage) = comp_h_max(ipage) + 1.0
               comp_h_tic_major(ipage) = 1.0
               comp_h_minor_tic_freq(ipage) = 1
            endif
c           Re-display tic/scale data
            if( .NOT. do_hlog(ipage) ) then
               call display_real (7 , 23, comp_h_tic_major(ipage), 1)
               call display_int  (8 , 27, comp_h_minor_tic_freq(ipage),
     1                            5)
            endif
            call display_int   (9 , 29, comp_h_exponent(ipage), 3)
            call display_string(5 , 63, str_habs, 3, 1)
            call display_real  (6 , 15, comp_h_min(ipage), 1)
            call display_real  (6 , 28, comp_h_max(ipage), 1)
c
         elseif ( item_number .eq. 9 ) then
            call get_string_with_prompt(comp_horiz_label(ipage),
     1              'Enter horizontal label (may use mathmode)')
            call display_string(7 , 54, comp_horiz_label(ipage), 27, 1)
c
         elseif ( item_number .eq. 10 ) then
            call get_string_with_prompt(comp_hunits(ipage),
     1              'Enter units (may use mathmode)')
            call display_string(8 , 54, comp_hunits(ipage), 20, 1)
c
         elseif ( item_number .eq. 11 ) then
c           Y axis min/max
            real_vec(1) = vmin_plot(ipage)
            real_vec(2) = vmax_plot(ipage)
            call get_reals_with_prompt(real_vec, 2,
     1              'Enter min and max values for vertical axis')
            if( real_vec(1) .ge. real_vec(2) ) then
               call display_message(
     1                 'Min value must be less than max')
            else if( do_vlog(ipage) .and.
     1         (real_vec(1) .le. 0.0 .or. real_vec(2) .le. 0.0) ) then
               call display_message(
     1                 'Log limits must be greater than zero')
            else
               vmin_plot(ipage) = real_vec(1)
               vmax_plot(ipage) = real_vec(2)
               if( do_draw_to_scale .and. .NOT. do_hlog(ipage) .and.
     1             .NOT. do_vlog(ipage) ) then
c                 Reset X plot bounds to maintain aspect ratio
                  call pm_get_compare_plot_ndc_space(xndc0, yndc0,
     1                    xndc1, yndc1)
                  ndc_aspect_ratio = (xndc1 - xndc0)/(yndc1 - yndc0)
                  if( ratio_y .eq. 0.0 ) then
c                    Don't have aspect ratio information yet,
c                    Use approximate ratio for PostScript device
                      aspect_ratio = ndc_aspect_ratio*1.30
                  else
                     aspect_ratio = ndc_aspect_ratio*(ratio_x/ratio_y)
                  endif
                  call pm_set_plot_bounds_to_scale(aspect_ratio,'Y',
     1                     comp_h_min(ipage),comp_h_max(ipage),
     2                     comp_h_tic_major(ipage),
     1                     comp_h_minor_tic_freq(ipage),
     3                     vmin_plot(ipage),vmax_plot(ipage),
     4                     v_tic_major(ipage),v_minor_tic_freq(ipage))
c                 Display X axis parameters
                  call display_real(6 ,15,comp_h_min(ipage), 1)
                  call display_real(6 ,28,comp_h_max(ipage), 1)
                  call display_real(7 ,23,comp_h_tic_major(ipage), 1)
                  call display_int (8 ,27,comp_h_minor_tic_freq(ipage),
     1                              5)
                  call display_int (9 ,29, comp_h_exponent(ipage), 3)
               endif
               call display_real(14, 15, vmin_plot(ipage), 1)
               call display_real(14, 28, vmax_plot(ipage), 1)
               if( .NOT. do_vlog(ipage) ) then
c                 Generate compatible tic mark information
                  call get_default_limits(
     1                 vmin_plot(ipage), vmax_plot(ipage),
     2                 rdum, rdum2,
     3                 v_tic_major(ipage), v_minor_tic_freq(ipage),
     4                 idum, .FALSE.)
                  v_tic_major(ipage) = v_tic_major(ipage)*10.0**idum
                  call display_real(15, 23, v_tic_major(ipage), 1)
                  call display_int (16, 27, v_minor_tic_freq(ipage),
     1                              5)
               endif
            endif
c
         elseif ( item_number .eq. 12 ) then
c           Major Y tic increment
            real_vec(1) = v_tic_major(ipage)
            if( do_vlog(ipage) ) then
               call display_message(
     1            'Not used when plotting on log scale')
            else
               call get_reals_with_prompt(real_vec, 1,
     1              'Enter major tic increment for vertical axis')
               delta = vmax_plot(ipage) - vmin_plot(ipage)
               if( real_vec(1) .gt. delta .or.
     1             real_vec(1) .le. 0.0 ) then
                   call display_message('Illegal value')
               else
                  v_tic_major(ipage) = real_vec(1)
                  call display_real(15, 23, v_tic_major(ipage), 1)
               endif
            endif
c
         elseif ( item_number .eq. 13 ) then
c           Minor tic frequency
            int_vec(ipage) = v_minor_tic_freq(ipage)
            if( do_vlog(ipage) ) then
               call display_message(
     1            'Not used when plotting on log scale')
            else
               call get_ints_with_prompt(int_vec, 1,
     1              'Enter minor tic frequency for vertical axis')
               if( int_vec(1) .lt. 1 ) then
                 call display_message('Illegal value')
               else
                  v_minor_tic_freq(ipage) = int_vec(1)
                  call display_int (16, 27, v_minor_tic_freq(ipage),
     1                              5)
               endif
            endif
c
         elseif ( item_number .eq. 14 ) then
c           Y axis scale factor
            if( do_vlog(ipage) ) then
               call display_message(
     1            'Cannot use scale factor on log plot')
            else
               int_vec(1) = v_exponent(ipage)
               call get_ints_with_prompt(int_vec, 1,
     1      'Enter integer power of ten that will scale (divide) data')
               if( int_vec(1) .gt. 38 .or.
     1            int_vec(1) .lt. -38 ) then
                  call display_message('Illegal value')
               else
                  iuser_vscale(ipage) = int_vec(1)
                  do_vuscale(ipage) = .TRUE.
                  h = 10.0**iuser_vscale(ipage)/10.0**v_exponent(ipage)
                  vmin_plot(ipage) = vmin_plot(ipage)/h
                  vmax_plot(ipage) = vmax_plot(ipage)/h
                  v_tic_major(ipage) = v_tic_major(ipage)/h
                  v_exponent(ipage) = iuser_vscale(ipage)
c                 Re-display tic/scale data
                  call display_real(14, 15, vmin_plot(ipage), 1)
                  call display_real(14, 28, vmax_plot(ipage), 1)
                  call display_int (17, 29, v_exponent(ipage), 3)
                  call display_real(15, 23, v_tic_major(ipage), 1)
                  call display_int (16, 27, v_minor_tic_freq(ipage), 5)
               endif
            endif
c
         elseif ( item_number .eq. 15 ) then
c           Logarithm
            do_vlog(ipage) = .NOT. do_vlog(ipage)
            if( do_vlog(ipage) ) then
               prev_vabs(ipage) = do_vabs(ipage)
               do_vabs(ipage)   = .TRUE.
               str_vlog = 'log'
               str_vabs = 'yes'
c              Save previous values
               prev_v_exponent(ipage) = v_exponent(ipage)
               prev_iuser_vscale(ipage) = iuser_vscale(ipage)
               prev_vmin_plot(ipage) = vmin_plot(ipage)
               prev_vmax_plot(ipage) = vmax_plot(ipage)
               prev_v_tic_major(ipage)  = v_tic_major(ipage)
               prev_v_minor_tic_freq(ipage) = v_minor_tic_freq(ipage)
               prev_v_saved(ipage) = .TRUE.
               iuser_vscale(ipage) = 0
c
               call recompute_minmax_comp(
     1                   vmin_full(ipage), vmax_full(ipage),
     2                   vmin_plot(ipage), vmax_plot(ipage),
     3                   v_tic_major(ipage), v_minor_tic_freq(ipage),
     4                   iuser_vscale(ipage), v_exponent(ipage),
     5                   do_vabs(ipage), do_vlog(ipage),
     6                   do_vuscale(ipage), ipage, work_vec(1,1),
     7                   plot_mode)
               if( iuser_vscale(ipage) .ne. 0 )
     1            v_exponent(ipage) = iuser_vscale(ipage)
c              Mark Y axis tics as not applicable
               call display_string(15, 23, '       NA', 9, 1)
               call display_string(16, 27, '   NA', 5, 1)
            else
               do_vabs(ipage) = prev_vabs(ipage)
               str_vlog = 'linear'
               if( do_vabs(ipage) ) then
                  str_vabs = 'yes'
               else
                  str_vabs = 'no'
               endif
c              Recompute vmin_full/vmax_full
               call recompute_minmax_comp(
     1                   vmin_full(ipage), vmax_full(ipage),
     2                   vmin_plot(ipage), vmax_plot(ipage),
     3                   v_tic_major(ipage), v_minor_tic_freq(ipage),
     4                   iuser_vscale(ipage), v_exponent(ipage),
     5                   do_vabs(ipage), do_vlog(ipage),
     6                   do_vuscale(ipage), ipage, work_vec(1,1),
     7                   plot_mode)
c              Restore previous axis values
               if( prev_v_saved(ipage) ) then
                  iuser_vscale(ipage) = prev_iuser_vscale(ipage)
                  v_exponent(ipage) = prev_v_exponent(ipage)
                  vmin_plot(ipage) = prev_vmin_plot(ipage)
                  vmax_plot(ipage) = prev_vmax_plot(ipage)
                  v_tic_major(ipage) = prev_v_tic_major(ipage)
                  v_minor_tic_freq(ipage) = prev_v_minor_tic_freq(ipage)
               endif
               if( iuser_vscale(ipage) .ne. 0 )
     1            v_exponent(ipage) = iuser_vscale(ipage)
c              Re-display tic/scale data
               call display_real  (15, 23, v_tic_major(ipage), 1)
               call display_int   (16, 27, v_minor_tic_freq(ipage), 5)
            endif
            call display_real  (12, 15, vmin_full(ipage), 1)
            call display_real  (12, 28, vmax_full(ipage), 1)
            call display_real  (14, 15, vmin_plot(ipage), 1)
            call display_real  (14, 28, vmax_plot(ipage), 1)
            call display_int   (17, 29, v_exponent(ipage), 3)
            call display_string(12, 59, str_vlog, 6, 1)
            call display_string(13, 63, str_vabs, 3, 1)
c
c           Absolute value
         elseif ( item_number .eq. 16 ) then
            if( do_vlog(ipage) .AND. do_vabs(ipage) ) then
               call display_message(
     1            'Absolute value required for log scale')
            else
               prev_vabs(ipage) = do_vabs(ipage)
               do_vabs(ipage) = .NOT. do_vabs(ipage)
               if( do_vabs(ipage) ) then
                  str_vabs = 'yes'
               else
                  str_vabs = 'no'
               endif
            endif
            call recompute_minmax_comp(
     1                vmin_full(ipage), vmax_full(ipage),
     2                vmin_plot(ipage), vmax_plot(ipage),
     3                v_tic_major(ipage), v_minor_tic_freq(ipage),
     4                iuser_vscale(ipage), v_exponent(ipage),
     5                do_vabs(ipage), do_vlog(ipage),
     6                do_vuscale(ipage), ipage, work_vec(1,1),
     7                plot_mode)
            if(do_vuscale(ipage)) v_exponent(ipage)=iuser_vscale(ipage)
c           Re-display tic/scale data
            call display_real  (15, 23, v_tic_major(ipage), 1)
            call display_int   (16, 27, v_minor_tic_freq(ipage), 5)
            call display_int   (17, 29, v_exponent(ipage), 3)
            call display_string(13, 63, str_vabs, 3, 1)
            call display_real  (12, 15, vmin_full(ipage), 1)
            call display_real  (12, 28, vmax_full(ipage), 1)
            call display_real  (14, 15, vmin_plot(ipage), 1)
            call display_real  (14, 28, vmax_plot(ipage), 1)
c
         elseif ( item_number .eq. 17 ) then
c           Autoscale vertical axis
            do_vuscale(ipage) = .FALSE.
            iuser_vscale(ipage) = 0
            call recompute_minmax_comp(
     1                rdum, rdum2,
     2                vmin_plot(ipage), vmax_plot(ipage),
     3                v_tic_major(ipage), v_minor_tic_freq(ipage),
     4                iuser_vscale(ipage), v_exponent(ipage),
     5                do_vabs(ipage), do_vlog(ipage), do_vuscale(ipage),
     6                ipage, work_vec(1,1), plot_mode)
c           Re-display tic/scale data
            if( .NOT. do_vlog(ipage) ) then
               call display_real  (15, 23, v_tic_major(ipage), 1)
               call display_int   (16, 27, v_minor_tic_freq(ipage), 5)
            endif
            call display_int   (17, 29, v_exponent(ipage), 3)
            call display_string(13, 63, str_vabs, 3, 1)
            call display_real  (14, 15, vmin_plot(ipage), 1)
            call display_real  (14, 28, vmax_plot(ipage), 1)
c
         elseif ( item_number .eq. 18 ) then
            call get_string_with_prompt(comp_vert_label(ipage),
     1              'Enter vertical label (may use mathmode)')
            call display_string(15, 54, comp_vert_label(ipage), 27, 1)
c
         elseif ( item_number .eq. 19 ) then
            call get_string_with_prompt(comp_vunits(ipage),
     1              'Enter units (may use mathmode)')
            call display_string(16, 54, comp_vunits(ipage), 20, 1)
c
         elseif ( item_number .eq. 20 ) then
            comp_nolabels = .NOT. comp_nolabels
            if( comp_nolabels ) then
               str_labels = 'yes'
            else
               str_labels = 'no '
            endif
            call display_string(18, 64, str_labels, 3, 1)
         else if( item_number .eq. 21 ) then ! Switch sets of axes
            if( ipage .eq. 1 ) then
               ipage = 2
            else
               ipage = 1
            endif
            goto 4000
         elseif ( item_number .eq. 22 ) then
            if( plot_mode .eq. 2 ) then
               call display_message(
     1            'Draw to scale not applicable to bar charts')
               goto 4000
            endif
            if( do_hlog(ipage) .or. do_vlog(ipage) ) then
               call display_message(
     1         'Draw to scale not applicable when log scales are used')
               goto 4000
            endif
            do_draw_to_scale = .NOT. do_draw_to_scale
            if( do_draw_to_scale ) then
              str_aspect = 'yes'
            else
              str_aspect = ' no'
            endif
            call display_string(20,64,str_aspect, 3, 1)
            if( do_draw_to_scale ) then
c              Reset plot bounds
               call pm_get_compare_plot_ndc_space(xndc0, yndc0,
     1                 xndc1, yndc1)
               ndc_aspect_ratio = (xndc1 - xndc0)/(yndc1 - yndc0)
               if( ratio_y .eq. 0.0 ) then
c                 Don't have aspect ratio information yet,
c                 Use approximate ratio for PostScript device
                  aspect_ratio = ndc_aspect_ratio*1.30
               else
                  aspect_ratio = ndc_aspect_ratio*(ratio_x/ratio_y)
               endif
               call pm_set_plot_bounds_to_scale(aspect_ratio,'E',
     1                  comp_h_min(ipage),comp_h_max(ipage),
     2                  comp_h_tic_major(ipage),
     1                  comp_h_minor_tic_freq(ipage),
     3                  vmin_plot(ipage),vmax_plot(ipage),
     4                  v_tic_major(ipage),v_minor_tic_freq(ipage))
c              Display X axis parameters
               call display_real (6 , 15, comp_h_min(ipage), 1)
               call display_real (6 , 28, comp_h_max(ipage), 1)
               call display_real (7 , 23, comp_h_tic_major(ipage), 1)
               call display_int  (8 , 27, comp_h_minor_tic_freq(ipage),
     1                            5)
               call display_int  (9 , 29, comp_h_exponent(ipage), 3)
c              Display Y axis parameters
               call display_real(14, 15, vmin_plot(ipage), 1)
               call display_real(14, 28, vmax_plot(ipage), 1)
               call display_real(15, 23, v_tic_major(ipage), 1)
               call display_int (16, 27, v_minor_tic_freq(ipage), 5)
               call display_int (17, 29, v_exponent(ipage), 3)
            endif
         else
            call display_message('Item number out of range')
         endif
         goto 2000
c
1000     call display_message('Illegal input')
         goto 2000
      endif
c
5000  continue
c
      return
      end
c
      subroutine do_multi_plot(data_vec,
     1     h_win_min, h_win_max,
     2     h_exponent, h_tic_major, h_minor_tic_freq,
     3     v_win_min, v_win_max,
     4     v_exponent, v_tic_major, v_minor_tic_freq,
     5     hunits, vunits, horiz_label, vert_label,
     6     h_log, v_log, h_axis_used, v_axis_used,
     7     do_header, plot_mode,
     8     add_notes_from_save, ws_window_open, add_notes_from_file,
     9     ierror)
c
c     Does the dirty work of calling GKS routines to plot:
c        1) plot frame
c        2) data with markers at data points
c        3) axis tic marks
c        4) axis numeric labels
c        5) axis alpha labels
c
c     NOTE: values in calls to tune_charsize were painfully picked to give
c     readable results with VTxxx terminals. If you change the
c     size of windows where text is written, expect the characters to be
c     screwed up.
c
c     NOTE: Some axis and data information is grabbed directly from
c     the 'sticky' information arrays (e.g. comp_...). This was done to
c     support the extensions used by CURV. At some later date,
c     we may want to pass this information via the subroutine parameters.
c
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'io.cmn'
      include 'defaults.cmn'
      include 'compare.cmn'
      include 'units.cmn'
c
      logical v_log(2), h_log(2), do_header, add_notes_from_save,
     1        ws_window_open, add_notes_from_file,
     2        h_axis_used(2), v_axis_used(2)
      character*(*) horiz_label(2), vert_label(2), hunits(2), vunits(2)
      real data_vec(COMP_ARRAYSIZE,MULTI_SIZE)
      real h_win_min(2), h_win_max(2), v_win_min(2), v_win_max(2)
      real h_tic_major(2), v_tic_major(2)
      integer plot_mode, ierror
      integer h_exponent(2), h_minor_tic_freq(2)
      integer v_exponent(2), v_minor_tic_freq(2)
c
c     locals
c
      logical exp_flag, lerr, is_vtxxx, lfound, x_char_data(2)
      integer i, j, ilen, ilen2, ilen3, string_len, window_id,
     1        return_stat
      integer get_new_window_id
      integer limit, slot, nactive, ixtag, iytag
      character*1 response
      character*12 num_label, num_label1, num_label2
      character*85 work_string
      real fuzz, delta, h_cur, v_cur, h_title, v_title,
     1     ticvalue, line_length, h_delta, dely1, dely2
      real temp_vec(COMP_ARRAYSIZE), temp_coord(COMP_ARRAYSIZE),
     1     temp_z(COMP_ARRAYSIZE)
      real h_scale(2), v_scale(2), h_max, v_max, h_minor, v_minor
c
      real hscl(2), vscl(2), xndc0, yndc0, xndc1, yndc1,
     1     xextent, yextent, xoffset, yoffset, h1, h2, ZEROO, ONEE
      real rleft, rrt, tic_start(2), short_tic_end(2), long_tic_end(2),
     1     short_log_tic_end(2), long_log_tic_end(2),
     2     tic_label_pos(2), log_tic_label_pos(2), axis_label_pos(2),
     3     axis_label_pos2(2), axis_pos(2), scaled_zero, ypos
      integer k, l, m, npts, ix, iy, iaxis, icount, ncurves, maxlen
c     Two parameters close to zero and one, but within the range (0,1)
      parameter( ZEROO = 0.000001 )
      parameter( ONEE  = 0.999999 )
c
      fuzz = 0.999
      axis_pos(1) = ZEROO
      axis_pos(2) = ONEE
c
c     Scale plot to 0.0 to 1.0 (avoid bug in GKS V5.0)
      do 50 ix=1,2
c        Compute horizontal and vertical scale factors
c        Meaning of {h,v}_scale is now inverse of pre 29-Apr-93 code
         if( h_axis_used(ix) ) then
            hscl(ix) = h_win_max(ix) - h_win_min(ix)
            if( hscl(ix) .le. 0.0 ) then
               ierror = PLT$INVALID_MINMAX
            return
            endif
            h_scale(ix) = 10.0**(-h_exponent(ix))
            hscl(ix) = 1.0/hscl(ix)
         endif
         if( v_axis_used(ix) ) then
            vscl(ix) = v_win_max(ix) - v_win_min(ix)
            if( vscl(ix) .le. 0.0 ) then
               ierror = PLT$INVALID_MINMAX
            return
            endif
c
            v_scale(ix) = 10.0**(-v_exponent(ix))
            vscl(ix) = 1.0/vscl(ix)
         endif
50    continue
c
c     open GKS
      window_id = get_new_window_id(PLT$COMPARE_WINDOW_ID)
      if( window_id .eq. -1 ) then
c        This should only happen when we are using the retained
c        style, and we run out of plot window ids.
         if( lbatch .or. lrestore ) then
            call plt$get_error_message(ierror,work_string)
            ilen = string_len(work_string)
            call plt$error_message('do_multi_plot', work_string(1:ilen))
         endif
         return
      endif
c
      call open_plot_device(window_id, plot_device, 'Postmini',
     1                      return_stat)
      if( return_stat .ne. PLT$NORMAL ) then
         if( lbatch .or. lrestore ) then
            call plt$get_error_message(return_stat,work_string)
            ilen = string_len(work_string)
            call plt$error_message('do_multi_plot', work_string(1:ilen))
         endif
         ierror = return_stat
         return
      endif
c
c     Activate plot device, clear it, start segment on VAXstations
      call activate_plot_device(window_id, .TRUE., .TRUE.)
      call pop_vs_window(window_id)
      call init_device_depend
c
c     Get plot area NDC space parameters
      call pm_get_compare_plot_ndc_space(xndc0, yndc0, xndc1, yndc1)
      xextent = (xndc1 - xndc0)
      yextent = (yndc1 - yndc0)
      xoffset = 0.5*(1.0 - xextent)
      yoffset = 0.5*(1.0 - yextent)
c
      call GKS$SET_WINDOW(1, 0.0,1.0,0.0,1.0)
c     create virtual plot window (NDC space)
      call GKS$SET_VIEWPORT(1, xndc0*ratio_x, xndc1*ratio_x,
     1                         yndc0*ratio_y, yndc1*ratio_y)
      call GKS$SET_CLIPPING(GKS$K_CLIP)
      call GKS$SELECT_XFORM(1)
c
c     plot data, set window to actual data range
c
      icount  = 0
      ncurves = 0
      x_char_data(1) = .FALSE.
      x_char_data(2) = .FALSE.
c     Count the number of active curves (for barchart plot)
c     Also, find out which datasets have X axis character data
c     associated with them.
      do i=1,n_cuts
         slot = slot_index(i)
         if( multi_1d(slot)%selected .eq. 1 ) then
            ncurves = ncurves + 1
            if( multi_1d(slot)%x_char_data ) then
               x_char_data(multi_1d(slot)%xaxis_tag) = .TRUE.
            endif
         endif
      end do
c
      do 100 i=1,n_cuts
         slot = slot_index(i)
         ixtag = multi_1d(slot)%xaxis_tag
         iytag = multi_1d(slot)%yaxis_tag
         if( multi_1d(slot)%selected .ne. 1 ) goto 100
         do 101 j=1,multi_1d(slot)%np
            temp_vec(j) = (data_vec(j,slot)*v_scale(iytag) -
     1                     v_win_min(iytag))*vscl(iytag)
            if( h_log(ixtag) ) then
               temp_coord(j) = (ALOG10(MAX(PM_MIN_FLOAT,
     1                ABS(multi_1d(slot)%coord(j)*pm_scale_cm2plot) ))
     2                 *h_scale(ixtag) - h_win_min(ixtag))*hscl(ixtag)
            else
               temp_coord(j) = (multi_1d(slot)%coord(j)*pm_scale_cm2plot
     1                   *h_scale(ixtag) - h_win_min(ixtag))*hscl(ixtag)
            endif
101      continue
         if( multi_1d(slot)%sort ) then
            call pm_sort(npts, temp_coord(k), temp_vec(k))
         endif
c        Underlying window system may not clip correctly when
c        zooming in on a small portion of data. This happens most often
c        in practice when zooming on absissa. Determine minimum points
c        needed for complete plot.
         rleft = 0.0
         k = multi_1d(slot)%np+1
         do 102 j=1,multi_1d(slot)%np
            if( temp_coord(j) .ge. rleft ) then
               k = MAX(1,j-1)
               goto 103
            endif
102      continue
103      continue
         rrt = 1.0
         l = 0
         do 104 j=multi_1d(slot)%np,1,-1
            if( temp_coord(j) .le. rrt ) then
               l = MIN(multi_1d(slot)%np,j+1)
               goto 105
            endif
104      continue
105      continue
         npts = l-k+1
c
         if( plot_mode .eq. 1 ) then
c        Plot polylines
         if( multi_1d(slot)%color .ne. 0 .and.
     1       multi_1d(slot)%linestyle .ne. 0 .and. npts .gt. 0 ) then
c           First, reset linetype to solid
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
            call GKS$SET_PLINE_COLOR_INDEX(1)
c           Use call to set_line_index so that different linestyles
c           will be used on monochrome output devices to "simulate"
c           different colors.
CC Is this really required? It causes a bug to occur in the exact case
CC it is trying to help.
CCc           Special case for Curv emulation and monochrome output device.
CC            if( .NOT. (l_curv .and. ncolors .eq. 1) )
CC     1         call set_line_index(multi_1d(slot)%color)
           call set_line_index(multi_1d(slot)%color)
           if( multi_1d(slot)%linestyle .ne. 1 ) then
               call set_pline_type(multi_1d(slot)%linestyle)
            endif
            if( multi_1d(slot)%linewidth .ne. 1.0 ) then
               if( l_curv ) then
                  call GKS$SET_PLINE_LINEWIDTH(multi_1d(slot)%linewidth)
               else
                  call GKS$SET_PLINE_LINEWIDTH(
     1                 multi_1d(slot)%linewidth*hw_tuned_width)
               endif
            endif
c           Handle case of "sweep" or "contour" data
            if( multi_1d(slot)%contour ) then
c              Note, for the case of "contour" data, we can't
c              do use the simple clipping trick above, so plot entire
c              data set.
               do m=1,multi_1d(slot)%np
                  temp_z(m) = multi_1d(slot)%zval(m)
               end do
               call pm_zpolyline(multi_1d(slot)%np,
     1              temp_coord, temp_vec, temp_z)
            else if( multi_1d(slot)%sweep ) then
               call swept_polyline(npts, temp_coord(k), temp_vec(k))
            else
               call GKS$POLYLINE(npts, temp_coord(k), temp_vec(k))
            endif
            if( multi_1d(slot)%linewidth .ne. 1.0 ) then
               call GKS$SET_PLINE_LINEWIDTH(hw_tuned_width)
            endif
         endif
c        Plot markers
         if( multi_1d(slot)%marker_color .ne. 0 .and.
     1       multi_1d(slot)%marker_type  .ne. 0 .and.
     1       npts .gt. 0) then
            call GKS$SET_PMARK_COLOR_INDEX(multi_1d(slot)%marker_color)
            call GKS$SET_PMARK_SIZE(multi_1d(slot)%marker_scale*
     1                              hw_tuned_pmark_size)
            call set_pmark_type(multi_1d(slot)%marker_type)
c           Plot markers
            if( multi_1d(slot)%contour ) then
               call GKS$POLYMARKER(multi_1d(slot)%np,
     1                             temp_coord, temp_vec)
            else
               call GKS$POLYMARKER(npts, temp_coord(k), temp_vec(k))
            endif
         endif
c        Plot error bars
         if( multi_1d(slot)%errorbars .and. npts .gt. 0) then
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c           Set error bar to color of marker
            if( multi_1d(slot)%color .eq. 0 ) then
               call GKS$SET_PLINE_COLOR_INDEX(
     1                  multi_1d(slot)%marker_color)
            else
               call set_line_index(multi_1d(slot)%color)
            endif
            do j=1,multi_1d(slot)%np
               if( v_log(iytag) ) then
c                 Need to "undo" log scaling.
c                 If we had the yabs flag passed in, we could just
c                 use values in multi_1d(slot)%data(j) and do the abs/log
c                 here... sigh.
                  dely1 = 10.0**data_vec(j,slot)
                  dely2 = dely1
c                 Back to log, making sure no negatives get in.
                  dely1 = MAX(PM_MIN_FLOAT,
     1                        dely1 - multi_1d(slot)%zval(j))
                  dely1 = ALOG10(dely1)
                  dely2 = MAX(PM_MIN_FLOAT,
     1                        dely2 + multi_1d(slot)%zval(j))
                  dely2 = ALOG10(dely2)
               else
                  dely1 = data_vec(j,slot) - multi_1d(slot)%zval(j)
                  dely2 = data_vec(j,slot) + multi_1d(slot)%zval(j)
               endif
c              Now scale to 0..1 window
               dely1 = (dely1*v_scale(iytag) -
     1                        v_win_min(iytag))*vscl(iytag)
               dely2 = (dely2*v_scale(iytag) -
     1                        v_win_min(iytag))*vscl(iytag)
c              The width of the error bar is fixed for now...
               call pm_draw_error_bar(temp_coord(j), dely1, dely2,
     1              0.025)
            end do
         endif
         else if( plot_mode .eq. 2 ) then
c           Plot bar chart
            icount = icount + 1
            scaled_zero = -v_win_min(iytag)*vscl(iytag)
            call pm_draw_barchart(npts, temp_coord(k), temp_vec(k),
     1                            multi_1d(slot)%color, ncurves,
     2                            icount, scaled_zero, ncolors)
         endif
100   continue
c
c     Draw horizontal line at Y=0 (bar chart plot only)
c     This assumes only one Y axis
      if( plot_mode .eq. 2 ) then
         call set_line_index(DEFAULT_LINE)
         temp_coord(1) = 0.0
         temp_coord(2) = 1.0
         temp_vec(1)   = scaled_zero
         temp_vec(2)   = scaled_zero
         call GKS$POLYLINE(2, temp_coord, temp_vec)
      endif
c
c     Draw in boundaries (suprem3 structure file only)
      call draw_1d_boundaries(h_scale(1), hscl(1),
     1                        h_win_min(1), h_log(1),
     2                        v_scale(1), vscl(1),
     3                        v_win_min(1), v_log(1), 2)
c
c     plot x-axis tic marks...
      call set_line_index(DEFAULT_LINE)
c
      do 110 iaxis=1,2
c        Skip entirely if axis is to be omitted
         if( comp_h_omit(iaxis) ) goto 110
c
         if( h_axis_used(iaxis) ) then
            ix = iaxis
         else
            if( h_axis_used(1) ) then
               ix = 1
            else
               ix = 2
            endif
         endif
c
c        Plot axis line
         call GKS$SET_WINDOW(1, 0.0,1.0,0.0,1.0)
         call GKS$SET_CLIPPING(GKS$K_NOCLIP)
         call GKS$SELECT_XFORM(1)
         call set_line_index(comp_h_color(ix))
         if( l_curv ) then
            call GKS$SET_PLINE_LINEWIDTH(comp_h_thickness(ix))
         else
            call GKS$SET_PLINE_LINEWIDTH(comp_h_thickness(ix)*
     1                                   hw_tuned_width)
         endif
         call move_abs(ZEROO, axis_pos(iaxis))
         call draw_abs(ONEE,  axis_pos(iaxis))
         call GKS$SET_CLIPPING(GKS$K_CLIP)
         call GKS$SELECT_XFORM(1)
c
c        Plot grid lines
         call pm_comp_plot_grid(h_win_min(ix), h_win_max(ix), 'X',
     1          h_tic_major(ix), h_minor_tic_freq(ix), h_log(ix),
     2          comp_h_major_grid(ix), comp_h_minor_grid(ix),
     3          comp_h_major_grid_linetype(ix),
     4          comp_h_minor_grid_linetype(ix))
c        Reset linetype to solid for tic marks
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c
c        Skip tic marks if barchart
         if( plot_mode .eq. 2 ) goto 110
c
c        Define tic mark start/end coordinates for each axis
c        1 == lower, 2 == upper
         tic_start(1)     = 0.0
         tic_start(2)     = 1.0
         short_tic_end(1) = comp_h_minortic_size(ix)
         short_tic_end(2) = 1.0 - comp_h_minortic_size(ix)
         long_tic_end(1)  = comp_h_majortic_factor(ix)*
     1                      comp_h_minortic_size(ix)
         long_tic_end(2)  = 1.0 - comp_h_majortic_factor(ix)*
     1                      comp_h_minortic_size(ix)
         short_log_tic_end(1) = comp_h_lminortic_size(ix)
         short_log_tic_end(2) = 1.0 - comp_h_lminortic_size(ix)
         long_log_tic_end(1)  = comp_h_majortic_factor(ix)*
     1                          comp_h_lminortic_size(ix)
         long_log_tic_end(2)  = 1.0 - comp_h_majortic_factor(ix)*
     1                          comp_h_lminortic_size(ix)
c
         if ( .not. h_log(ix) ) then
c           Find first 'nice' major tic mark below h_win_min
            call pm_get_axis_start(h_win_min(ix), h_win_max(ix),
     1                             h_tic_major(ix), h_cur)
c           Normalize to 0.0 to 1.0 scale
            h_cur = (h_cur - h_win_min(ix))/
     1              (h_win_max(ix) - h_win_min(ix))
            delta = h_tic_major(ix)/
     1      (FLOAT(h_minor_tic_freq(ix))*(h_win_max(ix)-h_win_min(ix)))
            i = 0
c           We use 1.0001 instead of 1.0 as a safeguard agains fp rounding
            do while(h_cur .le. 1.0001)
c              Plot tic marks only if within axis bounds
c              We use -0.0001 instead of 0.0 as a safeguard agains fp rounding
               if(h_cur .ge. -0.0001) then
                  call move_abs(h_cur,tic_start(iaxis))
c                 Determine where a big tic mark should go
                  if( MOD(i,h_minor_tic_freq(ix)) .eq. 0 ) then
                     call draw_abs(h_cur,long_tic_end(iaxis))
                  else
                     call draw_abs(h_cur,short_tic_end(iaxis))
                  endif
               endif
               h_cur = h_cur + delta
               i = i + 1
            enddo
         else
c           log plot; tics every decade
            call GKS$SET_WINDOW(1, h_win_min(ix),h_win_max(ix),0.0,1.0)
            call GKS$SELECT_XFORM(1)
            h_cur = FLOAT(INT(h_win_min(ix))) - 1.0
            h_max = FLOAT(INT(h_win_max(ix))) + 1.0
            limit = INT((h_max - h_cur)/fuzz)
            do i=1,limit
               if( h_cur .ge. h_win_min(ix) .or.
     1             h_cur .le. h_win_max(ix)) then
                  call move_abs(h_cur, tic_start(iaxis))
                  call draw_abs(h_cur, long_log_tic_end(iaxis))
               endif
c              minor tic marks...
               do j=2,9
                  h_minor = h_cur + ALOG10(FLOAT(j))
                  if( h_minor .ge. h_win_min(ix) .or.
     1                h_minor .le. h_win_max(ix)) then
                     call move_abs(h_minor,tic_start(iaxis))
                     call draw_abs(h_minor,short_log_tic_end(iaxis))
                  endif
               end do
               h_cur = h_cur + 1.0
            end do
            if( h_cur .le. h_win_max(ix)) then
               call move_abs(h_cur, tic_start(iaxis))
               call draw_abs(h_cur, long_log_tic_end(iaxis))
            endif
         endif
110   continue
c
c     plot y-axis tic marks...
      do 120 iaxis=1,2
c
c        Skip entirely if axis is to be omitted
         if( comp_v_omit(iaxis) ) goto 120
c
         if( v_axis_used(iaxis) ) then
            iy = iaxis
         else
            if( v_axis_used(1) ) then
               iy = 1
            else
               iy = 2
            endif
         endif
c
c        Plot axis line
         call GKS$SET_WINDOW(1, 0.0,1.0,0.0,1.0)
         call GKS$SET_CLIPPING(GKS$K_NOCLIP)
         call GKS$SELECT_XFORM(1)
         call set_line_index(comp_v_color(iy))
         if( l_curv ) then
            call GKS$SET_PLINE_LINEWIDTH(comp_v_thickness(iy))
         else
            call GKS$SET_PLINE_LINEWIDTH(comp_v_thickness(iy)*
     1                                   hw_tuned_width)
         endif
         call move_abs(axis_pos(iaxis), ZEROO)
         call draw_abs(axis_pos(iaxis), ONEE)
         call GKS$SET_CLIPPING(GKS$K_CLIP)
         call GKS$SELECT_XFORM(1)
c
         call pm_comp_plot_grid(v_win_min(iy), v_win_max(iy), 'Y',
     1          v_tic_major(iy), v_minor_tic_freq(iy), v_log(iy),
     2          comp_v_major_grid(iy), comp_v_minor_grid(iy),
     3          comp_v_major_grid_linetype(iy),
     4          comp_v_minor_grid_linetype(iy))
c        Reset linetype to solid for tic marks
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c
c        Define tic mark start/end coordinates for each axis
c        1 == left, 2 == right
         tic_start(1)     = 0.0
         tic_start(2)     = 1.0
         short_tic_end(1) = comp_v_minortic_size(iy)
         short_tic_end(2) = 1.0 - comp_v_minortic_size(iy)
         long_tic_end(1)  = comp_v_majortic_factor(iy)*
     1                      comp_v_minortic_size(iy)
         long_tic_end(2)  = 1.0 - comp_v_majortic_factor(iy)*
     1                      comp_v_minortic_size(iy)
         short_log_tic_end(1) = comp_v_lminortic_size(iy)
         short_log_tic_end(2) = 1.0 - comp_v_lminortic_size(iy)
         long_log_tic_end(1)  = comp_v_majortic_factor(iy)*
     1                          comp_v_lminortic_size(iy)
         long_log_tic_end(2)  = 1.0 - comp_v_majortic_factor(iy)*
     1                          comp_v_lminortic_size(iy)
c
         if ( .not. v_log(iy) ) then
c           Find first 'nice' major tic mark below v_win_min
            call pm_get_axis_start(v_win_min(iy), v_win_max(iy),
     1                             v_tic_major(iy), v_cur)
c           Normalize to 0.0 to 1.0 scale
            v_cur = (v_cur - v_win_min(iy))/
     1              (v_win_max(iy) - v_win_min(iy))
            delta = v_tic_major(iy)/
     1      (FLOAT(v_minor_tic_freq(iy))*(v_win_max(iy)-v_win_min(iy)))
            i = 0
c           We use 1.0001 instead of 1.0 as a safeguard agains fp rounding
            do while(v_cur .le. 1.0001)
c              Plot tic marks only if within axis bounds
c              We use -0.0001 instead of 0.0 as a safeguard agains fp rounding
               if( v_cur .ge. -0.0001 ) then
                  call move_abs(tic_start(iaxis),v_cur)
c                 Determine where a big tic mark should go
                  if( MOD(i,v_minor_tic_freq(iy)) .eq. 0 ) then
                     call draw_abs(long_tic_end(iaxis),v_cur)
                  else
                     call draw_abs(short_tic_end(iaxis),v_cur)
                  endif
               endif
               v_cur = v_cur + delta
               i = i + 1
            enddo
         else
c           log plot; tics every decade
            call GKS$SET_WINDOW(1, 0.0,1.0,v_win_min(iy),v_win_max(iy))
            call GKS$SELECT_XFORM(1)
c           left side of plot...
            v_cur = FLOAT(INT(v_win_min(iy))) - 1.0
            v_max = FLOAT(INT(v_win_max(iy))) + 1.0
            limit = INT((v_max - v_cur)/fuzz)
            do 620 i=1,limit
               if( v_cur .ge. v_win_min(iy) .or.
     1             v_cur .le. v_win_max(iy)) then
                  call move_abs(tic_start(iaxis), v_cur)
                  call draw_abs(long_log_tic_end(iaxis), v_cur)
               endif
c              minor tic marks...
               do 621 j=2,9
                  v_minor = v_cur + ALOG10(FLOAT(j))
                  if( v_minor .ge. v_win_min(iy) .or.
     1                v_minor .le. v_win_max(iy)) then
                     call move_abs(tic_start(iaxis),v_minor)
                     call draw_abs(short_log_tic_end(iaxis),v_minor)
                  endif
621            continue
               v_cur = v_cur + 1.0
620         continue
            if( v_cur .le. v_win_max(iy)) then
               call move_abs(tic_start(iaxis), v_cur)
               call draw_abs(long_log_tic_end(iaxis), v_cur)
            endif
         endif
120   continue
c
c     horizontal label...
c
c     Magic constants are amount window must be expanded from (0,1)
c     to keep in proportion with previous plotting window,
c     but with full viewport in x. They scale with the ratio EXTENT/100 =
c     for the plot window size used. EXTENT depends if we are using the
c     POSTMINI or CURV window extents
c        The left x window:  (0.0 - xndc0)/EXTENT
c        The right x window: (1.0 - (xndc0 + EXTENT))/EXTENT
      h1 = -xndc0/xextent
      h2 = 1.0 + (1.0 -(xndc0+xextent))/xextent
      call GKS$SET_WINDOW(1, h1, h2, 0.0, 0.2)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      tic_label_pos(1) = 0.12
      tic_label_pos(2) = 0.08
      log_tic_label_pos(1) = 0.11
      log_tic_label_pos(2) = 0.09
      axis_label_pos(1) = 0.01
      axis_label_pos(2) = 0.09
c
      do 130 ix=1,2
c
         if( h_axis_used(ix) .and. .NOT. comp_h_omit(ix) ) then
c
            if( ix .eq. 1 ) then
c              Bottom axis
               call GKS$SET_VIEWPORT(1, 0.0, ratio_x,
     1                              (yndc0-0.1)*ratio_y, yndc0*ratio_y)
            else
c              Top axis
               call GKS$SET_VIEWPORT(1, 0.0, ratio_x,
     1                               yndc1*ratio_y, 
     2                              (yndc1+0.10)*ratio_y)
            endif
            call GKS$SELECT_XFORM(1)
c
c           Skip over tic mark labeling, if asked to
            if(comp_h_nonumbering(ix)) goto 131
c
            if( l_curv ) then
               call tune_charsize(0.035, 0.06, -0.5, 0.03, 0.0458, 0.0)
            else
               call tune_charsize(0.035, 0.06, -0.5, 0.025, 0.0458, 0.0)
            endif
            if( .NOT. x_char_data(ix) ) then
            if ( h_log(ix) ) then
               call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                                 GKS$K_TEXT_VALIGN_HALF)
c              Log axis. Label each major tic mark, adjust left and
c              right label to align with plot frame.
               h_cur = FLOAT(INT(h_win_min(ix))) - 1
               h_minor = h_cur
c              Normalize to 0.0 to 1.0 scale
               h_cur = (h_cur - h_win_min(ix))/
     1                 (h_win_max(ix) - h_win_min(ix))
               delta = 1.0/(h_win_max(ix) - h_win_min(ix))
               do while( h_cur .le. 1.0001 )
                  if( h_cur .ge. -0.0001) then
c                    Format into mathmode string 10^{tic_value}
                     call format_mathmode_exponent(h_minor,num_label1)
                     call mathmode(ws_id, h_cur, log_tic_label_pos(ix),
     1                             num_label1, .FALSE., lerr)
                  endif
                  h_minor = h_minor + 1.0
                  h_cur = h_cur + delta
               end do
            else
               call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                                 GKS$K_TEXT_VALIGN_HALF)
c              Write out labels at major tic intervals
c              Find first 'nice' major tic mark below h_win_min
               call pm_get_axis_start(h_win_min(ix), h_win_max(ix),
     1                                h_tic_major(ix), h_cur)
c              Normalize to 0.0 to 1.0 scale
               h_cur = (h_cur - h_win_min(ix))/
     1                 (h_win_max(ix) - h_win_min(ix))
               delta = h_tic_major(ix)/
     1                 (FLOAT(h_minor_tic_freq(ix))*
     1                 (h_win_max(ix) - h_win_min(ix)))
               i = 0
c              We use 1.0001 instead of 1.0 as a safeguard against fp rounding
               do while(h_cur .le. 1.0001)
c                 Plot tic marks only if within axis bounds
c                 We use -0.0001 instead of 0.0 as a safeguard against
c                 fp rounding
                  if(h_cur .ge. -0.0001) then
c                    Determine when to plot a tic mark label
                     if( MOD(i,h_minor_tic_freq(ix)) .eq. 0 ) then
                        ticvalue=(h_cur*(h_win_max(ix) - h_win_min(ix))+
     1                            h_win_min(ix))
                        if(string_len(comp_h_format(ix)) .eq. 0) then
                           call format_axis_label(ticvalue,
     1                          h_tic_major(ix),num_label,exp_flag)
                        else
                           call cformat_axis_label(ticvalue,
     1                          comp_h_format(ix),num_label)
                        endif
                        ilen = string_len(num_label)
                        call GKS$TEXT(h_cur, tic_label_pos(ix),
     1                                num_label(1:ilen))
                     endif
                  endif
                  h_cur = h_cur + delta
                  i = i + 1
               enddo
            endif
            else
c              Bar chart with character data for X axis
               lfound = .FALSE.
c              Plot the first X axis character dataset found
               do i=1,n_cuts
                  slot = slot_index(i)
                  ixtag = multi_1d(slot)%xaxis_tag
                  if( multi_1d(slot)%selected .eq. 1 .and.
     1               ixtag .eq. ix .and. .NOT. lfound ) then
                     maxlen = 0
                     do j=1,multi_1d(slot)%np
                        ilen = string_len(multi_1d(slot)%char_data(j))
                        maxlen = MAX(ilen,maxlen)
                     end do
c                    Empirical formula to decide when to switch to vertical
c                    orientation of labels...
                     if( 80.0/(FLOAT(multi_1d(slot)%np)*FLOAT(maxlen))
     1                   .gt. 1.0 ) then
                       call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                                         GKS$K_TEXT_VALIGN_HALF)
                       call GKS$SET_TEXT_UPVEC(0.0,1.0)
                       ypos = tic_label_pos(ix)
                     else
                        if( ix .eq. 1 ) then
                        call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_RIGHT,
     1                                          GKS$K_TEXT_VALIGN_HALF)
                        ypos = tic_label_pos(ix) + .05
                        else
                        call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_LEFT,
     1                                          GKS$K_TEXT_VALIGN_HALF)
                        ypos = tic_label_pos(ix) - .05
                        endif
                        call GKS$SET_TEXT_UPVEC(-1.0,0.0)
                        call tune_charsize(0.035, 0.06, -0.5,
     1                                     0.025, 0.03, 0.0)
                     endif
                     do j=1,multi_1d(slot)%np
                        temp_coord(j) = (multi_1d(slot)%coord(j)
     1                  *pm_scale_cm2plot
     1                  *h_scale(ix) - h_win_min(ix))*hscl(ix)
                     end do
                     do k=1,multi_1d(slot)%np
                        if( temp_coord(k) .ge. -0.0001 .and.
     1                      temp_coord(k) .le.  1.0001 ) then
                        ilen = string_len(multi_1d(slot)%char_data(k))
                        if( ilen .gt. 0 ) then
                        call mathmode(ws_id, temp_coord(k), ypos,
     1                             multi_1d(slot)%char_data(k)(1:ilen),
     2                             .FALSE., lerr)
                        endif
                        endif
                     end do
                     lfound = .TRUE.
                  endif
               end do
               call GKS$SET_TEXT_UPVEC(0.0,1.0)
            endif
c
131         continue
c
c           Write horizontal label
            if( l_curv ) then
               call tune_charsize(0.035,0.06,-0.5, 0.0386,0.0629,0.0)
            else
               call tune_charsize(0.035, 0.06 , -0.5, 0.03, 0.05, 0.0)
            endif
            call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                              GKS$K_TEXT_VALIGN_BOTTOM)
            ilen = string_len(horiz_label(ix))
            ilen3 = string_len(hunits(ix))
            if( h_exponent(ix) .eq. 0 ) then
               if( ilen3 .eq. 0 ) then
                  work_string = horiz_label(ix)(1:ilen)
               else
                  work_string = horiz_label(ix)(1:ilen) //
     1               ' ' // left_unit_delim //
     2               hunits(ix)(1:ilen3) // right_unit_delim
                  ilen = ilen + ilen3 + 3
               endif
            else
               call format_mathmode_exponent(FLOAT(h_exponent(ix)),
     1                                       num_label1)
               ilen2 = string_len(num_label1)
               if( ilen3 .eq. 0 ) then
                  work_string =  horiz_label(ix)(1:ilen) //
     1                        ' ' // left_unit_delim //
     2                        num_label1(1:ilen2) // right_unit_delim
                  ilen = ilen + ilen2 + 3
               else
                  work_string =  horiz_label(ix)(1:ilen) //
     1                        ' ' // left_unit_delim //
     2                        num_label1(1:ilen2) // ' ' //
     3                        hunits(ix)(1:ilen3) // right_unit_delim
                  ilen = ilen + ilen2 + ilen3 + 4
               endif
            endif
            if( ilen .gt. 0 ) then
            call mathmode(ws_id, 0.5, axis_label_pos(ix),
     1                    work_string(1:ilen), .FALSE. , lerr)
            endif
         endif
130   continue
c
c     vertical tic mark labels...
c
c     Magic constants are amount window must be expanded from (0,1)
c     to keep in proportion with previous plotting window,
c     but with full viewport in y. They scale with the ratio EXTENT/100
c     for the current plot window used.
c        The lower y window:  (0.0 - yndc0)/EXTENT
c        The upper y window:  (1.0 - (yndc0 + EXTENT))/EXTENT
      h1 = -yndc0/yextent
      h2 = 1.0 + (1.0 - (yndc0 + yextent))/yextent
      call GKS$SET_WINDOW  (1, 0.0, 0.15, h1, h2)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      tic_label_pos(1) = 0.13
      tic_label_pos(2) = 0.02
      log_tic_label_pos(1) = 0.14
      log_tic_label_pos(2) = 0.01
      axis_label_pos(1)  = 0.035
      axis_label_pos(2)  = 0.115
C      axis_label_pos2(1) = 0.050
      axis_label_pos2(1) = 0.060  ! Mod value from Rafael's code
      axis_label_pos2(2) = 0.100
      do 140 iy=1,2
         exp_flag = .FALSE.
         if( v_axis_used(iy) .and. .NOT. comp_v_omit(iy) ) then
            if( iy .eq. 1 ) then
               call GKS$SET_VIEWPORT(1,
     1              (xndc0-xoffset)*ratio_x, xndc0*ratio_x,
     2              0.0*ratio_y, ratio_y)
               call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_RIGHT,
     1                                 GKS$K_TEXT_VALIGN_HALF)
            else
               call GKS$SET_VIEWPORT(1,
     1              xndc1*ratio_x, (xndc1+xoffset)*ratio_x,
     2              0.0*ratio_y, ratio_y)
               call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_LEFT,
     1                                 GKS$K_TEXT_VALIGN_HALF)
            endif
            call GKS$SELECT_XFORM(1)
c
c           Skip over tic mark labeling, if asked to
            if( comp_v_nonumbering(iy) ) goto 141
c
            if ( v_log(iy) ) then
c              Log axis. Label each major tic mark, adjust bottom and
c              top label to align with plot frame.
               if( l_curv ) then
                  call tune_charsize(0.015,0.0615,-0.2,0.0135,0.036,0.0)
               else
                  call tune_charsize(0.015,0.0615,-0.2,0.015,0.0346,0.0)
               endif
               v_cur = FLOAT(INT(v_win_min(iy))) - 1
               v_minor = v_cur
c              Normalize to 0.0 to 1.0 scale
               v_cur = (v_cur - v_win_min(iy))/
     1                 (v_win_max(iy) - v_win_min(iy))
               delta = 1.0/(v_win_max(iy) - v_win_min(iy))
               do while( v_cur .le. 1.0001 )
                  if( v_cur .ge. -0.0001) then
c                    Format into mathmode string 10^{tic_value}
                     call format_mathmode_exponent(v_minor,num_label1)
                     call mathmode(ws_id, log_tic_label_pos(iy), v_cur,
     1                             num_label1, .FALSE., lerr)
                  endif
                  v_minor = v_minor + 1.0
                  v_cur = v_cur + delta
               end do
            else
c              Linear axis
               if( l_curv ) then
                  call tune_charsize(0.015,0.0615,-0.2,0.0135,0.036,0.0)
               else
                  call tune_charsize(0.015,0.0615,-0.2,0.015,0.0346,0.0)
               endif
c              Write labels on major tic marks...
c              Find first 'nice' major tic mark below v_win_min
               call pm_get_axis_start(v_win_min(iy), v_win_max(iy),
     1                                v_tic_major(iy), v_cur)
c              Normalize to 0.0 to 1.0 scale
               v_cur = (v_cur - v_win_min(iy))/
     1                 (v_win_max(iy) - v_win_min(iy))
               delta = v_tic_major(iy)/
     1                 (FLOAT(v_minor_tic_freq(iy))*
     2                 (v_win_max(iy) - v_win_min(iy)))
               i = 0
c              We use 1.0001 instead of 1.0 as a safeguard against fp rounding
               do while(v_cur .le. 1.0001)
c                 Plot tic marks only if within axis bounds
c                 We use -0.0001 instead of 0.0 as a safeguard against
c                 fp rounding
                  if( v_cur .ge. -0.0001 ) then
c                    Determine where a big tic mark should go
                     if( MOD(i,v_minor_tic_freq(iy)) .eq. 0 ) then
                        ticvalue=(v_cur*(v_win_max(iy) - v_win_min(iy))
     1                             + v_win_min(iy))
                        if(string_len(comp_v_format(iy)) .eq. 0) then
                           call format_axis_label(ticvalue,
     1                          v_tic_major(iy),num_label2, exp_flag)
                        else
                           call cformat_axis_label(ticvalue,
     1                          comp_v_format(iy),num_label2)
                        endif
                        ilen = string_len(num_label2)
                        call GKS$TEXT(tic_label_pos(iy), v_cur,
     1                                num_label2(1:ilen))
                     endif
                  endif
                  v_cur = v_cur + delta
                  i = i + 1
               enddo
            endif
c
141         continue
c           vertical label
            call GKS$SET_TEXT_UPVEC(-1.0,0.0)
            if( l_curv ) then
               call tune_charsize(0.035, 0.015, 0.0, 0.050, 0.018, 0.0)
            else
c               call tune_charsize(0.037, 0.02, 0.0, 0.046, 0.02, 0.0)
               call tune_charsize(0.037, 0.02, 0.0, 0.0391, 0.017, 0.0)
            endif
            call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                              GKS$K_TEXT_VALIGN_HALF)
            ilen = string_len(vert_label(iy))
            ilen3 = string_len(vunits(iy))
            if( v_exponent(iy) .eq. 0 ) then
               if( ilen3 .eq. 0 ) then
                  work_string = vert_label(iy)(1:ilen)
               else
                  work_string = vert_label(iy)(1:ilen) //
     1               ' ' // left_unit_delim //
     2               vunits(iy)(1:ilen3) // right_unit_delim
                  ilen = ilen + ilen3 + 3
               endif
            else
               call format_mathmode_exponent(FLOAT(v_exponent(iy)),
     1                                       num_label1)
               ilen2 = string_len(num_label1)
               if( ilen3 .eq. 0 ) then
                  work_string =  vert_label(iy)(1:ilen) //
     1                        ' ' // left_unit_delim //
     2                        num_label1(1:ilen2) // right_unit_delim
                  ilen = ilen + ilen2 + 3
               else
                  work_string =  vert_label(iy)(1:ilen) //
     1                        ' ' // left_unit_delim //
     2                        num_label1(1:ilen2) // ' ' //
     3                        vunits(iy)(1:ilen3) // right_unit_delim
                  ilen = ilen + ilen2 + ilen3 + 4
               endif
            endif
            if( exp_flag ) then
               if( ilen .gt. 0 ) then
               call mathmode(ws_id, axis_label_pos(iy), 0.5,
     1                       work_string(1:ilen), .FALSE., lerr)
               endif
            else
               if( ilen .gt. 0 ) then
               call mathmode(ws_id, axis_label_pos2(iy), 0.5,
     1                       work_string(1:ilen), .FALSE., lerr)
               endif
            endif
         endif
         call GKS$SET_TEXT_UPVEC(0.0,1.0)
140   continue
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                        GKS$K_TEXT_VALIGN_NORMAL)
c
      call GKS$SET_PLINE_LINEWIDTH(hw_tuned_width)
c
c     titling...
c
      call compare_draw_titles(ws_id,ratio_x,ratio_y,l_curv,
     1                         h_axis_used(2), comp_nolabels)
c
c     plot labels...
      if( .NOT. comp_nolabels ) then
         call GKS$SET_WINDOW(1, 0.0, 1.0, 0.0, 0.2)
         call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 1.0*ratio_x,
     1                            0.0*ratio_y, 0.15*ratio_y)
         call GKS$SELECT_XFORM(1)
         h_title = 0.0
         line_length = 0.04
         delta = 0.05
         v_title = 0.15
         nactive = 0
         do i=1,n_cuts
            slot = slot_index(i)
            if( multi_1d(slot)%selected .ne. 0 .and.
     1          string_len(multi_1d(slot)%label) .gt. 0)
     2          nactive = nactive + 1
         end do
c        Gag: yet another set of magic constants :-) jvf
c        Set dependent on how many labels are to be plotted
         if( nactive .le. 12 ) then
            call tune_charsize(0.025, 0.04, -0.6, 0.013, 0.02, 0.0)
         else
            call tune_charsize(0.02, 0.032, -0.6, 0.0104, 0.016, 0.0)
         endif
c        Compute how much space each column of data labels is allowed
         h_delta = 1.0/(1.0 + INT((nactive-1)/4))

         call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                           GKS$K_TEXT_VALIGN_NORMAL)
         j = 0
         do 300 i=1,n_cuts
            slot = slot_index(i)
            if( multi_1d(slot)%selected .eq. 0 .or.
     1          string_len(multi_1d(slot)%label) .eq. 0 ) goto 300
c           Plot a line/marker to identify string with curve
            if( multi_1d(slot)%color .ne. 0 .and.
     1          multi_1d(slot)%linestyle .ne. 0 ) then
c              First, reset linetype to solid
               call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c              Use call to set_line_index so that different linestyles
c              will be used on monochrome output devices to "simulate"
c              different colors.
               call set_line_index(multi_1d(slot)%color)
               if( multi_1d(slot)%linestyle .ne. 1 ) then
                  call set_pline_type(multi_1d(slot)%linestyle)
               endif
               if( multi_1d(slot)%linewidth .ne. 1.0 ) then
                  call GKS$SET_PLINE_LINEWIDTH(
     1                 multi_1d(slot)%linewidth*hw_tuned_width)
               endif
               call move_abs(h_title, v_title)
               call draw_abs(h_title+line_length, v_title)
               if( multi_1d(slot)%linewidth .ne. 1.0 ) then
                  call GKS$SET_PLINE_LINEWIDTH(hw_tuned_width)
               endif
            endif
            if( multi_1d(slot)%marker_color .ne. 0 .and.
     1          multi_1d(slot)%marker_type .ne. 0 ) then
               call GKS$SET_PMARK_COLOR_INDEX(
     1                 multi_1d(slot)%marker_color)
               call set_pmark_type(multi_1d(slot)%marker_type)
c              Plot marker
               call GKS$POLYMARKER(1, h_title+0.5*line_length,
     1                             v_title)
            endif
            ilen = string_len(multi_1d(slot)%label)
c           Hack for VTxxx terminals. VTxxx only write out up to 35
c           characters for (legible) label, but stroke workstations
c           use a smaller character size, and can write longer labels.
            if( use_stroke_prec ) then
               if( ilen .gt. 0 ) then
               call mathmode(ws_id,h_title+delta, v_title,
     1                    multi_1d(slot)%label(1:ilen),.FALSE.,lerr)
               endif
            else
               ilen = MIN(35,ilen)
               if( ilen .gt. 0 ) then
               call mathmode(ws_id,h_title+delta, v_title,
     1                    multi_1d(slot)%label(1:ilen),.FALSE.,lerr)
               endif
            endif
            v_title = v_title - 0.04
            j = j + 1
c           switch to next column after every 4 plots
            if( MOD(j,4) .eq. 0 ) then
               h_title = h_title + h_delta
               v_title = 0.15
            endif
300      continue
      endif
c
      call set_line_index(DEFAULT_LINE)
c
c     Draw side labels (callable version only)
      call callable_draw_sidelabels
c
      call deactivate_plot_device
c
      if( plot_device .eq. hardcopy_device ) then
c        Add in annotations to hardcopy plot
         if( ws_window_open ) then
            call anno_replay_annotations(current_window_id, window_id,
     1                 .TRUE.)
         else
            if( add_notes_from_save ) then
               call anno_restore_notes_from_save(window_id, lerr)
            endif
         endif
         call close_plot_device(window_id)
      else
         if( add_notes_from_save ) then
            call anno_restore_notes_from_save(window_id, lerr)
         else if ( add_notes_from_file ) then
c           Restore any previous annotations
            call read_savefile('POSTMINI.NOTE',.TRUE.,lerr)
            if( .NOT. lerr )
     1         call anno_restore_notes_from_save(window_id, lerr)
         endif
         if( is_vtxxx(window_id) ) then
            call move_cursor_bottom_of_screen
            call beep_terminal
            call get_char(response,'Hit return to continue ',0)
c           Restore default color map
            call reset_vt_colormap
         endif
      endif
c
      ierror = PLT$NORMAL
c
      return
      end
c
      subroutine recompute_minmax_comp(dmin, dmax,
     1           dmin_scaled, dmax_scaled, majtic, ticfreq, iuser_scale,
     2           iaxis_scale, do_abs, do_log, do_uscale, iaxis, work,
     3           plot_mode)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'compare.cmn'
c
      integer ticfreq, iuser_scale, iaxis_scale, iaxis, plot_mode
      real work(COMP_ARRAYSIZE)
      real dmin, dmax, dmin_scaled, dmax_scaled, majtic
      real ymin, ymax, yplus, yminus
      logical do_abs, do_log, do_uscale
c
c     locals
      integer i, j, slot
      real user_scale, axis_scale, log_small
      logical no_data, do_errorbars
      real ONE
      parameter (ONE = 1)
c
      no_data = .TRUE.
      log_small = ALOG10(PM_MIN_FLOAT)
      if( do_log ) then
         dmin =  38
         dmax = -38
      else if( do_abs ) then
         dmin = 1.0E38
         dmax = 0.0
      else
         dmin =  1.0E38
         dmax = -1.0E38
      endif
      if( do_uscale ) then
         user_scale = 10.0**iuser_scale
      else
         user_scale = 1.0
      endif
c
      do 100 i=1,n_cuts
         slot = slot_index(i)
c        Only look at data which has been selected
         if( multi_1d(slot)%selected .eq. 1 .and.
     1       multi_1d(slot)%yaxis_tag .eq. iaxis ) then
            no_data = .FALSE.
            do_errorbars = multi_1d(slot)%errorbars
            do 110 j=1, multi_1d(slot)%np
               call pm_abs_log_scale_datum(multi_1d(slot)%data(j),
     1              do_abs,do_log,do_uscale,user_scale,ONE,log_small,
     2              work(j))
               if( do_errorbars ) then
                  call pm_abs_log_scale_datum(
     1               multi_1d(slot)%data(j)-multi_1d(slot)%zval(j),
     2               do_abs,do_log,do_uscale,user_scale,ONE,log_small,
     3               yminus)
                  call pm_abs_log_scale_datum(
     1               multi_1d(slot)%data(j)+multi_1d(slot)%zval(j),
     2               do_abs,do_log,do_uscale,user_scale,ONE,log_small,
     3               yplus)
               else
                  yplus  = work(j)
                  yminus = work(j)
               endif
c              Now compute min/max
               ymin = MIN(work(j),yminus)
               ymax = MAX(work(j),yplus)
               if( dmin .gt. ymin .and. ymin .ne. log_small)
     1             dmin = ymin
               if( dmax .lt. ymax ) dmax = ymax
110         continue
         endif
100   continue
c
      if( no_data ) then
c        Return dummy limits
         if( do_log ) then
            dmin = 1
            dmax = 2
         else
            dmin = 0.0
            dmax = 1.0
         endif
      endif
c
      if( plot_mode .eq. 2 ) then
c        Bar chart type plot
         if( dmin .gt. 0.0 ) dmin = 0.0
         if( dmax .lt. 0.0 ) dmax = 0.0
      endif
c
      call get_default_limits(dmin, dmax, dmin_scaled,
     1     dmax_scaled, majtic, ticfreq, iaxis_scale, do_log)
      if( do_log ) then
         axis_scale = 1.0
         dmin = 10.0**dmin
         dmax = 10.0**dmax
         dmin_scaled = 10.0**dmin_scaled
         dmax_scaled = 10.0**dmax_scaled
      else
         axis_scale  = 10.0**iaxis_scale
         if( do_uscale ) then
            majtic      = majtic*axis_scale
CC            dmin = dmin*axis_scale
CC            dmax = dmax*axis_scale
            dmin_scaled = dmin_scaled*axis_scale
            dmax_scaled = dmax_scaled*axis_scale
         endif
         dmin = dmin*user_scale
         dmax = dmax*user_scale
      endif
c
      return
      end
c
      subroutine recompute_minmax_coords(dmin, dmax,
     1           dmin_scaled, dmax_scaled, majtic, ticfreq, iuser_scale,
     2           iaxis_scale, do_abs, do_log, do_uscale, iaxis, work)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'units.cmn'
      include 'compare.cmn'
c
      integer ticfreq, iuser_scale, iaxis_scale, iaxis
      real work(COMP_ARRAYSIZE)
      real dmin, dmax, dmin_scaled, dmax_scaled, majtic
      logical do_abs, do_log, do_uscale
c
c     locals
      integer i, j, slot
      real user_scale, axis_scale, log_small
      logical no_data
c
      log_small = ALOG10(PM_MIN_FLOAT)
      no_data = .TRUE.
      if( do_log ) then
         dmin =  38
         dmax = -38
      else if( do_abs ) then
         dmin = 1.0E38
         dmax = 0.0
      else
         dmin =  1.0E38
         dmax = -1.0E38
      endif
      if( do_uscale ) then
         user_scale = 10.0**iuser_scale
      else
         user_scale = 1.0
      endif
c
      do 100 i=1,n_cuts
         slot = slot_index(i)
c        Only look at data which has been selected
         if( multi_1d(slot)%selected .eq. 1 .and.
     1       multi_1d(slot)%xaxis_tag .eq. iaxis ) then
            no_data = .FALSE.
            do 110 j=1, multi_1d(slot)%np
               call pm_abs_log_scale_datum(multi_1d(slot)%coord(j),
     1              do_abs,do_log,do_uscale,user_scale,pm_scale_cm2plot,
     2              log_small,work(j))
c              Now compute min/max
               if( dmin .gt. work(j) .and. work(j) .ne. log_small)
     1             dmin = work(j)
               if( dmax .lt. work(j) ) dmax = work(j)
110         continue
         endif
100   continue
c
      if( no_data ) then
c        Return dummy limits
         if( do_log ) then
            dmin = 1
            dmax = 2
         else
            dmin = 0.0
            dmax = 1.0
         endif
      endif
c
      call get_default_limits(dmin, dmax, dmin_scaled,
     1     dmax_scaled, majtic, ticfreq, iaxis_scale, do_log)
      if( do_log ) then
         axis_scale = 1.0
         dmin = 10.0**dmin
         dmax = 10.0**dmax
         dmin_scaled = 10.0**dmin_scaled
         dmax_scaled = 10.0**dmax_scaled
      else
         axis_scale  = 10.0**iaxis_scale
         if( do_uscale ) then
            majtic      = majtic*axis_scale
            dmin_scaled = dmin_scaled*axis_scale
            dmax_scaled = dmax_scaled*axis_scale
         endif
         dmin = dmin*user_scale
         dmax = dmax*user_scale
      endif
c
      return
      end
c
      subroutine swept_polyline(n, x, y)
      implicit none
c
      integer n
      real x(n), y(n)
c
      integer i, i0, k
      logical lforward
c
c     Scans x array and looks for places where the x coordinate
c     'retraces' backwards. At these points, output a polyline,
c     and skip to the next point.
c
c     Determine forward 'direction' from first two x coordinates
      if( n .ge. 2 ) then
         if( x(2) .ge. x(1) ) then
            lforward = .TRUE.
         else
            lforward = .FALSE.
         endif
      else
         lforward = .TRUE.
      endif
c
      k  = 1
      i0 = 1
      do 100 i=1,n-1
         if( (lforward .and. x(i+1) .lt. x(i)) .or.
     1       (.NOT. lforward .and. x(i+1) .gt. x(i)) ) then
c           Avoid any GKS errors if we have only one point in a sweep
            if( k .gt. 1 ) then
               call gks$polyline(k, x(i0), y(i0))
            endif
            k  = 1
            i0 = i+1
         else
            k  = k+1
         endif
100   continue
c     Output any remaining polyline
      if( k .ne. 1 ) call gks$polyline(k, x(i0), y(i0))
c
      return
      end
c
      subroutine compare_draw_titles(ws_id,ratio_x,ratio_y,
     1                               lcurv, lupper_axis, lno_labels)
      implicit none
c
c     Function: Draw titles on comparison plots
c
c     Programmer: John Faricelli
c     Date:    11-Mar-93
c
c     Input arguments:
c
c     ws_id     -- GKS workstation id of currently active workstation
c     ratio_x   -- X and Y aspect ratio of currently active workstation
c     ratio_y 
c     lcurv     -- curv style titles (if .TRUE.)
c                  center titles and use curv viewport
c     lupper_axis -- If .TRUE. then an upper X axis has been specified
c     lno_labels -- Labels have been suppressed, use larger viewport
c                   for titles
c
      logical lcurv, lupper_axis, lno_labels
      integer ws_id
      real ratio_x, ratio_y
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'gksdefs.cmn'
c
      include 'title.cmn'
c
      real ylabel, yheight, xpos, ypos, ydel
c
      integer i, ilen, string_len
      logical lerr
c
      if( lcurv ) then
         xpos = 0.5
         ylabel = 0.3
         if( lno_labels ) then
            yheight = 0.6
            call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 1.0*ratio_x,
     1                               0.70*ratio_y,1.0*ratio_y)
         else
            yheight = 0.3
            call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 1.0*ratio_x,
     1                               0.85*ratio_y,1.0*ratio_y)
         endif
         call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                           GKS$K_TEXT_VALIGN_HALF)
      else
         xpos = 0.0
         ylabel = 0.2
         if( lno_labels ) then
            yheight = 0.5
            call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 1.0*ratio_x,
     1                               0.75*ratio_y,1.0*ratio_y)
         else
            yheight = 0.2
            call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 1.0*ratio_x,
     1                               0.90*ratio_y,1.0*ratio_y)
         endif
         call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                           GKS$K_TEXT_VALIGN_HALF)
      endif
c     Try to keep Y world coordinate space in scale with viewport
c     so character size tuning less painful.
      call GKS$SET_WINDOW(1, 0.0,1.0,0.0,yheight)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      call GKS$SELECT_XFORM(1)
c
c     Code needs to be checked for case of VTxxx hardware fonts
c     Things look OK on a DECterm in ReGIS mode
c
      if( num_titles .le. 3 ) then
         call tune_charsize(0.02, 0.07, -0.3, 0.02, 0.05, 0.0)
         ydel = 0.08
      else
         call tune_charsize(0.01,0.0175,-0.02,0.014,0.035,0.0)
         ydel = 0.05
      endif
      if( lupper_axis ) then
         ypos = 0.66666*ylabel
      else
         ypos = ylabel/FLOAT(num_titles+1)
      endif
      do i=num_titles,1,-1
         ilen = string_len(title_dat(i))
         if( ilen .gt. 0 ) then
            call mathmode(ws_id, xpos, ypos, title_dat(i)(1:ilen),
     1                    .FALSE., lerr)
         endif
         ypos = ypos + ydel
      end do
c
      return
      end
c
      subroutine scale_like_zoom(old_low, old_high, new_low, new_high)
      implicit none
c
      real old_low(2), old_high(2), new_low(2), new_high(2)
c
      real h
c
c     Scale second axis limits in proportion to first axis limits
c     New axis limits are always within the old ones
c
c     Original 1st axis limits: old_low(1), old_high(1)
c     New      1st axis limits: new_low(1), new_high(1)
c     Original 2nd axis limits: old_low(2), old_high(2)
c
c     Outputs:
c     New      2nd axis limits: new_low(2), new_high(2)
c
      h = (new_low(1) - old_low(1))/(old_high(1) - old_low(1))
      new_low(2) = old_low(2) + (old_high(2) - old_low(2))*h
c
      h = (new_high(1) - old_low(1))/(old_high(1) - old_low(1))
      new_high(2) = old_low(2) + (old_high(2) - old_low(2))*h
c
      return
      end
c
      subroutine pm_get_compare_plot_ndc_space(xndc0, yndc0,
     1           xndc1, yndc1)
      implicit none
c
      real xndc0, yndc0, xndc1, yndc1
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'defaults.cmn'
      include 'compare.cmn'
c
c     This subroutine sets the NDC (normalized device coordinate)
c     space limits where the curves are plotted, exclusing of axis labels,
c     etc. for the comparison plot.
c
c     We support two spaces: the original POSTMINI drawing space
c     and the CURV drawing space.
c
c     Get plot area origin                  
      if( comp_nolabels ) then
         if( l_curv ) then
            xndc0 = COMP_ORIGIN_CURV_NO_LABELS_X
         else
            xndc0 = COMP_ORIGIN_NO_LABELS_X
         endif
         yndc0 = COMP_ORIGIN_NO_LABELS_Y
      else
         if( l_curv ) then
            xndc0 = COMP_ORIGIN_CURV_WITH_LABELS_X
         else
            xndc0 = COMP_ORIGIN_WITH_LABELS_X
         endif
         yndc0 = COMP_ORIGIN_WITH_LABELS_Y
      endif
c
      if( l_curv ) then
         xndc1 = xndc0 + COMP_CURV_EXTENT_X
         yndc1 = yndc0 + COMP_CURV_EXTENT_Y
      else
         xndc1 = xndc0 + COMP_EXTENT_X
         yndc1 = yndc0 + COMP_EXTENT_Y
      endif
c
      return
      end
c
      subroutine pm_comp_plot_grid(win_min, win_max, axis_type,
     1              tic_major, minor_tic_freq, llog,
     2              lmajor_grid, lminor_grid,
     3              major_grid_linetype, minor_grid_linetype)
      implicit none
c
      real win_min, win_max, tic_major
      integer minor_tic_freq, major_grid_linetype, minor_grid_linetype
      logical llog, lmajor_grid, lminor_grid
      character*1 axis_type
c
      include 'gksdefs.cmn'
c
c     Two parameters close to zero and one, but within the range (0,1)
      real ZEROO, ONEE
      parameter( ZEROO = 0.000001 )
      parameter( ONEE  = 0.999999 )
c
      real delta, cur, max, minor, fuzz, xx(2), yy(2)
      integer i, j, limit
c
      fuzz = 0.999
      xx(1) = ZEROO
      xx(2) = ONEE
      yy(1) = ZEROO
      yy(2) = ONEE
c
      if ( .not. llog ) then
c        ASSUMES that the GKS plot window is 0.0 to 1.0 in each direction
c        Find first 'nice' major tic mark below h_win_min
         call pm_get_axis_start(win_min, win_max, tic_major, cur)
c        Normalize to 0.0 to 1.0 scale
         cur = (cur - win_min)/(win_max - win_min)
         delta = tic_major/
     1   (FLOAT(minor_tic_freq)*(win_max-win_min))
         i = 0
c        We use 1.0001 instead of 1.0 as a safeguard agains fp rounding
         do while(cur .le. 1.0001)
c           Plot tic marks only if within axis bounds
c           We use -0.0001 instead of 0.0 as a safeguard agains fp rounding
            if(cur .ge. -0.0001) then
               if( axis_type .eq. 'X' ) then
                  xx(1) = cur
                  xx(2) = cur
               else
                  yy(1) = cur
                  yy(2) = cur
               endif
c              Determine where a major grid line should go
               if( MOD(i,minor_tic_freq) .eq. 0 ) then
                  if( lmajor_grid) then
                     call gks$set_pline_linetype(major_grid_linetype)
                     call gks$polyline(2,xx,yy)
                  endif
               else
                  if( lminor_grid) then
                     call gks$set_pline_linetype(minor_grid_linetype)
                     call gks$polyline(2,xx,yy)
                  endif
               endif
            endif
            cur = cur + delta
            i = i + 1
         enddo
      else
c
c        log plot, major tics every decade
         if( axis_type .eq. 'X' ) then
            call GKS$SET_WINDOW(1, win_min,win_max,0.0,1.0)
         else
            call GKS$SET_WINDOW(1, 0.0,1.0,win_min,win_max)
         endif
         call GKS$SELECT_XFORM(1)
         cur = FLOAT(INT(win_min)) - 1.0
         max = FLOAT(INT(win_max)) + 1.0
         limit = INT((max - cur)/fuzz)
         do i=1,limit
            if( axis_type .eq. 'X' ) then
               xx(1) = cur
               xx(2) = cur
            else
               yy(1) = cur
               yy(2) = cur
            endif
            if( cur .ge. win_min .or. cur .le. win_max) then
               if( lmajor_grid ) then
                  call gks$set_pline_linetype(major_grid_linetype)
                  call gks$polyline(2,xx,yy)
               endif
            endif
c           minor tic marks...
            do j=2,9
               minor = cur + ALOG10(FLOAT(j))
               if( axis_type .eq. 'X' ) then
                  xx(1) = minor
                  xx(2) = minor
               else
                  yy(1) = minor
                  yy(2) = minor
               endif
               if( minor .ge. win_min .or.
     1             minor .le. win_max) then
                  if( lminor_grid) then
                     call gks$set_pline_linetype(minor_grid_linetype)
                     call gks$polyline(2,xx,yy)
                  endif
               endif
            end do
            cur = cur + 1.0
         end do
         if( cur .le. win_max) then
            if( axis_type .eq. 'X' ) then
               xx(1) = cur
               xx(2) = cur
            else
               yy(1) = cur
               yy(2) = cur
            endif
            if( lmajor_grid) then
               call gks$set_pline_linetype(major_grid_linetype)
               call gks$polyline(2,xx,yy)
            endif
         endif
      endif
c
      end
c
      subroutine pm_sort(npts, x, y)
      implicit none
c
      integer npts
c
      real x(npts), y(npts)
c
      real tmp
      integer i, j
c
c     Simple (read slow) in place sort.
c     I'd rather use quicksort (qsort()), but I'd prefer to keep
c     any dupicate x values in original order. qsort() can't guarantee
c     that.
      do 100 j=npts-1,1,-1
         do 100 i=1,j
            if( x(i+1) .lt. x(i) ) then
               tmp    = x(i+1)
               x(i+1) = x(i)
               x(i)   = tmp
               tmp    = y(i+1)
               y(i+1) = y(i)
               y(i)   = tmp
            endif
100   continue
c
      return
      end
c
      subroutine pm_zpolyline(npts, x, y, z)
      implicit none
c
      integer npts
c
      real x(npts), y(npts), z(npts)
c
      real tmp
      integer i, j, istart, iend, n
c
c     Simple (read slow) in place sort.
c     I'd rather use quicksort (qsort()), but I NEED to keep
c     duplicate z values in original order. qsort() can't guarantee
c     that.
c
c     A double sort (sort on Z, then on X) would probably be OK, though,
c     even if the original order was not maintained.
c
      do 100 j=npts-1,1,-1
         do 100 i=1,j
            if( z(i) .lt. z(i+1) ) then
               tmp    = x(i+1)
               x(i+1) = x(i)
               x(i)   = tmp
               tmp    = y(i+1)
               y(i+1) = y(i)
               y(i)   = tmp
               tmp    = z(i+1)
               z(i+1) = z(i)
               z(i)   = tmp
            endif
100   continue
c
c     OK, now plot polylines for each constant value of Z
c
      istart = 1
      iend   = 1
      do 200 i=2,npts
         if( z(i-1) .ne. z(i) ) then
            n = iend - istart + 1
            if( n .gt. 1 ) call gks$polyline(n,x(istart),y(istart))
            istart = i
            iend   = i
         else
            iend = i
         endif
200   continue
c     Emit final polyline
      n = iend - istart + 1
      if( n .gt. 1 ) call gks$polyline(n,x(istart),y(istart))
c
      return
      end
c
      subroutine pm_get_axis_start(rmin,rmax,tic_major,start)
      implicit none
c
      real rmin, rmax, tic_major, start
c
c     Determine a 'nice' starting point for the major tics on an axis
c
      real delta, diff
c
c     Find first major tic mark below min
c     Special case: if max-min can be exactly divided by the major
c                   tic mark, then start at min.
      delta = rmax - rmin
      diff = ABS(delta - tic_major*NINT(delta/tic_major))
      if( diff .lt. 0.0001 ) then
         start = rmin
      else
         if( rmin .lt. 0.0 ) then
            start = FLOAT(INT(rmin/tic_major) - 1)*tic_major
         else
            start = FLOAT(INT(rmin/tic_major))*tic_major
         endif
      endif
c
      return
      end
c
      subroutine pm_draw_error_bar(x, y1, y2, size)
      implicit none
c
      real x, y1, y2, size
c
      real xx(2), yy(2)
c
      xx(1) = x
      xx(2) = x
      yy(1) = y1
      yy(2) = y2
      call GKS$POLYLINE(2,xx,yy)
      xx(1) = x - 0.5*size
      xx(2) = x + 0.5*size
      yy(1) = y1
      yy(2) = y1
      call GKS$POLYLINE(2,xx,yy)
      yy(1) = y2
      yy(2) = y2
      call GKS$POLYLINE(2,xx,yy)
c
      return
      end
c
      subroutine pm_abs_log_scale_datum(value,
     1              do_abs,do_log,do_uscale,user_scale,scale,
     2              log_small,result)
      implicit none
c
      real value, user_scale, log_small, scale, result
      logical do_abs, do_log, do_uscale
c
      real svalue
c
      svalue = scale*value
c
      if( do_abs .or. do_log ) then
         if( do_uscale ) then
            result = ABS(svalue/user_scale)
         else
            result = ABS(svalue)
         endif
         if( do_log ) then
c           check for log(0.0), fudge
            if( result .eq. 0.0 ) then
               result = log_small
            else
               result = ALOG10(result)
            endif
         endif
      else
         if( do_uscale ) then
            result = svalue/user_scale
         else
            result = svalue
         endif
      endif
c
      return
      end
c
      subroutine pm_draw_barchart(npts, coords, data, color,
     1                            num_datasets, idataset,
     2                            scaled_zero, ncolors)
      implicit none
c
      integer npts, color, num_datasets, idataset, ncolors
      real coords(npts), data(npts), scaled_zero
c
      include 'gksdefs.cmn'
c
      integer i
      real x(5), y(5), width, hwidth, offset
c
c     Code assumes uniform spacing of bar chart
c     Width of bar
      width = (coords(2) - coords(1))/FLOAT(num_datasets+1)
      hwidth = 0.5*width
      offset = - FLOAT(num_datasets - 1)*hwidth +
     1           width*FLOAT(idataset-1)
c
      if( ncolors .eq. 1 ) then
         if( color .eq. 1 ) then
            call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
            call GKS$SET_FILL_COLOR_INDEX(color)
         else
            call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_HATCH)
c           GKS HATCH values are negative (device specific) -1 to -9
            call GKS$SET_FILL_STYLE_INDEX(MAX(MIN(-1,-(color-1)),-9))
         endif
      else
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
         call GKS$SET_FILL_COLOR_INDEX(color)
      endif
c
      do i=1,npts
         x(1) = coords(i) - hwidth + offset
         x(2) = coords(i) + hwidth + offset
         x(3) = x(2)
         x(4) = x(1)
         y(1) = scaled_zero
         y(2) = scaled_zero
         y(3) = data(i)
         y(4) = data(i)
c        Workaround for GKS 6.5 bug: a box that lies outside the clipping
c        boundary but with one edge on it causes GKS to fault.
c        Note: the workaround assumes we have scaled all X data to 0..1
         if( x(1) .lt. 1.0 .and. x(2) .gt. 0.0) then
            call GKS$FILL_AREA(4,x,y)
         endif
         if( ncolors .eq. 1 ) then
c           Draw box around bar
            call GKS$SET_PLINE_COLOR_INDEX(1)
            call GKS$POLYLINE(5,x,y)
         endif
      end do
c
      return
      end
