1121 lines
39 KiB
Fortran
1121 lines
39 KiB
Fortran
!!******************************************************************************
|
|
!!
|
|
!! This file is part of the AMUN source code, a program to perform
|
|
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
|
!! adaptive mesh.
|
|
!!
|
|
!! Copyright (C) 2008-2024 Grzegorz Kowal <grzegorz@amuncode.org>
|
|
!!
|
|
!! This program is free software: you can redistribute it and/or modify
|
|
!! it under the terms of the GNU General Public License as published by
|
|
!! the Free Software Foundation, either version 3 of the License, or
|
|
!! (at your option) any later version.
|
|
!!
|
|
!! This program is distributed in the hope that it will be useful,
|
|
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
!! GNU General Public License for more details.
|
|
!!
|
|
!! You should have received a copy of the GNU General Public License
|
|
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
!!
|
|
!!*****************************************************************************
|
|
!!
|
|
!! module: OPERATORS
|
|
!!
|
|
!! This module provides differential operators like gradient, divergence, or
|
|
!! curl.
|
|
!!
|
|
!!*****************************************************************************
|
|
!
|
|
module operators
|
|
|
|
implicit none
|
|
|
|
! 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()
|
|
|
|
private
|
|
|
|
public :: initialize_operators, finalize_operators
|
|
public :: divergence, gradient, curl, laplace
|
|
public :: derivative_1st, derivative_2nd
|
|
|
|
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
!
|
|
contains
|
|
!
|
|
!===============================================================================
|
|
!!
|
|
!!*** PUBLIC SUBROUTINES *****************************************************
|
|
!!
|
|
!===============================================================================
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine INITIALIZE_OPERATORS:
|
|
! -------------------------------
|
|
!
|
|
! Subroutine initializes the module structures, pointers and variables.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! order - the order of the interpolation method;
|
|
! status - return flag of the procedure execution status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine initialize_operators(order, status)
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: order
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine initialize_operators
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine FINALIZE_OPERATORS:
|
|
! -----------------------------
|
|
!
|
|
! Subroutine releases the memory used by module variables and pointers.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - return flag of the procedure execution status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine finalize_operators(status)
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine finalize_operators
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine DIVERGENCE:
|
|
! ---------------------
|
|
!
|
|
! Subroutine calculates the cell centered divergence of the input vector
|
|
! field.
|
|
!
|
|
! div(U) = ∇.([Ux, Uy, Uz]) = ∂x Ux + ∂y Uy + ∂z Uz
|
|
!
|
|
! Arguments:
|
|
!
|
|
! dh - the spacial intervals in all direction;
|
|
! u - the input vector field;
|
|
! v - the output divergence field;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine divergence(dh, u, v)
|
|
|
|
implicit none
|
|
|
|
real(kind=8), dimension(3) , intent(in) :: dh
|
|
real(kind=8), dimension(:,:,:,:), intent(in) :: u
|
|
real(kind=8), dimension(:,:,:) , intent(out) :: v
|
|
|
|
integer :: dir
|
|
|
|
real(kind=8), dimension(size(u,2), size(u,3), size(u,4)) :: w
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
v(:,:,:) = 0.0d+00
|
|
|
|
! iterate over directions and update divergence with directional derivatives
|
|
!
|
|
do dir = 1, NDIMS
|
|
|
|
! calculate derivative along the current direction
|
|
!
|
|
call derivative_1st(dir, dh(dir), u(dir,:,:,:), w(:,:,:))
|
|
|
|
! update the divergence array
|
|
!
|
|
v(:,:,:) = v(:,:,:) + w(:,:,:)
|
|
|
|
end do ! directions
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine divergence
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine GRADIENT:
|
|
! -------------------
|
|
!
|
|
! Subroutine calculates the cell centered gradient of the input scalar field.
|
|
!
|
|
! grad(U) = ∇ U = [ ∂x U, ∂y U, ∂z U ]
|
|
!
|
|
! Arguments:
|
|
!
|
|
! dh - the spacial intervals in all direction;
|
|
! u - the input scalar field;
|
|
! v - the output gradient vector field;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine gradient(dh, u, v)
|
|
|
|
implicit none
|
|
|
|
real(kind=8), dimension(3) , intent(in) :: dh
|
|
real(kind=8), dimension(:,:,:) , intent(in) :: u
|
|
real(kind=8), dimension(:,:,:,:), intent(out) :: v
|
|
|
|
integer :: dir
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
v(:,:,:,:) = 0.0d+00
|
|
|
|
! iterate over directions and calculate gradient components
|
|
!
|
|
do dir = 1, NDIMS
|
|
|
|
! calculate derivative along the current direction
|
|
!
|
|
call derivative_1st(dir, dh(dir), u(:,:,:), v(dir,:,:,:))
|
|
|
|
end do ! directions
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine gradient
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine CURL:
|
|
! ---------------
|
|
!
|
|
! Subroutine calculates the cell centered curl of the input vector field.
|
|
!
|
|
! curl(U) = ∇x([Ux, Uy, Uz])
|
|
! = [∂y Uz - ∂z Uy, ∂z Ux - ∂x Uz, ∂x Uy - ∂y Ux]
|
|
!
|
|
! Arguments:
|
|
!
|
|
! dh - the spacial intervals in all direction;
|
|
! u - the input vector field;
|
|
! v - the output divergence field;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine curl(dh, u, v)
|
|
|
|
implicit none
|
|
|
|
real(kind=8), dimension(3) , intent(in) :: dh
|
|
real(kind=8), dimension(:,:,:,:), intent(in) :: u
|
|
real(kind=8), dimension(:,:,:,:), intent(out) :: v
|
|
|
|
real(kind=8), dimension(size(u,2), size(u,3), size(u,4)) :: w
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! === calculate Vx component ===
|
|
!
|
|
! contribution from the Y derivative of Uz
|
|
!
|
|
call derivative_1st(2, dh(2), u(3,:,:,:), w(:,:,:))
|
|
|
|
! update Vx
|
|
!
|
|
v(1,:,:,:) = w(:,:,:)
|
|
|
|
#if NDIMS == 3
|
|
! contribution from the Z derivative of Uy
|
|
!
|
|
call derivative_1st(3, dh(3), u(2,:,:,:), w(:,:,:))
|
|
|
|
! update Vx
|
|
!
|
|
v(1,:,:,:) = v(1,:,:,:) - w(:,:,:)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! === calculate Vy component ===
|
|
!
|
|
#if NDIMS == 3
|
|
! contribution from the Z derivative of Ux
|
|
!
|
|
call derivative_1st(3, dh(3), u(1,:,:,:), w(:,:,:))
|
|
|
|
! update Vy
|
|
!
|
|
v(2,:,:,:) = w(:,:,:)
|
|
|
|
! contribution from the X derivative of Uz
|
|
!
|
|
call derivative_1st(1, dh(1), u(3,:,:,:), w(:,:,:))
|
|
|
|
! update Vy
|
|
!
|
|
v(2,:,:,:) = v(2,:,:,:) - w(:,:,:)
|
|
#else /* NDIMS == 3 */
|
|
! contribution from the X derivative of Az
|
|
!
|
|
call derivative_1st(1, dh(1), u(3,:,:,:), w(:,:,:))
|
|
|
|
! update Vy
|
|
!
|
|
v(2,:,:,:) = - w(:,:,:)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! === calculate Vz component ===
|
|
!
|
|
! contribution from the X derivative of Uy
|
|
!
|
|
call derivative_1st(1, dh(1), u(2,:,:,:), w(:,:,:))
|
|
|
|
! update Vz
|
|
!
|
|
v(3,:,:,:) = w(:,:,:)
|
|
|
|
! contribution from the Y derivative of Ux
|
|
!
|
|
call derivative_1st(2, dh(2), u(1,:,:,:), w(:,:,:))
|
|
|
|
! update Vz
|
|
!
|
|
v(3,:,:,:) = v(3,:,:,:) - w(:,:,:)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine curl
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine LAPLACE:
|
|
! ------------------
|
|
!
|
|
! Subroutine calculates the Laplace operator of the input scalar field.
|
|
!
|
|
! laplace(U) = Δ(U) = ∂²x U + ∂²y U + ∂²z U
|
|
!
|
|
! Arguments:
|
|
!
|
|
! dh - the spacial intervals in all direction;
|
|
! u - the input scalar field U;
|
|
! v - the output field representing the laplacian of u;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine laplace(dh, u, v)
|
|
|
|
implicit none
|
|
|
|
real(kind=8), dimension(3) , intent(in) :: dh
|
|
real(kind=8), dimension(:,:,:), intent(in) :: u
|
|
real(kind=8), dimension(:,:,:), intent(out) :: v
|
|
|
|
real(kind=8), dimension(size(u,1), size(u,2), size(u,3)) :: w
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! calculate the second derivative of U along the X direction
|
|
!
|
|
call derivative_2nd(1, dh(1), u(:,:,:), w(:,:,:))
|
|
|
|
! update V
|
|
!
|
|
v(:,:,:) = w(:,:,:)
|
|
|
|
! calculate the second derivative of U along the X direction
|
|
!
|
|
call derivative_2nd(2, dh(2), u(:,:,:), w(:,:,:))
|
|
|
|
! update V
|
|
!
|
|
v(:,:,:) = v(:,:,:) + w(:,:,:)
|
|
|
|
#if NDIMS == 3
|
|
! calculate the second derivative of U along the Z direction
|
|
!
|
|
call derivative_2nd(3, dh(3), u(:,:,:), w(:,:,:))
|
|
|
|
! update V
|
|
!
|
|
v(:,:,:) = v(:,:,:) + w(:,:,:)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine laplace
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine DERIVATIVE_1ST_3O:
|
|
! ----------------------------
|
|
!
|
|
! Subroutine calculates the first 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 first derivative of the input field;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine derivative_1st_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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
m0 = size(u, d)
|
|
m1 = m0 - 1
|
|
m2 = m0 - 2
|
|
|
|
select case(d)
|
|
|
|
case(1)
|
|
|
|
v(1 ,:,:) = (u(2 ,:,:) - u(1 ,:,:)) / h
|
|
v(2:m1,:,:) = 5.0d-01 * (u(3:m0,:,:) - u(1:m2,:,:)) / h
|
|
v( m0,:,:) = (u( m0,:,:) - u( m1,:,:)) / h
|
|
|
|
case(2)
|
|
|
|
v(:,1 ,:) = (u(:,2 ,:) - u(:,1 ,:)) / h
|
|
v(:,2:m1,:) = 5.0d-01 * (u(:,3:m0,:) - u(:,1:m2,:)) / h
|
|
v(:, m0,:) = (u(:, m0,:) - u(:, m1,:)) / h
|
|
|
|
#if NDIMS == 3
|
|
case(3)
|
|
|
|
v(:,:,1 ) = (u(:,:,2 ) - u(:,:,1 )) / h
|
|
v(:,:,2:m1) = 5.0d-01 * (u(:,:,3:m0) - u(:,:,1:m2)) / h
|
|
v(:,:, m0) = (u(:,:, m0) - u(:,:, m1)) / h
|
|
#endif /* NDIMS == 3 */
|
|
|
|
end select
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine derivative_1st_3o
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine DERIVATIVE_1ST_5O:
|
|
! ----------------------------
|
|
!
|
|
! Subroutine calculates the first 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 first derivative of the input field;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine derivative_1st_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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
m0 = size(u, d)
|
|
m1 = m0 - 1
|
|
m2 = m0 - 2
|
|
|
|
h2 = h * h
|
|
|
|
select case(d)
|
|
|
|
case(1)
|
|
|
|
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
|
|
|
|
case(2)
|
|
|
|
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
|
|
|
|
#if NDIMS == 3
|
|
case(3)
|
|
|
|
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
|
|
#endif /* NDIMS == 3 */
|
|
|
|
end select
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine derivative_2nd_9o
|
|
|
|
!===============================================================================
|
|
!
|
|
end module operators
|