mhd3d.F90 – GVEC

!!matvec with matmul !#define __MATVEC_N(y,Mat,Vec) y=MATMUL(Mat,Vec) !#define __MATVEC_T(y,Mat,Vec) y=MATMUL(Vec,Mat) !#define __PMATVEC_N(fy,y,Mat,Vec) y=fyy+MATMUL(Mat,Vec) !#define __PMATVEC_T(fy,y,Mat,Vec) y=fyy+MATMUL(Vec,Mat) !#define __AMATVEC_N(y,fMat,Mat,Vec) y=fMatMATMUL(Mat,Vec) !#define __AMATVEC_T(y,fMat,Mat,Vec) y=fMatMATMUL(Vec,Mat) !#define __PAMATVEC_N(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Mat,Vec) !#define __PAMATVEC_T(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Vec,Mat)
!!#define __GENERICMATVEC(NT,fy,y,fMat,Mat,Vec) CALL DGEMV(NT,SIZE(Mat,1),SIZE(Mat,2),fMat,Mat,SIZE(Mat,1),Vec,1,fy,y,1)
!!matmat with matmul !#define __MATMAT_NN(Y,A,B) Y=MATMUL(A,B) !#define __MATMAT_TN(Y,A,B) Y=MATMUL(TRANSPOSE(A),B) !#define __MATMAT_NT(Y,A,B) Y=MATMUL(A,TRANSPOSE(B)) !#define __MATMAT_TT(Y,A,B) Y=TRANSPOSE(MATMUL(B,A))
!#define __PMATMAT_NN(fy,Y,A,B) Y=fyY+MATMUL(A,B) !#define __PMATMAT_TN(fy,Y,A,B) Y=fyY+MATMUL(TRANSPOSE(A),B) !#define __PMATMAT_NT(fy,Y,A,B) Y=fyY+MATMUL(A,TRANSPOSE(B)) !#define __PMATMAT_TT(fy,Y,A,B) Y=fyY+TRANSPOSE(MATMUL(B,A))
!#define __AMATMAT_NN(Y,fa,A,B) Y=faMATMUL(A,B) !#define __AMATMAT_TN(Y,fa,A,B) Y=faMATMUL(TRANSPOSE(A),B) !#define __AMATMAT_NT(Y,fa,A,B) Y=faMATMUL(A,TRANSPOSE(B)) !#define __AMATMAT_TT(Y,fa,A,B) Y=faTRANSPOSE(MATMUL(B,A))
!#define __PAMATMAT_NN(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,B) !#define __PAMATMAT_TN(fy,Y,fa,A,B) Y=fyY+faMATMUL(TRANSPOSE(A),B) !#define __PAMATMAT_NT(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,TRANSPOSE(B)) !#define __PAMATMAT_TT(fy,Y,fa,A,B) Y=fyY+faTRANSPOSE(MATMUL(B,A))
! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)
!#define __GENERICMATMAT_NN(fy,Y,fa,A,B) CALL DGEMM('N','N',SIZE(A,1),SIZE(B,2),SIZE(B,1),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TN(fy,Y,fa,A,B) CALL DGEMM('T','N',SIZE(A,2),SIZE(B,2),SIZE(B,1),fa,A,SIZE(B,1),B,SIZE(B,1),fy,Y,SIZE(A,2)) !#define __GENERICMATMAT_NT(fy,Y,fa,A,B) CALL DGEMM('N','T',SIZE(A,1),SIZE(B,1),SIZE(B,2),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TT(fy,Y,fa,A,B) CALL DGEMM('T','T',SIZE(A,2),SIZE(B,1),SIZE(B,2),fa,A,SIZE(B,2),B,SIZE(B,1),fy,Y,SIZE(A,2))
! SIMPLE INTERFACE FOR DGEMM, specifying nrows/ncols of mat A and nrows/ncols of mat B (for any transpose!) ! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)
Files dependent on this one

Help
Source Code

!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.FPP"

!===================================================================================================================================
!>
!!# Module **MHD3D**
!!
!! CONTAINS INITIALIZATION OF MHD 3D Energy functional that will be minimized
!!
!===================================================================================================================================
MODULE MODgvec_MHD3D
! MODULES
  USE MODgvec_Globals, ONLY:wp,abort,UNIT_stdOut,fmt_sep,MPIRoot,enter_subregion,exit_subregion
  USE MODgvec_MHD3D_minimize, ONLY: t_minimizer_mhd3d
  IMPLICIT NONE

  PRIVATE
  PUBLIC t_functional_mhd3d

!-----------------------------------------------------------------------------------------------------------------------------------
! TYPES
!-----------------------------------------------------------------------------------------------------------------------------------

  TYPE :: t_functional_mhd3d
    !-------------------------------------------------------------------------------------------------------------------------------
    LOGICAL :: initialized
    TYPE(t_minimizer_mhd3d), ALLOCATABLE :: minimizer
    !-------------------------------------------------------------------------------------------------------------------------------
    CONTAINS
      PROCEDURE :: init         => InitMHD3D
      PROCEDURE :: initSolution => InitSolutionMHD3D
      PROCEDURE :: minimize     => MinimizeMHD3D
      PROCEDURE :: free         => FinalizeMHD3D
  END TYPE t_functional_mhd3d



!===================================================================================================================================

CONTAINS

!===================================================================================================================================
!> Initialize Module
!!
!===================================================================================================================================
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

SUBROUTINE InitProfile(sf, var,var_profile)
  ! MODULES
  USE MODgvec_ReadInTools    , ONLY: GETSTR,GETLOGICAL,GETINT,GETINTARRAY,GETREAL,GETREALALLOCARRAY, GETREALARRAY
  USE MODgvec_rProfile_bspl  , ONLY: t_rProfile_bspl
  USE MODgvec_rProfile_poly  , ONLY: t_rProfile_poly
  USE MODgvec_cubic_spline   , ONLY: interpolate_cubic_spline
  USE MODgvec_rProfile_base, ONLY: c_rProfile
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
  CHARACTER(LEN=4), INTENT(IN) :: var
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(c_rProfile), ALLOCATABLE   :: var_profile
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER              :: n_profile_knots
  REAL(wp),ALLOCATABLE :: profile_knots(:)
  INTEGER              :: n_profile_coefs    !! number of polynomial/bspline coeffients for profile profile
  REAL(wp),ALLOCATABLE :: profile_coefs(:)   !! polynomial/bspline coefficients of the profile profile
  CHARACTER(LEN=20)    :: profile_type
  REAL(wp),ALLOCATABLE :: profile_rho2(:)
  REAL(wp),ALLOCATABLE :: profile_vals(:)
  INTEGER              :: n_profile_vals
  INTEGER              :: n_profile_rho2
  REAL(wp)             :: profile_BC_vals(1:2)
  CHARACTER(LEN=10)    :: profile_BC_type(1:2)
  REAL(wp)             :: profile_scale
  CHARACTER(LEN=10)    :: possible_BCs(0:2)
  INTEGER              :: BC(1:2),iBC,jBC
  !===================================================================================================================================
  CALL enter_subregion("get-profile-"//TRIM(var))
  possible_BCs(0)="not_a_knot"
  possible_BCs(1)="1st_deriv"
  possible_BCs(2)="2nd_deriv"

  profile_scale=GETREAL(var//"_scale",Proposal=1.0_wp)
  profile_type  = GETSTR(var//"_type") !make it mandatory
  IF (profile_type.EQ."polynomial") THEN
    CALL GETREALALLOCARRAY(var//"_coefs",profile_coefs,n_profile_coefs) !a+b*s+c*s^2...
    profile_coefs=profile_coefs*profile_scale
    var_profile = t_rProfile_poly(profile_coefs)
  ELSE IF (profile_type.EQ."bspline") THEN
    CALL GETREALALLOCARRAY(var//"_coefs",profile_coefs,n_profile_coefs)
    profile_coefs=profile_coefs*profile_scale
    CALL GETREALALLOCARRAY(var//"_knots",profile_knots,n_profile_knots)
    IF(ABS(profile_knots(1)).GT.1.0d-12) CALL abort(__STAMP__,&
        "First knot position must be =0 for bspline of "//TRIM(var)//" profile!",&
        TypeInfo="InvalidParameterError")
    IF(ABS(profile_knots(n_profile_knots)-1.0_wp).GT.1.0d-12) CALL abort(__STAMP__,&
        "Last knot position must be =1 for bspline of "//TRIM(var)//" profile!",&
        TypeInfo="InvalidParameterError")
    var_profile = t_rProfile_bspl(coefs=profile_coefs,knots=profile_knots)
  ELSE IF (profile_type.EQ."interpolation") THEN
    CALL GETREALALLOCARRAY(var//"_vals",profile_vals, n_profile_vals)
    CALL GETREALALLOCARRAY(var//"_rho2",profile_rho2, n_profile_rho2)
    IF(ABS(profile_rho2(1)).GT.1.0d-12) CALL abort(__STAMP__,&
      "First rho2 position must be =0 for interpolation of "//TRIM(var)//" profile!",&
        TypeInfo="InvalidParameterError")
    IF(ABS(profile_rho2(n_profile_rho2)-1.0_wp).GT.1.0d-12) CALL abort(__STAMP__,&
      "Last rho2 position must be =1 for interpolation of "//TRIM(var)//" profile!",&
        TypeInfo="InvalidParameterError")
    IF (n_profile_vals .NE. n_profile_rho2) THEN
      CALL abort(__STAMP__,&
      'Size of '//var//'_rho2 and '//var//'_vals must be equal!',&
        TypeInfo="InvalidParameterError")
    END IF
    profile_BC_type(1) = GETSTR(var//"_BC_type_axis",Proposal="not_a_knot")
    profile_BC_type(2) = GETSTR(var//"_BC_type_edge",Proposal="not_a_knot")
    BC=-1
    DO iBC=1,2
      DO jBC=0,2
        IF (INDEX(TRIM(profile_BC_type(iBC)),TRIM(possible_BCs(jBC)))>0) THEN
          BC(iBC)=jBC
          EXIT
        END IF
      END DO !jBC
    END DO !iBC
    IF(ANY(BC<0)) THEN
      CALL abort(__STAMP__,&
          "BC_type of profile must be 'not_a_knot', '1st_deriv' or '2nd_deriv' ... got '"//TRIM(profile_BC_type(1))//"' on axis and '"//TRIM(profile_BC_type(2))//"' on edge",&
        TypeInfo="InvalidParameterError")
    END IF

    IF(ANY(BC>0)) THEN
     profile_BC_vals = GETREALARRAY(var//"_BC_vals", 2, Proposal=(/0.0_wp, 0.0_wp/))
     CALL interpolate_cubic_spline(profile_rho2,profile_vals, profile_coefs, profile_knots, BC, profile_BC_vals)
    ELSE
     CALL interpolate_cubic_spline(profile_rho2,profile_vals, profile_coefs, profile_knots, BC)
    END IF

    profile_coefs=profile_coefs*profile_scale
    var_profile = t_rProfile_bspl(profile_knots,profile_coefs)
    SDEALLOCATE(profile_vals)
    SDEALLOCATE(profile_rho2)
  ELSE
    CALL abort(__STAMP__,&
         "Specified "//var//"_type unknown. Expecting 'polynomial', 'bspline' or 'interpolation' ... got '"//TRIM(profile_type)//"'",&
        TypeInfo="InvalidParameterError")
  END IF ! profile type

  SDEALLOCATE(profile_knots)
  SDEALLOCATE(profile_coefs)
  CALL exit_subregion("get-profile-"//TRIM(var))
END SUBROUTINE InitProfile

!===================================================================================================================================
!> Initialize Module
!!
!===================================================================================================================================
SUBROUTINE InitSolutionMHD3D(sf)
! MODULES
  USE MODgvec_MHD3D_Vars     , ONLY: which_init,init_LA,boundary_perturb,boundary_perturb_depth,boundary_perturb_type
  USE MODgvec_Restart_vars   , ONLY: doRestart,RestartFile
  USE MODgvec_Restart        , ONLY: RestartFromState
  USE MODgvec_Restart        , ONLY: WriteState
  USE MODgvec_MHD3D_EvalFunc , ONLY: InitProfilesGP,EvalEnergy,EvalForce
  USE MODgvec_Analyze        , ONLY: Analyze
  USE MODgvec_ReadInTools    , ONLY: GETLOGICAL
  USE MODgvec_MPI            , ONLY: par_Bcast,par_barrier
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER              :: JacCheck
!===================================================================================================================================
  CALL par_barrier(beforeScreenOut="    INITIALIZE SOLUTION...",afterScreenOut="                           ...")
  CALL enter_subregion("init-solution")
  ASSOCIATE(vars=>sf%minimizer%vars)
    IF(MPIroot) THEN
      IF(doRestart)THEN
        WRITE(UNIT_stdOut,'(4X,A)')'... restarting from file ... '
        CALL RestartFromState(RestartFile,vars%dofs(0))
        CALL InitSolution(vars%dofs(0),-1) ! (re-)apply BC and init LA (if init_LA is true)
      ELSE
        CALL InitSolution(vars%dofs(0),which_init)
      END IF
      IF(boundary_perturb)THEN
        CALL AddBoundaryPerturbation(vars%dofs(0),boundary_perturb_depth,boundary_perturb_type)
      END IF !boundary_perturb
    END IF !MPIroot
    CALL par_Bcast(vars%dofs(0)%X1,0)
    CALL par_Bcast(vars%dofs(0)%X2,0)
    CALL exit_subregion("init-solution")

    IF(init_LA) THEN
      CALL Init_LA_From_Solution(vars%dofs(0))  !BCast inside
    ELSE
      CALL par_Bcast(vars%dofs(0)%LA,0)
    END IF

    CALL vars%dofs(-1)%set_to(vars%dofs(0))
    CALL vars%dofs(-2)%set_to(vars%dofs(0))
    CALL vars%dofs(-3)%set_to(vars%dofs(0))

    JacCheck=2
    CALL InitProfilesGP() !evaluate profiles once at Gauss Points (on MPIroot + BCast)

    CALL enter_subregion("check-solution")
    vars%dofs(0)%W_MHD3D=EvalEnergy(vars%dofs(0),.TRUE.,JacCheck)
    IF(JacCheck.LE.-1)THEN
      CALL Analyze(0, vars%dofs(0), vars%force(0))
      SELECT CASE(JacCheck)
      CASE(-1)
        CALL abort(__STAMP__,&
          "JACOBIAN WITH SIGN CHANGE FOUND AFTER INITIALIZATION!", &
          TypeInfo="InitializationError")
      CASE(-2)
        CALL abort(__STAMP__,&
          "RELATIVE JACOBIAN TOO SMALL AFTER INITIALIZATION!", &
          TypeInfo="InitializationError")
      CASE(-3)
        CALL abort(__STAMP__,&
          "ALL JACOBIANS NEGATIVE AFTER INITIALIZATION! (LEFT-HANDED COORDINATES)", &
          TypeInfo="InitializationError")
      END SELECT
    END IF
    CALL WriteState(vars%dofs(0),0)
    CALL EvalForce(vars%dofs(0),.FALSE.,JacCheck, vars%force(0))
    SWRITE(UNIT_stdOut,'(8x,A,3E11.4)')'|Force|= ',SQRT(vars%force(0)%norm_2())
    CALL Analyze(0, vars%dofs(0), vars%force(0))
  END ASSOCIATE !vars
  CALL exit_subregion("check-solution")
  CALL par_barrier(afterScreenOut="    ...DONE")
  SWRITE(UNIT_stdOut,fmt_sep)
END SUBROUTINE InitSolutionMHD3D


!===================================================================================================================================
!> automatically build the string to be read from parameterfile, varname + m,n mode number, and then read it from parameterfile
!!
!===================================================================================================================================
FUNCTION get_iMode(varname_in,mn_in,nfp_in)
! MODULES
  USE MODgvec_ReadInTools    , ONLY: GETREAL,remove_blanks
!$ USE omp_lib
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
    CHARACTER(LEN=*),INTENT(IN) :: varname_in
    INTEGER         ,INTENT(IN) :: mn_in(2),nfp_in
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
    REAL(wp)                    :: get_iMode
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
    CHARACTER(LEN=100) :: varstr
!===================================================================================================================================
  SWRITE(varstr,'(A,"(",I4,";",I4,")")')TRIM(varname_in),mn_in(1),mn_in(2)/nfp_in
  varstr=remove_blanks(varstr)         !quiet on default=0.0
  get_iMode=GETREAL(TRIM(varstr),Proposal=0.0_wp,quiet_def_in=.TRUE.)
END FUNCTION get_iMode


!===================================================================================================================================
!> Overwrite axis with average axis by center of closed line of the boundary in each poloidal plane
!!
!===================================================================================================================================
SUBROUTINE InitAverageAxis()
! MODULES
  USE MODgvec_MHD3D_Vars   , ONLY:X1_base,X1_a,X1_b
  USE MODgvec_MHD3D_Vars   , ONLY:X2_base,X2_a,X2_b
  USE MODgvec_MHD3D_Vars   , ONLY:average_axis_move
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT Variables
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL Variables
  REAL(wp) :: X1_b_IP(X1_base%f%mn_nyq(1),X1_base%f%mn_nyq(2))
  REAL(wp) :: X2_b_IP(X2_base%f%mn_nyq(1),X2_base%f%mn_nyq(2))

  INTEGER  :: i_m,i_n
  REAL(wp) :: dl,lint,x1int,x2int
!-----------------------------------------------------------------------------------------------------------------------------------
  ASSOCIATE(m_nyq=>X1_base%f%mn_nyq(1),n_nyq=>X1_base%f%mn_nyq(2))
    X1_b_IP(:,:) = RESHAPE(X1_base%f%evalDOF_IP(0,X1_b),(/m_nyq,n_nyq/))
    X2_b_IP(:,:) = RESHAPE(X2_base%f%evalDOF_IP(0,X2_b),(/m_nyq,n_nyq/))
    DO i_n=1,n_nyq
      dl=X1_b_IP(m_nyq,i_n)*X2_b_IP(1,i_n)-X1_b_IP(1,i_n)*X2_b_IP(m_nyq,i_n)
      lint=dl
      x1int=(X1_b_IP(m_nyq,i_n)+X1_b_IP(1,i_n))*dl
      x2int=(X2_b_IP(m_nyq,i_n)+X2_b_IP(1,i_n))*dl
      DO i_m=2,m_nyq
        dl=SQRT((X1_b_IP(i_m,i_n)-X1_b_IP(i_m-1,i_n))**2+(X2_b_IP(i_m,i_n)-X2_b_IP(i_m-1,i_n))**2)
        dl=X1_b_IP(i_m-1,i_n)*X2_b_IP(i_m,i_n)-X1_b_IP(i_m,i_n)*X2_b_IP(i_m-1,i_n)
        lint=lint+dl
        x1int=x1int+(X1_b_IP(i_m-1,i_n)+X1_b_IP(i_m,i_n))*dl
        x2int=x2int+(X2_b_IP(i_m-1,i_n)+X2_b_IP(i_m,i_n))*dl
      END DO
      ! c_x= 1/(6A) sum_i (x_i-1+x_i)*(x_i-1*y_i - x_i*y_i-1), A=1/2 sum_i  (x_i-1*y_i - x_i*y_i-1)
      X1_b_IP(:,i_n) = x1int/(3.0_wp*lint) + average_axis_move(1)
      X2_b_IP(:,i_n) = x2int/(3.0_wp*lint) + average_axis_move(2)
    END DO
    X1_a = X1_base%f%initDOF(RESHAPE(X1_b_IP,(/X1_base%f%mn_IP/)))
    X2_a = X2_base%f%initDOF(RESHAPE(X2_b_IP,(/X2_base%f%mn_IP/)))
  END ASSOCIATE
END SUBROUTINE InitAverageAxis

!===================================================================================================================================
!> Initialize the solution with the given boundary condition
!!
!===================================================================================================================================
SUBROUTINE InitSolution(U_init,which_init_in)
! MODULES
  USE MODgvec_Globals,       ONLY:ProgressBar,getTime
  USE MODgvec_MHD3D_Vars   , ONLY:init_fromBConly,init_BC,init_average_axis,average_axis_move
  USE MODgvec_MHD3D_Vars   , ONLY:X1_base,X1_BC_Type,X1_a,X1_b
  USE MODgvec_MHD3D_Vars   , ONLY:X2_base,X2_BC_Type,X2_a,X2_b
  USE MODgvec_MHD3D_Vars   , ONLY:LA_base,LA_BC_Type,init_LA
  USE MODgvec_sol_var_MHD3D, ONLY:t_sol_var_mhd3d
  USE MODgvec_lambda_solve,  ONLY:lambda_solve
  USE MODgvec_VMEC_Vars,     ONLY:Rmnc_spl,Rmns_spl,Zmnc_spl,Zmns_spl
  USE MODgvec_VMEC_Vars,     ONLY:lmnc_spl,lmns_spl
  USE MODgvec_VMEC_Readin,   ONLY:lasym
  USE MODgvec_VMEC,          ONLY:VMEC_EvalSplMode
!$ USE omp_lib
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  INTEGER, INTENT(IN) :: which_init_in
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_sol_var_MHD3D), INTENT(INOUT) :: U_init
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER  :: iMode,i_m
  REAL(wp) :: BC_val(2)
  REAL(wp) :: rhopos
  REAL(wp) :: X1_gIP(1:X1_base%s%nBase,1:X1_base%f%modes)
  REAL(wp) :: X2_gIP(1:X2_base%s%nBase,1:X2_base%f%modes)
  REAL(wp) :: LA_gIP(1:LA_base%s%nBase,1:LA_base%f%modes)
!===================================================================================================================================
  IF(.NOT.MPIroot) CALL abort(__STAMP__, &
                       "InitSolution should only be called by MPIroot!")
  X1_gIP=0.0_wp; X2_gIP=0.0_wp; LA_gIP=0.0_wp

  SELECT CASE(which_init_in)
  CASE(-1) !restart
    X1_a(:)=U_init%X1(1,:)
    X2_a(:)=U_init%X2(1,:)
    X1_b(:)=U_init%X1(X1_base%s%nBase,:)
    X2_b(:)=U_init%X2(X2_base%s%nBase,:)
    IF (init_average_axis)THEN
      WRITE(UNIT_stdOut,'(A)') "WARNING: init_average_axis ignored due to restart!"
    END IF
  CASE(0)
    !X1_a,X2_a and X1_b,X2_b already filled from parameter file readin...
    IF(init_average_axis)THEN
      CALL InitAverageAxis()
    END IF !init_average_axis
  CASE(1) !VMEC
    IF((init_BC.EQ.-1).OR.(init_BC.EQ.1))THEN ! compute  axis from VMEC, else use the one defined in paramterfile
      rhopos=0.0_wp
      ASSOCIATE(sin_range    => X1_base%f%sin_range, cos_range    => X1_base%f%cos_range )
      DO imode=cos_range(1)+1,cos_range(2)
        X1_a(iMode:iMode)= VMEC_EvalSplMode(X1_base%f%Xmn(:,iMode),0,(/rhopos/),Rmnc_Spl)
      END DO
      IF(lasym)THEN
        DO imode=sin_range(1)+1,sin_range(2)
          X1_a(iMode:iMode)=X1_a(iMode:imode)  +VMEC_EvalSplMode(X1_base%f%Xmn(:,iMode),0,(/rhopos/),Rmns_Spl)
        END DO
      END IF !lasym
      END ASSOCIATE !X1
      ASSOCIATE(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:iMode)=VMEC_EvalSplMode(X2_base%f%Xmn(:,iMode),0,(/rhopos/),Zmns_Spl)
      END DO
      IF(lasym)THEN
        DO imode=cos_range(1)+1,cos_range(2)
          X2_a(iMode:iMode)=X2_a(iMode:iMode)  +VMEC_EvalSplMode(X2_base%f%Xmn(:,iMode),0,(/rhopos/),Zmnc_Spl)
        END DO
      END IF !lasym
      END ASSOCIATE !X2
    END IF
    IF((init_BC.EQ.-1).OR.(init_BC.EQ.0))THEN ! compute edge from VMEC, else use the one defined in paramterfile
      rhopos=1.0_wp
      ASSOCIATE(sin_range    => X1_base%f%sin_range, cos_range    => X1_base%f%cos_range )
      DO imode=cos_range(1)+1,cos_range(2)
        X1_b(iMode:iMode)=VMEC_EvalSplMode(X1_base%f%Xmn(:,iMode),0,(/rhopos/),Rmnc_Spl)
      END DO
      IF(lasym)THEN
        DO imode=sin_range(1)+1,sin_range(2)
          X1_b(iMode:iMode)=X1_b(iMode:iMode)  +VMEC_EvalSplMode(X1_base%f%Xmn(:,iMode),0,(/rhopos/),Rmns_Spl)
        END DO
      END IF !lasym
      END ASSOCIATE !X1
      ASSOCIATE(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:iMode)=VMEC_EvalSplMode(X2_base%f%Xmn(:,iMode),0,(/rhopos/),Zmns_Spl)
      END DO
      IF(lasym)THEN
        DO imode=cos_range(1)+1,cos_range(2)
          X2_b(iMode:iMode)=X2_b(iMode:iMode)  +VMEC_EvalSplMode(X2_base%f%Xmn(:,iMode),0,(/rhopos/),Zmnc_Spl)
        END DO
      END IF !lasym
      END ASSOCIATE !X2
    END IF
    IF(init_average_axis)THEN
      CALL InitAverageAxis()
    END IF !init_average_axis
    IF(.NOT.init_fromBConly)THEN !only boundary and axis from VMEC
      ASSOCIATE(s_IP         => X1_base%s%s_IP, &
                nBase        => X1_base%s%nBase, &
                sin_range    => X1_base%f%sin_range,&
                cos_range    => X1_base%f%cos_range )
      DO imode=cos_range(1)+1,cos_range(2)
        X1_gIP(:,iMode)  =VMEC_EvalSplMode(X1_base%f%Xmn(:,iMode),0,s_IP,Rmnc_Spl)
      END DO !imode=cos_range
      IF(lasym)THEN
        DO imode=sin_range(1)+1,sin_range(2)
          X1_gIP(:,iMode)=X1_gIP(:,iMode)  +VMEC_EvalSplMode(X1_base%f%Xmn(:,iMode),0,s_IP,Rmns_Spl)
        END DO !imode= sin_range
      END IF !lasym
      END ASSOCIATE !X1
      ASSOCIATE(s_IP         => X2_base%s%s_IP, &
                nBase        => X2_base%s%nBase, &
                sin_range    => X2_base%f%sin_range,&
                cos_range    => X2_base%f%cos_range )
      DO imode=sin_range(1)+1,sin_range(2)
        X2_gIP(:,iMode)=VMEC_EvalSplMode(X2_base%f%Xmn(:,iMode),0,s_IP,Zmns_Spl)
      END DO !imode=sin_range
      IF(lasym)THEN
        DO imode=cos_range(1)+1,cos_range(2)
          X2_gIP(:,iMode)=X2_gIP(:,iMode)  +VMEC_EvalSplMode(X2_base%f%Xmn(:,iMode),0,s_IP,Zmnc_Spl)
        END DO !imode= sin_range
      END IF !lasym
      END ASSOCIATE !X2
      ASSOCIATE(s_IP         => LA_base%s%s_IP, &
                nBase        => LA_base%s%nBase, &
                sin_range    => LA_base%f%sin_range,&
                cos_range    => LA_base%f%cos_range )
      DO imode=sin_range(1)+1,sin_range(2)
        LA_gIP(:,iMode)=VMEC_EvalSplMode(LA_base%f%Xmn(:,iMode),0,s_IP,lmns_Spl)
      END DO !imode= sin_range
      IF(lasym)THEN
        DO imode=cos_range(1)+1,cos_range(2)
          LA_gIP(:,iMode)=LA_gIP(:,iMode)  +VMEC_EvalSplMode(LA_base%f%Xmn(:,iMode),0,s_IP,lmnc_Spl)
        END DO !imode=cos_range
      END IF !lasym
      END ASSOCIATE !LA
     END IF !fullIntVmec
  END SELECT !which_init


  IF((which_init_in.NE.-1).AND.(init_fromBConly))THEN
    !no restart(=-1) and initialization only
    !smoothly interpolate between  edge and axis data
    ASSOCIATE(s_IP         =>X1_base%s%s_IP, &
              modes        =>X1_base%f%modes, &
              zero_odd_even=>X1_base%f%zero_odd_even)
    DO imode=1,modes
      SELECT CASE(zero_odd_even(iMode))
      CASE(MN_ZERO,M_ZERO) !X1_a only used here!!
        X1_gIP(:,iMode)=(1.0_wp-(s_IP(:)**2))*X1_a(iMode)+(s_IP(:)**2)*X1_b(iMode)  ! meet edge and axis, ~(1-s^2)
      CASE(M_ODD_FIRST)
        X1_gIP(:,iMode)=s_IP(:)*X1_b(iMode)      ! first odd mode ~s
      CASE DEFAULT ! ~s^m
        X1_gIP(:,iMode)=s_IP(:)*X1_b(iMode)
        DO i_m=2,X1_base%f%Xmn(1,iMode)
          X1_gIP(:,iMode)=X1_gIP(:,iMode)*s_IP(:)
        END DO
      END SELECT !X1(:,iMode) zero odd even
    END DO !iMode
    END ASSOCIATE

    ASSOCIATE(s_IP         =>X2_base%s%s_IP, &
              modes        =>X2_base%f%modes, &
              zero_odd_even=>X2_base%f%zero_odd_even)
    DO imode=1,modes
      SELECT CASE(zero_odd_even(iMode))
      CASE(MN_ZERO,M_ZERO) !X2_a only used here!!!
        X2_gIP(:,iMode)=(1.0_wp-(s_IP(:)**2))*X2_a(iMode)+(s_IP(:)**2)*X2_b(iMode) ! meet edge and axis, ~(1-s^2)
      CASE(M_ODD_FIRST)
        X2_gIP(:,iMode)=s_IP(:)*X2_b(iMode)      ! first odd mode ~s
      CASE DEFAULT ! ~s^m
        X2_gIP(:,iMode)=s_IP(:)*X2_b(iMode)
        DO i_m=2,X2_base%f%Xmn(1,iMode)
          X2_gIP(:,iMode)=X2_gIP(:,iMode)*s_IP(:)
        END DO
      END SELECT !X2(:,iMode) zero odd even
    END DO
    END ASSOCIATE
  END IF !init_fromBConly

  !apply strong boundary conditions
  ASSOCIATE(modes        =>X1_base%f%modes, &
            zero_odd_even=>X1_base%f%zero_odd_even)
  DO imode=1,modes
    SELECT CASE(zero_odd_even(iMode))
    CASE(MN_ZERO,M_ZERO)
      BC_val =(/ X1_a(iMode)    ,      X1_b(iMode)/)
    !CASE(M_ODD_FIRST,M_ODD,M_EVEN)
    CASE DEFAULT
      BC_val =(/          0.0_wp,      X1_b(iMode)/)
    END SELECT !X1(:,iMode) zero odd even
    IF(which_init_in.NE.-1) U_init%X1(:,iMode)=X1_base%s%initDOF( X1_gIP(:,iMode) )
    CALL X1_base%s%applyBCtoDOF(U_init%X1(:,iMode),X1_BC_type(:,iMode),BC_val)
  END DO
  END ASSOCIATE !X1

  ASSOCIATE(modes        =>X2_base%f%modes, &
            zero_odd_even=>X2_base%f%zero_odd_even)
  DO imode=1,modes
    SELECT CASE(zero_odd_even(iMode))
    CASE(MN_ZERO,M_ZERO)
      BC_val =(/     X2_a(iMode),      X2_b(iMode)/)
    !CASE(M_ODD_FIRST,M_ODD,M_EVEN)
    CASE DEFAULT
      BC_val =(/          0.0_wp,      X2_b(iMode)/)
    END SELECT !X1(:,iMode) zero odd even
    IF(which_init_in.NE.-1) U_init%X2(:,iMode)=X2_base%s%initDOF( X2_gIP(:,iMode) )
    CALL X2_base%s%applyBCtoDOF(U_init%X2(:,iMode),X2_BC_type(:,iMode),BC_val)
  END DO
  END ASSOCIATE !X2

  ASSOCIATE(modes        =>LA_base%f%modes, &
            zero_odd_even=>LA_base%f%zero_odd_even)
  IF((which_init_in.NE.-1).AND.(.NOT.init_LA)) THEN
    !lambda init might not be needed since it has no boundary condition and changes anyway after the update of the mapping...
    IF(.NOT.init_fromBConly)THEN
      WRITE(UNIT_stdOut,'(4X,A)') "... lambda initialized with VMEC ..."
    ELSE
      WRITE(UNIT_stdOut,'(4X,A)') "... initialize lambda =0 ..."
      LA_gIP=0.0_wp
    END IF
  END IF !!which_init_in>-1 and not init_LA
  !always apply strong BC
  DO imode=1,modes
    IF(zero_odd_even(iMode).EQ.MN_ZERO)THEN
      U_init%LA(:,iMode)=0.0_wp ! (0,0) mode should not be here, but must be zero if its used.
    ELSE
      BC_val =(/ 0.0_wp, 0.0_wp/)
      IF(which_init_in.NE.-1) U_init%LA(:,iMode)=LA_base%s%initDOF( LA_gIP(:,iMode) )
      CALL LA_base%s%applyBCtoDOF(U_init%LA(:,iMode),LA_BC_type(:,iMode),BC_val)
    END IF!iMode ~ MN_ZERO
  END DO !iMode

  END ASSOCIATE !LA
END SUBROUTINE InitSolution


!===================================================================================================================================
!> Initialize LAMBDA FROM U_init%X1,%X2 and iota profile, this computation is distributed over MPIranks
!!
!===================================================================================================================================
SUBROUTINE Init_LA_from_Solution(U_init)
! MODULES
  USE MODgvec_Globals,       ONLY:ProgressBar,getTime,myRank,nRanks
  USE MODgvec_MHD3D_Vars   , ONLY:X1_base,X2_base,LA_base,LA_BC_Type,hmap, Phi_profile, chi_profile
  USE MODgvec_sol_var_MHD3D, ONLY:t_sol_var_mhd3d
  USE MODgvec_lambda_solve,  ONLY:lambda_solve
  USE MODgvec_MPI           ,ONLY:par_reduce,par_BCast
  USE MODgvec_hmap          ,ONLY:hmap_new_auxvar,PP_T_HMAP_AUXVAR
!$ USE omp_lib
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_sol_var_MHD3D), INTENT(INOUT) :: U_init
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER  :: iMode,is,ns_str,ns_end,iRank,nBase
  REAL(wp) :: BC_val(2),rhopos
  REAL(wp) :: StartTime,EndTime
  REAL(wp),DIMENSION(1:LA_base%s%nBase):: PhiPrime,chiPrime
  REAL(wp) :: LA_gIP(1:LA_base%s%nBase,1:LA_base%f%modes)
#ifdef PP_WHICH_HMAP
  TYPE(PP_T_HMAP_AUXVAR),ALLOCATABLE  :: hmap_xv(:) !! auxiliary variables for hmap
#else
  CLASS(PP_T_HMAP_AUXVAR),ALLOCATABLE  :: hmap_xv(:) !! auxiliary variables for hmap
#endif
!===================================================================================================================================
  StartTime=GetTime()
  SWRITE(UNIT_stdOut,'(4X,A)') "... Initialize lambda from mapping ..."
  CALL enter_subregion("reinit-lambda")
  nBase        = LA_base%s%nBase
  ASSOCIATE(modes        => LA_base%f%modes, &
            s_IP         => LA_base%s%s_IP, &
            zero_odd_even=> LA_base%f%zero_odd_even, &
            modes_str    => LA_base%f%modes_str, &
            modes_end    => LA_base%f%modes_end, &
            offset_modes => LA_Base%f%offset_modes )
  !evaluate profiles only in MPIroot!
  IF(MPIroot)THEN
    DO is=1,nBase
      rhopos=MIN(1.0_wp-1.0e-12_wp,MAX(1.0e-4_wp,s_IP(is))) !exclude axis
      phiPrime(is)=Phi_profile%eval_at_rho(rhopos,deriv=1)
      chiPrime(is)=chi_profile%eval_at_rho(rhopos,deriv=1)
    END DO
  END IF !MPIroot
  CALL par_BCast(phiPrime,0)
  CALL par_BCast(chiPrime,0)
  !initialize Lambda, radially parallel
  ns_str = (nBase*(myRank  ))/nRanks+1
  ns_end = (nBase*(myRank+1))/nRanks
  LA_gIP=0.0_wp
  CALL ProgressBar(0,ns_end) !init

  CALL hmap_new_auxvar(hmap,X1_base%f%x_IP(2,:),hmap_xv,.FALSE.) !no 2nd derivative needed
  DO is=ns_str,ns_end
    rhopos=MIN(1.0_wp-1.0e-12_wp,MAX(1.0e-4_wp,s_IP(is))) !exclude axis
    CALL lambda_Solve(rhopos,hmap,hmap_xv,X1_base,X2_base,LA_base%f,U_init%X1,U_init%X2,LA_gIP(is,:),phiPrime(is),chiPrime(is))
    CALL ProgressBar(is,ns_end)
  END DO !is
  DEALLOCATE(hmap_xv)
!!!  CALL par_reduce(LA_gIP,'SUM',0)
!!!  IF(MPIroot)THEN
!!!    DO iMode=1,modes
!!!      IF(zero_odd_even(iMode).EQ.MN_ZERO)THEN
!!!        U_init%LA(:,iMode)=0.0_wp ! (0,0) mode should not be here, but must be zero if its used.
!!!      ELSE
!!!        U_init%LA(:,iMode)=LA_base%s%initDOF( LA_gIP(:,iMode) )
!!!      END IF!iMode ~ MN_ZERO
!!!      BC_val =(/ 0.0_wp, 0.0_wp/)
!!!      CALL LA_base%s%applyBCtoDOF(U_init%LA(:,iMode),LA_BC_type(:,iMode),BC_val)
!!!    END DO !iMode=1,modes
!!!  END IF
!!!  CALL par_BCast(U_init%LA,0)
  !reduce radially, different mode sets to different MPIranks (should be a gatherv)
  DO iRank=0,nRanks-1
    IF(offset_modes(iRank+1)-offset_modes(iRank).GT.0) &
      CALL par_Reduce(LA_gIP(1:nbase,offset_modes(iRank)+1:offset_modes(iRank+1)),'SUM',iRank)
  END DO
  DO iMode=modes_str,modes_end
    IF(zero_odd_even(iMode).EQ.MN_ZERO)THEN
      U_init%LA(1:nBase,iMode)=0.0_wp ! (0,0) mode should not be here, but must be zero if its used.
    ELSE
      U_init%LA(1:nBase,iMode)=LA_base%s%initDOF( LA_gIP(1:nBase,iMode) )
    END IF!iMode ~ MN_ZERO
    BC_val =(/ 0.0_wp, 0.0_wp/)
    CALL LA_base%s%applyBCtoDOF(U_init%LA(:,iMode),LA_BC_type(:,iMode),BC_val)
  END DO !iMode=modes_str, modes_end
  ! broadcast result: different mode ranges to different MPIranks
  DO iRank=0,nRanks-1
    IF(offset_modes(iRank+1)-offset_modes(iRank).GT.0) &
      CALL par_Bcast(U_init%LA(1:nBase,offset_modes(iRank)+1:offset_modes(iRank+1)),iRank)
  END DO
  END ASSOCIATE !LA
  EndTime=GetTime()
  CALL exit_subregion("reinit-lambda")
  SWRITE(UNIT_stdOut,'(4X,A,F9.2,A)') " init lambda took [ ",EndTime-StartTime," sec]"
END SUBROUTINE Init_LA_from_solution

!===================================================================================================================================
!> Add boundary perturbation
!!
!===================================================================================================================================
SUBROUTINE AddBoundaryPerturbation(U_init, depth, blend_type)
! MODULES
  USE MODgvec_MHD3D_Vars   , ONLY:X1_base,X1_BC_Type,X1_a,X1_b,X1pert_b
  USE MODgvec_MHD3D_Vars   , ONLY:X2_base,X2_BC_Type,X2_a,X2_b,X2pert_b
  USE MODgvec_sol_var_MHD3D, ONLY:t_sol_var_mhd3d
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  REAL(wp),INTENT(IN) :: depth ! depth of perturbation from boundary (0.1..0.3)
  INTEGER, INTENT(IN) :: blend_type ! 0/BLEND_LEGACY: legacy Gaussian, 1/BLEND_COSM: new cosine
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_sol_var_MHD3D), INTENT(INOUT) :: U_init
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER  :: iMode
  REAL(wp) :: BC_val(2)
  REAL(wp) :: X1pert_gIP(1:X1_base%s%nBase)
  REAL(wp) :: X2pert_gIP(1:X2_base%s%nBase)
!===================================================================================================================================
  IF(.NOT.MPIroot) CALL abort(__STAMP__, &
                       "AddBoundaryPerturbation should only be called by MPIroot!")
  WRITE(UNIT_stdOut,'(4X,A)') "ADD BOUNDARY PERTURBATION..."

  ASSOCIATE(s_IP         =>X1_base%s%s_IP, &
            modes        =>X1_base%f%modes )
  DO imode=1,modes

    X1_b(iMode)=X1_b(iMode)+X1pert_b(iMode)
    X1pert_gIP(:)=blend(s_IP, depth, X1_base%f%Xmn(1, iMode))*X1pert_b(iMode)
    U_init%X1(:,iMode)=U_init%X1(:,iMode) + X1_base%s%initDOF( X1pert_gIP(:) )
  END DO
  END ASSOCIATE

  ASSOCIATE(s_IP         =>X2_base%s%s_IP, &
            modes        =>X2_base%f%modes )
  DO imode=1,modes
    X2_b(iMode)=X2_b(iMode)+X2pert_b(iMode)
    X2pert_gIP(:)=blend(s_IP, depth, X2_base%f%Xmn(1, iMode))*X2pert_b(iMode)
    U_init%X2(:,iMode)=U_init%X2(:,iMode) + X2_base%s%initDOF( X2pert_gIP(:))
  END DO
  END ASSOCIATE

  !apply strong boundary conditions
  ASSOCIATE(modes        =>X1_base%f%modes, &
            zero_odd_even=>X1_base%f%zero_odd_even)
  DO imode=1,modes
    SELECT CASE(zero_odd_even(iMode))
    CASE(MN_ZERO,M_ZERO)
      BC_val =(/ X1_a(iMode)    ,      X1_b(iMode)/)
    !CASE(M_ODD_FIRST,M_ODD,M_EVEN)
    CASE DEFAULT
      BC_val =(/          0.0_wp,      X1_b(iMode)/)
    END SELECT !X1(:,iMode) zero odd even
    CALL X1_base%s%applyBCtoDOF(U_init%X1(:,iMode),X1_BC_type(:,iMode),BC_val)
  END DO
  END ASSOCIATE !X1

  ASSOCIATE(modes        =>X2_base%f%modes, &
            zero_odd_even=>X2_base%f%zero_odd_even)
  DO imode=1,modes
    SELECT CASE(zero_odd_even(iMode))
    CASE(MN_ZERO,M_ZERO)
      BC_val =(/     X2_a(iMode),      X2_b(iMode)/)
    !CASE(M_ODD_FIRST,M_ODD,M_EVEN)
    CASE DEFAULT
      BC_val =(/          0.0_wp,      X2_b(iMode)/)
    END SELECT !X1(:,iMode) zero odd even
    CALL X2_base%s%applyBCtoDOF(U_init%X2(:,iMode),X2_BC_type(:,iMode),BC_val)
  END DO
  END ASSOCIATE !X2

  WRITE(UNIT_stdOut,'(4X,A)') "... DONE."
  WRITE(UNIT_stdOut,fmt_sep)

  CONTAINS

  ELEMENTAL FUNCTION blend(s_in, depth, m)
    USE MODgvec_Globals, ONLY: wp, PI
    USE MODgvec_MHD3D_Vars, ONLY: BLEND_LEGACY
    REAL(wp),INTENT(IN) :: s_in !input coordinate [0,1]
    REAL(wp)            :: blend
    REAL(wp),INTENT(IN) :: depth
    INTEGER,INTENT(IN)  :: m     ! exponent for cosine blending (poloidal mode number)
    ASSOCIATE(shift => 1.0_wp - depth)
      IF (blend_type == BLEND_LEGACY) THEN
        blend = EXP(-4.0_wp * ((s_in - 1.0_wp) / depth)**2)
      ELSE IF (s_in .GE. shift) THEN
        blend = (COS(((s_in - shift) / (1.0_wp - shift) - 1.0_wp) * PI) / 2.0_wp + 0.5_wp)**m
      ELSE IF (m .EQ. 0) THEN
        blend = 1.0_wp
      ELSE
        blend = 0.0_wp
      END IF
    END ASSOCIATE
  END FUNCTION blend

END SUBROUTINE AddBoundaryPerturbation

!===================================================================================================================================
!> Compute Equilibrium, iteratively
!!
!===================================================================================================================================

SUBROUTINE MinimizeMHD3D(sf, dt_in)
! MODULES
  USE MODgvec_MHD3D_vars, ONLY: maxIter, minimize_tol, MinimizerType
  USE MODgvec_MHD3D_minimize, ONLY: t_minimizer_mhd3d, new_minimizer
  IMPLICIT NONE
  REAL(wp), OPTIONAL, INTENT(IN) :: dt_in
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
!===================================================================================================================================
  __PERFON('minimizer')
  CALL enter_subregion("Minimization")
  IF (PRESENT(dt_in)) THEN
      sf%minimizer%vars%dt = dt_in
  ENDIF
  CALL sf%minimizer%minimize(minimize_tol, maxIter)
  CALL exit_subregion("Minimization")
  __PERFOFF('minimizer')
END SUBROUTINE MinimizeMHD3D



!===================================================================================================================================
!> Finalize Module
!!
!===================================================================================================================================
SUBROUTINE FinalizeMHD3D(sf)
! MODULES
  USE MODgvec_MHD3D_Vars
  USE MODgvec_MHD3D_EvalFunc,ONLY:FinalizeMHD3D_EvalFunc
  USE MODgvec_VMEC,ONLY: FinalizeVMEC
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
!===================================================================================================================================
  CALL enter_subregion("finalize-MHD3D")
  IF(ALLOCATED(X1_base)) CALL X1_base%free()
  IF(ALLOCATED(X2_base)) CALL X2_base%free()
  IF(ALLOCATED(LA_base)) CALL LA_base%free()

  CALL sgrid%free()
  IF(ALLOCATED(BFF)) THEN
    CALL BFF%free()
    DEALLOCATE(BFF)
  END IF

  IF(ALLOCATED(sf%minimizer)) CALL sf%minimizer%free()

  SDEALLOCATE(X1_BC_type)
  SDEALLOCATE(X2_BC_type)
  SDEALLOCATE(LA_BC_type)
  SDEALLOCATE(X1_b)
  SDEALLOCATE(X2_b)
  SDEALLOCATE(X1pert_b)
  SDEALLOCATE(X2pert_b)
  SDEALLOCATE(LA_b)
  SDEALLOCATE(X1_a)
  SDEALLOCATE(X2_a)

  SDEALLOCATE(iota_profile)
  SDEALLOCATE(pres_profile)
  SDEALLOCATE(Phi_profile)
  SDEALLOCATE(chi_profile)

  CALL FinalizeMHD3D_EvalFunc()
  IF(which_init.EQ.1) CALL FinalizeVMEC()

  SDEALLOCATE(hmap)
  SDEALLOCATE(hmap_auxvar)

  SDEALLOCATE(X1_base)
  SDEALLOCATE(X2_base)
  SDEALLOCATE(LA_base)
  CALL exit_subregion("finalize-MHD3D")
END SUBROUTINE FinalizeMHD3D

END MODULE MODgvec_MHD3D
mhd3d.F90 Source File

Contents

Modules

Source Code

This file depends on

Files dependent on this one

Source Code
