module step ! ! (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. ! 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,& n_e1,u_e1,tpar_e1,tperp_e1,apar1, & n_e2,u_e2,tpar_e2,tperp_e2,apar2,phi2,delt) ! 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 itg_data use gryffin_layouts use linear use nonlinear use constants, only: ii, pi use bounds, only: bcon use fields use gryffin_grid, only: ld, ldb, kd, kdpass, md, l_left, l_right, & x0, r0, rkx, rkperp2, rky, rky2 complex, dimension(:,m_low:,n_low:,:) :: & density1, u_par1, t_par1, q_par1, & t_perp1, q_perp1, density2, u_par2, t_par2, q_par2, t_perp2, q_perp2 complex, dimension(:,m_low:,n_low:) :: & e_density1, e_p1, e_r1, e_t1, e_density2, e_p2, e_r2, e_t2, & n_e1, u_e1, tpar_e1, tperp_e1, apar1, & phi2, n_e2, u_e2, apar2, tpar_e2, tperp_e2 complex, dimension(ld, m_low:m_alloc, & n_low:n_alloc, nspecies) :: & kp_density2, kp_upar2, kp_tpar2, kp_qpar2, kp_tperp2, kp_qperp2, & kp_vphi, akp_qpar2, akp_qperp2, qparps, qperpps, & apar1_flr1, J0apar1 complex, dimension(ld, m_low:m_alloc, & n_low:n_alloc) :: kp_n_e12, kp_u_e2, akp_u_e2, kp2_n_e2, & ubar, teclose, dapardt1 complex, dimension(kdpass, m_low:m_alloc, & n_low:n_alloc) :: wd_t, awd_t real, dimension(m_low:) :: delt real nuur,nuui,vmom,nu_erk,nuky,znorm,dz0, omegad, a2, nu_ee real, dimension(md) :: switch integer i,l,m,n real, dimension(size(density1,1)) :: dzt complex :: tparfsa, tperpfsa logical skip_code !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Note: electron terms in subroutine 'finite_beta_terms' must be ! calculated first if beta_e > 0. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! electron dynamics for finite beta case ! if (beta_e > 0.0) call finite_beta_terms ! subroutine found below call ion_terms ! subroutine found below ! ! electron dynamics for finite beta case ! if (beta_e > 0.0) then ! zero out ubar ubar=0. do n = n_low, n_high do m = m_low, m_high do l=1,ld-1 ! ! ---------EM Ampere's Law (determines u_e)------------------ ! for now, use weights from poisson nbar for ampere ubar 4/97 ! do i=1,nspecies ubar(l,m,n)=ubar(l,m,n) & +(u_par1(l,m,n,i)*nwgt(l,m,n,i) & +q_perp1(l,m,n,i)*twgt(l,m,n,i))*g0wgtp(l,m,n,i) enddo u_e1(l,m,n)=ubar(l,m,n) & -rkperp2(l,m,n)*apar1(l,m,n)/(tiovte*beta_e) enddo enddo enddo endif ! ! bounce averaged electrons: ! if (epse > 0.0) call bounce_terms ! subroutine found below ! ! Make like ron's wedge ! skip_code=.true. if(skip_code) goto 33 do i=1,nspecies do n = n_low, n_high do m = m_low, m_high if ((abs(r0(m,n)) > pi+.001).or.(m == 1)) then density1(:,m,n,i)=cmplx(1.e-15,1.e-15) u_par1(:,m,n,i)=cmplx(1.e-15,1.e-15) t_par1(:,m,n,i)=cmplx(1.e-15,1.e-15) q_par1(:,m,n,i)=cmplx(1.e-15,1.e-15) t_perp1(:,m,n,i)=cmplx(1.e-15,1.e-15) q_perp1(:,m,n,i)=cmplx(1.e-15,1.e-15) 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 !--- EM electron boundary conditions ! probably not necessary for derived quantities like u_e, tpar_e, tperp_e if (beta_e > 0.) then call bcon(n_e1) call bcon(u_e1) call bcon(tpar_e1) call bcon(tperp_e1) call bcon(apar1) endif contains subroutine ion_terms 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 ! ! load kp_upar2,kp_vphi,etc. ! do l=1,ld-1 kp_density2(l,:,:,:)=density2(l,:,:,:)*bmaginv(l) kp_upar2(l,:,:,:)=u_par2(l,:,:,:)*bmaginv(l) kp_tpar2(l,:,:,:)=t_par2(l,:,:,:)*bmaginv(l) kp_vphi(l,:,:,:)=phi_u(l,:,:,:) enddo kp_density2(ld,:,:,:)=0. kp_upar2(ld,:,:,:)=0. kp_tpar2(ld,:,:,:)=0. kp_vphi(ld,:,:,:)=0. call kparfft(kp_density2, 0, 0.) call kparfft(kp_upar2, 0, 0.) call kparfft(kp_tpar2, 0, 0.) call kparfft(kp_vphi, 0, 0.) do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld omegad = omega_gb(l,m,n,i)+omega_kap(l,m,n,i) ! ! -----------------density equation---------------------- ! density1(l,m,n,i)=density(l,m,n,i)-delt(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*(phi_nd(l,m,n,i)+2.*density2(l,m,n,i)) & -ii*omegad*(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)-delt(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*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*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 if (beta_e > 0.) then do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld density1(l,m,n,i) = density1(l,m,n,i) + delt(m)*emnl_density1(l,m,n,i) u_par1(l,m,n,i)=u_par1(l,m,n,i)-delt(m)*( & -emnl_upar1(l,m,n,i)-emnl_upar2(l,m,n,i) & -ii*rky(m,n)/Ln(i)*vt(i)*apar_u(l,m,n,i) & +(J0apar1(l,m,n,i)-J0apar(l,m,n,i))/delt(m)*vt(i)*charge(i)/tau(i) & ) enddo enddo enddo enddo endif ! ! ---------three pole closure: parallel moments---------- ! if(nparmom == 3) q_par2 = rkkin * t_par2 ! ! load kp_qpar2 ! do l=1,ld-1 kp_qpar2(l,:,:,:)=q_par2(l,:,:,:)*bmaginv(l) enddo kp_qpar2(ld,:,:,:)=0.0 if (nparmom == 3) then call kparfft(kp_qpar2, 2, 0.) else call kparfft(kp_qpar2, 0, 0.) endif do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld omegad = omega_gb(l,m,n,i)+omega_kap(l,m,n,i) ! ! -------------parallel temperature equation----------- ! t_par1(l,m,n,i)=t_par(l,m,n,i)-delt(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*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) & ! CVS +(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) & ! CVS +(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*switch(m))*t_par2(l,m,n,i) & +2.*nu2r*abs(omegad*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 if (beta_e > 0.) then do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld t_par1(l,m,n,i) = t_par1(l,m,n,i) + delt(m)*emnl_tpar(l,m,n,i) enddo enddo enddo enddo endif if(nparmom == 4) then ! ! load kp_tpar2,akp_qpar2 ! qparps=0. qperpps=0. if (iflow == 2) then m = 1 do i=1,nspecies do n = n_low, n_high if(idx_local(m_lo, m, n)) then tparfsa=0.0 tperpfsa=0.0 do l=l_left,l_right tparfsa=tparfsa + t_par2(l,m,n,i)*dzt(l) tperpfsa=tperpfsa + t_perp2(l,m,n,i)*dzt(l) enddo if(igeo == 0) then qparps (:,m,n,i)=3*ii*rkx(:,m,n)*qsf/eps*bmaginv(:)*tparfsa qperpps(:,m,n,i)= ii*rkx(:,m,n)*qsf/eps*bmaginv(:)*tperpfsa else qparps (:,m,n,i)=3*ii*rkx(:,m,n)/eps_over_q*bmaginv(:)*tparfsa qperpps(:,m,n,i)= ii*rkx(:,m,n)/eps_over_q*bmaginv(:)*tperpfsa endif endif 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, 0.) call kparfft(akp_qpar2, 1, 0.) do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld omegad = omega_gb(l,m,n,i)+omega_kap(l,m,n,i) ! ! --------------parallel heat flux equation--------------- ! q_par1(l,m,n,i)=q_par(l,m,n,i)-delt(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* & ( (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*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 if (beta_e > 0.) then do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld q_par1(l,m,n,i) = q_par1(l,m,n,i) - delt(m)*( & -3.*ii*rky(m,n)/Ln(i)*vt(i)*apar_qpar(l,m,n,i)) enddo enddo enddo enddo endif endif ! ! ----------one pole closure: perpendicular moments------------ ! if(nperpmom == 1) then do i=1,nspecies q_perp2(:,:,:,i)= rkkperp*(t_perp2(:,:,:,i) & +charge(i)/tau(i)*phi_qperp(:,:,:,i)) enddo endif ! ! load kp_qperp2 ! do l=1,ld-1 kp_qperp2(l,:,:,:)=q_perp2(l,:,:,:)*bmaginv(l)**2 enddo kp_qperp2(ld,:,:,:)=0.0 if (nperpmom == 1) then call kparfft(kp_qperp2, 2, 0.) else call kparfft(kp_qperp2, 0, 0.) endif do i=1,nspecies do n=n_low, n_high do m=m_low, m_high do l=1,ld omegad = omega_gb(l,m,n,i)+omega_kap(l,m,n,i) ! ! ----------perpendicular temperature equation-------------- ! t_perp1(l,m,n,i)=t_perp(l,m,n,i)-delt(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* & ( 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) & ! CVS +(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*switch(m))*t_perp2(l,m,n,i) & +2.*nu3r*abs(omegad*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 ! BD 8.13.99 bug fix for nparmom = 3 (rarely used) if(nperpmom == 1) then q_perp2 = 0. end if if (beta_e > 0.) then do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld t_perp1(l,m,n,i) = t_perp1(l,m,n,i) - delt(m)*( & -emnl_tperp(l,m,n,i)) enddo enddo enddo enddo endif if(nperpmom >= 2) then ! ! load kp_tperp2,akp_qperp2 ! akp_qperp2 = (q_perp2-qperpps)*rkkperp2 akp_qperp2(ld,:,:,:) = 0. do i=1,nspecies kp_tperp2(:,:,:,i)=t_perp2(:,:,:,i)+charge(i)/tau(i)*phi_qperp(:,:,:,i) enddo kp_tperp2(ld,:,:,:)=0.0 call kparfft(akp_qperp2, 1, 0.) call kparfft(kp_tperp2, 0, 0.) do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld omegad = omega_gb(l,m,n,i)+omega_kap(l,m,n,i) ! ! -----------perpendicular heat flux equation---------------- ! q_perp1(l,m,n,i)=q_perp(l,m,n,i)-delt(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* & ( (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*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 if (beta_e > 0.) then do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld omegad = omega_gb(l,m,n,i)+omega_kap(l,m,n,i) q_perp1(l,m,n,i) = q_perp1(l,m,n,i) - delt(m)*( & -ii*rky(m,n)/Ln(i)*vt(i)*apar_qperp(l,m,n,i) & -ii*omegad*apar_flr1(l,m,n,i) & +(apar1_flr1(l,m,n,i)-apar_flr1(l,m,n,i))/delt(m) & *vt(i)*charge(i)/tau(i) & ) enddo enddo enddo enddo endif endif end subroutine ion_terms subroutine bounce_terms 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 = n_low, n_high do m = m_low, m_high 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)-delt(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+delt(m)*nu_erk) ! ! -------------electron pressure equation-------------- ! e_p1(l,m,n)=e_p(l,m,n)-delt(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+delt(m)*nu_erk) ! ! -------------electron r equation-------------- ! e_r1(l,m,n)=e_r(l,m,n)-delt(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+delt(m)*nu_erk) enddo enddo enddo ! ! -------------electron t equation-------------- ! if (nemom == 4) then do n=n_low, n_high do m=m_low, m_high do l=1,kd e_t1(l,m,n)=e_t(l,m,n)-delt(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+delt(m)*nu_erk) enddo enddo enddo endif ! Evaluate Collisions on trapped electrons: if(nueeff .ne. 0.0) call ecollis(e_density1,e_p1,e_r1,e_t1,delt) end subroutine bounce_terms subroutine finite_beta_terms ! initializations (probably unnecessary) a2=0. nu_ee=0. kp_u_e2 = 0. kp2_n_e2 = 0. teclose = 0. do n = n_low, n_high do m = m_low, m_high ! load kp_u_e (to be kpar*ue ) ! and kp2_n_e2 (to be kpar**2 * n_e) do l=1,ld-1 kp_u_e2(l,m,n)=u_e2(l,m,n) kp2_n_e2(l,m,n)=n_e2(l,m,n) ! T_e closure, contains linear term only for now (can add nonlinear term ! from previous timestep) teclose(l,m,n)=rky(m,n)*etae*apar2(l,m,n)*cs enddo teclose(ld,m,n)=0. kp_u_e2(ld,m,n)=0. kp2_n_e2(ld,m,n)=0. enddo enddo call kparffte(kp_u_e2, 0, 0.) call kparffte(teclose, 7, 1.) if (semi_imp > 0) call kparffte(kp2_n_e2, 6, 1.) do n = n_low, n_high do m = m_low, m_high do l=1,ld ! ! ---------EM electron density equation---------------------- ! ! calculate constant for semi-implicit method (usually only used for slab) if (semi_imp > 0) then if (pfilter2(l,m,n) == 0.) then a2=1./(tiovte*beta_e) else a2=rkperp2(l,m,n)* & (1./pfilter2(l,m,n)+1./tiovte)/ & (tiovte*beta_e) endif endif n_e1(l,m,n)=n_e(l,m,n) & -.5*dt0**2*(a2-a0**2)*kp2_n_e2(l,m,n) & -delt(m)*( & +emnl_n_e1(l,m,n)-emnl_n_e2(l,m,n) & +ii*kp_u_e2(l,m,n)*cs & +ii*rky(m,n)*phi2(l,m,n) & -u_e2(l,m,n)*cs*bgrad(l) & -ii*(omega_gb_e(l,m,n)+omega_kap_e(l,m,n))*( & (-2.)*tiovte*phi2(l,m,n) & +2.*n_e2(l,m,n) & +2.*teclose(l,m,n) & ! +tpar_e2(l,m,n) & ! +tperp_e2(l,m,n) & ) & ) enddo enddo enddo ! for now dividing whole RHS by the semi-implicit factor. pbs 8/97 if (a0 > 0.) call kparffte(n_e1, 3, .5*dt0**2*a0**2) ! now need fft half at old and half at new timestep for semi-implicit kp_n_e12 = 0. do n = n_low, n_high do m = m_low, m_high do l=1,ld-1 if(pfilter2(l,m,n) /= 0.) & kp_n_e12(l,m,n)=.5*(1./pfilter2(l,m,n)+1./tiovte)* & (n_e2(l,m,n)+n_e1(l,m,n)) & -(n_e2(l,m,n)/pfilter2(l,m,n)+phi2(l,m,n)) enddo kp_n_e12(ld,m,n)=0. enddo enddo call kparffte(kp_n_e12, 0, 0.) ! PS 10/98 calculate 2-moment Hammett-Perkins Landau damping term if ! meovmi>0 if (meovmi > 0.) then do n = n_low, n_high do m = m_low, m_high ! load akp_u_e2 (to be |kpar|*ue ) do l=1,ld-1 akp_u_e2(l,m,n)=u_e2(l,m,n) enddo akp_u_e2(ld,m,n)=0. enddo enddo call kparffte(akp_u_e2, 1, 0.) ! gives |kpar|*u_e2 else akp_u_e2 = 0. endif do n = n_low, n_high do m = m_low, m_high do l=1,ld ! ! ---------EM momentum eqn (apar evolution)------------------ ! dapardt1(l,m,n)= & emnl_apar(l,m,n) & +ii*kp_n_e12(l,m,n) & -ii*rky(m,n)*apar2(l,m,n)/tiovte & +nuei*(-1.*rkperp2(l,m,n)*apar2(l,m,n))/(tiovte*beta_e) & +sqrt(0.5*pi*meovmi/tiovte)*akp_u_e2(l,m,n) ! -(tpar_e2(l,m,n)-tperp_e2(l,m,n))* & ! bgrad(l)/tiovte if(pfilter2(l,m,n) == 0.) then apar1(l,m,n)=0. dapardt1(l,m,n)=0. ! apar1(l,m,n)=apar(l,m,n)-delt(m)*( & ! -(.5/tiovte)*ii*kp_n_e12(l,m,n) & ! +ii*rky(m,n)*apar2(l,m,n)/tiovte ) else apar1(l,m,n)=apar(l,m,n)-delt(m)*( & -dapardt1(l,m,n) & ) endif ! ! ---------EM tpar_e equation-------------------------------- ! ! tpar_e1(l,m,n)=tpar_e(l,m,n)-delt(m)*( & ! 2.*ii*kp_u_e2(l,m,n)*cs & ! +ii*rky(m,n)*etae*phi2(l,m,n) & ! -ii*omegad_e(l,m,n)*(6.*tpar_e2(l,m,n) & ! +2.*n_e2(l,m,n) & ! -2.*tiovte*phi2(l,m,n) ) & ! +(2./3.)*nu_ee*(tpar_e2(l,m,n)-tperp_e2(l,m,n)) ) tpar_e1(l,m,n)=teclose(l,m,n) ! ! ---------EM tperp_e equation-------------------------------- ! ! tperp_e1(l,m,n)=tperp_e(l,m,n)-delt(m)*( & ! -u_e2(l,m,n)*bgrad(l)*cs & ! +ii*rky(m,n)*etae*phi2(l,m,n) & ! -ii*omegad_e(l,m,n)*(4.*tperp_e2(l,m,n) & ! +n_e2(l,m,n) & ! -tiovte*phi2(l,m,n) ) & ! +nu_ee/3.*(tperp_e2(l,m,n)-tpar_e2(l,m,n)) ) tperp_e1(l,m,n)=teclose(l,m,n) enddo enddo enddo ! Calculate flr on new apar1 needed for d/dt in ion equations J0apar1 = spread(apar1,4,nspecies) * g0wgt apar1_flr1 = spread(apar1,4,nspecies) * flrwgt end subroutine finite_beta_terms end subroutine timestep subroutine ecollis(e_density1, e_p1, e_r1, e_t1, delt) use itg_data, only: epse, nueeff, tiovte, nemom use gryffin_layouts use gryffin_grid, only: kdpass, mrr, kd, nd use fields, only: phi_ba, e_density, e_p, e_r use bounce_avg, only: jac, diff use constants, only: pi ! ! 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(m_low:) :: delt complex, dimension(:,m_low:,n_low:) :: & e_density1, e_p1, e_r1, e_t1 real, dimension(kdpass, m_low:m_alloc) :: aec, bec, cec, & zbeta, zgamma ! Local Variables: real, dimension(m_low:m_alloc) :: 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 = delt*nueeff*epse/(8*epsb**2)*sqrt(2*pi)/16 ! if (nemom == 3) then if(idx_local(m_lo, 1, 1)) then if(mrr(1,1) .ne. 0) call aborter(6,'ecollis.f works only if the first m mode is m=0') endif ! solve for all (k_theta .ne. 0) modes: do n = n_low, n_high do m = m_low, m_high if(m == 1) cycle 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) m=1 do n = n_low, n_high if(idx_local(m_lo, m, n)) then do l=2,kdpass-1 ! ! -------------electron density equation----------------- ! e_density1(l,m,n)=e_density(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_p(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_r(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_density(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_p(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_r(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_density(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_p(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_r(l,m,n)-dtnu0(m)*0.3125*tiovte & *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))) endif 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=m_low, m_high 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=m_low, m_high 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=m_low, m_high do n = n_low, n_high 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=m_low, m_high if(m == 1) cycle do n=n_low, n_high 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=m_low, m_high if (m == 1) cycle do n=n_low, n_high 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: m=1 do l=2,kdpass do n=n_low, n_high if(idx_local(m_lo, m, n)) e_density1(l,m,n) = & (e_density1(l,m,n)-aec(l,m)*e_density1(l-1,m,n))/zbeta(l,m) enddo enddo m=1 do l=kdpass-1,1,-1 do n = n_low, n_high if(idx_local(m_lo, m, n)) 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 = m_low, m_high 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=m_low, m_high 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 = m_low, m_high do n = n_low, n_high 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=m_low, m_high if(m == 1) cycle do n=n_low, n_high 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=m_low, m_high if (m == 1) cycle do n = n_low, n_high 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: m=1 do l=2,kdpass do n = n_low, n_high if(idx_local(m_lo, m, n)) e_p1(l,m,n) = & (e_p1(l,m,n)-aec(l,m)*e_p1(l-1,m,n))/zbeta(l,m) enddo enddo m=1 do l=kdpass-1,1,-1 do n=1,nd if(idx_local(m_lo, m, n)) 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=m_low, m_high 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=m_low, m_high 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=m_low, m_high do n=n_low, n_high 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=m_low, m_high if(m == 1) cycle do n=n_low, n_high 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=m_low, m_high if (m == 1) cycle do n=n_low, n_high 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: m=1 do l=2,kdpass do n=n_low, n_high if(idx_local(m_lo, m, n)) e_r1(l,m,n) = & (e_r1(l,m,n)-aec(l,m)*e_r1(l-1,m,n))/zbeta(l,m) enddo enddo m=1 do l=kdpass-1,1,-1 do n=n_low, n_high if(idx_local(m_lo, m, n)) 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 subroutine courant use itg_data, only: beta_e, dtadj, dt0, dt, nspecies use fields, only: apar, potential use mp, only: proc0 use gryffin_layouts use linear, only: phi_u use nlps, only: maxv real :: vmax1, vmax2, vmax3 integer :: i, m ! ! ******************************************************************* ! Check Courant condition: ! vmax1 = 0. vmax2 = 0. vmax3 = 0. if (beta_e > 0.) then call maxv(apar,vmax2) call maxv(potential,vmax3) endif do i=1,nspecies call maxv(phi_u(:,:,:,i),vmax1) ! ! With our FLR models, the nonlinear FLR terms are always smaller than ! the main term, so the next two lines are unnecessary. This may not be ! true for such models as the Taylor series FLR with nonlinear FLR. ! ! call maxv(phi_flr1(:,:,:,i),vmax1) ! call maxv(phi_flr2(:,:,:,i),vmax1) enddo vmax1=max(dtadj/dt0,vmax1) if (vmax2 > vmax1) then if(proc0) write(*,*) 'Apar nonlinearity triggering courant ',vmax2,dtadj/vmax2 vmax1=vmax2 endif if (vmax3 > vmax1) then if(proc0) write(*,*) 'Phi nonlinearity triggering courant ',vmax3,dtadj/vmax3 vmax1=vmax3 endif do m=m_low, m_high dt(m)=min(dt0,dtadj/vmax1) enddo end subroutine courant end module step