
      program makebob
c
c Make unix byte file for use with Bob ("Brick Of Bytes") on Picard (SGI) 
c from timesgcm histories:
c
      include 'tigcm.h'
      dimension f(imx,kmx,jmx),f0(imx,kmx,jmx),f1(imx,kmx,jmx)
      logical istiegcm
      dimension iasf(13)
      data iasf/13*1/
      data luhist/8/, iframe/0/
c
c nfhist = 17: tn,un,vn,o2,o1,n4s,no,o+,n2d,ti,te,ne,o2+,w,z,pot,n2
      data fieldlab
     + /'TN      ','UN      ','VN      ','O2      ','O1      ',  
     +  'N4S     ','NO      ','O+      ','N2D     ','TI      ',
     +  'TE      ','NE      ','O2+     ','WN      ','Z       ',
     +  'POTEN   ','N2      '/
c
      call setgrid(gcmlat,glat1,dlat,jmx, gcmlon,glon1,dlon,imx,
     +             gcmzp,zp1,dzp,kmx)
      call getinp
      if (icontour.gt.0) then
        call opngks
        call gsclip(0)
        call gsasf(iasf)
        call gsfais(1)
        call cpsetup(cpspval)
        call cpseti('LLB',3)
        call cpseti('LLB',1)
        call cpseti('LLB',3)
        call cpseti('LBC',1)
        call setag
        call agseti('AXIS/LEFT/CONTROL.',1)
        call agseti('AXIS/BOTTOM/CONTROL.',1)
        ixzp = ixfind(gcmzp,kmx,szp,dzp)
        if (ixzp.le.0) then
          write(6,"('>>> bad szp=',f10.3)") szp
          stop 'szp'
        endif
      endif
      istart = 1
      ircp = 0
      totmin = 1.e36
      totmax = -1.e36
c
c Time loop: get ti(e)gcm data for current ut:
c
      it = 1
 200  continue
      if (ninterp.le.0.or.it.eq.1) then
        call getgcmf(histvol,nhvols,istart,mtimes(1,it),f,1,
     +    mkfields,imx,kmx,jmx,nfhist,iden,luhist,istart,ivol,
     +    istiegcm,ier)
        if (ier.eq.1) then
          write(6,"('tigcmbob: getgcm could not find history ',3i3,
     +      ' it=',i3)") (mtimes(i,it),i=1,3),it
          it = it+1
          goto 100
        elseif (ier.ne.0) then
          write(6,"('Error ',i3,' from getgcm: it=',i3,' mt=',3i3)")
     +      it,(mtimes(i,it),i=1,3)
          stop 'getgcm'
        endif
        ut(it) = mtimes(2,it)+mtimes(3,it)/60.
      endif
      istart = 0
      if (it.eq.ntms) ircp = 1
      iframe = iframe+1
      if (mkbob.gt.0) then
        write(6,"(' ')")
        write(6,"('Calling mkfile for ut=',f6.2,' it=',i3,' frame ',
     +    i4)") ut(it),it,iframe
        call mkfile(f,ircp,cmin,cmax)
      endif
      if (icontour.gt.0) call pltf(f,imx,kmx,jmx,ixzp,iframe,ut(it),
     +  cint,cmin,cmax,ilog(ixf),ncenslt,censlt,gcmlon,dlon)
      if (info.gt.0) call rptinfo(f,imx,kmx,jmx,totmin,totmax,iframe,
     +  ilog(ixf),ut(it))
      it = it+1
      if (it.gt.ntms) goto 300
      if (ninterp.gt.0) then
        f0(:,:,:) = f(:,:,:)
        call getgcmf(histvol,nhvols,istart,mtimes(1,it),f1,1,
     +    mkfields,imx,kmx,jmx,nfhist,iden,luhist,istart,ivol,
     +    istiegcm,ier)
        if (ier.eq.1) then
          write(6,"('tigcmbob: getgcm could not find history ',3i3,
     +      ' it=',i3)") (mtimes(i,it),i=1,3),it
          it = it+1
          goto 100
        elseif (ier.ne.0) then
          write(6,"('Error ',i3,' from getgcm: it=',i3,' mt=',3i3)")
     +      it,(mtimes(i,it),i=1,3)
          stop 'getgcm'
        endif
        ut(it) = mtimes(2,it)+mtimes(3,it)/60.
c
c Time interpolation(s) between f0 (at it-1) and f1 (at it):
c
        do in=1,ninterp
          if (ut(it).gt.ut(it-1)) then
            utx = ut(it-1)+float(in)/float(ninterp+1)*(ut(it)-ut(it-1))
            call tinterp(f0,f1,f,imx,kmx,jmx,ut(it-1),ut(it),utx)
          else
            utx = ut(it-1)+float(in)/float(ninterp+1)*
     +            (ut(it)+24.-ut(it-1))
            call tinterp(f0,f1,f,imx,kmx,jmx,ut(it-1),ut(it)+24.,utx)
          endif
          iframe = iframe+1
          if (mkbob.gt.0) then
            write(6,"(' ')")
            write(6,"('Calling mkfile for ut=',f6.2,' in=',i3,
     +        ' iframe ',i4)") utx,in,iframe
            call mkfile(f,ircp,cmin,cmax)
          endif
          if (icontour.gt.0) call pltf(f,imx,kmx,jmx,ixzp,iframe,utx,
     +      cint,cmin,cmax,ilog(ixf),ncenslt,censlt,gcmlon,dlon)
          if (info.gt.0) call rptinfo(f,imx,kmx,jmx,totmin,totmax,
     +      iframe,ilog(ixf),ut(it))
        enddo
        f(:,:,:) = f1(:,:,:) 
      endif
 100  continue
      if (it.le.ntms) goto 200
c
 300  continue
      if (icontour.gt.0) call clsgks
      stop 'success!'
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mkfile(f,ircp,rmin,rmax)
      include 'tigcm.h'
      dimension f(imx,kmx,jmx),glb(imx,kmx,jmx),ix(imx)
      character*16 flnm
      character*56 rdir
      data rdir
     + /'sting.hao.ucar.edu:/d/foster/tigcmbob                   '/
c
c Mkbyte writes byte file acceptable by Bob (on picard)
c If censlt was set from lexical reads, it will be local time at 
c   center of longitude axis, otherwise x-axis will be -180 to +180.
c
      if (ncenslt.le.0) then
        if (ilog(ixf).gt.0) call log10f(glb,imx*kmx*jmx,1.e-20,spval)
        call clearstr(flnm)
        write(flnm,"(a,'.dat')")fieldlab(ixf)(1:lenstr(fieldlab(ixf)))
        len = lenstr(flnm)
        ier = mkbyte(f,imx,kmx,jmx,flnm(1:len),rmin,rmax)
      else
        slt0 = censlt-12.
        if (slt0.lt.0.) slt0 = slt0-12.
        ix0 = ixslt(slt0,ut(it),xlon,gcmlon,imx,dlon)
        iflag = 0
        do i=1,imx
          ix(i) = mod(ix0+i-1,imx)
          if (ix(i).eq.0) ix(i) = imx
          if (i.gt.1.and.ix(i-1).gt.ix(i)) iflag = 1
          if (iflag.gt.0) ix(i) = ix(i)+1
          glb(i,:,:) = f(ix(i),:,:)
        enddo
        if (ilog(ixf).gt.0) call log10f(glb,imx*kmx*jmx,1.e-20,spval)
        call clearstr(flnm)
        write(flnm,"(a,'.dat')")fieldlab(ixf)(1:lenstr(fieldlab(ixf)))
        len = lenstr(flnm)
        ier = mkbyte(glb,imx,kmx,jmx,flnm(1:len),rmin,rmax)
      endif
      if (ircp.gt.0) call rcpfile(0,flnm,rdir)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine tinterp(f0,f1,f,imx,kmx,jmx,ut0,ut1,utx)
      dimension f0(imx,kmx,jmx),f1(imx,kmx,jmx),f(imx,kmx,jmx)
c
      frac = (utx-ut0) / (ut1-ut0)
      f(:,:,:) = frac*f1(:,:,:) + (1.-frac)*f0(:,:,:)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine rptinfo(f,imx,kmx,jmx,totmin,totmax,iframe,log,ut)
      dimension f(imx,kmx,jmx)
c
      if (log.gt.0) call log10f(f,imx*jmx*kmx,1.e-20,1.e36)
      call fminmax(f,imx*kmx*jmx,rmin,rmax,1.e36)
      if (rmin.lt.totmin) totmin = rmin
      if (rmax.gt.totmax) totmax = rmax
      write(6,"('info: frame ',i4,' ut=',f5.2,' rmin,max=',2e12.4,
     +  ' totmin,max=',2e12.4)") iframe,ut,rmin,rmax,totmin,totmax
      return
      end
