34 #include "sll_assert.h"
35 #include "sll_errors.h"
36 #include "sll_memory.h"
37 #include "sll_working_precision.h"
60 sll_f_cubic_spline_1d_get_x1_delta, &
71 sll_f_cubic_spline_2d_get_x1_min, &
72 sll_f_cubic_spline_2d_get_x1_max, &
73 sll_f_cubic_spline_2d_get_x1_delta, &
74 sll_f_cubic_spline_2d_get_x2_min, &
75 sll_f_cubic_spline_2d_get_x2_max, &
76 sll_f_cubic_spline_2d_get_x2_delta
82 sll_real64,
parameter ::
inv_6 = 1._f64/6._f64
89 sll_int32,
private :: n_points
90 sll_real64,
private :: delta
91 sll_real64,
private :: rdelta
92 sll_real64,
private :: xmin
93 sll_real64,
private :: xmax
94 sll_int32,
private :: bc_type
97 sll_real64,
dimension(:),
pointer,
private :: d => null()
99 sll_real64,
dimension(:),
pointer,
private :: coeffs => null()
101 sll_real64,
private :: slope_l
103 sll_real64,
private :: slope_r
105 logical,
private :: compute_slope_l
107 logical,
private :: compute_slope_r
111 logical,
private :: use_fast_algorithm
113 sll_real64,
dimension(:),
pointer,
private :: a => null()
115 sll_real64,
dimension(:),
pointer,
private :: cts => null()
117 sll_int32,
dimension(:),
pointer,
private :: ipiv => null()
119 sll_real64,
dimension(:),
pointer,
private :: f_aux => null()
128 sll_int32,
private :: num_pts_x1
129 sll_int32,
private :: num_pts_x2
130 sll_real64,
private :: x1_delta
131 sll_real64,
private :: x1_rdelta
132 sll_real64,
private :: x2_delta
133 sll_real64,
private :: x2_rdelta
134 sll_real64,
private :: x1_min
135 sll_real64,
private :: x1_max
136 sll_real64,
private :: x2_min
137 sll_real64,
private :: x2_max
138 sll_int32,
private :: x1_bc_type
139 sll_int32,
private :: x2_bc_type
141 sll_real64,
pointer,
private ::
data(:, :) => null()
142 sll_real64,
pointer,
private :: d1(:) => null()
145 sll_real64,
pointer,
private :: d2(:) => null()
146 sll_real64,
pointer,
private :: coeffs(:, :) => null()
147 sll_real64,
pointer,
private :: x1_min_slopes(:) => null()
148 sll_real64,
pointer,
private :: x1_max_slopes(:) => null()
149 sll_real64,
pointer,
private :: x2_min_slopes(:) => null()
150 sll_real64,
pointer,
private :: x2_max_slopes(:) => null()
151 sll_real64,
pointer,
private :: x1_min_slopes_coeffs(:) => null()
152 sll_real64,
pointer,
private :: x1_max_slopes_coeffs(:) => null()
153 sll_real64,
pointer,
private :: x2_min_slopes_coeffs(:) => null()
154 sll_real64,
pointer,
private :: x2_max_slopes_coeffs(:) => null()
155 logical,
private :: compute_slopes_x1_min
156 logical,
private :: compute_slopes_x1_max
157 logical,
private :: compute_slopes_x2_min
158 logical,
private :: compute_slopes_x2_max
165 #ifndef DOXYGEN_SHOULD_SKIP_THIS
176 #define FORWARD_FD_5PT( f, r_delta ) \
177 r_delta*(-(25.0_f64/12.0_f64)*f(1) + 4.0_f64*f(2) - 3.0_f64*f(3) + (4.0_f64/3.0_f64)*f(4) - 0.25_f64*f(5))
180 #define BACKWARD_FD_5PT( f, r_delta, np ) \
181 r_delta*(0.25_f64*f(np - 4) - (4.0_f64/3.0_f64)*f(np - 3) + 3.0_f64*f(np - 2) - 4.0_f64*f(np - 1) + (25.0_f64/12.0_f64)*f(np))
184 #define FORWARD_FD_3PT( f, r_delta ) \
185 r_delta*(-1.5_f64*f(1) + 2.0_f64*f(2) - 0.5_f64*f(3))
188 #define BACKWARD_FD_3PT( f, r_delta, np ) \
189 r_delta*(0.5_f64*f(np - 2) - 2.0_f64*f(np - 1) + 1.5_f64*f(np))
191 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
197 #ifndef DOXYGEN_SHOULD_SKIP_THIS
199 #define MAKE_GET_SLOT_FUNCTION( fname, datatype, slot, ret_type ) \
200 function fname(spline_obj)
result(val); \
201 type(datatype) :: spline_obj; \
203 val = spline_obj%slot; \
214 make_get_slot_function(sll_f_cubic_spline_2d_get_x1_delta,
sll_t_cubic_spline_2d, x1_delta, sll_real64)
215 make_get_slot_function(sll_f_cubic_spline_2d_get_x2_delta,
sll_t_cubic_spline_2d, x2_delta, sll_real64)
231 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
251 sll_int32,
intent(in) :: num_points
252 sll_real64,
intent(in) :: xmin
253 sll_real64,
intent(in) :: xmax
254 sll_int32,
intent(in) :: bc_type
255 sll_real64,
intent(in),
optional :: sl
256 sll_real64,
intent(in),
optional :: sr
257 logical,
intent(in),
optional :: fast_algorithm
261 self%n_points = num_points
264 self%delta = (xmax - xmin)/real((num_points - 1), f64)
265 self%rdelta = 1.0_f64/self%delta
266 self%bc_type = bc_type
267 if (num_points .lt. num_terms)
then
268 self%use_fast_algorithm = .false.
270 self%use_fast_algorithm = .true.
271 if (
present(fast_algorithm))
then
272 self%use_fast_algorithm = fast_algorithm
275 if (xmin .gt. xmax)
then
276 print *,
'ERROR, sll_s_cubic_spline_1d_init: xmin is greater than xmax, ', &
277 'this would cause all sorts of errors.'
282 select case (bc_type)
283 case (sll_p_periodic)
284 if (
present(sl) .or.
present(sr))
then
289 sll_warning(
'self',
'values of sl and sr are not taken into account')
293 self%slope_L = 0.0_f64
294 self%slope_R = 0.0_f64
296 if (self%use_fast_algorithm .eqv. .false.)
then
297 sll_allocate(self%a(3*(num_points - 1)), ierr)
298 sll_allocate(self%cts(7*(num_points - 1)), ierr)
299 sll_allocate(self%ipiv(num_points - 1), ierr)
303 self%a(1) = 1.0_f64/6.0_f64
304 self%a(2) = 4.0_f64/6.0_f64
305 self%a(3) = 1.0_f64/6.0_f64
306 self%a(3*num_points - 5) = 1.0_f64/6.0_f64
307 self%a(3*num_points - 4) = 4.0_f64/6.0_f64
308 self%a(3*num_points - 3) = 1.0_f64/6.0_f64
309 do i = 1, num_points - 3
310 self%a(3*i + 1) = 1.0_f64/6.0_f64
311 self%a(3*i + 2) = 4.0_f64/6.0_f64
312 self%a(3*i + 3) = 1.0_f64/6.0_f64
326 if (
present(sl))
then
328 self%compute_slope_L = .false.
330 self%slope_L = 0.0_f64
331 self%compute_slope_L = .true.
333 if (
present(sr))
then
335 self%compute_slope_R = .false.
337 self%slope_R = 0.0_f64
338 self%compute_slope_R = .true.
340 if (self%use_fast_algorithm .eqv. .false.)
then
341 sll_allocate(self%a(3*(num_points + 2)), ierr)
342 sll_allocate(self%cts(7*(num_points + 2)), ierr)
343 sll_allocate(self%ipiv(num_points + 2), ierr)
344 sll_allocate(self%f_aux(num_points + 2), ierr)
349 self%a(2) = 4.0_f64/6.0_f64
350 self%a(3) = 2.0_f64/6.0_f64
351 self%a(3*(num_points + 2) - 2) = 2.0_f64/6.0_f64
352 self%a(3*(num_points + 2) - 1) = 4.0_f64/6.0_f64
353 self%a(3*(num_points + 2)) = 0.0_f64
355 self%a(3*i + 1) = 1.0_f64/6.0_f64
356 self%a(3*i + 2) = 4.0_f64/6.0_f64
357 self%a(3*i + 3) = 1.0_f64/6.0_f64
371 print *,
'ERROR: sll_s_cubic_spline_1d_init(): not recognized boundary condition'
374 if (self%use_fast_algorithm .eqv. .true.)
then
375 sll_allocate(self%d(num_points), ierr)
383 sll_allocate(self%coeffs(0:num_points + 2), ierr)
499 sll_real64,
dimension(:),
intent(in) :: f
508 bc_type = spline%bc_type;
509 select case (bc_type)
510 case (sll_p_periodic)
515 print *,
'ERROR: compute_spline_1D(): not recognized boundary condition'
532 sll_real64,
dimension(:),
pointer :: f
533 sll_int32,
intent(in) :: num_pts
534 sll_real64,
dimension(:),
pointer :: d
535 sll_real64,
dimension(:),
pointer :: coeffs
536 sll_real64,
parameter :: a = sqrt((2.0_f64 + sqrt(3.0_f64))/6.0_f64)
537 sll_real64,
parameter :: r_a = 1.0_f64/a
538 sll_real64,
parameter :: b = sqrt((2.0_f64 - sqrt(3.0_f64))/6.0_f64)
539 sll_real64,
parameter :: b_a = b/a
540 sll_real64 :: coeff_tmp
544 sll_assert(
size(f) .ge. num_pts - 1)
545 sll_assert(
size(d) .ge. num_pts)
546 sll_assert(
size(coeffs) .ge. num_pts)
547 sll_assert(num_pts .gt. num_terms)
553 do i = 0, num_terms - 1
555 coeff_tmp = coeff_tmp*(-b_a)
556 d1 = d1 + coeff_tmp*f(np - 1 - i)
561 d(i) = r_a*(f(i) - b*d(i - 1))
567 coeff_tmp = coeff_tmp*(-b_a)
568 d1 = d1 + coeff_tmp*d(i)
575 coeffs(i + 1) = r_a*(d(i) - b*coeffs(i + 2))
577 coeffs(1) = coeffs(np)
578 coeffs(np + 1) = coeffs(2)
579 coeffs(np + 2) = coeffs(3)
580 coeffs(np + 3) = coeffs(4)
592 sll_real64,
dimension(:),
pointer :: f
593 sll_int32,
intent(in) :: num_pts
594 sll_real64,
dimension(:),
pointer :: d
595 sll_real64,
intent(in) :: slope_l
596 sll_real64,
intent(in) :: slope_r
597 sll_real64,
intent(in) :: delta
598 sll_real64,
dimension(:),
pointer :: coeffs
601 sll_real64,
parameter :: a = sqrt((2.0_f64 + sqrt(3.0_f64))/6.0_f64)
602 sll_real64,
parameter :: r_a = 1.0_f64/a
603 sll_real64,
parameter :: b = sqrt((2.0_f64 - sqrt(3.0_f64))/6.0_f64)
604 sll_real64,
parameter :: b_a = b/a
605 sll_real64,
parameter :: ralpha = sqrt(6.0_f64/sqrt(3.0_f64))
606 sll_real64 :: coeff_tmp
623 fnp = f(np) - delta*slope_r/3.0_f64
634 coeff_tmp = coeff_tmp*(-b_a)
635 d1 = d1 + coeff_tmp*(f(i) - 2.0_f64*slope_l*delta*real(i - 1, f64))
640 d(i) = r_a*(f(i) - b*d(i - 1))
642 d(np) = ralpha*(0.5_f64*fnp - b*d(np - 1))
645 coeffs(np + 1) = ralpha*d(np)
647 coeffs(i + 1) = r_a*(d(i) - b*coeffs(i + 2))
649 coeffs(1) = coeffs(3) - 2.0_f64*delta*slope_l
650 coeffs(np + 2) = coeffs(np) + 2.0_f64*delta*slope_r
651 coeffs(np + 3) = 0.0_f64
655 sll_real64,
dimension(:),
intent(in),
target :: f
657 sll_real64,
dimension(:),
pointer :: coeffs
659 sll_real64,
dimension(:),
pointer :: d
660 sll_real64,
dimension(:),
pointer :: fp
662 if (.not. (
size(f) .ge. spline%n_points - 1))
then
664 print *,
'ERROR: compute_spline_1D_periodic(): '
665 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
666 spline%n_points - 1,
' . Passed size: ',
size(f)
669 if (spline%use_fast_algorithm .eqv. .false.)
then
670 np = spline%n_points - 1
678 spline%coeffs(0) = spline%coeffs(np)
679 spline%coeffs(np + 1) = spline%coeffs(1)
680 spline%coeffs(np + 2) = spline%coeffs(2)
681 spline%coeffs(np + 3) = spline%coeffs(3)
686 coeffs => spline%coeffs(0:np + 2)
693 sll_real64,
dimension(:),
intent(in),
target :: f
695 sll_real64,
dimension(:),
pointer :: coeffs
697 sll_real64,
dimension(:),
pointer :: fp
698 sll_real64,
dimension(:),
pointer :: d
699 sll_real64 :: slope_l
700 sll_real64 :: slope_r
702 sll_real64 :: r_delta
705 if (.not. (
size(f) .ge. spline%n_points))
then
707 print *,
'ERROR: compute_spline_1D_hermite(): '
708 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
709 spline%n_points,
' . Passed size: ',
size(f)
715 coeffs => spline%coeffs(0:)
717 r_delta = 1.0_f64/delta
721 if (spline%compute_slope_L .eqv. .true.)
then
722 slope_l = forward_fd_5pt(f, r_delta)
724 slope_l = spline%slope_L
726 if (spline%compute_slope_R .eqv. .true.)
then
727 slope_r = backward_fd_5pt(f, r_delta, np)
729 slope_r = spline%slope_R
732 if (.not. spline%use_fast_algorithm)
then
734 spline%f_aux(2:np + 1) = fp(1:np)
736 spline%f_aux(1) = fp(1) + (1.0_f64/3.0_f64)*delta*slope_l
737 spline%f_aux(np + 2) = fp(np) - (1.0_f64/3.0_f64)*delta*slope_r
746 fp, np, d, slope_l, slope_r, delta, coeffs)
831 sll_real64,
intent(in) :: x
832 sll_real64,
intent(in) :: xmin
833 sll_real64,
dimension(:),
pointer :: coeffs
834 sll_real64,
intent(in) :: rh
850 dx = t0 - real(cell - 1, f64)
853 ci = coeffs(cell + 1)
854 cip1 = coeffs(cell + 2)
855 cip2 = coeffs(cell + 3)
860 t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
861 t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
876 intrinsic ::
associated, int, real
877 sll_real64,
intent(in) :: x
879 sll_real64,
dimension(:),
pointer :: coeffs
885 sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
888 coeffs => spline%coeffs(0:spline%n_points + 2)
904 intrinsic ::
associated, int, real
905 sll_int32,
intent(in) :: n
906 sll_real64,
dimension(1:n),
intent(in) :: a_in
907 sll_real64,
dimension(1:n),
intent(out) :: a_out
909 sll_real64,
dimension(:),
pointer :: coeffs
923 sll_int32 :: num_cells
927 num_cells = spline%n_points - 1
929 coeffs => spline%coeffs
933 if (.not. ((x .ge. spline%xmin) .and. (x .le. spline%xmax)))
then
934 print *,
'splines', x, spline%xmin, spline%xmax
937 sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
938 t0 = (x - spline%xmin)*rh
940 dx = t0 - real(cell - 1, f64)
943 cim1 = coeffs(cell - 1)
945 cip1 = coeffs(cell + 1)
946 cip2 = coeffs(cell + 2)
949 t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
950 t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
951 a_out(i) = (1.0_f64/6.0_f64)*(t2 + t4)
972 intrinsic ::
associated, int, real
973 sll_int32,
intent(in) :: n
974 sll_real64,
dimension(:),
pointer :: ptr_in
975 sll_real64,
dimension(:),
pointer :: ptr_out
977 sll_real64,
dimension(:),
pointer :: coeffs
991 sll_int32 :: num_cells
994 sll_assert(
associated(ptr_in))
995 sll_assert(
associated(ptr_out))
997 num_cells = spline%n_points - 1
999 coeffs => spline%coeffs
1004 sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
1005 t0 = (x - spline%xmin)*rh
1007 dx = t0 - real(cell - 1, f64)
1010 cim1 = coeffs(cell - 1)
1012 cip1 = coeffs(cell + 1)
1013 cip2 = coeffs(cell + 2)
1016 t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
1017 t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
1018 ptr_out(i) = (1.0_f64/6.0_f64)*(t2 + t4)
1027 intrinsic :: int, real
1028 sll_real64,
intent(in) :: x
1029 sll_real64,
intent(in) :: xmin
1030 sll_real64,
intent(in) :: rh
1031 sll_real64,
dimension(:),
pointer :: coeffs
1046 dx = t0 - real(cell - 1, f64)
1049 ci = coeffs(cell + 1)
1050 cip1 = coeffs(cell + 2)
1051 cip2 = coeffs(cell + 3)
1052 t1 = 2.0_f64*(cim1 - 2.0_f64*ci + cip1)
1053 t2 = -cim1 + 3.0_f64*(ci - cip1) + cip2
1067 intrinsic ::
associated
1068 sll_real64,
intent(in) :: x
1069 sll_real64,
dimension(:),
pointer :: coeffs
1074 sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
1075 coeffs => spline%coeffs(0:spline%n_points + 2)
1100 intrinsic ::
associated
1101 sll_real64,
dimension(:),
intent(in) :: array_in
1102 sll_int32,
intent(in) :: num_pts
1103 sll_real64,
dimension(:),
intent(out) :: array_out
1105 sll_real64,
dimension(:),
pointer :: coeffs
1108 sll_assert(num_pts .le.
size(array_in))
1110 coeffs => spline%coeffs(0:spline%n_points + 2)
1112 sll_assert((array_in(i) .ge. spline%xmin) .and. (array_in(i) .le. spline%xmax))
1114 array_in(i), spline%xmin, spline%rdelta, coeffs)
1134 intrinsic ::
associated
1135 sll_real64,
dimension(:),
pointer :: ptr_in
1136 sll_int32,
intent(in) :: num_pts
1137 sll_real64,
dimension(:),
pointer :: ptr_out
1139 sll_real64,
dimension(:),
pointer :: coeffs
1142 sll_assert(num_pts .le.
size(ptr_in))
1143 sll_assert(
associated(ptr_in))
1144 sll_assert(
associated(ptr_out))
1145 coeffs => spline%coeffs(0:spline%n_points + 2)
1148 sll_assert((ptr_in(i) .ge. spline%xmin) .and. (ptr_in(i) .le. spline%xmax))
1150 ptr_in(i), spline%xmin, spline%rdelta, coeffs)
1159 if (spline%use_fast_algorithm .eqv. .true.)
then
1160 sll_deallocate(spline%d, ierr)
1162 sll_deallocate(spline%coeffs, ierr)
1164 if (spline%use_fast_algorithm .eqv. .false.)
then
1165 sll_deallocate(spline%a, ierr)
1166 sll_deallocate(spline%cts, ierr)
1167 sll_deallocate(spline%ipiv, ierr)
1168 if (spline%bc_type == sll_p_hermite)
then
1169 sll_deallocate(spline%f_aux, ierr)
1238 const_slope_x1_min, &
1239 const_slope_x1_max, &
1240 const_slope_x2_min, &
1241 const_slope_x2_max, &
1245 x2_max_slopes)
result(this)
1248 sll_int32,
intent(in) :: num_pts_x1
1249 sll_int32,
intent(in) :: num_pts_x2
1250 sll_real64,
intent(in) :: x1_min
1251 sll_real64,
intent(in) :: x1_max
1252 sll_real64,
intent(in) :: x2_min
1253 sll_real64,
intent(in) :: x2_max
1254 sll_int32,
intent(in) :: x1_bc_type
1255 sll_int32,
intent(in) :: x2_bc_type
1256 sll_real64,
intent(in),
optional :: const_slope_x1_min
1257 sll_real64,
intent(in),
optional :: const_slope_x1_max
1258 sll_real64,
intent(in),
optional :: const_slope_x2_min
1259 sll_real64,
intent(in),
optional :: const_slope_x2_max
1260 sll_real64,
intent(in),
dimension(:),
optional :: x1_min_slopes
1261 sll_real64,
intent(in),
dimension(:),
optional :: x1_max_slopes
1262 sll_real64,
intent(in),
dimension(:),
optional :: x2_min_slopes
1263 sll_real64,
intent(in),
dimension(:),
optional :: x2_max_slopes
1276 const_slope_x1_min, &
1277 const_slope_x1_max, &
1278 const_slope_x2_min, &
1279 const_slope_x2_max, &
1297 const_slope_x1_min, &
1298 const_slope_x1_max, &
1299 const_slope_x2_min, &
1300 const_slope_x2_max, &
1307 sll_int32,
intent(in) :: num_pts_x1
1308 sll_int32,
intent(in) :: num_pts_x2
1309 sll_real64,
intent(in) :: x1_min
1310 sll_real64,
intent(in) :: x1_max
1311 sll_real64,
intent(in) :: x2_min
1312 sll_real64,
intent(in) :: x2_max
1313 sll_int32,
intent(in) :: x1_bc_type
1314 sll_int32,
intent(in) :: x2_bc_type
1315 sll_real64,
intent(in),
optional :: const_slope_x1_min
1316 sll_real64,
intent(in),
optional :: const_slope_x1_max
1317 sll_real64,
intent(in),
optional :: const_slope_x2_min
1318 sll_real64,
intent(in),
optional :: const_slope_x2_max
1319 sll_real64,
intent(in),
dimension(:),
optional :: x1_min_slopes
1320 sll_real64,
intent(in),
dimension(:),
optional :: x1_max_slopes
1321 sll_real64,
intent(in),
dimension(:),
optional :: x2_min_slopes
1322 sll_real64,
intent(in),
dimension(:),
optional :: x2_max_slopes
1324 sll_int32 :: bc_selector
1327 this%num_pts_x1 = num_pts_x1
1328 this%num_pts_x2 = num_pts_x2
1329 this%x1_min = x1_min
1330 this%x1_max = x1_max
1331 this%x2_min = x2_min
1332 this%x2_max = x2_max
1333 this%x1_delta = (x1_max - x1_min)/real((num_pts_x1 - 1), f64)
1334 this%x2_delta = (x2_max - x2_min)/real((num_pts_x2 - 1), f64)
1335 this%x1_rdelta = 1.0_f64/this%x1_delta
1336 this%x2_rdelta = 1.0_f64/this%x2_delta
1337 this%x1_bc_type = x1_bc_type
1338 this%x2_bc_type = x2_bc_type
1340 if ((num_pts_x1 .le. num_terms) .or. (num_pts_x2 .le. num_terms))
then
1341 sll_error(
"sll_s_cubic_spline_2d_init",
" Because of the algorithm used, this function is meant to be used with arrays that are at least of size = 28")
1343 if ((x1_min .gt. x1_max) .or. (x2_min .gt. x2_max))
then
1344 sll_error(
"sll_s_cubic_spline_2d_init",
" one of the xmin is greater than the corresponding xmax, this would cause all sorts of errors")
1350 if (
present(x1_min_slopes))
then
1351 sll_assert(
size(x1_min_slopes) .ge. num_pts_x2)
1353 if (
present(x1_max_slopes))
then
1354 sll_assert(
size(x1_max_slopes) .ge. num_pts_x2)
1356 if (
present(x2_min_slopes))
then
1357 sll_assert(
size(x2_min_slopes) .ge. num_pts_x1)
1359 if (
present(x2_max_slopes))
then
1360 sll_assert(
size(x2_max_slopes) .ge. num_pts_x1)
1363 sll_allocate(this%d1(num_pts_x1), ierr)
1364 sll_allocate(this%d2(num_pts_x2), ierr)
1371 if (x1_bc_type .eq. sll_p_periodic)
then
1372 bc_selector = bc_selector + 1
1374 if (x1_bc_type .eq. sll_p_hermite)
then
1375 bc_selector = bc_selector + 2
1377 if (x2_bc_type .eq. sll_p_periodic)
then
1378 bc_selector = bc_selector + 4
1380 if (x2_bc_type .eq. sll_p_hermite)
then
1381 bc_selector = bc_selector + 8
1384 select case (bc_selector)
1388 present(x1_min_slopes) .or.
present(x1_max_slopes) .or. &
1389 present(x2_min_slopes) .or.
present(x2_max_slopes) .or. &
1390 present(const_slope_x1_min) .or.
present(const_slope_x1_max) .or. &
1391 present(const_slope_x2_min) .or.
present(const_slope_x2_max))
then
1393 sll_warning(
'sll_s_cubic_spline_2d_init',
'values of slopes are not taken into account as we are in periodic-periodic')
1398 present(x2_min_slopes) .or.
present(x2_max_slopes) .or. &
1399 present(const_slope_x2_min) .or.
present(const_slope_x2_max))
then
1400 sll_warning(
'sll_s_cubic_spline_2d_init',
'values of slopes in x2 are not taken into account as we are periodic in x2')
1403 if (
present(const_slope_x1_min) .and.
present(x1_min_slopes))
then
1404 sll_error(
"sll_s_cubic_spline_2d_init",
" hermite-periodic-case, it is not allowed to specify simultaneously a constant value for the slopes at x1_min and an array-specified set of slopes")
1406 if (
present(const_slope_x1_max) .and.
present(x1_max_slopes))
then
1407 sll_error(
"sll_s_cubic_spline_2d_init",
"hermite-periodic-case, it is not allowed to specify simultaneously a constant value for the slopes at x1_max and an array-specified set of slopes")
1411 sll_allocate(this%x1_min_slopes(num_pts_x2), ierr)
1412 sll_allocate(this%x1_max_slopes(num_pts_x2), ierr)
1417 sll_allocate(this%x1_min_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1418 sll_allocate(this%x1_max_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1428 #ifndef DOXYGEN_SHOULD_SKIP_THIS
1430 #define FILL_SLOPES(const_opt, input_opt, numpts, output, slopes) \
1431 if (
present(input_opt)) then; \
1432 this%output(1:numpts) = input_opt(1:numpts); \
1433 this%slopes = .false.; \
1434 else if (
present(const_opt)) then; \
1435 this%output(1:numpts) = const_opt; \
1436 this%slopes = .false.; \
1438 this%slopes = .true.; \
1441 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
1444 fill_slopes(const_slope_x1_min, x1_min_slopes, num_pts_x2, x1_min_slopes, compute_slopes_x1_min)
1447 fill_slopes(const_slope_x1_max, x1_max_slopes, num_pts_x2, x1_max_slopes, compute_slopes_x1_max)
1452 present(x1_min_slopes) .or.
present(x1_max_slopes) .or. &
1453 present(const_slope_x1_min) .or.
present(const_slope_x1_max))
then
1454 sll_error(
"sll_s_cubic_spline_2d_init",
"periodic-hermite case, it is not allowed to specify the end slopes in the case of periodic boundary conditions.")
1455 sll_warning(
'sll_s_cubic_spline_2d_init',
'values of slopes in x1 are not taken into account as we are periodic in x1')
1458 if (
present(const_slope_x2_min) .and.
present(x2_min_slopes))
then
1459 sll_error(
"sll_s_cubic_spline_2d_init",
"periodic-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_min and an array-specified set of slopes")
1461 if (
present(const_slope_x2_max) .and.
present(x2_max_slopes))
then
1462 sll_error(
"sll_s_cubic_spline_2d_init",
"periodic-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_max and an array-specified set of slopes")
1466 sll_allocate(this%x2_min_slopes(num_pts_x1), ierr)
1467 sll_allocate(this%x2_max_slopes(num_pts_x1), ierr)
1471 fill_slopes(const_slope_x2_min, x2_min_slopes, num_pts_x1, x2_min_slopes, compute_slopes_x2_min)
1474 fill_slopes(const_slope_x2_max, x2_max_slopes, num_pts_x1, x2_max_slopes, compute_slopes_x2_max)
1478 if (
present(const_slope_x1_min) .and.
present(x1_min_slopes))
then
1479 sll_error(
"sll_s_cubic_spline_2d_init",
"hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x1_min and an array-specified set of slopes")
1481 if (
present(const_slope_x1_max) .and.
present(x1_max_slopes))
then
1482 sll_error(
"sll_s_cubic_spline_2d_init",
"hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_max and an array-specified set of slopes")
1484 if (
present(const_slope_x2_min) .and.
present(x2_min_slopes))
then
1485 sll_error(
"sll_s_cubic_spline_2d_init",
"hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_min and an array-specified set of slopes")
1487 if (
present(const_slope_x2_max) .and.
present(x2_max_slopes))
then
1488 sll_error(
"sll_s_cubic_spline_2d_init",
"hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_max and an array-specified set of slopes")
1492 sll_allocate(this%x1_min_slopes(num_pts_x2), ierr)
1493 sll_allocate(this%x1_max_slopes(num_pts_x2), ierr)
1494 sll_allocate(this%x2_min_slopes(num_pts_x1), ierr)
1495 sll_allocate(this%x2_max_slopes(num_pts_x1), ierr)
1500 sll_allocate(this%x1_min_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1501 sll_allocate(this%x1_max_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1504 fill_slopes(const_slope_x1_min, x1_min_slopes, num_pts_x2, x1_min_slopes, compute_slopes_x1_min)
1507 fill_slopes(const_slope_x1_max, x1_max_slopes, num_pts_x2, x1_max_slopes, compute_slopes_x1_max)
1510 fill_slopes(const_slope_x2_min, x2_min_slopes, num_pts_x1, x2_min_slopes, compute_slopes_x2_min)
1513 fill_slopes(const_slope_x2_max, x2_max_slopes, num_pts_x1, x2_max_slopes, compute_slopes_x2_max)
1516 sll_error(
"sll_s_cubic_spline_2d_init",
"did not recognize given boundary conditions")
1522 sll_allocate(this%coeffs(0:num_pts_x1 + 2, 0:num_pts_x2 + 2), ierr)
1527 sll_real64,
dimension(:, :),
intent(in),
target ::
data
1529 sll_real64,
dimension(:),
pointer :: coeffs
1532 sll_real64,
dimension(:),
pointer :: d1
1533 sll_real64,
dimension(:),
pointer :: d2
1534 sll_real64,
dimension(:),
pointer :: datap
1543 if ((
size(
data, 1) .lt. spline%num_pts_x1))
then
1545 print *,
'ERROR: compute_spline_2D_prdc_prdc(): '
1546 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1547 spline%num_pts_x1,
' . Passed size: ',
size(
data, 1)
1550 if ((
size(
data, 2) .lt. spline%num_pts_x2))
then
1552 print *,
'ERROR: compute_spline_2D_prdc_prdc(): '
1553 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1554 spline%num_pts_x2,
' . Passed size: ',
size(
data, 2)
1557 npx1 = spline%num_pts_x1
1558 npx2 = spline%num_pts_x2
1564 datap =>
data(i, 1:npx2)
1565 coeffs => spline%coeffs(i, 0:npx2 + 2)
1574 datap => spline%coeffs(1:npx1, j)
1577 coeffs => spline%coeffs(0:npx1 + 2, j)
1583 sll_real64,
dimension(:, :),
intent(in),
target ::
data
1585 sll_real64,
dimension(:),
pointer :: coeffs
1588 sll_real64,
dimension(:),
pointer :: d1
1589 sll_real64,
dimension(:),
pointer :: d2
1590 sll_real64,
dimension(:),
pointer :: datap
1591 sll_real64,
dimension(:),
pointer :: coeffs_ptr1
1592 sll_real64,
dimension(:),
pointer :: coeffs_ptr2
1593 sll_real64 :: min_slope
1594 sll_real64 :: max_slope
1597 sll_real64 :: r_x1_delta
1605 if ((
size(
data, 1) .lt. spline%num_pts_x1))
then
1607 print *,
'ERROR: compute_spline_2D_prdc_prdc(): '
1608 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1609 spline%num_pts_x1,
' . Passed size: ',
size(
data, 1)
1612 if ((
size(
data, 2) .lt. spline%num_pts_x2))
then
1614 print *,
'ERROR: compute_spline_2D_prdc_prdc(): '
1615 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1616 spline%num_pts_x2,
' . Passed size: ',
size(
data, 2)
1619 npx1 = spline%num_pts_x1
1620 npx2 = spline%num_pts_x2
1627 datap =>
data(i, 1:npx2)
1628 coeffs => spline%coeffs(i, 0:npx2 + 2)
1641 r_x1_delta = 1.0_f64/spline%x1_delta
1642 if (spline%compute_slopes_x1_min .eqv. .true.)
then
1646 spline%x1_min_slopes(j) = &
1647 r_x1_delta*(-(25.0_f64/12.0_f64)*
data(1, j) + &
1648 4.0_f64*
data(2, j) - &
1649 3.0_f64*
data(3, j) + &
1650 (4.0_f64/3.0_f64)*
data(4, j) - &
1651 0.25_f64*
data(5, j))
1654 if (spline%compute_slopes_x1_max .eqv. .true.)
then
1657 spline%x1_max_slopes(j) = &
1658 r_x1_delta*(+0.25_f64*
data(npx1 - 4, j) - &
1659 (4.0_f64/3.0_f64)*
data(npx1 - 3, j) + &
1660 3.0_f64*
data(npx1 - 2, j) - &
1661 4.0_f64*
data(npx1 - 1, j) + &
1662 (25.0_f64/12.0_f64)*
data(npx1, j))
1669 coeffs_ptr1 => spline%x1_min_slopes_coeffs(0:npx2 + 2)
1671 spline%x1_min_slopes, &
1675 coeffs_ptr2 => spline%x1_max_slopes_coeffs(0:npx2 + 2)
1677 spline%x1_max_slopes, &
1683 datap => spline%coeffs(1:npx1, j)
1684 coeffs => spline%coeffs(0:npx1 + 2, j)
1685 min_slope = spline%x1_min_slopes_coeffs(j)
1686 max_slope = spline%x1_max_slopes_coeffs(j)
1700 sll_real64,
dimension(:, :),
intent(in),
target ::
data
1702 sll_real64,
dimension(:),
pointer :: coeffs
1705 sll_real64,
dimension(:),
pointer :: d1
1706 sll_real64,
dimension(:),
pointer :: d2
1707 sll_real64,
dimension(:),
pointer :: datap
1708 sll_real64 :: min_slope
1709 sll_real64 :: max_slope
1712 sll_real64 :: r_x2_delta
1720 if ((
size(
data, 1) .lt. spline%num_pts_x1))
then
1722 print *,
'ERROR: compute_spline_2D_prdc_prdc(): '
1723 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1724 spline%num_pts_x1,
' . Passed size: ',
size(
data, 1)
1727 if ((
size(
data, 2) .lt. spline%num_pts_x2))
then
1729 print *,
'ERROR: compute_spline_2D_prdc_prdc(): '
1730 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1731 spline%num_pts_x2,
' . Passed size: ',
size(
data, 2)
1734 npx1 = spline%num_pts_x1
1735 npx2 = spline%num_pts_x2
1738 r_x2_delta = 1.0_f64/spline%x2_delta
1742 if (spline%compute_slopes_x2_min .eqv. .true.)
then
1746 spline%x2_min_slopes(i) = &
1747 r_x2_delta*(-(25.0_f64/12.0_f64)*
data(i, 1) + &
1748 4.0_f64*
data(i, 2) - &
1749 3.0_f64*
data(i, 3) + &
1750 (4.0_f64/3.0_f64)*
data(i, 4) - &
1751 0.25_f64*
data(i, 5))
1754 if (spline%compute_slopes_x2_max .eqv. .true.)
then
1758 spline%x2_max_slopes(i) = &
1759 r_x2_delta*(+0.25_f64*
data(i, npx2 - 4) - &
1760 (4.0_f64/3.0_f64)*
data(i, npx2 - 3) + &
1761 3.0_f64*
data(i, npx2 - 2) - &
1762 4.0_f64*
data(i, npx2 - 1) + &
1763 (25.0_f64/12.0_f64)*
data(i, npx2))
1770 datap =>
data(i, 1:npx2)
1771 coeffs => spline%coeffs(i, 0:npx2 + 2)
1772 min_slope = spline%x2_min_slopes(i)
1773 max_slope = spline%x2_max_slopes(i)
1789 datap => spline%coeffs(1:npx1, j)
1790 coeffs => spline%coeffs(0:npx1 + 2, j)
1796 sll_real64,
dimension(:, :),
intent(in),
target ::
data
1798 sll_real64,
dimension(:),
pointer :: coeffs
1801 sll_real64,
dimension(:),
pointer :: d1
1802 sll_real64,
dimension(:),
pointer :: d2
1803 sll_real64,
dimension(:),
pointer :: datap
1804 sll_real64 :: min_slope
1805 sll_real64 :: max_slope
1808 sll_real64 :: r_x1_delta
1809 sll_real64 :: r_x2_delta
1817 if ((
size(
data, 1) .lt. spline%num_pts_x1))
then
1819 print *,
'ERROR: compute_spline_2D_hrmt_hrmt(): '
1820 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1821 spline%num_pts_x1,
' . Passed size: ',
size(
data, 1)
1824 if ((
size(
data, 2) .lt. spline%num_pts_x2))
then
1826 print *,
'ERROR: compute_spline_2D_hrmt_hrmt(): '
1827 write (*,
'(a, i8, a, i8)')
'spline object needs data of size >= ', &
1828 spline%num_pts_x2,
' . Passed size: ',
size(
data, 2)
1831 npx1 = spline%num_pts_x1
1832 npx2 = spline%num_pts_x2
1835 r_x1_delta = 1.0_f64/spline%x1_delta
1836 r_x2_delta = 1.0_f64/spline%x2_delta
1839 if (spline%compute_slopes_x1_min .eqv. .true.)
then
1843 spline%x1_min_slopes(j) = &
1844 r_x1_delta*(-(25.0_f64/12.0_f64)*
data(1, j) + &
1845 4.0_f64*
data(2, j) - &
1846 3.0_f64*
data(3, j) + &
1847 (4.0_f64/3.0_f64)*
data(4, j) - &
1848 0.25_f64*
data(5, j))
1851 if (spline%compute_slopes_x1_max .eqv. .true.)
then
1855 spline%x1_max_slopes(j) = &
1856 r_x1_delta*(+0.25_f64*
data(npx1 - 4, j) - &
1857 (4.0_f64/3.0_f64)*
data(npx1 - 3, j) + &
1858 3.0_f64*
data(npx1 - 2, j) - &
1859 4.0_f64*
data(npx1 - 1, j) + &
1860 (25.0_f64/12.0_f64)*
data(npx1, j))
1864 if (spline%compute_slopes_x2_min .eqv. .true.)
then
1868 spline%x2_min_slopes(i) = &
1869 r_x2_delta*(-(25.0_f64/12.0_f64)*
data(i, 1) + &
1870 4.0_f64*
data(i, 2) - &
1871 3.0_f64*
data(i, 3) + &
1872 (4.0_f64/3.0_f64)*
data(i, 4) - &
1873 0.25_f64*
data(i, 5))
1876 if (spline%compute_slopes_x2_max .eqv. .true.)
then
1880 spline%x2_max_slopes(i) = &
1881 r_x2_delta*(+0.25_f64*
data(i, npx2 - 4) - &
1882 (4.0_f64/3.0_f64)*
data(i, npx2 - 3) + &
1883 3.0_f64*
data(i, npx2 - 2) - &
1884 4.0_f64*
data(i, npx2 - 1) + &
1885 (25.0_f64/12.0_f64)*
data(i, npx2))
1892 datap =>
data(i, 1:npx2)
1893 coeffs => spline%coeffs(i, 0:npx2 + 2)
1894 min_slope = spline%x2_min_slopes(i)
1895 max_slope = spline%x2_max_slopes(i)
1916 datap => spline%coeffs(1:npx1, 0)
1917 coeffs => spline%coeffs(0:npx1 + 2, 0)
1918 min_slope = spline%x1_min_slopes(2)
1919 max_slope = spline%x1_max_slopes(2)
1929 datap => spline%coeffs(1:npx1, npx2 + 1)
1930 coeffs => spline%coeffs(0:npx1 + 2, npx2 + 1)
1931 min_slope = spline%x1_min_slopes(npx2 - 1)
1932 max_slope = spline%x1_max_slopes(npx2 - 1)
1943 datap => spline%coeffs(1:npx1, j)
1944 coeffs => spline%coeffs(0:npx1 + 2, j)
1945 min_slope = spline%x1_min_slopes(j)
1946 max_slope = spline%x1_max_slopes(j)
1965 sll_real64,
dimension(:, :),
intent(in),
target ::
data
1969 sll_int32 :: bc_selector
1978 bc1 = spline%x1_bc_type
1979 bc2 = spline%x2_bc_type
1989 if (bc1 .eq. sll_p_periodic)
then
1990 bc_selector = bc_selector + 1
1993 if (bc1 .eq. sll_p_hermite)
then
1994 bc_selector = bc_selector + 2
1997 if (bc2 .eq. sll_p_periodic)
then
1998 bc_selector = bc_selector + 4
2001 if (bc2 .eq. sll_p_hermite)
then
2002 bc_selector = bc_selector + 8
2005 select case (bc_selector)
2019 print *,
'ERROR: compute_spline_2D(): ', &
2020 'did not recognize given boundary condition combination.'
2037 intrinsic :: real, int
2038 sll_real64,
dimension(1:, 1:),
intent(in) :: x1
2039 sll_real64,
dimension(1:, 1:),
intent(in) :: x2
2041 sll_real64,
dimension(:, :),
intent(out) :: a_out
2044 sll_real64 :: x1_min
2045 sll_real64 :: x2_min
2064 sll_real64 :: svalx1, svalx2, svalx3, svalx4
2065 sll_real64 :: svaly1, svaly2, svaly3, svaly4
2066 sll_real64 :: ax1, ax2, ax3, ay1, ay2, ay3
2068 sll_real64 :: t1, t2
2069 sll_int32 :: ipm1, ip, ipp1, ipp2
2070 sll_int32 :: jpm1, jp, jpp1, jpp2
2082 if ((
size(x1, 1) .ne. spline%num_pts_x1) .or. (
size(x1, 2) .ne. spline%num_pts_x2))
then
2084 print *,
'ERROR: sll_s_cubic_spline_2d_deposit_value(): '
2085 write (*,
'(a, i8, i8, a, i8, i8, a)')
'array of feets of characteristics needs data of size = (', &
2086 spline%num_pts_x1, spline%num_pts_x2,
') . Passed size: (',
size(x1, 1),
size(x1, 2),
')'
2089 if ((
size(x2, 1) .ne. spline%num_pts_x1) .or. (
size(x2, 2) .ne. spline%num_pts_x2))
then
2091 print *,
'ERROR: sll_s_cubic_spline_2d_deposit_value(): '
2092 write (*,
'(a, i8, i8, a, i8, i8, a)')
'array of feets of characteristics needs data of size = (', &
2093 spline%num_pts_x1, spline%num_pts_x2,
') . Passed size: (',
size(x2, 1),
size(x2, 2),
')'
2096 if ((
size(a_out, 1) .ne. spline%num_pts_x1) .or. (
size(a_out, 2) .ne. spline%num_pts_x2))
then
2098 print *,
'ERROR: sll_s_cubic_spline_2d_deposit_value(): '
2099 write (*,
'(a, i8, i8, a, i8, i8, a)')
'array of feets of characteristics needs data of size = (', &
2100 spline%num_pts_x1, spline%num_pts_x2,
') . Passed size: (',
size(a_out, 1),
size(a_out, 2),
')'
2104 bc1 = spline%x1_bc_type
2105 bc2 = spline%x2_bc_type
2107 x1_min = spline%x1_min
2108 x2_min = spline%x2_min
2109 rh1 = spline%x1_rdelta
2110 rh2 = spline%x2_rdelta
2112 n1 = spline%num_pts_x1
2113 n2 = spline%num_pts_x2
2115 if (bc1 .eq. sll_p_periodic)
then
2118 if (bc1 .eq. sll_p_hermite)
then
2121 if (bc2 .eq. sll_p_periodic)
then
2124 if (bc2 .eq. sll_p_hermite)
then
2134 t1 = (x1(i1, i2) - x1_min)*rh1
2135 cell1 = floor(t1) + 1
2136 dx1 = t1 - real(cell1 - 1, f64)
2137 cdx1 = 1.0_f64 - dx1
2140 t2 = (x2(i1, i2) - x2_min)*rh2
2141 cell2 = floor(t2) + 1
2142 dx2 = t2 - real(cell2 - 1, f64)
2143 cdx2 = 1.0_f64 - dx2
2148 if ((real(cell1 - 1, f64)/rh1 + x1_min > x1(i1, i2)) .or. (x1(i1, i2) > real(cell1, f64)/rh1 + x1_min))
then
2149 print *,
'problem with the localization of r', &
2150 real(cell1 - 1, f64)/rh1 + x1_min, x1(i1, i2),
real(cell1, f64)/rh1 + x1_min, dx1
2152 if ((real(cell2 - 1, f64)/rh2 + x2_min > x2(i1, i2)) .or. (x2(i1, i2) > real(cell2, f64)/rh2 + x2_min))
then
2153 print *,
'problem with the localization of theta', &
2154 real(cell2 - 1, f64)/rh2 + x2_min, x2(i1, i2),
real(cell2, f64)/rh2 + x2_min, dx2
2157 cij = spline%coeffs(i1, i2)
2160 if (bc1 .eq. sll_p_periodic)
then
2161 ipm1 = mod(cell1 + n1 - 3, n1 - 1) + 1
2162 ip = mod(cell1 + n1 - 2, n1 - 1) + 1
2163 ipp1 = mod(cell1 + n1 - 1, n1 - 1) + 1
2164 ipp2 = mod(cell1 + n1, n1 - 1) + 1
2166 if (bc1 .eq. sll_p_hermite)
then
2172 if (bc2 .eq. sll_p_periodic)
then
2173 jpm1 = mod(cell2 + n2 - 3, n2 - 1) + 1
2174 jp = mod(cell2 + n2 - 2, n2 - 1) + 1
2175 jpp1 = mod(cell2 + n2 - 1, n2 - 1) + 1
2176 jpp2 = mod(cell2 + n2, n2 - 1) + 1
2178 if (bc2 .eq. sll_p_hermite)
then
2194 svalx2 = (0.5_f64*(-ax3 + ax2 + ax1) + 1._f64/6._f64)
2195 svalx3 = (2._f64/3._f64 - ax2 + 0.5_f64*ax3)
2196 svalx4 = ((1._f64 - ax3)/6._f64 + 0.5_f64*(ax2 - ax1))
2199 svaly2 = 0.5_f64*(-ay3 + ay2 + ay1) + 1._f64/6._f64
2200 svaly3 = 2._f64/3.0_f64 - ay2 + 0.5_f64*ay3
2201 svaly4 = (1._f64 - ay3)/6._f64 + 0.5_f64*(ay2 - ay1)
2203 if (ipm1 .ge. 1)
then
2204 if (jpm1 .ge. 1)
then
2205 a_out(ipm1, jpm1) = a_out(ipm1, jpm1) + cij*svalx1*svaly1
2207 a_out(ipm1, jp) = a_out(ipm1, jp) + cij*svalx1*svaly2
2208 if (jpp1 .le. n2)
then
2209 a_out(ipm1, jpp1) = a_out(ipm1, jpp1) + cij*svalx1*svaly3
2211 if (jpp2 .le. n2)
then
2212 a_out(ipm1, jpp2) = a_out(ipm1, jpp2) + cij*svalx1*svaly4
2216 if (jpm1 .ge. 1)
then
2217 a_out(ip, jpm1) = a_out(ip, jpm1) + cij*svalx2*svaly1
2219 a_out(ip, jp) = a_out(ip, jp) + cij*svalx2*svaly2
2220 if (jpp1 .le. n2)
then
2221 a_out(ip, jpp1) = a_out(ip, jpp1) + cij*svalx2*svaly3
2223 if (jpp2 .le. n2)
then
2224 a_out(ip, jpp2) = a_out(ip, jpp2) + cij*svalx2*svaly4
2227 if (ipp1 .le. n1)
then
2228 if (jpm1 .ge. 1)
then
2229 a_out(ipp1, jpm1) = a_out(ipp1, jpm1) + cij*svalx3*svaly1
2231 a_out(ipp1, jp) = a_out(ipp1, jp) + cij*svalx3*svaly2
2232 if (jpp1 .le. n2)
then
2233 a_out(ipp1, jpp1) = a_out(ipp1, jpp1) + cij*svalx3*svaly3
2235 if (jpp2 .le. n2)
then
2236 a_out(ipp1, jpp2) = a_out(ipp1, jpp2) + cij*svalx3*svaly4
2240 if (ipp2 .le. n1)
then
2241 if (jpm1 .ge. 1)
then
2242 a_out(ipp2, jpm1) = a_out(ipp2, jpm1) + cij*svalx4*svaly1
2244 a_out(ipp2, jp) = a_out(ipp2, jp) + cij*svalx4*svaly2
2245 if (jpp1 .le. n2)
then
2246 a_out(ipp2, jpp1) = a_out(ipp2, jpp1) + cij*svalx4*svaly3
2248 if (jpp2 .le. n2)
then
2249 a_out(ipp2, jpp2) = a_out(ipp2, jpp2) + cij*svalx4*svaly4
2253 if (bc1 .eq. sll_p_hermite)
then
2255 if (jpm1 .ge. 1)
then
2256 a_out(1, jpm1) = a_out(1, jpm1) + spline%coeffs(0, i2)/6._f64*svaly1
2258 a_out(1, jp) = a_out(1, jp) + spline%coeffs(0, i2)/6._f64*svaly2
2259 if (jpp1 .le. n2)
then
2260 a_out(1, jpp1) = a_out(1, jpp1) + spline%coeffs(0, i2)/6._f64*svaly3
2262 if (jpp2 .le. n2)
then
2263 a_out(1, jpp2) = a_out(1, jpp2) + spline%coeffs(0, i2)/6._f64*svaly4
2267 if (i1 .eq. n1)
then
2268 if (jpm1 .ge. 1)
then
2269 a_out(n1, jpm1) = a_out(n1, jpm1) + spline%coeffs(n1 + 1, i2)/6._f64*svaly1
2271 a_out(n1, jp) = a_out(n1, jp) + spline%coeffs(n1 + 1, i2)/6._f64*svaly2
2272 if (jpp1 .le. n2)
then
2273 a_out(n1, jpp1) = a_out(n1, jpp1) + spline%coeffs(n1 + 1, i2)/6._f64*svaly3
2275 if (jpp2 .le. n2)
then
2276 a_out(n1, jpp2) = a_out(n1, jpp2) + spline%coeffs(n1 + 1, i2)/6._f64*svaly4
2281 if (bc2 .eq. sll_p_hermite)
then
2283 if (ipm1 .ge. 1)
then
2284 a_out(ipm1, 1) = a_out(ipm1, 1) + spline%coeffs(i1, 0)/6._f64*svalx1
2286 a_out(ip, 1) = a_out(ip, 1) + spline%coeffs(i1, 0)/6._f64*svalx2
2287 if (ipp1 .ne. n1)
then
2288 a_out(ipp1, 1) = a_out(ipp1, 1) + spline%coeffs(i1, 0)/6._f64*svalx3
2290 if (ipp2 .ne. n1)
then
2291 a_out(ipp2, 1) = a_out(ipp2, 1) + spline%coeffs(i1, 0)/6._f64*svalx4
2295 if (i2 .eq. n2)
then
2296 if (ipm1 .ge. 1)
then
2297 a_out(ipm1, n2) = a_out(ipm1, n2) + spline%coeffs(i1, n2 + 1)/6._f64*svalx1
2299 a_out(ip, n1) = a_out(ip, n1) + spline%coeffs(i1, n2 + 1)/6._f64*svalx2
2300 if (ipp1 .ne. n1)
then
2301 a_out(ipp1, n2) = a_out(ipp1, n2) + spline%coeffs(i1, n2 + 1)/6._f64*svalx3
2303 if (ipp2 .ne. n1)
then
2304 a_out(ipp2, n2) = a_out(ipp2, n2) + spline%coeffs(i1, n2 + 1)/6._f64*svalx4
2312 if (bc1 .eq. sll_p_periodic)
then
2313 a_out(n1, :) = a_out(1, :)
2315 if (bc2 .eq. sll_p_periodic)
then
2316 a_out(:, n2) = a_out(:, 1)
2330 intrinsic ::
associated, int, real
2331 sll_real64,
intent(in) :: x1
2332 sll_real64,
intent(in) :: x2
2348 sll_real64 :: x1_min
2349 sll_real64 :: x2_min
2350 sll_int32 :: num_pts_x1
2351 sll_int32 :: num_pts_x2
2352 sll_real64,
dimension(:),
pointer :: coeffs_line_jm1
2353 sll_real64,
dimension(:),
pointer :: coeffs_line_j
2354 sll_real64,
dimension(:),
pointer :: coeffs_line_jp1
2355 sll_real64,
dimension(:),
pointer :: coeffs_line_jp2
2358 sll_assert((x1 .ge. spline%x1_min) .and. (x1 .le. spline%x1_max))
2359 sll_assert((x2 .ge. spline%x2_min) .and. (x2 .le. spline%x2_max))
2361 x1_min = spline%x1_min
2362 x2_min = spline%x2_min
2363 num_pts_x1 = spline%num_pts_x1
2364 num_pts_x2 = spline%num_pts_x2
2365 rh1 = spline%x1_rdelta
2366 rh2 = spline%x2_rdelta
2368 t0 = (x2 - x2_min)*rh2
2370 dx = t0 - real(cell - 1, f64)
2379 coeffs_line_jm1 => spline%coeffs(0:num_pts_x1 + 2, cell - 1)
2380 coeffs_line_j => spline%coeffs(0:num_pts_x1 + 2, cell)
2381 coeffs_line_jp1 => spline%coeffs(0:num_pts_x1 + 2, cell + 1)
2382 coeffs_line_jp2 => spline%coeffs(0:num_pts_x1 + 2, cell + 2)
2389 t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
2390 t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
2414 intrinsic ::
associated, int, real
2415 sll_real64,
intent(in) :: x1
2416 sll_real64,
intent(in) :: x2
2432 sll_real64 :: x1_min
2433 sll_real64 :: x2_min
2434 sll_int32 :: num_pts_x1
2435 sll_int32 :: num_pts_x2
2436 sll_real64,
dimension(:),
pointer :: coeffs_line_jm1
2437 sll_real64,
dimension(:),
pointer :: coeffs_line_j
2438 sll_real64,
dimension(:),
pointer :: coeffs_line_jp1
2439 sll_real64,
dimension(:),
pointer :: coeffs_line_jp2
2442 sll_assert((x1 .ge. spline%x1_min) .and. (x1 .le. spline%x1_max))
2443 sll_assert((x2 .ge. spline%x2_min) .and. (x2 .le. spline%x2_max))
2445 x1_min = spline%x1_min
2446 x2_min = spline%x2_min
2447 num_pts_x1 = spline%num_pts_x1
2448 num_pts_x2 = spline%num_pts_x2
2449 rh1 = spline%x1_rdelta
2450 rh2 = spline%x2_rdelta
2452 t0 = (x2 - x2_min)*rh2
2454 dx = t0 - real(cell - 1, f64)
2463 coeffs_line_jm1 => spline%coeffs(0:num_pts_x1 + 2, cell - 1)
2464 coeffs_line_j => spline%coeffs(0:num_pts_x1 + 2, cell)
2465 coeffs_line_jp1 => spline%coeffs(0:num_pts_x1 + 2, cell + 1)
2466 coeffs_line_jp2 => spline%coeffs(0:num_pts_x1 + 2, cell + 2)
2473 t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
2474 t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
2500 intrinsic ::
associated, int, real
2501 sll_real64,
intent(in) :: x1
2502 sll_real64,
intent(in) :: x2
2517 sll_real64 :: x1_min
2518 sll_real64 :: x2_min
2519 sll_int32 :: num_pts_x1
2520 sll_int32 :: num_pts_x2
2521 sll_real64,
dimension(:),
pointer :: coeffs_line_jm1
2522 sll_real64,
dimension(:),
pointer :: coeffs_line_j
2523 sll_real64,
dimension(:),
pointer :: coeffs_line_jp1
2524 sll_real64,
dimension(:),
pointer :: coeffs_line_jp2
2527 sll_assert((x1 .ge. spline%x1_min) .and. (x1 .le. spline%x1_max))
2528 sll_assert((x2 .ge. spline%x2_min) .and. (x2 .le. spline%x2_max))
2530 x1_min = spline%x1_min
2531 x2_min = spline%x2_min
2532 num_pts_x1 = spline%num_pts_x1
2533 num_pts_x2 = spline%num_pts_x2
2534 rh1 = spline%x1_rdelta
2535 rh2 = spline%x2_rdelta
2537 t0 = (x2 - x2_min)*rh2
2539 dx = t0 - real(cell - 1, f64)
2548 coeffs_line_jm1 => spline%coeffs(0:num_pts_x1 + 2, cell - 1)
2549 coeffs_line_j => spline%coeffs(0:num_pts_x1 + 2, cell)
2550 coeffs_line_jp1 => spline%coeffs(0:num_pts_x1 + 2, cell + 1)
2551 coeffs_line_jp2 => spline%coeffs(0:num_pts_x1 + 2, cell + 2)
2556 t1 = 2.0_f64*(cim1 - 2.0_f64*ci + cip1)
2557 t2 = -cim1 + 3.0_f64*(ci - cip1) + cip2
2564 sll_real64,
dimension(:),
intent(out) :: coeff
2567 sll_int32 :: num_pts_x1
2568 sll_int32 :: num_pts_x2
2570 num_pts_x1 = spline%num_pts_x1
2571 num_pts_x2 = spline%num_pts_x2
2573 sll_assert(
size(coeff) >= (num_pts_x1 + 2)*(num_pts_x2 + 2))
2575 do j = 1, num_pts_x2 + 2
2576 do i = 1, num_pts_x1 + 2
2577 coeff(i + (num_pts_x1 + 2)*(j - 1)) = spline%coeffs(i - 1, j - 1)
2585 sll_deallocate(spline%d1, ierr)
2586 sll_deallocate(spline%d2, ierr)
2587 sll_deallocate(spline%coeffs, ierr)
2588 spline%data => null()
2589 if (
associated(spline%x1_min_slopes))
then
2590 sll_deallocate(spline%x1_min_slopes, ierr)
2592 if (
associated(spline%x1_max_slopes))
then
2593 sll_deallocate(spline%x1_max_slopes, ierr)
2595 if (
associated(spline%x2_min_slopes))
then
2596 sll_deallocate(spline%x2_min_slopes, ierr)
2598 if (
associated(spline%x2_max_slopes))
then
2599 sll_deallocate(spline%x2_max_slopes, ierr)
2601 if (
associated(spline%x1_min_slopes_coeffs))
then
2602 sll_deallocate(spline%x1_min_slopes_coeffs, ierr)
2604 if (
associated(spline%x1_max_slopes_coeffs))
then
2605 sll_deallocate(spline%x1_max_slopes_coeffs, ierr)
2607 if (
associated(spline%x2_min_slopes_coeffs))
then
2608 sll_deallocate(spline%x2_min_slopes_coeffs, ierr)
2610 if (
associated(spline%x2_max_slopes_coeffs))
then
2611 sll_deallocate(spline%x2_max_slopes_coeffs, ierr)
2618 sll_real64,
intent(in) :: alpha
2619 sll_real64,
intent(out) :: output_array(:)
2622 sll_real64 :: alpha0
2624 sll_real64 :: alpha1
2627 sll_int32 :: num_cells
2630 sll_assert(
size(output_array) == spline%n_points)
2632 alpha0 = alpha*spline%rdelta
2633 dcell = floor(alpha0)
2634 alpha1 = alpha0 - real(dcell, f64)
2635 select case (spline%bc_type)
2636 case (sll_p_periodic)
2637 num_cells = spline%n_points - 1
2638 case (sll_p_hermite)
2639 num_cells = spline%n_points
2642 do i = max(1, 1 - dcell), min(num_cells, num_cells - dcell)
2647 if (spline%bc_type == sll_p_periodic)
then
2648 do i = 1, max(1, 1 - dcell) - 1
2649 cell = modulo(i + dcell - 1, num_cells) + 1
2654 do i = min(num_cells, num_cells - dcell) + 1, num_cells
2655 cell = modulo(i + dcell - 1, num_cells) + 1
2661 output_array(spline%n_points) = output_array(1)
2669 output_array(i) = output_array(1)
2672 cell = spline%n_points
2675 do i = num_cells - dcell + 1, spline%n_points - 1
2676 output_array(i) = output_array(spline%n_points)
2687 sll_int32,
intent(in) :: cell
2688 sll_real64,
intent(in) :: dx
2689 sll_real64,
intent(out) :: out
2701 cim1 = spline%coeffs(cell - 1)
2702 ci = spline%coeffs(cell)
2703 cip1 = spline%coeffs(cell + 1)
2704 cip2 = spline%coeffs(cell + 2)
2707 t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
2708 t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
2709 out =
inv_6*(t2 + t4)
Describe different boundary conditions.
Provides capabilities for data and derivative interpolation with cubic B-splines and different bounda...
real(kind=f64), parameter inv_6
subroutine interpolate_pointer_values(ptr_in, ptr_out, n, spline)
returns the values of the images of a collection of abscissae,
subroutine, public sll_s_cubic_spline_2d_compute_interpolant(data, spline)
Computes the spline coefficients for the given data. The coeffcients are first computed in the second...
subroutine, public sll_s_cubic_spline_2d_deposit_value(x1, x2, spline, a_out)
Updated distribution function at time .
real(kind=f64) function, public sll_f_cubic_spline_1d_eval_deriv(x, spline)
returns the value of the derivative at the image of the abscissa 'x', using the spline coefficients s...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval_deriv_x2(x1, x2, spline)
Returns the interpolated value of the derivative.
subroutine compute_spline_2d_hrmt_hrmt(data, spline)
subroutine, public sll_s_cubic_spline_1d_compute_interpolant(f, spline)
Computes the cubic spline coefficients corresponding to the 1D data array 'f'.
subroutine, public sll_s_cubic_spline_1d_init(self, num_points, xmin, xmax, bc_type, sl, sr, fast_algorithm)
Returns a pointer to a heap-allocated cubic spline object.
subroutine compute_spline_1d_hermite_aux(f, num_pts, d, slope_l, slope_r, delta, coeffs)
subroutine, public sll_s_cubic_spline_2d_get_coeff(spline, coeff)
subroutine, public sll_s_cubic_spline_1d_eval_disp(spline, alpha, output_array)
Computes the interpolated values at each grid point replaced by alpha for the precomputed spline coef...
subroutine interpolate_pointer_derivatives(ptr_in, ptr_out, num_pts, spline)
analogous to the sll_s_cubic_spline_1d_eval_array_deriv() subroutine but its input and output arrays ...
real(kind=f64) function interpolate_derivative_aux(x, xmin, rh, coeffs)
subroutine spline_interpolate_from_interpolant_cell_dx(spline, cell, dx, out)
Helper function for sll_s_cubic_spline_1d_eval_disp: evaluate spline in given cell and normalized dis...
subroutine, public sll_s_cubic_spline_1d_eval_array(a_in, a_out, n, spline)
returns the values of the images of a collection of abscissae, represented by a 1D array in another o...
subroutine compute_spline_2d_hrmt_prdc(data, spline)
subroutine compute_spline_1d_periodic(f, spline)
subroutine compute_spline_1d_hermite(f, spline)
subroutine, public sll_s_cubic_spline_1d_eval_array_deriv(array_in, array_out, num_pts, spline)
returns the values of the derivatives evaluated at a collection of abscissae stored in the input 1D a...
subroutine compute_spline_1d_periodic_aux(f, num_pts, d, coeffs)
subroutine compute_spline_2d_prdc_prdc(data, spline)
real(kind=f64) function, public sll_f_cubic_spline_1d_eval(x, spline)
returns the value of the interpolated image of the abscissa x, using the spline coefficients stored i...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval_deriv_x1(x1, x2, spline)
Returns the interpolated value of the derivative in the x1 direction at the point (x1,...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval(x1, x2, spline)
Returns the interpolated value of the image of the point (x1,x2) using the spline decomposition store...
real(kind=f64) function interpolate_value_aux(x, xmin, rh, coeffs)
subroutine, public sll_s_cubic_spline_1d_free(spline)
deallocate the sll_t_cubic_spline_1d object
subroutine, public sll_s_cubic_spline_2d_init(this, num_pts_x1, num_pts_x2, x1_min, x1_max, x2_min, x2_max, x1_bc_type, x2_bc_type, const_slope_x1_min, const_slope_x1_max, const_slope_x2_min, const_slope_x2_max, x1_min_slopes, x1_max_slopes, x2_min_slopes, x2_max_slopes)
subroutine compute_spline_2d_prdc_hrmt(data, spline)
type(sll_t_cubic_spline_2d) function, pointer sll_f_new_cubic_spline_2d(num_pts_x1, num_pts_x2, x1_min, x1_max, x2_min, x2_max, x1_bc_type, x2_bc_type, const_slope_x1_min, const_slope_x1_max, const_slope_x2_min, const_slope_x2_max, x1_min_slopes, x1_max_slopes, x2_min_slopes, x2_max_slopes)
Returns a pointer to a heap-allocated 2D cubic spline object.
subroutine, public sll_s_cubic_spline_2d_free(spline)
Tridiagonal system solver.
subroutine, public sll_s_setup_cyclic_tridiag(a, n, cts, ipiv)
Give the factorization of the matrix in argument.
subroutine, public sll_s_solve_cyclic_tridiag_double(cts, ipiv, b, n, x)
Solves tridiagonal system.
basic type for one-dimensional cubic spline data.
Basic type for one-dimensional cubic spline data.