c c------------------------------------------------------------------ c subroutine pltlat(utc,fields) 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 'tigcmproc.h' character*56 toplab,fieldlab character*80 rec80 dimension plt(imx,kmx),viewport(4),yaxht(kmx), + plt0(imx,kmx),fields(imx,kmx,jmx,nfget),hts(imx,kmx) pointer(ppltht,pltht(1)) data viewport /.15,.89,.26,.91/ 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) c c Field loop: c write(6,"(' ')") do 100 ip=1,nftot if (ifplt(ip).le.0) goto 100 if (ip.eq.ixfof2.or.ip.eq.ixhmf2) goto 100 if (ip.eq.ixunvn.or.ip.eq.ixuivi) goto 100 c c ihtsc = 1,2 for pressure,height on y-axis: c ib = 0 ie = 0 if (nzprange.le.0.and.nhtscale.gt.0) then ib = 2 ie = 2 elseif (nzprange.gt.0.and.nhtscale.gt.0) then ib = 1 ie = 2 elseif (nzprange.gt.0.and.nhtscale.le.0) then ib = 1 ie = 1 endif if (ib.eq.0.or.ie.eq.0) then write(6,"('>>> pltlat: nzprange=',i3,' nhtscale=',i3, + ' need vertical scale for lon plots -- returning')") + nzprange,nhtscale return endif do 150 iht = ib,ie ihtsc = iht-1 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 alloc(ppltht,imx*nhtscale) endif c c Selected latitude loop: c do 200 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 200 endif rlat = gcmlat(ixlat) c c Get field and interpolate to height scale if necessary: c call getlat(fields,plt0,ixlat,ilog(ip),ip) if (ihtsc.le.0) then call getlat(fields,plt0,ixlat,ilog(ip),ip) 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(fields,plt0,ixlat,0,ip) call getlat(fields,hts,ixlat,0,ixz) call cuthtint(plt0,hts,imx,kmx,pltht,htscale,nhtscale, + ilog(ip),cpspval,ier,0) if (ilog(ip).gt.0) call log10f(pltht,imx*nhtscale,1.e-20, + 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 call clearstr(toplab) if (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)+.07,0.) fieldlab = toplab c c Lat, ut at top below field label: c call clearstr(toplab) 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 call clearstr(toplab) write(toplab,"('HISTORY=',a)") histvol(ivol) call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.15, + .012) 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) + fields(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 (write ascii data file if needed): c 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)', + histvol(ivol),fieldlab,iframe,rec80,'tigcmproc', + 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)',histvol(ivol), + fieldlab,iframe,rec80,'tigcmproc',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 200 continue c c Height/pressure: if (ihtsc.gt.0) call hpdeallc(ppltht,ier,1) 150 continue c c End field loop: ip=1,nftot 100 continue c return end