c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/pltlatht.f c------------------------------------------------------------------ c subroutine pltlatht(utc) c c Contour tgcm longitude slices (including dynamo fields): c (latitude on x-axis, height on y-axis) c include 'tgcmparam.h' include 'input.h' include 'tigcm.h' include 'tgcmhdr.h' include 'tigcmfld.h' include 'tiegcmfld.h' include 'selgrid.h' include 'tgcmlab.h' include 'cool.h' include 'plt.h' include 'color.h' dimension plotp(jmx,kmx) character*56 toplab,lab,lab1 c c Declarations for axes labeling: c (numy is determined in the code) c parameter (nx=7,ny=5) dimension numx(nx),numy(ny) data numx/-90,-60,-30,0,30,60,90/ character*8 fmtx,fmty character*16 labx,laby data fmtx/'(i3) '/, fmty/'(i3) '/ data nfmtx/3/, nfmty/3/, mnrx/3/, mnry/4/ data laby/' HEIGHT (KM) '/ data labx/' LATITUDE (DEG) '/ data xc/-87.5/,xd/87.5/ c c Set up conpack: c call cpseti('SET',0) call cpseti('MAP - mapping flag',0) c c Viewport may be changed here, e.g., for narrow height scales, c (like sme comparisons), we may want a shorter viewport, c so change slmapt and slmapb) (see also comments in tigcm.h) c slmapl = 0.10 slmapr = 0.93 c c slmapb,t below for typical 100-500km tigcm slmapt = 0.91 slmapb = 0.26 c c slmapb,t below for sme comparison (e.g., htscale=100->160km): c slmapt = 0.90 c slmapb = 0.50 c c Numeric labels for y-axis c (watch this when nhtscale,dht,ht1,htscale are changed) c write(6,"('pltlatht: ht1=',f12.8,' dht=',f12.8,' nhtscale=',i3, c + ' htscale=',/(6f12.8))") ht1,dht,nhtscale, c + (htscale(i),i=1,nhtscale) ii = 0 do 25 i=1,nhtscale,ny ii = ii+1 numy(ii) = nint(htscale(i)) c write(6,"('do 25: i=',i3,' ii=',i3,' htscale(i)=',f12.8, c + ' numy(ii)=',i5)") i,ii,htscale(i),numy(ii) 25 continue if (ii.ne.ny) then write(6,"(' >>> pltlatht warning: numy=',/(10i5))") + (numy(i),i=1,ny) write(6,"(' setting numy(ny)=htscale(nhtscale)')") numy(ny) = htscale(nhtscale) endif c rilx = 0.5*(slmapl+slmapr) rily = -.30 call cpsetr('ILX - info label x-coord',rilx) call cpsetr('ILY - info label y-coord',rily) call set(slmapl,slmapr,slmapb,slmapt, + -90.,90.,htscale(1),htscale(nhtscale),1) yc = htscale(1) yd = htscale(nhtscale) c c Take year-day from header: c This is not reliable -- removed c iyrday = (date(1)-1900.)*1000.+date(2) c c Field loop: c write(6,"(' ')") write(6,"('Longitude slices (lat on x-axis, height on y-axis):')") do 100 ip=1,ntigcmf+ntiegcmf ipp = ip-ntigcmf if (ip.le.ntigcmf.and.iptigcm(ip).le.0) goto 100 if (ip.gt.ntigcmf.and.(idyn.le.0.or.iptiegcm(ipp).le.0)) + goto 100 c c Height independent fields are skipped: if (ip.eq.ixfof2.or.ip.eq.ixhmf2) then write(6,"(' pltlatht: skipping ',a8,' because ', + 'this field is height independent')") labtigcm_short(ip) goto 100 endif c c Selected longitude loop: c do 105 i=1,nlon+nslt c c Find index to selected longitude, or local time: if (i.le.nlon) then if (slon(i).eq.r12flag) then ilon = ixslt(12.,utc,slon(i)) else ilon = islon(i) endif else ilon = ixslt(sslt(i-nlon),utc,dum) endif c c If we are doing zonal means, flag tgcmint to do extrapolation when c desired height > max height at grid point; otherwise let tgcmint c return cpspval when desired height > max height: c spvint = cpspval if (ilon.eq.izmflag) spvint = 0. c c Set up top label: if (ndays.eq.1) then if (i.le.nlon) then if (ilon.ne.izmflag) then write(toplab,"('UT = ',f5.2,' LON = ',f7.2, + ' DAY = ',i5,20x)") utc,slon(i),iyd(1) else write(toplab,"('UT = ',f5.2,' ZONAL MEAN', + ' DAY = ',i5,23x)") utc,iyd(1) endif else write(toplab,"('UT = ',f5.2,' SLT = ',f7.2, + ' DAY = ',i5,20x)") utc,sslt(i-nlon),iyd(1) endif else if (i.le.nlon) then if (ilon.ne.izmflag) then write(toplab,"('UT = ',f5.2,' LON = ',f7.2,32x)") + utc,slon(i) else write(toplab,"('UT = ',f5.2,' ZONAL MEAN',35x)") utc endif else write(toplab,"('UT = ',f5.2,' SLT = ',f7.2,32x)") + utc,sslt(i-nlon) endif endif c c Define plot array: c if (ilon.eq.izmflag) plotp(:,:) = 0. if (ip.le.ntigcmf) then if (ip.le.ngcmflds) then call tgcmint(pnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,htscale,nhtscale,igcmlog(ip),2, + finterp,imx,jmx,mxhtscal,kmx,gcmzp,ier,1,spvint,0) if (ilon.ne.izmflag) then do 110 k=1,nhtscale plotp(:,k) = finterp(ilon,:,k) 110 continue else do k=1,nhtscale do j=1,jmx plotp(j,k) = ave(finterp(1,j,k),imx-1,cpspval) enddo enddo endif elseif (ip.eq.ixrho) then call tgcmint(td,pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,htscale,nhtscale,igcmlog(ip),2, + finterp,imx,jmx,mxhtscal,kmx,gcmzp,ier,1,spvint,0) if (ilon.ne.izmflag) then do 115 k=1,nhtscale plotp(:,k) = finterp(ilon,:,k) 115 continue else do k=1,nhtscale do j=1,jmx plotp(j,k) = ave(finterp(1,j,k),imx-1,cpspval) enddo enddo endif c c o/(o2+n2): elseif (ip.eq.ixrat) then c c Temp store o2 in udat: call tgcmint(pnt(1,1,1,ixo2),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,htscale,nhtscale,igcmlog(ip),2, + finterp,imx,jmx,mxhtscal,kmx,gcmzp,ier,1,spvint,0) do 120 k=1,nhtscale do 120 j=1,jmx udat(j,k) = finterp(ilon,j,k) 120 continue c c Temp store n2 in vdat: call tgcmint(pnt(1,1,1,ixn2),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,htscale,nhtscale,igcmlog(ip),2, + finterp,imx,jmx,mxhtscal,kmx,gcmzp,ier,1,spvint,0) do 121 k=1,nhtscale do 121 j=1,jmx vdat(j,k) = finterp(ilon,j,k) 121 continue c c Temp store o1 in plotp: call tgcmint(pnt(1,1,1,ixo1),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,htscale,nhtscale,igcmlog(ip),2, + finterp,imx,jmx,mxhtscal,kmx,gcmzp,ier,1,spvint,0) do 122 k=1,nhtscale plotp(:,k) = finterp(ilon,:,k) 122 continue c c Find ratio and store in plotp: if (ilon.ne.izmflag) then do 123 k=1,nhtscale do 123 j=1,jmx plotp(j,k) = plotp(j,k) / (udat(j,k)+vdat(j,k)) 123 continue else write(6,"('Zonal mean not available for ratio')") goto 105 endif endif else c c Dynamo fields: call tgcmint(dynpnt(1,1,1,ipp),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,htscale,nhtscale,idynlog(ipp),2, + finterp,imx,jmx,mxhtscal,kmx,gcmzp,ier,1,spvint,0) if (ilon.ne.izmflag) then do 124 k=1,nhtscale plotp(:,k) = finterp(ilon,:,k) 124 continue else do k=1,nhtscale do j=1,jmx plotp(j,k) = ave(finterp(1,j,k),imx-1,cpspval) enddo enddo endif endif c c Always take logs of appropriate fields, since ht on y-axis: c if ((ip.le.ntigcmf.and.igcmlog(ip).gt.0).or. + (ip.gt.ntigcmf.and.idynlog(ipp).gt.0)) then do 125 k=1,kmx do 125 j=1,jmx if (plotp(j,k).eq.cpspval.or.plotp(j,k).le.1.e-30) + then plotp(j,k) = cpspval else plotp(j,k) = alog10(plotp(j,k)) endif 125 continue if (ip.le.ntigcmf) + write(lab,"('LOG10 ',a50)") labtigcm(ip)(1:50) if (ip.gt.ntigcmf) + write(lab,"('LOG10 ',a50)") labtiegcm(ipp)(1:50) else if (ip.le.ntigcmf) lab = labtigcm(ip) if (ip.gt.ntigcmf) lab = labtiegcm(ipp) endif c c Contour it: if (ip.le.ntigcmf) then cint = conints(ip) cmin = conmins(ip) cmax = conmaxs(ip) ipfm = ip else cint = cintdyn(ipp) cmin = cmindyn(ipp) cmax = cmaxdyn(ipp) ipfm = ip+ntigcmf endif if (iclrfill.le.0) then call bwcon(plotp,jmx,jmx,nhtscale,xc,xd,yc,yd, + cint,cmin,cmax, + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) else lbar = 2 call clrcon(plotp,jmx,jmx,nhtscale,xc,xd,yc,yd, + cint,cmin,cmax, + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,slmapl,slmapr,.09,.18,ipfm) endif c call box(1) c c Top labels: call top2lab(lab,0.11,toplab,0.05,.017) c c Bottom runlab label and history label: if (iclrfill.le.0) then ylab = 0.15 if (htscale(nhtscale).lt.200.) ylab = .35 call botlab(runlab(1:lnblnk2(runlab)),0.5*(slmapl+slmapr), + ylab) ylab = 0.11 if (htscale(nhtscale).lt.200.) ylab = .31 write(lab1,"(7x,'TGCM History = ',4a8,2x)") + (output(ii,1),ii=1,3),output(3,2) call botlab(lab1,0.5,ylab) else ylab = 0.05 if (htscale(nhtscale).lt.200.) ylab = .25 call botlab(runlab(1:lnblnk2(runlab)),0.5*(slmapl+slmapr), + ylab) ylab = 0.02 if (htscale(nhtscale).lt.200.) ylab = .21 write(lab1,"(7x,'TGCM History = ',4a8,2x)") + (output(ii,1),ii=1,3),output(3,2) call botlab(lab1,0.5,ylab) endif c c Finish up the frame: c call box(1) call frame iframe=iframe+1 if (i.le.nlon) then if (ilon.ne.izmflag) then write(6,"(' pltlatht frame ',i3,': ',a,' ut=',f4.1, + ' lon=',f6.1)") iframe,lab(1:36),utc,slon(i) else write(6,"(' pltlatht frame ',i3,': ',a,' ut=',f4.1, + ' Zonal Mean')") iframe,lab(1:36),utc endif else write(6,"(' pltlatht frame ',i3,': ',a,' ut=',f4.1, + ' slt=',f6.1)") iframe,lab(1:36),utc,sslt(i-nlon) endif c c End selected lon loop: 105 continue c c End field loop: 100 continue c c Not doing cooling terms with height interpolation: c if (iplcool.gt.0) then write(6,"(' Cooling terms with height on y-axis not yet ', + 'available')") endif c return end