c
      subroutine mklats(it)
c
c Prepare latitude slices and pass to conlat where contours are drawn:
c ccmf(ngcmlon,ngcmlat,nccmlev,nccmfld) (pointer pccmf is in common)
c fhist(imx,kmx,jmx,nfhist) (pointer pfhist is in common)
c
      include "flxproc.h"
      include "ccm.h"
      pointer (pplt,plt(1)), (pzpcol,zpcol(1)), (pht2d,ht2d(1))
      real slats(jmx),vp(4),zpranges(2,jmx)
      data vp /.13,.87,.28,.90/
c
      write(6,"(' ')")
      write(6,"('Latitude slices at ut = ',f4.1,'  mtime = ',3i3)")
     +  ut,(mtimes(i,it),i=1,3)
      if (iyaxright.gt.1) then	! allow room on right for 2 right y-axes
        vp(1) = .12
        vp(2) = .77
      endif
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
      call alloc(pht2d,imx*ncplzp)
c
c Field loop:
c
      do ip=1,mxflds
        if (ixfhist(ip).le.0) goto 100
        ilog = 0
        if (ilat_log10.gt.0.and.logplt(ip).gt.0) ilog = 1
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
        nzpranges = 0
        zpranges(:,:) = spval
        if (ixfcpl(ip).gt.0) then	! coupled field
          do j=1,jmx
            if (flat_zprange(1,j).ne.spval.and.
     +          flat_zprange(2,j).ne.spval) then
              izp0 = ixfind(cplzp,ncplzp,flat_zprange(1,j),dzp)
              izp1 = ixfind(cplzp,ncplzp,flat_zprange(2,j),dzp)
              if (izp0.le.0) then
                write(6,"('Reset bottom of zp axis from ',
     +            f5.1,' to bottom of coupled system at zp ',f5.1)")
     +            flat_zprange(1,j),cplzp(1)
                flat_zprange(1,j) = cplzp(1)
              endif
              if (izp1.le.0) then
                write(6,"('Reset top of zp axis from ',
     +            f5.1,' to top of coupled system at zp ',f5.1)")
     +            flat_zprange(2,j),cplzp(ncplzp)
                flat_zprange(2,j) = cplzp(ncplzp)
              endif
              nzpranges = nzpranges+1
              zpranges(1,nzpranges) = flat_zprange(1,j)
              zpranges(2,nzpranges) = flat_zprange(2,j)
            endif
          enddo
        else		! tgcm only field
          do j=1,jmx
            if (flat_zprange(1,j).ne.spval.and.
     +          flat_zprange(2,j).ne.spval) then
              izp0 = ixfind(gcmzp,kmx,flat_zprange(1,j),dzp)
              izp1 = ixfind(gcmzp,kmx,flat_zprange(2,j),dzp)
              if (izp0.le.0) then
                write(6,"('Reset bottom of zp axis from ',
     +            f5.1,' to bottom of tgcm at zp ',f5.1)")
     +            flat_zprange(1,j),gcmzp(1)
                flat_zprange(1,j) = gcmzp(1)
              endif
              if (izp1.le.0) then
                write(6,"('Reset top of zp axis from ',
     +            f5.1,' to top of tgcm at zp ',f5.1)")
     +            flat_zprange(2,j),gcmzp(kmx)
                flat_zprange(2,j) = gcmzp(kmx)
              endif
              nzpranges = nzpranges+1
              zpranges(1,nzpranges) = flat_zprange(1,j)
              zpranges(2,nzpranges) = flat_zprange(2,j)
            endif
          enddo
        endif
c
c Selected latitudes loop:
c
        do j=1,nlats
          ixj = ixfind(gcmlat,jmx,slats(j),dlat)
          if (ixj.le.0) then
            write(6,"('>>> mklats: j=',i2,' bad slats(j)=',f10.2)")
     +        j,slats(j)
            stop 'ixj' 
          endif
          slats(j) = gcmlat(ixj)	! nearest latitude grid point
c
c Get ht2d(imx,ncplzp) for yaxzpht:
c
          call xfercpl_lat(ht2d,ccmf,fhist,imx,jmx,kmx,nccmlev,
     +      ncplzp,nfhist,nccmfld,ixj,1,ncplzp,ncplzp,ixfhist(ixz),
     +      ixz2)
c
c Loop through specified zp ranges:
c
          do irange = 1,nzpranges
c
c Is a coupled field:
c
            if (ixfcpl(ip).gt.0) then	! coupled field
              izp0 = ixfind(cplzp,ncplzp,zpranges(1,irange),dzp)
              izp1 = ixfind(cplzp,ncplzp,zpranges(2,irange),dzp)
              ny = izp1-izp0+1
              call alloc(pplt,imx*ny)
              call alloc(pzpcol,ny)
              do k=1,ny
                zpcol(k) = cplzp(izp0+k-1)
              enddo
              call xfercpl_lat(plt,ccmf,fhist,imx,jmx,kmx,nccmlev,
     +          ncplzp,nfhist,nccmfld,ixj,izp0,izp1,ny,ixfhist(ip),
     +          ixfcpl(ip))
              if (ilog.gt.0) call log10f(plt,imx*ny,1.e-20,spval)
c
c Contour:
c     subroutine conlat(f,xx,yy,nx,ny,vp,cmin,cmax,cint,
c    +  scalefac,iht,icolor,ibox_clabs,gkswid,zmin,zmax,ciu)
c
              call conlat(plt,gcmlon,zpcol,imx,ny,vp,
     +          0.,0.,0.,1.,0,0,0,gkswid,zmin,zmax,ciu)
c
c Add local time axis:
c
              call sltxax(-180.,180.,ut,.10)
c
c Add extra right-hand axes:
c     subroutine yaxzpht(hts1d,hts2d,gcmzp,id1,kmx,dzp,yy,ny,p0,
c    +  iht,iyaxright)
c     real hts2d(id1,kmx), where ht2d=hts(imx,ncplzp), 
c          gcmzp(kmx),     where gcmzp=cplzp(ncplzp),
c          yy(ny),         where yy=zpcol(ny)
c
              if (iyaxright.gt.0)
     +          call yaxzpht(dum,ht2d,cplzp,imx,ncplzp,dzp,zpcol,ny,
     +            p0,0,iyaxright)
c
c Add info labels and advance frame:
c
              call lablat(1,slats(j),zmin,zmax,ciu,ilog,it,ip)
              call frame
              iframe = iframe+1
              write(6,"('Frame ',i3,' coupled field ',a,' lat=',f7.2,
     +          ' zp=',f6.1,' to ',f6.1)") iframe,flab8(ip),slats(j),
     +          zpcol(1),zpcol(ny)
c
c Is a non-coupled field (tgcm only):
c
            else				! tgcm field only
              izp0 = ixfind(gcmzp,kmx,zpranges(1,irange),dzp)
              izp1 = ixfind(gcmzp,kmx,zpranges(2,irange),dzp)
              ny = izp1-izp0+1
              call alloc(pplt,imx*ny)
              call alloc(pzpcol,ny)
              do k=1,ny
                zpcol(k) = gcmzp(izp0+k-1)
              enddo
              call xfergcm_lat(pltgcm,fhist,imx,jmx,kmx,nfhist,ixj,
     +          izp0,izp1,ny,ixfhist(ip))
              if (ilog.gt.0) call log10f(pltgcm,imx*ny,1.e-20,spval)
              call conlat(pltgcm,gcmlon,zpcol,imx,ny,vp,
     +          0.,0.,0.,1.,0,0,0,gkswid,zmin,zmax,ciu)
              call sltxax(-180.,180.,ut,.10)
              call lablat(0,slats(j),zmin,zmax,ciu,ilog,it,ip)
              call frame
              iframe = iframe+1
              write(6,"('Frame ',i3,' tgcm field ',a,' lat=',f7.2,
     +          ' min,max=',2e12.4,' zp=',f6.1,' to ',f6.1)") 
     +          iframe,flab8(ip),slats(j),zpcol(1),zpcol(ny)
            endif
            call hpdeallc(pplt,ier,1)
            call hpdeallc(pzpcol,ier,1)
 200        continue
          enddo		! irange=1,nzprange
        enddo		! j=1,nlats
 100    continue
      enddo		! ip=1,mxflds
      call hpdeallc(pht2d,ier,1)

      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine xfercpl_lat(pltcpl,fccm,fgcm,imx,jmx,kmx,nccmlev,
     +  ncplzp,nfgcm,nfccm,ixj,izp0,izp1,nzp,ipgcm,ipccm)
c
c Define pltcpl coupled latitude slice from fccm and fgcm at levels 
c corresponding to izp0 to izp1 (where izp0,izp1 are bottom and top
c of zprange within cplzp):
c
      real pltcpl(imx,nzp),fccm(imx,jmx,nccmlev,nfccm),
     +  fgcm(imx,kmx,jmx,nfgcm)
c
      kk = 0
      do k=1,nccmlev
        if (k.ge.izp0.and.k.le.izp1) then
          kk = kk+1
          pltcpl(:,kk) = fccm(:,ixj,k,ipccm)
        endif
      enddo
      do k=nccmlev+1,ncplzp
        if (k.ge.izp0.and.k.le.izp1) then
          kk = kk+1
          pltcpl(:,kk) = fgcm(:,k-nccmlev+1,ixj,ipgcm)
        endif
      enddo
      return
      end 
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine xfergcm_lat(pltgcm,fgcm,imx,jmx,kmx,nfgcm,ixj,
     +  ny,izp0,izp1,ipgcm)
      real pltgcm(imx,ny),fgcm(imx,kmx,jmx,nfgcm)
c
      kk = 0
      do k=1,kmx
        if (k.ge.izp0.and.k.le.izp1) then
          kk = kk+1
          pltgcm(:,kk) = fgcm(:,k,ixj,ipgcm)
        endif
      enddo
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine lablat(icpl,slat,zmin,zmax,ciu,ilog,it,ip)
c
c Add info labels to lat slice contour plot:
c
      include "flxproc.h"
      include "ccm.h"
      character*120 blab0,blab1,blab2,tlab0,tlab1
      data toffset/.025/, boffset/.20/
c
c ccm history (bottom of 3 bottom labels):
c
      call clearstr(blab0)
      write(blab0,"('CCM FILE ',a,' (DAY ',f7.3,')')") 
     +  ccmlsd(1:lenstr(ccmlsd)),ccmday
c
c tgcm history (middle of 3 bottom labels):
c
      call clearstr(blab1)
      write(blab1,"('TGCM HISTORY ',a,' (',i3,', ',i2,', ',i2,')')")
     +  histvols(ivol)(1:lenstr(histvols(ivol))),
     +  (mtimes(i,it),i=1,3)
c
c Min,max label (top of 3 bottom labels):
c
      call clearstr(blab2)
      write(blab2,"('MIN,MAX=',2(1pe12.4),' INTERVAL=',1pe12.4)")
     +  zmin,zmax,ciu
c     if (scalefac(ip).eq.1.) then
c       write(blab2,"('MIN,MAX=',2(1pe12.4),' INTERVAL=',1pe12.4)")
c    +    zmin,zmax,ciu
c     else
c       write(blab2,"('MIN,MAX=',2(1pe12.4),' INTERVAL=',1pe12.4,
c    +    ' (X',1pe8.2,')')") zmin,zmax,ciu,scalefac(ip)
c     endif
c
c Info label with time and selected latitude:
c (bottom of 2 top labels)
c
      call clearstr(tlab0)
      write(tlab0,"('UT=',f4.1,' LAT=',f7.2,' (DEG)')") ut,slat
c
c Field label:
c (top of 2 top labels)
c
      call clearstr(tlab1)
      if (icpl.gt.0) then
        if (ilog.le.0) then
          write(tlab1,"('FIELD ',a,' (COUPLED ',a,'/',a,')')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model)),
     +      ccmres(1:lenstr(ccmres))
        else
          write(tlab1,"('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(tlab1,"('FIELD ',a,' (',a,' ONLY)')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model))
        else
          write(tlab1,"('LOG10 ',a,' (',a,' ONLY)')") 
     +      flab8(ip)(1:lenstr(flab8(ip))),model(1:lenstr(model))
        endif
      endif
c
c Add the 5 labels to the plot:
c     subroutine wrlab5(tlab0,tlab1,toffset, blab0,blab1,blab2,
c    +  boffset, tlab0_sz, tlab1_sz, blab0_sz, blab1_sz, blab2_sz, 
c    +  ihq)
c
      call wrlab5(tlab0,tlab1,toffset, blab0,blab1,blab2,boffset,
     +  .015,.015,.014,.014,.014,0)
      return
      end