c
      subroutine getfld(fh,nlon,nzp,nlat,nflds,fout,id1,id2,type,
     +  glat,glon,zp,zpht,logint,ixhist,ip,ier)
c
c Return slice of fh in fout according to type:
c On input:
c   fh(nlon,nzp,nlat,nflds) contains fields on the 3-d global grid
c   ixhist(mxfproc) = field indices to fh
c   type = char slice type (GLOBAL, LONSLICE, LATSLICE, or VERTPROF)
c   glat,glon,zp = values of "fixed" dimension of fh, i.e.:
c     if type=GLOBAL  , then return fout(nlon,nlat) at zpht (zp or ht)
c     if type=LONSLICE, then return fout(nlat,nzp) at glon longitude
c     if type=LATSLICE, then return fout(nlon,nzp) at glat latitude
c     if type=VERTPROF, then return fout(nzp) at glat,glon location
c On output:
c   fout(id1,id2) is defined according to type
c   glat,glon rounded to nearest grid point
c   ier = 0 if no ier   
c
      include "flxproc.h"
      real fh(nlon,nzp,nlat,nflds),fout(id1,id2),fnzp(nzp),zp(nzp)
      integer ixhist(mxfproc)
      pointer(pfglb,fglb(nlon,nlat)),(pflatcol,flatcol(nlon,nzp))
      character*8 type
c
      if (type.ne.'GLOBAL  '.and.type.ne.'LONSLICE'.and.
     +    type.ne.'LATSLICE'.and.type.ne.'VERTPROF') then
        write(6,"('>>> WARNIING getfld: unknown type=',a)")
     +    type
        ier = 1
        return
      endif
c
c Return global fout(nlon,nlat) at zpht (pressure or height surface)
c (getglb does height interp if necessary)
c
      if (type.eq.'GLOBAL  ') then
        if (id1.ne.nlon.or.id2.ne.nlat) then
          write(6,"('>>> WARNING getfld (global): id1 should ',
     +      '=',i3,' id2 should=',i3)") nlon,nlat
          write(6,"('    id1=',i3,' id2=',i3)") id1,id2
          ier = 2
          return
        endif
        call getglb(fh,nlon,nzp,nlat,nflds,fout,zp,zpht,ixhist,ip)
        return
      endif
c
c Return lon slice fout(nlat,nzp) at glon (or zonal means if glon = zmflag)
c (also change glon to be nearest grid point longitude)
c
      if (type.eq.'LONSLICE') then
        if (glon.ne.zmflag) then
          do j=1,nlat        
            call getcol(fh,nlon,nzp,nlat,nflds,fnzp,gcmlat(j),glon,
     +        ixhist,ip)
            fout(j,:) = fnzp(:)
          enddo
        else	! zonal means
          call alloc(pflatcol,nlon*nzp)
          do j=1,nlat
            call getlatsli(fh,nlon,nzp,nlat,flatcol,gcmlat(j),
     +        ixhist,ip)
            do k=1,nzp
              fout(j,k) = fmean(flatcol(1,k),nlon,0,1.e-20,spval,0)
            enddo
          enddo
          call hpdeallc(pflatcol,ier,1)
        endif
        return
      endif
c
c Return lat slice fout(nlon,nzp) or fout(nlon,nhtscale) at glat
c (also change glat to be nearest grid point latitude)
c
      if (type.eq.'LATSLICE') then
        call getlatsli(fh,nlon,nzp,nlat,fout,glat,ixhist,ip)
        return
      endif
c
c Return vertical profile fout(nzp) at glat,glon:
c If glon == spval but glat != spval then return zonal means at glat
c If glon = glat = zmflag then return global means
c (glat,glon are from user input xylocs(2,mxlocs), which has been
c  verified in getinp; write warning one is bad anyway)
c
      if (type.eq.'VERTPROF') then
        if (glon.ne.zmflag) then		! check range of glon
          ixlon = ixfind(gcmlon,nlon,glon,dlon)
          if (ixlon.le.0)
     +      write(6,"('>>> getfld (VERT): bad glon=',f9.2)") glon
        endif
        if (glat.ne.zmflag) then		! check range of glat
          ixlat = ixfind(gcmlat,nlat,glat,dlat)
          if (ixlat.le.0)
     +      write(6,"('>>> getfld (VERT): bad glat=',f9.2)") glat
        endif
        if (glat.ne.zmflag.and.glon.ne.zmflag) then	! not global means
          glat = gcmlat(ixlat)
          glon = gcmlon(ixlon)
          call getcol(fh,nlon,nzp,nlat,nflds,fout(1,1),glat,glon,
     +      ixhist,ip)            
        else
          if (glat.ne.zmflag) then	! zonal means 
            glat = gcmlat(ixlat)
            call alloc(pflatcol,nlon*nzp)
            call getlatsli(fh,nlon,nzp,nlat,flatcol,glat,ixhist,ip)
            do k=1,nzp
              fout(k,1) = fmean(flatcol(1,k),nlon,0,1.e-20,spval,0)
            enddo 
            call hpdeallc(pflatcol,ier,1)
          else				! global means (fglbm is in util.a)
            call alloc(pfglb,nlon*nlat)
            do k=1,nzp
              call getglb(fh,nlon,nzp,nlat,nflds,fglb,zp,zp(k),
     +          ixhist,ip)
              fout(k,1) = fglbm(fglb,nlon,nlat,gcmlat,dlat,dlon,spval)
            enddo
            call hpdeallc(pfglb,ier,1)
          endif
        endif
      endif
      return
      end