source: LMDZ5/branches/testing/libf/misc/slopes_m.F90 @ 5161

Last change on this file since 5161 was 2839, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2785:2838 into testing branch
File size: 9.9 KB
RevLine 
[2839]1MODULE slopes_m
[2440]2
3  ! Author: Lionel GUEZ
[2839]4  ! Extension / factorisation: David CUGNET
[2440]5
[2839]6  IMPLICIT NONE
[2440]7
[2839]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.
[2440]12
[2839]13  ! The only difference between the specific functions is the rank
14  ! of the first argument and the equal rank of the result.
[2440]15
[2839]16  ! slope(ix,...) acts on ix th dimension.
[2440]17
[2839]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
[2440]24
[2839]25  PRIVATE
26  PUBLIC :: slopes
[2440]27
[2839]28CONTAINS
[2440]29
[2839]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
[2440]64
[2839]65END FUNCTION slopes1
66!
67!-------------------------------------------------------------------------------
[2440]68
69
[2839]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)
[2440]113
[2839]114END FUNCTION slopes2
115!
116!-------------------------------------------------------------------------------
[2440]117
118
[2839]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)
[2440]165
[2839]166END FUNCTION slopes3
167!
168!-------------------------------------------------------------------------------
[2440]169
170
[2839]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)
[2440]220
[2839]221END FUNCTION slopes4
222!
223!-------------------------------------------------------------------------------
[2440]224
225
[2839]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)
[2440]278
[2839]279END FUNCTION slopes5
280!
281!-------------------------------------------------------------------------------
[2440]282
[2839]283END MODULE slopes_m
Note: See TracBrowser for help on using the repository browser.