module fields_init use fields implicit none complex, dimension(:,:), allocatable, private :: ka contains subroutine init_fields use itg_data use gryffin_grid, only: ld, md, rkperp2, rky, rkx use gryffin_layouts use mp, only: broadcast, proc0 use io, only: wread, bwread, wreadnc integer :: ntim, ntimr, ndum1, ndum2 integer :: lin_sav, ne_sav, nfreq_sav, nstp_sav, ntotal_sav, ifilter_sav real :: dt1, tim, rmu1_sav real :: rkymin, rkymax character filename*80 integer :: l, m, n if(nread == 1) then ! ! Read previous *.ncp results file: ! ntim=0 ntimr=0 ne_sav=ne ! save some values ntotal_sav=ntotal nstp_sav=nstp nfreq_sav=nfreq lin_sav=lin ifilter_sav=ifilter rmu1_sav=rmu1 ndum1 = 1 ! These arguments are used only by postc. ndum2 = 1 ! These arguments are used only by postc. if(binary) then filename=runname(1:lrunname)//'.bresp' call bwread(filename,ntim,ntimr,tim,dt1,ndum1,ndum2) elseif(ascii) then filename=runname(1:lrunname)//'.resp' call wread(filename,ntim,ntimr,tim,dt1,ndum1,ndum2) elseif(cdf) then filename=runname(1:lrunname)//'.ncp' call wreadnc(filename,ntim,ntimr,tim,dt1,ndum1,ndum2) endif ! ! restore B variation of kperp rho which was removed for post ! do n = n_low, n_high do m = m_low, m_high do l=1,ld rkperp2(l,m,n)=rkperp2(l,m,n)*bmaginv(l)**2 enddo enddo enddo ! Reset some values just read from the previous results file, ! for example, restart time indexing from the beginning: ne=ne_sav ntotal=ntotal_sav nstp=nstp_sav nfreq_sav=nfreq lin=lin_sav ifilter=ifilter_sav rmu1=rmu1_sav else call pertb ! ! If icrit=1,2, advance each k_theta with a different dt. ! if(icrit /= 0) then if(idx_local(m_lo, 1, 1)) rkymin=rky(1,1) call broadcast(rkymin, proc_id(m_lo, 1, 1)) if(rkymin == 0.) then rkymin=1. if(md > 1 .and. idx_local(m_lo, 2, 1)) rkymin = rky(2,1) call broadcast(rkymin, proc_id(m_lo, 2, 1)) endif if(idx_local(m_lo, md, 1)) then rkymax=rky(md,1) if(rkymax == 0.) rkymax=1. endif call broadcast(rkymax, proc_id(m_lo, md, 1)) do n = n_low, n_high do m = m_low, m_high time(m)=0.0 dt(m)=rkymax/max(rkymin,rky(m,n))*dt0 enddo enddo else dt = dt0 endif endif end subroutine init_fields subroutine pertb ! ! (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. ! ! ! This subroutine inserts the perturbation. ! use itg_data use gryffin_layouts use mp use gryffin_grid, only: nd, nddp, nddm, ld, md, ldb, r, rkx, malias, nalias, nadd use constants, only: ii use nlps, only: reality use gryffin_redist use redistribute, only: gather, scatter type (gryffin_redist_type) :: redist real a,b,ax,g05caE,c1,c2,c3,c4,c5,c6,p00 integer i,l,m,n if(proc0) 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 = n_low, n_high do m = m_low, m_high 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=m_lo%nl_world, m_lo%nu_world do m=m_lo%ml_world, m_lo%mu_world a=g05caE(a) b=g05caE(b) if(idx_local(m_lo, m, n)) then do l=2,ldb potential(l,m,n)=cmplx(a*real(potential(l,m,n)),b*aimag(potential(l,m,n))) enddo endif 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 do n=n_low, n_high if (idx_local(m_lo, 1, n) .and. beta_e > 0) & potential(:,1,n) = 0.1*potential(:,1,n) do m=m_low, m_high density(:,m,n,i)=potential(:,m,n) enddo enddo enddo ! initialize n_e=phi for finite beta case PS 6/98 if (beta_e > 0.) then do n=n_low, n_high do m=m_low, m_high do l=1,ldb n_e(l,m,n)=potential(l,m,n) enddo enddo enddo endif ! 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) then density = 0. potential = 0. if(pmag .gt. 0.0) then do i=1,nspecies do n = n_low, n_high do m = m_low, m_high 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: ! ! not implemented for distributed kx's if(nd.ne.1) then call init_gryffin_redist (ldb, md, nd, malias, nalias, & nddp, nadd, redist) allocate (ka(nalias, & redist%x_lo%llim_proc:redist%x_lo%ulim_alloc)) call gather(redist%kx2x, potential, ka) call reality(redist, ka) call scatter(redist%kx2x, ka, potential) do i = 1, nspecies call gather(redist%kx2x, density(:,:,:,i), ka) call reality(redist, ka) call scatter(redist%kx2x, ka, density(:,:,:,i)) call gather(redist%kx2x, u_par(:,:,:,i), ka) call reality(redist, ka) call scatter(redist%kx2x, ka, u_par(:,:,:,i)) call gather(redist%kx2x, t_par(:,:,:,i), ka) call reality(redist, ka) call scatter(redist%kx2x, ka, t_par(:,:,:,i)) call gather(redist%kx2x, t_perp(:,:,:,i), ka) call reality(redist, ka) call scatter(redist%kx2x, ka, t_perp(:,:,:,i)) enddo deallocate (ka) call delete_gryffin_redist(redist) ! ! This could be set to any real number to satisfy the ! reality condition. ! if(idx_local(m_lo, 1, 1)) then 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 endif end subroutine pertb end module fields_init