FIXME : Add documentation
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| logical, | intent(in) | :: | write_text |
Write to text file |
subroutine do_write_gs(write_text)
use file_utils, only: open_output_file, close_output_file
use mp, only: proc0, sum_reduce, iproc, nproc
use gs2_transforms, only: transform2, inverse2
use fields_arrays, only: apar, phi
use kt_grids, only: naky, ntheta0, akx, aky, nx, ny
use theta_grid, only: ntgrid, delthet, jacob, gradpar, nperiod
use constants, only: pi, zi
use splines, only: fitp_surf1, fitp_surf2
implicit none
!> Write to text file
logical, intent(in) :: write_text
real, dimension (:), allocatable :: total
real, dimension (:,:,:), allocatable :: bxf, byf, vxf, vyf, bxfsavg, byfsavg
real, dimension (:,:,:), allocatable :: bxfs, byfs, vxfs, vyfs, rvx, rvy, rx, ry
real, dimension (:), allocatable :: xx4, yy4, dz
real, dimension (:,:), allocatable :: bxs, bys, vxs, vys
real, dimension (:,:), allocatable :: bxsavg, bysavg
complex, dimension (:,:,:), allocatable :: bx, by, vx, vy, vx2, vy2
real, dimension (:), allocatable :: stemp, zx1, zxm, zy1, zyn, xx, yy
real :: zxy11, zxym1, zxy1n, zxymn, rxt, ryt, bxt, byt, L_x, L_y
real, dimension (2*ntgrid) :: kpar
integer :: nnx, nny, nnx4, nny4, ulim, llim, iblock, i, unit
integer :: ik, it, ig, ierr
if (.not. write_text) return
nny = 2*ny
nnx = 2*nx
nnx4 = nnx+4
nny4 = nny+4
allocate (bx(-ntgrid:ntgrid,ntheta0,naky))
allocate (by(-ntgrid:ntgrid,ntheta0,naky))
allocate (vx(-ntgrid:ntgrid,ntheta0,naky))
allocate (vy(-ntgrid:ntgrid,ntheta0,naky))
do ik=1,naky
do it=1,ntheta0
do ig=-ntgrid, ntgrid
bx(ig,it,ik) = zi*aky(ik)*apar(ig,it,ik)
by(ig,it,ik) = -zi*akx(it)*apar(ig,it,ik)
vx(ig,it,ik) = -zi*aky(ik)*phi(ig,it,ik)
vy(ig,it,ik) = zi*akx(it)*phi(ig,it,ik)
end do
end do
end do
allocate (bxf(nnx,nny,-ntgrid:ntgrid))
allocate (byf(nnx,nny,-ntgrid:ntgrid))
allocate (vxf(nnx,nny,-ntgrid:ntgrid))
allocate (vyf(nnx,nny,-ntgrid:ntgrid))
call transform2 (bx, bxf, nny, nnx)
call transform2 (by, byf, nny, nnx)
call transform2 (vx, vxf, nny, nnx)
call transform2 (vy, vyf, nny, nnx)
! fields come out as (x, y, z)
deallocate (bx, by)
allocate ( bxfs(nnx4, nny4, -ntgrid:ntgrid))
allocate ( byfs(nnx4, nny4, -ntgrid:ntgrid))
allocate (bxfsavg(nnx4, nny4, -ntgrid:ntgrid))
allocate (byfsavg(nnx4, nny4, -ntgrid:ntgrid))
allocate ( vxfs(nnx4, nny4, -ntgrid:ntgrid))
allocate ( vyfs(nnx4, nny4, -ntgrid:ntgrid))
do ig=-ntgrid,ntgrid
do ik=1,2
do it=3,nnx4-2
bxfs(it,ik,ig) = bxf(it-2,nny-2+ik,ig)
byfs(it,ik,ig) = byf(it-2,nny-2+ik,ig)
vxfs(it,ik,ig) = vxf(it-2,nny-2+ik,ig)
vyfs(it,ik,ig) = vyf(it-2,nny-2+ik,ig)
bxfs(it,nny4-2+ik,ig) = bxf(it-2,ik,ig)
byfs(it,nny4-2+ik,ig) = byf(it-2,ik,ig)
vxfs(it,nny4-2+ik,ig) = vxf(it-2,ik,ig)
vyfs(it,nny4-2+ik,ig) = vyf(it-2,ik,ig)
end do
end do
do ik=3,nny4-2
do it=3,nnx4-2
bxfs(it,ik,ig) = bxf(it-2,ik-2,ig)
byfs(it,ik,ig) = byf(it-2,ik-2,ig)
vxfs(it,ik,ig) = vxf(it-2,ik-2,ig)
vyfs(it,ik,ig) = vyf(it-2,ik-2,ig)
end do
end do
do ik=1,nny4
do it=1,2
bxfs(it,ik,ig) = bxfs(nnx4-4+it,ik,ig)
byfs(it,ik,ig) = byfs(nnx4-4+it,ik,ig)
vxfs(it,ik,ig) = vxfs(nnx4-4+it,ik,ig)
vyfs(it,ik,ig) = vyfs(nnx4-4+it,ik,ig)
bxfs(nnx4-2+it,ik,ig) = bxfs(it+2,ik,ig)
byfs(nnx4-2+it,ik,ig) = byfs(it+2,ik,ig)
vxfs(nnx4-2+it,ik,ig) = vxfs(it+2,ik,ig)
vyfs(nnx4-2+it,ik,ig) = vyfs(it+2,ik,ig)
end do
end do
end do
deallocate (vxf, vyf)
allocate (xx4(nnx4), xx(nnx))
allocate (yy4(nny4), yy(nny))
L_x = 2.0*pi/akx(2)
L_y = 2.0*pi/aky(2)
do it = 1, nnx
xx4(it+2) = real(it-1)*L_x/real(nnx)
xx(it) = real(it-1)*L_x/real(nnx)
end do
do it=1,2
xx4(it) = xx4(nnx4-4+it)-L_x
xx4(nnx4-2+it) = xx4(it+2)+L_x
end do
do ik = 1, nny
yy4(ik+2) = real(ik-1)*L_y/real(nny)
yy(ik) = real(ik-1)*L_y/real(nny)
end do
do ik=1,2
yy4(ik) = yy4(nny4-4+ik)-L_y
yy4(nny4-2+ik) = yy4(ik+2)+L_y
end do
allocate (dz(-ntgrid:ntgrid))
dz = delthet*jacob
allocate ( bxs(3*nnx4*nny4,-ntgrid:ntgrid)) ; bxs = 0.
allocate ( bys(3*nnx4*nny4,-ntgrid:ntgrid)) ; bys = 0.
allocate ( vxs(3*nnx4*nny4,-ntgrid:ntgrid)) ; vxs = 0.
allocate ( vys(3*nnx4*nny4,-ntgrid:ntgrid)) ; vys = 0.
allocate (bxsavg(3*nnx4*nny4,-ntgrid:ntgrid))
allocate (bysavg(3*nnx4*nny4,-ntgrid:ntgrid))
allocate (stemp(nnx4+2*nny4))
allocate (zx1(nny4), zxm(nny4), zy1(nnx4), zyn(nnx4))
do ig=-ntgrid, ntgrid
call fitp_surf1(nnx4, nny4, xx4, yy4, bxfs(:,:,ig), &
nnx4, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, &
255, bxs(:,ig), stemp, 1., ierr)
call fitp_surf1(nnx4, nny4, xx4, yy4, byfs(:,:,ig), &
nnx4, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, &
255, bys(:,ig), stemp, 1., ierr)
call fitp_surf1(nnx4, nny4, xx4, yy4, vxfs(:,:,ig), &
nnx4, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, &
255, vxs(:,ig), stemp, 1., ierr)
call fitp_surf1(nnx4, nny4, xx4, yy4, vyfs(:,:,ig), &
nnx4, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, &
255, vys(:,ig), stemp, 1., ierr)
end do
deallocate (zx1, zxm, zy1, zyn, stemp)
do ig=-ntgrid, ntgrid-1
bxsavg(:,ig) = 0.5*(bxs(:,ig)+bxs(:,ig+1))
bysavg(:,ig) = 0.5*(bys(:,ig)+bys(:,ig+1))
bxfsavg(:,:,ig) = 0.5*(bxfs(:,:,ig)+bxfs(:,:,ig+1))
byfsavg(:,:,ig) = 0.5*(byfs(:,:,ig)+byfs(:,:,ig+1))
end do
! now, integrate to find a field line
allocate ( rx(nnx,nny,-ntgrid:ntgrid))
allocate ( ry(nnx,nny,-ntgrid:ntgrid))
allocate (rvx(nnx,nny,-ntgrid:ntgrid)) ; rvx = 0.
allocate (rvy(nnx,nny,-ntgrid:ntgrid)) ; rvy = 0.
do ik=1,nny
do it=1,nnx
rx(it,ik,-ntgrid) = xx(it)
ry(it,ik,-ntgrid) = yy(ik)
end do
end do
!Should these not come from layouts object?
iblock = (nnx*nny-1)/nproc + 1
llim = 1 + iblock * iproc
ulim = min(nnx*nny, llim+iblock-1)
do i=llim, ulim
it = 1 + mod(i-1, nnx)
ik = 1 + mod((i-1)/nnx, nny)
ig = -ntgrid
rxt = rx(it,ik,ig)
ryt = ry(it,ik,ig)
rvx(it,ik,ig) = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, vxfs(:,:,ig), nnx4, vxs(:,ig), 1.)
rvy(it,ik,ig) = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, vyfs(:,:,ig), nnx4, vys(:,ig), 1.)
do ig=-ntgrid,ntgrid-1
bxt = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, bxfs(:,:,ig), nnx4, bxs(:,ig), 1.)
byt = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, byfs(:,:,ig), nnx4, bys(:,ig), 1.)
rxt = rx(it,ik,ig) + 0.5*dz(ig)*bxt
ryt = ry(it,ik,ig) + 0.5*dz(ig)*byt
if (rxt > L_x) rxt = rxt - L_x
if (ryt > L_y) ryt = ryt - L_y
if (rxt < 0.) rxt = rxt + L_x
if (ryt < 0.) ryt = ryt + L_y
bxt = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, bxfsavg(:,:,ig), nnx4, bxsavg(:,ig), 1.)
byt = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, byfsavg(:,:,ig), nnx4, bysavg(:,ig), 1.)
rxt = rx(it,ik,ig) + dz(ig)*bxt
ryt = ry(it,ik,ig) + dz(ig)*byt
if (rxt > L_x) rxt = rxt - L_x
if (ryt > L_y) ryt = ryt - L_y
if (rxt < 0.) rxt = rxt + L_x
if (ryt < 0.) ryt = ryt + L_y
rx(it,ik,ig+1) = rxt
ry(it,ik,ig+1) = ryt
rvx(it,ik,ig+1) = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, vxfs(:,:,ig+1), nnx4, vxs(:,ig+1), 1.)
rvy(it,ik,ig+1) = fitp_surf2(rxt, ryt, nnx4, nny4, xx4, yy4, vyfs(:,:,ig+1), nnx4, vys(:,ig+1), 1.)
end do
end do
deallocate (bxfs, byfs, bxfsavg, byfsavg, vxfs, vyfs)
deallocate (rx, ry, bxs, bys, vxs, vys, bxsavg, bysavg)
allocate (total(2*nnx*nny*(2*ntgrid+1)))
i=1
do ig=-ntgrid,ntgrid
do ik=1,nny
do it=1,nnx
total(i) = rvx(it,ik,ig)
total(i+1) = rvy(it,ik,ig)
i = i + 2
end do
end do
end do
call sum_reduce(total, 0)
i=1
do ig=-ntgrid,ntgrid
do ik=1,nny
do it=1,nnx
rvx(it,ik,ig) = total(i)
rvy(it,ik,ig) = total(i+1)
i = i + 2
end do
end do
end do
call inverse2 (rvx, vx, nny, nnx)
call inverse2 (rvy, vy, nny, nnx)
allocate (vx2(-ntgrid:ntgrid,ntheta0,naky))
allocate (vy2(-ntgrid:ntgrid,ntheta0,naky))
call par_spectrum (vx, vx2)
call par_spectrum (vy, vy2)
if (.not. proc0) return
call open_output_file (unit, ".gs")
do ig = 1, ntgrid
kpar(ig) = (ig-1)*gradpar(ig)/real(2*nperiod-1)
kpar(2*ntgrid-ig+1)=-(ig)*gradpar(ig)/real(2*nperiod-1)
end do
do ik = 1, naky
do it = 1, ntheta0
do ig = ntgrid+1,2*ntgrid
write (unit, "(9(1x,e12.5))") &
kpar(ig), aky(ik), akx(it), &
real(vx2(ig-ntgrid-1,it,ik)), &
real(vy2(ig-ntgrid-1,it,ik))
end do
do ig = 1, ntgrid
write (unit, "(9(1x,e12.5))") &
kpar(ig), aky(ik), akx(it), &
real(vx2(ig-ntgrid-1,it,ik)), &
real(vy2(ig-ntgrid-1,it,ik))
end do
write (unit, "()")
end do
end do
call close_output_file (unit)
deallocate (vx2, vy2)
deallocate (bxf, byf, xx4, xx, yy4, yy, dz, total)
deallocate (vx, vy, rvx, rvy)
end subroutine do_write_gs