c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/dynglb.f c------------------------------------------------------------------ c subroutine dynglb(utc) c c Make global (cylindrical equidistant) contours 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,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/ data uvmxlat/50./ ! max lat to plot ui,vi vectors 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 ntiegcmf = 6: potential, ui, vi, wi, sumpot, hpot c write(6,"(' ')") write(6,"('Global contours (dynamo fields): npls=',i2, + ' nhts=',i2)") npls,nhts do 100 ip=1,ntiegcmf if (iptiegcm(ip).le.0) goto 100 c c Selected pressure and heights loop: c 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 6-9 are heelis potential, ui, vi, and wi: c if (ip.le.ntiegcmf-5) then plotp(:,:) = dynpnt(:,ispls(k),:,ip) elseif (ip.eq.ixspot) then ! sum potentials plotp(:,:) = dynpnt(:,ispls(k),:,ixpot) + hpoten(:,:) c c Height independent heelis potential and ion drifts: elseif (ip.gt.ixspot) then if (k.gt.1) goto 110 write(toplab,"('UT = ',f6.2,45x)") utc endif if (ip.eq.ixhpot) then plotp(:,:) = hpoten(:,:) elseif (ip.eq.ixhui) then plotp(:,:) = hui(:,:) elseif (ip.eq.ixhvi) then plotp(:,:) = hvi(:,:) elseif (ip.eq.ixhwi) then plotp(:,:) = hwi(:,:) endif c c Interpolate to desired height: else write(toplab,"('UT = ',f6.2,' HT = ',f6.1,33x)") + utc,shts(kk) c c Define plot array (interpolate to height surface): if (ip.le.ntiegcmf-5) then call tgcmint(dynpnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(ip),2,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) 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,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) plotp(:,:) = plotp(:,:) + hpoten(:,:) else ! height independent heelis potential write(6,"('pltglb skipping heelis because', + ' it is height independent: ip=',i3)") ip goto 110 endif endif c c Take logs if necessary: if (idynlog(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)") labtiegcm(ip)(1:50) else lab = labtiegcm(ip) if (ip.eq.ixpot.and.ipltuv.eq.-1) + write(lab,"('TIEGCM ION DRIFT (M/S)',34x)") endif c c Contour it, with axes and axis labels: c if (ip.eq.ixpot.and.ipltuv.eq.-1) goto 101 if (iclrfill.le.0) then call bwcon(plotp,imx,imx,jmx,xc,xd,yc,yd, + cintdyn(ip),cmindyn(ip),cmaxdyn(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, + cintdyn(ip),cmindyn(ip),cmaxdyn(ip), + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,pmapl,pmapr,.09,.19,ip+ntigcmf) endif 101 continue 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,"(8x,'TIGCM History = ',3a8,8x)") + (histvol(i,ivf),i=1,3) call botlab(lab1,0.5*(pmapl+pmapr),.20) endif c c Add tiegcm ion drift vectors if doing tiegcm potential: if ((ipltuv.eq.1.or.ipltuv.eq.-1).and. + (ip.eq.ixpot.or.ip.eq.ixhpot.or.ip.eq.ixspot)) then if (ip.eq.ixpot) then if (k.le.npls) then udat(:,:) = ui(:,ispls(k),:) vdat(:,:) = vi(:,ispls(k),:) else c c Do height interpolation of ui,vi: call tgcmint(dynpnt(1,1,1,ixui),pnt(1,1,1,ixz),gcmlat, + jmx,gcmlon,imx,shts(kk),1,igcmlog(ixui),2,udat,imx, + jmx,1,kmx,gcmzp,ier,1,cpspval,0) call tgcmint(dynpnt(1,1,1,ixvi),pnt(1,1,1,ixz),gcmlat, + jmx,gcmlon,imx,shts(kk),1,igcmlog(ixvi),2,vdat,imx, + jmx,1,kmx,gcmzp,ier,1,cpspval,0) 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 udat(:,:) = hui(:,:) vdat(:,:) = hvi(:,:) 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) endif endif rmin = 1.e36 rmax = -1.e36 do 95 j=1,jmx do 95 i=1,imx if ((udat(i,j).eq.0..and.vdat(i,j).eq.0.).or. + (abs(gcmlat(j)).gt.uvmxlat)) then udat(i,j) = cpspval vdat(i,j) = cpspval else rmin = amin1(rmin,sqrt(udat(i,j)**2+vdat(i,j)**2)) rmax = amax1(rmax,sqrt(udat(i,j)**2+vdat(i,j)**2)) endif 95 continue c uvmax = 0. uvmax = 300. m = 1 ispvec = 0 spvec = 1.e36 iskip = 0 call pltvec(udat,vdat,uvmax,jmx,iskip,gcmlat(1), + gcmlat(jmx),m,ispvec,spvec) endif 96 continue if (ip.eq.ixpot.and.ipltuv.eq.-1) then write(lab,"('MINIMUM=',f10.2,', MAXIMUM=',f10.2,18x)") + rmin,rmax call botlab(lab(1:37),0.5*(pmapl+pmapr),.15) call labaxes(nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) endif c c Finish up: call frame iframe=iframe+1 if (k.le.npls) then if (ip.ne.ixhpot) then write(6,"('dynglb frame ',i3,': ',a,' UT=',f4.1, + ' ZP=',f4.1)") iframe,lab(1:40),utc,spls(k) else write(6,"('dynglb frame ',i3,': ',a,' UT=',f4.1)") + iframe,lab(1:40),utc endif else write(6,"('dynglb frame ',i3,': ',a,' UT=',f4.1, + ' HT=',f6.1)") iframe,lab(1:40),utc,shts(kk) endif c c End pressure loop: 110 continue c c End field loop 100 continue c return end