source: LMDZ5/trunk/libf/misc/slopes_m.F90 @ 2554

Last change on this file since 2554 was 2440, checked in by lguez, 8 years ago

For read_climoz = 1 or 2, replaced first order conservative regridding
of ozone by second order conservative regridding, with Van Leer
slope-limiting. The replacement is done for both latitude and pressure
regridding. The replacement is beneficial if the resolution of the
input data is coarser than the resolution of LMDZ. If the resolution
of the input data is finer, then the replacement is neutral, it does
not change much.
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.