c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timesdif/pltpolht.f
c------------------------------------------------------------------
c
      subroutine pltpolht(utc,pert,cntr)
c
      include 'gettime.h'
      include 'timesdif.h'
      character*56 toplab,lab
      character*1 hem
      dimension viewport(4),pert(imx,kmx,jmx,nplfld),
     +  cntr(imx,kmx,jmx,ntimefld),pertglb(imx,jmx),cntrglb(imx,jmx),
     +  plt(imx,jmx)
      data viewport /.125,.875,.150,.900/
c
      xmid = 0.5*(viewport(1)+viewport(2))
      call cpseti('SET',0)
      call cpseti('MAP',1)
      call cpsetr('XC1',gcmlon(1))
      call cpsetr('XCM',gcmlon(imx))
      call cpsetr('ILX',0.5)
      call cpsetr('ILY',-.15)
      ixz = 0
      do ip=1,ntimefld
        if (ipltime(ip).gt.0) then
          ixz = ixz+1
          if (ip.eq.itxz) goto 99
        endif
      enddo
 99   continue
c
c Perimeter latitude loop:
c
      do ih=1,npollat
        nlats = (90.-abs(perimlat(ih))+dlat/2.)/dlat
        if (perimlat(ih).ge.0.) then 
          hem = 'N'
          yc1 = perimlat(ih)
          ycn = gcmlat(jmx)
          iskip = jmx-nlats
          rot = -utc*15.
          plat = 90.
        else              
          hem = 'S'
          yc1 = gcmlat(1)
          ycn = perimlat(ih)
          iskip = 0
          rot = utc*15.-180.
          plat = -90.
        endif
        call mkproj('ST',viewport,plat,0.,rot,perimlat(ih),dum)
        call cpsetr('YCN',ycn)
        call cpsetr('YC1',yc1)
c
c Field loop: 
c
        ixp = 0
        do 200 ip=1,ntimefld
          if (ipltime(ip).le.0) goto 200
          ixp = ixp+1
          if (ip.eq.itxz) goto 200
c
c Selected height loop:
c
          do 100 ihts = 1,nhts
c
c Interpolate perturbed to selected height:
            call glbhtint(pert(1,1,1,ixp),pert(1,1,1,ixz),imx,kmx,jmx,
     +        pertglb,shts(ihts),1,itimelog(ip),cpspval,ier,iprnt)
c
c Interpolate control to selected height:
            call glbhtint(cntr(1,1,1,ip),cntr(1,1,1,itxz),imx,kmx,jmx,
     +        cntrglb,shts(ihts),1,itimelog(ip),cpspval,ier,iprnt)
c      
c Take diffs of interpolated fields (results in pertglb),
c   and contour:
            call mkdifglb(pertglb,cntrglb,imx,jmx,wtime(ip),cpspval)
            do j=iskip+1,jmx
              plt(:,j-iskip) = pertglb(:,j)
            enddo
            if (icolor.le.0) then
              call contour(plt,imx,imx,nlats,cint(ip),cmin(ip),cmax(ip))
            else
              call conclr(plt,imx,imx,nlats,cint(ip),cmin(ip),cmax(ip))
            endif
            if (icont.gt.0) call maplot
            call labpolar(viewport(1),viewport(2),viewport(3),
     +        viewport(4),hem)
c
c Top label:
            if (wtime(ip).gt.0.) then
              write(toplab,"(a8,' (% DIFFERENCE) HT=',f5.1,' UT=',
     +          f5.2,' PERIMLAT=',f5.1)")
     +          timelab_short(ip),shts(ihts),utc,perimlat(ih)
            else
              write(toplab,"(a8,' (RAW DIFFERENCE) HT=',f5.1,' UT=',
     +          f5.2,' PERIM=',f5.1,1x)")
     +          timelab_short(ip),shts(ihts),utc,perimlat(ih)
            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)-.08,.012)
            call frame
            iframe = iframe+1
            write(6,"('pltpolht diffs: frame ',i4,' field ',a8,
     +        ' ht=',f5.1,' perimlat=',f6.2)") 
     +        iframe,timelab_short(ip),shts(ihts),perimlat(ih)
c
c End selected height loop: ihts=1,nhts
 100      continue
c
c End field loop: ip=1,ntimefld
 200    continue
c
c End hemisphere loop
      enddo
c
      return
      end