! ------------ ---------------------------------------------------------- ! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator ! http://lammps.sandia.gov, Sandia National Laboratories ! Steve Plimpton, sjplimp@sandia.gov ! ! Copyright (2003) Sandia Corporation. Under the terms of Contract ! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains ! certain rights in this software. This software is distributed under ! the GNU General Public License. ! ! See the README file in the top-level LAMMPS directory. ! ! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu !------------------------------------------------------------------------- module Spline2 !************************************************************************************ ! ! Two-dimensional cubic spline function. ! !--------------------------------------------------------------------------------------------------- ! ! Intel Fortran ! ! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017 ! !*************************************************************************************************** use Spline1 use iso_c_binding, only : c_int, c_double, c_char implicit none contains !****************************************************************************************** subroutine CreateSpline2 ( CL, CD, CR, CU, N1, N2, N, P1, P2, F, Fxx, Fyy, Fxxyy, FF, MM, DD, K0, K1, K2 ) integer(c_int), intent(in) :: CL, CD, CR, CU, N1, N2, N real(c_double), dimension(0:N1-1), intent(in) :: P1 real(c_double), dimension(0:N2-1), intent(in) :: P2 real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy real(c_double), dimension(0:N-1), intent(inout) :: FF, MM, DD, K0, K1, K2 integer(c_int) :: II !------------------------------------------------------------------------------------------- do II = 0, N2 - 1 FF(0:N1-1) = F(0:N1-1,II) MM(0) = Fxx(0,II) MM(N1-1) = Fxx(N1-1,II) call CreateSpline1 ( CL, CR, N1, P1, FF, MM, DD, K0, K1, K2 ) Fxx(0:N1-1,II) = MM(0:N1-1) end do do II = 0, N1 - 1 MM(0) = Fyy(II,0) MM(N-1) = Fyy(II,N2-1) FF(0:N2-1) = F(II,0:N2-1) call CreateSpline1 ( CD, CU, N2, P2, FF, MM, DD, K0, K1, K2 ) Fyy(II,0:N2-1) = MM(0:N2-1) end do FF(0:N1-1) = Fyy(0:N1-1,0 ) call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, K2 ) Fxxyy(0:N1-1,0) = MM(0:N1-1) FF(0:N1-1) = Fyy(0:N1-1,N2-1 ) call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, k2 ) Fxxyy(0:N1-1,N2-1) = MM(0:N1-1) do II = 1, N1 - 2 MM(0) = Fxxyy(II,0) MM(N-1) = Fxxyy(II,N2-1) FF(0:N2-1) = Fxx(II,0:N2-1) call CreateSpline1 ( 2 , 2, N2, P2, FF, MM, DD, K0, K1, K2 ) Fxxyy(II,0:N2-1) = MM(0:N2-1) end do end subroutine CreateSpline2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine CreateSpline2Ext ( CL, CD, CR, CU, N1, N1A, N2, N2A, N, P1, P2, F, Fxx, Fyy, Fxxyy, FF, MM, DD, K0, K1, K2 ) integer(c_int), intent(in) :: CL, CD, CR, CU, N1, N1A, N2, N2A, N real(c_double), dimension(0:N1-1), intent(in) :: P1 real(c_double), dimension(0:N2-1), intent(in) :: P2 real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy real(c_double), dimension(0:N-1), intent(inout) :: FF, MM, DD, K0, K1, K2 integer(c_int) :: II !------------------------------------------------------------------------------------------- Fxx = 0.0d+00 Fyy = 0.0d+00 Fxxyy = 0.0d+00 do II = 0, N2A FF(0:N1-1) = F(0:N1-1,II) MM(0) = Fxx(0,II) MM(N1-1) = Fxx(N1-1,II) call CreateSpline1 ( CL, CR, N1, P1, FF, MM, DD, K0, K1, K2 ) Fxx(0:N1-1,II) = MM(0:N1-1) end do do II = N2A + 1, N2 - 1 FF(0:N1-N1A-1) = F(N1A:N1-1,II) MM(0) = Fxx(N1A,II) MM(N1-N1A-1) = Fxx(N1-1,II) call CreateSpline1 ( CL, CR, N1 - N1A, P1, FF, MM, DD, K0, K1, K2 ) Fxx(N1A:N1-1,II) = MM(0:N1-N1A-1) end do do II = 0, N1A - 1 MM(0) = Fyy(II,0) MM(N2A) = Fyy(II,N2A) FF(0:N2A) = F(II,0:N2A) call CreateSpline1 ( CD, CU, N2A + 1, P2, FF, MM, DD, K0, K1, K2 ) Fyy(II,0:N2A) = MM(0:N2A) end do do II = N1A, N1 - 1 MM(0) = Fyy(II,0) MM(N-1) = Fyy(II,N2-1) FF(0:N2-1) = F(II,0:N2-1) call CreateSpline1 ( CD, CU, N2, P2, FF, MM, DD, K0, K1, K2 ) Fyy(II,0:N2-1) = MM(0:N2-1) end do FF(0:N1-1) = Fyy(0:N1-1,0) call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, K2 ) Fxxyy(0:N1-1,0) = MM(0:N1-1) FF(0:N1A) = Fyy(0:N1A,N2A) call CreateSpline1 ( 3, 3, N1A + 1, P1, FF, MM, DD, K0, K1, K2 ) Fxxyy(0:N1A,N2A) = MM(0:N1A) FF(0:N1-N1A-1) = Fyy(N1A:N1-1,N2-1 ) call CreateSpline1 ( 3, 3, N1-N1A, P1, FF, MM, DD, K0, K1, K2 ) Fxxyy(N1A:N1-1,N2-1) = MM(0:N1-N1A-1) do II = 1, N1A MM(0) = Fxxyy(II,0) MM(N2A) = Fxxyy(II,N2A) FF(0:N2A) = Fxx(II,0:N2A) call CreateSpline1 ( 2 , 2, N2A + 1, P2, FF, MM, DD, K0, K1, K2 ) Fxxyy(II,0:N2A) = MM(0:N2A) end do do II = N1A + 1, N1 - 2 MM(0) = Fxxyy(II,0) MM(N-1) = Fxxyy(II,N2-1) FF(0:N2-1) = Fxx(II,0:N2-1) call CreateSpline1 ( 2 , 2, N2, P2, FF, MM, DD, K0, K1, K2 ) Fxxyy(II,0:N2-1) = MM(0:N2-1) end do end subroutine CreateSpline2Ext !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real(c_double) function CalcSpline2_0 ( i, j, X, Y, N1, N2, P1, P2, F, Fxx, Fyy, Fxxyy ) !!! integer(c_int), intent(in) :: i, j, N1, N2 real(c_double), intent(in) :: X, Y real(c_double), dimension(0:N1-1), intent(in) :: P1 real(c_double), dimension(0:N2-1), intent(in) :: P2 real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy integer(c_int) :: i1, j1 real(c_double) :: T, Gy_0, Gy_1, Gxxy_0, Gxxy_1 !------------------------------------------------------------------------------------------- i1 = i - 1 j1 = j - 1 T = P2(j) - P2(j1) Gy_0 = ValueSpline1_0 ( Y, P2(j), P2(j1), F(i,j), F(i,j1), Fyy(i,j), Fyy(i,j1), T ) Gy_1 = ValueSpline1_0 ( Y, P2(j), P2(j1), F(i1,j), F(i1,j1), Fyy(i1,j), Fyy(i1,j1), T ) Gxxy_0 = ValueSpline1_0 ( Y, P2(j), P2(j1), Fxx(i,j), Fxx(i,j1), Fxxyy(i,j), Fxxyy(i,j1), T ) Gxxy_1 = ValueSpline1_0 ( Y, P2(j), P2(j1), Fxx(i1,j), Fxx(i1,j1), Fxxyy(i1,j), Fxxyy(i1,j1), T ) CalcSpline2_0 = ValueSpline1_0 ( X, P1(i), P1(i1), Gy_0, Gy_1,Gxxy_0, Gxxy_1, P1(i) - P1(i1) ) end function CalcSpline2_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine CalcSpline2_1 ( S, Sx1, Sy1, i, j, X, Y, N1, N2, P1, P2, F, Fxx, Fyy, Fxxyy ) !!! real(c_double), intent(out) :: S, Sx1, Sy1 integer(c_int), intent(in) :: i, j, N1, N2 real(c_double), intent(in) :: X, Y real(c_double), dimension(0:N1-1), intent(in) :: P1 real(c_double), dimension(0:N2-1), intent(in) :: P2 real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy integer(c_int) :: i1, j1 real(c_double) :: T, Gy_0, Gy_1, Gxxy_0, Gxxy_1 real(c_double) :: Gyy_0, Gyy_1, Gxxyy_0, Gxxyy_1 !------------------------------------------------------------------------------------------- i1 = i - 1 j1 = j - 1 T = P2(j) - P2(j1) call ValueSpline1_1 ( Gy_0, Gyy_0, Y, P2(j), P2(j1), F(i,j), F(i,j1), Fyy(i,j), Fyy(i,j1), T ) call ValueSpline1_1 ( Gy_1, Gyy_1, Y, P2(j), P2(j1), F(i1,j), F(i1,j1), Fyy(i1,j), Fyy(i1,j1), T ) call ValueSpline1_1 ( Gxxy_0, Gxxyy_0, Y, P2(j), P2(j1), Fxx(i,j), Fxx(i,j1), Fxxyy(i,j), Fxxyy(i,j1), T ) call ValueSpline1_1 ( Gxxy_1, Gxxyy_1, Y, P2(j), P2(j1), Fxx(i1,j), Fxx(i1,j1), Fxxyy(i1,j), Fxxyy(i1,j1), T ) call ValueSpline1_1 ( S, Sx1, X, P1(i), P1(i1), Gy_0, Gy_1,Gxxy_0, Gxxy_1, P1(i) - P1(i1) ) Sy1 = ValueSpline1_0 ( X, P1(i), P1(i1), Gyy_0, Gyy_1,Gxxyy_0, Gxxyy_1, P1(i) - P1(i1) ) end subroutine CalcSpline2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module Spline2 !********************************************************************************