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 contour plot utilities
c
      subroutine plot2d
      implicit none
c
c     Asks user what array to plot, then plot it!
c
c     Programmer: John Faricelli
c     Date: 18-Jun-85
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'plot.cmn'
      include 'io.cmn'
      include 'bincodes.cmn'
c
c     locals
c
      integer answer, string_len
      character*10 buffer
      character*80 errfmt
c
      errfmt = '('' '')'
c
c     Short circuit this menu if there is only one data item to select from
c
      if( max_menu_num .eq. 1 ) then
         if( is_data_1d(1) ) then
            call clear_screen
            call beep_terminal
            write(luttyo,100)
100         format(' ','Cannot plot contours of 1-D data'/
     1             ' ','Hit return to continue...')
            read(luttyi,106,end=300,err=300) buffer
         else
            call ship_contours(arrdat(arrind(menu_to_arr(1))),
     1                   ni,nj,
     1                   data_name(1), data_units(1),
     2                   data_nickname(1),
     3                   is_data_log(1), is_data_in_oxide(1),
     4                   .FALSE.)
         endif
         return
      endif
c
200   call clear_screen
      write(luttyo,errfmt)     ! print any error messages
      errfmt = '('' '')'  ! clear them
      call write_data_menu('Contour plot section:')
      write(luttyo,105)
105   format(/' ','Request? ',$)
      read(luttyi,106,err=300,end=300) buffer
106   format(a)
      call upcase(buffer)
      if( string_len(buffer) .eq. 0) then
         answer = EXIT_MENU
      else if( buffer(1:1) .eq. 'Q' .or.  buffer(1:1) .eq. 'E' ) then
         answer = EXIT_MENU
      else
         read(buffer,*,err=300) answer
      endif
c
c    Error check:
c       1) Exit menu
c       2) Answer in range
c       3) Is requested data in binary file?
c       4) Is data only 1D?
c
      if( answer .eq. EXIT_MENU ) then
         goto 1000
      else if ( answer .lt. 1 .or. answer .gt. max_menu_num ) then
         errfmt = '('' '',''Invalid input: Please try again...'')'
         call beep_terminal
      else if( .NOT. in_binary_file(answer) ) then
         errfmt = '('' '',''Data not available from input file'')'
         call beep_terminal
      else if( is_data_1d(answer) ) then
         errfmt = '('' '',''Cannot plot contours of 1D data'')'
         call beep_terminal
      else
         call ship_contours(arrdat(arrind(menu_to_arr(answer))),
     1                   ni,nj,
     1                   data_name(answer), data_units(answer),
     2                   data_nickname(answer),
     3                   is_data_log(answer), is_data_in_oxide(answer),
     4                   .FALSE.)
      endif
      goto 200
c
c     illegal input
300   errfmt = '('' '',''Invalid input: Please try again...'')'
      call beep_terminal
      goto 200
c
1000  continue
      return
      end
c
      subroutine ship_contours(data_array, isize, jsize, label,
     1                         units, nickname, log_data, val_in_oxide,
     2                         lsave)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
      include 'units.cmn'
      include 'defaults.cmn'
c
c     To get plot title
      include 'callable.cmn'
c
      integer isize, jsize
      logical log_data, val_in_oxide, lsave
      real data_array(isize, jsize)
      character*(*) label, units, nickname
c
c     locals
c
      logical do_plot_title, do_junctions
      logical contour_values_set, do_grid,
     1        fake_ln03, do_color, do_clabel, do_negate, do_yflip,
     2        do_solid, do_log_interp, do_draw_to_scale,
     3        lentire_matrix, do_colormat, do_log
      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
      integer item_number, ierror
      character*131 buffer
      character*80  string
      character*3 str_color, str_negate, str_label, str_grid, str_junct,
     1            str_yflip, str_solid, str_aspect, str_colormat,
     1            str_log
      character*6 str_log_interp
      logical eof, abort, is_graphics_terminal, lopen,
     1        delayed_message, 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_log = log_data
      if( binfiletype .eq. DT_FAIM .or.
     1    binfiletype .eq. DT_MCP2D ) then
         do_draw_to_scale = .TRUE.
      else
         do_draw_to_scale = .FALSE.
      endif
      cdum = ' '
c
      ulen = string_len(pm_length_unit)
c
      if( .not. lsave ) then
         plot_title = ' '
      endif
      if( callable_num_titles .eq. 1 ) then
         plot_title = callable_titles(1)
      endif
c
c     Find min/max of entire dataset
      lentire_matrix = .TRUE.
      if( data_is_trigrid ) then
         call minmax_triarray(data_array, tri_npts, isize*jsize,
     1                     cord,
     2                     rdum1, rdum2, rdum3, rdum4, lentire_matrix,
     3                     min_value, max_value)
      else
         call minmax_rectarray(data_array, isize, jsize, isize,
     1                     xyz(1,idisti), xyz(1,idistj),
     2                     rdum1, rdum2, rdum3, rdum4, lentire_matrix,
     3                     min_value, max_value)
       endif
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 = xyzlimits(1, idisti)
         twod_h_max = xyzlimits(2, idisti)
         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
         twod_v_min = xyzlimits(1, idistj)
         twod_v_max = xyzlimits(2, idistj)
         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( binfiletype .eq. DT_ASCII ) then
            horiz_label = 'X'
            horiz_units = ' '
         else if( binfiletype .eq. DT_ANALYTICAL ) then
            horiz_label = 'X'
            horiz_units = ' '
         else 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( binfiletype .eq. DT_ASCII .or.
     1       binfiletype .eq. DT_ANALYTICAL ) then
            vert_label = 'Y'
            vert_units = ' '
         else 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
c
c     Initialization for menu
      hmin_full = pm_scale_cm2plot*xyzlimits(1, idisti)
      hmax_full = pm_scale_cm2plot*xyzlimits(2, idisti)
      vmin_full = pm_scale_cm2plot*xyzlimits(1, idistj)
      vmax_full = pm_scale_cm2plot*xyzlimits(2, idistj)
      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.
      if( data_is_trigrid ) then
         call minmax_triarray(data_array, tri_npts, isize*jsize,
     1                     cord,
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
      else
         call minmax_rectarray(data_array, isize, jsize, isize,
     1                     xyz(1,idisti), xyz(1,idistj),
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
       endif
c
      if( lsave ) then
         do_color  = twod_color
         do_clabel = twod_clabel
         do_negate = twod_negate
         do_grid   = twod_grid
         do_junctions = twod_junctions
         do_yflip  = .NOT. twod_invert
         do_solid  = twod_solid
         do_log_interp = twod_log_interp
         do 110 i=1,twod_num_contours
            contour_values(i) = twod_contours(i)
110      continue
         do_colormat = twod_colormat
      else
         do_negate = .FALSE.
         twod_num_contours = MAX_CONTOURS
         if( data_is_trigrid ) then
            call get_default_contours_tri(twod_h_min, twod_h_max,
     1              twod_v_min, twod_v_max,
     2              data_array, tri_npts, do_log, do_negate,
     3              twod_num_contours, contour_values)
         else
            call get_default_contours(twod_h_min, twod_h_max,
     1              twod_v_min, twod_v_max,
     2              data_array, ni, nj, do_log, do_negate,
     3              twod_num_contours, contour_values)
         endif
c        Sort in ascending order (MUST BE DONE!!!)
         call sort_contours(contour_values,twod_num_contours)
         do_color  = .TRUE.
         do_clabel = .FALSE.
         do_grid   = .FALSE.
         if( binfiletype .eq. DT_FAIM  .or.
     1       binfiletype .eq. DT_MCP2D .or.
     1       binfiletype .eq. DT_ASCII .or.
     1       binfiletype .eq. DT_ANALYTICAL ) then
            do_yflip  = .FALSE.
         else
            do_yflip  = .TRUE.
         endif
         do_solid  = .FALSE.
         do_log_interp = do_log
         do_colormat = .FALSE.
         if( net_dop_array_index .ne. 0 ) then
            do_junctions  = .TRUE.
         else
            do_junctions  = .FALSE.
         endif
      endif
c
c     Sort in ascending order (MUST BE DONE!!!)
      call sort_contours(contour_values,twod_num_contours)
      do_plot_title = .TRUE.
c
      if( do_color ) then
         str_color = 'yes'
      else
         str_color = ' no'
      endif
      if( do_clabel ) then
         str_label = 'yes'
      else
         if( do_color ) then
            str_label = ' NA'
         else
            str_label = ' no'
         endif
      endif
      if( do_negate ) then
         str_negate = 'yes'
      else
         str_negate = ' no'
      endif
      if( do_grid ) then
         str_grid = 'yes'
      else
         str_grid = ' no'
      endif
      if( do_junctions ) then
         str_junct = 'yes'
      else
         str_junct = ' NA'
      endif
      if( do_yflip ) then
         str_yflip = 'yes'
      else
         str_yflip = ' no'
      endif
      if( do_draw_to_scale ) then
         str_aspect = 'yes'
      else
         str_aspect = ' no'
      endif
      if( do_solid ) then
         str_solid = 'yes'
      else
         if( do_color ) then
            str_solid = ' NA'
         else
            str_solid = ' no'
         endif
      endif
      if( do_log ) then
         str_log = 'yes'
      else
         str_log = ' no'
      endif
      if ( do_log_interp ) then
         str_log_interp = '   log'
      else
         str_log_interp = 'linear'
      endif
      if ( do_color ) then
         str_colormat = ' NA'
      else if( do_colormat ) then
         str_colormat = 'yes'
      else
         str_colormat = ' no'
      endif
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_2D_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_2d_to_save(contour_values, twod_num_contours,
     1     do_junctions, do_grid, do_color, do_negate, do_clabel,
     2     .NOT. do_yflip, do_solid, do_log_interp, do_colormat,
     3     nickname, label, units)
c
         call negate_array(work, data_array, ni, nj,
     1      tri_npts, data_is_trigrid, do_negate)
         if( do_log_interp ) then
            call log10_array(work, ni, nj, tri_npts, data_is_trigrid)
         endif
c
         call plot_contours(work,
     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     xyz(1,idisti), xyz(1,idistj), ni, nj,
     8     twod_num_contours, contour_values, label, units,
     9     do_plot_title, do_junctions, do_grid, do_color, do_clabel,
     A     do_yflip, do_solid, do_log_interp, do_colormat,
     B     val_in_oxide, 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 ,51,'CONTOUR PARAMETERS:',19,0)
      call display_string(4 ,42,'Quantity: ',10,0)
      call display_string(5 ,42,'Data range:            to',25,0)
      call display_string(6 ,42,'Plot range:            to',25,0)
      call display_string(8 ,42,'12) Contours',12,0)
      call display_string(8 ,56,'13) Color contours:',19,0)
      call display_string(9 ,56,'14) Label contours:',19,0)
      call display_string(10,56,'15) Junctions:',14,0)
      call display_string(11,56,'16) Grid:',9,0)
      call display_string(12,56,'17) Change sign:',16,0)
      call display_string(13,56,'18) Recompute contours',22,0)
      call display_string(14,56,'19) Invert plot:',16,0)
      call display_string(15,56,'20) Use black lines',20,0)
      call display_string(16,56,'21) Interpolation:',18,0)
      call display_string(17,56,'22) Draw to scale:',18,0)
      call display_string(18,56,'23) Color materials:',20,0)
      call display_string(19,56,'24) 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, plot_title, 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  (5 , 54, min_value, 1)
      call display_real  (5 , 68, max_value, 1)
      if( plot_min_value .gt. plot_max_value ) then
         call display_string(6 ,54,
     1                       'No grid points found   ',23,1)
         else
            call display_string(6 ,53,'            to  ',26,0)
            call display_real  (6 ,54, plot_min_value, 1)
            call display_real  (6 ,68, plot_max_value, 1)
      endif
c     Contours
      call display_string(4 ,52,label,29,0)
      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     Color contours
      call display_string(8 ,77,str_color, 3, 1)
c     Label contours
      call display_string(9 ,77,str_label, 3, 1)
c     Junctions
      call display_string(10,77,str_junct, 3, 1)
c     Grid
      call display_string(11,77,str_grid, 3, 1)
c     Change sign
      call display_string(12,77,str_negate, 3, 1)
c     Invert vertical axis
      call display_string(14,77,str_yflip, 3, 1)
c     Solid lines
      call display_string(15,77,str_solid, 3, 1)
c     Interpolation mode
      call display_string(16,75,str_log_interp, 6, 1)
c     Maintain aspect ratio flag
      call display_string(17,77,str_aspect, 3,1)
c     Color fill each material
      call display_string(18,77,str_colormat, 3,1)
c     Use log contours
      call display_string(19,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_2D_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_2d_to_save(contour_values, twod_num_contours,
     1     do_junctions, do_grid, do_color, do_negate, do_clabel,
     2     .NOT. do_yflip, do_solid, do_log_interp, do_colormat,
     3     nickname, label, units)
c
         call negate_array(work, data_array, ni, nj,
     1      tri_npts, data_is_trigrid, do_negate)
         if( do_log_interp ) then
            call log10_array(work, ni, nj, tri_npts, data_is_trigrid)
         endif
c
         call plot_contours(work,
     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     xyz(1,idisti), xyz(1,idistj), ni, nj,
     8     twod_num_contours, contour_values, label, units,
     9     do_plot_title, do_junctions, do_grid, do_color, do_clabel,
     A     do_yflip, do_solid, do_log_interp, do_colormat,
     B     val_in_oxide, 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 ) then
            if( is_graphics_terminal(current_window_id) ) 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_2D_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_2d_to_save(contour_values,twod_num_contours,
     1     do_junctions, do_grid, do_color, do_negate, do_clabel,
     2     .NOT. do_yflip, do_solid, do_log_interp, do_colormat,
     3     nickname, label, units)
c
                  call negate_array(work, data_array, ni, nj,
     1               tri_npts, data_is_trigrid, do_negate)
                  if( do_log_interp ) then
                     call log10_array(work, ni, nj,
     1                  tri_npts, data_is_trigrid)
                  endif
c
                  call plot_contours(work,
     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     xyz(1,idisti), xyz(1,idistj), ni, nj,
     8     twod_num_contours, contour_values, label, units,
     9     do_plot_title, do_junctions, do_grid, do_color, do_clabel,
     A     do_yflip, do_solid, do_log_interp, do_colormat,
     B     val_in_oxide, lsave, lopen, lannotate, ierror)
                  if( ierror .eq. PLT$NORMAL ) then
                     lopen = .TRUE.
                  endif
c                 Update dataset plot min/max limits
                  lentire_matrix = .FALSE.
                  if( data_is_trigrid ) then
                     call minmax_triarray(data_array, tri_npts,
     1                     isize*jsize, cord,
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
                  else
                     call minmax_rectarray(data_array, isize, jsize,
     1                     isize, xyz(1,idisti), xyz(1,idistj),
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
                  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_2D_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_2d_to_save(contour_values, twod_num_contours,
     1     do_junctions, do_grid, do_color, do_negate, do_clabel,
     2     .NOT. do_yflip, do_solid, do_log_interp, do_colormat,
     3     nickname, label, units)
c
         call negate_array(work, data_array, ni, nj,
     1               tri_npts, data_is_trigrid, do_negate)
         if( do_log_interp ) then
            call log10_array(work, ni, nj,
     1      tri_npts, data_is_trigrid)
         endif
c
         call plot_contours(work,
     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     xyz(1,idisti), xyz(1,idistj), ni, nj,
     8     twod_num_contours, contour_values, label, units,
     9     do_plot_title, do_junctions, do_grid, do_color, do_clabel,
     A     do_yflip, do_solid, do_log_interp, do_colormat,
     B     val_in_oxide, 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              Copy data values back to work, as they may be destroyed
c              by contour plot subroutine (if junctions are plotted)
               call negate_array(work, data_array, ni, nj,
     1                     tri_npts, data_is_trigrid, do_negate)
               if( do_log_interp ) then
                  call log10_array(work, ni, nj,
     1            tri_npts, data_is_trigrid)
               endif
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_interpolate( work, x1, y1,
     1            do_log_interp, data_is_trigrid, 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(plot_title, 
     1              'Enter title (may use mathmode)')
            call display_string(1 , 11, plot_title, 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.
               if( data_is_trigrid ) then
                  call minmax_triarray(data_array, tri_npts,
     1                     isize*jsize, cord,
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
               else
                  call minmax_rectarray(data_array, isize, jsize,
     1                     isize, xyz(1,idisti), xyz(1,idistj),
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
               endif
               if( plot_min_value .gt. plot_max_value ) then
                  call display_string(6 ,54,
     1               'No grid points found   ',23,1)
               else
                  call display_string(6 ,53,'            to  ',26,0)
                  call display_real  (6 ,54, plot_min_value, 1)
                  call display_real  (6 ,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.
               if( data_is_trigrid ) then
                  call minmax_triarray(data_array, tri_npts,
     1                     isize*jsize, cord,
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
               else
                  call minmax_rectarray(data_array, isize, jsize,
     1                     isize, xyz(1,idisti), xyz(1,idistj),
     2                     twod_h_min, twod_h_max,
     3                     twod_v_min, twod_v_max, lentire_matrix,
     4                     plot_min_value, plot_max_value)
               endif
               if( plot_min_value .gt. plot_max_value ) then
                  call display_string(6 ,54,
     1               'No grid points found   ',23,1)
               else
                  call display_string(6 ,53,'            to  ',26,0)
                  call display_real  (6 ,54, plot_min_value, 1)
                  call display_real  (6 ,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
            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 520 i=1,twod_num_contours
                     contour_values(i) = real_vec(i)
520               continue
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
               endif
            endif
         elseif ( item_number .eq. 13 ) then
c           Color contours
            do_color = .NOT. do_color
            if( do_color ) then
               str_color = 'yes'
               str_label = ' NA'
               do_clabel = .FALSE.
               do_solid  = .FALSE.
               str_solid = ' NA'
               str_colormat = ' NA'
            else
               str_color = ' no'
               str_label = ' no'
               str_solid = ' no'
               str_colormat = ' no'
               do_clabel = .FALSE.
               do_solid  = .FALSE.
               do_colormat = .FALSE.
            endif
            call display_string(8 ,77,str_color, 3, 1)
            call display_string(9 ,77,str_label, 3, 1)
            call display_string(15,77,str_solid, 3, 1)
            call display_string(18,77,str_colormat, 3, 1)
         elseif ( item_number .eq. 14 ) then
c           Labels
            if( do_color ) then
               call display_message(
     1            'Not applicable when using color contours')
            else
               do_clabel = .NOT. do_clabel
               if( do_clabel ) then
                  str_label = 'yes'
               else
                  str_label = ' no'
               endif
            endif
            call display_string(9 ,77,str_label, 3, 1)
         elseif ( item_number .eq. 15 ) then
c           Junctions
            if( net_dop_array_index .eq. 0 ) then
               call display_message('Not applicable, no doping')
            else
               do_junctions = .NOT. do_junctions
               if( do_junctions ) then
                  str_junct = 'yes'
               else
                  str_junct = ' no'
               endif
               call display_string(10,77,str_junct, 3, 1)
            endif
         elseif ( item_number .eq. 16 ) then
c           Grid
            do_grid = .NOT. do_grid
            if( do_grid ) then
               str_grid = 'yes'
            else
               str_grid = ' no'
            endif
            call display_string(11,77,str_grid, 3, 1)
         elseif ( item_number .eq. 17 ) then
c           Change sign
            do_negate = .NOT. do_negate
            if( do_negate ) then
               str_negate = 'yes'
            else
               str_negate = ' no'
            endif
            do i=1,twod_num_contours
               contour_values(i) = - contour_values(i)
            end do
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
            call display_string(12,77,str_negate, 3, 1)
         elseif ( item_number .eq. 18 .or. item_number .eq. 24 ) then
            if( item_number .eq. 24 ) then
c              Log contours
               do_log = .NOT. do_log
               if ( do_log ) then
                  str_log = 'yes'
               else
                  str_log = ' no'
               endif
               call display_string(19,77,str_log, 3, 1)
            endif
c           Recompute contours on current bounds
            twod_num_contours = 9
            if( data_is_trigrid ) then
            call get_default_contours_tri(twod_h_min, twod_h_max,
     1           twod_v_min, twod_v_max,
     2           data_array, tri_npts, do_log, do_negate,
     3           twod_num_contours, contour_values)
            else
            call get_default_contours(twod_h_min, twod_h_max,
     1           twod_v_min, twod_v_max,
     2           data_array, ni, nj, do_log, do_negate,
     3           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
         elseif ( item_number .eq. 19 ) then
            do_yflip = .NOT. do_yflip
c           Remember: do_yflip means plot 'normally' for device and
c           process simulators (positive Y down)
            if( do_yflip ) then
              str_yflip = 'yes'
            else
              str_yflip = ' no'
            endif
            call display_string(14,77,str_yflip, 3, 1)
         elseif ( item_number .eq. 20 ) then
c           Solid lines
            if( do_color ) then
               call display_message(
     1            'Not applicable when using color contours')
            else
               do_solid = .NOT. do_solid
               if( do_solid ) then
                  str_solid = 'yes'
               else
                  str_solid = ' no'
               endif
            endif
            call display_string(15,77,str_solid, 3, 1)
         elseif ( item_number .eq. 21 ) then
c           Interpolation mode
            do_log_interp = .NOT. do_log_interp
            if( do_log_interp ) then
               str_log_interp = '   log'
            else
               str_log_interp = 'linear'
            endif
            call display_string(16,75,str_log_interp, 6, 1)
         elseif ( item_number .eq. 22 ) 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(17,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. 23 ) then
c           Color materials
            if ( do_color ) then
               call display_message(
     1            'Not applicable when using color contours')
            else
               do_colormat = .NOT. do_colormat
               if ( do_colormat ) then
                  str_colormat = 'yes'
               else
                  str_colormat = ' no'
               endif
            endif
            call display_string(18,77,str_colormat, 3, 1)
         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_contours(data_array,
     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     h_vec, v_vec, num_h, num_v,
     6     num_contours, contour_values, label, units,
     7     do_header, do_junctions, do_grid, do_color, do_clabel,
     8     do_yflip, do_solid, do_log, do_colormat,
     9     val_in_oxide, add_notes_from_save, ws_window_open,
     A     add_notes_from_file, ierror)
c
      implicit none
c
c     Subroutine which actually calls graphics routines.
c
c     NOTE: Since we are working with MINIMOS 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_header, do_junctions, do_grid, do_color, do_clabel,
     1        val_in_oxide, add_notes_from_save, ws_window_open,
     2        add_notes_from_file, do_yflip, do_solid, do_log,
     3        do_colormat
      character*(*) label, units, h_label, v_label, h_units, v_units
      integer num_contours, num_h, num_v, ierror
      integer h_minor_tic_freq, v_minor_tic_freq
      real data_array(num_h,num_v), contour_values(1)
      real h_vec(1), v_vec(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 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'units.cmn'
      include 'defaults.cmn'
c
c     locals
c
      logical lerr, exp_flag, is_vtxxx
      integer i, ilen, string_len, window_id, return_stat,
     1        ilen2, itot
      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
      real tmp_contours(10)
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
      call find_max_in_array_subset(
     1        h_win_min, h_win_max, v_win_min, v_win_max,
     1        data_array, xyz(1,idisti), xyz(1,idistj),
     1        MAX(1,ni), MAX(1,nj), MAX(1,ni),
     1        data_is_trigrid, MAX(1,tri_npts),
     1        cord,
     1        do_log, rmax)
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
      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
c     Convert contour values to log if needed
c     Note: we've assumed max number of contours is 10. jvf
      do i=1,MIN(num_contours,10)
         tmp_contours(i) = contour_values(i)
      end do
      if( do_log ) then
         call log10_array(tmp_contours, 10, 1, num_contours, .TRUE.)
      endif
c
      if( do_color .and. can_do_color_contours ) then
c        Set colormap for color contours
         call set_color_contour_colormap
         if( data_is_trigrid ) then
         call color_contours_tri(data_array, cord, nop, tri_ne,
     1            tri2reg, imat,
     2            ARR2DSIZE,
     3            tmp_contours, num_contours,
     4            h_win_min, h_win_max, v_win_min, v_win_max,
     5            val_in_oxide,
     6            max_cont_colors, do_yflip, v_sum)
         else
            if( hidden_line_method .eq. HIDDEN_LINES_HISTOGRAM) then
               call color_contours_histo(data_array, tmp_contours,
     1            num_contours,
     2            h_vec, v_vec, num_h, num_v, num_h,
     3            h_win_min, h_win_max, v_win_min, v_win_max,
     4            max_cont_colors, do_yflip, v_sum)
            else
         call color_contours(data_array, tmp_contours, num_contours,
     1            h_vec, v_vec, num_h, num_v, num_h,
     2            h_win_min, h_win_max, v_win_min, v_win_max,
     3            max_cont_colors, do_yflip, v_sum)
            endif
         endif
c
      else
         if( data_is_trigrid ) then
         call contur_tri(data_array, cord, nop, tri_ne,
     1            tri2reg, imat,
     2            ARR2DSIZE,
     3            tmp_contours, num_contours,
     4            h_win_min, h_win_max, v_win_min, v_win_max,
     5            val_in_oxide,
     6            .FALSE., do_yflip, v_sum, do_clabel, do_solid,
     7            do_colormat)
         else
         call contur(data_array, tmp_contours, num_contours,
     1            h_vec, v_vec, num_h, num_v, num_h,
     2            h_win_min, h_win_max, v_win_min, v_win_max,
     3            .FALSE., do_yflip, v_sum, do_clabel, do_solid)
         endif
      endif
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!
      if ( do_junctions ) then
c        Negate_array call is used merely to copy doping into data_array
         call negate_array(data_array,
     1      arrdat(arrind(net_dop_array_index)), num_h, num_v,
     2      tri_npts, data_is_trigrid, .FALSE.)
         if( do_log ) then
            call log10_array(junction, 1, 1, 1, .TRUE.)
            call log10_array(data_array, num_h, num_v,
     1         tri_npts, .TRUE.)
         endif
         if( data_is_trigrid ) then
         call contur_tri(data_array,
     1            cord, nop, tri_ne,
     2            tri2reg, imat,
     3            ARR2DSIZE, junction, 1,
     4            h_win_min, h_win_max, v_win_min, v_win_max,
     5            .FALSE.,
     6            .TRUE., do_yflip, v_sum, .FALSE., do_solid, .FALSE.)
         else
         call contur(data_array, junction, 1,
     1            h_vec, v_vec, num_h, num_v, num_h,
     2            h_win_min, h_win_max, v_win_min, v_win_max,
     3            .TRUE., do_yflip, v_sum, .FALSE., do_solid)
         endif
      endif
      call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c
c     Draw in oxide and contact boundary information
      call draw_2d_boundaries(v_sum, h_win_min, h_win_max,
     1                        v_win_min, v_win_max, val_in_oxide,
     2                        do_yflip)
c
      if( do_grid ) then
         call plot_2d_grid(binfiletype,
     1           h_vec, v_vec, num_h, num_v,
     1           data_is_trigrid, cord, nop, tri_npts, tri_ne,
     1           v_sum, do_yflip)
      endif
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)
      if( do_header ) then
         call mathmode(ws_id, 0.0, 0.15, plot_title,
     1                 .FALSE., lerr)
      endif
c
c     Print auxiliary information
c
      call TIME(time_string)
      ilen = string_len(label)
      ilen2 = string_len(units)
      work_string = label(1:ilen) // ' ' // left_unit_delim //
     1         units(1:ilen2) // right_unit_delim // ' ' //
     1         date_string  // ' ' // time_string
      call mathmode(ws_id, 0.0, 0.10, work_string, .FALSE., lerr)
      ilen = string_len(binfilename)
      work_string =  'Data source: ' // binfilename(1:ilen)
      call GKS$TEXT(0.0, 0.05, work_string)
c
c     draw legend on side...
      call  plot_color_map_legend(contour_values, num_contours,
     1      rmax, do_color .and. can_do_color_contours, do_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 find_greatest_subset_of_array(hmesh, vmesh,
     1                          low_h, hi_h, low_v, hi_v,
     2                          nh, nv, ih1, ih2, iv1, iv2)
      implicit none
c
c     Find subset of array defined by the integers ih1, ih2, iy1, iv2
c     such that:
c        hmesh(ih1) is the largest  element .le. to low_h
c        hmesh(ih2) is the smallest element .ge. to hi_h
c        vmesh(iv1) is the largest  element .le. to low_v
c        vmesh(iv2) is the smallest element .ge. to hi_v
c
c     Mesh coordinates are stored in hmesh and vmesh
c     
c
      integer nh, nv
      real hmesh(nh), vmesh(nv)
      real low_h, hi_h, low_v, hi_v
      integer ih1, ih2, iv1, iv2
c
      integer getindex
c
c     getindex(a,b,c) returns the index i of b(i) such that
c     b(i) is the largest element of b .le. a
c     This is what we want for ih1 and iv1.
c     For ih2 and iv2, we want need to increase what getindex returns
c     by one to get the the smallest element of b .ge. a
c
      if( low_h .le. hmesh(1) ) then
         ih1 = 1
      else if ( low_h .ge. hmesh(nh) ) then
         ih1 = nh
      else
         ih1 = getindex(low_h, hmesh, nh)
      endif
c
      if( hi_h .le. hmesh(1) ) then
         ih2 = 1
      else if ( hi_h .ge. hmesh(nh) ) then
         ih2 = nh
      else
         ih2 = getindex(hi_h, hmesh, nh)
         if( ih2 .ne. nh  .and. hmesh(ih2) .ne. hi_h ) ih2 = ih2 + 1
      endif
c
      if( low_v .le. vmesh(1) ) then
         iv1 = 1
      else if ( low_v .ge. vmesh(nv) ) then
         iv1 = nv
      else
         iv1 = getindex(low_v, vmesh, nv)
      endif
c
      if( hi_v .le. vmesh(1) ) then
         iv2 = 1
      else if ( hi_v .ge. vmesh(nv) ) then
         iv2 = nv
      else
         iv2 = getindex(hi_v, vmesh, nv)
         if( iv2 .ne. nv  .and. vmesh(iv2) .ne. hi_v ) iv2 = iv2 + 1
      endif
      return
      end
c
      subroutine get_default_contours(hmin,hmax,vmin,vmax,
     1           data_array, isize, jsize, do_log, do_negate,
     2           num_contours, contour_values)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
c
      integer isize, jsize, num_contours
      logical do_log, do_negate
      real hmin, hmax, vmin, vmax, data_array(isize,jsize),
     1     contour_values(num_contours)
c
c     locals
      integer i, j, ih1, ih2, iv1, iv2
      real min_value, max_value
c
      call find_greatest_subset_of_array(xyz(1,idisti), xyz(1,idistj),
     1     hmin, hmax, vmin, vmax, ni, nj, ih1, ih2, iv1, iv2)
c
c     find min/max of data
      min_value = PM_MAX_FLOAT
      max_value = -PM_MAX_FLOAT
      do 100 j=iv1,iv2
         do 100 i=ih1,ih2
            if( min_value .gt. data_array(i,j) ) then
               min_value = data_array(i,j)
            endif
            if( max_value .lt. data_array(i,j) ) then
               max_value = data_array(i,j)
            endif
100   continue
c
      call compute_contours(min_value, max_value, do_log, do_negate,
     1                      .FALSE., num_contours, contour_values)
c
      return
      end
c
      subroutine compute_contours(rmin, rmax, do_log, do_negate,
     1            lexact, num_contours, contour_values)
      implicit none
c
      integer num_contours
      logical do_log, do_negate, lexact
      real rmin, rmax, contour_values(num_contours)
c
      real min_value, max_value, log_max, log_min, sign_max, sign_min,
     1     temp, temp2, delta
      real round_down_whole, round_up_whole
      integer i, imid
c
      min_value = rmin
      max_value = rmax
c
      if( do_negate ) then
         temp = max_value
         max_value = - min_value
         min_value = - temp
      endif
c
      if( do_log ) then
         if( max_value .lt. 0.0 ) then
            sign_max = -1.0
         else
            sign_max =  1.0
         endif
         if( min_value .lt. 0.0 ) then
            sign_min = -1.0
         else
            sign_min = 1.0
         endif
         min_value = ABS(min_value)
         max_value = ABS(max_value)
c        After absolute value, min and max may be reversed
         if( min_value .gt. max_value ) then
            temp = min_value
            min_value = max_value
            max_value = temp
            temp = sign_min
            sign_min = sign_max
            sign_max = temp
         endif
         if( min_value .gt. 0.0 ) then
c           Check for large spread of min/max
            if((ALOG10(max_value) - ALOG10(min_value)).gt. 10) then
c              Force steps of 10 in contours
               min_value = 0.0
            endif
         endif
         if( min_value .eq. 0.0 ) sign_min = sign_max
         if( max_value .eq. 0.0 ) sign_max = sign_min
         if( sign_min .ne. sign_max ) then
c           Signs are unequal. Count down powers of ten from max,
c           count up from min
            if( lexact ) then
               temp = sign_max*max_value
            else
               log_max = ALOG10(max_value)
               if( log_max .lt. 0.0 ) then
                     temp = sign_max*10.0**round_up_whole(log_max)
               else
                     temp = sign_max*10.0**round_down_whole(log_max)
               endif
            endif
            if( lexact ) then
               temp2 = sign_min*min_value
            else
               log_min = ALOG10(min_value)
               if( log_min .lt. 0.0 ) then
                  temp2 = sign_min*10.0**round_up_whole(log_min)
               else
                  temp2 = sign_min*10.0**round_down_whole(log_min)
               endif
            endif
c           Select an equal number of positive and negative contours
c           If number of contours is odd, then favor the larger
c           (in magnitude) extremum.
            imid = num_contours/2
            if( MOD(num_contours,2) .eq. 1 ) then
               imid = imid + 1
               if( max_value .lt. min_value ) imid = imid - 1
            endif
            do 200 i=1,imid
               contour_values(i) = temp
               temp = temp/10.0
200         continue
            do 201 i=num_contours,imid+1,-1
               contour_values(i) = temp2
               temp2 = temp2/10.0
201         continue
         else
c           Signs are equal. Spread contours evenly (log scale)
c           Protect against log(0.0) (input data was all zeroes)
            if( max_value .eq. 0.0 ) max_value = 1.0
            if( min_value .eq. 0.0 ) then
c              No lower bound -- step by 10.0
               if( lexact ) then
                  temp = sign_max*max_value
               else
c                 Round first contour to nearest power of ten
                  log_max = ALOG10(max_value)
                  if( log_max .lt. 0.0 ) then
                     temp = sign_max*10.0**round_up_whole(log_max)
                  else
                     temp = sign_max*10.0**round_down_whole(log_max)
                  endif
               endif
               do i=1,num_contours
                  contour_values(i) = temp
                  temp = temp/10.0
               end do
            else
c              Split min/max into equal logarithmic steps
               temp = sign_max*max_value
               if( lexact ) then
                  delta = 10.0**(ALOG10(max_value/min_value)/
     1                          (num_contours-1))
                  do i=1,num_contours
                     contour_values(i) = temp
                     temp = temp/delta
                  end do
               else
                  delta = 10.0**(ALOG10(max_value/min_value)/
     1                          (num_contours+1))
                  do i=1,num_contours
                     temp = temp/delta
                     contour_values(i) = temp
                  end do
               endif
            endif
         endif
      else
         if( lexact ) then
            do 300 i=1,num_contours
               contour_values(i) = min_value + (max_value - min_value)*
     1                             FLOAT(i-1)/FLOAT(num_contours-1)
300         continue
         else
            do 301 i=1,num_contours
               contour_values(i) = min_value + (max_value - min_value)*
     1                             FLOAT(i)/FLOAT(num_contours+1)
301         continue
         endif
      endif
c
      return
      end
c
      subroutine parse_contours(buffer, max_contours,
     1                          ncont, contours, lerror)
      implicit none
c
      character*(*) buffer
      integer max_contours, ncont
      logical lerror
      real contours(max_contours)
c
      character*50 item
      integer ipos, nsteps
      real delta, rmin, rmax
c
      lerror = .FALSE.
      ncont = 0
      ipos  = 1
100      call get_keyword(buffer,item,ipos)
         if( ipos .eq. 0 ) goto 200
         if( item(1:1) .eq. 'E' ) then
c           n equidistant contours between min and max
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) nsteps
            if( nsteps .le. 0 ) goto 1000
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmin
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmax
            if( rmin .ge. rmax ) goto 1000
            if( (nsteps + ncont) .gt. max_contours ) goto 1000
            call compute_contours(rmin, rmax, .FALSE., .FALSE.,
     1              .TRUE., nsteps, contours(ncont+1))
            ncont = ncont + nsteps
         else if( item(1:1) .eq. 'L' ) then
c           n log equidistant contours between min and max
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) nsteps
            if( nsteps .le. 0 ) goto 1000
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmin
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmax
            if( rmin .ge. rmax ) goto 1000
            if( rmin .eq. 0.0 .or. rmax .eq. 0.0 ) goto 1000
            if( (nsteps + ncont) .gt. max_contours ) goto 1000
            call compute_contours(rmin, rmax, .TRUE., .FALSE.,
     1              .TRUE., nsteps, contours(ncont+1))
            ncont = ncont + nsteps
         else if( item(1:1) .eq. 'S' ) then
c           step delta (additive) between min and max
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) delta
            if( delta .le. 0 ) goto 1000
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmin
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmax
            if( rmin .ge. rmax ) goto 1000
            nsteps = NINT((rmax - rmin)/delta + 1)
c           Must have at least two steps (contours)
            if( nsteps .eq. 1 ) goto 1000
            if( (nsteps + ncont) .gt. max_contours ) goto 1000
            call compute_contours(rmin, rmax, .FALSE., .FALSE.,
     1              .TRUE., nsteps, contours(ncont+1))
            ncont = ncont + nsteps
         else if( item(1:1) .eq. 'D' ) then
c            step delta (multiplicative) between min and max
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) delta
            if( delta .le. 0 ) goto 1000
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmin
            call get_keyword(buffer,item,ipos)
            if( ipos .eq. 0 ) goto 1000
            read(item,*,end=1000,err=1000) rmax
            if( rmin .ge. rmax ) goto 1000
            if( rmin .eq. 0.0 .or. rmax .eq. 0.0 ) goto 1000
c           Can't cross zero with a log step
            if( rmin .lt. 0.0 .and. rmax .gt. 0.0 ) goto 1000
            if( rmin .gt. 0.0 ) then
               nsteps = NINT( (ALOG10(rmax) - ALOG10(rmin))/
     1                         ALOG10(delta) ) + 1
            else
               nsteps = NINT( (ALOG10(ABS(rmin)) - ALOG10(ABS(rmax)))/
     1                         ALOG10(delta) ) + 1
            endif
            if( (nsteps + ncont) .gt. max_contours ) goto 1000
            call compute_contours(rmin, rmax, .TRUE., .FALSE.,
     1              .TRUE., nsteps, contours(ncont+1))
            ncont = ncont + nsteps
         else
            if( ncont+1 .gt. max_contours ) goto 1000
            ncont = ncont + 1
            read(item,*,end=1000,err=1000)
     1         contours(ncont)
         endif
         goto 100
200   continue
      if( ncont .eq. 0 ) lerror = .TRUE.
c
c       Branch over error handler
        goto 1001
c
c       Error handler
1000    continue
        lerror = .TRUE.
c
1001    continue
c
      return
      end
c
      subroutine draw_2d_boundaries(v_sum, hmin, hmax,
     1                              vmin, vmax, val_in_oxide, do_yflip)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'data.cmn'
      include 'mesh.cmn'
c
      real hmin, hmax, vmin, vmax, v_sum
      logical val_in_oxide, do_yflip
c
      if( binfiletype .eq. DT_MINIMOS ) then
         call draw_2d_bound_mini(v_sum, val_in_oxide)
      else if( binfiletype .eq. DT_VLSICAP ) then
         call draw_2d_bound_vlsi(v_sum, val_in_oxide)
      else if( binfiletype .eq. DT_PROMIS ) then
         call draw_2d_bound_prom(v_sum, val_in_oxide)
      else if( binfiletype .eq. DT_2DOP ) then
         call draw_2d_bound_2dop(v_sum, val_in_oxide)
      else if( binfiletype .eq. DT_PISCES .or.
     1         binfiletype .eq. DT_SUPREM4 .or.
     2         binfiletype .eq. DT_TIF .or.
     3         binfiletype .eq. DT_MEDICI .or.
     5        (binfiletype .eq. DT_ASCII  .and. data_is_trigrid)) then
         call draw_2d_bound_tri(cord, nop, imat, tri_nextel, tri2reg,
     1                          tri_npts, tri_ne, tri_nreg,
     1                          v_sum, hmin, hmax, vmin, vmax,
     1                          do_yflip)
      endif
c
      return
      end
c
      subroutine draw_2d_bound_mini(v_sum, val_in_oxide)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'defaults.cmn'
c
      real v_sum
      logical val_in_oxide
      integer i, nysem
      real fox_temp(2*ARRAYSIZE), dist_temp(2*ARRAYSIZE)
c
      call set_line_index(DEFAULT_LINE)
      nysem = nyoxm + ny
      if( nonplanar_boundary .or. three_d ) then
c        Use fill area to blot out any contours that "seep out"
c        over into areas where the function value is not defined
c        (e.g. electron density in oxide). Note, this will only work
c        on devices which support fill in the background color.
c        The alternative is to clip the contours at the (possibly
c        nonplanar) boundaries.
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c        foxo, foxu have already been set dependent on plane type
c        from 2D oxide/contact array in binary file
         if( plane_type .eq. XY_PLANE ) then
            call GKS$SET_FILL_COLOR_INDEX(0)
            do i=1,fox_nx
               fox_temp(i) = v_sum - foxu(i)
               dist_temp(i) = fox_x(i)
            end do
            do i=1,fox_nx
               fox_temp(i+fox_nx) = v_sum - foxo(fox_nx-i+1)
               dist_temp(i+fox_nx) = fox_x(fox_nx-i+1)
            end do
            if( .NOT. val_in_oxide .and. l_oxide_fill )
     1          call GKS$FILL_AREA(2*fox_nx,dist_temp,fox_temp)
c           Draw oxide boundary
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
            call GKS$POLYLINE(fox_nx,fox_x,fox_temp)
c           contact boundary
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
            do i=1,fox_nx
               fox_temp(i) = v_sum - foxo(i)
               dist_temp(i) = fox_x(i)
            end do
            do i=1,fox_nx
               fox_temp(i+fox_nx) = v_sum - xyz(1,2)
               dist_temp(i+fox_nx) = fox_x(fox_nx-i+1)
            end do
c           White out contact area if no values there (cheap workaround
c           for contours "seeping" into oxide area)
            if( .NOT. val_in_oxide .and. l_oxide_fill )
     1          call GKS$FILL_AREA(2*fox_nx,dist_temp,fox_temp)
c           If SOI, then fill in oxide film
            if( nysem .lt. nytot ) then
               fox_temp(1) = xyz(1,1)
               dist_temp(1) = v_sum - xyz(nysem,2)
               fox_temp(2) = xyz(nx,1)
               dist_temp(2) = v_sum - xyz(nysem,2)
               fox_temp(3) = xyz(nx,1)
               dist_temp(3) = v_sum - xyz(nytot,2)
               fox_temp(4) = xyz(1,1)
               dist_temp(4) = v_sum - xyz(nytot,2)
               if( .NOT. val_in_oxide .and. l_oxide_fill )
     1            call gks$fill_area(4,fox_temp,dist_temp)
c              Draw line for SOI interface
               call move_abs(xyz(1,1),  v_sum - xyz(nysem,2))
               call draw_abs(xyz(nx,1), v_sum - xyz(nysem,2))
            endif
c           Draw in contact outline
            call GKS$POLYLINE(fox_nx,fox_x,fox_temp)
            call GKS$SET_FILL_COLOR_INDEX(1)
         else if( plane_type .eq. YZ_PLANE ) then
            call GKS$SET_FILL_COLOR_INDEX(0)
            do i=1,fox_nz
               fox_temp(i) = v_sum - foxu(i)
               dist_temp(i) = fox_z(i)
            end do
            do i=1,fox_nz
               fox_temp(i+fox_nz)  = v_sum - foxo(fox_nz-i+1)
               dist_temp(i+fox_nz) = fox_z(fox_nz-i+1)
            end do
            if( .NOT. val_in_oxide .and. l_oxide_fill )
     1         call GKS$FILL_AREA(2*fox_nz,dist_temp,fox_temp)
c           Draw oxide boundary
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
            call GKS$POLYLINE(fox_nz,fox_z,fox_temp)
c           If SOI, then fill in oxide film
            if( nysem .lt. nytot ) then
               fox_temp(1) = xyz(1,3)
               dist_temp(1) = v_sum - xyz(nysem,2)
               fox_temp(2) = xyz(nz,3)
               dist_temp(2) = v_sum - xyz(nysem,2)
               fox_temp(3) = xyz(nz,3)
               dist_temp(3) = v_sum - xyz(nytot,2)
               fox_temp(4) = xyz(1,3)
               dist_temp(4) = v_sum - xyz(nytot,2)
               if( .NOT. val_in_oxide .and. l_oxide_fill )
     1            call gks$fill_area(4,fox_temp,dist_temp)
c              Draw line for SOI interface
               call move_abs(xyz(1,3),  v_sum - xyz(nysem,2))
               call draw_abs(xyz(nz,3), v_sum - xyz(nysem,2))
            endif
c           Draw contact boundary
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
            do i=1,fox_nz
               fox_temp(i) = v_sum - foxo(i)
            end do
            call GKS$POLYLINE(fox_nz,fox_z,fox_temp)
            call GKS$SET_FILL_COLOR_INDEX(1)
         else
c           Do nothing for top plane (XZ)
         endif
      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. val_in_oxide ) then
            call GKS$SET_FILL_COLOR_INDEX(0)
            fox_temp(1) = xyz(1,1)
            dist_temp(1) = v_sum
            fox_temp(2) = xyz(nx,1)
            dist_temp(2) = v_sum
            fox_temp(3) = xyz(nx,1)
            dist_temp(3) = v_sum - xyz(1,2)
            fox_temp(4) = xyz(1,1)
            dist_temp(4) = v_sum - xyz(1,2)
            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) = xyz(1,1)
               dist_temp(1) = v_sum - xyz(nysem,2)
               fox_temp(2) = xyz(nx,1)
               dist_temp(2) = v_sum - xyz(nysem,2)
               fox_temp(3) = xyz(nx,1)
               dist_temp(3) = v_sum - xyz(nytot,2)
               fox_temp(4) = xyz(1,1)
               dist_temp(4) = v_sum - xyz(nytot,2)
               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(xyz(1,1),  v_sum)
         call draw_abs(xyz(nx,1), v_sum)
c        Draw line for SOI interface
         if( nysem .lt. nytot ) then
            call move_abs(xyz(1,1),  v_sum - xyz(nysem,2))
            call draw_abs(xyz(nx,1), v_sum - xyz(nysem,2))
         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(xyz(nxsp-1,1), v_sum)
         call draw_abs(xyz(nxsp-1,1),(v_sum - xyz(1,2)))
c
         if( mini_smile ) then
            call move_abs(0.0,(v_sum - xyz(nytop,2)))
            do i=nxg1,nxg2
               call draw_abs(xyz(i,1),(v_sum - foxu(i)))
            end do
            call draw_abs(xl, (v_sum - xyz(nytop,2)))
            call draw_abs(0.0,(v_sum - xyz(nytop,2)))
         else
            call move_abs(0.0,(v_sum - xyz(nytop,2)))
            call draw_abs(0.0,(v_sum + tox))
            call draw_abs(xl, (v_sum + tox))
            call draw_abs(xl, (v_sum - xyz(nytop,2)))
            call draw_abs(0.0,(v_sum - xyz(nytop,2)))
         endif
c
         call move_abs(xyz(nxsoxm+1,1),  v_sum)
         call draw_abs(xyz(nxsoxm+1,1), (v_sum - xyz(1,2)))
      endif
c
      return
      end
c
      subroutine draw_2d_bound_vlsi(v_sum, val_in_oxide)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
c
      real v_sum
      logical val_in_oxide
c
c     locals
      integer idom, j, k, jpoly
      real xdum(2), ydum(2)
c
      if( .NOT. read_vlsi_inp_file ) return
c
      do idom=1,ndomain
c        Make electrodes a different color than other boundaries
         if( ABS(domaintype(idom)) .eq. 11 ) then
            call GKS$SET_PLINE_COLOR_INDEX(3)
         else
            call GKS$SET_PLINE_COLOR_INDEX(1)
         endif
         do j=1,MAX_POLYLIST
            if( polygon_list(idom,j) .eq. 0 ) goto 100
c           Ignore sense of polygon (pos or neg) since this is
c           immaterial when drawing a polyline
            jpoly = ABS(polygon_list(idom,j))
            do k=polygon2vertex(jpoly,1),
     1           polygon2vertex(jpoly,2)-1
               xdum(1) = vertex_coords(k,1)
               xdum(2) = vertex_coords(k+1,1)
               ydum(1) = v_sum - vertex_coords(k,2)
               ydum(2) = v_sum - vertex_coords(k+1,2)
               call gks$polyline(2,xdum,ydum)
            end do
         end do
100      continue
      end do
c
      return
      end
c
      subroutine draw_2d_bound_prom(v_sum, val_in_oxide)
      implicit none
c
c     Draw oxide body shape (into silicon)
c     Depends on having oxide shape in save file in location 101
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'defaults.cmn'
c
      real v_sum
      logical val_in_oxide
c
      integer i, j, iox
      logical do_oxide
      real fox_temp(2*ARRAYSIZE+2), dist_temp(2*ARRAYSIZE)
c
      call set_line_index(DEFAULT_LINE)
      do_oxide = .FALSE.
      do i=1,promis_nequ
         if( nelem(i) .eq. 101 ) then
            do_oxide = .TRUE.
            iox = i
         endif
      end do
c     Note: if data is on a triangular grid, then the oxide information
c     is (usually) placed in the arr2d arrays, with xyz containing the
c     lateral position.
      if( do_oxide .or. (data_is_trigrid .and. nx .gt. 1) ) then
c        Use fill area to blot out any contours that "seep out"
c        over into areas where the function value is not defined
c        (e.g. doping concentration in oxide). Note, this will only work
c        on devices which support fill in the background color.
c        The alternative is to clip the contours at the (possibly
c        nonplanar) boundaries.
         call GKS$SET_FILL_COLOR_INDEX(0)
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c        Lower oxide shape is given in microns (convert to cm)
         do i=1,ni
            fox_temp(i) = v_sum - 1.0E-4*arr2d(i,1)
            dist_temp(i) = xyz(i,idisti)
         end do
c        Upper oxide shape is given in microns (convert to cm)
c        Traverse 'backwards' to complete oxide polygon
         j = ni
         do i=ni,1,-1
            j = j + 1
            fox_temp(j) = v_sum - 1.0E-4*(arr2d(i,1)-arr2d(i,2))
            dist_temp(j) = xyz(i,idisti)
         end do
c        Draw oxide shape
         if( l_oxide_fill ) then
            call GKS$FILL_AREA(2*ni,dist_temp,fox_temp)
         endif
c        Draw oxide boundary
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
         call GKS$POLYLINE(ni,xyz(1,idisti),fox_temp)
c        Draw upper oxide shape (compute from oxide thickness)
c        Upper oxide shape is given in microns (convert to cm)
         do i=1,ni
            fox_temp(i) = v_sum - 1.0E-4*(arr2d(i,1)-arr2d(i,2))
            dist_temp(i) = xyz(i,idisti)
         end do
         call GKS$POLYLINE(ni,xyz(1,idisti),fox_temp)
      endif
      return
      end
c
      subroutine draw_2d_bound_2dop(v_sum, val_in_oxide)
      implicit none
c
c     Draw oxide body shape on doping file
c     Depends on having read associated interface file
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'defaults.cmn'
c
      real v_sum
      logical val_in_oxide
c
      integer i
      real fox_temp(2*ARRAYSIZE), dist_temp(2*ARRAYSIZE)
c
      call set_line_index(DEFAULT_LINE)
c
c     ASSUMES oxide shape is in first two 1D data positions
      if( in_binary_file(4) .and. in_binary_file(5) ) then
c        Use fill area to blot out any contours that "seep out"
c        over into areas where the function value is not defined
c        (e.g. doping concentration in oxide). Note, this will only work
c        on devices which support fill in the background color.
c        The alternative is to clip the contours at the (possibly
c        nonplanar) boundaries.
         call GKS$SET_FILL_COLOR_INDEX(0)
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c        Reminder: oxide shape has been negated on read in so that
c        it plots 'upright' on X-Y plots. Convert to cm.
         do i=1,ni
            fox_temp(i) = v_sum + 1.0E-4*arr2d(i,2)
            dist_temp(i) = xyz(i,idisti)
         end do
c        Complete polygon, draw oxide shape
         do i=1,ni
            fox_temp(i+ni) = v_sum + 1.0E-4*arr2d(ni-i+1,1)
            dist_temp(i+ni) = xyz(ni-i+1,idisti)
         end do
         if( l_oxide_fill) then
            call GKS$FILL_AREA(2*ni,dist_temp,fox_temp)
         endif
c        Draw oxide boundary
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_DASHED)
         call GKS$POLYLINE(ni,xyz(1,idisti),fox_temp)
         call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_SOLID)
c        Draw upper oxide shape, convert to cm
         do i=1,ni
            fox_temp(i) = v_sum + 1.0E-4*arr2d(i,1)
            dist_temp(i) = xyz(i,idisti)
         end do
         call GKS$POLYLINE(ni,xyz(1,idisti),fox_temp)
      endif
      return
      end
c
      subroutine find_max_in_array_subset(hmin, hmax, vmin, vmax,
     1           data_array, x, y, nx, ny, nsize,
     1           data_is_trigrid, tri_npts, cord, do_log, rmax)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
c
      integer tri_npts, nx, ny, nsize
      logical data_is_trigrid, do_log
      real hmin, hmax, vmin, vmax, data_array(*),
     1     rmax, x(nx), y(ny), cord(2,tri_npts)
c
c     locals
      integer i, j, ih1, ih2, iv1, iv2
c
      if( data_is_trigrid ) then
         rmax = -PM_MAX_FLOAT
         do 200 i=1,tri_npts
c           Is point in subset region?
            if( cord(1,i) .lt. hmin .or. cord(1,i) .gt. hmax .or.
     1          cord(2,i) .lt. vmin .or. cord(2,i) .gt. vmax) goto 200
            rmax = MAX(rmax,data_array(i))
200      continue
      else
         call find_greatest_subset_of_array(x, y,
     1        hmin, hmax, vmin, vmax, nx, ny, ih1, ih2, iv1, iv2)
c
c        find max of data
         rmax = -PM_MAX_FLOAT
         do 100 j=iv1,iv2
           do 100 i=ih1,ih2
               if( rmax .lt. data_array(i+(j-1)*nsize) ) then
                  rmax = data_array(i+(j-1)*nsize)
               endif
100      continue
      endif
c
      if( do_log ) then
         rmax = SIGN(10.0**MIN(38.0,ABS(rmax)),rmax)
      endif
      return
      end
c
      subroutine find_max_in_array_subset_old(hmin,hmax,vmin,vmax,
     1           data_array, nsize, data_is_trigrid, do_log, rmax)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
c
      integer nsize
      logical data_is_trigrid, do_log
      real hmin, hmax, vmin, vmax, data_array(nsize*nsize),
     1     rmax
c
c     locals
      integer i, j, ih1, ih2, iv1, iv2
c
      if( data_is_trigrid ) then
         rmax = -PM_MAX_FLOAT
         do 200 i=1,tri_npts
c           Is point in subset region?
            if( cord(1,i) .lt. hmin .or. cord(1,i) .gt. hmax .or.
     1          cord(2,i) .lt. vmin .or. cord(2,i) .gt. vmax) goto 200
            rmax = MAX(rmax,data_array(i))
200      continue
      else
         call find_greatest_subset_of_array(xyz(1,idisti),xyz(1,idistj),
     1        hmin, hmax, vmin, vmax, ni, nj, ih1, ih2, iv1, iv2)
c
c        find max of data
         rmax = -PM_MAX_FLOAT
         do 100 j=iv1,iv2
           do 100 i=ih1,ih2
               if( rmax .lt. data_array(i+(j-1)*nsize) ) then
                  rmax = data_array(i+(j-1)*nsize)
               endif
100      continue
      endif
c
      if( do_log ) then
         rmax = SIGN(10.0**MIN(38.0,ABS(rmax)),rmax)
      endif
      return
      end
c
      subroutine plot_2d_grid(type, h_vec, v_vec, num_h, num_v,
     1           data_is_trigrid, cord, nop, tri_npts, tri_ne,
     1           v_sum, do_yflip)
      implicit none
c
c     Plot grid lines over contours
c
      include 'const.h'
c
      logical data_is_trigrid, do_yflip
      integer type, num_h, num_v, tri_npts, tri_ne, nop(3,tri_ne)
      real h_vec(num_h), v_vec(num_v), cord(2,tri_npts), v_sum
c
      integer i,j
c
      call set_line_index(DEFAULT_LINE)
c
      if( type .eq. DT_VLSICAP) then
         call plot_vlsicap_mesh(v_sum)
      else if( type .eq. DT_SUPREM4 .or.
     1         type .eq. DT_PISCES  .or.
     2         type .eq. DT_TIF     .or.
     3         type .eq. DT_MCP2D   .or.
     4         type .eq. DT_MEDICI  .or.
     6        (type .eq. DT_ASCII  .and. data_is_trigrid) .or.
     7        (type .eq. DT_PROMIS .and. data_is_trigrid)
     8   ) then
         call plot_tri_mesh(cord, nop, tri_npts, tri_ne,
     1                      v_sum, do_yflip)
      else
         if( do_yflip ) then
            do 10 j=1,num_v
               call move_abs(h_vec(1), (v_sum - v_vec(j)))
               call draw_abs(h_vec(num_h),(v_sum - v_vec(j)))
10          continue
            do 11 i=1,num_h
               call move_abs(h_vec(i),(v_sum - v_vec(1)))
               call draw_abs(h_vec(i),(v_sum - v_vec(num_v)))
11          continue
         else
            do 12 j=1,num_v
               call move_abs(h_vec(1), v_vec(j))
               call draw_abs(h_vec(num_h),v_vec(j))
12          continue
            do 13 i=1,num_h
               call move_abs(h_vec(i),v_vec(1))
               call draw_abs(h_vec(i),v_vec(num_v))
13          continue
         endif
      endif
      return
      end
c
      subroutine plot_vlsicap_mesh(v_sum)
      implicit none
c
      real v_sum
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'data.cmn'
      include 'io.cmn'
c
c     Plot VLSICAP triangular mesh
c
c     locals
c
      character*132 lgf_file
      character*80 buffer
      integer i, j, itri, ioerr, ilen, string_len, nrecl, ierr
      real xdum(4), ydum(4)
c
c     Form LGF filespec (assumed to be in sibling directory 'lgf_vlsi')
      call get_parent_directory(binfilename,lgf_file)
      call add_subdirectory(lgf_file,'lgf_vlsi')
      ilen = string_len(lgf_file)
      lgf_file = lgf_file(1:ilen) // binfilestem
      call add_extension(lgf_file,'.lgf')
c
      write(luttyo,10)
10    format(' ','Reading LGF mesh file (wow, is this a long file...',$)
c
      nrecl=0
      call pm_open_file(binfil, lgf_file,
     1                 'OLD','FORMATTED','READ','FORTRAN',
     2                 'NONSHARED', nrecl, ioerr, ierr)
      if( ierr .ne. 0 ) goto 10010
c
c     Skip over previous triangular meshes
      do i=1,vlsi_cycle-1
         read(binfil,100,end=10020,err=10020) buffer
100      format(a)
         read(buffer(8:18),*,end=10020,err=10020) itri
         do j=1,2*itri
            read(binfil,*)
         end do
      end do
c
c     Get final mesh
      call gks$set_pline_color_index(1)
      read(binfil,100,end=10020,err=10020) buffer
      read(buffer(8:18),*,end=10020,err=10020) itri
      do j=1,itri
         read(binfil,*,end=10020,err=10020) (xdum(i),ydum(i),i=1,3)
         do i=1,3
c           Convert coordinates from meters to cm,
c           shift X,Y and flip Y coordinate.
            xdum(i) = 100.0*xdum(i) - vlsi_edge
            ydum(i) = v_sum - (vlsi_yint - 100.0*ydum(i))
         end do
         xdum(4) = xdum(1)
         ydum(4) = ydum(1)
         call gks$polyline(4,xdum,ydum)
      end do
c
      close(unit=binfil)
      goto 10000
c
10010 continue
      write(luttyo,10011) ioerr
10011 format(' ','FORTRAN error ',i5,' during open of lgf file')
      goto 10000
c
10020 continue
      write(luttyo,10021)
10021 format(' ','Error while reading lgf file')
      close(unit=binfil)
      goto 10000
c
10000 continue
c
      return
      end
c
      subroutine sort_contours(contour_values,num_contours)
      implicit none
c
      integer num_contours
      real contour_values(num_contours)
c
c     locals
      integer i,j
      real temp
c
      do 100 j=1,num_contours-1
         do 200 i=1,num_contours-j
            if( contour_values(i) .gt. contour_values(i+1) ) then
               temp = contour_values(i)
               contour_values(i) = contour_values(i+1)
               contour_values(i+1) = temp
            endif
200      continue
100   continue
c
      return
      end
c
      subroutine get_default_contours_tri(hmin,hmax,vmin,vmax,
     1           data_array, nsize, do_log, do_negate,
     2           num_contours, contour_values)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
c
      integer nsize, num_contours
      logical do_log, do_negate
      real hmin, hmax, vmin, vmax, data_array(nsize),
     1     contour_values(num_contours)
c
c     locals
      integer i
      real min_value, max_value
c     Statement function for test of point outide region
      real x, y, h1, h2, v1, v2
      logical loutside
      loutside(x,y,h1,h2,v1,v2) = x .lt. h1 .or. x .gt. h2 .or.
     1                            y .lt. v1 .or. y .gt. v2
c
c     find min/max of data
      min_value = PM_MAX_FLOAT
      max_value = -PM_MAX_FLOAT
      do 100 i=1,tri_npts
         if( loutside(cord(1,i),cord(2,i), hmin, hmax, vmin, vmax))
     1      goto 100
         if( min_value .gt. data_array(i) ) then
            min_value = data_array(i)
         endif
         if( max_value .lt. data_array(i) ) then
            max_value = data_array(i)
         endif
100   continue
c
c     Check to make sure we found some data points...
      if( min_value .eq. max_value ) then
         if( min_value .lt. 0.0 ) then
            min_value =  1.1*min_value
            max_value =  0.9*max_value
         else if ( min_value .gt. 0.0 ) then
            min_value =  0.9*min_value
            max_value =  1.1*max_value
         else
            min_value = -1.0
            max_value =  1.0
         endif
         call compute_contours(min_value, max_value, do_log, do_negate,
     1                      .FALSE., num_contours, contour_values)
      else if( min_value .lt. max_value ) then
         call compute_contours(min_value, max_value, do_log, do_negate,
     1                      .FALSE., num_contours, contour_values)
      else
         do i=1,num_contours
            contour_values(i) = 0.0
         end do
      endif
c
      return
      end
c
      subroutine plot_tri_mesh(cord, nop, tri_npts, tri_ne,
     1                         v_sum, do_yflip)
      implicit none
c
      integer tri_npts, tri_ne, nop(3,tri_ne)
      real cord(2,tri_npts)
      real v_sum
      logical do_yflip
c
      include 'const.h'
c
c     Plot triangular mesh
c
c     locals
c
      integer i, j
      real xdum(4), ydum(4)
c
      call set_line_index(DEFAULT_LINE)
c
      do j=1,tri_ne
         do i=1,3
c           Flip Y coordinates
            xdum(i) = cord(1,nop(i,j))
            if( do_yflip ) then
               ydum(i) = v_sum - cord(2,nop(i,j))
            else
               ydum(i) = cord(2,nop(i,j))
            endif
         end do
         xdum(4) = xdum(1)
         ydum(4) = ydum(1)
         call gks$polyline(4,xdum,ydum)
      end do
c
      return
      end
c
      subroutine draw_2d_bound_tri(cord, nop, imat, tri_nextel, tri2reg,
     1                             tri_npts, tri_ne, tri_nreg,
     1                             v_sum, hmin, hmax, vmin, vmax,
     1                             do_yflip)
      implicit none
c
      integer tri_npts, tri_ne, tri_nreg
      real    cord(2,tri_npts)
      integer nop(3,tri_ne), tri_nextel(3,tri_ne), tri2reg(tri_ne),
     1        imat(tri_nreg)
      real v_sum, hmin, hmax, vmin, vmax
      logical val_in_oxide, do_yflip
c
      include 'const.h'
      include 'sysdep.cmn'
c
      real xmax, xmin, ymax, ymin, xx(3), yy(3), xtemp(2), ytemp(2)
      integer ie, i, ntemp, node(3), nside(3), iside, ieladj,
     1        nodefr, nodeto
      logical ltest
      data nside/2,3,1/
c
      call set_line_index(DEFAULT_LINE)
c
c     Scan all elements
      do 100 ie=1,tri_ne
         xmin = PM_MAX_FLOAT
         xmax = - PM_MAX_FLOAT
         ymin = PM_MAX_FLOAT
         ymax = - PM_MAX_FLOAT
c        Look for boundaries at each element
         do 110 i=1,3
            ntemp = nop(i,ie)
            node(i) = ntemp
            xx(i)   = cord(1,ntemp)
            yy(i)   = cord(2,ntemp)
            xmin = MIN(xmin,xx(i))
            xmax = MAX(xmax,xx(i))
            ymin = MIN(ymin,yy(i))
            ymax = MAX(ymax,yy(i))
110      continue
c        Is element outside plot boundaries?
         if( xmax .lt. hmin .or. hmax .lt. xmin .or.
     1       ymax .lt. vmin .or. vmax .lt. ymin ) goto 100
c        Scan all three sides of element
         nodefr = 3
         do 120 nodeto=1,3
            iside = nside(nodeto)
            ieladj = tri_nextel(iside,ie)
c           Is this side a material boundary? Break test in two because
c           we can't guarantee order of evaluation in an if statement,
c           and we might run outside the tri2reg() array.
            if( ieladj .gt. 0 ) then
               ltest =  (imat(tri2reg(ie)) .ne. imat(tri2reg(ieladj)))
            else
               ltest = .FALSE.
            endif
c           Is this side an external boundary or material boundary, plot
            if( ieladj .le. 0 .or. ltest ) then
               xtemp(1) = xx(nodefr)
               xtemp(2) = xx(nodeto)
c              Flip Y coordintes
               if( do_yflip ) then
                  ytemp(1) = v_sum - yy(nodefr)
                  ytemp(2) = v_sum - yy(nodeto)
               else
                  ytemp(1) = yy(nodefr)
                  ytemp(2) = yy(nodeto)
               endif
               call gks$polyline(2,xtemp,ytemp)
            endif
            nodefr = nodeto
120      continue
100   continue
c
      return
      end
c
      subroutine pm_set_plot_bounds_to_scale(aspect_ratio, fixed_axis,
     1           xmin, xmax, x_tic_major, x_minor_tic_freq,
     2           ymin, ymax, y_tic_major, y_minor_tic_freq)
      implicit none
c
c     Reset X/Y plot parameters to scale with aspect_ratio
c     Aspect ratio is defined as x/y
c
c     If fixed_axis == 'X', keep X axis fixed, scaled Y
c     If fixed_axis == 'Y', keep Y axis fixed, scaled X
c     If fixed_axis == 'E' (either) scale both axes for best fit
c
      real aspect_ratio, xmin, xmax, x_tic_major,
     1     ymin, ymax, y_tic_major
      integer x_minor_tic_freq, y_minor_tic_freq
      character*1 fixed_axis
c
      real delx, dely, midx, midy, delta, rdum1, rdum2
      integer idum1
c
      delx = xmax - xmin
      dely = ymax - ymin
      midx = 0.5*(xmax + xmin)
      midy = 0.5*(ymax + ymin)
c
      if( fixed_axis .eq. 'X' ) then
c        Expand/contract Y bounds, centered on original y bounds
         delta = 0.5*delx/aspect_ratio
         ymin  = midy - delta
         ymax  = midy + delta
C
C         ymax = ymin + delx/aspect_ratio
      else if( fixed_axis .eq. 'Y' ) then
c        Expand/contract X bounds, centered on original x bounds
         delta = 0.5*dely*aspect_ratio
         xmin  = midx - delta
         xmax  = midx + delta
C
C         xmax = xmin + dely*aspect_ratio
      else
c        Try 'best fit'. Strategy is always to expand bounds
         if( delx .gt. aspect_ratio*dely ) then
            delta = 0.5*delx/aspect_ratio
            ymin  = midy - delta
            ymax  = midy + delta
C
C            ymax = ymin + delx/aspect_ratio
         else
            delta = 0.5*dely*aspect_ratio
            xmin  = midx - delta
            xmax  = midx + delta
C
C            xmax = xmin + dely*aspect_ratio
         endif
      endif
c
c     Recompute 'nice' tic marks
      call get_default_limits(xmin, xmax, rdum1, rdum2,
     1                 x_tic_major, x_minor_tic_freq,
     2                 idum1, .FALSE.)
      x_tic_major = x_tic_major*10.0**idum1
c
      call get_default_limits(ymin, ymax, rdum1, rdum2,
     1                 y_tic_major, y_minor_tic_freq,
     2                 idum1, .FALSE.)
      y_tic_major = y_tic_major*10.0**idum1
c
      return
      end
