program post c c Post-processor for the ITG program, makes plots, etc. c c ********* merge up to 'nfile' output (*.nc) files ******** c ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! c c use mp use itg_data use fields use diagnostics_arrays use graphics use plots use io use gryffin_layouts use command_line implicit none include 'netcdf.inc' c local post variables: character*80 str,filename,string character*80, dimension(:), allocatable :: filenames integer ipage,status real tim,dt1 integer net ! total # of time points from concatenated runs real, allocatable, dimension(:) :: wt1, wt2, delta, fphi, . aout, sdout, rnnum, rmnum, qlflux, dreta, drtor, . wupsi, dapar, dator, dtot real, allocatable, dimension(:,:) :: wt3, thing real, allocatable, dimension(:,:,:) :: wf2, whh complex omegaz,phi0,phi1 real, allocatable, dimension(:) :: omega, grate, aveth, avekp2, avekx2 real, allocatable, dimension(:) :: realp,imagp real norm,gdx,g05caE,fave,x1 ! real worky ! This was undefined in the previous version?? ! real*8, allocatable, dimension(:) :: xxx,yyy,work,xy,yy,xx,yx,workx real, allocatable, dimension(:) :: xxx,yyy,work,xy,yy,xx,yx,workx, worky real, allocatable, dimension(:) :: cxpr, cxpi, cypr, cypi, ctt real, allocatable, dimension(:) :: cxa, cya, cza, cxb, cyb real, allocatable, dimension(:) :: czb, cxc, cyc, czf, cze, czd real cxm,cym,czm,czcmin,czcmax real, allocatable, dimension(:,:,:) :: gamxk real, allocatable, dimension(:) :: gamxmax, gamxmin, t_fac real kxmax ! integer*4 ifail integer ifail integer, dimension(:), allocatable :: filenum, nt_end, ntimfp integer iharm,igrowth,igraph,idum,plot_drives integer kyspec logical file_exists logical skip_correlations logical skip_non_avg_correlations c very local variables: real a,ave,avpf,barwid,chiav,chimax,chimin,conplot,cx0 real cxmax,cxmin,cy0,cymax,cymin,cz0,czdmax,czdmin,czemax real czemin,czmax,czmin,dx,dxtot,dy,gmax,gmin real pfmax,pfmin,rldb,rmax,rp,rrat,sdv,sps,thetac real wmax,xcir0,xcir1,ycir0,ycir1 real xmax, xmin real xxmin, xxmax, yymin, yymax character*10 fx, fy integer i,j,l,m,mplot,mpts,m_min,n,ncnt integer ndpmax,ndpmaxerr integer nfile,nfile0,nox,noxm,noy,noym integer nplot,npts,nr,nx,nxdiv,ny,nydiv integer nsteps integer nff integer icount, ierr character command*80 interface subroutine mnmx(n,vector,vmin,vmax) integer :: n real, dimension(:) :: vector real :: vmin, vmax end subroutine mnmx end interface interface subroutine cplogp(wvec,wt,ne) real, dimension(:) :: wvec, wt integer :: ne end subroutine cplogp end interface interface subroutine nsumavp(wzs,aout,sdout,net,fave) real, dimension(:) :: aout, sdout real, dimension(:,:,:) :: wzs integer net real :: fave end subroutine nsumavp end interface interface subroutine timeavp(seq,av,sdv,net,fave,smin,smax) real, dimension(:) :: seq real :: av, sdv, fave, smin, smax integer net end subroutine timeavp end interface ! 'dt' cannot appear in the namelist since it is dynamically allocated. namelist/datlis/ld,md,nd,nffty,nfftx,ne,mfr,nfr,nspecies, * tim,shr,qsf,epsn,eps,alpha,etai,nuii,rmu1,etae,rmime, * dt0,x0,y0,z0,nfile,kyspec,conplot,rmass,charge,n_I,tau, * eta,eta_par,nsteps namelist/options/iharm,igrowth,fave,igraph,plot_drives,binary, * ascii c Fortran files used, and their associated unit numbers: c c input files (where * is the run name): c c 12 *.pin postc's namelist input file (separate from c itg's namelist input file *.in) c 21 *.nc the results file from an ITG run. c 22-29 *.res2 - *.res9 additional results files c to be concatenated to the end of *.res. c c output files: c c 9 *.post main output file from post_itg, ascii. c 2 *.m The NCAR/GKS graphics metafile (view with ctrans or ictrans) c 3 (Used by NCAR internal routines) c c On most unix computers, unit 5 is standard input, unit 6 is standard c output, and unit 0 is standard error, which all point to the terminal c unless redirected by the shell. It is thus best not to use units 0,5,6, c (nor units 100, 101, and 102 on the cray). c c (at one time there was also an option to write a y10res file which was c the concatenation of the results files y11res-y17res.) call init_mp ! initialize MPI c determine the name of the namelist input file: icount=iargc() if(icount >= 1) then call cl_getarg(1,runname,lrunname,ierr) else write(*,*) 'What is the run name?' read(*,*) runname lrunname=index(runname,' ')-1 endif c Check to make sure the postc input file is there: open(unit=12,file=runname(1:lrunname)//'.pin',status='old', > err=98) go to 99 98 write(*,*) 'Unable to find postc namelist file ' > //runname(1:lrunname)//'.pin' write(*,*) 'Will try the post.in default file instead.' command='cp post.in '//runname(1:lrunname)//'.pin'//char(0) call ishell(command) 99 close(unit=12) write(*,*) 'Running postc for ', runname(1:lrunname) c Run a shell script to strip out tab characters and comments, which c some compilers can't handle, from the namelist input file: command='gfbin/nmlstrip '//runname(1:lrunname)//'.pin '//char(0) call ishell(command) c default values for post.in namelist variables: mfr=2 nfr=1 nfile=1 kyspec=1 conplot=0 nsteps=1 ascii=.false. binary=.false. igraph=1 iharm=0 igrowth=0 fave=0.75 plot_drives=1 md = -1 nffty = -2 nd = -1 nfftx = -2 c Read the input namelist: open(unit=12,file=runname(1:lrunname)//'.pins',status='old') igraph=1 read(12,datlis) read(12,options) close(unit=12) open(unit=1,file=runname(1:lrunname)//'.freq') open(unit=9,file=runname(1:lrunname)//'.post') c open all of the existing results files (*.nc, and any additional c *.nc2 to *.nc9 results files to be concatenated on the end). c At least the first input file *.nc must exist: allocate (filenum(nfile), ntimfp(nfile), nt_end(nfile)) allocate (filenames(nfile)) c nfile0 counts the actual number of restart files concatenated. c nfile is the maximum number of restart files. nfile0=0 do i=1,nfile filenum(i)=20+i if(i == 1) then if(binary) then filename = runname(1:lrunname)//'.bres' elseif(ascii) then filename = runname(1:lrunname)//'.res' else filename = runname(1:lrunname)//'.nc' endif else if(binary) then write(filename,993) runname(1:lrunname),i 993 format(a,'.bres',i1) elseif(ascii) then write(filename,994) runname(1:lrunname),i 994 format(a,'.res',i1) else write(filename,995) runname(1:lrunname),i 995 format(a,'.nc',i1) endif endif filenames(i)=filename inquire(file=filenames(i),exist=file_exists) if(.not. file_exists) goto 7 nfile0=nfile0+1 enddo 7 continue if (nfile0 <= 0) call aborter(6, > 'postc error: Results file '//filenames(1)//' doesn''t exist') net=0 ncnt=0 do nff=1, nfile0 c c read a *.nc results files. The format of this file is given c by wpunchnc.f: c if(binary) then call bwread(filenames(nff), net, ncnt, tim, dt1, nsteps, nfile) elseif(ascii) then call wread(filenames(nff), net, ncnt, tim, dt1, nsteps, nfile) else call wreadnc(filenames(nff), net, ncnt, tim, dt1, nsteps, nfile) endif c Copy some things over from itg format to post-processor format: grtmx(net,:,:)=grtmx(net-1,:,:) ntimfp(nff)=ntimf nt_end(nff)=ncnt write(9,5)date enddo ! end of nff loop over all y*res input files allocate (aveth(md), avekp2(md), avekx2(md)) ! allocate (mrr(md, nd)) allocate(omega(ncnt), grate(ncnt)) allocate(realp(net), imagp(net)) c some arrays not allocated if electrons were not included... if(.not.allocated(fluxemn)) then allocate (fluxemn(net,md,nd),qfluxemn(net,md,nd)) fluxemn = 0. qfluxemn = 0. endif c write(*,datlis) if(md == 1) mfr=1 c define the lhistory variable as in itg: lhistory= float(ld)*.6+1 if(nd == 1) then nmin=0 nmax=0 nddm=1 nddp=1 allocate (nrr(1)) nrr=0 else allocate (nrr(nd)) nddm=nd/2 nddp=nd/2+1 nmin=-nd/2 nmax=nd/2 do n=1,nddp nrr(n)=n-1 enddo do n=1,nddm nrr(n+nddp)=n-nddm-1 enddo endif do n=1,nd do m=1,md mrr(m,n)=m+mmin(n)-1 if(mrr(m,n) == 0.and.nrr(n) == 0) then meq=m neq=n endif enddo enddo c Define default value for dt(m) and time(m) if they weren't set: if(icrit == 0) then dt = dt1 time = tim endif allocate (t_fac(md)) t_fac = dt/dt0 rldb=float(ldb) do l=1,ld r(l)=2.*x0*float(l-1)/rldb-x0 enddo c moved center of theta axis to zero, r0=theta0, mbeer do m=1,md do n=1,nd if (mrr(m,n) == 0.or.lin == 1) then if (nmax /= 0) then r0(m,n)=z0*float(nrr(n))/nmax else r0(m,n)=z0*float(nrr(n)) endif else if (nmax /= 0) then r0(m,n)=z0*float(nrr(n))/nmax/float(mrr(m,n)) else r0(m,n)=z0*float(nrr(n))/float(mrr(m,n)) endif endif enddo enddo do l=2,ldb dr(l)=r(l+1)-r(l-1) enddo dr(1)=2.*(r(2)-r(1)) dr(ld)=2.*(r(ld)-r(ldb)) allocate (qlflux(nspecies)) write(9,*) 'dt does not appear in the namelist because it is ', . 'dynamically allocated.' write(9,datlis) do i=1,nspecies qlflux(i)=-wpfx(net-1,i)/phirms(net-1)**2 enddo open(unit=19,file=runname(1:lrunname)//'.flux') write(19,*) net do i=1,net write(19,*) timo(i)*etai,-wpfx(i,1)/etai enddo close(19) write(9,3) ld, net write(9,4) tim, shr, qsf, epsn, alpha, etai, rmu1, etae, rmime, nspecies 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,'qsf = ',g11.3,/ * x,'epsn = ',g11.3,/ * x,'alpha = ',g11.3,/ * x,'etai = ',g11.3,/ * x,'rmu1 = ',g11.3,/ * x,'etae = ',g11.3,/ * x,'rmime = ',g11.3,/ * x,'nspecies = ',i5,/) 5 format(1x,a30) 51 format(1x,a80) write(9,101) (nrr(n),mmin(n),mmax(n),n=1,nd) 101 format(3x,'at n =',i4,5x,'mmin=',i4,3x,'mmax=',i4) write(9,3000) i=1+(1.-fave)*(net-1) write(9,1133) timo(i),timo(net) 1133 format(/7x,'time averages from ',f9.3,' to ',f9.3,/) do i=1,nspecies call timeavp(wpfx(:,i),avpf,sdv,net,fave,pfmin,pfmax) write(9,1103) i write(9,1104) avpf,sdv write(9,1107) pfmax,pfmin enddo 1103 format(/7x,' time-averaged thermal flux of species ',i1) 1104 format(7x,' avpf=',e16.8,' sdv of pf=',e16.8) 1107 format(7x,' max=',e16.8,' min=',e16.8) chimax=-pfmin/(etai+1.e-15) chiav=-avpf/(etai+1.e-15) chimin=-pfmax/(etai+1.e-15) write(9,1189) chimax,chiav,chimin 1189 format(7x,'chimax=',e16.8,2x,'chiav=',e16.8,2x, * 'chimin=',e16.8) call timeavp(wenx,avpf,sdv,net,fave,pfmin,pfmax) write(9,1105) write(9,1106) avpf,sdv 1105 format(/7x,' time-averaged total energy ') 1106 format(7x,' aven=',e16.8,' sdv of en=',e16.8) write(9,1108) pfmax,pfmin 1108 format(7x,' max=',e16.8,' min=',e16.8) call timeavp(phirms,avpf,sdv,net,fave,pfmin,pfmax) write(9,1325) write(9,1326) avpf,sdv 1325 format(/7x,' time-averaged Phi RMS ') 1326 format(7x,' phirms=',e16.8,' sdv of phirms=',e16.8) write(9,1328) pfmax,pfmin 1328 format(7x,' phimax=',e16.8,' phimin=',e16.8) call timeavp(wakx,avpf,sdv,net,fave,pfmin,pfmax) write(9,1125) write(9,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(9,1128) pfmax,pfmin 1128 format(7x,' max=',e16.8,' min=',e16.8) call timeavp(waky,avpf,sdv,net,fave,pfmin,pfmax) write(9,1135) write(9,1136) avpf,sdv 1135 format(/7x,' time-averaged y mode width ') 1136 format(7x,' =',e16.8,' sdv of =',e16.8) write(9,1138) pfmax,pfmin 1138 format(7x,' max=',e16.8,' min=',e16.8) call timeavp(wxsp,avpf,sdv,net,fave,pfmin,pfmax) write(9,1175) write(9,1176) avpf,sdv 1175 format(/7x,' time-averaged mode spread in x ') 1176 format(7x,' avxs=',e16.8,' sdv of xs=',e16.8) write(9,1178) pfmax,pfmin 1178 format(7x,' avxsmax=',e16.8,' avxsmin=',e16.8) c Quasilinear stuff write(9,*) write(9,*) 'Quasilinear stuff' write(9,*) ' norm. to vti/Lne, rhoi^2 vti/Lne^2 ne0, ', . 'rhoi^2 vti/Lne^2 pe0 ' write(9,*) n=nfr write(9,*) 'n= ',n allocate (gamxmax(md), gamxmin(md)) do m=2,md gamxmin(m)=mgamx(m,n,net) gamxmax(m)=mgamx(m,n,net) do i=net-10,net if (mgamx(m,n,i) > gamxmax(m)) gamxmax(m)=mgamx(m,n,i) if (mgamx(m,n,i) < gamxmin(m)) gamxmin(m)=mgamx(m,n,i) enddo write(9,3331) m,float(m-1)/y0,gamxmin(m),gamxmax(m) enddo 3331 format(i4,' ky=',f6.3,' gamma=',f8.4,' to ',f8.4) write(9,*) 'Particle fluxes: electrons, ion species...' do m=2,md write(9,3332) m,float(m-1)/y0,fluxemn(net,m,n)/phisq(net,m,n), . (fluximn(net,m,n,i)/phisq(net,m,n),i=1,nspecies) enddo 3332 format(i4,' ky=',f6.3,' flux='48f8.4) write(9,*) 'Heat fluxes: electrons, ion species...' do m=2,md write(9,3332) m,float(m-1)/y0,qfluxemn(net,m,n)/phisq(net,m,n), . (qfluximn(net,m,n,i)/phisq(net,m,n),i=1,nspecies) enddo write(9,*) 'Phi^2' do m=2,md write(9,3333) m,float(m-1)/y0,phisq(net,m,n) enddo 3333 format(i4,' ky=',f6.3,' phi^2=',e12.4) 3000 format( // ) allocate(wt1(net), wt2(net)) if(igraph == 1) then allocate (dash(max(md,nspecies,7))) c c Initialize the dash(n) array to provide the dash-line patterns which c label different m modes. (Do it here, and pass it through post.cmn c in order to get consistent line-types for a given m in different c plotting subroutines...) c c dash is a 16-bit pattern describing the dash pattern: dash(1)=65535 ! 2^16-1 is a solid line dash(2)=43690 dash(3)=52428 dash(4)=56079 dash(5)=56159 dash(6)=29013 c fill in the rest of the dash-line patterns randomly: call g05cbE(35) do m=7,max(md,nspecies,7) dash(m)=int(g05caE(a)*32768)+32768 c write(*,*) 'dash(',m,')= ',dash(m) enddo c Initialize the GKS graphics system, writing to a *.m meta file: call gopks(6,idum) filename = runname(1:lrunname)//'.m' call gesc(-1391,1,filename,0,idum,filename) call gopwk(1,2,1) ! 2=fortran unit # used for the gmeta file call gacwk(1) call reset ! should be called before each picture? 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 == 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 == 1) then call plchlq(.1,5.8,'Background pressure gradient maintained. ' * ,12.,0.,-1.) else call plchlq(.1,5.8,'Background pressure allowed to evolve. ' * ,12.,0.,-1.) endif write(str,111) 'Lr/2 = ',pi*y0/(abs(shr)*z0)*float(nd-1)/2.,' rho' 111 format(a11,f8.2,a4) call plchlq(.1,5.6,str,12.,0.,-1.) write(str,111) 'Ltheta/2 = ',y0*pi,' rho' call plchlq(.1,5.4,str,12.,0.,-1.) c if (iperiod == 2) then c write(str,111) 'Lz/2 = ',xp,' ' c else write(str,111) 'Lz/2 = ',x0,' ' c endif 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,i3) 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 == dt1*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.) 1061 format('r/R = ',f6.3,' nuii = ',f7.4) write(str,1061) eps,nuii call plchlq(.1,3.0,str,12.,0.,-1.) write(str,106) shr,qsf,epsn,alpha 106 format('shat=',f6.3,' q=',f6.3,' Ln/R=',f6.3,' alpha=',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 == 2) call plchlq(.25,2.6,'Multispecies',12.,0.,-1.) write(str,108) tiovte 108 format('Ti/Te = ',f4.1) call plchlq(.1,2.4,str,12.,0.,-1.) if (ifilter == 1) then write(str,109) rmu1 call plchlq(.1,2.2,str,12.,0.,-1.) endif 109 format('Perpendicular viscosities = ',f6.3) 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.) 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 == 2) then call plchlq(.5,4.6,'Multiple ion species:',12.,0.,-1.) call plchlq(.5,4.4,'See .post file',12.,0.,-1.) endif call finish(ipage, date) if (conplot == 1) then call cnplot0(potential,'Potential ',tim,ipage) call cnplot0(density(:,:,:,1),'Density ',tim,ipage) call cnplot0(u_par(:,:,:,1),'Parallel Velocity ',tim,ipage) call cnplot0(t_par(:,:,:,1),'Parallel Temp ',tim,ipage) if(nparmom == 4) then call cnplot0(q_par(:,:,:,1),'Par-Par Heat Flux ',tim,ipage) endif call cnplot0(t_perp(:,:,:,1),'Perp Temp ',tim,ipage) if(nperpmom == 2) then call cnplot0(q_perp(:,:,:,1),'Par-Perp Heat Flux',tim,ipage) endif endif c now do final time plots (was outp) c for now, only plot one mode if(mfr == 0) goto 712 mplot=mfr nplot=nfr mplot2=mplot nplot2=nplot ! used in hmplot, hmplotk call hmplot2(potential,'Potential',tim,ipage,mplot,nplot) call hmplot2(apar,'A_par',tim,ipage,mplot,nplot) c plot mode amp. vs. time (omega plot useless nonlinearly) do i=1,ncnt c GWH: ?? signs chosen to be backwards compatible with previous plots. c Should eventually change: omega(i)=-real(utim(i,mplot,nplot)) grate(i)=aimag(utim(i,mplot,nplot)) enddo gmin=0.0 gmax=0.0 call mnmx(ncnt/2,omega(ncnt/2:),gmin,gmax) call mnmx(ncnt/2,grate(ncnt/2:),gmin,gmax) gmax=1.2*max(abs(gmin),abs(gmax)) gmin=-gmax call setzer call plchlq(.5,.95,'Amplitude 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.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timf(1),timf(ncnt),gmin,gmax) c solid line call dashdb(65535) call curve(timf,omega,ncnt) c dashed line? call dashdb(29398) call curved(timf,grate,ncnt) call finish(ipage, date) endif !end of igraph=1 section c make a plot of omega for this mode mplot=mfr nplot=nfr do i=2,ncnt phi1=utim(i,mplot,nplot) phi0=utim(i-1,mplot,nplot) if(timf(i)-timf(i-1) > 1.e-10) then omegaz=ii*log(phi1/phi0)/(timf(i)-timf(i-1))/t_fac(mplot) else goto 712 endif omega(i)=-real(omegaz) grate(i)=aimag(omegaz) enddo omega(1)=omega(2) grate(1)=grate(2) gmin=0.0 gmax=0.0 call mnmx(ncnt/2,omega(ncnt/2:),gmin,gmax) call mnmx(ncnt/2,grate(ncnt/2:),gmin,gmax) gmax=1.2*max(abs(gmin),abs(gmax)) gmin=-gmax if(igraph == 1) then 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(ncnt),' + ',grate(ncnt),' i' endif !end of igraph=1 section c Write out the final frequencies to the file itg.freq: mplot=mfr nplot=nfr write(9,*) 'Frequencies: ' write(9,*) 'theta_0(mplot)=',r0(mplot,nplot) write(1,*) write(1,*) 'Frequencies: theta_0(mplot)= ',r0(mplot,nplot) i=ncnt m_min=2 if(md == 1) m_min=1 do m=m_min,md phi1=utim(i,m,nplot) phi0=utim(i-1,m,nplot) if(timf(i)-timf(i-1) > 1.e-10) then omegaz=ii*log(phi1/phi0)/(timf(i)-timf(i-1))/t_fac(m) else omegaz=666. endif write(9,1000) float(m-1)/y0,real(omegaz),aimag(omegaz) write(1,1000) float(m-1)/y0,-real(omegaz),aimag(omegaz) enddo 1000 format(3(f10.5)) if(igraph == 1) then call plchlq(.75,.05,string,15.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timf(1),timf(ncnt),gmin,gmax) c solid line call dashdb(65535) call curve(timf,omega,ncnt) c dashed line? call dashdb(29398) call curved(timf,grate,ncnt) call finish(ipage, date) endif ! end of igraph == 1 section 712 continue c now do all the hmplot calls (harmonics plotting) if(igraph == 1 .and. iharm == 1) then call hmplotcs(potential,'Potential',tim,ipage) if (epse /= 0.0) then call hmplotc(phi_bk,'Kappa Avg. of Bounce Avg. Phi',tim,ipage) call hmplotc(e_denk,'Kappa Avg. of n_e',tim,ipage) call hmplotk(phi_ba,'Bounce Avg. Phi',tim,ipage) c GWH: someday rename this to Trapped electron density: call hmplotk(e_density,'Electron Density',tim,ipage) call hmplotk(e_p,'Electron Pressure',tim,ipage) call hmplotk(e_r,'Electron R',tim,ipage) if (nemom == 4) then call hmplotk(e_t,'Electron T',tim,ipage) endif endif do i=1,nspecies call hmplotcs(density(:,:,:,i),'Density',tim,ipage) call hmplotcs(u_par(:,:,:,i),'Parallel Velocity',tim,ipage) call hmplotcs(t_par(:,:,:,i),'Parallel Temperature',tim,ipage) if(nparmom == 4) . call hmplotcs(q_par(:,:,:,i),'Par-Par Heat Flux',tim,ipage) call hmplotcs(t_perp(:,:,:,i),'Perp Temperature',tim,ipage) if(nperpmom >= 2) . call hmplotcs(q_perp(:,:,:,i),'Par-Perp Heat Flux',tim,ipage) c if(nperpmom >= 3)then c call hmplotcs(rp(:,:,:,i),'R Perp',tim,ipage) c endif c if(nperpmom >= 4)then c call hmplotcs(sps(:,:,:,i),'S Perp',tim,ipage) c endif enddo allocate(whh(ld,md,nd)) do n=1,nd do m=1,md do l=1,ld whh(l,m,n)=cabs(potential(l,m,n))**2 if(whh(l,m,n) < 1.e-30) whh(l,m,n) = 0.0 enddo enddo enddo call hmplot(whh,'Potential**2',tim,ipage) do n=1,nd do m=1,md do l=1,ld if (whh(l,m,n) >= 0.0) then whh(l,m,n)=sqrt(whh(l,m,n)) else whh(l,m,n)=0.0 endif if(whh(l,m,n) < 1.e-30) whh(l,m,n) = 0.0 enddo enddo enddo call hmplot(whh,'Abs(Potential)',tim,ipage) endif c c end of harmonics plotting c c plot growth rates vs. time c allocate(gamxk(net, md, nd)) if (igrowth == 1 .and. igraph == 1) then gmin=0.0 gmax=0.0 n=nfr do m=1,md gamxk(:,m,n)=mgamx(m,n,:) do i=net/2,net if(gmax < mgamx(m,n,i)) gmax=mgamx(m,n,i) if(gmin > mgamx(m,n,i)) gmin=mgamx(m,n,i) enddo enddo call setzer call plchlq(.5,.95,'Growth rates',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) c solid line call dashdb(65535) n=nfr do m=1,md do i=1,net gamx(i)=mgamx(m,n,i) if (gamx(i) > gmax) gamx(i)=gmax if (gamx(i) < gmin) gamx(i)=gmin enddo call dashdb(dash(m)) call curved(timo,gamx,net) enddo call finish(ipage, date) c c growth rate contour plot c do n=1,nd do m=1,md do i=1,net gamxk(i,m,n)=mgamx(m,n,i) enddo enddo enddo call cnplotk(gamxk,'Growth rates ',tim,ipage,net,fave) endif c c plot flows vs. time (linear) c if(igraph == 1) then gmin=0.0 gmax=0.0 call mnmx(net,utor,gmin,gmax) call mnmx(net,upol,gmin,gmax) call setzer call plchlq(.5,.95,'Poloidal and Toroidal flow (dashed)',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) c solid line call dashdb(21845) call curve(timo,upol,net) call curved(timo,utor,net) call finish(ipage, date) c c plot flows vs. time (log) c gmin=0.0 gmax=0.0 utor = alog(abs(utor)+1.e-15) upol = alog(abs(upol)+1.e-15) call mnmx(net,utor,gmin,gmax) call mnmx(net,upol,gmin,gmax) call setzer call plchlq(.5,.95, > 'Poloidal and Toroidal flow (dashed) log_e',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) c solid line call dashdb(21845) call curve(timo,upol,net) call curved(timo,utor,net) call finish(ipage, date) c FSA phi gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(phih0(i)) imagp(i)=aimag(phih0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave Phi',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA upar gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(uparh0(i)) imagp(i)=aimag(uparh0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave u_par',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA den gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(denh0(i)) imagp(i)=aimag(denh0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave n',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA tpar gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(tparh0(i)) imagp(i)=aimag(tparh0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave Tpar',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA tperp gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(tperph0(i)) imagp(i)=aimag(tperph0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave Tperp',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA qpar gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(qparh0(i)) imagp(i)=aimag(qparh0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave qpar',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA qperp gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(qperph0(i)) imagp(i)=aimag(qperph0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave qperp',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA v_E gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(veh0(i)) imagp(i)=aimag(veh0(i)) enddo open(unit=19,file=runname(1:lrunname)//'.flow') if(realp(1) /= 0.) then do i=1,net write(19,*) timo(i),realp(i)/realp(1) enddo close(19) endif call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave v_E',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c FSA uparp gmin=0.0 gmax=0.0 do i=1,net realp(i)=real(uparph0(i)) imagp(i)=aimag(uparph0(i)) enddo call mnmx(net,realp,gmin,gmax) call mnmx(net,imagp,gmin,gmax) call setzer call plchlq(.5,.95,'Flux Surf Ave u_par(eps/q)',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call dashdb(21845) call curve(timo,realp,net) call curved(timo,imagp,net) call finish(ipage, date) c plot potential energy vs. time (linear) gmin=0.0 gmax=0.0 do i=1,net n=nfr do m=1,md if(gmax < wtif(i,m,n)) gmax=wtif(i,m,n) if(gmin > wtif(i,m,n)) gmin=wtif(i,m,n) enddo enddo call setzer call plchlq(.5,.95,'Total Energy vs. mode',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) do m=1,md gamx(:)=wtif(:,m,nfr) call dashdb(dash(m)) call curved(timo,gamx,net) enddo call finish(ipage, date) c plot total energy vs. time (log) gmin=0.0 gmax=0.0 do i=1,net n=nfr do m=1,md gamx(i)=alog(abs(wtif(i,m,n))+1.e-30) if(gmax < gamx(i)) gmax=gamx(i) if(gmin > gamx(i)) gmin=gamx(i) enddo enddo call setzer call plchlq(.5,.95,'Total Energy (log_e) vs. mode',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) n=nfr do m=1,md do i=1,net gamx(i)=alog(abs(wtif(i,m,n))+1.e-30) enddo call dashdb(dash(m)) call curved(timo,gamx,net) enddo call finish(ipage, date) c c c plot kinetic energy vs. time (log) c gmin=0.0 gmax=0.0 do i=1,net n=nfr do m=1,md gamx(i)=alog(abs(wkif(i,m,n))+1.e-30) if(gmax < gamx(i)) gmax=gamx(i) if(gmin > gamx(i)) gmin=gamx(i) enddo enddo call setzer call plchlq(.5,.95,'Kinetic Energy (log_e) vs. mode',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) n=nfr do m=1,md do i=1,net gamx(i)=alog(abs(wkif(i,m,n))+1.e-30) enddo call dashdb(dash(m)) call curved(timo,gamx,net) enddo call finish(ipage, date) c c plot thermal energy vs. time (log) c gmin=0.0 gmax=0.0 do i=1,net n=nfr do m=1,md if(gmax < alog(wpif(i,m,n)+1.e-30)) . gmax=alog(wpif(i,m,n)+1.e-30) if(gmin > alog(wpif(i,m,n)+1.e-30)) . gmin=alog(wpif(i,m,n)+1.e-30) enddo enddo call setzer call plchlq(.5,.95,'Thermal Energy (log_e) vs. mode',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timo(1),timo(net),gmin,gmax) n=nfr do m=1,md gamx(:)=alog(wpif(:,m,n)+1.e-30) call dashdb(dash(m)) call curved(timo,gamx,net) enddo call finish(ipage, date) if(shr /= 0.0) then c c ballooning angle vs. m with circles, mbeer c can change to ovals by setting dx,dy c call setzer call plchlq(.5,.95,'Mode Amplitudes vs. m,n (phi)', > 20.,0.,0.) call plchlq(.5,.08,'ballooning coordinate',15.,0.,0.) call plchlq(.02,.5,'poloidal mode number',15.,90.,0.) gmin=mmin(1) gmax=mmax(1) do n=1,nd if (float(mmax(n)) > gmax) gmax=float(mmax(n)) if (float(mmin(n)) < gmin) gmin=float(mmin(n)) enddo gmax = gmax + 1. rmax=r(ld) c Choose instantaneous amplitude or growth rate for dy allocate(thing(md,nd)) if(lin == 1) then thing(:,:)=grtmx(net-1,:,:) else thing(:,:)=wtif(net-1,:,:) endif wmax=abs(thing(1,1)) do n=1,nd do m=1,md wmax=max(wmax,abs(thing(m,n))) rrat=r0(m,n) if(abs(rrat) > rmax) rmax = abs(rrat) enddo enddo call maps(-rmax,rmax,gmin,gmax) do n=1,nd do m=1,md rrat=r0(m,n) dy=thing(m,n)/abs(wmax) c c width isn't right. just make them pretty c dx=xp/10. xcir0=rrat+dx ycir0=float(mrr(m,n)) do l=0,10 thetac=float(l)/float(10)*6.2832 xcir1=rrat+dx*cos(thetac) ycir1=float(mrr(m,n))+dy*sin(thetac) if (dy > 0.) then call line(xcir0,ycir0,xcir1,ycir1) else call line(xcir1,ycir1,xcir1,ycir1) endif xcir0=xcir1 ycir0=ycir1 enddo enddo enddo call dashdb(21845) call lined(r(1),gmin,r(1),gmax) call lined(r(ld),gmin,r(ld),gmax) call line(-xp,gmin,-xp,gmax) call line(xp,gmin,xp,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, date) endif if(shr /= 0.0) then c c ballooning angle vs. m with circles, mbeer c can change to ovals by setting dx,dy c call setzer call plchlq(.5,.95,'Mode Amplitudes vs. m,n (pressure)', > 20.,0.,0.) call plchlq(.5,.08,'ballooning coordinate',15.,0.,0.) call plchlq(.02,.5,'poloidal mode number',15.,90.,0.) gmin=mmin(1) gmax=mmax(1) do n=1,nd if (float(mmax(n)) > gmax) gmax=float(mmax(n)) if (float(mmin(n)) < gmin) gmin=float(mmin(n)) enddo gmax = gmax + 1. rmax=r(ld) c Choose instantaneous amplitude or growth rate for dy do n=1,nd do m=1,md rrat=r0(m,n) if(abs(rrat) > rmax) rmax = abs(rrat) enddo enddo call maps(-rmax,rmax,gmin,gmax) do n=1,nd do m=1,md rrat=r0(m,n) dy=thing(m,n)/abs(wmax) c dx=wh(m,n) c c width isn't right. just make them pretty c dx=xp/10. xcir0=rrat+dx ycir0=float(mrr(m,n)) do l=0,10 thetac=float(l)/float(10)*6.2832 xcir1=rrat+dx*cos(thetac) ycir1=float(mrr(m,n))+dy*sin(thetac) if (dy > 0.) then call line(xcir0,ycir0,xcir1,ycir1) else call line(xcir1,ycir1,xcir1,ycir1) endif xcir0=xcir1 ycir0=ycir1 enddo enddo enddo 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, date) endif c c energy conservation if(icrit == 0) then gmin=0.0 gmax=0.0 call mnmx(net,pcerr,gmin,gmax) call setzer call plchlq(.5,.95,'Energy Conservation',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,pcerr,net-1) call finish(ipage, date) endif wt2 = wenx gmin=0.0 gmax=0.0 call mnmx(net,wt2,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Total Energy',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,wt2,net) call finish(ipage, date) call cplogp(wt2,wt1,net) gmin=-3.5 gmax=2.0 call mnmx(net,wt1,gmin,gmax) call setzer call plchlq(.5,.95,'Total 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 curve(timo,wt1,net) call finish(ipage, date) c c We have to be more careful with signs upon upgrading to multiple c ion species: c wt2 = 0. do j=1,nspecies wt2(:)=wt2(:)-wpfx(:,j) enddo wt2 = abs(wt2) call cplogp(wt2,wt1,net) gmin=0.0 gmax=0.0 call mnmx(net,wt1,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Thermal Flux',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, date) endif !end of igraph=1 section gmin=0.0 gmax=0.0 do n=1,net wt2(n)=0. enddo allocate(wt3(net,nspecies)) do i=1,nspecies wt2(:) = wt2(:) - wpfx(:,i) wt3(:,i) = -wpfx(:,i) call mnmx(net,wt3(:,i),gmin,gmax) enddo call mnmx(net,wt2,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Thermal Flux',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,wt2,net) do i=1,nspecies call dashdb(dash(i+1)) call curved(timo,wt3(1,i),net) enddo call finish(ipage, date) endif !end of igraph=1 section write(9,*) '---------' write(9,*) 'Flux vs. time: q_i,q_e,Gamma_e' do n=1,net write(9,136) timo(n),wt2(n),-qfluxe(n),-fluxe(n) enddo 136 format(4(e18.6,' ')) c c Particle fluxes c gmin=0.0 gmax=0.0 do n=1,net wt2(n)=0. enddo do i=1,nspecies wt2(:)=wt2(:)-fluxi(:,i) wt3(:,i)=-fluxi(:,i) call mnmx(net,wt3(:,i),gmin,gmax) enddo wt1 = -fluxe call mnmx(net,wt2,gmin,gmax) call mnmx(net,wt1,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Particle Fluxes',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,wt2,net) call dashdb(dash(2)) call curved(timo,wt1,net) do i=1,nspecies call dashdb(dash(i+2)) call curved(timo,wt3(1,i),net) enddo call finish(ipage, date) endif !end of igraph == 1 section c c Electron heat flux c gmin=0.0 gmax=0.0 do n=1,net wt2(n)=0. enddo do n=1,net wt2(n)=-qfluxe(n) enddo call mnmx(net,wt2,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Electron Heat Flux',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,wt2,net) call finish(ipage, date) c c do 1298 n=1,net c wt2(n)=-wt2(n) c 1298 continue c c Make a plot of the timestep vs. time c allocate(delta(ncnt)) m=1 do n=2,ncnt-1 delta(n)=timf(n+1)-timf(n) if(n >= nt_end(m)) m=m+1 delta(n)=delta(n)/float(ntimfp(m)) enddo delta(1)=delta(2) delta(ncnt)=delta(ncnt-1) gmin=0.0 gmax=1.1*dt0 call mnmx(ncnt,delta,gmin,gmax) call setzer call plchlq(.5,.95,'Time step',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call maps(timf(1),timf(ncnt),gmin,gmax) call curve(timf,delta,ncnt) call finish(ipage, date) endif ! end of igraph == 1 section allocate(wf2(net,md,nd)) wf2 = pi*y0*wtif/x0 allocate (aout(md), sdout(md), rmnum(md)) if(kyspec == 1) then call nsumavp(wf2,aout,sdout,net,fave) gmin=0.0 gmax=0.0 rmnum(:)=mrr(:,nd) call mnmx(md,aout,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin write(9,1676) (rmnum(m),aout(m),m=1,md) 1676 format(2x,3(' m=',f4.0,' ener=',e11.4)) c amplitude vs. mode number, mbeer if(igraph == 1) then 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 m=1,md call line(rmnum(m)-.3,0.0,rmnum(m)-.3,aout(m)) call line(rmnum(m)-.3,aout(m),rmnum(m)+.3,aout(m)) call line(rmnum(m)+.3,0.0,rmnum(m)+.3,aout(m)) enddo call finish(ipage, date) endif ! end of igraph == 1 section endif ! end of kyspec == 1 section if(igraph == 1) then c c Phi rms vs. time c gmin=0.0 gmax=0.0 call mnmx(net,phirms,gmin,gmax) call setzer call plchlq(.5,.95,'RMS Phi ',20. * ,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) c call plchlq(.02,.5,'',15.,90.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,phirms,net) call finish(ipage, date) xmin=gmin xmax=gmax gmin=0. gmax=0. wt1 = 0. do n=1,net do i=1,nspecies wt1(n)=wt1(n)-wpfx(n,i) enddo enddo call mnmx(net,wt1,gmin,gmax) call setzer call plchlq(.5,.95,'Flux vs. Phi ',20.,0.,0.) call plchlq(.5,.08,'Abs(Phi)',15.,0.,0.) c call plchlq(.02,.5,'',15.,90.,0.) call mapsm(xmin,xmax,gmin,gmax) call curve(phirms,wt1,net) call finish(ipage, date) c c Flux/Phi rms vs. time c gmin=0.0 gmax=0.0 allocate(fphi(net)) do n=1,net fphi(n)=0. do i=1,nspecies if(phirms(n) > 1.e-10) then fphi(n)=fphi(n)-wpfx(n,i)/phirms(n) else fphi(n)=0. endif enddo enddo call mnmx(net,fphi,gmin,gmax) call setzer call plchlq(.5,.95,'Flux/RMS Phi ',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) c call plchlq(.02,.5,'',15.,90.,0.) call mapsm(timo(1),timo(net),gmin,gmax) call curve(timo,fphi,net) call finish(ipage, date) c c new spectra c c instantaneous |phi(k_x)|^2 at z=0 (l=0) c c from phi* phi c c GWH: This has a BUG, because it doesn't average over positive and c negative k_x's, and the spectrum is not necessarily symmetric in k_x: c The same bug occurs in other k_x spectra plots also. This needs to c be fixed some day... ndpmax=(nd-1)/2+1 ! n corresponding to maximum positive k_x ndpmaxerr=(nd-1)/2 ! ?? reproduce previous erroneous value c ?? To remove this error later, replace all ndpmaxerr with ndpmax kxmax=abs(shr)*z0/y0 allocate (rnnum(nd)) if (nd /= 1) then gmin=0.0 gmax=0.0 do n=1,ndpmax rnnum(n)=float(n-1)/float((nd-1)/2)*z0*abs(shr)/y0 wt1(n)=psp(net,1,n) do m=2,md wt1(n)=wt1(n)+0.5*psp(net,m,n) enddo enddo call mnmx(ndpmaxerr,wt1,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Instantaneous |Phi|^2(k_x)' & ,20.,0.,0.) call plchlq(.02,.5,'|Phi|^2',15.,90.,0.) call plchlq(.5,.08,'k_x rho_i',15.,0.,0.) call mapsm(0.,kxmax,gmin,gmax) barwid=kxmax/3./float((nd-1)/2) do i=1,ndpmaxerr call line(rnnum(i)-barwid,0.0,rnnum(i)-barwid,wt1(i)) call line(rnnum(i)-barwid,wt1(i),rnnum(i)+barwid,wt1(i)) call line(rnnum(i)+barwid,0.0,rnnum(i)+barwid,wt1(i)) enddo call finish(ipage, date) endif c c instantaneous |phi(k_y)|^2 at z=0 (l=ldb/2) c c from phi* phi c if (md /= 1) then gmin=0.0 gmax=0.0 do m=1,md rmnum(m)=float(m-1)/y0 aout(m)=0. do n=1,nd aout(m)=aout(m)+psp(net,m,n) enddo if (m > 1) aout(m)=aout(m)/4. enddo call mnmx(md,aout,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin call setzer call plchlq(.5,.95,'Instantaneous |Phi|^2(k_y)' & ,20.,0.,0.) call plchlq(.02,.5,'|Phi|^2',15.,90.,0.) call plchlq(.5,.08,'k_y rho_i',15.,0.,0.) call mapsm(0.,float(md)/y0,gmin,gmax) barwid=float(md)/y0/3./float(md) do m=1,md call line(rmnum(m)-barwid,0.0,rmnum(m)-barwid,aout(m)) call line(rmnum(m)-barwid,aout(m),rmnum(m)+barwid,aout(m)) call line(rmnum(m)+barwid,0.0,rmnum(m)+barwid,aout(m)) enddo call finish(ipage, date) endif endif ! end of igraph == 1 section c c n_e spectra: this isn't quite n_e, but close, just phi w/o m=0 c c time averaged |n_e(k_x)|^2 at z=0 (l=ldb/2) c if (nd /= 1) then gmin=0.0 gmax=0.0 do n=1,ndpmax wt1(n)=0. do i=1,net wt2(i)=psp(i,1,n) enddo call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) c wt1(n)=wt1(n)+ave do m=2,md do i=1,net wt2(i)=0.5*psp(i,m,n) enddo call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) wt1(n)=wt1(n)+ave enddo enddo call mnmx(ndpmaxerr,wt1,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Time Averaged |n_e|^2(k_x)' & ,20.,0.,0.) call plchlq(.02,.5,'|n_e|^2',15.,90.,0.) call plchlq(.5,.08,'k_x rho_i',15.,0.,0.) call mapsm(0.,kxmax,gmin,gmax) barwid=kxmax/3./float((nd-1)/2) do n=1,ndpmaxerr call line(rnnum(n)-barwid,0.0,rnnum(n)-barwid,wt1(n)) call line(rnnum(n)-barwid,wt1(n),rnnum(n)+barwid,wt1(n)) call line(rnnum(n)+barwid,0.0,rnnum(n)+barwid,wt1(n)) enddo call finish(ipage, date) endif ! end of igraph == 1 section write(9,*) '----- Time averaged ne x-spectrum: ' do n=1,ndpmaxerr write(9,*) n,rnnum(n),wt1(n) enddo endif c c time averaged |n_e(k_y)|^2 at z=0 (l=ldb/2) c if (md /= 1) then gmin=0.0 gmax=0.0 aout(1)=0. do m=2,md aout(m)=0. do n=1,nd do i=1,net wt2(i)=psp(i,m,n) enddo call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) aout(m)=aout(m)+ave enddo if (m > 1) aout(m)=aout(m)/4. enddo call mnmx(md,aout,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Time Averaged |n_e|^2(k_y)' & ,20.,0.,0.) call plchlq(.02,.5,'|n_e|^2',15.,90.,0.) call plchlq(.5,.08,'k_y rho_i',15.,0.,0.) call mapsm(0.,float(md)/y0,gmin,gmax) barwid=float(md)/y0/3./float(md) do m=1,md call line(rmnum(m)-barwid,0.0,rmnum(m)-barwid,aout(m)) call line(rmnum(m)-barwid,aout(m),rmnum(m)+barwid,aout(m)) call line(rmnum(m)+barwid,0.0,rmnum(m)+barwid,aout(m)) enddo call finish(ipage, date) endif !end of igraph == 1 section write(9,*) '---- Time averaged ne y-spectrum' do m=1,md write(9,*) m,rmnum(m),aout(m) enddo endif c c new spectra c c time averaged |phi(k_x)|^2 at z=0 (l=ldb/2) c if (nd /= 1) then gmin=0.0 gmax=0.0 do n=1,ndpmax wt1(n)=0. do i=1,net wt2(i)=psp(i,1,n) enddo call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) wt1(n)=wt1(n)+ave do m=2,md do i=1,net wt2(i)=0.5*psp(i,m,n) enddo call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) wt1(n)=wt1(n)+ave enddo enddo call mnmx(ndpmaxerr,wt1,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Time Averaged |Phi|^2(k_x)' & ,20.,0.,0.) call plchlq(.02,.5,'|Phi|^2',15.,90.,0.) call plchlq(.5,.08,'k_x rho_i',15.,0.,0.) call mapsm(0.,kxmax,gmin,gmax) barwid=kxmax/3./float((nd-1)/2) do n=1,ndpmaxerr call line(rnnum(n)-barwid,0.0,rnnum(n)-barwid,wt1(n)) call line(rnnum(n)-barwid,wt1(n),rnnum(n)+barwid,wt1(n)) call line(rnnum(n)+barwid,0.0,rnnum(n)+barwid,wt1(n)) enddo call finish(ipage, date) endif !end of igraph == 1 section write(9,*) '----- Time averaged x-spectrum: ' do n=1,ndpmaxerr write(9,*) n,rnnum(n),wt1(n) enddo endif c c time averaged |phi(k_y)|^2 at z=0 (l=ldb/2) c if (md /= 1) then gmin=0.0 gmax=0.0 do m=1,md aout(m)=0. do n=1,nd do i=1,net wt2(i)=psp(i,m,n) enddo call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) aout(m)=aout(m)+ave enddo if (m > 1) aout(m)=aout(m)/4. enddo call mnmx(md,aout,gmin,gmax) gmax=1.1*gmax gmin=1.1*gmin if(igraph == 1) then call setzer call plchlq(.5,.95,'Time Averaged |Phi|^2(k_y)' & ,20.,0.,0.) call plchlq(.02,.5,'|Phi|^2',15.,90.,0.) call plchlq(.5,.08,'k_y rho_i',15.,0.,0.) call mapsm(0.,float(md)/y0,gmin,gmax) barwid=float(md)/y0/3./float(md) do m=1,md call line(rmnum(m)-barwid,0.0,rmnum(m)-barwid,aout(m)) call line(rmnum(m)-barwid,aout(m),rmnum(m)+barwid,aout(m)) call line(rmnum(m)+barwid,0.0,rmnum(m)+barwid,aout(m)) enddo call finish(ipage, date) endif ! end of igraph == 1 section write(9,*) '---- Time averaged y-spectrum' do m=1,md write(9,*) m,rmnum(m),aout(m) enddo endif allocate(xxx(ldb), yyy(ldb), work(ldb)) allocate (cxpr(nalias), cxpi(nalias)) allocate (cypr(2*malias), cypi(2*malias), worky(2*malias)) allocate (xx(nalias), yx(nalias), workx(nalias), cxa(nalias), cza(ld)) allocate(xy(2*malias), yy(2*malias), cya(2*malias)) allocate (cxb(4*nd), cyb(4*md), czb(ld)) allocate (cxc(4*nd), cyc(4*md), czf(ld), czd(ld), cze(ld)) skip_correlations=.true. if(skip_correlations) goto 175 !ignore correlation functions for now c c y-averaged correlation function C(x) c if (lin == 0) then npts=(nd-1)/2 xx(1)=wt1(1) yx(1)=0. do n=2,npts xx(n)=wt1(n) yx(n)=0. xx(nalias-n+2)=wt1(n) yx(nalias-n+2)=0. enddo xx(npts+1)=wt1(npts+1) yx(npts+1)=0. do n=npts+2,nalias-npts+1 xx(n)=0. yx(n)=0. enddo do n=1,nalias/2 cxa(n)=pi*float(n-1)/float(nalias/2-1)* & y0/(abs(shr)*z0)*float((nd-1)/2) enddo call c06fce(xx,yx,nalias,workx,ifail) ! fft do n=1,nalias cxpr(n)=xx(n)/xx(1) cxpi(n)=yx(n)/xx(1) enddo xmin=0. xmax=0. call mnmx(nalias/2,cxpr,xmin,xmax) call mnmx(nalias/2,cxpi,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'C_x, averaged in y ',20. * ,0.,0.) call plchlq(.5,.08,'x/rho_i',15.,0.,0.) call mapsm(cxa(1),cxa(nalias/2),xmin,xmax) call curve(cxa,cxpr,nalias/2) c call dashdb(21845) c call curved(cxa,cxpi,npts) call finish(ipage, date) endif !end of igraph == 1 section c c x-averaged correlation function C(y) c mpts=(md-1)*2 xy(1)=aout(1) yy(1)=0. do m=2,md xy(m)=aout(m) yy(m)=0. xy(2*malias-m+2)=aout(m) yy(2*malias-m+2)=0. enddo do m=md+1,2*malias-md+1 xy(m)=0. yy(m)=0. enddo do m=1,malias cya(m)=pi*float(m-1)/float(malias-1)*y0 enddo call c06fce(xy,yy,2*malias,worky,ifail) ! fft do m=1,malias*2 cypr(m)=xy(m)/xy(1) cypi(m)=yy(m)/xy(1) enddo xmin=0. xmax=0. call mnmx(malias,cypr,xmin,xmax) call mnmx(malias,cypi,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'C_y, averaged in x',20. * ,0.,0.) call plchlq(.5,.08,'y/rho_i',15.,0.,0.) call mapsm(cya(1),cya(malias),xmin,xmax) call curve(cya,cypr,malias) c call dashdb(21845) c call curved(cya,cypi,md) call finish(ipage, date) endif !end of igraph == 1 section c c Correlation functions c c c?a=?-axis, czd=z-axis for C(-z,z) c mpts=4*md npts=4*nd do i=1,npts c cxa(i)=2.0*pi*float(i-1-npts/2)/float(npts-1)* c & y0/(abs(shr)*z0)*float((nd-1)/2) cxb(i)=0. cxc(i)=0. enddo do i=1,mpts cya(i)=2.0*pi*float(i-1-mpts/2)/float(mpts-1)*y0 cyb(i)=0. cyc(i)=0. enddo do i=1,ld cza(i)=r(i) czb(i)=0. czf(i)=0. enddo do i=1,ldb/2 czd(i)=r(i+ldb/2) cze(i)=0. enddo c c The non-averaged correlation functions are not so interesting, c so skip this usually part. c skip_non_avg_correlations=.true. if(skip_non_avg_correlations) goto 34 c c subtract off means c c c?m=mean of phi(?) c do i=1,net cxm=0.0 do j=1,npts cxm=cxm+cx(i,j) enddo do j=1,npts cx(i,j)=cx(i,j)-cxm/npts enddo cym=0.0 do j=1,mpts cym=cym+cy(i,j) enddo do j=1,mpts cy(i,j)=cy(i,j)-cym/mpts enddo czm=0.0 do j=1,ld-1 czm=czm+cz(i,j) enddo do j=1,ld-1 cz(i,j)=cz(i,j)-czm/(ld-1) enddo enddo c c time average c c c?0=phi(0), c?c(i)=, c?b(i)= c cze= c do i=1,net cx0=cx(i,npts/2) do j=1,npts cxb(j)=cxb(j)+cx(i,j)*cx0 cxc(j)=cxc(j)+cx(i,j)*cx(i,j) enddo cy0=cy(i,mpts/2) do j=1,mpts cyb(j)=cyb(j)+cy(i,j)*cy0 cyc(j)=cyc(j)+cy(i,j)*cy(i,j) enddo cz0=cz(i,ldb/2+1) do j=1,ld-1 czb(j)=czb(j)+cz(i,j)*cz0 czf(j)=czf(j)+cz(i,j)*cz(i,j) enddo do j=1,ldb/2 cze(j)=cze(j)+cz(i,j+ldb/2)*cz(i,ldb/2+2-j) enddo enddo cxmax=0.0 cxmin=0.0 cymax=0.0 cymin=0.0 czmax=0.0 czmin=0.0 c c c?0= c cx0=cxb(npts/2) do j=1,npts cxb(j)=cxb(j)/sqrt(cx0)/sqrt(cxc(j)) if (cxb(j) > cxmax) cxmax=cxb(j) if (cxb(j) < cxmin) cxmin=cxb(j) enddo cy0=cyb(mpts/2) do j=1,mpts cyb(j)=cyb(j)/sqrt(cy0)/sqrt(cyc(j)) if (cyb(j) > cymax) cymax=cyb(j) if (cyb(j) < cymin) cymin=cyb(j) enddo cz0=czb(ldb/2+1) write(9,*) '---- correlation function C(z,0)' do j=1,ld-1 czb(j)=czb(j)/sqrt(cz0)/sqrt(czf(j)) if (czb(j) > czmax) czmax=czb(j) if (czb(j) < czmin) czmin=czb(j) write(9,*) cza(j),czb(j) enddo write(9,*) '---- correlation function C(z,-z)' do j=1,ldb/2 cze(j)=cze(j)/sqrt(czf(j+ldb/2))/sqrt(czf(ldb/2+2-j)) if (cze(j) > czemax) czemax=cze(j) if (cze(j) < czemin) czemin=cze(j) write(9,*) czd(j),cze(j) enddo write(9,*) '---- Phi^2 along field line' do j=1,ld-1 czf(j)=sqrt(czf(j)) if (czf(j) > czcmax) czcmax=czf(j) if (czf(j) < czcmin) czcmin=czf(j) write(9,*) cza(j),czf(j) enddo if(igraph == 1) then call setzer call plchlq(.5,.95,'C_x ',20. * ,0.,0.) call plchlq(.5,.08,'x',15.,0.,0.) call mapsm(cxa(1),cxa(npts),cxmin,cxmax) call curve(cxa,cxb,npts) call finish(ipage, date) call setzer call plchlq(.5,.95,'C_y ',20. * ,0.,0.) call plchlq(.5,.08,'y',15.,0.,0.) call mapsm(cya(1),cya(mpts),cymin,cymax) call curve(cya,cyb,mpts) call finish(ipage, date) call setzer call plchlq(.5,.95,'C_z ',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czmin,czmax) call curve(cza,czb,ld-1) call finish(ipage, date) call plchlq(.5,.95,'C(-z,z) ',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(czd(1),czd(ldb/2),czemin,czemax) call curve(czd,cze,ldb/2) call finish(ipage, date) call setzer call plchlq(.5,.95,'SQRT()',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czcmin,czcmax) call curve(cza,czf,ld-1) call finish(ipage, date) endif ! end of igraph == 1 section 34 continue c c parallel correlation function averaged in x and y, potential c czb = 0. cze = 0. do j=1,ld wt2(:)=czc(:,j) call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) czb(j)=ave wt2(:)=czcn(:,j) call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) cze(j)=ave enddo write(9,*) '---- Phi^2 along field line, averaged in x,y' write(9,*) ' and phi(z)phi(0)/sqrt(phi^2(0) phi^2(z))' write(9,*) ' and phi(z)phi(0)/sqrt(phi^2(0) phi^2(0))' do j=1,ldb czd(j)=czb(j)/cze(ldb/2+1) czb(j)=czb(j)/sqrt(cze(j))/sqrt(cze(ldb/2+1)) write(9,*) cza(j),sqrt(cze(j)),czb(j),czd(j) enddo czemin=0. czemax=0. czcmin=0. czcmax=0. czdmin=0. czdmax=0. do j=1,ldb cze(j)=sqrt(cze(j)) if (cze(j) > czemax) czemax=cze(j) if (cze(j) < czemin) czemin=cze(j) if (czb(j) > czcmax) czcmax=czb(j) if (czb(j) < czcmin) czcmin=czb(j) if (czd(j) > czdmax) czdmax=czd(j) if (czd(j) < czdmin) czdmin=czd(j) enddo if(igraph == 1) then call setzer call plchlq(.5,.95,'SQRT(), averaged in x,y, phi',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czemin,czemax) call curve(cza,cze,ld-1) call finish(ipage, date) call setzer call plchlq(.5,.95,'C_z, norm. to (0,z), ave. in x,y, phi',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czcmin,czcmax) call curve(cza,czb,ld-1) call finish(ipage, date) call plchlq(.5,.95,'C_z, norm. to (0,0), ave. in x,y, phi',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czdmin,czdmax) call curve(cza,czd,ld-1) call finish(ipage, date) endif ! end of igraph == 1 section c c x and y-averaged k_parallel spectrum, potential c do n=1,ldb/2 xxx(n)=czd(n+ldb/2) yyy(n)=0. xxx(ldb-n+1)=czd(ldb/2+1-n) yyy(ldb-n+1)=0. enddo do n=1,ldb/2 cza(n)=float(n-1)/float(ldb/2-1)*(ldb/2)*pi/x0 enddo call c06fce(xxx,yyy,ldb,work,ifail) ! fft do n=1,ldb/2 czd(n)=xxx(n)/xxx(1) cze(n)=yyy(n)/xxx(1) enddo xmin=0. xmax=0. call mnmx(ldb/2,czd,xmin,xmax) call mnmx(ldb/2,cze,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'k_par spectrum of phi, averaged in x,y ',20. * ,0.,0.) call plchlq(.5,.08,'k_par q R',15.,0.,0.) call mapsm(cza(1),cza(ldb/2),xmin,xmax) call curve(cza,czd,ldb/2) call dashdb(21845) call curved(cza,cze,ldb/2) call finish(ipage, date) endif ! end of igraph == 1 section c c parallel correlation function averaged in x and y, electron density c cza = r czb = 0. cze = 0. do j=1,ld wt2(:)=czn(:,j) call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) czb(j)=ave wt2(:)=cznn(:,j) call timeavp(wt2,ave,sdv,net,fave,pfmin,pfmax) cze(j)=ave enddo write(9,*) '---- N_e^2 along field line, averaged in x,y' write(9,*) ' and N_e(z)N_e(0)/sqrt(N_e^2(0) N_e^2(z))' write(9,*) ' and N_e(z)N_e(0)/sqrt(N_e^2(0) N_e^2(0))' do j=1,ldb czd(j)=czb(j)/cze(ldb/2+1) czb(j)=czb(j)/sqrt(cze(j))/sqrt(cze(ldb/2+1)) write(9,*) cza(j),sqrt(cze(j)),czb(j),czd(j) enddo czemin=0. czemax=0. czcmin=0. czcmax=0. czdmin=0. czdmax=0. do j=1,ldb cze(j)=sqrt(cze(j)) if (cze(j) > czemax) czemax=cze(j) if (cze(j) < czemin) czemin=cze(j) if (czb(j) > czcmax) czcmax=czb(j) if (czb(j) < czcmin) czcmin=czb(j) if (czd(j) > czdmax) czdmax=czd(j) if (czd(j) < czdmin) czdmin=czd(j) enddo if(igraph == 1) then call setzer call plchlq(.5,.95,'SQRT(), averaged in x,y',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czemin,czemax) call curve(cza,cze,ld-1) call finish(ipage, date) call setzer call plchlq(.5,.95,'C_z, norm. to (0,z), ave. in x,y, N_e',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czcmin,czcmax) call curve(cza,czb,ld-1) call finish(ipage, date) call plchlq(.5,.95,'C_z, norm. to (0,0), ave. in x,y, N_e',20. * ,0.,0.) call plchlq(.5,.08,'z',15.,0.,0.) call mapsm(cza(1),cza(ld-1),czdmin,czdmax) call curve(cza,czd,ld-1) call finish(ipage, date) endif ! end of igraph == 1 section c c x and y-averaged k_parallel spectrum, electron density c do n=1,ldb/2 xxx(n)=czd(n+ldb/2) yyy(n)=0. xxx(ldb-n+1)=czd(ldb/2+1-n) yyy(ldb-n+1)=0. enddo do n=1,ldb/2 cza(n)=float(n-1)/float(ldb/2-1)*(ldb/2)*pi/x0 enddo call c06fce(xxx,yyy,ldb,work,ifail) ! fft do n=1,ldb/2 czd(n)=xxx(n)/xxx(1) cze(n)=yyy(n)/xxx(1) enddo xmin=0. xmax=0. call mnmx(ldb/2,czd,xmin,xmax) call mnmx(ldb/2,cze,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'k_par spectrum of N_e, averaged in x,y ',20. * ,0.,0.) call plchlq(.5,.08,'k_par q R',15.,0.,0.) call mapsm(cza(1),cza(ldb/2),xmin,xmax) call curve(cza,czd,ldb/2) call dashdb(21845) call curved(cza,cze,ldb/2) call finish(ipage, date) endif ! end of igraph == 1 section endif c c time correlation function, not yet averaged over restarts c xmin=0. xmax=0. allocate (ctt(net)) do i=1,net if (mod(i-1,ninterv) == 0) norm=ct(i) ct(i)=ct(i)/norm if (ct(i) > xmax) xmax=ct(i) if (ct(i) < xmin) xmin=ct(i) ctt(i)=timo(i)-timo(1) enddo if(igraph == 1 .and. icrit == 0) then call setzer call plchlq(.5,.95,'C(t,z=0), averaged in x,y ',20. * ,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(ctt(1),ctt(net),xmin,xmax) call curve(ctt,ct,net) call finish(ipage, date) endif 175 continue c c Flux surface averaged potential c open(unit=35,file=runname(1:lrunname)//'.mov') write(35,*) net,4*nd do i=1,net npts=4*nd do j=1,npts cxb(j)=2.0*pi*float(j-1-npts/2)/float(npts-1)* & y0/(abs(shr)*z0)*float((nd-1)/2) x1=2.0*pi*float(j-1-npts/2)/float(npts-1) cxc(j)=0.0 do n=1,nd nr=nrr(n) cxc(j)=cxc(j)+real(phi00(i,n)*exp(ii*nr*x1)) enddo enddo write(35,*) (cxc(j),j=1,4*nd) enddo close(35) xmin=0. xmax=0. call mnmx(4*nd,cxc,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'Flux surface averaged Potential',20. * ,0.,0.) call plchlq(.5,.08,'x/rho_i',15.,0.,0.) call mapsm(cxb(1),cxb(4*nd),xmin,xmax) call curve(cxb,cxc,4*nd) call finish(ipage, date) endif c c Flux surface averaged parallel velocity c c open(unit=35,file=runname(1:lrunname)//'.mov') c write(35,*) net,4*nd do i=1,net npts=4*nd do j=1,npts cxb(j)=2.0*pi*float(j-1-npts/2)/float(npts-1)* & y0/(abs(shr)*z0)*float((nd-1)/2) x1=2.0*pi*float(j-1-npts/2)/float(npts-1) cxc(j)=0.0 do n=1,nd nr=nrr(n) cxc(j)=cxc(j)+real(upar00(i,n)*exp(ii*nr*x1)) enddo enddo c write(35,*) (cxc(j),j=1,4*nd) enddo c close(35) xmin=0. xmax=0. call mnmx(4*nd,cxc,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'Flux surface averaged parallel flow',20. * ,0.,0.) call plchlq(.5,.08,'x/rho_i',15.,0.,0.) call mapsm(cxb(1),cxb(4*nd),xmin,xmax) call curve(cxb,cxc,4*nd) call finish(ipage, date) endif c c Flux surface averaged total pressure (no FLR corrections) c i=net npts=4*nd do j=1,npts cxb(j)=2.0*pi*float(j-1-npts/2)/float(npts-1)* & y0/(abs(shr)*z0)*float((nd-1)/2) x1=2.0*pi*float(j-1-npts/2)/float(npts-1) cxc(j)=0.0 do n=1,nd nr=nrr(n) cxc(j)=cxc(j)+real((1.5*den00(i,n)+tpar00(i,n)/2+tperp00(i,n)) & *exp(ii*nr*x1)) enddo enddo xmin=0. xmax=0. call mnmx(4*nd,cxc,xmin,xmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'Flux surface averaged total pressure',20. * ,0.,0.) call plchlq(.5,.08,'x/rho_i',15.,0.,0.) call mapsm(cxb(1),cxb(4*nd),xmin,xmax) call curve(cxb,cxc,4*nd) call finish(ipage, date) endif c write(*,*) 'x/rho, ' c do i=1,4*nd c write(*,173) cxb(i),cxc(i) c enddo c c drive terms c xmin=0. xmax=0.0 do n=1,14 do i=1,net if (drt(n,i) > xmax) xmax=drt(n,i) if (drt(n,i) < xmin) xmin=drt(n,i) enddo enddo if(plot_drives == 1 .and. igraph == 1) then call setzer call plchlq(.5,.95,'Driving terms ',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),xmin,xmax) do n=1,14 do i=1,net wt1(i)=drt(n,i) enddo if (mod(n,2) == 0) then call curve(timo,wt1,net) else call dashdb(21845) call curved(timo,wt1,net) endif enddo xxmin=timo(1) xxmax=timo(net) yymin=xmin yymax=xmax if (yymin == yymax) yymax=yymax+1.0 nxdiv=2 nydiv=2 call nicer(xxmin,xxmax,nox,noxm,fx,nx) call nicer(yymin,yymax,noy,noym,fy,ny) xxmax=xxmax+(xxmax-xxmin)*.05 call set(0.2,0.95,0.15,0.9,xxmin,xxmax,yymin,yymax,1) c call set(0.2,0.95,0.15,0.9,timo(1),timo(net)/.95,xmin,xmax,1) do n=1,14 write(str,133) n c BD: array bound error here? call plchlq(timo(net),drt(n,net),str,12.,0.,-1.) enddo 133 format(i3) call finish(ipage, date) endif ! end of igraph == 1 section xmin=0. xmax=0. allocate(dreta(net),drtor(net),wupsi(net), . dapar(net),dator(net),dtot(net)) do i=1,net dreta(i)=drt(1,i)+drt(2,i)+drt(3,i) drtor(i)=drt(6,i)+drt(8,i)+drt(11,i) wupsi(i)=drt(14,i) dapar(i)=drt(4,i)+drt(5,i) dator(i)=drt(7,i)+drt(9,i)+drt(10,i)+drt(12,i)+drt(13,i) dtot(i)=dreta(i)+drtor(i)+wupsi(i)+dapar(i)+dator(i) enddo if(plot_drives == 1 .and. igraph == 1) then call mnmx(net,dreta,xmin,xmax) call mnmx(net,wupsi,xmin,xmax) call mnmx(net,drtor,xmin,xmax) call mnmx(net,dapar,xmin,xmax) call mnmx(net,dator,xmin,xmax) call mnmx(net,dtot,xmin,xmax) call setzer call plchlq(.5,.95,'Driving terms ',20.,0.,0.) call plchlq(.5,.08,'time',15.,0.,0.) call mapsm(timo(1),timo(net),xmin,xmax) call curve(timo,dreta,net) call curve(timo,drtor,net) call curve(timo,dtot,net) call curve(timo,wupsi,net) call curve(timo,dapar,net) call curve(timo,dator,net) xxmin=timo(1) xxmax=timo(net) yymin=xmin yymax=xmax if (yymin == yymax) yymax=yymax+1.0 nxdiv=2 nydiv=2 call nicer(xxmin,xxmax,nox,noxm,fx,nx) call nicer(yymin,yymax,noy,noym,fy,ny) xxmax=xxmax+(xxmax-xxmin)*.05 call set(0.2,0.95,0.15,0.9,xxmin,xxmax,yymin,yymax,1) 134 format(a5) write(str,134) 'Total' call plchlq(timo(net),dtot(net),str,10.,0.,-1.) write(str,134) 'Wupsi' call plchlq(timo(net),wupsi(net),str,10.,0.,-1.) write(str,134) 'DrEta' call plchlq(timo(net),dreta(net),str,10.,0.,-1.) write(str,134) 'DrTor' call plchlq(timo(net),drtor(net),str,10.,0.,-1.) write(str,134) 'DaPar' call plchlq(timo(net),dapar(net),str,10.,0.,-1.) write(str,134) 'DaTor' call plchlq(timo(net),dator(net),str,10.,0.,-1.) call finish(ipage, date) endif ! end of igraph == 1 section write(9,*) '--- Time, D_*, D_d+, D_d-, D_\|, W_u\Psi, T, pcerr' do i=1,net write(9,*) timo(i),dreta(i),drtor(i),dator(i),dapar(i), & wupsi(i),dtot(i),pcerr(i) enddo c c k-contour plots c if(plot_drives == 1 .and. igraph == 1) then call cnplotk(psp, '|Phi|^2 ',tim,ipage,net,fave) call cnplotk(wenrk,'Energy ',tim,ipage,net,fave) call cnplotk(dketa,'Slab drive ',tim,ipage,net,fave) call cnplotk(dkpar,'Parallel Damping',tim,ipage,net,fave) call cnplotk(dktor,'Toroidal drive ',tim,ipage,net,fave) call cnplotk(dktdp,'Toroidal damping',tim,ipage,net,fave) call cnplotk(wkups,'W_(u,Psi) ',tim,ipage,net,fave) call cnplotk(dktot,'Total drive ',tim,ipage,net,fave) endif c c linear stuff c if (lin == 1) then c growth rates vs. mode number, mbeer gmin=0. gmax=0. c do m=1,md c if (mgamx(m,1,net) > gmax) gmax=mgamx(m,1,net) c if (mgamx(m,1,net) < gmin) gmin=mgamx(m,1,net) c enddo call mnmx(md,mgamx(:,1,net),gmin,gmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'Growth rate spectrum',20.,0.,0.) call plchlq(.02,.5,'growth rate (v_ti/L_n)',15.,90.,0.) call plchlq(.5,.08,'k_theta rho_i',15.,0.,0.) c took out m=0 mode. dxtot=float(mrr(md,1))/y0 if(dxtot == 0.0) dxtot=1.0/y0 call mapsm(0.,dxtot,gmin,gmax) dy=(gmax-gmin)/50. dx=(float(mrr(md,1))/y0)/50. do m=1,md rmnum(m)=mrr(m,nd) c gmin=(rmnum(m)-.3)/y0 c gmax=(rmnum(m)+.3)/y0 gmin=rmnum(m)/y0+dx gmax=rmnum(m)/y0-dx call line(gmin,mgamx(m,1,net)+dy,gmin,mgamx(m,1,net)-dy) call line(gmax,mgamx(m,1,net)+dy,gmax,mgamx(m,1,net)-dy) call line(gmin,mgamx(m,1,net)+dy,gmax,mgamx(m,1,net)+dy) call line(gmin,mgamx(m,1,net)-dy,gmax,mgamx(m,1,net)-dy) c call line(gmin,0.0,gmin,mgamx(m,net)) c call line(gmin,aveth(m),gmax,mgamx(m,net)) c call line(gmax,0.0,gmax,mgamx(m,net)) c write(*,*) rmnum(m)/y0,mgamx(m,1,net) enddo call finish(ipage, date) endif ! end of igraph == 1 section c mode width vs. mode number, mbeer gmin=0. gmax=0. do m=1,md aveth(m)=0.0 avekp2(m)=0.0 norm=0.0 do l=1,ld aveth(m)=aveth(m)+(cabs(density(l,m,1,1))**2)*r(l)**2 avekp2(m)=avekp2(m)+(cabs(density(l,m,1,1))**2) > *rkperp2(l,m,1) norm=norm+(cabs(density(l,m,1,1))**2) c write(*,*) 'rkperp2',l,m,rkperp2(l,m,1) enddo if (norm /= 0.0) then aveth(m)=sqrt(aveth(m)/norm) avekp2(m)=avekp2(m)/norm else aveth(m)=0.0 avekp2(m)=(float(mrr(m,1))/y0)**2 endif enddo call mnmx(md,aveth,gmin,gmax) if(igraph == 1) then call setzer call plchlq(.5,.95,'Mode width spectrum',20.,0.,0.) call plchlq(.02,.5,'mode width ',15.,90.,0.) call plchlq(.5,.08,'k_theta rho_i',15.,0.,0.) c took out m=0 mode. call mapsm(0.,dxtot,gmin,gmax) dy=(gmax-gmin)/50. dx=(float(mrr(md,1))/y0)/50. do m=1,md c gmin=(rmnum(m)-.3)/y0 c gmax=(rmnum(m)+.3)/y0 gmin=rmnum(m)/y0+dx gmax=rmnum(m)/y0-dx call line(gmin,aveth(m)+dy,gmin,aveth(m)-dy) call line(gmax,aveth(m)+dy,gmax,aveth(m)-dy) call line(gmin,aveth(m)+dy,gmax,aveth(m)+dy) call line(gmin,aveth(m)-dy,gmax,aveth(m)-dy) c call line(gmin,0.0,gmin,aveth(m)) c call line(gmin,aveth(m),gmax,aveth(m)) c call line(gmax,0.0,gmax,aveth(m)) enddo call finish(ipage, date) endif ! end of igraph == 1 section c gamma delta_x **2 vs. mode number gmin=0. gmax=0. do m=1,md c BD: c I haven't defined carefully for igeo=1. c if ((mrr(m,1) /= 0).and.(aveth(m) /= 0.)) then gdx=mgamx(m,1,net)/(float(mrr(m,1))/y0*shr*aveth(m))**2 else gdx=0.0 endif gmax=max(gmax,gdx) gmin=min(gmin,gdx) if (gmin < 0.) gmin=0. enddo if(igraph == 1) then call setzer call plchlq(.5,.95,'gamma del_x**2 spectrum',20.,0.,0.) call plchlq(.02,.5,'gamma del_x**2',15.,90.,0.) call plchlq(.5,.08,'k_theta rho_i',15.,0.,0.) c took out m=0 mode. call mapsm(0.,dxtot,gmin,gmax) dy=(gmax-gmin)/50. dx=(float(mrr(md,1))/y0)/50. do m=1,md gmin=rmnum(m)/y0+dx gmax=rmnum(m)/y0-dx c gmin=(rmnum(m)-.3)/y0 c gmax=(rmnum(m)+.3)/y0 if ((mrr(m,1) /= 0).and.(aveth(m) /= 0)) then gdx=mgamx(m,1,net)/(float(mrr(m,1))/y0*shr*aveth(m))**2 else gdx=0.0 endif call line(gmin,gdx+dy,gmin,gdx-dy) call line(gmax,gdx+dy,gmax,gdx-dy) call line(gmin,gdx+dy,gmax,gdx+dy) call line(gmin,gdx-dy,gmax,gdx-dy) c call line(gmin,0.0,gmin,gdx) c call line(gmin,gdx,gmax,gdx) c call line(gmax,0.0,gmax,gdx) enddo call finish(ipage, date) endif ! end of igraph == 1 section c gamma delta **2 vs. mode number gmin=0. gmax=0. do m=1,md if(avekp2(m) /= 0) then c if (mrr(m,1) /= 0) then gdx=mgamx(m,1,net)/avekp2(m) c gdx=mgamx(m,1,net)/(float(mrr(m,1))/y0)**2/ c & (1.+(shr*aveth(m))**2) else c write(*,*) '< k_perp**2 > = 0 for m= ',mrr(m,1) gdx=0. endif gmax=max(gmax,gdx) gmin=min(gmin,gdx) if (gmin < 0.) gmin=0. enddo if(igraph == 1) then call setzer call plchlq(.5,.95,'gamma delta**2 spectrum' & ,20.,0.,0.) call plchlq(.02,.5,'gamma delta**2',15.,90.,0.) call plchlq(.5,.08,'k_theta rho_i',15.,0.,0.) c took out m=0 mode. call mapsm(0.,dxtot,gmin,gmax) dy=(gmax-gmin)/50. dx=(float(mrr(md,1))/y0)/50. do m=1,md gmin=rmnum(m)/y0+dx gmax=rmnum(m)/y0-dx c gmin=(rmnum(i)-.3)/y0 c gmax=(rmnum(i)+.3)/y0 if(avekp2(m) /= 0) then c if (mrr(m,1) /= 0) then gdx=mgamx(m,1,net)/avekp2(m) c gdx=mgamx(m,1,net)/(float(mrr(m,1))/y0)**2/ c & (1.+(shr*aveth(m))**2) else gdx=0.0 endif call line(gmin,gdx+dy,gmin,gdx-dy) call line(gmax,gdx+dy,gmax,gdx-dy) call line(gmin,gdx+dy,gmax,gdx+dy) call line(gmin,gdx-dy,gmax,gdx-dy) c call line(gmin,0.0,gmin,gdx) c call line(gmin,gdx,gmax,gdx) c call line(gmax,0.0,gmax,gdx) enddo call finish(ipage, date) endif !end of igraph=1 section write(1,*) write(9,*) 'm k_theta gamma D_x D' write(1,*) 'm k_theta gamma D_x D' c BD: c Makeshift definition for avekx2 if igeo=1 c do m=1,md avekx2(m)=(float(mrr(m,1))/y0*shr*aveth(m))**2 enddo do m=1,md if((avekp2(m) /= 0.).and.(avekx2(m) /= 0.)) then c if ((mrr(m,1) /= 0).and.(aveth(m) /= 0.)) then write(9,33) mrr(m,1),float(mrr(m,1))/y0,mgamx(m,1,net), & aveth(m),mgamx(m,1,net)/avekx2(m), & mgamx(m,1,net)/avekp2(m) write(1,33) mrr(m,1),float(mrr(m,1))/y0,mgamx(m,1,net), & aveth(m),mgamx(m,1,net)/avekx2(m), & mgamx(m,1,net)/avekp2(m) c write(9,33) mrr(m,1),float(mrr(m,1))/y0,mgamx(m,1,net), c & aveth(m), c & mgamx(m,1,net)/(float(mrr(m,1))/y0*shr*aveth(m))**2, c & mgamx(m,1,net)/(float(mrr(m,1))/y0)**2/ c & (1.+(shr*aveth(m))**2) else write(9,33) mrr(m,1),float(mrr(m,1))/y0,mgamx(m,1,net), & aveth(m),0.,0. write(1,33) mrr(m,1),float(mrr(m,1))/y0,mgamx(m,1,net), & aveth(m),0.,0. endif enddo 33 format(i5,f10.3,f10.3,f10.3,f10.3,f10.3) write(9,*) write(1,*) close(1) endif ! end linear stuff c c Close graphics output file: c if(igraph == 1) then CALL GDAWK (1) ! deactivate workstation CALL GCLWK (1) ! close workstation CALL GCLKS ! close GKS endif call finish_mp ! finish MPI stop end