SCALE-RM
scale_filter.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  use scale_prof
19  !-----------------------------------------------------------------------------
20  implicit none
21  private
22  !-----------------------------------------------------------------------------
23  !
24  !++ Public procedure
25  !
26  public :: filter_hyperdiff
27 
28  interface filter_hyperdiff
29  module procedure filter_hyperdiff_2d
30  module procedure filter_hyperdiff_3d
31  end interface filter_hyperdiff
32 
33  !-----------------------------------------------------------------------------
34  !
35  !++ Public parameters & variables
36  !
37  !-----------------------------------------------------------------------------
38  !
39  !++ Private procedure
40  !
41  !-----------------------------------------------------------------------------
42  !
43  !++ Private parameters & variables
44  !
45  !-----------------------------------------------------------------------------
46 contains
47  !-----------------------------------------------------------------------------
49  subroutine filter_hyperdiff_2d( &
50  IA, IS, IE, JA, JS, JE, &
51  data, order, nite, &
52  limiter_sign )
53  use scale_comm_cartesc, only: &
54  comm_vars8, &
55  comm_wait
56  use scale_prc, only: &
57  prc_abort
58  integer, intent(in) :: IA, IS, IE
59  integer, intent(in) :: JA, JS, JE
60 
61  real(RP), intent(inout) :: data(IA,JA)
62 
63  integer, intent(in) :: order
64  integer, intent(in) :: nite
65 
66  real(RP), intent(in), optional :: limiter_sign(IA,JA)
67 
68  real(RP), pointer :: p1(:,:)
69  real(RP), pointer :: p2(:,:)
70  real(RP), target :: work1(IA,JA)
71  real(RP), target :: work2(IA,JA)
72 
73  logical :: limiter
74 
75  integer :: i, j
76  integer :: ite, n
77 
78  call prof_rapstart('FILTER', 3)
79 
80  if ( mod(order,2) .ne. 0 ) then
81  log_error("FILTER_hyperdiff_2D", *) "order must be even"
82  call prc_abort
83  end if
84 
85 
86  limiter = present( limiter_sign )
87 
88  ! reduce grid-scale variation
89  do ite = 1, nite
90 
91  call comm_vars8( data(:,:), 1 )
92  call comm_wait ( data(:,:), 1, .false. )
93 
94  !$omp parallel do
95  do j = 1, ja
96  do i = 1, ia
97  work2(i,j) = data(i,j)
98  end do
99  end do
100 
101  p1 => work2
102  p2 => work1
103  do n = 1, order
104  !$omp parallel do
105  do j = max(js,2), min(je,ja-1)
106  do i = max(is,2), min(ie,ia-1)
107  p2(i,j) = ( - p1(i+1,j) + p1(i,j)*2.0_rp - p1(i-1,j) &
108  - p1(i,j+1) + p1(i,j)*2.0_rp - p1(i,j-1) ) / 8.0_rp
109  end do
110  end do
111  if ( js == 1 ) then
112  do i = max(is,2), min(ie,ia-1)
113  p2(i,js) = ( - p1(i+1,js) + p1(i,js)*2.0_rp - p1(i-1,js) &
114  - p1(i,js+1) + p1(i,js) ) / 6.0_rp
115  end do
116  else
117  do i = max(is,2), min(ie,ia-1)
118  p2(i,js-1) = p2(i,js)
119  end do
120  end if
121  if ( je == ja ) then
122  do i = max(is,2), min(ie,ia-1)
123  p2(i,je) = ( - p1(i+1,je) + p1(i,je)*2.0_rp - p1(i-1,je) &
124  + p1(i,je) - p1(i,je-1) ) / 6.0_rp
125  end do
126  else
127  do i = max(is,2), min(ie,ia-1)
128  p2(i,je+1) = p2(i,je)
129  end do
130  end if
131  if ( is == 1 ) then
132  do j = max(js,2), min(je,ja-1)
133  p2(is,j) = ( - p1(is+1,j) + p1(is,j) &
134  - p1(is,j+1) + p1(is,j)*2.0_rp - p1(is,j-1) ) / 6.0_rp
135  end do
136  if ( js == 1 ) then
137  p2(is,js) = ( - p1(is+1,js) + p1(is,js) &
138  - p1(is,js+1) + p1(is,js) ) / 4.0_rp
139  end if
140  if ( je == ja ) then
141  p2(is,je) = ( - p1(is+1,je) + p1(is,je) &
142  + p1(is,je) - p1(is,je-1) ) / 4.0_rp
143  end if
144  else
145  do j = max(js-1,1), min(je+1,ja)
146  p2(is-1,j) = p2(is,j)
147  end do
148  end if
149  if ( ie == ia ) then
150  do j = max(js,2), min(je,ja-1)
151  p2(ie,j) = ( + p1(ie,j) - p1(ie-1,j) &
152  - p1(ie,j+1) + p1(ie,j)*2.0_rp - p1(ie,j-1) ) / 6.0_rp
153  end do
154  if ( js == 1 ) then
155  p2(ie,js) = ( + p1(ie,js) - p1(ie-1,js) &
156  - p1(ie,js+1) + p1(ie,js) ) / 4.0_rp
157  end if
158  if ( je == ja ) then
159  p2(ie,je) = ( + p1(ie,je) - p1(ie-1,je) &
160  + p1(ie,je) - p1(ie,je-1) ) / 4.0_rp
161  end if
162  else
163  do j = max(js-1,1), min(je+1,ja)
164  p2(ie+1,j) = p2(ie,j)
165  end do
166  end if
167 
168  call comm_vars8( p2(:,:), 1 )
169  call comm_wait ( p2(:,:), 1, .false. )
170 
171  if ( mod(n,2) == 0 ) then
172  p1 => work2
173  p2 => work1
174  else
175  p1 => work1
176  p2 => work2
177  end if
178  end do
179 
180  !$omp parallel do
181  do j = js, je
182  do i = is, ie
183  data(i,j) = data(i,j) - p1(i,j)
184  end do
185  end do
186 
187  if ( limiter ) then
188  !$omp parallel do
189  do j = js, je
190  do i = is, ie
191  data(i,j) = sign( max( data(i,j) * limiter_sign(i,j), 0.0_rp ), limiter_sign(i,j) )
192  end do
193  end do
194  end if
195  end do
196 
197  call prof_rapend('FILTER', 3)
198 
199  return
200  end subroutine filter_hyperdiff_2d
201 
202  !-----------------------------------------------------------------------------
204  subroutine filter_hyperdiff_3d( &
205  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
206  data, order, nite, &
207  limiter_sign )
208  use scale_comm_cartesc, only: &
209  comm_vars8, &
210  comm_wait
211  use scale_prc, only: &
212  prc_abort
213  integer, intent(in) :: KA, KS, KE
214  integer, intent(in) :: IA, IS, IE
215  integer, intent(in) :: JA, JS, JE
216 
217  real(RP), intent(inout) :: data(KA,IA,JA)
218 
219  integer, intent(in) :: order
220  integer, intent(in) :: nite
221 
222  real(RP), intent(in), optional :: limiter_sign(KA,IA,JA)
223 
224  real(RP) :: work_data(IA,JA)
225  real(RP) :: work_sign(IA,JA)
226 
227  logical :: flag
228 
229  integer :: k, i, j
230 
231  if ( present(limiter_sign) ) then
232  flag = .true.
233  else
234  flag = .false.
235  end if
236 
237  do k = ks, ke
238 
239  !$omp parallel do
240  do j = js, je
241  do i = is, ie
242  work_data(i,j) = data(k,i,j)
243  end do
244  end do
245  if ( flag ) then
246  !$omp parallel do
247  do j = js, je
248  do i = is, ie
249  work_sign(i,j) = limiter_sign(k,i,j)
250  end do
251  end do
252 
253  call filter_hyperdiff_2d( ia, is, ie, ja, js, je, &
254  work_data(:,:), order, nite, &
255  limiter_sign = work_sign )
256 
257  else
258 
259  call filter_hyperdiff_2d( ia, is, ie, ja, js, je, &
260  work_data(:,:), order, nite )
261  end if
262 
263 
264  !$omp parallel do
265  do j = 1, ja
266  do i = 1, ia
267  data(k,i,j) = work_data(i,j)
268  end do
269  end do
270 
271  end do
272 
273 
274  return
275  end subroutine filter_hyperdiff_3d
276 
277 end module scale_filter
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:159
scale_filter::filter_hyperdiff_3d
subroutine filter_hyperdiff_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, data, order, nite, limiter_sign)
Hyper diffusion filter 3D.
Definition: scale_filter.F90:208
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_io
module STDIO
Definition: scale_io.F90:10
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_filter
module FILTER
Definition: scale_filter.F90:11
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217