InitMHD3D Subroutine

private subroutine InitMHD3D(sf)

Uses

  • proc~~initmhd3d~~UsesGraph proc~initmhd3d t_functional_mhd3d%InitMHD3D module~modgvec_base MODgvec_base proc~initmhd3d->module~modgvec_base module~modgvec_boundaryfromfile MODgvec_boundaryFromFile proc~initmhd3d->module~modgvec_boundaryfromfile module~modgvec_globals MODgvec_Globals proc~initmhd3d->module~modgvec_globals module~modgvec_hmap MODgvec_hmap proc~initmhd3d->module~modgvec_hmap module~modgvec_mhd3d_evalfunc MODgvec_MHD3D_evalFunc proc~initmhd3d->module~modgvec_mhd3d_evalfunc module~modgvec_mhd3d_minimize MODgvec_MHD3D_minimize proc~initmhd3d->module~modgvec_mhd3d_minimize module~modgvec_mhd3d_vars MODgvec_MHD3D_Vars proc~initmhd3d->module~modgvec_mhd3d_vars module~modgvec_mpi MODgvec_MPI proc~initmhd3d->module~modgvec_mpi module~modgvec_readintools MODgvec_ReadInTools proc~initmhd3d->module~modgvec_readintools module~modgvec_rprofile_poly MODgvec_rProfile_poly proc~initmhd3d->module~modgvec_rprofile_poly module~modgvec_sgrid MODgvec_sGrid proc~initmhd3d->module~modgvec_sgrid module~modgvec_vmec MODgvec_VMEC proc~initmhd3d->module~modgvec_vmec module~modgvec_vmec_readin MODgvec_VMEC_Readin proc~initmhd3d->module~modgvec_vmec_readin module~modgvec_vmec_vars MODgvec_VMEC_Vars proc~initmhd3d->module~modgvec_vmec_vars module~modgvec_base->module~modgvec_globals module~modgvec_base->module~modgvec_sgrid module~modgvec_fbase MODgvec_fBase module~modgvec_base->module~modgvec_fbase module~modgvec_sbase MODgvec_sBase module~modgvec_base->module~modgvec_sbase module~modgvec_boundaryfromfile->module~modgvec_globals module~modgvec_io_netcdf MODgvec_IO_NETCDF module~modgvec_boundaryfromfile->module~modgvec_io_netcdf iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env 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_mhd3d_evalfunc->module~modgvec_globals module~sll_m_spline_matrix sll_m_spline_matrix module~modgvec_mhd3d_evalfunc->module~sll_m_spline_matrix module~sll_m_spline_matrix_banded sll_m_spline_matrix_banded module~modgvec_mhd3d_evalfunc->module~sll_m_spline_matrix_banded module~modgvec_mhd3d_minimize->module~modgvec_globals module~modgvec_sol_var_mhd3d MODgvec_sol_var_MHD3D module~modgvec_mhd3d_minimize->module~modgvec_sol_var_mhd3d module~modgvec_mhd3d_vars->module~modgvec_base module~modgvec_mhd3d_vars->module~modgvec_boundaryfromfile module~modgvec_mhd3d_vars->module~modgvec_globals module~modgvec_mhd3d_vars->module~modgvec_hmap module~modgvec_mhd3d_vars->module~modgvec_sgrid module~modgvec_rprofile_base MODgvec_rProfile_base module~modgvec_mhd3d_vars->module~modgvec_rprofile_base module~modgvec_mhd3d_vars->module~modgvec_sol_var_mhd3d module~modgvec_readintools->module~modgvec_globals module~modgvec_rprofile_poly->module~modgvec_globals module~modgvec_rprofile_poly->module~modgvec_rprofile_base module~modgvec_sgrid->module~modgvec_globals module~modgvec_vmec->module~modgvec_globals module~modgvec_cubic_spline MODgvec_cubic_spline module~modgvec_vmec->module~modgvec_cubic_spline module~modgvec_vmec_readin->module~modgvec_globals module~modgvec_vmec_vars->module~modgvec_globals module~modgvec_vmec_vars->module~modgvec_cubic_spline module~modgvec_vmec_vars->module~modgvec_rprofile_base module~modgvec_c_hmap->module~modgvec_globals module~modgvec_cubic_spline->module~modgvec_globals module~sll_m_bsplines sll_m_bsplines module~modgvec_cubic_spline->module~sll_m_bsplines module~modgvec_fbase->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_globals module~modgvec_hmap_axisnb->module~modgvec_c_hmap module~modgvec_hmap_axisnb->module~modgvec_fbase 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~modgvec_io_netcdf->module~modgvec_globals module~modgvec_rprofile_base->module~modgvec_globals module~modgvec_sbase->module~modgvec_globals module~modgvec_sbase->module~modgvec_sgrid module~modgvec_sbase->module~sll_m_spline_matrix 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~modgvec_sol_var_mhd3d->module~modgvec_globals module~modgvec_c_sol_var MODgvec_c_sol_var module~modgvec_sol_var_mhd3d->module~modgvec_c_sol_var module~sll_m_spline_matrix->module~sll_m_spline_matrix_banded module~sll_m_errors sll_m_errors module~sll_m_spline_matrix->module~sll_m_errors 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_working_precision sll_m_working_precision module~sll_m_spline_matrix->module~sll_m_working_precision module~sll_m_spline_matrix_banded->iso_fortran_env module~sll_m_assert sll_m_assert 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~modgvec_c_sol_var->module~modgvec_globals module~sll_m_bsplines->module~sll_m_assert module~sll_m_bsplines->module~sll_m_errors module~sll_m_bsplines->module~sll_m_working_precision 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->iso_fortran_env module~sll_m_spline_interpolator_1d->module~sll_m_spline_matrix module~sll_m_spline_interpolator_1d->module~sll_m_assert module~sll_m_spline_interpolator_1d->module~sll_m_errors module~sll_m_spline_interpolator_1d->module~sll_m_working_precision 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_1d sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_spline_1d 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 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_working_precision module~sll_m_bsplines_non_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_uniform->module~sll_m_assert module~sll_m_bsplines_uniform->module~sll_m_errors module~sll_m_bsplines_uniform->module~sll_m_working_precision module~sll_m_bsplines_uniform->module~sll_m_bsplines_base module~sll_m_spline_1d->module~sll_m_assert module~sll_m_spline_1d->module~sll_m_working_precision module~sll_m_spline_1d->module~sll_m_bsplines_base

Initialize Module

for minimizer, accept step if dW<dW_allowed*W_MHD(iter=0) default +10e-10 choice phi=Phi_edge * s

Type Bound

t_functional_mhd3d

Arguments

Type IntentOptional Attributes Name
class(t_functional_mhd3d), intent(inout) :: sf

Calls

proc~~initmhd3d~~CallsGraph proc~initmhd3d t_functional_mhd3d%InitMHD3D antiderivative antiderivative proc~initmhd3d->antiderivative boundaryfromfile_new boundaryfromfile_new proc~initmhd3d->boundaryfromfile_new getint getint proc~initmhd3d->getint getintarray getintarray proc~initmhd3d->getintarray getlogical getlogical proc~initmhd3d->getlogical getreal getreal proc~initmhd3d->getreal getstr getstr proc~initmhd3d->getstr initvmec initvmec proc~initmhd3d->initvmec interface~enter_subregion enter_subregion proc~initmhd3d->interface~enter_subregion interface~exit_subregion exit_subregion proc~initmhd3d->interface~exit_subregion interface~par_bcast par_Bcast proc~initmhd3d->interface~par_bcast proc~base_new Base_new proc~initmhd3d->proc~base_new proc~bff_convert_to_modes t_boundaryFromFile%bff_convert_to_modes proc~initmhd3d->proc~bff_convert_to_modes proc~get_imode get_iMode proc~initmhd3d->proc~get_imode proc~getrealallocarray GETREALALLOCARRAY proc~initmhd3d->proc~getrealallocarray proc~hmap_new hmap_new proc~initmhd3d->proc~hmap_new proc~hmap_new_auxvar hmap_new_auxvar proc~initmhd3d->proc~hmap_new_auxvar proc~initializemhd3d_evalfunc InitializeMHD3D_evalFunc proc~initmhd3d->proc~initializemhd3d_evalfunc proc~initprofile InitProfile proc~initmhd3d->proc~initprofile proc~new_minimizer new_minimizer proc~initmhd3d->proc~new_minimizer proc~par_barrier par_Barrier proc~initmhd3d->proc~par_barrier proc~sgrid_init t_sGrid%sGrid_init proc~initmhd3d->proc~sgrid_init interface~enter_subregion->interface~enter_subregion interface~exit_subregion->interface~exit_subregion proc~par_bcast_array1d par_Bcast_array1D interface~par_bcast->proc~par_bcast_array1d proc~par_bcast_array1d_int par_Bcast_array1D_int interface~par_bcast->proc~par_bcast_array1d_int proc~par_bcast_array1d_str par_Bcast_array1D_str interface~par_bcast->proc~par_bcast_array1d_str proc~par_bcast_array2d par_Bcast_array2D interface~par_bcast->proc~par_bcast_array2d proc~par_bcast_scalar par_Bcast_scalar interface~par_bcast->proc~par_bcast_scalar proc~par_bcast_scalar_int par_Bcast_scalar_int interface~par_bcast->proc~par_bcast_scalar_int proc~par_bcast_scalar_str par_Bcast_scalar_str interface~par_bcast->proc~par_bcast_scalar_str proc~base_test Base_test proc~base_new->proc~base_test proc~sbase_new sBase_new proc~base_new->proc~sbase_new proc~fbase_change_base t_fBase%fBase_change_base proc~bff_convert_to_modes->proc~fbase_change_base proc~fbase_evaldof_xn t_fBase%fBase_evalDOF_xn proc~bff_convert_to_modes->proc~fbase_evaldof_xn proc~fbase_initdof t_fBase%fBase_initDOF proc~bff_convert_to_modes->proc~fbase_initdof proc~get_imode->getreal proc~remove_blanks remove_blanks proc~get_imode->proc~remove_blanks interface~findstr FindStr proc~getrealallocarray->interface~findstr proc~converttoproposalstr ConvertToProposalStr proc~getrealallocarray->proc~converttoproposalstr proc~count_sep count_sep proc~getrealallocarray->proc~count_sep proc~hmap_new->interface~enter_subregion proc~hmap_new->interface~exit_subregion proc~initializemhd3d_evalfunc->interface~enter_subregion proc~initializemhd3d_evalfunc->interface~exit_subregion init init proc~initializemhd3d_evalfunc->init proc~initprofile->getreal proc~initprofile->getstr proc~initprofile->interface~enter_subregion proc~initprofile->interface~exit_subregion proc~initprofile->proc~getrealallocarray getrealarray getrealarray proc~initprofile->getrealarray interface~interpolate_cubic_spline interpolate_cubic_spline proc~initprofile->interface~interpolate_cubic_spline proc~sol_var_mhd3d_copy t_sol_var_MHD3D%sol_var_MHD3D_copy proc~new_minimizer->proc~sol_var_mhd3d_copy proc~sol_var_mhd3d_init t_sol_var_MHD3D%sol_var_MHD3D_init proc~new_minimizer->proc~sol_var_mhd3d_init tau tau proc~new_minimizer->tau velocity velocity proc~new_minimizer->velocity proc~sgrid_test sGrid_test proc~sgrid_init->proc~sgrid_test interface~findstr->interface~findstr interface~interpolate_cubic_spline->interface~interpolate_cubic_spline proc~base_evaldof t_base%base_evalDOF proc~base_test->proc~base_evaldof proc~sbase_initdof t_sBase%sBase_initDOF proc~base_test->proc~sbase_initdof proc~converttoproposalstr->proc~remove_blanks proc~fbase_compare t_fBase%fBase_compare proc~fbase_change_base->proc~fbase_compare dgemv dgemv proc~fbase_evaldof_xn->dgemv proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_evaldof_xn->proc~fbase_eval_xn 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~sbase_init t_sBase%sBase_init proc~sbase_new->proc~sbase_init proc~sgrid_test->proc~sgrid_init proc~sgrid_compare t_sGrid%sGrid_compare proc~sgrid_test->proc~sgrid_compare proc~sgrid_find_elem t_sGrid%sGrid_find_elem proc~sgrid_test->proc~sgrid_find_elem proc~sol_var_mhd3d_copy->proc~sol_var_mhd3d_init none~set_to t_sol_var_MHD3D%set_to proc~sol_var_mhd3d_init->none~set_to proc~sol_var_mhd3d_test sol_var_MHD3D_test proc~sol_var_mhd3d_init->proc~sol_var_mhd3d_test proc~sol_var_mhd3d_set_to_scalar t_sol_var_MHD3D%sol_var_MHD3D_set_to_scalar none~set_to->proc~sol_var_mhd3d_set_to_scalar proc~sol_var_mhd3d_set_to_solvar t_sol_var_MHD3D%sol_var_MHD3D_set_to_solvar none~set_to->proc~sol_var_mhd3d_set_to_solvar proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens proc~base_evaldof->proc~fbase_evaldof_ip_tens dgemm dgemm proc~fbase_projectiptodof_tens->dgemm proc~fbase_projectxntodof->dgemv proc~fbase_projectxntodof->proc~fbase_eval_xn proc~sbase_init->init add_element add_element proc~sbase_init->add_element barycentricweights barycentricweights proc~sbase_init->barycentricweights dmatip dmatip proc~sbase_init->dmatip eval_basis eval_basis proc~sbase_init->eval_basis eval_basis_and_n_derivs eval_basis_and_n_derivs proc~sbase_init->eval_basis_and_n_derivs eval_deriv eval_deriv proc~sbase_init->eval_deriv factorize factorize proc~sbase_init->factorize get_interp_points get_interp_points proc~sbase_init->get_interp_points initializevandermonde initializevandermonde proc~sbase_init->initializevandermonde legendregaussnodesandweights legendregaussnodesandweights proc~sbase_init->legendregaussnodesandweights mthpolynomialderivativematrix mthpolynomialderivativematrix proc~sbase_init->mthpolynomialderivativematrix proc~getlu getLU proc~sbase_init->proc~getlu proc~inv INV proc~sbase_init->proc~inv proc~sbase_alloc sBase_alloc proc~sbase_init->proc~sbase_alloc proc~sbase_test sBase_test proc~sbase_init->proc~sbase_test proc~sll_s_bsplines_new sll_s_bsplines_new proc~sbase_init->proc~sll_s_bsplines_new proc~sll_s_spline_matrix_new sll_s_spline_matrix_new proc~sbase_init->proc~sll_s_spline_matrix_new proc~solvemat SOLVEMAT proc~sbase_init->proc~solvemat xiip xiip proc~sbase_init->xiip compute_interpolant compute_interpolant proc~sbase_initdof->compute_interpolant proc~sol_var_mhd3d_test->proc~sol_var_mhd3d_copy proc~sol_var_mhd3d_test->none~set_to proc~sol_var_mhd3d_axby t_sol_var_MHD3D%sol_var_MHD3D_AXBY proc~sol_var_mhd3d_test->proc~sol_var_mhd3d_axby proc~sol_var_mhd3d_norm_2 t_sol_var_MHD3D%sol_var_MHD3D_norm_2 proc~sol_var_mhd3d_test->proc~sol_var_mhd3d_norm_2 proc~fbase_evaldof_ip_tens->proc~fbase_evaldof_xn proc~fbase_evaldof_ip_tens->dgemm dgetrf dgetrf proc~getlu->dgetrf proc~inv->dgetrf dgetri dgetri proc~inv->dgetri proc~sbase_alloc->dmatip proc~sbase_alloc->xiip wbaryip wbaryip proc~sbase_alloc->wbaryip proc~sbase_test->proc~sbase_new proc~sbase_test->proc~sbase_initdof proc~sbase_applybctodof_lgm t_sBase%sBase_applyBCtoDOF_LGM proc~sbase_test->proc~sbase_applybctodof_lgm proc~sbase_change_base t_sBase%sBase_change_base proc~sbase_test->proc~sbase_change_base proc~sbase_compare t_sBase%sBase_compare proc~sbase_test->proc~sbase_compare proc~sbase_eval t_sBase%sBase_eval proc~sbase_test->proc~sbase_eval proc~sbase_evaldof_base t_sBase%sBase_evalDOF_base proc~sbase_test->proc~sbase_evaldof_base proc~sbase_evaldof_gp t_sBase%sBase_evalDOF_GP proc~sbase_test->proc~sbase_evaldof_gp proc~sbase_evaldof_s t_sBase%sBase_evalDOF_s proc~sbase_test->proc~sbase_evaldof_s proc~sll_s_bsplines_new->init proc~sll_s_spline_matrix_new->init proc~sll_s_error_handler sll_s_error_handler proc~sll_s_spline_matrix_new->proc~sll_s_error_handler proc~solvemat->dgetrf dgetrs dgetrs proc~solvemat->dgetrs proc~solve SOLVE proc~sbase_applybctodof_lgm->proc~solve proc~sbase_change_base->proc~sbase_initdof proc~sbase_change_base->proc~sbase_compare eval eval proc~sbase_change_base->eval evalDOF_base evalDOF_base proc~sbase_change_base->evalDOF_base proc~sbase_compare->proc~sgrid_compare proc~sbase_eval->proc~sgrid_find_elem proc~sbase_eval->eval_basis proc~sbase_eval->eval_basis_and_n_derivs lagrangeinterpolationpolys lagrangeinterpolationpolys proc~sbase_eval->lagrangeinterpolationpolys proc~sbase_evaldof_s->proc~sbase_eval proc~sbase_evaldof_s->proc~sbase_evaldof_base interface~c_abort~2 c_abort proc~sll_s_error_handler->interface~c_abort~2 proc~errout errout proc~sll_s_error_handler->proc~errout proc~solve->dgetrf proc~solve->dgetrs

Called by

proc~~initmhd3d~~CalledByGraph proc~initmhd3d t_functional_mhd3d%InitMHD3D proc~init Init proc~init->proc~initmhd3d proc~rungvec rungvec proc~rungvec->proc~initmhd3d program~gvec_post GVEC_POST program~gvec_post->proc~initmhd3d proc~start_rungvec start_rungvec proc~start_rungvec->proc~rungvec program~gvec GVEC program~gvec->proc~rungvec

Source Code

SUBROUTINE InitMHD3D(sf)
  ! MODULES
  USE MODgvec_MHD3D_Vars
  USE MODgvec_Globals        , ONLY: TWOPI
  USE MODgvec_sgrid          , ONLY: t_sgrid
  USE MODgvec_base           , ONLY: base_new
  USE MODgvec_boundaryFromFile, ONLY: boundaryFromFile_new
  USE MODgvec_hmap           , ONLY: hmap_new,hmap_new_auxvar
  USE MODgvec_VMEC           , ONLY: InitVMEC
  USE MODgvec_VMEC_vars      , ONLY: vmec_iota_profile,vmec_pres_profile
  USE MODgvec_VMEC_Readin    , ONLY: nfp,nFluxVMEC,Phi,xm,xn,lasym,mpol,ntor !<<< only exists on MPIroot!
  USE MODgvec_MHD3D_EvalFunc , ONLY: InitializeMHD3D_EvalFunc
  USE MODgvec_ReadInTools    , ONLY: GETSTR,GETLOGICAL,GETINT,GETINTARRAY,GETREAL,GETREALALLOCARRAY, GETREALARRAY
  USE MODgvec_MPI            , ONLY: par_BCast,par_barrier
  USE MODgvec_rProfile_poly  , ONLY: t_rProfile_poly
  USE MODgvec_MHD3D_minimize , ONLY: new_minimizer
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER          :: iMode,nElems,n_sgrid_rho
  INTEGER          :: grid_type
  REAL(wp),ALLOCATABLE :: sgrid_rho(:)
  INTEGER          :: X1X2_deg,X1X2_cont
  INTEGER          :: X1_mn_max(2),X2_mn_max(2)
  INTEGER          :: LA_deg,LA_cont,LA_mn_max(2)
  CHARACTER(LEN=8) :: X1_sin_cos
  CHARACTER(LEN=8) :: X2_sin_cos
  CHARACTER(LEN=8) :: LA_sin_cos
  INTEGER          :: degGP,mn_nyq(2),mn_nyq_min(2),fac_nyq
  INTEGER          :: nfp_loc
  INTEGER          :: X1X2_BCtype_axis(0:4),LA_BCtype_axis(0:4)
  INTEGER          :: proposal_mn_max(1:2)=(/2,0/) !!default proposals, changed for VMEC input to automatically match input!
  CHARACTER(LEN=8) :: proposal_X1_sin_cos="_cos_"  !!default proposals, changed for VMEC input to automatically match input!
  CHARACTER(LEN=8) :: proposal_X2_sin_cos="_sin_"  !!default proposals, changed for VMEC input to automatically match input!
  CHARACTER(LEN=8) :: proposal_LA_sin_cos="_sin_"  !!default proposals, changed for VMEC input to automatically match input!
  REAL(wp)         :: scale_minor_radius
  CHARACTER(LEN=255) ::boundary_filename
  CHARACTER(LEN=8) :: boundary_perturb_type_str !! readin variable for boundary_perturb_type: legacy, cosm
  INTEGER          :: varsize_for_dofs(6)
!===================================================================================================================================
  CALL enter_subregion("init-MHD3D")
  CALL par_Barrier(beforeScreenOut='INIT MHD3D ...')

  which_init = GETINT("whichInitEquilibrium",0)
  IF(which_init.EQ.1) CALL InitVMEC()

  !-----------MINIMIZER
  MinimizerType= GETINT("MinimizerType",Proposal=10)
  PrecondType  = GETINT("PrecondType",Proposal=1)

  dW_allowed=GETREAL("dW_allowed",Proposal=1.0e-10_wp) !! for minimizer, accept step if dW<dW_allowed*W_MHD(iter=0) default +10e-10
  maxIter   = GETINT("maxIter",Proposal=5000)
  outputIter= GETINT("outputIter",Proposal=maxIter)
  logIter   = GETINT("logIter",Proposal=maxIter)
  nlogScreen= GETINT("nLogScreen",Proposal=1)
  minimize_tol  =GETREAL("minimize_tol",Proposal=1.0e-12_wp)
  start_dt  =GETREAL("start_dt",Proposal=0.1_wp)
  doCheckDistance=GETLOGICAL("doCheckDistance",Proposal=.FALSE.)
  doCheckAxis=GETLOGICAL("doCheckAxis",Proposal=.TRUE.)
  !-----------

  fac_nyq = GETINT( "fac_nyq",Proposal=4)

  !constants
  mu_0    = 2.0e-07_wp*TWOPI
  gamm    = 0.0_wp  !fixed!
  sgammM1=1.0_wp/(gamm-1.0_wp)

  init_LA= GETLOGICAL("init_LA",Proposal=.TRUE.)

  SELECT CASE(which_init)
  CASE(0)
    init_fromBConly= .TRUE.
    init_BC        = 2

    !hmap
    which_hmap=GETINT("which_hmap",Proposal=1)
    CALL hmap_new(hmap,which_hmap)
    IF(hmap%nfp.NE.-1)THEN
      nfp_loc = hmap%nfp
    ELSE
      nfp_loc = GETINT("nfp")
    END IF

    Phi_edge   = GETREAL("PHIEDGE",Proposal=1.0_wp)
    Phi_edge   = Phi_edge/TWOPI !normalization like in VMEC!!!
    IF(MPIroot)THEN
      CALL InitProfile(sf,"iota",iota_profile)
      CALL InitProfile(sf,"pres",pres_profile)
    END IF

  CASE(1) !VMEC init
    init_fromBConly= GETLOGICAL("init_fromBConly",Proposal=.FALSE.)
    IF(init_fromBConly)THEN
      !=-1, keep vmec axis and boundary, =0: keep vmec boundary, overwrite axis, =1: keep vmec axis, overwrite boundary, =2: overwrite axis and boundary
      init_BC= GETINT("reinit_BC",Proposal=-1)
    ELSE
      init_BC=-1
    END IF

    IF(MPIroot)THEN
      init_with_profile_iota     = GETLOGICAL("init_with_profile_iota", Proposal=.FALSE.)
      IF(init_with_profile_iota)THEN
        CALL InitProfile(sf,"iota",iota_profile)
      ELSE
        iota_profile=vmec_iota_profile
      END IF ! iota from parameterfile
      init_with_profile_pressure = GETLOGICAL("init_with_profile_pressure", Proposal=.FALSE.)
      IF(init_with_profile_pressure)THEN
        CALL InitProfile(sf,"pres",pres_profile)
      ELSE
        pres_profile=vmec_pres_profile
      END IF ! pressure from parameterfile

      proposal_mn_max(:)=(/mpol-1,ntor/)
      IF(lasym)THEN !asymmetric
        proposal_X1_sin_cos="_sincos_"
        proposal_X2_sin_cos="_sincos_"
        proposal_LA_sin_cos="_sincos_"
      END IF
      nfp_loc = nfp
      Phi_edge = Phi(nFluxVMEC)
    END IF !MPIroot

    which_hmap=1 !hmap_RZ
    CALL hmap_new(hmap,which_hmap)
  END SELECT !which_init

  IF(MPIroot)THEN
    Phi_profile=t_rProfile_poly((/0.0_wp,Phi_edge/)) !! choice phi=Phi_edge * s
    !iota= (dChi/ds) / (dPhi/ds) = dchi_ds / Phi_edge  => chi = Phi_edge * int(iota ds)
    chi_profile=iota_profile%antiderivative()
    chi_profile%coefs=chi_profile%coefs*Phi_edge
  END IF !MPIroot

  getBoundaryFromFile=GETINT("getBoundaryFromFile",Proposal=-1)  ! =-1: OFF, get X1b and X2b from parameterfile. 1: get boundary from specific netcdf file
  SELECT CASE(getBoundaryFromFile)
  CASE(-1)
    !do nothing
  CASE(1)
    boundary_filename=GETSTR("boundary_filename")
    scale_minor_radius=GETREAL("scale_minor_radius",Proposal=1.0_wp)
    IF(MPIroot)THEN
      CALL boundaryFromFile_new(BFF,boundary_filename)
      IF(nfp_loc.NE.BFF%nfp) WRITE(UNIT_stdOut,'(6X,A,I4)')'INFO: changed to boundary file NFP= ',BFF%nfp
      nfp_loc=BFF%nfp
      proposal_mn_max(:)=(/BFF%m_max,BFF%n_max/)
      IF(BFF%lasym.EQ.1)THEN !asymmetric
        proposal_X1_sin_cos="_sincos_"
        proposal_X2_sin_cos="_sincos_"
        proposal_LA_sin_cos="_sincos_"
      END IF
    END IF
  END SELECT
  CALL par_BCast(proposal_mn_max,0)
  CALL par_BCast(proposal_X1_sin_cos,0)
  CALL par_BCast(proposal_X2_sin_cos,0)
  CALL par_BCast(proposal_LA_sin_cos,0)
  CALL par_BCast(nfp_loc,0)

  CALL enter_subregion("discretization")
  X1X2_deg     = GETINT(     "X1X2_deg")
  X1X2_cont    = GETINT(     "X1X2_continuity",Proposal=(X1X2_deg-1) )
  X1_mn_max    = GETINTARRAY("X1_mn_max"   ,2 ,Proposal=proposal_mn_max)
  X2_mn_max    = GETINTARRAY("X2_mn_max"   ,2 ,Proposal=proposal_mn_max)
  X1_sin_cos   = GETSTR(     "X1_sin_cos"     ,Proposal=proposal_X1_sin_cos)  !_sin_,_cos_,_sin_cos_
  X2_sin_cos   = GETSTR(     "X2_sin_cos"     ,Proposal=proposal_X2_sin_cos)


  LA_deg     = GETINT(     "LA_deg")
  LA_cont    = GETINT(     "LA_continuity",Proposal=(LA_deg-1))
  LA_mn_max  = GETINTARRAY("LA_mn_max", 2 ,Proposal=proposal_mn_max)
  LA_sin_cos = GETSTR(     "LA_sin_cos"   ,Proposal=proposal_LA_sin_cos)

  IF(fac_nyq.EQ.-1)THEN
    fac_nyq=2
    mn_nyq_min(1)=1+fac_nyq*MAXVAL((/X1_mn_max(1),X2_mn_max(1),LA_mn_max(1)/))
    mn_nyq_min(2)=1+fac_nyq*(MAXVAL((/X1_mn_max(2),X2_mn_max(2),LA_mn_max(2)/))+hmap%n_max)
    mn_nyq  = GETINTARRAY("mn_nyq",2)
    IF(mn_nyq(1).LT.mn_nyq_min(1))THEN
       SWRITE(UNIT_stdOut,'(A,I6)')'WARNING: mn_nyq(1) too small, should be >= ',mn_nyq_min(1)
    END IF
    IF(mn_nyq(2).LT.mn_nyq_min(2))THEN
       SWRITE(UNIT_stdOut,'(A,I6)') 'WARNING: mn_nyq(2) too small, should be >= ',mn_nyq_min(2)
    END IF
  ELSE
    mn_nyq(1)=1+fac_nyq*MAXVAL((/X1_mn_max(1),X2_mn_max(1),LA_mn_max(1)/))
    mn_nyq(2)=1+fac_nyq*(MAXVAL((/X1_mn_max(2),X2_mn_max(2),LA_mn_max(2)/))+hmap%n_max)
  END IF

  SWRITE(UNIT_stdOut,*)
  SWRITE(UNIT_stdOut,'(2(A,I4),A,I6," , ",I6,A)')'    fac_nyq = ', fac_nyq,' hmap%n_max = ',hmap%n_max,'  ==> interpolation points mn_nyq=( ',mn_nyq(:),' )'
  SWRITE(UNIT_stdOut,*)

  grid_type= GETINT("sgrid_grid_type",Proposal=0)
  degGP    = GETINT("degGP",Proposal=MAX(X1X2_deg,LA_deg)+2)

  !INITIALIZE GRID
  IF (grid_type.NE.GRID_TYPE_CUSTOM) THEN
    nElems   = GETINT("sgrid_nElems")
    CALL sgrid%init(nElems,grid_type)
  ELSE
    CALL GETREALALLOCARRAY("sgrid_rho", sgrid_rho, n_sgrid_rho)
    nElems = n_sgrid_rho - 1
    CALL sgrid%init(nElems,grid_type,sgrid_rho)
    SDEALLOCATE(sgrid_rho)
  END IF

  !INITIALIZE BASE        !sbase parameter                 !fbase parameter               ...exclude_mn_zero
  CALL base_new(X1_base  , X1X2_deg,X1X2_cont,sgrid,degGP , X1_mn_max,mn_nyq,nfp_loc,X1_sin_cos,.FALSE.)
  CALL base_new(X2_base  , X1X2_deg,X1X2_cont,sgrid,degGP , X2_mn_max,mn_nyq,nfp_loc,X2_sin_cos,.FALSE.)
  CALL base_new(LA_base  ,   LA_deg,  LA_cont,sgrid,degGP , LA_mn_max,mn_nyq,nfp_loc,LA_sin_cos,.TRUE. )

  CALL hmap_new_auxvar(hmap,X1_base%f%zeta_IP,hmap_auxvar,.FALSE.) !no second derivative needed!

  IF((which_init.EQ.1).AND.MPIroot) THEN !VMEC
    IF(lasym)THEN
      IF((X1_base%f%sin_cos.NE._SINCOS_).OR. &
         (X2_base%f%sin_cos.NE._SINCOS_).OR. &
         (LA_base%f%sin_cos.NE._SINCOS_) ) THEN
        WRITE(UNIT_stdOut,'(A)')'!!!!!!!! WARNING: !!!!!!!!!!!!!!!'
        WRITE(UNIT_stdOut,'(A)')'!!!!!!!!   ---->  VMEC was run asymmetric, you should use _sincos_ basis for all variables'
        WRITE(UNIT_stdOut,'(A)')'!!!!!!!! WARNING: !!!!!!!!!!!!!!!'
        !CALL abort(__STAMP__,&
        !    '!!!!  VMEC was run asymmetric, you should use _sincos_ basis for all variables')
      END IF
    END IF
    IF((MAXVAL(INT(xm(:))).GT.MINVAL((/X1_mn_max(1),X2_mn_max(1),LA_mn_max(1)/))).OR. &
       (MAXVAL(ABS(INT(xn(:))/nfp_loc)).GT.MINVAL((/X1_mn_max(2),X2_mn_max(2),LA_mn_max(2)/))))THEN
      WRITE(UNIT_stdOut,'(A)')    '!!!!!!!! WARNING: !!!!!!!!!!!!!!!'
      WRITE(UNIT_stdOut,'(A,2I6)')'!!!!!!!!   ---->  you use a lower mode number than the VMEC  run  ', &
                                    MAXVAL(INT(xm(:))),MAXVAL(ABS(INT(xn(:))/nfp_loc))
      WRITE(UNIT_stdOut,'(A)')    '!!!!!!!! WARNING: !!!!!!!!!!!!!!!'
      !  CALL abort(__STAMP__,&
      !'!!!!!  you use a lower mode number than the VMEC  run  (m,n)_max')
    END IF
  END IF

  nDOF_X1 = X1_base%s%nBase* X1_base%f%modes
  nDOF_X2 = X2_base%s%nBase* X2_base%f%modes
  nDOF_LA = LA_base%s%nBase* LA_base%f%modes



  !X1X2_BCtype_axis(MN_ZERO    )= GETINT("X1X2_BCtype_axis_mn_zero"    ,Proposal=0 ) !AUTOMATIC,m-dependent
  !X1X2_BCtype_axis(M_ZERO     )= GETINT("X1X2_BCtype_axis_m_zero"     ,Proposal=0 ) !AUTOMATIC,m-dependent
  !X1X2_BCtype_axis(M_ODD_FIRST)= GETINT("X1X2_BCtype_axis_m_odd_first",Proposal=0 ) !AUTOMATIC,m-dependent
  !X1X2_BCtype_axis(M_ODD      )= GETINT("X1X2_BCtype_axis_m_odd"      ,Proposal=0 ) !AUTOMATIC,m-dependent
  !X1X2_BCtype_axis(M_EVEN     )= GETINT("X1X2_BCtype_axis_m_even"     ,Proposal=0 ) !AUTOMATIC,m-dependent
  X1X2_BCtype_axis= 0 !fix to AUTOMATIC, m-dependent

  !boundary conditions (used in force, in init slightly changed)
  ASSOCIATE(modes        =>X1_base%f%modes, zero_odd_even=>X1_base%f%zero_odd_even)
  ALLOCATE(X1_BC_type(1:2,modes))
  X1_BC_type(BC_EDGE,:)=BC_TYPE_DIRICHLET
  DO imode=1,modes
    X1_BC_type(BC_AXIS,iMode)=X1X2_BCtype_axis(zero_odd_even(iMode))
    IF(X1_BC_type(BC_AXIS,iMode).EQ.0)THEN !AUTOMATIC, m-dependent BC, for m>deg, switch off all DOF up to deg+1
      X1_BC_type(BC_AXIS,iMode)=-1*MIN(X1_base%s%deg+1,X1_base%f%Xmn(1,iMode))
      IF(MPIroot.AND.(nElems.EQ.1).AND.(X1_base%f%Xmn(1,iMode).GT.X1_base%s%deg)) THEN
        IF(X1_base%f%Xmn(2,iMode).EQ.0) & !warning for all n-modes written once!
           WRITE(UNIT_stdOut,'(4X,A,I4,A)')'WARNING, 1-element spline with BC for m>deg, will ZERO edge coeff. X1_b(m=',&
                                          X1_base%f%Xmn(1,iMode),',n=-n_max,n_max)! (use 2elems instead)'
      END IF
    END IF
  END DO
  END ASSOCIATE !X1

  ASSOCIATE(modes        =>X2_base%f%modes, zero_odd_even=>X2_base%f%zero_odd_even)
  ALLOCATE(X2_BC_type(1:2,modes))
  X2_BC_type(BC_EDGE,:)=BC_TYPE_DIRICHLET
  DO imode=1,modes
    X2_BC_type(BC_AXIS,iMode)=X1X2_BCtype_axis(zero_odd_even(iMode))
    IF(X2_BC_type(BC_AXIS,iMode).EQ.0)THEN ! AUTOMATIC, m-dependent BC, for m>deg, switch off all DOF up to deg+1
      X2_BC_type(BC_AXIS,iMode)=-1*MIN(X2_base%s%deg+1,X2_base%f%Xmn(1,iMode))
      IF(MPIroot.AND.(nElems.EQ.1).AND.(X2_base%f%Xmn(1,iMode).GT.X2_base%s%deg)) THEN
        IF(X2_base%f%Xmn(2,iMode).EQ.0) & !warning for all n-modes written once!
          WRITE(UNIT_stdOut,'(4X,A,I4,A)')'WARNING, 1-element spline with BC for m>deg, will ZERO edge coeff. X2_b(m=',&
                                          X2_base%f%Xmn(1,iMode),',n=-n_max,n_max)! (use 2elems instead)'
      END IF
    END IF
  END DO
  END ASSOCIATE !X2

  !LA_BCtype_axis(MN_ZERO    )= GETINT("LA_BCtype_axis_mn_zero"    ,Proposal=BC_TYPE_SYMMZERO )
  !LA_BCtype_axis(M_ZERO     )= GETINT("LA_BCtype_axis_m_zero"     ,Proposal=BC_TYPE_SYMM     )
  !LA_BCtype_axis(M_ODD_FIRST)= GETINT("LA_BCtype_axis_m_odd_first",Proposal=0 ) !AUTOMATIC,m-dependent
  !LA_BCtype_axis(M_ODD      )= GETINT("LA_BCtype_axis_m_odd"      ,Proposal=0 ) !AUTOMATIC,m-dependent
  !LA_BCtype_axis(M_EVEN     )= GETINT("LA_BCtype_axis_m_even"     ,Proposal=0 ) !AUTOMATIC,m-dependent
  LA_BCtype_axis= 0 !fix to AUTOMATIC, m-dependent

  ASSOCIATE(modes        =>LA_base%f%modes, zero_odd_even=>LA_base%f%zero_odd_even)
  ALLOCATE(LA_BC_type(1:2,modes))
  LA_BC_type(BC_EDGE,:)=BC_TYPE_OPEN
  DO imode=1,modes
    LA_BC_type(BC_AXIS,iMode)=LA_BCtype_axis(zero_odd_even(iMode))
    IF(LA_BC_type(BC_AXIS,iMode).EQ.0)THEN ! AUTOMATIC, m-dependent BC, for m>deg, switch off all DOF up to deg+1
      LA_BC_type(BC_AXIS,iMode)=-1*MIN(LA_base%s%deg+1,LA_base%f%Xmn(1,iMode))
    END IF
  END DO
  END ASSOCIATE !LA

  CALL exit_subregion("discretization")
  CALL enter_subregion("boundary")
  !INITIALIZATION PARAMETERS (ONLY NECESSARY ON MPIroot)
  IF(MPIroot)THEN
    init_average_axis= GETLOGICAL("init_average_axis",Proposal=.FALSE.)
    IF(init_average_axis)THEN
      average_axis_move(1) = GETREAL("average_axis_move_X1",Proposal=0.0_wp)
      average_axis_move(2) = GETREAL("average_axis_move_X2",Proposal=0.0_wp)
    END IF
    ALLOCATE(X1_b(1:X1_base%f%modes) )
    ALLOCATE(X2_b(1:X2_base%f%modes) )
    ALLOCATE(LA_b(1:LA_base%f%modes) )
    ALLOCATE(X1_a(1:X1_base%f%modes) )
    ALLOCATE(X2_a(1:X2_base%f%modes) )
    X1_b=0.0_wp
    X2_b=0.0_wp
    LA_b=0.0_wp
    X1_a=0.0_wp
    X2_a=0.0_wp

    IF((init_BC.EQ.0).OR.(init_BC.EQ.2))THEN !READ axis values from input file
      WRITE(UNIT_stdOut,'(4X,A)')'... read axis data for X1:'
      ASSOCIATE(modes=>X1_base%f%modes,sin_range=>X1_base%f%sin_range,cos_range=>X1_base%f%cos_range)
      DO iMode=sin_range(1)+1,sin_range(2)
        X1_a(iMode)=get_iMode('X1_a_sin',X1_base%f%Xmn(:,iMode),X1_base%f%nfp)
      END DO !iMode
      DO iMode=cos_range(1)+1,cos_range(2)
        X1_a(iMode)=get_iMode('X1_a_cos',X1_base%f%Xmn(:,iMode),X1_base%f%nfp)
      END DO !iMode
      END ASSOCIATE
      WRITE(UNIT_stdOut,'(4X,A)')'... read axis data for X2:'
      ASSOCIATE(modes=>X2_base%f%modes,sin_range=>X2_base%f%sin_range,cos_range=>X2_base%f%cos_range)
      DO iMode=sin_range(1)+1,sin_range(2)
        X2_a(iMode)=get_iMode('X2_a_sin',X2_base%f%Xmn(:,iMode),X2_base%f%nfp)
      END DO !iMode
      DO iMode=cos_range(1)+1,cos_range(2)
        X2_a(iMode)=get_iMode('X2_a_cos',X2_base%f%Xmn(:,iMode),X2_base%f%nfp)
      END DO !iMode
      END ASSOCIATE
    END IF
    IF(getBoundaryFromFile.EQ.-1)THEN
      IF(((init_BC.EQ.1).OR.(init_BC.EQ.2)))THEN !READ edge values from input file
        WRITE(UNIT_stdOut,'(4X,A)')'... read edge boundary data for X1:'
        ASSOCIATE(modes=>X1_base%f%modes,sin_range=>X1_base%f%sin_range,cos_range=>X1_base%f%cos_range)
        DO iMode=sin_range(1)+1,sin_range(2)
          X1_b(iMode)=get_iMode('X1_b_sin',X1_base%f%Xmn(:,iMode),X1_base%f%nfp)
        END DO !iMode
        DO iMode=cos_range(1)+1,cos_range(2)
          X1_b(iMode)=get_iMode('X1_b_cos',X1_base%f%Xmn(:,iMode),X1_base%f%nfp)
        END DO !iMode
        END ASSOCIATE
        WRITE(UNIT_stdOut,'(4X,A)')'... read edge boundary data for X2:'
        ASSOCIATE(modes=>X2_base%f%modes,sin_range=>X2_base%f%sin_range,cos_range=>X2_base%f%cos_range)
        DO iMode=sin_range(1)+1,sin_range(2)
          X2_b(iMode)=get_iMode('X2_b_sin',X2_base%f%Xmn(:,iMode),X2_base%f%nfp)
        END DO !iMode
        DO iMode=cos_range(1)+1,cos_range(2)
          X2_b(iMode)=get_iMode('X2_b_cos',X2_base%f%Xmn(:,iMode),X2_base%f%nfp)
        END DO !iMode
        END ASSOCIATE
      END IF !init_BC
    ELSE !getBoundaryFromFile
      CALL BFF%convert_to_modes(X1_base%f,X2_base%f,X1_b,X2_b,scale_minor_radius)
    END IF
  END IF !MPIroot

  boundary_perturb = GETLOGICAL('boundary_perturb', Proposal=.FALSE.)
  boundary_perturb_type_str  = GETSTR("boundary_perturb_type", proposal="legacy")
  IF (boundary_perturb_type_str .EQ. "legacy") THEN
    boundary_perturb_type = BLEND_LEGACY
  ELSE IF (boundary_perturb_type_str .EQ. "cosm") THEN
    boundary_perturb_type = BLEND_COSM
  ELSE
    CALL abort(__STAMP__,&
    'boundary_perturb_type must be "legacy" or "cosm", found '//TRIM(boundary_perturb_type_str),&
    intInfo=boundary_perturb_type,&
    TypeInfo="InvalidParameterError")
  END IF
  boundary_perturb_depth = GETREAL('boundary_perturb_depth', proposal=0.6_wp)
  IF(boundary_perturb)THEN
    ALLOCATE(X1pert_b(1:X1_base%f%modes) )
    ALLOCATE(X2pert_b(1:X2_base%f%modes) )
    X1pert_b=0.0_wp
    X2pert_b=0.0_wp
    !READ boudnary values from input file
    ASSOCIATE(modes=>X1_base%f%modes,sin_range=>X1_base%f%sin_range,cos_range=>X1_base%f%cos_range)
    SWRITE(UNIT_stdOut,'(4X,A)')'... read data for X1pert:'
    DO iMode=sin_range(1)+1,sin_range(2)
      X1pert_b(iMode)=get_iMode('X1pert_b_sin',X1_base%f%Xmn(:,iMode),X1_base%f%nfp)
    END DO !iMode
    DO iMode=cos_range(1)+1,cos_range(2)
      X1pert_b(iMode)=get_iMode('X1pert_b_cos',X1_base%f%Xmn(:,iMode),X1_base%f%nfp)
    END DO !iMode
    END ASSOCIATE
    ASSOCIATE(modes=>X2_base%f%modes,sin_range=>X2_base%f%sin_range,cos_range=>X2_base%f%cos_range)
    SWRITE(UNIT_stdOut,'(4X,A)')'... read data for X2pert:'
    DO iMode=sin_range(1)+1,sin_range(2)
      X2pert_b(iMode)=get_iMode('X2pert_b_sin',X2_base%f%Xmn(:,iMode),X2_base%f%nfp)
    END DO !iMode
    DO iMode=cos_range(1)+1,cos_range(2)
      X2pert_b(iMode)=get_iMode('X2pert_b_cos',X2_base%f%Xmn(:,iMode),X2_base%f%nfp)
    END DO !iMode
    END ASSOCIATE
  END IF
  CALL exit_subregion("boundary")
  CALL enter_subregion("initialize minimizer")
  varsize_for_dofs = (/X1_base%s%nbase,X2_base%s%nbase,LA_base%s%nBase,&
    X1_base%f%modes,X2_base%f%modes,LA_base%f%modes/)
  CALL new_minimizer(&
    sf=sf%minimizer,varsize_in = varsize_for_dofs, &
    dt_initial=start_dt, MinType_in=MinimizerType, dW_allowed_in=dW_allowed,& ! Minimizer
    DoCheckDistance=DoCheckDistance, DoCheckAxis=DoCheckAxis,& !what to log
    outputIter=outputIter, nlogScreen=nlogScreen, logIter=logIter& !when to log
  )
  CALL exit_subregion("initialize minimizer")

  CALL InitializeMHD3D_EvalFunc()
  CALL exit_subregion("init-MHD3D")
  CALL par_barrier(afterScreenOut='...DONE')
  SWRITE(UNIT_stdOut,fmt_sep)

END SUBROUTINE InitMHD3D