c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcmdif/pltlon.f c------------------------------------------------------------------ c subroutine pltlon(utc,ihtsc,pert,cntr) c c Contour longitude slices (or zonal means): c (lat on x-axis, zp on y-axis if ihtsc <=0, c linear height scale if ihtsc > 0) c include 'tigcmdif.h' character*56 toplab,fieldlab character*24 flnm0,flnm1 character*80 rec80,histlab dimension pert(imx,kmx,jmx,nfget), cntr(imx,kmx,jmx,nfget) dimension viewport(4),yaxht(kmx),rimx(imx),plt(jmx,kmx), + plt0(jmx,kmx),hts(jmx,kmx) pointer(ppltht,pltht(1)), (ppltht0,pltht0(1)) data viewport /.15,.89,.20,.85/ c if (nlon+nslt.le.0) then write(6,"('pltlon: no selected longitudes -- returning')") return endif xmid = 0.5*(viewport(1)+viewport(2)) call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',gcmlat(1)) call cpsetr('XCM',gcmlat(jmx)) 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), + -90.,90.,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), + -90.,90.,htscale(1),htscale(nhtscale),1) call hpalloc(ppltht,jmx*nhtscale,ier,1) if (ier.ne.0) then write(6,"('pltlat: hpalloc error = ',i6,' for pltht')") ier stop 'ppltht' endif call hpalloc(ppltht0,jmx*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 100 ip=1,nftot if (ifplt(ip).le.0) goto 100 if (ip.eq.ixfof2.or.ip.eq.ixhmf2) goto 100 if (ip.eq.ixz.and.ihtsc.gt.0) goto 100 if (ip.eq.ixuivi.or.ip.eq.ixunvn) goto 100 ipercent = 0 if (idif.eq.1.and.iperc(ip).eq.1) ipercent = 1 c c Selected longitude/local time loop: c do 200 i=1,nlon+nslt if (i.le.nlon) then if (slon(i).ne.zmflag) then ixlon = ixfind(gcmlon,imx,slon(i),dlon) if (ixlon.le.0) then write(6,"('>>> pltlon: bad longitude=',f10.3,' i=',i3, + ' skipping this lon')") slon(i),i goto 200 endif rlon = gcmlon(ixlon) else ixlon = ifix(zmflag) rlon = slon(i) endif c c Solar local time: else islt = i-nlon ixlon = ixslt(sslt(islt),utc,rlon,gcmlon,imx,dlon) endif write(6,"('pltlon: rlon=',f7.2)") rlon c c Zp on y-axis (no ht interp): c Note getlon() handles zonal mean or selected longitude, and c log10 or no log10 of field. c if (ihtsc.le.0) then call getlon(pert,plt0,rlon,0,ip,1) call getlon(cntr,plt,rlon,0,ip,0) call mkdifs(plt0,plt,jmx*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 c c Ht on y-axis (do ht interp separately on pert and cntr before mkdifs): c (interp on pert uses hts from pert, and same for cntr) c else call getlon(pert,plt,rlon,0,ip,1) call getlon(pert,hts,rlon,0,ixz,1) call cuthtint(plt,hts,jmx,kmx,pltht,htscale,nhtscale, + ilog(ip),cpspval,ier,1) call getlon(cntr,plt0,rlon,0,ip,0) call getlon(cntr,hts,rlon,0,ixz,0) call cuthtint(plt0,hts,jmx,kmx,pltht0,htscale,nhtscale, + ilog(ip),cpspval,ier,1) call mkdifs(pltht,pltht0,jmx*nhtscale,ipercent,cpspval) ny = nhtscale endif c c Contour: c if (icolor.le.0) then if (ihtsc.le.0) then call contour(plt,jmx,jmx,ny,cint(ip),cmin(ip),cmax(ip)) else call contour(pltht,jmx,jmx,ny,cint(ip),cmin(ip),cmax(ip)) endif else if (ihtsc.le.0) then call conclr(plt,jmx,jmx,ny,cint(ip),cmin(ip),cmax(ip)) else call conclr(pltht,jmx,jmx,ny,cint(ip),cmin(ip),cmax(ip)) endif endif c c Add axes labels: c if (ihtsc.le.0) then call labrect(gcmlat,jmx,gcmzp(izprange(1)),ny,'LATITUDE', + 'ZP',0.) else call labrect(gcmlat,jmx,htscale(1),ny,'LATITUDE', + 'HEIGHT (KM)',0.) endif c c Field label at top: c call clearstr(toplab) if (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 lon, ut at top below field label: c call clearstr(toplab) if (i.le.nlon) then if (slon(i).ne.zmflag) then write(toplab,"('LON=',f8.2,' UT=',f6.2)") rlon,utc else write(toplab,"('ZONAL MEANS UT=',f6.2)") utc endif else write(toplab,"('SLT=',f5.1,' (LON=',f6.1,') UT=',f6.2)") + sslt(islt),rlon,utc endif call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.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) 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 latitudinally and pert/cntr 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 if (slon(i).ne.zmflag) then do j=1,jmx yaxht(izp) = yaxht(izp) + 0.5* + (pert(ixlon,k,j,ixfld(ixz))+cntr(ixlon,k,j,ixfld(ixz))) enddo else ! zonal means -- use global mean hts do j=1,jmx rimx(:) = 0.5*(pert(:,k,j,ixfld(ixz))+ + cntr(:,k,j,ixfld(ixz))) yaxht(izp) = yaxht(izp) + + calcmean(rimx,imx,0,1.e-20,cpspval) enddo endif yaxht(izp) = yaxht(izp) / float(jmx) 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,jmx,ny,gcmlat, + gcmzp(izprange(1)),'LATITUDE','LN(P0/P)', + histlab,fieldlab,iframe,rec80,'tigcmdif', + dirascii) endif if (i.le.nlon) then if (slon(i).ne.zmflag) then write(6,"('pltlon: frame ',i4,' field ',a8,' lon=',f8.2, + ' zprange=',2f8.2)") iframe,labshort(ip),rlon,zprange else write(6,"('pltlon: frame ',i4,' field ',a8, + ' (zonal means) zprange=',2f8.2)") + iframe,labshort(ip),zprange endif else write(6,"('pltlon: frame ',i4,' field ',a8,' slt=',f8.2, + ' zprange=',2f8.2)") iframe,labshort(ip),sslt(islt), + zprange endif else if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,pltht,jmx,ny,gcmlat,htscale, + 'LATITUDE','HEIGHT (KM)',histlab,fieldlab, + iframe,rec80,'tigcmdif',dirascii) endif if (i.le.nlon) then if (slon(i).ne.zmflag) then write(6,"('pltlon: frame ',i4,' field ',a8,' lon=',f8.2, + ' ht scale=',f7.2,' to ',f7.2)") iframe, + labshort(ip),rlon,htscale(1),htscale(nhtscale) else write(6,"('pltlon: frame ',i4,' field ',a8, + ' (zonal means) ht scale=',f7.2,' to ',f7.2)") + iframe,labshort(ip),htscale(1),htscale(nhtscale) endif else write(6,"('pltlon: frame ',i4,' field ',a8,' slt=',f8.2, + ' ht scale=',f7.2,' to ',f7.2)") iframe, + labshort(ip),sslt(islt),htscale(1),htscale(nhtscale) endif endif c c End selected longitude loop 200 continue c c End field loop: ip=1,ntimefld 100 continue if (ihtsc.gt.0) then call hpdeallc(ppltht,ier,1) call hpdeallc(ppltht0,ier,1) endif c return end