c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/dynpol.f c------------------------------------------------------------------ c subroutine dynpol(utc) c c Make polar (stereographic) contours of difference fields: c include 'tgcmparam.h' include 'input.h' include 'tigcm.h' include 'tgcmhdr.h' include 'tiegcmfld.h' include 'tigcmfld.h' include 'selgrid.h' include 'tgcmlab.h' include 'cool.h' include 'plt.h' include 'color.h' dimension plotp(imx,jmx) character*56 toplab,lab 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 dimension perimlat(nperim) data perimlat/27.5,-27.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 = .05 pmapr = .80 pmapb = .13 pmapt = .91 rilx = 0.5 rily = pmapb-0.25 call cpsetr('ILX - info label x-coord',rilx) call cpsetr('ILY - info label y-coord',rily) c c Hemisphere loop: c (ih=1 -> north, ih=2 -> south) c (iplpol=1 -> south only, iplpol=2 -> north only, iplpol=3 -> both) c if (iplpol.eq.1) then ! south only ifirst = 2 nhem = 2 elseif (iplpol.eq.2) then ! north only ifirst = 1 nhem = 1 elseif (iplpol.eq.3) then ! both hems ifirst = 1 nhem = 2 endif do 500 ih=ifirst,nhem hem = 'S' if (perimlat(ih).ge.0.) hem = 'N' nlats = (90.-abs(perimlat(ih))+dlat/2.)/dlat 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 = perimlat(ih) if (perimlat(ih).gt.0.) ycn = -perimlat(ih) elseif (hem.eq.'N') then iskip = jmx-nlats rot = -utc*15. plat = 90. yc1 = perimlat(ih) ycn = gcmlat(jmx) endif 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(perimlat(ih)) 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 (dynamo):')") if (ih.eq.2) + write(6,"('Polar (stereographic south) contours (dynamo):')") do 100 ip=1,ntiegcmf if (iptiegcm(ip).le.0) goto 100 c c Selected pressure loop: do 110 k=1,npls+nhts kk = k-npls if (k.le.npls) then write(toplab,"('UT = ',f6.2,' ZP = ',f4.1,35x)") utc,spls(k) c c Define plot array (selected pressure): do 120 j=iskip+1,jmx plotp(:,j-iskip) = dynpnt(:,ispls(k),j,ip) 120 continue else write(toplab,"('UT = ',f6.2,' HT = ',f6.1,33x)") + utc,shts(kk) c c Define plot array (selected height): call tgcmint(dynpnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(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 c c Take logs if necessary: if (idynlog(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)") labtiegcm(ip)(1:50) else lab = labtiegcm(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, + cintdyn(ip),cmindyn(ip),cmaxdyn(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, + cintdyn(ip),cmindyn(ip),cmaxdyn(ip), + 0,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,.83,.98,pmapb,pmapt,ip+ntigcmf) 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(labtiegcm(ip),0.11,toplab,0.07,.017) write(lab,"(12x,a,12x)") runlab call botlab(lab,0.45,.01) c c Add perimeter latitude label in upper right corner: write(latlab,"('Perim lat = ',f5.1,7x)") + perimlat(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 ion drift vectors if doing potential: c if (ip.eq.ixpot.and.ipltuv.eq.1) then if (k.le.npls) then do 150 j=iskip+1,jmx udat(:,j-iskip) = ui(:,ispls(k),j) vdat(:,j-iskip) = vi(:,ispls(k),j) 150 continue else c c Do height interpolation of ui,vi: c call tgcmint(dynpnt(1,1,1,ixui),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(ixui),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(dynpnt(1,1,1,ixvi),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(ixvi),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.le.npls) then write(6,"('dynpol frame ',i3,': ',a,' UT=',f4.1, + ' ZP=',f4.1)") iframe,labtiegcm(ip)(1:40),utc,spls(k) else write(6,"('dynpol frame ',i3,': ',a,' UT=',f4.1, + ' HT=',f6.1)") iframe,labtiegcm(ip)(1:40),utc,shts(kk) endif c c End pressure/height loop: 110 continue c c End field loop 100 continue c c End hem loop 500 continue c return end