AddBoundaryPerturbation Subroutine

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