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