c c------------------------------------------------------------------ c Begin file /home/sting/foster/timegcm/pltampha.f c------------------------------------------------------------------ c subroutine pltampha(fields,utc,ihtsc) c c Contour longitude slices (lat on x-axis, zp on y-axis if ihtsc > 0, c linear height scale if ihtsc <= 0) c include 'timesproc.h' parameter(jmxp1=jmx+1, nfreq=4) character*56 toplab,fieldlab character*80 rec80 dimension plt(jmx,kmx),viewport(4),yaxht(kmx),rimx(imx), + fields(imx,kmx,jmx,nfget),zmhts(jmx,kmx),plt0(jmx,kmx) pointer (ppltht,pltht(1)), (pfht,fht(1)), (pamp,amp(1)), + (pphase,phase(1)) c data viewport /.15,.89,.26,.91/ data viewport /.13,.89,.17,.91/ dimension flo(4),hi(4) data flo/0.,6.,8.,9./, hi/24.,18.,16.,15./ c c /conre1/ is in conrec.a (old modified conrec) c /trans/ is in libplt.a, transormation functions fx,fy in module velvct c (fx and fy are used only by velvct (via pltvec) and this routine; c conpack does not use fx and fy) c fx and fy (in libplt.a) were modified to allow iqzqz=2, so that they c are mapped to xmin,xmax,ymin,ymax and nset=1 is passed to conrec so c my set call is used. This way, I can do set call from -90 to +90 on c x-axis, but data will be plotted at -87.5 to +87.5. c common /conre1/ ioffp,spval common/trans/mm,nn,xmin,xmax,ymin,ymax,if,iqzqz c iqzqz = 2 mm = jmx xmin = gcmlat(1) xmax = gcmlat(jmx) spval = cpspval ioffp = 1 xmid = 0.5*(viewport(1)+viewport(2)) nzp = izprange(2)-izprange(1)+1 if (ihtsc.le.0) then call set(viewport(1),viewport(2),viewport(3),viewport(4), + -90.,90.,gcmzp(izprange(1)),gcmzp(izprange(2)),1) call hpalloc(pamp,jmx*kmx*nfreq,ier,1) call hpalloc(pphase,jmx*kmx*nfreq,ier,1) nn = nzp ymin = gcmzp(izprange(1)) ymax = gcmzp(izprange(2)) else call set(viewport(1),viewport(2),viewport(3),viewport(4), + -90.,90.,htscale(1),htscale(nhtscale),1) call hpalloc(pamp,jmx*nhtscale*nfreq,ier,1) call hpalloc(pphase,jmx*nhtscale*nfreq,ier,1) call hpalloc(ppltht,jmx*nhtscale,ier,1) call hpalloc(pfht,imx*nhtscale*jmx,ier,1) nn = nhtscale ymin = htscale(1) ymax = htscale(nhtscale) endif c c Height for extra right axis (use globally averaged heights): c if (ihtsc.le.0) then yaxht(:) = 0. do k=izprange(1),izprange(2) izp = k-izprange(1)+1 do j=1,jmx rimx(:) = fields(:,k,j,ifget(itxz)) yaxht(izp) = yaxht(izp) + calcmean(rimx,imx,0,1.e-20, + cpspval) enddo yaxht(izp) = yaxht(izp) / float(jmx) enddo rnd = 10. if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5. else c c Find zonal mean heights for use in interpolation: c do k=1,kmx do j=1,jmx rimx(:) = fields(:,k,j,ifget(itxz)) zmhts(j,k) = calcmean(rimx,imx,0,1.e-20,cpspval) enddo enddo endif c c Field loop: c do 200 ip=1,nfhist if (ifplt(ip).le.0) goto 200 if (ihtsc.le.0) then call mkampha(fields(1,1,1,ifget(ip)),amp,phase,imx,kmx,jmx, + utc,cpspval) if (ip.ne.itxn2) call fixphase(phase,jmx,kmx,nfreq,0, + cpspval) if (ip.eq.itxn2) call fixphase(phase,jmx,kmx,nfreq,1, + cpspval) else call glbhtint(fields(1,1,1,ifget(ip)), + fields(1,1,1,ifget(itxz)),imx,kmx,jmx,fht,htscale,nhtscale, + itimelog(ip),cpspval,ier,0) call mkampha0(fht,amp,phase,imx,nhtscale,jmx,utc,cpspval) if (ip.ne.itxn2) call fixphase(phase,jmx,nhtscale,nfreq,0, + cpspval) if (ip.eq.itxn2) call fixphase(phase,jmx,nhtscale,nfreq,1, + cpspval) endif c c Four frequencies: c do 100 nf = 1,nfreq if (nf.gt.iampha) goto 100 c c Plot amplitudes: c if (ihtsc.le.0) then call defamp(amp,jmx,kmx,nfreq,plt0,nf,cpspval) do k=izprange(1),izprange(2) kk = k-izprange(1)+1 plt(:,kk) = plt0(:,k) enddo call conrec(plt,jmx,jmx,nzp,0.,0.,0.,1,-1,-1430B,1) call labrect(gcmlat,jmx,gcmzp(izprange(1)),nzp,'LATITUDE', + 'ZP',0.) else call defamp(amp,jmx,nhtscale,nfreq,pltht,nf,cpspval) call conrec(pltht,jmx,jmx,nhtscale,0.,0.,0.,1,-1,-1430B,1) call labrect(gcmlat,jmx,htscale,nhtscale,'LATITUDE', + 'HEIGHT (KM)',0.) endif call clearstr(toplab) write(toplab,"(a8,' AMPLITUDE (',i1,') UT=',f6.2)") + flab(ip),nf,utc if (iwrascii.gt.0) then call clearstr(rec80) write(rec80,"(a)") toplab(1:lenstr(toplab)) call clearstr(fieldlab) write(fieldlab,"(a8)") flab(ip)(1:lenstr(flab(ip))) endif call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.05,0.) write(toplab,"('HISTORY=',a)") histvol(ivol) call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.11, + .012) if(ihtsc.le.0) call altyax(nzp,yaxht,gcmzp(izprange(1)), + rnd,6) call frame iframe = iframe+1 if (ihtsc.le.0) then if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,plt,jmx,nzp,gcmlat, + gcmzp(izprange(1)),'LATITUDE','LN(P0/P)', + histvol(ivol),fieldlab,iframe,rec80,'timesproc', + dirascii) endif write(6,"('pltampha: frame ',i4,' field ',a8,' amp(',i1, + ') zprange=',2f8.2)") iframe,flab(ip),nf,zprange else if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,pltht,jmx,nhtscale,gcmlat, + htscale,'LATITUDE','HEIGHT (KM)',histvol(ivol),fieldlab, + iframe,rec80,'timesproc',dirascii) endif write(6,"('pltampha: frame ',i4,' field ',a8,' amp(',i1, + ') htscale=',f6.2,' to ',f6.2,' by ',f5.1)") iframe, + flab(ip),nf,htscale(1),htscale(nhtscale),delht endif c c Plot phases: c fflo = flo(nf) hhi = hi(nf) fi = (hhi-fflo)/12. scale = (hhi-fflo)/(2.*pi) if (ihtsc.le.0) then call defphase(phase,jmx,kmx,nfreq,plt0,scale,nf,cpspval) do k=izprange(1),izprange(2) kk = k-izprange(1)+1 plt(:,kk) = plt0(:,k) enddo call conrec(plt,jmx,jmx,nzp,fflo,hhi,fi,1,-1,-1430B,2) call labrect(gcmlat,jmx,gcmzp(izprange(1)),nzp,'LATITUDE', + 'ZP',0.) else call defphase(phase,jmx,nhtscale,nfreq,pltht,scale,nf, + cpspval) call conrec(pltht,jmx,jmx,nhtscale,fflo,hhi,fi,1,-1, + -1430B,2) call labrect(gcmlat,jmx,htscale,nhtscale,'LATITUDE', + 'HEIGHT (KM)',0.) endif call clearstr(toplab) write(toplab,"(a8,' PHASE (',i1,') UT=',f6.2)") + flab(ip),nf,utc if (iwrascii.gt.0) then call clearstr(rec80) write(rec80,"(a)") toplab(1:lenstr(toplab)) call clearstr(fieldlab) write(fieldlab,"(a8)") flab(ip)(1:lenstr(flab(ip))) endif call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.05,0.) write(toplab,"('HISTORY=',a)") histvol(ivol) call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.11, + .012) if(ihtsc.le.0) call altyax(nzp,yaxht,gcmzp(izprange(1)), + rnd,6) call frame iframe = iframe+1 if (ihtsc.le.0) then if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,plt,jmx,nzp,gcmlat, + gcmzp(izprange(1)),'LATITUDE','LN(P0/P)', + histvol(ivol),fieldlab,iframe,rec80,'timesproc', + dirascii) endif write(6,"('pltampha: frame ',i4,' field ',a8,' phase(',i1, + ') zprange=',2f8.2)") iframe,flab(ip),nf,zprange else if (iwrascii.gt.0) then call wrascii(iwrascii,luascii,pltht,jmx,nhtscale,gcmlat, + htscale,'LATITUDE','HEIGHT (KM)',histvol(ivol),fieldlab, + iframe,rec80,'timesproc',dirascii) endif write(6,"('pltampha: frame ',i4,' field ',a8,' phase(',i1, + ') htscale=',f6.2,' to ',f6.2,' by ',f5.1)") iframe, + flab(ip),nf,htscale(1),htscale(nhtscale),delht endif c c Freq loop: 100 continue c c End field loop: ip=1,nftot 200 continue call hpdeallc(pphase,ier,1) call hpdeallc(pamp,ier,1) if (ihtsc.gt.0) then call hpdeallc(ppltht,ier,1) call hpdeallc(pfht,ier,1) endif c return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine fixphase(phase,id1,id2,id3,ixn2,spv) dimension phase(id1,id2,id3) do k=1,id2 if (phase(2,k,2).ne.spv) then phase(1,k,2) = phase(2,k,2) else phase(1,k,2) = spv endif if (phase(3,k,3).ne.spv) then phase(1,k,3) = phase(3,k,3) else phase(1,k,3) = spv endif if (phase(3,k,3).ne.spv) then phase(2,k,3) = phase(3,k,3) else phase(2,k,3) = spv endif if (phase(4,k,4).ne.spv) then phase(1,k,4) = phase(4,k,4) else phase(1,k,4) = spv endif if (phase(4,k,4).ne.spv) then phase(2,k,4) = phase(4,k,4) else phase(2,k,4) = spv endif if (phase(4,k,4).ne.spv) then phase(3,k,4) = phase(4,k,4) else phase(3,k,4) = spv endif if (phase(id1-1,k,2).ne.spv) then phase(id1,k,2) = phase(id1-1,k,2) else phase(id1,k,2) = spv endif if (phase(id1-2,k,3).ne.spv) then phase(id1,k,3) = phase(id1-2,k,3) else phase(id1,k,3) = spv endif if (phase(id1-2,k,3).ne.spv) then phase(id1-1,k,3) = phase(id1-2,k,3) else phase(id1-1,k,3) = spv endif if (phase(id1-3,k,4).ne.spv) then phase(id1,k,4) = phase(id1-3,k,4) else phase(id1,k,4) = spv endif if (phase(id1-3,k,4).ne.spv) then phase(id1-1,k,4) = phase(id1-3,k,4) else phase(id1-1,k,4) = spv endif if (phase(id1-3,k,4).ne.spv) then phase(id1-2,k,4) = phase(id1-3,k,4) else phase(id1-2,k,4) = spv endif enddo if (ixn2.le.0) then do nf=1,id3 do j=1,id1 if (phase(j,2,nf).ne.spv) then phase(j,1,nf) = phase(j,2,nf) else phase(j,1,nf) = spv endif enddo enddo endif return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine defamp(amp,id1,id2,id3,plt,nf,spv) dimension amp(id1,id2,id3),plt(id1,id2) c do k=1,id2 if (amp(1,k,nf).ne.spv) then plt(1,k) = amp(1,k,nf) else plt(1,k) = spv endif do j=2,id1 if (amp(j,k,nf).ne.spv.and.amp(j-1,k,nf).ne.spv) then plt(j,k) = .5*(amp(j,k,nf)+amp(j-1,k,nf)) else plt(j,k) = spv endif enddo enddo return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine defphase(phase,id1,id2,id3,plt,scale,nf,spv) parameter (pi = 3.14159) dimension phase(id1,id2,id3),plt(id1,id2) c do k=1,id2 if (phase(1,k,nf).ne.spv) then plt(1,k) = phase(1,k,nf)*scale+12. else plt(1,k) = spv endif do j=2,id1 if (phase(j,k,nf).ne.spv.and.phase(j-1,k,nf).ne.spv) then plt(j,k) = .5*(phase(j,k,nf)+phase(j-1,k,nf))*scale+12. if (abs(phase(j,k,nf)-phase(j-1,k,nf)).gt.pi) + plt(j,k) = plt(j,k)+pi*scale else plt(j,k) = spv endif enddo enddo return end