
      program tigcmpro
c
c ti(e)gcm processor:
c
      include 'tigcmproc.h'
      include 'tgcmhdr.h'
      dimension iasf(13)
      data iasf/13*1/
      data luhist/8/, lumag/99/, luascii/90/, ludipdec/20/
     +     iframe/0/, zmflag/999./, r12flag/999./, gmflag/999./
      pointer (pfields,fields(1))
c     pointer (phtscale,htscale(1))
      logical exists
      data labshort
     + /'TN      ','UN      ','VN      ','O2      ','O1      ',  
     +  'N4S     ','NO      ','O+      ','N2D     ','TI      ',  
     +  'TE      ','NE      ','O2+     ','W       ','Z       ',  
     +  'POT     ','N2      ','UI      ','VI      ','WI      ',
     +  'FOF2    ','HMF2    ','RHO     ','UN+VN   ','UI+VI   ',
     +  'O/O2    ','O/N2    ','O2/N2   ','O/O2+N2 '/
c
      data ilog
c          tn  un  vn  o2  o1  n4s  no  o+  n2d  ti  te  ne  o2+  w  z
     +    / 0,  0,  0,  1,  1,   1,  1,  1,   1,  0,  0,  1,  1,  0, 0,
c          pot  n2  ui  vi  wi fof2 hmf2 rho unvn uivi o/o2 o/n2 o2/n2
     +       0,  1,  0,  0,  0,   0,   0,  1,   0,   0,   0,   0,    0,
c          o/(o2+n2)
     +       0/
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
      if (iwrascii.gt.0) then
        inquire(file='tigcmproc.dat',exist=exists)
        if (exists) istat = ishell('rm tigcmproc.dat')
      endif
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',3)
      call cpseti('LLB',1)
      if (icolor.gt.0) then
        call cpseti('LLB',3)
        call cpseti('LBC',1)
      endif
c     call cpsetr('LLP',2)
c     call cpsetr('RC2',.5)
      call cpsetr('PC6',.5)
      call setag
      call agseti('AXIS/LEFT/CONTROL.',1)
      call agseti('AXIS/BOTTOM/CONTROL.',1)
c
c Allocate space for desired fields:
c
      call alloc(pfields,imx*kmx*jmx*nfget)
c
c Ut loop:
c
      istart = 1
      iprhdr = 1
      do it=1,ntms
        call getgcmf(histvol,nhvols,istart,mtimes(1,it),fields,nfget,
     +    ifget,imx,kmx,jmx,nfhist,iden,luhist,iprhdr,ivol,istiegcm,ier)   
        if (ier.ne.0) then
          write(6,"('tigcmproc: error ',i3,' from getgcmf: it=',i3,
     +      ' mtimes=',i2,':',i2,':',i2)") ier,it,(mtimes(i,it),i=1,3)
          goto 100
        else
          write(6,"('tigcmproc: getgcmf found history ',i2,':',i2,':',
     +      i2,' on vol ',a,' it=',i3,' istiegcm=',l1)") 
     +      (mtimes(i,it),i=1,3),histvol(ivol)(1:lenstr(histvol(ivol))),
     +      it,istiegcm
        endif
        if (it.eq.1) call setflab(istiegcm,iden,flab,nftot,ionvel)
        call getions(fields,it)
c
c Ascii files of each field at global grid:
c
        if (iwrgrid.gt.0) call wrgrid(fields,it,'tigcmproc')
c
c Contour CE and/or MO projections:
c
        if (nglb.gt.0) call pltglb('CE',ut(it),fields)
        if (nmolw.gt.0) call pltglb('MO',ut(it),fields)
c
c Polar stereographic plots:
c
        if (npollat.gt.0) call pltpol(ut(it),fields)
c
c Satellite view projection:
c
        if (nsatv.gt.0) call pltsatv(ut(it),fields)
c
c Longitude slices:
c
        if (nlon.gt.0.or.nslt.gt.0) call pltlon(ut(it),fields)
c
c Latitude slices:
c
        if (nlat.gt.0) call pltlat(ut(it),fields)
c
c Location line plots:
c
        if (nloc.gt.0) call pltloc(ut(it),mtimes(1,it),fields)
c
c Free ion drifts space if it was used in this iteration:
        if (ifplt(ixui).gt.0.or.ifplt(ixuivi).gt.0.or.
     +     (ipuv.gt.0.and.ifplt(ixpot).gt.0))call hpdeallc(pglbui,ier,1)
        if (ifplt(ixvi).gt.0.or.ifplt(ixuivi).gt.0.or.
     +     (ipuv.gt.0.and.ifplt(ixpot).gt.0))call hpdeallc(pglbvi,ier,1)
        if (ifplt(ixwi).gt.0) call hpdeallc(pglbwi,ier,1)
c
c End ut loop:
        istart = 0
        iprhdr = 0
      enddo
c
 100  continue
      call hpdeallc(pfields,ier,1)
      call clsgks
      if (iwrascii.eq.2.and.lenstr(dirascii).gt.0) 
     +  call rcpfile(luascii,'tigcmproc.dat',dirascii)
      if (lenstr(dircgm).ne.0) call rcpfile(0,'gmeta',dircgm)
      stop 'success!'
      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
      else
        call asign(flnm(1:lenflnm),lu,ier)
        read(lu) dipdec
        write(6,"('Successful read of dipdec from ',a)") flnm(1:lenflnm)
        close(lu)
      endif
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getions(fields,it)
      include 'tigcmproc.h'
      real fields(imx,kmx,jmx,nfget)
c
c Check for tigcm vs tiegcm history:
c
      if (.not.istiegcm.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
        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')") histvol(ivol)(1:lenstr(histvol(ivol)))
        ifplt(ixpot) = 0
        ifplt(ixui) = 0
        ifplt(ixvi) = 0
        ifplt(ixwi) = 0
        ifplt(ixuivi) = 0
        write(6,"(' ')")
        write(6,"('!!!!!!!!!!!!!!!!!!!! NOTICE !!!!!!!!!!!!!!!!!!!')")
        write(6,"(' ')")
      endif 
c
c Get ion drifts if needed (pointers to ui,vi,wi in common in tigcmproc.h):
c
      if (ifplt(ixui).gt.0.or.ifplt(ixuivi).gt.0.or.
     +   (ipuv.gt.0.and.ifplt(ixpot).gt.0)) then
        call alloc(pglbui,imx*kmx*jmx)
        call mkdrifts(ut(it),fields(1,1,1,ixfld(ixpot)),
     +    fields(1,1,1,ixfld(ixz)),gcmlat,glbui,imx,kmx,jmx,1,ionvel,
     +    0,lumag,ludipdec,fields(1,1,1,ixfld(ixu)),
     +    fields(1,1,1,ixfld(ixv)),fields(1,1,1,ixfld(ixw)),
     +    fields(1,1,1,ixfld(ixo2)),fields(1,1,1,ixfld(ixo1)),
     +    fields(1,1,1,ixfld(ixn2)),dum,dum,cpspval)
      endif
      if (ifplt(ixvi).gt.0.or.ifplt(ixuivi).gt.0.or.
     +   (ipuv.gt.0.and.ifplt(ixpot).gt.0)) then
        call alloc(pglbvi,imx*kmx*jmx)
        call mkdrifts(ut(it),fields(1,1,1,ixfld(ixpot)),
     +    fields(1,1,1,ixfld(ixz)),gcmlat,glbvi,imx,kmx,jmx,2,ionvel,
     +    0,lumag,ludipdec,fields(1,1,1,ixfld(ixu)),
     +    fields(1,1,1,ixfld(ixv)),fields(1,1,1,ixfld(ixw)),
     +    fields(1,1,1,ixfld(ixo2)),fields(1,1,1,ixfld(ixo1)),
     +    fields(1,1,1,ixfld(ixn2)),dum,dum,cpspval)
      endif
      if (ifplt(ixwi).gt.0) then
        call alloc(pglbwi,imx*kmx*jmx)
        call mkdrifts(ut(it),fields(1,1,1,ixfld(ixpot)),
     +    fields(1,1,1,ixfld(ixz)),gcmlat,glbwi,imx,kmx,jmx,3,ionvel,
     +    0,lumag,ludipdec,fields(1,1,1,ixfld(ixu)),
     +    fields(1,1,1,ixfld(ixv)),fields(1,1,1,ixfld(ixw)),
     +    fields(1,1,1,ixfld(ixo2)),fields(1,1,1,ixfld(ixo1)),
     +    fields(1,1,1,ixfld(ixn2)),dum,dum,cpspval)
      endif
      return
      end
