
      program tigcmloc
c
c ti(e)gcm location processor:
c
      include 'tigcmloc.h'
      include 'tgcmhdr.h'
      dimension iasf(13)
      data iasf/13*1/
      pointer (pfields,fields(1))
      pointer (putkmx,utkmx(1)), (phtutkmx,htutkmx(1)),
     +        (putjmx,utjmx(1))
      logical exists
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='tigcmloc.dat',exist=exists)
        if (exists) istat = ishell('rm tigcmloc.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',1)
      if (icolor.gt.0) then
        call cpseti('LBC',1)
        call cpseti('LLB',3)
      endif
      call setag
c
c Allocate space for needed fields and plot arrays:
c
c fields(imx,kmx,jmx,nfget) (will be deallocated after ut loop)
c utkmx(ntms,kmx,nfplt,nloc) (plot array for locations)
c htutkmx(ntms,kmx,nloc) = heights for extra right hand side axis (pltloc)
c utjmx(ntms,jmx,nfplt,nzphtlon) (plot array for ut vs lat)
c
c
      call alloc(pfields,imx*kmx*jmx*nfget)
      if (nloc.gt.0) then
        call alloc(putkmx,ntms*kmx*nfplt*nloc)
        call alloc(phtutkmx,ntms*kmx*nloc)
      endif
      if (nzphtlon.gt.0) call alloc(putjmx,ntms*jmx*nfplt*nzphtlon)
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,"('tigcmloc: error ',i3,' from getgcmf: it=',i3,
     +      ' mtimes=',i2,':',i2,':',i2)") ier,it,(mtimes(i,it),i=1,3)
          stop 'getgcmf'
        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),
     +      histvol(ivol))
        if (it.eq.1) then
          call setflab(istiegcm,iden,flab,nftot,ionvel)
          if (ionvel.eq.2.and.(ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.
     +                         ifplt(ixwi).gt.0)) then
            call alloc(pveutkmx,ntms*kmx*nloc)
            call alloc(pvuutkmx,ntms*kmx*nloc)
          endif 
        endif
c
c Get doppler fields if necessary:
c
        if (idop.gt.0) call mkdop(fields,it)
c
c Define plot arrays:
c
        call mkfields(fields,utkmx,htutkmx,utjmx,it)
c
c Make line plots at current time (field on x-axis, zp/ht on y-axis):
c
        if (nutline.gt.0) then
          nzprange = izprange(2)-izprange(1)
          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(fields,utkmx,htutkmx,it,0)
              if (nhtscale.gt.0) call lineut(fields,utkmx,htutkmx,it,1)
            endif
          enddo
        endif
c
c End ut loop:
        istart = 0
        iprhdr = 0
      enddo
      call hpdeallc(pfields,ier,1)
c
c Contour ut on x-axis, zp/ht on y-axis:
c
      if (nloc.gt.0) then
        nzprange = izprange(2)-izprange(1)
        do l=1,nloc      
          write(6,"(' ')")
          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) call linezpht(utkmx,htutkmx,l)
          if (idop.gt.0) call pltdop(l)
        enddo 
      endif
c
c Contour ut on x-axis, lat on y-axis:
c
      if (nzphtlon.gt.0) then
        do l=1,nzphtlon
          write(6,"(' ')")
          call pltutlat(utjmx,l)
        enddo
      endif
c
      if (nloc.gt.0) then
        call hpdeallc(putkmx,ier,1)
        call hpdeallc(phtutkmx,ier,1)
      endif
      if (ionvel.eq.2.and.(ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.
     +                     ifplt(ixwi).gt.0)) then
        call hpdeallc(pveutkmx,ier,1)
        call hpdeallc(pvuutkmx,ier,1)
      endif 
      call clsgks
      if (iwrascii.eq.2.and.lenstr(dirascii).gt.0)
     +  call rcpfile(luascii,'tigcmloc.dat',dirascii)
      if (lenstr(dircgm).ne.0) call rcpfile(0,'gmeta',dircgm)
      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 mkdop(fields,it)
c
c Define doppler fields for current ut:
c
      include 'tigcmloc.h'
      include 'tgcmhdr.h'
      include 'doppler.h'
      dimension fields(imx,kmx,jmx,nfget)
      dimension ht(kmx),tn(kmx),un(kmx),vn(kmx),te(kmx),xo2(kmx),
     +          xo1(kmx),xn2(kmx),xo2p(kmx),xne(kmx),verk(kmx)
c
c
      iyd = (date(1)-1900.)*1000+date(2) ! date(2) in tgcmhdr
      if (nloc.gt.0) then
        do l=1,nloc      
          if (gloc(1,l).eq.gmflag.or.gloc(2,l).eq.gmflag) then
            write(6,"('>>> mkdop: doppler fields not available ',
     +        'in zonal or global means')")
            doptn(it,l) = 0. 
            dopun(it,l) = 0. 
            dopvn(it,l) = 0. 
            ver(it,:,l) = 0. 
            ray(it,l) = 0.
          else
c
c     subroutine defutkmx(field,utkmx,glat,glon,it)
c
            call defutkmx(fields(1,1,1,ifget(ixt)),tn,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixu)),un,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixv)),vn,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixo2)),xo2,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixo1)),xo1,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixn2)),xn2,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixte)),te,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixz)),ht,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixo2p)),xo2p,gloc(1,l),
     +        gloc(2,l),it)
            call defutkmx(fields(1,1,1,ifget(ixne)),xne,gloc(1,l),
     +        gloc(2,l),it)
          endif
c
c     subroutine doptuv(luout,utc,iyd,f107d,f107a,glat,glon,kmx,
c    +  ht,tn,un,vn,te,xo2,xo1,xn2,xo2p,xne,doptn,dopun,dopvn,
c    +  ray,ver)
c
c         write(6,"('mkdop calling doptuv: it=',i2,' l=',i2,
c    +      ' ut(it)=',f5.2,' iyd=',i6,' f107d,a=',2f7.2,
c    +      ' glat,lon=',f6.2,' ',f7.2)")
c    +      it,l,ut(it),iyd,f107d,f107a,gloc(1,l),gloc(2,l)
          call doptuv(ludop,ut(it),iyd,f107d,f107a,gloc(1,l),
     +      gloc(2,l),kmx,ht,tn,un,vn,te,xo2,xo1,xn2,xo2p,xne,
     +      doptn(it,l),dopun(it,l),dopvn(it,l),ray(it,l),verk)
          ver(it,:,l) = verk(:)
c         write(6,"('mkdop after doptuv: ver=',/(6e12.4))") verk
        enddo ! l=1,nloc
      else
        write(6,"('>>> tigcmloc: nloc should be > 0 if doing ',
     +    'doppler fields')")
      endif
      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
      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
