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

!===================================================================================================================================
!>
!!# 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
  USE MODgvec_c_functional,   ONLY: t_functional
  IMPLICIT NONE

  PRIVATE
  PUBLIC t_functional_mhd3d

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

  TYPE,EXTENDS(t_functional) :: t_functional_mhd3d
    !-------------------------------------------------------------------------------------------------------------------------------
    LOGICAL :: initialized
    !-------------------------------------------------------------------------------------------------------------------------------
    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: t_boundaryFromFile,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
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER          :: i,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
  CLASS(t_boundaryFromFile),ALLOCATABLE:: BFF
  CHARACTER(LEN=255) ::boundary_filename
!===================================================================================================================================
  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)


  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=4
    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)/))
    mn_nyq  = GETINTARRAY("mn_nyq",2)
    IF(mn_nyq(1).LT.mn_nyq_min(1))THEN
       CALL abort(__STAMP__,&
                  'mn_nyq(1) too small, should be >= ',intInfo=mn_nyq_min(1))
    END IF
    IF(mn_nyq(2).LT.mn_nyq_min(2))THEN
       CALL abort(__STAMP__,&
                  'mn_nyq(2) too small, should be >= ',intInfo=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)/))
  END IF

  SWRITE(UNIT_stdOut,*)
  SWRITE(UNIT_stdOut,'(A,I4,A,I6," , ",I6,A)')'    fac_nyq = ', fac_nyq,'  ==> 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

  !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)
      CALL BFF%free()
    END IF
  END IF !MPIroot



  !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

  boundary_perturb=GETLOGICAL('boundary_perturb',Proposal=.FALSE.)
  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

  ! ALLOCATE DATA
  ALLOCATE(U(-3:1))
  CALL U(1)%init((/X1_base%s%nbase,X2_base%s%nbase,LA_base%s%nBase,  &
                   X1_base%f%modes,X2_base%f%modes,LA_base%f%modes/)  )
  DO i=-3,0
    CALL U(i)%copy(U(1))
  END DO
  ALLOCATE(F(-1:0))
  DO i=-1,0
    CALL F(i)%copy(U(1))
  END DO
  ALLOCATE(V(-1:1))
  DO i=-1,1
    CALL V(i)%copy(U(1))
  END DO
  ALLOCATE(P(-1:1))
  DO i=-1,1
    CALL P(i)%copy(U(1))
  END DO

  CALL InitializeMHD3D_EvalFunc()

  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
  !===================================================================================================================================
  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)
    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 (n_profile_vals .NE. n_profile_rho2) THEN
      CALL abort(__STAMP__,&
      'Size of '//var//'_rho2 and '//var//'_vals must be equal!')
    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 can only be 'not_a_knot', '1st_deriv' or '2nd_deriv'!")
    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. It must be either "polynomial", "bspline" or "interpolation".')
  END IF ! profile type

  SDEALLOCATE(profile_knots)
  SDEALLOCATE(profile_coefs)
END SUBROUTINE InitProfile

!===================================================================================================================================
!> Initialize Module
!!
!===================================================================================================================================
SUBROUTINE InitSolutionMHD3D(sf)
! MODULES
  USE MODgvec_MHD3D_Vars     , ONLY: which_init,U,F,init_LA,boundary_perturb
  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="                           ...")
  IF(MPIroot) THEN
    IF(doRestart)THEN
      WRITE(UNIT_stdOut,'(4X,A)')'... restarting from file ... '
      CALL RestartFromState(RestartFile,U(0))
      CALL InitSolution(U(0),-1) ! (re-)apply BC and init LA (if init_LA is true)
    ELSE
      CALL InitSolution(U(0),which_init)
    END IF
    IF(boundary_perturb)THEN
      CALL AddBoundaryPerturbation(U(0),0.3_wp)
    END IF !boundary_perturb
  END IF !MPIroot
  CALL par_Bcast(U(0)%X1,0)
  CALL par_Bcast(U(0)%X2,0)

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

  CALL U(-1)%set_to(U(0))

  JacCheck=2
  CALL InitProfilesGP() !evaluate profiles once at Gauss Points (on MPIroot + BCast)
  U(0)%W_MHD3D=EvalEnergy(U(0),.TRUE.,JacCheck)
  IF(JacCheck.EQ.-1)THEN
    CALL Analyze(0)
    CALL abort(__STAMP__,&
        "NEGATIVE JACOBIAN FOUND AFTER INITIALIZATION!")
  END IF
  CALL WriteState(U(0),0)
  CALL EvalForce(U(0),.FALSE.,JacCheck, F(0))
  SWRITE(UNIT_stdOut,'(8x,A,3E11.4)')'|Force|= ',SQRT(F(0)%norm_2())
  CALL Analyze(0)

  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 ..."
  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()
  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,h)
! 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) :: h ! depth of perturbation from boundary (0.1..0.3)
!-----------------------------------------------------------------------------------------------------------------------------------
! 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)*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)*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)
    REAL(wp),INTENT(IN) :: s_in !input coordinate [0,1]
    REAL(wp)            :: blend
    !blend= ( MIN(0., (s_in-(1.-h)) ) / h) **4
    !blend= 2. -2./(EXP(8.*(s_in-1.)/h) +1)
    !blend= EXP(-4.*((s_in-1.)/h)**2)
    !blend= s_in
    blend= EXP(-4.0_wp*((s_in-1.0_wp)/0.6_wp)**2)
  END FUNCTION blend

END SUBROUTINE AddBoundaryPerturbation

!===================================================================================================================================
!> Compute Equilibrium, iteratively
!!
!===================================================================================================================================
SUBROUTINE MinimizeMHD3D(sf)
! MODULES
  USE MODgvec_MHD3D_vars, ONLY: MinimizerType
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
!===================================================================================================================================
  __PERFON('minimizer')
  SELECT CASE(MinimizerType)
  CASE(0,10)
    CALL MinimizeMHD3D_descent(sf)
  CASE DEFAULT
    CALL abort(__STAMP__,&
        "Minimizertype does not exist",MinimizerType,-1.0_wp)
  END SELECT
  __PERFOFF('minimizer')
END SUBROUTINE MinimizeMHD3D

!===================================================================================================================================
!> Compute Equilibrium, iteratively
!!
!===================================================================================================================================
SUBROUTINE MinimizeMHD3D_descent(sf)
! MODULES
  USE MODgvec_MHD3D_Vars
  USE MODgvec_MHD3D_EvalFunc
  USE MODgvec_Analyze, ONLY:analyze
  USE MODgvec_Restart, ONLY:WriteState
  USE MODgvec_MHD3D_visu, ONLY:WriteSFLoutfile
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER   :: iter,nStepDecreased,nSkip_Jac,nSkip_dw
  INTEGER   :: JacCheck,lastoutputIter,StartTimeArray(8)
  REAL(wp)  :: dt,deltaW,absTol
  INTEGER,PARAMETER   :: ndamp=10
  REAL(wp)  :: tau(1:ndamp), tau_bar
  REAL(wp)  :: min_dt_out,max_dt_out,min_dw_out,max_dw_out,sum_dW_out,t_pseudo,Fnorm(3),Vnorm(3),Fnorm0(3),Fnorm_old(3),W_MHD3D_0
  INTEGER   :: logUnit !globally needed for logging
  INTEGER   :: logiter_ramp,logscreen
  LOGICAL   :: restart_iter
  LOGICAL   :: first_iter
!===================================================================================================================================
  SWRITE(UNIT_stdOut,'(A)') "MINIMIZE MHD3D FUNCTIONAL..."


  abstol=minimize_tol


  dt=start_dt
  nstepDecreased=0
  nSkip_Jac=0
  t_pseudo=0
  lastOutputIter=0
  iter=0
  Vnorm=0.0_wp
  logiter_ramp=1
  logscreen=1

  first_iter=.TRUE.
  restart_iter=.FALSE.

  CALL U(-3)%set_to(U(0)) !initial state, should remain unchanged

  DO WHILE(iter.LT.maxIter)
    IF((first_iter).OR.(restart_iter))THEN
      JacCheck=1 !abort if detJ<0
      CALL EvalAux(           U(0),JacCheck)
      U(0)%W_MHD3D=EvalEnergy(U(0),.FALSE.,JacCheck)
      W_MHD3D_0 = U(0)%W_MHD3D
      CALL EvalForce(         U(0),.FALSE.,JacCheck,F(0))
      Fnorm0=SQRT(F(0)%norm_2())
      Fnorm=Fnorm0
      Fnorm_old=1.1_wp*Fnorm0
      CALL U(-1)%set_to(U(0)) !last state
      CALL U(-2)%set_to(U(0)) !state at last logging interval
      !for hirshman method
      IF(MinimizerType.EQ.10)THEN
        CALL V(-1)%set_to(0.0_wp)
        CALL V( 0)%set_to(0.0_wp)
        tau(1:ndamp)=0.15_wp/dt
        tau_bar = 0.075_wp
      END IF
      min_dt_out=1.0e+30_wp
      max_dt_out=0.0_wp
      min_dW_out=1.0e+30_wp
      max_dW_out=-1.0e+30_wp
      sum_dW_out=0.0_wp
      nSkip_dW =0
      IF(restart_iter) restart_iter=.FALSE.
      IF(first_iter)THEN
        CALL StartLogging()
        first_iter=.FALSE.
      END IF
    END IF !before first iteration or after restart Jac<0

    !COMPUTE NEW SOLUTION P(1) as a prediction

    SELECT CASE(MinimizerType)
    CASE(0) !gradient descent, previously used for minimizerType=0
      CALL P(1)%AXBY(1.0_wp,U(0),dt,F(0)) !overwrites P(1), predicts solution U(1)
    CASE(10) !hirshman method
      !tau is damping parameter
      tau(1:ndamp-1) = tau(2:ndamp) !save old
      tau(ndamp)  = MIN(0.15_wp,ABS(LOG(SUM(Fnorm**2)/SUM(Fnorm_old**2))))/dt  !ln(|F_n|^2/|F_{n-1}|^2), Fnorm=|F_X1|,|F_X2|,|F_LA|
      tau_bar = 0.5_wp*dt*SUM(tau)/REAL(ndamp,wp)   !=1/2 * tauavg
      CALL V(1)%AXBY(((1.0_wp-tau_bar)/(1.0_wp+tau_bar)),V(0),(dt/(1.0_wp+tau_bar)),F(0)) !velocity V(1)
      CALL P(1)%AXBY(1.0_wp,U(0),dt,V(1)) !overwrites P(1), predicst solution U(1)
      Vnorm=SQRT(V(1)%norm_2())
    END SELECT


    JacCheck=2 !no abort,if detJ<0, JacCheck=-1
    P(1)%W_MHD3D=EvalEnergy(P(1),.TRUE.,JacCheck)
    IF(JacCheck.EQ.-1)THEN
      dt=0.9_wp*dt
      nstepDecreased=nStepDecreased+1
      nSkip_Jac=nSkip_Jac+1
      restart_iter=.TRUE.
      CALL U(0)%set_to(U(-3)) !reset to initial state
      SWRITE(UNIT_stdOut,'(8X,I8,A,E11.4,A)')iter,'...detJac<0, decrease stepsize to dt=',dt,  ' and RESTART simulation!!!!!!!'
    ELSE
      !detJ>0
      deltaW=P(1)%W_MHD3D-U(0)%W_MHD3D!should be <=0,
      IF(deltaW.LE.dW_allowed*W_MHD3D_0)THEN !valid step /hirshman method accept W increase!

        IF(ALL(Fnorm.LE.abstol))THEN
          CALL Logging(.FALSE.)
          SWRITE(UNIT_stdOut,'(4x,A)')'==>Iteration finished, |force| in relative tolerance'
          EXIT !DO LOOP
        END IF
        iter=iter+1
        t_pseudo=t_pseudo+dt
        ! for simple gradient & hirshman
        CALL U(-1)%set_to(U(0))
        CALL U(0)%set_to(P(1))
        ! for hirshman method
        IF(MinimizerType.EQ.10)THEN
          CALL V(-1)%set_to(V(0))
          CALL V(0)%set_to(V(1))
        END IF

        CALL EvalForce(P(1),.FALSE.,JacCheck,F(0)) !evalAux was already called on P(1)=U(0), so that its set false here.
        Fnorm_old=Fnorm
        Fnorm=SQRT(F(0)%norm_2())

        nstepDecreased=0
        min_dt_out=MIN(min_dt_out,dt)
        max_dt_out=MAX(max_dt_out,dt)
        min_dW_out=MIN(min_dW_out,deltaW)
        max_dW_out=MAX(max_dW_out,deltaW)
        sum_dW_out=sum_dW_out+deltaW
        IF(MOD(iter,logIter_ramp).EQ.0)THEN

          CALL Logging(.NOT.((logIter_ramp.GE.logIter).AND.(MOD(logscreen,nLogScreen).EQ.0)))
          IF(.NOT.(logIter_ramp.LT.logIter))THEN !only reset for logIter
            logscreen=logscreen+1
            min_dt_out=1.0e+30_wp
            max_dt_out=0.0_wp
            min_dW_out=1.0e+30_wp
            max_dW_out=-1.0e+30_wp
            sum_dW_out=0.0_wp
            nSkip_dW =0
          END IF
          logIter_ramp=MIN(logIter,logIter_ramp*2)
        END IF

      ELSE !not a valid step, decrease timestep and skip P(1)
        dt=0.9_wp*dt
        nstepDecreased=nStepDecreased+1
        nSkip_dW=nSkip_dW+1
        !CALL U(0)%set_to(U(-2))
        restart_iter=.TRUE.
        SWRITE(UNIT_stdOut,'(8X,I8,A,E8.1,A,E11.4)')iter,'...deltaW>',dW_allowed,'*W_MHD3D_0, skip step and decrease stepsize to dt=',dt
      END IF
    END IF !JacCheck

    IF(nStepDecreased.GT.130) THEN ! 0.9^130 ~10^-6
      SWRITE(UNIT_stdOut,'(A,E21.11)')'Iteration stopped since timestep has been decreased by 0.9^130: ', dt
      SWRITE(UNIT_stdOut,fmt_sep)
      RETURN
    END IF
    IF((MOD(iter,outputIter).EQ.0).AND.(lastoutputIter.NE.iter))THEN
      __PERFON('output')
      SWRITE(UNIT_stdOut,'(A)')'##########################  OUTPUT ##################################'
      CALL Analyze(iter)
      CALL WriteState(U(0),iter)
      SWRITE(UNIT_stdOut,'(A)')'#####################################################################'
      lastOutputIter=iter
      __PERFOFF('output')
    END IF
  END DO !iter
  IF(iter.GE.MaxIter)THEN
    SWRITE(UNIT_stdOut,'(A,E21.11)')"maximum iteration count exceeded"
  END IF
  SWRITE(UNIT_stdOut,'(A)') "... DONE."
  SWRITE(UNIT_stdOut,fmt_sep)
  IF(lastoutputIter.NE.iter)THEN
    CALL Analyze(MIN(iter,MaxIter))
    CALL WriteState(U(0),MIN(iter,MaxIter))
  END IF
  CALL FinishLogging()
  CALL writeSFLoutfile(U(0),MIN(iter,MaxIter))


CONTAINS

  !=================================================================================================================================
  !> all screen and logfile tasks, can use all variables from subroutine above
  !!
  !=================================================================================================================================
  SUBROUTINE StartLogging()
  USE MODgvec_Globals,     ONLY: GETFREEUNIT
  USE MODgvec_Output_Vars, ONLY: ProjectName,outputLevel
  USE MODgvec_MHD3D_visu,  ONLY: checkAxis
  IMPLICIT NONE
  !---------------------------------------------------------------------------------------------------------------------------------
  CHARACTER(LEN=255)  :: fileString
  INTEGER             :: TimeArray(8),iLogDat
  REAL(wp)            :: AxisPos(2,2),W_MHD3D
  INTEGER,PARAMETER   :: nLogDat=20
  REAL(wp)            :: LogDat(1:nLogDat)
  !=================================================================================================================================
  IF(.NOT.MPIroot) RETURN
  __PERFON('log_output')
  W_MHD3D=U(0)%W_MHD3D
  CALL DATE_AND_TIME(values=TimeArray) ! get System time
  WRITE(UNIT_stdOut,'(A,E11.4,A)')'%%%%%%%%%%  START ITERATION, dt= ',dt, '  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
  WRITE(UNIT_stdOut,'(A,I4.2,"-",I2.2,"-",I2.2,1X,I2.2,":",I2.2,":",I2.2)') &
                 '%%% Sys date : ',timeArray(1:3),timeArray(5:7)
  WRITE(UNIT_stdOut,'(A,3E21.14)') &
          '%%% dU = |Force|= ',Fnorm(1:3)
  IF(MinimizerType.EQ.10) THEN
    WRITE(UNIT_stdOut,'(A,E11.4,A,3E11.4)') &
          '%%% accel.GD: tau= ',tau_bar,' |vel|= ',Vnorm(1:3)
  END IF

  WRITE(UNIT_stdOut,'(40(" -"))')
  !------------------------------------
  StartTimeArray=TimeArray !save first time stamp

  logUnit=GETFREEUNIT()
  WRITE(FileString,'("logMinimizer_",A,"_",I4.4,".csv")')TRIM(ProjectName),outputLevel
  OPEN(UNIT     = logUnit       ,&
     FILE     = TRIM(FileString) ,&
     STATUS   = 'REPLACE'   ,&
     ACCESS   = 'SEQUENTIAL' )
  !header
  iLogDat=0
  WRITE(logUnit,'(A)',ADVANCE="NO")'"#iterations","runtime(s)","min_dt","max_dt"'
  WRITE(logUnit,'(A)',ADVANCE="NO")',"W_MHD3D","min_dW","max_dW","sum_dW"'
  WRITE(logUnit,'(A)',ADVANCE="NO")',"normF_X1","normF_X2","normF_LA"'
  LogDat(ilogDat+1:iLogDat+11)=(/0.0_wp,0.0_wp,dt,dt,W_MHD3D,0.0_wp,0.0_wp,0.0_wp,Fnorm(1:3)/)
  iLogDat=11
  IF(MinimizerType.EQ.10) THEN
    WRITE(logUnit,'(A)',ADVANCE="NO")',"tau","normV_X1","normV_X2","normV_LA"'
    LogDat(ilogDat+1:iLogDat+4)=(/tau_bar,Vnorm(1:3)/)
    iLogDat=iLogDat+4
  END IF
  IF(doCheckDistance) THEN
    WRITE(logUnit,'(A)',ADVANCE="NO")',"max_Dist","avg_Dist"'
    LogDat(iLogDat+1:iLogDat+2)=(/0.0_wp,0.0_wp/)
    iLogDat=iLogDat+2
  END IF!doCheckDistance
  IF(doCheckAxis) THEN
    WRITE(logUnit,'(A)',ADVANCE="NO")',"X1_axis_0","X2_axis_0","X1_axis_1","X2_axis_1"'
    CALL CheckAxis(U(0),2,AxisPos)
    LogDat(iLogDat+1:iLogDat+4)=RESHAPE(AxisPos,(/4/))
    iLogDat=iLogDat+4
  END IF!doCheckAxis
  WRITE(logUnit,'(A)')' '
  !first data line
  WRITE(logUnit,'(*(e23.15,:,","))') logDat(1:iLogDat)
  __PERFOFF('log_output')
  END SUBROUTINE StartLogging

  !=================================================================================================================================
  !> all screen and logfile tasks, can use all variables from subroutine above
  !!
  !=================================================================================================================================
  SUBROUTINE Logging(quiet)
  USE MODgvec_MHD3D_visu, ONLY: checkDistance
  USE MODgvec_MHD3D_visu, ONLY: checkAxis
  IMPLICIT NONE
  LOGICAL, INTENT(IN) :: quiet !! True: no screen output
  !---------------------------------------------------------------------------------------------------------------------------------
  INTEGER             :: TimeArray(8),runtime_ms,iLogDat
  REAL(wp)            :: AxisPos(2,2),maxDist,avgDist,W_MHD3D
  INTEGER,PARAMETER   :: nLogDat=20
  REAL(wp)            :: LogDat(1:nLogDat)
  !=================================================================================================================================
  IF(.NOT.MPIroot) RETURN
  __PERFON('log_output')
  CALL DATE_AND_TIME(values=TimeArray) ! get System time
  W_MHD3D=U(0)%W_MHD3D
  IF(.NOT.quiet)THEN
    WRITE(UNIT_stdOut,'(80("%"))')
    WRITE(UNIT_stdOut,'(A,I4.2,"-",I2.2,"-",I2.2,1X,I2.2,":",I2.2,":",I2.2)') &
                      '%%% Sys date : ',timeArray(1:3),timeArray(5:7)
    WRITE(UNIT_stdOut,'(A,I8,A,2I8,A,E11.4,A,2E11.4,A,E21.14,A,3E12.4)') &
                      '%%% #ITERATIONS= ',iter,', #skippedIter (Jac/dW)= ',nSkip_Jac,nSkip_dW, &
              '\n%%% t_pseudo= ',t_pseudo,', min/max dt= ',min_dt_out,max_dt_out, &
              '\n%%% W_MHD3D= ',W_MHD3D,', min/max/sum deltaW= ' , min_dW_out,max_dW_out,sum_dW_out
    WRITE(UNIT_stdOut,'(A,3E21.14)') &
                '%%% dU = |Force|= ',Fnorm(1:3)
    !------------------------------------
  END IF!.NOT.quiet
  iLogDat=0
  runtime_ms=MAX(0,SUM((timeArray(5:8)-StartTimearray(5:8))*(/360000,6000,100,1/)))
  LogDat(ilogDat+1:iLogDat+11)=(/REAL(iter,wp),REAL(runtime_ms,wp)/100.0_wp, &
                                min_dt_out,max_dt_out,W_MHD3D,min_dW_out,max_dW_out,sum_dW_out, &
                                Fnorm(1:3)/)
  iLogDat=11
  IF(MinimizerType.EQ.10) THEN
    IF(.NOT.quiet)THEN
      WRITE(UNIT_stdOut,'(A,E11.4,A,3E11.4)') &
            '%%% accel.GD: tau= ',tau_bar,' |vel|= ',Vnorm(1:3)
    END IF!.NOT.quiet
    LogDat(ilogDat+1:iLogDat+4)=(/tau_bar,Vnorm(1:3)/)
    iLogDat=iLogDat+4
  END IF
  IF(doCheckDistance) THEN
    CALL CheckDistance(U(0),U(-2),maxDist,avgDist)
    CALL U(-2)%set_to(U(0))
    IF(.NOT.quiet)THEN
      WRITE(UNIT_stdOut,'(A,2E11.4)') &
      '               %%% Dist to last log (max/avg) : ',maxDist,avgDist
    END IF!.NOT.quiet
    LogDat(iLogDat+1:iLogDat+2)=(/maxDist,avgDist/)
    iLogDat=iLogDat+2
  END IF!doCheckDistance
  IF(doCheckAxis) THEN
    CALL CheckAxis(U(0),2,AxisPos)
    IF(.NOT.quiet)THEN
      WRITE(UNIT_stdOut,'(2(A,2E22.14))') &
        '%%% axis position (X1,X2,zeta=0     ): ',AxisPos(1:2,1), &
      '\n%%% axis position (X1,X2,zeta=pi/nfp): ',AxisPos(1:2,2)
    END IF!.NOT.quiet
    LogDat(iLogDat+1:iLogDat+4)=RESHAPE(AxisPos,(/4/))
    iLogDat=iLogDat+4
  END IF !doCheckAxis
  IF(.NOT.quiet)THEN
    WRITE(UNIT_stdOut,'(40(" -"))')
  END IF!.NOT.quiet
  WRITE(logUnit,'(*(e23.15,:,","))') logDat(1:iLogDat)
  __PERFOFF('log_output')
  END SUBROUTINE Logging

  !=================================================================================================================================
  !>
  !!
  !=================================================================================================================================
  SUBROUTINE FinishLogging()
  IMPLICIT NONE
  !---------------------------------------------------------------------------------------------------------------------------------
  CLOSE(logUnit)
  END SUBROUTINE FinishLogging

END SUBROUTINE MinimizeMHD3D_descent


!===================================================================================================================================
!> 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
  INTEGER :: i
!===================================================================================================================================
  CALL X1_base%free()
  CALL X2_base%free()
  CALL LA_base%free()

  DO i=-1,1
    CALL U(i)%free()
    CALL P(i)%free()
    CALL V(i)%free()
  END DO
  DO i=-1,0
    CALL F(i)%free()
  END DO
  CALL sgrid%free()

  SDEALLOCATE(U)
  SDEALLOCATE(P)
  SDEALLOCATE(V)
  SDEALLOCATE(F)
  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)
END SUBROUTINE FinalizeMHD3D

END MODULE MODgvec_MHD3D
