c
      subroutine getglb(fields,plt,ixpl,ht,ip,ipert)
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 'tigcmdif.h'
      dimension plt(imx,jmx), fields(imx,kmx,jmx,nfget),plt0(imx,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 (ipert.gt.0) 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
        else
          if (ixpl.gt.0) then
            if (ip.eq.ixui) plt(:,:) = glbui0(:,ixpl,:)
            if (ip.eq.ixvi) plt(:,:) = glbvi0(:,ixpl,:)
            if (ip.eq.ixwi) plt(:,:) = glbwi0(:,ixpl,:)
          else
            if (ip.eq.ixui) 
     +        call glbhtint(glbui0,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +          plt,ht,1,ilog(ip),cpspval,ier,0)
            if (ip.eq.ixvi) 
     +        call glbhtint(glbvi0,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +          plt,ht,1,ilog(ip),cpspval,ier,0)
            if (ip.eq.ixwi) 
     +        call glbhtint(glbwi0,fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +          plt,ht,1,ilog(ip),cpspval,ier,0)
          endif
        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 vpar (neutral wind parallel to magnetic field):
c vmag = u*sin(dec)+v*cos(dec) and vpar = vmag*cos(dip)
c (for ht interp, u and v are interpolated before vpar is calculated, 
c  and vertical variation of dip and dec are ignored)
c
      elseif (ip.eq.ixvpar) then
        if (ixpl.gt.0) then
          plt(:,:) = (fields(:,ixpl,:,ixfld(ixu))*sin(dipdec(:,:,2))+
     +                fields(:,ixpl,:,ixfld(ixv))*cos(dipdec(:,:,2)))*
     +               cos(dipdec(:,:,1))
        else
          call glbhtint(fields(1,1,1,ixfld(ixu)),
     +                  fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt,ht,1,ilog(ip),cpspval,ier,0)
          call glbhtint(fields(1,1,1,ixfld(ixv)),
     +                  fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt0,ht,1,ilog(ip),cpspval,ier,0)
          plt(:,:) = (plt(:,:)*sin(dipdec(:,:,2))+
     +               plt0(:,:)*cos(dipdec(:,:,2)))*cos(dipdec(:,:,1))
        endif
c
c wpar (vertical component of vpar):
c vmag = u*sin(dec)+v*cos(dec) and wpar = -vmag*sin(dip)*cos(dip)
c (for ht interp, u and v are interpolated before vpar is calculated, 
c  and vertical variation of dip and dec are ignored)
c 
      elseif (ip.eq.ixwpar) then
        if (ixpl.gt.0) then
          plt(:,:) =-((fields(:,ixpl,:,ixfld(ixu))*sin(dipdec(:,:,2))+
     +                fields(:,ixpl,:,ixfld(ixv))*cos(dipdec(:,:,2)))*
     +                sin(dipdec(:,:,1))*cos(dipdec(:,:,1)))
        else
          call glbhtint(fields(1,1,1,ixfld(ixu)),
     +                  fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt,ht,1,ilog(ip),cpspval,ier,0)
          call glbhtint(fields(1,1,1,ixfld(ixv)),
     +                  fields(1,1,1,ixfld(ixz)),
     +      imx,kmx,jmx,plt0,ht,1,ilog(ip),cpspval,ier,0)
          plt(:,:) = -((plt(:,:)*sin(dipdec(:,:,2))+
     +               plt0(:,:)*cos(dipdec(:,:,2)))*sin(dipdec(:,:,1))*
     +               cos(dipdec(:,:,1)))
        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,glon,log,ip,ipert)
      include 'tigcmdif.h'
      dimension plt(jmx,kmx),rimx(imx),fields(imx,kmx,jmx,nfget),
     +  rimx0(imx),rkmx0(kmx),rkmx1(kmx),rkmx2(kmx)
      dimension xne(kmx),hts(kmx)
c
      if (ip.le.0) then
        write(6,"('>>> getlon: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      endif
      ixlon = ifix(zmflag)
      if (glon.ne.zmflag) then
        ixlon = ixfind(gcmlon,imx,glon,dlon)
        if (ixlon.le.0) then
          write(6,"('>>> getlon: bad longitude=',f10.3,
     +      ' -- returning zeroes')") glon
          plt(:,:) = 0.
          return
        endif
      endif
      if (ip.le.nfhist) then
        if (ixlon.ne.ifix(zmflag)) then
          if (ibilin.le.0) then
            do k=1,kmx
              do j=1,jmx
                plt(j,k) = fields(ixlon,k,j,ixfld(ip))
              enddo
            enddo
          else
c           write(6,"('getlon bilin: glon=',f7.1,' ip=',i2,
c    +        ' ipert=',i2)") glon,ip,ipert
            do j=1,jmx
              call locbilin(fields(1,1,1,ixfld(ip)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
              plt(j,:) = rkmx0(:)
            enddo
          endif
          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 (ipert.gt.0) then
          if (ixlon.ne.ifix(zmflag)) then
            if (ibilin.le.0) 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 bilinear interp to glon at each grid lat:
              do j=1,jmx
                if (ip.eq.ixui) call locbilin(glbui,imx,kmx,jmx,
     +              gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
                if (ip.eq.ixvi) call locbilin(glbvi,imx,kmx,jmx,
     +              gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
                if (ip.eq.ixwi) call locbilin(glbwi,imx,kmx,jmx,
     +              gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
                plt(j,:) = rkmx0(:)
              enddo
            endif
          else ! zonal means of perturbed case
            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
        else ! control
          if (ixlon.ne.ifix(zmflag)) then
            if (ibilin.le.0) then
              do k=1,kmx
                do j=1,jmx
                  if (ip.eq.ixui) plt(j,k) = glbui0(ixlon,k,j) 
                  if (ip.eq.ixvi) plt(j,k) = glbvi0(ixlon,k,j) 
                  if (ip.eq.ixwi) plt(j,k) = glbwi0(ixlon,k,j) 
                enddo
              enddo
            else ! do bilinear interp to glon at each grid lat:
              do j=1,jmx
                if (ip.eq.ixui) call locbilin(glbui0,imx,kmx,jmx,
     +              gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
                if (ip.eq.ixvi) call locbilin(glbvi0,imx,kmx,jmx,
     +              gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
                if (ip.eq.ixwi) call locbilin(glbwi0,imx,kmx,jmx,
     +              gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
                plt(j,:) = rkmx0(:)
              enddo
            endif
          else ! zonal means of base case
            do k=1,kmx
              do j=1,jmx
                if (ip.eq.ixui) rimx(:) = glbui0(:,k,j)
                if (ip.eq.ixvi) rimx(:) = glbvi0(:,k,j)
                if (ip.eq.ixwi) rimx(:) = glbwi0(:,k,j)
                plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0)
              enddo
            enddo
          endif
        endif
c
c fof2 and hmf2 are height independent -- return redundant values in vertical:
c
      elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
c       write(6,"('>>> getlon: ip=',i3,' ixfof2 and ixhmf2 are ',
c    +    'height independent -- returning zeroes')") ip
c       plt(:,:) = 0.
c       return
        if (ixlon.ne.ifix(zmflag)) then
          if (ibilin.le.0) then
            do j=1,jmx
              xne(:) = fields(ixlon,:,j,ixfld(ixne))
              hts(:) = fields(ixlon,:,j,ixfld(ixz))
              call fof2int(xne,hts,kmx,hmf2out,fof2out,0,ixlon,j)
              do k=1,kmx ! redundant in vertical
                if (ip.eq.ixfof2) plt(j,k) = fof2out
                if (ip.eq.ixhmf2) plt(j,k) = hmf2out
              enddo
            enddo
          else
            do j=1,jmx
              call locbilin(fields(1,1,1,ixfld(ixne)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,xne,cpspval,0)
              call locbilin(fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,hts,cpspval,0)
              call fof2int(xne,hts,kmx,hmf2out,fof2out,0,ixlon,j)
              do k=1,kmx ! redundant in vertical
                if (ip.eq.ixfof2) plt(j,k) = fof2out
                if (ip.eq.ixhmf2) plt(j,k) = hmf2out
              enddo
            enddo
          endif
        else ! zonal means of fof2,hmf2:
          do j=1,jmx
            do k=1,kmx
              rimx(:) = fields(:,k,j,ixfld(ixne))
              xne(k) = calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ixfld(ixz))
              hts(k) = calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
            enddo
            call fof2int(xne,hts,kmx,hmf2out,fof2out,0,ixlon,j)
            do k=1,kmx
              if (ip.eq.ixfof2) plt(j,k) = fof2out
              if (ip.eq.ixhmf2) plt(j,k) = hmf2out
            enddo
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
c
c rho (o2+o+n2):
c
      elseif (ip.eq.ixrho) then
        if (ixlon.ne.ifix(zmflag)) then
          if (ibilin.le.0) then
            do k=1,kmx
              plt(:,k) = fields(ixlon,k,:,ixfld(ixo1)) +
     +                   fields(ixlon,k,:,ixfld(ixo2)) +
     +                   fields(ixlon,k,:,ixfld(ixn2))
            enddo
          else
            do j=1,jmx
              call locbilin(fields(1,1,1,ixfld(ixo1)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
              call locbilin(fields(1,1,1,ixfld(ixo2)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx1,cpspval,0)
              call locbilin(fields(1,1,1,ixfld(ixn2)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx2,cpspval,0)
              plt(j,:) = rkmx0(:)+rkmx1(:)+rkmx2(:)
            enddo
          endif
        else ! zmeans of rho
          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 vpar (neutral wind parallel to magnetic field):
c vmag = u*sin(dec)+v*cos(dec) and vpar = vmag*cos(dip)
c (for zonal means, u, v, dip, and dec are averaged before vpar is calculated)
c (vertical variation of dip and dec are ignored)
c
      elseif (ip.eq.ixvpar) then
        if (ixlon.ne.ifix(zmflag)) then
          if (ibilin.le.0) then
            do k=1,kmx
              plt(:,k) = (fields(ixlon,k,:,ixfld(ixu))*
     +          sin(dipdec(ixlon,:,2))+fields(ixlon,k,:,ixfld(ixv))*
     +          cos(dipdec(ixlon,:,2)))*cos(dipdec(ixlon,:,1))
            enddo
          else ! bilin to glon at grid lat (dont bother to interp dipdec)
            do j=1,jmx
              call locbilin(fields(1,1,1,ixfld(ixu)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
              call locbilin(fields(1,1,1,ixfld(ixv)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx1,cpspval,0)
              plt(j,:) = (rkmx0(:)*sin(dipdec(ixlon,j,2))+rkmx1(:)*
     +          cos(dipdec(ixlon,j,2)))*cos(dipdec(ixlon,j,1))
            enddo
          endif
        else ! zonal means of vpar
          do k=1,kmx
            do j=1,jmx
              uzm = calcmean(fields(1,k,j,ixfld(ixu)),imx-1,0,1.e-20,
     +          cpspval,0)
              vzm = calcmean(fields(1,k,j,ixfld(ixv)),imx-1,0,1.e-20,
     +          cpspval,0)
              if (k.eq.1) then
                dipzm = calcmean(dipdec(1,j,1),imx-1,0,1.e-20,cpspval,0)
                deczm = calcmean(dipdec(1,j,2),imx-1,0,1.e-20,cpspval,0)
              endif
              plt(j,k) = (uzm*sin(deczm)+vzm*cos(deczm))*cos(dipzm)
            enddo
          enddo
        endif
c
c wpar (vertical component of vpar):
c vmag = u*sin(dec)+v*cos(dec) and wpar = -vmag*sin(dip)*cos(dip)
c (for zonal means, u, v, dip, and dec are averaged before wpar is calculated)
c (vertical variation of dip and dec are ignored)
c
      elseif (ip.eq.ixwpar) then
        if (ixlon.ne.ifix(zmflag)) then
          if (ibilin.le.0) then
            do k=1,kmx
              plt(:,k) = -((fields(ixlon,k,:,ixfld(ixu))*
     +          sin(dipdec(ixlon,:,2))+fields(ixlon,k,:,ixfld(ixv))*
     +          cos(dipdec(ixlon,:,2)))*sin(dipdec(ixlon,:,1))*
     +          cos(dipdec(ixlon,:,1)))
            enddo
          else
            do j=1,jmx
              call locbilin(fields(1,1,1,ixfld(ixu)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx0,cpspval,0)
              call locbilin(fields(1,1,1,ixfld(ixv)),imx,kmx,jmx,
     +          gcmlat,gcmlon,gcmlat(j),glon,rkmx1,cpspval,0)
              plt(j,:) = -((rkmx0(:)*sin(dipdec(ixlon,j,2))+rkmx1(:)*
     +          cos(dipdec(ixlon,j,2)))*sin(dipdec(ixlon,j,1))*
     +          cos(dipdec(ixlon,j,1)))
            enddo
          endif
        else ! zonal means of wpar
          do k=1,kmx
            do j=1,jmx
              uzm = calcmean(fields(1,k,j,ixfld(ixu)),imx-1,0,1.e-20,
     +          cpspval,0)
              vzm = calcmean(fields(1,k,j,ixfld(ixv)),imx-1,0,1.e-20,
     +          cpspval,0)
              if (k.eq.1) then
                dipzm = calcmean(dipdec(1,j,1),imx-1,0,1.e-20,cpspval,0)
                deczm = calcmean(dipdec(1,j,2),imx-1,0,1.e-20,cpspval,0)
              endif
              plt(j,k) = -((uzm*sin(deczm)+vzm*cos(deczm))*sin(dipzm)*
     +                    cos(dipzm))
            enddo
          enddo
        endif
      else
        write(6,"('>>> getlon: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getlat(fields,plt,glat,log,ip,ipert)
      include 'tigcmdif.h'
      dimension plt(imx,kmx),fields(imx,kmx,jmx,nfget),rkmx0(kmx),
     +  rkmx1(kmx),rkmx2(kmx)
c
      if (ip.le.0) then
        write(6,"('>>> getlat: unsupported ip=',i3)") ip
        plt(:,:) = 0.
        return
      endif
      ixlat = ixfind(gcmlat,jmx,glat,dlat)
      if (ixlat.le.0) then
        write(6,"('>>> getlat: bad latitude=',f10.3,
     +    ' -- returning zeroes')") glat
        plt(:,:) = 0.
        return
      endif
      if (ip.le.nfhist) then
        if (ibilin.le.0) then
          plt(:,:) = fields(:,:,ixlat,ixfld(ip))
        else
          do i=1,imx
            call locbilin(fields(1,1,1,ixfld(ip)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
            plt(i,:) = rkmx0(:) 
          enddo
        endif
      elseif (ip.ge.ixui.and.ip.le.ixwi) then
        if (ibilin.le.0) then
          if (ipert.gt.0) then
            if (ip.eq.ixui) plt(:,:) = glbui(:,:,ixlat)
            if (ip.eq.ixvi) plt(:,:) = glbvi(:,:,ixlat)
            if (ip.eq.ixwi) plt(:,:) = glbwi(:,:,ixlat)
          else
            if (ip.eq.ixui) plt(:,:) = glbui0(:,:,ixlat)
            if (ip.eq.ixvi) plt(:,:) = glbvi0(:,:,ixlat)
            if (ip.eq.ixwi) plt(:,:) = glbwi0(:,:,ixlat)
          endif
        else
          if (ipert.gt.0) then
            do i=1,imx
              if (ip.eq.ixui) call locbilin(glbui,imx,kmx,jmx,
     +          gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
              if (ip.eq.ixvi) call locbilin(glbvi,imx,kmx,jmx,
     +          gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
              if (ip.eq.ixwi) call locbilin(glbwi,imx,kmx,jmx,
     +          gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
              plt(i,:) = rkmx0(:)
            enddo
          else
            do i=1,imx
              if (ip.eq.ixui) call locbilin(glbui0,imx,kmx,jmx,
     +          gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
              if (ip.eq.ixvi) call locbilin(glbvi0,imx,kmx,jmx,
     +          gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
              if (ip.eq.ixwi) call locbilin(glbwi0,imx,kmx,jmx,
     +          gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
              plt(i,:) = rkmx0(:)
            enddo
          endif
        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,"('>>> 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
        if (ibilin.le.0) then
          plt(:,:) = fields(:,:,ixlat,ixfld(ixo1)) +
     +               fields(:,:,ixlat,ixfld(ixo2)) +
     +               fields(:,:,ixlat,ixfld(ixn2))
        else
          do i=1,imx
            call locbilin(fields(1,1,1,ixfld(ixo1)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixo2)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx1,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixn2)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx2,cpspval,0)
            plt(i,:) = rkmx0(:)+rkmx1(:)+rkmx2(:)
          enddo
        endif
c
c vpar:
c
      elseif (ip.eq.ixvpar) then
        if (ibilin.le.0) then
          do k=1,kmx
            plt(:,k) = (fields(:,k,ixlat,ixfld(ixu))*
     +        cos(dipdec(:,ixlat,2))+fields(:,k,ixlat,ixfld(ixv))*
     +        sin(dipdec(:,ixlat,2)))*cos(dipdec(:,ixlat,1))
          enddo
        else
          do i=1,imx
            call locbilin(fields(1,1,1,ixfld(ixu)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixv)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx1,cpspval,0)
            plt(i,:) = (rkmx0(:)*cos(dipdec(i,ixlat,2))+rkmx1(:)*
     +        sin(dipdec(i,ixlat,2)))*cos(dipdec(i,ixlat,1))
          enddo
        endif
c
c wpar:
c
      elseif (ip.eq.ixwpar) then
        if (ibilin.le.0) then
          do k=1,kmx
            plt(:,k) = -((fields(:,k,ixlat,ixfld(ixu))*
     +        cos(dipdec(:,ixlat,2))+fields(:,k,ixlat,ixfld(ixv))*
     +        sin(dipdec(:,ixlat,2)))*sin(dipdec(:,ixlat,1))*
     +        cos(dipdec(:,ixlat,1)))
          enddo
        else
          do i=1,imx
            call locbilin(fields(1,1,1,ixfld(ixu)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixv)),imx,kmx,jmx,
     +        gcmlat,gcmlon(i),glat,glon,rkmx1,cpspval,0)
            plt(i,:) = -((rkmx0(:)*cos(dipdec(i,ixlat,2))+rkmx1(:)*
     +        sin(dipdec(i,ixlat,2)))*sin(dipdec(i,ixlat,1))*
     +        cos(dipdec(i,ixlat,1)))
          enddo
        endif
c
c Bad ip:
c
      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,glat,glon,log,ip,ipert)
      include 'tigcmdif.h'
      dimension plt(kmx),fields(imx,kmx,jmx,nfget),glb(imx,jmx)
      dimension rkmx0(kmx),rkmx1(kmx)
c
      if (ip.le.0.or.ip.gt.ixwpar) then
        write(6,"('>>> getloc: unsupported ip=',i3,' -- returning ',
     +    'zeroes')") ip
        plt(:) = 0.
        return
      endif
c
c At specified location:
c (do bilinear interpolation of ibilin > 0):
c
      if (glat.ne.gmflag.and.glon.ne.gmflag) then
        ixlat = ixfind(gcmlat,jmx,glat,dlat)
        if (ixlat.le.0) then
          write(6,"('>>> getloc: bad lat=',f10.3,' -- returning ',
     +      'zeroes')") glat 
          plt(:) = 0.
          return
        endif
        ixlon = ixfind(gcmlon,imx,glon,dlon)
        if (ixlon.le.0) then
          write(6,"('>>> getloc: bad lon=',f10.3,' -- returning ', 
     +      'zeroes')") glon
          plt(:) = 0.
          return
        endif
        if (ip.le.nfhist) then
          if (ibilin.le.0) then
            plt(:) = fields(ixlon,:,ixlat,ixfld(ip))
          else
            call locbilin(fields(1,1,1,ixfld(ip)),imx,kmx,jmx,gcmlat,
     +        gcmlon,glat,glon,plt,cpspval,0)
          endif
        elseif (ip.ge.ixui.and.ip.le.ixwi) then
          if (ibilin.le.0) then
            if (ipert.gt.0) 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)
            else
              if (ip.eq.ixui) plt(:) = glbui0(ixlon,:,ixlat)
              if (ip.eq.ixvi) plt(:) = glbvi0(ixlon,:,ixlat)
              if (ip.eq.ixwi) plt(:) = glbwi0(ixlon,:,ixlat)
            endif
          else
            if (ipert.gt.0) then
              if (ip.eq.ixui) call locbilin(glbui,imx,kmx,jmx,gcmlat,
     +          gcmlon,glat,glon,plt,cpspval,0)
              if (ip.eq.ixvi) call locbilin(glbvi,imx,kmx,jmx,gcmlat,
     +          gcmlon,glat,glon,plt,cpspval,0)
              if (ip.eq.ixwi) call locbilin(glbwi,imx,kmx,jmx,gcmlat,
     +          gcmlon,glat,glon,plt,cpspval,0)
            else
              if (ip.eq.ixui) call locbilin(glbui0,imx,kmx,jmx,gcmlat,
     +          gcmlon,glat,glon,plt,cpspval,0)
              if (ip.eq.ixvi) call locbilin(glbvi0,imx,kmx,jmx,gcmlat,
     +          gcmlon,glat,glon,plt,cpspval,0)
              if (ip.eq.ixwi) call locbilin(glbwi0,imx,kmx,jmx,gcmlat,
     +          gcmlon,glat,glon,plt,cpspval,0)
            endif
          endif
        elseif (ip.eq.ixrho) then
          if (ibilin.le.0) then
            plt(:) = fields(ixlon,:,ixlat,ixfld(ixo1)) +
     +               fields(ixlon,:,ixlat,ixfld(ixo2)) +
     +               fields(ixlon,:,ixlat,ixfld(ixn2))
          else ! (o2 in rkmx0, n2 in rkmx1)
            call locbilin(fields(1,1,1,ixfld(ixo1)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,plt,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixo2)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixn2)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx1,cpspval,0)
            plt(:) = plt(:) + rkmx0(:) + rkmx1(:)
          endif
        elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
          if (ibilin.le.0) then
            rkmx0(:) = fields(ixlon,:,ixlat,ixfld(ixne))
            rkmx1(:) = fields(ixlon,:,ixlat,ixfld(ixz))
            call fof2int(rkmx0,rkmx1,kmx,hmf2out,fof2out,0,ixlon,ixlat)
            if (ip.eq.ixfof2) plt(:) = fof2out
            if (ip.eq.ixhmf2) plt(:) = hmf2out
          else
            call locbilin(fields(1,1,1,ixfld(ixne)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixz)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx1,cpspval,0)
            call fof2int(rkmx0,rkmx1,kmx,hmf2out,fof2out,0,ixlon,ixlat)
            if (ip.eq.ixfof2) plt(:) = fof2out
            if (ip.eq.ixhmf2) plt(:) = hmf2out
          endif
        elseif (ip.eq.ixvpar) then
          if (ibilin.le.0) then
            plt(:) = (fields(ixlon,:,ixlat,ixfld(ixu))*
     +        sin(dipdec(ixlon,ixlat,2))+
     +        fields(ixlon,:,ixlat,ixfld(ixv))*
     +        cos(dipdec(ixlon,ixlat,2)))*cos(dipdec(ixlon,ixlat,1))
          else ! un in rkmx0, vn in rkmx1
            call locbilin(fields(1,1,1,ixfld(ixu)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixv)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx1,cpspval,0)
            plt(:) = (rkmx0(:)*sin(dipdec(ixlon,ixlat,2))+rkmx1(:)*
     +        cos(dipdec(ixlon,ixlat,2)))*cos(dipdec(ixlon,ixlat,1))
          endif
        elseif (ip.eq.ixwpar) then
          if (ibilin.le.0) then
            plt(:) = -((fields(ixlon,:,ixlat,ixfld(ixu))*
     +        sin(dipdec(ixlon,ixlat,2))+
     +        fields(ixlon,:,ixlat,ixfld(ixv))*
     +        cos(dipdec(ixlon,ixlat,2)))*cos(dipdec(ixlon,ixlat,1))*
     +        sin(dipdec(ixlon,ixlat,1)))
          else
            call locbilin(fields(1,1,1,ixfld(ixu)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx0,cpspval,0)
            call locbilin(fields(1,1,1,ixfld(ixv)),imx,kmx,jmx,
     +        gcmlat,gcmlon,glat,glon,rkmx1,cpspval,0)
            plt(:) = -((rkmx0(:)*sin(dipdec(ixlon,ixlat,2))+rkmx1(:)*
     +        cos(dipdec(ixlon,ixlat,2)))*cos(dipdec(ixlon,ixlat,1))*
     +        sin(dipdec(ixlon,ixlat,1)))
          endif
        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
          if (ipert.gt.0) 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
          else
            do k=1,kmx
              if (ip.eq.ixui) glb(:,:) = glbui0(:,k,:)
              if (ip.eq.ixvi) glb(:,:) = glbvi0(:,k,:)
              if (ip.eq.ixwi) glb(:,:) = glbwi0(:,k,:)
              plt(k) = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            enddo
          endif
        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.ixvpar) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixu))
            uglbm = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ixfld(ixv))
            vglbm = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            dipglbm = calcmean(dipdec,imx*jmx,0,1.e-20,cpspval)
            decglbm = calcmean(dipdec(1,1,2),imx*jmx,0,1.e-20,cpspval)
            plt(k) = (uglbm*sin(decglbm)+vglbm*cos(decglbm))*
     +        cos(dipglbm) 
          enddo
        elseif (ip.eq.ixwpar) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ixfld(ixu))
            uglbm = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ixfld(ixv))
            vglbm = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            dipglbm = calcmean(dipdec,imx*jmx,0,1.e-20,cpspval)
            decglbm = calcmean(dipdec(1,1,2),imx*jmx,0,1.e-20,cpspval)
            plt(k) = -((uglbm*sin(decglbm)+vglbm*cos(decglbm))*
     +        sin(dipglbm)*cos(dipglbm)) 
          enddo
        endif
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mkdrifts(fields,ivel,vel,ier)
      include 'tigcmdif.h'
      dimension fields(imx,kmx,jmx,nfget),vel(imx,kmx,jmx)
      character*32 magfile
      data ncalls/0/
      data magfile /'/FOSTER/magfield/magfield.cray  '/
      save ncalls
c
      ncalls = ncalls+1
c
c Use mkuivi (~/getgcm/lib) to make ui, vi, and/or wi, defining
c appropriate index in fields. (fields has already been allocated,
c including potential and height, which are necessary to calculate
c the drift velocities)
c
      ier = 0
c
c Get mag field strengths from mss:
c
      if (ncalls.eq.1) then
        write(6,"(' ')")
        write(6,"('mkdrifts: getting file ',a,' from mss...')")
     +    magfile
        call getmss(magfile,lumag)
        write(6,"(' ')")
c
c Check indices to height and potential (both needed for ion drifts)
c
        if (ixfld(ixz).le.0.or.ixfld(ixpot).le.0) then
          if (ixfld(ixz).le.0) write(6,"('>>> mkdrifts: need heights')")
          if (ixfld(ixpot).le.0) 
     +      write(6,"('>>> mkdrifts: need potential')")
          ier = 1
          return
        endif 
      else
        rewind lumag
      endif
c
c Get desired component:
c
      call mkuivi(fields(1,1,1,ixfld(ixpot)),fields(1,1,1,ixfld(ixz)),
     +  gcmlat,vel,ivel,lumag)
c
      return
      end