c
      subroutine mklons(mtime)
c
c Prepare longitude slices and pass to pltlon where contours are drawn:
c fccm(ngcmlon,ngcmlat,nccmlev,nccmfld) (pointer pfccm is in common)
c fhist(imx,kmx,jmx,nfhist) (pointer pfhist is in common)
c
      include "flxproc.h"
      include "ccm.h"
      include "pltopt.h"
      pointer (pplt,plt(1)), (pzpcol,zpcol(1)), (pht2d,ht2d(1)),
     +  (pzp,zp(1)), (pf2d,f2d(1)), (phtscale,htscale(1))
      character*96 tlabs(3),blabs(3)
      character*80 msgout
      real slons(imx),sslts(imx),vplon(4),zpranges(2,mxvrange),
     +  blabchsz(3),tlabchsz(3),vp(4)
      real htscales(3,mxvrange)
      integer mtime(3)
      data vplon /.13,.87,.23,.88/
      data tlabchsz/.020,.020,.020/, blabchsz/.02,.02,.02/
      data toffset /.03/, boffset/.10/
c
      write(6,"(/'Longitude slices at ut = ',f4.1,'  mtime = ',3i3,
     +  ' ccmday = ',f10.4)") ut,(mtime(i),i=1,3),ccmday
      if (multiadvfr.gt.0) nppf = 0
c
c Determine number of selected lons, and place in local array slons:
c
      nlons = 0
      nslts = 0
      do i=1,imx
        if (flons(i).ne.spval) then
          nlons = nlons+1
          slons(nlons) = flons(i)
        endif  
        if (fslts(i).ne.spval) then
          nslts = nslts+1
          sslts(nslts) = fslts(i)
        endif
      enddo
      if (nlons.le.0.and.nslts.le.0) then
        write(6,"('mklons: no lons or slts to plot -- returning')")
        return
      endif
c
c Field loop:
c (non-zero values of ifproc are numbered in user defined order)
c (common nfproc is number of fields to be plotted (max of ifproc values),
c  and nplt is current field in user defined order)
c
      nplt = 1
 101  continue
      do ip=1,mxfproc
        if (ifproc(ip).ne.nplt) goto 100
        if (ip.eq.ixunvn.or.ip.eq.ixuivi) goto 100
        if (ihtindep(ip).gt.0) goto 100
        ilog = 0
        if (ilon_log10.gt.0.and.logplt(ip).gt.0) ilog = 1
c
c Define plot options (common in pltopt.h, included by pltlon)
c
        pltmin = cmin(ip) 
        pltmax = cmax(ip) 
        conint = cint(ip) 
        scfac = scalefac(ip)
        labelv = ivec_label
	pltchsize = 0.
c
c Set up coupled vs tgcm-only field:
c (zp(nzp) = full zp range available, zpcol(ny) = range actually plotted)
c (also allocate heights array ht2d)
c
        if (ifcpl(ip).gt.0) then
          nzp = ncplzp
          call alloc(pzp,nzp)
          do k=1,nzp
            zp(k) = cplzp(k)
          enddo
        else
          nzp = kmx
          call alloc(pzp,nzp)
          do k=1,nzp
            zp(k) = gcmzp(k)
          enddo
        endif
        call alloc(pht2d,jmx*nzp)	! heights for yaxzpht
c
c Determine number of and validate zp ranges requested and transfer 
c   to local array. This must be done in field loop, since valid zpranges
c   are different for coupled vs non-coupled fields.
c
        call setyverts(flon_zprange,zpranges,nzpranges,
     +                 flon_htscale,htscales,nhtscales,
     +                 mxvrange,zp,nzp,dzp,spval,0)
        if (nzpranges.le.0.and.nhtscales.le.0) goto 100
c
c Lon(or zm)/slt loop:
c
        do ilonslt = 1,nlons+nslts
c
c At selected longitudes or zonal means:
c
          if (ilonslt.le.nlons) then
            glon = slons(ilonslt)
            slt = spval
            if (glon.ne.zmflag) slt = fslt(dum,ut,glon,1)
c
c At selected local times:
c
          else
            glon = fslt(sslts(ilonslt-nlons),ut,dum,3)
            slt = sslts(ilonslt-nlons)
          endif
c
c Get ht2d(jmx,nzp) for right hand yaxis (yaxzpht):
c
          call getfldlon(ifcpl(ip),ht2d,imx,jmx,1,nzp,nzp,mxfproc,
     +      glon,logplt(ixz),ixz,fhist,kmx,nfhist,gcmzp,ixfhist,
     +      fccm,nccmlev,nccmfld,cplzp,ixfccm)
c
c Loop through specified zp ranges:
c
          if (nzpranges.le.0) goto 50
          if (iyaxright.gt.1) then ! allow room on right for 2 right y-axes
            vp(1) = .12
            vp(2) = .77
          endif
          do irange = 1,nzpranges
c
c Is a coupled field:
c
            izp0 = ixfind(zp,nzp,zpranges(1,irange),dzp)
            izp1 = ixfind(zp,nzp,zpranges(2,irange),dzp)
            ny = izp1-izp0+1
            call alloc(pplt,jmx*ny)
            call alloc(pzpcol,ny)
            do k=1,ny
              zpcol(k) = zp(izp0+k-1)
            enddo
c
c Transfer to local array and take log10:
c
            call getfldlon(ifcpl(ip),plt,imx,jmx,izp0,izp1,ny,mxfproc,
     +        glon,logplt(ip),ip,fhist,kmx,nfhist,gcmzp,ixfhist,
     +        fccm,nccmlev,nccmfld,cplzp,ixfccm)
            if (ilog.gt.0) call log10f(plt,jmx*ny,1.e-20,spval)
c
c Contour:
c
            if (iyaxright.gt.1) then
              if (multiplt.le.0) vplon(1) = .10
              vplon(2) = .76
            endif
            call setmultivp(vplon,iadvfr,nppf+1,multiplt,ipltrowcol,
     +        vp)
            if (mkplt.gt.0)
     +        call pltlon(plt,ht2d,zp,gcmlat,zpcol,jmx,ny,nzp,glon,
     +          iyaxright,p0,ut,vp,spval)
c
c Info labels, advance frame, report to stdout:
c
            call mklonlab(ifcpl(ip),glon,zmin,zmax,ciu,ilog,mtime,
     +        ip,blabs,tlabs)
            if (mkplt.gt.0)
     +        call wrlab6(tlabs,toffset,tlabchsz,
     +                    blabs,boffset,blabchsz,vp,0)
            nppf = nppf+1               ! number of plots per frame
            call clearstr(msgout)
            if (glon.ne.zmflag) then
              write(msgout,"('field ',a,' lon=',f7.2,
     +          ' zp=',f6.1,' to ',f6.1)") flab8(ip),glon,
     +          zpcol(1),zpcol(ny)
            else
              write(msgout,"('field ',a,
     +          ' (zonal means)  zp=',f6.1,' to ',f6.1)") 
     +          flab8(ip),zpcol(1),zpcol(ny)
            endif
c
c Write to output xdr data file:
c
            if (iwrxdr.gt.0) then
              iframe_xdr = iframe_xdr+1
!             call wrxdr(flnm_xdr(1:lenstr(flnm_xdr)),plt,jmx,ny,
!    +          gcmlat,zpcol,'LATITUDE (DEG)','ZP (LN(P0/P))',
!    +          tlabs(2),tlabs(1),blabs(2),blabs(1))
              call wrxdr(trim(flnm_xdr),plt,jmx,ny,
     +          gcmlat,zpcol,'LATITUDE (DEG)','ZP (LN(P0/P))',
     +          tlabs(2),tlabs(1),blabs(2),blabs(1),mtime,0)
            endif       ! iwrxdr > 0
c                
c Advance frame and report to stdout:
c
            call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps,
     +        multiplt,iadvfr,nppf,msgout,'longitude slices',iframe_plt)
            call hpdeallc(pplt,ier,1)
            call hpdeallc(pzpcol,ier,1)
          enddo		! irange=1,nzprange
c
c Plot at linear height scales:
c
  50      continue
          if (nhtscales.le.0) goto 150
          if (iyaxright.ge.1) then
            vp(1) = .13
            vp(2) = .87
          endif
          call alloc(pf2d,jmx*nzp)
c
c Get 2d field at full zp:
c
          call getfldlon(ifcpl(ip),f2d,imx,jmx,1,nzp,nzp,mxfproc,
     +      glon,logplt(ip),ip,fhist,kmx,nfhist,gcmzp,ixfhist,
     +      fccm,nccmlev,nccmfld,cplzp,ixfccm)
c
c Loop through desired height scales:
c
          do iscale = 1,nhtscales
c
c Set up current height scale:
c
            nhtscale = ifix((htscales(2,iscale)-htscales(1,iscale)) /
     +        htscales(3,iscale) + 1.0000001)
            call alloc(phtscale,nhtscale)
            do k=1,nhtscale
              htscale(k) = htscales(1,iscale) + (k-1)*htscales(3,iscale)
            enddo
            call alloc(pplt,jmx*nhtscale)
c
c Interpolate to current linear height scale and take log10 if
c necessary:
c
            call cuthtint(f2d,ht2d,jmx,nzp,plt,htscale,nhtscale,
     +        logplt(ip),spval,ier,1)
            if (ilog.gt.0) call log10f(plt,jmx*nhtscale,1.e-20,spval)
c
c Contour:
c
            if (iyaxright.gt.1) then
              if (multiplt.le.0) vplon(1) = .10
              vplon(2) = .80
            endif
            call setmultivp(vplon,iadvfr,nppf+1,multiplt,ipltrowcol,
     +        vp)
            if (mkplt.gt.0)
     +        call pltlon(plt,ht2d,zp,gcmlat,htscale,jmx,nhtscale,nzp,
     +          glon,10+iyaxright,p0,ut,vp,spval)
c
c Info labels, advance frame, report to stdout:
c
            call mklonlab(ifcpl(ip),glon,zmin,zmax,ciu,ilog,mtime,
     +        ip,blabs,tlabs)
            nppf = nppf+1               ! number of plots per frame
            call clearstr(msgout)
            if(mkplt.gt.0)
     +        call wrlab6(tlabs,toffset,tlabchsz,
     +                    blabs,boffset,blabchsz,vp,0)
            if (glon.ne.zmflag) then
              write(msgout,"('field ',a,' lon=',f7.2,' slt=',f5.1,
     +          ' ht=',f6.1,' to ',f6.1)") flab8(ip),glon,slt,
     +          htscale(1),htscale(nhtscale)
            else
              write(msgout,"('field ',a,
     +          ' (zonal means)  ht=',f6.1,' to ',f6.1)") 
     +          flab8(ip),htscale(1),htscale(nhtscale)
            endif
c
c Write to output xdr data file:
c
            if (iwrxdr.gt.0) then
              iframe_xdr = iframe_xdr+1
              call fminmax(plt,jmx*nhtscale,rmin,rmax,spval)
c             write(6,"('mklons: rmin,rmax=',2e12.4)") rmin,rmax
!             call wrxdr(flnm_xdr,plt,jmx,nhtscale,gcmlat,htscale,
!    +          'LATITUDE (DEG)','HEIGHT (KM)',
!    +          tlabs(2),tlabs(1),blabs(2),blabs(1))
              call wrxdr(trim(flnm_xdr),plt,jmx,nhtscale,gcmlat,htscale,
     +          'LATITUDE (DEG)','HEIGHT (KM)',
     +          tlabs(2),tlabs(1),blabs(2),blabs(1),mtime,0)
            endif       ! iwrxdr > 0
c
c Advance frame and report to stdout:
c
            call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps,
     +        multiplt,iadvfr,nppf,msgout,'longitude slices',iframe_plt)
            call hpdeallc(pplt,ier,1)
            call hpdeallc(phtscale,ier,1)
          enddo		! iscale=1,nhtscales
          call hpdeallc(pf2d,ier,1)
 150      continue
        enddo		! j=1,nlons
        call hpdeallc(pht2d,ier,1)
        call hpdeallc(pzp,ier,1)
c
c Increment field count:
c
        nplt = nplt+1
        if (nplt.gt.nfproc) then        ! done with all fields
          goto 102
        else
          goto 101                      ! go back for next field
        endif
 100    continue
      enddo     ! ip=1,mxfproc
 102  continue  ! nplt = 1,nfproc
c
c If doing multiplt and multiadvfr > 0, advance frame if page is only
c partially full (was advanced above if page is full)
c                  
      if (multiplt.gt.0.and.multiadvfr.eq.1.and.iadvfr.le.0)
     +  call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps,
     +    multiplt,1,nppf,' ','longitude slices',iframe_plt)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mklonlab(icpl,slon,zmin,zmax,ciu,ilog,mtime,ip,
     +  blabs,tlabs)
c
c Add info labels to lon slice contour plot:
c
      include "flxproc.h"
      include "ccm.h"
      include "tgcmhdr.h"
      character*(*) blabs(3),tlabs(3)
      character*8 datestr
      integer mtime(3)
c
c-------------------------- top labels --------------------------
c
c Info label with date, time and selected longitude:
c (bottom of 3 top labels)
c
      call clearstr(tlabs(1))
c
c date(2) in header is unreliable, so for now use model-day (nday from
c history header) as day-of-year (1-366) and ignore year:
c
c     iyd = (date(1)-(date(1)/100*100))*1000+date(2)
      iyd = 96000 + nday
      call iyd2date(iyd,imo,ida,iyr)
      write(datestr,"(i2.2,'/',i2.2)") imo,ida
      if (slon.ne.zmflag) then
        slt = fslt(dum,ut,slon,1)
c
c Temp: Write ut and slt with i2.2 for animxdr, so animation label 
c       does not jump around when ut and slt go from 1 to 2-digit values.
c
c       write(tlabs(1),"(a,' UT=',f4.1,' SLT=',f4.1,' LON=',f7.2)")
c    +    datestr(1:lenstr(datestr)),ut,slt,slon
        write(tlabs(1),"(a,' UT=',i2.2,' SLT=',i2.2,' LON=',f6.1)")
     +    datestr(1:lenstr(datestr)),ifix(ut),ifix(slt),slon
c
      else
c       write(tlabs(1),"('UT=',f4.1,'  ZONAL MEANS')") ut
        write(tlabs(1),"(a,' UT=',f4.1,'  ZONAL MEANS')") 
     +    datestr(1:lenstr(datestr)),ut
      endif
c     write(6,"('mklonlab: tlabs(1)=',a)") tlabs(1)(1:lenstr(tlabs(1)))
c
c Field label:
c (middle of 3 top labels)
c
      call clearstr(tlabs(2))
      if (icpl.gt.0) then
        if (ilog.le.0) then
c         write(tlabs(2),"('FIELD ',a,' (COUPLED ',a,'/',a,')')") 
c    +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model)),
c    +      ccmres(1:lenstr(ccmres))
          write(tlabs(2),"(a,' (',a,'/',a,')')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model)),
     +      ccmres(1:lenstr(ccmres))
        else
          write(tlabs(2),"('LOG10 ',a,' (COUPLED ',a,'/',a,')')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model)),
     +      ccmres(1:lenstr(ccmres))
        endif
      else
        if (ilog.le.0) then
          write(tlabs(2),"('FIELD ',a,' (',a,' ONLY)')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model))
        else
          write(tlabs(2),"('LOG10 ',a,' (',a,' ONLY)')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model))
        endif
      endif
c
c Top (third) of 3 top labels is optional 80-char annotation:
c
      call clearstr(tlabs(3))
      if (lenstr(lon_top_anno).gt.0) then
        write(tlabs(3),"(a)") lon_top_anno
      endif
c
c-------------------------- bottom labels --------------------------
c
c ccm history (bottom of 3 bottom labels):
c
      call clearstr(blabs(1))
      write(blabs(1),"('CCM FILE ',a,' (DAY ',f7.3,')')") 
     +  ccmlsd(1:lenstr(ccmlsd)),ccmday
c
c tgcm history (middle of 3 bottom labels):
c
      call clearstr(blabs(2))
      write(blabs(2),"('TGCM HISTORY ',a,
     +  ' (',i3.3,', ',i2.2,', ',i2.2,')')")
     +  histvols(ivol)(1:lenstr(histvols(ivol))),(mtime(i),i=1,3)
c
c Min,max label (top of 3 bottom labels):
c
      call clearstr(blabs(3))
      write(blabs(3),"('MIN,MAX=',2(1pe12.4),' INTERVAL=',1pe12.4)")
     +  zmin,zmax,ciu
      if (scalefac(ip).eq.1.) then
        write(blabs(3),"('MIN,MAX=',2(1pe12.4),' INTERVAL=',1pe12.4)")
     +    zmin,zmax,ciu
      else
        write(blabs(3),"('MIN,MAX=',2(1pe12.4),' INTERVAL=',1pe12.4,
     +    ' (X',1pe8.2,')')") zmin,zmax,ciu,scalefac(ip)
      endif
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getfldlon(icpl,fout,imx,jmx,izp0,izp1,nzpout,mxfproc,
     +  glon,log,ip,
     +  fgcm,ngcmzp,nfgcm,gcmzp,ixfgcm,
     +  fccm,nccmzp,nfccm,ccmzp,ixfccm)
c
c Call getfld to return tgcm lon slice in fout(jmx,nzpout) if icpl <= 0, 
c   or coupled lon slice if icpl > 0. Bottom and top of zp scale is from
c   izp0 to izp1
c
      real fout(jmx,nzpout),fccmsli(jmx,nccmzp),fgcmsli(jmx,ngcmzp)
      real fgcm(imx,ngcmzp,jmx,nfgcm),gcmzp(ngcmzp)
      real fccm(imx,nccmzp,jmx,nfccm),ccmzp(nccmzp)
      integer ixfgcm(mxfproc),ixfccm(mxfproc)
c
c Coupled slice:
c
      if (icpl.gt.0) then
        call getfld(fgcm,imx,ngcmzp,jmx,nfgcm,fgcmsli,jmx,ngcmzp,
     +    'LONSLICE',dum,glon,gcmzp,dum,log,ixfgcm,ip,ier)
        call getfld(fccm,imx,nccmzp,jmx,nfccm,fccmsli,jmx,nccmzp,
     +    'LONSLICE',dum,glon,ccmzp,dum,log,ixfccm,ip,ier)
        do k=izp0,izp1
          if (k.le.nccmzp) then
            fout(:,k-izp0+1) = fccmsli(:,k)
             
c           call fminmax(fccmsli(1,k),jmx,fmin,fmax,1.e36)
c           write(6,"('getfldlon: k=',i2,' izp0,1=',2i3,' nccmzp=',
c    +        i2,' fccmsli min,max=',2e12.4)") k,izp0,izp1,nccmzp,
c    +        fmin,fmax

          else
            fout(:,k-izp0+1) = fgcmsli(:,k-nccmzp+1)
          endif
        enddo
c
c Tgcm only slice:
c
      else
        call getfld(fgcm,imx,ngcmzp,jmx,nfgcm,fgcmsli,jmx,ngcmzp,
     +    'LONSLICE',dum,glon,gcmzp,dum,log,ixfgcm,ip,ier)
        do k=izp0,izp1
          fout(:,k-izp0+1) = fgcmsli(:,k)
        enddo
      endif
      return
      end