gvec_to_jorek_prepare Subroutine

public subroutine gvec_to_jorek_prepare(X1_base_in, X1_in, X2_base_in, X2_in, LG_base_in, LG_in)

Uses

  • proc~~gvec_to_jorek_prepare~~UsesGraph proc~gvec_to_jorek_prepare gvec_to_jorek_prepare module~modgvec_base MODgvec_base proc~gvec_to_jorek_prepare->module~modgvec_base module~modgvec_fbase MODgvec_fBase proc~gvec_to_jorek_prepare->module~modgvec_fbase module~modgvec_globals MODgvec_Globals proc~gvec_to_jorek_prepare->module~modgvec_globals module~modgvec_gvec_to_jorek_vars MODgvec_gvec_to_jorek_Vars proc~gvec_to_jorek_prepare->module~modgvec_gvec_to_jorek_vars module~modgvec_readstate_vars MODgvec_ReadState_Vars proc~gvec_to_jorek_prepare->module~modgvec_readstate_vars module~modgvec_base->module~modgvec_fbase module~modgvec_base->module~modgvec_globals module~modgvec_sbase MODgvec_sBase module~modgvec_base->module~modgvec_sbase module~modgvec_sgrid MODgvec_sGrid module~modgvec_base->module~modgvec_sgrid module~modgvec_fbase->module~modgvec_globals iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env module~modgvec_gvec_to_jorek_vars->module~modgvec_base module~modgvec_gvec_to_jorek_vars->module~modgvec_fbase module~modgvec_gvec_to_jorek_vars->module~modgvec_globals module~modgvec_readstate_vars->module~modgvec_base module~modgvec_readstate_vars->module~modgvec_globals module~modgvec_hmap MODgvec_hmap module~modgvec_readstate_vars->module~modgvec_hmap module~modgvec_readstate_vars->module~modgvec_sbase module~modgvec_readstate_vars->module~modgvec_sgrid module~modgvec_c_hmap MODgvec_c_hmap module~modgvec_hmap->module~modgvec_c_hmap module~modgvec_hmap_axisnb MODgvec_hmap_axisNB module~modgvec_hmap->module~modgvec_hmap_axisnb module~modgvec_hmap_cyl MODgvec_hmap_cyl module~modgvec_hmap->module~modgvec_hmap_cyl module~modgvec_hmap_frenet MODgvec_hmap_frenet module~modgvec_hmap->module~modgvec_hmap_frenet module~modgvec_hmap_knot MODgvec_hmap_knot module~modgvec_hmap->module~modgvec_hmap_knot module~modgvec_hmap_rz MODgvec_hmap_RZ module~modgvec_hmap->module~modgvec_hmap_rz module~modgvec_sbase->module~modgvec_globals module~modgvec_sbase->module~modgvec_sgrid module~sll_m_bsplines sll_m_bsplines module~modgvec_sbase->module~sll_m_bsplines module~sll_m_spline_interpolator_1d sll_m_spline_interpolator_1d module~modgvec_sbase->module~sll_m_spline_interpolator_1d module~sll_m_spline_matrix sll_m_spline_matrix module~modgvec_sbase->module~sll_m_spline_matrix module~modgvec_sgrid->module~modgvec_globals module~modgvec_c_hmap->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_fbase module~modgvec_hmap_axisnb->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_c_hmap module~modgvec_io_netcdf MODgvec_IO_NETCDF module~modgvec_hmap_axisnb->module~modgvec_io_netcdf module~modgvec_hmap_cyl->module~modgvec_globals module~modgvec_hmap_cyl->module~modgvec_c_hmap module~modgvec_hmap_frenet->module~modgvec_globals module~modgvec_hmap_frenet->module~modgvec_c_hmap module~modgvec_hmap_knot->module~modgvec_globals module~modgvec_hmap_knot->module~modgvec_c_hmap module~modgvec_hmap_rz->module~modgvec_globals module~modgvec_hmap_rz->module~modgvec_c_hmap module~sll_m_assert sll_m_assert module~sll_m_bsplines->module~sll_m_assert module~sll_m_bsplines_base sll_m_bsplines_base module~sll_m_bsplines->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform sll_m_bsplines_non_uniform module~sll_m_bsplines->module~sll_m_bsplines_non_uniform module~sll_m_bsplines_uniform sll_m_bsplines_uniform module~sll_m_bsplines->module~sll_m_bsplines_uniform module~sll_m_errors sll_m_errors module~sll_m_bsplines->module~sll_m_errors module~sll_m_working_precision sll_m_working_precision module~sll_m_bsplines->module~sll_m_working_precision module~sll_m_spline_interpolator_1d->module~sll_m_spline_matrix module~sll_m_spline_interpolator_1d->module~sll_m_assert module~sll_m_boundary_condition_descriptors sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_bsplines_base module~sll_m_spline_interpolator_1d->module~sll_m_errors module~sll_m_spline_1d sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_working_precision module~sll_m_spline_matrix->module~sll_m_errors module~sll_m_spline_matrix_banded sll_m_spline_matrix_banded module~sll_m_spline_matrix->module~sll_m_spline_matrix_banded module~sll_m_spline_matrix_base sll_m_spline_matrix_base module~sll_m_spline_matrix->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_working_precision module~modgvec_io_netcdf->module~modgvec_globals module~sll_m_boundary_condition_descriptors->module~sll_m_working_precision module~sll_m_bsplines_base->module~sll_m_assert module~sll_m_bsplines_base->module~sll_m_working_precision module~sll_m_bsplines_non_uniform->module~sll_m_assert module~sll_m_bsplines_non_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform->module~sll_m_working_precision module~sll_m_bsplines_uniform->module~sll_m_assert module~sll_m_bsplines_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_uniform->module~sll_m_errors module~sll_m_bsplines_uniform->module~sll_m_working_precision module~sll_m_errors->iso_fortran_env module~sll_m_spline_1d->module~sll_m_assert module~sll_m_spline_1d->module~sll_m_bsplines_base module~sll_m_spline_1d->module~sll_m_working_precision module~sll_m_spline_matrix_banded->iso_fortran_env module~sll_m_spline_matrix_banded->module~sll_m_assert module~sll_m_spline_matrix_banded->module~sll_m_errors module~sll_m_spline_matrix_banded->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_banded->module~sll_m_working_precision module~sll_m_spline_matrix_base->module~sll_m_working_precision module~sll_m_spline_matrix_dense->iso_fortran_env module~sll_m_spline_matrix_dense->module~sll_m_assert module~sll_m_spline_matrix_dense->module~sll_m_errors module~sll_m_spline_matrix_dense->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense->module~sll_m_working_precision

prepare all data to be written

spos = s_pos(i_s) SCALE DOMAIN TO s=s_logicals_max ! SCALED DOMAIN: NEW LOGICAL DOMAIN REMAINS s_logical=[0,1], ! --> position to evaluated was s = s_logicals_max , d/ds_logical = ds/ds_logical * d/ds = s_max *d/ds

! data_scalar3D(ithet,izeta,i_s, S__) = spos =s_new, scale s-derivative with new domain size d/ds_new=d/ds*ds/ds_new scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size scale s-derivative with new domain size

Arguments

Type IntentOptional Attributes Name
class(t_base), intent(in) :: X1_base_in
real(kind=wp), intent(in) :: X1_in(1:X1_base_in%s%nBase,1:X1_base_in%f%modes)
class(t_base), intent(in) :: X2_base_in
real(kind=wp), intent(in) :: X2_in(1:X2_base_in%s%nBase,1:X2_base_in%f%modes)
class(t_base), intent(in) :: LG_base_in
real(kind=wp), intent(in) :: LG_in(1:LG_base_in%s%nBase,1:LG_base_in%f%modes)

Calls

proc~~gvec_to_jorek_prepare~~CallsGraph proc~gvec_to_jorek_prepare gvec_to_jorek_prepare eval_dxdq eval_dxdq proc~gvec_to_jorek_prepare->eval_dxdq interface~cross CROSS proc~gvec_to_jorek_prepare->interface~cross interface~progressbar ProgressBar proc~gvec_to_jorek_prepare->interface~progressbar proc~fbase_evaldof_x t_fBase%fBase_evalDOF_x proc~gvec_to_jorek_prepare->proc~fbase_evaldof_x proc~fbase_initdof t_fBase%fBase_initDOF proc~gvec_to_jorek_prepare->proc~fbase_initdof proc~get_field Get_Field proc~gvec_to_jorek_prepare->proc~get_field proc~sbase_evaldof2d_s t_sBase%sBase_evalDOF2D_s proc~gvec_to_jorek_prepare->proc~sbase_evaldof2d_s proc~sbase_evaldof_s t_sBase%sBase_evalDOF_s proc~gvec_to_jorek_prepare->proc~sbase_evaldof_s interface~cross->interface~cross interface~progressbar->interface~progressbar proc~fbase_eval t_fBase%fBase_eval proc~fbase_evaldof_x->proc~fbase_eval proc~fbase_projectiptodof_tens t_fBase%fBase_projectIPtoDOF_tens proc~fbase_initdof->proc~fbase_projectiptodof_tens proc~fbase_projectxntodof t_fBase%fBase_projectxntoDOF proc~fbase_initdof->proc~fbase_projectxntodof proc~get_field->eval_dxdq proc~get_field->interface~cross proc~get_field->proc~fbase_evaldof_x proc~get_field->proc~fbase_initdof proc~get_field->proc~sbase_evaldof2d_s proc~get_field->proc~sbase_evaldof_s proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens proc~get_field->proc~fbase_evaldof_ip_tens proc~sbase_applybctodof_lgm t_sBase%sBase_applyBCtoDOF_LGM proc~get_field->proc~sbase_applybctodof_lgm dgemv dgemv proc~sbase_evaldof2d_s->dgemv proc~sbase_eval t_sBase%sBase_eval proc~sbase_evaldof2d_s->proc~sbase_eval proc~sbase_evaldof_s->proc~sbase_eval proc~sbase_evaldof_base t_sBase%sBase_evalDOF_base proc~sbase_evaldof_s->proc~sbase_evaldof_base proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_eval->proc~fbase_eval_xn dgemm dgemm proc~fbase_evaldof_ip_tens->dgemm proc~fbase_evaldof_xn t_fBase%fBase_evalDOF_xn proc~fbase_evaldof_ip_tens->proc~fbase_evaldof_xn proc~fbase_projectiptodof_tens->dgemm proc~fbase_projectxntodof->dgemv proc~fbase_projectxntodof->proc~fbase_eval_xn proc~solve SOLVE proc~sbase_applybctodof_lgm->proc~solve eval_basis eval_basis proc~sbase_eval->eval_basis eval_basis_and_n_derivs eval_basis_and_n_derivs proc~sbase_eval->eval_basis_and_n_derivs lagrangeinterpolationpolys lagrangeinterpolationpolys proc~sbase_eval->lagrangeinterpolationpolys proc~sgrid_find_elem t_sGrid%sGrid_find_elem proc~sbase_eval->proc~sgrid_find_elem proc~fbase_evaldof_xn->dgemv proc~fbase_evaldof_xn->proc~fbase_eval_xn dgetrf dgetrf proc~solve->dgetrf dgetrs dgetrs proc~solve->dgetrs

Source Code

SUBROUTINE gvec_to_jorek_prepare(X1_base_in,X1_in,X2_base_in,X2_in,LG_base_in,LG_in)
! MODULES
USE MODgvec_gvec_to_jorek_Vars
USE MODgvec_Globals,        ONLY: CROSS,TWOPI,ProgressBar
USE MODgvec_ReadState_Vars, ONLY: profiles_1d,hmap_r,sbase_prof !for profiles
USE MODgvec_Base,           ONLY: t_base
USE MODgvec_fBase,          ONLY: t_fbase


IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
CLASS(t_base) ,INTENT(IN) :: X1_base_in,X2_base_in,LG_base_in
REAL(wp)      ,INTENT(IN) :: X1_in(1:X1_base_in%s%nBase,1:X1_base_in%f%modes)
REAL(wp)      ,INTENT(IN) :: X2_in(1:X2_base_in%s%nBase,1:X2_base_in%f%modes)
REAL(wp)      ,INTENT(IN) :: LG_in(1:LG_base_in%s%nBase,1:LG_base_in%f%modes) ! is either LA if SFLcoord=0/1 or G if SFLcoord=2
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER                                  :: i_s,ithet,izeta,iVar,jVar                                        ! Indices for enumeration
REAL(wp)                                 :: spos,xp(2),sqrtG                                                 ! Current s, theta, zeta, position and metric tensor
REAL(wp)                                 :: dX1ds,dX1dthet,dX1dzeta,d2X1dsdthet                              ! X coordinate and derivatives
REAL(wp)                                 :: dX2ds,dX2dthet,dX2dzeta,d2X2dsdthet                              ! Z coordinate and derivatives
REAL(wp)                                 :: dLAdthet,dLAdzeta                                                ! SFL transformation
REAL(wp)                                 :: Phi_int,dPhids_int,Chi_int,dChids_int,iota_int                   ! Flux variables and derivatives
REAL(wp)                                 :: P_int,dPds_int                                                   ! Pressure and derivatives
REAL(wp)                                 :: A_R_int, dA_Rds_int, dA_Rdthet_int, d2A_Rdsdthet_int             ! Vector potential components and derivatives
REAL(wp)                                 :: A_Z_int, dA_Zds_int, dA_Zdthet_int, d2A_Zdsdthet_int
REAL(wp)                                 :: A_phi_int, dA_phids_int, dA_phidthet_int, d2A_phidsdthet_int
REAL(wp)                                 :: B_R_int, dB_Rds_int, dB_Rdthet_int, d2B_Rdsdthet_int             ! Magnetic field components and derivatives
REAL(wp)                                 :: B_Z_int, dB_Zds_int, dB_Zdthet_int, d2B_Zdsdthet_int
REAL(wp)                                 :: B_phi_int, dB_phids_int, dB_phidthet_int, d2B_phidsdthet_int
REAL(wp)                                 :: J_R_int, dJ_Rds_int, dJ_Rdthet_int, d2J_Rdsdthet_int             ! Current density components and derivatives
REAL(wp)                                 :: J_Z_int, dJ_Zds_int, dJ_Zdthet_int, d2J_Zdsdthet_int
REAL(wp)                                 :: J_phi_int, dJ_phids_int, dJ_phidthet_int, d2J_phidsdthet_int
REAL(wp)                                 :: X1_int,X2_int,G_int,dGds,dGdthet,dGdzeta
REAL(wp)                                 :: AR_diff,AZ_diff,Aphi_diff,BR_diff,BZ_diff,Bphi_diff              ! Diagnostics for the convergence of the generated field representation
REAL(wp)                                 :: AR_diff_max,AZ_diff_max,Aphi_diff_max,BR_diff_max,BZ_diff_max,Bphi_diff_max
REAL(wp),DIMENSION(3)                    :: qvec,e_s,e_thet,e_zeta                                           ! Vectors for local covariant coordinate system
REAL(wp),DIMENSION(3)                    :: Acart,A_orig,Bcart,B_orig                                        ! Test variables for magnetic field and position
REAL(wp)                                 :: Bthet, Bzeta
REAL(wp),DIMENSION(3)                    :: grad_s,grad_thet,grad_zeta,grad_R,grad_Z                         ! Contravariant coordinates

! 2D theta/zeta fourier representation of variables in GVEC
REAL(wp),DIMENSION(1:X1_base_in%f%modes) :: X1_s,dX1ds_s
REAL(wp),DIMENSION(1:X2_base_in%f%modes) :: X2_s,dX2ds_s
REAL(wp),DIMENSION(1:out_base%f%modes)   :: A_R_s, A_Rds_int, A_Z_s, A_Zds_int, A_phi_s, A_phids_int
REAL(wp),DIMENSION(1:out_base%f%modes)   :: B_R_s, B_Rds_int, B_Z_s, B_Zds_int, B_phi_s, B_phids_int
REAL(wp),DIMENSION(1:out_base%f%modes)   :: J_R_s, J_Rds_int, J_Z_s, J_Zds_int, J_phi_s, J_phids_int
REAL(wp),DIMENSION(1:LG_base_in%f%modes) :: LG_s


! Variables for new GVEC fourier representations needed for JOREK inputs
REAL(wp),DIMENSION(out_base%s%nBase,out_base%f%modes):: A_R, A_Z, A_phi !< data (1:nBase,1:modes) of A_* in GVEC coords
REAL(wp),DIMENSION(out_base%s%nBase,out_base%f%modes):: B_R, B_Z, B_phi !< data (1:nBase,1:modes) of B_* in GVEC coords
REAL(wp),DIMENSION(out_base%s%nBase,out_base%f%modes):: J_R, J_Z, J_phi !< data (1:nBase,1:modes) of J_* in GVEC coords

!===================================================================================================================================
! ------------------------------------------------------------------------------------------------------------
! ------- CALCULATE 2D FOURIER REPRESENTATION OF Vector Potential, Magnetic Field, and Current Density -------
! ------------------------------------------------------------------------------------------------------------
SWRITE(UNIT_stdOut,'(A)')'PREPARE FIELD BASES ...'
CALL get_field(2,4,A_R)
CALL get_field(2,3,A_Z)
CALL get_field(2,5,A_phi)
CALL get_field(1,4,B_R)
CALL get_field(1,3,B_Z)
CALL get_field(1,5,B_phi)
CALL get_field(3,4,J_R)
CALL get_field(3,3,J_Z)
CALL get_field(3,5,J_phi)

! -----------------------------------------------------------------------------
! ---------- GENERATE 3D POINTS FROM GVEC REPRESENTATION ----------------------
! -----------------------------------------------------------------------------
SWRITE(UNIT_stdOut,'(A)')'PREPARE 3D DATA FOR GVEC-TO-JOREK ...'
CALL ProgressBar(0,Ns_out) !init
DO i_s=1,Ns_out
  !!spos  = s_pos(i_s)
  !! SCALE DOMAIN TO s=s_logical*s_max
  spos = s_pos(i_s)*s_max

  Phi_int     = sbase_prof%evalDOF_s(spos,       0 ,profiles_1d(:,1))
  dPhids_int  = sbase_prof%evalDOF_s(spos, DERIV_S ,profiles_1d(:,1))
  Chi_int     = sbase_prof%evalDOF_s(spos,       0 ,profiles_1d(:,2)) !Chi representation inconsistent - see ReadState routine
  dChids_int  = sbase_prof%evalDOF_s(spos, DERIV_S ,profiles_1d(:,2)) !Chi representation inconsistent - see ReadState routine
  iota_int    = sbase_prof%evalDOF_s(spos,       0 ,profiles_1d(:,3))
  P_int       = sbase_prof%evalDOF_s(spos,       0 ,profiles_1d(:,4))
  dPds_int    = sbase_prof%evalDOF_s(spos, DERIV_S ,profiles_1d(:,4))

  dLAdthet= 0.0_wp !only changed for SFLcoord=0
  dLAdzeta= 0.0_wp !only changed for SFLcoord=0
  G_int   = 0.0_wp !only changed for SFLcoords=2
  dGds    = 0.0_wp !only changed for SFLcoords=2
  dGdthet = 0.0_wp !only changed for SFLcoords=2
  dGdzeta = 0.0_wp !only changed for SFLcoords=2

  !interpolate radially
  X1_s(   :) = X1_base_in%s%evalDOF2D_s(spos,X1_base_in%f%modes,       0,X1_in(:,:))
  dX1ds_s(:) = X1_base_in%s%evalDOF2D_s(spos,X1_base_in%f%modes, DERIV_S,X1_in(:,:))

  X2_s(   :) = X2_base_in%s%evalDOF2D_s(spos,X2_base_in%f%modes,       0,X2_in(:,:))
  dX2ds_s(:) = X2_base_in%s%evalDOF2D_s(spos,X2_base_in%f%modes, DERIV_S,X2_in(:,:))
  IF(SFLcoord.EQ.0)THEN !GVEC coordinates
    LG_s(  :) = LG_base_in%s%evalDOF2D_s(spos,LG_base_in%f%modes,       0,LG_in(:,:))
  ELSE
    SWRITE(UNIT_StdOut,*)'This SFLcoord is not valid',SFLcoord
    STOP
  END IF

  ! Generate representations of the fields in (R,Z,phi)
  A_R_s(:)       =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,         0, A_R(:,:))
  A_Rds_int(:)   =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,   DERIV_S, A_R(:,:))
  A_Z_s(:)       =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,         0, A_Z(:,:))
  A_Zds_int(:)   =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,   DERIV_S, A_Z(:,:))
  A_phi_s(:)     =   out_base%s%evalDOF2D_s(spos, out_base%f%modes,         0, A_phi(:,:))
  A_phids_int(:) =   out_base%s%evalDOF2D_s(spos, out_base%f%modes,   DERIV_S, A_phi(:,:))
  B_R_s(:)       =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,         0, B_R(:,:))
  B_Rds_int(:)   =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,   DERIV_S, B_R(:,:))
  B_Z_s(:)       =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,         0, B_Z(:,:))
  B_Zds_int(:)   =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,   DERIV_S, B_Z(:,:))
  B_phi_s(:)     =   out_base%s%evalDOF2D_s(spos, out_base%f%modes,         0, B_phi(:,:))
  B_phids_int(:) =   out_base%s%evalDOF2D_s(spos, out_base%f%modes,   DERIV_S, B_phi(:,:))
  J_R_s(:)       =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,         0, J_R(:,:))
  J_Rds_int(:)   =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,   DERIV_S, J_R(:,:))
  J_Z_s(:)       =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,         0, J_Z(:,:))
  J_Zds_int(:)   =   out_base%s%evalDOF2D_s(spos,   out_base%f%modes,   DERIV_S, J_Z(:,:))
  J_phi_s(:)     =   out_base%s%evalDOF2D_s(spos, out_base%f%modes,         0, J_phi(:,:))
  J_phids_int(:) =   out_base%s%evalDOF2D_s(spos, out_base%f%modes,   DERIV_S, J_phi(:,:))

  BR_diff=0.0_wp; BZ_diff=0.0_wp; Bphi_diff=0.0_wp; AR_diff=0.0_wp; AZ_diff=0.0_wp; Aphi_diff=0.0_wp
  BR_diff_max=0.0_wp; BZ_diff_max=0.0_wp; Bphi_diff_max=0.0_wp; AR_diff_max=0.0_wp; AZ_diff_max=0.0_wp; Aphi_diff_max=0.0_wp
!$OMP PARALLEL DO  SCHEDULE(STATIC) DEFAULT(NONE) COLLAPSE(2)                                                                  &
!$OMP   REDUCTION(+:BR_diff, BZ_diff, Bphi_diff, AR_diff, AZ_diff, Aphi_diff)                                                  &
!$OMP   REDUCTION(max: BR_diff_max, BZ_diff_max, Bphi_diff_max, AR_diff_max,AZ_diff_max,Aphi_diff_max)                         &
!$OMP   PRIVATE(izeta,ithet,X1_int,dX1ds,dX1dthet,dX1dzeta,d2X1dsdthet, X2_int,dX2ds,dX2dthet,dX2dzeta, d2X2dsdthet,           &
!$OMP           xp,qvec,e_s,e_thet,e_zeta,sqrtG,                                                                               &
!$OMP           A_R_int,dA_Rds_int,dA_Rdthet_int,d2A_Rdsdthet_int,                                                             &
!$OMP           A_Z_int,dA_Zds_int,dA_Zdthet_int,d2A_Zdsdthet_int,                                                             &
!$OMP           A_phi_int,dA_phids_int,dA_phidthet_int,d2A_phidsdthet_int,                                                     &
!$OMP           B_R_int,dB_Rds_int,dB_Rdthet_int,d2B_Rdsdthet_int,                                                             &
!$OMP           B_Z_int,dB_Zds_int,dB_Zdthet_int,d2B_Zdsdthet_int,                                                             &
!$OMP           B_phi_int,dB_phids_int,dB_phidthet_int,d2B_phidsdthet_int,                                                     &
!$OMP           J_R_int,dJ_Rds_int,dJ_Rdthet_int,d2J_Rdsdthet_int,                                                             &
!$OMP           J_Z_int,dJ_Zds_int,dJ_Zdthet_int,d2J_Zdsdthet_int,                                                             &
!$OMP           J_phi_int,dJ_phids_int,dJ_phidthet_int,d2J_phidsdthet_int,                                                     &
!$OMP           Acart,A_orig,grad_s,grad_thet,grad_zeta,grad_R,grad_Z,                                                   &
!$OMP           Bthet,Bzeta,Bcart,B_orig)                                                                                      &
!$OMP   FIRSTPRIVATE(dLAdthet,dLAdzeta,G_int,dGds,dGdthet,dGdzeta)                                                             &
!$OMP   SHARED(i_s,Nzeta_out,Nthet_out,spos,s_pos,s_max, thet_pos,zeta_pos,X1_base_in,X2_base_in,LG_base_in,out_base,               &
!$OMP          hmap_r,X1_s,dX1ds_s,X2_s,dX2ds_s,LG_s,SFLcoord, Phi_int, dPhids_int, iota_int, Chi_int, dChids_int, P_int, dPds_int, &
!$OMP          A_R_s, A_Rds_int,A_Z_s, A_Zds_int,A_phi_s, A_phids_int, B_R_s, B_Rds_int,B_Z_s, B_Zds_int,B_phi_s, B_phids_int,      &
!$OMP          J_R_s, J_Rds_int,J_Z_s, J_Zds_int,J_phi_s, J_phids_int,                                                              &
!$OMP          data_scalar3D)
  !interpolate in the angles
  DO izeta=1,Nzeta_out; DO ithet=1,Nthet_out
    xp=(/TWOPI * thet_pos(ithet),zeta_pos(izeta)/)

    X1_int      = X1_base_in%f%evalDOF_x(xp,          0, X1_s  )
    dX1ds       = X1_base_in%f%evalDOF_x(xp,          0,dX1ds_s)
    dX1dthet    = X1_base_in%f%evalDOF_x(xp, DERIV_THET, X1_s  )
    d2X1dsdthet = X1_base_in%f%evalDOF_x(xp, DERIV_THET,dX1ds_s)
    dX1dzeta    = X1_base_in%f%evalDOF_x(xp, DERIV_ZETA, X1_s  )

    X2_int      = X2_base_in%f%evalDOF_x(xp,          0, X2_s  )
    dX2ds       = X2_base_in%f%evalDOF_x(xp,          0,dX2ds_s)
    dX2dthet    = X2_base_in%f%evalDOF_x(xp, DERIV_THET, X2_s  )
    d2X2dsdthet = X2_base_in%f%evalDOF_x(xp, DERIV_THET,dX2ds_s)
    dX2dzeta    = X2_base_in%f%evalDOF_x(xp, DERIV_ZETA, X2_s  )

    ! Get A components
    A_R_int            = out_base%f%evalDOF_x(xp,          0,   A_R_s)
    dA_Rds_int         = out_base%f%evalDOF_x(xp,          0, A_Rds_int)
    dA_Rdthet_int      = out_base%f%evalDOF_x(xp, DERIV_THET,   A_R_s)
    d2A_Rdsdthet_int   = out_base%f%evalDOF_x(xp, DERIV_THET, A_Rds_int)
    A_Z_int            = out_base%f%evalDOF_x(xp,          0,   A_Z_s)
    dA_Zds_int         = out_base%f%evalDOF_x(xp,          0, A_Zds_int)
    dA_Zdthet_int      = out_base%f%evalDOF_x(xp, DERIV_THET,   A_Z_s)
    d2A_Zdsdthet_int   = out_base%f%evalDOF_x(xp, DERIV_THET, A_Zds_int)
    A_phi_int          = out_base%f%evalDOF_x(xp,          0,   A_phi_s)
    dA_phids_int       = out_base%f%evalDOF_x(xp,          0, A_phids_int)
    dA_phidthet_int    = out_base%f%evalDOF_x(xp, DERIV_THET,   A_phi_s)
    d2A_phidsdthet_int = out_base%f%evalDOF_x(xp, DERIV_THET, A_phids_int)

    ! Get B components
    B_R_int            = out_base%f%evalDOF_x(xp,          0,   B_R_s)
    dB_Rds_int         = out_base%f%evalDOF_x(xp,          0, B_Rds_int)
    dB_Rdthet_int      = out_base%f%evalDOF_x(xp, DERIV_THET,   B_R_s)
    d2B_Rdsdthet_int   = out_base%f%evalDOF_x(xp, DERIV_THET, B_Rds_int)
    B_Z_int            = out_base%f%evalDOF_x(xp,          0,   B_Z_s)
    dB_Zds_int         = out_base%f%evalDOF_x(xp,          0, B_Zds_int)
    dB_Zdthet_int      = out_base%f%evalDOF_x(xp, DERIV_THET,   B_Z_s)
    d2B_Zdsdthet_int   = out_base%f%evalDOF_x(xp, DERIV_THET, B_Zds_int)
    B_phi_int          = out_base%f%evalDOF_x(xp,          0,   B_phi_s)
    dB_phids_int       = out_base%f%evalDOF_x(xp,          0, B_phids_int)
    dB_phidthet_int    = out_base%f%evalDOF_x(xp, DERIV_THET,   B_phi_s)
    d2B_phidsdthet_int = out_base%f%evalDOF_x(xp, DERIV_THET, B_phids_int)

    ! Get J components
    J_R_int            = out_base%f%evalDOF_x(xp,          0,   J_R_s)
    dJ_Rds_int         = out_base%f%evalDOF_x(xp,          0, J_Rds_int)
    dJ_Rdthet_int      = out_base%f%evalDOF_x(xp, DERIV_THET,   J_R_s)
    d2J_Rdsdthet_int   = out_base%f%evalDOF_x(xp, DERIV_THET, J_Rds_int)
    J_Z_int            = out_base%f%evalDOF_x(xp,          0,   J_Z_s)
    dJ_Zds_int         = out_base%f%evalDOF_x(xp,          0, J_Zds_int)
    dJ_Zdthet_int      = out_base%f%evalDOF_x(xp, DERIV_THET,   J_Z_s)
    d2J_Zdsdthet_int   = out_base%f%evalDOF_x(xp, DERIV_THET, J_Zds_int)
    J_phi_int          = out_base%f%evalDOF_x(xp,          0,   J_phi_s)
    dJ_phids_int       = out_base%f%evalDOF_x(xp,          0, J_phids_int)
    dJ_phidthet_int    = out_base%f%evalDOF_x(xp, DERIV_THET,   J_phi_s)
    d2J_phidsdthet_int = out_base%f%evalDOF_x(xp, DERIV_THET, J_phids_int)

    ! Get straight field line transformation
    IF(SFLcoord.EQ.0)THEN !GVEC coordinates (else=0)
      G_int    = LG_base_in%f%evalDOF_x(xp, 0, LG_s)
      dLAdthet = LG_base_in%f%evalDOF_x(xp, DERIV_THET, LG_s)
      dLAdzeta = LG_base_in%f%evalDOF_x(xp, DERIV_ZETA, LG_s)
    END IF

    ! --- Compare generated field representations of A and B with calculations from
    ! --- the original representation to ensure convergence
    ! Get the covariant basis vectors
    qvec     = (/ X1_int, X2_int, xp(2) /) !(X1,X2,zeta)
    e_s      = hmap_r%eval_dxdq(qvec,(/dX1ds   ,dX2ds   , 0.0    /)) !dxvec/ds
    e_thet   = hmap_r%eval_dxdq(qvec,(/dX1dthet,dX2dthet, 0.0    /)) !dxvec/dthet
    e_zeta   = hmap_r%eval_dxdq(qvec,(/dX1dzeta,dX2dzeta, 1.0_wp /)) !dxvec/dzeta
    sqrtG    = SUM(e_s*(CROSS(e_thet,e_zeta)))

    ! Get contravarian basis vectors
    grad_s    = CROSS(e_thet,e_zeta) /sqrtG
    grad_thet = CROSS(e_zeta,e_s   ) /sqrtG
    grad_zeta = CROSS(e_s   ,e_thet) /sqrtG

    ! Get Grad R and Grad Z pol - WARNING: this implementation only works for PEST coordinates
    grad_R = dX1ds * grad_s + dX1dthet * grad_thet + dX1dzeta * grad_zeta
    grad_Z = dX2ds * grad_s + dX2dthet * grad_thet + dX2dzeta * grad_zeta

    ! Calculate ave/max error in (R,Z, phi) magnetic field
    Bthet = ((iota_int-dLAdzeta )*dPhids_int)   !/sqrtG
    Bzeta = ((1.0_wp  +dLAdthet )*dPhids_int)   !/sqrtG
    Bcart(:) =  ( e_thet(:)*Bthet+e_zeta(:)*Bzeta) /sqrtG
    B_orig(1) =  Bcart(1) * grad_R(1)    + Bcart(2) * grad_R(2)    + Bcart(3) * grad_R(3)
    B_orig(2) =  Bcart(1) * grad_Z(1)    + Bcart(2) * grad_Z(2)    + Bcart(3) * grad_Z(3)
    B_orig(3) =  X1_int * (Bcart(1) * grad_zeta(1) + Bcart(2) * grad_zeta(2) + Bcart(3) * grad_zeta(3))
    BR_diff   = BR_diff   + ABS((B_orig(1) - B_R_int));   BR_diff_max   = MAX(BR_diff_max,   ABS(B_orig(1) - B_R_int))
    BZ_diff   = BZ_diff   + ABS((B_orig(2) - B_Z_int));   BZ_diff_max   = MAX(BZ_diff_max,   ABS(B_orig(2) - B_Z_int))
    Bphi_diff = Bphi_diff + ABS((B_orig(3) - B_phi_int)); Bphi_diff_max = MAX(Bphi_diff_max, ABS(B_orig(3) - B_phi_int))

    ! Calculate ave/max error in (R,Z,phi) vector potential
    Acart(:)  = (Phi_int * grad_thet(:) - (G_int * dPhids_int) * grad_s(:) - chi_int * grad_zeta(:))
    A_orig(1) =  Acart(1) * grad_R(1)    + Acart(2) * grad_R(2)    + Acart(3) * grad_R(3)
    A_orig(2) =  Acart(1) * grad_Z(1)    + Acart(2) * grad_Z(2)    + Acart(3) * grad_Z(3)
    A_orig(3) =  X1_int * (Acart(1) * grad_zeta(1) + Acart(2) * grad_zeta(2) + Acart(3) * grad_zeta(3))
    AR_diff   = AR_diff   + ABS((A_orig(1) - A_R_int));   AR_diff_max   = MAX(AR_diff_max,   ABS(A_orig(1) - A_R_int))
    AZ_diff   = AZ_diff   + ABS((A_orig(2) - A_Z_int));   AZ_diff_max   = MAX(AZ_diff_max,   ABS(A_orig(2) - A_Z_int))
    Aphi_diff = Aphi_diff + ABS((A_orig(3) - A_phi_int)); Aphi_diff_max = MAX(Aphi_diff_max, ABS(A_orig(3) - A_phi_int))

    !==========
    ! save data

    !!! SCALED DOMAIN: NEW LOGICAL DOMAIN REMAINS s_logical=[0,1],
    !!!                -->  position to evaluated was s = s_logical*s_max ,  d/ds_logical = ds/ds_logical * d/ds = s_max *d/ds

    !!! data_scalar3D(ithet,izeta,i_s, S__)          = spos
    data_scalar3D(ithet,izeta,i_s, S__)          = s_pos(i_s) !! =s_new,
    data_scalar3D(ithet,izeta,i_s, THET__)       = thet_pos(ithet)
    data_scalar3D(ithet,izeta,i_s, ZETA__)       = zeta_pos(izeta)
    data_scalar3D(ithet,izeta,i_s, X1__)         = X1_int
    data_scalar3D(ithet,izeta,i_s, X1_S__)       = dX1ds *s_max !! scale s-derivative with new domain size d/ds_new=d/ds*ds/ds_new
    data_scalar3D(ithet,izeta,i_s, X1_T__)       = dX1dthet
    data_scalar3D(ithet,izeta,i_s, X1_ST__)      = d2X1dsdthet *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, X2__)         = X2_int
    data_scalar3D(ithet,izeta,i_s, X2_S__)       = dX2ds *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, X2_T__)       = dX2dthet
    data_scalar3D(ithet,izeta,i_s, X2_ST__)      = d2X2dsdthet *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, P__)          = P_int
    data_scalar3D(ithet,izeta,i_s, P_S__)        = dPds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, A_R__)        = A_R_int
    data_scalar3D(ithet,izeta,i_s, A_R_S__)      = dA_Rds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, A_R_T__)      = dA_Rdthet_int
    data_scalar3D(ithet,izeta,i_s, A_R_ST__)     = d2A_Rdsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, A_Z__)        = A_Z_int
    data_scalar3D(ithet,izeta,i_s, A_Z_S__)      = dA_Zds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, A_Z_T__)      = dA_Zdthet_int
    data_scalar3D(ithet,izeta,i_s, A_Z_ST__)     = d2A_Zdsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, A_phi__)      = A_phi_int
    data_scalar3D(ithet,izeta,i_s, A_phi_S__)    = dA_phids_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, A_phi_T__)    = dA_phidthet_int
    data_scalar3D(ithet,izeta,i_s, A_phi_ST__)   = d2A_phidsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, B_R__)        = B_R_int
    data_scalar3D(ithet,izeta,i_s, B_R_S__)      = dB_Rds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, B_R_T__)      = dB_Rdthet_int
    data_scalar3D(ithet,izeta,i_s, B_R_ST__)     = d2B_Rdsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, B_Z__)        = B_Z_int
    data_scalar3D(ithet,izeta,i_s, B_Z_S__)      = dB_Zds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, B_Z_T__)      = dB_Zdthet_int
    data_scalar3D(ithet,izeta,i_s, B_Z_ST__)     = d2B_Zdsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, B_phi__)      = B_phi_int
    data_scalar3D(ithet,izeta,i_s, B_phi_S__)    = dB_phids_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, B_phi_T__)    = dB_phidthet_int
    data_scalar3D(ithet,izeta,i_s, B_phi_ST__)   = d2B_phidsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, J_R__)        = J_R_int
    data_scalar3D(ithet,izeta,i_s, J_R_S__)      = dJ_Rds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, J_R_T__)      = dJ_Rdthet_int
    data_scalar3D(ithet,izeta,i_s, J_R_ST__)     = d2J_Rdsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, J_Z__)        = J_Z_int
    data_scalar3D(ithet,izeta,i_s, J_Z_S__)      = dJ_Zds_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, J_Z_T__)      = dJ_Zdthet_int
    data_scalar3D(ithet,izeta,i_s, J_Z_ST__)     = d2J_Zdsdthet_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, J_phi__)      = J_phi_int
    data_scalar3D(ithet,izeta,i_s, J_phi_S__)    = dJ_phids_int *s_max !! scale s-derivative with new domain size
    data_scalar3D(ithet,izeta,i_s, J_phi_T__)    = dJ_phidthet_int
    data_scalar3D(ithet,izeta,i_s, J_phi_ST__)   = d2J_phidsdthet_int *s_max !! scale s-derivative with new domain size
    !==========

  END DO ; END DO !izeta,ithet
!$OMP END PARALLEL DO
  CALL ProgressBar(i_s,Ns_out)
END DO !i_s=1,Ns_out
BR_diff = BR_diff / REAL(Ns_out*Nthet_out*Nzeta_out,wp);BZ_diff = BZ_diff / REAL(Ns_out*Nthet_out*Nzeta_out,wp);Bphi_diff = Bphi_diff / REAL(Ns_out*Nthet_out*Nzeta_out,wp)
AR_diff = AR_diff / REAL(Ns_out*Nthet_out*Nzeta_out,wp);AZ_diff = AZ_diff / REAL(Ns_out*Nthet_out*Nzeta_out,wp);Aphi_diff = Aphi_diff / REAL(Ns_out*Nthet_out*Nzeta_out,wp)
SWRITE(UNIT_stdOut,'(A,6E16.7)')'AVE/MAX diff in B(R, Z, phi) :', BR_diff, BR_diff_max, BZ_diff, BZ_diff_max, Bphi_diff, Bphi_diff_max
SWRITE(UNIT_stdOut,'(A,6E16.7)')'MIN/MAX B(R, Z, phi)         :',MINVAL(data_scalar3D(:,:,:,B_R__)),MAXVAL(data_scalar3D(:,:,:,B_R__)),&
                                                                 MINVAL(data_scalar3D(:,:,:,B_Z__)),MAXVAL(data_scalar3D(:,:,:,B_Z__)),&
                                                                 MINVAL(data_scalar3D(:,:,:,B_phi__)),MAXVAL(data_scalar3D(:,:,:,B_phi__))
SWRITE(UNIT_stdOut,'(A,6E16.7)')'AVE/MAX diff in A(R, Z, phi) :', AR_diff,AR_diff_max,AZ_diff,AZ_diff_max,Aphi_diff,Aphi_diff_max
SWRITE(UNIT_stdOut,'(A,6E16.7)')'MIN/MAX A(R, Z, phi)         :',MINVAL(data_scalar3D(:,:,:,A_R__)),MAXVAL(data_scalar3D(:,:,:,A_R__)),&
                                                                 MINVAL(data_scalar3D(:,:,:,A_Z__)),MAXVAL(data_scalar3D(:,:,:,A_Z__)),&
                                                                 MINVAL(data_scalar3D(:,:,:,A_phi__)),MAXVAL(data_scalar3D(:,:,:,A_phi__))

! -----------------------------------------------------------------------------
! ---------------- CONVERT TO 1D TOROIDAL REPRESENTATION ----------------------
! -----------------------------------------------------------------------------
SWRITE(Unit_stdOut, *)  "Writing 3D variables to toroidal fourier representation..."
map_vars_3D_2D(R__     )      = X1__
map_vars_3D_2D(R_S__   )      = X1_S__
map_vars_3D_2D(R_T__   )      = X1_T__
map_vars_3D_2D(R_ST__  )      = X1_ST__
map_vars_3D_2D(Z__     )      = X2__
map_vars_3D_2D(Z_S__   )      = X2_S__
map_vars_3D_2D(Z_T__   )      = X2_T__
map_vars_3D_2D(Z_ST__  )      = X2_ST__
map_vars_3D_2D(P2D__   )      = P__
map_vars_3D_2D(P2D_S__ )      = P_S__
map_vars_3D_2D(P2D_T__ )      = -1 !leave zero
map_vars_3D_2D(P2D_ST__)      = -1 !leave zero
map_vars_3D_2D(A_R2D__ )      =A_R__
map_vars_3D_2D(A_R2D_S__)     =A_R_S__
map_vars_3D_2D(A_R2D_T__)     =A_R_T__
map_vars_3D_2D(A_R2D_ST__)    =A_R_ST__
map_vars_3D_2D(A_Z2D__)       =A_Z__
map_vars_3D_2D(A_Z2D_S__)     =A_Z_S__
map_vars_3D_2D(A_Z2D_T__)     =A_Z_T__
map_vars_3D_2D(A_Z2D_ST__)    =A_Z_ST__
map_vars_3D_2D(A_phi2D__)     =A_phi__
map_vars_3D_2D(A_phi2D_S__)   =A_phi_S__
map_vars_3D_2D(A_phi2D_T__)   =A_phi_T__
map_vars_3D_2D(A_phi2D_ST__)  =A_phi_ST__
map_vars_3D_2D(B_R2D__)       =B_R__
map_vars_3D_2D(B_R2D_S__)     =B_R_S__
map_vars_3D_2D(B_R2D_T__)     =B_R_T__
map_vars_3D_2D(B_R2D_ST__)    =B_R_ST__
map_vars_3D_2D(B_Z2D__)       =B_Z__
map_vars_3D_2D(B_Z2D_S__)     =B_Z_S__
map_vars_3D_2D(B_Z2D_T__)     =B_Z_T__
map_vars_3D_2D(B_Z2D_ST__)    =B_Z_ST__
map_vars_3D_2D(B_phi2D__)     =B_phi__
map_vars_3D_2D(B_phi2D_S__)   =B_phi_S__
map_vars_3D_2D(B_phi2D_T__)   =B_phi_T__
map_vars_3D_2D(B_phi2D_ST__)  =B_phi_ST__
map_vars_3D_2D(J_R2D__)       =J_R__
map_vars_3D_2D(J_R2D_S__)     =J_R_S__
map_vars_3D_2D(J_R2D_T__)     =J_R_T__
map_vars_3D_2D(J_R2D_ST__)    =J_R_ST__
map_vars_3D_2D(J_Z2D__)       =J_Z__
map_vars_3D_2D(J_Z2D_S__)     =J_Z_S__
map_vars_3D_2D(J_Z2D_T__)     =J_Z_T__
map_vars_3D_2D(J_Z2D_ST__)    =J_Z_ST__
map_vars_3D_2D(J_phi2D__)     =J_phi__
map_vars_3D_2D(J_phi2D_S__)   =J_phi_S__
map_vars_3D_2D(J_phi2D_T__)   =J_phi_T__
map_vars_3D_2D(J_phi2D_ST__)  =J_phi_ST__

data_scalar2D=0.0_wp
DO iVar=1,nVarScalar2D
  jVar=map_vars_3D_2D(iVar)
  IF(jVar.LE.0)CYCLE !do not set these
  DO ithet=1, Nthet_out
    DO i_s=1, Ns_out
      ! Store DOFs for writing to file
      data_scalar2D(ithet, i_s, 1:n_modes,iVar) = fbase_zeta%initDOF(data_scalar3D(ithet,:,i_s,jVar))
    END DO ! ithet=1, Nthet_out
  END DO ! i_s=1, Ns_out
END DO !iVar=1,nVarScalar2D


SWRITE(UNIT_stdOut,'(A)')'... DONE'
SWRITE(UNIT_stdOut,fmt_sep)

END SUBROUTINE gvec_to_jorek_prepare