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     3D surface plots
c
      subroutine plot3d
      implicit none
c
c     Asks user what array to plot, then plot it!
c
c     Programmer: John Faricelli
c     Date: 14-Sep-87
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 surface of 1-D data'/
     1             ' ','Hit return to continue...')
            read(luttyi,106,end=300,err=300) buffer
         else if( data_is_trigrid .and. default_hidden_method .eq.
     1            HIDDEN_LINES_HISTOGRAM ) then
            call clear_screen
            call beep_terminal
            write(luttyo,101)
101         format(' ','3D histogram on triangular grid not supported ',
     1       'yet'/' ','Hit return to continue...')
            read(luttyi,106,end=300,err=300) buffer
         else
            call ship_surface(arrdat(arrind(menu_to_arr(1))), ni, nj,
     1                     data_name(1), data_units(1),
     2                     data_nickname(1), is_data_log(1),
     3                     .FALSE.)
         endif
         return
      endif
c
200   call clear_screen
      write(luttyo,errfmt)     ! print any error messages
      errfmt = '('' '')'  ! clear them
      call write_data_menu('3D surface 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 binary file'')'
         call beep_terminal
      else if( data_is_trigrid .and. default_hidden_method .eq.
     1         HIDDEN_LINES_HISTOGRAM ) then
         errfmt =
     1   '('' '',''3D histogram on triangular grid not supported'')'
         call beep_terminal
      else if( is_data_1d(answer) ) then
         errfmt = '('' '',''Cannot perform surface plot of 1D data'')'
         call beep_terminal
      else
         call ship_surface(arrdat(arrind(menu_to_arr(answer))), ni, nj,
     1                     data_name(answer), data_units(answer),
     2                     data_nickname(answer), is_data_log(answer),
     3                     .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_surface(data_array, isize, jsize, label,
     1                         units, nickname, llog, 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
      logical lsave, llog
      integer isize, jsize
      real data_array(isize, jsize)
      character*(*) label, units, nickname
c
c     locals
c
      logical do_plot_title, do_abs, do_log, fake_ln03
      logical loutside
      character*50 horiz_label, vert_label, xyz_labels(3)
      character*20 horiz_units, vert_units, xyz_units(3)
      integer ilen, ulen, string_len
      integer ni_subset, nj_subset
      integer z_minor_tic_freq
      integer exponent(3)
      integer idum1, idum2
      real real_vec(10)
      real good_min(3), good_max(3), minor_tic(3), major_tic_freq(3)
      real delta
      real meshh(COMP_ARRAYSIZE), meshv(COMP_ARRAYSIZE)
      real rdum1, rdum2, rdum3
      logical tri_clip(MAX_TRI)
CC
CC BEGIN MENU
      real hmin_full, hmax_full, vmin_full, vmax_full,
     1     zmin_full, zmax_full
      real hmin_plot, hmax_plot, vmin_plot, vmax_plot,
     1     zmin_plot, zmax_plot
      real prev_zmin_plot, prev_zmax_plot, prev_z_tic_major
      real h_tic_major, v_tic_major, z_tic_major
      real h1, h2, zclip1, zclip2
      integer z_exponent, prev_z_minor_tic_freq
      integer iuser_scale, prev_iuser_scale,
     1        prev_z_exponent
      integer item_number, int_vec(2), ierror
      character*80 string
      character*50 z_label
      character*20 z_units
      character*3 str_abs, str_swap
      character*6 str_log
      logical eof, abort, is_graphics_terminal, prev_abs, do_uscale,
     1        prev_set, lopen, delayed_message, lannotate
CC END MENU
c
      character*1 cdum
CC
CC    BEGIN ORBIT
      real x1, y1, arg
      logical do_orbit
CC END ORBIT
c
      logical lfirst
      data lfirst/.TRUE./
      save lfirst
c
      lopen = .FALSE.
      lannotate = .FALSE.
      delayed_message = .FALSE.
      cdum = ' '
      ulen = string_len(pm_length_unit)
      prev_set = .FALSE.
c
      if( .not. lsave ) then
         plot_title = ' '
      endif
c
      if( .NOT. lsave .and.
     1    (threed_lfirst .or. l_autoscale .or.
     2    (prev_plane_type .ne. plane_type) ) ) then
c        Actually, we should reset bounds every time we read in a new
c        binary file.
         threed_lfirst = .FALSE.
         threed_h_min = xyzlimits(1,idisti)
         threed_h_max = xyzlimits(2,idisti)
         call get_default_limits(threed_h_min, threed_h_max,
     1      rdum1, rdum2, threed_h_tic_major, threed_h_minor_tic_freq,
     2      idum1, .FALSE.)
         threed_h_tic_major = threed_h_tic_major*10.0**idum1
         threed_v_min = xyzlimits(1,idistj)
         threed_v_max = xyzlimits(2,idistj)
         call get_default_limits(threed_v_min, threed_v_max,
     1      rdum1, rdum2, threed_v_tic_major, threed_v_minor_tic_freq,
     2      idum1, .FALSE.)
         threed_v_tic_major = threed_v_tic_major*10.0**idum1
c        Now account for user defined min/max
         if( l_x_axis_min ) threed_h_min = x_axis_min*pm_scale_plot2cm
         if( l_x_axis_max ) threed_h_max = x_axis_max*pm_scale_plot2cm
         if( l_y_axis_min ) threed_v_min = y_axis_min*pm_scale_plot2cm
         if( l_y_axis_max ) threed_v_max = y_axis_max*pm_scale_plot2cm
c
         if( binfiletype .eq. DT_FAIM .or.
     1       binfiletype .eq. DT_ASCII ) then
            threed_swapxy  = .FALSE.
         else
            threed_swapxy  = .TRUE.
         endif
c
         if( l_threed_view ) then
            threed_eye_position(1) = threed_view(1)
            threed_eye_position(2) = threed_view(2)
            threed_eye_position(3) = threed_view(3)
         else
            threed_eye_position(1) =   8.0
            threed_eye_position(2) =  45.0
            if( threed_swapxy ) then
               threed_eye_position(3) = -45.0
            else
               threed_eye_position(3) = -135.0
            endif
         endif
      endif
c
c     Define axis labeling
      if( lsave ) then
            horiz_label = threed_h_label
            horiz_units = threed_h_units
            vert_label  = threed_v_label
            vert_units  = threed_v_units
      else if (lfirst .or. prev_plane_type .ne. plane_type) then
         lfirst = .FALSE.
         if( binfiletype .eq. DT_ASCII .or.
     1       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
      if( lsave ) then
         z_label = threed_z_label
         z_units = threed_z_units
      else
         z_label = label
         z_units = units
      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*threed_h_min
      hmax_plot = pm_scale_cm2plot*threed_h_max
      vmin_plot = pm_scale_cm2plot*threed_v_min
      vmax_plot = pm_scale_cm2plot*threed_v_max
      h_tic_major = pm_scale_cm2plot*threed_h_tic_major
      v_tic_major = pm_scale_cm2plot*threed_v_tic_major
c
      zclip1 = -PM_MAX_FLOAT
      zclip2 =  PM_MAX_FLOAT
c     Find subset of array which user specified.
c     Interpolate, if necessary, to get exact user specified bounds.
c     Swap i and j indices, if requested (many simulator coordinate systems
c     are left handed).
c     It also returns interpolated mesh arrays, scaled by pm_scale_cm2plot.
      if( data_is_trigrid ) then
         call compute_clip_of_triarray(data_array, work,
     1                             cord, nop,
     1                             meshh, meshv, tri_clip,
     1                             threed_h_min, threed_h_max,
     1                             threed_v_min, threed_v_max,
     1                             tri_npts, tri_ne,
     1                             threed_swapxy,
     1                             loutside)
      else
         call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
      endif
      if( loutside ) then
         string = 
     1      'No data selected: plot bounds are outside all data values'
         delayed_message = .TRUE.
      endif
      if( lsave ) then
         zmin_plot        = threed_z_min
         zmax_plot        = threed_z_max
         z_tic_major      = threed_z_tic_major
         z_minor_tic_freq = threed_z_minor_tic_freq
         z_exponent       = threed_z_exponent
         do_abs = threed_z_abs
         do_log = threed_z_log
         do_uscale = .FALSE.
         iuser_scale = 0
         if( data_is_trigrid ) then
            call abs_and_log_tri3d(work, tri_npts, do_log, do_abs)
            call recompute_minmax_tri3d(work, tri_npts,
     1           do_log, do_abs, iuser_scale, do_uscale,
     1           zmin_full, zmax_full, rdum1, rdum2,
     1           rdum3, idum1, idum2)
         else
            call abs_and_log_3d(work, ni_subset, nj_subset,
     1              do_log, do_abs)
            call recompute_minmax_3d(work, ni_subset, nj_subset,
     1           do_log, do_abs, iuser_scale, do_uscale,
     2           zmin_full, zmax_full, rdum1, rdum2,
     3           rdum3, idum1, idum2)
         endif
         if( idum2 .ne. threed_z_exponent ) then
            do_uscale = .TRUE.
            iuser_scale = threed_z_exponent
         endif
      else
         if( binfiletype .eq. DT_FAIM .or.
     1       binfiletype .eq. DT_ASCII ) then
            threed_swapxy  = .FALSE.
         else
            threed_swapxy  = .TRUE.
         endif
         iuser_scale = 0
         do_abs  = .FALSE.
         if( llog ) then
            do_abs  = .TRUE.
            do_log  = .TRUE.
         else
            do_log  = .FALSE.
         endif
         do_uscale = .FALSE.
         if( data_is_trigrid ) then
            call abs_and_log_tri3d(work, tri_npts, do_log, do_abs)
            call recompute_minmax_tri3d(work, tri_npts,
     1           do_log, do_abs, iuser_scale, do_uscale,
     1           zmin_full, zmax_full, zmin_plot, zmax_plot,
     1           z_tic_major, z_minor_tic_freq, z_exponent)
         else
            call abs_and_log_3d(work, ni_subset, nj_subset,
     1           do_log, do_abs)
            call recompute_minmax_3d(work, ni_subset, nj_subset,
     1           do_log, do_abs, iuser_scale, do_uscale,
     2           zmin_full, zmax_full, zmin_plot, zmax_plot,
     3           z_tic_major, z_minor_tic_freq, z_exponent)
         endif
         if( do_log ) then
            z_minor_tic_freq = 1
            z_tic_major = 1.0
         endif
      endif
c
      prev_iuser_scale = 0
      prev_abs = .FALSE.
c
      if( do_abs ) then
         str_abs = 'yes'
      else
         str_abs = ' no'
      endif
      if( do_log ) then
         str_log = '   log'
      else
         str_log = 'linear'
      endif
      if( threed_swapxy ) then
         str_swap = 'yes'
      else
         str_swap = 'no '
      endif
      do_plot_title = .TRUE.
c
c     If batch or quiet processing, perform plot 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
         call copy_axis_to_save(PLOT_3D_COMMAND,
     1        hmin_plot, hmax_plot,
     2        0, h_tic_major,
     3        threed_h_minor_tic_freq, horiz_units,
     4        horiz_label, .FALSE., .FALSE., .TRUE.,
     5        vmin_plot, vmax_plot,
     6        0, v_tic_major,
     7        threed_v_minor_tic_freq, vert_units,
     8        vert_label, .FALSE., .FALSE., .TRUE.,
     9        zmin_plot, zmax_plot,
     A        z_exponent, z_tic_major,
     B        z_minor_tic_freq, z_units,
     C        z_label, do_log, do_abs, .TRUE.)
c
c        Stuff data information into save structures
         call copy_3d_to_save(nickname)
c
c        Set up stuff for calling SURF routines.
         if( do_log ) then
            h1 = ALOG10(zmin_plot)
            h2 = ALOG10(zmax_plot)
         else
            h1 = zmin_plot
            h2 = zmax_plot
         endif
         call conv_pm3d_axes_to_surf(hmin_plot, vmin_plot, h1,
     1        hmax_plot, vmax_plot, h2,
     2        h_tic_major, v_tic_major, z_tic_major,
     3        threed_h_minor_tic_freq, threed_v_minor_tic_freq,
     4        z_minor_tic_freq,
     5        0, 0, z_exponent,
     6        horiz_label, horiz_units, vert_label, vert_units,
     7        z_label, z_units,
     8        good_min, good_max, minor_tic, major_tic_freq, exponent,
     9        xyz_labels, xyz_units,
     A        threed_swapxy)
c
c        Clip data to Z axis limits
         if( do_log ) then
            zclip1 = ALOG10(zmin_plot)
            zclip2 = ALOG10(zmax_plot)
         elseif( do_abs ) then
            zclip1 = 0.0
            zclip2 = ABS(zmax_plot)*10.0**z_exponent
         else
            zclip1 = zmin_plot*10.0**z_exponent
            zclip2 = zmax_plot*10.0**z_exponent
         endif
         if( data_is_trigrid ) then
            call clip_and_normalize_array_tri3d(work, tri_npts,
     1                      zclip1, zclip2, z_exponent)
         else
            call clip_and_normalize_array_3d(work, ni_subset, nj_subset,
     1                      zclip1, zclip2, z_exponent)
         endif
c
         call plot_surface(work, ni_subset, nj_subset,
     1        meshh, meshv,
     1        cord, nop, tri_clip, data_is_trigrid, tri_npts, tri_ne,
     1        threed_eye_position(1), threed_eye_position(2),
     1        threed_eye_position(3),
     1        good_min, good_max, minor_tic, major_tic_freq, exponent,
     1        xyz_labels, xyz_units,
     1        do_log, do_plot_title, lsave, lopen, lannotate, ierror)
         plot_device = previous_plot_device
         return
      endif
c
c     Main return point for menu
c
4000  continue
c     Write out menu form:
      call clear_screen
      call display_string(1,1 ,'1)  Title:',10, 0)
c
      call display_string(3 ,1 ,'HORIZONTAL AXIS PARAMETERS:',27,0)
      call display_string(4 ,1 ,'Data range:',11,0)
      call display_string(4 ,24,' to ',4,0)
      call display_string(6 ,1 ,'2)  Plot range:',15,0)
      call display_string(6 ,24,' to ',4,0)
      call display_string(7 ,1 ,'3)  Major tics at:',18,0)
      call display_string(8 ,1 ,'4)  Minor tic freq:',19,0)
      call display_string(9 ,1 ,'5)  Label: ',11, 0)
      call display_string(10,1 ,'6)  Units: ',11, 0)
c
      call display_string(12,1 ,'VERTICAL AXIS PARAMETERS:',25,0)
      call display_string(13,1 ,'Data range:',11,0)
      call display_string(13,24,' to',3,0)
      call display_string(15,1 ,'7)  Plot range:',15,0)
      call display_string(15,24,' to',3,0)
      call display_string(16,1 ,'8)  Major tics at:',18,0)
      call display_string(17,1 ,'9)  Minor tic freq:',19,0)
      call display_string(18,1 ,'10) Label: ',11, 0)
      call display_string(19,1 ,'11) Units: ',11, 0)
c
      call display_string(3 ,42,'Z AXIS PARAMETERS:',18,0)
      call display_string(4 ,42,'Data range:               to',28,0)
      call display_string(6 ,42,'12) Plot range:',15,0)
      call display_string(6 ,24,' to',3,0)
      call display_string(7 ,42,'13) Major tics at:',18,0)
      call display_string(8 ,42,'14) Minor tic freq:',19,0)
      call display_string(9 ,42,'15) Scale factor:     10**',26,0)
      call display_string(10,42,'16) Linear/log:',15,0)
      call display_string(11,42,'17) Absolute value:',19,0)
      call display_string(12,42,
     1     '18) Recompute default axis parameters',37,0)
      call display_string(13,42,'19) Label: ',11, 0)
      call display_string(14,42,'20) Units: ',11, 0)
c
      call display_string(16,42,'21) Eye position: ',18, 0)
      call display_string(16,61,'22) Swap X/Y: ',14, 0)
      call display_string(16,61,'22) Flip Y axis: ',17, 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, threed_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, 20, 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, threed_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, 20, 1)
c
c     Z data min/max
      call display_real  (4 , 58, zmin_full, 1)
      call display_real  (4 , 71, zmax_full, 1)
      call display_real  (6 , 58, zmin_plot, 1)
      call display_real  (6 , 71, zmax_plot, 1)
      if( do_log ) then
         call display_string(7, 62, '       NA', 9, 1)
         call display_string(8, 66, '   NA', 5, 1)
      else
         call display_real  (7 , 62, z_tic_major, 1)
         call display_int   (8 , 66, z_minor_tic_freq, 5)
      endif
      call display_int   (9 , 68, z_exponent, 3)
      call display_string(10 , 65, str_log, 6, 1)
      call display_string(11, 68, str_abs, 3, 1)
      call display_string(13, 53, z_label, 29, 1)
      call display_string(14, 53, z_units, 20, 1)
c
      call display_real  (17, 46, threed_eye_position(1), 1)
      call display_real  (18, 46, threed_eye_position(2), 1)
      call display_real  (19, 46, threed_eye_position(3), 1)
c
      call display_string(16, 78, str_swap, 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 number of item, or '//
     2 'Annotate Plot Hardcopy Keepwin Orbit Repaint 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
         call copy_axis_to_save(PLOT_3D_COMMAND,
     1        hmin_plot, hmax_plot,
     2        0, h_tic_major,
     3        threed_h_minor_tic_freq, horiz_units,
     4        horiz_label, .FALSE., .FALSE., .TRUE.,
     5        vmin_plot, vmax_plot,
     6        0, v_tic_major,
     7        threed_v_minor_tic_freq, vert_units,
     8        vert_label, .FALSE., .FALSE., .TRUE.,
     9        zmin_plot, zmax_plot,
     A        z_exponent, z_tic_major,
     B        z_minor_tic_freq, z_units,
     C        z_label, do_log, do_abs, .TRUE.)
c
c        Stuff data information into save structures
         call copy_3d_to_save(nickname)
c
c        Set up stuff for calling SURF routines.
         if( do_log ) then
            h1 = ALOG10(zmin_plot)
            h2 = ALOG10(zmax_plot)
         else
            h1 = zmin_plot
            h2 = zmax_plot
         endif
         call conv_pm3d_axes_to_surf(hmin_plot, vmin_plot, h1,
     1        hmax_plot, vmax_plot, h2,
     2        h_tic_major, v_tic_major, z_tic_major,
     3        threed_h_minor_tic_freq, threed_v_minor_tic_freq,
     4        z_minor_tic_freq,
     5        0, 0, z_exponent,
     6        horiz_label, horiz_units, vert_label, vert_units,
     7        z_label, z_units,
     8        good_min, good_max, minor_tic, major_tic_freq, exponent,
     9        xyz_labels, xyz_units,
     A        threed_swapxy)
c
c        Clip data to Z axis limits
         if( do_log ) then
            zclip1 = ALOG10(zmin_plot)
            zclip2 = ALOG10(zmax_plot)
         elseif( do_abs) then
            zclip1 = 0.0
            zclip2 = ABS(zmax_plot)*10.0**z_exponent
         else
            zclip1 = zmin_plot*10.0**z_exponent
            zclip2 = zmax_plot*10.0**z_exponent
         endif
         if( data_is_trigrid ) then
            call compute_clip_of_triarray(data_array, work,
     1                             cord, nop,
     1                             meshh, meshv, tri_clip,
     1                             threed_h_min, threed_h_max,
     1                             threed_v_min, threed_v_max,
     1                             tri_npts, tri_ne,
     1                             threed_swapxy,
     1                             loutside)
            call abs_and_log_tri3d(work, tri_npts,
     1                       do_log, do_abs)
            call clip_and_normalize_array_tri3d(work, tri_npts,
     1                      zclip1, zclip2, z_exponent)
         else
            call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
            call abs_and_log_3d(work, ni_subset, nj_subset,
     1                       do_log, do_abs)
            call clip_and_normalize_array_3d(work, ni_subset, nj_subset,
     1                      zclip1, zclip2, z_exponent)
         endif
c
         call plot_surface(work, ni_subset, nj_subset,
     1        meshh, meshv,
     1        cord, nop, tri_clip, data_is_trigrid, tri_npts, tri_ne,
     1        threed_eye_position(1), threed_eye_position(2),
     1        threed_eye_position(3),
     1        good_min, good_max, minor_tic, major_tic_freq, exponent,
     1        xyz_labels, xyz_units,
     1        do_log, do_plot_title, lsave, lopen, lannotate, ierror)
         if( ierror .eq. PLT$NORMAL ) then
            lopen = .TRUE.
            if( loutside ) then
               call display_message(
     1      'No data selected: plot bounds are outside all data values')
            endif
         else
            call plt$get_error_message(ierror,string)
            delayed_message = .TRUE.
         endif
         goto 4000
c
      else if( string(1:1) .eq. 'H' ) then
         fake_ln03 = .TRUE.
         call inquire_hardcopy_menu(fake_ln03, abort)
         if( abort ) goto 4000
         call display_prompt('Making hardcopy')
c
c        BRANCH IN POINT FOR BATCH PROCESSING
c
2500     continue
         previous_plot_device = plot_device
         plot_device = hardcopy_device
c
         call copy_axis_to_save(PLOT_3D_COMMAND,
     1        hmin_plot, hmax_plot,
     2        0, h_tic_major,
     3        threed_h_minor_tic_freq, horiz_units,
     4        horiz_label, .FALSE., .FALSE., .TRUE.,
     5        vmin_plot, vmax_plot,
     6        0, v_tic_major,
     7        threed_v_minor_tic_freq, vert_units,
     8        vert_label, .FALSE., .FALSE., .TRUE.,
     9        zmin_plot, zmax_plot,
     A        z_exponent, z_tic_major,
     B        z_minor_tic_freq, z_units,
     C        z_label, do_log, do_abs, .TRUE.)
c
c        Stuff data information into save structures
         call copy_3d_to_save(nickname)
c
c        Set up stuff for calling SURF routines.
         if( do_log ) then
            h1 = ALOG10(zmin_plot)
            h2 = ALOG10(zmax_plot)
         else
            h1 = zmin_plot
            h2 = zmax_plot
         endif
         call conv_pm3d_axes_to_surf(hmin_plot, vmin_plot, h1,
     1        hmax_plot, vmax_plot, h2,
     2        h_tic_major, v_tic_major, z_tic_major,
     3        threed_h_minor_tic_freq, threed_v_minor_tic_freq,
     4        z_minor_tic_freq,
     5        0, 0, z_exponent,
     6        horiz_label, horiz_units, vert_label, vert_units,
     7        z_label, z_units,
     8        good_min, good_max, minor_tic, major_tic_freq, exponent,
     9        xyz_labels, xyz_units,
     A        threed_swapxy)
c
c        Clip data to Z axis limits
         if( do_log ) then
            zclip1 = ALOG10(zmin_plot)
            zclip2 = ALOG10(zmax_plot)
         elseif( do_abs ) then
            zclip1 = 0.0
            zclip2 = ABS(zmax_plot)*10.0**z_exponent
         else
            zclip1 = zmin_plot*10.0**z_exponent
            zclip2 = zmax_plot*10.0**z_exponent
         endif
         if( data_is_trigrid ) then
            call compute_clip_of_triarray(data_array, work,
     1                             cord, nop,
     1                             meshh, meshv, tri_clip,
     1                             threed_h_min, threed_h_max,
     1                             threed_v_min, threed_v_max,
     1                             tri_npts, tri_ne,
     1                             threed_swapxy,
     1                             loutside)
            call abs_and_log_tri3d(work, tri_npts,
     1                       do_log, do_abs)
            call clip_and_normalize_array_tri3d(work, tri_npts,
     1                      zclip1, zclip2, z_exponent)
         else
            call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
            call abs_and_log_3d(work, ni_subset, nj_subset,
     1                       do_log, do_abs)
            call clip_and_normalize_array_3d(work, ni_subset, nj_subset,
     1                      zclip1, zclip2, z_exponent)
         endif
c
         call plot_surface(work, ni_subset, nj_subset,
     1        meshh, meshv,
     1        cord, nop, tri_clip, data_is_trigrid, tri_npts, tri_ne,
     1        threed_eye_position(1), threed_eye_position(2),
     1        threed_eye_position(3),
     1        good_min, good_max, minor_tic, major_tic_freq, exponent,
     1        xyz_labels, xyz_units,
     1        do_log, do_plot_title, lsave, lopen, lannotate, ierror)
         plot_device = previous_plot_device
         if( ierror .eq. PLT$NORMAL ) then
            ilen = string_len(plot_file_name)
            write(string,200,err=201) plot_file_name(1:ilen)
200         format('Plot file ',a' created')
201         continue
         else 
            call plt$get_error_message(ierror,string)
         endif
         delayed_message = .TRUE.
         goto 4000
c
      else if( string(1:1) .eq. 'O' ) then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: cannot orbit on this plot device')
               goto 2000
            else
c              Orbit loop:
c              Draw initial bounding box
               call surf_boundingbox(current_window_id,
     1                      threed_eye_position(1),
     2                      threed_eye_position(2),
     3                      threed_eye_position(3),
     4                      .TRUE.)
               call display_prompt(
     1            'Click on window with left button to rotate ' //
     1            'bounding box. '// PM_CANCEL_BUTTON //
     1            ' button re-plots.')
3000           continue
               call get_orbit_locator(current_window_id,
     1               x1, y1, arg, do_orbit)
               if( do_orbit ) then
                  call orbit_modify_eye(x1, y1, arg,
     1               threed_eye_position)
c                 Draw bounding box
                  call surf_boundingbox(current_window_id,
     1                      threed_eye_position(1),
     2                      threed_eye_position(2),
     3                      threed_eye_position(3),
     4                      .FALSE.)
c                 Redisplay eye position
                  call display_real  (17, 46, threed_eye_position(1), 1)
                  call display_real  (18, 46, threed_eye_position(2), 1)
                  call display_real  (19, 46, threed_eye_position(3), 1)
                  goto 3000
               endif
c
                  call copy_axis_to_save(PLOT_3D_COMMAND,
     1               hmin_plot, hmax_plot,
     2               0, h_tic_major,
     3               threed_h_minor_tic_freq, horiz_units,
     4               horiz_label, .FALSE., .FALSE., .TRUE.,
     5               vmin_plot, vmax_plot,
     6               0, v_tic_major,
     7               threed_v_minor_tic_freq, vert_units,
     8               vert_label, .FALSE., .FALSE., .TRUE.,
     9               zmin_plot, zmax_plot,
     A               z_exponent, z_tic_major,
     B               z_minor_tic_freq, z_units,
     C               z_label, do_log, do_abs, .TRUE.)
c
c                 Stuff data information into save structures
                  call copy_3d_to_save(nickname)
c
c                 Set up stuff for calling SURF routines.
                  if( do_log ) then
                     h1 = ALOG10(zmin_plot)
                     h2 = ALOG10(zmax_plot)
                  else
                     h1 = zmin_plot
                     h2 = zmax_plot
                  endif
                  call conv_pm3d_axes_to_surf(hmin_plot, vmin_plot, h1,
     1               hmax_plot, vmax_plot, h2,
     2               h_tic_major, v_tic_major, z_tic_major,
     3               threed_h_minor_tic_freq, threed_v_minor_tic_freq,
     4               z_minor_tic_freq,
     5               0, 0, z_exponent,
     6               horiz_label, horiz_units, vert_label, vert_units,
     7               z_label, z_units,
     8               good_min, good_max, minor_tic, major_tic_freq,
     9               exponent, xyz_labels, xyz_units,
     A               threed_swapxy)
c
c                 Clip data to Z axis limits
                  if( do_log ) then
                     zclip1 = ALOG10(zmin_plot)
                     zclip2 = ALOG10(zmax_plot)
                  elseif( do_abs) then
                     zclip1 = 0.0
                     zclip2 = ABS(zmax_plot)*10.0**z_exponent
                  else
                     zclip1 = zmin_plot*10.0**z_exponent
                     zclip2 = zmax_plot*10.0**z_exponent
                  endif
                  if( data_is_trigrid ) then
                     call compute_clip_of_triarray(data_array, work,
     1                             cord, nop,
     1                             meshh, meshv, tri_clip,
     1                             threed_h_min, threed_h_max,
     1                             threed_v_min, threed_v_max,
     1                             tri_npts, tri_ne,
     1                             threed_swapxy,
     1                             loutside)
                     call abs_and_log_tri3d(work, tri_npts,
     1                       do_log, do_abs)
                     call clip_and_normalize_array_tri3d(work,
     1                      tri_npts, zclip1, zclip2, z_exponent)
                  else
                     call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                     call abs_and_log_3d(work, ni_subset, nj_subset,
     1                       do_log, do_abs)
                     call clip_and_normalize_array_3d(work,
     1                      ni_subset, nj_subset,
     2                      zclip1, zclip2, z_exponent)
                  endif
c
                  call plot_surface(work, ni_subset, nj_subset,
     1               meshh, meshv,
     1               cord, nop, tri_clip, data_is_trigrid,
     1               tri_npts, tri_ne,
     1               threed_eye_position(1), threed_eye_position(2),
     1               threed_eye_position(3),
     1               good_min, good_max,
     1               minor_tic, major_tic_freq, exponent,
     1               xyz_labels, xyz_units,
     1               do_log, do_plot_title, lsave, lopen, lannotate,
     1               ierror)
                  if( ierror .eq. PLT$NORMAL ) then
                     lopen = .TRUE.
                  endif
            endif
         else
            call display_message(
     1         'Plot must be displayed before you can rotate it')
         endif
         goto 2000
c
      else if( string(1:1) .eq. 'A' ) then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: cannot annotate on this graphics device')
            else
               call anno_annotate_graph(current_window_id)
               lannotate = .TRUE.
               goto 4000
            endif
         else
            call display_message(
     1         'Plot must be displayed before you can annotate it')
         endif
         goto 2000
c     Special case to catch E or Q (which are "legal" floating point
c     numbers!)
      else if( string(1:1) .eq. 'E' .or. string(1:1) .eq. 'Q') then
         goto 5000
      else if( string(1:1) .eq. 'K') then
         if( lopen ) then
            if( is_graphics_terminal(current_window_id) ) then
               call display_message(
     1           'Sorry: no window functions on this graphics device')
               goto 2000
            else
c              Set "save this window" flag
c              (only for workstations and when using retained window style)
               call retain_current_window
            endif
         else
            call display_message(
     1      'Plot must be displayed in a window before you can keep it')
         endif
         goto 2000
      else if( string(1:1) .eq. '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 if( MAX(real_vec(1),real_vec(2)) .lt.
     1               pm_scale_cm2plot*xyzlimits(1,idisti) .or.
     2               MIN(real_vec(1),real_vec(2)) .gt.
     3               pm_scale_cm2plot*xyzlimits(2,idisti))then
               call display_message(
     1            'Requested bounds are outside all data values')
            else
               hmin_plot = real_vec(1)
               hmax_plot = real_vec(2)
               threed_h_min = pm_scale_plot2cm*hmin_plot
               threed_h_max = pm_scale_plot2cm*hmax_plot
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
               endif
               call display_real(6 , 15, hmin_plot, 1)
               call display_real(6 , 28, hmax_plot, 1)
c              Check tic mark parameters for consistency with new min/max
               call get_default_limits(threed_h_min, threed_h_max,
     1            rdum1, rdum2, threed_h_tic_major,
     2            threed_h_minor_tic_freq, idum1, .FALSE.)
                  threed_h_tic_major = threed_h_tic_major*10.0**idum1
               h_tic_major = pm_scale_cm2plot*threed_h_tic_major
               call display_real( 7, 23, h_tic_major, 1)
               call display_int ( 8, 27, threed_h_minor_tic_freq, 5)
            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 .or.
     2          real_vec(1) .lt. delta/100 ) then
               call display_message('Illegal value')
            else
               h_tic_major = real_vec(1)
               threed_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) = threed_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
               threed_h_minor_tic_freq = int_vec(1)
               call display_int (8 , 27, threed_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 if( MAX(real_vec(1),real_vec(2)) .lt.
     1               pm_scale_cm2plot*xyzlimits(1,idistj) .or.
     2               MIN(real_vec(1),real_vec(2)) .gt.
     3               pm_scale_cm2plot*xyzlimits(2,idistj))then
               call display_message(
     1            'Requested bounds are outside all data values')
            else
               vmin_plot = real_vec(1)
               vmax_plot = real_vec(2)
               threed_v_min = pm_scale_plot2cm*vmin_plot
               threed_v_max = pm_scale_plot2cm*vmax_plot
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
              endif
c
               call display_real(15, 15, vmin_plot, 1)
               call display_real(15, 28, vmax_plot, 1)
c              Check tic mark parameters for consistency with new min/max
               call get_default_limits(threed_v_min, threed_v_max,
     1            rdum1, rdum2, threed_v_tic_major,
     2            threed_v_minor_tic_freq, idum1, .FALSE.)
                  threed_v_tic_major = threed_v_tic_major*10.0**idum1
               v_tic_major = pm_scale_cm2plot*threed_v_tic_major
               call display_real(16, 23, v_tic_major, 1)
               call display_int (17, 27, threed_v_minor_tic_freq, 5)
            endif
c
         elseif ( item_number .eq. 8 ) then
c           Major axis 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 .or.
     2          real_vec(1) .lt. delta/100 ) then
                call display_message('Illegal value')
            else
               v_tic_major = real_vec(1)
               threed_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) = threed_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
               threed_v_minor_tic_freq = int_vec(1)
               call display_int (17, 27, threed_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           Z axis min/max
            real_vec(1) = zmin_plot
            real_vec(2) = zmax_plot
            call get_reals_with_prompt(real_vec, 2,
     1              'Enter min and max values for Z axis')
            if( real_vec(1) .ge. real_vec(2) ) then
               call display_message(
     1                 'Min value must be less than max')
            else if( do_log .and.
     1         (real_vec(1) .le. 0.0 .or. real_vec(2) .le. 0.0) ) then
               call display_message(
     1                 'Log limits must be greater than zero')
            else if( do_log .and.
     1         (real_vec(2)/MAX(1.0E-30,real_vec(1)) .lt. 10.0 ) ) then
               call display_message(
     1                 'Log limits must span an order of magnitude')
            else
               zmin_plot = real_vec(1)
               zmax_plot = real_vec(2)
               call display_real  (6 , 58, zmin_plot, 1)
               call display_real  (6 , 71, zmax_plot, 1)
               if( .NOT. do_log ) then
c                 Check tic mark parameters for consistency with new min/max
                  call get_default_limits(zmin_plot, zmax_plot,
     1               rdum1, rdum2, z_tic_major,
     2               z_minor_tic_freq, idum1, .FALSE.)
                  z_tic_major = z_tic_major*10.0**idum1
                  call display_real(7 , 62, z_tic_major, 1)
                  call display_int (8 , 66, z_minor_tic_freq, 5)
               endif
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
               endif
            endif
         elseif ( item_number .eq. 13 ) then
c           Major Z tic increment
            real_vec(1) = z_tic_major
            if( do_log ) then
               call display_message(
     1            'Not used when plotting on log scale')
            else
               call get_reals_with_prompt(real_vec, 1,
     1              'Enter major tic increment for Z axis')
               delta = zmax_plot - zmin_plot
               if( real_vec(1) .gt. delta .or.
     1             real_vec(1) .le. 0.0 .or.
     2             real_vec(1) .lt. delta/100 ) then
                   call display_message('Illegal value')
               else
                  z_tic_major = real_vec(1)
                  call display_real  (7 , 62, z_tic_major, 1)
               endif
            endif
         elseif ( item_number .eq. 14 ) then
c           Minor tic frequency
            int_vec(1) = z_minor_tic_freq
            if( do_log ) then
               call display_message(
     1            'Not used when plotting on log scale')
            else
               call get_ints_with_prompt(int_vec, 1,
     1              'Enter minor tic frequency for Z axis')
               if( int_vec(1) .lt. 1 ) then
                 call display_message('Illegal value')
               else
                  z_minor_tic_freq = int_vec(1)
                  call display_int (8 , 66, z_minor_tic_freq, 5)
               endif
            endif
         elseif ( item_number .eq. 15 ) then
c           Z axis scale factor
            int_vec(1) = z_exponent
            call get_ints_with_prompt(int_vec, 1,
     1   'Enter integer power of ten that will scale (divide) data')
            if( int_vec(1) .gt. 38 .or.
     1          int_vec(1) .lt. -38 ) then
               call display_message('Illegal value')
            else
               iuser_scale = int_vec(1)
               do_uscale = .TRUE.
               h1 = 10.0**iuser_scale/10.0**z_exponent
               zmin_plot = zmin_plot/h1
               zmax_plot = zmax_plot/h1
               if( .NOT. do_log ) then
                  z_tic_major = z_tic_major/h1
               endif
               z_exponent = iuser_scale
c              Re-display tic/scale data
               call display_real(6 , 58, zmin_plot, 1)
               call display_real(6 , 71, zmax_plot, 1)
               call display_int (9 , 68, z_exponent, 3)
               if( .NOT. do_log ) then
                  call display_real(7 , 62, z_tic_major, 1)
                  call display_int (8 , 66, z_minor_tic_freq, 5)
               endif
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
               endif
            endif
c
         elseif ( item_number .eq. 16 ) then
c           Logarithm
            do_log = .NOT. do_log
            if( do_log ) then
               prev_abs = do_abs
               do_abs = .TRUE.
               str_log = '   log'
               str_abs = 'yes'
c              Save previous values
               prev_z_exponent = z_exponent
               prev_iuser_scale = iuser_scale
               prev_zmin_plot = zmin_plot
               prev_zmax_plot = zmax_plot
               prev_z_tic_major  = z_tic_major
               prev_z_minor_tic_freq = z_minor_tic_freq
               prev_set = .TRUE.
               iuser_scale = 0
c
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
                  call recompute_minmax_tri3d(work, tri_npts,
     1               do_log, do_abs, iuser_scale, do_uscale,
     1               zmin_full, zmax_full, zmin_plot, zmax_plot,
     1               z_tic_major, z_minor_tic_freq, z_exponent)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
                  call recompute_minmax_3d(work, ni_subset, nj_subset,
     1               do_log, do_abs, iuser_scale, do_uscale,
     2               zmin_full, zmax_full, zmin_plot, zmax_plot,
     3               z_tic_major, z_minor_tic_freq, z_exponent)
               endif
               if( iuser_scale .ne. 0 ) z_exponent = iuser_scale
c              Mark Z axis tics as not applicable
               call display_string(7 , 62, '       NA', 9, 1)
               call display_string(8 , 66, '   NA', 5, 1)
               z_minor_tic_freq = 1
               z_tic_major = 1.0
            else
               do_abs = prev_abs
               str_log = 'linear'
               if( do_abs ) then
                  str_abs = 'yes'
               else
                  str_abs = ' no'
               endif
c              Recompute zmin_full/zmax_full
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
                  call recompute_minmax_tri3d(work, tri_npts,
     1               do_log, do_abs, iuser_scale, do_uscale,
     1               zmin_full, zmax_full, zmin_plot, zmax_plot,
     1               z_tic_major, z_minor_tic_freq, z_exponent)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
                  call recompute_minmax_3d(work, ni_subset, nj_subset,
     1               do_log, do_abs, iuser_scale, do_uscale,
     2               zmin_full, zmax_full, zmin_plot, zmax_plot,
     3               z_tic_major, z_minor_tic_freq, z_exponent)
               endif
c              Restore previous axis values
               if( prev_set ) then
                  iuser_scale = prev_iuser_scale
                  z_exponent = prev_z_exponent
                  zmin_plot = prev_zmin_plot
                  zmax_plot = prev_zmax_plot
                  z_tic_major = prev_z_tic_major
                  z_minor_tic_freq = prev_z_minor_tic_freq
               endif
               if( iuser_scale .ne. 0 ) z_exponent = iuser_scale
c              Re-display tic/scale data
               call display_real(7 ,62, z_tic_major, 1)
               call display_int (8 ,66, z_minor_tic_freq, 5)
            endif
            call display_real  (4 ,58, zmin_full, 1)
            call display_real  (4 ,71, zmax_full, 1)
            call display_real  (6 ,58, zmin_plot, 1)
            call display_real  (6 ,71, zmax_plot, 1)
            call display_int   (9 ,68, z_exponent, 3)
            call display_string(10,65, str_log, 6, 1)
            call display_string(11,68, str_abs, 3, 1)
c
         elseif ( item_number .eq. 17 ) then
            if( do_log .AND. do_abs ) then
               call display_message(
     1            'Absolute value required for log scale')
            else
               prev_abs = do_abs
               do_abs = .NOT. do_abs
               if( do_abs ) then
                  str_abs = 'yes'
               else
                  str_abs = ' no'
               endif
c              Recompute zmin_full/zmax_full
               if( data_is_trigrid ) then
                  call compute_clip_of_triarray(data_array, work,
     1                          cord, nop,
     1                          meshh, meshv, tri_clip,
     1                          threed_h_min, threed_h_max,
     1                          threed_v_min, threed_v_max,
     1                          tri_npts, tri_ne,
     1                          threed_swapxy,
     1                          loutside)
                  call abs_and_log_tri3d(work, tri_npts,
     1                          do_log, do_abs)
                  call recompute_minmax_tri3d(work, tri_npts,
     1               do_log, do_abs, iuser_scale, do_uscale,
     1               zmin_full, zmax_full, zmin_plot, zmax_plot,
     1               z_tic_major, z_minor_tic_freq, z_exponent)
               else
                  call interp_subset_of_array(data_array, work,
     1                             xyz(1,idisti), xyz(1,idistj),
     2                             meshh, meshv,
     3                             threed_h_min, threed_h_max,
     4                             threed_v_min, threed_v_max,
     5                             ni, nj, ni_subset, nj_subset,
     6                             COMP_ARRAYSIZE, threed_swapxy,
     7                             loutside)
                  call abs_and_log_3d(work, ni_subset, nj_subset,
     1                             do_log, do_abs)
                  call recompute_minmax_3d(work, ni_subset, nj_subset,
     1               do_log, do_abs, iuser_scale, do_uscale,
     2               zmin_full, zmax_full, zmin_plot, zmax_plot,
     3               z_tic_major, z_minor_tic_freq, z_exponent)
               endif
               if( do_uscale ) z_exponent = iuser_scale
            endif
c           Re-display tic/scale data
            call display_real(4 ,58, zmin_full, 1)
            call display_real(4 ,71, zmax_full, 1)
            call display_real(6 ,58, zmin_plot, 1)
            call display_real(6 ,71, zmax_plot, 1)
            call display_real(7 ,62, z_tic_major, 1)
            call display_int (8 ,66, z_minor_tic_freq, 5)
            call display_int (9 ,68, z_exponent, 3)
            call display_string(11,68, str_abs, 3, 1)
c
         elseif ( item_number .eq. 18 ) then
            call display_message('Not implemented yet')
         elseif ( item_number .eq. 19 ) then
            call get_string_with_prompt(z_label,
     1              'Enter label (may use mathmode)')
            call display_string(13, 53, z_label, 29, 1)
         elseif ( item_number .eq. 20 ) then
            call get_string_with_prompt(z_units,
     1              'Enter units (may use mathmode)')
            call display_string(14, 53, z_units, 20, 1)
         elseif ( item_number .eq. 21 ) then
            real_vec(1) = threed_eye_position(1)
            real_vec(2) = threed_eye_position(2)
            real_vec(3) = threed_eye_position(3)
            call get_reals_with_prompt(real_vec, 3,
     1      'Enter eye position (radius, polar angle, rotation angle)')
            threed_eye_position(1) = real_vec(1)
            threed_eye_position(2) = real_vec(2)
            threed_eye_position(3) = real_vec(3)
            call display_real  (17, 46, threed_eye_position(1), 1)
            call display_real  (18, 46, threed_eye_position(2), 1)
            call display_real  (19, 46, threed_eye_position(3), 1)
         elseif ( item_number .eq. 22 ) then
            threed_swapxy = .NOT. threed_swapxy
            if( threed_swapxy ) then
               str_swap = 'yes'
            else
               str_swap = 'no '
            endif
            call display_string(16, 78, str_swap, 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_surface(data_array, isize, jsize,
     1     h_vec, v_vec,
     1     cord, nop, tri_clip, data_is_trigrid, tri_npts, tri_ne,
     1     eye_radius, eye_phi, eye_theta,
     1     good_min, good_max, minor_tic, major_tic_freq, exponent,
     1     labels, units,
     1     do_log, do_header, add_notes_from_save, ws_window_open,
     1     add_notes_from_file, ierror)
c
      implicit none
c
c     Subroutine which actually calls graphics routines.
c
      logical do_header, do_log, add_notes_from_save,
     1        ws_window_open, add_notes_from_file
      character*(*) labels(3), units(3)
      integer isize, jsize, ierror
      integer exponent(3)
      real data_array(isize,jsize)
      real h_vec(isize), v_vec(jsize)
      real good_min(3), good_max(3), minor_tic(3), major_tic_freq(3)
      real eye_radius, eye_phi, eye_theta
c
      integer tri_npts, tri_ne
      integer nop(3,tri_ne)
      real    cord(2,tri_npts)
      logical data_is_trigrid, tri_clip(tri_ne)
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'io.cmn'
      include 'units.cmn'
      include 'defaults.cmn'
c
c     locals
c
      logical lerr, is_vtxxx
      logical minimum_specified(3), maximum_specified(3),
     1        axis_log_flags(3), axis_omit_flags(3),
     2        axis_nonumbering_flags(3)
      logical modify_eye, modify_axis
      integer axis_colors(3), selected_axis(3)
      integer ilen, ilen2, ilen3, itot, string_len, window_id,
     1        return_stat
      integer get_new_window_id
      real    axis_shortic_spacing(3), axis_longtic_frequency(3)
      real    axis_minimum(3), axis_maximum(3)
      real    theta, phi
      character*80 title1, title2, title3, label3, tmpstr
      character*80 label1, label2
      character*50 expstr
      character*1 response
      character*80 axis_label1, axis_label2
      character*20 ctmp
      character*8 time_string
c
c     open GKS
c
      window_id = get_new_window_id(PLT$3D_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 windows
         ierror = PLT$TOO_MANY_WINDOWS
         if( lbatch .or. lrestore ) then
            call plt$get_error_message(ierror,tmpstr)
            ilen = string_len(tmpstr)
            call plt$error_message('plot_surface', tmpstr(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,tmpstr)
            ilen = string_len(tmpstr)
            call plt$error_message('plot_surface', tmpstr(1:ilen))
         endif
         ierror = return_stat
         return
      endif
c
c     Activate plot device, clear it, start segment on VAXstations
      call activate_plot_device(window_id, .TRUE., .TRUE.)
      call pop_vs_window(window_id)
      call init_device_depend
c
c     Reset character sizing to defaults
      call GKS$SET_TEXT_EXPFAC(1.0)
      call GKS$SET_TEXT_SPACING(0.0)
c
c     Setup for call to SURF
c
      modify_eye  = .TRUE.
      modify_axis = .TRUE.
c
      axis_label1 = labels(1)
      axis_label2 = labels(2)
c
      if( l_x_axis_label ) axis_label1 = x_axis_label
      if( l_y_axis_label ) axis_label2 = y_axis_label
c
c     Plot axis such that it will always appear close to eye position
c
c     Reduce phi to range of 0 <= phi <= 180
      phi = MOD(eye_phi,360.0)
      if( phi .lt. 0.0 ) phi = phi + 360.0
      if( phi .gt. 180.0 ) then
         phi = 360.0 - phi
         theta = eye_theta + 180.0
      else
         theta = eye_theta
      endif
c     Reduce theta to range of 0 <= theta <= 360
      theta = MOD(theta,360.0)
      if( theta .lt. 0.0 ) theta = theta + 360.0
      if( theta .ge. 0.0 .and. theta .le. 90.0 ) then
         selected_axis(1) = 2
         selected_axis(2) = 2
         selected_axis(3) = 2
      else if( theta .gt. 90.0 .and. theta .le. 180.0 ) then
         selected_axis(1) = 2
         selected_axis(2) = 1
         selected_axis(3) = 3
      else if( theta .gt. 180.0 .and. theta .le. 270.0 ) then
         selected_axis(1) = 1
         selected_axis(2) = 1
         selected_axis(3) = 4
      else
         selected_axis(1) = 1
         selected_axis(2) = 2
         selected_axis(3) = 1
      endif    
      axis_log_flags(1) = .FALSE.
      axis_log_flags(2) = .FALSE.
      axis_log_flags(3) = .FALSE.
      axis_omit_flags(1) = .FALSE.
      axis_omit_flags(2) = .FALSE.
      axis_omit_flags(3) = .FALSE.
      axis_nonumbering_flags(1) = .FALSE.
      axis_nonumbering_flags(2) = .FALSE.
      axis_nonumbering_flags(3) = .FALSE.
c
c     X axis defaults
      minimum_specified(1) = .TRUE.
      maximum_specified(1) = .TRUE.
      axis_minimum(1)      = good_min(1)
      axis_maximum(1)      = good_max(1)
      axis_shortic_spacing(1)   = minor_tic(1)
      axis_longtic_frequency(1) = major_tic_freq(1)
      ilen  = string_len(axis_label1)
      ilen2 = string_len(units(1))
      label1 = ' '
      itot = 0
      if( ilen .gt. 0 ) then
         itot = ilen
         label1 = axis_label1(1:ilen)
      endif
      if( exponent(1) .ne. 0 ) then
         write(ctmp,200) exponent(1)
200      format('10**',i3,' ')
         if( ilen2 .gt. 0 ) then
            itot = ilen + ilen2 + 11
            label1(ilen+1:itot) = ' ' // left_unit_delim //
     1         ctmp(1:8) // units(1)(1:ilen2) // right_unit_delim
         else
            itot = ilen + 11
            label1(ilen+1:itot) = ' ' // left_unit_delim //
     1          ctmp(1:8) // right_unit_delim
         endif
      else
         if( ilen2 .gt. 0 ) then
            itot = ilen + ilen2 + 3
            label1(ilen+1:itot) = ' ' // left_unit_delim //
     1         units(1)(1:ilen2) // right_unit_delim
         endif
      endif
c
c     Y axis defaults
      minimum_specified(2) = .TRUE.
      maximum_specified(2) = .TRUE.
      axis_minimum(2)      = good_min(2)
      axis_maximum(2)      = good_max(2)
      axis_shortic_spacing(2)   = minor_tic(2)
      axis_longtic_frequency(2) = major_tic_freq(2)
      ilen  = string_len(axis_label2)
      ilen2 = string_len(units(2))
      itot = 0
      label2 = ' '
      if( ilen .gt. 0 ) then
         itot = ilen
         label2 = axis_label2(1:ilen)
      endif
      if( exponent(2) .ne. 0 ) then
         write(ctmp,200) exponent(2)
         if( ilen2 .gt. 0 ) then
            itot = ilen + ilen2 + 11
            label2(ilen+1:itot) = ' ' // left_unit_delim // ctmp(1:8) //
     1         units(2)(1:ilen2) // right_unit_delim
         else
            itot = ilen + 11
            label2(ilen+1:itot) = ' ' // left_unit_delim // ctmp(1:8) //
     1                            right_unit_delim
         endif
      else
         if( ilen2 .gt. 0 ) then
            itot = ilen + ilen2 + 3
            label2(ilen+1:itot) = ' ' // left_unit_delim //
     1         units(2)(1:ilen2) // right_unit_delim
         endif
      endif
c
c     Z axis
      minimum_specified(3) = .TRUE.
      maximum_specified(3) = .TRUE.
      axis_minimum(3)      = good_min(3)
      axis_maximum(3)      = good_max(3)
      axis_shortic_spacing(3)   = minor_tic(3)
      axis_longtic_frequency(3) = major_tic_freq(3)
c     Add remote exponent to label, if necessary
      if( l_z_axis_label ) then
         tmpstr = z_axis_label
      else
         if( do_log ) then
            tmpstr = 'Log ' // labels(3)
         else
            tmpstr = labels(3)
         endif
      endif
      ilen = string_len(tmpstr)
      ilen2 = string_len(units(3))     
      if( exponent(3) .ne. 0 ) then
         if( ilen2 .eq. 0 ) then
            write(expstr,220) exponent(3)
220         format(' (10^{',i3,'} )')
         else
            write(expstr,221) exponent(3), units(3)(1:ilen2)
221         format(' (10^{',i3,'}',1x,a,')')
         endif
         ilen3 = string_len(expstr)
         label3 = tmpstr(1:ilen) // expstr(1:ilen3)
      else
         if( ilen2 .eq. 0 ) then
            label3 = tmpstr(1:ilen)
         else
            label3 = tmpstr(1:ilen) // ' ' // left_unit_delim //
     1               units(3)(1:ilen2) // right_unit_delim
         endif
      endif
c
      axis_colors(1) = 1
      axis_colors(2) = 1
      axis_colors(3) = 1
c
c     Print auxiliary information
c
      if( do_header ) then
         title1 = plot_title
      else
         title1 = ' '
      endif
      call TIME(time_string)
      ilen = string_len(units(3))
      ilen2 = string_len(labels(3))
      if( ilen .eq. 0 ) then
         title2 = labels(3)(1:ilen2)
      else
         title2 = labels(3)(1:ilen2) // ' ' // left_unit_delim //
     1            units(3)(1:ilen) // right_unit_delim
      endif
      ilen = string_len(binfilename)
      title3 =  'Data source: ' // binfilename(1:ilen)
c
      if( .NOT. data_is_trigrid ) then
c        'Safe' values so that RTL will not complain of a zero row
c        index for adjustable size arrays.
         tri_npts = 1
         tri_ne   = 1
      endif
c
      call surf_surf(ws_id, title1, title2, title3,
     1             label1, label2, label3,
     2             h_vec, v_vec, data_array,
     2             isize, jsize, isize, jsize,
     2             data_is_trigrid, nop, tri_clip,
     2             tri_npts, tri_ne,
     2             minimum_specified, axis_minimum,
     2             maximum_specified, axis_maximum,
     3             modify_eye, eye_radius, phi, theta,
     4             modify_axis,
     5             selected_axis, axis_log_flags, axis_omit_flags,
     6             axis_nonumbering_flags, axis_shortic_spacing,
     7             axis_longtic_frequency, axis_colors,
     8             hidden_line_method, exponent(3))
c
c     Make sure all output gets onto screen...
      call deactivate_plot_device
c
      if( plot_device .eq. hardcopy_device ) then
c        Add in annotations to hardcopy plot
         if( ws_window_open ) then
            call anno_replay_annotations(current_window_id, window_id,
     1                 .TRUE.)
         else
            if( add_notes_from_save ) then
               call anno_restore_notes_from_save(window_id, lerr)
            endif
         endif
         call close_plot_device(window_id)
      else
         if( add_notes_from_save ) then
            call anno_restore_notes_from_save(window_id, lerr)
         else if ( add_notes_from_file ) then
c           Restore any previous annotations
            call read_savefile('POSTMINI.NOTE',.TRUE.,lerr)
            if( .NOT. lerr )
     1         call anno_restore_notes_from_save(window_id, lerr)
         endif
         if( is_vtxxx(window_id) ) then
            call move_cursor_bottom_of_screen
            call beep_terminal
            call get_char(response,'Hit return to continue ',0)
c           Restore default color map
            call reset_vt_colormap
         endif
      endif
c
      ierror = PLT$NORMAL
c
      end
c
      subroutine find_subset_of_array(hmesh, vmesh,
     1                          low_h, hi_h, low_v, hi_v,
     2                          ni, nj, ih1, ih2, iv1, iv2)
      implicit none
c
c     Find subset of array that lies within coordinates
c     low_h <= h <= hi_h  and low_v <= v <= hi_v
c     Mesh coordinates are stored in hmesh and vmesh
c     
c
      integer ni, nj
      real hmesh(ni), vmesh(nj)
      real low_h, hi_h, low_v, hi_v
      integer ih1, ih2, iv1, iv2
c
      real d1, d2
      integer getindex
c
      if( low_h .le. hmesh(1) ) then
         ih1 = 1
      else if ( low_h .ge. hmesh(ni) ) then
         ih1 = ni
      else
         ih1 = getindex(low_h, hmesh, ni)
c        Look at hmesh(ih1) and hmesh(ih1+1) to see which
c        is closer to low_h.
         if( ih1 .ne. ni ) then
            d1 = low_h - hmesh(ih1)
            d2 = hmesh(ih1+1) - low_h
            if( d1 .gt. d2 ) ih1 = ih1 + 1
         endif
      endif
c
      if( hi_h .le. hmesh(1) ) then
         ih2 = 1
      else if ( hi_h .ge. hmesh(ni) ) then
         ih2 = ni
      else
         ih2 = getindex(hi_h, hmesh, ni)
c        Look at hmesh(ih2) and hmesh(ih2+1) to see which
c        is closer to hi_h.
         if( ih2 .ne. ni ) then
            d1 = hi_h - hmesh(ih2)
            d2 = hmesh(ih2+1) - hi_h
            if( d1 .gt. d2 ) ih2 = ih2 + 1
         endif
      endif
c
      if( low_v .le. vmesh(1) ) then
         iv1 = 1
      else if ( low_v .ge. vmesh(nj) ) then
         iv1 = nj
      else
         iv1 = getindex(low_v, vmesh, nj)
c        Look at vmesh(iv1) and vmesh(iv1+1) to see which
c        is closer to low_v.
         if( iv1 .ne. nj ) then
            d1 = low_v - vmesh(iv1)
            d2 = vmesh(iv1+1) - low_v
            if( d1 .gt. d2 ) iv1 = iv1 + 1
         endif
      endif
      if( hi_v .le. vmesh(1) ) then
         iv2 = 1
      else if ( hi_v .ge. vmesh(nj) ) then
         iv2 = nj
      else
         iv2 = getindex(hi_v, vmesh, nj)
c        Look at vmesh(iv1) and vmesh(iv1+1) to see which
c        is closer to low_v.
         if( iv2 .ne. nj ) then
            d1 = hi_v - vmesh(iv2)
            d2 = vmesh(iv2+1) - hi_v
            if( d1 .gt. d2 ) iv2 = iv2 + 1
         endif
c
      endif
      return
      end
c
      
      subroutine get_surf_axis_limits(data_min, data_max,
     1           good_min, good_max, minor_tic, major_tic_freq,
     1           exponent, log_data)
      implicit none
c
      include 'sysdep.cmn'
c
c     Get limits for z axis for Nadim Khalil's 3D SURF subroutine
c
      real data_min, data_max, good_min, good_max
      real    minor_tic, major_tic_freq
      integer exponent
      logical log_data
c
c     locals
c
      real temp_max, log_small
      real tic_major
      integer minor_tic_freq
      real round_up_whole, round_down_whole
c
      log_small = ALOG10(PM_MIN_FLOAT)
c
c     Try to choose good default min/max values for plot
c     from data min/max.
c
c     First, check for funny limits. Try to avoid good_min = good_max.
      temp_max = data_max
      if( data_min .eq. data_max ) then
         if( data_max .eq. 0.0 ) then
            temp_max = data_max + 1.0
         else
            temp_max = 1.001*data_max
         endif
      endif
c
      if( log_data ) then
         if(  data_min .le. log_small) then
            good_min = -38.0
         else
            good_min = round_down_whole(data_min)
         endif
         good_max = round_up_whole(temp_max)
c        Set default plot limits depending on data range
         if( (good_max - good_min) .le. 6.0 ) then
            minor_tic      = 1.0
            major_tic_freq = 1.0
         else if( (good_max - good_min) .lt. 15.0 ) then
c           Round limits to power of two
            good_min = good_min - MOD(good_min,2.0)
            good_max = good_max + MOD(good_max,2.0)
            minor_tic      = 1.0
            major_tic_freq = 2.0
         else
c           Round minimum down to even power of 5.
c           Plot major tic marks every five decades.
c           Leave maximum value alone. Alternatively, we could
c           round maximum up to nearest 5th power, but this would
c           compress the plot more than necessary.
c           Also, many device quantities are near 10^{20}.
            good_min = good_min - MOD(good_min,5.0)
            minor_tic      = 1.0
            major_tic_freq = 5.0
         endif
      else
         call linear_autoscale(data_min, data_max, 1,
     1        good_min, good_max, exponent, tic_major, minor_tic_freq)
c        Convert autoscale's output to what SURF expects...
         minor_tic      = tic_major/FLOAT(minor_tic_freq)
         major_tic_freq = FLOAT(minor_tic_freq)
      endif
      return
      end
c
      subroutine interp_subset_of_array(data_array, work,
     1                             orig_meshh, orig_meshv,
     2                             meshh, meshv,
     3                             low_h, hi_h, low_v, hi_v,
     4                             ni, nj, ni_subset, nj_subset,
     5                             isize, swapxy, loutside)
      implicit none
c
      integer ni, nj, ni_subset, nj_subset, isize
      real    data_array(ni,nj), work(ni*nj)
      real    orig_meshh(isize), meshh(isize)
      real    orig_meshv(isize), meshv(isize)
      real    low_h, hi_h, low_v, hi_v
      logical swapxy, loutside
c
      include 'const.h'
      include 'units.cmn'
c
      real delx, dely
      integer i, j, ih1, ih2, iv1, iv2, kl, kt, itemp
      logical exact_left, exact_right, exact_top, exact_bottom
c
c     First, check for bounding box completely outside mesh
      if( MAX(low_h,hi_h) .lt. orig_meshh(1)  .or.
     1    MIN(low_h,hi_h) .gt. orig_meshh(ni) .or.
     1    MAX(low_v,hi_v) .lt. orig_meshv(1)  .or.
     1    MIN(low_v,hi_v) .gt. orig_meshv(nj) ) then
         loutside = .TRUE.
c        Return some innocuous values
         ni_subset = 1
         nj_subset = 1
         meshh(1) = 0.0
         meshv(1) = 0.0
         work(1)  = 0.0
         return
      else
         loutside = .FALSE.
      endif
c
      exact_left   = .FALSE.
      exact_right  = .FALSE.
      exact_top    = .FALSE.
      exact_bottom = .FALSE.
c
      ih1 = 1
      do 100 i=1,ni
         if( low_h .le. orig_meshh(i) ) then
            if( low_h .eq. orig_meshh(i) .or.
     1          i .eq. 1) exact_left = .TRUE.
            ih1 = i
            goto 101
         endif
100   continue
101   continue
c
      ih2 = ni
      do 102 i=ni,1,-1
         if( hi_h .ge. orig_meshh(i) ) then
            if( hi_h .eq. orig_meshh(i) .or.
     1          i .eq. ni ) exact_right = .TRUE.
            ih2 = i
            goto 103
         endif
102   continue
103   continue
c
      iv1 = 1
      do 104 i=1,nj
         if( low_v .le. orig_meshv(i) ) then
            if( low_v .eq. orig_meshv(i) .or.
     1          i .eq. 1 ) exact_top = .TRUE.
            iv1 = i
            goto 105
         endif
104   continue
105   continue
c
      iv2 = nj
      do 106 i=nj,1,-1
         if( hi_v .ge. orig_meshv(i) ) then
            if( hi_v .eq. orig_meshv(i) .or.
     1          i .eq. nj ) exact_bottom = .TRUE.
            iv2 = i
            goto 107
         endif
106   continue
107   continue
c
      ni_subset = ih2 - ih1 + 1
      if( .NOT. exact_left  ) then
         ni_subset = ni_subset + 1
         kl = 2
      else
         kl = 1
      endif
      if( .NOT. exact_right ) ni_subset = ni_subset + 1
c
      nj_subset = iv2 - iv1 + 1
      if( .NOT. exact_top    ) then
         nj_subset = nj_subset + 1
         kt = 2
      else
         kt = 1
      endif
      if( .NOT. exact_bottom ) nj_subset = nj_subset + 1
c
c     Copy subset arrays of mesh, data into interpolation storage.
c
      if( swapxy ) then
c
c        Swap x and y axes -- many simulators (e.g. MINIMOS) have positive
c        Y axis going down into depth (e.g. left handed coordinate system)
c
         meshv(1) = low_h*pm_scale_cm2plot
         meshv(ni_subset) = hi_h*pm_scale_cm2plot
         do i=ih1,ih2
            meshv(i-ih1+kl) = orig_meshh(i)*pm_scale_cm2plot
         end do
         meshh(1) = low_v*pm_scale_cm2plot
         meshh(nj_subset) = hi_v*pm_scale_cm2plot
         do j=iv1,iv2
            meshh(j-iv1+kt) = orig_meshv(j)*pm_scale_cm2plot
         end do
c
c        (linear) interpolate data onto work array.
c        SURF's coordinate system is right handed; MINIMOS's is left handed.
c        Exchange matrix indices to compensate.
c
c        Inner core
         do j=iv1,iv2
            do i=ih1,ih2
               work(j-iv1+kt+(i-ih1+kl-1)*nj_subset) = data_array(i,j)
            end do
         end do
c
c        Left
         if( .NOT. exact_left ) then
            delx = (low_h - orig_meshh(ih1-1))/
     1            (orig_meshh(ih1) - orig_meshh(ih1-1))
            do j=iv1,iv2
               work(j-iv1+kt) = data_array(ih1-1,j)  + 
     1            delx*(data_array(ih1,j) - data_array(ih1-1,j))
            end do
         endif
c
c        Right
         if( .NOT. exact_right ) then
            delx = (hi_h - orig_meshh(ih2))/
     1             (orig_meshh(ih2+1) - orig_meshh(ih2))
            do j=iv1,iv2
               work(j-iv1+kt+(ni_subset-1)*nj_subset) =
     1            data_array(ih2,j)  + 
     2            delx*(data_array(ih2+1,j) - data_array(ih2,j))
            end do
         endif
c
c        Top
         if( .NOT. exact_top ) then
            dely = (low_v - orig_meshv(iv1-1))/
     1            (orig_meshv(iv1) - orig_meshv(iv1-1))
            do j=ih1,ih2
               work(1+(j-ih1+kl-1)*nj_subset) = data_array(j,iv1-1)  + 
     1            dely*(data_array(j,iv1) - data_array(j,iv1-1))
            end do
         endif
c
c        Bottom
         if( .NOT. exact_bottom ) then
            dely = (hi_v - orig_meshv(iv2))/
     1            (orig_meshv(iv2+1) - orig_meshv(iv2))
            do j=ih1,ih2
               work(nj_subset+(j-ih1+kl-1)*nj_subset) =
     1            data_array(j,iv2)  + 
     2            dely*(data_array(j,iv2+1) - data_array(j,iv2))
            end do
         endif
c
c        Corners
         if( .NOT. (exact_left .or. exact_top) ) then
            delx = (low_h - orig_meshh(ih1-1))/
     1            (orig_meshh(ih1) - orig_meshh(ih1-1))
            dely = (low_v - orig_meshv(iv1-1))/
     1            (orig_meshv(iv1) - orig_meshv(iv1-1))
            work(1) = data_array(ih1,iv1)
     1          + delx*(data_array(ih1,iv1) - data_array(ih1-1,iv1))
     2          + dely*(data_array(ih1,iv1) - data_array(ih1,iv1-1))
         endif
         if( .NOT. (exact_left .or. exact_bottom) ) then
            delx = (low_h - orig_meshh(ih1-1))/
     1            (orig_meshh(ih1) - orig_meshh(ih1-1))
            dely = (hi_v - orig_meshv(iv2))/
     1            (orig_meshv(iv2+1) - orig_meshv(iv2))
            work(nj_subset) = data_array(ih1,iv2)
     1          + delx*(data_array(ih1,iv2) - data_array(ih1-1,iv2))
     2          + dely*(data_array(ih1,iv2+1) - data_array(ih1,iv2))
         endif
         if( .NOT. (exact_right .or. exact_top) ) then
            delx = (hi_h - orig_meshh(ih2))/
     1             (orig_meshh(ih2+1) - orig_meshh(ih2))
           dely = (low_v - orig_meshv(iv1-1))/
     1             (orig_meshv(iv1) - orig_meshv(iv1-1))
            work(1+(ni_subset-1)*nj_subset) = data_array(ih2,iv1)
     1          + delx*(data_array(ih2+1,iv1) - data_array(ih2,iv1))
     2          + dely*(data_array(ih2,iv1) - data_array(ih2,iv1-1))
         endif
         if( .NOT. (exact_right .or. exact_bottom) ) then
            delx = (hi_h - orig_meshh(ih2))/
     1            (orig_meshh(ih2+1) - orig_meshh(ih2))
            dely = (hi_v - orig_meshv(iv2))/
     1            (orig_meshv(iv2+1) - orig_meshv(iv2))
            work(nj_subset+(ni_subset-1)*nj_subset) =
     1          data_array(ih2,iv2)
     2          + delx*(data_array(ih2+1,iv2) - data_array(ih2,iv2))
     3          + dely*(data_array(ih2,iv2+1) - data_array(ih2,iv2))
         endif
c
         itemp = nj_subset
         nj_subset = ni_subset
         ni_subset = itemp
      else
c
c        Normal X and Y coordinate system
c
         meshh(1) = low_h*pm_scale_cm2plot
         meshh(ni_subset) = hi_h*pm_scale_cm2plot
         do i=ih1,ih2
            meshh(i-ih1+kl) = orig_meshh(i)*pm_scale_cm2plot
         end do
         meshv(1) = low_v*pm_scale_cm2plot
         meshv(nj_subset) = hi_v*pm_scale_cm2plot
         do j=iv1,iv2
            meshv(j-iv1+kt) = orig_meshv(j)*pm_scale_cm2plot
         end do
c
c        Inner core
         do j=iv1,iv2
            do i=ih1,ih2
               work(i-ih1+kl+(j-iv1+kt-1)*ni_subset) = data_array(i,j)
            end do
         end do
c
c        Left
         if( .NOT. exact_left ) then
            delx = (low_h - orig_meshh(ih1-1))/
     1            (orig_meshh(ih1) - orig_meshh(ih1-1))
            do j=iv1,iv2
               work(1+(j-iv1+kt-1)*ni_subset) = data_array(ih1-1,j)  + 
     1            delx*(data_array(ih1,j) - data_array(ih1-1,j))
            end do
         endif
c
c        Right
         if( .NOT. exact_right ) then
            delx = (hi_h - orig_meshh(ih2))/
     1             (orig_meshh(ih2+1) - orig_meshh(ih2))
            do j=iv1,iv2
               work(ni_subset+(j-iv1+kt-1)*ni_subset) =
     1            data_array(ih2,j)  + 
     2            delx*(data_array(ih2+1,j) - data_array(ih2,j))
            end do
         endif
c
c        Top
         if( .NOT. exact_top ) then
            dely = (low_v - orig_meshv(iv1-1))/
     1            (orig_meshv(iv1) - orig_meshv(iv1-1))
            do j=ih1,ih2
               work(j-ih1+kl) = data_array(j,iv1-1)  + 
     1            dely*(data_array(j,iv1) - data_array(j,iv1-1))
            end do
         endif
c
c        Bottom
         if( .NOT. exact_bottom ) then
            dely = (hi_v - orig_meshv(iv2))/
     1            (orig_meshv(iv2+1) - orig_meshv(iv2))
            do j=ih1,ih2
               work(j-ih1+kl+(nj_subset-1)*ni_subset) =
     1            data_array(j,iv2)  + 
     2            dely*(data_array(j,iv2+1) - data_array(j,iv2))
            end do
         endif
c
c        Corners
         if( .NOT. (exact_left .or. exact_top) ) then
            delx = (low_h - orig_meshh(ih1-1))/
     1            (orig_meshh(ih1) - orig_meshh(ih1-1))
            dely = (low_v - orig_meshv(iv1-1))/
     1            (orig_meshv(iv1) - orig_meshv(iv1-1))
             work(1) = data_array(ih1,iv1)
     1          + delx*(data_array(ih1,iv1) - data_array(ih1-1,iv1))
     2          + dely*(data_array(ih1,iv1) - data_array(ih1,iv1-1))
         endif
         if( .NOT. (exact_left .or. exact_bottom) ) then
            delx = (low_h - orig_meshh(ih1-1))/
     1            (orig_meshh(ih1) - orig_meshh(ih1-1))
            dely = (hi_v - orig_meshv(iv2))/
     1            (orig_meshv(iv2+1) - orig_meshv(iv2))
            work(1+(nj_subset-1)*ni_subset) = data_array(ih1,iv2)
     1          + delx*(data_array(ih1,iv2) - data_array(ih1-1,iv2))
     2          + dely*(data_array(ih1,iv2+1) - data_array(ih1,iv2))
         endif
         if( .NOT. (exact_right .or. exact_top) ) then
            delx = (hi_h - orig_meshh(ih2))/
     1             (orig_meshh(ih2+1) - orig_meshh(ih2))
            dely = (low_v - orig_meshv(iv1-1))/
     1             (orig_meshv(iv1) - orig_meshv(iv1-1))
            work(ni_subset) = data_array(ih2,iv1)
     1          + delx*(data_array(ih2+1,iv1) - data_array(ih2,iv1))
     2          + dely*(data_array(ih2,iv1) - data_array(ih2,iv1-1))
         endif
         if( .NOT. (exact_right .or. exact_bottom) ) then
            delx = (hi_h - orig_meshh(ih2))/
     1            (orig_meshh(ih2+1) - orig_meshh(ih2))
            dely = (hi_v - orig_meshv(iv2))/
     1            (orig_meshv(iv2+1) - orig_meshv(iv2))
            work(ni_subset+(nj_subset-1)*ni_subset) =
     1         data_array(ih2,iv2)
     2         + delx*(data_array(ih2+1,iv2) - data_array(ih2,iv2))
     3         + dely*(data_array(ih2,iv2+1) - data_array(ih2,iv2))
         endif
      endif
c
      return
      end
c
      subroutine abs_and_log_3d(work, ni, nj, do_log, do_abs)
      implicit none
c
      include 'sysdep.cmn'
c
      logical do_log, do_abs
      integer ni, nj
      real work(ni, nj)
c
c     locals
      integer i,j
c
c     Take log/absolute value if necessary...
      do 100 j=1,nj
         do 100  i=1,ni
            if( do_abs ) then
               work(i,j) = ABS(work(i,j))
               if( do_log ) then
c                 Protect against log 0.0
                  work(i,j) = ALOG10(MAX(work(i,j),PM_MIN_FLOAT))
               endif
            endif
100   continue
c
      return
      end
c
      subroutine clip_and_normalize_array_3d(work, ni, nj,
     1           floor, ceil, iexp)
      implicit none
c
      integer ni, nj, iexp
      real work(ni,nj), floor, ceil
c
c     Locals
      integer i,j
      real scale
c
      scale = 10.0**iexp
      do 100 j=1,nj
         do 100  i=1,ni
c           Rudimentary clipping
            if ( work(i,j) .lt. floor ) work(i,j) = floor
            if ( work(i,j) .gt. ceil ) work(i,j) = ceil
            work(i,j) = work(i,j)/scale
100   continue
c
      return
      end
c
      subroutine recompute_minmax_3d(work, ni, nj,
     1           do_log, do_abs, iuser_scale, do_uscale,
     2           dmin, dmax, dmin_scaled, dmax_scaled, majtic, ticfreq,
     3           iaxis_scale)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'units.cmn'
c
      logical do_log, do_abs, do_uscale
      integer ni, nj, iuser_scale, iaxis_scale
      real work(ni,nj)
      real dmin, dmax, dmin_scaled, dmax_scaled, majtic
      integer ticfreq
c
c     locals
      integer i,j
      real user_scale, axis_scale, alog1, h
c
      if( do_uscale ) then
         user_scale = 10.0**iuser_scale
      else
         user_scale = 1.0
      endif
c
c     Get Z axis plot factors
      if( do_abs .and. .NOT. do_log ) then
         dmin = PM_MAX_FLOAT
         dmax = PM_MIN_FLOAT
      else if( do_log ) then
         dmin =  38.0
         dmax = -38.0
      else
         dmin =  PM_MAX_FLOAT
         dmax = -PM_MAX_FLOAT
      endif
      alog1 = ALOG10(PM_MIN_FLOAT)
      do 110 j=1,nj
         do 110 i=1,ni
            if( do_uscale ) then
               if( do_log ) then
                  h = work(i,j) - FLOAT(iuser_scale)
               else
                  h = work(i,j)/user_scale
               endif
            else
               h = work(i,j)
            endif
            if( dmin .gt. h ) then
c              Avoid taking "zero" values as minimum on log scale
               if( .not. do_log .or. h .ne. alog1) then
                  dmin = h
               endif
            endif
            if( dmax .lt. h ) dmax = h
110   continue
c
      call get_default_limits(dmin, dmax,
     1     dmin_scaled, dmax_scaled, majtic,
     2     ticfreq, iaxis_scale, do_log)
      if( do_log ) then
         axis_scale = 1.0
         dmin = 10.0**dmin
         dmax = 10.0**dmax
         dmin_scaled = 10.0**dmin_scaled
         dmax_scaled = 10.0**dmax_scaled
      else
         axis_scale  = 10.0**iaxis_scale
         if( do_uscale ) then
            majtic      = majtic*axis_scale
            dmin_scaled = dmin_scaled*axis_scale
            dmax_scaled = dmax_scaled*axis_scale
         endif
         dmin = dmin*user_scale
         dmax = dmax*user_scale
      endif
c
      return
      end
c
      subroutine conv_pm3d_axes_to_surf(hmin_plot, vmin_plot, zmin_plot,
     1     hmax_plot, vmax_plot, zmax_plot,
     2     h_tic_major, v_tic_major, z_tic_major,
     3     h_minor_tic_freq, v_minor_tic_freq, z_minor_tic_freq,
     5     h_exponent, v_exponent, z_exponent,
     6     horiz_label, horiz_units, vert_label, vert_units,
     7     z_label, z_units,
     8     good_min, good_max, minor_tic, major_tic_freq, exponent,
     9     xyz_labels, xyz_units,
     A     swapxy)
      implicit none
c
      real hmin_plot, vmin_plot, zmin_plot,
     1     hmax_plot, vmax_plot, zmax_plot,
     2     h_tic_major, v_tic_major, z_tic_major
      integer h_minor_tic_freq, v_minor_tic_freq, z_minor_tic_freq
      integer h_exponent, v_exponent, z_exponent
      integer exponent(3)
      real good_min(3), good_max(3), minor_tic(3), major_tic_freq(3)
      character*(*) horiz_label, horiz_units, vert_label, vert_units,
     1              z_label, z_units
      character*(*) xyz_labels(3), xyz_units(3)
      logical swapxy
c
      if( swapxy ) then
         good_min(1) = vmin_plot
         good_min(2) = hmin_plot
         good_max(1) = vmax_plot
         good_max(2) = hmax_plot
         minor_tic(1) = v_tic_major/v_minor_tic_freq
         minor_tic(2) = h_tic_major/h_minor_tic_freq
         major_tic_freq(1) = v_minor_tic_freq
         major_tic_freq(2) = h_minor_tic_freq
         exponent(1) = v_exponent
         exponent(2) = h_exponent
         xyz_labels(1) = vert_label
         xyz_labels(2) = horiz_label
         xyz_units(1) = vert_units
         xyz_units(2) = horiz_units
      else
         good_min(1) = hmin_plot
         good_min(2) = vmin_plot
         good_max(1) = hmax_plot
         good_max(2) = vmax_plot
         minor_tic(1) = h_tic_major/h_minor_tic_freq
         minor_tic(2) = v_tic_major/v_minor_tic_freq
         major_tic_freq(1) = h_minor_tic_freq
         major_tic_freq(2) = v_minor_tic_freq
         exponent(1) = h_exponent
         exponent(2) = v_exponent
         xyz_labels(1) = horiz_label
         xyz_labels(2) = vert_label
         xyz_units(1) = horiz_units
         xyz_units(2) = vert_units
      endif
c
c     Z axis is left alone in any case
      good_min(3) = zmin_plot
      good_max(3) = zmax_plot
      minor_tic(3) = z_tic_major/z_minor_tic_freq
      major_tic_freq(3) = z_minor_tic_freq
      exponent(3) = z_exponent
      xyz_labels(3) = z_label
      xyz_units(3)  = z_units
c
      return
      end
c
      subroutine get_orbit_locator(window_id,
     1               x1, y1, arg, do_orbit)
      implicit none
c
      integer window_id
      real x1, y1, arg
      logical do_orbit
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
c
      integer lookup_ws_id, work_ws_type
      integer work_id,
     1        xform, pet_type, dev_num, data_record(2),
     2        record_buffer_length, input_status, err_stat
      real    echo_area(4)
c
      if( .NOT. gks_is_open ) return
c
      work_id = lookup_ws_id(window_id)
      if( work_id .eq. PLT$EMPTY_WINDOW_SLOT ) then
c        Bad window id
         call plt$error_message('get_orbit_locator','Bad window id')
         do_orbit = .FALSE.
         return
      endif
c
      work_ws_type = ws_struct(work_id)%type
      if( work_ws_type .ne. GKS$K_VSII .and.
     1    work_ws_type .ne. GKS$K_VSII_GPX .and.
     2    work_ws_type .ne. GKS$K_DECWINDOWS .and.
     3    work_ws_type .ne. GKS$K_MOTIF ) then
         call plt$error_message('get_orbit_locator','Not a workstation')
         do_orbit = .FALSE.
         return
      endif
c
      call anno_setup_ws_for_locator(window_id, echo_area, err_stat)
c
c     Put locator at initial point
c     Use echo type 2 to get cross hairs
      dev_num = 1
      xform = 1
      pet_type = 2
c     gks$init_locator modifies record_buffer_length.
c     Reset it to expected value.
      record_buffer_length = 0
      call gks$init_locator(work_id, dev_num, 0.5, 0.5,
     1     xform, pet_type, echo_area,
     2     data_record, record_buffer_length)
      call gks$request_locator(work_id, dev_num, input_status,
     1     xform, x1, y1)
      if( input_status .ne. 1 ) then
c        Input was cancelled
         do_orbit = .FALSE.
      else
c        Decode position of locator into orbit action
         do_orbit = .TRUE.
c        Get angle of locator with respect to mid screen
         arg = ATAN2(y1 - 0.5,x1 - 0.5)
c        Convert to degrees
         arg = arg*180.0/3.14159246
c        Convert to 0 to 360 degrees
         if( arg .lt. 0.0 ) arg = arg + 360.0
      endif
c     Deactivate window
      call deactivate_plot_device
c
      return
      end
c
      subroutine orbit_modify_eye(x1, y1, arg,
     1                            threed_eye_position)
      implicit none
c
      real x1, y1, arg
      real threed_eye_position(3)
c
c     x1, y1  -- position of locator
c     arg     -- angle of vector from mid point of screen to locator
c                in degrees from 0 to 360, counter clockwise from (1,0)
c
      if( arg .le. 45.0 .or. arg .gt. 315.0 ) then
c        Move eye position 10 degrees to right
         threed_eye_position(3) = threed_eye_position(3) + 10
      else if( arg .gt. 45.0 .and. arg .lt. 135.0 ) then
c        Move eye position 10 degrees up
         threed_eye_position(2) = threed_eye_position(2) - 10
      else if( arg .gt. 135.0 .and. arg .lt. 225.0 ) then
c        Move eye position 10 degrees to left
         threed_eye_position(3) = threed_eye_position(3) - 10
      else
c        Move eye position 10 degrees down
         threed_eye_position(2) = threed_eye_position(2) + 10
      endif
c
      return
      end

c
      subroutine compute_clip_of_triarray(data, work,
     1                             cord, nop,
     1                             xvals, yvals, tri_clip,
     1                             min_x, max_x, min_y, max_y,
     1                             tri_npts, tri_ne,
     1                             swapxy, loutside)
      implicit none
c
      integer tri_npts, tri_ne
      integer nop(3,tri_ne)
      real    cord(2,tri_npts), data(tri_npts), work(tri_npts)
      real    xvals(tri_npts), yvals(tri_npts)
      real    min_x, max_x, min_y, max_y
      logical swapxy, loutside, tri_clip(tri_ne)
c
      include 'const.h'
      include 'units.cmn'
c
      real x1,x2,x3, y1,y2,y3, minx, maxx, miny, maxy
      integer i
c
      loutside = .TRUE.
      do i=1,tri_ne
c        First, check to see if triangle is totally outside, or clips
c        bounds.
         x1 = cord(1,nop(1,i))
         y1 = cord(2,nop(1,i))
         x2 = cord(1,nop(2,i))
         y2 = cord(2,nop(2,i))
         x3 = cord(1,nop(3,i))
         y3 = cord(2,nop(3,i))
         minx = MIN(x1,x2,x3)
         maxx = MAX(x1,x2,x3)
         miny = MIN(y1,y2,y3)
         maxy = MAX(y1,y2,y3)
         if( maxx .le. min_x .or. minx .ge. max_x .or.
     1       maxy .le. min_y .or. miny .ge. max_y ) then
c           Triangle totally outside clipping region
            tri_clip(i) = .TRUE.
         else if( maxx .le. max_x .and. minx .ge. min_x .and.
     1            maxy .le. max_y .and. miny .ge. min_y ) then
c           Triangle totally inside clipping region
            tri_clip(i) = .FALSE.
            loutside    = .FALSE.
         else
c           Triangle intersects clipping region
            tri_clip(i) = .TRUE.
            loutside    = .FALSE.
         endif
      end do
c
      if( swapxy ) then
c
c        Swap x and y axes -- many simulators (e.g. MINIMOS) have positive
c        Y axis going down into depth (e.g. left handed coordinate system)
         do i=1,tri_npts
           xvals(i) = pm_scale_cm2plot*cord(2,i)
           yvals(i) = pm_scale_cm2plot*cord(1,i)
           work(i) = data(i)
         end do
c
      else
c
c        Normal X and Y coordinate system
         do i=1,tri_npts
           xvals(i) = pm_scale_cm2plot*cord(1,i)
           yvals(i) = pm_scale_cm2plot*cord(2,i)
           work(i) = data(i)
         end do
c
      endif
c
      return
      end
c
      subroutine abs_and_log_tri3d(work, tri_npts, do_log, do_abs)
      implicit none
c
      include 'sysdep.cmn'
c
      logical do_log, do_abs
      integer tri_npts
      real work(tri_npts)
c
c     locals
      integer i
c
c     Take log/absolute value if necessary...
      do 100 i=1,tri_npts
         if( do_abs ) then
            work(i) = ABS(work(i))
            if( do_log ) then
c              Protect against log 0.0
               work(i) = ALOG10(MAX(work(i),PM_MIN_FLOAT))
            endif
         endif
100   continue
c
      return
      end
c
      subroutine recompute_minmax_tri3d(work, tri_npts,
     1           do_log, do_abs, iuser_scale, do_uscale,
     2           dmin, dmax, dmin_scaled, dmax_scaled, majtic, ticfreq,
     3           iaxis_scale)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'units.cmn'
c
      logical do_log, do_abs, do_uscale
      integer tri_npts, iuser_scale, iaxis_scale
      real work(tri_npts)
      real dmin, dmax, dmin_scaled, dmax_scaled, majtic
      integer ticfreq
c
c     locals
      integer i
      real user_scale, axis_scale, alog1, h
c
      if( do_uscale ) then
         user_scale = 10.0**iuser_scale
      else
         user_scale = 1.0
      endif
c
c     Get Z axis plot factors
      if( do_abs .and. .NOT. do_log ) then
         dmin = PM_MAX_FLOAT
         dmax = PM_MIN_FLOAT
      else if( do_log ) then
         dmin =  38.0
         dmax = -38.0
      else
         dmin =  PM_MAX_FLOAT
         dmax = -PM_MAX_FLOAT
      endif
      alog1 = ALOG10(PM_MIN_FLOAT)
      do 110 i=1,tri_npts
         if( do_uscale ) then
            if( do_log ) then
               h = work(i) - FLOAT(iuser_scale)
            else
               h = work(i)/user_scale
            endif
         else
            h = work(i)
         endif
         if( dmin .gt. h ) then
c           Avoid taking "zero" values as minimum on log scale
            if( .not. do_log .or. h .ne. alog1) then
               dmin = h
            endif
         endif
         if( dmax .lt. h ) dmax = h
110   continue
c
      call get_default_limits(dmin, dmax,
     1     dmin_scaled, dmax_scaled, majtic,
     2     ticfreq, iaxis_scale, do_log)
      if( do_log ) then
         axis_scale = 1.0
         dmin = 10.0**dmin
         dmax = 10.0**dmax
         dmin_scaled = 10.0**dmin_scaled
         dmax_scaled = 10.0**dmax_scaled
      else
         axis_scale  = 10.0**iaxis_scale
         if( do_uscale ) then
            majtic      = majtic*axis_scale
            dmin_scaled = dmin_scaled*axis_scale
            dmax_scaled = dmax_scaled*axis_scale
         endif
         dmin = dmin*user_scale
         dmax = dmax*user_scale
      endif
c
      return
      end
c
      subroutine clip_and_normalize_array_tri3d(work, tri_npts,
     1           floor, ceil, iexp)
      implicit none
c
      integer tri_npts, iexp
      real work(tri_npts), floor, ceil
c
c     Locals
      integer i
      real scale
c
      scale = 10.0**iexp
      do 100  i=1,tri_npts
c        Rudimentary clipping
         if ( work(i) .lt. floor ) work(i) = floor
         if ( work(i) .gt. ceil ) work(i) = ceil
         work(i) = work(i)/scale
100   continue
c
      return
      end
