c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcmdif/pltpol.f c------------------------------------------------------------------ c subroutine pltpol(utc,pert,cntr) c include 'tigcmdif.h' character*56 toplab,fieldlab character*80 rec80,histlab character*1 hem character*24 flnm0,flnm1 dimension plt(imx,jmx),vp(4),udat(imx,jmx), + vdat(imx,jmx),plt0(imx,jmx),ylat(jmx),incx(imx),incy(jmx) dimension pert(imx,kmx,jmx,nfget),cntr(imx,kmx,jmx,nfget) data vp /.150,.850,.153,.853/, dum/0./ data iblack/0/, iwhite/1/ data incx/imx*3/ c xmid = 0.5*(vp(1)+vp(2)) call cpseti('SET',0) call cpseti('MAP',1) call cpsetr('XC1',gcmlon(1)) call cpsetr('XCM',gcmlon(imx)) call cpsetr('ILX',xmid) call cpsetr('ILY',-.12) call cpsetr('ILS',.016) call cpseti('ILP',0) c c Field loop: write(6,"(' ')") do 100 ip=1,nftot if (ifplt(ip).le.0) goto 100 ipercent = 0 if (idif.eq.1.and.iperc(ip).eq.1) ipercent = 1 c c Selected pressure/height loop: do 200 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,"('>>> pltpol: bad selected pressure=',f10.3, + ' izp=',i3)") spls(izp),izp goto 200 endif else ht = shts(izp-npls) endif if (izp.gt.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) goto 100 if (ip.eq.ixz.and.izp.gt.npls) goto 100 if (ip.ne.ixunvn.and.ip.ne.ixuivi) then call getglb(pert,plt0,ixpl,ht,ip,1) call getglb(cntr,plt,ixpl,ht,ip,0) call mkdifs(plt0,plt,imx*jmx,ipercent,cpspval) endif c c Perimeter latitude loop: c do 300 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. do i=1,jmx-4 incy(i) = 2 enddo incy(jmx-3) = 4 incy(jmx-2) = 6 incy(jmx-1) = 8 incy(jmx) = 24 else hem = 'S' yc1 = gcmlat(1) ycn = perimlat(ih) iskip = 0 rot = utc*15.-180. plat = -90. incy(1) = 24 incy(2) = 8 incy(3) = 6 incy(4) = 4 do i=5,jmx incy(i) = 2 enddo endif call mkproj('ST',vp,plat,0.,rot,perimlat(ih),dum) call cpsetr('YC1',yc1) call cpsetr('YCN',ycn) do j=iskip+1,jmx plt(:,j-iskip) = plt0(:,j) ylat(j-iskip) = gcmlat(j) enddo if (ip.eq.ixunvn.or.ip.eq.ixuivi) then if (icont.gt.0) call maplot call labpolar(vp(1),vp(2),vp(3),vp(4),hem) if (ip.eq.ixunvn) goto 205 if (ip.eq.ixuivi) goto 210 endif c c Draw contours and map outlines: c if (icolor.le.0) then call polyclr(iwhite) if (icont.gt.0) call maplot call contour(plt,imx,imx,nlats,cint(ip),cmin(ip),cmax(ip)) else call polyclr(iblack) call conclr(plt,imx,imx,nlats,cint(ip),cmin(ip),cmax(ip)) if (icont.gt.0) call maplot endif if (icolor.gt.0) call polyclr(iwhite) call labpolar(vp(1),vp(2),vp(3),vp(4),hem) c c Add un+vn vectors if plotting 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,plt0,ixpl,ht,ixu,1) call getglb(cntr,plt,ixpl,ht,ixu,0) call mkdifs(plt0,plt,imx*jmx,ipercen,cpspval) do j=iskip+1,jmx udat(:,j-iskip) = plt0(:,j) enddo call getglb(pert,plt0,ixpl,ht,ixv,1) call getglb(cntr,plt,ixpl,ht,ixv,0) call mkdifs(plt0,plt,imx*jmx,ipercen,cpspval) do j=iskip+1,jmx vdat(:,j-iskip) = plt0(:,j) enddo do i=1,imx do j=1,nlats 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 = 2 if (hem.eq.'N') m = 3 ispvec = 0 spvec = 1.e36 ixcolor = 1 if (icolor.gt.0) then call polyclr(iblack) ixcolor = 0 endif scarrow = 120. call getset(vl,vr,vb,vt,wl,wr,wb,wt,l) call pltvec(udat,vdat,uvmax,imx,jmx,nlats,iskip,yc1,ycn, + m,ispvec,spvec,ixcolor,scarrow) call set(vl,vr,vb,vt,wl,wr,wb,wt,l) c c subroutine pltvect(u,v,id1,id2,x0,x1,y0,y1,incx,incy,ylaboff, c + vlc,vhc,spv) c c vlaboff = -.15 c call pltvect(udat,vdat,imx,nlats,gcmlon,gcmlon(imx), c + ylat,ylat(nlats),incx,incy,ylaboff,0.,0.,cpspval) if (iwrascii.gt.0.and.ip.eq.ixunvn) then do j=iskip+1,jmx do i=1,imx if (udat(i,j-iskip).ne.cpspval.and. + vdat(i,j-iskip).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 ui+vi vectors if plotting poten (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,plt0,ixpl,ht,ixui,1) call getglb(cntr,plt,ixpl,ht,ixui,0) call mkdifs(plt0,plt,imx*jmx,ipercen,cpspval) do j=iskip+1,jmx udat(:,j-iskip) = plt0(:,j) enddo call getglb(pert,plt0,ixpl,ht,ixvi,1) call getglb(cntr,plt,ixpl,ht,ixvi,0) call mkdifs(plt0,plt,imx*jmx,ipercen,cpspval) do j=iskip+1,jmx vdat(:,j-iskip) = plt0(:,j) enddo do i=1,imx do j=1,nlats 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 = 2 if (hem.eq.'N') m = 3 ispvec = 0 spvec = 1.e36 ixcolor = 1 if (icolor.gt.0) then call polyclr(iblack) ixcolor = 0 endif call getset(vl,vr,vb,vt,wl,wr,wb,wt,l) call pltvec(udat,vdat,uvmax,imx,jmx,nlats,iskip,yc1,ycn, + m,ispvec,spvec,ixcolor) call set(vl,vr,vb,vt,wl,wr,wb,wt,l) if (iwrascii.gt.0.and.ip.eq.ixuivi) then do j=iskip+1,jmx do i=1,imx if (udat(i,j-iskip).ne.cpspval.and. + vdat(i,j-iskip).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 Field label at top: c if (icolor.gt.0) call polyclr(iwhite) 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,vp(4)+.13,0.) fieldlab = 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,vp(4)+.09,0.) c c ut, zp, perimlat 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,' PERIMLAT=', + f5.1)") spls(izp),utc,perimlat(ih) else write(toplab,"('UT=',f6.2,' PERIMLAT=',f5.1)") + utc,perimlat(ih) endif else if (ip.ne.ixfof2.and.ip.ne.ixhmf2) then write(toplab,"('HT=',f5.1,' UT=',f6.2,' PERIMLAT=', + f5.1)") ht,utc,perimlat(ih) else write(toplab,"('UT=',f6.2,' PERIMILAT=',f5.1)") + utc,perimlat(ih) endif endif call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.05,0.) call clearstr(rec80) write(rec80,"(a)") toplab(1:lenstr(toplab)) c c History volume label 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 FROM ',a,' MINUS ',a)") + flnm0(1:lenpert),flnm1(1:lencntr) call wrlab(toplab(1:lenstr(toplab)),xmid,vp(3)-0.14,.012) endif call clearstr(histlab) write(histlab,"(a)") toplab(1:lenstr(toplab)) 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,"('pltpol: frame ',i4,' field ',a8,' zp=',f5.1, + ' perimlat=',f6.1)") + iframe,labshort(ip),spls(izp),perimlat(ih) else write(6,"('pltpol: frame ',i4,' field ',a8,' perimlat=', + f6.1)") iframe,labshort(ip),perimlat(ih) endif else if (ip.ne.ixfof2.and.ip.ne.ixhmf2) then write(6,"('pltpol: frame ',i4,' field ',a8,' ht=',f5.1, + ' perimlat=',f6.1)") iframe,labshort(ip),ht, + perimlat(ih) else write(6,"('pltpol: frame ',i4,' field ',a8,f5.1, + ' perimlat=',f6.1)") iframe,labshort(ip),perimlat(ih) endif endif c c Write ascii data file if desired: if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,plt,imx,nlats,gcmlon, + ylat,'LONGITUDE','LATITUDE',histlab,fieldlab, + iframe,rec80,'tigcmdif',dirascii) endif c c End perimeter latitude loop 300 continue c c End selected pressure/height loop: izp=1,npls+nhts 200 continue c c End field loop 100 continue return end