Merge branch 'master' into reconnection
This commit is contained in:
commit
7939c98cf5
@ -46,6 +46,28 @@ module operators
|
|||||||
integer , save :: imi, imd, img, imc, iml, im1, im2
|
integer , save :: imi, imd, img, imc, iml, im1, im2
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
! interfaces for procedure pointers
|
||||||
|
!
|
||||||
|
abstract interface
|
||||||
|
subroutine derivative_1st_iface(d, h, u, v)
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
end subroutine
|
||||||
|
subroutine derivative_2nd_iface(d, h, u, v)
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
end subroutine
|
||||||
|
end interface
|
||||||
|
|
||||||
|
! procedure pointers for derivatives
|
||||||
|
!
|
||||||
|
procedure(derivative_1st_iface), pointer, save :: derivative_1st => null()
|
||||||
|
procedure(derivative_2nd_iface), pointer, save :: derivative_2nd => null()
|
||||||
|
|
||||||
! by default everything is public
|
! by default everything is public
|
||||||
!
|
!
|
||||||
private
|
private
|
||||||
@ -81,19 +103,15 @@ module operators
|
|||||||
!
|
!
|
||||||
subroutine initialize_operators(status)
|
subroutine initialize_operators(status)
|
||||||
|
|
||||||
! local variables are not implicit by default
|
use interpolations, only : order
|
||||||
!
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
! subroutine arguments
|
|
||||||
!
|
|
||||||
integer, intent(out) :: status
|
integer, intent(out) :: status
|
||||||
!
|
!
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
! set timer descriptions
|
|
||||||
!
|
|
||||||
call set_timer('operators:: initialization', imi)
|
call set_timer('operators:: initialization', imi)
|
||||||
call set_timer('operators:: divergence' , imd)
|
call set_timer('operators:: divergence' , imd)
|
||||||
call set_timer('operators:: gradient' , img)
|
call set_timer('operators:: gradient' , img)
|
||||||
@ -102,18 +120,30 @@ module operators
|
|||||||
call set_timer('operators:: 1st derivative', im1)
|
call set_timer('operators:: 1st derivative', im1)
|
||||||
call set_timer('operators:: 2nd derivative', im2)
|
call set_timer('operators:: 2nd derivative', im2)
|
||||||
|
|
||||||
! start accounting time for the module initialization/finalization
|
|
||||||
!
|
|
||||||
call start_timer(imi)
|
call start_timer(imi)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
! reset the status flag
|
select case(order)
|
||||||
!
|
case(9)
|
||||||
|
derivative_1st => derivative_1st_9o
|
||||||
|
derivative_2nd => derivative_2nd_9o
|
||||||
|
|
||||||
|
case(7)
|
||||||
|
derivative_1st => derivative_1st_7o
|
||||||
|
derivative_2nd => derivative_2nd_7o
|
||||||
|
|
||||||
|
case(5)
|
||||||
|
derivative_1st => derivative_1st_5o
|
||||||
|
derivative_2nd => derivative_2nd_5o
|
||||||
|
|
||||||
|
case default
|
||||||
|
derivative_1st => derivative_1st_3o
|
||||||
|
derivative_2nd => derivative_2nd_3o
|
||||||
|
end select
|
||||||
|
|
||||||
status = 0
|
status = 0
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
! stop accounting time for the module initialization/finalization
|
|
||||||
!
|
|
||||||
call stop_timer(imi)
|
call stop_timer(imi)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
@ -502,195 +532,782 @@ module operators
|
|||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
! subroutine DERIVATIVE_1ST:
|
! subroutine DERIVATIVE_1ST_3O:
|
||||||
! -------------------------
|
! ----------------------------
|
||||||
!
|
!
|
||||||
! Subroutine calculates the first order derivative of the input scalar field
|
! Subroutine calculates the first order derivative of the input scalar field
|
||||||
! along a given direction.
|
! along the given direction with the 3rd order approximation.
|
||||||
!
|
!
|
||||||
! Arguments:
|
! Arguments:
|
||||||
!
|
!
|
||||||
! dir - the direction of derivative;
|
! d - the direction of derivative;
|
||||||
! dh - the spacial interval;
|
! h - the spacial interval;
|
||||||
! u - the input scalar field;
|
! u - the input scalar field;
|
||||||
! v - the first derivative of the input field;
|
! v - the first derivative of the input field;
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
subroutine derivative_1st(dir, dh, u, v)
|
subroutine derivative_1st_3o(d, h, u, v)
|
||||||
|
|
||||||
! local variables are not implicit by default
|
|
||||||
!
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
! subroutine arguments
|
integer , intent(in) :: d
|
||||||
!
|
real(kind=8) , intent(in) :: h
|
||||||
integer , intent(in) :: dir
|
|
||||||
real(kind=8) , intent(in) :: dh
|
|
||||||
real(kind=8), dimension(:,:,:), intent(in) :: u
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
real(kind=8), dimension(:,:,:), intent(out) :: v
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
! local variables
|
|
||||||
!
|
|
||||||
integer :: m0, m1, m2
|
integer :: m0, m1, m2
|
||||||
!
|
!
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
#ifdef DEBUG
|
|
||||||
! return if the direction is wrong
|
|
||||||
!
|
|
||||||
if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return
|
|
||||||
#endif /* DEBUG */
|
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
! start accounting time for the 1st derivative calculation
|
|
||||||
!
|
|
||||||
call start_timer(im1)
|
call start_timer(im1)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
! prepare index limits
|
m0 = size(u, d)
|
||||||
!
|
|
||||||
m0 = size(u, dir)
|
|
||||||
m1 = m0 - 1
|
m1 = m0 - 1
|
||||||
m2 = m0 - 2
|
m2 = m0 - 2
|
||||||
|
|
||||||
! select the direction
|
select case(d)
|
||||||
!
|
|
||||||
select case(dir)
|
|
||||||
|
|
||||||
! derivative along the X direction
|
|
||||||
!
|
|
||||||
case(1)
|
case(1)
|
||||||
|
|
||||||
v(2:m1,:,:) = 0.5d+00 * (u(3:m0,:,:) - u(1:m2,:,:)) / dh
|
v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h
|
||||||
v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / dh
|
v(2:m1,:,:) = 5.0d-01 * (u(3:m0,:,:) - u(1:m2,:,:)) / h
|
||||||
v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / dh
|
v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h
|
||||||
|
|
||||||
! derivative along the Y direction
|
|
||||||
!
|
|
||||||
case(2)
|
case(2)
|
||||||
|
|
||||||
v(:,2:m1,:) = 0.5d+00 * (u(:,3:m0,:) - u(:,1:m2,:)) / dh
|
v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h
|
||||||
v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / dh
|
v(:,2:m1,:) = 5.0d-01 * (u(:,3:m0,:) - u(:,1:m2,:)) / h
|
||||||
v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / dh
|
v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h
|
||||||
|
|
||||||
#if NDIMS == 3
|
#if NDIMS == 3
|
||||||
! derivative along the Z direction
|
|
||||||
!
|
|
||||||
case(3)
|
case(3)
|
||||||
|
|
||||||
v(:,:,2:m1) = 0.5d+00 * (u(:,:,3:m0) - u(:,:,1:m2)) / dh
|
v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h
|
||||||
v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / dh
|
v(:,:,2:m1) = 5.0d-01 * (u(:,:,3:m0) - u(:,:,1:m2)) / h
|
||||||
v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / dh
|
v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h
|
||||||
#endif /* NDIMS == 3 */
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
end select
|
end select
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
! stop accounting time for the 1st derivative calculation
|
|
||||||
!
|
|
||||||
call stop_timer(im1)
|
call stop_timer(im1)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
end subroutine derivative_1st
|
end subroutine derivative_1st_3o
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
! subroutine DERIVATIVE_2ND:
|
! subroutine DERIVATIVE_1ST_5O:
|
||||||
! -------------------------
|
! ----------------------------
|
||||||
!
|
!
|
||||||
! Subroutine calculates the second order derivative of the input scalar field
|
! Subroutine calculates the first order derivative of the input scalar field
|
||||||
! along a given direction.
|
! along the given direction with the 5th order approximation.
|
||||||
!
|
!
|
||||||
! Arguments:
|
! Arguments:
|
||||||
!
|
!
|
||||||
! dir - the direction of derivative;
|
! d - the direction of derivative;
|
||||||
! dh - the spacial interval;
|
! h - the spacial interval;
|
||||||
! u - the input scalar field;
|
! u - the input scalar field;
|
||||||
! v - the output scalar field representing the second derivative of u;
|
! v - the first derivative of the input field;
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
subroutine derivative_2nd(dir, dh, u, v)
|
subroutine derivative_1st_5o(d, h, u, v)
|
||||||
|
|
||||||
! local variables are not implicit by default
|
|
||||||
!
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
! subroutine arguments
|
integer , intent(in) :: d
|
||||||
!
|
real(kind=8) , intent(in) :: h
|
||||||
integer , intent(in) :: dir
|
|
||||||
real(kind=8) , intent(in) :: dh
|
|
||||||
real(kind=8), dimension(:,:,:), intent(in) :: u
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
real(kind=8), dimension(:,:,:), intent(out) :: v
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
! local variables
|
integer :: m0, m1, m2, m3, m4
|
||||||
!
|
|
||||||
integer :: m0, m1, m2
|
|
||||||
!
|
!
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
#ifdef DEBUG
|
#ifdef PROFILE
|
||||||
! return if the direction is wrong
|
call start_timer(im1)
|
||||||
!
|
#endif /* PROFILE */
|
||||||
if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return
|
|
||||||
#endif /* DEBUG */
|
m0 = size(u, d)
|
||||||
|
m1 = m0 - 1
|
||||||
|
m2 = m0 - 2
|
||||||
|
m3 = m0 - 3
|
||||||
|
m4 = m0 - 4
|
||||||
|
|
||||||
|
select case(d)
|
||||||
|
|
||||||
|
case(1)
|
||||||
|
|
||||||
|
v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h
|
||||||
|
v(2 ,:,:) = 5.0d-01 * (u(3 ,:,:) - u(1 ,:,:)) / h
|
||||||
|
v(3:m2,:,:) = (8.0d+00 * (u(4:m1,:,:) - u(2:m3,:,:)) &
|
||||||
|
- (u(5:m0,:,:) - u(1:m4,:,:))) / (1.2d+01 * h)
|
||||||
|
v( m1,:,:) = 5.0d-01 * (u( m0,:,:) - u( m2,:,:)) / h
|
||||||
|
v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
|
||||||
|
v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h
|
||||||
|
v(:,2 ,:) = 5.0d-01 * (u(:,3 ,:) - u(:,1 ,:)) / h
|
||||||
|
v(:,3:m2,:) = (8.0d+00 * (u(:,4:m1,:) - u(:,2:m3,:)) &
|
||||||
|
- (u(:,5:m0,:) - u(:,1:m4,:))) / (1.2d+01 * h)
|
||||||
|
v(:, m1,:) = 5.0d-01 * (u(:, m0,:) - u(:, m2,:)) / h
|
||||||
|
v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
case(3)
|
||||||
|
|
||||||
|
v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h
|
||||||
|
v(:,:,2 ) = 5.0d-01 * (u(:,:,3 ) - u(:,:,1 )) / h
|
||||||
|
v(:,:,3:m2) = (8.0d+00 * (u(:,:,4:m1) - u(:,:,2:m3)) &
|
||||||
|
- (u(:,:,5:m0) - u(:,:,1:m4))) / (1.2d+01 * h)
|
||||||
|
v(:,:, m1) = 5.0d-01 * (u(:,:, m0) - u(:,:, m2)) / h
|
||||||
|
v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
! start accounting time for the 2nd derivative calculation
|
call stop_timer(im1)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
|
end subroutine derivative_1st_5o
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine DERIVATIVE_1ST_7O:
|
||||||
|
! ----------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculates the first order derivative of the input scalar field
|
||||||
|
! along the given direction with the 7th order approximation.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! d - the direction of derivative;
|
||||||
|
! h - the spacial interval;
|
||||||
|
! u - the input scalar field;
|
||||||
|
! v - the first derivative of the input field;
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine derivative_1st_7o(d, h, u, v)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
|
integer :: m0, m1, m2, m3, m4, m5, m6
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
#ifdef PROFILE
|
||||||
|
call start_timer(im1)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
m0 = size(u, d)
|
||||||
|
m1 = m0 - 1
|
||||||
|
m2 = m0 - 2
|
||||||
|
m3 = m0 - 3
|
||||||
|
m4 = m0 - 4
|
||||||
|
m5 = m0 - 5
|
||||||
|
m6 = m0 - 6
|
||||||
|
|
||||||
|
select case(d)
|
||||||
|
|
||||||
|
case(1)
|
||||||
|
|
||||||
|
v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h
|
||||||
|
v(2 ,:,:) = 5.0d-01 * (u(3 ,:,:) - u(1 ,:,:)) / h
|
||||||
|
v(3 ,:,:) = (8.0d+00 * (u(4 ,:,:) - u(2 ,:,:)) &
|
||||||
|
- (u(5 ,:,:) - u(1 ,:,:))) / (1.2d+01 * h)
|
||||||
|
v(4:m3,:,:) = (4.5d+01 * (u(5:m2,:,:) - u(3:m4,:,:)) &
|
||||||
|
- 9.0d+00 * (u(6:m1,:,:) - u(2:m5,:,:)) &
|
||||||
|
+ (u(7:m0,:,:) - u(1:m6,:,:))) / (6.0d+01 * h)
|
||||||
|
v( m2,:,:) = (8.0d+00 * (u( m1,:,:) - u( m3,:,:)) &
|
||||||
|
- (u( m0,:,:) - u( m4,:,:))) / (1.2d+01 * h)
|
||||||
|
v( m1,:,:) = 5.0d-01 * (u( m0,:,:) - u( m2,:,:)) / h
|
||||||
|
v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
|
||||||
|
v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h
|
||||||
|
v(:,2 ,:) = 5.0d-01 * (u(:,3 ,:) - u(:,1 ,:)) / h
|
||||||
|
v(:,3 ,:) = (8.0d+00 * (u(:,4 ,:) - u(:,2 ,:)) &
|
||||||
|
- (u(:,5 ,:) - u(:,1 ,:))) / (1.2d+01 * h)
|
||||||
|
v(:,4:m3,:) = (4.5d+01 * (u(:,5:m2,:) - u(:,3:m4,:)) &
|
||||||
|
- 9.0d+00 * (u(:,6:m1,:) - u(:,2:m5,:)) &
|
||||||
|
+ (u(:,7:m0,:) - u(:,1:m6,:))) / (6.0d+01 * h)
|
||||||
|
v(:, m2,:) = (8.0d+00 * (u(:, m1,:) - u(:, m3,:)) &
|
||||||
|
- (u(:, m0,:) - u(:, m4,:))) / (1.2d+01 * h)
|
||||||
|
v(:, m1,:) = 5.0d-01 * (u(:, m0,:) - u(:, m2,:)) / h
|
||||||
|
v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
case(3)
|
||||||
|
|
||||||
|
v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h
|
||||||
|
v(:,:,2 ) = 5.0d-01 * (u(:,:,3 ) - u(:,:,1 )) / h
|
||||||
|
v(:,:,3 ) = (8.0d+00 * (u(:,:,4 ) - u(:,:,2 )) &
|
||||||
|
- (u(:,:,5 ) - u(:,:,1 ))) / (1.2d+01 * h)
|
||||||
|
v(:,:,4:m3) = (4.5d+01 * (u(:,:,5:m2) - u(:,:,3:m4)) &
|
||||||
|
- 9.0d+00 * (u(:,:,6:m1) - u(:,:,2:m5)) &
|
||||||
|
+ (u(:,:,7:m0) - u(:,:,1:m6))) / (6.0d+01 * h)
|
||||||
|
v(:,:, m2) = (8.0d+00 * (u(:,:, m1) - u(:,:, m3)) &
|
||||||
|
- (u(:,:, m0) - u(:,:, m4))) / (1.2d+01 * h)
|
||||||
|
v(:,:, m1) = 5.0d-01 * (u(:,:, m0) - u(:,:, m2)) / h
|
||||||
|
v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
#ifdef PROFILE
|
||||||
|
call stop_timer(im1)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine derivative_1st_7o
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine DERIVATIVE_1ST_9O:
|
||||||
|
! ----------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculates the first order derivative of the input scalar field
|
||||||
|
! along the given direction with the 9th order approximation.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! d - the direction of derivative;
|
||||||
|
! h - the spacial interval;
|
||||||
|
! u - the input scalar field;
|
||||||
|
! v - the first derivative of the input field;
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine derivative_1st_9o(d, h, u, v)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
|
integer :: m0, m1, m2, m3, m4, m5, m6, m7, m8
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
#ifdef PROFILE
|
||||||
|
call start_timer(im1)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
m0 = size(u, d)
|
||||||
|
m1 = m0 - 1
|
||||||
|
m2 = m0 - 2
|
||||||
|
m3 = m0 - 3
|
||||||
|
m4 = m0 - 4
|
||||||
|
m5 = m0 - 5
|
||||||
|
m6 = m0 - 6
|
||||||
|
m7 = m0 - 7
|
||||||
|
m8 = m0 - 8
|
||||||
|
|
||||||
|
select case(d)
|
||||||
|
|
||||||
|
case(1)
|
||||||
|
|
||||||
|
v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h
|
||||||
|
v(2 ,:,:) = 5.00d-01 * (u(3 ,:,:) - u(1 ,:,:)) / h
|
||||||
|
v(3 ,:,:) = (8.00d+00 * (u(4 ,:,:) - u(2 ,:,:)) &
|
||||||
|
- (u(5 ,:,:) - u(1 ,:,:))) / (1.2d+01 * h)
|
||||||
|
v(4 ,:,:) = (4.50d+01 * (u(5 ,:,:) - u(3 ,:,:)) &
|
||||||
|
- 9.00d+00 * (u(6 ,:,:) - u(2 ,:,:)) &
|
||||||
|
+ (u(7 ,:,:) - u(1 ,:,:))) / (6.0d+01 * h)
|
||||||
|
v(5:m4,:,:) = (6.72d+02 * (u(6:m3,:,:) - u(4:m5,:,:)) &
|
||||||
|
- 1.68d+02 * (u(7:m2,:,:) - u(3:m6,:,:)) &
|
||||||
|
+ 3.20d+01 * (u(8:m1,:,:) - u(2:m7,:,:)) &
|
||||||
|
- 3.00d+00 * (u(9:m0,:,:) - u(1:m8,:,:))) / (8.4d+02 * h)
|
||||||
|
v( m3,:,:) = (4.50d+01 * (u( m2,:,:) - u( m4,:,:)) &
|
||||||
|
- 9.00d+00 * (u( m1,:,:) - u( m5,:,:)) &
|
||||||
|
+ (u( m0,:,:) - u( m6,:,:))) / (6.0d+01 * h)
|
||||||
|
v( m2,:,:) = (8.00d+00 * (u( m1,:,:) - u( m3,:,:)) &
|
||||||
|
- (u( m0,:,:) - u( m4,:,:))) / (1.2d+01 * h)
|
||||||
|
v( m1,:,:) = 5.00d-01 * (u( m0,:,:) - u( m2,:,:)) / h
|
||||||
|
v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
|
||||||
|
v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h
|
||||||
|
v(:,2 ,:) = 5.00d-01 * (u(:,3 ,:) - u(:,1 ,:)) / h
|
||||||
|
v(:,3 ,:) = (8.00d+00 * (u(:,4 ,:) - u(:,2 ,:)) &
|
||||||
|
- (u(:,5 ,:) - u(:,1 ,:))) / (1.2d+01 * h)
|
||||||
|
v(:,4 ,:) = (4.50d+01 * (u(:,5 ,:) - u(:,3 ,:)) &
|
||||||
|
- 9.00d+00 * (u(:,6 ,:) - u(:,2 ,:)) &
|
||||||
|
+ (u(:,7 ,:) - u(:,1 ,:))) / (6.0d+01 * h)
|
||||||
|
v(:,5:m4,:) = (6.72d+02 * (u(:,6:m3,:) - u(:,4:m5,:)) &
|
||||||
|
- 1.68d+02 * (u(:,7:m2,:) - u(:,3:m6,:)) &
|
||||||
|
+ 3.20d+01 * (u(:,8:m1,:) - u(:,2:m7,:)) &
|
||||||
|
- 3.00d+00 * (u(:,9:m0,:) - u(:,1:m8,:))) / (8.4d+02 * h)
|
||||||
|
v(:, m3,:) = (4.50d+01 * (u(:, m2,:) - u(:, m4,:)) &
|
||||||
|
- 9.00d+00 * (u(:, m1,:) - u(:, m5,:)) &
|
||||||
|
+ (u(:, m0,:) - u(:, m6,:))) / (6.0d+01 * h)
|
||||||
|
v(:, m2,:) = (8.00d+00 * (u(:, m1,:) - u(:, m3,:)) &
|
||||||
|
- (u(:, m0,:) - u(:, m4,:))) / (1.2d+01 * h)
|
||||||
|
v(:, m1,:) = 5.00d-01 * (u(:, m0,:) - u(:, m2,:)) / h
|
||||||
|
v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
case(3)
|
||||||
|
|
||||||
|
v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h
|
||||||
|
v(:,:,2 ) = 5.00d-01 * (u(:,:,3 ) - u(:,:,1 )) / h
|
||||||
|
v(:,:,3 ) = (8.00d+00 * (u(:,:,4 ) - u(:,:,2 )) &
|
||||||
|
- (u(:,:,5 ) - u(:,:,1 ))) / (1.2d+01 * h)
|
||||||
|
v(:,:,4 ) = (4.50d+01 * (u(:,:,5 ) - u(:,:,3 )) &
|
||||||
|
- 9.00d+00 * (u(:,:,6 ) - u(:,:,2 )) &
|
||||||
|
+ (u(:,:,7 ) - u(:,:,1 ))) / (6.0d+01 * h)
|
||||||
|
v(:,:,5:m4) = (6.72d+02 * (u(:,:,6:m3) - u(:,:,4:m5)) &
|
||||||
|
- 1.68d+02 * (u(:,:,7:m2) - u(:,:,3:m6)) &
|
||||||
|
+ 3.20d+01 * (u(:,:,8:m1) - u(:,:,2:m7)) &
|
||||||
|
- 3.00d+00 * (u(:,:,9:m0) - u(:,:,1:m8))) / (8.4d+02 * h)
|
||||||
|
v(:,:, m3) = (4.50d+01 * (u(:,:, m2) - u(:,:, m4)) &
|
||||||
|
- 9.00d+00 * (u(:,:, m1) - u(:,:, m5)) &
|
||||||
|
+ (u(:,:, m0) - u(:,:, m6))) / (6.0d+01 * h)
|
||||||
|
v(:,:, m2) = (8.00d+00 * (u(:,:, m1) - u(:,:, m3)) &
|
||||||
|
- (u(:,:, m0) - u(:,:, m4))) / (1.2d+01 * h)
|
||||||
|
v(:,:, m1) = 5.00d-01 * (u(:,:, m0) - u(:,:, m2)) / h
|
||||||
|
v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
#ifdef PROFILE
|
||||||
|
call stop_timer(im1)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine derivative_1st_9o
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine DERIVATIVE_2ND_3O:
|
||||||
|
! ----------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculates the second order derivative of the input scalar field
|
||||||
|
! along the given direction with the 3rd order approximation.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! d - the direction of derivative;
|
||||||
|
! h - the spacial interval;
|
||||||
|
! u - the input scalar field;
|
||||||
|
! v - the output scalar field representing the second derivative of u;
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine derivative_2nd_3o(d, h, u, v)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
|
integer :: m0, m1, m2
|
||||||
|
real(kind=8) :: h2
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
#ifdef PROFILE
|
||||||
call start_timer(im2)
|
call start_timer(im2)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
! prepare index limits
|
m0 = size(u, d)
|
||||||
!
|
|
||||||
m0 = size(u, dir)
|
|
||||||
m1 = m0 - 1
|
m1 = m0 - 1
|
||||||
m2 = m0 - 2
|
m2 = m0 - 2
|
||||||
|
|
||||||
! select the direction
|
h2 = h * h
|
||||||
!
|
|
||||||
select case(dir)
|
select case(d)
|
||||||
|
|
||||||
! derivative along the X direction
|
|
||||||
!
|
|
||||||
case(1)
|
case(1)
|
||||||
|
|
||||||
v(2:m1,:,:) = ((u(3:m0,:,:) + u(1:m2,:,:)) - 2.0d+00 * u(2:m1,:,:)) / dh**2
|
|
||||||
v(1 ,:,:) = 0.0d+00
|
v(1 ,:,:) = 0.0d+00
|
||||||
|
v(2:m1,:,:) = ((u(3:m0,:,:) + u(1:m2,:,:)) - 2.0d+00 * u(2:m1,:,:)) / h2
|
||||||
v( m0,:,:) = 0.0d+00
|
v( m0,:,:) = 0.0d+00
|
||||||
|
|
||||||
! derivative along the Y direction
|
|
||||||
!
|
|
||||||
case(2)
|
case(2)
|
||||||
|
|
||||||
v(:,2:m1,:) = ((u(:,3:m0,:) + u(:,1:m2,:)) - 2.0d+00 * u(:,2:m1,:)) / dh**2
|
|
||||||
v(:,1 ,:) = 0.0d+00
|
v(:,1 ,:) = 0.0d+00
|
||||||
|
v(:,2:m1,:) = ((u(:,3:m0,:) + u(:,1:m2,:)) - 2.0d+00 * u(:,2:m1,:)) / h2
|
||||||
v(:, m0,:) = 0.0d+00
|
v(:, m0,:) = 0.0d+00
|
||||||
|
|
||||||
#if NDIMS == 3
|
#if NDIMS == 3
|
||||||
! derivative along the Z direction
|
|
||||||
!
|
|
||||||
case(3)
|
case(3)
|
||||||
|
|
||||||
v(:,:,2:m1) = ((u(:,:,3:m0) + u(:,:,1:m2)) - 2.0d+00 * u(:,:,2:m1)) / dh**2
|
|
||||||
v(:,:,1 ) = 0.0d+00
|
v(:,:,1 ) = 0.0d+00
|
||||||
|
v(:,:,2:m1) = ((u(:,:,3:m0) + u(:,:,1:m2)) - 2.0d+00 * u(:,:,2:m1)) / h2
|
||||||
v(:,:, m0) = 0.0d+00
|
v(:,:, m0) = 0.0d+00
|
||||||
#endif /* NDIMS == 3 */
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
end select
|
end select
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
! stop accounting time for the 2nd derivative calculation
|
|
||||||
!
|
|
||||||
call stop_timer(im2)
|
call stop_timer(im2)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
end subroutine derivative_2nd
|
end subroutine derivative_2nd_3o
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine DERIVATIVE_2ND_5O:
|
||||||
|
! ----------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculates the second order derivative of the input scalar field
|
||||||
|
! along the given direction with the 5th order approximation.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! d - the direction of derivative;
|
||||||
|
! h - the spacial interval;
|
||||||
|
! u - the input scalar field;
|
||||||
|
! v - the output scalar field representing the second derivative of u;
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine derivative_2nd_5o(d, h, u, v)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
|
integer :: m0, m1, m2, m3, m4
|
||||||
|
real(kind=8) :: h2
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
#ifdef PROFILE
|
||||||
|
call start_timer(im2)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
m0 = size(u, d)
|
||||||
|
m1 = m0 - 1
|
||||||
|
m2 = m0 - 2
|
||||||
|
m3 = m0 - 3
|
||||||
|
m4 = m0 - 4
|
||||||
|
|
||||||
|
h2 = h * h
|
||||||
|
|
||||||
|
select case(d)
|
||||||
|
|
||||||
|
case(1)
|
||||||
|
|
||||||
|
v(1 ,:,:) = 0.0d+00
|
||||||
|
v(2 ,:,:) = ((u(3 ,:,:) + u(1 ,:,:)) - 2.0d+00 * u(2 ,:,:)) / h2
|
||||||
|
v(3:m2,:,:) = (1.6d+01 * (u(4:m1,:,:) + u(2:m3,:,:)) &
|
||||||
|
- (u(5:m0,:,:) + u(1:m4,:,:)) &
|
||||||
|
- 3.0d+01 * u(3:m2,:,:)) / (1.2d+01 * h2)
|
||||||
|
v( m1,:,:) = ((u( m0,:,:) + u( m2,:,:)) - 2.0d+00 * u( m1,:,:)) / h2
|
||||||
|
v( m0,:,:) = 0.0d+00
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
|
||||||
|
v(:,1 ,:) = 0.0d+00
|
||||||
|
v(:,2 ,:) = ((u(:,3 ,:) + u(:,1 ,:)) - 2.0d+00 * u(:,2 ,:)) / h2
|
||||||
|
v(:,3:m2,:) = (1.6d+01 * (u(:,4:m1,:) + u(:,2:m3,:)) &
|
||||||
|
- (u(:,5:m0,:) + u(:,1:m4,:)) &
|
||||||
|
- 3.0d+01 * u(:,3:m2,:)) / (1.2d+01 * h2)
|
||||||
|
v(:, m1,:) = ((u(:, m0,:) + u(:, m2,:)) - 2.0d+00 * u(:, m1,:)) / h2
|
||||||
|
v(:, m0,:) = 0.0d+00
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
case(3)
|
||||||
|
|
||||||
|
v(:,:,1 ) = 0.0d+00
|
||||||
|
v(:,:,2 ) = ((u(:,:,3 ) + u(:,:,1 )) - 2.0d+00 * u(:,:,2 )) / h2
|
||||||
|
v(:,:,3:m2) = (1.6d+01 * (u(:,:,4:m1) + u(:,:,2:m3)) &
|
||||||
|
- (u(:,:,5:m0) + u(:,:,1:m4)) &
|
||||||
|
- 3.0d+01 * u(:,:,3:m2)) / (1.2d+01 * h2)
|
||||||
|
v(:,:, m1) = ((u(:,:, m0) + u(:,:, m2)) - 2.0d+00 * u(:,:, m1)) / h2
|
||||||
|
v(:,:, m0) = 0.0d+00
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
#ifdef PROFILE
|
||||||
|
call stop_timer(im2)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine derivative_2nd_5o
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine DERIVATIVE_2ND_7O:
|
||||||
|
! ----------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculates the second order derivative of the input scalar field
|
||||||
|
! along the given direction with the 7th order approximation.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! d - the direction of derivative;
|
||||||
|
! h - the spacial interval;
|
||||||
|
! u - the input scalar field;
|
||||||
|
! v - the output scalar field representing the second derivative of u;
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine derivative_2nd_7o(d, h, u, v)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
|
integer :: m0, m1, m2, m3, m4, m5, m6
|
||||||
|
real(kind=8) :: h2
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
#ifdef PROFILE
|
||||||
|
call start_timer(im2)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
m0 = size(u, d)
|
||||||
|
m1 = m0 - 1
|
||||||
|
m2 = m0 - 2
|
||||||
|
m3 = m0 - 3
|
||||||
|
m4 = m0 - 4
|
||||||
|
m5 = m0 - 5
|
||||||
|
m6 = m0 - 6
|
||||||
|
|
||||||
|
h2 = h * h
|
||||||
|
|
||||||
|
select case(d)
|
||||||
|
|
||||||
|
case(1)
|
||||||
|
|
||||||
|
v(1 ,:,:) = 0.0d+00
|
||||||
|
v(2 ,:,:) = ((u(3 ,:,:) + u(1 ,:,:)) - 2.0d+00 * u(2 ,:,:)) / h2
|
||||||
|
v(3 ,:,:) = (1.60d+01 * (u(4 ,:,:) + u(2 ,:,:)) &
|
||||||
|
- (u(5 ,:,:) + u(1 ,:,:)) &
|
||||||
|
- 3.00d+01 * u(3 ,:,:)) / (1.2d+01 * h2)
|
||||||
|
v(4:m3,:,:) = (1.35d+02 * (u(5:m2,:,:) + u(3:m4,:,:)) &
|
||||||
|
- 1.35d+01 * (u(6:m1,:,:) + u(2:m5,:,:)) &
|
||||||
|
+ (u(7:m0,:,:) + u(1:m6,:,:)) &
|
||||||
|
- 2.45d+02 * u(4:m3,:,:)) / (9.0d+01 * h2)
|
||||||
|
v( m2,:,:) = (1.60d+01 * (u( m1,:,:) + u( m3,:,:)) &
|
||||||
|
- (u( m0,:,:) + u( m4,:,:)) &
|
||||||
|
- 3.00d+01 * u( m2,:,:)) / (1.2d+01 * h2)
|
||||||
|
v( m1,:,:) = ((u( m0,:,:) + u( m2,:,:)) - 2.0d+00 * u( m1,:,:)) / h2
|
||||||
|
v( m0,:,:) = 0.0d+00
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
|
||||||
|
v(:,1 ,:) = 0.0d+00
|
||||||
|
v(:,2 ,:) = ((u(:,3 ,:) + u(:,1 ,:)) - 2.0d+00 * u(:,2 ,:)) / h2
|
||||||
|
v(:,3 ,:) = (1.60d+01 * (u(:,4 ,:) + u(:,2 ,:)) &
|
||||||
|
- (u(:,5 ,:) + u(:,1 ,:)) &
|
||||||
|
- 3.00d+01 * u(:,3 ,:)) / (1.2d+01 * h2)
|
||||||
|
v(:,4:m3,:) = (1.35d+02 * (u(:,5:m2,:) + u(:,3:m4,:)) &
|
||||||
|
- 1.35d+01 * (u(:,6:m1,:) + u(:,2:m5,:)) &
|
||||||
|
+ (u(:,7:m0,:) + u(:,1:m6,:)) &
|
||||||
|
- 2.45d+02 * u(:,4:m3,:)) / (9.0d+01 * h2)
|
||||||
|
v(:, m2,:) = (1.60d+01 * (u(:, m1,:) + u(:, m3,:)) &
|
||||||
|
- (u(:, m0,:) + u(:, m4,:)) &
|
||||||
|
- 3.00d+01 * u(:, m2,:)) / (1.2d+01 * h2)
|
||||||
|
v(:, m1,:) = ((u(:, m0,:) + u(:, m2,:)) - 2.0d+00 * u(:, m1,:)) / h2
|
||||||
|
v(:, m0,:) = 0.0d+00
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
case(3)
|
||||||
|
|
||||||
|
v(:,:,1 ) = 0.0d+00
|
||||||
|
v(:,:,2 ) = ((u(:,:,3 ) + u(:,:,1 )) - 2.0d+00 * u(:,:,2 )) / h2
|
||||||
|
v(:,:,3 ) = (1.60d+01 * (u(:,:,4 ) + u(:,:,2 )) &
|
||||||
|
- (u(:,:,5 ) + u(:,:,1 )) &
|
||||||
|
- 3.00d+01 * u(:,:,3 )) / (1.2d+01 * h2)
|
||||||
|
v(:,:,4:m3) = (1.35d+02 * (u(:,:,5:m2) + u(:,:,3:m4)) &
|
||||||
|
- 1.35d+01 * (u(:,:,6:m1) + u(:,:,2:m5)) &
|
||||||
|
+ (u(:,:,7:m0) + u(:,:,1:m6)) &
|
||||||
|
- 2.45d+02 * u(:,:,4:m3)) / (9.0d+01 * h2)
|
||||||
|
v(:,:, m2) = (1.60d+01 * (u(:,:, m1) + u(:,:, m3)) &
|
||||||
|
- (u(:,:, m0) + u(:,:, m4)) &
|
||||||
|
- 3.00d+01 * u(:,:, m2)) / (1.2d+01 * h2)
|
||||||
|
v(:,:, m1) = ((u(:,:, m0) + u(:,:, m2)) - 2.0d+00 * u(:,:, m1)) / h2
|
||||||
|
v(:,:, m0) = 0.0d+00
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
#ifdef PROFILE
|
||||||
|
call stop_timer(im2)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine derivative_2nd_7o
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine DERIVATIVE_2ND_9O:
|
||||||
|
! ----------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculates the second order derivative of the input scalar field
|
||||||
|
! along the given direction with the 9th order approximation.
|
||||||
|
!
|
||||||
|
! Arguments:
|
||||||
|
!
|
||||||
|
! d - the direction of derivative;
|
||||||
|
! h - the spacial interval;
|
||||||
|
! u - the input scalar field;
|
||||||
|
! v - the output scalar field representing the second derivative of u;
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine derivative_2nd_9o(d, h, u, v)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer , intent(in) :: d
|
||||||
|
real(kind=8) , intent(in) :: h
|
||||||
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
||||||
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
||||||
|
|
||||||
|
integer :: m0, m1, m2, m3, m4, m5, m6, m7, m8
|
||||||
|
real(kind=8) :: h2
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
#ifdef PROFILE
|
||||||
|
call start_timer(im2)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
m0 = size(u, d)
|
||||||
|
m1 = m0 - 1
|
||||||
|
m2 = m0 - 2
|
||||||
|
m3 = m0 - 3
|
||||||
|
m4 = m0 - 4
|
||||||
|
m5 = m0 - 5
|
||||||
|
m6 = m0 - 6
|
||||||
|
m7 = m0 - 7
|
||||||
|
m8 = m0 - 8
|
||||||
|
|
||||||
|
h2 = h * h
|
||||||
|
|
||||||
|
select case(d)
|
||||||
|
|
||||||
|
case(1)
|
||||||
|
|
||||||
|
v(1 ,:,:) = 0.0d+00
|
||||||
|
v(2 ,:,:) = ((u(3 ,:,:) + u(1 ,:,:)) - 2.0d+00 * u(2 ,:,:)) / h2
|
||||||
|
v(3 ,:,:) = (1.600d+01 * (u(4 ,:,:) + u(2 ,:,:)) &
|
||||||
|
- (u(5 ,:,:) + u(1 ,:,:)) &
|
||||||
|
- 3.000d+01 * u(3 ,:,:)) / (1.2d+01 * h2)
|
||||||
|
v(4 ,:,:) = (1.350d+02 * (u(5 ,:,:) + u(3 ,:,:)) &
|
||||||
|
- 1.350d+01 * (u(6 ,:,:) + u(2 ,:,:)) &
|
||||||
|
+ (u(7 ,:,:) + u(1 ,:,:)) &
|
||||||
|
- 2.450d+02 * u(4 ,:,:)) / (9.0d+01 * h2)
|
||||||
|
v(5:m4,:,:) = (8.064d+03 * (u(6:m3,:,:) + u(4:m5,:,:)) &
|
||||||
|
- 1.008d+03 * (u(7:m2,:,:) + u(3:m6,:,:)) &
|
||||||
|
+ 1.280d+02 * (u(8:m1,:,:) + u(2:m7,:,:)) &
|
||||||
|
- 9.000d+00 * (u(9:m0,:,:) + u(1:m8,:,:)) &
|
||||||
|
- 1.435d+04 * u(5:m4,:,:)) / (5.04d+03 * h2)
|
||||||
|
v( m3,:,:) = (1.350d+02 * (u( m2,:,:) + u( m4,:,:)) &
|
||||||
|
- 1.350d+01 * (u( m1,:,:) + u( m5,:,:)) &
|
||||||
|
+ (u( m0,:,:) + u( m6,:,:)) &
|
||||||
|
- 2.450d+02 * u( m3,:,:)) / (9.0d+01 * h2)
|
||||||
|
v( m2,:,:) = (1.600d+01 * (u( m1,:,:) + u( m3,:,:)) &
|
||||||
|
- (u( m0,:,:) + u( m4,:,:)) &
|
||||||
|
- 3.000d+01 * u( m2,:,:)) / (1.2d+01 * h2)
|
||||||
|
v( m1,:,:) = ((u( m0,:,:) + u( m2,:,:)) - 2.0d+00 * u( m1,:,:)) / h2
|
||||||
|
v( m0,:,:) = 0.0d+00
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
|
||||||
|
v(:,1 ,:) = 0.0d+00
|
||||||
|
v(:,2 ,:) = ((u(:,3 ,:) + u(:,1 ,:)) - 2.0d+00 * u(:,2 ,:)) / h2
|
||||||
|
v(:,3 ,:) = (1.600d+01 * (u(:,4 ,:) + u(:,2 ,:)) &
|
||||||
|
- (u(:,5 ,:) + u(:,1 ,:)) &
|
||||||
|
- 3.000d+01 * u(:,3 ,:)) / (1.2d+01 * h2)
|
||||||
|
v(:,4 ,:) = (1.350d+02 * (u(:,5 ,:) + u(:,3 ,:)) &
|
||||||
|
- 1.350d+01 * (u(:,6 ,:) + u(:,2 ,:)) &
|
||||||
|
+ (u(:,7 ,:) + u(:,1 ,:)) &
|
||||||
|
- 2.450d+02 * u(:,4 ,:)) / (9.0d+01 * h2)
|
||||||
|
v(:,5:m4,:) = (8.064d+03 * (u(:,6:m3,:) + u(:,4:m5,:)) &
|
||||||
|
- 1.008d+03 * (u(:,7:m2,:) + u(:,3:m6,:)) &
|
||||||
|
+ 1.280d+02 * (u(:,8:m1,:) + u(:,2:m7,:)) &
|
||||||
|
- 9.000d+00 * (u(:,9:m0,:) + u(:,1:m8,:)) &
|
||||||
|
- 1.435d+04 * u(:,5:m4,:)) / (5.04d+03 * h2)
|
||||||
|
v(:, m3,:) = (1.350d+02 * (u(:, m2,:) + u(:, m4,:)) &
|
||||||
|
- 1.350d+01 * (u(:, m1,:) + u(:, m5,:)) &
|
||||||
|
+ (u(:, m0,:) + u(:, m6,:)) &
|
||||||
|
- 2.450d+02 * u(:, m3,:)) / (9.0d+01 * h2)
|
||||||
|
v(:, m2,:) = (1.600d+01 * (u(:, m1,:) + u(:, m3,:)) &
|
||||||
|
- (u(:, m0,:) + u(:, m4,:)) &
|
||||||
|
- 3.000d+01 * u(:, m2,:)) / (1.2d+01 * h2)
|
||||||
|
v(:, m1,:) = ((u(:, m0,:) + u(:, m2,:)) - 2.0d+00 * u(:, m1,:)) / h2
|
||||||
|
v(:, m0,:) = 0.0d+00
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
case(3)
|
||||||
|
|
||||||
|
v(:,:,1 ) = 0.0d+00
|
||||||
|
v(:,:,2 ) = ((u(:,:,3 ) + u(:,:,1 )) - 2.0d+00 * u(:,:,2 )) / h2
|
||||||
|
v(:,:,3 ) = (1.600d+01 * (u(:,:,4 ) + u(:,:,2 )) &
|
||||||
|
- (u(:,:,5 ) + u(:,:,1 )) &
|
||||||
|
- 3.000d+01 * u(:,:,3 )) / (1.2d+01 * h2)
|
||||||
|
v(:,:,4 ) = (1.350d+02 * (u(:,:,5 ) + u(:,:,3 )) &
|
||||||
|
- 1.350d+01 * (u(:,:,6 ) + u(:,:,2 )) &
|
||||||
|
+ (u(:,:,7 ) + u(:,:,1 )) &
|
||||||
|
- 2.450d+02 * u(:,:,4 )) / (9.0d+01 * h2)
|
||||||
|
v(:,:,5:m4) = (8.064d+03 * (u(:,:,6:m3) + u(:,:,4:m5)) &
|
||||||
|
- 1.008d+03 * (u(:,:,7:m2) + u(:,:,3:m6)) &
|
||||||
|
+ 1.280d+02 * (u(:,:,8:m1) + u(:,:,2:m7)) &
|
||||||
|
- 9.000d+00 * (u(:,:,9:m0) + u(:,:,1:m8)) &
|
||||||
|
- 1.435d+04 * u(:,:,5:m4)) / (5.04d+03 * h2)
|
||||||
|
v(:,:, m3) = (1.350d+02 * (u(:,:, m2) + u(:,:, m4)) &
|
||||||
|
- 1.350d+01 * (u(:,:, m1) + u(:,:, m5)) &
|
||||||
|
+ (u(:,:, m0) + u(:,:, m6)) &
|
||||||
|
- 2.450d+02 * u(:,:, m3)) / (9.0d+01 * h2)
|
||||||
|
v(:,:, m2) = (1.600d+01 * (u(:,:, m1) + u(:,:, m3)) &
|
||||||
|
- (u(:,:, m0) + u(:,:, m4)) &
|
||||||
|
- 3.000d+01 * u(:,:, m2)) / (1.2d+01 * h2)
|
||||||
|
v(:,:, m1) = ((u(:,:, m0) + u(:,:, m2)) - 2.0d+00 * u(:,:, m1)) / h2
|
||||||
|
v(:,:, m0) = 0.0d+00
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
#ifdef PROFILE
|
||||||
|
call stop_timer(im2)
|
||||||
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine derivative_2nd_9o
|
||||||
|
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user