subroutine timestep(wc1,ws1,vc1,vs1,tparc1,tpars1,qparc1, * qpars1,tc1,ts1,qc1,qs1,rc1,rs1,sc1,ss1,wc2,ws2,vc2,vs2, * tparc2,tpars2,prc2,prs2,qparc2,qpars2,tc2,ts2,qc2,qs2,rc2, * rs2,sc2,ss2,dt,unc,uns,uvc,uvs,utparc,utpars,utc,uts,uqc, * uqs,urc,urs) c c mbeer -- switched to Cowley coordinates, May 8, 1992 c c $Header: /afs/pppl.gov/u/qpliu/cvsroot/codes/gryffin/src/old/timesteprw.f,v 4.0 1995/05/24 04:46:58 mbeer Exp $ c c Revision 3.0 1994/02/03 22:17:18 bdorland c No change from version 2.0. c c Advance equations one time step, not including collisional c dissipation c implicit none include 'itg.par' include 'itg.cmn' real ws1(lz,mz,nz,nspecz),wc1(lz,mz,nz,nspecz), * vs1(lz,mz,nz,nspecz),vc1(lz,mz,nz,nspecz), * tpars1(lz,mz,nz,nspecz),tparc1(lz,mz,nz,nspecz), * ts1(lz,mz,nz,nspecz),tc1(lz,mz,nz,nspecz), * qs1(lz,mz,nz,nspecz),qs2(lz,mz,nz,nspecz), * qc1(lz,mz,nz,nspecz),qc2(lz,mz,nz,nspecz), * qpars1(lz,mz,nz,nspecz),qparc1(lz,mz,nz,nspecz) real ws2(lz,mz,nz,nspecz),wc2(lz,mz,nz,nspecz), * vs2(lz,mz,nz,nspecz),vc2(lz,mz,nz,nspecz), * tpars2(lz,mz,nz,nspecz),tparc2(lz,mz,nz,nspecz), * ts2(lz,mz,nz,nspecz),tc2(lz,mz,nz,nspecz),dt, * uqs(lz,mz,nz,nspecz),uqc(lz,mz,nz,nspecz), * uns(lz,mz,nz,nspecz),unc(lz,mz,nz,nspecz) real uvs(lz,mz,nz,nspecz),uvc(lz,mz,nz,nspecz), * utpars(lz,mz,nz,nspecz),utparc(lz,mz,nz,nspecz), * uts(lz,mz,nz,nspecz),utc(lz,mz,nz,nspecz), * ky,kz,qpars2(lz,mz,nz,nspecz),qparc2(lz,mz,nz,nspecz), * prs2(lz,mz,nz,nspecz),prc2(lz,mz,nz,nspecz) real rc1(lz,mz,nz,nspecz),rs1(lz,mz,nz,nspecz), * sc1(lz,mz,nz,nspecz),ss1(lz,mz,nz,nspecz) real rc2(lz,mz,nz,nspecz),rs2(lz,mz,nz,nspecz), * sc2(lz,mz,nz,nspecz),ss2(lz,mz,nz,nspecz) real urc(lz,mz,nz,nspecz),urs(lz,mz,nz,nspecz) real kpart,kdopp,nuur,nuui,vmom real kpvc2(lz,mz,nz,nspecz),kpvs2(lz,mz,nz,nspecz) real kpuvc(lz,mz,nz,nspecz),kpuvs(lz,mz,nz,nspecz) real kpprc(lz,mz,nz,nspecz),kpprs(lz,mz,nz,nspecz) real kpqparc2(lz,mz,nz,nspecz),kpqpars2(lz,mz,nz,nspecz) real akpqparc2(lz,mz,nz,nspecz),akpqpars2(lz,mz,nz,nspecz) real kptparc2(lz,mz,nz,nspecz),kptpars2(lz,mz,nz,nspecz) real kpqc2(lz,mz,nz,nspecz),kpqs2(lz,mz,nz,nspecz) real akpqc2(lz,mz,nz,nspecz),akpqs2(lz,mz,nz,nspecz) real kptc2(lz,mz,nz,nspecz),kpts2(lz,mz,nz,nspecz) integer l,m,n,mr c if (nparmom.eq.3) then nuur=nu5r nuui=nu5i vmom=0. else nuur=0.0 nuui=0.0 vmom=1. endif c c cleaning up. at some time, may want to use kdopp. mbeer c kz=0.0 kpart=0.0 kdopp=0.0 do n=1,nd do m=1,md c c load kpvc2,kpvs2,kpuvc,kpuvs c do l=1,ld-1 kpvc2(l,m,n,1)=vc2(l,m,n,1)*bmaginv(l) kpvs2(l,m,n,1)=vs2(l,m,n,1)*bmaginv(l) kpuvc(l,m,n,1)=uvc(l,m,n,1) kpuvs(l,m,n,1)=uvs(l,m,n,1) kpprc(l,m,n,1)=prc2(l,m,n,1)*bmaginv(l) kpprs(l,m,n,1)=prs2(l,m,n,1)*bmaginv(l) enddo kpvc2(ld,m,n,1)=0.0 kpvs2(ld,m,n,1)=0.0 kpuvc(ld,m,n,1)=0.0 kpuvs(ld,m,n,1)=0.0 kpprc(ld,m,n,1)=0.0 kpprs(ld,m,n,1)=0.0 enddo enddo call kparfft(kpvc2,kpvs2,0) call kparfft(kpuvc,kpuvs,0) call kparfft(kpprc,kpprs,0) c do n=1,nd do m=1,md do l=1,ld c c -----------------density equation---------------------- c ws1(l,m,n,1)=ws(l,m,n,1)+dt*(wpns(l,m,n,1)+wgts(l,m,n,1) cpm * +lambda2/4.*pmptns(l,m,n,1)+nu2/4.*pmt5s(l,m,n,1) * +rky(m,n)*unc(l,m,n,1)-kpvs2(l,m,n,1)*bmag(l)+kdopp*wc2(l,m,n,1) * +omegad(l,m,n)*( uqc(l,m,n,1) * +(prc2(l,m,n,1)+tc2(l,m,n,1)+wc2(l,m,n,1))/2.)) c wc1(l,m,n,1)=wc(l,m,n,1)+dt*(wpnc(l,m,n,1)+wgtc(l,m,n,1) cpm * +lambda2/4.*pmptnc(l,m,n,1)+nu2/4.*pmt5c(l,m,n,1) * -rky(m,n)*uns(l,m,n,1)-kpvc2(l,m,n,1)*bmag(l)-kdopp*ws2(l,m,n,1) * -omegad(l,m,n)*( uqs(l,m,n,1) * +(prs2(l,m,n,1)+ts2(l,m,n,1)+ws2(l,m,n,1))/2.)) c c -------------parallel velocity equation---------------- c vs1(l,m,n,1)=vs(l,m,n,1)+dt*(wpvs(l,m,n,1)+wgqs(l,m,n,1) cpm * -lambda1*pmpvs(l,m,n,1)+nu1*wgvs(l,m,n,1) * -kpuvs(l,m,n,1)-kpprs(l,m,n,1)*bmag(l) * -rky(m,n)*vz*uvc(l,m,n,1)+kdopp*vc2(l,m,n,1) * +omegad(l,m,n)*((2.+nuui)*vc2(l,m,n,1) * +0.5*vmom*(qparc2(l,m,n,1)+qc2(l,m,n,1)) ) * -nuur*abs(omegad(l,m,n))*vs2(l,m,n,1) * -(ts2(l,m,n,1)+ws2(l,m,n,1)) *bgrad(l) ) c added dissipative and reactive omegad piece c vc1(l,m,n,1)=vc(l,m,n,1)+dt*(wpvc(l,m,n,1)+wgqc(l,m,n,1) cpm * -lambda1*pmpvc(l,m,n,1)+nu1*wgvc(l,m,n,1) * -kpuvc(l,m,n,1)-kpprc(l,m,n,1)*bmag(l) * +rky(m,n)*vz*uvs(l,m,n,1)-kdopp*vs2(l,m,n,1) * -omegad(l,m,n)*((2.+nuui)*vs2(l,m,n,1) * +0.5*vmom*(qpars2(l,m,n,1)+qs2(l,m,n,1)) ) * -nuur*abs(omegad(l,m,n))*vc2(l,m,n,1) * -(tc2(l,m,n,1)+wc2(l,m,n,1)) *bgrad(l) ) c added dissipative and reactive omegad piece enddo enddo enddo c call bcon(wc1,ws1,1,0) call bcon(vc1,vs1,1,0) c if(nparmom.eq.3) then do n=1,nd do m=1,md do l=1,ld c c ---------three pole closure: parallel moments---------- c qpars2(l,m,n,1)=rkkin*tpars2(l,m,n,1) qparc2(l,m,n,1)=rkkin*tparc2(l,m,n,1) enddo enddo enddo endif c do n=1,nd do m=1,md c c load kpqparc2,kpqpars2 c do l=1,ld-1 kpqparc2(l,m,n,1)=qparc2(l,m,n,1)*bmaginv(l) kpqpars2(l,m,n,1)=qpars2(l,m,n,1)*bmaginv(l) enddo kpqparc2(ld,m,n,1)=0.0 kpqpars2(ld,m,n,1)=0.0 enddo enddo if (nparmom.eq.3) then call kparfft(kpqparc2,kpqpars2,1) else call kparfft(kpqparc2,kpqpars2,0) endif c do n=1,nd do m=1,md do l=1,ld c c -------------parallel temperature equation----------- c tpars1(l,m,n,1)=tpars(l,m,n,1)+dt*(wptps(l,m,n,1) cpm * -lambda1*pmptps(l,m,n,1)+nu1*wgtps(l,m,n,1) * -2.*kpvs2(l,m,n,1)*bmag(l)+rky(m,n)*utparc(l,m,n,1) * -kpqpars2(l,m,n,1)*bmag(l)+kdopp*tparc2(l,m,n,1) * +omegad(l,m,n)*(1.5*prc2(l,m,n,1) * +(1.5+nu1i)*tparc2(l,m,n,1) * +(0.5+nu2i)*tc2(l,m,n,1) * -0.5*(tc2(l,m,n,1)+wc2(l,m,n,1)) * +uqc(l,m,n,1) ) * -nu1r*abs(omegad(l,m,n))*tpars2(l,m,n,1) * -nu2r*abs(omegad(l,m,n))*ts2(l,m,n,1) * -(2.*vs2(l,m,n,1)+2.*qs2(l,m,n,1))*bgrad(l) * -2./3.*nuii*(tpars2(l,m,n,1)-ts2(l,m,n,1)) ) c tparc1(l,m,n,1)=tparc(l,m,n,1)+dt*(wptpc(l,m,n,1) cpm * -lambda1*pmptpc(l,m,n,1)+nu1*wgtpc(l,m,n,1) * -(2.*kpvc2(l,m,n,1)*bmag(l)+rky(m,n)*utpars(l,m,n,1)) * -kpqparc2(l,m,n,1)*bmag(l)-kdopp*tpars2(l,m,n,1) * -omegad(l,m,n)*(1.5*prs2(l,m,n,1) * +(1.5+nu1i)*tpars2(l,m,n,1) * +(0.5+nu2i)*ts2(l,m,n,1) * -0.5*(ts2(l,m,n,1)+ws2(l,m,n,1)) * +uqs(l,m,n,1) ) * -nu1r*abs(omegad(l,m,n))*tparc2(l,m,n,1) * -nu2r*abs(omegad(l,m,n))*tc2(l,m,n,1) * -(2.*vc2(l,m,n,1)+2.*qc2(l,m,n,1))*bgrad(l) * -2./3.*nuii*(tparc2(l,m,n,1)-tc2(l,m,n,1)) ) c enddo enddo enddo c call bcon(tparc1,tpars1,1,0) c if(nparmom.eq.4) then do n=1,nd do m=1,md c c load kptparc2,kptpars2,akpqparc2,akpqpars2 c do l=1,ld-1 kptparc2(l,m,n,1)=tparc2(l,m,n,1) kptpars2(l,m,n,1)=tpars2(l,m,n,1) akpqparc2(l,m,n,1)=qparc2(l,m,n,1)*bmaginv(l) akpqpars2(l,m,n,1)=qpars2(l,m,n,1)*bmaginv(l) enddo kptparc2(ld,m,n,1)=0.0 kptpars2(ld,m,n,1)=0.0 akpqparc2(ld,m,n,1)=0.0 akpqpars2(ld,m,n,1)=0.0 enddo enddo call kparfft(kptparc2,kptpars2,0) call kparfft(akpqparc2,akpqpars2,1) c do n=1,nd do m=1,md do l=1,ld c c --------------parallel heat flux equation--------------- c qpars1(l,m,n,1)=qpars(l,m,n,1)+dt*(wpqps(l,m,n,1) cpm * -lambda1*pmpqps(l,m,n,1)+nu1*wgqps(l,m,n,1) * -zBq1*kptpars2(l,m,n,1)-zDq1*akpqpars2(l,m,n,1)*bmag(l) * +kdopp*qparc2(l,m,n,1) * +0.5*omegad(l,m,n)*((nu5i+6.)*vc2(l,m,n,1) * +(nu6i-3.)*qparc2(l,m,n,1)+(nu7i-3.)*qc2(l,m,n,1) * +nu11i*tparc2(l,m,n,1)+nu12i*tc2(l,m,n,1) ) * -0.5*abs(omegad(l,m,n))*(nu5r*vs2(l,m,n,1) * +nu6r*qpars2(l,m,n,1)+nu7r*qs2(l,m,n,1) * +nu11r*tpars2(l,m,n,1)+nu12r*ts2(l,m,n,1) ) c * -zBq1*tpars2(l,m,n,1)*bgrad(l) * -nuii*qpars2(l,m,n,1) ) c qparc1(l,m,n,1)=qparc(l,m,n,1)+dt*(wpqpc(l,m,n,1) cpm * -lambda1*pmpqpc(l,m,n,1)+nu1*wgqpc(l,m,n,1) * -zBq1*kptparc2(l,m,n,1)-zDq1*akpqparc2(l,m,n,1)*bmag(l) * -kdopp*qpars2(l,m,n,1) * -0.5*omegad(l,m,n)*((nu5i+6.)*vs2(l,m,n,1) * +(nu6i-3.)*qpars2(l,m,n,1)+(nu7i-3.)*qs2(l,m,n,1) * +nu11i*tpars2(l,m,n,1)+nu12i*ts2(l,m,n,1) ) * -0.5*abs(omegad(l,m,n))*(nu5r*vc2(l,m,n,1) * +nu6r*qparc2(l,m,n,1)+nu7r*qc2(l,m,n,1) * +nu11r*tparc2(l,m,n,1)+nu12r*tc2(l,m,n,1) ) c * -zBq1*tparc2(l,m,n,1)*bgrad(l) * -nuii*qparc2(l,m,n,1) ) c enddo enddo enddo c call bcon(qparc1,qpars1,1,0) endif c if(nperpmom.eq.4) goto 10 do n=1,nd do m=1,md do l=1,ld if(nperpmom.eq.1) then c c ----------one pole closure: perpendicular moments------------ c c c qs2(l,m,n,1)=rkkperp*(ts2(l,m,n,1)+uqs(l,m,n,1)) c qc2(l,m,n,1)=rkkperp*(tc2(l,m,n,1)+uqc(l,m,n,1)) qs2(l,m,n,1)=rkkperp*(ts2(l,m,n,1)) qc2(l,m,n,1)=rkkperp*(tc2(l,m,n,1)) else if(nperpmom.eq.3) then c c ---------three pole closure: perpendicular moments----------- c c this is not used, only good for 1 or 2 perp. moments, mbeer c c ss2(l,m,n,1)= sign(1.0,kpar)*rkkin*(rc2(l,m,n,1)-tc2(l,m,n,1)) c sc2(l,m,n,1)=-sign(1.0,kpar)*rkkin*(rs2(l,m,n,1)-ts2(l,m,n,1)) endif enddo enddo enddo 10 continue c do n=1,nd do m=1,md c c load kpqc2,kpqs2 c do l=1,ld-1 kpqc2(l,m,n,1)=qc2(l,m,n,1)*bmaginv(l)**2 kpqs2(l,m,n,1)=qs2(l,m,n,1)*bmaginv(l)**2 enddo kpqc2(ld,m,n,1)=0.0 kpqs2(ld,m,n,1)=0.0 enddo enddo if (nperpmom.eq.1) then call kparfft(kpqc2,kpqs2,1) else call kparfft(kpqc2,kpqs2,0) endif c do n=1,nd do m=1,md do l=1,ld c c ----------perpendicular temperature equation-------------- c ts1(l,m,n,1)=ts(l,m,n,1)+dt*(wpts(l,m,n,1) * +(wdts(l,m,n,1)+wgns(l,m,n,1)) cpm * -(lambda2*pmpts(l,m,n,1)-nu2*wgts(l,m,n,1)) * +rky(m,n)*utc(l,m,n,1)-kpqs2(l,m,n,1)*bmag(l)**2 * +kdopp*tc2(l,m,n,1) * +omegad(l,m,n)*(tc2(l,m,n,1)+wc2(l,m,n,1) * +(0.5+nu3i)*tparc2(l,m,n,1) * +(1.+nu4i)*tc2(l,m,n,1) * -0.5*prc2(l,m,n,1) * +urc(l,m,n,1)) * -nu4r*abs(omegad(l,m,n))*ts2(l,m,n,1) * -nu3r*abs(omegad(l,m,n))*tpars2(l,m,n,1) * +vs2(l,m,n,1)*bgrad(l) * +nuii/3.*(tpars2(l,m,n,1)-ts2(l,m,n,1)) ) c tc1(l,m,n,1)=tc(l,m,n,1)+dt*(wptc(l,m,n,1) * +(wdtc(l,m,n,1)+wgnc(l,m,n,1)) cpm * -(lambda2*pmptc(l,m,n,1)-nu2*wgtc(l,m,n,1)) * -rky(m,n)*uts(l,m,n,1)-kpqc2(l,m,n,1)*bmag(l)**2 * -kdopp*ts2(l,m,n,1) * -omegad(l,m,n)*(ts2(l,m,n,1)+ws2(l,m,n,1) * +(0.5+nu3i)*tpars2(l,m,n,1) * +(1.+nu4i)*ts2(l,m,n,1) * -0.5*prs2(l,m,n,1) * +urs(l,m,n,1)) * -nu4r*abs(omegad(l,m,n))*tc2(l,m,n,1) * -nu3r*abs(omegad(l,m,n))*tparc2(l,m,n,1) * +vc2(l,m,n,1)*bgrad(l) * +nuii/3.*(tparc2(l,m,n,1)-tc2(l,m,n,1)) ) enddo enddo enddo c call bcon(tc1,ts1,1,0) c if(nperpmom.ge.2) then do n=1,nd do m=1,md c c load kptc2,kpts2,akpqc2,akpqs2 c do l=1,ld-1 akpqc2(l,m,n,1)=qc2(l,m,n,1)*rkkperp2*bmaginv(l)**2 akpqs2(l,m,n,1)=qs2(l,m,n,1)*rkkperp2*bmaginv(l)**2 kptc2(l,m,n,1)=tc2(l,m,n,1)+uqc(l,m,n,1) kpts2(l,m,n,1)=ts2(l,m,n,1)+uqs(l,m,n,1) enddo akpqc2(ld,m,n,1)=0.0 akpqs2(ld,m,n,1)=0.0 kptc2(ld,m,n,1)=0.0 kpts2(ld,m,n,1)=0.0 enddo enddo call kparfft(akpqc2,akpqs2,1) call kparfft(kptc2,kpts2,0) c do n=1,nd do m=1,md do l=1,ld c c -----------perpendicular heat flux equation---------------- c qs1(l,m,n,1)=qs(l,m,n,1)+dt*(wpqs(l,m,n,1) * +wgvs(l,m,n,1) * -(kpts2(l,m,n,1)+akpqs2(l,m,n,1)*bmag(l)**2)+kdopp*qc2(l,m,n,1) * +0.5*omegad(l,m,n)*((nu8i+1.)*vc2(l,m,n,1) * +(nu9i-1.)*qparc2(l,m,n,1)+(nu10i-1.)*qc2(l,m,n,1) * +nu13i*tparc2(l,m,n,1)+nu14i*tc2(l,m,n,1) ) * -0.5*abs(omegad(l,m,n))*(nu8r*vs2(l,m,n,1) * +nu9r*qpars2(l,m,n,1)+nu10r*qs2(l,m,n,1) * +nu13r*tpars2(l,m,n,1)+nu14r*ts2(l,m,n,1) ) * -(ts2(l,m,n,1)-tpars2(l,m,n,1))*bgrad(l) * -nuii*qs2(l,m,n,1) ) c qc1(l,m,n,1)=qc(l,m,n,1)+dt*(wpqc(l,m,n,1) * +wgvc(l,m,n,1) * -(kptc2(l,m,n,1)+akpqc2(l,m,n,1)*bmag(l)**2)-kdopp*qs2(l,m,n,1) * -0.5*omegad(l,m,n)*((nu8i+1.)*vs2(l,m,n,1) * +(nu9i-1.)*qpars2(l,m,n,1)+(nu10i-1.)*qs2(l,m,n,1) * +nu13i*tpars2(l,m,n,1)+nu14i*ts2(l,m,n,1) ) * -0.5*abs(omegad(l,m,n))*(nu8r*vc2(l,m,n,1) * +nu9r*qparc2(l,m,n,1)+nu10r*qc2(l,m,n,1) * +nu13r*tparc2(l,m,n,1)+nu14r*tc2(l,m,n,1) ) * -(tc2(l,m,n,1)-tparc2(l,m,n,1))*bgrad(l) * -nuii*qc2(l,m,n,1) ) enddo enddo enddo call bcon(qc1,qs1,1,0) endif c if(nperpmom.ge.3) then do n=1,nd do m=1,md do l=1,ld c c -------------------r_perp equation-------------------- c c only good for 1 or 2 perp. moments, mbeer c rs1(l,m,n,1)=rps(l,m,n,1)+dt*(rky(m,n)*urc(l,m,n,1) * +kpart*(3.*qc2(l,m,n,1)+sc2(l,m,n,1))+kdopp*rc2(l,m,n,1)) c rc1(l,m,n,1)=rpc(l,m,n,1)+dt*(-rky(m,n)*urs(l,m,n,1) * -kpart*(3.*qs2(l,m,n,1)+ss2(l,m,n,1))-kdopp*rs2(l,m,n,1)) enddo enddo enddo c call bcon(rc1,rs1,1,0) endif c if(nperpmom.ge.4) then do n=1,nd do m=1,md do l=1,ld c c -------------------------s_perp equation------------------------ c ss1(l,m,n,1)=sps(l,m,n,1)+dt*( * +zBq1*kpart*(rc2(l,m,n,1)-tc2(l,m,n,1)) * -zDq1*abs(kpart)*ss2(l,m,n,1)+kdopp*sc2(l,m,n,1)) c sc1(l,m,n,1)=spc(l,m,n,1)+dt*( * -zBq1*kpart*(rs2(l,m,n,1)-ts2(l,m,n,1)) * -zDq1*abs(kpart)*sc2(l,m,n,1)-kdopp*ss2(l,m,n,1)) enddo enddo enddo call bcon(sc1,ss1,1,0) endif c c For now, specialize second species to electrons: c if(nspecies.eq.1) goto 100 do n=1,nd do m=1,md do l=1,ld c c -------------electron density equation----------------- c ws1(l,m,n,2)=ws(l,m,n,2)+dt*(wpns(l,m,n,2) * +rky(m,n)*unc(l,m,n,2)+kpart*vc2(l,m,n,2)+kdopp*wc2(l,m,n,2)) c wc1(l,m,n,2)=wc(l,m,n,2)+dt*(wpnc(l,m,n,2) * -rky(m,n)*uns(l,m,n,2)-kpart*vs2(l,m,n,2)-kdopp*ws2(l,m,n,2)) c c ------electron parallel velocity equation-------------- c vs1(l,m,n,2)=vs(l,m,n,2)+dt*(wpvs(l,m,n,2) * +rmime*kpart*(prc2(l,m,n,2)-uvc(l,m,n,2)) * -rmime*rky(m,n)*vz*uvc(l,m,n,2)+kdopp*vc2(l,m,n,2)) c vc1(l,m,n,2)=vc(l,m,n,2)+dt*(wpvc(l,m,n,2) * -rmime*kpart*(prs2(l,m,n,2)-uvs(l,m,n,2)) * +rmime*rky(m,n)*vz*uvs(l,m,n,2)-kdopp*vs2(l,m,n,2)) enddo enddo enddo c call bcon(wc1,ws1,2,0) call bcon(vc1,vs1,2,0) c if(nparmom.eq.3) then do n=1,nd do m=1,md do l=1,ld c c ---------three pole closure: parallel moments---------- c qpars2(l,m,n,2)=sign(1.0,kpart)*rkkin*tparc2(l,m,n,2) * *sqrt(rmime) qparc2(l,m,n,2)=-sign(1.0,kpart)*rkkin*tpars2(l,m,n,2) * *sqrt(rmime) enddo enddo enddo endif c do n=1,nd do m=1,md do l=1,ld c c ------electron parallel temperature equation-------- c tpars1(l,m,n,2)=tpars(l,m,n,2)+dt*(wptps(l,m,n,2) * +2.*kpart*vc2(l,m,n,2)+rky(m,n)*utparc(l,m,n,2) * +kpart*qparc2(l,m,n,2)+kdopp*tparc2(l,m,n,2)) c tparc1(l,m,n,2)=tparc(l,m,n,2)+dt*(wptpc(l,m,n,2) * -2.*kpart*vs2(l,m,n,2)-rky(m,n)*utpars(l,m,n,2) * -kpart*qpars2(l,m,n,2)-kdopp*tpars2(l,m,n,2)) c enddo enddo enddo c call bcon(tparc1,tpars1,2,0) c if(nparmom.eq.4) then do n=1,nd do m=1,md do l=1,ld c c --------electron parallel heat flux equation------------- c qpars1(l,m,n,2)=qpars(l,m,n,2)+dt*(wpqps(l,m,n,2) * +zEq1*rmime*kpart*tparc2(l,m,n,2)-sqrt(rmime)* * zDq1*abs(kpart)*qpars2(l,m,n,2)+kdopp*qparc2(l,m,n,2)) c qparc1(l,m,n,2)=qparc(l,m,n,2)+dt*(wpqpc(l,m,n,2) * -zEq1*rmime*kpart*tpars2(l,m,n,2)-sqrt(rmime)* * zDq1*abs(kpart)*qparc2(l,m,n,2)-kdopp*qpars2(l,m,n,2)) c enddo enddo enddo c call bcon(qparc1,qpars1,2,0) endif 100 continue c c Make like ron's wedge c goto 33 do m=1,md do n=1,nd if ((abs(r0(m,n)).gt.pi+.001).or.(m.eq.1)) then do l=1,ld ws1(l,m,n,1)=1.e-15 wc1(l,m,n,1)=1.e-15 vs1(l,m,n,1)=1.e-15 vc1(l,m,n,1)=1.e-15 tpars1(l,m,n,1)=1.e-15 tparc1(l,m,n,1)=1.e-15 qpars1(l,m,n,1)=1.e-15 qparc1(l,m,n,1)=1.e-15 ts1(l,m,n,1)=1.e-15 tc1(l,m,n,1)=1.e-15 qs1(l,m,n,1)=1.e-15 qc1(l,m,n,1)=1.e-15 enddo endif enddo enddo 33 continue return end