c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/pltglb.f c------------------------------------------------------------------ c subroutine pltglb(utc) c c Make global (cylindrical equidistant) contours c include 'tgcmparam.h' include 'input.h' include 'tigcm.h' include 'tgcmhdr.h' include 'tigcmfld.h' include 'selgrid.h' include 'tgcmlab.h' include 'cool.h' include 'plt.h' include 'color.h' dimension plotp(imx,jmx),tncol(kmx),htcol(kmx) character*56 toplab,lab,lab1 c c Declarations for axes labeling: c parameter (nx=7,ny=7) dimension numx(nx),numy(ny) data numx/-180,-120,-60,0,60,120,180/ data numy/-90,-60,-30,0,30,60,90/ character*8 fmtx,fmty character*16 labx,laby data fmtx/'(i4) '/, fmty/'(i3) '/ data nfmtx/4/, nfmty/3/, mnrx/3/, mnry/3/ data laby/' LATITUDE (DEG) '/ data labx/' LONGITUDE (DEG)'/ data xc/-180./,xd/180./ data yc/-87.5/,yd/87.5/ c c ispvec = special value flag for velvct c spvec = specil value for velvct data ispvec/0/, spvec/1.e36/ c c Temp: data maps/1/, iccont/1/ c c Set up conpack: c call cpseti('SET',0) pmapl = 0.10 pmapr = 0.95 pmapb = 0.455 pmapt = 0.880 wl = -180. wr = +180. wb = -90. wt = +90. rilx = 0.5 rily = 0.08 rily = -((pmapb-rily) / (pmapt-pmapb)) call cpsetr('ILX - info label x-coord',rilx) call cpsetr('ILY - info label y-coord',rily) iel = 0 c c Set up ezmap, if doing continental outlines: c if (maps.le.0) then call set(pmapl,pmapr,pmapb,pmapt,wl,wr,wb,wt,1) else proj = 'CE ' plat = 0. plon = 0. call cpseti('MAP - mapping flag',1) rot = 0. call mappos(pmapl,pmapr,pmapb,pmapt) call mapsti('EL',iel) call mapstc('OU','CO') call mapsti('DO',1) call maproj(proj,plat,plon,rot) call mapset('MA',0.,0.,0.,0.) call mapsti('G1',1) call mapsti('G2',2) call mapsti('VS',0) call mapint endif c c Field loop: c write(6,"(' ')") write(6,"('Global (cylindrical equidistant) contours:')") do 100 ip=1,ntigcmf if (iptigcm(ip).le.0) goto 100 c c Selected pressure and heights loop: c do 110 k=1,npls+nhts kk = k-npls c c Height independent fields: if (ip.eq.ixfof2.or.ip.eq.ixhmf2) then if (k.eq.1) then write(toplab,"('UT = ',f6.2,45x)") utc if (ip.eq.ixfof2) plotp(:,:) = fof2(:,:) if (ip.eq.ixhmf2) plotp(:,:) = hmf2(:,:) else goto 110 endif c c Height dependent fields: else if (k.le.npls) then write(toplab,"('UT = ',f6.2,' ZP = ',f4.1,35x)") + utc,spls(k) else write(toplab,"('UT = ',f6.2,' HT = ',f6.1,33x)") + utc,shts(kk) endif if (ip.le.ngcmflds) then ! use pnt c c At constant pressure surface: if (k.le.npls) then plotp(:,:) = pnt(:,ispls(k),:,ip) c c Interpolate to constant height surface: else call tgcmint(pnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ip),2,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) endif c c Total density: elseif (ip.eq.ixrho) then if (k.le.npls) then plotp(:,:) = td(:,ispls(k),:) else call tgcmint(td,pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ip),2,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) endif c c Atomic/Mol ratio: elseif (ip.eq.ixrat) then if (k.le.npls) then plotp(:,:) = pnt(:,ispls(k),:,ixo1) / + (pnt(:,ispls(k),:,ixo2) + pnt(:,ispls(k),:,ixn2)) c c Height interp for ratio (use udat,vdat for temp storage of c interpolated o2 and n2, and use plotp for o1, then put c ratio in plotp: c else call tgcmint(pnt(1,1,1,ixo1),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ip),2,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) call tgcmint(pnt(1,1,1,ixo2),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ip),2,udat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) call tgcmint(pnt(1,1,1,ixn2),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ip),2,vdat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) plotp(:,:) = plotp(:,:) / (udat(:,:) + vdat(:,:)) endif endif endif c c Take logs if necessary: if (igcmlog(ip).gt.0.and.log10pl.gt.0) then do 125 j=1,jmx do 125 i=1,imx if (plotp(i,j).eq.cpspval.or.plotp(i,j).le.1.e-30) + then plotp(i,j) = cpspval else plotp(i,j) = alog10(plotp(i,j)) endif 125 continue write(lab,"('LOG10 ',a50)") labtigcm(ip)(1:50) else lab = labtigcm(ip) endif c c Contour it, with axes and axis labels: c if (iclrfill.le.0) then call bwcon(plotp,imx,imx,jmx,xc,xd,yc,yd, + conints(ip),conmins(ip),conmaxs(ip), + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) else lbar = 2 call clrcon(plotp,imx,imx,jmx,xc,xd,yc,yd, + conints(ip),conmins(ip),conmaxs(ip), + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,pmapl,pmapr,.09,.19,ip) endif c call box(0) c c Add local time axis at bottom: call sltxax(utc) c c Draw continental outlines, if desired: c (white for b/w plots, iccont for color fill) c if (maps.gt.0) then call gqplci(ierr,isplci) call gqpmci(ierr,ispmci) if (iclrfill.le.0) then call gsplci(iwhite) call gspmci(iwhite) else call gsplci(iccont) call gspmci(iccont) endif call maplot call gsplci(isplci) call gspmci(ispmci) endif c c Add labels: c call top2lab(lab,0.15,toplab,0.07,.017) call botlab(runlab(1:lnblnk2(runlab)),0.5*(pmapl+pmapr),0.25) c c Bottom history label: if (iclrfill.le.0) then write(lab1,"(7x,'TGCM History = ',4a8,2x)") + (output(ii,1),ii=1,3),output(3,2) call botlab(lab1,0.5*(pmapl+pmapr),.20) endif c c Add neutral wind vectors if doing temp: c if (ip.eq.ixt.and.ipltuv.eq.1) then if (k.le.npls) then udat(:,:) = pnt(:,ispls(k),:,ixu) vdat(:,:) = pnt(:,ispls(k),:,ixv) else c c Do height interpolation of un,vn: c call tgcmint(pnt(1,1,1,ixu),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ixu),2,udat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) call tgcmint(pnt(1,1,1,ixv),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ixv),2,vdat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) endif do 95 j=1,jmx do 95 i=1,imx if (udat(i,j).eq.0..and.vdat(i,j).eq.0.) then udat(i,j) = cpspval vdat(i,j) = cpspval endif 95 continue uvmax = 0. m = 1 ispvec = 0 spvec = 1.e36 iskip = 0 call pltvec(udat,vdat,uvmax,jmx,iskip,gcmlat(1), + gcmlat(jmx),m,ispvec,spvec) endif c c Finish up: call frame iframe=iframe+1 if (k.eq.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) then write(6,"('pltglb frame ',i3,': ',a,' UT=',f4.1)") + iframe,lab(1:40),utc else if (k.le.npls) then write(6,"('pltglb frame ',i3,': ',a,' UT=',f4.1, + ' ZP=',f4.1)") iframe,lab(1:40),utc,spls(k) else write(6,"('pltglb frame ',i3,': ',a,' UT=',f4.1, + ' HT=',f6.1)") iframe,lab(1:40),utc,shts(kk) endif endif c c End selected pressure loop: 110 continue c c End field loop (ip=1,ngcmflds) 100 continue c c Height integrated cooling/heating terms: c (there are ncoolt terms: c no heat, no cool, no heat-cool, co2 cool, o3p cool, irt) c if (iplcool.gt.0) then do 200 ip=1,ncoolt write(toplab,"('UT = ',f6.2,' (Height integration) ',23x)") + utc do 210 j=1,jmx do 210 i=1,imx plotp(i,j) = 0. tncol(:) = pnt(i,:,j,ixt) htcol(:) = pnt(i,:,j,ixz) do 205 k=2,kmx c c at k: tn = pnt(i,k,j,ixt) xo2 = pnt(i,k,j,ixo2) xo = pnt(i,k,j,ixo1) xn4s = pnt(i,k,j,ixn4s) xno = pnt(i,k,j,ixno) xn2d = pnt(i,k,j,ixn2d) xn2 = pnt(i,k,j,ixn2) xne = pnt(i,k,j,ixne) c c at k-1: tnm1 = pnt(i,k-1,j,ixt) xo2m1= pnt(i,k-1,j,ixo2) xom1 = pnt(i,k-1,j,ixo1) xn4sm1 = pnt(i,k-1,j,ixn4s) xnom1 = pnt(i,k-1,j,ixno) xn2dm1 = pnt(i,k-1,j,ixn2d) xn2m1 = pnt(i,k-1,j,ixn2) xnem1 = pnt(i,k-1,j,ixne) c c do integration: rz=1.66e-24*(16.*xo + 32.*xo2 + 28.*xn2) rzm=1.66e-24*(16.*xom1 + 32.*xo2m1 + 28.*xn2m1) coolk = cool(tn,xo2,xo,xn4s,xno,xn2d,xn2,xne,k, + tncol,htcol,ip) coolkm1 = cool(tnm1,xo2m1,xom1,xn4sm1,xnom1,xn2dm1, + xn2m1,xnem1,k-1,tncol,htcol,ip) rw = 0.5*(pnt(i,k,j,ixz) - pnt(i,k-1,j,ixz))*1.e+5 plotp(i,j) = plotp(i,j) + (coolk*rz + coolkm1*rzm) * rw c c End cool integrate k loop: 205 continue c c End cool i,j loop: 210 continue call bwcon(plotp,imx,imx,jmx,xc,xd,yc,yd,0.,1.,0., + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) call sltxax(utc) if (maps.gt.0) then call gqplci(ierr,isplci) call gqpmci(ierr,ispmci) if (iclrfill.le.0) then call gsplci(iwhite) call gspmci(iwhite) else call gsplci(iccont) call gspmci(iccont) endif call maplot call gsplci(isplci) call gspmci(ispmci) endif call top2lab(labcooli(ip),0.15,toplab,0.07,.017) call botlab(runlab(1:lnblnk2(runlab)),0.5*(pmapl+pmapr),0.11) call frame iframe=iframe+1 write(6,"('pltglb frame ',i3,': ',a,' (ht integrated)')") + iframe,labcooli(ip)(1:43) c c End field loop for cooling terms: 200 continue c c End cooling terms conditional: endif c return end