c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcmdif/pltglb.f
c------------------------------------------------------------------
c
      subroutine pltglb(proj,utc,pert,cntr)
c
      include 'tigcmdif.h'
      character*56 toplab
      character*80 rec80,histlab,fieldlab
      character*24 flnm0,flnm1
      character*2 proj
      dimension plt(imx,jmx),viewport(4),pert(imx,kmx,jmx,nfget),
     +  udat(imx,jmx),vdat(imx,jmx),lonlab(9),latlab(7),
     +  cntr(imx,kmx,jmx,nfget),glb(imx,jmx)
      data viewport /.14,.92,.32,.82/, dum/0./, dgrid/30./
      data latlab/-90,-60,-30,0,30,60,90/
      data iblack/0/, iwhite/1/
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('YC1',gcmlat(1))
      call cpsetr('YCN',gcmlat(jmx))
      call cpsetr('ILX',xmid)
      call cpsetr('ILY',-.48)
      call cpsetr('ILS',.015)
      call cpseti('ILP',0)
      call mapstr('GR',dgrid)
      vlaboff = -.85
c
c Field loop: 
c (height interp, and log10 is done separately on pert and cntr
c  in getglb if necessary, before mkdifs)
c
      write(6,"(' ')")
      do 200 ip=1,nftot
        if (ifplt(ip).le.0) goto 200
c
c Selected pressure/height loop:
        do 100 izp = 1,npls+nhts
          ixpl = 0
          ht = 0.
          if (izp.le.npls) then
            ixpl = ixfind(gcmzp,kmx,spls(izp),dzp)
            if (ixpl.le.0) then
              write(6,"('>>> pltglb: bad selected pressure=',f10.3,
     +          ' izp=',i3)") spls(izp),izp
              goto 100
            endif
          else
            ht = shts(izp-npls)
          endif
          if (izp.gt.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) goto 200
          if (ip.eq.ixz.and.izp.gt.npls) goto 200
          ipercent = 0
          if (idif.eq.1.and.iperc(ip).eq.1) ipercent = 1
          if (ip.ne.ixunvn.and.ip.ne.ixuivi) then
            call getglb(pert,plt,ixpl,ht,ip,1)
            call getglb(cntr,vdat,ixpl,ht,ip,0)
            call mkdifs(plt,vdat,imx*jmx,ipercent,cpspval)
          endif
c
c Projection loop (projection orientations):
c (ixslt() changes plon if necessary from r12flag to local noon longitude)
c
          nproj = nglb
          if (proj.eq.'MO') nproj = nmolw
          do 150 iproj=1,nproj
            if (proj.eq.'CE') then
              plat = cenglb(1,iproj)
              plon = cenglb(2,iproj)
            else
              plat = cenmolw(1,iproj)
              plon = cenmolw(2,iproj)
            endif
            if (plon.eq.r12flag)
     +        ixlon = ixslt(12.,utc,plon,gcmlon,imx,dlon)
            call mkproj(proj,viewport,plat,plon,0.,dum,dum)
c
c Draw map and contours: 
c If monochrome, add contours after map so continental outlines and
c   map grid will not go through line labels. If color, however, contours
c   must be last so they do not cover the map.
c
            if (ip.eq.ixunvn.or.ip.eq.ixuivi) then
              if (icont.gt.0) call maplot
              if (plat.ne.0..or.proj.eq.'MO') then
                call mapgrd
                call maplbl
              endif
              if (ip.eq.ixunvn) goto 205
              if (ip.eq.ixuivi) goto 210
            endif
            if (icolor.le.0) then
              call polyclr(iwhite)
              if (icont.gt.0) call maplot
              if (plat.ne.0..or.proj.eq.'MO') then
                call mapgrd
                call maplbl
              endif
              call contour(plt,imx,imx,jmx,cint(ip),cmin(ip),cmax(ip))
            else
              call polyclr(iblack)
              call conclr(plt,imx,imx,jmx,cint(ip),cmin(ip),cmax(ip))
              if (icont.gt.0) call maplot
              if (plat.ne.0..or.proj.eq.'MO') then
                call mapgrd
                call maplbl
              endif
            endif
c
c Plot un+vn vectors if contouring tn (black if color):
c
 205        continue
            if (ip.eq.ixunvn.or.(ip.eq.ixt.and.ituv.gt.0)) then
              ipercen = 0
              if (idif.eq.1.and.iperc(ixu).eq.1.and.iperc(ixv).eq.1) 
     +          ipercen = 1
              call getglb(pert,udat,ixpl,ht,ixu,1)
              call getglb(cntr,vdat,ixpl,ht,ixu,0)
              call mkdifs(udat,vdat,imx*jmx,ipercen,cpspval)
              call getglb(pert,vdat,ixpl,ht,ixv,1)
              call getglb(cntr,glb,ixpl,ht,ixv,0)
              call mkdifs(vdat,glb,imx*jmx,ipercen,cpspval)
              do i=1,imx
                do j=1,jmx
                  if (udat(i,j).le.1.e-20.and.udat(i,j).ge.-1.e-20.and.
     +                vdat(i,j).le.1.e-20.and.vdat(i,j).ge.-1.e-20) then
                    udat(i,j) = cpspval
                    vdat(i,j) = cpspval
                  endif
                enddo
              enddo
c             uvmax = 0.
              m = 1
              ispvec = 1
              spvec = cpspval
              iskip = 0
              ixcolor = iwhite
              if (icolor.gt.0.and.ip.ne.ixunvn) then
                call polyclr(iblack)
                ixcolor = iblack
              endif
              if (ip.eq.ixunvn) then
                if (icolor.le.0) call polyclr(iwhite)
                if (icont.gt.0) call maplot
              endif
c             call getset(vl,vr,vb,vt,wl,wr,wb,wt,l)
c             call pltvec(udat,vdat,uvmax,imx,jmx,jmx,iskip,gcmlat(1),
c    +          gcmlat(jmx),m,ispvec,spvec,ixcolor)
c             call set(vl,vr,vb,vt,wl,wr,wb,wt,l)
              incx = 3
              incy = 3
              vlc = 0.
              vhc = uvmax
              call pltvect(udat,vdat,imx,jmx,gcmlon(1),gcmlon(imx),
     +          gcmlat(1),gcmlat(jmx),incx,incy,vlaboff,vlc,vhc,cpspval)
              if (iwrascii.gt.0.and.ip.eq.ixunvn) then
                do j=1,jmx
                  do i=1,imx
                    if (udat(i,j).ne.cpspval.and.vdat(i,j).ne.cpspval)
     +                then
                      plt(i,j) = sqrt(udat(i,j)**2+vdat(i,j)**2)
                    else
                      plt(i,j) = cpspval
                    endif
                  enddo
                enddo
              endif
            endif
c
c Plot ui+vi vectors if contouring potential (black if color):
c
 210        continue
            if (ip.eq.ixuivi.or.(ip.eq.ixpot.and.ipuv.gt.0)) then
              ipercen = 0
              if (idif.eq.1.and.iperc(ixui).eq.1.and.iperc(ixvi).eq.1) 
     +          ipercen = 1
              call getglb(pert,udat,ixpl,ht,ixui,1)
              call getglb(cntr,vdat,ixpl,ht,ixui,0)
              call mkdifs(udat,vdat,imx*jmx,ipercen,cpspval)
              call getglb(pert,vdat,ixpl,ht,ixvi,1)
              call getglb(cntr,glb,ixpl,ht,ixvi,0)
              call mkdifs(vdat,glb,imx*jmx,ipercen,cpspval)
              do i=1,imx
                do j=1,jmx
                  if (udat(i,j).le.1.e-20.and.udat(i,j).ge.-1.e-20.and.
     +                vdat(i,j).le.1.e-20.and.vdat(i,j).ge.-1.e-20) then
                    udat(i,j) = cpspval
                    vdat(i,j) = cpspval
                  endif
                enddo
              enddo
c             uvmax = 0.
              m = 1
              ispvec = 1
              spvec = cpspval
              iskip = 0
              ixcolor = iwhite
              if (icolor.gt.0) then
                call polyclr(iblack)
                ixcolor = iblack
              endif
c             call getset(vl,vr,vb,vt,wl,wr,wb,wt,l)
c             call pltvec(udat,vdat,uvmax,imx,jmx,jmx,iskip,gcmlat(1),
c    +          gcmlat(jmx),m,ispvec,spvec,ixcolor)
c             call set(vl,vr,vb,vt,wl,wr,wb,wt,l)
              incx = 3
              incy = 3
              vlc = 0.
              vhc = uvmax
              call pltvect(udat,vdat,imx,jmx,gcmlon(1),gcmlon(imx),
     +          gcmlat(1),gcmlat(jmx),incx,incy,vlaboff,vlc,vhc,cpspval)
              if (iwrascii.gt.0.and.ip.eq.ixuivi) then
                do j=1,jmx
                  do i=1,imx
                    if (udat(i,j).ne.cpspval.and.vdat(i,j).ne.cpspval)
     +                then
                      plt(i,j) = sqrt(udat(i,j)**2+vdat(i,j)**2)
                    else
                      plt(i,j) = cpspval
                    endif
                  enddo
                enddo
              endif
            endif
c
c Add local time axis:
c (fake out sltxax with local time if plon.ne.0. (it thinks center lon=0.)
c
            if (icolor.gt.0) call polyclr(iwhite)
            if (proj.ne.'MO') then
              if (plat.eq.0..and.plon.eq.0.) then
                call sltxax(utc)
              elseif (plat.eq.0.) then
                rslt = utc + plon/15.
                if (rslt.lt.0.) rslt = rslt+24.
                call sltxax(rslt)
              endif
            endif
c
c Label axes:
c
            if (plat.eq.0..and.proj.ne.'MO') then
              if (plon.eq.0.) then
                call labrect(gcmlon,imx,gcmlat,jmx,'LONGITUDE',
     +            'LATITUDE',.04)
              else
                rlon1= plon
                if (rlon1.le.-180.) rlon1 = rlon1+360.
                lonlab(5) = ifix(rlon1)
                rlon = rlon1
                do i=1,4
                  rlon = rlon-45.
                  if (rlon.le.-180.) rlon = rlon+360.
                  lonlab(5-i) = ifix(rlon)
                enddo
                rlon = rlon1
                do i=1,4
                  rlon = rlon+45.
                  if (rlon.gt.180.) rlon = rlon-360.
                  lonlab(5+i) = ifix(rlon)
                enddo
                call labaxes(9,lonlab,3,'LONGITUDE',
     +                       7,latlab,3,'LATITUDE', 0.)
              endif
            endif
c
c Field label at top:
c
            call clearstr(toplab)
            if (log10map.gt.0.and.ilog(ip).gt.0) then
              len = lenstr(flab(ip))
              if (len+6.le.56) then
                write(toplab,"('LOG10 ',a)") flab(ip)
              endif
            else
              toplab = flab(ip)
            endif
            call wrlab(toplab(1:lenstr(toplab)),xmid,
     +        viewport(4)+.11,0.)
            call clearstr(fieldlab)
            write(fieldlab,"(a)") toplab(1:lenstr(toplab))
c
c Difference label:
c
            call clearstr(toplab)
            if (ipercent.gt.0) then
              write(toplab,"('PERCENT DIFFERENCE FIELD')")
            else
              write(toplab,"('RAW DIFFERENCE FIELD')")
            endif
            call wrlab(toplab(1:lenstr(toplab)),xmid,
     +        viewport(4)+.07,0.)
c
c Label ut, zp, ht at top below field label:
c
            call clearstr(toplab)
            if (izp.le.npls) then
              if (ip.ne.ixfof2.and.ip.ne.ixhmf2) then
                write(toplab,"('ZP=',f5.1,' UT=',f6.2)")
     +            spls(izp),utc
              else
                write(toplab,"('UT=',f6.2)") utc
              endif
            else
              if (ip.ne.ixfof2.and.ip.ne.ixhmf2) then
                write(toplab,"('HT=',f5.1,' UT=',f6.2)") ht,utc
              else
                write(toplab,"('UT=',f6.2)") utc
              endif
            endif
            call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+.03,0.)
            call clearstr(rec80)
            write(rec80,"(a)") toplab(1:lenstr(toplab))
c
c Label history volume and center of projection at bottom:
c
            call clearstr(toplab)
            call tail(pertvol(ipvol),flnm0)
            lenpert = lenstr(flnm0)
            call tail(cntrvol(icvol),flnm1)
            lencntr = lenstr(flnm1)
            if (lenpert+lencntr+18.le.56) then
              write(toplab,"('DIFFS OF ',a,' MINUS ',a)") 
     +          flnm0(1:lenpert),flnm1(1:lencntr)
              call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.17,
     +          .012)
            endif
            call clearstr(histlab)
            write(histlab,"(a)") toplab(1:lenstr(toplab))
            call clearstr(toplab)
            write(toplab,"('CENTER OF PROJECTION = ',f5.1,',',f6.1)") 
     +        plat,plon
            call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.20,
     +        .012)
c
c Wrap it up:
c
            call frame
            iframe = iframe+1
            if (izp.le.npls) then
              if (ip.ne.ixfof2.and.ip.ne.ixhmf2) then
                write(6,"('pltglb: frame ',i4,' field ',a8,' zp=',f5.1,
     +            ' proj=',a,' cenproj=',f6.2,',',f7.2)") iframe,
     +            labshort(ip),spls(izp),proj,plat,plon
              else
                write(6,"('pltglb: frame ',i4,' field ',a8,
     +            ' proj=',a,' cenproj=',f6.2,',',f7.2)")
     +            iframe,labshort(ip),proj,plat,plon
              endif
            else
              if (ip.ne.ixfof2.and.ip.ne.ixhmf2) then
                write(6,"('pltglb: frame ',i4,' field ',a8,' ht=',f5.1,
     +            ' proj=',a,' cenproj=',f6.2,',',f7.2)")
     +            iframe,labshort(ip),ht,proj,plat,plon
              else
                write(6,"('pltglb: frame ',i4,' field ',a8,
     +            'proj=',a,' cenproj=',f6.2,',',f7.2)")
     +            iframe,labshort(ip),proj,plat,plon
              endif
            endif
c
c Write ascii data file if desired:
            if (iwrascii.gt.0) then
              call wrascii(iwrascii,luascii,plt,imx,jmx,gcmlon,gcmlat,
     +          'LONGITUDE','LATITUDE',histlab,fieldlab,
     +          iframe,rec80,'tigcmdif',dirascii)
            endif
c
c End projection loop
 150      continue
c
c End selected pressure/height loop: izp=1,npls+nhts
 100    continue
c
c End field loop
 200  continue
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine polyclr(iclr)
      call sflush
      call gsplci(iclr)
      call gspmci(iclr)
      return
      end