c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcm/tigcm.f
c------------------------------------------------------------------
c
      program tigcmp
c
c Plot tigcm fields 
c
      include 'tgcmparam.h'
      include 'input.h'
      include 'tigcmfld.h'
      include 'tiegcmfld.h'
      include 'selgrid.h'
      include 'tigcm.h'
      include 'color.h'
c
      dimension iasf(13)
      logical found
      integer getmss
      data iasf/13*1/
      data runlab/'                                '/
      data ifok/1/, iprnt/0/
c
c Define tgcm grid, set up proc defaults (including conpack),
c   set up labels, and initialize gks:
c
      call setgrid 
      call setdef
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
c
c Case loop:
c
      call lexcon(ncase   ,8HNCASE   ,icase)
      call readlx(5,ier)
      if (icase.le.0) ncase = 1
      do 200 ic=1,ncase 
        write(6,"(' ')")
        write(6,"('-------------------------------------------',
     +            '-----------------------------')")
        write(6,"('tigcm proc: start case ',i3,' of ',i3)") 
     +    ic,ncase
        write(6,"('-------------------------------------------',
     +            '-----------------------------')")
c
c Get user input for this case:
c
        call getinp
        call mklabti(iden)
c
c Define linear height scale:
c (ht1, dht, and nhtscale are defined in tigcm.h, and
c  are in common /linht/ with htscale and finterp)
c
      do 25 ih=1,nhtscale
 25     htscale(ih) = ht1 + (ih-1)*dht
c     write(6,"('tigcm: ht1=',f13.8,' dht=',f13.8,' nhtscale=',i3,
c    +  ' htscale=',/(5f13.8))") ht1,dht,nhtscale,
c    +  (htscale(i),i=1,nhtscale)
c
c Define color table:
        if (iclrfill.gt.0) call defclrs(0)
c
c Get magnetic field data from mss, if this is a dynamo case:
c
        lumag = 99
        if (idyn.gt.0) then
          call mklabtie
          istat = getmss('/FOSTER/magfield',lumag)
          if (istat.ne.ifok) then
            write(6,"('Tigcmf: error ',i3,' from getmss for magfield: ',
     +        'lumag=',i3)") istat,lumag
            stop 'getmss'
          endif
        endif
c
c ut loop:
c
        istart = 1
        mdp = md(1)
        id = 1
        do 100 it=1,ntms
          write(6,"(' ')")
          write(6,"('-------------------------------------------',
     +              '-----------------------------')")
          call gethist(histvol,nhvols,md(it),mh(it),mm(it),
     +      istart,found,ivf)
          istart = 0
          if (.not.found) then
            write(6,"('tigcm proc: history ',i2,':',i2,':',i2,
     +        ' was not found')")
     +        md(it),mh(it),mm(it)
            stop 'hist'
          endif
          if (md(it).gt.mdp) then
            id = id+1
            if (it.eq.ntms) iyd(id) = iyd(id-1)+1 
          endif
          if (found.and.it.eq.1) then
            write(6,"(' ')")
            write(6,"('TIGCM Header:')")
            call prgcmhdr
          endif
c
c Find ion drifts from potential:
c
          if (idyn.gt.0) then
            call mkuivi(poten,pnt(1,1,1,ixz),ui,vi,wi,lumag)
          endif
c
c Write out sample fields:
c
          if (iprnt.gt.0) then
            write(6,"(' ')")
            write(6,"('Loc: samples of history:')")
            do 110 ip=1,ngcmflds 
              do 112 j=1,jmx,6 
                write(6,"('ip=',i2,' j=',i2,' i= 1 k=1->kmx:',/
     +            (5e13.5))") ip,j,(pnt(1,k,j,ip),k=1,kmx)
                write(6,"('ip=',i2,' j=',i2,' i=',i2,' k=1->kmx:',/
     +            (5e13.5))") ip,j,imx,(pnt(imx,k,j,ip),k=1,kmx)
 112          continue
 110        continue
          endif
c
c Make global (cyl equidist) plots at current ut:
          if (iplglb.gt.0) call pltglb(ut(it))
          if (idyn.gt.0.and.iplglb.gt.0) call dynglb(ut(it))
c
c Make polar (stereographic) plots at current ut:
          if (iplpol.gt.0) call pltpol(ut(it))
          if (idyn.gt.0.and.iplpol.gt.0) call dynpol(ut(it))
c
c Make satellite view projection plots at current ut:
          if (iplsatv.gt.0) call pltsatv(ut(it),id)
          if (idyn.gt.0.and.iplsatv.gt.0) call dynsatv(ut(it),id)
c
c Make longitude slices (lat on x-axis, zp on y-axis):
          if (ipllatzp.gt.0) call pltlatzp(ut(it))
c
c Make longitude slices (lat on x-axis, ht on y-axis):
          if (ipllatht.gt.0) call pltlatht(ut(it))
c
c Update lat vs ut scratch files 
c (pltutlat will plot if doing last ut):
          if (iplutlat.gt.0) call pltutlat(ic,it)
c
c Update zp vs ut global mean scratch file 
c (pltutzpm will plot if doing last ut):
          if (iplutzpm.gt.0) call pltutzpm(ic,it)
c
c Update ht vs ut global mean scratch file 
c (pltuthtm will plot if doing last ut):
          if (ipluthtm.gt.0) call pltuthtm(ic,it)
c
c End ut loop:
        mdp = md(it)
 100    continue
c
c End case loop:
 200  continue
c
c Close gks:
      call clsgks
      stop
      end