c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcm/pltutzpm.f
c------------------------------------------------------------------
c
      subroutine pltutzpm(ic,it)
c
c Contour global means with ut on x-axis, zp on y-axis
c

      include 'tgcmparam.h'
      include 'tgcmhdr.h'
      include 'input.h'
      include 'tigcmfld.h'
      include 'tgcmlab.h'
      include 'tigcm.h'
      dimension plotp(mxtms,kmx),fut(kmx,ntigcmf)
      character*56 toplab,lab,lab1
c
c Declarations for axes labeling:
c
      parameter (ny=13)
      dimension numx(mxtms),numy(ny)
      data numy/-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5/
      character*8 fmtx,fmty
      character*16 labx,laby
      data fmty/'(i2)    '/
      data nfmty/2/, mnry/2/
      data laby/'    -LN P/P0    '/
      data icp/1/
c
c For finding global means:
      data re /6371.e+5/
      data ncalls/0/
      save ncalls,sud,sum,dtor
c
      ncalls = ncalls+1
      if (ncalls.eq.1) then
        sud = re*re*dlon*dlat*dtr*dtr
        sum = 0.
        do 60 i=1,imx
        do 60 j=1,jmx
 60       sum = sum + sud*cos(gcmlat(j)*dtr)
      endif
c
c Write global means to scratch file for current ut:
c
      fut(:,:) = 0.     
      lu = luglbm1
      do 70 ip=1,ntigcmf
        if (iptigcm(ip).le.0) goto 70
c
c Some fields not yet available or inappropriate for this
c type of display (e.g., fof2 and hmf2 are height-independent
c so are not plotted zp vs ut):
c NOTE: changes to this conditional should also be made to
c       conditional in do 100 read loop below
c
        if (ip.eq.ixrat.or.ip.eq.ixfof2.or.ip.eq.ixhmf2) goto 70
        if (ip.eq.ixrho) then
          do 80 k=1,kmx
          do 80 j=1,jmx
          do 80 i=1,imx
            abb = cos(gcmlat(j)*dtr)*sud/sum
            fut(k,ip) = fut(k,ip) + td(i,k,j) * abb 
 80       continue
        else
          do 82 k=1,kmx
          do 82 j=1,jmx
          do 82 i=1,imx
            abb = cos(gcmlat(j)*dtr)*sud/sum
            fut(k,ip) = fut(k,ip) + pnt(i,k,j,ip) * abb 
 82       continue
        endif 
 70   continue
      if (ic.gt.icp) rewind lu
      write(lu) fut
c
c If just wrote last ut, jump down to make plots:
      if (it.eq.ntms) goto 500
c
c Otherwise exit, since not doing plots yet:
      goto 501
c
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
c If doing last ut, read array back from scratch file, and contour:
c
 500  continue
      rewind lu
      do 200 iit=1,ntms
        read(lu) fut  
        do 205 k=1,kmx
        do ip=1,ntigcmf
          gcmuts(iit,k,ip) = fut(k,ip) 
        enddo
 205    continue
 200  continue
c
c Set up conpack:
c
      call cpseti('SET',0)
      call cpseti('MAP - mapping flag',0)
      utmapl = 0.09
      utmapr = 0.83
      utmapb = 0.15
      utmapt = 0.85
      rilx = 0.5*(utmapl+utmapr)
      rily = utmapb - 0.35
      call cpsetr('ILX - info label x-coord',rilx)
      call cpsetr('ILY - info label y-coord',rily)
      utinc1 = float(mhut(1))
      utincn = utinc1 + float(ntms-1)
      xc = utinc1
      xd = utincn
      call set(utmapl,utmapr,utmapb,utmapt,
     +         utinc1,utincn,gcmzp(1),gcmzp(kmx),1)
      yc = gcmzp(1)
      yd = gcmzp(kmx)
      write(toplab,"('Global Mean  ',a32,11x)") runlab
c
c Set up x-axis:
c
      nx = 0
      if (ndays.eq.1) then
        mnrx = 3
        write(labx,"(' UT (day ',i5,') ')") iyd(1)
        do 45 i=1,ntms
          if (mod(mhut(i),6).eq.0) then
            nx = nx+1
            numx(nx) = mhut(i)
          endif
 45     continue
      elseif (ndays.le.3) then
        mnrx = 3
        write(labx,"('days ',i5,'-',i5)") iyd(1),iyd(ndays)
        do 50 i=1,ntms
          if (mod(mhut(i),6).eq.0) then
            nx = nx+1
            numx(nx) = mhut(i)
          endif
 50     continue
      else
        mnrx = 4
        do 55 i=1,ndays
          write(labx,"('  DAY OF ',i4,3x)") 1900+iyd(i)/1000
          nx = nx+1
          numx(nx) = iyd(i)-iyd(i)/1000*1000
 55     continue
        nx = nx+1
        iydlast = iyd(ndays)+1
        numx(nx) = iydlast-iydlast/1000*1000
      endif
      fmtx = '(i2)    '
      nfmtx = 2
      if (numx(1).gt.99) then
        fmtx = '(i3)    '
        nfmtx = 3
      endif
c
c Field loop:
      write(6,"(' ')")
      do 100 ip=1,ntigcmf
        if (iptigcm(ip).le.0) goto 100
        if (ip.eq.ixrat) then
          write(6,"(' ')")
          write(6,"('>>> pltutzpm: plots for ',a,' not yet',
     +      ' available -- skipping this field')") 
     +      labtigcm_short(ip)
          goto 100
        elseif (ip.eq.ixfof2.or.ip.eq.ixhmf2) then
          write(6,"(' ')")
          write(6,"('>>> pltutzpm: ',a,' is height independent',
     +      ' and therefore inappropriate for zp vs ut plot --',
     +      /' skipping this field')") labtigcm_short(ip)
          goto 100
        else
c
c Always take logs of species densities since zp on y-axis:
          if (igcmlog(ip).gt.0) then
            do 105 k=1,kmx
            do 105 iit=1,ntms
              if (gcmuts(iit,k,ip).le.1.e-20) plotp(iit,k) = cpspval
              if (gcmuts(iit,k,ip).gt.1.e-20) 
     +          plotp(iit,k) = alog10(gcmuts(iit,k,ip))
 105        continue
            write(lab1,"('LOG10 ',a50)") labtigcm(ip)(1:50)
          else 
            do 110 k=1,kmx
              plotp(:,k) = gcmuts(:,k,ip)
 110        continue
            write(lab1,"(a)") labtigcm(ip)
          endif
          call bwcon(plotp,mxtms,ntms,kmx,xc,xd,yc,yd,
     +                 conints(ip),conmins(ip),conmaxs(ip),
     +                 1,nx,numx,mnrx,labx,fmtx,nfmtx,
     +                 ny,numy,mnry,laby,fmty,nfmty)
          call top2lab(lab1,0.10,toplab,0.05,.017)
c
c Bottom history label:
          write(lab,"(7x,'TGCM History = ',4a8,2x)") 
     +      (output(ii,1),ii=1,3),output(3,2)
          call botlab(lab,0.5*(utmapl+utmapr),.055)
c
c Wrap it up:
          call frame
          iframe=iframe+1
          write(6,"('Pltutzpm frame ',i3,': ',a)") iframe,lab1(1:43)
        endif
c
c End fields loop:
 100  continue
c
 501  continue
      icp = ic
      return
      end