
      subroutine getglb(fields,plt,ixpl,ht,ip)
c
c Return desired global field (at ip) in plt(imx,jmx).
c ixpl = index to desired pressure. If ixpl=0, assume interpolation to
c        constant height surface ht
c
      include 'tigcmproc.h'
      dimension plt(imx,jmx), fields(imx,kmx,jmx,nfget),plt0(imx,jmx),
     +  plt1(imx,kmx,jmx)
      dimension xne(kmx),hts(kmx)
c
      if (ixpl.le.0.and.ixfld(ixz).le.0) then
        write(6,"('>>> getglb: need heights for interpolation: ',
     +    ' ixz=',i3,' ifget(ixz)=',i3,' ixfld(ixz)=',i3)") 
     +    ixz,ifget(ixz),ixfld(ixz)
      endif
c
c Field from history:
c
      if (ip.le.0) then
        write(6,"('>>> getglb: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      elseif (ip.le.nfhist) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ixfld(ip))
        else
          call glbhtint(fields(1,1,1,ixfld(ip)),
     +                  fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt,ht,1,ilog(ip),cpspval,ier,0)
        endif
c
c ion drifts:
      elseif (ip.ge.ixui.and.ip.le.ixwi) then
        if (ixpl.gt.0) then
          if (ip.eq.ixui) plt(:,:) = glbui(:,ixpl,:)
          if (ip.eq.ixvi) plt(:,:) = glbvi(:,ixpl,:)
          if (ip.eq.ixwi) plt(:,:) = glbwi(:,ixpl,:)
        else
          if (ip.eq.ixui) call glbhtint(glbui,fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt,ht,1,ilog(ip),cpspval,ier,0)
          if (ip.eq.ixvi) call glbhtint(glbvi,fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt,ht,1,ilog(ip),cpspval,ier,0)
          if (ip.eq.ixwi) call glbhtint(glbwi,fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt,ht,1,ilog(ip),cpspval,ier,0)
        endif
c
c fof2:
      elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
        do i=1,imx
          do j=1,jmx
            xne(:) = fields(i,:,j,ixfld(ixne))
            hts(:) = fields(i,:,j,ixfld(ixz))
            call fof2int(xne,hts,kmx,hmf2out,fof2out,0,i,j)
            if (ip.eq.ixfof2) plt(i,j) = fof2out
            if (ip.eq.ixhmf2) plt(i,j) = hmf2out
          enddo
        enddo
c
c rho (o2+o+n2):
      elseif (ip.eq.ixrho) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ixfld(ixo1)) +
     +               fields(:,ixpl,:,ixfld(ixo2)) + 
     +               fields(:,ixpl,:,ixfld(ixn2))
        else
          call glbhtint(fields(1,1,1,ixfld(ixo1)),
     +      fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,plt,ht,1,ilog(ip),
     +      cpspval,ier,0)
          call glbhtint(fields(1,1,1,ixfld(ixo2)),
     +      fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,plt0,ht,1,ilog(ip),
     +      cpspval,ier,0)
          do j=1,jmx
            do i=1,imx
              if (plt0(i,j).ne.cpspval.and.plt(i,j).ne.cpspval) then
                plt(i,j) = plt(i,j) + plt0(i,j)
              else
                plt(i,j) = cpspval
              endif
            enddo
          enddo
          call glbhtint(fields(1,1,1,ixfld(ixn2)),
     +      fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,plt0,ht,1,ilog(ip),
     +      cpspval,ier,0)
          do j=1,jmx
            do i=1,imx
              if (plt0(i,j).ne.cpspval.and.plt(i,j).ne.cpspval) then
                plt(i,j) = plt(i,j) + plt0(i,j)
              else
                plt(i,j) = cpspval
              endif
            enddo
          enddo
        endif
c
c o/o2:
      elseif (ip.eq.ixroo2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ixfld(ixo1)) /
     +               fields(:,ixpl,:,ixfld(ixo2))
        else
          plt1(:,:,:) = fields(:,:,:,ixfld(ixo1)) /
     +                  fields(:,:,:,ixfld(ixo2))
          call glbhtint(plt1,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +      plt,ht,1,ilog(ip),cpspval,ier,0)
        endif 
c
c o/n2:
      elseif (ip.eq.ixron2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ixfld(ixo1)) /
     +               fields(:,ixpl,:,ixfld(ixn2))
        else
          plt1(:,:,:) = fields(:,:,:,ixfld(ixo1)) /
     +                  fields(:,:,:,ixfld(ixn2))
          call glbhtint(plt1,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +      plt,ht,1,ilog(ip),cpspval,ier,0)
        endif 
c
c o2/n2:
      elseif (ip.eq.ixro2n2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ixfld(ixo2)) /
     +               fields(:,ixpl,:,ixfld(ixn2))
        else
          plt1(:,:,:) = fields(:,:,:,ixfld(ixo2)) /
     +                  fields(:,:,:,ixfld(ixn2))
          call glbhtint(plt1,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +      plt,ht,1,ilog(ip),cpspval,ier,0)
        endif 
c
c o/(o2+n2):
      elseif (ip.eq.ixroo2n2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ixfld(ixo1)) /
     +      (fields(:,ixpl,:,ixfld(ixo2))+fields(:,ixpl,:,ixfld(ixn2)))
        else
          plt1(:,:,:) = fields(:,:,:,ixfld(ixo1)) /
     +      (fields(:,:,:,ixfld(ixo2)) + fields(:,:,:,ixfld(ixn2)))
          call glbhtint(plt1,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +      plt,ht,1,ilog(ip),cpspval,ier,0)
        endif 
c
c Bad ip:
      else
        write(6,"('>>> getglb: unsupported ip=',i3)") ip
      endif
c
c Take log10 if desired:
      if (log10map.gt.0.and.ilog(ip).gt.0) 
     +  call log10f(plt,imx*jmx,1.e-20,cpspval)
c
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getlon(fields,plt,ixlon,log,ip)
      include 'tigcmproc.h'
      dimension plt(jmx,kmx),rimx(imx),fields(imx,kmx,jmx,nfget)
c
      if (ip.le.0) then
        write(6,"('>>> getlon: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      elseif (ip.le.nfhist) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ixfld(ip))
            enddo
          enddo
          if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ixfld(ip))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0)
            enddo
          enddo
        endif
c
c ion drifts:
      elseif (ip.ge.ixui.and.ip.le.ixwi) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              if (ip.eq.ixui) plt(j,k) = glbui(ixlon,k,j) 
              if (ip.eq.ixvi) plt(j,k) = glbvi(ixlon,k,j) 
              if (ip.eq.ixwi) plt(j,k) = glbwi(ixlon,k,j) 
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              if (ip.eq.ixui) rimx(:) = glbui(:,k,j)
              if (ip.eq.ixvi) rimx(:) = glbvi(:,k,j)
              if (ip.eq.ixwi) rimx(:) = glbwi(:,k,j)
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0)
            enddo
          enddo
        endif
c
c fof2 and hmf2 are height independent, so not returned here:
c
      elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
        write(6,"('>>> getlon: ip=',i3,' ixfof2 and ixhmf2 are ',
     +    'height independent -- returning zeroes')")
        plt(:,:) = 0.
        return
c
c rho (o2+o+n2):
c
      elseif (ip.eq.ixrho) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ixfld(ixo1)) +
     +                   fields(ixlon,k,j,ixfld(ixo2)) +
     +                   fields(ixlon,k,j,ixfld(ixn2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ixfld(ixo1))
              plt(j,k) = calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ixfld(ixo2))
              plt(j,k) = plt(j,k) + calcmean(rimx,imx-1,0,1.e-20,
     +          cpspval,0)
              rimx(:) = fields(:,k,j,ixfld(ixn2))
              plt(j,k) = plt(j,k) + calcmean(rimx,imx-1,0,1.e-20,
     +          cpspval,0)
            enddo
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
c
c o/o2:
      elseif (ip.eq.ixroo2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ixfld(ixo1)) /
     +                   fields(ixlon,k,j,ixfld(ixo2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ixfld(ixo1)) /
     +                  fields(:,k,j,ixfld(ixo2))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0) 
            enddo
          enddo
        endif
c
c o/n2:
      elseif (ip.eq.ixron2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ixfld(ixo1)) /
     +                   fields(ixlon,k,j,ixfld(ixn2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ixfld(ixo1)) /
     +                  fields(:,k,j,ixfld(ixn2))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0) 
            enddo
          enddo
        endif
c
c o2/n2:
      elseif (ip.eq.ixro2n2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ixfld(ixo2)) /
     +                   fields(ixlon,k,j,ixfld(ixn2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ixfld(ixo2)) /
     +                  fields(:,k,j,ixfld(ixn2))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0) 
            enddo
          enddo
        endif
      elseif (ip.eq.ixroo2n2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ixfld(ixo1)) /
     +          (fields(ixlon,k,j,ixfld(ixo2)) + 
     +           fields(ixlon,k,j,ixfld(ixn2)))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ixfld(ixo1)) /
     +          (fields(:,k,j,ixfld(ixo2)) +
     +           fields(:,k,j,ixfld(ixn2)))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0) 
            enddo
          enddo
        endif

c
c Unknown field:
      else
        write(6,"('>>> getlon: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getlat(fields,plt,ixlat,log,ip)
      include 'tigcmproc.h'
      dimension plt(imx,kmx),fields(imx,kmx,jmx,nfget)
c
      if (ip.le.0) then
        write(6,"('>>> getlat: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      elseif (ip.le.nfhist) then
        plt(:,:) = fields(:,:,ixlat,ixfld(ip))
      elseif (ip.ge.ixui.and.ip.le.ixwi) then
        if (ip.eq.ixui) plt(:,:) = glbui(:,:,ixlat)
        if (ip.eq.ixvi) plt(:,:) = glbvi(:,:,ixlat)
        if (ip.eq.ixwi) plt(:,:) = glbwi(:,:,ixlat)
c
c fof2 and hmf2 are height independent, so not returned here:
c
      elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
        write(6,"('>>> getlat: ip=',i3,' ixfof2 and ixhmf2 are ',
     +    'height independent -- returning zeroes')")
        plt(:,:) = 0.
        return
c
c rho (o2+o+n2):
c
      elseif (ip.eq.ixrho) then
        plt(:,:) = fields(:,:,ixlat,ixfld(ixo1)) +
     +             fields(:,:,ixlat,ixfld(ixo2)) +
     +             fields(:,:,ixlat,ixfld(ixn2))
c
c Ratios o/o2, o/n2, o2/n2, and o/(o2+n2):
      elseif (ip.eq.ixroo2) then
        plt(:,:) = fields(:,:,ixlat,ixfld(ixo1)) /
     +             fields(:,:,ixlat,ixfld(ixo2))
      elseif (ip.eq.ixron2) then
        plt(:,:) = fields(:,:,ixlat,ixfld(ixo1)) /
     +             fields(:,:,ixlat,ixfld(ixn2))
      elseif (ip.eq.ixro2n2) then
        plt(:,:) = fields(:,:,ixlat,ixfld(ixo2)) /
     +             fields(:,:,ixlat,ixfld(ixn2))
      elseif (ip.eq.ixroo2n2) then
        plt(:,:) = fields(:,:,ixlat,ixfld(ixo1)) /
     +    (fields(:,:,ixlat,ixfld(ixo2))+fields(:,:,ixlat,ixfld(ixn2)))
c
c Unknown field:
      else
        write(6,"('>>> getlat: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      endif
      if (log.gt.0) call log10f(plt,imx*kmx,1.e-20,cpspval)
      return 
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getloc(fields,plt,ixlat,ixlon,log,ip)
      include 'tigcmproc.h'
      dimension plt(kmx),fields(imx,kmx,jmx,nfget),glb(imx,jmx),
     +  rimx(imx)
c
c fof2 and hmf2 are height independent, so not returned here:
c
      if (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
        write(6,"('>>> getloc: ip=',i3,' ixfof2 and ixhmf2 are ',
     +    'height independent -- returning zeroes')") ip
        plt(:) = 0.
        return
      elseif (ip.le.0.or.ip.gt.nftot) then
        write(6,"('>>> getloc: unsupported ip=',i3,' -- returning ',
     +    'zeroes')") ip
        plt(:) = 0.
        return
      endif
c
c At specified location:
c
      if (ixlat.gt.0.and.ixlon.gt.0) then
        if (ip.le.nfhist) then
          plt(:) = fields(ixlon,:,ixlat,ixfld(ip))
        elseif (ip.ge.ixui.and.ip.le.ixwi) then
          if (ip.eq.ixui) plt(:) = glbui(ixlon,:,ixlat)
          if (ip.eq.ixvi) plt(:) = glbvi(ixlon,:,ixlat)
          if (ip.eq.ixwi) plt(:) = glbwi(ixlon,:,ixlat)
        elseif (ip.eq.ixrho) then
          plt(:) = fields(ixlon,:,ixlat,ixfld(ixo1)) +
     +             fields(ixlon,:,ixlat,ixfld(ixo2)) +
     +             fields(ixlon,:,ixlat,ixfld(ixn2))
        elseif (ip.eq.ixroo2) then	! o/o2
          plt(:) = fields(ixlon,:,ixlat,ixfld(ixo1)) /
     +             fields(ixlon,:,ixlat,ixfld(ixo2))
        elseif (ip.eq.ixron2) then	! o/n2
          plt(:) = fields(ixlon,:,ixlat,ixfld(ixo1)) /
     +             fields(ixlon,:,ixlat,ixfld(ixn2))
        elseif (ip.eq.ixro2n2) then	! o2/n2
          plt(:) = fields(ixlon,:,ixlat,ixfld(ixo2)) /
     +             fields(ixlon,:,ixlat,ixfld(ixn2))
        elseif (ip.eq.ixroo2n2) then	! o/(o2+n2)
          plt(:) = fields(ixlon,:,ixlat,ixfld(ixo1)) /
     +             (fields(ixlon,:,ixlat,ixfld(ixo2))+
     +              fields(ixlon,:,ixlat,ixfld(ixn2)))
        endif
c
c Zonal means at a selected latitude:
c
      elseif (ixlat.gt.0.and.ixlon.eq.0) then
        if (ip.le.nfhist) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ixfld(ip))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.ge.ixui.and.ip.le.ixwi) then
          if (ip.eq.ixui) then
            do k=1,kmx
              rimx(:) = glbui(:,k,ixlat)
              plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          endif
          if (ip.eq.ixvi) then
            do k=1,kmx
              rimx(:) = glbvi(:,k,ixlat)
              plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          endif
          if (ip.eq.ixwi) then
            do k=1,kmx
              rimx(:) = glbwi(:,k,ixlat)
              plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          endif
        elseif (ip.eq.ixrho) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ixfld(ixo1)) +
     +                fields(:,k,ixlat,ixfld(ixo2)) +
     +                fields(:,k,ixlat,ixfld(ixn2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixroo2) then	! o/o2
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ixfld(ixo1)) /
     +                fields(:,k,ixlat,ixfld(ixo2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixron2) then	! o/n2
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ixfld(ixo1)) /
     +                fields(:,k,ixlat,ixfld(ixn2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixro2n2) then	! o2/n2
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ixfld(ixo2)) /
     +                fields(:,k,ixlat,ixfld(ixn2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixroo2n2) then	! o/(o2+n2)
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ixfld(ixo1)) /
     +                (fields(:,k,ixlat,ixfld(ixo2))+
     +                 fields(:,k,ixlat,ixfld(ixn2)))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        endif
c
c Global means:
c
      else
        if (ip.le.nfhist) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ip))
            plt(k) = calcmean(glb,imx*jmx,log,1.e-20,cpspval)
          enddo
        elseif (ip.ge.ixui.and.ip.le.ixwi) then
          do k=1,kmx
            if (ip.eq.ixui) glb(:,:) = glbui(:,k,:)
            if (ip.eq.ixvi) glb(:,:) = glbvi(:,k,:)
            if (ip.eq.ixwi) glb(:,:) = glbwi(:,k,:)
            plt(k) = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixrho) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixo1))
            plt(k) = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ixfld(ixo2))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ixfld(ixn2))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixroo2) then	! o/o2
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixo1)) /
     +                 fields(:,k,:,ixfld(ixo2))
            plt(k) = calcmean(glb,imx*jmx,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixron2) then	! o/n2
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixo1)) /
     +                 fields(:,k,:,ixfld(ixn2))
            plt(k) = calcmean(glb,imx*jmx,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixro2n2) then	! o2/n2
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixo2)) /
     +                 fields(:,k,:,ixfld(ixn2))
            plt(k) = calcmean(glb,imx*jmx,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.ixroo2n2) then	! o/(o2+n2)
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixo1)) /
     +        (fields(:,k,:,ixfld(ixo2))+fields(:,k,:,ixfld(ixn2)))
            plt(k) = calcmean(glb,imx*jmx,log,1.e-20,cpspval)
          enddo
        endif
      endif
      return
      end
