c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine defutkmx(field,utkmx,glat,glon,it)
      include 'tigcmloc.h'
      dimension field(imx,kmx,jmx), utkmx(kmx),glb(imx,jmx),
     +  rimx(imx),rk(kmx) 
c
c Define utkmx for current time (it). If ibilin > 0, then interpolate
c to desired location (no interp is done if doing global means. If
c doing zonal means, interp is done to desired latitude) 
c
      call getlocix(glat,glon,ilat,ilon)
c
c Change glat,glon to reflect actual grid point used, if not doing
c interpolation to location:
c
      if (it.eq.1.and.ibilin.le.0) then
        if (glat.ne.gmflag) glat = gcmlat(ilat)
        if (glon.ne.gmflag) glon = gcmlon(ilon)
      endif
c
c Want global means:
c
      if (ilat.eq.ifix(gmflag).and.ilon.eq.ifix(gmflag)) then
        do k=1,kmx
          glb(:,:) = field(:,k,:)
          utkmx(k) = fglbm(glb,imx,jmx,gcmlat,dlat,dlon,cpspval)
        enddo
c
c Want zonal means at selected latitude:
c (if ibilin > 0, do linear interpolation to desired latitude)
c
      elseif (ilon.eq.ifix(gmflag).and.ilat.ne.ifix(gmflag)) then
        if (ibilin.le.0) then
          do k=1,kmx
            rimx(:) = field(:,k,ilat) 
            utkmx(k) = calcmean(rimx,imx,0,1.e-20,cpspval,0)
          enddo
        else
          utkmx(:) = 0.
          do i=1,imx-1
            call locbilin(field,imx,kmx,jmx,gcmlat,gcmlon,glat,
     +        gcmlon(i),rk,cpspval,1)
            utkmx(:) = utkmx(:) + rk(:)
          enddo
          utkmx(:) = utkmx(:) / float(imx-1)
        endif
c
c Want lat,lon location:
c (if ibilin > 0, do linear interpolation to desired location)
c
      else
        if (ibilin.le.0) then
          utkmx(:) = field(ilon,:,ilat)
        else
          call locbilin(field,imx,kmx,jmx,gcmlat,gcmlon,glat,glon,
     +      rk,cpspval,1)
          utkmx(:) = rk(:)
        endif
      endif
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine defutjmx(field,hts,utjmx,zpht,ixzp,slon,ixlon,it,ip)
c
c Define utjmx(it,jmx) where "it" is index at current ut
c (if ibilin > 0, interpolate to desired longitude, otherwise use
c  ixlon, which is index to closest grid longitude)
c On input:
c   fields(imx,kmx,jmx) = original at global grid and current ut
c   hts(imx,kmx,jmx) = gcm heights at global grid and current ut
c   zpht = selected pressure or height
c   ixzp = index to zpht if its a pressure; ixzp=0 if zpht is a height
c   slon = selected longitude (if = gmflag, then get zonal means)
c   ixlon = index to slon; ixlon = ifix(gmflag) if want zonal means
c   it = time loop index
c   ip = current field index
c On output:
c   utjmx is defined, nothing else is changed
c
      include 'tigcmloc.h'
      dimension field(imx,kmx,jmx), utjmx(ntms,jmx), hts(imx,kmx,jmx),
     +  rimx(imx), fjkmx(jmx,kmx), hjkmx(jmx,kmx),rjmx(jmx), rkmx(kmx)
c
      igmflag = ifix(gmflag)
c
c At selected zp and longitude:
        if (ixzp.gt.0.and.ixlon.ne.igmflag) then
          if (ibilin.le.0) then
            utjmx(it,:) = field(ixlon,ixzp,:)
          else
            do j=1,jmx
              call locbilin(field,imx,kmx,jmx,gcmlat,gcmlon,
     +          gcmlat(j),slon,rkmx,cpspval,1)
              utjmx(it,j) = rkmx(ixzp)
            enddo
          endif
          if (ip.eq.ixfof2.or.ip.eq.ixhmf2) return
c
c At selected zp and zonal means:
        elseif (ixzp.gt.0.and.ixlon.eq.igmflag) then
          do j=1,jmx
            rimx(:) = field(:,ixzp,j)
            utjmx(it,j) = calcmean(rimx,imx,0,1.e-20,cpspval,0)
          enddo
          if (ip.eq.ixfof2.or.ip.eq.ixhmf2) return
c
c At selected height and longitude (if ibilin > 0, then interpolate
c to desired longitude before doing height interpolation):
c
        elseif (ixzp.le.0.and.ixlon.ne.igmflag.and.ip.ne.ixfof2.and.
     +          ip.ne.ixhmf2) then
          if (ibilin.le.0) then
            do k=1,kmx
              do j=1,jmx
                hjkmx(j,k) = hts(ixlon,k,j)
                fjkmx(j,k) = field(ixlon,k,j)
              enddo
            enddo
          else
            do j=1,jmx
              call locbilin(hts,imx,kmx,jmx,gcmlat,gcmlon,
     +          gcmlat(j),slon,rkmx,cpspval,1)
              hjkmx(j,:) = rkmx(:)     
              call locbilin(field,imx,kmx,jmx,gcmlat,gcmlon,
     +          gcmlat(j),slon,rkmx,cpspval,1)
              fjkmx(j,:) = rkmx(:)     
            enddo
          endif
          call cuthtint(fjkmx,hjkmx,jmx,kmx,rjmx,zpht,1,logterp(ip),
     +      cpspval,ier,0)
          utjmx(it,:) = rjmx(:)
c
c At selected height and zonal means (interpolate first, then take means)
        elseif (ixzp.le.0.and.ixlon.eq.igmflag.and.ip.ne.ixfof2.and.
     +          ip.ne.ixhmf2) then
          utjmx(it,:) = 0.
          do i=1,imx-1
            do k=1,kmx
              do j=1,jmx
                hjkmx(j,k) = hts(i,k,j)
                fjkmx(j,k) = field(i,k,j)
              enddo
            enddo
            call cuthtint(fjkmx,hjkmx,jmx,kmx,rjmx,zpht,1,logterp(ip),
     +        cpspval,ier,0)
            utjmx(it,:) = utjmx(it,:) + rjmx(:) 
          enddo
          utjmx(it,:) = utjmx(it,:) / float(imx-1)
        endif
c
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getlocix(glat,glon,ilat,ilon)
      include 'tigcmloc.h'
c
c Selected lat,lon location:
      if (glat.ne.gmflag.and.glon.ne.gmflag) then
        ilat = ixfind(gcmlat,jmx,glat,dlat)
        ilon = ixfind(gcmlon,imx,glon,dlon)
c
c Zonal means at selected latitude:
      elseif (glat.ne.gmflag.and.glon.eq.gmflag) then
        ilat = ixfind(gcmlat,jmx,glat,dlat)
        ilon = ifix(gmflag)
c
c Global means:
      elseif (glat.eq.gmflag.and.glon.eq.gmflag) then
        ilat = ifix(gmflag)
        ilon = ifix(gmflag)
      endif
c
c Bad glat or glon:
      if (ilat.le.0.or.ilon.le.0) then
        write(6,"('>>> getlocix: bad loc? glat,glon=',2f12.3)")
     +    glat,glon
      endif
      return
      end