AddBoundaryPerturbation Subroutine

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