
      subroutine mkfields(fields,utkmx,htutkmx,utkmxlt,hutkmxlt,
     +  utjmx,utjmxlt,it)
c
c Define plot array utkmx(ntms,kmx,nfplt,nloc) (output)
c Use fields, and heights in htutkmx(ntms,kmx,nloc) (both input)
c
      include 'timesloc.h'
      include 'tgcmhdr.h'
      dimension fields(imx,kmx,jmx,nfget), utkmx(ntms,kmx,nfplt,nloc),
     +  htutkmx(ntms,kmx,nloc),utjmx(ntms,jmx,nfplt,nzphtlon),
     +  hts(kmx),fjmx8(jmx,8),fion3(3),utjmxlt(ntms,jmx,nfplt,nzphtslt),
     +  utkmxlt(ntms,kmx,nfplt,nloclt),hutkmxlt(ntms,kmx,nloclt)
      logical didions
c
c Define utkmx at desired locations:
c
      if (nloc.gt.0) then
        do l=1,nloc
          glat = gloc(1,l)
          glon = gloc(2,l)
          call defutkmx(fields(1,1,1,ifget(ixz)),htutkmx(1,1,l),
     +      glat,glon,it)
          call mkutkmx(fields,htutkmx(1,1,l),utkmx(1,1,1,l),glat,
     +      glon,it,l,1)
        enddo
      endif
c
c Define utkmx at desired latitudes and fixed local times:
c (doppler fields not available at fixed local times)
c
      if (nloclt.gt.0) then
        do l=1,nloclt
          glat = gloclt(1,l)
          glon = fslt(gloclt(2,l),ut(it),0.,3) ! lon at required slt
          call defutkmx(fields(1,1,1,ifget(ixz)),hutkmxlt(1,1,l),
     +      glat,glon,it)
          call mkutkmx(fields,hutkmxlt(1,1,l),utkmxlt(1,1,1,l),glat,
     +      glon,it,l,0)
        enddo
      endif
c
c Define plot array utjmx(ntms,jmx,nfplt,nzphtlon)
c For ut vs jmx contours:
c
      if (nzphtlon.le.0) goto 2000 
      do l=1,nzphtlon
        didions = .false.
        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
c Ion drifts (any necessary ht interp is done in defjmx):
c (zonal means of ion drifts not available)
            else
              if (ip.ge.ixui.and.ip.le.ixwi) then
                if (.not.didions) then
                  if (ionvel.eq.2.and.zphtlon(2,l).eq.zmflag) then
                    goto 200
                    didions = .true.
                  endif
                  call defjmx(fields(1,1,1,ifget(ixu)),
     +              fields(1,1,1,ifget(ixz)),fjmx8,zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixv)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixw)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixo2)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixo1)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,5),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixn2)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,6),zphtlon(1,l),
     +              izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                  if (ionvel.eq.2) then
                    altgrf = 300.
                    yr = float(date(1))
                    if (yr.gt.1985.) yr = 1985.
                  endif
                  do j=1,jmx
                    if (ionvel.eq.2) call igrf(gcmlat(j),zphtlon(2,l),
     +                altgrf,yr,gmlat,gmlon,rdip,rdec)
                    rdip = rdip*dtr
                    rdec = rdec*dtr
                    call mkionvel(fion3,1,gcmlat(j),zphtlon(2,l),
     +                gmlat,gmlon,rdip,rdec,float(date(2)),ut(it),
     +                fjmx8(j,1),fjmx8(j,2),fjmx8(j,3),fjmx8(j,4),
     +                fjmx8(j,5),fjmx8(j,6),ionvel-1,cpspval)
                    utjmx(it,j,ifplt(ixui),l) = fion3(1)
                    utjmx(it,j,ifplt(ixvi),l) = fion3(2)
                    utjmx(it,j,ifplt(ixwi),l) = fion3(3)
                  enddo
                  didions = .true.
                endif
c
c E5577 (greenline volume emission, see ~/lib/mke5577):
c
              elseif (ip.eq.ixe5577) then
                call defjmx(fields(1,1,1,ifget(ixt)),
     +            fields(1,1,1,ifget(ixz)),fjmx8,zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo1)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixn2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call mke5577(fjmx8,fjmx8(1,2),fjmx8(1,3),fjmx8(1,4),
     +            jmx,fjmx8(1,6),cpspval)
                utjmx(it,:,ifplt(ixe5577),l) = fjmx8(:,6)
c
c E6300 (redline volume emission, see ~/lib/mke6300):
c
              elseif (ip.eq.ixe6300) then
                call defjmx(fields(1,1,1,ifget(ixt)),
     +            fields(1,1,1,ifget(ixz)),fjmx8,zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixte)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo1)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixn2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,5),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2p)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,6),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixne)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,7),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call mke6300(fjmx8,fjmx8(1,2),fjmx8(1,3),fjmx8(1,4),
     +            fjmx8(1,5),fjmx8(1,6),fjmx8(1,7),jmx,fjmx8(1,8),
     +            cpspval)
                utjmx(it,:,ifplt(ixe6300),l) = fjmx8(:,8)
c
c Eo200 (o2 (0-0) atmos band emission, see ~/lib/mkeo200):
c
              elseif (ip.eq.ixeo200) then
                call defjmx(fields(1,1,1,ifget(ixt)),
     +            fields(1,1,1,ifget(ixz)),fjmx8,zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo1)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixn2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtlon(1,l),
     +            izphtlon(1,l),zphtlon(2,l),izphtlon(2,l),it,ip)
                call mkeo200(fjmx8,fjmx8(1,2),fjmx8(1,3),fjmx8(1,4),
     +            jmx,fjmx8(1,6),cpspval)
                utjmx(it,:,ifplt(ixeo200),l) = fjmx8(:,6)
              else
                write(6,"('mkfields: unknown ip=',i3)") ip
              endif
            endif
          endif
 200      continue
        enddo ! fields ip loop
      enddo
 2000 continue
c
c Define plot array utjmxlt(ntms,jmx,nfplt,nzphtslt)
c For ut vs jmx contours:
c
      if (nzphtslt.le.0) goto 3000 
      do l=1,nzphtslt
        didions = .false.
        sltlon = fslt(zphtslt(2,l),ut(it),0.,3) ! lon at required slt
        izphtslt(2,l) = ixfind(gcmlon,imx,sltlon,dlon) ! index to sltlon
c       write(6,"('mkfields: l=',i2,' zphtslt=',2f8.2,' sltlon=',f8.2,
c    +    ' izphtslt=',2i4)") l,(zphtslt(ii,l),ii=1,2),sltlon,
c    +    (izphtslt(ii,l),ii=1,2)
        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)),utjmxlt(1,1,ifplt(ip),l),
     +          zphtslt(1,l),izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
c
c Field not on history -- must be calculated from history fields:
c
c Ion drifts (any necessary ht interp is done in defjmx):
c (zonal means of ion drifts not available)
            else
              if (ip.ge.ixui.and.ip.le.ixwi) then
                if (.not.didions) then
                  call defjmx(fields(1,1,1,ifget(ixu)),
     +              fields(1,1,1,ifget(ixz)),fjmx8,zphtslt(1,l),
     +              izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixv)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtslt(1,l),
     +              izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixw)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtslt(1,l),
     +              izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixo2)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtslt(1,l),
     +              izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixo1)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,5),zphtslt(1,l),
     +              izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                  call defjmx(fields(1,1,1,ifget(ixn2)),
     +              fields(1,1,1,ifget(ixz)),fjmx8(1,6),zphtslt(1,l),
     +              izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                  if (ionvel.eq.2) then
                    altgrf = 300.
                    yr = float(date(1))
                    if (yr.gt.1985.) yr = 1985.
                  endif
                  do j=1,jmx
                    if (ionvel.eq.2) call igrf(gcmlat(j),zphtslt(2,l),
     +                altgrf,yr,gmlat,gmlon,rdip,rdec)
                    rdip = rdip*dtr
                    rdec = rdec*dtr
                    call mkionvel(fion3,1,gcmlat(j),zphtslt(2,l),
     +                gmlat,gmlon,rdip,rdec,float(date(2)),ut(it),
     +                fjmx8(j,1),fjmx8(j,2),fjmx8(j,3),fjmx8(j,4),
     +                fjmx8(j,5),fjmx8(j,6),ionvel-1,cpspval)
                    utjmxlt(it,j,ifplt(ixui),l) = fion3(1)
                    utjmxlt(it,j,ifplt(ixvi),l) = fion3(2)
                    utjmxlt(it,j,ifplt(ixwi),l) = fion3(3)
                  enddo
                  didions = .true.
                endif
c
c E5577 (greenline volume emission, see ~/lib/mke5577):
c
              elseif (ip.eq.ixe5577) then
                call defjmx(fields(1,1,1,ifget(ixt)),
     +            fields(1,1,1,ifget(ixz)),fjmx8,zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo1)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixn2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call mke5577(fjmx8,fjmx8(1,2),fjmx8(1,3),fjmx8(1,4),
     +            jmx,fjmx8(1,6),cpspval)
                utjmxlt(it,:,ifplt(ixe5577),l) = fjmx8(:,6)
c
c E6300 (redline volume emission, see ~/lib/mke6300):
c
              elseif (ip.eq.ixe6300) then
                call defjmx(fields(1,1,1,ifget(ixt)),
     +            fields(1,1,1,ifget(ixz)),fjmx8,zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixte)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo1)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixn2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,5),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2p)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,6),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixne)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,7),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call mke6300(fjmx8,fjmx8(1,2),fjmx8(1,3),fjmx8(1,4),
     +            fjmx8(1,5),fjmx8(1,6),fjmx8(1,7),jmx,fjmx8(1,8),
     +            cpspval)
                utjmxlt(it,:,ifplt(ixe6300),l) = fjmx8(:,8)
c
c Eo200 (o2 (0-0) atmos band emission, see ~/lib/mkeo200):
c
              elseif (ip.eq.ixeo200) then
                call defjmx(fields(1,1,1,ifget(ixt)),
     +            fields(1,1,1,ifget(ixz)),fjmx8,zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,2),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixo1)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,3),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call defjmx(fields(1,1,1,ifget(ixn2)),
     +            fields(1,1,1,ifget(ixz)),fjmx8(1,4),zphtslt(1,l),
     +            izphtslt(1,l),sltlon,izphtslt(2,l),it,ip)
                call mkeo200(fjmx8,fjmx8(1,2),fjmx8(1,3),fjmx8(1,4),
     +            jmx,fjmx8(1,6),cpspval)
                utjmxlt(it,:,ifplt(ixeo200),l) = fjmx8(:,6)
              else
                write(6,"('mkfields: unknown ip=',i3)") ip
              endif
            endif
          endif
 300      continue
        enddo ! fields ip loop
      enddo
 3000 continue
      return
      end
