c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcmdif/pltsatv.f c------------------------------------------------------------------ c subroutine pltsatv(utc,pert,cntr) c c Contour on sat view projections at selected zp and/or hts: c include 'tigcmdif.h' character*56 toplab,fieldlab character*80 rec80,histlab character*2 outline character*24 flnm0,flnm1 dimension plt(imx,jmx),viewport(4),udat(imx,jmx),vdat(imx,jmx) dimension pert(imx,kmx,jmx,nfget),cntr(imx,kmx,jmx,nfget) data viewport /.15,.85,.15,.85/, dum/0./, eradii/6.631/, + dgrid/30./ 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',-.07) call cpsetr('ILS',.018) call cpseti('ILP',0) call mapstr('GR',dgrid) c c Field loop: c Note getglb handles selected pressure vs selected constant height c and does interpolation in case of the latter. c 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,"('>>> pltsatv: 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,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 (note ixslt changes plon) c do 300 iproj=1,nsatv plat = censatv(1,iproj) plon = censatv(2,iproj) if (plon.eq.r12flag) + ixlon = ixslt(12.,utc,plon,gcmlon,imx,dlon) call mkproj('SV',viewport,plat,plon,0.,dum,eradii) if (ip.eq.ixunvn.or.ip.eq.ixuivi) then if (icont.gt.0) call maplot call mapgrd call maplbl if (ip.eq.ixunvn) goto 205 if (ip.eq.ixuivi) goto 210 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 first so they do not cover the map. c if (icolor.le.0) then call polyclr(iwhite) if (icont.gt.0) call maplot call mapgrd call maplbl 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 call mapgrd call maplbl endif 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,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,plt,ixpl,ht,ixv,0) call mkdifs(vdat,plt,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 = 0 spvec = 1.e36 iskip = 0 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,jmx,iskip,gcmlat(1), + gcmlat(jmx),m,ispvec,spvec,ixcolor) call set(vl,vr,vb,vt,wl,wr,wb,wt,l) 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 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,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,plt,ixpl,ht,ixvi,0) call mkdifs(vdat,plt,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 = 0 spvec = 1.e36 iskip = 0 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,jmx,iskip,gcmlat(1), + gcmlat(jmx),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=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 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, + viewport(4)+.11,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, + 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 FROM ',a,' MINUS ',a)") + flnm0(1:lenpert),flnm1(1:lencntr) call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.09, + .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.13, + .012) c c Wrap it up: c call frame iframe = iframe+1 if (izp.le.npls) then write(6,"('pltsatv: frame ',i4,' field ',a8,' zp=',f5.1, + ' censatv=',2f9.2)") iframe,labshort(ip),spls(izp), + plat,plon else write(6,"('pltsatv: frame ',i4,' field ',a8,' ht=',f5.1, + ' censatv=',2f9.2)") iframe,labshort(ip),ht,plat, + plon 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: iproj=1,nsatv 300 continue c c End selected pressure/height loop: izp=1,npls+nhts 200 continue c c End field loop: ip=1,nftot 100 continue c return end