c
      subroutine mklats(mtime)
c
c Prepare latitude slices and pass to conlat 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))
      real slats(jmx),vp(4),vplat(4),zpranges(2,mxvrange),tlabchsz(3),
     +  blabchsz(3),htscales(3,mxvrange)
      integer mtime(3)
      character*96 tlabs(3),blabs(3)
      character*80 msgout
      data vplat /.13,.87,.27,.90/
      data tlabchsz/.020,.020,.020/, blabchsz/.02,.02,.02/
      data toffset /.03/, boffset/.26/
c
      write(6,"(/'Latitude slices at ut = ',f4.1,'  tgcm mtime = ',3i3,
     +  ' ccmday = ',f10.4)") ut,(mtime(i),i=1,3),ccmday
      if (multiadvfr.gt.0) nppf = 0
c
c Determine number of selected lats, and place in local array slats:
c
      nlats = 0
      do j=1,jmx
        if (flats(j).ne.spval) then
          nlats = nlats+1
          slats(nlats) = flats(j)
        endif  
      enddo
      if (nlats.le.0) return
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 (ilat_log10.gt.0.and.logplt(ip).gt.0) ilog = 1
c
c Define plot options (common in pltopt.h, included by pltlat)
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,imx*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 Setyverts also checks linear height scales.
c
        call setyverts(flat_zprange,zpranges,nzpranges,
     +                 flat_htscale,htscales,nhtscales,
     +                 mxvrange,zp,nzp,dzp,spval,0)
        if (nzpranges.le.0.and.nhtscales.le.0) goto 100
c
c Selected latitudes loop:
c
        do j=1,nlats
c
c Get heights ht2d(imx,nzp) at current lat:
c
          call getfldlat(ifcpl(ip),ht2d,imx,jmx,1,nzp,nzp,mxfproc,
     +      slats(j),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
            vplat(1) = .12
            vplat(2) = .77
          endif
          do irange = 1,nzpranges
            izp0 = ixfind(zp,nzp,zpranges(1,irange),dzp)
            izp1 = ixfind(zp,nzp,zpranges(2,irange),dzp)
            ny = izp1-izp0+1
            call alloc(pplt,imx*ny)
            call alloc(pzpcol,ny)
            do k=1,ny
              zpcol(k) = zp(izp0+k-1)
            enddo
            call getfldlat(ifcpl(ip),plt,imx,jmx,izp0,izp1,ny,mxfproc,
     +        slats(j),logplt(ip),ip,fhist,kmx,nfhist,gcmzp,ixfhist,
     +        fccm,nccmlev,nccmfld,cplzp,ixfccm)
            if (ilog.gt.0) call log10f(plt,imx*ny,1.e-20,spval)
c
c Contour (pltlat contours, adds slt x-axis and extra right-hand yaxes, 
c   and saves zmin,zmax,ciu in pltopt.h):
c     subroutine pltlat(f,fht,zp,xx,yy,nx,ny,nzp,glat,iyax,p0,ut,vp,
c    +  spval)
c
            call setmultivp(vplat,iadvfr,nppf+1,multiplt,ipltrowcol,vp)
            if (mkplt.gt.0) then
              call pltlat(plt,ht2d,zp,gcmlon,zpcol,imx,ny,nzp,
     +          slats(j),iyaxright,p0,ut,vp,spval)
            endif
c
c Add info labels and advance frame:
c
            call mklatlab(ifcpl(ip),slats(j),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
            call clearstr(msgout)
            write(msgout,"('field ',a,' lat=',f7.2,
     +        ' zp=',f6.1,' to ',f6.1)") flab8(ip),slats(j),
     +        zpcol(1),zpcol(ny)
c
c Save to xdr file:
c
            if (iwrxdr.gt.0) then
              iframe_xdr = iframe_xdr+1
!             call wrxdr(flnm_xdr(1:lenstr(flnm_xdr)),plt,imx,ny,
!    +          gcmlon,zpcol,'LONGITUDE (DEG)','ZP (LN(P0/P))',
!    +          tlabs(2),tlabs(1),blabs(2),blabs(1))
              call wrxdr(trim(flnm_xdr),plt,imx,ny,
     +          gcmlon,zpcol,'LONGITUDE (DEG)','ZP (LN(P0/P))',
     +          tlabs(2),tlabs(1),blabs(2),blabs(1),mtime,0)
            endif
c
c Advance frame and release heap space:
c
            call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps,
     +        multiplt,iadvfr,nppf,msgout,'latitude slices',iframe_plt)
            call hpdeallc(pplt,ier,1)
            call hpdeallc(pzpcol,ier,1)
 200        continue
          enddo		! irange=1,nzpranges
c
c Plot at linear height scales:
c
 50       continue
          if (nhtscales.le.0) goto 150
          if (iyaxright.ge.1) then
            vplat(1) = .12
            vplat(2) = .80
          endif
          call alloc(pf2d,imx*nzp)
c
c Get 2d field at full zp:
c
          call getfldlat(ifcpl(ip),f2d,imx,jmx,1,nzp,nzp,mxfproc,
     +      slats(j),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,imx*nhtscale)
c
c Interpolate to current linear height scale and take log10 if
c necessary:
c
            call cuthtint(f2d,ht2d,imx,nzp,plt,htscale,nhtscale,
     +        logplt(ip),spval,ier,1)
            if (ilog.gt.0) call log10f(plt,imx*nhtscale,1.e-20,spval)
c
c Contour, labels, etc:
c
            call setmultivp(vplat,iadvfr,nppf+1,multiplt,ipltrowcol,vp)
            if (mkplt.gt.0)
c     subroutine pltlat(f,fht,zp,xx,yy,nx,ny,nzp,glat,iyax,p0,ut,vp,
c    +  spval)
     +        call pltlat(plt,ht2d,zp,gcmlon,htscale,imx,nhtscale,nzp,
     +          slats(j),10+iyaxright,p0,ut,vp,spval)
            call mklatlab(ifcpl(ip),slats(j),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
            call clearstr(msgout)
            write(msgout,"('field ',a,' lat=',f7.2,
     +        ' ht=',f6.1,' to ',f6.1)") flab8(ip),slats(j),
     +        htscale(1),htscale(nhtscale)
c
c Save to xdr file:
c     subroutine wrxdr(flnm,f,nx,ny,xx,yy,xlab,ylab,lab1,
c    +  lab2,lab3,lab4)
c
            if (iwrxdr.gt.0) then
              iframe_xdr = iframe_xdr+1
!             call wrxdr(flnm_xdr(1:lenstr(flnm_xdr)),plt,imx,nhtscale,
!    +          gcmlon,htscale,'LONGITUDE (DEG)','HEIGHT (KM)',
!    +          tlabs(2),tlabs(1),blabs(2),blabs(1))
              call wrxdr(trim(flnm_xdr),plt,imx,nhtscale,
     +          gcmlon,htscale,'LONGITUDE (DEG)','HEIGHT (KM)',
     +          tlabs(2),tlabs(1),blabs(2),blabs(1),mtime,0)
            endif
c
c Advance frame and relase heap space:
c
            call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps,
     +        multiplt,iadvfr,nppf,msgout,'latitude slices',iframe_plt)
            call hpdeallc(phtscale,ier,1)
            call hpdeallc(pplt,ier,1)
          enddo		! iscale=1,nhtscales
          call hpdeallc(pf2d,ier,1)
 150      continue
        enddo		! j=1,nlats
        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
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,' ','latitude slices',iframe_plt)
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mklatlab(icpl,slat,zmin,zmax,ciu,ilog,mtime,ip,
     +  blabs,tlabs)
c
c Add info labels to lat slice contour:
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 time and selected latitude:
c (bottom of 3 top labels)
c
      call clearstr(tlabs(1))
c
c Temp: put integer ut on label so animxdr animation labels don't jump
c       when ut goes from single to double digit or double to single digit.
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     write(tlabs(1),"('UT=',f4.1,' LAT=',f7.2,' (DEG)')") ut,slat
c     iyd = (date(1)-(date(1)/100*100))*1000+date(2)
c     call iyd2date(iyd,imo,ida,iyr)
c     write(datestr,"(i2.2,'/',i2.2)") imo,ida
      iyd = 96000 + nday
      call iyd2date(iyd,imo,ida,iyr)
      write(datestr,"(i2.2,'/',i2.2)") imo,ida
      write(tlabs(1),"(a,' UT=',i2.2,' LAT=',f7.2,' (DEG)')") 
     +  datestr(1:lenstr(datestr)),ifix(ut),slat
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,' (COUPLED ',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
c         write(tlabs(2),"('FIELD ',a,' (',a,' ONLY)')") 
c    +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model))
          write(tlabs(2),"(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(lat_top_anno).gt.0) then
        write(tlabs(3),"(a)") lat_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 getfldlat(icpl,fout,imx,jmx,izp0,izp1,nzpout,mxfproc,
     +  glat,log,ip,
     +  fgcm,ngcmzp,nfgcm,gcmzp,ixfgcm,
     +  fccm,nccmzp,nfccm,ccmzp,ixfccm)
c
c Call getfld to return tgcm lat slice in fout(imx,nzpout) if icpl <= 0, 
c   or coupled lat slice if icpl > 0. Bottom and top of zp scale is from
c   izp0 to izp1
c
      real fout(imx,nzpout),fccmsli(imx,nccmzp),fgcmsli(imx,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,imx,ngcmzp,
     +    'LATSLICE',glat,dum,gcmzp,dum,log,ixfgcm,ip,ier)
        call getfld(fccm,imx,nccmzp,jmx,nfccm,fccmsli,imx,nccmzp,
     +    'LATSLICE',glat,dum,ccmzp,dum,log,ixfccm,ip,ier)
        do k=izp0,izp1
          if (k.le.nccmzp) then
            fout(:,k-izp0+1) = fccmsli(:,k)
          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,imx,ngcmzp,
     +    'LATSLICE',glat,dum,gcmzp,dum,log,ixfgcm,ip,ier)
        do k=izp0,izp1
          fout(:,k-izp0+1) = fgcmsli(:,k)
        enddo
      endif
      return
      end