c
c------------------------------------------------------------------
c
      subroutine pltut(gcmut,nut,iut0,iut1)
c
c Plot interpolated gcm data with sat ut on x-axis, linear height
c   on y-axis:
c
      include 'tracsat.h'      
      common/cpmpcom/ pltl,pltr,npoints,pltsatut(mxpts)
      dimension gcmut(nut,nhtscale,nfplt),vp(4),
     +  plot(nut,nhtscale),sltvp(4)
c     data vp /.17,.87,.61,.96/, chsize/.040/
      data vp /.17,.87,.61,.94/, chsize/.040/
      character*32 labx
      character*56 histlab
      character*80 lab80
      common /agx/ xutlab,nlab,iget
      data dxlab/0.5/
      save vp
c
c Satellite ut for orbit is on x-axis:
c (it may need to altered to increase across a day boundary)
c
      npoints = iut1-iut0+1
      inc = 0
      do i=1,npoints-1
        pltsatut(i) = satut(iut0+i-1) ! for cpmpxy
        if (satut(iut0+i-1).gt.satut(iut0+i)) inc = i+1
      enddo
      pltsatut(npoints) = satut(iut0+npoints-1) ! for cpmpxy
      if (inc.gt.0) then
c       write(6,"('pltut: need to make satut increasing. npoints=',
c    +    i3,' Orig satut=',
c    +    /(9f7.2))") npoints,(satut(i),i=iut0,iut0+npoints-1)
        do i=inc,npoints
          satut(iut0+i-1) = satut(iut0+i-1)+24.     
          pltsatut(i) = pltsatut(i)+24.
        enddo  
c       write(6,"('Increasing satut=',/(9f7.2))")
c    +    (satut(i),i=iut0,iut0+npoints-1)
      endif
      npts = npoints
      if (npts.gt.mxpts.or.npts.le.0) then
        write(6,"('>>> pltut: bad npts=',i4,' iut0,1=',2i4)")
     +    npts,iut0,iut1
        return
      endif
      if (iut0.le.0.or.iut1.le.0.or.iut1.gt.nut) then
        write(6,"('>>> pltut: bad iut0 or iut1: iut0,1=',2i4)")
     +    iut0,iut1
        return
      endif
      write(6,"(' ')")
      write(6,"('pltut: npts=',i5,' iut0,1=',2i5)") npts,iut0,iut1
      write(6,"('       satut=',f6.3,' to ',f6.3)") satut(iut0),
     +  satut(iut1)
      write(6,"('       satlat,lon=',f6.2,',',f7.2,' to ',f6.2,',',
     +  f7.2)") satlat(iut0),satlon(iut0),satlat(iut1),satlon(iut1)
      xc = satut(iut0)
      xd = satut(iut1)
      pltl = xc ! for cpmpxy
      pltr = xd ! for cpmpxy
c
c Set up conpack for plotting with cylindrical equidistant maps:
      call cpsetr('XC1',xc)
      call cpsetr('XCM',xd)
      call cpsetr('YC1',htscale(1))
      call cpsetr('YCN',htscale(nhtscale))
      call cpseti('SET',0)
      call cpseti('MAP',0)
      xmid = 0.5*(vp(1)+vp(2))
      call cpsetr('ILX',xmid)
      call cpsetr('ILY',-.2)
      call cpseti('ILP',0)
      call cpseti('RWC',300)
      call clearstr(histlab)
      if (inttime.gt.0) then
        write(histlab,"('HISTORY ',a,' (with time interpolation)')") 
     +    histvol(ivol)(1:lenstr(histvol(ivol)))
      else
        write(histlab,"('HISTORY = ',a,' (no time interpolation)')") 
     +    histvol(ivol)(1:lenstr(histvol(ivol)))
      endif
      nf = 0
c
c Field loop (don't need to plot height):
      do ip=1,nftot
        if (ifplt(ip).le.0) goto 100
        nf = nf+1
        if (ip.eq.ixz) goto 100
        call cpseti('MAP',3) ! for irregular ut on x-axis (see cpmpxy)
        call set(vp(1),vp(2),vp(3),vp(4),xc,xd,htscale(1),
     +    htscale(nhtscale),1)
        do i=1,npts
          plot(i,:) = gcmut(iut0+i-1,:,nf)
        enddo
c
c Take log10 if necessary:
        if (logint(ip).gt.0) 
     +    call log10f(plot,nut*nhtscale,1.e-20,spval)
c
c Draw contours:
        call contour(plot,nut,npts,nhtscale,0,
     +    cint(ip),cmin(ip),cmax(ip))
c
c Define top label:
        call cpgetr('ZMN',zmn)
        call cpgetr('ZMX',zmx)
        call cpgetr('CIU',ciu)
        lenflab = lenstr(flab(ip))
        call clearstr(lab80)
        if (logint(ip).le.0) then
          write(lab80,"(a,' ',a,' (MIN ',e10.3,', MAX ',e10.3,
     +      ', INTERVAL ',f9.2,')')") model(1:lenstr(model)),
     +      flab(ip)(1:lenflab),zmn,zmx,ciu
        else
          write(lab80,"(a,' LOG10(',a,') (',e9.3,' to ',e9.3,' by ',
     +      e9.3,')')") model(1:lenstr(model)),
     +      flab(ip)(1:lenflab),zmn,zmx,ciu
        endif
c
c Draw labels:
        call wrlab(lab80(1:lenstr(lab80)),xmid,vp(4)+.05,.012)
        call wrlab(histlab(1:lenstr(histlab)),xmid,vp(4)+.02,.011)
        call clearstr(labx)
        write(labx,"('UT ',f5.2,'-',f5.2,' (DAY ',i5,')')") 
     +    satut(iut0),satut(iut0+npoints-1),iydgcm
        iget = 1
        call labrect(satut(iut0),npts,htscale,nhtscale,
     +    labx(1:lenstr(labx)),"HEIGHT (KM)",chsize)
        iget = 0
c
c Local time plot in center:
c
        sltvp(1) = 0.17
        sltvp(2) = 0.87
        sltvp(3) = 0.45 
        sltvp(4) = 0.55
        call pltslt(satut(iut0),satlon(iut0),npts,sltvp,chsize) 
c
c Show sat orbit in line plot over CE map at bottom:
c     subroutine maptrac(satut,satlat,satlon,npts,chsize,chsat)
c
        rvl = 0.17
        rvr = 0.87
        rvb = 0.06
        rvt = 0.41
        call set(rvl,rvr,rvb,rvt,0.,1.,0.,1.,1)
        call maptrac(satut(iut0),satlat(iut0),satlon(iut0),iydsat,
     +    npts,chsize,chsat)
        call wrlab(chsat(1:lenstr(chsat))//' ORBIT TRACK',
     +    xmid,rvb+.06,.011)
        call tail(satfile,labx)
        call wrlab('('//labx(1:lenstr(labx))//')',xmid,rvb+.03,.011)
c       call box(1)
        call frame
        write(6,"('Contoured field ',a,' (min,max,int=',3e12.4,
     +    ')')") flab(ip)(1:lenflab),zmn,zmx,ciu
 100    continue
      enddo
      return
      end