SCALE-RM
scale_dft.F90
Go to the documentation of this file.
1 #include "scalelib.h"
2 ! special 2D-DFT routines for low wavenumber spectral transform
3 module scale_dft
4  use mpi
5  use scale_precision
6  use scale_prc, only: &
7  prc_local_comm_world, &
8  prc_myrank
9  use scale_prc_cartesc, only: &
10  prc_num_x, &
11  prc_num_y, &
12  prc_2drank
13  use scale_comm_cartesc, only: &
14  comm_datatype
15  use scale_const, only: &
16  pi => const_pi
17  implicit none
18 
19  public :: dft_setup
20  public :: dft_g2g
21  public :: dft_g2g_divfree
22 
23  integer, private :: IMAX, JMAX ! number of total grids for spectral transform
24  integer, private :: IGS, JGS ! global index
25  integer, private :: LMM, MMM ! maximum truncation wavenumber
26 
27  ! table for trigonometric function
28  real(RP), private, allocatable :: table_x(:,:), table_y(:,:)
29  real(RP), private, allocatable :: table_l(:), table_m(:)
30  ! work array
31  real(RP), private, allocatable :: work(:,:,:)
32 
33 contains
34 
35  subroutine dft_setup(KA,KS,KE,IA,IS,IE,JA,JS,JE,LM,MM)
36  implicit none
37  integer,intent(in) :: ka, ia, ja
38  integer,intent(in) :: ks, ke, is, ie, js, je
39  integer,intent(in) :: lm, mm
40  real(rp) :: x, y
41  integer :: i, j, k, l, m
42 
43  imax = (ie-is+1)*prc_num_x
44  jmax = (je-js+1)*prc_num_y
45 
46  igs = (ie-is+1)*prc_2drank(prc_myrank,1)+1
47  jgs = (je-js+1)*prc_2drank(prc_myrank,2)+1
48 
49  lmm = lm
50  mmm = mm
51 
52  allocate( table_x(is:ie,0:2*lm), table_y(js:je,0:2*mm) )
53  allocate( table_l(0:2*lm), table_m(0:2*mm) )
54  allocate( work(ka,0:2*lm,ja) )
55 
56  do i = is, ie
57  x = 2*pi/imax*(i-is+igs)
58  table_x(i,0) = 1
59  enddo
60 
61  do l = 1, lm
62  do i = is, ie
63  x = 2*pi/imax*(i-is+igs)
64  table_x(i,2*l-1) = cos(l*x)
65  table_x(i,2*l) = -sin(l*x)
66  enddo
67  enddo
68 
69  do j = js, je
70  y = 2*pi/jmax*(j-js+jgs)
71  table_y(j,0) = 1
72  enddo
73 
74  do m = 1, mm
75  do j = js, je
76  y = 2*pi/jmax*(j-js+jgs)
77  table_y(j,2*m-1) = cos(m*y)
78  table_y(j,2*m) = -sin(m*y)
79  enddo
80  enddo
81 
82  table_l(0) = 1
83  do l = 1, lm
84  table_l(2*l-1) = cos(pi*l/imax)
85  table_l(2*l) = sin(pi*l/imax)
86  enddo
87 
88  table_m(0) = 1
89  do m = 1, mm
90  table_m(2*m-1) = cos(pi*m/jmax)
91  table_m(2*m) = sin(pi*m/jmax)
92  enddo
93 
94  end subroutine dft_setup
95 
96  subroutine dft_g2s(KA,KS,KE,IA,IS,IE,JA,JS,JE,LM,MM,f,s)
97  ! Grid to spectral transformation
98  integer,intent(in) :: KA, IA, JA
99  integer,intent(in) :: KS, KE, IS, IE, JS, JE
100  integer,intent(in) :: LM, MM
101  real(RP), intent(in) :: f(KA, IA, JA) ! x: full level, y: full level
102  real(RP), intent(out) :: s(KA,0:2*LM,0:2*MM)
103  real(RP) :: work_s(KA,0:2*LM,0:2*MM)
104  real(RP) :: c, tb
105  integer :: i, j, k, l, m
106  integer :: ierr
107 
108  c = 1.0_rp/imax
109  do j = js, je
110  do l = 0, 2*lm
111  do k = ks, ke
112  work(k,l,j) = 0
113  enddo
114  do i = is, ie
115  tb = table_x(i,l)*c
116  do k = ks, ke
117  work(k,l,j) = work(k,l,j) + f(k,i,j)*tb
118  enddo
119  enddo
120  enddo
121  enddo
122 
123  do m = 0, 2*mm
124  do l = 0, 2*lm
125  do k = ks, ke
126  work_s(k,l,m) = 0
127  enddo
128  enddo
129  enddo
130 
131  c = 1.0_rp/jmax
132  do m = 0, 2*mm
133  do j = js, je
134  tb = table_y(j,m)*c
135  do l = 0, 2*lm
136  do k = ks, ke
137  work_s(k,l,m) = work_s(k,l,m) + work(k,l,j)*tb
138  enddo
139  enddo
140  enddo
141  enddo
142 
143  call mpi_allreduce(work_s, s, ka*(2*lm+1)*(2*mm+1), comm_datatype, mpi_sum, prc_local_comm_world, ierr)
144 
145  end subroutine dft_g2s
146 
147  subroutine dft_s2g(KA,KS,KE,IA,IS,IE,JA,JS,JE,LM,MM,s,f)
148  ! Grid to spectral transformation
149  integer,intent(in) :: KA, IA, JA
150  integer,intent(in) :: KS, KE, IS, IE, JS, JE
151  integer,intent(in) :: LM, MM
152  real(RP), intent(in) :: s(KA,0:2*LM,0:2*MM)
153  real(RP), intent(out) :: f(KA, IA, JA) ! x: full level, y: full level
154  real(RP) :: c, tb
155  integer :: i, j, k, l, m
156 
157  do j = js, je
158  do i = is, ie
159  do k = ks, ke
160  f(k,i,j) = 0
161  enddo
162  enddo
163  enddo
164 
165  do j = js, je
166  do l = 0, 2*lm
167  do k = ks, ke
168  work(k,l,j) = 0
169  enddo
170  enddo
171  enddo
172 
173  do m = 0, 2*mm
174  if( m == 0 ) then
175  c = 1
176  else
177  c = 2
178  endif
179  do j = js, je
180  tb = table_y(j,m)*c
181  do l = 0, 2*lm
182  do k = ks, ke
183  work(k,l,j) = work(k,l,j) + s(k,l,m)*tb
184  enddo
185  enddo
186  enddo
187  enddo
188 
189  do j = js, je
190  do i = is, ie
191  do k = ks, ke
192  f(k,i,j) = 0
193  enddo
194  enddo
195  do l = 0, 2*lm
196  if( l == 0 ) then
197  c = 1
198  else
199  c = 2
200  endif
201  do i = is, ie
202  tb = table_x(i,l)*c
203  do k = ks, ke
204  f(k,i,j) = f(k,i,j) + work(k,l,j)*tb
205  enddo
206  enddo
207  enddo
208  enddo
209 
210  end subroutine dft_s2g
211 
212  ! Grid to Grid transformation
213  subroutine dft_g2g(KA,KS,KE,IA,IS,IE,JA,JS,JE,LM,MM,f)
214  integer,intent(in) :: ka, ia, ja
215  integer,intent(in) :: ks, ke, is, ie, js, je
216  integer,intent(in) :: lm, mm
217  real(rp), intent(inout) :: f(ka, ia, ja) ! x,y: full or half level
218  real(rp) :: s(ka,0:2*lm,0:2*mm)
219 
220  call dft_g2s(ka,ks,ke,ia,is,ie,ja,js,je,lm,mm,f,s)
221  call dft_s2g(ka,ks,ke,ia,is,ie,ja,js,je,lm,mm,s,f)
222 
223  end subroutine dft_g2g
224 
225  ! Grid to Grid transformation of vector field with horizontal divergence filter
226  subroutine dft_g2g_divfree(KA,KS,KE,IA,IS,IE,JA,JS,JE,LM,MM,u,v)
227  integer,intent(in) :: ka, ia, ja
228  integer,intent(in) :: ks, ke, is, ie, js, je
229  integer,intent(in) :: lm, mm
230  real(rp), intent(inout) :: u(ka, ia, ja) ! x: half level, y: full level
231  real(rp), intent(inout) :: v(ka, ia, ja) ! x: full level, y: half level
232  integer :: k, l, m
233  real(rp) :: s1(ka,0:2*lm,0:2*mm)
234  real(rp) :: s2(ka,0:2*lm,0:2*mm)
235  real(rp) :: s3(ka,0:2*lm,0:2*mm)
236  real(rp) :: a, b, fac
237 
238  call dft_g2s(ka,ks,ke,ia,is,ie,ja,js,je,lm,mm,u,s1)
239 
240  ! phase shift
241  do m = 0, 2*mm
242  do l = 1, lm
243  do k = ks, ke
244  a = s1(k,2*l-1,m)
245  b = s1(k,2*l,m)
246  s1(k,2*l-1,m) = a*table_l(2*l-1) + b*table_l(2*l)
247  s1(k,2*l,m) = b*table_l(2*l-1) - a*table_l(2*l)
248  enddo
249  enddo
250  enddo
251 
252  call dft_g2s(ka,ks,ke,ia,is,ie,ja,js,je,lm,mm,v,s2)
253 
254  ! phase shift
255  do m = 1, mm
256  do l = 0, 2*lm
257  do k = ks, ke
258  a = s2(k,l,2*m-1)
259  b = s2(k,l,2*m)
260  s2(k,l,2*m-1) = a*table_m(2*m-1) + b*table_m(2*m)
261  s2(k,l,2*m) = b*table_m(2*m-1) - a*table_m(2*m)
262  enddo
263  enddo
264  enddo
265 
266  ! rotation
267  do m = 0, 2*mm
268  do l = 0, 2*lm
269  do k = ks, ke
270  s3(k,l,m) = 0
271  enddo
272  enddo
273  enddo
274 
275  ! ∂v/∂x
276  do m = 0, 2*mm
277  do l = 1, lm
278  do k = ks, ke
279  s3(k,2*l-1,m) = -l*s2(k,2*l,m)
280  s3(k,2*l,m) = l*s2(k,2*l-1,m)
281  enddo
282  enddo
283  enddo
284 
285  ! + (- ∂u/∂y)
286  do m = 1, mm
287  do l = 0, 2*lm
288  do k = ks, ke
289  s3(k,l,2*m-1) = s3(k,l,2*m-1) + m*s1(k,l,2*m)
290  s3(k,l,2*m) = s3(k,l,2*m) - m*s1(k,l,2*m-1)
291  enddo
292  enddo
293  enddo
294 
295  ! minus inverse laplacian ( stream function on model plane )
296  do k = ks, ke
297  s3(k,0,0) = 0
298  enddo
299 
300  do l = 1, lm
301  fac = 1.0_rp/ (l*l)
302  do k = ks, ke
303  s3(k,2*l-1,0) = s3(k,2*l-1,0)*fac
304  s3(k,2*l,0) = s3(k,2*l,0)*fac
305  enddo
306  enddo
307 
308  do m = 1, mm
309  fac = 1.0_rp/ (m*m)
310  do k = ks, ke
311  s3(k,0,2*m-1) = s3(k,0,2*m-1)*fac
312  s3(k,0,2*m) = s3(k,0,2*m)*fac
313  enddo
314  enddo
315 
316  do m = 1, mm
317  do l = 1, lm
318  fac = 1.0_rp/ (l*l+m*m)
319  do k = ks, ke
320  s3(k,2*l-1,2*m-1) = s3(k,2*l-1,2*m-1)*fac
321  s3(k,2*l-1,2*m) = s3(k,2*l-1,2*m)*fac
322  s3(k,2*l,2*m-1) = s3(k,2*l,2*m-1)*fac
323  s3(k,2*l,2*m) = s3(k,2*l,2*m)*fac
324  enddo
325  enddo
326  enddo
327 
328  ! divergence free 2D vector
329 
330  ! ∂ψ/∂y
331  do l = 1, 2*lm
332  do k = ks, ke
333  s1(k,l,0) = 0
334  enddo
335  enddo
336  do m = 1, mm
337  do l = 0, 2*lm
338  do k = ks, ke
339  s1(k,l,2*m-1) = -m*s3(k,l,2*m)
340  s1(k,l,2*m) = m*s3(k,l,2*m-1)
341  enddo
342  enddo
343  enddo
344 
345  ! -∂ψ/∂x
346  do m = 1, 2*mm
347  do k = ks, ke
348  s2(k,0,m) = 0
349  enddo
350  enddo
351  do m = 0, 2*mm
352  do l = 1, lm
353  do k = ks, ke
354  s2(k,2*l-1,m) = l*s3(k,2*l,m)
355  s2(k,2*l,m) = -l*s3(k,2*l-1,m)
356  enddo
357  enddo
358  enddo
359 
360  ! phase shift
361  do m = 0, 2*mm
362  do l = 1, lm
363  do k = ks, ke
364  a = s1(k,2*l-1,m)
365  b = s1(k,2*l,m)
366  s1(k,2*l-1,m) = a*table_l(2*l-1) - b*table_l(2*l)
367  s1(k,2*l,m) = b*table_l(2*l-1) + a*table_l(2*l)
368  enddo
369  enddo
370  enddo
371 
372  call dft_s2g(ka,ks,ke,ia,is,ie,ja,js,je,lm,mm,s1,u)
373 
374  ! phase shift
375  do m = 1, mm
376  do l = 0, 2*lm
377  do k = ks, ke
378  a = s2(k,l,2*m-1)
379  b = s2(k,l,2*m)
380  s2(k,l,2*m-1) = a*table_m(2*m-1) - b*table_m(2*m)
381  s2(k,l,2*m) = b*table_m(2*m-1) + a*table_m(2*m)
382  enddo
383  enddo
384  enddo
385 
386  call dft_s2g(ka,ks,ke,ia,is,ie,ja,js,je,lm,mm,s2,v)
387 
388  end subroutine dft_g2g_divfree
389 
390 end module scale_dft
scale_comm_cartesc::comm_datatype
integer, public comm_datatype
datatype of variable
Definition: scale_comm_cartesC.F90:98
scale_dft::dft_g2s
subroutine dft_g2s(KA, KS, KE, IA, IS, IE, JA, JS, JE, LM, MM, f, s)
Definition: scale_dft.F90:97
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_prc::prc_local_comm_world
integer, public prc_local_comm_world
local communicator
Definition: scale_prc.F90:88
scale_dft::dft_g2g_divfree
subroutine, public dft_g2g_divfree(KA, KS, KE, IA, IS, IE, JA, JS, JE, LM, MM, u, v)
Definition: scale_dft.F90:227
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:90
scale_dft
Definition: scale_dft.F90:3
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_prc_cartesc::prc_2drank
integer, dimension(:,:), allocatable, public prc_2drank
node index in 2D topology
Definition: scale_prc_cartesC.F90:44
scale_prc_cartesc
module process / cartesC
Definition: scale_prc_cartesC.F90:11
scale_dft::dft_setup
subroutine, public dft_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, LM, MM)
Definition: scale_dft.F90:36
scale_dft::dft_g2g
subroutine, public dft_g2g(KA, KS, KE, IA, IS, IE, JA, JS, JE, LM, MM, f)
Definition: scale_dft.F90:214
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_prc_cartesc::prc_num_y
integer, public prc_num_y
y length of 2D processor topology
Definition: scale_prc_cartesC.F90:42
scale_prc_cartesc::prc_num_x
integer, public prc_num_x
x length of 2D processor topology
Definition: scale_prc_cartesC.F90:41
scale_dft::dft_s2g
subroutine dft_s2g(KA, KS, KE, IA, IS, IE, JA, JS, JE, LM, MM, s, f)
Definition: scale_dft.F90:148
Generated by   doxygen 1.8.17