
      subroutine labutxy(mt,ntms,yy,ny,ylab,charsize,isltax,rlon)
c
c Label time x-axis; also label y-axis;
c Add local time axis below time axis if isltax > 0
c Assumes constant delta time along x-axis
c
      dimension mt(3,ntms),yy(ny)
      character*(*) ylab
      data chsize/.015/
      cmt2dhrs(md,mh,mm) = float(md)*24.+float(mh)+float(mm)/60.
      cmt2dday(md,mh,mm) = float(md)+float(mh)/24.+float(mm)/(24.*60.)
c
      dday0 = cmt2dday(mt(1,1),mt(2,1),mt(3,1))
      dday1 = cmt2dday(mt(1,ntms),mt(2,ntms),mt(3,ntms))
      call getset(vl,vr,vb,vt, wl,wr,wb,wt, log)
      if (ny.gt.0) then
        call set(vl,vr,vb,vt,dday0,dday1,yy(1),yy(ny),1)
        call labrect(dum,0,yy,ny,' ',ylab,chsize+.01)
      else
        call set(vl,vr,vb,vt,dday0,dday1,0.,1.,1)
      endif
      if (charsize.gt.0.) chsize = charsize
c
c Let autograph label y-axis:
c
      dhrs = cmt2dhrs(mt(1,ntms),mt(2,ntms),mt(3,ntms))-
     +       cmt2dhrs(mt(1,1),mt(2,1),mt(3,1))
c
c time = 0->0.25 hrs (0-15 min); label every 3 minutes:
      if (dhrs.le.0.25) then
        call labutax(mt,ntms,3.,3,3,chsize)
        if (isltax.gt.0) call labsltax(rlon,.05,3,chsize)
c
c time = .25->0.75 hrs (15-45 min); label every 5 minutes:
      elseif (dhrs.le.0.75) then
        call labutax(mt,ntms,5.,3,5,chsize)
        if (isltax.gt.0) call labsltax(rlon,.25,5,chsize)
c
c time = 0.75->1.9 hrs; label every 10 minutes:
      elseif (dhrs.le.1.9) then
        call labutax(mt,ntms,10.,3,5,chsize)
        if (isltax.gt.0) call labsltax(rlon,.25,5,chsize)
c
c time = 2.0->5.0 hrs; label every 30 minutes:
      elseif (dhrs.le.5.0) then
        call labutax(mt,ntms,30.,3,3,chsize)
        if (isltax.gt.0) call labsltax(rlon,.75,3,chsize)
c
c time = 5->12 hrs; label every hour:
      elseif (dhrs.le.12.0) then
        call labutax(mt,ntms,1.,2,4,chsize)
        if (isltax.gt.0) call labsltax(rlon,2.,4,chsize)
c
c time = 12->24 hrs (.5-1 days); label every 3 hours:
      elseif (dhrs.le.24.0) then
        call labutax(mt,ntms,3.,2,3,chsize)
        if (isltax.gt.0) call labsltax(rlon,4.,4,chsize)
c
c time = 24->48 hrs (1-2 days); label every 4 hours:
      elseif (dhrs.le.48.0) then
        call labutax(mt,ntms,4.,2,4,chsize)
        if (isltax.gt.0) call labsltax(rlon,6.,3,chsize)
c
c time = 48->72 hrs (2-3 days); label every 8 hours:
      elseif (dhrs.le.72.0) then
        call labutax(mt,ntms,8.,2,4,chsize)
        if (isltax.gt.0) call labsltax(rlon,12.,6,chsize)
c
c Otherwise label every day:
      else
        call labutax(mt,ntms,1.,1,4,chsize)
        if (isltax.gt.0) call labsltax(rlon,24.,4,chsize)
      endif
c
      call set(vl,vr,vb,vt, wl,wr,wb,wt, log)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine labutax(mt,ntms,del,itime,nmnr,charsize)
c
c Plot numeric labels and tics on x-axis, at intervals of del time
c If itime = 1,2,3 -> del is in days, hours, or minutes, and labels
c                     should be in days, hours, or minutes
c (Also plot appropriate text label)
c nmnr = number of minor spaces between each pair of labeled tics
c charsize = character size (has default if 0 on input)
c NOTE: on input, set call has been made, with wl,wr in decimal days
c

      dimension mt(3,ntms)
      data yoff/.03/, chsize/.02/, ticlen/.014/ 
      character*16 str
      character*40 mtimelab
      character*56 xlab
c
      cmt2dmin(md,mh,mm) = float(md)*24.*60.+float(mh)*60.+float(mm)
      cmt2dhrs(md,mh,mm) = float(md)*24.+float(mh)+float(mm)/60.
      cmt2dday(md,mh,mm) = float(md)+float(mh)/24.+float(mm)/(24.*60.)
c
      call getset(vl,vr,vb,vt, wl,wr,wb,wt, log)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      ypos = vb-yoff
      rminday = 24.*60.
      idel = int(del)
      if (charsize.gt.0.) chsize = charsize
      write(mtimelab,"(' (MTIMES ',i2.2,' ',i2.2,':',i2.2,
     +  ' TO ',i2.2,' ',i2.2,':',i2.2,') ')") (mt(i,1),i=1,3),
     +  (mt(i,ntms),i=1,3)
c
c Delta is in days:
c
      if (itime.eq.1) then
        call clearstr(xlab)
        write(xlab,"('UT (DAYS)',a)") mtimelab
        rlab = float(int(wl)/idel*idel)
        if (rlab.lt.wl) rlab = rlab+del
        xpos0 = vl + (rlab-wl)*(vr-vl)/(wr-wl)
        xpos1 = vl + (rlab+del-wl)*(vr-vl)/(wr-wl)
        dmnr = 1./nmnr * (xpos1-xpos0)
        do i=1,1000
          xpos = vl + (rlab-wl)*(vr-vl)/(wr-wl)
          if (xpos.lt.vl.or.xpos.gt.vr) goto 100
          do ii=1,nmnr-1
            xmnr = xpos-ii*dmnr
            if (xmnr.gt.vl) then
              call line(xmnr,vb,xmnr,vb+.5*ticlen)
              call line(xmnr,vt,xmnr,vt-.5*ticlen)
            endif
          enddo         
          call clearstr(str)
          call ddaytomt(rlab,md,mh,mm)
          write(str,"(i2)") md
          call rmblnk(str)
          len = lenstr(str)
c         write(6,"('i=',i2,' rlab=',f9.4,' str=',a,' xpos=',f9.4)")
c    +      i,rlab,str,xpos
          call plchmq(xpos,ypos,str(1:len),chsize,0.,0.)
          call line(xpos,vb,xpos,vb+ticlen)
          call line(xpos,vt,xpos,vt-ticlen)
          rlab = rlab+del
          xposp = xpos
        enddo
 100    continue
        do ii=1,nmnr-1
          xmnr = xpos-ii*dmnr
          if (xmnr.ge.vl.and.xmnr.le.vr) then
            call line(xmnr,vb,xmnr,vb+.5*ticlen)
            call line(xmnr,vt,xmnr,vt-.5*ticlen)
          endif
        enddo         
c
c Delta is in hours:
c
      elseif (itime.eq.2) then
        write(xlab,"('UT (HOURS)',a)") mtimelab
        wlhrs = cmt2dhrs(mt(1,1),mt(2,1),mt(3,1))
        wrhrs = cmt2dhrs(mt(1,ntms),mt(2,ntms),mt(3,ntms))
        rlab = float(int(wlhrs)/idel*idel)
        if (rlab.lt.wlhrs) rlab = rlab+del
        xpos0 = vl + (rlab-wlhrs)*(vr-vl)/(wrhrs-wlhrs)
        xpos1 = vl + (rlab+del-wlhrs)*(vr-vl)/(wrhrs-wlhrs)
        dmnr = 1./nmnr * (xpos1-xpos0)
        do i=1,1000
          xpos = vl + (rlab-wlhrs)*(vr-vl)/(wrhrs-wlhrs)
          if (xpos.lt.vl.or.xpos.gt.vr) goto 200
          do ii=1,nmnr-1
            xmnr = xpos-ii*dmnr
            if (xmnr.gt.vl) then
              call line(xmnr,vb,xmnr,vb+.5*ticlen)
              call line(xmnr,vt,xmnr,vt-.5*ticlen)
            endif
          enddo         
          call clearstr(str)
          call dhrstomt(rlab,md,mh,mm)
          write(str,"(i2)") mh
          call rmblnk(str)
          len = lenstr(str)
c         write(6,"('i=',i2,' rlab=',f9.4,' str=',a,' xpos=',f9.4)")
c    +      i,rlab,str,xpos
          call plchmq(xpos,ypos,str(1:len),chsize,0.,0.)
          call line(xpos,vb,xpos,vb+ticlen)
          call line(xpos,vt,xpos,vt-ticlen)
          rlab = rlab+del
          xposp = xpos
        enddo
 200    continue
        do ii=1,nmnr-1
          xmnr = xpos-ii*dmnr
          if (xmnr.ge.vl.and.xmnr.le.vr) then 
            call line(xmnr,vb,xmnr,vb+.5*ticlen)
            call line(xmnr,vt,xmnr,vt-.5*ticlen)
          endif
        enddo         
c
c Delta is in minutes
c
      elseif (itime.eq.3) then
        write(xlab,"('UT (MINUTES)',a)") mtimelab
        wlmin = cmt2dmin(mt(1,1),mt(2,1),mt(3,1))
        wrmin = cmt2dmin(mt(1,ntms),mt(2,ntms),mt(3,ntms))
        rlab = float(int(wlmin)/idel*idel)
        if (rlab.lt.wlmin) rlab = rlab+del
        xpos0 = vl + (rlab-wlmin)*(vr-vl)/(wrmin-wlmin)
        xpos1 = vl + (rlab+del-wlmin)*(vr-vl)/(wrmin-wlmin)
        dmnr = 1./nmnr * (xpos1-xpos0)
        do i=1,1000
          xpos = vl + (rlab-wlmin)*(vr-vl)/(wrmin-wlmin)
          if (xpos.lt.vl.or.xpos.gt.vr) goto 300
          do ii=1,nmnr-1
            xmnr = xpos-ii*dmnr
            if (xmnr.gt.vl) then
              call line(xmnr,vb,xmnr,vb+.5*ticlen)
              call line(xmnr,vt,xmnr,vt-.5*ticlen)
            endif
          enddo         
          call clearstr(str)
          call dmintomt(rlab,md,mh,mm)
          write(str,"(i2)") mm
          call rmblnk(str)
          len = lenstr(str)
c         write(6,"('i=',i2,' rlab=',f9.4,' str=',a,' xpos=',f9.4)")
c    +      i,rlab,str,xpos
          call plchmq(xpos,ypos,str(1:len),chsize,0.,0.)
          call line(xpos,vb,xpos,vb+ticlen)
          call line(xpos,vt,xpos,vt-ticlen)
          rlab = rlab+del
          xposp = xpos
        enddo
 300    continue
        do ii=1,nmnr-1
          xmnr = xpos-ii*dmnr
          if (xmnr.ge.vl.and.xmnr.le.vr) then
            call line(xmnr,vb,xmnr,vb+.5*ticlen)
            call line(xmnr,vt,xmnr,vt-.5*ticlen)
          endif
        enddo         
      endif
c
c Add time label:
c
      call plchmq(0.5*(vl+vr),vb-2.2*yoff,xlab(1:lenstr(xlab)),
     +  chsize,0.,0.)
      call set(vl,vr,vb,vt, wl,wr,wb,wt, log)
c
      return
      end 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine labsltax(rlon,del,nmnr,charsize)
c
c Draw slt axis (hrs) below existing ut axis 
c (set call has been done with wl,wr = decimal days)
c
      character*8 str
      character*32 xlab
      data chsize/.015/
      data yoff/.12/, ticlen/.014/, ynumoff/.03/
      cmt2dhrs(md,mh,mm) = float(md)*24.+float(mh)+float(mm)/60.

      call getset(vl,vr,vb,vt, wl,wr,wb,wt, log)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      yposax = vb-yoff
      if (charsize.gt.0) chsize = charsize
      idel = int(del)
      if (idel.le.0) idel = 1
      call clearstr(xlab)
      write(xlab,"('LOCAL TIME (HRS) (LON=',f7.2,')')") rlon
c
c Draw the axis:
c
      call line(vl,yposax,vr,yposax)
c
c Convert to decimal slt (including model day):
c
      call ddaytomt(wl,mdl,mhl,mml)
      wlhrs = cmt2dhrs(mdl,mhl,mml)
      wlslt = dut2dslt(wlhrs,rlon)
      call ddaytomt(wr,mdr,mhr,mmr)
      wrhrs = cmt2dhrs(mdr,mhr,mmr)
      wrslt = wlslt + (wrhrs-wlhrs)
      sltdhrs = float(int(wlslt))
      sltdhrs = rnd(sltdhrs,del)
      slt = dhrs2ut(sltdhrs)
      if (sltdhrs.lt.wlslt) then
        do i=1,100
          sltdhrs = sltdhrs+del
          if (sltdhrs.ge.wlslt) goto 50
        enddo
      endif
 50   continue
      ypos = vb-yoff-ynumoff
      xpos0 = vl + (sltdhrs-wlslt)*(vr-vl)/(wrslt-wlslt)
      xpos1 = vl + (sltdhrs+del-wlslt)*(vr-vl)/(wrslt-wlslt)
      dmnr = 1./nmnr * (xpos1-xpos0)
c
c Loop on labeled intervals:
c
      do i=1,20
        xpos = vl + (sltdhrs-wlslt)*(vr-vl)/(wrslt-wlslt)
        if (xpos.gt.vr) goto 100
        if (xpos.ge.vl) then
          do ii=1,nmnr-1
            xmnr = xpos-ii*dmnr
            if (xmnr.ge.vl.and.xmnr.le.vr) 
     +        call line(xmnr,yposax,xmnr,yposax+.5*ticlen)
          enddo         
          slt = dhrs2ut(sltdhrs)
          call clearstr(str)
          write(str,"(f5.2)") slt
          call rmblnk(str)
          len = lenstr(str)
c         write(6,"('i=',i2,' slt=',f9.4,' sltdhrs=',f12.5,' str=',a,
c    +      ' xpos=',f9.4)") i,slt,sltdhrs,str,xpos
          call plchmq(xpos,ypos,str(1:len),chsize,0.,0.)
          call line(xpos,yposax,xpos,yposax+ticlen)
        endif
        sltdhrs = sltdhrs+del
      enddo
 100  continue
      do ii=1,nmnr-1
        xmnr = xpos-ii*dmnr
        if (xmnr.ge.vl.and.xmnr.le.vr) 
     +    call line(xmnr,yposax,xmnr,yposax+.5*ticlen)
      enddo         
      call plchmq(0.5*(vl+vr),vb-yoff-2.2*ynumoff,xlab,chsize,0.,0.)
      call set(vl,vr,vb,vt, wl,wr,wb,wt, log)
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine rmblnk(str)
c
c Remove any blanks from str
c
      character*(*) str
      character*80 rstr
c
      mxlen = 80
      length = len(str)
      if (length.gt.80) then
        write(6,"('>>> rmblnk: str too big=',i5,' need <= ',i2)") mxlen
        return
      endif
      ipos = 0
      call clearstr(rstr)
      do i=1,length
        if (str(i:i).ne.' ') then
          ipos = ipos+1
          rstr(ipos:ipos) = str(i:i)
        endif
      enddo
      call clearstr(str)
      if (ipos.gt.0) then
        do i=1,ipos
          str(i:i) = rstr(i:i)
        enddo
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function rnd(val,rndval)
      ival = (val + (rndval-.0001)) / amax1(rndval,1.)
      rnd = ival * amax1(rndval,1.)
      return
      end 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine ddaytomt(decday,md,mh,mm)
c
c Return model time (md,mh,mm), given time in decimal days (decday):
c
      md = int(decday)
      mh = (decday-float(md))*24.
      mm = (decday*(24.*60.)-float(md)*(24.*60.)-float(mh)*60.)+.5
      if (mm.eq.60) then
        mm = 0
        mh = mh+1
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine dhrstomt(dechrs,md,mh,mm)
c
c Return model time (md,mh,mm), given time in decimal hours (dechrs):
c
      rminday = 24.*60.
      md = dechrs/24.
      mh = dechrs-float(md)*24.
      if (mh.eq.24) then
        mh = 0
        md = md+1
      endif
      mm = dechrs*60.-float(md)*rminday-float(mh)*60.+.5
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine dmintomt(decmin,md,mh,mm)
c
c Return model time (md,mh,mm), given time in decimal minutes (decmin):
c
      rminday = 24.*60.
      md = decmin/rminday+1.e-10
      mh = (decmin-float(md)*rminday)/60.
      mm = decmin-float(md)*rminday-float(mh)*60.
      if (mm.eq.60) mm = 0
      if (mh.eq.24) mh = 0
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function dhrs2ut(dechrs)
c
c Return ut, given time in dechrs (dechrs includes model day)
c
      call dhrstomt(dechrs,md,mh,mm)
      dhrs2ut = float(mh)+float(mm)/60.
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function dut2dslt(dechrs,rlon)
c
c Convert decimal ut (includes model day) to decimal slt (including 
c model day)
c
      call dhrstomt(dechrs,md,mh,mm)
      ut = float(mh)+float(mm)/60.
      slt = fslt(0.,ut,rlon,1)
      dut2dslt = float(md)*24.+slt 
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function fslt(slt,ut,glon,iflag)
c
c If iflag = 1 -> return slt, given ut and glon
c If iflag = 2 -> return ut, given slt and glon
c If iflag = 3 -> return glon, given slt and ut
c (slt and ut always in decimal hours, glon in -180 -> 180)
c
      if (iflag.eq.1) then
        fslt = ut + glon/15.
        if (fslt.lt.0.) fslt = fslt+24.
        if (fslt.ge.24.) fslt = fslt-24.
      elseif (iflag.eq.2) then
        fslt = slt - glon/15.  
        if (fslt.lt.0.) fslt = fslt+24.
        if (fslt.gt.24.) fslt = fslt-24.
      elseif (iflag.eq.3) then
        fslt = (slt-ut) * 15.
        if (fslt.ge.180.) fslt = fslt-360.
        if (fslt.lt.-180.) fslt = fslt+360.
      else
        write(6,"('>>> fslt: bad iflag=',i3)") iflag
        fslt = 0.
      endif
      return
      end
