
      subroutine getinp
      include 'tigcmproc.h'
      dimension htsc(3)
c
      call lexcon(histvol ,8HHISTVOL ,nhvols)
      call lexcon(mtimes  ,8HDDHHMM  ,ntimes)
      call lexcon(iden    ,8HIDEN    ,niden)
      call lexcon(modelhts,8HMODELHTS,nmodhts)
      call lexcon(icolor  ,8HCOLOR   ,ncolor)
      call lexcon(cenglb  ,8HCENGLB  ,nglb)
      call lexcon(cenmolw ,8HCENMOLW ,nmolw)
      call lexcon(censatv ,8HCENSATV ,nsatv)
      call lexcon(perimlat,8HPOLLAT  ,npollat)
      call lexcon(icont   ,8HICONT   ,ncont)
      call lexcon(ituv    ,8HTUV     ,ntuv)
      call lexcon(ipuv    ,8HPUV     ,npuv)
      call lexcon(uvmax   ,8HUVMAX   ,nuvmax)
      call lexcon(log10map,8HLOG10   ,nlog10)
      call lexcon(spls    ,8HZP      ,npls)
      call lexcon(shts    ,8HHTS     ,nhts)
      call lexcon(slon    ,8HLON     ,nlon)
      call lexcon(sslt    ,8HSLT     ,nslt)
      call lexcon(slat    ,8HLAT     ,nlat)
      call lexcon(zprange ,8HZP_RANGE,nzprange)
      call lexcon(htsc    ,8HHT_SCALE,nhtscale)
      call lexcon(gloc    ,8HLOC     ,nloc)
      call lexcon(ifplt   ,8HFIELDS  ,npltact)
      call lexcon(ionvel  ,8HIONVEL  ,nionvel)
      call lexcon(iwrgrid ,8HWRGRID  ,nwrgrid)
      call lexcon(iwrascii,8HWRASCII ,nwrascii)
      call lexcon(mssascii,8HMSSASCII,nmsascii)
      call lexcon(dirascii,8HDIRASCII,ndirdat)
      call lexcon(dircgm  ,8HDIRCGM  ,ndircgm)
      call lexcon(cmin    ,8HCMIN    ,ncmin)
      call lexcon(cmax    ,8HCMAX    ,ncmax)
      call lexcon(cint    ,8HCINT    ,ncint)
      call readlx(5,ierr)
      if (ierr.ne.1) then
        write(6,"(' ')")
        write(6,"('>>> getinp readlx ierr=',i10)") ierr
        write(6,"(' ')")
        stop 'readlx'
      endif
c
c Iden = 0,1,2 for mass mixing ratios, number densities, or mixing
c ratios of number densities
c
      if (niden.lt.0.or.iden.gt.2) then
        write(6,"(' >>> getinp: no iden -- will default to cm-3')")
        iden = 1
      endif
c
c modelhts = 1 -> use model heights rather than newhts (see gettimes.a)
c
      if (nmodhts.le.0) modelhts = 0
c
c Color option:
c
      if (ncolor.le.0) icolor = 0
c
c Histvol(s)
c
      if (mod(nhvols,3).ne.0) then
        write(6,"('>>> getinp: history volumes must be 24-chars each',
     +    /' nhvols=',i3)") nhvols
        stop 'nhvols'
      endif
      nhvols = nhvols/3
      if (nhvols.gt.mxvols) then
        write(6,"(' >>> getinp: nhvols.gt.mxvols: nhvols=',I3,
     +    ' mxvols=',i3)") nhvols,mxvols
        stop 'nhvols'
      endif
c        
c Number of desired times:
c        
      if (ntimes.le.0) then
        write(6,"('>>> getinp: need model times (DDHHMM)')")
        stop 'mtimes'
      endif
      if (mod(ntimes,3).ne.0) then
        write(6,"('>> getinp: model times should come in triplets ',
     +    '(dd,hh,mm): ntimes=',i3)") ntimes
        stop 'ntimes'
      endif
      ntms = ntimes/3
      if (ntms.gt.mxtms) then
        write(6,"('>>> getinp: too many model times (DDHHMM);',
     +    ' ntms=',i3,' mxtms=',i3)") ntms,mxtms
        stop 'mxtms'
      endif
c
c Ut loop:
c (mh,mm=99 means search only)
c
      do it=1,ntms
        if (mtimes(2,it).ne.99.and.
     +     (mtimes(2,it).lt.0.or.mtimes(2,it).gt.23)) then
          write(6,"('>>> getinp: bad model hour: ',i3,' it=',i3)")
     +      mtimes(2,it),it
          stop 'mhour'
        endif
        if (mtimes(3,it).ne.99.and.
     +     (mtimes(3,it).lt.0.or.mtimes(3,it).gt.59)) then
          write(6,"('>>> getinp: bad model minute: ',i3,' it=',i3)")
     +      mtimes(3,it),it
          stop 'mm'
        endif
        ut(it) = float(mtimes(2,it))+float(mtimes(3,it))/60.
      enddo
c
c Ion drifts flag:
      if (nionvel.le.0) ionvel = 1
c
c ifplt(nftot) = desired fields to plot
c ifget(nftot) = fields needed (ifplt plus others as needed below):
c Thus, alloc for fields and call to getgcmf use ifget, whereas plot
c   routines use ifplt.
c
      if (npltact.ne.nftot) then
        write(6,"('>>> getinp: need ',i3,' values for FIELDS',
     +    ' -- got ',i3)") nftot,npltact
        stop 'nftot'
      endif
      do i=1,nfhist
        ifget(i) = ifplt(i)
      enddo 
c
c Need potential and heights for ion drifts:
      if (ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0.or.
     +    ifplt(ixuivi).gt.0) then
        ifget(ixpot) = 1
        ifget(ixz) = 1
c
c Need several fields if ionvel=2 (adding effect of neutral atmos):
c
        if (ionvel.eq.2) then
          ifget(ixu) = 1
          ifget(ixv) = 1
          ifget(ixw) = 1
          ifget(ixo1) = 1
          ifget(ixo2) = 1
          ifget(ixn2) = 1
        endif
      endif
c
c Need ne and hts to get fof2 and/or hmf2:
      if (ifplt(ixfof2).gt.0.or.ifplt(ixhmf2).gt.0) then
        ifget(ixne) = 1
        ifget(ixz) = 1
      endif
c
c Need o2,o,n2 for total density (rho):
      if (ifplt(ixrho).gt.0) then
        ifget(ixo1) = 1
        ifget(ixo2) = 1
        ifget(ixn2) = 1
      endif
c
c Need o2,o, and/or n2 for the various ratios (o/o2, o/n2, o2/n2, o/(o2+n2)):
      if (ifplt(ixroo2).gt.0) then
        ifget(ixo1) = 1
        ifget(ixo2) = 1
      endif
      if (ifplt(ixron2).gt.0) then
        ifget(ixo1) = 1
        ifget(ixn2) = 1
      endif
      if (ifplt(ixro2n2).gt.0) then
        ifget(ixo2) = 1
        ifget(ixn2) = 1
      endif
      if (ifplt(ixroo2n2).gt.0) then
        ifget(ixo1) = 1
        ifget(ixo2) = 1
        ifget(ixn2) = 1
      endif
c
c------------------ begin mapped projection inputs ------------------
c
      if (nglb.le.0.and.nmolw.le.0.and.nsatv.le.0.and.npollat.le.0) 
     +  goto 100
c
c Continental outlines flag (default is to draw continental outlines):
      if (ncont.le.0) icont = 1
c
c Need heights for height interpolation on mapped projections:
      if (nhts.gt.0.and.ifget(ixz).le.0) ifget(ixz) = 1
c
c Plot un+vn vectors over tn if ituv > 0: 
      if (ntuv.le.0) ituv = 1
      if (ifplt(ixunvn).gt.0.or.(ituv.gt.0.and.ifplt(ixt).gt.0)) then
        ifget(ixu) = 1
        ifget(ixv) = 1
      endif
c
c Plot ui+vi vectors over poten if ipuv > 0 (will need ht to find ui,vi):
      if (npuv.le.0) ipuv = 1
      if (ifplt(ixuivi).gt.0.or.(ifplt(ixpot).gt.0.and.ipuv.gt.0)) 
     +  ifget(ixz) = 1
c
c Max magnitude for vector scaling:
      if (nuvmax.le.0) uvmax = 0.
c
c Plot log10 on mapped projections if log10map > 0:
      if (nlog10.le.0) log10map = 0
c
c CE projection(s) (nglb may be = 0):
c (cenglb(2,i) = r12flag means use local noon longitude)
c
      if (nglb.gt.0)
     +  call checkcen(nglb,mxmap,cenglb,gcmlon,imx,r12flag,'CENGLB')
c
c MO projection(s) (nmolw may be = 0):
c (cenmolw(2,i) = r12flag means use local noon longitude)
c
      if (nmolw.gt.0)
     +  call checkcen(nmolw,mxmap,cenmolw,gcmlon,imx,r12flag,'CENMOLW')
c
c Sat view projection(s) (nsatv may be = 0):
c (censatv(2,i) = r12flag means use local noon longitude)
c
      if (nsatv.gt.0)
     +  call checkcen(nsatv,mxmap,censatv,gcmlon,imx,r12flag,'CENSATV')
c
c Perimeter latitudes (for polars)
c
      if (npollat.gt.mxmap) then
        write(6,"('>>> getinp: too many perim lats (POLLAT)=',i3
     +    ' will default to ',i3)") npollat,mxmap
        npollat = mxmap
      endif
      if (npollat.gt.0) then
        do i=1,npollat
          ixlat = ixfind(gcmlat,jmx,perimlat(i),dlat)
          if (ixlat.le.0) then
            write(6,"('>>> getinp: bad perimlat (POLLAT)=',f10.3,
     +        ' i=',i3)") perimlat(i),i
            stop 'perimlat'
          endif
          perimlat(i) = gcmlat(ixlat)
        enddo
      endif
c
c Selected pressures for mapped projections:
      if (npls.gt.kmx) then
        write(6,"('>>> getinp: too many spls: npls=',i3,
     +    ' will default to ',i3)") npls,kmx
        npls = kmx
      endif
c
c Selected heights for mapped projections:
      if (nhts.gt.kmx) then
        write(6,"('>>> getinp: too many hts: nhts=',i3,
     +    ' will default to ',i3)") nhts,kmx
        nhts = kmx
      endif
 100  continue
c
c-------------------- end mapped projection inputs ------------------
c
c-------------- begin vertical slice and loc plots inputs -----------
c
      if (nlon.le.0.and.nslt.le.0.and.nlat.le.0.and.nloc.le.0) goto 200
c
c Need heights for extra side axis and/or interp to linear ht scale:
      if (ifget(ixz).le.0) ifget(ixz) = 1
c
c Locations for line plots:
c
      if (nloc.gt.0) then
        if (mod(nloc,2).ne.0) then
          write(6,"('>>> getinp: need even number of locs (lat,lon)',
     +      ' nloc=',i3)") nloc
          stop 'nloc'
        else
          nloc = nloc/2
        endif
        if (nloc.gt.mxloc) then
          write(6,"('>>> getinp: too many locs: nloc=',i3,
     +      ' will default to ',i3)") nloc,mxloc
          stop 'nloc'
        endif
        write(6,"('getinp: nloc=',i2)") nloc
        do i=1,nloc
          if (gloc(1,i).ne.gmflag) then
            ilat = ixfind(gcmlat,jmx,gloc(1,i),dlat)
            if (ilat.le.0) then
              write(6,"('>>> getinp: bad latitude in LOC at i=',i2,
     +          ' read lat=',f10.3)") i,gloc(1,i)
              stop 'lon loc'
            endif
          endif
          if (gloc(2,i).ne.gmflag) then
            ilon = ixfind(gcmlon,imx,gloc(2,i),dlon)
            if (ilon.le.0) then
              write(6,"('>>> getinp: bad longitude in LOC at i=',i2,
     +          ' read lon=',f10.3)") i,gloc(2,i)
              stop 'lat loc'
            endif
          endif
          if (gloc(1,i).eq.gmflag.and.gloc(2,i).ne.gmflag) then
            write(6,"('>>> getinp: bad LOC ',i2,' read lat,lon=',
     +        2f10.3,/' (use lat,999. for zonal means, or 999.,999.',
     +        ' for global means')") i,(gloc(ii,i),ii=1,2)
            stop 'loc'
          endif
          if (gloc(1,i).ne.gmflag.and.gloc(2,i).ne.gmflag) then
            write(6,"('LOC ',i2,' = ',2f9.2)") i,(gloc(ii,i),ii=1,2)
          elseif (gloc(1,i).ne.gmflag) then
            write(6,"('LOC ',i2,': lat = ',f9.2,' (ZONAL MEANS)')") 
     +        i,gloc(1,i)
          else
            write(6,"('LOC ',i2,' (GLOBAL MEANS)')") i
          endif 
        enddo
      endif
c
c Selected longitudes (for lon slices made in pltlon):
      if (nlon.gt.imx) then
        write(6,"('>>> getinp: too many lons: nlon=',i3,
     +    ' will default to ',i3)") nlon,imx
        nlon = imx
      endif
c
c Selected local times (for lon slices made in pltlon):
      if (nslt.gt.imx) then
        write(6,"('>>> getinp: too many slt: nslt=',i3,
     +    ' will default to ',i3)") nslt,imx
        nslt = imx
      endif
      if (nslt.gt.0) then
        do i=1,nslt
          if (sslt(i).lt.0..or.sslt(i).gt.23.99999) then
            write(6,"('>>> getinp: bad SLT=',f10.2)") sslt(i)
            stop 'SLT'
          endif
        enddo
      endif
c
c Selected latitudes (for lat slices made in pltlat):
      if (nlat.gt.jmx) then
        write(6,"('>>> getinp: too many lats: nlat=',i3,
     +    ' will default to ',i3)") nlat,jmx
        nlat = jmx
      endif
c
c Range of zp for y-axis for vertical slices:
      if (nzprange.gt.0) then
        if (zprange(1).lt.gcmzp(1).or.zprange(2).gt.gcmzp(kmx).or.
     +      zprange(1).ge.zprange(2)) then
          write(6,"('>>> BAD ZP_RANGE = ',2f10.2,' gcmzp=',/(6e12.4))")
     +      zprange,gcmzp
          stop 'zprange'
        endif
        izprange(1) = ixfind(gcmzp,kmx,zprange(1),dzp)
        if (izprange(1).lt.0) then
          write(6,"('>>> error from ixfind for zprange(1): zprange(1)=',
     +      f10.2,' kmx=',i3,' gcmzp=',/(6e12.4))") zprange(1),kmx,gcmzp
          stop 'izprange'
        endif
        izprange(2) = ixfind(gcmzp,kmx,zprange(2),dzp)
        if (izprange(2).lt.0) then
          write(6,"('>>> error from ixfind for zprange(2): zprange(2)=',
     +      f10.2,' kmx=',i3,' gcmzp=',/(6e12.4))") zprange(2),kmx,gcmzp
          stop 'izprange'
        endif
        nzprange = izprange(2)-izprange(1)+1
        write(6,"('getinp: zprange=',2f8.2,' izprange=',2i3,
     +    ' nzprange=',i3)") zprange,izprange,nzprange
      else 
        izprange(1) = 0
        izprange(2) = 0
      endif
c
c Linear height scale (for pltlat,pltlon):
c User specifies bottom, top, and delta in that order:
c
      if (nhtscale.gt.0.and.nhtscale.ne.3) then
        write(6,"('>>> getinp: HTSCALE should be = htbot,httop,delht ',
     +    'e.g., 50.,500.,25.')")
        stop 'htscale'
      endif
      if (nhtscale.gt.0) then
        if (htsc(1).ge.htsc(2)) then
          write(6,"('>>> getinp: HTSCALE: top must be > bot: bot=',
     +      f10.2,' top=',f10.2)") htsc(1),htsc(2)
          stop 'htscale'
        endif
        if (htsc(3).le.0.) then
          write(6,"('>>> getinp: HTSCALE: bad delta ht=',f10.2)")htsc(3)
          stop 'htscale'
        endif
        htbot = htsc(1)
        httop = htsc(2)
        delht = htsc(3)
      else
        delht = 0.
      endif
c
c Set up linear height scale, if needed:
c
      nhtscale = 0
      if (delht.gt.0.) then
        nhtscale = ifix((httop-htbot) / delht + 1.0000001)
        call alloc(phtscale,nhtscale)
        do i=1,nhtscale
          htscale(i) = htbot + (i-1)*delht
        enddo
        write(6,"('getinp: nhtscale=',i2,' linear height scale = ',
     +    /(6e10.3))") nhtscale,(htscale(i),i=1,nhtscale)
      endif
 200  continue
c
c---------------- end vertical slice and loc plots inputs -----------
c
c Fields to get (nfget = number of fields to get (sum of ifget > 0)
c (max nfget = nfhist)
c
      nfget = 0
      do ip=1,nfhist
        if (ifget(ip).gt.0) nfget = nfget+1
      enddo 
      ixfld(:) = 0
      ii = 0
      do i=1,nfhist
        if (ifget(i).gt.0) then
          ii = ii+1
          ixfld(i) = ii
        endif
      enddo
c     write(6,"('getinp: ifget=',/(20i3))") ifget
c     write(6,"('getinp: ixfld=',/(20i3))") ixfld
c     write(6,"('getinp: ifplt=',/(20i3))") ifplt
c
c Contour intervals:
c
      if (ncint.le.0) then
        do ip=1,nftot
          cint(ip) = 0.
        enddo
      elseif (ncint.ne.nftot) then
        write(6,"('>>> Need ',i3,' contour intervals (CINT),'
     +    ' or none: got ',i3,' -- stop ')")
     +    nftot,ncint
        stop 'CINT'
      endif
c
c Contour minimums:
c
      if (ncmin.le.0) then
        do ip=1,nftot
          cmin(ip) = 1.
        enddo
      elseif (ncmin.ne.nftot) then
        write(6,"('>>> Need ',i3,' contour minimums (CMIN),'
     +    ' or none: got ',i3,' -- stop ')")
     +    nftot,ncmin
        stop 'CMIN'
      endif
c
c Contour maximums:
c
      if (ncmax.le.0) then
        do ip=1,nftot
          cmax(ip) = 0.
        enddo
      elseif (ncmax.ne.nftot) then
        write(6,"('>>> Need ',i3,' contour maximums (CMAX),'
     +    ' or none: got ',i3,' -- stop ')")
     +    nftot,ncmax
        stop 'CMAX'
      endif
c
c iwrgrid = 1 -> write ascii file of entire grid for each field
c                and each history
c
      if (nwrgrid.le.0) iwrgrid = 0
c
c iwrascii = 0 -> do not write ascii data file
c iwrascii = 1 -> write ascii data file per frame
c iwrascii = 2 -> write single ascii data file (append each frame)
c dirascii = remote dir for rcp of ascii files 
c
      if (nwrascii.le.0) iwrascii = 0
      if (ndirdat.gt.0) then
        if (ndirdat.ne.7) then
          write(6,"('>>> DIRASCII must be 56 characters: ',
     +      'ndirdat=',i3)") ndirdat
          stop 'DIRASCII'
        endif
      else
        call clearstr(dirascii)
      endif
c
c mssascii = mass store dir to put ascii files:
c (currently only used for WRGRID, i.e., global grid ascii files)
c
      if (nmsascii.gt.0) then
        if (nmsascii.ne.7) then
          write(6,"('>>> MSSASCII must be 56 characters: ',
     +      'nmsascii=',i3)") nmsascii
          stop 'MSSASCII'
        endif
      else
        call clearstr(nmsascii)
      endif
c
c dircgm = remote dir for rcp of gmeta
c
      if (ndircgm.gt.0) then
        if (ndircgm.ne.7) then
          write(6,"('>>> DIRCGM must be 56 characters: ',
     +      'ndircgm=',i3)") ndircgm
          stop 'DIRCGM'
        endif
      else
        call clearstr(dircgm)
      endif
c
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine checkcen(ncen,mxmap,cen,gcmlon,imx,glonflag,lab)
      dimension cen(2,mxmap), gcmlon(imx)
      character*(*) lab
c
      if (mod(ncen,2).ne.0) then
        write(6,"('>>> checkcen ',a,': projection center must be in ',
     +    'pairs (lat,lon) of center of projection: ncen=',i3)") 
     +  lab,ncen
        stop 'ncen'
      endif
      ncen = ncen/2
      if (ncen.gt.mxmap) then
        write(6,"('>>> checkcen warning ',a,': too many proj=',i3,
     +    ' will default to ',i3)") lab,ncen,mxmap
        ncen = mxmap
      endif
      do i=1,ncen
        if (cen(1,i).lt.-90..or.cen(1,i).gt.90.) then
          write(6,"('>>> checkcen ',a,': bad lat in cen=',f10.2,' i=',
     +      i2)") lab,cen(1,i),i
          stop 'cen'
        endif
        if (cen(2,i).ne.glonflag.and.(cen(2,i).lt.gcmlon(1).or.
     +      cen(2,i).gt.gcmlon(imx))) then
          write(6,"('>>> checkcen ',a,': bad lon in cen=',f10.2,
     +      ' i=',i2)") lab,cen(2,i),i
          stop 'cen'
        endif
      enddo
c
      return
      end
