c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcmproc/pltglb.f c------------------------------------------------------------------ c subroutine pltglb(proj,utc,fields) c c Contour on cylindrical equidistant projection with longitude on c x-axis, latitude on y-axis, at selected zp/ht: c include 'tigcmproc.h' character*56 toplab,fieldlab character*80 rec80 character*2 proj dimension plt(imx,jmx),viewport(4),fields(imx,kmx,jmx,nfget), + udat(imx,jmx),vdat(imx,jmx),lonlab(9),latlab(7) 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',-.45) call cpsetr('ILS',.016) call cpseti('ILP',0) call mapstr('GR',dgrid) vlaboff = -.8 c if (isltxax.le.0) then c rily = -.25 c vlaboff = -.46 c endif c c Field loop: write(6,"(' ')") do 100 ip=1,nftot if (ifplt(ip).le.0) goto 100 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 Selected pressure/height loop: do 200 izp = 1,npls+nhts if (izp.gt.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) goto 100 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 200 endif else ht = shts(izp-npls) endif c c Set up map and get data: c if (ip.eq.ixunvn.or.ip.eq.ixuivi) then if (icolor.le.0) call polyclr(iwhite) if (plat.ne.0..or.proj.eq.'MO') then call mapgrd call maplbl endif if (icont.gt.0) call maplot if (ip.eq.ixunvn) goto 105 if (ip.eq.ixuivi) goto 110 else call getglb(fields,plt,ixpl,ht,ip) endif 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 (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 (if doing ascii files for unvn, save vector sums) c 105 continue if (ip.eq.ixunvn.or.(ip.eq.ixt.and.ituv.gt.0)) then call getglb(fields,udat,ixpl,ht,ixu) call getglb(fields,vdat,ixpl,ht,ixv) m = 1 ispvec = 0 spvec = 1.e36 iskip = 0 ixcolor = iwhite if (icolor.gt.0.and.ip.ne.ixunvn) then call polyclr(iblack) ixcolor = iblack endif 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 (if doing ascii files for uivi, save vector sums) c 110 continue if (ip.eq.ixuivi.or.(ip.eq.ixpot.and.ipuv.gt.0)) then call getglb(fields,udat,ixpl,ht,ixui) call getglb(fields,vdat,ixpl,ht,ixvi) m = 1 ispvec = 0 spvec = 1.e36 iskip = 0 ixcolor = iwhite if (icolor.gt.0.and.ip.ne.ixuivi) then call polyclr(iblack) ixcolor = iblack endif 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 write(toplab,"('LOG10 ',a)") flab(ip)(1:lenstr(flab(ip))) else toplab = flab(ip) endif call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+.08,0.) fieldlab = toplab 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) write(toplab,"('HISTORY=',a)") histvol(ivol) call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.15, + .012) call clearstr(toplab) write(toplab,"('CENTER OF PROJECTION = ',f5.1,',',f6.1)") + plat,plon call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.18, + .012) c c Wrap it up: c c call box(1) 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',histvol(ivol),fieldlab, + iframe,rec80,'tigcmproc',dirascii) endif c c End selected pressure/height loop: izp=1,npls+nhts 200 continue c c End projection loop 150 continue c c End field loop 100 continue c return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine polyclr(iclr) call sflush call gsplci(iclr) call gspmci(iclr) return end