c
c------------------------------------------------------------------
c
      program timespro
c
c Timesgcm snapshot processor:
c
      include 'timesproc.h'
      dimension iasf(13)
      data iasf/13*1/
      logical exists,issech
      pointer (pfields,fields(1))
      dimension cpcit(10),icplit(10),plt(imx,jmx) 
      character*56 toplab,proclab
      data cpcit /1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 7.5/
      data icplit/  2,   2,   2,   2,   2,   2,   2,   2,   2,   2/
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)
      do k=1,kmx
        pmb(k) = p0 * exp(-gcmzp(k))	! pressure (mb)
      enddo
      call getinp
      call setag
      call agseti('AXIS/LEFT/CONTROL.',1)
      call agseti('AXIS/BOTTOM/CONTROL.',1)
      call pcsetc('FC','$')
c
c If labels=1, will do only simple labelling and contouring:
c
      if (labels.gt.1) then
        call cpsetup(cpspval)
        if (icolor.gt.0) then
          call cpseti('LLP',2)
          call cpseti('LLB',0)
c         call cpseti('LLB',3)
c         call cpseti('LBC',1)
c         call cpsetr('LLS',.012)
          call cpsetr('LLS',.016)
          call cpseti('LIS',2)
        else
          call cpseti('LLP',3)
c         call cpsetr('LLS',.021) ! temporary
c         call cpsetr('PC6',.70) ! temporary
        endif
        call cpsetr('PC6',.7)
      else
        call pcseti('QU',2) ! use low quality chars
        call cpsetr('SPV',cpspval)
        call cpsetc('ILT','MIN $ZMN$, MAX $ZMX$, INT $CIU$')
        call cpsetr('ORV',1.e12)
        call cpseti('LLP',2)
        call cpsetr('LLS',.017)
        call cpsetr('RC2',.5)
        finc = -10. ! number of contour levels in cpcnrc (labels=1)
        call cpseti('LLO -- line label orientation',1)
        do i=1,10
          call cpseti('PAI -- parameter array index',i)
          call cpsetr('CIT -- contour interval table',cpcit(i))
          call cpseti('LIT -- label interval table',icplit(i))
        enddo
      endif
c
c Delete and ascii data files left over from previous runs:
c
      if (iwrascii.gt.0) then
        inquire(file='timesproc.dat',exist=exists)
        if (exists) istat = unlink('timesproc.dat')
      endif
c
c Set up gks:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
c
c Set up linear height scale, if needed:
c
      nzprange = izprange(2)-izprange(1)
      nhtscale = 0
      if (delht.gt.0.) then
        nhtscale = ifix((httop-htbot) / delht + 1.0000001)
        call alloc(phtscale,nhtscale)
        do i=1,nhtscale
          htscale(i) = htbot + (i-1)*delht
        enddo
        if (ipllon.gt.0.or.ipllat.gt.0.or.iplloc.gt.0.or.iplglbm.gt.0)
     +  write(6,"('timesproc: nhtscale=',i3,' linear height scale = ',
     +    /(6e10.3))") nhtscale,(htscale(i),i=1,nhtscale)
      endif
c
c Loop through requested times, reading histories and plotting:
c 10/94: gettgcm now accepts timesgcm secondary histories, which
c        contain an arbitrary number of fields (max 29)
c
      istart = 1
      iprhdr = 1
      do it=1,ntms
        write(6,"(' ')")
        write(6,"('Read timesgcm hist: it=',i2,' ntms=',i2,
     +    ' md:mh:mm=',i3,':',i2,':',i2)") 
     +    it,ntms,(mtimes(i,it),i=1,3)
c       call gettgcm(histvol,nhvols,istart,mtimes(1,it),pfields,
c    +    ifget,nfhist,imx,kmx,jmx,nflds,iden,luhist,iprhdr,ivol,
c    +    isdyn,istimes,ier)
c
c     subroutine rdtgcm(vol,nvol,tempdir,istart,mt,pfields,ifields,
c    +  mxfields,nlon,nzp,nlat,nflds,iden,ionvel,lu,iprhdr,ivol,isdyn,
c    +  istimes,issech,ier)
c
        call rdtgcm(histvol,nhvols,'/usr/tmp/TIGCM',istart,
     +    mtimes(1,it),pfields,ifget,nfhist,imx,nzp,jmx,nflds,iden,
     +    ionvel,luhist,iprhdr,ivol,isdyn,istimes,issech,ier)
        if (nzp.ne.kmx) then
          write(6,"('>>> rdtgcm returned nzp=',i2,' .ne. kmx=',i2)")
     +      nzp,kmx
          write(6,"('Use this proc for timesgcm kmx=42 only')")
          stop 'kmx'
        endif

        if (ier.eq.1) then
          write(6,"('>>> Stopping because proc could not find ',
     +      'model time ',i3,':',i2,':',i2)") (mtimes(i,it),i=1,3)
          write(6,"('    History volume(s) searched = ')")
          do j=1,nhvols
            write(6,"(a)") histvol(j)
          enddo
          stop 'histtime'
        elseif (ier.ne.0) then
          write(6,"('>>> Error ',i3,' from rdtgcm')") ier
          stop 'rdtgcm'
        endif
        istart = 0
        iprhdr = 0
c
c Turn off any unavailable fields, and any requested diagnostic fields 
c that depend on history fields that are unavailable:
c
        call chflds(ifget,nfhist,ifplt,nftot,ifdep,nftot-nfhist,0,
     +    flab,1)
c
c Cannot plot vectors if winds/drifts not available:
c
        if (ituv.gt.0.and.ifplt(itxt).gt.0.and.
     +     (ifget(itxu).le.0.or.ifget(itxv).le.0)) then
          write(6,"('>>> Warning: un & vn not available',
     +      ' (cannot plot neutral winds over tn)')")
          ituv = 0
        endif
        if (izuv.gt.0.and.ifplt(itxz).gt.0.and.
     +     (ifget(itxu).le.0.or.ifget(itxv).le.0)) then
          write(6,"('>>> Warning: un & vn not available',
     +      ' (cannot plot neutral winds over height)')")
          izuv = 0
        endif
        if (ipuv.gt.0.and.ifplt(itxpot).gt.0.and.
     +     (ifget(itxui).le.0.or.ifget(itxvi).le.0)) then
          write(6,"('>>> Warning: ui & vi not available',
     +      ' (cannot plot ion drifts over e-poten)')")
          ipuv = 0
        endif
c
c Set up ifget as indices to fields array:
c
        call incind(ifget,nfhist,nfget)	! nfget is returned 
        nfplt = 0
        do i=1,nftot
          if (ifplt(i).gt.0) nfplt = nfplt+1
        enddo
c
c Report to stdout if this is first time:
c
        if (it.eq.1) then
          write(6,"(' ')")
          write(6,"('Will plot ',i2,' fields:')") nfplt
          write(6,"('    tn    un    vn    o2    ox   n4s',
     +      '   noz    co   co2   h2o',/'    h2   hox    o+   ch4',
     +      '  o21d   no2    no    o3    o1    oh',/,'   ho2     h',
     +      '   n2d    ti    te    ne   o2+     w     z  epot')")
          write(6,"('    ui    vi    wi    n2  o/o2  o/n2 o2/n2   rho',
     +      '  htot e5577 e6300 eo200')")
          write(6,"(3(10i6/),12i6))") ifplt 
          write(6,"(' ')")
        endif
c
c Get ion drifts: 
c
c       call getions(fields,it)
c
c Ascii files of each field at global grid:
c
        if (iwrgrid.gt.0) then
          ifr = 0
          do ip=1,nftot
            if (ifplt(ip).gt.0) then
              call clearstr(proclab)
              write(proclab,"('timesproc.',a,'.',3i2.2)")
     +          flab(ip)(1:lenstr(flab(ip))),
     +          (mtimes(i,it),i=1,3)
              lenproc = lenstr(proclab)
              iwr = 2
              do k=1,kmx
                call getglb(fields,plt,k,0.,ip)
                write(toplab,"('FIELD ',a,' ZP=',f5.1,' MTIME=',3i3)")
     +            flab(ip),gcmzp(k),(mtimes(i,it),i=1,3)
                ifr = ifr+1
                call wrascii(iwr,luascii,plt,imx,jmx,gcmlon,gcmlat,
     +           'LONGITUDE','LATITUDE',histvol(ivol),flab(ip),
     +            ifr,toplab(1:lenstr(toplab)),
     +            proclab(1:lenproc),dirascii)
              enddo
              proclab = proclab(1:lenproc)//'.dat'
              lenproc = lenstr(proclab)
              if (lenstr(dirascii).gt.0) then
                call rcpfile(luascii,proclab(1:lenproc),dirascii)
              else
                close(luascii)
              endif
            endif
          enddo
        endif
c
c Global cylindrical equidistant plots:
c
        if (iplglb.gt.0) call pltglb(fields,mtimes(1,it))
        if (ipsltglb.gt.0) call sltglb(fields,ut(it))
c
c Polar stereographic plots:
c
        if (npollat.gt.0) call pltpol(fields,ut(it))
c
c Satellite view projection:
c
        if (nsatv.gt.0) call pltsatv(fields,ut(it))
c
c Longitude slices (lat on x-axis, pressure and/or height on y-axis):
c
        if (ipllon.gt.0) then
          istream = 0
          nz = 0
          do i=1,nlon
            if (slon(i).eq.zmflag) nz = 1 
          enddo
          if (ifplt(itxv).gt.0.and.nz.gt.0) istream = 1
          if (nzprange.gt.0) then
            call pltlon(fields,ut(it),0)
            if (istream.gt.0) call pltstrm(fields,ut(it),0)
          endif
          if (nhtscale.gt.0) then
            call pltlon(fields,ut(it),1)
            if (istream.gt.0) call pltstrm(fields,ut(it),1)
          endif
        endif
c
c Latitude slices (lon on x-axis, pressure and/or height on y-axis):
c
        if (ipllat.gt.0) then
          if (nzprange.gt.0) call pltlat(fields,ut(it),0)
          if (nhtscale.gt.0) call pltlat(fields,ut(it),1)
        endif
c
c Calculate and plot amplitude and phases (fourier transform)
c
        if (iampha.gt.0) then
          if (nzprange.gt.0) call pltampha(fields,ut(it),0)
          if (nhtscale.gt.0) call pltampha(fields,ut(it),1)
        endif
c
c Make line plots of zp/ht profiles at selected locations:
c
        if (iplloc.gt.0.or.iplglbm.gt.0) call setag
        if (iplloc.gt.0) then
          if (nzprange.gt.0) call pltloc(fields,ut(it),0,0)
          if (nhtscale.gt.0) call pltloc(fields,ut(it),0,1)
        endif
c
c Make line plots of zp/ht profiles at global means:
c
        if (iplglbm.gt.0) then
          if (nzprange.gt.0) call pltloc(fields,ut(it),1,0)
          if (nhtscale.gt.0) call pltloc(fields,ut(it),1,1)
        endif
c
c Release space that was allocated off the heap for this model time:
c (pfields was allocated by gettgcm)
c
        call hpdeallc(pfields,ier,1)
        if (ifplt(itxui).gt.0.or.(ipuv.gt.0.and.ifplt(itxpot).gt.0))
     +    call hpdeallc(pglbui,ier,1)
        if (ifplt(itxvi).gt.0.or.(ipuv.gt.0.and.ifplt(itxpot).gt.0))
     +    call hpdeallc(pglbvi,ier,1)
        if (ifplt(itxwi).gt.0) call hpdeallc(pglbwi,ier,1)
c
c End ut loop:
      enddo
      if (delht.gt.0.) then
        call hpdeallc(phtscale,ier,1)
        if (ier.ne.0) write(6,"('>>> warning: hpdeallc error = ',i6,
     +      ' for phtscale')") ier
      endif
c
      call clsgks
      if (iwrascii.eq.2.and.lenstr(dirascii).gt.0)
     +  call rcpfile(luascii,'timesproc.dat',dirascii)
      if (lenstr(dircgm).ne.0) call rcpfile(0,'gmeta',dircgm)
      stop 'done'
      end
c
c------------------------------------------------------------------
c
      subroutine getions(fields,it)
      include 'timesproc.h'
      real fields(imx,kmx,jmx,nfget)
c
c Get ion drifts if needed (pointers to ui,vi,wi in common in timesproc.h):
c
      if (ifplt(itxui).gt.0.or.(ipuv.gt.0.and.ifplt(itxpot).gt.0))then
        call alloc(pglbui,imx*kmx*jmx)
        call mkdrifts(ut(it),fields(1,1,1,ifget(itxpot)),
     +    fields(1,1,1,ifget(itxz)),gcmlat,glbui,imx,kmx,jmx,1,ionvel,
     +    0,lumag,ludipdec,fields(1,1,1,ifget(itxu)),
     +    fields(1,1,1,ifget(itxv)),fields(1,1,1,ifget(itxw)),
     +    fields(1,1,1,ifget(itxo2)),fields(1,1,1,ifget(itxo1)),
     +    fields(1,1,1,ifget(itxn2)),dum,dum,cpspval)
      endif
      if (ifplt(itxvi).gt.0.or.(ipuv.gt.0.and.ifplt(itxpot).gt.0))then
        call alloc(pglbvi,imx*kmx*jmx)
        call mkdrifts(ut(it),fields(1,1,1,ifget(itxpot)),
     +    fields(1,1,1,ifget(itxz)),gcmlat,glbvi,imx,kmx,jmx,2,ionvel,
     +    0,lumag,ludipdec,fields(1,1,1,ifget(itxu)),
     +    fields(1,1,1,ifget(itxv)),fields(1,1,1,ifget(itxw)),
     +    fields(1,1,1,ifget(itxo2)),fields(1,1,1,ifget(itxo1)),
     +    fields(1,1,1,ifget(itxn2)),dum,dum,cpspval)
      endif
      if (ifplt(itxwi).gt.0) then
        call alloc(pglbwi,imx*kmx*jmx)
        call mkdrifts(ut(it),fields(1,1,1,ifget(itxpot)),
     +    fields(1,1,1,ifget(itxz)),gcmlat,glbwi,imx,kmx,jmx,3,ionvel,
     +    0,lumag,ludipdec,fields(1,1,1,ifget(itxu)),
     +    fields(1,1,1,ifget(itxv)),fields(1,1,1,ifget(itxw)),
     +    fields(1,1,1,ifget(itxo2)),fields(1,1,1,ifget(itxo1)),
     +    fields(1,1,1,ifget(itxn2)),dum,dum,cpspval)
      endif
      return
      end