c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timesdif/timedif.f
c------------------------------------------------------------------
c
      program timesdif
c
c Difference processor for timegcm model 
c (-13 to +5 column with 31 fields including n2):
c
      include 'gettime.h'
      include 'timesdif.h'
      logical found
      dimension iasf(13),mt(3)
      pointer (ppnt,pnt(imx,kmx,jmx,ntimefld)),(pfld,fld(1))
      data iasf/13*1/, luhist/8/, iprhdr/1/, iframe/0/
c
c Initialization: set up tgcm grid, get user input, set
c                 some graphics default parameters:
c
      call setgrid(gcmlat,glat1,dlat,jmx, gcmlon,glon1,dlon,imx,
     +             gcmzp,zp1,dzp,kmx)
      call getinp
      call cpsetup(cpspval)	! in libplt.a
      do ip=1,nflds
        ixtime(ip) = ip
      enddo
c
c Allocate space for global array with only desired fields:
c nplfld = number of desired fields
c
      call hpalloc(pfld,imx*kmx*jmx*nplfld,ier,1)
      if (ier.ne.0) then
        write(6,"('hpalloc error = ',i6,' for fld')") ier
        stop 'fld'
      endif
c
c Set up gks:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
c
c UT loop:
c
      istart = 1
      do it=1,ntms
c
c Allocate space for all global fields at single ut:
c
        call hpalloc(ppnt,imx*kmx*jmx*ntimefld,ier,1)
        if (ier.ne.0) then
          write(6,"('hpalloc error = ',i6,' for pnt')") ier
          stop 'pnt'
        endif
c
c Get perturbed history:
c
        mt(1) = mdpert(it)
        mt(2) = mhpert(it)
        mt(3) = mmpert(it)
        call gettimes(pertvol,npvol,istart,mt,pnt,imx,kmx,jmx,
     +    ntimefld,iden,luhist,iprhdr,ivol,ier)
        if (ier.ne.0) then
          write(6,"('>>> timedif: error ',i3,' from gettimes for time ',
     +      3i4)") ier,(mt(i),i=1,3)
          stop 'gettimes'
        endif
        lenpertv = lenstr(pertvol(ivol))
        write(6,"('timedif: obtained timegcm for perturbed time ',3i4,
     +    ' from vol ',a)") (mt(i),i=1,3), pertvol(ivol)(1:lenpertv)
        call tail(pertvol(ivol)(1:lenpertv),cpertvol)
c
c Save perturbed fields in fld (only fields to be plotted):
c fld(imx,kmx,jmx,nplfld)
c
        ixp = 0
        do ip=1,ntimefld
          if (ipltime(ip).gt.0) then
            ixp = ixp+1
            do j=1,jmx
              do k=1,kmx
                do i=1,imx
                  fld((ixp-1)*jmx*kmx*imx+(j-1)*kmx*imx+(k-1)*imx+i) =
     +              pnt(i,k,j,ip)
                enddo
              enddo
            enddo
          endif
        enddo
c
c Get control history:
c
        mt(1) = mdcntr(it)
        mt(2) = mhcntr(it)
        mt(3) = mmcntr(it)
        call gettimes(cntrvol,ncvol,istart,mt,pnt,imx,kmx,jmx,
     +    ntimefld,iden,luhist,iprhdr,ivol,ier)
        if (ier.ne.0) then
          write(6,"('>>> timedif: error ',i3,' from gettimes for time ',
     +      3i4)") ier,(mt(i),i=1,3)
          stop 'gettimes'
        endif
        lencntrv = lenstr(cntrvol(ivol))
        write(6,"('timedif: obtained timegcm for control time ',3i4,
     +    ' from vol ',a)") (mt(i),i=1,3),cntrvol(ivol)(1:lencntrv)
        call tail(cntrvol(ivol)(1:lencntrv),ccntrvol)
        istart = 0
        iprhdr = 0
c
c Make diffs (fld = fld-pnt):
c
        write(6,"(' ')")
        write(6,"('Making diffs...')")
        call mkdifs(fld,pnt)
c
c If not plotting at selected heights, we can release pnt:
c
        if (nhts.le.0) then
          call hpdeallc(ppnt,ier,1)
          if (ier.ne.0) then
            write(6,"('hpdeallc error = ',i6,' for ppnt')") ier
            stop 'ppnt'
          endif
        endif
c
c Global cylindrical equidistant plots at selected zp's:
c
        if (iplglb.gt.0.and.npls.gt.0) call pltglb(ut(it),fld)
c
c Polar stereographic plots at selected zp's:
c
        if (npollat.gt.0.and.npls.gt.0) call pltpol(ut(it),fld)
c
c Latitude slices:
c
        if (ipllat.gt.0) call pltlat(ut(it),fld,0)
c
c Longitude slices (zp on y-axis):
c
        if (ipllon.gt.0) call pltlon(ut(it),fld,0)
c
c Make line plots of zp profiles at selected locations:
c
        if (iplloc.gt.0) call pltloc(0,fld)
c
c Make line plots of zp profiles at global means:
c
        if (iplglbm.gt.0) call pltloc(1,fld)
c
c If plotting at selected heights, first recover perturbed fields:
c (At this point, fld contains diffs, pnt contains control;
c  after recpert, fld will contain pert and pnt will still contain control)
c
        if (nhts.gt.0.or.ipllonht.gt.0) then
          call recpert(fld,pnt)
c
c Cyl Equid plots at selected heights:
c
          if (iplglb.gt.0.and.nhts.gt.0) call pltglbht(ut(it),fld,pnt)
          if (npollat.gt.0.and.nhts.gt.0) call pltpolht(ut(it),fld,pnt)
c
c Longitude slices (ht on y-axis):
c
          if (ipllonht.gt.0) then
c
c Now interpolate both pert and cntr to htscale 
c (use 2nd dim of fld and pnt):
c
            call htint(fld,pnt,htscale,nhtscale)
c
c Make differences at selected heights (middle dimension beyond nhtscale
c   should henceforth be ignored):
c
            write(6,"(' ')")
            write(6,"('Making diffs for linear height scale...')")
            call mkdifs(fld,pnt)
            if (ht0.eq.0..and.ht1.eq.0..and.dht.eq.0.) then
              write(6,"('timedif: defaulting height scale')")
              ht0 = 30.
              ht1 = 500.
              dht = 20.
            endif
            nhtscale = 0
            do i=1,500
              if (ht0 + (i-1)*dht .gt. ht1) goto 500
              nhtscale = nhtscale+1
            enddo
 500        continue
            call hpalloc(phtscale,nhtscale,ier,1)
            if (ier.ne.0) then
              write(6,"('hpalloc error = ',i6,' for htscale')") ier
              stop 'htscale'
            endif
            do i=1,nhtscale
              htscale(i) = ht0 + (i-1)*dht
            enddo
            call pltlat(ut(it),fld,htscale,nhtscale,1)
          endif
c
c Free pnt:
          call hpdeallc(ppnt,ier,1)
          if (ier.ne.0) then
            write(6,"('hpdeallc error = ',i6,' for ppnt')") ier
            stop 'ppnt'
          endif

        endif
c
c Time loop:
      enddo
c
c Release selected fields array:
c
      call hpdeallc(pfld,ier,1)
      if (ier.ne.0) then
        write(6,"('hpdeallc error = ',i6,' for pfld')") ier
        stop 'pfld'
      endif
      call clsgks
c
      stop
      end