c
c------------------------------------------------------------------
c
      subroutine pltlon(fields,utc,ihtsc)
c
c Contour longitude slices (lat on x-axis, zp or ht on y-axis,
c   according to ihtsc)
c Also plot with x-axis latitudes over-the-pole, according to
c   ovpollat() if nvpollat > 0
c
      include 'timesproc.h'
      character*56 toplab,fieldlab
      character*80 rec80
      character*32 xlab
      dimension plt(jmx,kmx),vp(4),yaxht(kmx),rimx(imx),
     +  plt0(jmx,kmx),fields(imx,kmx,jmx,nfget),hts(jmx,kmx),
     +  glats(jmx),fjmx0(jmx),fjmx1(jmx)
      pointer(ppltht,pltht(1)), (pplthto,plthto(1)),
     +  (ppmbht,pmbht(1))
      common/islat/ iplat
      data vp /.13,.87,.26,.91/, iplat/0/
c
c For pressure (mb) right hand y-axis: note strings pfcode (with
c function codes) must match values of presslab (see yaxright.f)
c
      parameter(npresslab=10,nhtlab=6)
      dimension presslab(npresslab),htlab(nhtlab)
      character*8 pfcode(npresslab)
      data presslab
     +  / 1.e+1,     1.e+0,     1.e-1,     1.e-2,     1.e-3, 
     +    1.e-4,     1.e-5,     1.e-6,     1.e-7,     1.e-8/
      data pfcode  
     +  /'10$S$1  ','10$S$0  ','10$S2$-1','10$S2$-2','10$S2$-3',
     +   '10$S2$-4','10$S2$-5','10$S2$-6','10$S2$-7','10$S2$-8'/
c
c Begin execution:
c
      if (nlon+nslt.le.0) then
        write(6,"('pltlon: no selected longitudes or local times',
     +    ' -- returning')")
        return
      endif
      if (iyax_r.gt.1.and.ihtsc.le.0) then
        vp(1) = .1
        vp(2) = .75
      else
        vp(1) = .13
        vp(2) = .87
      endif
c
c Set up conpack for y-axis:
c
      xmid = 0.5*(vp(1)+vp(2))
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('ILX',xmid)
      call cpseti('ILP',0)
      call cpsetr('ILS',.016)
      if (labels.le.1) call cpsetr('ILS',.02)
      if (ihtsc.le.0) then
        call cpsetr('YC1',gcmzp(izprange(1)))
        call cpsetr('YCN',gcmzp(izprange(2)))
      else
        call cpsetr('YC1',htscale(1))
        call cpsetr('YCN',htscale(nhtscale))
        call hpalloc(ppltht,jmx*nhtscale,ier,1)
        if (ier.ne.0) then
          write(6,"('pltlon: hpalloc error = ',i6,' for pltht')") ier
          stop 'ppltht'
        endif
        if (nvpollat.gt.0) then
          call hpalloc(pplthto,jmx*nhtscale,ier,1)
          if (ier.ne.0) then
            write(6,"('pltlon: hpalloc error = ',i6,' for plthto')") ier
            stop 'pplthto'
          endif
        endif
      endif
c
c Field loop: 
c
      do 200 ip=1,nftot
        if (ifplt(ip).le.0.or.(ihtsc.gt.0.and.ip.eq.itxz)) goto 200
        ilog = 0
        if (itimelog(ip).gt.0.and.logvert.gt.0) ilog = 1
        nplat = 0
        iplat = 0
c
c after 105, nplat may be > 0 (over-the-pole latitude on x-axis):
c Set up x-axis accordingly:
c
 105    continue
        if (nplat.le.0) then
          call cpsetr('XC1',gcmlat(1))
          call cpsetr('XCM',gcmlat(jmx))
          if (ihtsc.le.0) then
            call set(vp(1),vp(2),vp(3),vp(4),
     +               -90.,90.,gcmzp(izprange(1)),gcmzp(izprange(2)),1)
          else
            call set(vp(1),vp(2),vp(3),vp(4),
     +               -90.,90.,htscale(1),htscale(nhtscale),1)
          endif
        endif
c
c Selected longitude loop:
c
        do 100 i=1,nlon+nslt
          call cpsetr('ILY',-.16)
          if (i.le.nlon) then
            if (slon(i).ne.zmflag) then
              ixlon = ixfind(gcmlon,imx,slon(i),dlon)
              if (ixlon.le.0) then
                write(6,"('>>> pltlon: bad longitude=',f10.3,' i=',i3,
     +            ' skipping this lon')") slon(i),i
                goto 100
              endif
              rlon = gcmlon(ixlon)
              if (nplat.gt.0) then
                iplat = 1
                rlono = rlon+180.
                if (rlono.ge.180.) rlono = rlono-360.
                ixlono = ixfind(gcmlon,imx,rlono,dlon)
                if (ixlono.le.0) then
                  write(6,"('>>> pltlon: bad opposite longitude=',
     +              f10.3,' i=',i3,' skipping this lon')") rlono,i
                  goto 100
                endif
                rlono = gcmlon(ixlono)
              endif
            else ! zonal means
              ixlon = ifix(zmflag)
              rlon = slon(i)
              ixlono = ixlon
              rlono = rlon
            endif
c
c Solar local time:
          else
            islt = i-nlon
            ixlon = ixslt(sslt(islt),utc,rlon,gcmlon,imx,dlon)
            if (nplat.gt.0) then
              iplat = 1
              rlono = rlon+180.
              if (rlono.ge.180.) rlono = rlono-360.
              ixlono = ixfind(gcmlon,imx,rlono,dlon)
              if (ixlono.le.0) then
                write(6,"('>>> pltlon: bad opposite longitude=',
     +            f10.3,' i=',i3,' skipping this lon')") rlono,i
                goto 100
              endif
              rlono = gcmlon(ixlono) ! opposite longitude
              slto = fslt(dum,utc,rlono,1) ! opposite local time 
            endif
          endif
c
c Get field along desired longitude (also along opposite longitde
c   if doing over-the-pole x-axis)
c Note getlon() handles zonal mean or selected longitude, and
c   log10 or no log10 of field.
c
          if (ihtsc.le.0) then
            call getlon(fields,plt0,ixlon,ilog,ip)
            do k=izprange(1),izprange(2)
              kk = k-izprange(1)+1
              plt(:,kk) = plt0(:,k)
            enddo
c
c Store opposite longitude in plt0 for use when nplat > 0 (use hts for temp)
            if (nplat.gt.0.and.rlon.ne.zmflag) then
              call getlon(fields,hts,ixlono,ilog,ip)
              do k=izprange(1),izprange(2)
                kk = k-izprange(1)+1
                plt0(:,kk) = hts(:,k)
              enddo
            endif
            ny = izprange(2)-izprange(1)+1
c
c Interpolate to linear height scale (take logs after interpolation):
c (log flag in interp routine is for linear or log interpolation)
          else 
            if (ip.lt.ixe5577.or.ip.gt.ixeo200) then
              call getlon(fields,plt0,ixlon,0,ip)
              call getlon(fields,hts,ixlon,0,itxz)
              call cuthtint(plt0,hts,jmx,kmx,pltht,htscale,nhtscale,
     +          ilog,cpspval,ier,1)
              if (ilog.gt.0)
     +          call log10f(pltht,jmx*nhtscale,1.e-20,cpspval)
c
c Interpolate along opposite longitude for over-the-pole x-axis (plt0):
              if (nplat.gt.0.and.rlon.ne.zmflag) then
                call getlon(fields,plt0,ixlono,0,ip)
                call cuthtint(plt0,hts,jmx,kmx,plthto,htscale,nhtscale,
     +            ilog,cpspval,ier,1)
                if (ilog.gt.0) call log10f(plthto,jmx*nhtscale,
     +            1.e-20,cpspval)
              endif
            else ! get interpolated emissions
              call getelonh(fields,pltht,ixlon,ip)
              if (ilog.gt.0)
     +          call log10f(pltht,jmx*nhtscale,1.e-20,cpspval)
              if (nplat.gt.0.and.rlon.ne.zmflag) then
                call getelonh(fields,plthto,ixlono,ip)
                if (ilog.gt.0) call log10f(plthto,jmx*nhtscale,
     +            1.e-20,cpspval)
              endif 
            endif
            ny = nhtscale 
          endif
c
c Specified lat for x-axis over the pole: 
c ovpollat(nvpollat) = lats at which to plot x-axis over the pole (in common)
c mkpollat changes plt to be field at new latitudes (ovpollat->pole->ovpollat)
c   (new latitudes are returned in glats(nlats))
c subroutine mkpollat(plt,jmx,ny,pollat,gcmlat,dlat,glats,nlats)
c Also do new set call and inform conpack re new x-axis
c
          nx = jmx
          glats(:) = gcmlat(:)
          write(xlab,"('LATITUDE')")
          if (nplat.gt.0) then
            izm = 0
            if (ixlono.eq.ifix(zmflag)) izm = 1
            if (ihtsc.le.0) then
              call mkpollat(plt,plt0,jmx,ny,ovpollat(nplat),gcmlat,
     +          dlat,glats,nx,izm)
            else
              call mkpollat(pltht,plthto,jmx,ny,ovpollat(nplat),gcmlat,
     +          dlat,glats,nx,izm)
            endif
            if (nx.le.0) then
              write(6,"('Error from mkpollat -- skipping')")
              goto 101
            endif
            call cpsetr('XC1',glats(1))
            call cpsetr('XCM',glats(nx))
            if (ihtsc.le.0) then
              call set(vp(1),vp(2),vp(3),vp(4),glats(1),glats(nx),
     +          gcmzp(izprange(1)),gcmzp(izprange(2)),1)
            else
              call set(vp(1),vp(2),vp(3),vp(4),glats(1),glats(nx),
     +          htscale(1),htscale(nhtscale),1)
            endif
            write(xlab,"('LATITUDE (',f5.1,' DEG AT ENDS)')") glats(1)
            if (izm.le.0) call cpsetr('ILY',-.30)
          endif
c
c Contour:
c
          if (labels.gt.1) then
            if (icolor.le.0) then
              if (ihtsc.le.0) then
                call contour(plt,jmx,nx,ny,cint(ip),cmin(ip),cmax(ip))
              else
                call contour(pltht,jmx,nx,ny,cint(ip),cmin(ip),cmax(ip))
              endif
            else
              if (ihtsc.le.0) then
                call conclr(plt,jmx,nx,ny,cint(ip),cmin(ip),cmax(ip))
              else
                call conclr(pltht,jmx,nx,ny,cint(ip),cmin(ip),cmax(ip))
              endif
            endif
            lenxlab = lenstr(xlab)
            if (ihtsc.le.0) then
c
c Temporary: try pressure (mb) on y-axis:
c
              call labrect(glats,nx,gcmzp(izprange(1)),ny,
     +          xlab(1:lenxlab),'ZP',0.)
            else
              call labrect(glats,nx,htscale(1),ny,xlab(1:lenxlab),
     +          'HEIGHT (KM)',0.)
            endif
c
c Simple contour and labels:
          else
            if (ihtsc.le.0) then
              call cpcnrc(plt,jmx,nx,ny,0.,0.,finc,1,-1,-1634B)
              call labrect(glats,nx,gcmzp(izprange(1)),ny,'LAT',
     +          'ZP',0.)
            else
              call cpcnrc(pltht,jmx,nx,ny,0.,0.,finc,1,-1,-1634B)
              call labrect(glats,nx,htscale,ny,'LAT','HT',0.)
            endif
          endif
c
c Top label with field name, lon/slt, ut, etc:
c
          call clearstr(toplab)
          if (i.le.nlon) then ! selected lon
            if (slon(i).ne.zmflag) then ! is not zonal means
              if (ilog.le.0) then
                if (nplat.le.0) then
                  write(toplab,"(a8,' LON=',f6.1,' UT=',f6.2)")
     +              flab(ip),rlon,utc
                else
                  write(toplab,"(a8,' LON=(',f6.1,',',f6.1,') UT=',
     +              f6.2)") flab(ip),rlon,rlono,utc
                endif
              else
                if (nplat.le.0) then
                  write(toplab,"('LOG10 ',a8,' LON=',f6.1,' UT=',f6.2)")
     +              flab(ip),rlon,utc
                else
                  write(toplab,"('LOG10 ',a8,' LON=(',f6.1,',',f6.1,
     +              ') UT=',f6.2)") 
     +            flab(ip),rlon,rlono,utc
                endif
              endif
            else ! is zonal means
              if (ilog.le.0) then
                write(toplab,"(a8,' ZONAL MEANS  UT=',f6.2)")
     +            flab(ip),utc
              else
                write(toplab,"('LOG10 ',a8,' ZONAL MEANS  UT=',f6.2)")
     +            flab(ip),utc
              endif
            endif
          else ! selected slt
            if (ilog.le.0) then
              if (nplat.le.0) then
                write(toplab,"(a8,' SLT=',f5.1,' UT=',f6.2)")
     +            flab(ip),sslt(islt),utc
              else
                write(toplab,"(a8,' SLT=(',f5.1,',',f5.1,') UT=',f6.2)") 
     +            flab(ip),sslt(islt),slto,utc
              endif
            else
              if (nplat.le.0) then
                write(toplab,"('LOG10 ',a8,' SLT=',f5.1,' UT=',f6.2)")
     +            flab(ip),sslt(islt),utc
              else
                write(toplab,"('LOG10 ',a8,' SLT=(',f5.1,',',f5.1,
     +            ') UT=',f6.2)")
     +            flab(ip),sslt(islt),slto,utc
              endif
            endif
          endif
          call clearstr(rec80)
          write(rec80,"(a)") toplab(1:lenstr(toplab))
          call clearstr(fieldlab)
          write(fieldlab,"(a)") flab(ip)
          if (labels.gt.1) then
            call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+0.05,
     +        0.)
            write(toplab,"('HISTORY=',a)") histvol(ivol)
            call wrlab(toplab(1:lenstr(toplab)),xmid,vp(3)-0.15,.012)
          else
            call wrlablq(toplab(1:lenstr(toplab)),xmid,vp(4)+0.05,.02)
          endif
c
c Put height axis on right (use latitudinally averaged heights):
c (if doing over-the-pole lat x-axis, then heights are averaged from
c  lats across the pole):
c
          yaxht(:) = 0.
          if (ihtsc.le.0) then
            izp0 = izprange(1)
            izp1 = izprange(2)
          else
            izp0 = 1
            izp1 = kmx
          endif 
          do k=izp0,izp1
            izp = k-izp0+1
            if (slon(i).ne.zmflag) then
              if (nplat.le.0) then
                do j=1,jmx
                  yaxht(izp) = yaxht(izp)+fields(ixlon,k,j,ifget(itxz))
                enddo
              else
                fjmx0(:) = fields(ixlon,k,:,ifget(itxz))
                fjmx1(:) = fields(ixlono,k,:,ifget(itxz))
                call mkpollat(fjmx0,fjmx1,jmx,1,ovpollat(nplat),
     +            gcmlat,dlat,dum,999,0)
                do j=1,nx
                  yaxht(izp) = yaxht(izp) + fjmx0(j)
                enddo
              endif
            else		! zonal means -- use global mean hts
              if (nplat.le.0) then
                do j=1,jmx
	          rimx(:) = fields(:,k,j,ifget(itxz))
                  yaxht(izp) = yaxht(izp) + 
     +              calcmean(rimx,imx,0,1.e-20,cpspval)
                enddo
              else
                do j=1,jmx
	          rimx(:) = fields(:,k,j,ifget(itxz))
                  fjmx0(j) = calcmean(rimx,imx,0,1.e-20,cpspval)
                enddo
                call mkpollat(fjmx0,dum,jmx,1,ovpollat(nplat),
     +            gcmlat,dlat,dum,999,1)
                do j=1,nx
                  yaxht(izp) = yaxht(izp) + fjmx0(j)
                enddo
              endif
            endif
            yaxht(izp) = yaxht(izp) / float(nx)
          enddo 
c
c If doing zp on left axis, add right axis in height and/or pmb if desired:
c
          if (ihtsc.le.0.and.iyax_r.ge.1) then
            rndval = 10.
            if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rndval = 5.
            botrnd = rnd(yaxht(1),rndval)
            toprnd = rnd(yaxht(ny),rndval)
            del = (toprnd-botrnd)/nhtlab
            delrnd = rnd(del,rndval)
            do k=1,nhtlab
              htlab(k) = botrnd + (k-1)*delrnd
            enddo
            call yaxright(gcmzp(izprange(1)),yaxht,ny,
     +        htlab,nhtlab,'(I3)',1,pfcode,'HEIGHT (km)',0.,0.)
            if (iyax_r.eq.2) then
              call yaxright(gcmzp(izprange(1)),pmb(izprange(1)),ny,
     +          presslab,npresslab,' ',0,pfcode,'PRESSURE (mb)',.13,0.)
            endif
          endif
c
c If left axis is linear height, and we are plotting a right hand axis
c in press (mb), then we need pmbht(nhtscale), i.e., pressure in mb
c at each height in htscale(nhtscale)
c     subroutine vecterp(fy,f,ny,fylin,flin,nlin)
c
          if (ihtsc.gt.0.and.iyax_r.ge.1) then
            call hpalloc(ppmbht,nhtscale,ier,1)
            call vecterp(yaxht,gcmzp,kmx,htscale,pmbht,nhtscale,
     +        cpspval,0)
            do k=1,nhtscale
              if (pmbht(k).ne.cpspval)
     +          pmbht(k) = p0 * exp(-pmbht(k)) ! pressure(mb) at each htscale
            enddo
            call yaxright(htscale,pmbht,ny,presslab,npresslab,' ',0,
     +        pfcode,'PRESSURE (mb)',0.,0.)
          endif
c
c If doing over the pole x-axis, show extra x-axis with the 2 lons,
c or both lons and slts if doing selected slt (no extra x-axis if
c doing zonal means):
c
          if (nplat.gt.0.and.slon(i).ne.zmflag) then 
            if (i.le.nlon) then
              call lab2lon(rlon,rlono,utc,.12,0.,1)
            else
              call lab2lon(rlon,rlono,utc,.12,0.,3)
            endif
          endif
c
c Wrap it up, including ascii data file, if needed:
c
          call frame
          iframe = iframe+1
c
c Report to stdout:
c Pressure on y-axis:
          if (ihtsc.le.0) then
            if (iwrascii.gt.0) then
              call wrascii(iwrascii,luascii,plt,nx,ny,glats,
     +          gcmzp(izprange(1)),'LATITUDE','LN(P0/P)',
     +          histvol(ivol),fieldlab,iframe,rec80,'timesproc',
     +          dirascii)
            endif
            if (i.le.nlon) then
              if (slon(i).ne.zmflag) then
                if (nplat.le.0) then
                  write(6,"('pltlon: frame ',i4,' field ',a8,' lon=',
     +              f8.2,' zp=',f5.1,',',f5.1)") 
     +            iframe,flab(ip),rlon,zprange
                else
                  write(6,"('pltlon: frame ',i4,' field ',a8,' lon=',
     +              f6.1,',',f6.1,' (ovplat=',f5.1,') zp=',
     +              f5.1,',',f5.1)") 
     +            iframe,flab(ip),rlon,rlono,glats(1),zprange
                endif
              else ! zonal means
                if (nplat.le.0) then
                  write(6,"('pltlon: frame ',i4,' field ',a8,
     +              ' (zonal means) zp=',f5.1,',',f5.1)") 
     +              iframe,flab(ip),zprange
                else
                  write(6,"('pltlon: frame ',i4,' field ',a8,
     +              ' (zonal means) zp=',f5.1,',',f5.1,
     +              ' ovplat=',f5.1)") 
     +              iframe,flab(ip),zprange,glats(1)
                endif
              endif
            else ! local time
              if (nplat.le.0) then
                write(6,"('pltlon: frame ',i4,' field ',a8,' slt=',f5.1,
     +            ' zp=',f5.1,',',f5.1)") iframe,flab(ip),
     +            sslt(islt),zprange
              else
                write(6,"('pltlon: frame ',i4,' field ',a8,' slt=',f5.1,
     +            ',',f5.1,' (ovplat=',f5.1,') zp=',f5.1,',',f5.1)") 
     +            iframe,flab(ip),sslt(islt),slto,glats(1),
     +            zprange
              endif
            endif
c
c Height on y-axis:
          else
            if (iwrascii.gt.0) then
              call wrascii(iwrascii,luascii,pltht,nx,ny,glats,
     +          htscale,'LATITUDE','HEIGHT (KM)',histvol(ivol),fieldlab,
     +            iframe,rec80,'timesproc',dirascii)
            endif
            if (i.le.nlon) then
              if (slon(i).ne.zmflag) then
                if (nplat.le.0) then
                  write(6,"('pltlon: frame ',i4,' field ',a8,' lon=',
     +              f6.1,' ht=',f7.2,',',f7.2)") iframe,
     +              flab(ip),rlon,htscale(1),htscale(nhtscale)
                else
                  write(6,"('pltlon: frame ',i4,' field ',a8,' ovplat=',
     +              f5.1,' ht=',f7.2,',',f7.2)") iframe,
     +              flab(ip),glats(1),htscale(1),
     +              htscale(nhtscale)
                endif
              else
                write(6,"('pltlon: frame ',i4,' field ',a8,
     +            ' (zonal means) ht=',f6.1,',',f6.1)")
     +            iframe,flab(ip),htscale(1),htscale(nhtscale)
              endif
            else ! local time
              if (nplat.le.0) then
                write(6,"('pltlon: frame ',i4,' field ',a8,' slt=',f5.1,
     +            ' ht=',f6.1,',',f6.1)") iframe,
     +            flab(ip),sslt(islt),htscale(1),
     +            htscale(nhtscale)
              else
                write(6,"('pltlon: frame ',i4,' field ',a8,' slt=',f5.1,
     +            ',',f5.1,' (ovplat=',f5.1,') ht=',2f6.1)") iframe,
     +            flab(ip),sslt(islt),slto,glats(1),htscale(1),
     +            htscale(nhtscale)
              endif
            endif
          endif
c
c End selected longitude/slt loop
 100    continue
c
c Redo with "over the pole" lat x-axis:
 101    continue
        if (nvpollat.gt.0) then
          nplat = nplat+1
          if (nplat.le.nvpollat) goto 105
        endif
c
c End field loop: ip=1,nftot
 200  continue
c
      if (ihtsc.gt.0) call hpdeallc(ppltht,ier,0)
      if (ihtsc.gt.0.and.nvpollat.gt.0) call hpdeallc(pplthto,ier,0)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mkpollat(plt,plto,jmx,ny,plat,gcmlat,dlat,
     +  glats,nlats,izm)
c
c Change plt to contain field at lats plat->pole->plat:
c On input: 
c   plt(jmx,ny) contains all latitudes at selected lon (slon)
c   plto(jmx,ny) contains all latitudes at lon opposite to slon
c   plat = "perimeter" lat
c   gcmlat(jmx) = model latitudes, dlat = delta lat grid
c   izm=0/1 for isnot/is a zonal means case (if is a zonal means case,
c     then field will be symmetrical across the pole, plto not used)
c On output: 
c   plt(jmx,ny) is changed to contain field over the pole
c   if nlats.ne.999 on input, then glats(nlats) will contains lats 
c     plat->pol->plat increasing over the poles 
c     (i.e., will be > 90 or < -90 on far side of pole)
c   nlats = number of latitudes in glats (if not 999 on input)
c
      dimension plt(jmx,ny),plto(jmx,ny),fjmx(jmx),gcmlat(jmx),
     +  glats(jmx)
c
      ixplat = ixfind(gcmlat,jmx,plat,dlat)
      if (ixplat.le.0) then
        write(6,"('>>> mkpollat: bad plat=',f9.2)") plat
        nlats = 0
        return
      endif 
      inc = 1
      ipole = jmx
      if (plat.lt.0.) then ! south pole
        inc = -1
        ipole = 1
      endif
      do k=1,ny
        idolat = 0
        if (k.eq.1.and.nlats.ne.999) idolat = 1
        jj = 0
        fjmx(:) = plt(:,k)
c
c plat to pole along selected lon:
        do j=ixplat,ipole,inc
          jj = jj+1
          if (idolat.gt.0) glats(jj) = gcmlat(j)
          plt(jj,k) = fjmx(j) 
        enddo
c
c pole to plat along opposite lon (or if zonal means, take from plt):
        if (izm.le.0) then
          do j=ipole,ixplat,-inc
            jj = jj+1
            if (idolat.gt.0) glats(jj) = gcmlat(j)
            plt(jj,k) = plto(j,k)
          enddo
        else
          do j=ipole,ixplat,-inc
            jj = jj+1
            if (idolat.gt.0) glats(jj) = gcmlat(j)
            plt(jj,k) = fjmx(j)
          enddo
        endif
        if (idolat.gt.0) nlats = jj
      enddo
c
c Make latitudes increasing across the poles:
      if (nlats.ne.999) then
        do j=1,nlats
          if (j.gt.1) then
            if (inc.gt.0.and.glats(j).le.glats(j-1).or.
     +          inc.lt.0.and.glats(j).ge.glats(j-1)) then
              glats(j)=(180.-abs(glats(j)))*inc
            endif
          endif
        enddo
      endif
      return
      end 
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine lab2lon(rlon0,rlon1,ut,yoff,charsize,iflag)
      character*24 lab0,lab1
c
c Input:
c   rlon0,rlon1 = longitudes for left half and right half of current frame
c   ut = current ut (decimal hours)
c   charsize = desired character size
c   iflag = 0 -> no labels on the axis
c   iflag = 1 -> label left half of axis with lon rlon0, right half rlon1
c   iflag = 2 -> label left half with slt for rlon0, right half with slt
c              for rlon1
c   iflag = 3 -> label with both lons and slt on both sides of axis
c Output:
c   Draw extra x-axis below existing x-axis, with tics at each end, 
c     and one in the middle (yoff = y-offset below x-axis), labelled 
c     as per iflag
c   Inputs are unchanged
c
      chsize = 0.012
      if (charsize.ne.0.) chsize=charsize
      call getset(vl,vr,vb,vt,wl,wr,wb,wt,ll)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1) 
      xhalf = 0.5*(vl+vr)
      call line(vl,vb-yoff,vr,vb-yoff)
      call line(vl,vb-yoff+.015,vl,vb-yoff-.015) ! left tic
      call line(xhalf,vb-yoff+.015,xhalf,vb-yoff-.015) ! middle tic
      call line(vr,vb-yoff+.015,vr,vb-yoff-.015) ! right tic
      if (iflag.gt.0) then
        call clearstr(lab0)
        call clearstr(lab1)
        if (iflag.eq.1) then
          write(lab0,"('LON=',f6.1)") rlon0
          write(lab1,"('LON=',f6.1)") rlon1
        else
          slt0 = fslt(dum,ut,rlon0,1)
          slt1 = fslt(dum,ut,rlon1,1)
          if (iflag.eq.2) then
            write(lab0,"('SLT=',f4.1,' (HRS)')") slt0
            write(lab1,"('SLT=',f4.1,' (HRS)')") slt1
          elseif (iflag.eq.3) then
            write(lab0,"('LON=',f6.1,' SLT=',f4.1)") rlon0,slt0
            write(lab1,"('LON=',f6.1,' SLT=',f4.1)") rlon1,slt1
          endif
        endif
        call plchhq(vl+.25*(vr-vl),vb-yoff-.02,lab0(1:lenstr(lab0)),
     +    chsize,0.,0.)
        call plchhq(vl+.75*(vr-vl),vb-yoff-.02,lab1(1:lenstr(lab1)),
     +    chsize,0.,0.)
      endif
      call set(vl,vr,vb,vt,wl,wr,wb,wt,ll)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getelonh(fields,pltht,ixlon,ip)
c
c Return volume emission in pltht, interpolated to linear height scale
c Interpolate individual required fields before calculating emission
c
      include 'timesproc.h'
      dimension pltht(jmx,nhtscale),fields(imx,kmx,jmx,nfget),
     +  rkmx(kmx),hts(kmx)
      pointer(pfht,fht(1))
c
      if (ip.eq.ixe5577) then ! need t,o2,o,n2
        call hpalloc(pfht,nhtscale*5,ier,1)
        do j=1,jmx
          call getloc(fields,hts,j,ixlon,0,itxz)
          call getloc(fields,rkmx,j,ixlon,0,itxt)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,0,fht,1,1,
     +      nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo2)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo1)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*2+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxn2)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*3+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call mke5577(fht,fht(nhtscale+1),fht(nhtscale*2+1),
     +      fht(nhtscale*3+1),nhtscale,fht(nhtscale*4+1),cpspval)
          do k=1,nhtscale
            pltht(j,k) = fht(nhtscale*4+k)
          enddo
        enddo
        call hpdeallc(pfht,ier,1)
      elseif (ip.eq.ixe6300) then
        call hpalloc(pfht,nhtscale*8,ier,1)
        do j=1,jmx
          call getloc(fields,hts,j,ixlon,0,itxz)
          call getloc(fields,rkmx,j,ixlon,0,itxt)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,0,fht,1,1,
     +      nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxte)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,0,fht(nhtscale+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo2)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*2+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo1)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*3+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxn2)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*4+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo2p)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*5+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxne)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*6+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call mke6300(fht,fht(nhtscale+1),fht(nhtscale*2+1),
     +      fht(nhtscale*3+1),fht(nhtscale*4+1),fht(nhtscale*5+1),
     +      fht(nhtscale*6+1),nhtscale,fht(nhtscale*7+1),
     +      cpspval)
          do k=1,nhtscale
            pltht(j,k) = fht(nhtscale*7+k)
          enddo
        enddo
        call hpdeallc(pfht,ier,1)
      elseif (ip.eq.ixeo200) then
        call hpalloc(pfht,nhtscale*5,ier,1)
        do j=1,jmx
          call getloc(fields,hts,j,ixlon,0,itxz)
          call getloc(fields,rkmx,j,ixlon,0,itxt)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,0,fht,1,1,
     +      nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo2)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxo1)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*2+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call getloc(fields,rkmx,j,ixlon,0,itxn2)
          call intloc(rkmx,hts,kmx,htscale,nhtscale,1,fht(nhtscale*3+1),
     +      1,1,nhtscale,ier,cpspval,0)
          call mkeo200(fht,fht(nhtscale+1),fht(nhtscale*2+1),
     +      fht(nhtscale*3+1),nhtscale,fht(nhtscale*4+1),cpspval)
          do k=1,nhtscale
            pltht(j,k) = fht(nhtscale*4+k)
          enddo
        enddo
        call hpdeallc(pfht,ier,1)
      else
        write(6,"('getelonh: unknown ip=',i3)") ip
      endif
      return
      end