
      program timesloc
c
c ti(e)gcm location processor:
c
      include 'timesloc.h'
      dimension iasf(13)
      logical isdyn
      data iasf/13*1/
      data luhist/8/, lumag/99/, iframe/0/, zmflag/999./, r12flag/999./
      data gmflag/999./, luascii/90/
      pointer (pfields,fields(1))
      pointer (putkmx,utkmx(1)), (phtutkmx,htutkmx(1)),
     +        (putjmx,utjmx(1)), (putjmxlt,utjmxlt(1)),
     +        (putkmxlt,utkmxlt(1)), (phutkmxlt,hutkmxlt(1))
      logical exists
      dimension cpcit(10),icplit(10) 
      data cpcit /1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 7.5/
      data icplit/  3,   3,   3,   3,   3,   3,   3,   3,   3,   3/
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='timesloc.dat',exist=exists)
        if (exists) istat = ishell('rm timesloc.dat')
      endif
c
c Set up gks, conpack, and autograph:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
      if (labels.gt.1) then
        call cpsetup(cpspval)
      else
        finc = -10. ! number of contour levels in cpcnrc (labels=1)
        call cpsetr('RC2',.5)
        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('LLO -- line label orientation',1)
      endif
      call cpseti('LLP',2)
      call cpsetr('LLS',.017)
      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
      call cpsetr('RC2',.5)
c     call cpseti('LLB',3)
c     if (icolor.gt.0) call cpseti('LBC',1)
      call setag
c
c Allocate space for needed fields and plot arrays:
c
c htutkmx(ntms,kmx,nloc) = heights for extra right hand side axis (pltloc)
c
      if (nloc.gt.0) then
        call alloc(phtutkmx,ntms*kmx*nloc)
        write(6,"('Allocated phtutkmx: ntms=',i2,' kmx=',i2,
     +    ' nloc=',i2)") ntms,kmx,nloc
      endif
c
c hutkmxlt(ntms,kmx,nloclt) = heights for extra right hand side axis (pltloclt)
c
      if (nloclt.gt.0) then
        call alloc(phutkmxlt,ntms*kmx*nloclt)
        write(6,"('Allocated phutkmxlt: ntms=',i2,' kmx=',i2,
     +    ' nloclt=',i2)") ntms,kmx,nloclt
      endif
c
c Ut loop:
c
      istart = 1
      iprhdr = 1
      do it=1,ntms
        call gettgcm(histvol,nhvols,istart,mtimes(1,it),pfields,ifget,
     +    nfhist,imx,kmx,jmx,nflds,iden,-luhist,iprhdr,ivol,isdyn,
     +    istimes,ier)   
        if (ier.ne.0) then
          write(6,"('timesloc: error ',i3,' from gettgcm: it=',i3,
     +      ' mtimes=',i2,':',i2,':',i2)") ier,it,(mtimes(i,it),i=1,3)
          stop 'getimesf'
        endif
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,
     +    labshort,1)
c
c Set up ifget as indices to fields array:
c
        call incind(ifget,nfhist,nfget)	! nfget is returned 
        call incind(ifplt,nftot,nfplt)	! nfplt is returned 
c
c Report to stdout if this is first time:
c
        if (it.eq.1) then
          write(6,"(' ')")
          write(6,"('Need ',i2,' fields from history:')") nfget
          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    n2')") 
          write(6,"(10i6,/10i6,/11i6)") ifget
          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    n2')") 
          write(6,"('    ui    vi    wi e5577 e6300 eo200')")
          write(6,"(10i6,/10i6,/11i6,/6i6)") ifplt 
          write(6,"(' ')")
        endif
c
c utkmx(ntms,kmx,nfplt,nloc) (plot array for locations)
c utjmx(ntms,jmx,nfplt,nzphtlon) (plot array for ut vs lat)
c utjmxlt(ntms,jmx,nfplt,nzphtslt) (for ut vs lat at fixed local time)
c
        if (it.eq.1) then
          if (nloc.gt.0) then
            call alloc(putkmx,ntms*kmx*nfplt*nloc)
            write(6,"('Allocated putkmx: ntms=',i2,' kmx=',i2,
     +        ' nfplt=',i2,' nloc=',i2)") ntms,kmx,nfplt,nloc
          endif
          if (nloclt.gt.0) then
            call alloc(putkmxlt,ntms*kmx*nfplt*nloclt)
            write(6,"('Allocated putkmxlt: ntms=',i2,' kmx=',i2,
     +        ' nfplt=',i2,' nloclt=',i2)") ntms,kmx,nfplt,nloclt
          endif
          if (nzphtlon.gt.0) then
            call alloc(putjmx,ntms*jmx*nfplt*nzphtlon)
            write(6,"('Allocated putjmx: ntms=',i2,' jmx=',i2,
     +        ' nfplt=',i2,' nloc=',i2)") ntms,jmx,nfplt,nloc
          endif
          if (nzphtslt.gt.0) then
            call alloc(putjmxlt,ntms*jmx*nfplt*nzphtslt)
            write(6,"('Allocated putjmxlt: ntms=',i2,' jmx=',i2,
     +        ' nfplt=',i2,' nloc=',i2)") ntms,jmx,nfplt,nloc
          endif
        endif
c
        if (it.eq.1) call setflab(iden,flab,nftot)
        call mkfields(fields,utkmx,htutkmx,utjmx,utjmxlt,it)
c
c Make line plots at current time (field on x-axis, zp/ht on y-axis):
c
        if (nutline.gt.0) then
          do i=1,nutline
            if (mutline(1,i).eq.mtimes(1,it).and.
     +          mutline(2,i).eq.mtimes(2,it).and.
     +          mutline(3,i).eq.mtimes(3,it)) then
              if (nzprange.gt.0) call lineut(utkmx,htutkmx,it,0)
              if (nhtscale.gt.0) call lineut(utkmx,htutkmx,it,1)
            endif
          enddo
        endif
c
c End ut loop:
        istart = 0
        iprhdr = 0
        call hpdeallc(pfields,ier,1)	! was alloc'd by gettgcm
      enddo
c
c Contour ut on x-axis, zp/ht on y-axis:
c
      if (nloc.gt.0) then
        do l=1,nloc      
          write(6,"(' ')")
          call agseti('AXIS/LEFT/CONTROL.',1)
          call agsetc ('LABEL/NAME.','L')
          call agseti ('LABEL/CONTROL.',1)
          if (iutvert.gt.0) then
            if (nzprange.gt.0) call pltloc(utkmx,htutkmx,0,l)
            if (nhtscale.gt.0) call pltloc(utkmx,htutkmx,1,l)
          endif
c
c Make line plots with ut on x-axis at selected zp and/or ht
c
          call setag
          call agseti('AXIS/LEFT/CONTROL.',1)
          call agseti('AXIS/BOTTOM/CONTROL.',1)
          if (nzphtlin.gt.0) then
            call agsetc ('LABEL/NAME.','L')
            call agseti ('LABEL/CONTROL.',2)
            call linezpht(utkmx,htutkmx,l)
          endif
c
c Add doppler plots from e5577 emission:
c
          if (idop5577.gt.0.and.ifplt(ixe5577).gt.0.and.
     +        nhtscale.gt.0) 
     +      call pltdop(dop5577,l,'E5577')
c
c Add doppler plots from e6300 emission:
c
          if (idop6300.gt.0.and.ifplt(ixe6300).gt.0.and.
     +        nhtscale.gt.0) 
     +      call pltdop(dop6300,l,'E6300')
c
c Add doppler plots from eo200 emission:
c
          if (idopo200.gt.0.and.ifplt(ixeo200).gt.0.and.
     +        nhtscale.gt.0) 
     +      call pltdop(dopo200,l,'EO200')
        enddo 
      endif
c
c zp/ht vs ut at selected latitudes and fixed local times:
c
      if (nloclt.gt.0.and.iutvert.gt.0) then
        do l=1,nloclt
          if (nzprange.gt.0) call pltloclt(utkmxlt,hutkmx,0,l)
          if (nhtscale.gt.0) call pltloclt(utkmxlt,hutkmx,1,l)
        enddo
      endif
c
c Make line plots with field on x-axis, zp/ht on y-axis, with
c each ut line on same frame for each field:
c
      if (nutlines.gt.0) then
        if (nzprange.gt.0) call lineuts(utkmx,htutkmx,0)
        if (nhtscale.gt.0) call lineuts(utkmx,htutkmx,1)
      endif
c
c Contour ut on x-axis, lat on y-axis at selected longitudes:
c
      if (nzphtlon.gt.0) then
        do l=1,nzphtlon
          write(6,"(' ')")
          call pltutlat(utjmx,l)
        enddo
      endif
c
c Contour ut on x-axis, lat on y-axis at fixed local times:
c
      if (nzphtslt.gt.0) then
        do l=1,nzphtslt
          write(6,"(' ')")
          call putlatlt(utjmxlt,l)
        enddo
      endif
c
      if (nloc.gt.0) call hpdeallc(putkmx,ier,1)
      if (nzphtlon.gt.0) call hpdeallc(putjmx,ier,1)
      if (nzphtslt.gt.0) call hpdeallc(putjmxlt,ier,1)
      call clsgks
      if (iwrascii.eq.2.and.lenstr(dirascii).gt.0)
     +  call rcpfile(luascii,'timesloc.dat',dirascii)
      if (lenstr(dircgm).ne.0) call rcpfile(0,'gmeta',dircgm)
      stop 'success!'
      end
c
