module step use itg_data implicit none contains subroutine timestep(density1,u_par1,t_par1,q_par1,t_perp1,q_perp1, & density2,u_par2,t_par2,q_par2,t_perp2,q_perp2, & e_density1,e_p1,e_r1,e_t1, & e_density2,e_p2,e_r2,e_t2,dt) ! ! (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, ! and Gregory W. Hammett. ALL RIGHTS RESERVED. ! ! mbeer -- switched to Cowley coordinates, May 8, 1992 ! ! Revision 3.0 1994/02/03 22:53:53 bdorland ! Upgraded to complex notation. ! ! Advance equations one time step, not including collisional ! dissipation ! use flr_terms, only: phi_n, phi_u, phi_tpar, phi_tperp, phi_qperp, phi_rperp, & phi_flr1, phi_flr2, phi_flr3, phi_nd, phi_tpard, phi_tperpd, phi_qperpb complex, dimension(:,:,:,:) :: density1, u_par1, t_par1, q_par1, & t_perp1, q_perp1, density2, u_par2, t_par2, q_par2, t_perp2, q_perp2 complex, dimension(:,:,:) :: e_density1, e_p1, e_r1, e_t1, & e_density2, e_p2, e_r2, e_t2 complex, dimension(size(density1, 1), size(density1, 2), & size(density1, 3), size(density1, 4)) :: kp_density2, kp_upar2, & kp_tpar2, kp_qpar2, kp_tperp2, kp_qperp2, kp_vphi, akp_qpar2, & akp_qperp2, qparps, qperpps complex, dimension(size(e_p1, 1),size(e_p1, 2),size(e_p1, 3)) :: wd_t, awd_t real, dimension(:) :: dt ! real, dimension(size(density1, 1), size(density1, 2), & ! size(density1, 3), size(density1, 4)) :: dt4 ! real, dimension(size(e_p1, 1), size(e_p1, 2), size(e_p1, 3)) :: dt3 interface subroutine bcon(a) complex, dimension(:,:,:) :: a end subroutine bcon end interface interface subroutine kparfft(a, i) complex, dimension(:,:,:,:) :: a integer :: i end subroutine kparfft end interface real nuur,nuui,vmom,nu_erk,nuky,znorm,dz0 real, dimension(size(dt)) :: switch integer i,l,m,n real, dimension(size(density1,1)) :: dzt complex, dimension(size(density1,3), size(density1,4)) :: tparfsa,tperpfsa dz0=2.*x0/ldb do l=l_left,l_right dzt(l)=dz0*jacobian(l) enddo znorm=0.0 do l=l_left,l_right znorm=znorm+dzt(l) enddo do l=l_left,l_right dzt(l)=dzt(l)/znorm enddo ! ! turn off omegad dissipative terms for m=0 if (iflow == 1) then switch(1)=1.0 else switch(1)=0.0 endif do m=2,md switch(m)=1.0 enddo if (nparmom == 3) then nuur=nu5r nuui=nu5i vmom=0. else nuur=0.0 nuui=0.0 vmom=1. endif ! if(lin == 0) dt4 = spread(spread(spread(dt,1,ld),3,nd),4,nspecies) do n=1,nd do i=1,nspecies do m=1,md ! ! load kp_upar2,kp_vphi,etc. ! do l=1,ld-1 kp_density2(l,m,n,i)=density2(l,m,n,i)*bmaginv(l) kp_upar2(l,m,n,i)=u_par2(l,m,n,i)*bmaginv(l) kp_tpar2(l,m,n,i)=t_par2(l,m,n,i)*bmaginv(l) kp_vphi(l,m,n,i)=phi_u(l,m,n,i) enddo kp_density2(ld,m,n,i)=0. kp_upar2(ld,m,n,i)=0. kp_tpar2(ld,m,n,i)=0. kp_vphi(ld,m,n,i)=0. enddo enddo enddo call kparfft(kp_density2,0) call kparfft(kp_upar2,0) call kparfft(kp_tpar2,0) call kparfft(kp_vphi,0) do n=1,nd do i=1,nspecies do m=1,md do l=1,ld ! ! -----------------density equation---------------------- ! density1(l,m,n,i)=density(l,m,n,i)-dt(m)*( & nl_density1(l,m,n,i)+nl_density2(l,m,n,i) & +ii*rky(m,n)/Ln(i)*phi_n(l,m,n,i) & +ii*kp_upar2(l,m,n,i)*b_mag(l)*vt(i) & -ii*omegad(l,m,n,i)*(phi_nd(l,m,n,i)+2.*density2(l,m,n,i)) & -ii*omegad(l,m,n,i)*(t_par2(l,m,n,i)+t_perp2(l,m,n,i))) ! ! -------------parallel momentum equation---------------- ! u_par1(l,m,n,i)=u_par(l,m,n,i)-dt(m)*( & nl_upar1(l,m,n,i)+nl_upar2(l,m,n,i) & +ii*kp_vphi(l,m,n,i)*vt(i)*charge(i)/tau(i) & +ii*(kp_tpar2(l,m,n,i)+kp_density2(l,m,n,i))*b_mag(l)*vt(i) & -ii*omegad(l,m,n,i)*2.*((2.+nuui*switch(m))*u_par2(l,m,n,i) & +0.5*vmom*(q_par2(l,m,n,i)+q_perp2(l,m,n,i))) & +2.*nuur*abs(omegad(l,m,n,i)*switch(m))*u_par2(l,m,n,i) & +(t_perp2(l,m,n,i)+density2(l,m,n,i)+phi_qperp(l,m,n,i) & *charge(i)/tau(i))*bgrad(l)*vt(i)) enddo enddo enddo enddo ! ! --------------nonlinear terms: density and upar -------- ! ! if(lin == 0) then ! density1 = density1 - dt4*(nl_density1+nl_density2) ! u_par1 = u_par1 - dt4*(nl_upar1+nl_upar2) ! endif ! ! ---------three pole closure: parallel moments---------- ! if(nparmom == 3) q_par2 = rkkin * t_par2 ! ! load kp_qpar2 ! do n=1,nd do i=1,nspecies do m=1,md do l=1,ld-1 kp_qpar2(l,m,n,i)=q_par2(l,m,n,i)*bmaginv(l) enddo kp_qpar2(ld,m,n,i)=0.0 enddo enddo enddo if (nparmom == 3) then call kparfft(kp_qpar2,2) else call kparfft(kp_qpar2,0) endif do n=1,nd do i=1,nspecies do m=1,md do l=1,ld ! ! -------------parallel temperature equation----------- ! t_par1(l,m,n,i)=t_par(l,m,n,i)-dt(m)*( & nl_tpar(l,m,n,i) & +2.*ii*kp_upar2(l,m,n,i)*b_mag(l)*vt(i) & +ii*rky(m,n)/Ln(i)*phi_tpar(l,m,n,i) & +ii*kp_qpar2(l,m,n,i)*b_mag(l)*vt(i) & -ii*omegad(l,m,n,i)*2.* & ( 1.5*(density2(l,m,n,i)+t_par2(l,m,n,i)) & ! orig +(1.5+nu1i*switch(m))*t_par2(l,m,n,i) & ! figs +(1.5-3+(nu1i+3)*switch(m))*t_par2(l,m,n,i) & +(1.5+1+(nu1i-1)*switch(m))*t_par2(l,m,n,i) & ! +(1.5-2)*t_par2(l,m,n,i) & ! orig +(0.5+nu2i*switch(m))*t_perp2(l,m,n,i) & ! figs +(0.5+1+(nu2i-1)*switch(m))*t_perp2(l,m,n,i) & +(0.5-1+(nu2i+1)*switch(m))*t_perp2(l,m,n,i) & -0.5*(t_perp2(l,m,n,i)+density2(l,m,n,i)) & +0.5*phi_tpard(l,m,n,i) ) & +2.*nu1r*abs(omegad(l,m,n,i)*switch(m))*t_par2(l,m,n,i) & +2.*nu2r*abs(omegad(l,m,n,i)*switch(m))*t_perp2(l,m,n,i) & +2.*(u_par2(l,m,n,i)+q_perp2(l,m,n,i))*bgrad(l)*vt(i) & +2./3.*nu_ii(i)*(t_par2(l,m,n,i)-t_perp2(l,m,n,i)) ) enddo enddo enddo enddo ! ! -------nonlinear t_par term -------------- ! ! if(lin == 0) t_par1 = t_par1 - dt4*nl_tpar if(nparmom == 4) then ! ! load kp_tpar2,akp_qpar2 ! qparps=0. qperpps=0. if (iflow == 2) then do i=1,nspecies do n=1,nd tparfsa(n,i)=0.0 tperpfsa(n,i)=0.0 do l=l_left,l_right tparfsa(n,i)=t_par2(l,1,n,i)*dzt(l) tperpfsa(n,i)=t_perp2(l,1,n,i)*dzt(l) enddo do l=1,ld qparps(l,1,n,i)=3*ii*rkx(l,1,n)*qsf/eps*bmaginv(l) & *tparfsa(n,i) qperpps(l,1,n,i)=ii*rkx(l,1,n)*qsf/eps*bmaginv(l) & *tperpfsa(n,i) enddo enddo enddo endif kp_tpar2 = t_par2 akp_qpar2 = q_par2-qparps kp_tpar2(ld,:,:,:) = 0. akp_qpar2(ld,:,:,:) = 0. call kparfft(kp_tpar2,0) call kparfft(akp_qpar2,1) do n=1,nd do i=1,nspecies do m=1,md do l=1,ld ! ! --------------parallel heat flux equation--------------- ! q_par1(l,m,n,i)=q_par(l,m,n,i)-dt(m)*( & nl_qpar(l,m,n,i) & +zBq1*ii*kp_tpar2(l,m,n,i)*vt(i) & +zDq1*akp_qpar2(l,m,n,i)*vt(i) & -ii*omegad(l,m,n,i)* & ( (nu5i*switch(m)+6.)*u_par2(l,m,n,i) & + (nu6i*switch(m)-3.)*q_par2(l,m,n,i) & + (nu7i*switch(m)-3.)*q_perp2(l,m,n,i) ) & +abs(omegad(l,m,n,i)*switch(m))* & ( nu5r*u_par2(l,m,n,i) & + nu6r*q_par2(l,m,n,i) & + nu7r*q_perp2(l,m,n,i) ) & +nu_ii(i)*q_par2(l,m,n,i) ) enddo enddo enddo enddo ! ! ------- nonlinear q_par term -------------------- ! ! if(lin == 0) q_par1 = q_par1 - dt4*nl_qpar endif ! ! ----------one pole closure: perpendicular moments------------ ! if(nperpmom == 1) then do n=1,nd do i=1,nspecies do m=1,md do l=1,ld q_perp2(l,m,n,i)= rkkperp*(t_perp2(l,m,n,i) & +charge(i)/tau(i)*phi_qperp(l,m,n,i)) enddo enddo enddo enddo endif ! ! load kp_qperp2 ! do n=1,nd do i=1,nspecies do m=1,md do l=1,ld-1 kp_qperp2(l,m,n,i)=q_perp2(l,m,n,i)*bmaginv(l)**2 enddo kp_qperp2(ld,m,n,i)=0.0 enddo enddo enddo if (nperpmom == 1) then call kparfft(kp_qperp2,2) else call kparfft(kp_qperp2,0) endif do n=1,nd do i=1,nspecies do m=1,md do l=1,ld ! ! ----------perpendicular temperature equation-------------- ! t_perp1(l,m,n,i)=t_perp(l,m,n,i)-dt(m)*( & nl_tperp1(l,m,n,i)+nl_tperp2(l,m,n,i)+nl_tperp3(l,m,n,i) & +ii*rky(m,n)/Ln(i)*phi_tperp(l,m,n,i) & +ii*kp_qperp2(l,m,n,i)*b_mag(l)**2*vt(i) & -ii*2.*omegad(l,m,n,i)* & ( t_perp2(l,m,n,i)+density2(l,m,n,i) & +(0.5+nu3i*switch(m))*t_par2(l,m,n,i) & ! orig +(1.+nu4i*switch(m))*t_perp2(l,m,n,i) & ! figs +(1.-1.5+(nu4i+1.5)*switch(m))*t_perp2(l,m,n,i) & +(1.+(nu4i)*switch(m))*t_perp2(l,m,n,i) & ! +(1.-1.5)*t_perp2(l,m,n,i) & -0.5*(density2(l,m,n,i)+t_par2(l,m,n,i)) & +0.5*phi_tperpd(l,m,n,i) ) & +2.*nu4r*abs(omegad(l,m,n,i)*switch(m))*t_perp2(l,m,n,i) & +2.*nu3r*abs(omegad(l,m,n,i)*switch(m))*t_par2(l,m,n,i) & -u_par2(l,m,n,i)*bgrad(l)*vt(i) & +nu_ii(i)/3.*(t_perp2(l,m,n,i)-t_par2(l,m,n,i)) ) enddo enddo enddo enddo ! ! --------- t_perp nonlinear terms ----------------------- ! ! if(lin == 0) t_perp1 = t_perp1 - dt4*(nl_tperp1+nl_tperp2+nl_tperp3) if(nperpmom >= 2) then ! ! load kp_tperp2,akp_qperp2 ! do n=1,nd do i=1,nspecies do m=1,md do l=1,ld-1 akp_qperp2(l,m,n,i)=(q_perp2(l,m,n,i)-qperpps(l,m,n,i)) & *rkkperp2 kp_tperp2(l,m,n,i)=t_perp2(l,m,n,i) & +charge(i)/tau(i)*phi_qperp(l,m,n,i) enddo akp_qperp2(ld,m,n,i)=0.0 kp_tperp2(ld,m,n,i)=0.0 enddo enddo enddo call kparfft(akp_qperp2,1) call kparfft(kp_tperp2,0) do n=1,nd do i=1,nspecies do m=1,md do l=1,ld ! ! -----------perpendicular heat flux equation---------------- ! q_perp1(l,m,n,i)=q_perp(l,m,n,i)-dt(m)*( & nl_qperp1(l,m,n,i)+nl_qperp2(l,m,n,i) & +nl_qperp3(l,m,n,i) & +ii*kp_tperp2(l,m,n,i)*vt(i) & +akp_qperp2(l,m,n,i)*vt(i) & -ii*omegad(l,m,n,i)* & ( (nu8i*switch(m)+1.)*u_par2(l,m,n,i) & + (nu9i*switch(m)-1.)*q_par2(l,m,n,i) & + (nu10i*switch(m)-1.)*q_perp2(l,m,n,i) ) & +abs(omegad(l,m,n,i)*switch(m))* & ( nu8r*u_par2(l,m,n,i) & + nu9r*q_par2(l,m,n,i) & + nu10r*q_perp2(l,m,n,i) ) & +(t_perp2(l,m,n,i)-t_par2(l,m,n,i)+phi_qperpb(l,m,n,i) & *charge(i)/tau(i))*bgrad(l)*vt(i) & +nu_ii(i)*q_perp2(l,m,n,i) ) enddo enddo enddo enddo ! ! ------- nonlinear q_perp terms --------------------- ! ! if(lin == 0) q_perp1 = q_perp1 - dt4*(nl_qperp1+nl_qperp2+nl_qperp3) endif ! ! bounce averaged electrons: ! if (epse > 0.0) then if (nemom == 3) then wd_t = nuai*e_density2 + nubi*e_p2 + nuci*e_r2 awd_t = nuar*e_density2 + nubr*e_p2 + nucr*e_r2 else wd_t = e_t2 awd_t = 0. endif ! nuky=rmu1 for hyperviscosity in electron eqns nuky=0 ! if (ncycle >= 4000) nuky=nuky/10. do n=1,nd do m=1,md do l=1,kd ! Try Rewoldt's version of a Krook model collision operator, as given ! in Eq. 4 and 6 of Kotschenreuther, Rewoldt, Tang, "Comparison of Initial ! Value and Eigenvalue Codes for Kinetic Toroidal Plasma Instabilities", ! PPPL-2986 (1994) (to be publ. in Comp. Phys. Comm.?). ! ! One of the interesting features of this model is that the loss rate ! goes to infinity near the trapped-passing boundary, so we will treat ! this term implicitly below. For simplicity for now, we will ignore ! the energy dependence of the collision frequency so the collision ! operator is almost diagonal (except for the phi_ba term...). ! (Also, take the delta=0 limit of Rewoldt's full formula.) ! ! nueeff = nu_e_effective = nu_e / (r/R) ! ! nu_erk=nueeff/4.0*1.31/(1-kap(l)**2)**2 ! For now, turn off the Krook model, and try out the full Lorentz ! collisions: ! ! add damping proportional to ky^4 to electron moments nu_erk=nuky*rky2(m,n)**2 ! ! GWH 9/12/94, it would seem to me that a better formula would be: ! nu_erk=nu_e/2.0/(1-kap(l))**2 ! based on Dimensional arguments and D=dx**2/dt. This would also ! give more naturally a sqrt(nu) dependence at low collisionality. ! ! ! -------------electron density equation----------------- ! e_density1(l,m,n)=e_density(l,m,n)-dt(m)*( & nl_edensity(l,m,n)+ & ii*omegade(l,m,n)*1.5* & (e_p2(l,m,n)/tiovte-phi_ba(l,m,n)) & +ii*rky(m,n)*phi_ba(l,m,n) & -nu_erk*phi_ba(l,m,n)) e_density1(l,m,n)=e_density1(l,m,n)/(1+dt(m)*nu_erk) ! ! -------------electron pressure equation-------------- ! e_p1(l,m,n)=e_p(l,m,n)-dt(m)*( & nl_ep(l,m,n)+ & ii*omegade(l,m,n)*2.5* & (e_r2(l,m,n)/tiovte-phi_ba(l,m,n)) & +ii*rky(m,n)*phi_ba(l,m,n)*(1.+etae) & -nu_erk*phi_ba(l,m,n)) e_p1(l,m,n)=e_p1(l,m,n)/(1+dt(m)*nu_erk) ! ! -------------electron r equation-------------- ! e_r1(l,m,n)=e_r(l,m,n)-dt(m)*( & nl_er(l,m,n)+ & ii*omegade(l,m,n)*3.5* & (wd_t(l,m,n)/tiovte-phi_ba(l,m,n)) & +ii*rky(m,n)*phi_ba(l,m,n)*(1.+2.*etae) & +3.5*abs(omegade(l,m,n))*awd_t(l,m,n)/tiovte & -nu_erk*phi_ba(l,m,n)) e_r1(l,m,n)=e_r1(l,m,n)/(1+dt(m)*nu_erk) enddo enddo enddo ! ! ------ nonlinear electron terms ---------------------- ! if (lin == 0) then ! dt3 = spread(spread(dt,1,ld),3,nd) ! e_density1 = e_density1 - dt3*nl_edensity ! e_p1 = e_p1 - dt3*nl_ep ! e_r1 = e_r1 - dt3*nl_er endif ! ! -------------electron t equation-------------- ! if (nemom == 4) then do n=1,nd do m=1,md do l=1,kd e_t1(l,m,n)=e_t(l,m,n)-dt(m)*( & nl_et(l,m,n)+ & ii*omegade(l,m,n)*4.5* & ((nuai*e_density2(l,m,n)+nubi*e_p2(l,m,n) & +nuci*e_r2(l,m,n)+nudi*e_t2(l,m,n))/tiovte & -phi_ba(l,m,n)) & +ii*rky(m,n)*phi_ba(l,m,n)*(1.+3.*etae) & +4.5*abs(omegade(l,m,n))* & (nuar*e_density2(l,m,n)+nubr*e_p2(l,m,n) & +nucr*e_r2(l,m,n)+nudr*e_t2(l,m,n))/tiovte) e_t1(l,m,n)=e_t1(l,m,n)/(1+dt(m)*nu_erk) enddo enddo enddo ! ! -------- nonlinear e_t term ------------------- ! ! if(lin == 0) e_t1 = e_t1 - dt3*nl_et endif ! Evaluate Collisions on trapped electrons: if(nueeff .ne. 0.0) call ecollis(e_density1,e_p1,e_r1,e_t1,dt) endif ! ! Make like ron's wedge ! goto 33 do n=1,nd do i=1,nspecies do m=1,md if ((abs(r0(m,n)) > pi+.001).or.(m == 1)) then do l=1,ld density1(l,m,n,i)=cmplx(1.e-15,1.e-15) u_par1(l,m,n,i)=cmplx(1.e-15,1.e-15) t_par1(l,m,n,i)=cmplx(1.e-15,1.e-15) q_par1(l,m,n,i)=cmplx(1.e-15,1.e-15) t_perp1(l,m,n,i)=cmplx(1.e-15,1.e-15) q_perp1(l,m,n,i)=cmplx(1.e-15,1.e-15) enddo endif enddo enddo enddo 33 continue do i=1,nspecies call bcon(density1(:,:,:,i)) call bcon(u_par1(:,:,:,i)) call bcon(t_par1(:,:,:,i)) if(nparmom >= 4) call bcon(q_par1(:,:,:,i)) call bcon(t_perp1(:,:,:,i)) if(nperpmom >= 2) call bcon(q_perp1(:,:,:,i)) enddo end subroutine timestep subroutine ecollis(e_density1,e_p1,e_r1,e_t1,dt) ! ! bounce-averaged Lorentz collision operator for trapped electrons ! gwh 9/11/94 ! ! Use implicit time-advancement, for stability and robustness, ! employing standard tri-diagonal techniques. ! ! real, dimension(:) :: dt complex, dimension(:,:,:) :: e_density1, e_p1, e_r1, e_t1 real, dimension(:, :), pointer, save :: aec, bec, cec, zbeta, zgamma ! Local Variables: real, dimension(size(dt)) :: dtnu0 real :: epsb integer :: l, m, n !************************************************************************ ! For simplicity for now, we will ignore the energy dependence of the ! collision frequency so the collision operator is almost diagonal ! (except for the phi_ba term, which is treated explicitly, which ! appears to work fine, particularly given the way Mike set up the ! Poisson equation solver in its depedence on <_b>_kappa from ! the past time step). ! ! Use the definitions: ! ! nueeff= nu_e_eff L_n/v_t, effective electron collision frequency ! MB: clarification: nueeff= nu_e_eff L_ne/v_ti ! ! nu_e_eff = nu_e/(r/R) ! = nustare * (1+1/Zeff) *3 *sqrt(pi*r/R) *Ln/R/q *sqrt(mi/me) ! MB: slight error in previous line, it's really nueeff, and missing Ti/Te ! nueeff = nustare * (1+1/Zeff) *3*sqrt(pi*r/R/2)*Lne/R/q *sqrt(mi*Te/me*Ti) ! where nustare is the usual neoclassical/snap collision ! frequency parameter (ignoring e-e collisions), ! nustare = nu_ei*Zeff/(r/R)/omega_b ! nu_ei = 1/tau_e = Braginskii's definition ! omega_b = sqrt(2*r/R*Te/me)/(q R) = bounce frequency ! MB: TRANSP,SNAP use: omega_b = sqrt(r/R*Te/me)/(q R), ! nueeff above reflects this sqrt(2) change ! epsb=epse/(1+epse) dtnu0 = dt*nueeff*epse/(8*epsb**2)*sqrt(2*pi)/16 if(.not.associated(aec)) allocate (aec(kdpass, md), bec(kdpass, md), cec(kdpass, md), & zbeta(kdpass, md), zgamma(kdpass, md)) ! if (nemom == 3) then if(mrr(1,1) .ne. 0) call aborter(6,'ecollis.f works only if the first m mode is m=0') ! solve for all (k_theta .ne. 0) modes: do n=1,nd do m=2,md do l=2,kd-1 ! ! -------------electron density equation----------------- ! e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) ! ! -------------electron pressure equation-------------- ! e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) ! ! -------------electron r equation-------------- ! e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) enddo ! handle the boundary at l=kd (everything vanishes for l>kd): l=kd e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte & *jac(l)*(diff(l)*( -phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte & *jac(l)*(diff(l)*( -phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*( -phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) ! handle the boundary at the origin at l=1 (no flux across kappa=0 ! because of symmetry...): l=1 e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) enddo enddo ! ! -------------electron t equation-------------- ! ! Collision operator hasn't been added to the t-equation yet... if (nemom == 4) then call aborter(6,"Collision operator not yet valid for nemom=4") endif !*********************************************************** ! solve for all (k_theta = 0) modes: ! (most of the below is a slightly modified copy of the section ! above for k_theta .ne. 0 modes) do n=1,nd m=1 do l=2,kdpass-1 ! ! -------------electron density equation----------------- ! e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) ! ! -------------electron pressure equation-------------- ! e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) ! ! -------------electron r equation-------------- ! e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) enddo ! handle the boundary at l=kdpass (no flux beyond max kappa): ! (These expressions could be simplified by combining into the above ! do loops if phi_ba(kdpass+1) was defined, since diff(kdpass)=0): l=kdpass e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte & *jac(l)*(diff(l)*( -phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte & *jac(l)*(diff(l)*( -phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*( -phi_ba(l,m,n)) & -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n))) ! handle the boundary at the origin at l=1 (no flux across kappa=0 ! because of symmetry...): l=1 e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) enddo !********************************************************************** ! ! Set up coefficients of the tridiagonal equation to be solved, ! of the form: ! ! a(l)*den(l-1) + b(l)*den(l) + c(l)*den(l+1) = den_old(l) ! ! ! For dt(k_theta) this gets a little more complicated: ! ! density moment ! do m=1,md do l=1,kdpass aec(l,m)=-dtnu0(m)*jac(l)*diff(l-1) cec(l,m)=-dtnu0(m)*jac(l)*diff(l) bec(l,m)=1-aec(l,m)-cec(l,m) enddo aec(1,m)=0.0 cec(kdpass,m)=0.0 enddo ! set up coefficients of the inversion formula: do m=1,md zbeta(1,m)=bec(1,m) do l=2,kdpass zgamma(l,m)=cec(l-1,m)/zbeta(l-1,m) zbeta(l,m)=bec(l,m)-aec(l,m)*zgamma(l,m) enddo enddo ! recursive-solution to tridagonal equations doesn't vectorize so ! put it as the outer loop. (When optimizing for multitasking, ! may want to multitask over m (and over different fluid moments), ! and vectorize over n.) ! ! Try to speed up by combining the various loops... do m=1,md do n=1,nd e_density1(1,m,n)= e_density1(1,m,n)/zbeta(1,m) enddo enddo ! solve for (k_theta .ne. 0) modes: do l=2,kd do m=2,md do n=1,nd e_density1(l,m,n)=(e_density1(l,m,n)-aec(l,m)*e_density1(l-1,m,n))/zbeta(l,m) enddo enddo enddo do l=kd-1,1,-1 do m=2,md do n=1,nd e_density1(l,m,n)=e_density1(l,m,n)-zgamma(l+1,m)*e_density1(l+1,m,n) enddo enddo enddo ! solve for (k_theta == 0) modes: do l=2,kdpass m=1 do n=1,nd e_density1(l,m,n)=(e_density1(l,m,n)-aec(l,m)*e_density1(l-1,m,n))/zbeta(l,m) enddo enddo do l=kdpass-1,1,-1 m=1 do n=1,nd e_density1(l,m,n)=e_density1(l,m,n)-zgamma(l+1,m)*e_density1(l+1,m,n) enddo enddo ! ! p moment ! do m=1,md do l=1,kdpass aec(l,m)=-dtnu0(m)/2*jac(l)*diff(l-1) cec(l,m)=-dtnu0(m)/2*jac(l)*diff(l) bec(l,m)=1-aec(l,m)-cec(l,m) enddo aec(1,m)=0.0 cec(kdpass,m)=0.0 enddo ! set up coefficients of the inversion formula: do m=1,md zbeta(1,m)=bec(1,m) do l=2,kdpass zgamma(l,m)=cec(l-1,m)/zbeta(l-1,m) zbeta(l,m)=bec(l,m)-aec(l,m)*zgamma(l,m) enddo enddo ! recursive-solution to tridagonal equations doesn't vectorize so ! put it as the outer loop. (When optimizing for multitasking, ! may want to multitask over m (and over different fluid moments), ! and vectorize over n.) ! ! Try to speed up by combining the various loops... do m=1,md do n=1,nd e_p1(1,m,n)= e_p1(1,m,n)/zbeta(1,m) enddo enddo ! solve for (k_theta .ne. 0) modes: do l=2,kd do m=2,md do n=1,nd e_p1(l,m,n)=(e_p1(l,m,n)-aec(l,m)*e_p1(l-1,m,n))/zbeta(l,m) enddo enddo enddo do l=kd-1,1,-1 do m=2,md do n=1,nd e_p1(l,m,n)=e_p1(l,m,n)-zgamma(l+1,m)*e_p1(l+1,m,n) enddo enddo enddo ! solve for (k_theta == 0) modes: do l=2,kdpass m=1 do n=1,nd e_p1(l,m,n)=(e_p1(l,m,n)-aec(l,m)*e_p1(l-1,m,n))/zbeta(l,m) enddo enddo do l=kdpass-1,1,-1 m=1 do n=1,nd e_p1(l,m,n)=e_p1(l,m,n)-zgamma(l+1,m)*e_p1(l+1,m,n) enddo enddo ! ! r moment ! do m=1,md do l=1,kdpass aec(l,m)=-dtnu0(m)*0.3125*jac(l)*diff(l-1) cec(l,m)=-dtnu0(m)*0.3125*jac(l)*diff(l) bec(l,m)=1-aec(l,m)-cec(l,m) enddo aec(1,m)=0.0 cec(kdpass,m)=0.0 enddo ! set up coefficients of the inversion formula: do m=1,md zbeta(1,m)=bec(1,m) do l=2,kdpass zgamma(l,m)=cec(l-1,m)/zbeta(l-1,m) zbeta(l,m)=bec(l,m)-aec(l,m)*zgamma(l,m) enddo enddo ! recursive-solution to tridagonal equations doesn't vectorize so ! put it as the outer loop. (When optimizing for multitasking, ! may want to multitask over m (and over different fluid moments), ! and vectorize over n.) ! ! Try to speed up by combining the various loops... do m=1,md do n=1,nd e_r1(1,m,n)= e_r1(1,m,n)/zbeta(1,m) enddo enddo ! solve for (k_theta .ne. 0) modes: do l=2,kd do m=2,md do n=1,nd e_r1(l,m,n)=(e_r1(l,m,n)-aec(l,m)*e_r1(l-1,m,n))/zbeta(l,m) enddo enddo enddo do l=kd-1,1,-1 do m=2,md do n=1,nd e_r1(l,m,n)=e_r1(l,m,n)-zgamma(l+1,m)*e_r1(l+1,m,n) enddo enddo enddo ! solve for (k_theta == 0) modes: do l=2,kdpass m=1 do n=1,nd e_r1(l,m,n)=(e_r1(l,m,n)-aec(l,m)*e_r1(l-1,m,n))/zbeta(l,m) enddo enddo do l=kdpass-1,1,-1 m=1 do n=1,nd e_r1(l,m,n)=e_r1(l,m,n)-zgamma(l+1,m)*e_r1(l+1,m,n) enddo enddo end subroutine ecollis end module step