FIXME : Add documentation
subroutine write_output
use splines, only: inter_cspl
use gs2_io_grids, only: nc_write_grid_sizes, nc_write_grid_scalars, nc_write_grids
use gs2_io_grids, only: nc_write_lambda_grid
use gs2_io_grids, only: nc_grid_file_open, nc_grid_file_close
use constants, only: run_name_size
use mp, only: proc0
implicit none
integer :: unit
integer :: i
! real, allocatable, dimension (:) :: bmaginaux, bmagsmaux, tmp
real, allocatable, dimension (:) :: thetagridout
real, allocatable, dimension (:) :: gbdriftout, gradparout, grhoout
real, allocatable, dimension (:) :: cvdriftout, gds2out, bmaggridout
real, allocatable, dimension (:) :: gds21out, gds22out
real, allocatable, dimension (:) :: cvdrift0out, gbdrift0out
real, allocatable, dimension (:) :: Rplotout, Rprimeout
real, allocatable, dimension (:) :: Zplotout, Zprimeout
real, allocatable, dimension (:) :: aplotout, aprimeout
! real :: th, bmin
real :: bmin
! integer :: ierr
integer, dimension (1) :: minloca
character(len=run_name_size) :: ncout
call open_output_file (unit, ".input.out")
write (unit=unit, fmt="('#',i5)") nthetain+1
do i = 1, nthetain+1
write (unit=unit, fmt="(3(1x,g19.12))") thetain(i), bmagin(i), bmagsm(i)
enddo
call close_output_file (unit)
! allocate (bmaginaux(nthetain+1), bmagsmaux(nthetain+1), tmp(2*nthetain))
! call fitp_curvp1 (nthetain,thetain,bmagin,twopi,bmaginaux,tmp,1.0,ierr)
! call fitp_curvp1 (nthetain,thetain,bmagsm,twopi,bmagsmaux,tmp,1.0,ierr)
! call open_output_file (unit, ".input.fine")
! do i = -nfinegrid, nfinegrid
! th = pi*real(i)/real(nfinegrid)
! write (unit=unit, fmt="(3(x,g19.12))") th, &
! fitp_curvp2(th,nthetain,thetain,bmagin,twopi,bmaginaux,1.0), &
! fitp_curvp2(th,nthetain,thetain,bmagsm,twopi,bmagsmaux,1.0)
! end do
! call close_output_file (unit)
! deallocate (bmaginaux, bmagsmaux, tmp)
call open_output_file (unit, ".bmag.out")
write (unit=unit, fmt="('#',i5)") ntheta
do i = 1, ntheta
write (unit=unit, fmt="(2(1x,g19.12))") thetaout(i), bmagout(i)
enddo
call close_output_file (unit)
call open_output_file (unit, ".lambda.out")
write (unit=unit, fmt="('#',i5)") nlambda
do i = 1, nlambda
write (unit=unit, fmt=*) alambdaout(i)
enddo
call close_output_file (unit)
ntgrid = ntheta/2 + (nperiodout-1)*ntheta
allocate (thetagridout(-ntgrid:ntgrid))
allocate (gbdriftout(-ntgrid:ntgrid))
allocate (gradparout(-ntgrid:ntgrid))
allocate (grhoout(-ntgrid:ntgrid))
allocate (cvdriftout(-ntgrid:ntgrid))
allocate (gds2out(-ntgrid:ntgrid))
allocate (bmaggridout(-ntgrid:ntgrid))
allocate (gds21out(-ntgrid:ntgrid))
allocate (gds22out(-ntgrid:ntgrid))
allocate (cvdrift0out(-ntgrid:ntgrid))
allocate (gbdrift0out(-ntgrid:ntgrid))
allocate (Rplotout(-ntgrid:ntgrid))
allocate (Rprimeout(-ntgrid:ntgrid))
allocate (Zplotout(-ntgrid:ntgrid))
allocate (Zprimeout(-ntgrid:ntgrid))
allocate (aplotout(-ntgrid:ntgrid))
allocate (aprimeout(-ntgrid:ntgrid))
thetagridout(-ntheta/2:ntheta/2-1) = thetaout(:ntheta)
bmaggridout(-ntheta/2:ntheta/2-1) = bmagout(:ntheta)
thetagridout(ntheta/2) = thetagridout(-ntheta/2) + twopi*iperiod
bmaggridout(ntheta/2) = bmaggridout(-ntheta/2)
do i = 1, nperiodout - 1
thetagridout(-ntheta/2-ntheta*i:ntheta/2-ntheta*i-1) &
= thetagridout(-ntheta/2:ntheta/2-1) - twopi*i
thetagridout(-ntheta/2+ntheta*i+1:ntheta/2+ntheta*i) &
= thetagridout(-ntheta/2+1:ntheta/2) + twopi*i
bmaggridout(-ntheta/2-ntheta*i:ntheta/2-ntheta*i-1) &
= bmaggridout(-ntheta/2:ntheta/2-1)
bmaggridout(-ntheta/2+ntheta*i+1:ntheta/2+ntheta*i) &
= bmaggridout(-ntheta/2+1:ntheta/2)
end do
call inter_cspl(thetagrid, gbdrift, thetagridout, gbdriftout)
call inter_cspl(thetagrid, gradpar, thetagridout, gradparout)
call inter_cspl(thetagrid, grho, thetagridout, grhoout)
call inter_cspl(thetagrid, cvdrift, thetagridout, cvdriftout)
call inter_cspl(thetagrid, gds2, thetagridout, gds2out)
call inter_cspl(thetagrid, bmag, thetagridout, bmagout)
call inter_cspl(thetagrid, gds21, thetagridout, gds21out)
call inter_cspl(thetagrid, gds22, thetagridout, gds22out)
call inter_cspl(thetagrid, cvdrift0, thetagridout, cvdrift0out)
call inter_cspl(thetagrid, gbdrift0, thetagridout, gbdrift0out)
call inter_cspl(thetagrid, Rplot, thetagridout, Rplotout)
call inter_cspl(thetagrid, Rprime, thetagridout, Rprimeout)
call inter_cspl(thetagrid, Zplot, thetagridout, Zplotout)
call inter_cspl(thetagrid, Zprime, thetagridout, Zprimeout)
call inter_cspl(thetagrid, aplot, thetagridout, aplotout)
call inter_cspl(thetagrid, aprime, thetagridout, aprimeout)
call open_output_file (unit, ".out")
write (unit=unit, fmt=*) 'nlambda'
write (unit=unit, fmt=*) nlambda
write (unit=unit, fmt=*) 'lambda'
do i = 1, nlambda
write (unit=unit, fmt=*) alambdaout(i)
enddo
write (unit=unit, fmt=*) 'ntgrid nperiod ntheta drhodpsi rmaj shat kxfac q'
write (unit=unit, fmt="(3i4,5(1x,g19.10))") ntgrid, nperiodout, ntheta, drhodpsi, rmaj, shat, kxfac, qval
write (unit=unit, fmt=*) 'gbdrift gradpar grho tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
gbdriftout(i),gradparout(i),grhoout(i),thetagridout(i)
end do
write (unit=unit, fmt=*) 'cvdrift gds2 bmag tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
cvdriftout(i),gds2out(i),bmaggridout(i),thetagridout(i)
end do
write (unit=unit, fmt=*) 'gds21 gds22 tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
gds21out(i),gds22out(i),thetagridout(i)
end do
write (unit=unit, fmt=*) 'cvdrift0 gbdrift0 tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
cvdrift0out(i),gbdrift0out(i),thetagridout(i)
end do
write (unit=unit, fmt=*) 'Rplot Rprime tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
Rplotout(i),Rprimeout(i),thetagridout(i)
end do
write (unit=unit, fmt=*) 'Zplot Zprime tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
Zplotout(i),Zprimeout(i),thetagridout(i)
end do
write (unit=unit, fmt=*) 'aplot aprime tgrid'
do i = -ntgrid,ntgrid
write (unit=unit, fmt="(8(1x,g19.10))") &
aplotout(i),aprimeout(i),thetagridout(i)
end do
call close_output_file (unit)
if (screenout) then
write (*, *) 'cvdrift gds2 bmag theta'
do i = -ntheta/2,ntheta/2
write (unit=*, fmt="(4(1x,g13.5))") &
cvdriftout(i),gds2out(i),bmaggridout(i),thetagridout(i)
end do
end if
minloca = minloc(bmagout(:ntheta))
bmin = bmagout(minloca(1))
print *, "theta,bmag minimum: ", thetaout(minloca(1)), bmin
minloca = minloc(bmagin)
print *, "theta,bmag input minimum: ", &
thetain(minloca(1)), bmagin(minloca(1))
if (bmin < bmagin(minloca(1))) print *, "WARNING: interpolated new minimum"
if (proc0) then
ncout = trim(run_name)//".out.nc"
call nc_grid_file_open(ncout, "w")
call nc_write_grid_sizes(ntheta, ntgrid, nperiodout)
call nc_write_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)
call nc_write_grids(ntgrid, bmaggridout, gradparout, gbdriftout, gbdrift0out, &
cvdriftout, cvdrift0out, gds2out, gds21out, gds22out, grhoout, thetagridout, &
Rplot=Rplotout, Rprime=Rprimeout, Zplot=Zplotout, Zprime=Zprimeout, &
aplot=aplotout, aprime=aprimeout, no_end_point_in=no_end_point)
call nc_write_lambda_grid(nlambda, alambdaout(:nlambda))
call nc_grid_file_close()
end if
end subroutine write_output