c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/pltpol.f c------------------------------------------------------------------ c subroutine pltpol(utc) c c Make polar (stereographic) contours of difference fields: 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 character*24 latlab c c Below are essentially dummies for bwcon and clrcon calls: 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)'/ c c Temp: parameter(nperim=2) character*1 hem c dimension perimlat(nperim) c data perimlat/27.5,-27.5/ c data perimlat/47.5,-47.5/ data maps/1/, iwhite/1/, iccont/1/ c c ispvec = special value flag for velvct c spvec = specil value for velvct data ispvec/0/, spvec/1.e36/ c c Tell conpack we will use cpmpxy for transformations: c call cpseti('MAP - mapping flag',1) call cpseti('SET - do set call flag',0) c c Using stereographic (elliptic) projection c pmapl,r,b,t = left,right,bottom,top map position c rilx,y are coords of conpack info label c proj = 'ST ' iel = 1 pmapl = .04 pmapr = .79 pmapb = .13 pmapt = .91 if (iclrfill.le.0) then pmapl = .125 pmapr = .875 pmapb = .15 pmapt = .90 endif rilx = 0.5 rily = pmapb-0.28 call cpsetr('ILX - info label x-coord',rilx) call cpsetr('ILY - info label y-coord',rily) c c Loop through selected perimeter latitudes (polperim): c do 500 ih=1,iplpol hem = 'S' if (polperim(ih).ge.0.) hem = 'N' nlats = nint((90.-abs(polperim(ih))+dlat/2.)/dlat) write(6,"('pltpol: polperim=',f10.3,' nlats=',i3)") + polperim(ih),nlats c c Set up rotation, etc: c xc1 = gcmlon(1) xcm = gcmlon(imx) plon = 0. if (hem.eq.'S') then iskip = 0 rot = utc*15.-180. plat = -90. yc1 = gcmlat(1) ycn = polperim(ih) if (polperim(ih).gt.0.) ycn = -polperim(ih) elseif (hem.eq.'N') then iskip = jmx-nlats rot = -utc*15. plat = 90. yc1 = polperim(ih) ycn = gcmlat(jmx) endif write(6,"('pltpol: ih=',i3,' polperim=',f9.2,' nlats=',i3, + ' iskip=',i3)") ih,polperim(ih),nlats,iskip C C Set up EZMAP: C call mappos(pmapl,pmapr,pmapb,pmapt) call mapsti('EL',iel) call mapstc('OU','CO') call mapsti('DO',1) call maproj(proj,plat,plon,rot) plat = 90.-abs(polperim(ih)) - .01 call mapset('AN',plat,plat,plat,plat) call mapsti('G1',1) call mapsti('G2',2) call mapsti('VS',0) call mapint call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) c c Field loop: c write(6,"(' ')") if (ih.eq.1) + write(6,"('Polar (stereographic north) contours:')") if (ih.eq.2) + write(6,"('Polar (stereographic south) contours:')") do 100 ip=1,ntigcmf if (iptigcm(ip).le.0) goto 100 c c Selected pressure 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) then do 120 j=iskip+1,jmx plotp(:,j-iskip) = fof2(:,j) 120 continue elseif (ip.eq.ixhmf2) then do 125 j=iskip+1,jmx plotp(:,j-iskip) = hmf2(:,j) 125 continue endif else goto 110 endif else c c Height dependent fields: 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 c c At constant pressure surface: if (k.le.npls) then do 130 j=iskip+1,jmx plotp(:,j-iskip) = pnt(:,ispls(k),j,ip) 130 continue 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,udat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) do 132 j=iskip+1,jmx plotp(:,j-iskip) = udat(:,j) 132 continue endif elseif (ip.eq.ixrho) then c c Total density: if (k.le.npls) then do 135 j=iskip+1,jmx plotp(:,j-iskip) = td(:,ispls(k),j) 135 continue else call tgcmint(td,pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ip),2,udat,imx,jmx, + 1,jmx,gcmzp,ier,1,cpspval,0) do 137 j=iskip+1,jmx plotp(:,j-iskip) = udat(:,j) 137 continue endif c c Atomic/Mol ratio: elseif (ip.eq.ixrat) then if (k.le.npls) then do 140 j=iskip+1,jmx plotp(:,j-iskip) = pnt(:,ispls(k),j,ixo1) / + (pnt(:,ispls(k),j,ixo2) + pnt(:,ispls(k),j,ixn2)) 140 continue 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,jmx,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,jmx,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,jmx,gcmzp,ier,1,cpspval,0) do 142 j=iskip+1,jmx plotp(:,j-iskip) = plotp(:,j) / (udat(:,j) + vdat(:,j)) 142 continue endif endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Temp mods: hardwire some stuff for Ray: c (also see commented out labels below) c c if (ip.eq.ixo2) then c plotp(:,:) = plotp(:,:) / 1.e+10 c write(toplab, c + "('TIGCM O2 NUMBER DENSITY x 1.E10 (CM-3)',18x)") c elseif (ip.eq.ixo1) then c plotp(:,:) = plotp(:,:) / 1.e+11 c write(toplab, c + "('TIGCM O NUMBER DENSITY x 1.E11 (CM-3)',19x)") c elseif (ip.eq.ixno) then c plotp(:,:) = plotp(:,:) / 1.e+07 c write(toplab, c + "('TIGCM NO NUMBER DENSITY x 1.E7 (CM-3)',19x)") c else c toplab = labtigcm(ip) c endif c toplab = labtigcm(ip) c xmid = 0.5*(pmapl+pmapr) c ypos = pmapt+.08 c call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) c call set(0.,1.,0.,1.,0.,1.,0.,1.,1) c call plchhq(xmid,ypos,toplab(1:lnblnk2(toplab)), c + .015,0.,0.) c call set(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) c c End temp mods: cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Take logs if necessary: if (igcmlog(ip).gt.0.and.log10pl.gt.0) then do 145 j=iskip+1,jmx do 145 i=1,imx if (plotp(i,j-iskip).eq.cpspval.or. + plotp(i,j-iskip).le.1.e-30) then plotp(i,j-iskip) = cpspval else plotp(i,j-iskip) = alog10(plotp(i,j-iskip)) endif 145 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,nlats,xc1,xcm,yc1,ycn, + conints(ip),conmins(ip),conmaxs(ip), + 0,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) else lbar = 1 call clrcon(plotp,imx,imx,nlats,xc1,xcm,yc1,ycn, + conints(ip),conmins(ip),conmaxs(ip), + 0,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,.83,.98,pmapb,pmapt,ip) endif 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(labtigcm(ip),0.11,toplab,0.07,.017) write(lab,"(a)") runlab call botlab(runlab(1:lenstr(runlab)),0.45,.01) c c Bottom history label: if (iclrfill.le.0) then write(lab1,"(8x,'TIGCM History = ',3a8,8x)") + (histvol(i,ivf),i=1,3) call botlab(lab1,0.5*(pmapl+pmapr),.04) endif c c Add perimeter latitude label in upper right corner: write(latlab,"('Perim lat = ',f5.1,7x)") + polperim(ih) szfs = .017*(vr-vl) call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) ypos = vt + 0.11*(vt-vb) call plchhq(0.96,ypos,latlab(1:17),szfs,0.,1.) call set(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) c c Add local time labels: c call labpolar(pmapl,pmapr,pmapb,pmapt,hem) c c Add neutral wind vectors if doing temp: c if (ip.eq.ixt.and.ipltuv.eq.1) then if (k.le.npls) then do 150 j=iskip+1,jmx udat(:,j-iskip) = pnt(:,ispls(k),j,ixu) vdat(:,j-iskip) = pnt(:,ispls(k),j,ixv) 150 continue 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,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) do 151 j=iskip+1,jmx udat(:,j-iskip) = plotp(:,j) 151 continue call tgcmint(pnt(1,1,1,ixv),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,igcmlog(ixv),2,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) do 152 j=iskip+1,jmx vdat(:,j-iskip) = plotp(:,j) 152 continue 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 = 2 if (hem.eq.'N') m = 3 call pltvec(udat,vdat,uvmax,nlats,iskip,yc1,ycn,m, + ispvec,spvec) endif c c Complete the frame: c call frame iframe=iframe+1 if (k.eq.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) then write(6,"('pltpol frame ',i3,': ',a,' UT=',f4.1)") + iframe,labtigcm(ip)(1:40),utc else if (k.le.npls) then write(6,"('pltpol frame ',i3,': ',a,' UT=',f4.1, + ' ZP=',f4.1)") + iframe,labtigcm(ip)(1:40),utc,spls(k) else write(6,"('pltpol frame ',i3,': ',a,' UT=',f4.1, + ' HT=',f6.1)") + iframe,labtigcm(ip)(1:40),utc,shts(kk) endif endif c c End selected pressure/height loop 110 continue c c End field loop 100 continue c c Height integrate cooling terms (use udat as temp storage): 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=iskip+1,jmx do 210 i=1,imx plotp(i,j-iskip) = 0. tncol(:) = pnt(i,:,j,ixt) htcol(:) = pnt(i,:,j,ixz) c c cool integrate k loop: 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-iskip) = plotp(i,j-iskip) + + (coolk*rz + coolkm1*rzm) * rw c c End cool integrate k loop: 205 continue c c End cool i,j loop: 210 continue c c Contour integrated cooling term: call bwcon(plotp,imx,imx,nlats,xc1,xcm,yc1,ycn,0.,1.,0., + 0,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) 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.11,toplab,0.07,.017) write(lab,"(12x,a,12x)") runlab call botlab(lab,0.45,.01) write(latlab,"('Perim lat = ',f5.1,7x)") + polperim(ih) szfs = .017*(vr-vl) call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) ypos = vt + 0.11*(vt-vb) call plchhq(0.96,ypos,latlab(1:17),szfs,0.,1.) call set(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) call labpolar(pmapl,pmapr,pmapb,pmapt,hem) call frame iframe=iframe+1 write(6,"('pltpol frame ',i3,': ',a, + ' (height integrated)')") iframe,labcooli(ip)(1:40) c c End field loop for cooling terms: 200 continue c c End cooling terms conditional: endif c c End hem loop 500 continue c return end