c
c------------------------------------------------------------------
c Begin file /home/sting/foster/vid/vid.f
c------------------------------------------------------------------
c
      program vid
      include "vid.h"
      dimension glb(imx,jmx),iasf(13),plot(imx,jmx),plot0(imx,jmx),
     +  plot1(imx,jmx),mtime(3)
      data lu/8/, lukp/20/, luflash/90/, idflash/20/,
     +  lutitles/80/
      data ifbbeg/1/, ifbend/2/, ifbkp/3/
      data iasf/13*1/
      data iframe/0/
      data spv/1.e36/
c
c Fields array f(ntms,imx,jmx,nfhist)
      pointer (pf,f(1))
c
      call setgrid(gcmlat,glat1,dlat,jmx, gcmlon,glon1,dlon,imx,
     +             gcmzp,gzp1,dzp,kmx)
      call getinp
      write(6,"('vid: nvol=',i3,' histvol=')") nvol
      do i=1,nvol
        write(6,"(a)") histvol(i)(1:lenstr(histvol(i)))
      enddo
      write(6,"('vid: ntms=',i3,' mtimes=',/(5(i2,':',i2,':',i2,
     +  ',')))") ntms,((mtimes(i,it),i=1,3),it=1,ntms) 
c
c Read kp from disk file (file was made from ~/lexinp/lexplt 
c kp originally from geophysical indices database)
c
      if (kpbar.gt.0.or.kpgraph.gt.0) then
        ngpi = 0
        do i=1,mxgpi
          read(lukp,"(20x,i5,5x,f5.2,5x,e12.4)",end=101)
     +      igpiday(i),gpiut(i),gpikp(i)
          ngpi = ngpi+1
        enddo
        write(6,"('>>> vid: reached 200 reading kp from lukp=',
     +    i3)") lukp
 101    continue
        if (ngpi.le.1) then
          write(6,"('>>> vid: need kp file: ngpi=',i3)") ngpi
          write(6,"('Kp file may be in sting:/d/foster/vid/mar79kp.dat',
     +      /'or on mss as /FOSTER/vid/mar79kp.dat')")
          write(6,"('(assign to fort.',i2,')')") lukp
          stop 'kp'
        endif
      endif
c
c Allocate fields array:
      call hpalloc(pf,imx*jmx*nfhist,ier,1)
      if (ier.ne.0) then
        write(6,"('vid: hpalloc error = ',i6,' for f')") ier
        stop 'f alloc'
      endif
c
c Set up gks:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
c     call gslwsc(2.0)
      call gslwsc(1.0)
c
c Fields loop:
c
      nfplt = 0
      iip = 0
      do ip=1,nflds
        if (ipltgcm(ip).gt.0) nfplt = nfplt+1
      enddo
      do 100 ip=1,nflds
        if (ipltgcm(ip).le.0) goto 100
        iip = iip+1
c
c Add titles if desired:
c (this must be done before gopwk for freeze frames and kp -- see below)
c
c Title page(s) for beginning of video:
c
        if (ititles(1).ne.0.and.iip.eq.1) then
          if (ititles(1).lt.0) then
            itest = 1
            write(6,"('vid calling ftitle ',i2,
     +        ' times for beginning (test run)')") abs(ititles(1))
          else
            itest = 0
            write(6,"('vid calling ftitle ',i2,
     +        ' times for beginning (production run!)')") ititles(1)
          endif
          call slsetr('FIN',1.5)
          call slsetr('FOU',1.5)
          call slseti('LOG',lutitles)
          do i=1,abs(ititles(1))
            call ftitle(itest)
          enddo
        endif
c
c Title page(s) before each field segment:
c
        if (ititles(2).ne.0) then
          if (ititles(2).lt.0) then
            itest = 1
            write(6,"('vid calling ftitle ',i2,' times for iip=',i2,
     +        ' (test run)')") abs(ititles(2)),iip
          else
            itest = 0
            write(6,"('vid calling ftitle ',i2,' times for iip=',i2,
     +        ' (production run!)')") ititles(2),iip
          endif
          do i=1,abs(ititles(2))
            call ftitle(itest)
          enddo
        endif
c
c Open workstation WISS for freeze frames and/or kp graph:
c
        if (nfbegend.gt.0.or.kpgraph.gt.0) then
          call gopwk(idflash,luflash,3)
        endif
c
c Make kp line graph for later use:
c
        if (kpgraph.gt.0) then
          kpgraph = ifbkp
          call mkkpline(iproj(iip),jday1)
        endif
c
c Time loop:
c
        totmin = 1.e36
        totmax = -1.e36
        init = 1
        it = 1
 500    continue	! time loop
          ut = float(mtimes(2,it)) + float(mtimes(3,it))/60.
          utday = float(mtimes(1,it))*24. + float(mtimes(2,it)) +
     +            float(mtimes(3,it))/60.
c
c slt > 0 -> make frames at desired solar local time:
c
          slon = 0.
          if (slt.ge.0.) islt = ixslt(slt,ut,slon,gcmlon,imx,dlon)
          istart = 0
          if (it.eq.1) istart = 1
          if (it.eq.1.or.ninterp.le.0) then
            call getimesz(histvol,nvol,istart,slev,f,mtimes(1,it),
     +        imx,jmx,nfhist,iden,ipltgcm,lu,0,ivol,spv,ln,ier)
            if (ier.ne.0) then
              write(6,"('>>> vid: error ',i3,' from getimesz: it=',i3)") 
     +          ier,it
              stop 'getimesz'
            endif
            call glbxfer(f,imx,jmx,nfhist,plot,ip)
          endif
          if (it.gt.1) init=0
c
c If doing first or last frame of a field segment, save frame in
c flash buffer and repeat frame nfbegend times:
c
          ib = 0
          if (mkframes.gt.0.and.nfbegend.gt.0.and.
     +       (it.eq.1.or.it.eq.ntms)) then
            ib = ifbbeg
            if (it.eq.ntms) ib = ifbend
            write(6,"('Vid calling gflas1: it=',i3,
     +        ' will save next frame in flash buffer ',i3)") it,ib
            call gflas1(ib)
          endif
c
c Plot current time (at time it):
c
          if (iproj(iip).eq.1) then
            call pltglb(plot,imx,jmx,gcmlon,gcmlat,init,ut,
     +        cint(ip),cmin(ip),cmax(ip),slon,slt,slev,fldlab(ip),
     +        mtimes(1,it),info,totmin,totmax,icolor,labform(ip),
     +        icella,kpbar,kpgraph,jday1,mkframes,gcmzp(kmx),nfixclr)
          elseif (iproj(iip).eq.2) then
            call pltsat(plot,imx,jmx,gcmlon,gcmlat,init,ut,
     +        cint(ip),cmin(ip),cmax(ip),slon,slt,slev,fldlab(ip),
     +        mtimes(1,it),info,totmin,totmax,icolor,labform(ip),
     +        eradii,censatv,icella,kpbar,kpgraph,jday1,mkframes,
     +        gcmzp(kmx),nfixclr)
          endif
c
c Not making freeze frames:
          if (ib.eq.0) then
            if (kpgraph.gt.0) call gflas3(ifbkp)
            write(6,"('vid: frame ',i4,' it=',i3,' ut=',f6.3,' jday=',
     +        i5,' mt=',i2,':',i2,':',i2)") iframe+1,it,ut,
     +        jday1+mtimes(1,it)-1,(mtimes(i,it),i=1,3) 
            write(6,"(' ')")
            if (mkframes.gt.0) call frame
            iframe = iframe+1
c
c Make freeze frames (flash beg or end frame nfbegend times):
c
          else
            write(6,"('Vid calling gflas2: it=',i3,' (flashing ',i3,
     +        ' copies of frame at time ',3i3,' from buffer ',i3)") 
     +        it,nfbegend,(mtimes(i,it),i=1,3),ib
            write(6,"('     these will be frames ',i4,' to ',i4)")
     +        iframe+1,iframe+nfbegend+1
            call gflas2(ib)
            do i=1,nfbegend+1
              call gflas3(ib)
              if (kpgraph.gt.0) call gflas3(ifbkp)
              call frame
              iframe = iframe+1
            enddo
          endif
          it = it+1
c
c Interp loop:
c (plot0 is data at time it, plot1 is data at time it+1, plot is
c  interpolated data between plot0 and plot1, done by interp())
c
          if (ninterp.gt.0.and.it.le.ntms) then
            init=0
            plot0(:,:) = plot(:,:)
            call getimesz(histvol,nvol,0,slev,f,mtimes(1,it),
     +        imx,jmx,nfhist,iden,ipltgcm,lu,0,ivol,spv,ln,ier)
            if (ier.ne.0) then
              write(6,"('>>> vid: error ',i3,' from getimesz: it=',i3)") 
     +          ier,it
              stop 'getimesz'
            endif
            call glbxfer(f,imx,jmx,nfhist,plot1,ip)
            utday0 = utday
            utday1 = float(mtimes(1,it))*24. + float(mtimes(2,it)) +
     +               float(mtimes(3,it)) / 60.
            do iit=1,ninterp
              frac = float(iit) / float(ninterp+1)
              utday = frac * (utday1-utday0) + utday0
              call ut2mt(utday,mtime)      
              ut = float(mtime(2)) + float(mtime(3)) / 60.
              slon = 0.
              if (slt.ge.0.) islt = ixslt(slt,ut,slon,gcmlon,imx,dlon)
              call interp(plot0,plot1,plot,imx,jmx,frac,spv)
              if (iproj(iip).eq.1) then
                call pltglb(plot,imx,jmx,gcmlon,gcmlat,init,ut,
     +            cint(ip),cmin(ip),cmax(ip),slon,slt,slev,fldlab(ip),
     +            mtime,info,totmin,totmax,icolor,labform(ip),icella,
     +            kpbar,kpgraph,jday1,mkframes,gcmzp(kmx),nfixclr)
              elseif (iproj(iip).eq.2) then
                call pltsat(plot,imx,jmx,gcmlon,gcmlat,init,ut,
     +            cint(ip),cmin(ip),cmax(ip),slon,slt,slev,fldlab(ip),
     +            mtime,info,totmin,totmax,icolor,labform(ip),
     +            eradii,censatv,icella,kpbar,kpgraph,jday1,mkframes,
     +            gcmzp(kmx),nfixclr)
              endif
              if (kpgraph.gt.0) call gflas3(ifbkp)
              if (mkframes.gt.0) call frame
              iframe = iframe+1
              write(6,"('vid: frame ',i4,' it=',i3,' iit=',i3,' ut=',
     +          f6.3,' jday=',i5,' mt=',i2,':',i2,':',i2)") 
     +          iframe,it,iit,ut,jday1+mtime(1)-1,mtime
              write(6,"(' ')")
            enddo	! iit = 1,ninterp
            plot(:,:) = plot1(:,:) 
          endif		! if interp loop
        if (it.le.ntms) goto 500	! time loop
c
c Add blank frames between field segments:
c
        if (nfbetseg.gt.0.and.nfplt.gt.1.and.iip.lt.nfplt) then
          write(6,"('vid: making ',i3,' blank frames (frames ',i4,
     +      ' to ',i4,')')") nfbetseg,iframe+1,iframe+nfbetseg
          do i=1,nfbetseg
            call frame 
            iframe = iframe+1
          enddo
        endif
c
c Close workstation WISS used for freeze frames and/or kpgraph:
c (this must be closed before calling ftitle below)
c
        if (nfbegend.gt.0.or.kpgraph.gt.0) call gclwk(idflash)
c
c Add end titles if desired:
c
        if (ititles(3).ne.0.and.iip.eq.nfplt) then
          if (ititles(3).lt.0) then
            write(6,"('vid calling ftitle ',i2,' times for end',
     +        ' (test run)')") abs(ititles(3))
            itest = 0
          else
            write(6,"('vid calling ftitle ',i2,' times for end',
     +        ' (production run!)')") abs(ititles(3))
            itest = 1
          endif
          do i=1,abs(ititles(3))
            call ftitle(itest)
          enddo
        endif
 100  continue	! ip=1,nflds
c
c Clean up:
c      
      call hpdeallc(pf,ier,1)
      if (ier.ne.0) then
        write(6,"('hpdeallc error = ',i6,' for pf')") ier
        stop 'pf'
      endif
      call clsgks
c
      stop
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine glbxfer(f,imx,jmx,nfhist,plot,ip)
      dimension f(imx,jmx,nfhist),plot(imx,jmx)
      common/timeix/itxt,itxu,itxv,itxo2,itxox,itxn4s,itxnoz,itxco,
     +              itxco2,itxh2o,itxh2,itxhox,itxop,itxch4,itxo21d,
     +              itxno2,itxno,itxo3,itxo1,itxoh,itxho2,itxh,
     +              itxn2d,itxti,itxte,itxne,itxo2p,itxw,itxz
c
c Define plot from f at time it, field ip.  
c Note nfhist = fields from getimesz (31 hist fields + N2 + FOF2)
c
      itxn2 = itxz+1
      if (ip.le.nfhist) then
        plot(:,:) = f(:,:,ip)
c
c rho (total density):
c (will need to check for spval after adding height interp to code)
c
      elseif (ip.eq.nfhist+1) then
        itxn2 = nfhist
        plot(:,:) = f(:,:,itxo2) + f(:,:,itxo1) + f(:,:,itxn2)
      else
        write(6,"('>>> glbxfer: unknown ip=',i3,' nfhist=',i3,
     +    ' returning zero field')")
     +    ip,nfhist
        plot(:,:) = 0.
      endif
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine ut2mt(ut,mt)
      dimension mt(3)
c
      mt(1) = ut / 24.
      mt(2) = ut - mt(1)*24.
      mt(3) = nint(60. * (ut - ifix(ut)))
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine interp(plot0,plot1,plot,imx,jmx,frac,spv)
c
      dimension plot0(imx,jmx),plot1(imx,jmx),plot(imx,jmx)
c
      rmin0 = 1.e36
      rmax0 = -1.e36
      do j=1,jmx
        do i=1,imx
          if (plot0(i,j).lt.rmin0) then
            rmin0 = plot0(i,j)
            imin0 = i
            jmin0 = j 
          endif
          if (plot0(i,j).gt.rmax0) then
            rmax0 = plot0(i,j)
            imax0 = i
            jmax0 = j
          endif
        enddo
      enddo
      rmin1 = 1.e36
      rmax1 = -1.e36
      do j=1,jmx
        do i=1,imx
          if (plot1(i,j).lt.rmin1) then
            rmin1 = plot1(i,j)
            imin1 = i
            jmin1 = j 
          endif
          if (plot1(i,j).gt.rmax1) then
            rmax1 = plot1(i,j)
            imax1 = i
            jmax1 = j
          endif
        enddo
      enddo
      write(6,"(' ')")
      write(6,"('interp: frac=',f6.3)") frac
      write(6,"('rmin0=',e12.4,' imin0,jmin0=',2i3)") rmin0,imin0,jmin0
      write(6,"('rmax0=',e12.4,' imax0,jmax0=',2i3)") rmax0,imax0,jmax0
      write(6,"('rmin1=',e12.4,' imin1,jmin1=',2i3)") rmin1,imin1,jmin1
      write(6,"('rmax1=',e12.4,' imax1,jmax1=',2i3)") rmax1,imax1,jmax1
c
      do j=1,jmx
        do i=1,imx
          plot(i,j) = spv
          if (plot0(i,j).ne.spv.and.plot1(i,j).ne.spv) 
     +      plot(i,j) = frac * (plot1(i,j) - plot0(i,j)) + plot0(i,j)
          if ((i.eq.imin0.and.j.eq.jmin0).or.(i.eq.imax0.and.j.eq.jmax0)
     +    .or.(i.eq.imin1.and.j.eq.jmin1).or.(i.eq.imax1.and.j.eq.jmax1)
     +    )
     +      write(6,"('i j=',2i3,' plot0,1=',2e12.4,' plot=',e12.4)")
     +        i,j,plot0(i,j),plot1(i,j),plot(i,j)
        enddo
      enddo

      rmin = 1.e36
      rmax = -1.e36
      do j=1,jmx
        do i=1,imx
          if (plot(i,j).lt.rmin) then
            rmin = plot(i,j)
            imin = i
            jmin = j 
          endif
          if (plot(i,j).gt.rmax) then
            rmax = plot(i,j)
            imax = i
            jmax = j
          endif
        enddo
      enddo
      write(6,"('rmin (of plot)=',e12.4,' imin,jmin=',2i3)") 
     +  rmin,imin,jmin
      write(6,"('rmax (of plot)=',e12.4,' imax,jmax=',2i3)") 
     +  rmax,imax,jmax
      write(6,"('plot0,1 (at imin,jmin)=',2e12.4,' plot=',e12.4)") 
     +  plot0(imin,jmin),plot1(imin,jmin),plot(imin,jmin)
      write(6,"('plot0,1 (at imax,jmax)=',2e12.4,' plot=',e12.4)") 
     +  plot0(imax,jmax),plot1(imax,jmax),plot(imax,jmax)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mkkpline(iproj,jday1)
c
c Store line plot of kp in flash buffer ib for later use by pltglb
c and pltsat, etc.
c
      parameter (mxgpi=200,ndays=10,ntmstot=ndays*24+1)
      common/kp/ kpbar,kpgraph,ngpi,gpikp(mxgpi),igpiday(mxgpi),
     +  gpiut(mxgpi),nkp,vpkpgce(4),vpkpgsv(4),rkp(ntmstot)
      dimension view(4),xx(ntmstot)
      data nkp/81/ ! through ut 1.5 day 79090
c
c Note view left and right should match viewport left and right in pltglb
      data xlaboff/.03/, ylaboff/.03/, ticlen/.014/
c
      if (iproj.eq.1) then
        view(:) = vpkpgce(:)
      elseif (iproj.eq.2) then
        view(:) = vpkpgsv(:)
      else
        write(6,"('>>> mkkpline: unknown iproj=',i3)") iproj
        return 
      endif
      charsize = .03*(view(2)-view(1))
      uthr = -1.
      iday = 1.
      do it=1,ntmstot
        xx(it) = float(it)
        uthr = uthr+1.
        if (uthr.gt.23.) then
          uthr = 0.
          iday = iday+1 
        endif
        utinc = (iday-1)*24.+uthr
        do i=1,ngpi-1
          utinc0 = float(igpiday(i)-jday1)*24.+gpiut(i)
          utinc1 = float(igpiday(i+1)-jday1)*24.+gpiut(i+1)
          if (utinc.lt.utinc0) then
            rkp(it) = gpikp(i)
            write(6,"('utinc < utinc0: it=',i3,' rkp=',f7.2)") 
     +        it,rkp(it)
            goto 101
          endif
          if (utinc.ge.utinc0.and.utinc.le.utinc1) then
            rkp0 = gpikp(i)
            rkp1 = gpikp(i+1)
            goto 100
          endif
        enddo
        write(6,"('mkkpline: could not bracket utinc=',f10.2,' iday=',
     +    i3,' it=',i4,' utinc0,1=',2f7.2,' jday1=',i5)") 
     +    utinc,iday,it,utinc0,utinc1
        stop 'kp'
 100    continue
        rkp(it) = (utinc-utinc0)*(rkp1-rkp0) / (utinc1-utinc0) + rkp0
c       write(6,"('mkkpline: it=',i3,' utinc=',f7.2,' utinc0,1=',
c    +    2f8.2,' rkp0,1=',2f7.2,' rkp(it)=',f7.2)")
c    +    it,utinc,utinc0,utinc1,rkp0,rkp1,rkp(it)
 101    continue
      enddo
      call gflas1(kpgraph)
      call set(view(1),view(2),view(3),view(4),
     +  1.,float(ntmstot),0.,7.,1)
      call agsetf('SET.',4.)
      call agsetf('FRA.',2.)
      call agsetf('LEFT/TYPE.',0.)
      call agsetf('RIGHT/TYPE.',0.)
      call agsetf('BOTTOM/TYPE.',0.)
      call agsetf('TOP/TYPE.',0.)
      call agsetf('LABEL/CONTROL.',0.)
      call agsetf('LEFT/MAJOR.',0.)
      call agsetf('RIGHT/MAJOR.',0.)
      call agsetf('BOTTOM/MAJOR.',0.)
      call agsetf('TOP/MAJOR.',0.)
      call ezxy(xx,rkp,ntmstot,' ')
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      do i=1,9
        xpos = view(1)+float(i)/10.*(view(2)-view(1))
        call line(xpos,view(4),xpos,view(4)-ticlen)
        call line(xpos,view(3),xpos,view(3)+ticlen)
      enddo
      call plchhq(view(1)-xlaboff,view(3),'0',charsize,0.,1.)
      call plchhq(view(1)-xlaboff,view(4),'7',charsize,0.,1.)
      call plchhq(view(1)-xlaboff,0.5*(view(4)+view(3)),'KP',
     +  charsize,0.,1.)
      call plchhq(view(1),view(3)-ylaboff,'79',charsize,0.,0.)
      call plchhq(view(2),view(3)-ylaboff,'89',charsize,0.,0.)
      call plchhq(0.5*(view(1)+view(2)),view(3)-ylaboff,'DAY',
     +  charsize,0.,0.)
      write(6,"('mkkpline calling gflas2(',i3,')')") kpgraph
      call gflas2
      return
      end