program post c c Post-processor for the ITG program. c c c $Header: /afs/pppl.gov/u/qpliu/cvsroot/codes/gryffin/src/old94/post.f,v 1.1 1992/06/02 20:02:08 mbeer Exp $ c $Log: post.f,v $ c Revision 1.1 1992/06/02 20:02:08 mbeer c Initial revision c c Revision 3.20 1992/04/04 07:03:51 bdorland c mfr=0 option added. c (no frequency plot) c c Revision 3.19 1992/04/03 06:45:20 bdorland c *** empty log message *** c c Revision 3.18 1992/04/03 05:56:30 bdorland c ifilter message updated. c c Revision 3.17 1992/04/03 00:19:25 bdorland c multiple files accepted again. c c Revision 3.16 1992/04/02 21:01:36 bdorland c n's reindexed the way they are in itg -- long overdue label fix c c Revision 3.15 1992/04/01 23:36:32 bdorland c mfr -> 1 if md=1. c c Revision 3.14 1992/04/01 19:19:20 bdorland c mapsm added. c c include 'itg.par' c c ********* merge up to 7 output (-res) files ******** c c.... if newg=1, then introduce newgrid size and mode numbers. c Right now, only decrease mode numbers. c character*80 str,string,note integer ipage,dash(10) real wt1(7*nez),wt2(7*nez),pfdx(7*nez),delta(7*nez), * wf(7*nez,mz,nz),wg(mz,nz),wh(mz,nz),wt3(7*nez), * wt4(7*nez),wf2(7*nez,mz,nz),wd(mz,nz),wa(7*nez,15), * wb(7*nez,15),aout(mz),sdout(mz),e(lz),ek(lz), * akxa(lz),eec(lz),rmnum(mz),p2c(mz,nz),p2s(mz,nz), * p2(mz,nz),whh(lz,mz,mz,nspecz),thing(mz,nz) complex omegaz,phi0,phi1 real omega(7*nez),grate(7*nez) c include 'post.cmn' namelist/datlis/ld,ne,mfr,nfr,nspecies *,tim,shr,etai,rmu1,rmu2,rkap,etae,rmime *,dt0,dt,x0,y0,z0,nfile,kyspec,conplot,tkap,qkap,qkappar namelist/wmod1/newg,ld1,nmin1,nmax1,mmin1,mmax1 dash(1)=65535 dash(2)=43690 dash(3)=52428 dash(4)=56079 dash(5)=56159 dash(6)=29013 open(unit=6,file='postout',status='new') open(unit=7,file='y10res',status='new',form='unformatted') open(unit=5,file='y11res',status='old',form='formatted',err=7) c open all of the existing resuts files (y11res-y17res), otherwise c skip down to line 7. (At least y11 must exist): open(unit=8,file='y12res',status='old',form='formatted',err=7) open(unit=9,file='y13res',status='old',form='formatted',err=7) open(unit=10,file='y14res',status='old',form='formatted',err=7) open(unit=11,file='y15res',status='old',form='formatted',err=7) open(unit=13,file='y16res',status='old',form='formatted',err=7) open(unit=14,file='y17res',status='old',form='formatted',err=7) 7 open(unit=12,file='post.in',status='old') c input files: c c 12 post.in main namelist input file to control the run. c 5 y11res a copy of itg.res, the results file from an ITG run. c 8-11,13-14 y12res-y17res additional results files to be concatenated c to the end of y11res. c c output files: c c 6 postout main output file from post_itg, ascii. c 7 y10res a concatenation of the results files y11res-y17res c read(12,datlis) read(12,wmod1) c nd1=nmax1-nmin1+1 ldb1=ld1-1 c nff=0 net=0 ncnt=0 c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(5,1)ld,ldb,nmin,nmax,nd,lin read(5,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(5,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(5,1)(mmin(n),n=1,nd) read(5,1)(mmax(n),n=1,nd) read(5,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(5,2) dt0,dt,x0,y0,z0 n2=nspecies net=net+ne read(5,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(5,2)(timo(i),i=netb,net) read(5,2)(wpfx(i),i=netb,net) read(5,2)(gamx(i),i=netb,net) read(5,2)(wrhx(i),i=netb,net) read(5,2)(wrhy(i),i=netb,net) read(5,2)(wdux(i),i=netb,net) read(5,2)(wu2x(i),i=netb,net) read(5,2)(wp2x(i),i=netb,net) read(5,2)(wv2x(i),i=netb,net) read(5,2)(wlux(i),i=netb,net) read(5,2)(wdpx(i),i=netb,net) read(5,2)(wdvx(i),i=netb,net) read(5,2)(wenx(i),i=netb,net) read(5,2)(wdfx(i),i=netb,net) read(5,2)(pcerr(i),i=netb,net) read(5,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(5,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(5,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(5,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(5,2)(wakx(i),i=netb,net) read(5,2)(waky(i),i=netb,net) read(5,2)(therm(i),i=1,ld) read(5,2)(wxsp(i),i=netb,net) read(5,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(5,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(5,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(5,2)(timf(i),i=1,ntotal) read(5,2)avgflag read(5,5)date read(5,51)note write(6,5)date netm=net-1 do 11 n=1,nd do 11 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 11 continue c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(8,1)ld,ldb,nmin,nmax,nd,lin read(8,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(8,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(8,1)(mmin(n),n=1,nd) read(8,1)(mmax(n),n=1,nd) read(8,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(8,2)dt0,dt,x0,y0,z0 net=net+ne read(8,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(8,2)(timo(i),i=netb,net) read(8,2)(wpfx(i),i=netb,net) read(8,2)(gamx(i),i=netb,net) read(8,2)(wrhx(i),i=netb,net) read(8,2)(wrhy(i),i=netb,net) read(8,2)(wdux(i),i=netb,net) read(8,2)(wu2x(i),i=netb,net) read(8,2)(wp2x(i),i=netb,net) read(8,2)(wv2x(i),i=netb,net) read(8,2)(wlux(i),i=netb,net) read(8,2)(wdpx(i),i=netb,net) read(8,2)(wdvx(i),i=netb,net) read(8,2)(wenx(i),i=netb,net) read(8,2)(wdfx(i),i=netb,net) read(8,2)(pcerr(i),i=netb,net) read(8,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(8,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(8,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(8,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(8,2)(wakx(i),i=netb,net) read(8,2)(waky(i),i=netb,net) read(8,2)(therm(i),i=1,ld) read(8,2)(wxsp(i),i=netb,net) read(8,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(8,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(8,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(8,2)(timf(i),i=1,ntotal) read(8,2)x read(8,5)date read(8,51)note write(6,5)date netm=net-1 do 12 n=1,nd do 12 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 12 continue c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(9,1)ld,ldb,nmin,nmax,nd,lin read(9,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(9,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(9,1)(mmin(n),n=1,nd) read(9,1)(mmax(n),n=1,nd) read(9,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(9,2)dt0,dt,x0,y0,z0 net=net+ne read(9,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(9,2)(timo(i),i=netb,net) read(9,2)(wpfx(i),i=netb,net) read(9,2)(gamx(i),i=netb,net) read(9,2)(wrhx(i),i=netb,net) read(9,2)(wrhy(i),i=netb,net) read(9,2)(wdux(i),i=netb,net) read(9,2)(wu2x(i),i=netb,net) read(9,2)(wp2x(i),i=netb,net) read(9,2)(wv2x(i),i=netb,net) read(9,2)(wlux(i),i=netb,net) read(9,2)(wdpx(i),i=netb,net) read(9,2)(wdvx(i),i=netb,net) read(9,2)(wenx(i),i=netb,net) read(9,2)(wdfx(i),i=netb,net) read(9,2)(pcerr(i),i=netb,net) read(9,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(9,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(9,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(9,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(9,2)(wakx(i),i=netb,net) read(9,2)(waky(i),i=netb,net) read(9,2)(therm(i),i=1,ld) read(9,2)(wxsp(i),i=netb,net) read(9,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(9,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(9,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(9,2)(timf(i),i=1,ntotal) read(9,2)x read(9,5)date read(9,51)note write(6,5)date netm=net-1 do 13 n=1,nd do 13 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 13 continue c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(10,1)ld,ldb,nmin,nmax,nd,lin read(10,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(10,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(10,1)(mmin(n),n=1,nd) read(10,1)(mmax(n),n=1,nd) read(10,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(10,2)dt0,dt,x0,y0,z0 net=net+ne read(10,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(10,2)(timo(i),i=netb,net) read(10,2)(wpfx(i),i=netb,net) read(10,2)(gamx(i),i=netb,net) read(10,2)(wrhx(i),i=netb,net) read(10,2)(wrhy(i),i=netb,net) read(10,2)(wdux(i),i=netb,net) read(10,2)(wu2x(i),i=netb,net) read(10,2)(wp2x(i),i=netb,net) read(10,2)(wv2x(i),i=netb,net) read(10,2)(wlux(i),i=netb,net) read(10,2)(wdpx(i),i=netb,net) read(10,2)(wdvx(i),i=netb,net) read(10,2)(wenx(i),i=netb,net) read(10,2)(wdfx(i),i=netb,net) read(10,2)(pcerr(i),i=netb,net) read(10,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(10,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(10,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(10,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(10,2)(wakx(i),i=netb,net) read(10,2)(waky(i),i=netb,net) read(10,2)(therm(i),i=1,ld) read(10,2)(wxsp(i),i=netb,net) read(10,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(10,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(10,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(10,2)(timf(i),i=1,ntotal) read(10,2)x read(10,5)date read(10,51)note write(6,5)date netm=net-1 do 14 n=1,nd do 14 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 14 continue c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(11,1)ld,ldb,nmin,nmax,nd,lin read(11,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(11,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(11,1)(mmin(n),n=1,nd) read(11,1)(mmax(n),n=1,nd) read(11,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(11,2)dt0,dt,x0,y0,z0 net=net+ne read(11,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(11,2)(timo(i),i=netb,net) read(11,2)(wpfx(i),i=netb,net) read(11,2)(gamx(i),i=netb,net) read(11,2)(wrhx(i),i=netb,net) read(11,2)(wrhy(i),i=netb,net) read(11,2)(wdux(i),i=netb,net) read(11,2)(wu2x(i),i=netb,net) read(11,2)(wp2x(i),i=netb,net) read(11,2)(wv2x(i),i=netb,net) read(11,2)(wlux(i),i=netb,net) read(11,2)(wdpx(i),i=netb,net) read(11,2)(wdvx(i),i=netb,net) read(11,2)(wenx(i),i=netb,net) read(11,2)(wdfx(i),i=netb,net) read(11,2)(pcerr(i),i=netb,net) read(11,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(11,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(11,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(11,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(11,2)(wakx(i),i=netb,net) read(11,2)(waky(i),i=netb,net) read(11,2)(therm(i),i=1,ld) read(11,2)(wxsp(i),i=netb,net) read(11,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(11,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(11,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(11,2)(timf(i),i=1,ntotal) read(11,2)x read(11,5)date read(11,51)note write(6,5)date netm=net-1 do 15 n=1,nd do 15 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 15 continue c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(13,1)ld,ldb,nmin,nmax,nd,lin read(13,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(13,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(13,1)(mmin(n),n=1,nd) read(13,1)(mmax(n),n=1,nd) read(13,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(13,2)dt0,dt,x0,y0,z0 net=net+ne read(13,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(13,2)(timo(i),i=netb,net) read(13,2)(wpfx(i),i=netb,net) read(13,2)(gamx(i),i=netb,net) read(13,2)(wrhx(i),i=netb,net) read(13,2)(wrhy(i),i=netb,net) read(13,2)(wdux(i),i=netb,net) read(13,2)(wu2x(i),i=netb,net) read(13,2)(wp2x(i),i=netb,net) read(13,2)(wv2x(i),i=netb,net) read(13,2)(wlux(i),i=netb,net) read(13,2)(wdpx(i),i=netb,net) read(13,2)(wdvx(i),i=netb,net) read(13,2)(wenx(i),i=netb,net) read(13,2)(wdfx(i),i=netb,net) read(13,2)(pcerr(i),i=netb,net) read(13,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(13,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(13,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(13,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(13,2)(wakx(i),i=netb,net) read(13,2)(waky(i),i=netb,net) read(13,2)(therm(i),i=1,ld) read(13,2)(wxsp(i),i=netb,net) read(13,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(13,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(13,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(13,2)(timf(i),i=1,ntotal) read(13,2)x read(13,5)date read(13,51)note write(6,5)date netm=net-1 do 16 n=1,nd do 16 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 16 continue c c nff=nff+1 if(nff.gt.nfile) goto 1000 netb=net+1 ncntb=ncnt+1 read(14,1)ld,ldb,nmin,nmax,nd,lin read(14,1)nsteps,ne,nfreq,ntotal,md,nspecies,ikx read(14,1)nparmom,nperpmom,iodd,iflr,iphi00,ifilter,iexp,igradon read(14,1)(mmin(n),n=1,nd) read(14,1)(mmax(n),n=1,nd) read(14,2)tim,shr,etai,rmu1,rmu2,rkap,tiovte,vy,vz,rmime,etae read(14,2)dt0,dt,x0,y0,z0 net=net+ne read(14,2)((((ws(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((wc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((us(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((uc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((tpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((tparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((prs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((prc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((vs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((vc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((ts(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((tc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((qs(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((qc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((qpars(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((qparc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((rps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((rpc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((sps(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)((((spc(l,m,n,i),l=1,ld),m=1,md),n=1,nd),i=1,n2) read(14,2)(timo(i),i=netb,net) read(14,2)(wpfx(i),i=netb,net) read(14,2)(gamx(i),i=netb,net) read(14,2)(wrhx(i),i=netb,net) read(14,2)(wrhy(i),i=netb,net) read(14,2)(wdux(i),i=netb,net) read(14,2)(wu2x(i),i=netb,net) read(14,2)(wp2x(i),i=netb,net) read(14,2)(wv2x(i),i=netb,net) read(14,2)(wlux(i),i=netb,net) read(14,2)(wdpx(i),i=netb,net) read(14,2)(wdvx(i),i=netb,net) read(14,2)(wenx(i),i=netb,net) read(14,2)(wdfx(i),i=netb,net) read(14,2)(pcerr(i),i=netb,net) read(14,2)(((wkif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(14,2)(((wpif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(14,2)(((wvif(i,m,n),i=netb,net),m=1,md),n=1,nd) read(14,2)(((wplx(i,m,n),i=netb,net),m=1,md),n=1,nd) read(14,2)(wakx(i),i=netb,net) read(14,2)(waky(i),i=netb,net) read(14,2)(therm(i),i=1,ld) read(14,2)(wxsp(i),i=netb,net) read(14,2)(((grtmx(i,m,n),i=netb,net),m=1,md),n=1,nd) ncnt=ncnt+ntotal read(14,2)(((uctim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(14,2)(((ustim(i,m,n),i=ncntb,ncnt),m=1,md),n=1,nd) read(14,2)(timf(i),i=1,ntotal) read(14,2)x read(14,5)date read(14,51)note write(6,5)date netm=net-1 do 17 n=1,nd do 17 m=1,md grtmx(net,m,n)=grtmx(netm,m,n) 17 continue c 1000 continue c if(md.eq.1) mfr=1 c c c define the lhistory variable as in itg: c if(shr .ne. 0.0) then lhistory=sqrt(1./shr)*ld/x0 else lhistory=ld/4 endif if(iodd .eq. 2) lhistory=0. if(lhistory .gt. ld/4) lhistory=ld/4 lhistory=lhistory + ld/2+1 c if(nd.eq.1) then nmin=0 nmax=0 nddm=1 nddp=1 nrr(1)=0 else nddm=nd/2 nddp=nd/2+1 nmin=-nd/2 nmax=nd/2 c do 888 n=1,nddp nrr(n)=n-1 888 continue do 889 n=1,nddm nrr(n+nddp)=n-nddm-1 889 continue c endif c do 89 n=1,nd do 89 m=1,md mrr(m,n)=m+mmin(n)-1 c if(mrr(m,n).eq.0.and.nrr(n).eq.0) then meq=m neq=n endif c 89 continue c pi=3.141592653589793 rldb=float(ldb) c do 90 l=1,ld r(l)=x0*float(l-1)/rldb 90 continue c r0=r(ld)/2.0 c do 91 l=1,ld bdel(l)=shr*(r(l)-r0) 91 continue c do 95 l=2,ldb 95 dr(l)=r(l+1)-r(l-1) dr(1)=2.*(r(2)-r(1)) dr(ld)=2.*(r(ld)-r(ldb)) c write(6,datlis) write(6,3)ld,net write(6,4)tim,shr,etai,rmu1,rmu2,rkap,etae,rmime,nspecies c 1 format(1x,3i10) 2 format(1x,3e22.14) 3 format(x,'ld = ',i5,/ * x,'ne = ',i5) 4 format(x,'tim = ',g11.3,/ * x,'shr = ',g11.3,/ * x,'etai = ',g11.3,/ * x,'rmu1 = ',g11.3,/ * x,'rmu2 = ',g11.3,/ * x,'rkap2 = ',g11.3,/ * x,'etae = ',g11.3,/ * x,'rmime = ',g11.3,/ * x,'nspecies = ',i5,/) 5 format(1x,a30) 51 format(1x,a80) c write(6,101) (nrr(n),mmin(n),mmax(n),n=1,nd) 101 format(3x,'at n =',i4,5x,'mmin=',i4,3x,'mmax=',i4) c write(6,3000) c call width(us,uc,avdx,avdy,akx2,aky2,wg,meq,neq) write(6,1109) write(6,1101) avdx,avdy call width(prs,prc,avdx,avdy,akx2,aky2,wh,meq,neq) write(6,1102) write(6,1101) avdx,avdy c 1101 format(7x,' = ',e16.8,' = ',e16.8) 1102 format(/7x,' mode width = ( < p^2 > / < (dp/dx)^2 > )**.5 ') 1109 format(/7x,' mode width = ( < u^2 > / < (du/dx)^2 > )**.5 ') c call timeavp(wpfx,avpf,sdv,net) c write(6,1103) write(6,1104) avpf,sdv 1103 format(/7x,' time-averaged thermal flux ' ) 1104 format(7x,' avpf=',e16.8,' sdv of pf=',e16.8) pfmax=avpf+sdv pfmin=avpf-sdv write(6,1107) pfmax,pfmin 1107 format(7x,' avpfmax=',e16.8,' avpfmin=',e16.8) c c call timeavp(wenx,avpf,sdv,net) c write(6,1105) write(6,1106) avpf,sdv 1105 format(/7x,' time-averaged total energy ') 1106 format(7x,' aven=',e16.8,' sdv of en=',e16.8) pfmax=avpf+sdv pfmin=avpf-sdv write(6,1108) pfmax,pfmin 1108 format(7x,' avenmax=',e16.8,' avenmin=',e16.8) c c call timeavp(wakx,avpf,sdv,net) c write(6,1125) write(6,1126) avpf,sdv 1125 format(/7x,' time-averaged radial mode width ') 1126 format(7x,' =',e16.8,' sdv of =',e16.8) pfmax=avpf+sdv pfmin=avpf-sdv write(6,1128) pfmax,pfmin 1128 format(7x,' max=',e16.8,' min=',e16.8) c c call timeavp(waky,avpf,sdv,net) c write(6,1135) write(6,1136) avpf,sdv 1135 format(/7x,' time-averaged y mode width ') 1136 format(7x,' =',e16.8,' sdv of =',e16.8) pfmax=avpf+sdv pfmin=avpf-sdv write(6,1138) pfmax,pfmin 1138 format(7x,' max=',e16.8,' min=',e16.8) c call timeavp(wxsp,avpf,sdv,net) c write(6,1175) write(6,1176) avpf,sdv 1175 format(/7x,' time-averaged mode spread in x ') 1176 format(7x,' avxs=',e16.8,' sdv of xs=',e16.8) pfmax=avpf+sdv pfmin=avpf-sdv write(6,1178) pfmax,pfmin 1178 format(7x,' avxsmax=',e16.8,' avxsmin=',e16.8) c c do 40 i=1,net absv=abs(wxsp(i)) if(absv.gt.1.0e-33) then pfdx(i)=(pi*y0)*wpfx(i)/wxsp(i) else pfdx(i)=0.0 endif 40 continue c c call timeavp(pfdx,avpf,sdv,net) c write(6,1185) write(6,1186) avpf,sdv 1185 format(7x,' pflux / mode spread in x ') 1186 format(7x,' pfdx=',e16.8,' sdv of xs=',e16.8) pfmax=avpf+sdv pfmin=avpf-sdv write(6,1188) pfmax,pfmin 1188 format(7x,' pfdxmax=',e16.8,' pfdxmin=',e16.8) c chimax=-pfmin/(etai+1.) chiav=-avpf/(etai+1.) chimin=-pfmax/(etai+1.) write(6,1189) chimax,chiav,chimin 1189 format(7x,'chimax=',e16.8,2x,'chiav=',e16.8,2x, * 'chimin=',e16.8) c 3000 format( // ) c call wpunchp(tim,dt,net,newg) c call opngks call reset ipage=1 call set(0.0,1.0,0.0,1.0,0.0,1.0,0.0,7.5,1) call plchlq(.5,7.2,'Gyrofluid ITG Simulation',25.,0.,0.) if (lin.eq.0) then call plchlq(.1,6.0,'E x B nonlinearities included. ',12.,0.,-1.) else call plchlq(.1,6.0,'This is a linear run. ',12.,0.,-1.) endif if (igradon.eq.1) then call plchlq(.1,5.8,'Background pressure gradient maintained. ' * ,12.,0.,-1.) else call plchlq(.1,5.8,'Background pressure allowed to relax. ' * ,12.,0.,-1.) endif write(str,111) 'Lx = ',x0,' rho' 111 format(a5,f6.2,a4) call plchlq(.1,5.6,str,12.,0.,-1.) write(str,111) 'Ly = ',y0*2.0*pi,' rho' call plchlq(.1,5.4,str,12.,0.,-1.) write(str,111) 'Lz = ',z0*2.0*pi,' Ln ' call plchlq(.1,5.2,str,12.,0.,-1.) write(str,102) ld 102 format('Radial grid points = ',i3) call plchlq(.1,5.0,str,12.,0.,-1.) call plchlq(.1,4.8,'Poloidal modes: ',12.,0.,-1.) write(str,103) 'Mmin = ',mmin(nd) 103 format(a7,i2) call plchlq(.2,4.6,str,12.,0.,-1.) write(str,103) 'Mmax = ',mmax(nd) call plchlq(.2,4.4,str,12.,0.,-1.) call plchlq(.1,4.2,'Toroidal modes: ',12.,0.,-1.) write(str,103) 'Nmin = ',nmin call plchlq(.2,4.0,str,12.,0.,-1.) write(str,103) 'Nmax = ',nmax call plchlq(.2,3.8,str,12.,0.,-1.) if(tim.eq.dt*float(ne)) then write(str,104) dt0 else write(str,9104) endif 104 format('Time step = ',f5.2) 9104 format('Time step adjustable') call plchlq(.1,3.4,str,12.,0.,-1.) write(str,105) tim 105 format('Time at end of run = ',f7.2) call plchlq(.1,3.2,str,12.,0.,-1.) write(str,106) shr 106 format('Ln / Ls = ',f6.3) call plchlq(.1,2.8,str,12.,0.,-1.) write(str,107) etai 107 format('ETAi = ',f4.1) call plchlq(.1,2.6,str,12.,0.,-1.) if(nspecies.eq.2) then write(str,1122) 'Etae = ',etae call plchlq(.25,2.6,str,12.,0.,-1.) 1122 format(a7,f4.1) endif write(str,108) tiovte 108 format('Ti/Te = ',f4.1) call plchlq(.1,2.4,str,12.,0.,-1.) if (ifilter.eq.1) then write(str,109) rmu1 109 format('Perpendicular viscosities = ',f6.3) call plchlq(.1,2.2,str,12.,0.,-1.) endif write(str,115) nparmom call plchlq (.1,2.0,str,12.,0.,-1.) 115 format('Number of parallel moments: ',i2) write(str,112) nperpmom call plchlq (.1,1.8,str,12.,0.,-1.) 112 format('Number of perpendicular moments: ',i2) call plchlq (.1,1.6,note,12.,0.,-1.) c c second column of inputs c if (vy.eq.0. .and. vz.eq.0.) then call plchlq(.5,6.0,'Sheared velocity flow not included', * 12.,0.,-1.) else call plchlq(.5,6.0,'Sheared velocity flow included',12.,0.,-1.) write(str,1013) 'vy = ',vy call plchlq(.75,5.8,str,12.,0.,-1.) write(str,1013) 'vz = ',vz call plchlq(.75,5.6,str,12.,0.,-1.) 1013 format(a5,f12.6) endif write(str,103) 'iodd = ',iodd call plchlq(.5,5.4,str,12.,0.,-1.) write(str,103) 'iflr = ',iflr call plchlq(.5,5.2,str,12.,0.,-1.) write(str,203) 'iphi00 = ',iphi00 203 format(a9,i2) call plchlq(.5,5.0,str,12.,0.,-1.) write(str,204) 'ifilter = ',ifilter call plchlq(.5,4.8,str,12.,0.,-1.) 204 format(a10,i2) write(str,2204) 'ikx = ',ikx 2204 format(a6,i2) call plchlq(.75,4.8,str,12.,0.,-1.) if(nspecies.eq.2) then call plchlq(.5,4.6,'Non-adiabatic electrons included', * 12.,0.,-1.) write(str,1121) 'Mi / me = ',rmime call plchlq(.5,4.4,str,12.,0.,-1.) 1121 format(a10,f7.2) endif call finish(ipage) c if (conplot .eq. 1) then call cnplot0(uc,us,'Potential',tim,ipage) call cnplot0(wc,ws,'Density',tim,ipage) call cnplot0(vc,vs,'Parallel Velocity',tim,ipage) call cnplot0(tparc,tpars,'Parallel Temperature',tim,ipage) call cnplot0(prc,prs,'Parallel Pressure',tim,ipage) if(nparmom.eq.4) then call cnplot0(qparc,qpars,'Par-Par Heat Flux',tim,ipage) endif call cnplot0(tc,ts,'Perpendicular Temperature',tim,ipage) if(nperpmom.eq.2) then call cnplot0(qc,qs,'Par-Perp Heat Flux',tim,ipage) endif endif c c now do final time plots (was outp) c c for now, only plot one mode if(mfr.eq.0) goto 712 mplot=mfr nplot=nfr call hmplot2(uc,us,'Potential',tim,ipage,mplot,nplot) c c make a plot of omega for this mode c do 92 i=2,ntotal phi1=cmplx(uctim(i,mplot,nplot),ustim(i,mplot,nplot)) phi0=cmplx(uctim(i-1,mplot,nplot),ustim(i-1,mplot,nplot)) if(timf(i)-timf(i-1).gt.1.e-30) then omegaz=log(phi1/phi0)/(timf(i)-timf(i-1))/(0.0,-1.0) else goto 712 endif omega(i)=-real(omegaz) grate(i)=aimag(omegaz) 92 continue omega(1)=omega(2) grate(1)=grate(2) gmin=0.0 gmax=0.0 call mnmx(ntotal/2,omega(ntotal/2),gmin,gmax) call mnmx(ntotal/2,grate(ntotal/2),gmin,gmax) gmax=1.2*max(abs(gmin),abs(gmax)) gmin=-gmax call setzer call plchlq(.5,.95,'Frequency vs. Time',20.,0.,0.) write(string,'(a,i2,a,i2,a,f6.1)') 'for the m=', * mrr(mplot,nplot), ' n=',nrr(nplot),' mode at x=',r(lhistory) call plchlq(.85,.92,string,15.,0.,0.) write(string,'(a,1pe9.2,a,1pe9.2,a)') 'final frequency =', * omega(ntotal),' + ',grate(ntotal),' i' c c Write out the final frequency to the file itg.freq: c open(unit=35,file='itg.freq',status='new',form='formatted') write(35,711)omega(ntotal),grate(ntotal) 711 format(1x,3e22.14) call plchlq(.75,.05,string,15.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timf(1),timf(ntotal),gmin,gmax) c solid line call dashdb(65535) call curve(timf,omega,ntotal) c dashed line? call dashdb(29398) call curved(timf,grate,ntotal) call finish(ipage) 712 continue c c now do all the hmplot calls (harmonics plotting) call hmplot(ws,'Density (sine)',tim,ipage) call hmplot(wc,'Density (cosine)',tim,ipage) call hmplot(us,'Potential (sine)',tim,ipage) call hmplot(uc,'Potential (cosine)',tim,ipage) call hmplot(vs,'Parallel Velocity (sine)',tim,ipage) call hmplot(vc,'Parallel Velocity (cosine)',tim,ipage) call hmplot(tpars,'Parallel Temperature (sine)',tim,ipage) call hmplot(tparc,'Parallel Temperature (cosine)',tim,ipage) call hmplot(prs,'Parallel Pressure (sine)',tim,ipage) call hmplot(prc,'Parallel Pressure (cosine)',tim,ipage) if(nparmom.eq.4)then call hmplot(qpars,'Par-Par Heat Flux (sine)',tim,ipage) call hmplot(qparc,'Par-Par Heat Flux (cosine)',tim,ipage) endif call hmplot(ts,'Perpendicular Temperature (sine)',tim,ipage) call hmplot(tc,'Perpendicular Temperature (cosine)',tim,ipage) if(nperpmom.ge.2)then call hmplot(qs,'Par-Perp Heat Flux (sine)',tim,ipage) call hmplot(qc,'Par-Perp Heat Flux (cosine)',tim,ipage) endif if(nperpmom.ge.3)then call hmplot(rps,'R Perp (sine)',tim,ipage) call hmplot(rpc,'R Perp (cosine)',tim,ipage) endif if(nperpmom.ge.4)then call hmplot(sps,'S Perp (sine)',tim,ipage) call hmplot(spc,'S Perp (cosine)',tim,ipage) endif c do 10 n=1,nd do 10 m=1,md do 10 l=1,ld whh(l,m,n,1)=us(l,m,n,1)*us(l,m,n,1)+uc(l,m,n,1)*uc(l,m,n,1) if(whh(l,m,n,1).lt.1.e-33) whh(l,m,n,1) = 0.0 10 continue c call hmplot(whh,'Potential**2',tim,ipage) c end of harmonics plotting gmin=0.0 gmax=0.0 c do 25 n=1,net if(gmax.lt.gamx(n)) gmax=gamx(n) if(gmin.gt.gamx(n)) gmin=gamx(n) 25 continue c gmax=1.1*gmax c gmin=1.1*gmin c call setzer call plchlq(.5,.95,'Growth rate',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) c solid line call dashdb(65535) call curve(timo,gamx,net) call finish(ipage) c if(shr .ne. 0.0) then c c rational surfaces vs. m with circles, mbeer c can change to ovals by setting dx,dy c call setzer call plchlq(.5,.95,'Rational Surfaces (phi)',20.,0.,0.) call plchlq(.5,.08,'radial position',15.,0.,0.) call plchlq(.02,.5,'poloidal mode number',15.,90.,0.) gmin=mmin(1) gmax=mmax(1) do 82344 n=1,nd if (float(mmax(n)).gt.gmax) gmax=float(mmax(n)) if (float(mmin(n)).lt.gmin) gmin=float(mmin(n)) 82344 continue gmax = gmax + 1. rmax=(r(ld)-r(1))/2. c Choose instantaneous amplitude or growth rate for dy if(lin.eq.1) then do 8220 n=1,nd do 8220 m=1,md thing(m,n)=grtmx(netm,m,n) 8220 continue else do 8221 n=1,nd do 8221 m=1,md thing(m,n)=alog(wkif(netm,m,n)+1.e-15) 8221 continue endif wmax=abs(thing(1,1)) wmaxnl=thing(1,1) wminnl=thing(1,1) do 8222 n=1,nd do 8222 m=1,md wmax=max(wmax,abs(thing(m,n))) wmaxnl=max(wmaxnl,abs(thing(m,n))) wminnl=max(wminnl,abs(thing(m,n))) if(mrr(m,n).ne.0) then rrat=(float(nrr(n))*y0)/(float(mrr(m,n))*shr*z0) endif if(abs(rrat).gt.rmax) rmax = abs(rrat) 8222 continue c c shift to make thing positive, mbeer c if (lin.eq.0) then do 8227 n=1,nd do 8227 m=1,md thing(m,n)=thing(m,n)+wminnl+.1 8227 continue wmax=(wminnl+wmaxnl+.01)/2. endif c c shift to make thing positive, mbeer c rmax=1.1*rmax call maps(-rmax+r0,rmax+r0,gmin,gmax) do 80517 n=1,nd do 80517 m=1,md if ((mrr(m,n).eq.0).and.(nrr(n).ne.0)) goto 80517 if (mrr(m,n).ne.0) & rrat=(float(nrr(n))*y0)/(float(mrr(m,n))*shr*z0) if (mrr(m,n).eq.0) rrat=0. dy=thing(m,n)/(2.*abs(wmax)) dx=wg(m,n) xcir0=r0+rrat+dx ycir0=float(mrr(m,n)) do 82345 l=0,40 thetac=float(l)/float(40)*6.2832 xcir1=r0+rrat+dx*cos(thetac) ycir1=float(mrr(m,n))+dy*sin(thetac) if (dy.gt.0.) then call line(xcir0,ycir0,xcir1,ycir1) else call line(xcir1,ycir1,xcir1,ycir1) endif xcir0=xcir1 ycir0=ycir1 82345 continue 80517 continue call line(r(1),gmin,r(1),gmax) call line(r(ld),gmin,r(ld),gmax) c following line is hard in NCAR c CALL YGRAXS(gmin/y0,(gmax-gmin)/(y0*4.25),gmax/y0, c & 6.,'ky rho $',-100,9.,0.) call finish(ipage) endif c if(shr .ne. 0.0) then c c rational surfaces vs. m with circles, mbeer c can change to ovals by setting dx,dy c call setzer call plchlq(.5,.95,'Rational Surfaces (pressure)',20.,0.,0.) call plchlq(.5,.08,'radial position',15.,0.,0.) call plchlq(.02,.5,'poloidal mode number',15.,90.,0.) gmin=mmin(1) gmax=mmax(1) do 8100 n=1,nd if (float(mmax(n)).gt.gmax) gmax=float(mmax(n)) if (float(mmin(n)).lt.gmin) gmin=float(mmin(n)) 8100 continue gmax = gmax + 1. rmax=(r(ld)-r(1))/2. c Choose instantaneous amplitude or growth rate for dy do 8101 n=1,nd do 8101 m=1,md if(mrr(m,n).ne.0) then rrat=(float(nrr(n))*y0)/(float(mrr(m,n))*shr*z0) endif if(abs(rrat).gt.rmax) rmax = abs(rrat) 8101 continue rmax=1.1*rmax call maps(-rmax+r0,rmax+r0,gmin,gmax) do 8102 n=1,nd do 8102 m=1,md if ((mrr(m,n).eq.0).and.(nrr(n).ne.0)) goto 8102 if (mrr(m,n).ne.0) & rrat=(float(nrr(n))*y0)/(float(mrr(m,n))*shr*z0) if (mrr(m,n).eq.0) rrat=0. dy=thing(m,n)/(2.*abs(wmax)) dx=wh(m,n) xcir0=r0+rrat+dx ycir0=float(mrr(m,n)) do 8103 l=0,40 thetac=float(l)/float(40)*6.2832 xcir1=r0+rrat+dx*cos(thetac) ycir1=float(mrr(m,n))+dy*sin(thetac) if (dy.gt.0.) then call line(xcir0,ycir0,xcir1,ycir1) else call line(xcir1,ycir1,xcir1,ycir1) endif xcir0=xcir1 ycir0=ycir1 8103 continue 8102 continue call line(r(1),gmin,r(1),gmax) call line(r(ld),gmin,r(ld),gmax) c CALL YGRAXS(gmin/y0,(gmax-gmin)/(y0*4.25),gmax/y0, c & 6.,'ky rho $',-100,9.,0.) call finish(ipage) endif c do 1570 i=1,net wt3(i)=pi*y0*wdux(i)/x0 wt4(i)=pi*y0*wp2x(i)/x0 1570 continue c call cplogp(wt3,wt1,net) call cplogp(wt4,wt2,net) c gmin=-10.0 gmax=3.0 do 171 n=1,net if(gmax.lt.wt1(n)) gmax=wt1(n) if(gmin.gt.wt1(n)) gmin=wt1(n) if(gmax.lt.wt2(n)) gmax=wt2(n) if(gmin.gt.wt2(n)) gmin=wt2(n) 171 continue call setzer call plchlq(.5,.95,'Perpendicular kinetic and thermal energy',20. * ,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call plchlq(.02,.5,'log base ten',15.,90.,0.) call maps(timo(1),timo(net),gmin,gmax) call dashdb(65535) call curve(timo,wt1,net) call dashdb(21845) call curve(timo,wt2,net) call set(0.1,0.9,0.15,0.9,0.,1.,0.,1.,1) call line(.05,.1,.1,.1) call plchlq(.12,.1,'kinetic',12.,0.,-1.) call lined(.05,.05,.1,.05) call plchlq(.12,.05,'thermal',12.,0.,-1.) call finish(ipage) c end of modifications (mbeer) c do 1573 i=1,net wt2(i)=pi*y0*wenx(i)/x0 1573 continue c gmin=0.0 gmax=0.0 do 172 n=1,net if(gmax.lt.wt2(n)) gmax=wt2(n) if(gmin.gt.wt2(n)) gmin=wt2(n) 172 continue gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Total Energy / x0',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) call curve(timo,wt2,net) call finish(ipage) c call cplogp(wt2,wt1,net) gmin=-3.5 gmax=2.0 do 173 n=1,net if(gmax.lt.wt1(n)) gmax=wt1(n) if(gmin.gt.wt1(n)) gmin=wt1(n) 173 continue call setzer call plchlq(.5,.95,'Total Energy / x0',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call plchlq(.02,.5,'log base ten',15.,90.,0.) call maps(timo(1),timo(net),gmin,gmax) call curve(timo,wt1,net) call finish(ipage) c do 1574 i=1,net wt2(i)=pi*y0*wpfx(i)/x0 1574 continue c call cplogp(wt2,wt1,net) gmin=0.0 gmax=0.0 do 174 n=1,net if(gmax.lt.wt1(n)) gmax=wt1(n) if(gmin.gt.wt1(n)) gmin=wt1(n) 174 continue gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Thermal Flux / x0',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call plchlq(.02,.5,'log base ten',15.,90.,0.) call maps(timo(1),timo(net),gmin,gmax) call curve(timo,wt1,net) call finish(ipage) c gmin=0.0 gmax=0.0 do 9174 n=1,net wt2(n)=-wt2(n) if(gmax.lt.wt2(n)) gmax=wt2(n) if(gmin.gt.wt2(n)) gmin=wt2(n) 9174 continue gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Thermal Flux / x0',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) call curve(timo,wt2,net) call finish(ipage) c do 1298 n=1,net wt2(n)=-wt2(n) 1298 continue c c Make a plot of the timestep vs. time c do 1300 n=2,net-1 delta(n)=timo(n+1)-timo(n) delta(n)=delta(n)*float(net)/float(nsteps) 1300 continue delta(1)=delta(2) delta(net)=delta(net-1) c gmin=0.0 gmax=1.1*dt0 call setzer call plchlq(.5,.95,'Time step',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) call curve(timo,delta,net) call finish(ipage) c do 1576 n=1,nd do 1576 m=1,md do 1576 i=1,net wf2(i,m,n)=pi*y0*wkif(i,m,n)/x0 1576 continue c do 20 n=1,nd do 20 m=1,md call cplogp(wf2(1,m,n),wf(1,m,n),net) 20 continue c ndex=1 call smplot(gmin,gmax,wf,net,ndex) c if(kyspec.eq.1) then call nsumavp(wf2,aout,sdout,net) c gmin=0.0 gmax=0.0 do 169 m=1,md c rmnum(m)=mrr(m,nd) if(gmax.lt.aout(m)) gmax=aout(m) if(gmin.gt.aout(m)) gmin=aout(m) 169 continue gmax=1.1*gmax gmin=1.1*gmin write(6,1676) (rmnum(m),aout(m),m=1,md) 1676 format(2x,3(' m=',f4.0,' ener=',e11.4)) c amplitude vs. mode number, mbeer call setzer call plchlq(.5,.95,'Time averaged perp. energy spectrum' &,20.,0.,0.) call plchlq(.02,.5,'mode amplitude',15.,90.,0.) call plchlq(.5,.08,'poloidal mode number',15.,0.,0.) call mapsm(-1.,float(md),gmin,gmax) do i=1,md call line(rmnum(i)-.3,0.0,rmnum(i)-.3,aout(i)) call line(rmnum(i)-.3,aout(i),rmnum(i)+.3,aout(i)) call line(rmnum(i)+.3,0.0,rmnum(i)+.3,aout(i)) enddo call finish(ipage) endif c c********************************************************** c call penrgyp(uc,phi2c,p2c) call penrgyp(us,phi2s,p2s) gmin=0.0 gmax=0.0 do 9169 n=1,nd do 9169 m=1,md rmnum(m)=mrr(m,nd) p2(m,n)=p2c(m,n)+p2s(m,n) p2(m,n)=abs(p2(m,n))+1.e-15 p2(m,n)=alog10(p2(m,n)) if(gmax.lt.p2(m,n)) gmax=p2(m,n) if(gmin.gt.p2(m,n)) gmin=p2(m,n) 9169 continue gmax=1.1*gmax gmin=1.1*gmin c CALL PAGE (11.0,8.5) c CALL AREA2D (9.,6.) c CALL HEADIN ('Phi**2 vs. ky$',100,1.3,1) c CALL XTICKS (5) c CALL XREVTK c CALL YTICKS (5) c CALL YREVTK c CALL XNAME ('mode number $',100) c CALL YNAME ('log base ten $',100) c CALL DFRAME c CALL GRAF (-1.,'scale',float(md),gmin,'scale',gmax) c do 468 n=1,nd c CALL CURVE (rmnum,p2(1,n),md,-1) c 468 continue c CALL DENDPL(-1) c do 1577 n=1,nd do 1577 m=1,md do 1577 i=1,net wf2(i,m,n)=pi*y0*wpif(i,m,n)/x0 1577 continue do 30 n=1,nd do 30 m=1,md call cplogp(wf2(1,m,n),wf(1,m,n),net) 30 continue c ndex=2 call smplot(gmin,gmax,wf,net,ndex) c call cplogp(pfdx,wt1,net) gmin=0.0 gmax=0.0 do 175 n=1,net if(gmax.lt.wt1(n)) gmax=wt1(n) if(gmin.gt.wt1(n)) gmin=wt1(n) 175 continue gmax=1.1*gmax gmin=1.1*gmin c CALL PAGE (11.0,8.5) c CALL AREA2D (9.,6.) c CALL HEADIN ('Adjusted thermal flux / x0 $',100,1.3,1) c CALL XTICKS (5) c CALL XREVTK c CALL YTICKS (5) c CALL YREVTK c CALL XNAME ('time $',100) c CALL YNAME ('log base ten $',100) c CALL DFRAME c CALL GRAF (timo(1),'scale',timo(net),gmin,'scale',gmax) c CALL CURVE (timo,wt1,net,0) c CALL DENDPL(-1) c do 1579 i=1,net wt2(i)=pfdx(i)/(etai+1.) 1579 continue c call cplogp(wt2,wt1,net) gmin=0.0 gmax=0.0 do 176 n=1,net if(gmax.lt.wt1(n)) gmax=wt1(n) if(gmin.gt.wt1(n)) gmin=wt1(n) 176 continue gmax=1.1*gmax gmin=1.1*gmin c CALL PAGE (11.0,8.5) c CALL AREA2D (9.,6.) c CALL HEADIN ('Chi $',100,1.3,1) c CALL XTICKS (5) c CALL XREVTK c CALL YTICKS (5) c CALL YREVTK c CALL XNAME ('time $',100) c CALL YNAME ('log base ten $',100) c CALL DFRAME c CALL GRAF (timo(1),'scale',timo(net),gmin,'scale',gmax) c CALL CURVE (timo,wt1,net,0) c CALL DENDPL(-1) c c mm=min0(7,nd) c mcent=meq+3 c c do 1610 m=1,mm c do 1610 i=1,net c nmode=2*m-1 c wa(i,m)=pi*y0*wkif(i,mcent,nmode)/x0 c 1610 continue c nmode=2*mm-1 c c do 1611 m=1,mm c call cplogp(wa(1,m),wb(1,m),net) c 1611 continue c c gmin=-11.0 c gmax=0.0 c c call saplot(gmin,gmax,timo,wb,net,mm) c call setch(70.,2.,1,1,1,0) c write(100,37) c call setch(2.,60.,1,1,1,1) c write(100,38) c call setch(1.,1.,1.,0,2,0) c write(100,1612) mrr(mcent,neq),nrr(1),nrr(nmode) c call frame c 1612 format('PERP.K.E. OF THE M=',I2,',N= FROM',I3,' TO',I3, c * ' MODES') c c call hmplot(prc,6,tim) c c call cnplot0(uc,us,2,tim) c call clsgks call exit end