
      subroutine mkfields(fields,utkmx,htutkmx,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),
     +  fkmx10(kmx,10),hts(kmx),fjmx8(jmx,8),fion3(3),
     +  utjmxlt(ntms,jmx,nfplt,nzphtslt)
      pointer(pfhts,fhts(1))
      logical didions
c
      if (nloc.le.0) goto 1000
      do l=1,nloc
        glat = gloc(1,l)
        if (locslt.le.0) then
          glon = gloc(2,l)
        else ! requested fixed local time
          glon = fslt(gloc(2,l),ut(it),0.,3) ! lon at required slt
        endif
        call defutkmx(fields(1,1,1,ifget(ixz)),htutkmx(1,1,l),
     +    glat,glon,it)
        didions = .false.
        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),glat,glon,it)
            else
c
c Desired field is not on history, and must be calculated:
c
c Ion drifts:
c (use fkmx10(kmx,10) for t,u,v,w,o2,o,n2)
c (igrf (from ~/lib/igrf) is used to determine dec and dip at loc)
c (note igrf extrapolates beyond 1985 -- avoid this)
c Zonal or global means of ion drifts with efield not available
c   (this would require finding average rdip,rdec by multiple calls
c    to igrf, and calling efield over all lons or entire globe, so 
c    is not worth it at this time)
c
              if (ip.ge.ixui.and.ip.le.ixwi) then
                if (.not.didions) then
                  if (ionvel.eq.2.and.(glat.eq.zmflag.or.
     +                glon.eq.zmflag)) then
                    didions = .true.
                    goto 100
                  endif ! means with efield not available
                  altgrf = 300.
                  yr = float(date(1))
                  if (yr.gt.1985.) yr = 1985.
                  call igrf(glat,glon,altgrf,yr,gmlat,gmlon,rdip,rdec)
                  rdip = rdip*dtr
                  rdec = rdec*dtr
                  call defkmx(fields(1,1,1,ifget(ixu)),fkmx10,
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,2),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixw)),fkmx10(1,3),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,4),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,5),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,6),
     +              glat,glon,it)
c
c     subroutine mkionvel(fout,kmx,glat,glon,gmlat,gmlon,rdip,rdec,
c    +  day,ut,u,v,w,xo2,xo,xn2,iefield,spv)
c
                  call mkionvel(fkmx10(1,8),kmx,glat,glon,gmlat,
     +              gmlon,rdip,rdec,float(date(2)),ut(it),fkmx10,
     +              fkmx10(1,2),fkmx10(1,3),fkmx10(1,4),fkmx10(1,5),
     +              fkmx10(1,6),ionvel-1,cpspval)
                  utkmx(it,:,ifplt(ixui),l) = fkmx10(:,8)
                  utkmx(it,:,ifplt(ixvi),l) = fkmx10(:,9)
                  utkmx(it,:,ifplt(ixwi),l) = fkmx10(:,10)
c
c If doing linear height scale on y-axis, we need to interpolate
c the 6 needed neutral fields first, then call mkionvel to calculate 
c drifts from the interpolated fields; this defines the array 
c htions(ntms,nhtscale,3,nloc) (phtions pointer in common, allocated 
c in getinp), which is plotted from pltloc: 
c
                  if (nhtscale.gt.0) then
                    call alloc(pfhts,nhtscale*10)
                    hts(:) = htutkmx(it,:,l)
                    do i=1,6
                      call intloc(fkmx10(1,i),hts,kmx,
     +                  htscale,nhtscale,0,fhts(nhtscale*(i-1)+1),
     +                  1,1,nhtscale,ier,cpspval,0)
                    enddo
                    call mkionvel(fhts(nhtscale*7+1),nhtscale,glat,glon,
     +                gmlat,gmlon,rdip,rdec,float(date(2)),ut(it),
     +                fhts,fhts(nhtscale+1),fhts(nhtscale*2+1),
     +                fhts(nhtscale*3+1),fhts(nhtscale*4+1),
     +                fhts(nhtscale*5+1),ionvel-1,cpspval)
c
c pointers: htions(ntms,nhtscale,3,nloc) and fhts(nhtscale,10),
c           with ui,vi,wi in fhts(:,8-10)
c
                    do i=1,3
                      do k=1,nhtscale
                        ix0 = (l-1)*3*nhtscale*ntms+(i-1)*nhtscale*
     +                        ntms+(k-1)*ntms+it
                        ix1 = (i+6)*nhtscale+k
                        htions(ix0) = fhts(ix1)
                      enddo
                    enddo 
                    call hpdeallc(pfhts,ier,1)
                  endif
                  didions = .true.
                endif
c
c E5577 (greenline volume emission, see mke5577):
c (note defkmx returns loc, global, or zonal means as needed)
c
              elseif (ip.eq.ixe5577) then
                call defkmx(fields(1,1,1,ifget(ixt)),fkmx10,
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,2),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,3),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,4),
     +            glat,glon,it)
                if (idop5577.gt.0) then ! put un,vn in fkmx10 5,6
                  call defkmx(fields(1,1,1,ifget(ixu)),fkmx10(1,5),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,6),
     +              glat,glon,it)
                  sigma = 17924. ! for doppler from greenline
                endif
                call mke5577(fkmx10,fkmx10(1,2),fkmx10(1,3),
     +            fkmx10(1,4),kmx,fkmx10(1,10),cpspval)
                utkmx(it,:,ifplt(ixe5577),l) = fkmx10(:,10)
                if (nhtscale.gt.0) then
                  call alloc(pfhts,nhtscale*7)
                  hts(:) = htutkmx(it,:,l)
                  do i=1,4 ! interpolate tn,o,o2,n2 to htscale
                    lg = 0
                    if (i.gt.1) lg = 1
                    call intloc(fkmx10(1,i),hts,kmx,htscale,nhtscale,
     +                lg,fhts(nhtscale*(i-1)+1),1,1,nhtscale,ier,
     +                cpspval,0)
                  enddo                  
                  call mke5577(fhts,fhts(nhtscale+1),
     +              fhts(nhtscale*2+1),fhts(nhtscale*3+1),
     +              nhtscale,fhts(nhtscale*6+1),cpspval)
c
c pointers: hte5577(ntms,nhtscale,nloc), fhts(nhtscale,7), with
c           e5577 currently in fhts(:,7) 
c
                  do k=1,nhtscale
                    ix0 = (l-1)*nhtscale*ntms + (k-1)*ntms + it
                    hte5577(ix0) = fhts(nhtscale*6+k)
                  enddo
c
c If doing doppler (i.e., doing e5577, nhtscale > 0),
c   then define 4 doppler fields: tn,un,vn,brightness in dop5577(ntms,nloc,4) 
c   (dop5577 is pdop pointee, pdop5577 is in common) 
c (doppler and quadrat routine in ~/lib/doppler.f)
c     SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma,
c    +  VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ)
c doppler takes ht,t,u,v, and ver as input (..HM,TM,ZM,VM,VER..
c (t,u,v currently in fkmx10, slots 1,5,6; put ver in fkmx10(:,10):
c
                  if (idop5577.gt.0) then ! interpolate u,v -> fhts 5,6
                    call mke5577(fkmx10,fkmx10(1,2),fkmx10(1,3),
     +                fkmx10(1,4),kmx,fkmx10(1,10),cpspval)
                    hts(:) = htutkmx(it,:,l)
                    call doppler(kmx,hts,fkmx10,fkmx10(1,5),
     +                fkmx10(1,6),fkmx10(1,10),sigma,dopvn,dopun,doptn,
     +                dattpv,attpz,dattpz)
                    hts(:) = htutkmx(it,:,l)*1.e5
                    dopb = quadrat(kmx,hts,fkmx10(1,10))*1.e-6
                    do ipp=1,4
                      ix = (ipp-1)*nloc*ntms+(l-1)*ntms+it
                      if (ipp.eq.1) dop5577(ix) = doptn
                      if (ipp.eq.2) dop5577(ix) = dopun
                      if (ipp.eq.3) dop5577(ix) = dopvn
                      if (ipp.eq.4) dop5577(ix) = dopb
                    enddo
                  endif
                  call hpdeallc(pfhts,ier,1)
                endif ! nhtscale > 0
c
c E6300 (redline volume emission, see mke6300):
c (note defkmx returns loc, global, or zonal means as needed)
c
              elseif (ip.eq.ixe6300) then
                call defkmx(fields(1,1,1,ifget(ixt)),fkmx10,
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixte)),fkmx10(1,2),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,3),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,4),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,5),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo2p)),fkmx10(1,6),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixne)),fkmx10(1,7),
     +            glat,glon,it)
                if (idop6300.gt.0) then ! put un,vn in fkmx10 8,9
                  call defkmx(fields(1,1,1,ifget(ixu)),fkmx10(1,8),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,9),
     +              glat,glon,it)
                  sigma = 15867.852 ! for doppler from redline
                endif
c pass in t,te,o2,o,n2,o2+, get out redline ver in slot 10 of fkmx10:
c     subroutine mke6300(tn,te,xo2,xo,xn2,xo2p,xne,id,fout,spval)
                call mke6300(fkmx10,fkmx10(1,2),fkmx10(1,3),
     +            fkmx10(1,4),fkmx10(1,5),fkmx10(1,6),fkmx10(1,7),kmx,
     +            fkmx10(1,10),cpspval)
                utkmx(it,:,ifplt(ixe6300),l) = fkmx10(:,10)
                if (nhtscale.gt.0) then
                  call alloc(pfhts,nhtscale*10)
                  hts(:) = htutkmx(it,:,l)
                  do i=1,7 ! interpolate tn,te,o,o2,n2,o2+,ne to htscale
                    lg = 0
                    if (i.gt.2) lg = 1
                    call intloc(fkmx10(1,i),hts,kmx,htscale,nhtscale,
     +                lg,fhts(nhtscale*(i-1)+1),1,1,nhtscale,ier,
     +                cpspval,0)
                  enddo                  
c ht interpolated ver will be in fhts(:,10):
                  call mke6300(fhts,fhts(nhtscale+1),
     +              fhts(nhtscale*2+1),fhts(nhtscale*3+1),
     +              fhts(nhtscale*4+1),fhts(nhtscale*5+1),
     +              fhts(nhtscale*6+1),
     +              nhtscale,fhts(nhtscale*9+1),cpspval)
c
c pointers: hte6300(ntms,nhtscale,nloc), fhts(nhtscale,10), with
c           e6300 currently in fhts(:,10) 
c
                  do k=1,nhtscale
                    ix0 = (l-1)*nhtscale*ntms + (k-1)*ntms + it
                    hte6300(ix0) = fhts(nhtscale*9+k)
                  enddo
c
c If doing doppler (i.e., doing e6300, nhtscale > 0)
c   then define 4 doppler fields: tn,un,vn,brightness in dop6300(ntms,nloc,4) 
c   (dop6300 is pdop pointee, pdop6300 is in common) 
c (doppler and quadrat routine in ~/lib/doppler.f)
c     SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma,
c    +  VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ)
c doppler takes ht,t,u,v, and ver as input (..HM,TM,ZM,VM,VER..
c (t,u,v currently in fkmx10, slots 1,8,9; put ver in fkmx10(:,10):
c
                  if (idop6300.gt.0) then
                    call mke6300(fkmx10,fkmx10(1,2),fkmx10(1,3),
     +                fkmx10(1,4),fkmx10(1,5),fkmx10(1,6),fkmx10(1,7),
     +                kmx,fkmx10(1,10),cpspval)
                    hts(:) = htutkmx(it,:,l)
                    call doppler(kmx,hts,fkmx10,fkmx10(1,8),
     +                fkmx10(1,9),fkmx10(1,10),sigma,dopvn,dopun,doptn,
     +                dattpv,attpz,dattpz)
                    hts(:) = htutkmx(it,:,l)*1.e5
                    dopb = quadrat(kmx,hts,fkmx10(1,10))*1.e-6
                    do ipp=1,4
                      ix = (ipp-1)*nloc*ntms+(l-1)*ntms+it
                      if (ipp.eq.1) dop6300(ix) = doptn
                      if (ipp.eq.2) dop6300(ix) = dopun
                      if (ipp.eq.3) dop6300(ix) = dopvn
                      if (ipp.eq.4) dop6300(ix) = dopb
                    enddo
                  endif
                  call hpdeallc(pfhts,ier,1)
                endif ! nhtscale > 0
c
c Eo200 (O2 (0-0) atmos band emission, see ~/lib/mkeo200):
c (note defkmx returns loc, global, or zonal means as needed)
c
              elseif (ip.eq.ixeo200) then
                call defkmx(fields(1,1,1,ifget(ixt)),fkmx10,
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,2),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,3),
     +            glat,glon,it)
                call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,4),
     +            glat,glon,it)
                if (idopo200.gt.0) then ! put un,vn in fkmx10 5,6
                  call defkmx(fields(1,1,1,ifget(ixu)),fkmx10(1,5),
     +              glat,glon,it)
                  call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,6),
     +              glat,glon,it)
                  sigma = 17924. ! for doppler from greenline
                endif
                call mkeo200(fkmx10,fkmx10(1,2),fkmx10(1,3),
     +            fkmx10(1,4),kmx,fkmx10(1,10),cpspval)
                utkmx(it,:,ifplt(ixeo200),l) = fkmx10(:,10)
                if (nhtscale.gt.0) then
                  call alloc(pfhts,nhtscale*7)
                  hts(:) = htutkmx(it,:,l)
                  do i=1,4 ! interpolate tn,o,o2,n2 to htscale
                    lg = 0
                    if (i.gt.1) lg = 1
                    call intloc(fkmx10(1,i),hts,kmx,htscale,nhtscale,
     +                lg,fhts(nhtscale*(i-1)+1),1,1,nhtscale,ier,
     +                cpspval,0)
                  enddo                  
                  call mkeo200(fhts,fhts(nhtscale+1),
     +              fhts(nhtscale*2+1),fhts(nhtscale*3+1),
     +              nhtscale,fhts(nhtscale*6+1),cpspval)
c
c pointers: hteo200(ntms,nhtscale,nloc), fhts(nhtscale,7), with
c           eo200 currently in fhts(:,7) 
c
                  do k=1,nhtscale
                    ix0 = (l-1)*nhtscale*ntms + (k-1)*ntms + it
                    hteo200(ix0) = fhts(nhtscale*6+k)
                  enddo
c
c If doing doppler (i.e., doing eo200, nhtscale > 0),
c   then define 4 doppler fields: tn,un,vn,brightness in dopo200(ntms,nloc,4) 
c   (dopo200 is pdop pointee, pdopo200 is in common) 
c (doppler and quadrat routine in ~/lib/doppler.f)
c     SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma,
c    +  VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ)
c doppler takes ht,t,u,v, and ver as input (..HM,TM,ZM,VM,VER..
c (t,u,v currently in fkmx10, slots 1,5,6; put ver in fkmx10(:,10):
c
                  if (idopo200.gt.0) then ! interpolate u,v -> fhts 5,6
                    call mkeo200(fkmx10,fkmx10(1,2),fkmx10(1,3),
     +                fkmx10(1,4),kmx,fkmx10(1,10),cpspval)
                    hts(:) = htutkmx(it,:,l)
                    call doppler(kmx,hts,fkmx10,fkmx10(1,5),
     +                fkmx10(1,6),fkmx10(1,10),sigma,dopvn,dopun,doptn,
     +                dattpv,attpz,dattpz)
                    hts(:) = htutkmx(it,:,l)*1.e5
                    dopb = quadrat(kmx,hts,fkmx10(1,10))*1.e-6
                    do ipp=1,4
                      ix = (ipp-1)*nloc*ntms+(l-1)*ntms+it
                      if (ipp.eq.1) dopo200(ix) = doptn
                      if (ipp.eq.2) dopo200(ix) = dopun
                      if (ipp.eq.3) dopo200(ix) = dopvn
                      if (ipp.eq.4) dopo200(ix) = dopb
                    enddo
                  endif
                  call hpdeallc(pfhts,ier,1)
                endif ! nhtscale > 0
              else
                write(6,"('mkfields: unknown ip=',i3)") ip
              endif
            endif
          endif
 100      continue
          enddo ! fields loop
        enddo
 1000 continue
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
