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) c data perimlat/27.5,-27.5/ c data perimlat/42.5,-42.5/ 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 c pmapl = .05 c pmapr = .80 pmapl = .125 pmapr = .875 pmapb = .15 pmapt = .90 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 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 = nint((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): c (first 4 dynamo fields are in dynpnt, 5th is sum of potentials, c and 6th is height independent heelis potential: c if (ip.le.ntiegcmf-5) then do j=iskip+1,jmx plotp(:,j-iskip) = dynpnt(:,ispls(k),j,ip) enddo elseif (ip.eq.ixspot) then ! sum potentials do j=iskip+1,jmx plotp(:,j-iskip) = dynpnt(:,ispls(k),j,ixpot) + + hpoten(:,j) enddo elseif (ip.gt.ixspot) then if (k.gt.1) goto 110 write(toplab,"('UT = ',f6.2,45x)") utc endif if (ip.eq.ixhpot) then ! height independent heelis potential: do j=iskip+1,jmx plotp(:,j-iskip) = hpoten(:,j) enddo elseif (ip.eq.ixhui) then do j=iskip+1,jmx plotp(:,j-iskip) = hui(:,j) enddo elseif (ip.eq.ixhvi) then do j=iskip+1,jmx plotp(:,j-iskip) = hvi(:,j) enddo elseif (ip.eq.ixhwi) then do j=iskip+1,jmx plotp(:,j-iskip) = hwi(:,j) enddo endif else c c Interpolate to desired height: write(toplab,"('UT = ',f6.2,' HT = ',f6.1,33x)") + utc,shts(kk) c c Define plot array (selected height): c SUBROUTINE TGCMINT(GCMIN,GCMHT,RLAT,NLAT,RLON,NLON,HTS,NHTS, c + LOGHT,IXFLAG,GCMOUT,IDIM1,IDIM2,IDIM3, c + NPRES,PRES,IER,IGRID,SPVAL,IPRNT) c if (ip.le.ntiegcmf-2) then 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 j=iskip+1,jmx plotp(:,j-iskip) = udat(:,j) enddo elseif (ip.eq.ixspot) then ! sum potentials call tgcmint(dynpnt(1,1,1,ixpot),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 j=iskip+1,jmx plotp(:,j-iskip) = udat(:,j) + hpoten(:,j) enddo else ! height independent heelis potential write(6,"('pltpol skipping heelis potential because', + ' it is height independent')") goto 110 endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Temp mods: hardwire some stuff for Ray: c (also see commented out labels below) c c toplab = labtiegcm(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 (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,"(a)") runlab call botlab(lab(1:lenstr(lab)),0.45,.01) c c Bottom history label: if (iclrfill.le.0) then write(lab,"('TIGCM History = ',3a8,16x)") + (histvol(i,ivf),i=1,3) call botlab(lab(1:lenstr(lab)),0.45,.035) endif 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 Overlay ion drift velocity vectors: if (ipltuv.eq.1.and. + (ip.eq.ixpot.or.ip.eq.ixhpot.or.ip.eq.ixspot)) then c c Ion drift ui,vi from dynamo potential: if (ip.eq.ixpot) 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 c c Add heelis ion drift vectors if doing heelis potential c (height independent): elseif (ip.eq.ixhpot) then if (k.gt.1) goto 96 do j=iskip+1,jmx udat(:,j-iskip) = hui(:,j) vdat(:,j-iskip) = hvi(:,j) enddo c c Add ion drifts calculated from sum of dynamo and heelis potentials: c (Cannot yet do this at interpolated heights) c elseif (ip.eq.ixspot) then if (k.gt.npls) then write(6,"('dynglb: ion drifts for sum of ', + 'potentials at interpolated heights not yet', + ' available')") goto 96 else call mkuivis(k,udat,vdat,pnt(1,1,1,ixz),lumag) plotp(:,:) = udat(:,:) do j=iskip+1,jmx udat(:,j-iskip) = plotp(:,j) enddo plotp(:,:) = vdat(:,:) do j=iskip+1,jmx vdat(:,j-iskip) = plotp(:,j) enddo endif endif c c Check vectors for zeros and plot: 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 96 continue c c Complete the frame: c c call box(1) 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