subroutine itg(icontrol) ! ! (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. ! ! See itg.hlp for more info. ! use mp use itg_data use gryffin_grid, only: md, rkperp2 use fields use linear, only: flr, flrinit, emflr, filter use step, only: timestep, courant use diagnostics, only: prstep, volsq, posten, preenr use diagnostics_arrays use nonlinear use bounce_avg use field_eq use gryffin_layouts use dmix implicit none complex, dimension(:,:,:,:), allocatable :: old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp complex, dimension(:,:,:), allocatable :: old_pot, save_phi, & old_e_density, old_e_p, old_e_r, old_e_t, phi_bat, old_n_e, old_u_e, & old_tpar_e, old_tperp_e, old_apar real, dimension(:), allocatable :: hdt real g_target,error,d_target real t_start integer m,n,ncyclea,ncycleb,icontrol(2) integer n_time_diag ! ! Read namelists and set defaults call input ! Set up arrays, constants, and ic's call init ! Allocate local arrays call allocate_local icontrol(2)=1 t_start = 0. if(ifilter == -1 .and. idx_local(m_lo, md, 1)) t_start=time(md) call broadcast(t_start, proc_id(m_lo, md, 1)) ! ! **************************************************** ! Top of the loop ! **************************************************** 100 continue ncycle=ncycle+1 ! ! Find L_T_crit: ! if(icrit /= 0 .and. ncycle >= 26) call critical(g_target, d_target, icontrol(2)) ! ! The code advances the fields in time using a simple explicit 2-step ! Runge-Kutta algorithm that is second order accurate in time. ! ! Take the first step with hdt=0.5*dt ! ! Calculate FLR corrections to the potential ! call flr(potential) if (beta_e > 0.0) call emflr(apar) if (epse > 0.0) call baphi(potential) ! Evaluate non-linear terms if needed if(lin /= 1) then call nonlin(density, u_par, t_par, q_par, t_perp, q_perp, potential) if (beta_e > 0.) call nonlin_em(density, u_par, t_par, q_par, t_perp, & q_perp, potential, n_e, u_e, tpar_e, tperp_e, apar) if (epse > 0.0) call nonline(e_density, e_p, e_r, phi_ba) endif ! ! Simple-minded way to adjust the time step for nonlinear runs to the ! radial Courant condition ! if(lin /= 1 .and. idt == 1) then if (.not.((nread == 1) .and. (ncycle < 3))) call courant endif ! time = time + dt hdt =0.5*dt ! ! Take a half time step. ! 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, & old_n_e, old_u_e, old_tpar_e, old_tperp_e, old_apar, & n_e, u_e, tpar_e, tperp_e, apar, potential, hdt) ! ! Solve Poisson equation ! if (epse > 0.0) call kapav(phi_ba, old_e_density) call poisson(old_pot, old_n_e, old_n, old_tperp) ! ! If ifilter=1, include (hyper-)viscosity-like terms in each equation. ! (n,v,tpar,qpar,T,q) ! if (iabs(ifilter) == 1) then call filter(old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp, hdt) endif ! !***************************************************************** ! Take the second step with a full dt to find values at (n+1)*dt !***************************************************************** ! ! Evaluate various FLR corrections to phi. ! call flr(old_pot) if (beta_e > 0.0) call emflr(old_apar) if (epse > 0.0) call baphi(old_pot) ! ! Evaluate non-linear terms if needed. ! if(lin /= 1) then call nonlin (old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp, old_pot) if (beta_e > 0.) call nonlin_em (old_n, old_upar, old_tpar, & old_qpar, old_tperp, old_qperp, old_pot, & old_n_e, old_u_e, old_tpar_e, old_tperp_e, old_apar) if (epse > 0.0) call nonline(old_e_density, old_e_p, old_e_r, phi_ba) endif ! ! Take a full time step. ! 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, n_e, & u_e, tpar_e, tperp_e, apar, old_n_e, old_u_e, old_tpar_e, & old_tperp_e, old_apar, old_pot, dt) ! ! check particle conservation ! ! if (mod(ncycle,nprnt) == 0) then ! ! save old phi_ba, so it's not overwritten ! ! phi_bat = phi_ba ! call baphi(potential) ! call pcons ! ! reload old phi_ba ! ! phi_ba = phi_bat ! endif ! ! Solve Poisson equation: ! if (epse > 0.0) call kapav(phi_ba,e_density) call poisson(potential, n_e, density, t_perp) ! ! If ifilter = 1, include (hyper-)viscosity-like terms in each equation. ! (n,v,tpar,qpar,T,q) ! if (iabs(ifilter) == 1) then call filter(density, u_par, t_par, q_par, t_perp, q_perp, dt) endif ! ! Keep time history of the potential. ! if(mod(ncycle, nfreq) == 0) then nfreqcnt=nfreqcnt+1 utim(nfreqcnt,:,:)=potential(lhistory,:,:) if(idx_local(m_lo, md, 1)) timf(nfreqcnt)=time(md) call broadcast(timf(nfreqcnt), proc_id(m_lo, md, 1)) endif ! ! Check the energy balance. ! ncyclea=ncycle+1 ncycleb=ncycle-1 if(mod(ncyclea,nprnt) == 0) call preenr if(ncycle /= 1 .and. mod(ncycleb,nprnt) == 0) call posten ! ! Save some results for printing later. ! if(mod(ncycle,nprnt) == 0 .or. (ncycle == 1)) then ne = ne + 1 call prstep(t_start) if(icrit /= 0 .and. ncycle >= 26) & call tc_out(g_target, d_target, error, icontrol, save_phi) ! ! if ifilter = -1, use time-dependent viscosity-like terms in each equation. ! (n,v,tpar,qpar,T,q) proportional to sqrt(int(dt chi)) ! if(ifilter == -1) call rdiffset endif if(mod(ncycle+1, nprnt) == 0 .and. icrit /= 0) save_phi=potential if(ntrace == 0) goto 1000 if(mod(ncycle,nprnt) == 0) then if(ntrace == 1) then if (epse > 0.0) then if(idx_local(m_lo, md, 1)) & write(*,2000) ncycle,dt(md),wpfx(ne,1),qfluxe(ne), & wpfx(ne,1)/qfluxe(ne),fluxe(ne),fluxi(ne,1) else if(ifilter == -1) then if(idx_local(m_lo, 1, 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) endif else if (iphi00 == 3) then if(idx_local(m_lo, md, 1)) write(*,2003) ncycle,dt(md),wpfx(ne,1) else if(idx_local(m_lo, md, 1)) write(*,2001) ncycle,dt(md),wpfx(ne,1) endif endif endif else if(idx_local(m_lo, md, 1)) then if (epse > 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 endif 2000 format(' n: ',i5,' 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: ',i5,' dt: ',f6.4,' Q_i: ',1pe10.2,' adiabatic electrons') 2002 format(' n: ',i5,' D: ',1pe10.2,' chi: ',1pe10.2) 2003 format(' n: ',i5,' dt: ',f6.4,' Q_i: ',1pe10.2,' adiabatic ions') 1000 continue if(icontrol(1) == 1) goto 200 if(ncycle < nstp) goto 100 200 continue ! ! take out B variation of kperp rho for post gamma/k^2 calculation ! do n=n_low, n_high do m=m_low, m_high rkperp2(:,m,n)=rkperp2(:,m,n)/bmaginv(:)**2 enddo enddo contains subroutine init ! ! Set up arrays and define values of constants. ! use itg_data, only: alloc_itg, init_itg, geo, s_alpha use linear, only: init_kparfft, flrinit use nlps, only: init_nlps use bounds, only: init_theta_fft use fields, only: alloc_fields use fields_init, only: init_fields use gryffin_grid, only: init_grid use bounce_avg, only: init_bounce_avg use nonlinear, only: alloc_nonlinear use diagnostics, only: init_diagnostics use diagnostics_arrays, only: alloc_diagnostics ! Initialize grid quantities call init_grid (epse, iperiod, lin, linlay) ! Allocate arrays (mostly diagnostics) call alloc_itg ! Allocate main 3-D arrays call alloc_fields ! Allocate arrays for nonlinear terms call alloc_nonlinear ! Initialize fft's for nonlinear terms if(lin == 0) call init_nlps ! Initialize remaining physics and control parameters call init_itg ! Do periodicity-dependent initialization for kpar call init_theta_fft ! in module bounds ! Do remaining initialization for kpar call init_kparfft ! in module linear ! Initialize geometry-dependent coefficients if(igeo == 1) then call geo else call s_alpha endif ! allocate some diagnostic arrays ! the logic here is that info is saved (printed) every nprnt steps, ! and also on the first step, so the total # of saved time points (n_time_diag) ! in the diagnostic arrays is determined by the next two lines: n_time_diag=1+nstp/nprnt if(nprnt == 1) n_time_diag=nstp call alloc_diagnostics(n_time_diag, nstp/nfreq) ! Initialize the main fields call init_fields ! in module fields_init ! Initialize diagnostics call init_diagnostics ! Initialize bounce-averaged electron model if (epse > 0.) call init_bounce_avg ! Initialize weights for FLR terms call flrinit(1) end subroutine init subroutine allocate_local use gryffin_grid, only: kdpass, ld allocate (hdt(m_low:m_alloc)) hdt = 0. allocate (old_n(ld, m_low:m_alloc,n_low:n_alloc,nspecies), & old_upar(ld, m_low:m_alloc,n_low:n_alloc,nspecies), & old_tpar(ld, m_low:m_alloc,n_low:n_alloc,nspecies), & old_qpar(ld, m_low:m_alloc,n_low:n_alloc,nspecies), & old_tperp(ld,m_low:m_alloc,n_low:n_alloc,nspecies), & old_qperp(ld,m_low:m_alloc,n_low:n_alloc,nspecies)) old_n = 0.; old_upar = 0.; old_tpar = 0.; old_qpar = 0.; old_tperp = 0.; old_qperp = 0. allocate (old_pot(ld,m_low:m_alloc,n_low:n_alloc), & save_phi(ld,m_low:m_alloc,n_low:n_alloc)) old_pot = 0.; save_phi = 0. if(epse > 0.) then allocate (old_e_density(kdpass, m_low:m_alloc, n_low:n_alloc), & old_e_p(kdpass, m_low:m_alloc, n_low:n_alloc), & old_e_r(kdpass, m_low:m_alloc, n_low:n_alloc), & old_e_t(kdpass, m_low:m_alloc, n_low:n_alloc), & phi_bat(kdpass, m_low:m_alloc, n_low:n_alloc)) old_e_density = 0. old_e_p = 0. old_e_r = 0. old_e_t = 0. phi_bat = 0. endif if (beta_e > 0.) then allocate (old_n_e(ld, m_low:m_alloc, n_low:n_alloc), & old_u_e(ld, m_low:m_alloc, n_low:n_alloc), & old_tpar_e(ld, m_low:m_alloc, n_low:n_alloc), & old_tperp_e(ld, m_low:m_alloc, n_low:n_alloc), & old_apar(ld, m_low:m_alloc, n_low:n_alloc)) old_n_e = 0. old_u_e = 0. old_tpar_e = 0. old_tperp_e = 0. old_apar = 0. endif end subroutine allocate_local subroutine rdiffset ! ! This routine calculates the diffusion associated with particle noise ! in GK PIC simulations. The code evaluates an analytical derivation ! by Greg Hammett. The parameters ifilter, N_p, n_zpic, and bmax are ! described in the input file (itg.in). ! use constants, only: pi integer m, n real tmp, fa_par, fac fa_par=2*pi*qsf/epsn/float(n_zpic) ! a_par/L_n do n = n_low, n_high do m = m_low, m_high tmp=0.5*sqrt(pi/2)*eta(1)**2/float(n_p)*fa_par*chi_int*bmax/(2+bmax) ! first time step hack if(rdiff(m,n) == 0.) then rdiff(m,n)=tmp else fac=1+sqrt(2/pi)/fa_par/(rdiff(m,n)-wpfx(ne,1)/eta(1)/1.5) & *2/bmax*(2+bmax)/(1+bmax) rdiff(m,n)=tmp*alog(fac) endif enddo enddo end subroutine rdiffset end subroutine itg