
      program tigcmdif
c
c ti(e)gcm difference processor:
c
      include 'tigcmdif.h'
      include 'tgcmhdr.h'
      dimension iasf(13),fkmx0(kmx),fkmx1(kmx)
      character*80 datflnm
      logical exists
      data iasf/13*1/
c
c fields() will be perturbed, fields0() will be control. Differences
c will be stored in fields: fields = fields - fields0
c
      pointer (pfields,fields(1)), (pfields0,fields0(1))
      pointer (putloc,futloc(1)), (putlocht,futlocht(1))
      pointer (phtutpert,htutpert(1)), (phtutcntr,htutcntr(1))
      pointer (putlat,futlat(1))
c
c ixgcm = initialize field indices
c
      do i=1,nftot
        ixgcm(i) = i
      enddo
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
c
c Set up gks, conpack, and autograph:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
      call cpsetup(cpspval)
      call cpseti('LLB',1)
      if (icolor.gt.0) then
        call cpseti('LLB',3)
        call cpseti('LBC',1)
      endif
      call setag
      call agseti('AXIS/LEFT/CONTROL.',1)
      call agseti('AXIS/BOTTOM/CONTROL.',1)
c
c Allocate space for perturbed and control fields:
c
      call hpalloc(pfields,imx*kmx*jmx*nfget,ier,1)
      call hpalloc(pfields0,imx*kmx*jmx*nfget,ier,1)
c
c Allocate for ut vs zp contours (futloc(ntms,kmx,nfplt,nloc)):
c
      if (iutloc.gt.0.or.iutline.gt.0) then
        if ((iutloc.gt.0.and.nzprange.gt.0).or.
     +      (iutline.gt.0.and.npls.gt.0))
     +    call hpalloc(putloc,ntms*kmx*nfplt*nloc,ier,1)
        if ((iutloc.gt.0.and.nhtscale.gt.0).or.
     +      (iutline.gt.0.and.nhts.gt.0))
     +    call hpalloc(putlocht,ntms*nhtscale*nfplt*nloc,ier,1)
        call hpalloc(phtutpert,ntms*kmx*nloc,ier,1)
        call hpalloc(phtutcntr,ntms*kmx*nloc,ier,1)
      endif
c
c Allocate for ut vs latitude contours (futlat(ntms,jmx,nfplt,nzphtlon)):
c
      if (iutlat.gt.0) then
        call hpalloc(putlat,ntms*jmx*nfplt*nzphtlon,ier,0)
        if (ier.ne.0) then
          write(6,"('>>> Error from hpalloc for utlat: increase ',
     +      'memory allocations')")
          stop 'mem'
        endif
      endif 
c
c Ut loop:
c
      istartp = 1
      istartc = 1
      ipvolp = 1
      icvolp = 1
      iprhdr = 1
      do it=1,ntms
c
c Get perturbed fields:
c
        call getgcmf(pertvol,npvols,istartp,mtpert(1,it),fields,nfget,
     +    ifget,imx,kmx,jmx,nfhist,iden,lupert,iprhdr,ipvol,
     +    istiegcm,ier)   
        if (ier.ne.0) then
          write(6,"('tigcmdif: error ',i3,' from getgcmf for perturbed',
     +      ' time: it=',i3,' mtpert=',i2,':',i2,':',i2)") 
     +    ier,it,(mtpert(i,it),i=1,3)
          stop 'getpert'
c       else
c         write(6,"('tigcmdif: found perturbed history ',
c    +      i2,':',i2,':',i2,' on vol ',a,' it=',i3,' istiegcm=',l1)") 
c    +      (mtpert(i,it),i=1,3),
c    +      pertvol(ipvol)(1:lenstr(pertvol(ipvol))),it,istiegcm
        endif
        if (.not.istiegcm.and.(ifplt(ixpot).gt.0.or.ifplt(ixui).gt.0.or.
     +      ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0)) 
     +    call dynoff(ifplt(ixpot),ifplt(ixui),ifplt(ixvi),ifplt(ixwi),
     +      pertvol(ipvol))
c
c Get control fields:
c Note if pert vol has been incremented (ipvol > ipvolp), we want search
c  for unpert history to continue from previous vol, so pass it istart < 0
c  with abs(istart) as desired volume index (see 9/93 change to getgcmf.f
c  re istart and iv1)
c
        if (ipvol.gt.ipvolp) istartc = -(ipvol-1)
        call getgcmf(cntrvol,ncvols,istartc,mtcntr(1,it),fields0,nfget,
     +    ifget,imx,kmx,jmx,nfhist,iden,lucntr,iprhdr,icvol,
     +    istiegcm,ier)   
        if (ier.ne.0) then
          write(6,"('tigcmdif: error ',i3,' from getgcmf for control',
     +      ' time: it=',i3,' mtcntr=',i2,':',i2,':',i2)") 
     +    ier,it,(mtcntr(i,it),i=1,3)
          stop 'getcntr'
c       else
c         write(6,"('tigcmdif: found control history ',
c    +      i2,':',i2,':',i2,' on vol ',a,' it=',i3,' istiegcm=',l1)") 
c    +      (mtcntr(i,it),i=1,3),
c    +      cntrvol(icvol)(1:lenstr(cntrvol(icvol))),it,istiegcm
        endif
        if (.not.istiegcm.and.(ifplt(ixpot).gt.0.or.ifplt(ixui).gt.0.or.
     +      ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0)) 
     +    call dynoff(ifplt(ixpot),ifplt(ixui),ifplt(ixvi),ifplt(ixwi),
     +      cntrvol(icvol))
c
c Get ion drifts if needed (diffs will be found later)
c (pointers to ui,vi,wi in common in tigcmdif.h):
c
        if ((imap.gt.0.and.ipuv.gt.0.and.ifplt(ixpot).gt.0).or.
     +    ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0.or.
     +    ifplt(ixuivi).gt.0) then
          if (ifplt(ixuivi).gt.0.or.ifplt(ixui).gt.0.or.
     +       (ipuv.gt.0.and.ifplt(ixpot).gt.0)) then
            call hpalloc(pglbui,imx*kmx*jmx,iier,1)
            call mkdrifts(fields,1,glbui,ier)
            call hpalloc(pglbui0,imx*kmx*jmx,iier,1)
            call mkdrifts(fields0,1,glbui0,ier)
            if (iier.ne.0) then
              write(6,"('>>> tigcmdif: error ',i3,' from hpalloc ',
     +          'for ui')") iier
              stop 'memory'
            endif
            if (ier.ne.0) then
              write(6,"('>>> error ',i2,' from mkdrifts for ui <<<')") 
     +          ier
              glbui(:,:,:) = 0.
              glbui0(:,:,:) = 0.
            endif
          endif
          if (ifplt(ixuivi).gt.0.or.ifplt(ixvi).gt.0.or.
     +       (ipuv.gt.0.and.ifplt(ixpot).gt.0)) then
            call hpalloc(pglbvi,imx*kmx*jmx,iier,1)
            call mkdrifts(fields,2,glbvi,ier)
            call hpalloc(pglbvi0,imx*kmx*jmx,iier,1)
            call mkdrifts(fields0,2,glbvi0,ier)
            if (iier.ne.0) then
              write(6,"('>>> tigcmdif: error ',i3,' from hpalloc ',
     +          'for vi')") iier
              stop 'memory'
            endif
            if (ier.ne.0) then
              write(6,"('>>> error ',i2,' from mkdrifts for vi <<<')") 
     +          ier
              glbvi(:,:,:) = 0.
              glbvi0(:,:,:) = 0.
            endif
          endif
          if (ifplt(ixwi).gt.0) then
            call hpalloc(pglbwi,imx*kmx*jmx,iier,1)
            call mkdrifts(fields,3,glbwi,ier)
            call hpalloc(pglbwi0,imx*kmx*jmx,iier,1)
            call mkdrifts(fields0,3,glbwi0,ier)
            if (iier.ne.0) then
              write(6,"('>>> tigcmdif: error ',i3,' from hpalloc ',
     +          'for wi')") iier
              stop 'memory'
            endif
            if (ier.eq.0) then
c             write(6,"('tigcmdif: got vertical ion velocity')")
            else
              write(6,"('>>> error ',i2,' from mkdrifts for wi <<<')") 
     +          ier
              glbwi(:,:,:) = 0.
              glbwi0(:,:,:) = 0.
            endif
          endif
        endif
c
c If doing vpar or wpar, need global dip and declination:
c
        if (it.eq.1.and.(ifplt(ixvpar).gt.0.or.ifplt(ixwpar).gt.0)) 
     +    call getdipdec(dipdec,imx,jmx,date(1),ludipdec,gcmlat,gcmlon)
c
c Set up labels:
c
        if (it.eq.1) call setflab(istiegcm,iden,flab,nftot)
c
c Contour CE projections:
c
        if (nglb.gt.0) call pltglb('CE',ut(it),fields,fields0)
        if (nmolw.gt.0) call pltglb('MO',ut(it),fields,fields0)
c
c Polar stereographic plots:
c
        if (npollat.gt.0) call pltpol(ut(it),fields,fields0)
c
c Satellite view projection:
c
        if (nsatv.gt.0) call pltsatv(ut(it),fields,fields0)
c
c Longitude slices:
c
        if ((nlon.gt.0.or.nslt.gt.0).and.nzprange.gt.0)
     +    call pltlon(ut(it),0,fields,fields0)
        if ((nlon.gt.0.or.nslt.gt.0).and.nhtscale.gt.0)
     +    call pltlon(ut(it),1,fields,fields0)
c
c Latitude slices:
c
        if (nlat.gt.0.and.nzprange.gt.0) 
     +    call pltlat(ut(it),0,fields,fields0)
        if (nlat.gt.0.and.nhtscale.gt.0) 
     +    call pltlat(ut(it),1,fields,fields0)
c
c Location line plots:
c
        if (nloc.gt.0) call setag
        if (nloc.gt.0.and.nzprange.gt.0.and.lineloc.gt.0) 
     +    call pltloc(ut(it),0,fields,fields0)
        if (nloc.gt.0.and.nhtscale.gt.0.and.lineloc.gt.0) 
     +    call pltloc(ut(it),1,fields,fields0)
c
c Update diffs for location contours (ut vs zp) and/or 
c line plots with ut on x-axis (at selected zp's):
c
        if (iutloc.gt.0.or.iutline.gt.0)
     +    call getutloc(it,fields,fields0,futloc,futlocht,
     +      htutpert,htutcntr)
c
c Update diffs for ut vs lat contours:
c
        if (iutlat.gt.0) 
     +    call getutlat(it,fields,fields0,futlat)
c
c End ut loop:
c
        istartp = 0
        istartc = 0
        ipvolp = ipvol
        icvolp = icvol
        iprhdr = 0
        if ((imap.gt.0.and.ipuv.gt.0.and.ifplt(ixpot).gt.0).or.
     +    ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0) then
          if (ifplt(ixui).gt.0.or.(ipuv.gt.0.and.ifplt(ixpot).gt.0))then 
            call hpdeallc(pglbui,ier,1)
            call hpdeallc(pglbui0,ier,1)
          endif
          if (ifplt(ixvi).gt.0.or.(ipuv.gt.0.and.ifplt(ixpot).gt.0))then 
            call hpdeallc(pglbvi,ier,1)
            call hpdeallc(pglbvi0,ier,1)
          endif
          if (ifplt(ixwi).gt.0) then
            call hpdeallc(pglbwi,ier,1)
            call hpdeallc(pglbwi0,ier,1)
          endif
          if (ier.ne.0) 
     +      write(6,"('>>> tigcmdif: error ',i3,' from hpdeallc for ',
     +                'ion drifts')") ier
        endif
      enddo
c
      call hpdeallc(pfields,ier,1)
      call hpdeallc(pfields0,ier,1)
c
c Make ut vs zp loc contours and/or ut vs field line plots:
c 
      if (iutloc.gt.0.or.iutline.gt.0) then
        if (iutloc.gt.0) then
          if (nzprange.gt.0) 
     +      call pltutloc(futloc,futlocht,htutpert,htutcntr,0)
          if (nhtscale.gt.0) 
     +      call pltutloc(futloc,futlocht,htutpert,htutcntr,1)
        endif
        if (iutline.gt.0.and.(npls.gt.0.or.nhts.gt.0)) 
     +      call pltutlin(futloc,futlocht,htutpert,htutcntr)
c
c Release space used by ut plots:
c
        if ((iutloc.gt.0.and.nzprange.gt.0).or.
     +      (iutline.gt.0.and.npls.gt.0)) call hpdeallc(putloc,ier,1)
        if ((iutloc.gt.0.and.nhtscale.gt.0).or.
     +      (iutline.gt.0.and.nhts.gt.0)) call hpdeallc(putlocht,ier,1)
        call hpdeallc(phtutpert,ier,1)
        call hpdeallc(phtutcntr,ier,1)
      endif
c
c Make ut vs latitude contours:
c
      if (iutlat.gt.0) call pltutlat(futlat)
c
c Return gmeta and ascii data files to remote unix directories:
c
      call clsgks
      if (iwrascii.eq.2.and.lenstr(dirascii).gt.0)
     +  call rcpfile(luascii,'tigcmdif.dat',dirascii)
      if (lenstr(dircgm).ne.0) call rcpfile(0,'gmeta',dircgm)
c
c Send ascii data file(s) to mss if desired:
c
      lenmss = lenstr(mssascii)
      if (iwrascii.gt.0.and.lenmss.gt.0) then
        if (iwrascii.eq.1) then
          do i=1,iframe
            write(datflnm,"('tigcmdat.',i3.3,'.dat')") i
            call putmss(datflnm,mssascii(1:lenmss)//'/'//datflnm)
          enddo  
        elseif (iwrascii.eq.2) then
          call putmss('tigcmdif.dat',
     +      mssascii(1:lenmss)//'/tigcmdif.dat')
        endif
      endif
      stop 'success!'
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine dynoff(i1,i2,i3,i4,vol)
      character*(*) vol
      write(6,"(' ')")
      write(6,"('!!!!!!!!!!!!!!!!!!!! NOTICE !!!!!!!!!!!!!!!!!!!')")
      write(6,"(' ')")
      write(6,"('HISTORY VOLUME ',a,' IS NOT A TIEGCM',
     +  ' VOLUME ',/,'I AM TURNING OFF POTENTIAL AND/OR ION ',
     +  'DRIFT PLOTS')") vol(1:lenstr(vol))
      i1 = 0
      i2 = 0
      i3 = 0
      i4 = 0
      write(6,"(' ')")
      write(6,"('!!!!!!!!!!!!!!!!!!!! NOTICE !!!!!!!!!!!!!!!!!!!')")
      write(6,"(' ')")
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getdipdec(dipdec,imx,jmx,iyr,lu,gcmlat,gcmlon)
      dimension dipdec(imx,jmx,2),gcmlat(jmx),gcmlon(imx)
      character*40 flnm
      data flnm/'/u0/foster/bin/dipdec.gcmgrid.bin       '/
      parameter (pi = 3.14159)
      parameter (dtr = pi/180.)
      logical exists
c
      lenflnm = lenstr(flnm)
      inquire(file=flnm(1:lenflnm),exist=exists)
      if (.not.exists) then
        write(6,"('>>> Cannot find file ',a,' -- will call igrf')") 
     +    flnm(1:lenflnm)
        yr = float(iyr)
        if (yr.gt.1985.) yr = 1985.
        write(6,"('Calling igrf globally to obtain dip,dec',
     +    ' for calculation of vpar and/or wpar: yr=',f7.1)") yr
        do j=1,jmx
          do i=1,imx
            call igrf(gcmlat(j),gcmlon(i),300.,yr,gmlat,gmlon,
     +        dipdec(i,j,1),dipdec(i,j,2))
            dipdec(i,j,1) = dipdec(i,j,1)*dtr
            dipdec(i,j,2) = dipdec(i,j,2)*dtr
          enddo
        enddo
c       call pltce(dipdec,'INCLINATION')
c       call pltce(dipdec(1,1,2),'DECLINATION')
      else
        call asign(flnm(1:lenflnm),lu,ier)
        read(lu) dipdec
        write(6,"('Successful read of dipdec from ',a)") flnm(1:lenflnm)
        close(lu)
c       call pltce(dipdec,'INCLINATION')
c       call pltce(dipdec(1,1,2),'DECLINATION')
      endif
      return
      end
