c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timesdif/pltlatht.f
c------------------------------------------------------------------
c
      subroutine pltlatht(utc,fld,htscale,nhtscale)
c
c Contour longitude slices (lat on x-axis, zp on y-axis)
c
      include 'gettime.h'
      include 'timesdif.h'
      character*56 toplab,lab
      dimension plt(jmx,nhtscale),viewport(4),yaxht(kmx),rimx(imx),
     +  fld(imx,kmx,jmx,nplfld),htscale(nhtscale)
      data viewport /.15,.89,.26,.91/
c
      if (nlon.le.0) then
        write(6,"('pltlatht: no selected longitudes -- returning')")
        return
      endif
      write(6,"('pltlatht: nhtscale=',i3,' Using height scale = ',
     +  /(6f10.2))") nhtscale,htscale
      xmid = 0.5*(viewport(1)+viewport(2))
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('XC1',gcmlat(1))
      call cpsetr('XCM',gcmlat(jmx))
      call cpsetr('YC1',htscale(1))
      call cpsetr('YCN',htscale(nhtscale))
      call cpsetr('ILX',xmid)
      call cpsetr('ILY',-.20)
      call cpseti('ILP',0)
      call set(viewport(1),viewport(2),viewport(3),viewport(4),
     +         -90.,90.,htscale(1),htscale(nhtscale),1)
c
c Field loop: 
      ixp = 0
      do 200 ip=1,ntimefld
        if (ipltime(ip).le.0) goto 200
        ixp = ixp+1
c
c Selected longitude loop:
c
        do 100 i=1,nlon
          if (slon(i).ne.zmlon) then
            ixlon = ixfind(gcmlon,imx,slon(i),dlon)
            if (ixlon.le.0) then
              write(6,"('>>> pltlatht: bad longitude=',f10.3,' i=',i3,
     +          ' skipping this lon')") slon(i),i
              goto 100
            endif
            rlon = gcmlon(ixlon)
          else
            ixlon = ifix(zmlon)
            rlon = slon(i)
          endif
c
          do k=1,nhtscale
            do j=1,jmx
              if (slon(i).ne.zmlon) then 
                plt(j,k) = fld(ixlon,k,j,ixp)
              else                              ! zonal means
                rimx(:) = fld(:,k,j,ixp)
                plt(j,k) = calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              endif
            enddo
          enddo
          if (icolor.le.0) then
            call contour(plt,jmx,jmx,nhtscale,cint(ip),cmin(ip),
     +        cmax(ip))
          else
            call conclr(plt,jmx,jmx,nhtscale,cint(ip),cmin(ip),cmax(ip))
          endif
          call labrect(gcmlat,jmx,htscale(nhtscale),nhtscale,'LATITUDE',
     +      'HT',0.)
c
c Top label:
          if (slon(i).ne.zmlon) then
            if (wtime(ip).gt.0.) then
              write(toplab,"(a8,' (% DIFFERENCE)  LON=',f8.2,
     +          ' UT=',f6.2,9x)") timelab_short(ip),rlon,utc
            else
              write(toplab,"(a8,' (RAW DIFFERENCE)  LON=',f8.2,
     +          ' UT=',f6.2,7x)") timelab_short(ip),rlon,utc
            endif
          else
            if (wtime(ip).gt.0.) then
              write(toplab,"(a8,' (% DIFFERENCE)  ZONAL MEANS  UT=',
     +          f6.2,9x)") timelab_short(ip),utc
            else
              write(toplab,"(a8,' (RAW DIFFERENCE)  ZONAL MEANS  UT=',
     +          f6.2,7x)") timelab_short(ip),utc
            endif
          endif
          call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.05,0.)
c
c History vols label:
          write(lab,"(2x,'Difference of TIGCM histories ',a8,
     +      ' minus ',a8,1x)") cpertvol,ccntrvol
          call wrlab(lab(1:lenstr(lab)),xmid,viewport(3)-.09,.012)
          call frame
          iframe = iframe+1
          if (slon(i).ne.zmlon) then
            write(6,"('pltlatht: frame ',i4,' field ',a8,' lon=',f8.2,
     +        )") iframe,timelab_short(ip),rlon
          else
            write(6,"('pltlatht: frame ',i4,' field ',a8,
     +        ' (zonal means)')")
     +        iframe,timelab_short(ip)
          endif
c
c End selected longitude loop
 100    continue
c
c End field loop: ip=1,ntimefld
 200  continue
c
      return
      end