c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcmdif/pltlat.f c------------------------------------------------------------------ c subroutine pltlat(utc,ihtsc,pert,cntr) c c Contour latitude slices: c (lon on x-axis, zp on y-axis if ihtsc <=0, c linear height scale if ihtsc > 0) c include 'tigcmdif.h' character*56 toplab character*80 rec80,histlab,fieldlab character*24 flnm0,flnm1 dimension plt(imx,kmx),viewport(4),yaxht(kmx),plt0(imx,kmx), + hts(imx,kmx) dimension pert(imx,kmx,jmx,nfget),cntr(imx,kmx,jmx,nfget) pointer(ppltht,pltht(1)), (ppltht0,pltht0(1)) data viewport /.15,.89,.20,.85/ c if (nlat.le.0) then write(6,"('pltlat: no selected latitudes -- returning')") return endif xmid = 0.5*(viewport(1)+viewport(2)) call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',gcmlon(1)) call cpsetr('XCM',gcmlon(imx)) call cpsetr('ILX',xmid) call cpsetr('ILY',-.16) call cpseti('ILP',0) call cpsetr('ILS',.016) if (ihtsc.le.0) then call cpsetr('YC1',gcmzp(izprange(1))) call cpsetr('YCN',gcmzp(izprange(2))) call set(viewport(1),viewport(2),viewport(3),viewport(4), + gcmlon(1),gcmlon(imx),gcmzp(izprange(1)),gcmzp(izprange(2)),1) else call cpsetr('YC1',htscale(1)) call cpsetr('YCN',htscale(nhtscale)) call set(viewport(1),viewport(2),viewport(3),viewport(4), + gcmlon(1),gcmlon(imx),htscale(1),htscale(nhtscale),1) call hpalloc(ppltht,imx*nhtscale,ier,1) if (ier.ne.0) then write(6,"('pltlat: hpalloc error = ',i6,' for pltht')") ier stop 'ppltht' endif call hpalloc(ppltht0,imx*nhtscale,ier,1) if (ier.ne.0) then write(6,"('pltlat: hpalloc error = ',i6,' for pltht0')") ier stop 'ppltht0' endif endif c c Field loop: c write(6,"(' ')") do 200 ip=1,nftot if (ifplt(ip).le.0) goto 200 if (ip.eq.ixfof2.or.ip.eq.ixhmf2) goto 200 if (ip.eq.ixuivi.or.ip.eq.ixunvn) goto 200 ipercent = 0 if (idif.eq.1.and.iperc(ip).eq.1) ipercent = 1 c c Selected latitude loop: c do 100 j=1,nlat ixlat = ixfind(gcmlat,jmx,slat(j),dlat) if (ixlat.le.0) then write(6,"('>>> pltlat: bad latitude=',f10.3,' j=',i3, + ' skipping this lat')") slat(j),j goto 100 endif rlat = gcmlat(ixlat) if (ihtsc.le.0) then call getlat(pert,plt0,rlat,ilog(ip),ip,1) call getlat(cntr,plt,rlat,ilog(ip),ip,0) call mkdifs(plt0,plt,imx*kmx,ipercent,cpspval) do k=izprange(1),izprange(2) kk = k-izprange(1)+1 plt(:,kk) = plt0(:,k) enddo ny = izprange(2)-izprange(1)+1 else call getlat(pert,plt,rlat,ilog(ip),ip,1) call getlat(pert,hts,rlat,0,ixz,1) call cuthtint(plt,hts,imx,kmx,pltht,htscale,nhtscale, + ilog(ip),cpspval,ier,1) call getlat(cntr,plt0,rlat,ilog(ip),ip,0) call getlat(cntr,hts,rlat,0,ixz,0) call cuthtint(plt0,hts,imx,kmx,pltht0,htscale,nhtscale, + ilog(ip),cpspval,ier,1) call mkdifs(pltht,pltht0,imx*nhtscale,ipercent,cpspval) ny = nhtscale endif c c Contour: c if (icolor.le.0) then if (ihtsc.le.0) then call contour(plt,imx,imx,ny,cint(ip),cmin(ip),cmax(ip)) else call contour(pltht,imx,imx,ny,cint(ip),cmin(ip),cmax(ip)) endif else if (ihtsc.le.0) then call conclr(plt,imx,imx,ny,cint(ip),cmin(ip),cmax(ip)) else call conclr(pltht,imx,imx,ny,cint(ip),cmin(ip),cmax(ip)) endif endif c c Add axes labels: c if (ihtsc.le.0) then call labrect(gcmlon,imx,gcmzp(izprange(1)),ny,'LONGITUDE', + 'ZP',0.) else call labrect(gcmlon,imx,htscale(1),ny,'LONGITUDE', + 'HEIGHT (KM)',0.) endif c c Field label at top: c do ii=1,56 toplab(ii:ii) = ' ' enddo if (ilog(ip).gt.0) then len = lenstr(flab(ip)) if (len+6.le.56) write(toplab,"('LOG10 ',a)") flab(ip) 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 do ii=1,56 toplab(ii:ii) = ' ' enddo 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 Lat, ut at top below field label: c do ii=1,56 toplab(ii:ii) = ' ' enddo write(toplab,"('LAT=',f8.2,' UT=',f6.2)") rlat,utc call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.03,0.) call clearstr(rec80) write(rec80,"(a)") toplab(1:lenstr(toplab)) c c History label at bottom: c do ii=1,56 toplab(ii:ii) = ' ' enddo call tail(pertvol(ipvol),flnm0) lenpert = lenstr(flnm0) call tail(cntrvol(icvol),flnm1) lencntr = lenstr(flnm1) write(toplab,"('DIFFS FROM ',a,' MINUS ',a)") + flnm0(1:lenpert),flnm1(1:lencntr) call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.15, + .012) call clearstr(histlab) write(histlab,"(a)") toplab(1:lenstr(toplab)) c c Put height axis on right if pressure on y-axis c (use zonally averaged heights): c if (ixfld(ixz).le.0.or.ihtsc.gt.0) goto 300 yaxht(:) = 0. do k=izprange(1),izprange(2) izp = k-izprange(1)+1 do i=1,imx-1 yaxht(izp) = yaxht(izp) + .5*(pert(i,k,ixlat,ixfld(ixz))+ + cntr(i,k,ixlat,ixfld(ixz))) enddo yaxht(izp) = yaxht(izp) / float(imx-1) enddo rnd = 10. if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5. call altyax(ny,yaxht,gcmzp(izprange(1)),rnd,6) c c Wrap it up: 300 continue call frame iframe = iframe+1 if (ihtsc.le.0) then if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,plt,imx,ny,gcmlon, + gcmzp(izprange(1)),'LONGITUDE','LN(P0/P)', + histlab,fieldlab,iframe,rec80,'tigcmdif', + dirascii) endif write(6,"('pltlat: frame ',i4,' field ',a8,' lat=',f8.2, + ' zprange=',2f8.2)") iframe,labshort(ip),rlat,zprange else if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,pltht,imx,ny,gcmlon, + htscale,'LONGITUDE','HEIGHT (KM)',histlab, + fieldlab,iframe,rec80,'tigcmdif',dirascii) endif write(6,"('pltlat: frame ',i4,' field ',a8,' lat=',f8.2, + ' ht scale=',f7.2,' to ',f7.2)") iframe, + labshort(ip),rlat,htscale(1),htscale(nhtscale) endif c c End selected latitude loop 100 continue c c End field loop: ip=1,nftot 200 continue if (ihtsc.gt.0) then call hpdeallc(ppltht,ier,1) call hpdeallc(ppltht0,ier,1) endif c return end