source: LMDZ6/branches/Amaury_dev/libf/misc/lmdz_slopes.f90 @ 5122

Last change on this file since 5122 was 5119, checked in by abarral, 4 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way
File size: 9.9 KB
Line 
1MODULE lmdz_slopes
2
3  ! Author: Lionel GUEZ
4  ! Extension / factorisation: David CUGNET
5
6  IMPLICIT NONE; PRIVATE
7  PUBLIC :: slopes
8
9  ! Those generic function computes second order slopes with Van
10  ! Leer slope-limiting, given cell averages. Reference: Dukowicz,
11  ! 1987, SIAM Journal on Scientific and Statistical Computing, 8,
12  ! 305.
13
14  ! The only difference between the specific functions is the rank
15  ! of the first argument and the equal rank of the result.
16
17  ! slope(ix,...) acts on ix th dimension.
18
19  ! REAL, INTENT(IN), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1
20  ! REAL, INTENT(IN):: x(:) ! (n + 1) cell edges
21  ! real slopes, same shape as f ! (n, ...)
22  INTERFACE slopes
23     MODULE procedure slopes1, slopes2, slopes3, slopes4, slopes5
24  END INTERFACE
25
26
27
28CONTAINS
29
30!-------------------------------------------------------------------------------
31
32PURE 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
65END FUNCTION slopes1
66
67!-------------------------------------------------------------------------------
68
69
70!-------------------------------------------------------------------------------
71
72PURE 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
114END FUNCTION slopes2
115
116!-------------------------------------------------------------------------------
117
118
119!-------------------------------------------------------------------------------
120
121PURE 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
166END FUNCTION slopes3
167
168!-------------------------------------------------------------------------------
169
170
171!-------------------------------------------------------------------------------
172
173PURE 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
221END FUNCTION slopes4
222
223!-------------------------------------------------------------------------------
224
225
226!-------------------------------------------------------------------------------
227
228PURE 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
279END FUNCTION slopes5
280
281!-------------------------------------------------------------------------------
282
283END MODULE lmdz_slopes
Note: See TracBrowser for help on using the repository browser.