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     2D comparison contour plot utilities
c
c
      subroutine compare_plot_2d(lsave,lplotfirst)
c
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
      include 'compare_2d.cmn'
      include 'units.cmn'
      include 'defaults.cmn'
      include 'title.cmn'
      include 'callable.cmn'
c
      logical lsave, lplotfirst
c
c     locals
c
      logical contour_values_set, fake_ln03,
     2        do_draw_to_scale, do_uniform_contours, do_log, do_negate,
     3        do_yflip, lentire_matrix
      character*80 label, units
      character*50 horiz_label, vert_label
      character*20 horiz_units, vert_units
      integer i, int_vec(10), ilen, ulen, string_len
      integer idum1
      real real_vec(MAX_CONTOURS)
      real contour_values(MAX_CONTOURS)
      real min_value, max_value
      real delta
      real rdum1, rdum2, rdum3, rdum4
      save horiz_label, vert_label, horiz_units, vert_units
CC
CC BEGIN MENU
      real hmin_full, hmax_full, vmin_full, vmax_full
      real hmin_plot, hmax_plot, vmin_plot, vmax_plot
      real plot_min_value, plot_max_value
      real h_tic_major, v_tic_major
      real aspect_ratio
      real overall_hmin, overall_hmax, overall_vmin, overall_vmax,
     1     hmin, hmax, vmin, vmax
      integer item_number, ierror, slot, entry, probe_dataset
      character*131 buffer
      character*80  string
      character*1  answer
      character*3 str_aspect, str_uniform_contours, str_yflip,
     1            str_log
      logical eof, abort, is_graphics_terminal, lopen,
     1        delayed_message, data_selected, lcontour_error, lannotate
CC END MENU
CC
CC BEGIN ZOOM
       real x1, x2, y1, y2
       logical do_zoom
CC END ZOOM
CC BEGIN PROBE
      integer ierr
      real probe_value, lastx, lasty, delx, dely
      logical do_probe, first_probe
CC END PROBE
c
      real rdum
      integer idum
      character*1 cdum
c
      logical lfirst
      data lfirst/.TRUE./
      save lfirst
c
      lopen = .FALSE.
      lannotate = .FALSE.
      delayed_message = .FALSE.
      contour_values_set = .FALSE.
      do_draw_to_scale = .FALSE.
      cdum = ' '
      do_negate = .FALSE.
      do_log    = .FALSE.
c
      ulen = string_len(pm_length_unit)
c
c     Check for data selected
      data_selected = .FALSE.
      do i=1,n_data_2d
         slot = slot_index_2d(i)
         if( multi_2d(slot)%selected ) data_selected = .TRUE.
      end do
      if( .NOT. data_selected ) then
         if( lbatch ) then
            call plt$error_message('compare_plot_2d',
     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
      if( .not. lsave ) then
         num_titles = 1
         title_dat(1) = ' '
      endif
      if( callable_num_titles .ge. 1 ) then
         num_titles = callable_num_titles
         do i=1,num_titles
            title_dat(i) = callable_titles(i)
         end do
      endif
c
c     Find min/max of entire dataset(s)
      lentire_matrix = .TRUE.
      call pm_minmax_multi_2d(rdum1, rdum2, rdum3, rdum4,
     1                        lentire_matrix,
     1                        min_value, max_value)
c
       overall_hmin =  PM_MAX_FLOAT
       overall_hmax = -PM_MAX_FLOAT
       overall_vmin =  PM_MAX_FLOAT
       overall_vmax = -PM_MAX_FLOAT
       do i=1,n_data_2d
         call pm_get_2d_coord_limits(i, hmin, hmax, vmin, vmax)
         overall_hmin = MIN(overall_hmin,hmin)
         overall_hmax = MAX(overall_hmax,hmax)
         overall_vmin = MIN(overall_vmin,vmin)
         overall_vmax = MAX(overall_vmax,vmax)
      end do
c
      if( .NOT. lsave .and.
     1    (twod_lfirst .or. l_autoscale .or.
     2    (prev_plane_type .ne. plane_type)) ) then
         twod_lfirst = .FALSE.
         twod_h_min =  overall_hmin
         twod_h_max =  overall_hmax
         twod_v_min =  overall_vmin
         twod_v_max =  overall_vmax
         if( twod_h_min .eq. twod_h_max ) then
            if( twod_h_min .eq. 0.0 ) then
               twod_h_min = -0.1
               twod_h_max =  0.1
            else if ( twod_h_min .lt. 0.0 ) then
               twod_h_min = 1.1*twod_h_min
               twod_h_max = 0.9*twod_h_min
            else
               twod_h_min = 0.9*twod_h_min
               twod_h_max = 1.1*twod_h_min
            endif
         endif
         if( twod_v_min .eq. twod_v_max ) then
            if( twod_v_min .eq. 0.0 ) then
               twod_v_min = -0.1
               twod_v_max =  0.1
            else if ( twod_v_min .lt. 0.0 ) then
              twod_v_min = 1.1*twod_v_min
               twod_v_max = 0.9*twod_v_min
            else
               twod_v_min = 0.9*twod_v_min
               twod_v_max = 1.1*twod_v_min
            endif
         endif
         call get_default_limits(twod_h_min, twod_h_max, rdum1, rdum2,
     1      twod_h_tic_major, twod_h_minor_tic_freq, idum1, .FALSE.)
         twod_h_tic_major = twod_h_tic_major*10.0**idum1
         call get_default_limits(twod_v_min, twod_v_max, rdum1, rdum2,
     1      twod_v_tic_major, twod_v_minor_tic_freq, idum1, .FALSE.)
         twod_v_tic_major = twod_v_tic_major*10.0**idum1
         if( do_draw_to_scale ) then
c           Reset plot bounds
c           Magic number 0.9286 is aspect ratio of NDC plot space
c           See plot_contours(), call to gks$set_viewport()
            if( ratio_y .eq. 0.0 ) then
c              Don't have aspect ratio information yet,
c              Use approximate ratio for PostScript device
               aspect_ratio = 0.9286*1.30
            else
               aspect_ratio = 0.9286*(ratio_x/ratio_y)
            endif
            hmin_plot = pm_scale_cm2plot*twod_h_min
            hmax_plot = pm_scale_cm2plot*twod_h_max
            vmin_plot = pm_scale_cm2plot*twod_v_min
            vmax_plot = pm_scale_cm2plot*twod_v_max
            h_tic_major = pm_scale_cm2plot*twod_h_tic_major
            v_tic_major = pm_scale_cm2plot*twod_v_tic_major
            call pm_set_plot_bounds_to_scale(aspect_ratio,'E',
     1               hmin_plot,hmax_plot,
     2               h_tic_major,twod_h_minor_tic_freq,
     3               vmin_plot,vmax_plot,
     4               v_tic_major,twod_v_minor_tic_freq)
            twod_h_min = pm_scale_plot2cm*hmin_plot
            twod_h_max = pm_scale_plot2cm*hmax_plot
            twod_h_tic_major = pm_scale_plot2cm*h_tic_major
            twod_v_min = pm_scale_plot2cm*vmin_plot
            twod_v_max = pm_scale_plot2cm*vmax_plot
            twod_v_tic_major = pm_scale_plot2cm*v_tic_major
         endif
      endif
c
c     Define axis labeling
      if (lsave) then
         horiz_label = twod_h_label
         horiz_units = twod_h_units
         vert_label  = twod_v_label
         vert_units  = twod_v_units
      else if( lfirst .or. prev_plane_type .ne. plane_type) then
         lfirst = .FALSE.
         if( plane_type .eq. XY_PLANE ) then
            horiz_label = 'Length'
            horiz_units = pm_length_unit
         else if( plane_type .eq. YZ_PLANE ) then
            horiz_label = 'Width'
            horiz_units = pm_length_unit
         else if( plane_type .eq. XZ_PLANE ) then
            horiz_label = 'Length'
            horiz_units = pm_length_unit
         else if( plane_type .eq. XT_PLANE ) then
            horiz_label = 'Time'
            horiz_units = 's'
         else if( plane_type .eq. YT_PLANE ) then
            horiz_label = 'Time'
            horiz_units = 's'
         endif
         if( plane_type .eq. XY_PLANE ) then
            vert_label = 'Depth'
            vert_units = pm_length_unit
         elseif( plane_type .eq. YZ_PLANE ) then
            vert_label = 'Depth'
            vert_units = pm_length_unit
         elseif( plane_type .eq. XZ_PLANE ) then
            vert_label = 'Width'
            vert_units = pm_length_unit
         elseif( plane_type .eq. XT_PLANE ) then
            vert_label = 'Length'
            vert_units = pm_length_unit
         elseif( plane_type .eq. YT_PLANE ) then
            vert_label = 'Depth'
            vert_units = pm_length_unit
         endif
      endif
      label = ' '
      units = ' '
      if( n_data_2d .eq. 1 ) then
         slot = slot_index_2d(1)
         label = multi_2d(slot)%quantity
         units = multi_2d(slot)%units
      endif
c
c     Initialization for menu
      hmin_full = pm_scale_cm2plot*overall_hmin
      hmax_full = pm_scale_cm2plot*overall_hmax
      vmin_full = pm_scale_cm2plot*overall_vmin
      vmax_full = pm_scale_cm2plot*overall_vmax
      hmin_plot = pm_scale_cm2plot*twod_h_min
      hmax_plot = pm_scale_cm2plot*twod_h_max
      vmin_plot = pm_scale_cm2plot*twod_v_min
      vmax_plot = pm_scale_cm2plot*twod_v_max
      h_tic_major = pm_scale_cm2plot*twod_h_tic_major
      v_tic_major = pm_scale_cm2plot*twod_v_tic_major
      lentire_matrix = .FALSE.
      call pm_minmax_multi_2d(twod_h_min, twod_h_max,
     1                        twod_v_min, twod_v_max,
     1                        lentire_matrix,
     1                        plot_min_value, plot_max_value)
c
      if( lsave ) then
         do_yflip  = .NOT. twod_invert
         do 110 i=1,twod_num_contours
            contour_values(i) = twod_contours(i)
110      continue
      else
c        Compute contours for overall dataset
         twod_num_contours = MAX_CONTOURS
         if( plot_min_value .gt. plot_max_value ) then
c           Current (X,Y) bounds are outide entire dataset
            call compute_contours(min_value, max_value,
     1           do_log, do_negate, .FALSE.,
     2           twod_num_contours, contour_values)
         else
            call compute_contours(plot_min_value, plot_max_value,
     1           do_log, do_negate, .FALSE.,
     2           twod_num_contours, contour_values)
         endif
c
c        Default for flipping Y axis: only if all datasets are ASCII...
         do_yflip  = .FALSE.
         do i=1,n_data_2d
            slot = slot_index_2d(i)
            if( multi_2d(slot)%datatype .eq. DT_FAIM .or.
     1          multi_2d(slot)%datatype .eq. DT_MCP2D .or.
     2          multi_2d(slot)%datatype .eq. DT_ASCII ) then
c              Don't do anything
            else
               do_yflip  = .TRUE.
            endif
         end do
      endif
c     Sort in ascending order (MUST BE DONE!!!)
      call sort_contours(contour_values,twod_num_contours)
c
      do_uniform_contours = .FALSE.
c
      if( do_uniform_contours ) then
         str_uniform_contours = 'yes'
      else
         str_uniform_contours = ' no'
      endif
c
      if( do_draw_to_scale ) then
         str_aspect = 'yes'
      else
         str_aspect = ' no'
      endif
c
      if( do_yflip ) then
         str_yflip = 'yes'
      else
         str_yflip = ' no'
      endif
c
      if( do_log ) then
         str_log = 'yes'
      else
         str_log = ' no'
      endif
c
      probe_dataset = 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_OVERLAY_COMMAND,
     1     pm_scale_cm2plot*twod_h_min, pm_scale_cm2plot*twod_h_max,
     2     idum,
     3     pm_scale_cm2plot*twod_h_tic_major, twod_h_minor_tic_freq,
     4     horiz_units, horiz_label, .FALSE., .FALSE., .TRUE.,
     5     pm_scale_cm2plot*twod_v_min, pm_scale_cm2plot*twod_v_max,
     6     idum,
     7     pm_scale_cm2plot*twod_v_tic_major, twod_v_minor_tic_freq,
     8     vert_units, vert_label, .FALSE., .FALSE., .TRUE.,
     9     rdum, rdum,
     A     idum, rdum, idum, cdum,
     B     cdum, .FALSE., .FALSE., .FALSE.)
c
c        Stuff data information into save structures
         call copy_multi_2d_to_save(contour_values, twod_num_contours,
     1     .NOT. do_yflip, do_uniform_contours)
c
         call plot_multi_contours(
     1     twod_h_min, twod_h_max,
     2     twod_h_tic_major, twod_h_minor_tic_freq,
     3     horiz_label, horiz_units,
     4     twod_v_min, twod_v_max,
     5     twod_v_tic_major, twod_v_minor_tic_freq,
     6     vert_label, vert_units,
     7     twod_num_contours, contour_values, do_uniform_contours,
     8     probe_dataset, do_yflip, lsave, lopen, lannotate, ierror)
         plot_device = previous_plot_device
         return
      endif
c
c     Main return point for menu
c
4000  continue
c     Write out menu form:
      call clear_screen
      call display_string(1,1 ,'1)  Title:',10, 0)
c
      call display_string(3 ,1 ,'HORIZONTAL AXIS PARAMETERS:',27,0)
      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)  Label: ',11, 0)
      call display_string(10,1 ,'6)  Units: ',11, 0)
c
      call display_string(12,1 ,'VERTICAL AXIS PARAMETERS:',25,0)
      call display_string(13,1 ,'Data range:',11,0)
      call display_string(13,24,' to',3,0)
      call display_string(15,1 ,'7)  Plot range:',15,0)
      call display_string(15,24,' to',3,0)
      call display_string(16,1 ,'8)  Major tics at:',18,0)
      call display_string(17,1 ,'9)  Minor tic freq:',19,0)
      call display_string(18,1 ,'10) Label: ',11, 0)
      call display_string(19,1 ,'11) Units: ',11, 0)
c
      call display_string(3 ,42,'Data range:            to',25,0)
      call display_string(4 ,42,'Plot range:            to',25,0)
      call display_string(5 ,42,'Sample quantity: ',17,0)
      call display_string(7 ,42,'12) Global',10,0)
      call display_string(8 ,42,'    contours',12,0)
      call display_string(7 ,56,'13) Use global',14,0)
      call display_string(8 ,56,'    contours:',13,0)
      call display_string(9 ,56,'14) Sample dataset:',20,0)
      call display_string(10,56,'15) Draw to scale:',18,0)
      call display_string(11,56,'16) Positive y down:',20,0)
      call display_string(12,56,'17) Recompute contours:',23,0)
      call display_string(13,56,'18) Use log contours:',21,0)
c
      call display_string(23,1 ,'Messages:',9,0)
c
c     Cut line
c     Title
      call display_string(1 , 11, title_dat(1), 69, 1)
c
c     X Range
      call display_real  (4 , 15, hmin_full, 1)
      call display_real  (4 , 28, hmax_full, 1)
      call display_real  (6 , 15, hmin_plot, 1)
      call display_real  (6 , 28, hmax_plot, 1)
c     X major tics
      call display_real  (7 , 23, h_tic_major, 1)
c     X minor tic freq
      call display_int   (8 , 27, twod_h_minor_tic_freq, 5)
c     Horizontal label
      call display_string(9 , 12, horiz_label, 27, 1)
c     Horizontal Units
      call display_string(10, 12, horiz_units, 12, 1)
c
c     Y Range
      call display_real  (13, 15, vmin_full, 1)
      call display_real  (13, 28, vmax_full, 1)
      call display_real  (15, 15, vmin_plot, 1)
      call display_real  (15, 28, vmax_plot, 1)
c     Y major tics
      call display_real  (16, 23, v_tic_major, 1)
c     Y minor tic freq
      call display_int   (17, 27, twod_v_minor_tic_freq, 5)
c     Vertical label
      call display_string(18, 12, vert_label, 27, 1)
c     Vertical Units
      call display_string(19, 12, vert_units, 12, 1)
c
c     Data min/max
      call display_real  (3 , 54, min_value, 1)
      call display_real  (3 , 68, max_value, 1)
      if( plot_min_value .gt. plot_max_value ) then
         call display_string(4 ,54,
     1                       'No grid points found   ',23,1)
         else
            call display_string(4 ,53,'            to  ',26,0)
            call display_real  (4 ,54, plot_min_value, 1)
            call display_real  (4 ,68, plot_max_value, 1)
      endif
c     Probe dataset
      slot = slot_index_2d(probe_dataset)
      ilen = string_len(multi_2d(slot)%quantity)
      call display_string(5,60,multi_2d(slot)%quantity(1:ilen),ilen,1)
c
c     Contours
      do 200 i=1,twod_num_contours
         call display_real(8+i, 43, contour_values(i), 1)
200   continue
      do 210 i=twod_num_contours+1,10
         call display_string(8+i, 43,'         ', 9, 0)
210   continue
c
c     Uniform contours
      call display_string(7 ,77,str_uniform_contours, 3,1)
c     Probe dataset (entry) number
      call display_int(9 ,77,probe_dataset,3)
c     Maintain aspect ratio flag
      call display_string(10,77,str_aspect, 3,1)
c     Positive Y axis down
      call display_string(11,77,str_yflip, 3,1)
c     Use log contours
      call display_string(13,77,str_log, 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')
      call read_from_input_area(string, eof)
      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_OVERLAY_COMMAND,
     1     pm_scale_cm2plot*twod_h_min, pm_scale_cm2plot*twod_h_max,
     2     idum,
     3     pm_scale_cm2plot*twod_h_tic_major, twod_h_minor_tic_freq,
     4     horiz_units, horiz_label, .FALSE., .FALSE., .TRUE.,
     5     pm_scale_cm2plot*twod_v_min, pm_scale_cm2plot*twod_v_max,
     6     idum,
     7     pm_scale_cm2plot*twod_v_tic_major, twod_v_minor_tic_freq,
     8     vert_units, vert_label, .FALSE., .FALSE., .TRUE.,
     9     rdum, rdum,
     A     idum, rdum, idum, cdum,
     B     cdum, .FALSE., .FALSE., .FALSE.)
c
c        Stuff data information into save structures
         call copy_multi_2d_to_save(contour_values, twod_num_contours,
     1     .NOT. do_yflip, do_uniform_contours)
c
         call plot_multi_contours(
     1     twod_h_min, twod_h_max,
     2     twod_h_tic_major, twod_h_minor_tic_freq,
     3     horiz_label, horiz_units,
     4     twod_v_min, twod_v_max,
     5     twod_v_tic_major, twod_v_minor_tic_freq,
     6     vert_label, vert_units,
     7     twod_num_contours, contour_values, do_uniform_contours,
     8     probe_dataset, do_yflip, lsave, lopen, lannotate, ierror)
         if( ierror .eq. PLT$NORMAL ) then
            lopen = .TRUE.
         else
            call plt$get_error_message(ierror,string)
            delayed_message = .TRUE.
         endif
         goto 4000
c
      else if( string(1:1) .eq. 'Z' .or. string(1:1) .eq. 'U') then
         if( lopen .or. string(1:1) .eq. 'U' ) then
            if( is_graphics_terminal(current_window_id) .and.
     1          string(1:1) .eq. 'Z' ) then
               call display_message(
     1           'Sorry: cannot zoom on this plot device')
               goto 2000
            else
               if( string(1:1) .eq. 'Z' ) then
                  if( string(2:2) .eq. '2' ) then
                     delta = 0.25*(twod_h_max - twod_h_min)
                     x1 = twod_h_min + delta
                     x2 = twod_h_max - delta
                     delta = 0.25*(twod_v_max - twod_v_min)
                     y1 = twod_v_min + delta
                     y2 = twod_v_max - delta
                     do_zoom = .TRUE.
                  else
                  call display_prompt(
     1     'Enter lower left, upper right. ' // PM_CANCEL_BUTTON //
     1     ' button cancels input')
c                 Magic constants are values to GKS$SET_VIEWPORT in
c                 subroutine plot_contours.
                  call window_zoom(current_window_id,
     1               twod_h_min, twod_h_max, twod_v_min, twod_v_max,
     2               0.175, 0.825, 0.10, 0.80, do_yflip,
     3               x1, x2, y1, y2, do_zoom)
                  endif
               else
                  if( string(2:2) .eq. '2' ) then
                     delta = 0.5*(twod_h_max - twod_h_min)
                     x1 = twod_h_min - delta
                     x2 = twod_h_max + delta
                     delta = 0.5*(twod_v_max - twod_v_min)
                     y1 = twod_v_min - delta
                     y2 = twod_v_max + delta
                  else
                     x1 = pm_scale_plot2cm*hmin_full
                     x2 = pm_scale_plot2cm*hmax_full
                     y1 = pm_scale_plot2cm*vmin_full
                     y2 = pm_scale_plot2cm*vmax_full
                  endif
                  do_zoom = .TRUE.
               endif
               if( do_zoom ) then
                  twod_h_min = x1
                  twod_h_max = x2
                  twod_v_min = y1
                  twod_v_max = y2
                  hmin_plot  = pm_scale_cm2plot*twod_h_min
                  hmax_plot  = pm_scale_cm2plot*twod_h_max
                  vmin_plot  = pm_scale_cm2plot*twod_v_min
                  vmax_plot  = pm_scale_cm2plot*twod_v_max
                  if( do_draw_to_scale ) then
c                    Magic number 0.9286 is aspect ratio of NDC plot space
c                    See plot_contours(), call to gks$set_viewport()
                     if( ratio_y .eq. 0.0 ) then
c                       Don't have aspect ratio information yet,
c                       Use approximate ratio for PostScript device
                        aspect_ratio = 0.9286*1.30
                     else
                        aspect_ratio = 0.9286*(ratio_x/ratio_y)
                     endif
                     call pm_set_plot_bounds_to_scale(aspect_ratio,'E',
     1                  hmin_plot,hmax_plot,
     2                  h_tic_major,twod_h_minor_tic_freq,
     3                  vmin_plot,vmax_plot,
     4                  v_tic_major,twod_v_minor_tic_freq)
                     twod_h_min = pm_scale_plot2cm*hmin_plot
                     twod_h_max = pm_scale_plot2cm*hmax_plot
                     twod_h_tic_major = pm_scale_plot2cm*h_tic_major
                     twod_v_min = pm_scale_plot2cm*vmin_plot
                     twod_v_max = pm_scale_plot2cm*vmax_plot
                     twod_v_tic_major = pm_scale_plot2cm*v_tic_major
                  endif
c
c                 Get tic mark parameters
CC                  call get_default_limits(twod_h_min, twod_h_max,
CC     1                 rdum1, rdum2,
CC     2                 twod_h_tic_major, twod_h_minor_tic_freq,
CC     3                 idum1, .FALSE.)
CC                  twod_h_tic_major = twod_h_tic_major*10.0**idum1
CC                  h_tic_major = pm_scale_cm2plot*twod_h_tic_major
CC  POSSIBLE FIX TO GET BETTER TIC MARKS FOR SOME AXIS SETTINGS
                  call get_default_limits(hmin_plot, hmax_plot,
     1                 rdum1, rdum2,
     2                 h_tic_major, twod_h_minor_tic_freq,
     3                 idum1, .FALSE.)
                  h_tic_major = h_tic_major*10.0**idum1
                  twod_h_tic_major = pm_scale_plot2cm*h_tic_major
                  call get_default_limits(twod_v_min, twod_v_max,
     1                 rdum1, rdum2,
     2                 twod_v_tic_major, twod_v_minor_tic_freq,
     3                 idum1, .FALSE.)
                  twod_v_tic_major = twod_v_tic_major*10.0**idum1
                  v_tic_major = pm_scale_cm2plot*twod_v_tic_major
c                 Stuff axis information into save structures
                  call copy_axis_to_save(PLOT_OVERLAY_COMMAND,
     1     pm_scale_cm2plot*twod_h_min, pm_scale_cm2plot*twod_h_max,
     2     idum,
     3     pm_scale_cm2plot*twod_h_tic_major, twod_h_minor_tic_freq,
     4     horiz_units, horiz_label, .FALSE., .FALSE., .TRUE.,
     5     pm_scale_cm2plot*twod_v_min, pm_scale_cm2plot*twod_v_max,
     6     idum,
     7     pm_scale_cm2plot*twod_v_tic_major, twod_v_minor_tic_freq,
     8     vert_units, vert_label, .FALSE., .FALSE., .TRUE.,
     9     rdum, rdum,
     A     idum, rdum, idum, cdum,
     B     cdum, .FALSE., .FALSE., .FALSE.)
c
c                 Stuff data information into save structures
                  call copy_multi_2d_to_save(contour_values,
     1            twod_num_contours,.NOT. do_yflip,do_uniform_contours)
c
                  call plot_multi_contours(
     1            twod_h_min, twod_h_max,
     2            twod_h_tic_major, twod_h_minor_tic_freq,
     3            horiz_label, horiz_units,
     4            twod_v_min, twod_v_max,
     5            twod_v_tic_major, twod_v_minor_tic_freq,
     6            vert_label, vert_units,
     7            twod_num_contours,contour_values,do_uniform_contours,
     8            probe_dataset, do_yflip,
     9            lsave, lopen, lannotate, ierror)
                  if( ierror .eq. PLT$NORMAL ) then
                     lopen = .TRUE.
                  endif
c                 Update dataset plot min/max limits
                  lentire_matrix = .FALSE.
                  call pm_minmax_multi_2d(twod_h_min, twod_h_max,
     1                                 twod_v_min, twod_v_max,
     1                                 lentire_matrix,
     1                                 plot_min_value, plot_max_value)
               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_OVERLAY_COMMAND,
     1     pm_scale_cm2plot*twod_h_min, pm_scale_cm2plot*twod_h_max,
     2     idum,
     3     pm_scale_cm2plot*twod_h_tic_major, twod_h_minor_tic_freq,
     4     horiz_units, horiz_label, .FALSE., .FALSE., .TRUE.,
     5     pm_scale_cm2plot*twod_v_min, pm_scale_cm2plot*twod_v_max,
     6     idum,
     7     pm_scale_cm2plot*twod_v_tic_major, twod_v_minor_tic_freq,
     8     vert_units, vert_label, .FALSE., .FALSE., .TRUE.,
     9     rdum, rdum,
     A     idum, rdum, idum, cdum,
     B     cdum, .FALSE., .FALSE., .FALSE.)
c
c        Stuff data information into save structures
         call copy_multi_2d_to_save(contour_values, twod_num_contours,
     1     .NOT. do_yflip, do_uniform_contours)
c
         call plot_multi_contours(
     1     twod_h_min, twod_h_max,
     2     twod_h_tic_major, twod_h_minor_tic_freq,
     3     horiz_label, horiz_units,
     4     twod_v_min, twod_v_max,
     5     twod_v_tic_major, twod_v_minor_tic_freq,
     6     vert_label, vert_units,
     7     twod_num_contours, contour_values, do_uniform_contours,
     8     probe_dataset, do_yflip, lsave, lopen, lannotate, ierror)
         plot_device = previous_plot_device
         if( ierr .eq. PLT$NORMAL ) then
            ilen = string_len(plot_file_name)
            write(string,600,err=601) plot_file_name(1:ilen)
600         format('Plot file ',a' created')
601         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 plot 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. 'K') then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: no window functions on this plot 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 plot 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              Magic constants are values to GKS$SET_VIEWPORT in
c              subroutine plot_contours.
c              Top of sample loop...
720            continue
               call get_probe_location(current_window_id,
     1            twod_h_min, twod_h_max, twod_v_min, twod_v_max,
     2            0.175, 0.825, 0.10, 0.80, do_yflip,
     3            x1, y1, do_probe)
               if( .NOT. do_probe ) goto 710
               call probe_multi_interpolate(probe_dataset, x1, y1,
     1                                      probe_value, ierr)
               if( ierr .eq. 0 ) then
                  if( first_probe ) then
                  write(string,700,err=701) probe_value,
     1               pm_scale_cm2plot*x1, pm_scale_cm2plot*y1
700               format('Value: ',1pe11.4,' at (',1pe11.4,
     1                   ',',1pe11.4,')')
701               continue
                  first_probe = .FALSE.
                  else
                  delx = x1 - lastx
                  dely = y1 - lasty
                  write(string,702,err=703) probe_value,
     1               pm_scale_cm2plot*x1, pm_scale_cm2plot*y1,
     2               pm_scale_cm2plot*delx, pm_scale_cm2plot*dely
702               format('Value: ',1pe11.4,' at (',1pe11.4,',',
     1                   1pe11.4,') Change(',1pe11.4,',',1pe11.4,')')
703               continue
                  endif
                  lastx = x1
                  lasty = y1
                  call display_message_quietly(string)
               else
                  call display_message(
     1               'Data probe outside array bounds')
               endif
               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 if( string(1:1) .eq. 'R') then ! Repaint screen
         goto 4000
      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)')
            call display_string(1 , 11, title_dat(1), 69, 1)
c
         elseif ( item_number .eq. 2 ) then
c           X axis min/max
            real_vec(1) = hmin_plot
            real_vec(2) = hmax_plot
            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
               hmin_plot = real_vec(1)
               hmax_plot = real_vec(2)
               if( do_draw_to_scale ) then
c                 Magic number 0.9286 is aspect ratio of NDC plot space
c                 See plot_contours(), call to gks$set_viewport()
                  if( ratio_y .eq. 0.0 ) then
c                    Don't have aspect ratio information yet,
c                    Use approximate ratio for PostScript device
                     aspect_ratio = 0.9286*1.30
                  else
                     aspect_ratio = 0.9286*(ratio_x/ratio_y)
                  endif
                  call pm_set_plot_bounds_to_scale(aspect_ratio,'X',
     1               hmin_plot,hmax_plot,
     2               h_tic_major,twod_h_minor_tic_freq,
     3               vmin_plot,vmax_plot,
     4               v_tic_major,twod_v_minor_tic_freq)
c                 Display Y axis parameters
                  twod_v_min = pm_scale_plot2cm*vmin_plot
                  twod_v_max = pm_scale_plot2cm*vmax_plot
                  call display_real(15, 15, vmin_plot, 1)
                  call display_real(15, 28, vmax_plot, 1)
                  twod_v_tic_major = pm_scale_plot2cm*v_tic_major
                  call display_real(16, 23, v_tic_major, 1)
                  call display_int (17, 27, twod_v_minor_tic_freq, 5)
               else
c                 Generate compatible tic mark parameters
CC                  call get_default_limits(twod_h_min, twod_h_max,
CC     1              rdum1, rdum2,
CC     2              twod_h_tic_major, twod_h_minor_tic_freq,
CC     3              idum1, .FALSE.)
CC               twod_h_tic_major = twod_h_tic_major*10.0**idum1
CC               h_tic_major = pm_scale_cm2plot*twod_h_tic_major
CC     POSSIBLE FIX TO GET BETTER TIC MARKS FOR SOME AXIS SETTINGS
                  call get_default_limits(hmin_plot, hmax_plot,
     1                 rdum1, rdum2,
     2                 h_tic_major, twod_h_minor_tic_freq,
     3                 idum1, .FALSE.)
                  h_tic_major = h_tic_major*10.0**idum1
               endif
c              Display X axis parameters (for both cases above)
               twod_h_min = pm_scale_plot2cm*hmin_plot
               twod_h_max = pm_scale_plot2cm*hmax_plot
               call display_real(6 , 15, hmin_plot, 1)
               call display_real(6 , 28, hmax_plot, 1)
               twod_h_tic_major = pm_scale_plot2cm*h_tic_major
               call display_real(7 , 23, h_tic_major, 1)
               call display_int (8 , 27, twod_h_minor_tic_freq, 5)
c              Update dataset plot min/max limits
               lentire_matrix = .FALSE.
               call pm_minmax_multi_2d(twod_h_min, twod_h_max,
     1                                 twod_v_min, twod_v_max,
     1                                 lentire_matrix,
     1                                 plot_min_value, plot_max_value)
               if( plot_min_value .gt. plot_max_value ) then
                  call display_string(4 ,54,
     1               'No grid points found   ',23,1)
               else
                  call display_string(4 ,53,'            to  ',26,0)
                  call display_real  (4 ,54, plot_min_value, 1)
                  call display_real  (4 ,68, plot_max_value, 1)
               endif
            endif
c
         elseif ( item_number .eq. 3 ) then
c           X axis major tic increment
            real_vec(1) = h_tic_major
            call get_reals_with_prompt(real_vec, 1,
     1              'Enter major tic increment for horizontal axis')
            delta = hmax_plot - hmin_plot
            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
               h_tic_major = real_vec(1)
               twod_h_tic_major = pm_scale_plot2cm*h_tic_major
               call display_real  (7 , 23, h_tic_major, 1)
            endif
c
         elseif ( item_number .eq. 4 ) then
c           Minor tic frequency
            int_vec(1) = twod_h_minor_tic_freq
            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
               twod_h_minor_tic_freq = int_vec(1)
               call display_int (8 , 27, twod_h_minor_tic_freq, 5)
            endif
c
         elseif ( item_number .eq. 5 ) then
            call get_string_with_prompt(horiz_label,
     1              'Enter horizontal label (may use mathmode)')
            call display_string( 9, 12, horiz_label, 27, 1)
c
         elseif ( item_number .eq. 6 ) then
            call get_string_with_prompt(horiz_units,
     1              'Enter units (may use mathmode)')
            call display_string(10, 12, horiz_units, 12, 1)
c
         elseif ( item_number .eq. 7 ) then
c           Y axis min/max
            real_vec(1) = vmin_plot
            real_vec(2) = vmax_plot
            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
               vmin_plot = real_vec(1)
               vmax_plot = real_vec(2)
               if( do_draw_to_scale ) then
c                 Magic number 0.9286 is aspect ratio of NDC plot space
c                 See plot_contours(), call to gks$set_viewport()
                  if( ratio_y .eq. 0.0 ) then
c                    Don't have aspect ratio information yet,
c                    Use approximate ratio for PostScript device
                     aspect_ratio = 0.9286*1.30
                  else
                     aspect_ratio = 0.9286*(ratio_x/ratio_y)
                  endif
                  call pm_set_plot_bounds_to_scale(aspect_ratio,'Y',
     1               hmin_plot,hmax_plot,
     2               h_tic_major,twod_h_minor_tic_freq,
     3               vmin_plot,vmax_plot,
     4               v_tic_major,twod_v_minor_tic_freq)
c                 Display X axis parameters
                  twod_h_min = pm_scale_plot2cm*hmin_plot
                  twod_h_max = pm_scale_plot2cm*hmax_plot
                  call display_real(6 , 15, hmin_plot, 1)
                  call display_real(6 , 28, hmax_plot, 1)
                  twod_h_tic_major = pm_scale_plot2cm*h_tic_major
                  call display_real(7 , 23, h_tic_major, 1)
                  call display_int (8 , 27, twod_h_minor_tic_freq, 5)
               else
c                 Generate compatible tic mark parameters
C                 POSSIBLE FIX TO GET BETTER TIC MARKS FOR SOME AXIS SETTINGS
                  call get_default_limits(vmin_plot, vmax_plot,
     1                 rdum1, rdum2,
     2                 v_tic_major, twod_v_minor_tic_freq,
     3                 idum1, .FALSE.)
                  v_tic_major = v_tic_major*10.0**idum1
               endif
c              Display Y axis parameters (for both cases above)
               twod_v_min = pm_scale_plot2cm*vmin_plot
               twod_v_max = pm_scale_plot2cm*vmax_plot
               call display_real(15, 15, vmin_plot, 1)
               call display_real(15, 28, vmax_plot, 1)
               twod_v_tic_major = pm_scale_plot2cm*v_tic_major
               call display_real(16, 23, v_tic_major, 1)
               call display_int (17, 27, twod_v_minor_tic_freq, 5)
c              Update dataset plot min/max limits
               lentire_matrix = .FALSE.
               call pm_minmax_multi_2d(twod_h_min, twod_h_max,
     1                                 twod_v_min, twod_v_max,
     1                                 lentire_matrix,
     1                                 plot_min_value, plot_max_value)
               if( plot_min_value .gt. plot_max_value ) then
                  call display_string(4 ,54,
     1               'No grid points found   ',23,1)
               else
                  call display_string(4 ,53,'            to  ',26,0)
                  call display_real  (4 ,54, plot_min_value, 1)
                  call display_real  (4 ,68, plot_max_value, 1)
               endif
            endif
c
         elseif ( item_number .eq. 8 ) then
c           Major Y tic increment
            real_vec(1) = v_tic_major
            call get_reals_with_prompt(real_vec, 1,
     1           'Enter major tic increment for vertical axis')
            delta = vmax_plot - vmin_plot
            if( real_vec(1) .gt. delta .or.
     1          real_vec(1) .le. 0.0 ) then
                call display_message('Illegal value')
            else
               v_tic_major = real_vec(1)
               twod_v_tic_major = pm_scale_plot2cm*v_tic_major
               call display_real  (16, 23, v_tic_major, 1)
            endif
c
         elseif ( item_number .eq. 9 ) then
c           Minor tic frequency
            int_vec(1) = twod_v_minor_tic_freq
            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
               twod_v_minor_tic_freq = int_vec(1)
               call display_int (17, 27, twod_v_minor_tic_freq, 5)
            endif
c
         elseif ( item_number .eq. 10 ) then
            call get_string_with_prompt(vert_label,
     1              'Enter vertical label (may use mathmode)')
            call display_string(18, 12, vert_label, 27, 1)
c
         elseif ( item_number .eq. 11 ) then
            call get_string_with_prompt(vert_units,
     1              'Enter units (may use mathmode)')
            call display_string(19, 12, vert_units, 12, 1)
c
         elseif ( item_number .eq. 12 ) then
c           Contours
            if( do_uniform_contours ) then
               call display_prompt('Enter up to 9 contour values')
               call read_from_input_area(buffer, eof)
               if( .NOT. eof ) then
                  call upcase(buffer)
                  call parse_contours(buffer, MAX_CONTOURS,
     1                 int_vec(1), real_vec, lcontour_error)
                  if( lcontour_error ) then
                     call display_message(
     1               'Invalid contour specification')
                  else
                     twod_num_contours = int_vec(1)
                     do i=1,twod_num_contours
                        contour_values(i) = real_vec(i)
                     end do
c                    Sort in ascending order (MUST BE DONE!!!)
                     call sort_contours(contour_values,
     1                                  twod_num_contours)
                     do i=1,twod_num_contours
                        call display_real(8+i, 43, contour_values(i), 1)
                     end do
                     do i=twod_num_contours+1,10
                       call display_string(8+i, 43,'         ', 9, 0)
                     end do
                  endif
               endif
            else
               call display_message('Individual contours are enabled')
            endif
         elseif ( item_number .eq. 13 ) then
c           Same contours on all plots
            do_uniform_contours = .NOT. do_uniform_contours
            if( do_uniform_contours ) then
               str_uniform_contours = 'yes'
            else
               str_uniform_contours = ' no'
            endif
            call display_string(7 ,77,str_uniform_contours, 3,1)
         elseif ( item_number .eq. 14 ) then
c           Probe dataset
            int_vec(1) = probe_dataset
            call get_ints_with_prompt(int_vec, 1,
     1 'Enter dataset to use for sample function and plotting sidebar')
            if( int_vec(1) .lt. 1 .or. int_vec(1) .gt. n_data_2d ) then
              call display_message('Illegal value')
            else
               probe_dataset = int_vec(1)
               call display_int(9,77,probe_dataset,3)
            endif
            slot = slot_index_2d(probe_dataset)
            ilen = string_len(multi_2d(slot)%quantity)
            call display_string(5,60,multi_2d(slot)%quantity(1:ilen),
     1                          ilen,1)
         elseif ( item_number .eq. 15 ) then
            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(10,77,str_aspect, 3, 1)
            if( do_draw_to_scale ) then
c              Reset plot bounds
c              Magic number 0.9286 is aspect ratio of NDC plot space
c              See plot_contours(), call to gks$set_viewport()
               if( ratio_y .eq. 0.0 ) then
c                 Don't have aspect ratio information yet,
c                 Use approximate ratio for PostScript device
                  aspect_ratio = 0.9286*1.30
               else
                  aspect_ratio = 0.9286*(ratio_x/ratio_y)
               endif
               call pm_set_plot_bounds_to_scale(aspect_ratio,'E',
     1                  hmin_plot,hmax_plot,
     2                  h_tic_major,twod_h_minor_tic_freq,
     3                  vmin_plot,vmax_plot,
     4                  v_tic_major,twod_v_minor_tic_freq)
c              Display X axis parameters
               twod_h_min = pm_scale_plot2cm*hmin_plot
               twod_h_max = pm_scale_plot2cm*hmax_plot
               call display_real(6 , 15, hmin_plot, 1)
               call display_real(6 , 28, hmax_plot, 1)
               twod_h_tic_major = pm_scale_plot2cm*h_tic_major
               call display_real(7 , 23, h_tic_major, 1)
               call display_int (8 , 27, twod_h_minor_tic_freq, 5)
c              Display Y axis parameters
               twod_v_min = pm_scale_plot2cm*vmin_plot
               twod_v_max = pm_scale_plot2cm*vmax_plot
               call display_real(15, 15, vmin_plot, 1)
               call display_real(15, 28, vmax_plot, 1)
               twod_v_tic_major = pm_scale_plot2cm*v_tic_major
               call display_real(16, 23, v_tic_major, 1)
               call display_int (17, 27, twod_v_minor_tic_freq, 5)
            endif
         elseif ( item_number .eq. 16 ) then
            do_yflip = .NOT. do_yflip
            if( do_yflip ) then
              str_yflip = 'yes'
            else
              str_yflip = ' no'
            endif
            call display_string(11,77,str_yflip, 3, 1)
         elseif ( item_number .eq. 17 .or. item_number .eq. 18) then
c           Recompute contours on current bounds
            if( do_uniform_contours ) then
c              Switch log contour flag, then recompute contours
               if( item_number .eq. 18 ) then
                  do_log = .NOT. do_log
                  if ( do_log ) then
                     str_log = 'yes'
                  else
                     str_log = ' no'
                  endif
                  call display_string(13,77,str_log, 3,1)
               endif
               twod_num_contours = 9
               do_negate = .FALSE.
               twod_num_contours = MAX_CONTOURS
               if( plot_min_value .gt. plot_max_value ) then
c                 Current (X,Y) bounds are outide entire dataset
                  call compute_contours(min_value, max_value,
     1               do_log, do_negate, .FALSE.,
     2               twod_num_contours, contour_values)
               else
                  call compute_contours(plot_min_value, plot_max_value,
     1               do_log, do_negate, .FALSE.,
     2               twod_num_contours, contour_values)
               endif
c              Sort in ascending order (MUST BE DONE!!!)
               call sort_contours(contour_values,twod_num_contours)
               do i=1,twod_num_contours
                  call display_real(8+i, 43, contour_values(i), 1)
               end do
               do i=twod_num_contours+1,10
                  call display_string(8+i, 43,'         ', 9, 0)
               end do
            else
               call display_message('Individual contours are enabled')
            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
      return
      end
c
      subroutine plot_multi_contours(
     1     h_win_min, h_win_max, h_tic_major, h_minor_tic_freq,
     2     h_label, h_units,
     3     v_win_min, v_win_max, v_tic_major, v_minor_tic_freq,
     4     v_label, v_units,
     5     num_contours, contour_values, do_uniform_contours,
     6     probe_dataset, do_yflip,
     7     add_notes_from_save, ws_window_open,
     8     add_notes_from_file, ierror)
c
      implicit none
c
c     Subroutine which actually calls graphics routines.
c
c     NOTE: Since many simulators use inverted vertical coordinates,
c     we must flip the vertical coordinate value of anything pertaining to
c     the plot window. The transformation is (v_sum - actual_v_value).
c     We now allow a natural vertical coordinate, via the do_yflip
c     parameter.
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, clipped outside the plot window, have terrible
c     aspect ratios, etc. etc. etc.
c
      logical do_uniform_contours, do_yflip,
     1        add_notes_from_save, ws_window_open,
     2        add_notes_from_file
      character*(*) h_label, v_label, h_units, v_units
      integer num_contours, ierror, probe_dataset
      integer h_minor_tic_freq, v_minor_tic_freq
      real contour_values(1)
      real h_win_min, h_win_max, h_tic_major
      real v_win_min, v_win_max, v_tic_major
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'io.cmn'
      include 'units.cmn'
      include 'defaults.cmn'
      include 'compare_2d.cmn'
      include 'title.cmn'
c
c     locals
c
      real tmp_contours(MAX_CONTOURS)
      logical lerr, exp_flag, is_vtxxx
      integer i, j, ilen, string_len, window_id, return_stat,
     1        ilen2, itot, slot, tmp_num_contours, do_log
      integer get_new_window_id
      character*1 response
      character*8 time_string
      character*12 num_label
      character*80 work_string
      real junction(1), v_sum, h_cur, v_cur, delta, rmax, ticvalue
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
      if( do_yflip ) then
         v_sum = (v_win_max+v_win_min)
      else
c        ??? Probably not right, but "OK" for use with non-simulator
c        data. jvf
         v_sum = v_win_min
      endif
c
c     open GKS
c
      window_id = get_new_window_id(PLT$2D_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
         ierror = PLT$TOO_MANY_WINDOWS
         if( lbatch .or. lrestore ) then
            call plt$get_error_message(ierror,work_string)
            ilen = string_len(work_string)
            call plt$error_message('plot_contours', 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('plot_contours', work_string(1:ilen))
         endif
         ierror = return_stat
         return
      endif
c
c     Find maximum in plot area for side bar label later on
      do_log = .FALSE.
      call find_max_in_multi_subset(probe_dataset,
     1        h_win_min, h_win_max, v_win_min, v_win_max,
     1        do_log, rmax)
c
c     Activate plot device, clear it, start segment on workstations
      call activate_plot_device(window_id, .TRUE., .TRUE.)
      call pop_vs_window(window_id)
      call init_device_depend
c
      call GKS$SET_WINDOW(1,h_win_min,h_win_max,v_win_min,v_win_max)
c     create virtual plot window at 65% of the screen (NDC space)
      call GKS$SET_VIEWPORT(1, 0.175*ratio_x, 0.825*ratio_x,
     1                         0.1*ratio_y, 0.80*ratio_y)
      call GKS$SET_CLIPPING(GKS$K_CLIP)
      call GKS$SELECT_XFORM(1)
c
      do i=1,n_data_2d
         slot = slot_index_2d(i)
         if( multi_2d(slot).selected ) then
         if( do_uniform_contours ) then
            tmp_num_contours = num_contours
            do j=1,tmp_num_contours
               tmp_contours(j) = contour_values(j)
            end do
         else
            tmp_num_contours = multi_2d(slot)%num_contours
            do j=1,tmp_num_contours
               tmp_contours(j) = multi_2d(slot)%contours(j)
            end do
         endif
         if( multi_2d(slot)%color .and. can_do_color_contours ) then
c           Set colormap for color contours
            call set_color_contour_colormap
            if( multi_2d(slot)%data_is_trigrid ) then
               call color_contours_tri(multi_2d(slot)%data,
     1            multi_2d(slot)%cord, multi_2d(slot)%nop,
     1            multi_2d(slot)%tri_ne, multi_2d(slot)%tri2reg,
     1            multi_2d(slot)%imat,
     1            multi_2d(slot)%tri_npts,
     1            tmp_contours, tmp_num_contours,
     1            h_win_min, h_win_max, v_win_min, v_win_max,
     1            multi_2d(slot)%val_in_oxide,
     1            max_cont_colors, do_yflip, v_sum)
            else
               if( hidden_line_method .eq. HIDDEN_LINES_HISTOGRAM) then
                  call color_contours_histo(multi_2d(slot)%data,
     1               tmp_contours, tmp_num_contours,
     1               multi_2d(slot)%x, multi_2d(slot)%y,
     1               multi_2d(slot)%nx, multi_2d(slot)%ny,
     1               multi_2d(slot)%nx,
     1               h_win_min, h_win_max, v_win_min, v_win_max,
     1               max_cont_colors, do_yflip, v_sum)
               else
                  call color_contours(multi_2d(slot)%data,
     1               tmp_contours, tmp_num_contours,
     1               multi_2d(slot)%x, multi_2d(slot)%y,
     1               multi_2d(slot)%nx, multi_2d(slot)%ny,
     1               multi_2d(slot)%nx,
     1               h_win_min, h_win_max, v_win_min, v_win_max,
     1               max_cont_colors, do_yflip, v_sum)
               endif
            endif
         else
            if( multi_2d(slot)%data_is_trigrid ) then
               call contur_tri(multi_2d(slot)%data,
     1            multi_2d(slot)%cord, multi_2d(slot)%nop,
     1            multi_2d(slot)%tri_ne,
     1            multi_2d(slot)%tri2reg, multi_2d(slot)%imat,
     1            multi_2d(slot)%tri_npts,
     1            tmp_contours, tmp_num_contours,
     1            h_win_min, h_win_max, v_win_min, v_win_max,
     1            multi_2d(slot)%val_in_oxide,
     1            .FALSE., do_yflip, v_sum,
     1            multi_2d(slot)%clabel, multi_2d(slot)%solid,
     1            multi_2d(slot)%colormat)
            else
               call contur(multi_2d(slot)%data,
     1            tmp_contours, tmp_num_contours,
     1            multi_2d(slot)%x, multi_2d(slot)%y,
     1            multi_2d(slot)%nx, multi_2d(slot)%ny,
     1            multi_2d(slot)%nx,
     1            h_win_min, h_win_max, v_win_min, v_win_max,
     1            .FALSE., do_yflip, v_sum,
     1            multi_2d(slot)%clabel, multi_2d(slot)%solid)
            endif
         endif
         endif
      end do
c
c     Draw in junctions, tell contour to use default line type
c     with dash and over-ride set_line_index (a hack for this situation).
      junction(1) = 1.0e5  ! not zero, otherwise we go crazy in the oxide!
      do i=1,n_data_2d
         slot = slot_index_2d(i)
         if( multi_2d(slot).selected .and.
     1       multi_2d(slot)%junctions ) then
            if( multi_2d(slot)%data_is_trigrid ) then
            call contur_tri(multi_2d(slot)%net_doping,
     1            multi_2d(slot)%cord, multi_2d(slot)%nop,
     1            multi_2d(slot)%tri_ne,
     1            multi_2d(slot)%tri2reg, multi_2d(slot)%imat,
     1            multi_2d(slot)%tri_npts, junction, 1,
     1            h_win_min, h_win_max, v_win_min, v_win_max,
     1            .FALSE.,
     1            .TRUE., do_yflip, v_sum, .FALSE.,
     1            multi_2d(slot)%solid, .FALSE.)
            else
            call contur(multi_2d(slot)%net_doping, junction, 1,
     1            multi_2d(slot)%x, multi_2d(slot)%y,
     1            multi_2d(slot)%nx, multi_2d(slot)%ny,
     1            multi_2d(slot)%nx,
     1            h_win_min, h_win_max, v_win_min, v_win_max,
     1            .TRUE., do_yflip, v_sum, .FALSE., 
     1            multi_2d(slot)%solid)
            endif
         endif
      end do
      call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c
c     Draw in boundary and grid
      do i=1,n_data_2d
         slot = slot_index_2d(i)
         if( multi_2d(slot).selected ) then
            call draw_multi_boundaries(i,
     1        h_win_min, h_win_max, v_win_min, v_win_max,
     2        v_sum, do_yflip)
            if( multi_2d(slot).grid ) then
               call plot_multi_grid(i, v_sum, do_yflip)
            endif
         endif
      end do
c
c     draw box around contours
      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(DEFAULT_LINE)
      call move_abs(ZEROO, ZEROO)
      call draw_abs(ONEE,  ZEROO)
      call draw_abs(ONEE,  ONEE)
      call draw_abs(ZEROO, ONEE)
      call draw_abs(ZEROO, ZEROO)
c
c     horizontal label...
c
c     Magic constants are amount window must be expanded on each side
c     to keep in proportion with previous plotting window,
c     but with full viewport in x. They scale with the ratio 13/20 =
c     (0.825 - 0.175)/1.0 for the current plot window setting.
c        The left x window:  13/20 = (0.0 - 0.175)/x  x = -7/26
c        The right x window: 13/20 = (1.0 - 0.825)/x  x =  7/26
      call GKS$SET_WINDOW(1, -0.269, 1.269, 0.0, 0.2)
      call GKS$SET_VIEWPORT(1, 0.0, ratio_x,
     1                         0.0*ratio_y,  0.1*ratio_y)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      call GKS$SELECT_XFORM(1)
c
c     Plot tic marks
c     Find first major tic mark below h_win_min
      if( h_win_min .lt. 0.0 ) then
         h_cur = FLOAT(INT(h_win_min/h_tic_major) - 1)*h_tic_major
      else
         h_cur = FLOAT(INT(h_win_min/h_tic_major))*h_tic_major
      endif
c     Normalize to 0.0 to 1.0 scale
      h_cur = (h_cur - h_win_min)/(h_win_max - h_win_min)
      delta = h_tic_major/
     1        (FLOAT(h_minor_tic_freq)*(h_win_max - h_win_min))
      i = 0
      call tune_charsize(0.035, 0.06, -0.5, 0.025, 0.0458, 0.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                        GKS$K_TEXT_VALIGN_HALF)
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,0.2)
c           Determine where a big tic mark should go
            if( MOD(i,h_minor_tic_freq) .eq. 0 ) then
               call draw_abs(h_cur,0.15)
c              Plot tic label
               ticvalue = pm_scale_cm2plot*
     1             (h_cur*(h_win_max - h_win_min) + h_win_min)
               call format_axis_label(ticvalue, 
     1            pm_scale_cm2plot*h_tic_major, num_label, exp_flag)
               ilen = string_len(num_label)
               call GKS$TEXT(h_cur, 0.12, num_label(1:ilen))
            else
               call draw_abs(h_cur,0.175)
            endif
         endif
         h_cur = h_cur + delta
         i = i + 1
      enddo
c
c     Horizontal axis label
      call tune_charsize(0.035, 0.06 , -0.5, 0.03, 0.05, 0.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                        GKS$K_TEXT_VALIGN_BOTTOM)
      ilen = string_len(h_label)
      ilen2 = string_len(h_units)
      itot = 0
      if( ilen .gt. 0 ) then
         itot = ilen
         work_string = h_label(1:ilen)
      endif
      if( ilen2 .gt. 0 ) then
         itot = ilen + ilen2 + 3
         work_string(ilen+1:itot) = ' ' // left_unit_delim //
     1      h_units(1:ilen2) // right_unit_delim
      endif
      if( itot .gt. 0 ) then
CC         call GKS$TEXT(0.5,0.02, work_string(1:itot) )
         call mathmode(ws_id,0.5,0.02, work_string(1:itot),.FALSE.,lerr)
      endif
c
c     vertical label...
c
c     Magic constants are amount window must be expanded on each side
c     to keep in proportion with previous plotting window,
c     but with full viewport in y. They scale with the ratio 7/10 =
c     (0.8 - 0.1)/1.0 for the current plot window setting.
c        The lower y window:  7/10 = (0.0 - 0.1)/x  x = -1/7
c        The upper y window:  7/10 = (1.0 - 0.8)/x  x =  2/7
      call GKS$SET_WINDOW  (1, 0.0, 0.15, -0.1429, 1.2857)
      call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 0.175*ratio_x,
     1                         0.0*ratio_y, ratio_y)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      call GKS$SELECT_XFORM(1)
c
c     Plot tic marks
c     Find first tic mark
      if( v_win_min .lt. 0.0 ) then
         v_cur = FLOAT(INT(v_win_min/v_tic_major) - 1)*v_tic_major
      else
         v_cur = FLOAT(INT(v_win_min/v_tic_major))*v_tic_major
      endif
c     Normalize to 0.0 to 1.0 scale
      if( do_yflip ) then
         v_cur = 1.0 -
     1           (v_cur - v_win_min)/(v_win_max - v_win_min)
      else
         v_cur = (v_cur - v_win_min)/(v_win_max - v_win_min)
      endif
      delta = v_tic_major/
     1        (FLOAT(v_minor_tic_freq)*(v_win_max - v_win_min))
      i = 0
c      call tune_charsize(0.015, 0.04, -0.2, 0.015, 0.025, 0.0)
      call tune_charsize(0.015, 0.057, -0.2, 0.015, 0.035, 0.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_RIGHT,
     1                        GKS$K_TEXT_VALIGN_HALF)
      if( do_yflip ) then
c        We use -0.0001 instead of 0.0 as a safeguard agains fp rounding
         do while(v_cur .ge. -0.0001)
c           Plot tic marks only if within axis bounds
c           We use 1.0001 instead of 1.0 as a safeguard agains fp rounding
            if( v_cur .le. 1.0001 ) then
               call move_abs(0.15,v_cur)
c              Determine where a big tic mark should go
               if( MOD(i,v_minor_tic_freq) .eq. 0 ) then
                  call draw_abs(0.125,v_cur)
c                 Plot tic label
                  ticvalue = pm_scale_cm2plot*
     1                (v_cur*(v_win_min - v_win_max) + v_win_max)
                  call format_axis_label(ticvalue, 
     1               pm_scale_cm2plot*v_tic_major, num_label, exp_flag)
                  ilen = string_len(num_label)
                  call GKS$TEXT(0.12, v_cur, num_label(1:ilen))
               else
                  call draw_abs(0.1375,v_cur)
               endif
            endif
            v_cur = v_cur - delta
            i = i + 1
         enddo
      else
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(0.15,v_cur)
c              Determine where a big tic mark should go
               if( MOD(i,v_minor_tic_freq) .eq. 0 ) then
                  call draw_abs(0.125,v_cur)
c                 Plot tic label
                  ticvalue = pm_scale_cm2plot*
     1                (v_cur*(v_win_max - v_win_min) + v_win_min)
                  call format_axis_label(ticvalue, 
     1               pm_scale_cm2plot*v_tic_major, num_label, exp_flag)
                  ilen = string_len(num_label)
                  call GKS$TEXT(0.12, v_cur, num_label(1:ilen))
               else
                  call draw_abs(0.1375,v_cur)
               endif
            endif
            v_cur = v_cur + delta
            i = i + 1
         enddo
      endif
c
c     Vertical axis label
c      call tune_charsize(0.024, 0.02, 0.0, 0.03, 0.02, 0.0)
      call tune_charsize(0.034, 0.02, 0.0, 0.043, 0.02, 0.0)
      call GKS$SET_TEXT_UPVEC(-1.0,0.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                        GKS$K_TEXT_VALIGN_HALF)
      ilen = string_len(v_label)
      ilen2 = string_len(v_units)
      itot = 0
      if( ilen .gt. 0 ) then
         itot = ilen
         work_string = v_label(1:ilen)
      endif
      if( ilen2 .gt. 0 ) then
         itot = ilen + ilen2 + 3
         work_string(ilen+1:itot) = ' ' // left_unit_delim //
     1      v_units(1:ilen2) // right_unit_delim
      endif
      if( itot .gt. 0 ) then
CC         call GKS$TEXT(0.03, 0.45, work_string(1:itot) )
         call mathmode(ws_id,0.03,0.5,work_string(1:itot),.FALSE.,lerr)
      endif
      call GKS$SET_TEXT_UPVEC(0.0,1.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                        GKS$K_TEXT_VALIGN_NORMAL)
c
c     titling...
c
      call GKS$SET_WINDOW(1, 0.0,1.0,0.0,.2)
      call GKS$SET_VIEWPORT(1, 0.0*ratio_x, 1.0*ratio_x,
     1                         0.80*ratio_y, 1.0*ratio_y)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      call GKS$SELECT_XFORM(1)
c     Gag: yet another set of magic constants :-) jvf
      call tune_charsize(0.025, 0.04, -0.5, 0.02, 0.025, 0.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                        GKS$K_TEXT_VALIGN_NORMAL)
      v_cur = 0.15
      do i=1,MIN(3,num_titles)
         ilen = string_len(title_dat(i))
         if( ilen .gt. 0 ) then
            call mathmode(ws_id,0.0,v_cur,title_dat(i),.FALSE.,lerr)
         endif
         v_cur = v_cur - 0.05
      end do
c
c     draw legend on side...
      if( do_uniform_contours ) then
         tmp_num_contours = num_contours
         do i=1,tmp_num_contours
            tmp_contours(i) = contour_values(i)
         end do
      else
         slot = slot_index_2d(probe_dataset)
         tmp_num_contours = multi_2d(slot)%num_contours
         do i=1,tmp_num_contours
            tmp_contours(i) = multi_2d(slot)%contours(i)
         end do
      end if
      call  plot_color_map_legend(tmp_contours, tmp_num_contours,
     1      rmax, multi_2d(slot)%color .and. can_do_color_contours,
     1      multi_2d(slot)%solid)
c
      call deactivate_plot_device
c
c     Plot annotations from the save area, if specified
c
      if( plot_device .eq. hardcopy_device ) then
         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
      end
c
      subroutine pm_minmax_multi_2d(hmin, hmax, vmin, vmax,
     1                              entire_matrix,
     2                              overall_min, overall_max)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'compare_2d.cmn'
c
      real hmin, hmax, vmin, vmax, overall_min, overall_max
      logical entire_matrix
c
c     Locals
c
      real temp_min, temp_max
      integer i, slot
c
      overall_min =  PM_MAX_FLOAT
      overall_max = -PM_MAX_FLOAT
c
      do i=1,n_data_2d
         slot = slot_index_2d(i)
         if( multi_2d(slot)%data_is_trigrid ) then
            call minmax_triarray(multi_2d(slot)%data,
     1                     multi_2d(slot)%tri_npts,
     1                     multi_2d(slot)%tri_npts,
     1                     multi_2d(slot)%cord,
     2                     hmin, hmax, vmin, vmax, entire_matrix,
     4                     temp_min, temp_max)
         else
            call minmax_rectarray(multi_2d(slot)%data,
     1                     multi_2d(slot)%nx, multi_2d(slot)%ny,
     1                     multi_2d(slot)%nx,
     1                     multi_2d(slot)%x, multi_2d(slot)%y,
     2                     hmin, hmax, vmin, vmax, entire_matrix,
     4                     temp_min, temp_max)
         endif
         overall_min = MIN(overall_min, temp_min)
         overall_max = MAX(overall_max, temp_max)
      end do
c
      end
c
      subroutine pm_get_2d_coord_limits(entry, hmin, hmax, vmin, vmax)
      implicit none
c
      integer entry
      real hmin, hmax, vmin, vmax
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'compare_2d.cmn'
c
      integer slot, i
c
      hmin = 0.0
      hmax = 0.0
      vmin = 0.0
      vmax = 0.0
c
      slot = slot_index_2d(entry)
      if( slot .eq. -1 ) then
         write(luttyo,*) '(pm_get_2d_coord_limits) Bad entry passed in'
         return
      endif
c
      if( multi_2d(slot).data_is_trigrid ) then
         hmin =   PM_MAX_FLOAT
         hmax = - PM_MAX_FLOAT
         vmin =   PM_MAX_FLOAT
         vmax = - PM_MAX_FLOAT
         do i=1,multi_2d(slot).tri_npts
            hmin = MIN(hmin,multi_2d(slot).cord(1,i))
            hmax = MAX(hmax,multi_2d(slot).cord(1,i))
            vmin = MIN(vmin,multi_2d(slot).cord(2,i))
            vmax = MAX(vmax,multi_2d(slot).cord(2,i))
         end do
      else
         hmin = multi_2d(slot).x(1)
         hmax = multi_2d(slot).x(multi_2d(slot).nx)
         vmin = multi_2d(slot).y(1)
         vmax = multi_2d(slot).y(multi_2d(slot).ny)
      endif
c
      return
      end
c
      subroutine find_max_in_multi_subset(entry,
     1                                    hmin, hmax, vmin, vmax,
     1                                    do_log, rmax)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'compare_2d.cmn'
c

      logical do_log
      integer entry
      real hmin, hmax, vmin, vmax, rmax
c
c     locals
      integer slot, i, j, ih1, ih2, iv1, iv2, off
c
      slot = slot_index_2d(entry)
      if( slot .eq. -1 ) return
c
      if( multi_2d(slot)%data_is_trigrid ) then
         rmax = -PM_MAX_FLOAT
         do 200 i=1,multi_2d(slot)%tri_npts
c           Is point in subset region?
            if( multi_2d(slot)%cord(1,i) .lt. hmin .or.
     1          multi_2d(slot)%cord(1,i) .gt. hmax .or.
     1          multi_2d(slot)%cord(2,i) .lt. vmin .or.
     1          multi_2d(slot)%cord(2,i) .gt. vmax) goto 200
            rmax = MAX(rmax,multi_2d(slot)%data(i))
200      continue
      else
         call find_greatest_subset_of_array(
     1           multi_2d(slot)%x, multi_2d(slot)%y,
     1           hmin, hmax, vmin, vmax,
     1           multi_2d(slot)%nx, multi_2d(slot)%ny,
     1           ih1, ih2, iv1, iv2)
c
c        find max of data
         rmax = -PM_MAX_FLOAT
         do j=iv1,iv2
            off = (j-1)*multi_2d(slot)%nx
            do i=ih1,ih2
               if( rmax .lt. multi_2d(slot)%data(i+off) ) then
                  rmax = multi_2d(slot)%data(i+off)
               endif
            end do
         end do
      endif
c
      if( do_log ) then
         rmax = SIGN(10.0**MIN(38.0,ABS(rmax)),rmax)
      endif
      return
      end
c
      subroutine plot_multi_grid(entry, v_sum, do_yflip)
      implicit none
c
c     Plot grid lines over contours
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'compare_2d.cmn'
c
      integer entry
      logical do_yflip
      real v_sum
c
      integer i, j, slot, nx, ny
c
      slot = slot_index_2d(entry)
      if( slot .eq. -1 ) return
c
      call set_line_index(DEFAULT_LINE)
c
      if( multi_2d(slot)%datatype .eq. DT_VLSICAP ) then
c        VLSICAP mesh is not supported. 7-Apr-99 jvf
      else if( multi_2d(slot)%datatype .eq. DT_SUPREM4 .or.
     1    multi_2d(slot)%datatype .eq. DT_PISCES  .or.
     1    multi_2d(slot)%datatype .eq. DT_TIF     .or.
     1    multi_2d(slot)%datatype .eq. DT_MCP2D   .or.
     1    multi_2d(slot)%datatype .eq. DT_MEDICI  .or.
     1    (multi_2d(slot)%datatype .eq. DT_ASCII  .and.
     1     multi_2d(slot)%data_is_trigrid) .or.
     1    (multi_2d(slot)%datatype .eq. DT_PROMIS .and.
     1     multi_2d(slot)%data_is_trigrid)
     1   ) then
         call plot_tri_mesh(multi_2d(slot)%cord, multi_2d(slot)%nop,
     1            multi_2d(slot)%tri_npts, multi_2d(slot)%tri_ne,
     1            v_sum, do_yflip)
      else
         nx = multi_2d(slot)%nx
         ny = multi_2d(slot)%ny
         if( do_yflip ) then
            do j=1,ny
               call move_abs(multi_2d(slot)%x(1),
     1                      (v_sum - multi_2d(slot)%y(j)))
               call draw_abs(multi_2d(slot)%x(nx),
     1                      (v_sum - multi_2d(slot)%y(j)))
            end do
            do i=1,multi_2d(slot)%nx
               call move_abs(multi_2d(slot)%x(i),
     1                      (v_sum - multi_2d(slot)%y(1)))
               call draw_abs(multi_2d(slot)%x(i),
     1                      (v_sum - multi_2d(slot)%y(ny)))
            end do
         else
            do j=1,ny
               call move_abs(multi_2d(slot)%x(1), multi_2d(slot)%y(j))
               call draw_abs(multi_2d(slot)%x(nx),multi_2d(slot)%y(j))
            end do
            do i=1,nx
               call move_abs(multi_2d(slot)%x(i),multi_2d(slot)%y(1))
               call draw_abs(multi_2d(slot)%x(i),multi_2d(slot)%y(ny))
            end do
         endif
      endif
c
      return
      end
c
      subroutine draw_multi_boundaries(entry, hmin, hmax, vmin, vmax,
     1                                 v_sum, do_yflip)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'compare_2d.cmn'
c
      integer entry
      real hmin, hmax, vmin, vmax, v_sum
      logical do_yflip
c
c     Locals
      integer slot, n
      real x(5), y(5)
c
      slot = slot_index_2d(entry)
      if( slot .eq. -1 ) return
c
      if( multi_2d(slot)%datatype .eq. DT_MINIMOS ) then
         call draw_multi_bound_mini(entry, v_sum)
      else if( multi_2d(slot)%datatype .eq. DT_VLSICAP ) then
c        VLSICAP boundary is not supported. 7-Apr-99 jvf
      else if( multi_2d(slot)%datatype .eq. DT_PROMIS ) then
c        PROMIS boundary is not supported. 7-Apr-99 jvf
      else if( multi_2d(slot)%datatype .eq. DT_2DOP ) then
c        2DOP boundary minimal support
         call set_line_index(DEFAULT_LINE)
c
         x(1) = multi_2d(slot)%x(1)
         x(2) = multi_2d(slot)%x(multi_2d(slot)%nx)
         x(3) = x(2)
         x(4) = x(1)
         x(5) = x(1)
c
         y(1) = multi_2d(slot)%y(1)
         y(2) = y(1)
         y(3) = multi_2d(slot)%y(multi_2d(slot)%ny)
         y(4) = y(3)
         y(5) = y(1)
         call GKS$POLYLINE(n, x, y)
      else if( multi_2d(slot)%datatype .eq. DT_PISCES .or.
     1         multi_2d(slot)%datatype .eq. DT_SUPREM4 .or.
     2         multi_2d(slot)%datatype .eq. DT_TIF .or.
     3         multi_2d(slot)%datatype .eq. DT_MEDICI .or.
     5        (multi_2d(slot)%datatype .eq. DT_ASCII  .and.
     6         multi_2d(slot)%data_is_trigrid)) then
         call draw_2d_bound_tri(multi_2d(slot)%cord,
     1                          multi_2d(slot)%nop,
     1                          multi_2d(slot)%imat,
     1                          multi_2d(slot)%tri_nextel,
     1                          multi_2d(slot)%tri2reg,
     1                          multi_2d(slot)%tri_npts,
     1                          multi_2d(slot)%tri_ne,
     1                          multi_2d(slot)%tri_nreg,
     1                          v_sum, hmin, hmax, vmin, vmax,
     1                          do_yflip)
      else if( multi_2d(slot)%datatype .eq. DT_ASCII .and.
     1         .NOT. multi_2d(slot)%data_is_trigrid ) then
         call set_line_index(DEFAULT_LINE)
c
         n = 5
         x(1) = multi_2d(slot)%x(1)
         x(2) = multi_2d(slot)%x(multi_2d(slot)%nx)
         x(3) = x(2)
         x(4) = x(1)
         x(5) = x(1)
c
         y(1) = multi_2d(slot)%y(1)
         y(2) = y(1)
         y(3) = multi_2d(slot)%y(multi_2d(slot)%ny)
         y(4) = y(3)
         y(5) = y(1)
         call GKS$POLYLINE(n, x, y)
      endif
c
      return
      end
c
      subroutine draw_multi_bound_mini(entry, v_sum)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'defaults.cmn'
      include 'compare_2d.cmn'
      include 'units.cmn'
c
      real v_sum
      integer entry
c
c     Locals
      integer slot
      integer nx, ny, nxsp, nxsoxm, nysem, nyoxm, nytop, nytot
      integer string_len
      logical nonplanar_boundary, three_d
      real    fox_temp(2*ARRAYSIZE), dist_temp(2*ARRAYSIZE), xl, tox,
     1        xzero, yzero
c
      slot = slot_index_2d(entry)
      if( slot .eq. -1 ) return
c
      nx      = multi_2d(slot)%nx
      ny      = multi_2d(slot)%ny
      xl      = multi_2d(slot)%misc_real(8)
      tox     = multi_2d(slot)%misc_real(9)
      nxsp    = multi_2d(slot)%misc_int(1)
      nxsoxm  = multi_2d(slot)%misc_int(2)
      ny      = multi_2d(slot)%misc_int(3)
      nyoxm   = multi_2d(slot)%misc_int(4)
      nytop   = multi_2d(slot)%misc_int(5)
      nonplanar_boundary = multi_2d(slot)%misc_int(6)
      three_d = multi_2d(slot)%misc_int(7)
c
      if( string_len(multi_2d(slot)%x_expression) .ne. 0 ) then
c        do something here eventually...
      else
         xl = multi_2d(slot).mx*xl + pm_scale_plot2cm*multi_2d(slot).ax
         xzero = pm_scale_plot2cm*multi_2d(slot).ax
      endif
c
      nysem   = nyoxm + ny
      nytot   = multi_2d(slot).ny
c
      call set_line_index(DEFAULT_LINE)
      if( nonplanar_boundary .or. three_d ) then
c        We don't do anything (yet) for 3D MINIMOS
c        We'd need to save foxo,foxu arrays
      else
c        Standard rectangular boundary
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c        White out oxide region if there is no data there
         if( .NOT. multi_2d(slot)%val_in_oxide ) then
            call GKS$SET_FILL_COLOR_INDEX(0)
            fox_temp(1)  = multi_2d(slot)%x(1)
            dist_temp(1) = v_sum
            fox_temp(2)  = multi_2d(slot)%x(nx)
            dist_temp(2) = v_sum
            fox_temp(3)  = multi_2d(slot)%x(nx)
            dist_temp(3) = v_sum - multi_2d(slot)%y(1)
            fox_temp(4)  = multi_2d(slot)%x(1)
            dist_temp(4) = v_sum - multi_2d(slot)%y(1)
            if( l_oxide_fill )
     1         call gks$fill_area(4,fox_temp,dist_temp)
c           If SOI, then fill in oxide film
            if( nysem .lt. nytot ) then
               fox_temp(1)  = multi_2d(slot)%x(1)
               dist_temp(1) = v_sum - multi_2d(slot)%y(nysem)
               fox_temp(2)  = multi_2d(slot)%x(nx)
               dist_temp(2) = v_sum - multi_2d(slot)%y(nysem)
               fox_temp(3)  = multi_2d(slot)%x(nx)
               dist_temp(3) = v_sum - multi_2d(slot)%y(nytot)
               fox_temp(4)  = multi_2d(slot)%x(1)
               dist_temp(4) = v_sum - multi_2d(slot)%y(nytot)
               if( l_oxide_fill )
     1            call gks$fill_area(4,fox_temp,dist_temp)
            endif
            call GKS$SET_FILL_COLOR_INDEX(1)
         endif
c
c        Draw line for oxide
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c        Flip y coordinate with v_sum
         call move_abs(multi_2d(slot)%x(1),  v_sum)
         call draw_abs(multi_2d(slot)%x(nx), v_sum)
c        Draw line for SOI interface
         if( nysem .lt. nytot ) then
            call move_abs(multi_2d(slot)%x(1),
     1                    v_sum - multi_2d(slot)%y(nysem))
            call draw_abs(multi_2d(slot)%x(nx),
     1                    v_sum - multi_2d(slot)%y(nysem))
         endif
c        Draw box for gate region, S/D contacts
c        Source side is at 0.0, drain side is at xl
c        Remember to flip y, too...
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
         call move_abs(multi_2d(slot)%x(nxsp-1), v_sum)
         call draw_abs(multi_2d(slot)%x(nxsp-1),
     1                 (v_sum - multi_2d(slot)%y(1)))
c
         call move_abs(xzero,(v_sum - multi_2d(slot)%y(nytop)))
         call draw_abs(xzero,(v_sum + tox))
         call draw_abs(xl, (v_sum + tox))
         call draw_abs(xl, (v_sum - multi_2d(slot)%y(nytop)))
         call draw_abs(xzero,(v_sum - multi_2d(slot)%y(nytop)))
c
         call move_abs(multi_2d(slot)%x(nxsoxm+1),  v_sum)
         call draw_abs(multi_2d(slot)%x(nxsoxm+1),
     1                 (v_sum - multi_2d(slot)%y(1)))
      endif
c
      return
      end
