!
      module input
!
! Read and validate user input parameters.
! (also define several variables in proc)
!
      use proc
      use flist,only: deredund
      use mk_v5d
      use mk_hvols
      use fields,only: ne5577,e5577lab
      implicit none
!
! User namelist inputs:
!
! histvols(mxhvols) = history volume mss paths requested by user
! mtimes(3,mxtms)   = model times requested by user (day,hr,min)
! cfields(mxfproc)  = char*8 names of requested fields
! modelhts          = if==1, get heights from history, otherwise calculate
! iden		    = density conversion flag:
!   iden = 0 -> leave species as on history (most are mass mixing ratio)
!   iden = 1 -> convert species to number densities (cm3)
!   iden = 2 -> convert species to number density mixing ratios
!   iden = 3 -> convert species to mass density (gm/cm3)
!   iden = 4 -> convert species to parts per billion by volume (PPBV)
! For special Bougher format data:
!   aux_fields(mxfproc) = char*8 names of requested auxillary fields
!   bf_format           = format statement describing Bougher output format
!
      character(len=lenhvol) :: histvols(mxhvols) = " "	! history volumes
      character(len=lenhvol) :: histvols_cntr(mxhvols) = " "
      character(len=8)       :: cfields(mxfproc) = " "	! requested fields
      character(len=16)      :: cunits(2,mxfproc) = " "	! optional units
      character(len=80)      :: tmpdir = ' '		! temp directory
      character(len=120)     :: rec80			! temp char var
      character(len=8)       :: aux_fields(mxfproc) = " " ! auxillary fields
      character(len=132)     :: bf_format=" "             ! Bougher format statement
      character(len=80)      :: grid_levels(2,mxfproc)=" "  ! midpoints or interfaces
      integer,private :: i
      integer :: 
     +  mtimes(3,mxtms)=ispval,	! requested model times
     +  mtimes_cntr(3,mxtms)=ispval,	! requested model times for control
     +  idifpercent=ispval,	! difference fields flag
     +  isearch=0,		! history search flag
! 12/8/06 btf: upped default maxdiskvols from 50 to 500
     +  maxdiskvols=500,	! default max number of hist vols allowed on disk
     +  modelhts = 1,		! heights flag
     +  iden = 1,		! density conversion flag
     +  ionvel = 1,		! ion velocity flag
     +  ie6300 = 0,		! E6300 emission flag
     +  ie5577(ne5577)=(/1,0,0,0,0/),! E5577 emission flags
				     ! (see also e5577lab in fields.f)
     +  ieo200(3) = (/1,0,0/),	! flags for O2 00 calculation
     +  ibohv_watts = 0,	! units flag for ht-integ OH
     +  integht = 1,		! integrate certain fields in height
     +  iprint_fknown=0,	! print list of all valid field names
     |  ilab_hq=0,              ! see label_quality
     |  iplot = 1,              ! whether or not to make plots
     |  ismooth = 1             ! smooth height-interpolated diffs (mklons,mkxyloc)
!
! lat slices:
      integer :: 
     +  ipltlat=0,		! plot flag for latitude slices
     +  ilat_log10=1,		! plot log10 sp if ilat_log10 > 0
     +  ilat_yaxright=2		! extra right-hand y-axes
      character(len=8) ::
     +  clat_xaxes(2)=' '	! desired x axes ('LON' and/or 'LT')
      real(kind=8) :: 
     +  flats(mxslice)=spval 	! selected latitudes
      real :: 
     +  flat_zprange(2)=spval,	! log pressure bottom,top (y-axis)
     +  flat_htscale(3)=spval	! height bottom,top,delta (y-axis)
      real ::
     |  flat_censlt=spval       ! center x-axis of lat slices at desired local time
!
! lon slices:
      integer ::
     +  ipltlon=0,		! plot flag for longitude slices
     +  ilon_log10=1,		! plot log10 sp if ilon_log10 > 0
     +  ilon_yaxright=2,	! extra right-hand y-axes
     +  amphase=0,		! Wave number for amplitude/phase
     +  istream=0,		! Stream functions for V and/or W
     +  ilon_epvdiv_yz=0,	! if > 0, add epvy+epvz vectors to 
                                !   epvdiv contours
     +  ilon_epvyzmag_yz=0,	! if > 0, add epvy+epvz vectors to 
                                !   epvyzmag contours
     |  ilon_timeave=0,         ! if > 0, save time-averaged data
     |  ilon_timeave_only=0     ! if > 0, plot time-averaged data only
      real(kind=8) ::
     +  flons(mxslice)=spval	! selected longitudes for lon slices
      real ::
     +  fslts(mxslice)=spval,	! selected local times for lon slices
     +  flon_zprange(2)=spval,	! log pressure bottom,top (y-axis)
     +  flon_htscale(3)=spval,	! height bottom,top,delta (y-axis)
     +  flon_xlatrange(2)=spval	! latitude range for x-axis of lon slices
!
! ipltzmslt:
! lon slices of zm at selected slt: must be a time series. Calculate the
! the average across all ut's of longitudes corresponding to the selected slt,
! and plot lat vs zp/ht. See correspondence w/ Steve Bougher in January, 2011:
!
! 2-D  (latitude vs ht/zp)
! Field at a specific LT that has been averaged over all longitudes:
! An entire diurnal cycle of histories (say at 1-hour intervals) would
! be needed in order to capture all longitudes (and UTs) corresponding
! to a specific local time.  The averaging of this field over all these
! longitudes would yield spatial averages at a given local time.
! The 2-D plot most useful would indeed be latitude vs ht/zp. 
!
      integer ::
     |  ipltzmslt=0, 
     |  izmslt_yaxright=2,
     |  izmslt_log10=1
      real ::
     |  zmslt(mxslice) = spval             ! requested local time at which to plot zonal means
      real ::
     +  zmslt_zprange(2)=spval,	! range in zp (bot,top) for y-axis
     +  zmslt_htscale(3)=spval	! height scale (bot,top,delta)
!
! maps:
      integer ::
     +  ipltmaps=0,		! plot flag for maps
     +  map_global=0,		! CE projections (cyl equidistant)
     +    map_global_cenlon=ispval,	!   center x-axis longitude
     +    map_global_censlt=ispval,	!   center x-axis local time
     +  map_polar=0,		! ST projections (polar stereographic)
     +  map_satview=0,		! SV projections (sat view)
     +  map_mollweide=0,	! MO projections (mollweide)
     +  map_continents=0,	! continental outlines flag
     +  map_tn_unvn=0,		! if > 0, add un+vn vectors to tn contours
     +  map_ht_unvn=0,		! if > 0, add un+vn vectors to ht contours
     +  map_ep_uivi=0,		! if > 0, add ui+vi vectors to ep contours
     +  map_hmf2=0,		! if > 0, add hmf2 heights to utvert plots
!
! If map_log10==1, take log10 of density fields. 
! If map_log10==2, force log10 of all fields processed. 
! (default is map_log10=0 -> do not take log10)
     +  map_log10=0
!
! These (and others) are typed 8-bytes so they can pass through fixnamelist
      real(kind=8) ::
     |  fmap_polar_perimlat(mxperimlat)=spval
      real(kind=8) ::
     |  fmap_zpht(mxzpht)=spval             ! selected zp/ht
!
      real ::
     +  fmap_satview_eradii=6.631,		  ! earth radii from satv 
     +  fmap_satview_latlon(2)=(/spval,spval/),	  ! lat-lon center of satv
     +  fmap_satview_latslt(2)=(/spval,spval/),	  ! lat-slt center of satv
     +  fmap_mollweide_latlon(2)=(/spval,spval/), ! lat-lon center of moll
     +  fmap_mollweide_latslt(2)=(/spval,spval/), ! lat-slt center of moll
     |  fmap_latminmax(2)=(/spval,spval/)         ! min,max latitude to plot
!
! xy vertical profiles at selected locations:
      integer ::
     +  ipltxyloc=0,		! flag for vertical profile xy plots
     +  ixyloc_log10=1,		! flag for log10 of vertical profiles
     +  ixyloc_yaxright=1	! flag for extra y-axes of ipltxyloc
      real(kind=8) ::
     +  xylocs(2,mxloc)=spval	! lat,lon locations for ipltxyloc
      real ::
     +  xyloc_zprange(2)=spval,	! zp range bottom,top for ipltxyloc
     +  xyloc_htscale(3)=spval	! ht scale for ipltxyloc
      character(len=56) ::
     +  xyloc_locname(mxloc)=' '	! optional names of locations
!
! ut vs zp/ht contours:
      integer ::
     +  ipltutvert=0,		! flag for ut vs vertical contours
     +  iutvert_log10=1,	! log10 flag
     +  iutvert_yaxright=2	! extra right-hand y-axes
      real(kind=8) ::
     +  utvert_locs(2,mxloc)=spval	! lat,lon pairs of locations
      real ::
     +  utvert_zprange(2)=spval,	! range in zp (bot,top) for y-axis
     +  utvert_htscale(3)=spval		! height scale (bot,top,delta)
      character(len=56) ::
     +  utvert_locname(mxloc)=' '	! optional names of locations
!
! ut vs zp/ht (satellite track) contours:
      integer ::
     +  ipltsatut=0,                    ! flag for ut vs vertical (satellite) contours
     +  isatut_log10=1, ! log10 flag
     +  isatut_yaxright=2               ! extra right-hand y-axes
      real(kind=8) ::
     +  satut_locs(2,mxloc)=spval       ! lat,lon pairs of locations
      real ::
     +  satut_zprange(2)=spval,         ! range in zp (bot,top) for y-axis
     +  satut_htscale(3)=spval          ! height scale (bot,top,delta)
      character(len=56) ::
     +  satut_locname(mxloc)=' '        ! optional names of locations
      character(len=120) :: satfile	! file with sat orbit points
!
! ut vs lat contours:
      integer ::
     +  ipltutlat=0,		! flag for ut vs latitude contours
     +  iutlat_log10=1,		! log10 flag
     +  iutlat_sltmaps=0	! maps at constant slt's
      real(kind=8) ::
     +  utlat_zphtlon(2,mxslice)=spval	! zp/ht, lon/slt/zm pairs
      real :: utlat_latminmax(2)=spval  ! min,max latitude to plot
!
! ut vs lon contours:
      integer ::
     +  ipltutlon=0,            ! flag for ut vs longitude contours
     +  iutlon_yslt,
     +  iutlon_log10=1
      real(kind=8) ::
     +  utlon_zphtlat(2,mxslice)=spval  ! zp/ht, lat pairs
!
! ut vs field xy line plots:
      integer ::
     +  ipltxyut=0,		! flag for ut vs field line plots
     +  ixyut_log10=1,		! log10 flag
     +  ixyut_doppler=0		! if > 0, make doppler t,u,v
      real(kind=8) ::
     +  xyut_locs(2,mxloc)=spval	! lat,lon pairs of locations
      real ::
     +  xyut_zpht(mxzpht)=spval		! selected zp and/or hts
      character(len=56) ::
     +  xyut_locname(mxloc)=' '	! optional names of locations
!
! For vis5d files:
      real :: v5d_zprange(2)=(/spval,spval/),	! zp range (bot,top)
     +  v5d_slt=spval			! fixed local time
      character(len=8) :: 
     +  v5d_vtype = ' '			! 'reg_zp' or 'ireg_mb'
!
! msc plotting:
      character(len=8) :: outplt(3) = 
     +  (/'cgm     ','        ','        '/)
      character(len=8) :: 
     |  psmode = 'port    ',
     |  label_quality = 'medium  '! 'low','medium','high' for plchlq, mq, or hq
      integer :: 
     +  icolor=0,		! color fill contour flag
     +  igreyscale=0,		! color fill using greyscale contour flag
     +  ibox_clabs=0,		! box contour line labels
     +  iboxplt=0,		! put box around plot if iboxplt > 0
     +  ishadeneg=0, 		! shade neg contour areas (monochrome)
     +  multiplt=0, 		! flag for multiple plots per frame
     +  multiadvfr=0, 		! frame advance flag for multiplt
     +  iprint_fldminmax=0,	! print 3d min,max of fields every hist
     +  ipltrowcol(2)=(/2,2/),	! number of rows,cols of plots per frame
     +  idoppler=0,		! sets ixyut_doppler (compat w/ old proc)
     +  ivec_label,		! for compatability w/ old proc inputs
     +  iwr_label(6)=1		! array to set which labels are written
!
! 1st dim of fmnmxint and cmnmxint is 'fldname',cmin,cmax,cint)
      real(kind=8) ::
     +  fmnmxint(4,mxfproc)=spval ! optional min,max,cint control
      character(len=16) :: cmnmxint(4,mxfproc)=' ' ! min,max,cint for AIX
      character(len=16) :: cscale(2,mxfproc)=' '   ! scale factor for AIX
      real ::
     +  fscale(2,mxfproc)=spval,  ! optional scale factor for each field
     +  vn_scale(3)=(/0.,0.,0./), ! low,high,scale for neutral vectors
     +  vi_scale(3)=(/0.,0.,0./), ! low,high,scale for ion drifts
     +  vep_scale(3)=(/0.,0.,0./),! low,high,scale for epflux vectors
     +  vepz_viewfac=1.,	  ! scale factor for epz vector
     +  vmag_len=0.		  ! vector arrow length	 	
!
! cdf_htscale(3) for height interpolation of fields saved to netcdf
! output file (sendcdf) (see mkcdf.f)
      real :: cdf_htscale(3)=(/spval,spval,spval/) ! optional height scale
      real :: cdf_zprange(2)=(/spval,spval/)       ! optional pressure scale
      real :: cdf_zonalmean=spval                  ! optional zonal means
      integer :: ncid_out=0       ! file unit of netcdf output file (see sub mkcdf)
      integer :: mtime_del_mins   ! delta time between histories (mins)
!
! Machine:path directories or files to which output files are to be rcp'd:
      character(len=120) :: 
     +  sendcgm = ' ',		! remote path for cgm
     +  sendps  = ' ',		! remote path for ps
     +  sendcdf = ' ',		! remote path for netcdf
     +  senddat = ' ',		! remote path for ascii data file
     +  senddat_bf = ' ',       ! remote path for ascii data file in Bougher format
     +  sendxdr = ' ',		! remote path for xdr data file
     +  sendv5d = ' ',		! remote path for vis5d data file
     +  sendsat = ' ',		! remote path for satellite netcdf
     +  sendcdf_lons = ' ',	! remote path for lon slice netcdf file 
     +  sendcdf_xyut = ' ',	! remote path for ut vs field netcdf file 
     +  sendcdf_utvert = ' ',	! remote path for ut vs vertical netcdf file 
!
! Mss paths to dispose output files:
     +  sendms_v5d = ' ',	! remote path for vis5d data file
     +  sendms_xdr = ' '	! mss remote path for xdr data file
!
! Local disk file names for output files:
! These files are built on the disk by the processor, and scp'd at the
!   end of the run (by scpfile) to the users requested machine:path
!   (namelist read parameters sendcgm, sendps, senddat, etc).
! (flnm_v5d (vis5d output file name) is declared in main tgcmproc, 
!  and defined by mkv5d))
!
      character(len=120) ::
     |  flnm_cgm = "tgcmproc.cgm",	 ! cgm metacode
     |  flnm_ps  = "tgcmproc.ps",	 ! postscript
     |  flnm_dat = " ",			 ! ascii data
     |  flnm_xdr = "tgcmproc.xdr",	 ! xdr data
     |  flnm_bf_dat = "tgcmproc_bf.dat", ! ascii data (Bougher)
     |  flnm_cdf = " ",                  ! netcdf output file name
     |  flnm_cdf_xyut = " "              ! netcdf xyut output file name
!
! Path to NASA AMES GCM lower boundaries file for mtgcm (mars tgcm):
! (mgcm_fspath is path on file server, e.g., /fs/othrorgs/home0/foster/...)
!
      character(len=80) :: mgcm_fspath = ' '
!
! NASA AMES MGCM lower boundaries for mars tgcm (read from mgcm_fspath):
      real :: zkmbot_mgcm(nlat), tbot_mgcm(nlat)
!
! Path to N2 file for vtgcm (venus tgcm):
! (vtgcm_fspath is path on file server like mgcm_fspath)
!
      character(len=80) :: vtgcm_fspath = ' '
!
! Path to magfield data and dip/declination data files (see mkdrifts.F):
      character(len=1024) :: ncfile_dipdec=' ', ncfile_magfield=' '
!
! Sattrack settings
      integer :: isattrack = 0,   ! set to 1 for sattrack runs
     +           idetrend = 0,    ! set to 1 to detrend up/down legs
     +           isatplt=1,       ! 1=up, 2=down, 3=average
     +           iwrsat=1,        ! 1=write sattrack only, 2=write all fields
     +           ibaseday=0       ! baseline day for precession calculation
      real ::    prec_rate=0.,    ! precession rate in mins/day, E+, W-
     +           latmax=90.,      ! maximum latitude of satellite path
     +           latmin=-90.,     ! minimum latitude of satellite path
     +           loctime0=0.      ! local time of latmax
!
! For compatability w/ old processor:
      integer ::
     +  iproc_sut,iproc_mut,	! for compatability w/ old proc inputs
     +  iemis_integ(6),		! for compatability w/ old proc inputs
     +  ieohv(10),		! for compatability w/ old proc inputs
     +  ibohv(60)		! for compatability w/ old proc inputs
      contains
!===================================================================
      subroutine getinp
!
! Read and verify namelist inputs.
!
! Namelist inputs (from above):
      namelist/proc_input/ histvols,histvols_cntr,mtimes,mtimes_cntr,
     +  isearch,maxdiskvols,cfields,cunits,modelhts,iden,ionvel,tmpdir,
     +  ie6300,ie5577,ieo200,ibohv_watts,integht,idifpercent,
     +  iprint_fknown,iplot,ismooth,
! Bougher format:
     |  aux_fields,bf_format,
! midpoints vs interfaces:
     |  grid_levels,
! lat slices:
     +  ipltlat,flats,flat_zprange,flat_htscale,ilat_log10,
     +  ilat_yaxright,clat_xaxes,flat_censlt,
! lon slices:
     +  ipltlon,flons,fslts,flon_zprange,flon_htscale,ilon_log10,
     +  ilon_yaxright,flon_xlatrange,amphase,istream,ilon_epvdiv_yz,
     +  ilon_epvyzmag_yz,ilon_timeave,ilon_timeave_only,sendcdf_lons,
! zmslt lon slices:
     |  ipltzmslt,zmslt,zmslt_zprange,zmslt_htscale,izmslt_yaxright,
     |  izmslt_log10,
! xy vertical profiles at selected locations:
     +  ipltxyloc,ixyloc_log10,ixyloc_yaxright,xylocs,xyloc_zprange,
     +  xyloc_htscale,xyloc_locname,
! maps:
     +  ipltmaps,map_global,map_global_cenlon,map_global_censlt,
     +  map_polar,map_satview,map_mollweide,fmap_zpht,map_continents,
     +  fmap_polar_perimlat,fmap_satview_eradii,fmap_satview_latlon,
     +  fmap_satview_latslt,fmap_mollweide_latlon,fmap_mollweide_latslt,
     +  map_tn_unvn,map_ht_unvn,map_ep_uivi,map_hmf2,map_log10,
     |  fmap_latminmax,
! utvert:
     +  ipltutvert,iutvert_log10,utvert_locs,utvert_zprange,
     +  utvert_htscale,utvert_locname,iutvert_yaxright,sendcdf_utvert,
! satut:
     +  ipltsatut,isatut_log10,satut_locs,satut_zprange,
     +  satut_htscale,satut_locname,satfile,isatut_yaxright,
! utlat:
     +  ipltutlat,iutlat_log10,utlat_zphtlon,iutlat_sltmaps,
     |  utlat_latminmax,
! utlon:
     +  ipltutlon,iutlon_yslt,iutlon_log10,utlon_zphtlat,
! xyut (ut on x-axis, field on y-axis at selected locations and zpht):
     +  ipltxyut,ixyut_log10,xyut_locs,xyut_zpht,xyut_locname,
     +  ixyut_doppler,sendcdf_xyut,
! msc plotting/processing:
     +  icolor,igreyscale,ibox_clabs,iboxplt,ishadeneg,outplt,psmode,
     +  multiplt,
     +  multiadvfr,ipltrowcol,fmnmxint,cmnmxint,fscale,iprint_fldminmax,
     +  idoppler,ivec_label,vn_scale,vi_scale,vep_scale,vepz_viewfac,
     +  vmag_len,label_quality,iwr_label,cscale,
! zp range bottom,top for vis5d files:
     +  v5d_zprange,v5d_slt,v5d_vtype,
! height scale for netcdf output file:
     |  cdf_htscale,
! range of zp for netcdf output file:
     |  cdf_zprange,
! zonal mean switch for netcdf output file (0 or 1):
     |  cdf_zonalmean,
! output files:
     +  sendcgm,sendps,sendcdf,senddat,sendxdr,sendv5d,sendms_v5d,
     +  sendms_xdr,sendsat,senddat_bf,
! output disk files
     |  flnm_cgm,flnm_ps,flnm_dat,flnm_xdr,flnm_bf_dat,flnm_cdf,
     |  flnm_cdf_xyut,
! input file for mtgcm:
     +  mgcm_fspath,
! input file for vtgcm:
     +  vtgcm_fspath,
! sattrack settings
     +  isattrack,idetrend,isatplt,iwrsat,
     +  ibaseday,latmax,latmin,prec_rate,loctime0,
! compatability w/ old processor inputs (not used by this proc):
     +  iproc_sut,iproc_mut,iemis_integ,ieohv,ibohv,
! data files used in UI,VI calculations:
     |  ncfile_dipdec, ncfile_magfield
! End namelist inputs.
!
! Locals:
      integer ios,i,ii,n,nfinp,lenf,luin,nzpht,ntimes_cntr,nzprange,
     +  nhtscale,nlons,nlats,nslts,idt,nfinp_aux
      character(len=lenhvol) :: hvols(mxhvols) = " "
      real :: slt,fmin,fmax,fcint
!
! Externals:
      integer,external :: nextlu,iunlink,fseries,mtime_to_mins,
     |  fixnamelist
      character(len=16),external :: float_to_str,float8_to_str
      logical,external :: isslt
!
! Remove comments from input file:
      write(6,"(/20('-'),' Input file: ',20('-'))")
      luin = rmcomments(5,';',1)
      write(6,"(18('-'),' End input file, begin input verification ',
     +  18('-'))")
!
! IBM namelist does not allow character strings values to be input
! to floats, e.g., xlocs = 60.,'slt12.', 'zm','zm' (This is allowed
! under SGI, and is detected by subs chlons, chlats with use of
! function float_is_string (util.F). So if we are under AIX, call
! sub fixnamelist. This sub checks for, e.g., 'zm', and if found
! replaces it with zmflag for the namelist read.
!
#if defined(AIX) || defined(SUN) || defined(LINUX)
      write(6,"('before fixnamelist: luin=',i3)") luin
      luin = fixnamelist(luin)
      write(6,"('after fixnamelist: luin=',i3)") luin
#endif
!
! Namelist read:
      read(luin, nml=proc_input,iostat=ios,err=900)
      if (ios.ne.0) then
        write(6,"('>>> getinp: error from namelist read: ios=',i5)")
     +    ios
        stop 'INPUT'
      endif

!     write(6,"('input: xylocs=')")
!     do i=1,mxloc
!       write(6,"('  i=',i3,' xylocs(:,i)=',2f10.5)") 
!    +    i,xylocs(:,i)
!       if (trim(float_to_str8(xylocs(1,i))) == 'zm'.or.
!    +      trim(float_to_str8(xylocs(1,i))) == 'ZM')
!    +    write(6,"('(zm lat)')")
!       if (trim(float_to_str8(xylocs(2,i))) == 'zm'.or.
!    +      trim(float_to_str8(xylocs(2,i))) == 'ZM')
!    +    write(6,"('(zm lon)')")
!     enddo

      close(luin)
      write(6,"(' ')")
!
! tmpdir:
! 6/22/04: If tmpdir is set to cwd, ignore it (set it to blanks).
!          This is to avoid linking problems between cwd and tmpdir
!          when they are the same.
!
      if (len_trim(tmpdir) > 0) then
        if (trim(tmpdir)=='.'.or.trim(tmpdir)==trim(cwd)) then
          write(6,"(/,'NOTE: Temporary scratch directory (',a,')')")
     |      trim(cwd)
          write(6,"('      is same as execution directory -- will',
     |      ' use only cwd.')")
          tmpdir = ' '
        else
          write(6,"('Temporary scratch directory (tmpdir) = ',a)")
     +      trim(tmpdir)
        endif
      else
        write(6,"(/,'NOTE: Temporary scratch directory (tmpdir)',
     +    ' was NOT specified.')")
        write(6,"(6x,'Processor will search only in the cwd')")
      endif
!
! Expand any env vars in history file paths:
      do i=1,mxhvols
        if (len_trim(histvols(i))>0) then 
          call expand_path(histvols(i))
        endif
      enddo
!
! History volumes:
! (mkhvols either echoes histvols to hvols, or if histvols(2)=='to',
!  then it expands histvols from 'volfirst','to','vollast','by','n'
!  to hvols)
!
      nvols = mkhvols(histvols,hvols,mxhvols)
      if (nvols.eq.0) then
        write(6,"(/'>>> Need history volume paths (histvols)')")
        stop 'histvols'
      endif
      histvols = hvols
      write(6,"('There are ',i3,' history volumes (histvols):')") nvols
      do i=1,mxhvols
        if (len_trim(histvols(i))>0) write(6,"(a)") trim(histvols(i))
      enddo
!
! Model times:
      ntimes = chmtimes(mtimes,mxtms)
      if (ntimes.le.0) then
        write(6,"('>>> Need model times (mtimes)')")
        stop 'mtimes'
      endif
      write(6,"('Model times (mtimes): '/
     +  7(i3,':',i2.2,':',i2.2,' '))") 
     +  ((mtimes(ii,i),ii=1,3),i=1,ntimes)
!
! Must have ntimes > 1 to plot ut on x-axis:
      if (ntimes==1) then
        if (ipltutvert > 0) then
          write(6,"('>>> WARNING: ntimes=',i4,' but ipltutvert=',
     |      i3,': need ntimes > 1 to put ut on x-axis',/,
     |      13x,'Turning ipltutvert off.')") ntimes,ipltutvert
          ipltutvert = 0
        endif
        if (ipltutlat > 0) then
          write(6,"('>>> WARNING: ntimes=',i4,' but ipltutlat=',
     |      i3,': need ntimes > 1 to put ut on x-axis',/,
     |      13x,'Turning ipltutlat off.')") ntimes,ipltutlat
          ipltutlat = 0
        endif
        if (ipltutlon > 0) then
          write(6,"('>>> WARNING: ntimes=',i4,' but ipltutlon=',
     |      i3,': need ntimes > 1 to put ut on x-axis',/,
     |      13x,'Turning ipltutlon off.')") ntimes,ipltutlon
          ipltutlon = 0
        endif
        if (ipltxyut > 0) then
          write(6,"('>>> WARNING: ntimes=',i4,' but ipltxyut=',
     |      i3,': need ntimes > 1 to put ut on x-axis',/,
     |      13x,'Turning ipltxyut off.')") ntimes,ipltxyut
          ipltxyut = 0
        endif
      endif ! ntimes==1
!
! Get delta time between data points (mins):
      if (ntimes > 1.and.
     |  (ipltutvert > 0.or.ipltutlat>0.or.ipltutlon>0.or.ipltxyut>0.or.
     |   ipltsatut > 0)) then
        i=0
        do n=2,ntimes
          mtime_del_mins = mtime_to_mins(mtimes(:,n))-
     |                     mtime_to_mins(mtimes(:,n-1))
          if (n > 2 .and. i /= mtime_del_mins) write(6,"('WARNING ',
     |      'input: delta time is changing: n=',i4,' mtime_del_mins=',
     |      i4,' previous del=',i4)") n,mtime_del_mins,i
          i = mtime_del_mins
        enddo ! n=1,ntimes
        write(6,"('input: mtime_del_mins=',i4)") mtime_del_mins
      endif
!
! histvols_cntr and mtimes_cntr for diff fields:
      diffs = .false.
      nvols_cntr = mkhvols(histvols_cntr,hvols,mxhvols)
      histvols_cntr = hvols
      ntimes_cntr = chmtimes(mtimes_cntr,mxtms)
      if (ntimes_cntr>0.and.nvols_cntr==0) then
        write(6,"('>>> input: need histvols_cntr (got mtimes_cntr)')")
        stop 'nvols_cntr'
      elseif (ntimes_cntr==0.and.nvols_cntr>0) then
        write(6,"('>>> input: need mtimes_cntr (got histvols_cntr)')")
        stop 'mtimes_cntr'
      elseif (ntimes_cntr>0.and.nvols_cntr>0) then
        diffs = .true.
        write(6,"(/'This will be a DIFFERENCE FIELDS run:')")
        write(6,"('  Control history volumes (histvols_cntr):')")
        do i=1,mxhvols
          if (len_trim(histvols_cntr(i)) > 0) then
            write(6,"('  ',a,'  ')",advance="NO") trim(histvols_cntr(i))
            if (mod(i,3).eq.0) write(6,"(' ')")
          endif
        enddo
        write(6,"(' ')")
        write(6,"('  Control model times (mtimes_cntr): '/,'  ',
     +    7(i3,':',i2.2,':',i2.2,' '))") 
     +    ((mtimes_cntr(ii,i),ii=1,3),i=1,ntimes_cntr)
        if (ntimes/=ntimes_cntr) then
          write(6,"('>>> input: for diffs, must have same number',
     +      ' of mtimes and mtimes_cntr: ntimes=',i4,' ntimes_cntr=',
     +      i3)") ntimes,ntimes_cntr
          stop 'mtimes_cntr'
        endif
        do i=1,ntimes
          if (mtimes(2,i)/=mtimes_cntr(2,i).or.
     +        mtimes(3,i)/=mtimes_cntr(3,i)) then
            write(6,"('>>> input: ut times must be same for ',
     +        'mtimes and mtimes_cntr')")
            write(6,"('    (model days may be different)')")
            write(6,"('    i=',i3,' mtimes(:,i)=',3i4,
     +        ' mtimes_cntr(:,i)=',3i4)") i,mtimes(:,i),mtimes_cntr(:,i)
            stop 'mtimes_cntr'
          endif
        enddo
!
! idifpercent==ispval (user did not provide override) -> 
!   take percent diffs of species, raw diffs of all other fields.
! idifpercent==0 (user overrode w/ idifpercent<=0) ->
!   take raw diffs of all fields.
! idifpercent==1 (user overrode w/ idifpercent>0) ->
!   take percent diffs of all fields.
!
        if (idifpercent==ispval) then
          write(6,"('  Diffs will be percent of densities, and raw ',
     +      ' diffs of all other fields.')")
        elseif (idifpercent<=0) then
          idifpercent = 0
          write(6,"('  idifpercent=',i2,' Diffs will be raw diffs ',
     +      'of all fields.')") idifpercent
        else
          idifpercent = 1
          write(6,"('  idifpercent=',i2,' Diffs will be percent diffs ',
     +      'of all fields.')") idifpercent
        endif
      endif ! diffs
!
! isearch=0 -> stop if a history not found
! isearch=1 -> continue search for next history if a history is not found
      if (isearch == 0) then
        write(6,"('isearch = ',i2,' (stop if a history is not found)')")
     +    isearch
      else
        write(6,"('isearch = ',i2,' (continue to next history if ',
     +    'current history is not found)')") isearch
      endif
!
! maxdiskvols (mxdiskvols in proc.f) = maximum number of history volumes
! to remain on the disk during execution.
      if (maxdiskvols <= 0) maxdiskvols = 1
      mxdiskvols = maxdiskvols
      write(6,"('maxdiskvols = ',i3,' (max number of history volume',
     +  ' files allowed on the disk)')") maxdiskvols
!
! Requested fields:
      nfinp = deredund(cfields)	! remove redundancy in field names
      if (nfinp.eq.0) then
        write(6,"('>>> Need field names (cfields)')")
        stop 'cfields'
      else
        write(6,"('Requested fields (cfields):',/(8a9))")
     +    cfields(1:nfinp)
      endif
!
! Auxillary fields:
      nfinp_aux = deredund(aux_fields)  ! remove redundancy in auxillary names
      if (nfinp_aux/=0) then
        write(6,"('Auxillary fields (aux_fields):',/(8a9))")
     +    aux_fields(1:nfinp_aux)
      endif
!
! Optional specification of units:
      nfinp = 0
      do i=1,mxfproc
        if (len_trim(cunits(1,i)) > 0) then
          if (len_trim(cunits(2,i)) > 0) then
            nfinp = nfinp+1 
          else
            write(6,"('>>> INPUT: cunits must be specified in character'
     +        ,' pairs: ''FIELD1'',''UNITS1, ''FIELD2'',''UNITS2''',
     +        ', etc')") 
            stop 'cunits'
          endif
        endif
      enddo
      if (nfinp > 0) then
        write(6,"('User specified field units:')")
        do i=1,mxfproc
          if (len_trim(cunits(1,i)) > 0) 
     +      write(6,"('  Field: ',a,' Units: ',a)")
     +        cunits(1,i),cunits(2,i)
        enddo
      endif
!
! Grid levels (interfaces or midpoints in vertical dimension):
!   grid_levels = 'INTERFACES': process all fields on interfaces
!   grid_levels = 'MIDPOINTS' : process all fields on midpoints
!   grid_levels = 'HISTORY'   : fields as on the history
!   grid_levels = 'HIST+TUV'  : fields as on the history, except T,U,V on interfaces
!   Default: T,U,V are interpolated to interfaces, with TLBC,ULBC,VLBC
!     as the bottom interface boundary. All other fields processed as
!     on the history (lev or ilev).
!   Default can be overriden on a per field basis, e.g.:
!   grid_levels = 'TN','MID', 'O2','INT'...
!
      call verify_grid_levels(grid_levels) 
!
! Emissions fields options:
      if (any(cfields=='E6300   ')) then
        if (ie6300 > 0) then
          write(6,"('ie6300=',i2,': will add SR63 to E6300 ',
     +      'calculation.')") ie6300
        else
          write(6,"('ie6300=',i2,': will NOT add SR63 to E6300 ',
     +      'calculation.')") ie6300
        endif 
      endif
      if (any(cfields=='E5577   ')) then
        write(6,"('Components for E5577 emission: ie5577(',i2,')=',
     +    5i3)") ne5577,ie5577
        do i=1,ne5577
          if (ie5577(i) > 0) write(6,"(a,',')",advance='no') 
     +      trim(e5577lab(i))
        enddo 
        write(6,"(' ')")
      endif
!
! Heights flag:
      if (modelhts > 0) then
        write(6,"('modelhts=',i2,': Heights will be obtained from ',
     +    'the model histories.')") modelhts
      else
        write(6,"('modelhts=',i2,': Heights will be calculated from',
     +    ' tn and mean mass.')") modelhts
      endif
!
! Density conversion flag:
      select case (iden)
      case (0)
        write(6,"('iden = 0: leave species as on history ',
     +    '(most are mass mixing ratio)')")
      case (1)
        write(6,"('iden = 1: convert species to number densities',
     +    ' (cm3)')")
      case (2)
        write(6,"('iden = 2: convert species to number density ',
     +    'mixing ratios')")
      case (3)
        write(6,"('iden = 3: convert species to mass density ',
     +    '(gm/cm3)')")
      case (4)
        write(6,"('iden = 4: convert species to parts per billion',
     |    ' by volume (ppbv)')")
      case default
        write(6,"(/'>>> Bad iden=',i3,' (iden must be 0,1,2, or 3',
     +    ' as follows:')")
        write(6,"('iden = 0 -> leave species as on history ',
     +    '(most are mass mixing ratio)')")
        write(6,"('iden = 1 -> convert species to number densities',
     +    ' (cm3)')")
        write(6,"('iden = 2 -> convert species to number density ',
     +    'mixing ratios')")
        write(6,"('iden = 3 -> convert species to mass density ',
     +    '(gm/cm3)')")
        write(6,"('iden = 4 -> convert species to ppbv')")
        write(6,"('Will default to iden=1 (cm3)')")
        iden = 1
      end select
!
!   ionvel = flag for whether or not to add neutral atmos effects:
!     ionvel = 0 -> use UI,VI from the history (do not derive)
!     ionvel = 1 -> ExB only
!     ionvel = 2 -> ExB+unvn  (total ion velocity)
!     ionvel = 3 -> ExB+unvn, with E=0 (set electric potential = 0)
!     ionvel = 4 -> ExB+unvn, with unvnwn = 0.
!
      if (ionvel < 0 .or. ionvel > 5) then
        write(6,"('>>> Bad ionvel=',i2,' ionvel is single integer ',
     +    'flag for ion velocities as follows:')") ionvel
        write(6,"('ionvel = 0 -> use ion drifts from history',
     |    ' (do not derive)')")
        write(6,"('ionvel = 1 -> ExB only')")
        write(6,"('ionvel = 2 -> ExB+unvn  (total ion velocity)')")
        write(6,"('ionvel = 3 -> ExB+unvn, with E = 0')")
        write(6,"('ionvel = 4 -> ExB+unvn, with un,vn,wn = 0.')")
        write(6,"('ionvel = 5 -> WI due to ion diffusion')")
        stop 'ionvel'
      else
        select case (ionvel)
          case (1)
            write(6,"('ionvel = 1: ion velocities will be ExB only')")
          case (2)
            write(6,"('ionvel = 2: ion velocities will be ExB+unvn',
     +        ' (total ion velocity)')")
          case (3)
            write(6,"('ionvel = 3: ion velocities will be ExB+unvn',
     +        ' with E=0 (set electric potential = 0)')")
          case (4)
            write(6,"('ionvel = 4: ion velocities will be ExB+unvn',
     +        ' with un,vn,wn = 0')")
          case (5)
            write(6,"('ionvel = 5: ion velocities will be ',
     +        ' vertical ion drift from ion diffusion')")
        end select
      endif
!
! Validate lat slices:
      if (ipltlat > 0) then
        write(6,"('ipltlat = ',i1,' (will make latitude slices)')")
     +    ipltlat
        nlats = fseries(flats,mxslice,glat1,glat1+(nlat-1)*dlat)
        if (nlats>0) then
          call chlatsli
          if (ilat_log10>0.and.diffs) then
            write(6,"('  >>> NOTE: ilat_log10=',i2,' diffs=',l1,
     +        ' I am turning off ilat_log10')") ilat_log10,diffs
            ilat_log10 = 0
          endif
        else
          write(6,"('>>> WARNING INPUT: no valid latitudes for ',
     |      'latitude slices -- setting ipltlat off')") 
          ipltlat = 0
        endif
!
! clat_xaxes = 'LON' and/or 'LT' for longitude and/or local time x-axes
!   (default is both):
!   if clat_xaxes(1:2) == 'LON','LT', then draw lon x-axis and 
!     an extra slt axis below the lon axis (this is the default)
!   if clat_xaxes(1:2) == 'LON',' ' , then draw only lon x-axis
!   if clat_xaxes(1:2) == 'LT' ,' ' , then draw only slt x-axis
!   if clat_xaxes(1:2) == 'LT' ,'LON' , then this is same as 'LON','LT'
!     (i.e., cannot plot slt on x-axis w/ extra axis in lon)
!
        if (len_trim(clat_xaxes(1))==0.and.len_trim(clat_xaxes(2))==0) 
     +    then
          clat_xaxes = (/'LON     ','LT      '/)	! default
        endif
        do i=1,2
          if (len_trim(clat_xaxes(i)) > 0) then
            if (trim(clat_xaxes(i))=='GLON') clat_xaxes(i)='LON     '
            if (trim(clat_xaxes(i))=='glon') clat_xaxes(i)='LON     '
            if (trim(clat_xaxes(i))=='lon')  clat_xaxes(i)='LON     '
            if (trim(clat_xaxes(i))=='SLT')  clat_xaxes(i)='LT      '
            if (trim(clat_xaxes(i))=='slt')  clat_xaxes(i)='LT      '
            if (trim(clat_xaxes(i))=='lt')   clat_xaxes(i)='LT      '
            if (trim(clat_xaxes(i)) /= 'LON'.and.
     +          trim(clat_xaxes(i)) /= 'LT') then
              write(6,"(/'>>> input: unrecognized clat_xaxes=',a)")
     +          clat_xaxes
              write(6,"(10x,'clat_xaxes(1:2) is char*8 input ',
     +          'specifying LON and/or LT x-axes for lat slices.')") 
              stop 'clat_xaxes'
            endif
          endif
        enddo
        if (clat_xaxes(1)==clat_xaxes(2)) clat_xaxes(2) = '        '
        write(6,"('  clat_xaxes(1:2)=',a)") clat_xaxes
        if (flat_censlt /= spval)
     |    write(6,"('  flat_censlt=',f6.2)") flat_censlt
      endif	! ipltlat > 0
!
! Lon slices:
      if (ipltlon > 0) then
        write(6,"('ipltlon = ',i1,' (will make longitude slices)')")
     +    ipltlon
!
! Let chlons call fseries if necessary:
!       nlons = fseries(flons,mxslice,glon1,glon1+(nlon-1)*dlon)
!       nslts = fseries(fslts,mxslice,0.,23.99999)
        nlons = 0 ; nslts = 0
        do i=1,mxslice
          if (flons(i) /= spval) nlons = nlons+1
          if (fslts(i) /= spval) nslts = nslts+1 
        enddo
        if (nlons==0.and.nslts==0) then
          write(6,"('>>> WARNING INPUT: no valid longitudes or  local',
     |      ' times for longitude slices -- setting ipltlon off')") 
          ipltlon = 0
        endif
        if (nlons>0) then
          call chlons(flons,mxslice,1) 
          if (ilon_log10>0.and.diffs) then
            write(6,"('  >>> NOTE: ilon_log10=',i2,' diffs=',l1,
     +        ' I am turning off ilon_log10')") ilon_log10,diffs
            ilon_log10 = 0
          endif
          write(6,"('  flons = ')",advance="no")
          ii = 0 
          do i=1,mxslice
            if (flons(i)/=spval) then
              ii = ii+1
              if (flons(i)==zmflag) then
                write(6,"(' (Zonal Means), ')",advance="no")
              else if (isslt(flons(i),slt)) then
                write(6,"(2x,a10,2x)",advance="no") flons(i)
              else
                write(6,"(f8.2)",advance="no") flons(i)
              endif
              if (mod(ii,9)==0) write(6,"(/)")
            endif
          enddo
          if (len_trim(sendcdf_lons) > 0) 
     |      write(6,"(' Will make netcdf file of longitude slices.')")
        endif ! ipltlon > 0
!
! Check zprange:
! Note this will have to be moved to after history is read, because
!   pressure scale is unkown until after history is read.
!
        if (flon_zprange(1) /= spval .and. flon_zprange(2) /= spval)
     +    nzprange = 1	! temporary
!
! Check htscale:
        if (flon_htscale(1) /= spval .and. flon_htscale(2) /= spval.and.
     +      flon_htscale(3) /= spval) call chhtscale(flon_htscale)
        if (flon_htscale(1) /= spval .and. flon_htscale(2) /= spval.and.
     +      flon_htscale(3) /= spval) nhtscale = 1
!
! One or both of flon_zprange and flon_htscale must be set:
        if (nzprange == 0 .and. nhtscale == 0) then
          write(6,"('>>> Need valid flon_zprange and/or flon_htscale',
     +      ' -- turning ipltlon off')")
          ipltlon = 0
        endif
!
! Check for time-averaging:
! (delta time must be constant)
        if (ilon_timeave > 0) then
          if (ntimes <= 1) then
            write(6,"(/,'>>> INPUT: Must have more than one model time',
     |        ' to do time-averaging: ntimes=',i4)") ntimes
            write(6,"('Setting ILON_TIMEAVE=0')")
            ilon_timeave = 0
          else  
            idt = mtime_to_mins(mtimes(:,2))-mtime_to_mins(mtimes(:,1))
            do i=3,ntimes
              if ((mtime_to_mins(mtimes(:,i))-
     |             mtime_to_mins(mtimes(:,i-1))) /= idt) then
                write(6,"(/,'>>> INPUT: I will not time-average uneven',
     |            ' time intervals: mtimes = ',/,(5(i4,2i3,',')))")
     |            (mtimes(:,ii),ii=1,ntimes)
                stop 'timeav'
              endif
            enddo
            write(6,"('ilon_timeave=',i3,' -- will save time-averaged ',
     |        'longitude slices.')") ilon_timeave
          endif
        endif
        write(6,"(' ')")
      endif ! ipltlon > 0
!
! Amplitudes and phases:
      if (amphase>mxwave) then
        write(6,"('>>> amphase=',i2,' -- can only plot up to ',
     +    'wave number ',i2,' -- resetting amphase to ',i2)")
     +    mxwave,mxwave
        amphase = 4
      endif
!     if (amphase>0.and.diffs) then
!       write(6,"('>>> note: amplitudes/phases not available',
!    +    ' as difference fields -- setting amphase 0.')")
!       amphase = 0
!     endif
!
! lon slices of zm at selected slt: must be a time series. Calculate the
! the average across all ut's of longitudes corresponding to the selected slt,
! and plot latitude vs zp/ht:
!
      if (ipltzmslt > 0) then
        write(6,"('ipltzmslt = ',i1,' (will make zonal mean slices ',
     |    'at constant local time)')") ipltzmslt
        nslts = 0
        do i=1,mxslice
          if (zmslt(i) /= spval) then
            nslts = nslts+1
            if (zmslt(i) < 0..or.zmslt(i) > 24.) then
              write(6,"('>>> Bad local time zmslt=',e12.4,' (must be',
     |          ' between 0 and 24 (inclusive).')") zmslt(i)
              stop 'zmslt'
            endif
          endif
        enddo
        if (nslts==0) then
          write(6,"('>>> WARNING INPUT: am setting default zmslt=12.',
     |      ' for zonal mean local time slices.')")
          zmslt(1) = 12.
        endif
      endif ! ipltzmslt > 0
!
! Vertical profiles at selected locations (xyloc):
      if (ipltxyloc > 0) then
        write(6,"('ipltxyloc = ',i2,' (will make vertical profile',
     +    ' line plots at selected locations)')") ipltxyloc
	if(mod(count(xylocs /= spval),2) /= 0) then
          write(6,"('>>> input: incorrect locs count from xylocs')")
          stop 'INPUT'
	endif
        call chlats(xylocs(1,:),mxloc,1)
        call chlons(xylocs(2,:),mxloc,1) 
!       write(6,"('input: xylocs=')")
!       do i=1,mxloc
!         if (xylocs(1,i)/=spval.and.xylocs(2,i)/=spval)
!    +      write(6,"('lat,lon=',2f10.2)") xylocs(:,i)
!       enddo
        if (ixyloc_log10>0.and.diffs) then
          write(6,"('  >>> NOTE: ixyloc_log10=',i2,' diffs=',l1,
     +      ' I am turning off ixyloc_log10')") ixyloc_log10,diffs
          ixyloc_log10 = 0
        endif
      endif
!
! Validate map options:
!     if (ipltmaps==0.and.(map_global>0.or.map_polar>0.or.
!    +                     map_satview>0.or.map_mollweide>0))
!    +  ipltmaps = 1
      if (ipltmaps > 0) then
        write(6,"('ipltmaps=',i2)") ipltmaps
        write(6,"('  map_global    =',i2)") map_global
        write(6,"('  map_polar     =',i2)") map_polar
        if (map_polar > 0) then
          call chlats(fmap_polar_perimlat,mxperimlat,0)
        endif
        write(6,"('  map_satview   =',i2)") map_satview
        write(6,"('  map_mollweide =',i2)") map_mollweide
        write(6,"('  map_continents=',i2)") map_continents
        write(6,"('  map_log10=',i2)") map_log10

        if (fmap_latminmax(1)/=spval.or.fmap_latminmax(2)/=spval) then
          write(6,"(/,'>>> Sorry, fmap_latminmax input option is',
     |      ' not working at this time -- please remove this from',
     |      ' your namelist input file.',/)")
          call shutdown('fmap_latminmax not working')
        endif
!
! Maps of fields at hmf2 heights not available (this can be done
! for xyut line plots). Currently no practical limits on zpht:
!
        nzpht = fseries(fmap_zpht,mxzpht,-spval,spval)
        if (nzpht == 0) then
          write(6,"('>>> input: ipltmaps=',i2,' but no valid',
     +      ' fmap_zpht values were read -- turning maps off.')") 
     +      ipltmaps
          ipltmaps = 0
        endif
        write(6,"('  nzpht=',i3,' fmap_zpht=',/(8f9.3))") 
     +    nzpht,fmap_zpht(1:nzpht)
        write(6,"('  map_tn_unvn = ',i3)") map_tn_unvn
        write(6,"('  map_ht_unvn = ',i3)") map_ht_unvn
        if ((map_tn_unvn>0.or.map_ht_unvn>0.or.any(cfields=='UNVN'))
     +    .and.any(vn_scale>0.))
     +    write(6,"('  vn_scale (vector low,high,scale=',3f10.2)")
     +      vn_scale
        write(6,"('  map_ep_uivi = ',i3)") map_ep_uivi
        if ((map_ep_uivi>0.or.any(cfields=='UIVI'))
     +    .and.any(vi_scale>0.))
     +    write(6,"('  vi_scale (vector low,high,scale=',3f10.2)")
     +      vi_scale
      endif
!
! utvert contours:
      if (ipltutvert > 0) then
	if(mod(count(utvert_locs /= spval),2) /= 0) then
          write(6,"('>>> input: incorrect locs count from utvert')")
          stop 'INPUT'
	endif
        call chlats(utvert_locs(1,:),mxloc,1)
        call chlons(utvert_locs(2,:),mxloc,1) 
        if (iutvert_log10>0.and.diffs) then
          write(6,"('  >>> NOTE: iutvert_log10=',i2,' diffs=',l1,
     +      ' I am turning off iutvert_log10')") iutvert_log10,diffs
          iutvert_log10 = 0
        endif
        write(6,"('ipltutvert=',i2,' (ut vs vertical at selected',
     +    ' locations.')") ipltutvert
        if (len_trim(sendcdf_utvert) > 0) 
     |    write(6,"(' Will make netcdf file of ut vs vert slices.')")
      endif
!
! satut contours:
      if (ipltsatut > 0) then
        write(6,"('ipltsatut=',i2,' (ut vs vertical along ',
     +    ' satellite track)')") ipltsatut
      endif
!
! utlat contours:
      if (ipltutlat > 0) then
        call chlons(utlat_zphtlon(2,:),mxslice,1)
        if (iutlat_log10>0.and.diffs) then
          write(6,"('  >>> NOTE: iutlat_log10=',i2,' diffs=',l1,
     +      ' I am turning off iutlat_log10')") iutlat_log10,diffs
          iutlat_log10 = 0
        endif
        write(6,"('ipltutlat=',i2,' (ut vs latitude at selected',
     +    ' pressure (or height) and longitude (or slt or zm).')") 
     +    ipltutlat
      endif
!
! utlon contours:
      if (ipltutlon > 0) then
        call chlats(utlon_zphtlat(2,:),mxslice,1)
        if (iutlon_log10>0.and.diffs) then
          write(6,"('  >>> NOTE: iutlon_log10=',i2,' diffs=',l1,
     +      ' I am turning off iutlon_log10')") iutlon_log10,diffs
          iutlon_log10 = 0
        endif
        write(6,"('ipltutlon=',i2,' (ut vs longitude at selected',
     +    ' pressure (or height) and latitude.')")
     +    ipltutlon
      endif
!
! xyut line plots:
      if (ipltxyut > 0) then
        if (ntimes<=1) then
          write(6,"('>>> input: ntimes=',i4,': must have at least',
     +      ' 2 times to make xyut plots (ipltxyut)')") ntimes
          ipltxyut = 0
        else
	  if(mod(count(xyut_locs /= spval),2) /= 0) then
            write(6,"('>>> input: incorrect locs count from xyut')")
            stop 'INPUT'
	  endif
          call chlats(xyut_locs(1,:),mxloc,1)
          call chlons(xyut_locs(2,:),mxloc,1) 
          if (ixyut_log10>0.and.diffs) then
            write(6,"('  >>> NOTE: ixyut_log10=',i2,' diffs=',l1,
     +        ' I am turning off ixyut_log10')") ixyut_log10,diffs
            ixyut_log10 = 0
          endif
!
! Get doppler fields if either ixyut_doppler > 0 or idoppler > 0
! (idoppler was an old proc input parameter)
          write(6,"('ipltxyut=',i2,' (ut vs field at selected',
     +      ' locations and vertical levels.')") ipltxyut
          if (ixyut_doppler > 0 .or. idoppler > 0) then
            write(6,"('  ixyut_doppler=',i2,' idoppler=',i2,
     +      ' (will plot doppler t,u,v)')") ixyut_doppler,idoppler
            if (idoppler <= 0) idoppler = 1
            if (ixyut_doppler <= 0) ixyut_doppler = 1
          endif
          if (len_trim(sendcdf_xyut) > 0) 
     |      write(6,"(' Will make netcdf file of xyut data.')")
        endif
      endif
!
! Optional fixing of min,max and contour intervals:
! (see also setmnmxint, called by main tgcmproc after getflds)
      floop: do i=1,mxfproc
        n = 0      
        if (fmnmxint(1,i)/=spval) then	! user set
!
! fmnmxint must be in groups of 4 values ('field',cmin,cmax,cint):
          n = n+1   
          do ii=2,4
            if (fmnmxint(ii,i)/=spval) n = n+1
          enddo
          if (n /= 4) then
            write(6,"('>>> WARNING: fmnmxint must be given in ',
     +        'groups of 4 values:',/4x,'''FIELD'',cmin,cmax,cint',
     +        ' (i=',i3,')')") i
            fmnmxint(:,i) = spval
            cycle floop
          endif
!
! Fixed contour intervals must be >= 0.:
          if (fmnmxint(4,i) < 0.) then
            write(6,"('>>> WARNING: illegal contour interval=',
     +        e12.4,' specified for field ',a)") fmnmxint(4,i),
     +        fmnmxint(1,i)
            write(6,"(4x,'Contour intervals (fmnmxint(4,i)) must be ',
     +        '>= 0.',/4x,'(if zero, then graphics will choose).')")
            fmnmxint(:,i) = spval
            cycle floop
          endif
          write(6,"('Fixed min,max,int for field ',a,': ',
     +      3e12.4)") fmnmxint(1,i),fmnmxint(2:4,i)
        endif
      enddo floop
!
! cmnmxint same as fmnmxint, except it is char instead of real.
! This is mainly for AIX runs. See also sub setfmnmxint in tgcmproc.f.
      cloop: do i=1,mxfproc
        n = 0      
        if (len_trim(cmnmxint(1,i))>0) then	! user set
!
! cmnmxint must be in groups of 4 values ('field','cmin','cmax','cint'):
          n = n+1
          do ii=2,4
            if (len_trim(cmnmxint(ii,i))>0) n = n+1
          enddo
          if (n /= 4) then
            write(6,"('>>> WARNING: cmnmxint must be given in ',
     |        'groups of 4 strings:',/4x,'''FIELD'',''cmin'',''cmax'',
     |        ''cint'' (i=',i3,')')") i
            cmnmxint(:,i) = ' '
            cycle cloop
          endif
          read(cmnmxint(2,i),fmt = *) fmin
          read(cmnmxint(3,i),fmt = *) fmax
          read(cmnmxint(4,i),fmt = *) fcint
!
! Fixed contour intervals must be >= 0.:
          if (fcint < 0.) then
            write(6,"('>>> WARNING: illegal contour interval=',
     |        e12.4,' specified for field ',a)") fcint,cmnmxint(1,i)
            write(6,"(4x,'Contour intervals (fmnmxint(4,i)) must be ',
     |        '>= 0.',/,4x,'(if zero, then graphics will choose).')")
            cmnmxint(:,i) = ' '
            cycle cloop
          endif
          write(6,"('Fixed min,max,int for field ',a,': ',3e12.4)") 
     |      cmnmxint(1,i),fmin,fmax,fcint
        endif
      enddo cloop
!
! Scale factors:
! fscale must be given in pairs ('field',scale):
! If set by user, these are transferred to f(ixf)%scalefac by
!   setmnmxint, which is called by tgcmproc after getflds.
!
      floop1: do i=1,mxfproc
        if (fscale(1,i)/=spval) then	! user set
          if (fscale(2,i)==spval) then
            write(6,"('>>> WARNING: fscale must be given in ',
     +        'pairs (scale_factor is a float):',/4x,
     +        '''FIELD'',scale_factor')")
            write(6,"('             (i=',i3,')')") i
            fscale(:,i) = spval
            cycle floop1
          else
            write(6,"('Will set scale factor for field ',a,
     +        ' to ',1pe12.4)")  trim(float_to_str(fscale(1,i))),
     +        fscale(2,i)
          endif
        endif
      enddo floop1
!
! cscale:
      floop2: do i=1,mxfproc
        if (len_trim(cscale(1,i))>0) then	! user set
          if (len_trim(cscale(2,i))==0) then
            write(6,"('>>> WARNING: cscale must be given in ',
     +        'quoted pairs (char strings): field,scale_factor')")
            cscale(:,i) = ' '
            cycle floop2
          else
            read(cscale(2,i),fmt=*) fmin
            write(6,"('Will set scale factor for field ',a,
     +        ' to ',1pe12.4)")  trim(cscale(1,i)),fmin
          endif
        endif
      enddo floop2
!
! If making vis5d output file, define v5d (type(vis5d) in module mk_v5d) 
!   from defaults and user input:
!   (v5d%zprange is set in setflev (getflds.f) when zp of hist is known)
!
      v5d%sendv5d = ' '
      v5d%sendmsv5d =  ' '
      if (len_trim(sendv5d)>0.or.len_trim(sendms_v5d)>0) then
#ifdef SUN
        write(6,"('>>> v5d not available on Sun.')")
        stop 'v5d on Sun'
#else
        v5d%flnm = ' '
        v5d%sendv5d = sendv5D
        v5d%sendmsv5d = sendms_v5D
        v5d%vtype = 'reg_zp  '
        if (len_trim(v5d_vtype)>0) then
          if(trim(v5d_vtype)=='reg_zp'.or.trim(v5d_vtype)=='ireg_mb'.or.
     +       trim(v5d_vtype)=='REG_ZP'.or.trim(v5d_vtype)=='IREG_MB')
     +      then
            v5d%vtype=v5d_vtype
          else
            write(6,"('>>> input: bad v5d_vtype: ',a)")
            write(6,"('    v5d may be ''reg_zp'' (regular ZP) ',
     +        'or ''ireg_mb'' (irregular millibars)')")
            stop 'v5d_vtype'
          endif
        endif
!
! Irregular mb vertical coord system does not work (pressures get too 
! small) -- must use generic regular (zp):
        if (trim(v5d%vtype)=='ireg_mb'.or.trim(v5d%vtype)=='IREG_MB')
     +    then
          write(6,"('>>> input: v5d_vtype=ireg_mb does not work.')")
          write(6,"('    changing v5d_type to reg_zp')")
          v5d%vtype = 'reg_zp  '
        endif
!
        v5d%fix_slt = spval
        if (v5d_slt/=spval) then
          if (v5d_slt>=0..and.v5d_slt<=24.) then
            v5d%fix_slt = v5d_slt
          else
            write(6,"('>>> input: bad v5d_slt =',f10.2,
     +        '(must be >= 0. and <= 24.')") v5d_slt
            stop 'v5d_slt'
          endif
        endif
        v5d%zprange = spval	! will be set in setflev (getflds.f)
#endif
      endif      
      if (len_trim(sendv5d) > 0.or.len_trim(sendms_v5d) > 0) then
        write(6,"('sendv5d = ',a)") trim(sendv5d)
        write(6,"('  sendms_v5d = ',a)") trim(sendms_v5d)
        if (v5d_zprange(1)/=spval.and.v5d_zprange(2)/=spval)
     +    write(6,"('  v5d_zprange=',2f8.2)") v5d_zprange
        if (v5d_slt/=spval) then
          write(6,"('  v5d_slt=',f8.2)") v5d_slt
        else
          write(6,"('  v5d_slt not set (data will NOT be fixed ',
     +      'by local time)')")
        endif
      endif
!
      if (iplot <= 0) then
        write(6,"('  iplot = ',i3,' --> plots will NOT be made.')")iplot
        write(6,"('  NOTE: values of sendcgm, sendps, etc. will be ',
     |    'IGNORED because iplot=',i3)") iplot
      else
        write(6,"('  iplot = ',i3,' --> plots WILL be made.')") iplot
      endif
!
! Verify outplt:
!
      do i=1,3
        if (len_trim(outplt(i)) > 0 .and.
     |    (trim(outplt(i)) /= 'cgm' .and. trim(outplt(i)) /= 'ps')) then
          write(6,"('>>> Input: unknown value for OUTPLT = ',a)")
     |      outplt(i)
          write(6,"('    OUTPLT must be one or both of cgm or ps')")
          stop 'OUTPLT'
        endif
      enddo
!
! Remote paths for output:
! If a sendcgm path was given but outplt does not include 'cgm', then
!   assume user wanted cgm and add it to outplt. Do same for sendps.
!
      if (len_trim(sendcgm) > 0) then
        write(6,"('sendcgm = ',a)") trim(sendcgm)
        if (.not.any(outplt=='cgm     ').and..not.any(outplt=='CGM')) 
     +    then
          loop_outplt1: do i=1,3
            if (len_trim(outplt(i))==0) then
              outplt(i) = 'cgm     '
              write(6,"('  (added cgm to outplt)')")
              exit loop_outplt1
            endif
          enddo loop_outplt1
        endif
      endif
      if (len_trim(sendps ) > 0) then
        write(6,"('sendps  = ',a)") trim(sendps)
        if (.not.any(outplt=='ps      ').and..not.any(outplt=='PS')) 
     +    then
          loop_outplt2: do i=1,3
            if (len_trim(outplt(i))==0) then
              outplt(i) = 'ps      '
              write(6,"('  (added ps to outplt)')")
              exit loop_outplt2
            endif
          enddo loop_outplt2
        endif
      endif
      if (len_trim(sendcdf) > 0.or.len_trim(flnm_cdf) > 0) then
        write(6,"('sendcdf = ',a)") trim(sendcdf)
        if (any(cdf_htscale /= spval)) then
          if (cdf_htscale(1)==cdf_htscale(2)) then ! 1 height is allowed
            cdf_htscale(3) = 1.
          else
            call chhtscale(cdf_htscale)
            if (any(cdf_htscale==spval)) then
              write(6,"(/,'>>> INPUT: bad cdf_htscale')")
              stop 'cdf_htscale'
            endif
          endif
          write(6,"('  cdf_htscale = ',3f10.2)") cdf_htscale
        endif
!
! Note cdf_zprange will be validated by sub setflev in getflds.F.
        if (any(cdf_zprange /= spval)) then
          if (.not.all(cdf_zprange /= spval)) then
            write(6,"(/,'>>> INPUT: bad cdf_zprange=',2e12.4)")
     |        cdf_zprange
            stop 'cdf_zprange' 
          endif
        endif
!
! Cannot do both zp and ht netcdf file:
        if (any(cdf_htscale /= spval).and.
     |      any(cdf_zprange /= spval)) then
          write(6,"(/,'>>> Input: cannot write both pressure and ',
     |      'height-interpolated fields',/,'to the same netcdf file',
     |      ' in the same run.')")
          write(6,"('Please provide either cdf_htscale or cdf_zprange',
     |      ' but not both.')")
          stop 'netcdf vertical'
        endif
!
! Zonal means of fields on netcdf file:
        if (cdf_zonalmean /= spval) then
          if (cdf_zonalmean > 0.) then
            cdf_zonalmean = 1. 
          else
            cdf_zonalmean = 0.
          endif 
          write(6,"('input: set cdf_zonalmean = ',f4.1)") cdf_zonalmean
        else                   ! default is no zonal means
          cdf_zonalmean = 0.
        endif
      endif
      if (len_trim(sendsat) > 0) then
        write(6,"('sendsat = ',a)") trim(sendsat)
      endif
      iwrdat = 0
!
! Build ascii data file:
! Default disk file name flnm_dat and senddat are blank.
!
      if (len_trim(senddat) > 0 .or. len_trim(flnm_dat) > 0) then
        if (len_trim(senddat) > 0) write(6,"('senddat = ',a)") 
     |    trim(senddat)
        if (len_trim(flnm_dat) > 0) write(6,"('flnm_dat = ',a)") 
     |    trim(flnm_dat)
! btf 11/29/07: Eliminate append to dat file option (iwrdat=2)
        iwrdat = 1 
!
! Set ludat (proc.f) to avoid possible conflict with crayopen:
!       ludat = nextlu()
        ludat = 98 
!
! If senddat was given, but not flnm_dat, set flnm_dat to default:
        if (len_trim(flnm_dat) == 0) flnm_dat = 'tgcmproc.dat'
!
! Remove any pre-existing local ascii data disk file 
! (this will exist only when proc is run interactive):
!       ios = iunlink(flnm_dat,1)
!       if (ios==0) then
!         write(6,"(/'Removed pre-existing data output file ',
!    +      a,' from the disk'/)") trim(flnm_dat)
!       endif
!
        write(6,"('input: flnm_dat    =',a)") trim(flnm_dat)
        write(6,"('       senddat     =',a)") trim(senddat)
        write(6,"('       iwrdat      =',i3)") iwrdat
!
! Open data file:
        open(ludat,file=flnm_dat,status='REPLACE')
        write(6,"('input: opened ascii data file ',a,' ludat=',i3)")
     |    trim(flnm_dat),ludat
      endif ! senddat and/or flnm_dat are non-blank
!
! Build netcdf data file:
      if (len_trim(sendcdf) > 0 .or. len_trim(flnm_cdf) > 0) then

        if (len_trim(sendcdf) > 0) write(6,"('sendcdf = ',a)") 
     |    trim(sendcdf)
        if (len_trim(flnm_cdf) > 0) write(6,"('flnm_cdf = ',a)") 
     |    trim(flnm_cdf)

        if (len_trim(flnm_cdf) == 0) flnm_cdf = 'tgcmproc.nc'

      endif
!
! Build xdr data:
      iwrxdr = 0
      if (len_trim(sendxdr) > 0.or.len_trim(sendms_xdr) > 0) then
        write(6,"('sendxdr = ',a)") trim(sendxdr)
        write(6,"('sendms_xdr = ',a)") trim(sendms_xdr)
        iwrxdr = 1
        luxdr = nextlu()
!
! Remove any pre-existing local xdr disk file 
! (this will exist only when proc is run interactive):
        ios = iunlink(flnm_xdr,1)
        if (ios==0) then
          write(6,"(/'Removed pre-existing xdr output file ',
     +      a,' from the disk'/)") flnm_xdr
        endif
      endif
!
! Bougher format:
      iwrdat_bf = 0
      if (len_trim(senddat_bf) > 0) then
        write(6,"('senddat_bf = ',a)") trim(senddat_bf)
        iwrdat_bf = 2
        ludat_bf = nextlu()     ! nextlu may not be working properly
!
! Remove any pre-existing local ascii data disk file (Bougher format)
! (this will exist only when proc is run interactive):
        ios = iunlink(flnm_bf_dat,0)
        if (ios==0) then
          write(6,"(/,'Removed pre-existing data output file ',
     +      a,' from the disk',/)") trim(flnm_bf_dat)
        endif
      endif
!
!     iplot = 1
!
! If sendcgm and sendps are not provided, then do not make plots.
! If none of sendcgm, sendps, senddat, sendxdr, sendcdf,, sendsat, or sendv5d 
!   are provided, stop.
!
!     if (len_trim(sendcgm)==0.and.len_trim(sendps)==0) then
!       if (iwrxdr==0.and.iwrdat==0.and.len_trim(sendcdf)==0.and.
!    |      len_trim(sendv5d)==0.and.len_trim(sendsat)==0.and.
!    |      len_trim(sendcdf_lons)==0.and.len_trim(sendcdf_xyut)==0)then
!         write(6,"(/'>>> PLEASE SPECIFY PLOT OR DATA OUTPUT PATHS:')")
!         write(6,"('Need one or more of the following to be set:')")
!         write(6,"('  sendcgm = machine:path for cgm graphics',
!    +      ' output file.')")
!         write(6,"('  sendps  = machine:path for ps graphics output',
!    +      ' file.')")
!         write(6,"('  senddat = machine:path for ascii data output',
!    +      ' file.')")
!         write(6,"('  sendxdr = machine:path for xdr data output',
!    +      ' file.')")
!         write(6,"('  sendcdf = machine:path for netcdf data output',
!    +      ' file.')")
!         write(6,"('  sendsat = machine:path for netcdf satellite',
!    +      ' file.')")
!         write(6,"('  sendv5d = machine:path for vis5d data output',
!    +      ' file.')")
!         stop
!       else
!         write(6,"(/'NOTE: NO PLOTS WILL BE PRODUCED THIS RUN.')")
!         write(6,"('      (only output data files will be made.)')")
!         write(6,"('      This is because neither sendcgm nor sendps ',
!    +      'were provided.'/)") 
!         iplot = 0
!       endif
!     endif
!
! File server path to NASA AMES lower boundaries for mars tgcm:
      if (len_trim(mgcm_fspath) > 0) then
        write(6,"('mgcm_fspath = ',a)") trim(mgcm_fspath)
      endif
!
! File server path to vtgcm N2:
      if (len_trim(vtgcm_fspath) > 0) then
        write(6,"('vtgcm_fspath = ',a)") trim(vtgcm_fspath)
      endif
!
! Get list of requested locations (reqlocs(2,mxloc)):
      n = 0
      do i=1,mxloc
        if (ipltxyloc > 0) then
          if (xylocs(1,i)/=spval.and.xylocs(2,i)/=spval) then
            if (.not.findloc(reqlocs,mxloc,xylocs(1,i),xylocs(2,i)))then 
              n = n+1
              reqlocs(:,n) = xylocs(:,i) 
            endif
          endif
        endif
        if (ipltutvert > 0) then
          if (utvert_locs(1,i)/=spval.and.utvert_locs(2,i)/=spval)then
            if (.not.findloc(reqlocs,mxloc,
     +        utvert_locs(1,i),utvert_locs(2,i))) then 
              n = n+1
              reqlocs(:,n) = utvert_locs(:,i) 
            endif
          endif
        endif
        if (ipltxyut > 0) then
          if (xyut_locs(1,i)/=spval.and.xyut_locs(2,i)/=spval) then
            if (.not.findloc(reqlocs,mxloc,
     +        xyut_locs(1,i),xyut_locs(2,i))) then 
              n = n+1
              reqlocs(:,n) = xyut_locs(:,i) 
            endif
          endif
        endif
      enddo 
!
! If iohglb>0, then mkderived will calculate oh fields at global grid, 
!   otherwise, it will calculate oh fields only at reqlocs locations. 
!   Note iohglb has been set to 1 if either of the following are true:
! 1. if any requested locations are zonal means, global means, or local time
! 2. if ipltmaps>0, ipltlon>0, or ipltlat>0
!
      iohglb = 0
      if (n > 0) then
        write(6,"('The following ',i3,' locations have been requested:'
     +    )") n
        do i=1,n
          if (isslt(reqlocs(2,i),slt)) then
            write(6,"('Loc ',i3,': lat=',f10.2,' (lon at local time ',
     +        f5.1,')')") i,reqlocs(1,i),slt
            iohglb = 1
          elseif (reqlocs(1,i)/=zmflag.and.reqlocs(2,i)/=zmflag) then
            write(6,"('Loc ',i3,': lat,lon=',2f10.2)") i,reqlocs(:,i)
          elseif (reqlocs(1,i)/=zmflag.and.reqlocs(2,i)==zmflag)then
            write(6,"('Loc ',i3,': lat,lon=',2f10.2,' (zonal means)')") 
     +        i,reqlocs(:,i)
            iohglb = 1
          elseif (reqlocs(1,i)==zmflag.and.reqlocs(2,i)==zmflag) then
            write(6,"('Loc ',i3,': lat,lon=',2f10.2,' (global means)')") 
     +        i,reqlocs(:,i)
            iohglb = 1
          endif
        enddo
      endif
      if (ipltmaps>0.or.ipltlat>0.or.ipltlon>0.or.
     +    ipltutlat>0.or.ipltutlon>0) iohglb=1
      write(6,"('iohglb=',i2)") iohglb
!
      if (ilab_hq > 0) label_quality = 'high    '
      if (label_quality(1:3) /= 'low'.and.
     |    label_quality(1:3) /= 'med'.and.
     |    label_quality(1:3) /= 'hig') then
        write(6,"(/,'>>> Bad label_quality = ',a)") label_quality
        write(6,"(/,'>>> label_quality must be low, medium, or high.')")
        stop 'label_quality'
      endif
!
      write(6,"(20('-'),' End input verification ',20('-'))")
      return
 900  continue
!
! Error in namelist read:
! (on the Crays, ios=1324 means unrecognized keyword in namelist input)
!
      write(6,"(/72('>')/'ERROR in namelist read of user inputs: ',
     +  ' lu=',i2,' ios=',i5/)") luin,ios
      backspace(luin)
      read(luin,"(a)") rec80
      write(6,"('This might be the result of an unrecognized ',
     +  'or misspelled keyword in the input file.')")
      write(6,"('Please check your input file in the vicinity ',
     +  'of the following line:')")
      write(6,"(/a/)") rec80
      write(6,"(72('<')/)")
      close(luin)
      stop 'input'
      end subroutine getinp
!-------------------------------------------------------------------
      integer function mtime_to_min(iday,ihr,min)
      integer,intent(in) :: iday,ihr,min
      mtime_to_min = iday*1440 + ihr*60 + min
      return
      end function mtime_to_min
!-------------------------------------------------------------------
      integer function chmtimes(mtimes,mxtms)
      use proc,only: mtimeflag
!     
! Args:
      integer,intent(in) :: mxtms
      integer,intent(inout) :: mtimes(3,mxtms)
!
! Locals:
      integer :: idmin,ntms,i,min,maxmin
      logical isshort
      logical,external :: int_is_str
!
! Check mtimes as given by namelist.
! If isshort=T, this means mtimes has been given in the form:
!   "id0,ih0,im0,'to',id1,ih1,im1,'by',idelmin", where:
!   id0,ih0,im0 = beginning model time, id1,ih1,im1 = ending model time, and
!   idelmin is delta time in minutes.
!   e.g:  mtimes = 0,0,0,'to',1,0,0,'by',60
! In this case (isshort), expand to proper mtimes
! If not isshort, just check validity of each component of each model time
! Return chmtimes as number of valid times.
!
      isshort = .false.
      if ((int_is_str(mtimes(1,2),'to').or.
     +     int_is_str(mtimes(1,2),'TO')).and.
     +    (int_is_str(mtimes(2,3),'by').or.
     +     int_is_str(mtimes(2,3),'BY'))) isshort = .true.
      if (mtimes(1,2)==mtimeflag.and.mtimes(2,3)==mtimeflag)
     |  isshort = .true.
      if (isshort) then
        if (mtimes(1,1) < 0 .or.mtimes(2,2) < 0) then
          write(6,"('WARNING: Model day < 0 not allowed in implied ',
     +      'times loop (mtimes).')")
          if (mtimes(1,1) < 0) mtimes(1,1) = -mtimes(1,1)
          if (mtimes(2,2) < 0) mtimes(2,2) = -mtimes(2,2)
        endif
        idmin = mtimes(3,3) 
        maxmin = mtime_to_min(mtimes(2,2),mtimes(3,2),mtimes(1,3))
        ntms = (maxmin - 
     +    mtime_to_min(mtimes(1,1),mtimes(2,1),mtimes(3,1))) / idmin + 1
        if (ntms.gt.mxtms) then
          write(6,"(/'>>> WARNING: too many times. ntms=',i3,
     +      ' mxtms=',i3,' -- number of times will be limited to ',
     +      'mxtms=',i3)") ntms,mxtms,mxtms
        endif
        do i=1,mxtms
          if (i.gt.1) then
            min=mtime_to_min(mtimes(1,i-1),mtimes(2,i-1),mtimes(3,i-1))+
     +        idmin
            call min_to_mtime(min,mtimes(1,i),mtimes(2,i),mtimes(3,i))
            if (mtime_to_min(mtimes(1,i),mtimes(2,i),mtimes(3,i)).gt.
     +          maxmin.or.i.eq.mxtms) then
              chmtimes = i-1
              if (i.eq.mxtms) chmtimes = i
              return
            endif
          endif
        enddo
      else
        do i=1,mxtms
          if (mtimes(1,i).ne.ispval.and.mtimes(2,i).ne.ispval.and.
     +        mtimes(3,i).ne.ispval) then
!
! Allow mtimes(1,i) < 0 to uncondionally accept next history:
!           if (mtimes(1,i).lt.0.or.mtimes(1,i).gt.367) then
            if (mtimes(1,i).gt.367) then
              write(6,"('>>> Warning: bad model day=',i5)") mtimes(1,i)
            endif 
            if (mtimes(2,i).lt.0.or.mtimes(2,i).gt.23) then
              write(6,"('>>> Warning: bad model hour=',i5)") mtimes(2,i)
            endif
            if (mtimes(3,i).lt.0.or.mtimes(3,i).gt.59) then
              write(6,"('>>> Warning: bad model minute=',i5)") 
     +          mtimes(3,i)
            endif
          endif
        enddo
        chmtimes = 0
        do i=1,mxtms
          if (mtimes(1,i).ne.ispval.and.mtimes(2,i).ne.ispval.and.
     +        mtimes(3,i).ne.ispval) chmtimes = chmtimes+1
        enddo
      endif
      return
      end function chmtimes
!-------------------------------------------------------------------
      integer function rmcomments(lu,comcharin,echo)
      implicit none
!
! Read input lines from unit lu. If current line contains the comment
!   character comcharin, strip the line from position of comchar to end,
!   and write any remaining line to a new unit. If no comment in current 
!   line, write entire line to new unit. 
! Return new unit, rewound (e.g., ready to be read by namelist).
! If echo > 0, echo output lines to stdout.
! If comcharin is ' ', then default comment char is ';'
!
! Some inputs may be changed for compatability w/ old proc and 
!   f90 namelist restrictions, e.g., beginning flag for namelist,
!   quotes around cfields field names.
!
! Args:
      integer,intent(in) :: lu,echo
      character(len=1),intent(in) :: comcharin
! Local:
      character(len=1) :: comchar
      logical isopen
      integer :: i,lens,ios,compos,nline,nlcfields
      character(len=10000) :: line
      character(len=64) :: newcfields(30)
! Externals:
      integer,external :: nextlu
!     logical,external :: iscfields
!
      if (lu <= 0) then
        write(6,"('>>> rmcomments: bad input lu=',i5)") lu
        rmcomments = -1
        return
      endif
      if (len_trim(comcharin) > 0) then
        comchar = comcharin
      else
        comchar = ';'
        write(6,"('rmcomments: using default semicolon as ',
     +    'comment character.')")
      endif
      inquire(unit=lu,opened=isopen)
      if (.not.isopen) then
        open(unit=lu,iostat=ios)
        if (ios /= 0) then
          write(6,"('>>> WARNING rmcomments: error opening input',
     +      ' file with unit lu=',i2,' ios=',i5)") lu,ios
          rmcomments = -1
          return
        endif
      endif
      rmcomments = nextlu()
!     rewind lu
      nline = 0
      read_loop: do
        line = ' '
        read(lu,"(a)",iostat=ios) line
        if (ios > 0) then
          write(6,"('>>> rmcomments: error reading from input',
     +      ' unit lu=',i3,' at line ',i5)") lu,nline
          return
        endif
        if (ios < 0) exit read_loop	! eof
        nline = nline+1
!
! Remove line if it has only "E" in column 1 (this was an
! old "Echo" directive from f77/cray namelist):
!
        if (line(1:1)=='E'.and.trim(line)=='E') cycle read_loop
!
! Use only non-commented part of line:
!
        compos = index(line,comchar)
        if (compos == 1) cycle read_loop
        if (compos > 0) line = line(1:compos-1)
        if (len_trim(adjustl(line))==0) cycle read_loop
!
! Check for old proc style input group and if found, use new input group.
! This enables new proc to read old proc input files.
!
        if (index(line,'&input_sut')>0.or.index(line,'&input_mut')>0)
     +    write(line,"('&proc_input')")
!
! If current line is cfields, add quotes around field names:
!
        if (iscfields(line,nline)) then
          call fixcfields(line,newcfields,nlcfields)
          do i=1,nlcfields
            write(rmcomments,"(a)") trim(newcfields(i))
            if (echo > 0) write(6,"(a)") 
     +        newcfields(i)(1:len_trim(newcfields(i)))
          enddo
        else
          write(rmcomments,"(a)") trim(line)
          if (echo > 0) write(6,"(a)") line(1:len_trim(line))
        endif
      enddo read_loop  
      rewind rmcomments
      return      
      end function rmcomments
!-------------------------------------------------------------------
      logical function iscfields(line,nline)
      character(len=*),intent(in) :: line
      integer,intent(in) :: nline
      logical,save :: wascfields=.false.
!
      iscfields = .false.
      if (index(line,'cfields') > 0) then
        iscfields = .true.
        wascfields = .true.
        return
      elseif (wascfields) then
        if (index(line,'=')==0.and.index(line,'&')==0) then
          iscfields = .true.
        else
          wascfields = .false.
        endif
      endif
      end function iscfields
!-------------------------------------------------------------------
      subroutine fixcfields(line,newlines,nlout)
      implicit none
!
! Given cfields input line (or continuation of a cfields line), add 
!   quotes around the field names if necessary. This is for 
!   compatability w/ oldproc f77 namelist, which allowed unquoted 
!   char input. New f90 namelist does not allow unquoted char input.
! (line is inout only so any tabs can be replaced w/ blanks, and
!  leading blanks can be removed)
!
! Args:
      character(len=*),intent(inout) :: line
      integer,intent(out) :: nlout
      character(len=*),intent(out) :: newlines(30)
!
! Locals:
      integer :: icfpos,lens,lenf,icomma1,icomma2,nf,i,ii,itab
      character(len=16) :: fnames(30)
      character(len=1) :: comma = ','
!
! Replace any tabs with blanks:
      do
        itab = index(line,'	')
        if (itab > 0) line(itab:itab)=' '
        if (itab==0) exit
      enddo
      line = adjustl(line)
!
! icfpos == 0 when line is a continuation line of cfields, otherwise
!   line is first cfields line.
!
      icfpos = index(line,'cfields')
      lens = len_trim(line)
      icomma1 = 0 ; icomma2 = 0 ; nf = 0
!
! Build list of fnames which are separated by commas: 
!
      do i=icfpos+1,lens
        if (line(i:i)==comma) then
          if (icomma1==0) then
            icomma1 = i
            if (icfpos > 0) then
              do ii=i,1,-1
                if (line(ii:ii)=='=') then
                  nf = nf+1
                  fnames(nf) = ' '
                  lenf = i-ii-1
                  if (lenf > 16) lenf = 16
                  fnames(nf)(1:lenf) = line(ii+1:i-1)
                  fnames(nf) = adjustl(fnames(nf))
                endif
              enddo
            else
              nf = nf+1
              ii_loop: do ii=1,lens
                if (line(ii:ii)/=' ') exit ii_loop
              enddo ii_loop
              fnames(nf) = line(ii:i-1)
              fnames(nf) = adjustl(fnames(nf))
            endif
          else
            icomma2 = i
            nf = nf+1
            fnames(nf) = ' '
            lenf = icomma2-icomma1-1
            if (lenf > 16) lenf = 16
            fnames(nf)(1:lenf) = line(icomma1+1:icomma2-1)
            fnames(nf) = adjustl(fnames(nf))
            icomma1 = i; icomma2 = 0
          endif
        endif	! comma in position i
      enddo
!
! No commas were found -> there is only one field name
      if (icomma1==0) then
        if (icfpos > 0) then
          do ii=lens,1,-1
            if (line(ii:ii)=='=') then
              nf = nf+1
              fnames(nf) = ' '
              lenf = lens-ii
              if (lenf > 16) lenf = 16
              fnames(nf)(1:lenf) = line(ii+1:lens)
              fnames(nf) = adjustl(fnames(nf))
            endif
          enddo
        else
          nf = nf+1
          fnames(nf) = line(1:lens)
          fnames(nf) = adjustl(fnames(nf))
        endif
      endif
!
! No trailing comma on last field name:
      if (icomma1 > 0 .and. icomma1 < lens) then
        nf = nf+1
        fnames(nf) = ' '
        fnames(nf) = line(icomma1+1:lens) 
        fnames(nf) = adjustl(fnames(nf))
      endif
!
! Insure single quotes around each fname:
      do i=1,nf
        if (fnames(i)(1:1) /= '''') then 
          if (fnames(i)(1:1)=='"') then
            fnames(i)(1:1) = ''''
          else
            fnames(i) = ''''//trim(fnames(i))
          endif
        endif
        lenf = len_trim(fnames(i))
        if (fnames(i)(lenf:lenf) /= '''') then
          if (fnames(i)(lenf:lenf) == '"') then
            fnames(i)(lenf:lenf) = ''''
          else  
            fnames(i) = fnames(i)(1:lenf)//''''
          endif
        endif
      enddo
!
! Write new line(s):
      newlines = ' '
      nlout = 1
      if (icfpos > 0) newlines(nlout) = 'cfields='
      do i=1,nf
        if (len_trim(newlines(nlout))+len_trim(fnames(i))+2 > 
     +    len(newlines(nlout))) then
          nlout = nlout+1
        endif
        if (len_trim(newlines(nlout))==0.and.(nlout>1.or.icfpos==0)) 
     +    then
          newlines(nlout) = '  '//trim(newlines(nlout))//trim(fnames(i))
        else
          newlines(nlout) = trim(newlines(nlout))//trim(fnames(i))
        endif
        newlines(nlout) = trim(newlines(nlout))//comma
!       write(6,"('fixcfields: i=',i2,' nlout=',i2,' newlines(nlout)=',
!    +    /2x,'$'a,'$')") i,nlout,newlines(nlout)   
      enddo 
      end subroutine fixcfields
!-------------------------------------------------------------------
      subroutine chlatsli
!
! Check lat slice inputs:
!
! Locals:
      integer :: nzprange=0, nhtscale=0 
!
! Check selected latitudes:
      call chlats(flats,mxslice,1)
!
! Check zprange:
! Note this will have to be moved to after history is read, because
!   pressure scale is unkown until after history is read.
!
      if (flat_zprange(1) /= spval .and. flat_zprange(2) /= spval)
     +  nzprange = 1	! temporary
!     if (flat_zprange(1) /= spval .and. flat_zprange(2) /= spval) then
!       call chzp(flat_zprange,2)
!       if (flat_zprange(1) >= flat_zprange(2)) then
!         write(6,"('>>> bad flat_zprange: ',2f8.2)")
!    +      flat_zprange
!         write(6,"('    (flat_zprange(1) must be < flat_zprange(2))')") 
!         flat_zprange = spval
!       else
!         nzprange = 1
!       endif
!     endif
!
! Check htscale:
      if (flat_htscale(1) /= spval .and. flat_htscale(2) /= spval .and.
     +    flat_htscale(3) /= spval) call chhtscale(flat_htscale)
      if (flat_htscale(1) /= spval .and. flat_htscale(2) /= spval .and.
     +    flat_htscale(3) /= spval) nhtscale = 1
!
! One or both of flat_zprange and flat_htscale must be set:
      if (nzprange == 0 .and. nhtscale == 0) then
        write(6,"('>>> Need valid flat_zprange and/or flat_htscale',
     +    ' -- turning ipltlat off')")
        ipltlat = 0
      endif
!
      if (flat_censlt /= spval) then
        if (flat_censlt < 0..or.flat_censlt > 24.) then
          write(6,"('chlatsli: bad flat_censlt=',f6.2,' (must be ',
     |      'within 0. and 24.')") flat_censlt
        else
!         write(6,"('chlatsli: flat_censlt=',f6.2)") flat_censlt
        endif
      endif
      return
      end subroutine chlatsli
!-------------------------------------------------------------------
      subroutine chhtscale(htscale)
!
! Validate height scale (bottom, top, delta):
!
      real,intent(inout) :: htscale(3)
!
      if (htscale(1) < 0. .or. htscale(2) < 0. .or. htscale(3) < 0.) 
     +  then
        write(6,"('>>> bad htscale (bottom,top,delta)=',3f9.2,
     +   ' (all 3 values must be > 0.)')") htscale
        htscale = spval
      endif
      if (htscale(1) >= htscale(2)) then
        write(6,"('>>> bad htscale (bottom,top,delta)=',3f9.2,
     +    ' (htscale(1) must be < htscale(2))')") htscale
        htscale = spval
      endif 
      if (htscale(3) > htscale(2)) then
        write(6,"('>>> bad htscale (bottom,top,delta)=',3f9.2,
     +    ' (htscale(3) must be < htscale(2))')") htscale
        htscale = spval
      endif
      return
      end subroutine chhtscale
!-------------------------------------------------------------------
      subroutine chzp(zp,nzp)
      use proc,only: gcmlev,npress
!
! Validate selected pressures zp -- set to nearest
! log pressure grid point.
!
      integer,intent(in) :: nzp
      real,intent(inout) :: zp(nzp)
      integer :: k,kk
      integer,external :: ixfind
!
      lev_loop: do k=1,nzp
        do kk=1,npress
          if (zp(k) == gcmlev(kk)) cycle lev_loop
        enddo
        if (zp(k) < gcmlev(1)) then
          write(6,"('Note: changed selected pressure ',f8.2,
     +      ' to ',f8.2)") zp(k),gcmlev(1)
          zp(k) = gcmlev(1)
        elseif (zp(k) > gcmlev(npress)) then
          write(6,"('Note: changed selected pressure ',f8.2,
     +      ' to ',f8.2)") zp(k),gcmlev(npress)
          zp(k) = gcmlev(npress)
        else
          kk = ixfind(gcmlev,npress,zp(k),dlev)
          write(6,"('Note: changing selected pressure level ',
     +      f9.2,' to nearest tgcm log pressure grid point ',f9.2)")
     +      zp(k),gcmlev(kk)
          zp(k) = gcmlev(kk)
        endif
      enddo lev_loop
      return
      end subroutine chzp      
!-------------------------------------------------------------------
      subroutine chlons(glons,nlons,ichange)
      use proc,only: gcmlon,nlon
!
! Validate selected longitudes glons -- if ichange > 0 then change 
!   lons if necessary to nearest tgcm grid longitude.
!
! Args:
      integer,intent(in) :: nlons,ichange
      real(kind=8),intent(inout) :: glons(nlons)
!
! Locals:
      integer :: i,ii,nseries
      real :: slt, glons_tmp
      logical :: isshort
!
! Externals:
      character(len=16),external :: float_to_str
      integer,external :: ixfind,fseries
      logical,external :: isslt,float_is_str
      real,external :: convlon
!
! Check for "short form" and expand w/ fseries if necessary:
      isshort = .false.
      if ((float_is_str(glons(2),'to').or.
     |     float_is_str(glons(2),'TO')).and.
     |    (float_is_str(glons(4),'by').or.
     |     float_is_str(glons(4),'BY'))) isshort = .true.
      if (isshort) nseries = fseries(glons,nlons,gcmlon(1),gcmlon(nlon))
!
      slice_loop: do i=1,nlons
        if (glons(i) == spval) cycle slice_loop
!
! Check for exact grid point:
        do ii=1,nlon
          if (glons(i) == gcmlon(ii)) cycle slice_loop
        enddo
!
! Check for zonal means:
        if (trim(float_to_str(glons(i))) == 'zm'.or.
     +      trim(float_to_str(glons(i))) == 'ZM') then
          glons(i) = zmflag

          write(6,"('chlons: i=',i2,' set glons(i)=zmflag=',e12.4)")
     +      i,zmflag

          cycle slice_loop
        elseif (glons(i)==zmflag) then
          cycle slice_loop
        endif
!
! Check for local times:
        if (isslt(glons(i),slt)) then
          if (slt < 0. .or. slt > 24.) then
            write(6,"('>>> Local time ',f10.2,' out of range.')")
     +        slt
            glons(i) = spval
          endif
          cycle slice_loop
        endif
!
! Check for outside grid: 
! (allow 0-360 lons -- convlon will convert them to -180->+180)
        glons_tmp = glons(i)     ! switch to kind=4 if default
        glons_tmp = convlon(glons_tmp,180)

        write(6,"('chlons: i=',i3,' glons(i)=',f8.2,' glons_tmp=',
     |    f8.2)") i,glons(i),glons_tmp

	glons(i) = glons_tmp
        if (glons(i) < gcmlon(1)) then
          write(6,"('Note: selected longitude ',f9.2,' is ',
     +      ' < western most longitude (',f9.2,') -- will ',
     +      'reset it to ',f9.2)") glons(i),gcmlon(1),gcmlon(1)
          glons(i) = gcmlon(1)
          cycle slice_loop
        elseif (glons(i) > gcmlon(nlon)) then
          write(6,"('Note: selected longitude ',f9.2,' is ',
     +      ' > eastern most longitude (',f9.2,') -- will ',
     +      'reset it to ',f9.2)") glons(i),gcmlon(nlon),gcmlon(nlon)
          glons(i) = gcmlon(nlon)
          cycle slice_loop
!
! Set to nearest grid point:
        elseif (ichange > 0) then
	  glons_tmp = glons(i)    ! switch to kind=4 if default
          ii = ixfind(gcmlon,nlon,glons_tmp,dlon)
          write(6,"('Note: changing selected longitude ',f7.2,' to ',
     +      'nearest tgcm grid longitude of ',f7.2)")
     +      glons(i),gcmlon(ii)
          glons(i) = gcmlon(ii)
        endif
      enddo slice_loop
      return
      end subroutine chlons

!-------------------------------------------------------------------
      subroutine chlats(glats,nlats,ichange)
      use proc,only: gcmlat,nlat
!
! Validate selected latitudes glats -- if ichange > 0 then change 
!   lats if necessary to nearest tgcm grid latitude.
!
! Args:
      integer,intent(in) :: nlats,ichange
      real(kind=8),intent(inout) :: glats(nlats)
!
! Locals:
      integer :: j,jj
      real :: glats_tmp
!
! Externals:
      integer,external :: ixfind
      character(len=16),external :: float_to_str
!
      slice_loop: do j=1,nlats
        if (glats(j) == spval) cycle slice_loop
        do jj=1,nlat
          if (glats(j) == gcmlat(jj)) cycle slice_loop
        enddo
!
! Sometimes lats are allowed to be 'zm', e.g. xylocs(:,i)='zm' to
! indicate global means:
        if (trim(float_to_str(glats(j))) == 'zm'.or.
     +      trim(float_to_str(glats(j))) == 'ZM') then
          glats(j) = zmflag
          cycle slice_loop
        elseif (glats(j)==zmflag) then
          cycle slice_loop
        endif
!
! Check for out-of-range:
        if (glats(j) < gcmlat(1)) then
          write(6,"('Note: selected latitude ',f9.2,' is ',
     +      ' < southern most latitude (',f9.2,') -- will ',
     +      'reset it to ',f9.2)") glats(j),gcmlat(1),gcmlat(1)
          glats(j) = gcmlat(1)
          cycle slice_loop
        elseif (glats(j) > gcmlat(nlat)) then
          write(6,"('Note: selected latitude ',f9.2,' is ',
     +      ' > northern most latitude (',f9.2,') -- will ',
     +      'reset it to ',f9.2)") glats(j),gcmlat(nlat),gcmlat(nlat)
          glats(j) = gcmlat(nlat)
          cycle slice_loop
        elseif (ichange > 0) then
	  glats_tmp = glats(j)   ! switch to kind=4 if default
          jj = ixfind(gcmlat,nlat,glats_tmp,dlat)
          write(6,"('Note: changing selected latitude ',f6.2,' to ',
     +      'nearest tgcm grid latitude of ',f6.2)")
     +      glats(j),gcmlat(jj)
          glats(j) = gcmlat(jj)
        endif
      enddo slice_loop
      return
      end subroutine chlats
!-------------------------------------------------------------------
      logical function findloc(glocs,nlocs,glat,glon)
!
! Return true if (glat,glon) is in glocs(2,nlocs), false otherwise:
!
      integer,intent(in) :: nlocs
      real(kind=8),intent(in) :: glocs(2,nlocs),glat,glon
      integer :: i
      findloc = .false.
      do i=1,nlocs
        if (glocs(1,i)==glat.and.glocs(2,i)==glon) findloc = .true. 
      enddo
      end function findloc
!-------------------------------------------------------------------
      subroutine verify_grid_levels(grid_levels) 
      use fields
      character(len=80),intent(inout) :: grid_levels(2,mxfproc)
!
! grid_levels(1,1)=='INTERFACES'-> plot all fields on interfaces (default)
! grid_levels(1,1)=='MIDPOINTS' -> plot all fields on midpoints
! grid_levels(1,1)=='HISTORY'   -> plot fields as on the history
! grid_levels(1,1)=='HIST+TUV'  -> as on the history, except T,U,V on interfaces
!
! Otherwise grid_levels(1,i) is field name and grid_levels(2,i)
!  is either INTERFACES or MIDPOINTS.
!
! If grid_levels is all blank, then set default 'HIST+TUV'
!
! Sub setflev (getflds.F) transfers grid_levels to f(ix)%zptype_req, 
! which is later compared with f(ix)%zptype to determine if conversion 
! is necessary. Conversions are done in proclat (T,U,V), and sub 
! set_zptype (getflds.F).
!
      integer :: i
      logical :: blank
      blank = .true.
      do i=1,mxfproc
        if (len_trim(grid_levels(1,i))>0.or.  
     |      len_trim(grid_levels(2,i))>0) blank=.false.
      enddo
      if (blank) then
        write(6,"('verify_grid_levels: set default grid_levels ',
     |    'to HIST+TUV (as on the history, except T,U,V on ',
     |    'interfaces)')")
        grid_levels(1,1) = "HIST+TUV"   ! as on history, except T,U,V to i-faces
      endif
      end subroutine verify_grid_levels
!-------------------------------------------------------------------
      end module input
