1 |
|
|
MODULE slopes_m |
2 |
|
|
|
3 |
|
|
! Author: Lionel GUEZ |
4 |
|
|
! Extension / factorisation: David CUGNET |
5 |
|
|
|
6 |
|
|
IMPLICIT NONE |
7 |
|
|
|
8 |
|
|
! Those generic function computes second order slopes with Van |
9 |
|
|
! Leer slope-limiting, given cell averages. Reference: Dukowicz, |
10 |
|
|
! 1987, SIAM Journal on Scientific and Statistical Computing, 8, |
11 |
|
|
! 305. |
12 |
|
|
|
13 |
|
|
! The only difference between the specific functions is the rank |
14 |
|
|
! of the first argument and the equal rank of the result. |
15 |
|
|
|
16 |
|
|
! slope(ix,...) acts on ix th dimension. |
17 |
|
|
|
18 |
|
|
! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1 |
19 |
|
|
! real, intent(in):: x(:) ! (n + 1) cell edges |
20 |
|
|
! real slopes, same shape as f ! (n, ...) |
21 |
|
|
INTERFACE slopes |
22 |
|
|
MODULE procedure slopes1, slopes2, slopes3, slopes4, slopes5 |
23 |
|
|
END INTERFACE |
24 |
|
|
|
25 |
|
|
PRIVATE |
26 |
|
|
PUBLIC :: slopes |
27 |
|
|
|
28 |
|
|
CONTAINS |
29 |
|
|
|
30 |
|
|
!------------------------------------------------------------------------------- |
31 |
|
|
! |
32 |
|
|
PURE FUNCTION slopes1(ix, f, x) |
33 |
|
|
! |
34 |
|
|
!------------------------------------------------------------------------------- |
35 |
|
|
! Arguments: |
36 |
|
|
INTEGER, INTENT(IN) :: ix |
37 |
|
|
REAL, INTENT(IN) :: f(:) |
38 |
|
|
REAL, INTENT(IN) :: x(:) |
39 |
|
|
REAL :: slopes1(SIZE(f,1)) |
40 |
|
|
!------------------------------------------------------------------------------- |
41 |
|
|
! Local: |
42 |
|
|
INTEGER :: n, i, j, sta(2), sto(2) |
43 |
|
|
REAL :: xc(SIZE(f,1)) ! (n) cell centers |
44 |
|
|
REAL :: h(2:SIZE(f,1)-1), delta_xc(2:SIZE(f,1)-1) ! (2:n-1) |
45 |
|
|
REAL :: fm, ff, fp, dx |
46 |
|
|
!------------------------------------------------------------------------------- |
47 |
|
|
n=SIZE(f,ix) |
48 |
|
|
FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. |
49 |
|
|
FORALL(i=2:n-1) |
50 |
|
|
h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) |
51 |
|
|
END FORALL |
52 |
|
|
slopes1(:)=0. |
53 |
|
|
DO i=2,n-1 |
54 |
|
|
ff=f(i); fm=f(i-1); fp=f(i+1) |
55 |
|
|
IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN |
56 |
|
|
slopes1(i)=0.; CYCLE !--- Local extremum |
57 |
|
|
!--- 2nd order slope ; (fm, ff, fp) strictly monotonous |
58 |
|
|
slopes1(i)=(fp-fm)/delta_xc(i) |
59 |
|
|
!--- Slope limitation |
60 |
|
|
slopes1(i) = SIGN(MIN(ABS(slopes1(i)), & |
61 |
|
|
ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes1(i) ) |
62 |
|
|
END IF |
63 |
|
|
END DO |
64 |
|
|
|
65 |
|
|
END FUNCTION slopes1 |
66 |
|
|
! |
67 |
|
|
!------------------------------------------------------------------------------- |
68 |
|
|
|
69 |
|
|
|
70 |
|
|
!------------------------------------------------------------------------------- |
71 |
|
|
! |
72 |
|
|
PURE FUNCTION slopes2(ix, f, x) |
73 |
|
|
! |
74 |
|
|
!------------------------------------------------------------------------------- |
75 |
|
|
! Arguments: |
76 |
|
|
INTEGER, INTENT(IN) :: ix |
77 |
|
|
REAL, INTENT(IN) :: f(:, :) |
78 |
|
|
REAL, INTENT(IN) :: x(:) |
79 |
|
|
REAL :: slopes2(SIZE(f,1),SIZE(f,2)) |
80 |
|
|
!------------------------------------------------------------------------------- |
81 |
|
|
! Local: |
82 |
|
|
INTEGER :: n, i, j, sta(2), sto(2) |
83 |
|
|
REAL, ALLOCATABLE :: xc(:) ! (n) cell centers |
84 |
|
|
REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) |
85 |
|
|
REAL :: fm, ff, fp, dx |
86 |
|
|
!------------------------------------------------------------------------------- |
87 |
|
|
n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) |
88 |
|
|
FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. |
89 |
|
|
FORALL(i=2:n-1) |
90 |
|
|
h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) |
91 |
|
|
END FORALL |
92 |
|
|
slopes2(:,:)=0. |
93 |
|
|
sta=[1,1]; sta(ix)=2 |
94 |
|
|
sto=SHAPE(f); sto(ix)=n-1 |
95 |
|
|
DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) |
96 |
|
|
DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) |
97 |
|
|
ff=f(i,j) |
98 |
|
|
SELECT CASE(ix) |
99 |
|
|
CASE(1); fm=f(i-1,j); fp=f(i+1,j) |
100 |
|
|
CASE(2); fm=f(i,j-1); fp=f(i,j+1) |
101 |
|
|
END SELECT |
102 |
|
|
IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN |
103 |
|
|
slopes2(i,j)=0.; CYCLE !--- Local extremum |
104 |
|
|
!--- 2nd order slope ; (fm, ff, fp) strictly monotonous |
105 |
|
|
slopes2(i,j)=(fp-fm)/dx |
106 |
|
|
!--- Slope limitation |
107 |
|
|
slopes2(i,j) = SIGN(MIN(ABS(slopes2(i,j)), & |
108 |
|
|
ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes2(i,j) ) |
109 |
|
|
END IF |
110 |
|
|
END DO |
111 |
|
|
END DO |
112 |
|
|
DEALLOCATE(xc,h,delta_xc) |
113 |
|
|
|
114 |
|
|
END FUNCTION slopes2 |
115 |
|
|
! |
116 |
|
|
!------------------------------------------------------------------------------- |
117 |
|
|
|
118 |
|
|
|
119 |
|
|
!------------------------------------------------------------------------------- |
120 |
|
|
! |
121 |
|
|
PURE FUNCTION slopes3(ix, f, x) |
122 |
|
|
! |
123 |
|
|
!------------------------------------------------------------------------------- |
124 |
|
|
! Arguments: |
125 |
|
|
INTEGER, INTENT(IN) :: ix |
126 |
|
|
REAL, INTENT(IN) :: f(:, :, :) |
127 |
|
|
REAL, INTENT(IN) :: x(:) |
128 |
|
|
REAL :: slopes3(SIZE(f,1),SIZE(f,2),SIZE(f,3)) |
129 |
|
|
!------------------------------------------------------------------------------- |
130 |
|
|
! Local: |
131 |
|
|
INTEGER :: n, i, j, k, sta(3), sto(3) |
132 |
|
|
REAL, ALLOCATABLE :: xc(:) ! (n) cell centers |
133 |
|
|
REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) |
134 |
|
|
REAL :: fm, ff, fp, dx |
135 |
|
|
!------------------------------------------------------------------------------- |
136 |
|
|
n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) |
137 |
|
|
FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. |
138 |
|
|
FORALL(i=2:n-1) |
139 |
|
|
h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) |
140 |
|
|
END FORALL |
141 |
|
|
slopes3(:,:,:)=0. |
142 |
|
|
sta=[1,1,1]; sta(ix)=2 |
143 |
|
|
sto=SHAPE(f); sto(ix)=n-1 |
144 |
|
|
DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) |
145 |
|
|
DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) |
146 |
|
|
DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) |
147 |
|
|
ff=f(i,j,k) |
148 |
|
|
SELECT CASE(ix) |
149 |
|
|
CASE(1); fm=f(i-1,j,k); fp=f(i+1,j,k) |
150 |
|
|
CASE(2); fm=f(i,j-1,k); fp=f(i,j+1,k) |
151 |
|
|
CASE(3); fm=f(i,j,k-1); fp=f(i,j,k+1) |
152 |
|
|
END SELECT |
153 |
|
|
IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN |
154 |
|
|
slopes3(i,j,k)=0.; CYCLE !--- Local extremum |
155 |
|
|
!--- 2nd order slope ; (fm, ff, fp) strictly monotonous |
156 |
|
|
slopes3(i,j,k)=(fp-fm)/dx |
157 |
|
|
!--- Slope limitation |
158 |
|
|
slopes3(i,j,k) = SIGN(MIN(ABS(slopes3(i,j,k)), & |
159 |
|
|
ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes3(i,j,k) ) |
160 |
|
|
END IF |
161 |
|
|
END DO |
162 |
|
|
END DO |
163 |
|
|
END DO |
164 |
|
|
DEALLOCATE(xc,h,delta_xc) |
165 |
|
|
|
166 |
|
|
END FUNCTION slopes3 |
167 |
|
|
! |
168 |
|
|
!------------------------------------------------------------------------------- |
169 |
|
|
|
170 |
|
|
|
171 |
|
|
!------------------------------------------------------------------------------- |
172 |
|
|
! |
173 |
|
|
PURE FUNCTION slopes4(ix, f, x) |
174 |
|
|
! |
175 |
|
|
!------------------------------------------------------------------------------- |
176 |
|
|
! Arguments: |
177 |
|
|
INTEGER, INTENT(IN) :: ix |
178 |
|
|
REAL, INTENT(IN) :: f(:, :, :, :) |
179 |
|
|
REAL, INTENT(IN) :: x(:) |
180 |
|
|
REAL :: slopes4(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4)) |
181 |
|
|
!------------------------------------------------------------------------------- |
182 |
|
|
! Local: |
183 |
|
|
INTEGER :: n, i, j, k, l, m, sta(4), sto(4) |
184 |
|
|
REAL, ALLOCATABLE :: xc(:) ! (n) cell centers |
185 |
|
|
REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) |
186 |
|
|
REAL :: fm, ff, fp, dx |
187 |
|
|
!------------------------------------------------------------------------------- |
188 |
|
|
n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) |
189 |
|
|
FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. |
190 |
|
|
FORALL(i=2:n-1) |
191 |
|
|
h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) |
192 |
|
|
END FORALL |
193 |
|
|
slopes4(:,:,:,:)=0. |
194 |
|
|
sta=[1,1,1,1]; sta(ix)=2 |
195 |
|
|
sto=SHAPE(f); sto(ix)=n-1 |
196 |
|
|
DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l) |
197 |
|
|
DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) |
198 |
|
|
DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) |
199 |
|
|
DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) |
200 |
|
|
ff=f(i,j,k,l) |
201 |
|
|
SELECT CASE(ix) |
202 |
|
|
CASE(1); fm=f(i-1,j,k,l); fp=f(i+1,j,k,l) |
203 |
|
|
CASE(2); fm=f(i,j-1,k,l); fp=f(i,j+1,k,l) |
204 |
|
|
CASE(3); fm=f(i,j,k-1,l); fp=f(i,j,k+1,l) |
205 |
|
|
CASE(4); fm=f(i,j,k,l-1); fp=f(i,j,k,l+1) |
206 |
|
|
END SELECT |
207 |
|
|
IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN |
208 |
|
|
slopes4(i,j,k,l)=0.; CYCLE !--- Local extremum |
209 |
|
|
!--- 2nd order slope ; (fm, ff, fp) strictly monotonous |
210 |
|
|
slopes4(i,j,k,l)=(fp-fm)/dx |
211 |
|
|
!--- Slope limitation |
212 |
|
|
slopes4(i,j,k,l) = SIGN(MIN(ABS(slopes4(i,j,k,l)), & |
213 |
|
|
ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes4(i,j,k,l) ) |
214 |
|
|
END IF |
215 |
|
|
END DO |
216 |
|
|
END DO |
217 |
|
|
END DO |
218 |
|
|
END DO |
219 |
|
|
DEALLOCATE(xc,h,delta_xc) |
220 |
|
|
|
221 |
|
|
END FUNCTION slopes4 |
222 |
|
|
! |
223 |
|
|
!------------------------------------------------------------------------------- |
224 |
|
|
|
225 |
|
|
|
226 |
|
|
!------------------------------------------------------------------------------- |
227 |
|
|
! |
228 |
|
|
PURE FUNCTION slopes5(ix, f, x) |
229 |
|
|
! |
230 |
|
|
!------------------------------------------------------------------------------- |
231 |
|
|
! Arguments: |
232 |
|
|
INTEGER, INTENT(IN) :: ix |
233 |
|
|
REAL, INTENT(IN) :: f(:, :, :, :, :) |
234 |
|
|
REAL, INTENT(IN) :: x(:) |
235 |
|
|
REAL :: slopes5(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4),SIZE(f,5)) |
236 |
|
|
!------------------------------------------------------------------------------- |
237 |
|
|
! Local: |
238 |
|
|
INTEGER :: n, i, j, k, l, m, sta(5), sto(5) |
239 |
|
|
REAL, ALLOCATABLE :: xc(:) ! (n) cell centers |
240 |
|
|
REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) |
241 |
|
|
REAL :: fm, ff, fp, dx |
242 |
|
|
!------------------------------------------------------------------------------- |
243 |
|
|
n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) |
244 |
|
|
FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. |
245 |
|
|
FORALL(i=2:n-1) |
246 |
|
|
h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) |
247 |
|
|
END FORALL |
248 |
|
|
slopes5(:,:,:,:,:)=0. |
249 |
|
|
sta=[1,1,1,1,1]; sta(ix)=2 |
250 |
|
|
sto=SHAPE(f); sto(ix)=n-1 |
251 |
|
|
DO m=sta(5),sto(5); IF(ix==5) dx=delta_xc(m) |
252 |
|
|
DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l) |
253 |
|
|
DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) |
254 |
|
|
DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) |
255 |
|
|
DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) |
256 |
|
|
ff=f(i,j,k,l,m) |
257 |
|
|
SELECT CASE(ix) |
258 |
|
|
CASE(1); fm=f(i-1,j,k,l,m); fp=f(i+1,j,k,l,m) |
259 |
|
|
CASE(2); fm=f(i,j-1,k,l,m); fp=f(i,j+1,k,l,m) |
260 |
|
|
CASE(3); fm=f(i,j,k-1,l,m); fp=f(i,j,k+1,l,m) |
261 |
|
|
CASE(4); fm=f(i,j,k,l-1,m); fp=f(i,j,k,l+1,m) |
262 |
|
|
CASE(5); fm=f(i,j,k,l,m-1); fp=f(i,j,k,l,m+1) |
263 |
|
|
END SELECT |
264 |
|
|
IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN |
265 |
|
|
slopes5(i,j,k,l,m)=0.; CYCLE !--- Local extremum |
266 |
|
|
!--- 2nd order slope ; (fm, ff, fp) strictly monotonous |
267 |
|
|
slopes5(i,j,k,l,m)=(fp-fm)/dx |
268 |
|
|
!--- Slope limitation |
269 |
|
|
slopes5(i,j,k,l,m) = SIGN(MIN(ABS(slopes5(i,j,k,l,m)), & |
270 |
|
|
ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes5(i,j,k,l,m) ) |
271 |
|
|
END IF |
272 |
|
|
END DO |
273 |
|
|
END DO |
274 |
|
|
END DO |
275 |
|
|
END DO |
276 |
|
|
END DO |
277 |
|
|
DEALLOCATE(xc,h,delta_xc) |
278 |
|
|
|
279 |
|
|
END FUNCTION slopes5 |
280 |
|
|
! |
281 |
|
|
!------------------------------------------------------------------------------- |
282 |
|
|
|
283 |
|
|
END MODULE slopes_m |