Newer
Older
Antoine Cyril David Hoffmann
committed
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
MODULE calculus
! Routine to evaluate gradients, interpolation schemes and integrals
USE basic
USE prec_const
USE grid
IMPLICIT NONE
REAL(dp), dimension(-2:2) :: dz_usu = &
(/ onetwelfth, -twothird, &
0._dp, &
twothird, -onetwelfth /) ! fd4 centered stencil
REAL(dp), dimension(-2:1) :: dz_o2e = &
(/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
REAL(dp), dimension(-1:2) :: dz_e2o = &
(/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
PUBLIC :: simpson_rule_z, interp_z, grad_z
CONTAINS
SUBROUTINE grad_z(target,f,ddzf)
implicit none
! Compute the periodic boundary condition 4 points centered finite differences
! formula among staggered grid or not.
! not staggered : the derivative results must be on the same grid as the field
! staggered : the derivative is computed from a grid to the other
INTEGER, INTENT(IN) :: target
COMPLEX(dp),dimension( izs:ize ), intent(in) :: f
COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzf
IF(SG) THEN
IF(TARGET .EQ. 0) THEN
CALL grad_z_o2e(f,ddzf)
ELSE
CALL grad_z_e2o(f,ddzf)
ENDIF
ELSE ! No staggered grid -> usual centered finite differences
DO iz = 3,Nz-2
ddzf(iz) = dz_usu(-2)*f(iz-2) + dz_usu(-1)*f(iz-1) &
+dz_usu( 1)*f(iz+1) + dz_usu( 2)*f(iz+2)
ENDDO
! Periodic boundary conditions
ddzf(Nz-1) = dz_usu(-2)*f(Nz-3) + dz_usu(-1)*f(Nz-2) &
+dz_usu(+1)*f(Nz ) + dz_usu(+2)*f(1 )
ddzf(Nz ) = dz_usu(-2)*f(Nz-2) + dz_usu(-1)*f(Nz-1) &
+dz_usu(+1)*f(1 ) + dz_usu(+2)*f(2 )
ddzf(1 ) = dz_usu(-2)*f(Nz-1) + dz_usu(-1)*f(Nz ) &
+dz_usu(+1)*f(2 ) + dz_usu(+2)*f(3 )
ddzf(2 ) = dz_usu(-2)*f(Nz ) + dz_usu(-1)*f(1 ) &
+dz_usu(+1)*f(3 ) + dz_usu(+2)*f(4)
ENDIF
END SUBROUTINE grad_z
SUBROUTINE grad_z_o2e(fo,dzfe) ! Paruta 2018 eq (27)
! gives the value of a field from the odd grid to the even one
implicit none
COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
COMPLEX(dp),dimension(1:Nz), intent(out) :: dzfe !
DO iz = 3,Nz-1
dzfe(iz) = dz_o2e(-2)*fo(iz-2) + dz_o2e(-1)*fo(iz-1) &
+dz_o2e( 0)*fo(iz ) + dz_o2e( 1)*fo(iz+1)
ENDDO
! Periodic boundary conditions
dzfe(Nz) = dz_o2e(-2)*fo(Nz-2) + dz_o2e(-1)*fo(Nz-1) &
+dz_o2e( 0)*fo(Nz ) + dz_o2e( 1)*fo(1)
dzfe(1 ) = dz_o2e(-2)*fo(Nz-1) + dz_o2e(-1)*fo(Nz ) &
+dz_o2e( 0)*fo(1 ) + dz_o2e( 1)*fo(2)
dzfe(2 ) = dz_o2e(-2)*fo(Nz ) + dz_o2e(-1)*fo(1 ) &
+dz_o2e( 0)*fo(2 ) + dz_o2e( 1)*fo(3)
END SUBROUTINE grad_z_o2e
SUBROUTINE grad_z_e2o(fe,dzfo)
! gives the value of a field from the even grid to the odd one
implicit none
COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
COMPLEX(dp),dimension(1:Nz), intent(out) :: dzfo
DO iz = 2,Nz-2
dzfo(iz) = dz_e2o(-1)*fe(iz-1) + dz_e2o(0)*fe(iz ) &
+dz_e2o( 1)*fe(iz+1) + dz_e2o(2)*fe(iz+2)
ENDDO
! Periodic boundary conditions
dzfo(Nz-1) = dz_e2o(-1)*fe(Nz-2) + dz_e2o(0)*fe(Nz-1) &
+dz_e2o( 1)*fe(Nz ) + dz_e2o(2)*fe(1 )
dzfo(Nz ) = dz_e2o(-1)*fe(Nz-1) + dz_e2o(0)*fe(Nz ) &
+dz_e2o( 1)*fe(1 ) + dz_e2o(2)*fe(2 )
dzfo(1 ) = dz_e2o(-1)*fe(Nz ) + dz_e2o(0)*fe(1 ) &
+dz_e2o( 1)*fe(2 ) + dz_e2o(2)*fe(3 )
END SUBROUTINE grad_z_e2o
SUBROUTINE interp_z(target,f_in,f_out)
! Function meant to interpolate one field defined on a even/odd z into
! the other odd/even z grid.
! If Staggered Grid flag (SG) is false, returns identity
implicit none
INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd
COMPLEX(dp),dimension(1:Nz), intent(in) :: f_in
COMPLEX(dp),dimension(1:Nz), intent(out) :: f_out !
IF(SG) THEN
IF(TARGET .EQ. 0) THEN
CALL interp_o2e_z(f_in,f_out)
ELSE
CALL interp_e2o_z(f_in,f_out)
ENDIF
ELSE ! No staggered grid -> identity
f_out(:) = f_in(:)
ENDIF
END SUBROUTINE interp_z
SUBROUTINE interp_o2e_z(fo,fe)
! gives the value of a field from the odd grid to the even one
implicit none
COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
COMPLEX(dp),dimension(1:Nz), intent(out) :: fe !
DO iz = 2,Nz
fe(iz) = 0.5_dp*(fo(iz)+fo(iz-1))
ENDDO
! Periodic boundary conditions
fe(1) = 0.5_dp*(fo(1) + fo(Nz))
END SUBROUTINE interp_o2e_z
SUBROUTINE interp_e2o_z(fe,fo)
! gives the value of a field from the even grid to the odd one
implicit none
COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
COMPLEX(dp),dimension(1:Nz), intent(out) :: fo
DO iz = 1,Nz-1
fo(iz) = 0.5_dp*(fe(iz+1)+fe(iz))
ENDDO
! Periodic boundary conditions
fo(Nz) = 0.5_dp*(fe(1) + fe(Nz))
END SUBROUTINE interp_e2o_z
SUBROUTINE simpson_rule_z(f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei
implicit none
complex(dp),dimension(izs:ize), intent(in) :: f
COMPLEX(dp), intent(out) :: intf
COMPLEX(dp) :: buff_
IF(Nz .GT. 1) THEN
IF(mod(Nz,2) .ne. 0 ) THEN
ERROR STOP 'Simpson rule: Nz must be an even number !!!!'
ENDIF
buff_ = 0._dp
DO iz = izs, Nz/2
IF(iz .eq. Nz/2) THEN ! ... iz = ize
buff_ = buff_ + (f(izs) + 4._dp*f(ize) + f(ize-1 ))
ELSE
buff_ = buff_ + (f(2*iz+1) + 4._dp*f(2*iz) + f(2*iz-1 ))
ENDIF
ENDDO
intf = buff_*deltaz/3._dp
ELSE
intf = f(izs)
ENDIF
END SUBROUTINE simpson_rule_z
END MODULE calculus