
      subroutine mkfields(fields,utkmx,htutkmx,utjmx,it)
      include 'tigcmloc.h'
      dimension fields(imx,kmx,jmx,nfget), utkmx(ntms,kmx,nfplt,nloc),
     +  htutkmx(ntms,kmx,nloc),utjmx(ntms,jmx,nfplt,nzphtlon),
     +  glb(imx,kmx,jmx),xne(kmx),hts(kmx)
c
c Define plot array utkmx(ntms,kmx,nfget,nloc) 
c Will also need heights in htutkmx(ntms,kmx,nloc)
c
      if (nloc.gt.0) then
        do l=1,nloc
          call defutkmx(fields(1,1,1,ifget(ixz)),htutkmx(1,1,l),
     +      gloc(1,l),gloc(2,l),it)
          do ip=1,nftot 
            if (ifplt(ip).gt.0) then
              if (ip.le.nfhist) then
                call defutkmx(fields(1,1,1,ifget(ip)),
     +            utkmx(1,1,ifplt(ip),l),gloc(1,l),gloc(2,l),it)
              else
c
c Desired field is not on history, and must be calculated
c (fof2 and hmf2 will be used only by linezpht)
c
                if (ip.ge.ixui.and.ip.le.ixwi) then
                  if (ip.eq.ixui) iflag = 1
                  if (ip.eq.ixvi) iflag = 2
                  if (ip.eq.ixwi) iflag = 3
                  call mkdrifts(fields,iflag,glb,ier)
                  call defutkmx(glb,utkmx(1,1,ifplt(ip),l),gloc(1,l),
     +                          gloc(2,l),it)
                elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
                  do i=1,imx
                    do j=1,jmx
                      xne(:) = fields(i,:,j,ifget(ixne))
                      hts(:) = fields(i,:,j,ifget(ixz))
                      call fof2int(xne,hts,kmx,hmf2out,fof2out,0,i,j)
                      if (ip.eq.ixfof2) glb(i,:,j) = fof2out
                      if (ip.eq.ixhmf2) glb(i,:,j) = hmf2out
                    enddo
                  enddo
                  call defutkmx(glb,utkmx(1,1,ifplt(ip),l),gloc(1,l),
     +                          gloc(2,l),it)
                elseif (ip.eq.ixrho) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) + 
     +                         fields(:,:,:,ifget(ixo2)) + 
     +                         fields(:,:,:,ifget(ixn2))
                  call defutkmx(glb,utkmx(1,1,ifplt(ip),l),gloc(1,l),
     +                          gloc(2,l),it)
                elseif (ip.eq.ixamr) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) /
     +             (fields(:,:,:,ifget(ixo2))+fields(:,:,:,ifget(ixn2)))
                  call defutkmx(glb,utkmx(1,1,ifplt(ip),l),gloc(1,l),
     +                          gloc(2,l),it)
                elseif (ip.eq.ixon2r) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) /
     +              fields(:,:,:,ifget(ixn2))
                  call defutkmx(glb,utkmx(1,1,ifplt(ip),l),gloc(1,l),
     +                          gloc(2,l),it)
                endif
              endif
            endif
          enddo
        enddo
      endif
c
c Define plot array utjmx(ntms,jmx,nfplt,nzphtlon)
c
      if (nzphtlon.gt.0) then 
c       izfield = (ifget(ixz)-1)*imx*kmx*jmx + 1
        do l=1,nzphtlon
          do ip=1,nftot
            if (ifplt(ip).gt.0) then
              if (ip.le.nfhist) then
                call defutjmx(fields(1,1,1,ifget(ip)),
     +            fields(1,1,1,ifget(ixz)),utjmx(1,1,ifplt(ip),l),
     +            zphtlon(1,l),izphtlon(1,l),zphtlon(2,l),
     +            izphtlon(2,l),it,ip)
c
c Field not on history -- must be calculated from history fields:
c
              else
                if (ip.ge.ixui.and.ip.le.ixwi) then
                  if (ip.eq.ixui) iflag = 1
                  if (ip.eq.ixvi) iflag = 2
                  if (ip.eq.ixwi) iflag = 3
                  call mkdrifts(fields,iflag,glb,ier)
                  call defutjmx(glb,fields(1,1,1,ifget(ixz)),
     +              utjmx(1,1,ifplt(ip),l),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
                  do i=1,imx
                    do j=1,jmx
                      xne(:) = fields(i,:,j,ifget(ixne))
                      hts(:) = fields(i,:,j,ifget(ixz))
                      call fof2int(xne,hts,kmx,hmf2out,fof2out,0,i,j)
                      if (ip.eq.ixfof2) glb(i,:,j) = fof2out
                      if (ip.eq.ixhmf2) glb(i,:,j) = hmf2out
                    enddo
                  enddo
                  call defutjmx(glb,fields(1,1,1,ifget(ixz)),
     +              utjmx(1,1,ifplt(ip),l),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                elseif (ip.eq.ixrho) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) + 
     +                         fields(:,:,:,ifget(ixo2)) + 
     +                         fields(:,:,:,ifget(ixn2))
                  call defutjmx(glb,fields(1,1,1,ifget(ixz)),
     +              utjmx(1,1,ifplt(ip),l),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                elseif (ip.eq.ixamr) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) /
     +             (fields(:,:,:,ifget(ixo2))+fields(:,:,:,ifget(ixn2)))
                  call defutjmx(glb,fields(1,1,1,ifget(ixz)),
     +              utjmx(1,1,ifplt(ip),l),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                elseif (ip.eq.ixon2r) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) /
     +              fields(:,:,:,ifget(ixo1))
                  call defutjmx(glb,fields(1,1,1,ifget(ixz)),
     +              utjmx(1,1,ifplt(ip),l),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                endif
              endif
            endif
          enddo
        enddo
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mkdrifts(fields,ivel,vel,ier)
      include 'tigcmloc.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 (ifget(ixz).le.0.or.ifget(ixpot).le.0) then
          if (ifget(ixz).le.0) write(6,"('>>> mkdrifts: need heights')")
          if (ifget(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,ifget(ixpot)),fields(1,1,1,ifget(ixz)),
     +  gcmlat,vel,ivel,lumag)
c     
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getloc(fields,plt,ixlat,ixlon,log,ip)
      include 'tigcmloc.h'
      dimension plt(kmx),fields(imx,kmx,jmx,nfget),glb(imx,jmx)
c
c Return kmx column in plt from fields at lat index ixlat, lon index
c ixlon.  If ixlon=int(gmflag), return zonal means; 
c if ixlon=ixlat=int(gmflag) return global means
c
      igmflag = ifix(gmflag)
      if (ixlat.ne.igmflag.and.ixlon.ne.igmflag) then
        plt(:) = fields(ixlon,:,ixlat,ifget(ip))  
      elseif (ixlat.ne.igmflag.and.ixlon.eq.igmflag) then
        do k=1,kmx
          plt(k) = calcmean(fields(1,k,ixlat,ifget(ip)),log,1.e-20,
     +      cpspval)
        enddo
      else
        do k=1,kmx
          glb(:,:) = fields(:,k,:,ifget(ip))
          plt(k) = calcmean(glb,imx*jmx,log,1.e-20,cpspval)
        enddo
      endif
      return
      end
