subroutine itg(icontrol,time,dt) c c (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, c and Gregory W. Hammett. ALL RIGHTS RESERVED. c c See itg.hlp for more info. c implicit none include 'itg.par' include 'itg.cmn' complex old_n(lz,mz,nz,nspecz),old_upar(lz,mz,nz,nspecz), * old_tpar(lz,mz,nz,nspecz),old_qpar(lz,mz,nz,nspecz), * old_tperp(lz,mz,nz,nspecz),old_qperp(lz,mz,nz,nspecz), * old_pot(lz,mz,nz),save_phi(lz,mz,nz) complex old_e_density(kz,mz,nz),old_e_p(kz,mz,nz), * old_e_r(kz,mz,nz),old_e_t(kz,mz,nz) complex phi_bat(kz,mz,nz) real time(mz),dt(mz),hdt(mz),g_target,error,d_target real t_start, D_noise integer i,j,l,m,n,ncyclea,ncycleb,icontrol(2) c input files: c c 12 RUNNAME.in main namelist input file to control the run. c 5 RUNNAME.resp a copy of itg.res to continue from a previous run. c c output files: c c 7 RUNNAME.res a complete dump of the common blocks. c 9 RUNNAME.out printout of various results c 20 RUNNAME.fields phi(x,y) for particle code. c c ************************************************************** c c Read namelists and set defaults call input(dt) c Set up arrays, constants, and ic's call init(old_qpar,old_qperp,old_e_density,time,dt) c Write initial conditions call echoin(time) icontrol(2)=1 if(ifilter.eq.-1) t_start=time(md) c c **************************************************** c Top of the loop c **************************************************** 100 continue c ncycle=ncycle+1 c c Find L_T_crit: c if(icrit.ne.0 .and. ncycle.ge.26) . call critical(g_target,d_target,time,icontrol(2)) c c The code advances the fields in time using a simple explicit 2-step c Runge-Kutta algorithm that is second order accurate in time. c c Take the first step with hdt=0.5*dt c c Calculate FLR corrections to the potential c call flr(potential) if (epse.gt.0.0) call baphi(potential) c c Evaluate non-linear terms if needed c if(lin.ne.1) then call nonlin(density,u_par,t_par,q_par,t_perp,q_perp,potential) if (epse.gt.0.0) call nonline(e_density,e_p,e_r,phi_ba) endif c c Simple-minded way to adjust the time step for nonlinear runs to the c radial Courant condition c if(lin.ne.1 .and. idt.eq.1) then if (.not.((nread.eq.1).and.(ncycle.lt.3))) then call courant(dt) endif endif c do m=1,md time(m)=time(m)+dt(m) hdt(m)=0.5*dt(m) enddo c c Take a half time step. c call timestep(old_n,old_upar,old_tpar,old_qpar,old_tperp, * old_qperp, * density,u_par,t_par,q_par,t_perp,q_perp, * old_e_density,old_e_p,old_e_r,old_e_t, * e_density,e_p,e_r,e_t,hdt) c c Electron collisions c c call ecollis(old_e_density,old_e_p,old_e_r,old_e_t, c . e_density,e_p,e_r,e_t,hdt) c c Solve Poisson equation c if (epse.gt.0.0) call kapav(phi_ba,old_e_density) call poisson(old_pot,old_n,old_tperp) c c If ifilter=1, include (hyper-)viscosity-like terms in each equation. c (n,v,tpar,qpar,T,q) c if (iabs(ifilter).eq.1) then call filter(old_n,old_upar,old_tpar,old_qpar,old_tperp, . old_qperp,hdt) endif c c***************************************************************** c Take the second step with a full dt to find values at (n+1)*dt c***************************************************************** c c Evaluate various FLR corrections to phi. c call flr(old_pot) if (epse.gt.0.0) call baphi(old_pot) c c Evaluate non-linear terms if needed. c if(lin.ne.1) then call nonlin(old_n,old_upar,old_tpar,old_qpar,old_tperp, * old_qperp,old_pot) if (epse.gt.0.0) call nonline(old_e_density,old_e_p,old_e_r,phi_ba) endif c c Take a full time step. c call timestep(density,u_par,t_par,q_par,t_perp,q_perp, * old_n,old_upar,old_tpar,old_qpar,old_tperp,old_qperp, * e_density,e_p,e_r,e_t, * old_e_density,old_e_p,old_e_r,old_e_t,dt) c c check particle conservation c c if (mod(ncycle,nprnt).eq.0) then c c save old phi_ba, so it's not overwritten c c do n=1,nd c do m=1,md c do j=1,kdpass c phi_bat(j,m,n)=phi_ba(j,m,n) c enddo c enddo c enddo c call baphi(potential) c call pcons c c reload old phi_ba c c do n=1,nd c do m=1,md c do j=1,kdpass c phi_ba(j,m,n)=phi_bat(j,m,n) c enddo c enddo c enddo c endif c c Solve Poisson equation: c if (epse.gt.0.0) call kapav(phi_ba,e_density) call poisson(potential,density,t_perp) c c If ifilter = 1, include (hyper-)viscosity-like terms in each equation. c (n,v,tpar,qpar,T,q) c if (iabs(ifilter).eq.1) then call filter(density,u_par,t_par,q_par,t_perp,q_perp,dt) endif c c Keep time history of the potential. c if(mod(ncycle,nfreq).eq.0) then timf(nfreqcnt)=time(md) if(nfreqcnt.eq.1) ntimf=ncycle do n=1,nd do m=1,md utim(nfreqcnt,m,n)=potential(lhistory,m,n) enddo enddo nfreqcnt=nfreqcnt+1 ntotal=nfreqcnt-1 endif c c Check the energy balance. c ncyclea=ncycle+1 ncycleb=ncycle-1 c if(mod(ncyclea,nprnt).eq.0) then call preenr endif c if(ncycle.ne.1.and.mod(ncycleb,nprnt).eq.0) then call posten(dt) endif c c Save some results for printing later. c if(mod(ncycle,nprnt).eq.0 .or. (ncycle.eq.1)) then ne=ne+1 call prstep(t_start,dt,time) if(icrit.ne.0 .and. ncycle.ge.26) . call tc_out(save_phi,g_target,d_target, . error,time,dt,icontrol) c c if ifilter = -1, use time-dependent viscosity-like terms in each equation. c (n,v,tpar,qpar,T,q) proportional to sqrt(int(dt chi)) c if(ifilter.eq.-1) call rdiffset endif if(mod(ncycle+1,nprnt).eq.0 .and. icrit.ne.0) then do n=1,nd do m=1,md do l=1,ldb save_phi(l,m,n)=potential(l,m,n) enddo enddo enddo endif c c Make a movie frame of the potential if desired. c if(ncycle.gt.movieon) then if(mod(ncycle,ninterv).eq.0 .or. ncycle.eq.1) then call wframe endif endif c if(ntrace.eq.0) goto 1000 if(mod(ncycle,nprnt).eq.0) then if(ntrace.eq.1) then if (epse.gt.0.0) then write(*,2000) ncycle,dt(md),wpfx(ne,1),qfluxe(ne), & wpfx(ne,1)/qfluxe(ne),fluxe(ne),fluxi(ne,1) else if(ifilter.eq.-1) then write(*,2002) ncycle,rdiff(1,1),-wpfx(ne,1)/eta(1) write(123,*) timo(ne),rdiff(1,1),-wpfx(ne,1)/eta(1), & chi_int,1.5*rdiff(1,1)-wpfx(ne,1)/eta(1) else write(*,2001) ncycle,dt(md),wpfx(ne,1) endif endif else if (epse.gt.0.0) then write(59,2000) ncycle,dt(md),wpfx(ne,1),qfluxe(ne), & wpfx(ne,1)/qfluxe(ne),fluxe(ne),fluxi(ne,1) else write(59,2001) ncycle,dt(md),wpfx(ne,1) endif endif endif 2000 format(' n: ',i4,' dt: ',f6.4,' Q_i: ',1pe10.2,' Q_e: ',1pe10.2, & ' Q_i/Q_e: ',1pe10.2, ' Gamma_e: ',1pe10.2,' Gamma_i: ',1pe10.2) 2001 format(' n: ',i4,' dt: ',f6.4,' Q_i: ',1pe10.2,' adiabatic electrons') 2002 format(' n: ',i4,' D: ',1pe10.2,' chi: ',1pe10.2) 1000 continue if(icontrol(1).eq.1) goto 200 if(ncycle.lt.nstp) goto 100 200 continue c c take out B variation of kperp rho for post gamma/k^2 calculation c do l=1,ld do m=1,md do n=1,nd rkperp2(l,m,n)=rkperp2(l,m,n)/bmaginv(l)**2 enddo enddo enddo return end