1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include "cctk.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
subroutine deriv_gf_2_1 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar )
implicit none
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
CCTK_REAL, parameter :: zero = 0.0
integer, parameter :: wp = kind(zero)
CCTK_INT, intent(IN) :: ni, nj, nk
CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var
CCTK_INT, intent(IN) :: dir
CCTK_INT, intent(IN) :: bb(2)
CCTK_INT, intent(IN) :: gsize
CCTK_REAL, intent(IN) :: delta
CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar
CCTK_REAL, dimension(1), save :: a
CCTK_REAL, dimension(3,2), save :: q
CCTK_REAL :: idel
CCTK_INT :: il, ir, jl, jr, kl, kr
logical, save :: first = .true.
if ( first ) then
a(1) = 1.0_wp/2.0_wp
q(1,1) = -1.0_wp; q(2,1) = 1.0_wp; q(3,1) = zero
q(1,2) = -1.0_wp/2.0_wp; q(2,2) = zero; q(3,2) = 1.0_wp/2.0_wp
first = .false.
end if
idel = 1.0_wp / delta
direction: select case (dir)
case (0) direction
if ( bb(1) == 0 ) then
il = 1 + gsize
else
dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) ) * idel
dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) ) * idel
il = 3
end if
if ( bb(2) == 0 ) then
ir = ni - gsize
else
dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + &
q(3,2) * var(ni-2,:,:) ) * idel
dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + &
q(2,1) * var(ni-1,:,:) ) * idel
ir = ni - 2
end if
dvar(il:ir,:,:) = a(1) * ( var(il+1:ir+1,:,:) - var(il-1:ir-1,:,:) ) * idel
case (1) direction
if ( zero_derivs_y /= 0 ) then
dvar = zero
else
if ( bb(1) == 0 ) then
jl = 1 + gsize
else
dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) ) * idel
dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) ) * idel
jl = 3
end if
if ( bb(2) == 0 ) then
jr = nj - gsize
else
dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + &
q(3,2) * var(:,nj-2,:) ) * idel
dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + &
q(2,1) * var(:,nj-1,:) ) * idel
jr = nj - 2
end if
dvar(:,jl:jr,:) = a(1) * ( var(:,jl+1:jr+1,:) - &
var(:,jl-1:jr-1,:) ) * idel
end if
case (2) direction
if ( zero_derivs_z /= 0 ) then
dvar = zero
else
if ( bb(1) == 0 ) then
kl = 1 + gsize
else
dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) ) * idel
dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) ) * idel
kl = 3
end if
if ( bb(2) == 0 ) then
kr = nk - gsize
else
dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + &
q(3,2) * var(:,:,nk-2) ) * idel
dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + &
q(2,1) * var(:,:,nk-1) ) * idel
kr = nk - 2
end if
dvar(:,:,kl:kr) = a(1) * ( var(:,:,kl+1:kr+1) - &
var(:,:,kl-1:kr-1) ) * idel
end if
end select direction
end subroutine deriv_gf_2_1
|
