source: LMDZ5/branches/IPSLCM6.0.8/libf/misc/slopes_m.F90

Last change on this file was 2471, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2434:2457 into testing branch
File size: 6.3 KB
Line 
1module slopes_m
2
3  ! Author: Lionel GUEZ
4
5  implicit none
6
7  interface slopes
8     ! This 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     ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1
17     ! real, intent(in):: x(:) ! (n + 1) cell edges
18     ! real slopes, same shape as f ! (n, ...)
19
20     module procedure slopes1, slopes2, slopes3, slopes4
21  end interface
22
23  private
24  public slopes
25
26contains
27
28  pure function slopes1(f, x)
29
30    real, intent(in):: f(:)
31    real, intent(in):: x(:)
32    real slopes1(size(f))
33
34    ! Local:
35    integer n, i
36    real xc(size(f)) ! (n) cell centers
37    real h
38
39    !------------------------------------------------------
40
41    n = size(f)
42    forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
43    slopes1(1) = 0.
44    slopes1(n) = 0.
45
46    do i = 2, n - 1
47       if (f(i) >= max(f(i - 1), f(i + 1)) .or. f(i) &
48            <= min(f(i - 1), f(i + 1))) then
49          ! Local extremum
50          slopes1(i) = 0.
51       else
52          ! (f(i - 1), f(i), f(i + 1)) strictly monotonous
53
54          ! Second order slope:
55          slopes1(i) = (f(i + 1) - f(i - 1)) / (xc(i + 1) - xc(i - 1))
56
57          ! Slope limitation:
58          h = abs(x(i + 1) - xc(i))
59          slopes1(i) = sign(min(abs(slopes1(i)), abs(f(i + 1) - f(i)) / h, &
60               abs(f(i) - f(i - 1)) / h), slopes1(i))
61       end if
62    end do
63
64  end function slopes1
65
66  !*************************************************************
67
68  pure function slopes2(f, x)
69
70    real, intent(in):: f(:, :)
71    real, intent(in):: x(:)
72    real slopes2(size(f, 1), size(f, 2))
73
74    ! Local:
75    integer n, i, j
76    real xc(size(f, 1)) ! (n) cell centers
77    real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
78
79    !------------------------------------------------------
80
81    n = size(f, 1)
82    forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
83
84    forall (i = 2:n - 1)
85       h(i) = abs(x(i + 1) - xc(i))
86       delta_xc(i) = xc(i + 1) - xc(i - 1)
87    end forall
88
89    do j = 1, size(f, 2)
90       slopes2(1, j) = 0.
91
92       do i = 2, n - 1
93          if (f(i, j) >= max(f(i - 1, j), f(i + 1, j)) .or. &
94               f(i, j) <= min(f(i - 1, j), f(i + 1, j))) then
95             ! Local extremum
96             slopes2(i, j) = 0.
97          else
98             ! (f(i - 1, j), f(i, j), f(i + 1, j))
99             ! strictly monotonous
100
101             ! Second order slope:
102             slopes2(i, j) = (f(i + 1, j) - f(i - 1, j)) / delta_xc(i)
103
104             ! Slope limitation:
105             slopes2(i, j) = sign(min(abs(slopes2(i, j)), &
106                  abs(f(i + 1, j) - f(i, j)) / h(i), &
107                  abs(f(i, j) - f(i - 1, j)) / h(i)), slopes2(i, j))
108          end if
109       end do
110
111       slopes2(n, j) = 0.
112    end do
113
114  end function slopes2
115
116  !*************************************************************
117
118  pure function slopes3(f, x)
119
120    real, intent(in):: f(:, :, :)
121    real, intent(in):: x(:)
122    real slopes3(size(f, 1), size(f, 2), size(f, 3))
123
124    ! Local:
125    integer n, i, j, k
126    real xc(size(f, 1)) ! (n) cell centers
127    real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
128
129    !------------------------------------------------------
130
131    n = size(f, 1)
132    forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
133
134    forall (i = 2:n - 1)
135       h(i) = abs(x(i + 1) - xc(i))
136       delta_xc(i) = xc(i + 1) - xc(i - 1)
137    end forall
138
139    do k = 1, size(f, 3)
140       do j = 1, size(f, 2)
141          slopes3(1, j, k) = 0.
142
143          do i = 2, n - 1
144             if (f(i, j, k) >= max(f(i - 1, j, k), f(i + 1, j, k)) .or. &
145                  f(i, j, k) <= min(f(i - 1, j, k), f(i + 1, j, k))) then
146                ! Local extremum
147                slopes3(i, j, k) = 0.
148             else
149                ! (f(i - 1, j, k), f(i, j, k), f(i + 1, j, k))
150                ! strictly monotonous
151
152                ! Second order slope:
153                slopes3(i, j, k) = (f(i + 1, j, k) - f(i - 1, j, k)) &
154                     / delta_xc(i)
155
156                ! Slope limitation:
157                slopes3(i, j, k) = sign(min(abs(slopes3(i, j, k)), &
158                     abs(f(i + 1, j, k) - f(i, j, k)) / h(i), &
159                     abs(f(i, j, k) - f(i - 1, j, k)) / h(i)), slopes3(i, j, k))
160             end if
161          end do
162
163          slopes3(n, j, k) = 0.
164       end do
165    end do
166
167  end function slopes3
168
169  !*************************************************************
170
171  pure function slopes4(f, x)
172
173    real, intent(in):: f(:, :, :, :)
174    real, intent(in):: x(:)
175    real slopes4(size(f, 1), size(f, 2), size(f, 3), size(f, 4))
176
177    ! Local:
178    integer n, i, j, k, l
179    real xc(size(f, 1)) ! (n) cell centers
180    real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
181
182    !------------------------------------------------------
183
184    n = size(f, 1)
185    forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
186
187    forall (i = 2:n - 1)
188       h(i) = abs(x(i + 1) - xc(i))
189       delta_xc(i) = xc(i + 1) - xc(i - 1)
190    end forall
191
192    do l = 1, size(f, 4)
193       do k = 1, size(f, 3)
194          do j = 1, size(f, 2)
195             slopes4(1, j, k, l) = 0.
196
197             do i = 2, n - 1
198                if (f(i, j, k, l) >= max(f(i - 1, j, k, l), f(i + 1, j, k, l)) &
199                     .or. f(i, j, k, l) &
200                     <= min(f(i - 1, j, k, l), f(i + 1, j, k, l))) then
201                   ! Local extremum
202                   slopes4(i, j, k, l) = 0.
203                else
204                   ! (f(i - 1, j, k, l), f(i, j, k, l), f(i + 1, j, k, l))
205                   ! strictly monotonous
206
207                   ! Second order slope:
208                   slopes4(i, j, k, l) = (f(i + 1, j, k, l) &
209                        - f(i - 1, j, k, l)) / delta_xc(i)
210
211                   ! Slope limitation:
212                   slopes4(i, j, k, l) = sign(min(abs(slopes4(i, j, k, l)), &
213                        abs(f(i + 1, j, k, l) - f(i, j, k, l)) / h(i), &
214                        abs(f(i, j, k, l) - f(i - 1, j, k, l)) / h(i)), &
215                        slopes4(i, j, k, l))
216                end if
217             end do
218
219             slopes4(n, j, k, l) = 0.
220          end do
221       end do
222    end do
223
224  end function slopes4
225
226end module slopes_m
Note: See TracBrowser for help on using the repository browser.