c
      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),fkmx(kmx)
      pointer (pglb0,glb0(1)), (pglb1,glb1(1)), 
     +        (pvelion,velion(imx,kmx,jmx,3))
c
c Define plot array utkmx(ntms,kmx,nfplt,nloc) 
c Will also need heights in htutkmx(ntms,kmx,nloc)
c If ion drifts are requested, allocate space at pvelion, and for ve,vu: 
c
      if (nloc.gt.0) then
        if (ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0.) 
     +    then
          call alloc(pvelion,imx*kmx*jmx*3)
          if (ionvel.eq.2) then
            call alloc(pglb0,imx*kmx*jmx)
            call alloc(pglb1,imx*kmx*jmx)
            ivevu = 1
          else
            ivevu = 0
          endif
c
c Get global ion drifts:
c (note if ionvel < 2, fields at u,v,w,o2,o1,n2 may be undefined,
c  but will not be used by mkdrifts)
c
          do i=ixui,ixwi
            if (ifplt(i).gt.0) then
              if (i.eq.ixui) ix = 1
              if (i.eq.ixvi) ix = 2
              if (i.eq.ixwi) ix = 3
              call mkdrifts(ut(it),fields(1,1,1,ifget(ixpot)),
     +          fields(1,1,1,ifget(ixz)),gcmlat,velion(1,1,1,ix),
     +          imx,kmx,jmx,ix,ionvel,ivevu,lumag,ludipdec,
     +          fields(1,1,1,ifget(ixu)),fields(1,1,1,ifget(ixv)),
     +          fields(1,1,1,ifget(ixw)),fields(1,1,1,ifget(ixo2)),
     +          fields(1,1,1,ifget(ixo1)),fields(1,1,1,ifget(ixn2)),
     +          glb0,glb1,cpspval)
            endif
          enddo
        endif
        do l=1,nloc
          call defutkmx(fields(1,1,1,ifget(ixz)),fkmx,
     +      gloc(1,l),gloc(2,l),it)
          htutkmx(it,:,l) = fkmx(:)
          do ip=1,nftot 
            if (ifplt(ip).gt.0) then
              if (ip.le.nfhist) then
                call defutkmx(fields(1,1,1,ifget(ip)),fkmx,
     +            gloc(1,l),gloc(2,l),it)
                utkmx(it,:,ifplt(ip),l) = fkmx(:)
              else
c
c Desired field is not on history, and must be calculated
c (fof2 and hmf2 will be used only by linezpht)
c (ve and vu are found by mkuivi (which is called by mkdrifts), only
c  if ionvel=2. ve and vu will be plotted like 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 defutkmx(velion(1,1,1,iflag),fkmx,gloc(1,l),
     +              gloc(2,l),it)
                  utkmx(it,:,ifplt(ip),l) = fkmx(:)
                  if (ionvel.eq.2) then
                    call defutkmx(glb0,fkmx,gloc(1,l),gloc(2,l),it)
                    call xfervevu(veutkmx,ntms,kmx,nloc,fkmx,it,l)
                    call defutkmx(glb1,fkmx,gloc(1,l),gloc(2,l),it)
                    call xfervevu(vuutkmx,ntms,kmx,nloc,fkmx,it,l)
                  endif
                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,fkmx,gloc(1,l),gloc(2,l),it)
                  utkmx(it,:,ifplt(ip),l) = fkmx(:)
                elseif (ip.eq.ixrho) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) + 
     +                         fields(:,:,:,ifget(ixo2)) + 
     +                         fields(:,:,:,ifget(ixn2))
                  call defutkmx(glb,fkmx,gloc(1,l),gloc(2,l),it)
                  utkmx(it,:,ifplt(ip),l) = fkmx(:)
                elseif (ip.eq.ixamr) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) /
     +             (fields(:,:,:,ifget(ixo2))+fields(:,:,:,ifget(ixn2)))
                  call defutkmx(glb,fkmx,gloc(1,l),gloc(2,l),it)
                  utkmx(it,:,ifplt(ip),l) = fkmx(:)
                elseif (ip.eq.ixon2r) then
                  glb(:,:,:) = fields(:,:,:,ifget(ixo1)) /
     +              fields(:,:,:,ifget(ixn2))
                  call defutkmx(glb,fkmx,gloc(1,l),gloc(2,l),it)
                  utkmx(it,:,ifplt(ip),l) = fkmx(:)
                endif
              endif
            endif
          enddo
        enddo
        if (ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0.) 
     +    then
          call hpdeallc(pvelion,ier,1)
          if (ionvel.eq.2) then
            call hpdeallc(pglb0,ier,1) ! for ve
            call hpdeallc(pglb1,ier,1) ! for vu
          endif
        endif
      endif
c
c Define plot array utjmx(ntms,jmx,nfplt,nzphtlon)
c
      if (nzphtlon.gt.0) then 
        if (ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0.) 
     +    then
          ivevu = 0
          call alloc(pvelion,imx*kmx*jmx*3)
          do i=ixui,ixwi
            if (ifplt(i).gt.0) then
              if (i.eq.ixui) ix = 1
              if (i.eq.ixvi) ix = 2
              if (i.eq.ixwi) ix = 3
              call mkdrifts(ut(it),fields(1,1,1,ifget(ixpot)),
     +          fields(1,1,1,ifget(ixz)),gcmlat,velion(1,1,1,ix),
     +          imx,kmx,jmx,ix,ionvel,ivevu,lumag,ludipdec,
     +          fields(1,1,1,ifget(ixu)),fields(1,1,1,ifget(ixv)),
     +          fields(1,1,1,ifget(ixw)),fields(1,1,1,ifget(ixo2)),
     +          fields(1,1,1,ifget(ixo1)),fields(1,1,1,ifget(ixn2)),
     +          glb0,glb1,cpspval)
            endif
          enddo
        endif
        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 defutjmx(velion(1,1,1,iflag),
     +              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
        if (ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.ifplt(ixwi).gt.0.) 
     +    call hpdeallc(pvelion,ier,1)
      endif
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine xfervevu(v,ntms,kmx,nloc,fkmx,it,l)
      dimension v(ntms,kmx,nloc), fkmx(kmx)
c
      v(it,:,l) = fkmx(:)
      return
      end