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

Last change on this file since 5127 was 5119, checked in by abarral, 11 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.