module_tools.f90
Go to the documentation of this file.
1 !*************************************************************************
2 ! This file is part of HERCULES.
3 !
4 ! HERCULES is free software; you can redistribute it and/or modify
5 ! it under the terms of the GNU General Public License as published by
6 ! the Free Software Foundation; either version 3 of the License, or
7 ! (at your option) any later version.
8 !
9 ! HERCULES is distributed in the hope that it will be useful,
10 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ! GNU General Public License for more details.
13 !
14 ! You should have received a copy of the GNU General Public License
15 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
16 !
17 ! Copyright 2015, Ping He, MEAS, NCSU
18 !*************************************************************************
19 
20 
21 
23 
25  implicit none
26 
27 
28 contains
29 
30 
31 
32  ! this function allocates arrays, initiates fields etc
33  subroutine initiate_fields
34 
36 
37  implicit none
38 
39  ! allocate arrays, by default, all the variables are on z pencil
40  ! field variables
41  allocate( u(nia3,nja3,nka3) )
42  allocate( v(nia3,nja3,nka3) )
43  allocate( w(nia3,nja3,nka3) )
44  allocate( p(nia3,nja3,nka3) )
45  allocate( phi(nia3,nja3,nka3) )
46  allocate( t(nia3,nja3,nka3) )
47 
48  ! convective terms
49  allocate( ff(nia3,nja3,nka3) )
50  allocate( gg(nia3,nja3,nka3) )
51  allocate( hh(nia3,nja3,nka3) )
52  allocate( ss(nia3,nja3,nka3) )
53 
54  if (cds .ne. 3) then ! for spectral method, we don't need these vars
55  ! source terms
56  allocate( qu(nia3,nja3,nka3) )
57  allocate( qv(nia3,nja3,nka3) )
58  allocate( qw(nia3,nja3,nka3) )
59  ! convection terms at t=n-1
60  allocate( ff0(nia3,nja3,nka3) )
61  allocate( gg0(nia3,nja3,nka3) )
62  allocate( hh0(nia3,nja3,nka3) )
63  allocate( ss0(nia3,nja3,nka3) )
64  ! interpolated variables
65  allocate( uhx(nia3,nja3,nka3) )
66  allocate( uhy(nia3,nja3,nka3) )
67  allocate( uhz(nia3,nja3,nka3) )
68  allocate( vhx(nia3,nja3,nka3) )
69  allocate( vhy(nia3,nja3,nka3) )
70  allocate( vhz(nia3,nja3,nka3) )
71  allocate( whx(nia3,nja3,nka3) )
72  allocate( why(nia3,nja3,nka3) )
73  allocate( thx(nia3,nja3,nka3) )
74  allocate( thy(nia3,nja3,nka3) )
75  allocate( thz(nia3,nja3,nka3) )
76  endif
77 
78  ! mean and turbulence statistics
79  allocate( mean_u(nia3,nja3,nka3) )
80  allocate( mean_v(nia3,nja3,nka3) )
81  allocate( mean_w(nia3,nja3,nka3) )
82  allocate( mean_t(nia3,nja3,nka3) )
83  allocate( mean_uu(nia3,nja3,nka3) )
84  allocate( mean_vv(nia3,nja3,nka3) )
85  allocate( mean_ww(nia3,nja3,nka3) )
86  allocate( mean_uw(nia3,nja3,nka3) )
87  allocate( mean_vw(nia3,nja3,nka3) )
88  allocate( mean_tt(nia3,nja3,nka3) )
89  allocate( mean_tw(nia3,nja3,nka3) )
90 
91  ! coefficients for Thomas algorithm,
92  ! aa, bb, cc are also on z pencil
93  allocate( aa(zsize(1),zsize(2),zsize(3)+2*ghst) )
94  allocate( bb(zsize(1),zsize(2),zsize(3)+2*ghst) )
95  allocate( cc(zsize(1),zsize(2),zsize(3)+2*ghst) )
96  allocate( aa_u( zsize(3)+2*ghst) )
97  allocate( bb_u(zsize(1),zsize(2),zsize(3)+2*ghst) )
98  allocate( cc_u( zsize(3)+2*ghst) )
99  allocate( aa_w( zsize(3)+2*ghst) )
100  allocate( bb_w(zsize(1),zsize(2),zsize(3)+2*ghst) )
101  allocate( cc_w( zsize(3)+2*ghst) )
102  allocate( aa_p( zsize(3)+2*ghst) )
103  allocate( bb_p(zsize(1),zsize(2),zsize(3)+2*ghst) )
104  allocate( cc_p( zsize(3)+2*ghst) )
105  allocate( aa_t( zsize(3)+2*ghst) )
106  allocate( bb_t(zsize(1),zsize(2),zsize(3)+2*ghst) )
107  allocate( cc_t( zsize(3)+2*ghst) )
108 
109  ! temporary arrays
110  ! c_src, and c_var are on z pencil
111  allocate( c_src(zsize(1),zsize(2),zsize(3)) )
112  allocate( c_var(zsize(1),zsize(2),zsize(3)) )
113  ! c_tmp1, and c_tmp2 are on x and y pencil
114  allocate( c_tmp1(xsize(1),xsize(2),xsize(3)) )
115  allocate( c_tmp2(ysize(1),ysize(2),ysize(3)) )
116 
117  if (cds .eq. 3) then
118  !!! variables for spectral method
119  ! Fourier coefficients of the fundamental variables
120  allocate( c_u(zsize(1),zsize(2),zsize(3)) )
121  allocate( c_v(zsize(1),zsize(2),zsize(3)) )
122  allocate( c_w(zsize(1),zsize(2),zsize(3)) )
123  allocate( c_p(zsize(1),zsize(2),zsize(3)) )
124  allocate(c_phi(zsize(1),zsize(2),zsize(3)) )
125  allocate( c_t(zsize(1),zsize(2),zsize(3)) )
126 
127  allocate( c_qu(zsize(1),zsize(2),zsize(3)) )
128  allocate( c_qv(zsize(1),zsize(2),zsize(3)) )
129  allocate( c_qw(zsize(1),zsize(2),zsize(3)) )
130  allocate( c_ff(zsize(1),zsize(2),zsize(3)) )
131  allocate( c_gg(zsize(1),zsize(2),zsize(3)) )
132  allocate( c_hh(zsize(1),zsize(2),zsize(3)) )
133  allocate( c_ss(zsize(1),zsize(2),zsize(3)) )
134  endif
135 
136  ! initiate the fields with zeros
137  ! field variables
138  u(:,:,:)=0.d0
139  v(:,:,:)=0.d0
140  w(:,:,:)=0.d0
141  t(:,:,:)=0.d0
142  p(:,:,:)=0.d0
143  phi(:,:,:)=0.d0
144 
145  ! convective terms
146  ff(:,:,:)=0.d0
147  gg(:,:,:)=0.d0
148  hh(:,:,:)=0.d0
149  ss(:,:,:)=0.d0
150 
151  if (cds .ne. 3) then
152  ! source terms
153  qu(:,:,:)=0.d0
154  qv(:,:,:)=0.d0
155  qw(:,:,:)=0.d0
156  ! convection terms at t=n-1
157  ff0(:,:,:)=0.d0
158  gg0(:,:,:)=0.d0
159  hh0(:,:,:)=0.d0
160  ss0(:,:,:)=0.d0
161  ! these are interpolated variables for u,v,w,t
162  uhx(:,:,:)=0.d0
163  uhy(:,:,:)=0.d0
164  uhz(:,:,:)=0.d0
165  vhx(:,:,:)=0.d0
166  vhy(:,:,:)=0.d0
167  vhz(:,:,:)=0.d0
168  whx(:,:,:)=0.d0
169  why(:,:,:)=0.d0
170  thx(:,:,:)=0.d0
171  thy(:,:,:)=0.d0
172  thz(:,:,:)=0.d0
173  endif
174 
175  ! Thomas algorithm coefficients
176  aa(:,:,:)=0.d0
177  bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
178  cc(:,:,:)=0.d0
179  aa_u(:)=0.d0
180  bb_u(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
181  cc_u(:)=0.d0
182  aa_w(:)=0.d0
183  bb_w(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
184  cc_w(:)=0.d0
185  aa_p(:)=0.d0
186  bb_p(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
187  cc_p(:)=0.d0
188  aa_t(:)=0.d0
189  bb_t(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
190  cc_t(:)=0.d0
191 
192  ! mean and turbulence statistics
193  mean_u(:,:,:)=0.d0
194  mean_v(:,:,:)=0.d0
195  mean_w(:,:,:)=0.d0
196  mean_t(:,:,:)=0.d0
197  mean_uu(:,:,:)=0.d0
198  mean_vv(:,:,:)=0.d0
199  mean_ww(:,:,:)=0.d0
200  mean_uw(:,:,:)=0.d0
201  mean_vw(:,:,:)=0.d0
202  mean_tt(:,:,:)=0.d0
203  mean_tw(:,:,:)=0.d0
204 
205  ! these are temporary arrays
206  c_src(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
207  c_var(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
208  c_tmp1(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
209  c_tmp2(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
210 
211  if (cds .eq. 3) then
212  ! Fourier coefficients
213  c_u(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
214  c_v(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
215  c_w(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
216  c_t(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
217  c_p(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
218  c_phi(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
219 
220  c_qu(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
221  c_qv(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
222  c_qw(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
223  c_ff(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
224  c_gg(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
225  c_hh(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
226  c_ss(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
227  endif
228 
229  if (restart .eq. 1) then
230 
231  ! read initial fields from files
232  call read_init_fields(u,v,w,p,t,ff,gg,hh,ss)
233 
234  ! add some noise to the initial fields?
235  if (isnoise .eq. 1) then
236  call add_noise2fields(u,v,w,p,t)
237  endif
238 
239  elseif (restart .eq. 0) then
240 
241  ! use default initial values (all zeros)
242 
243  ! add some noise to the initial fields?
244  if (isnoise .eq. 1) then
245  call add_noise2fields(u,v,w,p,t)
246  endif
247 
248  else
249 
250  print *, 'restart error!'
251 
252  endif
253 
254  return
255  end subroutine
256 
257 
258 
259  ! this function reads initial fields from files
260  subroutine read_init_fields(u,v,w,p,t,ff,gg,hh,ss)
262  use decomp_2d_io
263 
264  implicit none
265 
266  real(mytype),dimension(:,:,:),intent(out) :: u, v, w, t, p
267  real(mytype),dimension(:,:,:),intent(out) :: ff,gg,hh,ss
268 
269  if (myid .eq. 0) then
270  print *, 'Reading Initial Fields.'
271  endif
272 
273  call decomp_2d_read_one(3,u(istr3:iend3,jstr3:jend3,kstr3:kend3), &
274  'results/init_u.dat')
275 
276  call decomp_2d_read_one(3,v(istr3:iend3,jstr3:jend3,kstr3:kend3), &
277  'results/init_v.dat')
278 
279  call decomp_2d_read_one(3,w(istr3:iend3,jstr3:jend3,kstr3:kend3), &
280  'results/init_w.dat')
281 
282  call decomp_2d_read_one(3,p(istr3:iend3,jstr3:jend3,kstr3:kend3), &
283  'results/init_p.dat')
284 
285  call decomp_2d_read_one(3,t(istr3:iend3,jstr3:jend3,kstr3:kend3), &
286  'results/init_t.dat')
287 
288  ! subtract the moving frame velocity u_mrf from u
291 
292  ! for RK3, we don't need to read f,g,h,s from the previous steps
293  if (dts .ne. 2) then
294 
295  call decomp_2d_read_one(3,ff(istr3:iend3,jstr3:jend3,kstr3:kend3), &
296  'results/init_ff.dat')
297 
298  call decomp_2d_read_one(3,gg(istr3:iend3,jstr3:jend3,kstr3:kend3), &
299  'results/init_gg.dat')
300 
301  call decomp_2d_read_one(3,hh(istr3:iend3,jstr3:jend3,kstr3:kend3), &
302  'results/init_hh.dat')
303 
304  call decomp_2d_read_one(3,ss(istr3:iend3,jstr3:jend3,kstr3:kend3), &
305  'results/init_ss.dat')
306 
307  endif
308 
309 
310  if (myid .eq. 0) then
311  print *, 'Reading Initial Fields Finished!'
312  endif
313 
314  return
315 
316  end subroutine
317 
318 
319 
320  ! this function adds some noise to the fields
321  subroutine add_noise2fields(u,v,w,p,t)
323  implicit none
324 
325  real(mytype),dimension(:,:,:),intent(inout) :: u, v, w, t, p
326 
327  integer :: i,j,k
328  real(mytype) :: tmp
329 
330  if (myid .eq. 0) then
331  print *,'Adding some noise to the initial fields'
332  endif
333 
334  call random_seed
335 
336  do k=kstr3,kend3,1
337  do j=jstr3,jend3,1
338  do i=istr3,iend3,1
339 
340  call random_number(tmp)
341  u(i,j,k)=u(i,j,k)+u(i,j,k)*tmp*noise_mag
342  call random_number(tmp)
343  v(i,j,k)=v(i,j,k)+v(i,j,k)*tmp*noise_mag
344  call random_number(tmp)
345  w(i,j,k)=w(i,j,k)+w(i,j,k)*tmp*noise_mag
346  call random_number(tmp)
347  p(i,j,k)=p(i,j,k)+p(i,j,k)*tmp*noise_mag
348  call random_number(tmp)
349  t(i,j,k)=t(i,j,k)+t(i,j,k)*tmp*noise_mag
350 
351  enddo
352  enddo
353  enddo
354 
355  return
356 
357  end subroutine
358 
359 
360 
361  ! this function calculates mean fields and turbulence statistics
362  subroutine calculate_mean_fields &
363  (u,v,w,t,uhx,vhy,mean_u,mean_v,mean_w,mean_t, mean_uu,mean_vv, &
364  mean_ww,mean_uw,mean_vw,mean_tt,mean_tw)
366  implicit none
367 
368  real(mytype),dimension(:,:,:),intent(in) :: u,v,w,t,uhx, vhy
369  real(mytype),dimension(:,:,:),intent(out):: &
370  mean_u,mean_v,mean_w,mean_t, mean_uu,mean_vv,mean_ww,mean_uw, mean_vw, &
371  mean_tt,mean_tw
372 
373  integer :: i,j,k
374 
375  do k=kstr3,kend3,1
376  do j=jstr3,jend3,1
377  do i=istr3,iend3,1
378 
379  ! note that here we need to add the speed of moving reference frame
380  mean_uu(i,j,k)=mean_uu(i,j,k)+(u(i,j,k)+u_mrf)**2.d0
381  mean_vv(i,j,k)=mean_vv(i,j,k)+v(i,j,k)**2.d0
382  mean_ww(i,j,k)=mean_ww(i,j,k)+w(i,j,k)**2.d0
383  mean_uw(i,j,k)=mean_uw(i,j,k)+0.5d0*( w(i,j,k)+w(i,j,k-1) ) * &
384  ( uhx(i,j,k)+u_mrf )
385  mean_vw(i,j,k)=mean_vw(i,j,k)+0.5d0*( w(i,j,k)+w(i,j,k-1) ) * vhy(i,j,k)
386 
387  mean_tt(i,j,k)=mean_tt(i,j,k)+t(i,j,k)**2.d0
388  mean_tw(i,j,k)=mean_tw(i,j,k)+t(i,j,k) * 0.5d0 * &
389  ( w(i,j,k)+w(i,j,k-1) )
390 
391  ! note that here we need to add the speed of moving reference frame
392  mean_u(i,j,k)=mean_u(i,j,k) + u(i,j,k)+u_mrf
393  mean_v(i,j,k)=mean_v(i,j,k) + v(i,j,k)
394  mean_w(i,j,k)=mean_w(i,j,k) + w(i,j,k)
395  mean_t(i,j,k)=mean_t(i,j,k) + t(i,j,k)
396 
397  enddo
398  enddo
399  enddo
400 
401 
402  return
403 
404  end subroutine
405 
406 
407 
408  ! this function calculates the interpolated fields for u,v,w
409  ! we don't need to call this function for spectral method
410  subroutine get_interp_fields_uvwhxyz(u,v,w,uhx,uhy,uhz,vhx,vhy,vhz,whx,why)
412  implicit none
413 
414  real(mytype),dimension(:,:,:),intent(in) :: u,v,w
415  real(mytype),dimension(:,:,:),intent(out):: uhx,uhy,uhz,vhx,vhy,vhz,&
416  whx,why
417 
418  integer :: i,j,k
419 
420  ! calculate interpolation variables
421  do k=kstr3,kend3,1
422  do j=jstr3,jend3,1
423  do i=istr3,iend3,1
424 
425  if (cds .eq. 1) then ! linear interpolation
426 
427  uhx(i,j,k)=0.5d0*( u(i,j,k)+u(i-1,j,k) )
428  uhy(i,j,k)=0.5d0*( u(i,j+1,k)+u(i,j,k) )
429  vhx(i,j,k)=0.5d0*( v(i+1,j,k)+v(i,j,k) )
430  vhy(i,j,k)=0.5d0*( v(i,j,k)+v(i,j-1,k) )
431  whx(i,j,k)=0.5d0*( w(i+1,j,k)+w(i,j,k) )
432  why(i,j,k)=0.5d0*( w(i,j+1,k)+w(i,j,k) )
433 
434  elseif (cds .eq. 2) then ! 4th order interpolation
435 
436  uhx(i,j,k)=( -u(i+1,j,k)+9.d0*u(i,j,k)+ &
437  9.d0*u(i-1,j,k)-u(i-2,j,k) ) / 16.d0
438 
439  uhy(i,j,k)=( -u(i,j+2,k)+9.d0*u(i,j+1,k)+ &
440  9.d0*u(i,j,k)-u(i,j-1,k) ) / 16.d0
441 
442  vhx(i,j,k)=( -v(i+2,j,k)+9.d0*v(i+1,j,k)+ &
443  9.d0*v(i,j,k)-v(i-1,j,k) ) / 16.d0
444 
445  vhy(i,j,k)=( -v(i,j+1,k)+9.d0*v(i,j,k)+ &
446  9.d0*v(i,j-1,k)-v(i,j-2,k) ) / 16.d0
447 
448  whx(i,j,k)=( -w(i+2,j,k)+9.d0*w(i+1,j,k)+ &
449  9.d0*w(i,j,k)-w(i-1,j,k) ) / 16.d0
450 
451  why(i,j,k)=( -w(i,j+2,k)+9.d0*w(i,j+1,k)+ &
452  9.d0*w(i,j,k)-w(i,j-1,k) ) / 16.d0
453 
454  endif
455 
456  ! this is always 2nd order interpolation
457  uhz(i,j,k)=u(i,j,k+1)*l_t(k)+u(i,j,k)*(1.d0-l_t(k))
458  ! this is always 2nd order interpolation
459  vhz(i,j,k)=v(i,j,k+1)*l_t(k)+v(i,j,k)*(1.d0-l_t(k))
460 
461  enddo
462  enddo
463  enddo
464 
465 
466  return
467 
468  end subroutine
469 
470 
471  ! this function calculates the interpolated fields for T
472  ! we don't need to call this function for spectral method
473  subroutine get_interp_fields_thxyz(t,thx,thy,thz)
475  implicit none
476 
477  real(mytype),dimension(:,:,:),intent(in) :: t
478  real(mytype),dimension(:,:,:),intent(out):: thx,thy,thz
479 
480  integer :: i,j,k
481 
482  ! calculate interpolation variables
483  do k=kstr3,kend3,1
484  do j=jstr3,jend3,1
485  do i=istr3,iend3,1
486 
487  if (cds .eq. 1) then
488 
489  thx(i,j,k)=0.5d0*( t(i+1,j,k)+t(i,j,k) )
490  thy(i,j,k)=0.5d0*( t(i,j+1,k)+t(i,j,k) )
491 
492  elseif (cds .eq. 2) then
493 
494  thx(i,j,k)=( -t(i+2,j,k)+9.d0*t(i+1,j,k)+ &
495  9.d0*t(i,j,k)-t(i-1,j,k) ) / 16.d0
496 
497  thy(i,j,k)=( -t(i,j+2,k)+9.d0*t(i,j+1,k)+ &
498  9.d0*t(i,j,k)-t(i,j-1,k) ) / 16.d0
499 
500  endif
501 
502  ! this is always 2nd order interpolation
503  thz(i,j,k)=t(i,j,k+1)*l_t(k)+t(i,j,k)*(1.d0-l_t(k))
504 
505  enddo
506  enddo
507  enddo
508 
509 
510  return
511 
512  end subroutine
513 
514 
515 
516  ! this function conducts 2d fft or ifft
517  ! we use 2d decomposition lib for array transposition between processors
518  subroutine mpi_2d_fft(c_inout3,direction)
520  ! C standard for using fftw
521  use, intrinsic :: iso_c_binding
522  ! this is a temporary array, it is on x/y pencil, it is used for 2d fft
523  use module_variables, only : c_tmp1,c_tmp2
524 
525  implicit none
526 
527  ! for using fftw
528  include 'include/fftw3.f03'
529 
530  ! flag: -1=fft; 1=ifft
531  integer, intent(in) :: direction
532  ! note that c_inout3 is on z pencil
533  complex(mytype),dimension(:,:,:),intent(inout) :: c_inout3
534 
535  integer :: i,j,k,l,m
536 
537 
538  ! we always remove the Nyquist wavenumber
539  if (direction .eq. 1) then
540 
541  do k=1,zsize(3),1
542  do j=1,zsize(2),1
543  do i=1,zsize(1),1
544 
545  ! deal with index
546  ! convert i to l
547  if ( (i+i_offset) .le. nx/2) then
548  l=i+i_offset-1
549  else
550  l=i+i_offset-nx-1
551  endif
552  ! convert j to m
553  if ( (j+j_offset) .le. ny/2) then
554  m=j+j_offset-1
555  else
556  m=j+j_offset-ny-1
557  endif
558 
559  ! remove Nyquist wavenumber
560  if (l .eq. -nx/2) then
561  c_inout3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
562  endif
563 
564  if (m .eq. -ny/2) then
565  c_inout3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
566  endif
567 
568  enddo
569  enddo
570  enddo
571 
572  endif
573 
574 
575  ! mpi 2d fft
576  ! first transpose from z to y
577  call transpose_z_to_y(c_inout3,c_tmp2)
578 
579  ! check it is a forward (direction=-1) or a backward (direction=-1) fft
580  ! note that we are on y pencil now!!
581  if (direction .eq. -1) then
582 
583  ! do 1d fft for c_tmp2, we are on y pencil
584  do k=1,ysize(3),1
585  do i=1,ysize(1),1 ! y pencil index
586 
587  call fftw_execute_dft(fft_plan2,c_tmp2(i,:,k), &
588  c_tmp2(i,:,k))
589 
590  enddo
591  enddo
592  ! ***NOTE*** since fftw is not normalized,
593  ! we need to normalize it after the forward fft
594  c_tmp2(:,:,:)=c_tmp2(:,:,:)/ny
595 
596  elseif (direction .eq. 1) then
597 
598  ! do 1d ifft for c_tmp2, we are on y pencil
599  do k=1,ysize(3),1
600  do i=1,ysize(1),1 ! y pencil index
601 
602  call fftw_execute_dft(ifft_plan2,c_tmp2(i,:,k), &
603  c_tmp2(i,:,k))
604 
605  enddo
606  enddo
607 
608  else
609 
610  print *, "direction is either -1 (forward) or 1 (backward)"
611  stop
612 
613  endif
614 
615  ! transpose from y to x
616  call transpose_y_to_x(c_tmp2,c_tmp1)
617 
618  ! check if it is a forward (direction=-1) or a backward (direction=-1) fft
619  ! NOTE we are on x pencil now!!!
620  if (direction .eq. -1) then
621 
622  ! do 1d fft for c_tmp1, we are on x pencil
623  do k=1,xsize(3),1
624  do j=1,xsize(2),1 ! x pencil index
625 
626  call fftw_execute_dft(fft_plan1,c_tmp1(:,j,k), &
627  c_tmp1(:,j,k))
628 
629  enddo
630  enddo
631  ! ***NOTE*** since fftw is not normalized,
632  ! we need to normalize it after the forward fft
633  c_tmp1(:,:,:)=c_tmp1(:,:,:)/nx
634 
635  elseif (direction .eq. 1) then
636 
637  ! do 1d ifft for c_tmp1, we are on x pencil
638  do k=1,xsize(3),1
639  do j=1,xsize(2),1 ! x pencil index
640 
641  call fftw_execute_dft(ifft_plan1,c_tmp1(:,j,k), &
642  c_tmp1(:,j,k))
643 
644  enddo
645  enddo
646 
647  else
648 
649  print *, "direction is either -1 (forward) or 1 (backward)"
650  stop
651 
652  endif
653 
654  ! transpose back from x pencil to z pencil
655  ! the output is c_inout3
656  call transpose_x_to_y(c_tmp1,c_tmp2)
657  call transpose_y_to_z(c_tmp2,c_inout3)
658 
659  ! we always remove the Nyquist wave number
660  if (direction .eq. -1) then
661 
662  do k=1,zsize(3),1
663  do j=1,zsize(2),1
664  do i=1,zsize(1),1
665 
666  ! deal with index
667  ! convert i to l
668  if ( (i+i_offset) .le. nx/2) then
669  l=i+i_offset-1
670  else
671  l=i+i_offset-nx-1
672  endif
673  ! convert j to m
674  if ( (j+j_offset) .le. ny/2) then
675  m=j+j_offset-1
676  else
677  m=j+j_offset-ny-1
678  endif
679 
680  ! remove Nyquist wavenumber
681  if (l .eq. -nx/2) then
682  c_inout3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
683  endif
684 
685  if (m .eq. -ny/2) then
686  c_inout3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
687  endif
688 
689  enddo
690  enddo
691  enddo
692 
693  endif
694 
695 
696  ! mpi 2d fft finished
697 
698  return
699 
700  end subroutine
701 
702 
703 
704  ! Rayleigh damping layer at the top of the boundary layer
705  subroutine rayleigh_damping(varin)
707  use mpi
708 
709  implicit none
710 
711  real(mytype),dimension(:,:,:),intent(inout) :: varin
712 
713  real(mytype) :: zd,sigz,mean_local,mean_all
714  integer :: i,j,k
715 
716 
717  ! do the damping for the topmost nzdamp points
718  do k=kend3-nzdamp,kend3,1
719 
720  ! first calculate the mean of the variable on the local processor
721  mean_local=0.d0
722  do j=jstr3,jend3,1
723  do i=istr3,iend3,1
724  mean_local=mean_local+varin(i,j,k)
725  enddo
726  enddo
727  mean_local=mean_local/zsize(1)/zsize(2)
728 
729  ! exchange mean_local between processors
730  ! sum up all the mean_local between all processors
731  call mpi_allreduce(mean_local,mean_all,1,real_type,mpi_sum, &
732  mpi_comm_world,ierr)
733  mean_all=mean_all/nprc ! calculate mean all
734 
735  ! do the damping
736  zd=(k-kend3+nzdamp)*1.d0/nzdamp
737  sigz=cdamp*0.5d0*( 1.d0-cos(pi*zd) )
738  do j=jstr3,jend3,1
739  do i=istr3,iend3,1
740  varin(i,j,k)=varin(i,j,k)-sigz*( varin(i,j,k)-mean_all )
741  enddo
742  enddo
743 
744  enddo
745 
746 
747  return
748 
749  end subroutine
750 
751 
752 
753 
754 
755 
756 end module
757 
758 
759 
real(mytype), dimension(:,:,:), allocatable, save mean_uw
real(mytype), dimension(:), allocatable, save aa_w
subroutine initiate_fields
integer, save, protected myid
real(mytype), dimension(:), allocatable, save aa_t
complex(mytype), dimension(:,:,:), allocatable, save c_qu
integer, save, protected nprc
integer, save, protected isnoise
integer, save, protected ny
subroutine add_noise2fields(u, v, w, p, t)
subroutine mpi_2d_fft(c_inout3, direction)
real(mytype), dimension(:,:,:), allocatable, save vhy
complex(mytype), dimension(:,:,:), allocatable, save bb_u
real(mytype), dimension(:,:,:), allocatable, save thx
real(mytype), dimension(:,:,:), allocatable, save ss
complex(mytype), dimension(:,:,:), allocatable, save c_u
complex(mytype), dimension(:,:,:), allocatable, save c_qv
real(mytype), dimension(:,:,:), allocatable, save uhz
complex(mytype), dimension(:,:,:), allocatable, save c_src
real(mytype), parameter pi
real(mytype), dimension(:), allocatable, save aa_p
integer, save, protected dts
complex(mytype), dimension(:,:,:), allocatable, save c_ss
real(mytype), dimension(:,:,:), allocatable, save mean_w
complex(mytype), dimension(:,:,:), allocatable, save c_hh
real(mytype), dimension(:,:,:), allocatable, save ff0
integer, save, protected jend3
real(mytype), dimension(:,:,:), allocatable, save mean_u
real(mytype), dimension(:), allocatable, save cc_p
complex(mytype), dimension(:,:,:), allocatable, save c_p
subroutine read_init_fields(u, v, w, p, t, ff, gg, hh, ss)
real(mytype), dimension(:), allocatable, save cc_w
real(mytype), dimension(:,:,:), allocatable, save ss0
real(mytype), dimension(:,:,:), allocatable, save mean_v
real(mytype), dimension(:,:,:), allocatable, save cc
real(mytype), dimension(:,:,:), allocatable, save p
real(mytype), dimension(:,:,:), allocatable, save thz
real(mytype), dimension(:,:,:), allocatable, save mean_vw
integer, save, protected nia3
type(c_ptr), save, protected fft_plan1
real(mytype), dimension(:,:,:), allocatable, save aa
complex(mytype), dimension(:,:,:), allocatable, save c_tmp2
complex(mytype), dimension(:,:,:), allocatable, save c_phi
complex(mytype), dimension(:,:,:), allocatable, save bb_t
real(mytype), dimension(:), allocatable, save cc_u
real(mytype), dimension(:,:,:), allocatable, save mean_tt
real(mytype), save, protected cdamp
complex(mytype), dimension(:,:,:), allocatable, save bb_w
complex(mytype), dimension(:,:,:), allocatable, save c_gg
real(mytype), dimension(:,:,:), allocatable, save gg
subroutine rayleigh_damping(varin)
real(mytype), dimension(:,:,:), allocatable, save vhx
integer, save, protected istr3
real(mytype), dimension(:,:,:), allocatable, save hh0
complex(mytype), dimension(:,:,:), allocatable, save c_tmp1
real(mytype), dimension(:,:,:), allocatable, save mean_ww
integer, save, protected jstr3
real(mytype), dimension(:,:,:), allocatable, save mean_t
real(mytype), dimension(:,:,:), allocatable, save w
real(mytype), dimension(:,:,:), allocatable, save u
type(c_ptr), save, protected ifft_plan1
integer, save, protected i_offset
integer, save, protected nja3
type(c_ptr), save, protected ifft_plan2
integer, save, protected j_offset
real(mytype), save, protected u_mrf
real(mytype), dimension(:,:,:), allocatable, save hh
complex(mytype), dimension(:,:,:), allocatable, save c_v
complex(mytype), dimension(:,:,:), allocatable, save c_ff
complex(mytype), dimension(:,:,:), allocatable, save c_qw
integer, parameter ghst
integer, save, protected cds
real(mytype), dimension(:,:,:), allocatable, save mean_vv
integer, save, protected restart
subroutine get_interp_fields_thxyz(t, thx, thy, thz)
real(mytype), dimension(:), allocatable, save aa_u
complex(mytype), dimension(:,:,:), allocatable, save c_w
real(mytype), dimension(:,:,:), allocatable, save t
integer, save, protected nzdamp
real(mytype), save, protected noise_mag
real(mytype), dimension(:,:,:), allocatable, save qw
real(mytype), dimension(:,:,:), allocatable, save whx
real(mytype), dimension(:,:,:), allocatable, save uhx
complex(mytype), dimension(:,:,:), allocatable, save c_t
real(mytype), dimension(:,:,:), allocatable, save why
integer, save, protected nx
real(mytype), dimension(:,:,:), allocatable, save mean_tw
real(mytype), dimension(:), allocatable, save cc_t
real(mytype), dimension(:,:,:), allocatable, save ff
real(mytype), dimension(:,:,:), allocatable, save mean_uu
real(mytype), dimension(:,:,:), allocatable, save v
type(c_ptr), save, protected fft_plan2
real(mytype), dimension(:,:,:), allocatable, save gg0
real(mytype), dimension(:,:,:), allocatable, save qu
integer, save, protected iend3
real(mytype), dimension(:), allocatable, save, protected l_t
real(mytype), dimension(:,:,:), allocatable, save phi
subroutine get_interp_fields_uvwhxyz(u, v, w, uhx, uhy, uhz, vhx, vhy, vhz, whx, why)
complex(mytype), dimension(:,:,:), allocatable, save bb
subroutine calculate_mean_fields(u, v, w, t, uhx, vhy, mean_u, mean_v, mean_w, mean_t, mean_uu, mean_vv, mean_ww, mean_uw, mean_vw, mean_tt, mean_tw)
real(mytype), dimension(:,:,:), allocatable, save thy
real(mytype), dimension(:,:,:), allocatable, save qv
complex(mytype), dimension(:,:,:), allocatable, save c_var
real(mytype), dimension(:,:,:), allocatable, save uhy
complex(mytype), dimension(:,:,:), allocatable, save bb_p
integer, save, protected kstr3
integer, save, protected nka3
real(mytype), dimension(:,:,:), allocatable, save vhz
integer, save, protected kend3