subroutine pertb ! ! This subroutine inserts the perturbation. ! use itg_data implicit none real a,b,ax,g05caE,c1,c2,c3,c4,c5,c6,p00 integer i,l,m,n write(9,302) ncycle 302 format(i7,'***** perturbation inserted ***************') potential = 0. density = 0. u_par = 0. t_par = 0. t_perp = 0. phiav_old = 0. ! ! Choose parity of run. ! a=1. b=1. if(iodd.eq.1) a=0. !iodd = 1 -> odd modes only if(iodd.eq.2) b=0. !iodd = 2 -> even modes only do n=1,nd do m=1,md do l=2,ldb p00=pmag*(b*r(l)+a)*exp(-r(l)**2) potential(l,m,n)=cmplx(p00,-p00) enddo enddo enddo ! ! randomize: don't want to multitask this section ! !CMIC$ GUARD call g05cbE(iseed) !random ic's, initialize seed do n=1,nd do m=1,md a=g05caE(a) b=g05caE(b) do l=2,ldb potential(l,m,n)=cmplx(a*real(potential(l,m,n)),b*aimag(potential(l,m,n))) enddo enddo enddo !CMIC$ END GUARD ! Does f90 understand this directive in the .f format? ! ! gaussian blob in center of box ! ! a=pi*y0/4.*pmag ! do l=1,ldb ! do n=1,nd ! us(l,1,n,1)=0.0 ! uc(l,1,n,1)=a*exp(-(a**2/2.*rkx2(ldb/2,1,n))) ! do m=2,md ! us(l,m,n,1)=0.0 ! uc(l,m,n,1)=2.*a*exp(-(a**2/2.*(rkx2(ldb/2,m,n) ! & +rky2(m,n)))) ! enddo ! enddo ! enddo ! ! the following code for initializing the density needs to be upgraded ! to use the new FLR approximations: ! ! ! bill suggests not doing this do i=1,nspecies density(:,:,:,i)=potential(:,:,:) enddo ! iphi00=9 = set phi=0 for all Fourier components, not just (m=0,n=0). ! (no electric field, for neoclassical comparisons): ! Also, for this neoclassical test case, initialize a cos(k_r*r) ! perturbation, which is constant along the field line, either in ! the density (pmag > 0) or the temperatures (pmag < 0). ! if((iphi00 .eq. 9).or.((md.eq.1).and.(mrr(1,1).eq.0))) then density = 0. potential = 0. if(pmag .gt. 0.0) then do i=1,nspecies do n=1,nd do m=1,md do l=1,ld ! bumpy equilibria c1=1.0 c2=0. c3=8*c2-2*c1 c4=0. c5=0. c6=0. ! density(l,m,n,i)=rkkperp2*eps*( (c1-2*c2)*sin(r(l))-eps/2*(c1-c2)* ! & sin(2*r(l)) )+eps*(c5-c4)*cos(r(l))+c6 ! u_par(l,m,n,i)=c1*b_mag(l) ! q_perp(l,m,n,i)=-c1*b_mag(l)+c2*b_mag(l)**2 ! q_par(l,m,n,i)=c3*b_mag(l)-2*c2*b_mag(l)**2 ! t_par(l,m,n,i)=-rkkperp2*eps*( (c1-2*c2)*sin(r(l))-eps/2*(c1-c2)* ! & sin(2*r(l))) + c4 ! t_perp(l,m,n,i)=-rkkperp2*eps*( (c1-2*c2)*sin(r(l))-eps/2*(c1-3*c2)* ! & sin(2*r(l)))-c4*eps*(cos(r(l))+eps)+c5*bmaginv(l) ! potential(l,m,n)=0.0 ! ! good flow initialization: density(l,m,n,i)=pmag*rkx(l,m,n)*(1.-ii) potential(l,m,n)=pmag/rkx(l,m,n)*(1.-ii) ! u_par(l,m,n,i)=(pmag*(1-ii) ! & -2*ii*rkx(l,m,n)*qsf*cos(r(l)) ! & *potential(l,m,n))*b_mag(l) enddo enddo enddo enddo else t_perp = abs(pmag) t_par = abs(pmag) endif endif ! ! Enforce the reality condition for the ky=0 components: ! if(nd.ne.1) then do i=1,nspecies do n=1,nddm do l=1,ldb potential(l,1,nddp+n)=conjg(potential(l,1,nddp+1-n)) density(l,1,nddp+n,i)=conjg(density(l,1,nddp+1-n,i)) u_par(l,1,nddp+n,i)=conjg(u_par(l,1,nddp+1-n,i)) t_par(l,1,nddp+n,i)=conjg(t_par(l,1,nddp+1-n,i)) t_perp(l,1,nddp+n,i)=conjg(t_perp(l,1,nddp+1-n,i)) enddo enddo enddo ! ! This could be set to any real number to satisfy the ! reality condition. ! do i=1,nspecies do l=2,ldb potential(l,1,1)=0. density(l,1,1,i)=0. u_par(l,1,1,i)=0. t_par(l,1,1,i)=0. t_perp(l,1,1,i)=0. enddo enddo endif end subroutine pertb