module_spectral.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 
29 contains
30 
31 
32 
33  ! this function calculates spatial derivatives and the
34  ! right hand side terms for momentum equations
35  subroutine get_momentum_cnvdiff_spec &
36  (u,v,w,p,t,c_u,c_v,c_w,c_p,c_ff,c_gg,c_hh,c_qu,c_qv,c_qw)
37 
38  ! these are temporary arrays, their dimensions are: nx/p_row,ny/p_col,nz
39  ! they are on z pencil
40  use module_variables, only : c_src,c_var
41  use module_tools
42 
43  implicit none
44 
45  real(mytype),dimension(:,:,:),intent(in) :: u,v,w,p,t
46  complex(mytype),dimension(:,:,:),intent(in) :: c_u,c_v,c_w,c_p
47  complex(mytype),dimension(:,:,:),intent(out) :: c_ff,c_gg,c_hh
48  complex(mytype),dimension(:,:,:),intent(out) :: c_qu,c_qv,c_qw
49 
50  integer :: i,j,k
51  real(mytype) :: tmp1,tmp2,a1,a2
52  real(mytype) :: coe_sk
53 
54  ! check if the skew-symmetric form is used
55  if (issk .eq. 1) then
56  coe_sk = 0.5d0
57  else
58  coe_sk = 1.d0
59  endif
60 
61  ! first calculate dpx_drive
62  call update_dpx_drive
63 
64  !!!!!!!!!! U equation !!!!!!!!!!!
65  !!!!!!!!!! now start building the right hand side !!!!!!!!!!
66  ! qu=u^n+dt*(rk_a*ff^n+rk_b*ff^n-1) + &
67  ! rk_c*0.5*dt*nu*(dudx2+dudy^2+dudz^2)^n-rk_c*dt*dpdx^n
68 
69  ! add u
70  c_qu(:,:,:)=c_u(:,:,:)
71 
72  ! add dudx2
73  call dev_fft_linear_cmplx(c_u,c_var,2,1)
74  c_qu(:,:,:)=c_qu(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
75  ! add dudy2
76  call dev_fft_linear_cmplx(c_u,c_var,2,2)
77  c_qu(:,:,:)=c_qu(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
78  ! add dudz2, we store dudz2 to tmp1 in physical space
79  do k=kstr3,kend3,1
80  do j=jstr3,jend3,1
81  do i=istr3,iend3,1
82 
83  tmp1= ( (u(i,j,k+1)-u(i,j,k))/dz_t(k) - &
84  (u(i,j,k)-u(i,j,k-1))/dz_b(k) ) / dz(k)
85  c_var(i-ghst,j-ghst,k-ghst)=cmplx(tmp1, 0.d0, kind=mytype)
86 
87  enddo
88  enddo
89  enddo
90  call mpi_2d_fft(c_var,-1)
91  c_qu(:,:,:)=c_qu(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
92 
93  ! add dpdx
94  call dev_fft_linear_cmplx(c_p,c_var,1,1)
95  c_qu(:,:,:)=c_qu(:,:,:)-rk_c*dt*c_var(:,:,:)
96 
97  ! add ff^n-1
98  c_qu(:,:,:)=c_qu(:,:,:)+rk_b*dt*c_ff(:,:,:)
99 
100  !!!!!!!!!!!!! calculate ff !!!!!!!!!!
101  ! ff=-0.5*(duudx+duvdy+duwdz)^n - &
102  ! 0.5*(u*dudx+v*dudy+w*dudz)^n+dpx_drive-u_mrf*dudx^n
103  c_ff(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
104  ! add u_mrf*dudx
105  call dev_fft_linear_cmplx(c_u,c_var,1,1)
106  c_ff(:,:,:)=c_ff(:,:,:)-u_mrf*c_var(:,:,:)
107  ! add duudx
108  call dev_fft_nonlinear_cmplx(u,u,c_var,1,0)
109  c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*c_var(:,:,:)
110  ! add duvdy
111  call dev_fft_nonlinear_cmplx(u,v,c_var,2,0)
112  c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*c_var(:,:,:)
113  ! add duwdz and dpx_drive, we store duwdz to tmp1 in physical space
114  do k=kstr3,kend3,1
115  do j=jstr3,jend3,1
116  do i=istr3,iend3,1
117 
118  a1=u(i,j,k+1)*l_t(k)+u(i,j,k)*(1.d0-l_t(k))
119  a2=u(i,j,k-1)*l_b(k)+u(i,j,k)*(1.d0-l_b(k))
120  tmp1 = ( a1 * w(i,j,k) - &
121  a2 * w(i,j,k-1) ) / dz(k)
122 
123  if (dp_opt .eq. 3) then
124  tmp2 = fc*( v(i,j,k) - vg )
125  else
126  tmp2 = dpx_drive
127  endif
128 
129  ! note that here we use 1/coe_sk*tmp2
130  ! since we will do a coe_sk*c_var later
131  c_var(i-ghst,j-ghst,k-ghst)= cmplx(1.d0/coe_sk*tmp2 - tmp1, &
132  0.d0, kind=mytype )
133 
134  enddo
135  enddo
136  enddo
137  call mpi_2d_fft(c_var,-1)
138  c_ff(:,:,:)=c_ff(:,:,:)+coe_sk*c_var(:,:,:)
139 
140  ! check if the skew-symmetric form is used,
141  ! if yes we need to add the other terms
142  if (issk .eq. 1) then
143 
144  ! add u*dudx
145  call dev_fft_linear_cmplx(c_u,c_src,1,1) ! c_dudx
146  call mpi_2d_fft(c_src,1)
147  do k=kstr3,kend3,1
148  do j=jstr3,jend3,1
149  do i=istr3,iend3,1
150 
151  ! multiply u and dudx in physical space
152  c_var(i-ghst,j-ghst,k-ghst)= &
153  cmplx( u(i,j,k) * real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
154  0.d0, kind=mytype )
155 
156  enddo
157  enddo
158  enddo
159  call mpi_2d_fft(c_var,-1)
160  c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*c_var(:,:,:)
161 
162  ! add v*dudy
163  call dev_fft_linear_cmplx(c_u,c_src,1,2) ! c_dudy
164  call mpi_2d_fft(c_src,1)
165  do k=kstr3,kend3,1
166  do j=jstr3,jend3,1
167  do i=istr3,iend3,1
168 
169  ! multiply u and dudx in physical space
170  c_var(i-ghst,j-ghst,k-ghst)= cmplx( v(i,j,k) * &
171  real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
172  0.d0, kind=mytype )
173 
174  enddo
175  enddo
176  enddo
177  call mpi_2d_fft(c_var,-1)
178  c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*c_var(:,:,:)
179 
180  ! add w*dudz
181  do k=kstr3,kend3,1
182  do j=jstr3,jend3,1
183  do i=istr3,iend3,1
184 
185  a1=u(i,j,k+1)*l_t(k)+u(i,j,k)*(1.d0-l_t(k))
186  a2=u(i,j,k-1)*l_b(k)+u(i,j,k)*(1.d0-l_b(k))
187  tmp1 = ( a1 - a2 ) / dz(k)
188  tmp2 = 0.5d0*( w(i,j,k) + w(i,j,k-1) )
189  c_var(i-ghst,j-ghst,k-ghst)= cmplx( tmp1*tmp2,0.d0, kind=mytype )
190 
191  enddo
192  enddo
193  enddo
194  call mpi_2d_fft(c_var,-1)
195  c_ff(:,:,:)=c_ff(:,:,:)-coe_sk*c_var(:,:,:)
196 
197  endif
198  !!!!!!!!!!!! finish ff !!!!!!!!!
199 
200  ! add ff to qu
201  c_qu(:,:,:)=c_qu(:,:,:)+rk_a*dt*c_ff(:,:,:)
202  !!!!!!!!!!!!!!! U equation finish !!!!!!
203 
204 
205 
206  !!!!!!!!!! V equation !!!!!!!!!!!
207  !!!!!!!!!! now start building the right hand side !!!!!!!!!!
208  ! qv=v^n+dt*(rk_a*gg^n+rk_b*gg^n-1) + &
209  ! rk_c*0.5*dt*nu*(dvdx2+dvdy^2+dvdz^2)^n-rk_c*dt*dpdy^n
210 
211  ! add v
212  c_qv(:,:,:)=c_v(:,:,:)
213 
214  ! add dvdx2
215  call dev_fft_linear_cmplx(c_v,c_var,2,1)
216  c_qv(:,:,:)=c_qv(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
217  ! add dvdy2
218  call dev_fft_linear_cmplx(c_v,c_var,2,2)
219  c_qv(:,:,:)=c_qv(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
220  ! add dvdz2, we store dvdz2 to tmp1 in physical space
221  do k=kstr3,kend3,1
222  do j=jstr3,jend3,1
223  do i=istr3,iend3,1
224 
225  tmp1 = ( (v(i,j,k+1)-v(i,j,k))/dz_t(k) - &
226  (v(i,j,k)-v(i,j,k-1))/dz_b(k) ) / dz(k)
227  c_var(i-ghst,j-ghst,k-ghst)=cmplx(tmp1, 0.d0, kind=mytype)
228 
229  enddo
230  enddo
231  enddo
232  call mpi_2d_fft(c_var,-1)
233  c_qv(:,:,:)=c_qv(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
234 
235  ! add dpdy
236  call dev_fft_linear_cmplx(c_p,c_var,1,2)
237  c_qv(:,:,:)=c_qv(:,:,:)-rk_c*dt*c_var(:,:,:)
238 
239  ! add gg0
240  c_qv(:,:,:)=c_qv(:,:,:)+rk_b*dt*c_gg(:,:,:)
241 
242  !!!!!!!!!!!!! calculate gg !!!!!!!!!!
243  ! gg=-0.5*(duvdx+dvvdy+dvwdz)^n - &
244  ! 0.5*(u*dvdx+v*dvdy+w*dvdz)^n-u_mrf*dvdx^n
245  c_gg(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
246  ! add u_mrf*dvdx
247  call dev_fft_linear_cmplx(c_v,c_var,1,1)
248  c_gg(:,:,:)=c_gg(:,:,:)-u_mrf*c_var(:,:,:)
249  ! add duvdx
250  call dev_fft_nonlinear_cmplx(u,v,c_var,1,0)
251  c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*c_var(:,:,:)
252  ! add dvvdy
253  call dev_fft_nonlinear_cmplx(v,v,c_var,2,0)
254  c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*c_var(:,:,:)
255  ! add dvwdz, we store dvwdz to tmp1 in physical space
256  do k=kstr3,kend3,1
257  do j=jstr3,jend3,1
258  do i=istr3,iend3,1
259 
260  a1=v(i,j,k+1)*l_t(k)+v(i,j,k)*(1.d0-l_t(k))
261  a2=v(i,j,k-1)*l_b(k)+v(i,j,k)*(1.d0-l_b(k))
262  tmp1 = ( a1 * w(i,j,k) - &
263  a2 * w(i,j,k-1) ) / dz(k)
264 
265  if (dp_opt .eq. 3) then
266  tmp2 = fc*( ug- u(i,j,k) - u_mrf )
267  else
268  tmp2 = 0.d0
269  endif
270 
271  c_var(i-ghst,j-ghst,k-ghst)= &
272  cmplx( 1.d0/coe_sk*tmp2-tmp1,0.d0, kind=mytype )
273 
274  enddo
275  enddo
276  enddo
277  call mpi_2d_fft(c_var,-1)
278  c_gg(:,:,:)=c_gg(:,:,:)+coe_sk*c_var(:,:,:)
279 
280 
281  ! check if the skew-symmetric form is used,
282  ! if yes we need to add the other terms
283  if (issk .eq. 1) then
284 
285  ! add u*dvdx
286  call dev_fft_linear_cmplx(c_v,c_src,1,1) ! c_dvdx
287  call mpi_2d_fft(c_src,1)
288  do k=kstr3,kend3,1
289  do j=jstr3,jend3,1
290  do i=istr3,iend3,1
291 
292  ! multiply u and dvdx in physical space
293  c_var(i-ghst,j-ghst,k-ghst)= cmplx( u(i,j,k) * &
294  real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
295  0.d0, kind=mytype )
296 
297  enddo
298  enddo
299  enddo
300  call mpi_2d_fft(c_var,-1)
301  c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*c_var(:,:,:)
302 
303  ! add v*dvdy
304  call dev_fft_linear_cmplx(c_v,c_src,1,2) ! c_dvdy
305  call mpi_2d_fft(c_src,1)
306  do k=kstr3,kend3,1
307  do j=jstr3,jend3,1
308  do i=istr3,iend3,1
309 
310  ! multiply v and dvdy in physical space
311  c_var(i-ghst,j-ghst,k-ghst)= cmplx( v(i,j,k) * &
312  real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
313  0.d0, kind=mytype )
314 
315  enddo
316  enddo
317  enddo
318  call mpi_2d_fft(c_var,-1)
319  c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*c_var(:,:,:)
320 
321  ! add w*dvdz
322  do k=kstr3,kend3,1
323  do j=jstr3,jend3,1
324  do i=istr3,iend3,1
325 
326  a1=v(i,j,k+1)*l_t(k)+v(i,j,k)*(1.d0-l_t(k))
327  a2=v(i,j,k-1)*l_b(k)+v(i,j,k)*(1.d0-l_b(k))
328  tmp1 = ( a1 - a2 ) / dz(k)
329  tmp2 = 0.5d0*( w(i,j,k) + w(i,j,k-1) )
330  c_var(i-ghst,j-ghst,k-ghst)= cmplx( tmp1*tmp2,0.d0, kind=mytype )
331 
332  enddo
333  enddo
334  enddo
335  call mpi_2d_fft(c_var,-1)
336  c_gg(:,:,:)=c_gg(:,:,:)-coe_sk*c_var(:,:,:)
337 
338  endif
339  !!!!!!!!!!!! finish gg !!!!!!!!!
340 
341  ! add gg to qv
342  c_qv(:,:,:)=c_qv(:,:,:)+rk_a*dt*c_gg(:,:,:)
343  !!!!!!!!!!!!!!! V equation finish !!!!!!
344 
345 
346 
347  !!!!!!!!!! W equation !!!!!!!!!!! ------ note that w is staggered
348  !!!!!!!!!! now start building the right hand side !!!!!!!!!!
349  ! qw=w^n+dt*(rk_a*hh^n+rk_b*hh^n-1) + &
350  ! rk_c*0.5*dt*nu*(dwdx2+dwdy^2+dwdz^2)^n-rk_c*dt*dpdz^n
351 
352  ! add w
353  c_qw(:,:,:)=c_w(:,:,:)
354 
355  ! add dwdx2
356  call dev_fft_linear_cmplx(c_w,c_var,2,1)
357  c_qw(:,:,:)=c_qw(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
358  ! add dwdy2
359  call dev_fft_linear_cmplx(c_w,c_var,2,2)
360  c_qw(:,:,:)=c_qw(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
361  ! add dwdz2, we store dwdz2 to tmp1 in physical space
362  do k=kstr3,kend3,1 !!!!! we might need a different range here
363  do j=jstr3,jend3,1
364  do i=istr3,iend3,1
365 
366  tmp1 = ( (w(i,j,k+1)-w(i,j,k))/dz(k+1) - &
367  (w(i,j,k)-w(i,j,k-1))/dz(k) ) / dz_t(k)
368 
369  c_var(i-ghst,j-ghst,k-ghst)=cmplx(tmp1, 0.d0, kind=mytype)
370 
371  enddo
372  enddo
373  enddo
374  call mpi_2d_fft(c_var,-1)
375  c_qw(:,:,:)=c_qw(:,:,:)+rk_c*0.5d0*dt*nu*c_var(:,:,:)
376 
377  ! add dpdz
378  do k=kstr3,kend3,1 !!!!! we might need a different range here
379  do j=jstr3,jend3,1
380  do i=istr3,iend3,1
381 
382  tmp1= ( p(i,j,k+1)-p(i,j,k) ) / dz_t(k)
383 
384  c_var(i-ghst,j-ghst,k-ghst)=cmplx(tmp1, 0.d0, kind=mytype)
385 
386  enddo
387  enddo
388  enddo
389  call mpi_2d_fft(c_var,-1)
390  c_qw(:,:,:)=c_qw(:,:,:)-rk_c*dt*c_var(:,:,:)
391 
392  ! add hh0
393  c_qw(:,:,:)=c_qw(:,:,:)+rk_b*dt*c_hh(:,:,:)
394 
395  !!!!!!!!!!!!! calculate hh !!!!!!!!!!
396  ! hh=-0.5*(duwdx+dvwdy+dwwdz)^n - &
397  ! 0.5*(u*dwdx+v*dwdy+w*dwdz)^n-u_mrf*dwdx^n + g*(T-T0)/T0
398  c_hh(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
399  ! add u_mrf*dwdx
400  call dev_fft_linear_cmplx(c_w,c_var,1,1)
401  c_hh(:,:,:)=c_hh(:,:,:)-u_mrf*c_var(:,:,:)
402  ! add duwdx
403  call dev_fft_nonlinear_cmplx(u,w,c_var,1,1)
404  c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*c_var(:,:,:)
405  ! add dvwdy
406  call dev_fft_nonlinear_cmplx(v,w,c_var,2,1)
407  c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*c_var(:,:,:)
408  ! add dwwdz and the buoyancy term, we store dwwdz to tmp1 in physical space
409  do k=kstr3,kend3,1 ! range
410  do j=jstr3,jend3,1
411  do i=istr3,iend3,1
412 
413  !dwwdz
414  a1 = ( 0.5d0*(w(i,j,k)+w(i,j,k+1)) )**2.d0
415  a2 = ( 0.5d0*(w(i,j,k)+w(i,j,k-1)) )**2.d0
416  tmp1= (a1-a2) / dz_t(k)
417 
418  if (isscalar .eq. 1) then
419  ! buoyancy
420  tmp2=got*( (1.d0-l_t(k))*t(i,j,k)+l_t(k)*t(i,j,k+1) -t_ref )
421  ! note that here we use 1/coe_sk*buoyancy
422  ! since we will do a coe_sk*c_var later
423  c_var(i-ghst,j-ghst,k-ghst) = &
424  cmplx( 1.d0/coe_sk*tmp2 - tmp1,0.d0, kind=mytype )
425  else
426  c_var(i-ghst,j-ghst,k-ghst)= cmplx( -tmp1,0.d0, kind=mytype )
427  endif
428 
429  enddo
430  enddo
431  enddo
432  call mpi_2d_fft(c_var,-1)
433  c_hh(:,:,:)=c_hh(:,:,:)+coe_sk*c_var(:,:,:)
434 
435  ! check if the skew-symmetric form is used,
436  ! if yes we need to add the other terms
437  if (issk .eq. 1) then
438 
439  ! add u*dwdx
440  call dev_fft_linear_cmplx(c_w,c_src,1,1) ! c_dwdx
441  call mpi_2d_fft(c_src,1)
442  do k=kstr3,kend3,1
443  do j=jstr3,jend3,1
444  do i=istr3,iend3,1
445 
446  ! multiply u and dwdx in physical space
447  tmp1= l_t(k)*u(i,j,k+1)+(1.d0-l_t(k))*u(i,j,k)
448  c_var(i-ghst,j-ghst,k-ghst)= cmplx( tmp1 * &
449  real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
450  0.d0, kind=mytype )
451 
452  enddo
453  enddo
454  enddo
455  call mpi_2d_fft(c_var,-1)
456  c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*c_var(:,:,:)
457 
458  ! add v*dwdy
459  call dev_fft_linear_cmplx(c_w,c_src,1,2) ! c_dwdy
460  call mpi_2d_fft(c_src,1)
461  do k=kstr3,kend3,1
462  do j=jstr3,jend3,1
463  do i=istr3,iend3,1
464 
465  ! multiply v and dwdy in physical space
466  tmp1= l_t(k)*v(i,j,k+1)+(1.d0-l_t(k))*v(i,j,k)
467  c_var(i-ghst,j-ghst,k-ghst)= cmplx( tmp1 * &
468  real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype) , &
469  0.d0, kind=mytype )
470 
471  enddo
472  enddo
473  enddo
474  call mpi_2d_fft(c_var,-1)
475  c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*c_var(:,:,:)
476 
477  ! add w*dwdz
478  do k=kstr3,kend3,1
479  do j=jstr3,jend3,1
480  do i=istr3,iend3,1
481 
482  a1 = 0.5d0*(w(i,j,k)+w(i,j,k+1))
483  a2 = 0.5d0*(w(i,j,k)+w(i,j,k-1))
484  tmp1= (a1-a2) / dz_t(k)
485  c_var(i-ghst,j-ghst,k-ghst) = &
486  cmplx( w(i,j,k)*tmp1,0.d0, kind=mytype )
487 
488  enddo
489  enddo
490  enddo
491  call mpi_2d_fft(c_var,-1)
492  c_hh(:,:,:)=c_hh(:,:,:)-coe_sk*c_var(:,:,:)
493 
494 
495  endif
496  !!!!!!!!!!!! finish hh !!!!!!!!!
497 
498  ! add hh to qw
499  c_qw(:,:,:)=c_qw(:,:,:)+rk_a*dt*c_hh(:,:,:)
500  !!!!!!!!!!!!!!! W equation finish !!!!!!
501 
502 
503  return
504  end subroutine
505 
506 
507  ! this function calculates spatial derivatives and the
508  ! right hand side terms for the T equation
509  subroutine get_temperature_cnvdiff_spec &
510  (u,v,w,t,c_t,c_ss,c_qt)
511 
512  ! these are temporary arrays, their dimensions are: nx/p_row,ny/p_col,nz
513  ! they are on z pencil
514  use module_variables, only : c_var
515  use module_tools
516 
517  implicit none
518 
519  real(mytype),dimension(:,:,:),intent(in) :: u,v,w,t
520  complex(mytype),dimension(:,:,:),intent(in) :: c_t
521  complex(mytype),dimension(:,:,:),intent(out) :: c_ss
522  complex(mytype),dimension(:,:,:),intent(out) :: c_qt
523 
524  integer :: i,j,k
525  real(mytype) :: tmp1,a1,a2
526 
527 
528  !!!!!!!!!! T equation !!!!!!!!!!!
529  !!!!!!!!!! now start building the right hand side !!!!!!!!!!
530  ! qt=T^n+dt*(rk_a*ss^n+rk_b*ss^n-1) + &
531  ! rk_c*0.5*dt*alpha*(dTdx2+dTdy^2+dTdz^2)^n
532 
533  ! add T
534  c_qt(:,:,:)=c_t(:,:,:)
535 
536  ! add dTdx2
537  call dev_fft_linear_cmplx(c_t,c_var,2,1)
538  c_qt(:,:,:)=c_qt(:,:,:)+rk_c*0.5d0*dt*alpha*c_var(:,:,:)
539  ! add dTdy2
540  call dev_fft_linear_cmplx(c_t,c_var,2,2)
541  c_qt(:,:,:)=c_qt(:,:,:)+rk_c*0.5d0*dt*alpha*c_var(:,:,:)
542  ! add dTdz2, we store dTdz2 to tmp1 in physical space
543  do k=kstr3,kend3,1
544  do j=jstr3,jend3,1
545  do i=istr3,iend3,1
546 
547  tmp1= ( (t(i,j,k+1)-t(i,j,k))/dz_t(k) - &
548  (t(i,j,k)-t(i,j,k-1))/dz_b(k) ) / dz(k)
549  c_var(i-ghst,j-ghst,k-ghst)=cmplx(tmp1, 0.d0, kind=mytype)
550 
551  enddo
552  enddo
553  enddo
554  call mpi_2d_fft(c_var,-1)
555  c_qt(:,:,:)=c_qt(:,:,:)+rk_c*0.5d0*dt*alpha*c_var(:,:,:)
556 
557  ! add ss^n-1
558  c_qt(:,:,:)=c_qt(:,:,:)+rk_b*dt*c_ss(:,:,:)
559 
560  !!!!!!!!!!!!! calculate ss !!!!!!!!!!
561  ! ss= - (dTudx+dTvdy+dTwdz)^n - u_mrf*dTdx
562  c_ss(:,:,:)=cmplx( 0.d0,0.d0,kind=mytype )
563  ! add u_mrf*dTdx
564  call dev_fft_linear_cmplx(c_t,c_var,1,1)
565  c_ss(:,:,:)=c_ss(:,:,:)-u_mrf*c_var(:,:,:)
566  ! add dTudx
567  call dev_fft_nonlinear_cmplx(u,t,c_var,1,0)
568  c_ss(:,:,:)=c_ss(:,:,:)-c_var(:,:,:)
569  ! add dTvdy
570  call dev_fft_nonlinear_cmplx(t,v,c_var,2,0)
571  c_ss(:,:,:)=c_ss(:,:,:)-c_var(:,:,:)
572  ! add dTwdz, we store dTwdz to tmp1 in physical space
573  do k=kstr3,kend3,1
574  do j=jstr3,jend3,1
575  do i=istr3,iend3,1
576 
577  a1=t(i,j,k+1)*l_t(k)+t(i,j,k)*(1.d0-l_t(k))
578  a2=t(i,j,k-1)*l_b(k)+t(i,j,k)*(1.d0-l_b(k))
579  tmp1 = ( a1 * w(i,j,k) - &
580  a2 * w(i,j,k-1) ) / dz(k)
581  c_var(i-ghst,j-ghst,k-ghst)= cmplx(tmp1,0.d0, kind=mytype )
582 
583  enddo
584  enddo
585  enddo
586  call mpi_2d_fft(c_var,-1)
587  c_ss(:,:,:)=c_ss(:,:,:)-c_var(:,:,:)
588 
589  ! add ss to qt
590  c_qt(:,:,:)=c_qt(:,:,:)+rk_a*dt*c_ss(:,:,:)
591  !!!!!!!!!!!!!!! U equation finish !!!!!!
592 
593 
594  return
595  end subroutine
596 
597 
598 
599  ! Poisson equation fft solver for spectral method,
600  ! we call it Poisson equation solver,
601  ! although we also use it to solve Helmholtz equation for U,V,W,T
602  subroutine poisson_solver_fft_spec &
603  (aa,bb,cc,c_src_in,c_var_out,zstagger,delnull,tbnd,ubnd)
605  ! C standard, for using fftw
606  use, intrinsic :: iso_c_binding
607  use module_tools
608  ! these are temporary arrays, their dimensions are: nx/p_row,ny/p_col,nz
609  ! they are on z pencil
610  use module_variables, only : c_src
611 
612  implicit none
613 
614  ! for using fftw
615  include 'include/fftw3.f03'
616 
617  ! these are Thomas algorithm coefficients
618  real(mytype), dimension(:,:,:),intent(inout) :: aa,cc
619  complex(mytype),dimension(:,:,:),intent(inout) :: bb
620  ! 3d real array, input from source term
621  complex(mytype), dimension(:,:,:), intent(in) :: c_src_in
622  ! these are flags
623  ! zstagger: if is staggered grid in z direction.
624  ! set it to 1 for w equation only
625  ! delnull: whether to delete null space, set it to 1 for p equation only
626  ! tbnd: whether to treat the Dirichlet boundary condition,
627  ! set it to 1 for temperature equation only
628  ! ubnd: special treatment for u boundary condition (for u_mrf),
629  ! this is because in a moving coordinate, the bottom wall u is not zero!!
630  ! set it to 1 for u equation only
631  integer, intent(in) :: zstagger,delnull,tbnd,ubnd
632  ! 3d real array, output for the Poisson fft solution
633  complex(mytype),dimension(:,:,:), intent(out) :: c_var_out
634 
635  ! temporary variable used in Thomas algorithm
636  complex(mytype) :: tmp
637  ! temporary variable used in T and u eqn boundary condition
638  real(mytype) :: aa_t_zstr,cc_t_zend
639  real(mytype) :: aa_u_zstr,cc_u_zend
640  integer :: i,j,k,kmax3
641 
642  ! initiate variables
643  c_var_out(:,:,:)=0.d0
644  c_src(:,:,:)=c_src_in(:,:,:)
645 
646  ! check some flags
647  ! check if it is w equation, the z range will be different
648  if (zstagger .eq. 1) then
649  kmax3=kend3-1
650  elseif (zstagger .eq. 0) then
651  kmax3=kend3
652  else
653  print *, 'Invalid zstagger flag!'
654  stop
655  endif
656  ! check if we need to remove the null space for pressure Poisson equation
657  ! we set the mean component of aa,cc,c_src to 0,
658  ! and therefore the mean pressure at the 1st grid level is 0
659  if (delnull .eq. 1) then
660  ! remove null space
661  if (myid .eq. 0) then
662  aa(1,1,kstr3)=0.d0
663  cc(1,1,kstr3)=0.d0
664  c_src(1,1,1)=cmplx(0.d0,0.d0 , kind=mytype )
665  endif
666  elseif (delnull .eq. 0) then
667 
668  else
669  print *, 'delnull flag error!'
670  stop
671  endif
672  ! check if it is t equation.
673  if (tbnd .eq. 1) then
674  ! add the mising terms to qt for the Dirichlet condition
675  if (myid .eq. 0) then
676  ! applied boundary condition for temperature
677  aa_t_zstr=-0.5d0*rk_c*dt*alpha/dz(kstr3)/dz_b(kstr3)
678  c_src(1,1,1)=c_src(1,1,1)-2.d0*aa_t_zstr*tbot
679  cc_t_zend=-0.5d0*rk_c*dt*alpha/dz(kend3)/dz_t(kend3)
680  c_src(1,1,nz)=c_src(1,1,nz)-2.d0*cc_t_zend*ttop
681  endif
682 
683  elseif (tbnd .eq. 0) then
684 
685  else
686  print *, 'tbnd flag error!'
687  stop
688  endif
689  ! check if it is u equation. We need some treatment for
690  ! the moving reference frame
691  if (ubnd .eq. 1) then
692  ! add the missing terms to qu for the Dirichlet condition
693  if (myid .eq. 0) then
694  ! applied boundary condition for u_mrf
695  aa_u_zstr=-0.5d0*rk_c*dt*nu/dz(kstr3)/dz_b(kstr3)
696  c_src(1,1,1)=c_src(1,1,1)+2.d0*aa_u_zstr*u_mrf
697  if (ochannel .eq. 0) then
698  cc_u_zend=-0.5d0*rk_c*dt*nu/dz(kend3)/dz_t(kend3)
699  c_src(1,1,nz)=c_src(1,1,nz)+2.d0*cc_u_zend*u_mrf
700  endif
701  endif
702 
703  elseif (ubnd .eq. 0) then
704 
705  else
706  print *, 'ubnd flag error!'
707  stop
708  endif
709 
710  ! Thomas algorithm to calculate c_var_out
711  ! forward loop
712  do k=kstr3+1,kmax3,1
713  do j=1,zsize(2),1
714  do i=1,zsize(1),1
715 
716  tmp=aa(i,j,k)/bb(i,j,k-1)
717  bb(i,j,k)=bb(i,j,k)-cc(i,j,k-1)*tmp
718  c_src(i,j,k-ghst)=c_src(i,j,k-ghst)-c_src(i,j,k-ghst-1)*tmp
719 
720  enddo
721  enddo
722  enddo
723  ! backward loop
724  do j=1,zsize(2),1
725  do i=1,zsize(1),1
726 
727  c_var_out(i,j,kmax3-ghst)=c_src(i,j,kmax3-ghst)/bb(i,j,kmax3)
728 
729  enddo
730  enddo
731 
732  do k=kmax3-1,kstr3,-1
733  do j=1,zsize(2),1
734  do i=1,zsize(1),1
735 
736  c_var_out(i,j,k-ghst)=( c_src(i,j,k-ghst)- &
737  cc(i,j,k)*c_var_out(i,j,k-ghst+1) ) / bb(i,j,k)
738 
739  enddo
740  enddo
741  enddo
742 
743 
744  return
745 
746  end subroutine
747 
748 
749 
750 
751  ! this function calculates the divergence
752  subroutine get_div_spec(c_u,c_v,w,c_div)
753 
754  ! these are temporary arrays, their dimensions are: nx/p_row,ny/p_col,nz
755  ! they are on z pencil
756  use module_variables, only : c_var
757  use mpi
758  use module_tools
759 
760  implicit none
761 
762  complex(mytype),dimension(:,:,:),intent(in) :: c_u,c_v
763  real(mytype),dimension(:,:,:),intent(in) :: w
764  complex(mytype),dimension(:,:,:),intent(out) :: c_div
765 
766  integer :: i,j,k
767  real(mytype) :: dwdz
768 
769 
770  ! add dudx
771  call dev_fft_linear_cmplx(c_u,c_var,1,1)
772  c_div(:,:,:)=c_var(:,:,:)
773 
774  ! add dvdy
775  call dev_fft_linear_cmplx(c_v,c_var,1,2)
776  c_div(:,:,:)=c_div(:,:,:)+c_var(:,:,:)
777 
778  ! add dwdz
779  do k=kstr3,kend3,1
780  do j=jstr3,jend3,1
781  do i=istr3,iend3,1
782 
783  dwdz = ( w(i,j,k)-w(i,j,k-1) )/dz(k)
784  c_var(i-ghst,j-ghst,k-ghst)= cmplx( dwdz,0.d0, kind=mytype )
785 
786  enddo
787  enddo
788  enddo
789  call mpi_2d_fft(c_var,-1)
790  c_div(:,:,:)=c_div(:,:,:)+c_var(:,:,:)
791 
792 
793  return
794 
795  end subroutine
796 
797 
798 
799 
800  ! calculate u^(n+1) using phi
801  subroutine calculate_uvw_spec(phi,c_phi,c_u,c_v,c_w)
802 
803  ! temporary complex array to store derivatives
804  use module_variables, only : c_var
805  use module_tools
806 
807  implicit none
808 
809  real(mytype),dimension(:,:,:),intent(in) :: phi
810  complex(mytype),dimension(:,:,:),intent(in) :: c_phi
811  complex(mytype),dimension(:,:,:),intent(inout) :: c_u,c_v,c_w
812 
813  integer :: i,j,k
814  real(mytype) :: dphidz
815 
816  ! calculate new u
817  call dev_fft_linear_cmplx(c_phi,c_var,1,1)
818  c_u(:,:,:)=c_u(:,:,:) - rk_c*dt*c_var(:,:,:)
819 
820  ! calculate new v
821  call dev_fft_linear_cmplx(c_phi,c_var,1,2)
822  c_v(:,:,:)=c_v(:,:,:) - rk_c*dt*c_var(:,:,:)
823 
824  ! calculate new w
825  do k=kstr3,kend3-1,1 ! note that we don't need to update w at nz
826  do j=jstr3,jend3,1
827  do i=istr3,iend3,1
828 
829  dphidz= ( phi(i,j,k+1) - phi(i,j,k) ) / dz_t(k)
830  c_var(i-ghst,j-ghst,k-ghst)= cmplx( dphidz,0.d0, kind=mytype )
831 
832  enddo
833  enddo
834  enddo
835  call mpi_2d_fft(c_var,-1)
836 
837  ! we don't update w(nz), we will assign it to zeros when applying
838  ! boundary conidition
839  do k=1,zsize(3)-1,1 ! notice this range!
840  do j=1,zsize(2),1
841  do i=1,zsize(1),1
842  c_w(i,j,k)=c_w(i,j,k) - rk_c*dt*c_var(i,j,k)
843  enddo
844  enddo
845  enddo
846 
847  return
848 
849  end subroutine
850 
851 
852 
853  ! calculate p^(n+1) using p^n and phi
854  subroutine calculate_new_p_spec(phi,c_phi,c_p)
855 
856  ! temporary complex array to store derivatives
857  use module_variables, only : c_var
858  use module_tools
859 
860  implicit none
861 
862  real(mytype),dimension(:,:,:),intent(in) :: phi
863  complex(mytype),dimension(:,:,:),intent(in) :: c_phi
864  complex(mytype),dimension(:,:,:),intent(inout) :: c_p
865 
866  integer :: i,j,k
867  real(mytype) :: dphidz2
868 
869  ! add phi
870  c_p(:,:,:)=c_p(:,:,:) + c_phi(:,:,:)
871  ! add dphidx2
872  call dev_fft_linear_cmplx(c_phi,c_var,2,1)
873  c_p(:,:,:)=c_p(:,:,:) - rk_c*0.5d0*dt*nu*c_var(:,:,:)
874  ! add dphidy2
875  call dev_fft_linear_cmplx(c_phi,c_var,2,2)
876  c_p(:,:,:)=c_p(:,:,:) - rk_c*0.5d0*dt*nu*c_var(:,:,:)
877  ! add dphidz2
878  do k=kstr3,kend3,1
879  do j=jstr3,jend3,1
880  do i=istr3,iend3,1
881 
882  dphidz2= ( (phi(i,j,k+1)-phi(i,j,k))/dz_t(k) - &
883  (phi(i,j,k)-phi(i,j,k-1))/dz_b(k) ) / dz(k)
884  c_var(i-ghst,j-ghst,k-ghst)= cmplx( dphidz2,0.d0, kind=mytype )
885 
886  enddo
887  enddo
888  enddo
889  call mpi_2d_fft(c_var,-1)
890  c_p(:,:,:)=c_p(:,:,:) - rk_c*0.5d0*dt*nu*c_var(:,:,:)
891 
892  return
893 
894  end subroutine
895 
896 
897 
898  ! print the div info to screen
899  subroutine screen_div_error_spec(u,v,w,c_u,c_v,isfinal,ishistout)
901  use mpi
902  use module_tools
903  ! temporary complex array to store derivatives
904  use module_variables, only : c_src
905 
906  implicit none
907 
908  real(mytype),dimension(:,:,:),intent(in) :: u,v,w
909  complex(mytype),dimension(:,:,:),intent(in) :: c_u,c_v
910 
911  integer,intent(in) :: isfinal,ishistout
912 
913  integer :: i,j,k
914  real(mytype) :: div_avg, div_max
915  real(mytype) :: div_avg_all,div_max_all
916  real(mytype) :: div_hist_max_all,tmp1
917 
918  call get_div_spec(c_u,c_v,w,c_src)
919 
920  call mpi_2d_fft(c_src,1)
921 
922  ! calculate continuity error
923  div_avg=0.d0
924  div_max=0.d0
925  do k=kstr3,kend3,1
926  do j=jstr3,jend3,1
927  do i=istr3,iend3,1
928 
929  tmp1=real(c_src(i-ghst,j-ghst,k-ghst),kind=mytype)
930 
931  div_avg=div_avg+abs(tmp1)
932 
933  if ( abs(tmp1) .gt. div_max ) then
934  div_max=abs(tmp1)
935  endif
936 
937  if ( u(i,j,k) .ne. u(i,j,k) .or. &
938  v(i,j,k) .ne. v(i,j,k) .or. &
939  w(i,j,k) .ne. w(i,j,k) ) then
940  print *, 'NaN error! ',i,' ',j,' ',k
941  stop
942  endif
943 
944  if (isfinal .eq. 1) then
945  if ( abs(tmp1) .gt. div_hist_max ) then
946  div_hist_max=abs(tmp1)
947  endif
948  endif
949 
950  enddo
951  enddo
952  enddo
953  div_avg=div_avg/nx/ny/nz
954 
955  ! exchange div_avg div_max between processors
956  ! sum up all the div_avg between all processors
957  call mpi_reduce(div_avg,div_avg_all,1,real_type,mpi_sum,0, &
958  mpi_comm_world,ierr)
959  ! find the max div_max between all processors
960  call mpi_reduce(div_max,div_max_all,1,real_type,mpi_max,0, &
961  mpi_comm_world,ierr)
962  ! find the max div_hist_max between all processors
963  call mpi_reduce(div_hist_max,div_hist_max_all,1,real_type,mpi_max,0, &
964  mpi_comm_world,ierr)
965 
966  ! print div to screen
967  if (myid .eq. 0) then
968 
969  if ( isfinal .eq. 0) then
970 
971  write (*,104) div_avg_all,div_max_all
972  104 format(' DIV-ERR INIT : DIV_AVG = ',es11.4,'; DIV_MAX = ',es11.4)
973 
974  elseif ( isfinal .eq. 1 .and. ishistout .eq. 0) then
975 
976  write (*,105) div_avg_all,div_max_all
977  105 format(' DIV-ERR FINAL: DIV_AVG = ',es11.4,'; DIV_MAX = ',es11.4)
978 
979  elseif ( isfinal .eq. 1 .and. ishistout .eq. 1) then
980 
981  write (*,115) div_avg_all,div_max_all
982  115 format(' DIV-ERR FINAL: DIV_AVG = ',es11.4,'; DIV_MAX = ',es11.4)
983  write (*,106) div_hist_max_all
984  106 format(' DIV-ERR HIST_MAX = ', es11.4)
985 
986  ! print ri if is_ri_var=1
987  if (is_ri_var .eq. 1) then
988  write (*,258) got
989  258 format(' Ri_b = ', es12.5)
990  endif
991  ! print bulk velocity if dp_opt=2
992  if (dp_opt .eq. 2) then
993  write (*,259) ubar+u_mrf
994  259 format(' Ubar = ', es12.5)
995  endif
996 
997  else
998  print *,'div out flag error!'
999  stop
1000  endif
1001 
1002  endif
1003 
1004 
1005  return
1006 
1007  end subroutine
1008 
1009 
1010 
1011  ! calculate mean fields
1012  subroutine calculate_mean_fields_spec &
1013  (u,v,w,t,mean_u,mean_v,mean_w,mean_t,mean_uu,mean_vv,mean_ww, &
1014  mean_uw,mean_vw,mean_tt,mean_tw)
1016  implicit none
1017 
1018  real(mytype),dimension(:,:,:),intent(in) :: u,v,w,t
1019  real(mytype),dimension(:,:,:),intent(out) :: mean_u,mean_v,mean_w,mean_t, &
1020  mean_uu,mean_vv,mean_ww,mean_uw, mean_vw, mean_tt,mean_tw
1021 
1022  integer :: i,j,k
1023 
1024  do k=kstr3,kend3,1
1025  do j=jstr3,jend3,1
1026  do i=istr3,iend3,1
1027 
1028  ! note that here we need to add the speed of moving reference frame
1029  mean_uu(i,j,k)=mean_uu(i,j,k)+(u(i,j,k)+u_mrf)**2.d0
1030  mean_vv(i,j,k)=mean_vv(i,j,k)+v(i,j,k)**2.d0
1031  mean_ww(i,j,k)=mean_ww(i,j,k)+w(i,j,k)**2.d0
1032  mean_uw(i,j,k)=mean_uw(i,j,k)+0.5d0*( w(i,j,k)+w(i,j,k-1) ) * &
1033  ( u(i,j,k)+u_mrf )
1034  mean_vw(i,j,k)=mean_vw(i,j,k)+0.5d0*( w(i,j,k)+w(i,j,k-1) ) * v(i,j,k)
1035 
1036  mean_tt(i,j,k)=mean_tt(i,j,k)+t(i,j,k)**2.d0
1037  mean_tw(i,j,k)=mean_tw(i,j,k)+t(i,j,k) * 0.5d0 * ( w(i,j,k)+w(i,j,k-1) )
1038 
1039  ! note that here we need to add the speed of moving reference frame
1040  mean_u(i,j,k)=mean_u(i,j,k) + u(i,j,k)+u_mrf
1041  mean_v(i,j,k)=mean_v(i,j,k) + v(i,j,k)
1042  mean_w(i,j,k)=mean_w(i,j,k) + w(i,j,k)
1043  mean_t(i,j,k)=mean_t(i,j,k) + t(i,j,k)
1044 
1045  enddo
1046  enddo
1047  enddo
1048 
1049 
1050  return
1051 
1052  end subroutine
1053 
1054 
1055 
1056  ! this function calculates derivatives for non-linear terms using spectral method
1057  ! dir=1: derivative in x direction
1058  ! dir=2: derivative in y direction
1059  ! weqn=1: for w equation. since w is staggered, we need special treatment
1060  ! if weqn=1, make sure to put w as r_in2
1061  subroutine dev_fft_nonlinear_cmplx(r_in1,r_in2,c_out3,dir,weqn)
1063  use module_tools
1064 
1065  implicit none
1066 
1067  real(mytype),dimension(:,:,:),intent(in) :: r_in1,r_in2
1068  complex(mytype),dimension(:,:,:),intent(out) :: c_out3
1069  integer,intent(in) :: dir,weqn
1070 
1071  integer :: i,j,k, l, m
1072  real(mytype) :: tmp1
1073 
1074  ! multiply these truncated variables
1075  do k=1,zsize(3),1
1076  do j=1,zsize(2),1
1077  do i=1,zsize(1),1
1078 
1079  if (weqn .eq. 1) then ! we need to interpolate u or v to w level
1080 
1081  tmp1=r_in1(i+ghst,j+ghst,k+1+ghst)*l_t(k+ghst) + &
1082  r_in1(i+ghst,j+ghst,k+ghst)*(1.d0-l_t(k+ghst))
1083 
1084  c_out3(i,j,k)=cmplx( tmp1 * &
1085  r_in2(i+ghst,j+ghst,k+ghst),0.d0,kind=mytype )
1086 
1087  else
1088 
1089  c_out3(i,j,k)=cmplx( r_in1(i+ghst,j+ghst,k+ghst) * &
1090  r_in2(i+ghst,j+ghst,k+ghst),0.d0,kind=mytype )
1091 
1092  endif
1093 
1094 
1095 
1096  enddo
1097  enddo
1098  enddo
1099 
1100  ! mpi 2d fft for c_out3
1101  call mpi_2d_fft(c_out3,-1)
1102 
1103  ! calculate derivative of the non-linear term
1104  do k=1,zsize(3),1
1105  do j=1,zsize(2),1
1106  do i=1,zsize(1),1
1107 
1108  ! deal with index
1109  ! convert i to l
1110  if ( (i+i_offset) .le. nx/2) then
1111  l=i+i_offset-1
1112  else
1113  l=i+i_offset-nx-1
1114  endif
1115  ! convert j to m
1116  if ( (j+j_offset) .le. ny/2) then
1117  m=j+j_offset-1
1118  else
1119  m=j+j_offset-ny-1
1120  endif
1121 
1122  if (dir .eq. 1) then ! d?dx
1123 
1124  c_out3(i,j,k) = c_out3(i,j,k)*iii*2.d0*pi/lx*l
1125 
1126  elseif (dir .eq. 2) then ! d?dy
1127 
1128  c_out3(i,j,k) = c_out3(i,j,k)*iii*2.d0*pi/ly*m
1129 
1130  else
1131 
1132  print *, 'flag error'
1133  stop
1134 
1135  endif
1136 
1137  ! remove Nyquist wavenumber
1138  if (l .eq. -nx/2 .or. m .eq. -ny/2) then
1139  c_out3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
1140  endif
1141 
1142 
1143  enddo
1144  enddo
1145  enddo
1146 
1147 
1148  return
1149 
1150  end subroutine
1151 
1152 
1153 
1154 
1155 
1156  ! do fft and calculate Fourier coefficients for a variable
1157  subroutine real_2_cmplx(r_in,c_out)
1159  use module_tools
1160 
1161  implicit none
1162 
1163  real(mytype),dimension(:,:,:),intent(in) :: r_in
1164  complex(mytype),dimension(:,:,:),intent(out) :: c_out
1165 
1166  integer :: i,j,k
1167 
1168  do k=1,zsize(3),1
1169  do j=1,zsize(2),1
1170  do i=1,zsize(1),1
1171 
1172  c_out(i,j,k)=cmplx(r_in(i+ghst,j+ghst,k+ghst),0.d0,kind=mytype)
1173 
1174  enddo
1175  enddo
1176  enddo
1177 
1178  call mpi_2d_fft(c_out,-1)
1179 
1180  return
1181 
1182  end subroutine
1183 
1184 
1185 
1186 
1187  ! do ifft and calculate a variable based on the Fourier coefficients
1188  subroutine cmplx_2_real(c_in,r_out)
1190  ! temporary complex array, their dimension is: nx/p_row,ny/p_col,nz
1191  use module_variables, only : c_var
1192  use module_tools
1193 
1194  implicit none
1195 
1196  real(mytype),dimension(:,:,:),intent(out) :: r_out
1197  complex(mytype),dimension(:,:,:),intent(in) :: c_in
1198 
1199  integer :: i,j,k
1200 
1201  c_var(:,:,:)=c_in(:,:,:)
1202 
1203  call mpi_2d_fft(c_var,1)
1204 
1205  do k=1,zsize(3),1
1206  do j=1,zsize(2),1
1207  do i=1,zsize(1),1
1208 
1209  r_out(i+ghst,j+ghst,k+ghst)=real(c_var(i,j,k),kind=mytype)
1210 
1211  enddo
1212  enddo
1213  enddo
1214 
1215  return
1216 
1217  end subroutine
1218 
1219 
1220 
1221  ! truncate the Fourier coefficients using the 2/3 rule
1222  subroutine spectral_truncation(c_inout3)
1224  implicit none
1225 
1226  complex(mytype),dimension(:,:,:),intent(inout) :: c_inout3
1227 
1228  integer :: i,j,k,l,m
1229 
1230  do k=1,zsize(3),1
1231  do j=1,zsize(2),1
1232  do i=1,zsize(1),1
1233 
1234  ! deal with index
1235  ! convert i to l
1236  if ( (i+i_offset) .le. nx/2) then
1237  l=i+i_offset-1
1238  else
1239  l=i+i_offset-nx-1
1240  endif
1241  ! convert j to m
1242  if ( (j+j_offset) .le. ny/2) then
1243  m=j+j_offset-1
1244  else
1245  m=j+j_offset-ny-1
1246  endif
1247 
1248  if (l .ge. int(nx/3.d0) .or. l .le. -int(nx/3.d0) ) then
1249 
1250  c_inout3(i,j,k) = cmplx(0.d0,0.d0,kind=mytype)
1251 
1252  endif
1253 
1254 
1255  if (m .ge. int(ny/3.d0) .or. m .le. -int(ny/3.d0) ) then
1256 
1257  c_inout3(i,j,k) = cmplx(0.d0,0.d0,kind=mytype)
1258 
1259  endif
1260 
1261 
1262  enddo
1263  enddo
1264  enddo
1265 
1266  return
1267 
1268  end subroutine
1269 
1270 
1271 
1272  ! this function calculate derivatives for linear terms using spectral method
1273  ! ord: 1st or 2nd order derivatives; dir: direction? x or y?
1274  subroutine dev_fft_linear_cmplx(c_in3,c_out3,ord,dir)
1276  implicit none
1277 
1278  complex(mytype),dimension(:,:,:),intent(in) :: c_in3
1279  complex(mytype),dimension(:,:,:),intent(out) :: c_out3
1280  integer,intent(in) :: ord,dir
1281 
1282  integer :: i,j,k,l,m
1283 
1284  c_out3(:,:,:)=c_in3(:,:,:)
1285 
1286  ! calculate derivative
1287  do k=1,zsize(3),1
1288  do j=1,zsize(2),1
1289  do i=1,zsize(1),1
1290 
1291  ! deal with index
1292  ! convert i to l
1293  if ( (i+i_offset) .le. nx/2) then
1294  l=i+i_offset-1
1295  else
1296  l=i+i_offset-nx-1
1297  endif
1298  ! convert j to m
1299  if ( (j+j_offset) .le. ny/2) then
1300  m=j+j_offset-1
1301  else
1302  m=j+j_offset-ny-1
1303  endif
1304 
1305  if (ord .eq. 1 .and. dir .eq. 1) then ! d?dx
1306 
1307  c_out3(i,j,k) = c_out3(i,j,k)*iii*2.d0*pi/lx*l
1308 
1309  elseif (ord .eq. 1 .and. dir .eq. 2) then ! d?dy
1310 
1311  c_out3(i,j,k) = c_out3(i,j,k)*iii*2.d0*pi/ly*m
1312 
1313  elseif (ord .eq. 2 .and. dir .eq. 1) then ! d?dx2
1314 
1315  c_out3(i,j,k) = -c_out3(i,j,k)*(2.d0*pi/lx*l)**2.d0
1316 
1317  elseif (ord .eq. 2 .and. dir .eq. 2) then ! d?dy2
1318 
1319  c_out3(i,j,k) = -c_out3(i,j,k)*(2.d0*pi/ly*m)**2.d0
1320 
1321  else
1322 
1323  print *, 'flag error'
1324  stop
1325 
1326  endif
1327 
1328  ! remove Nyquist wavenumber
1329  if (l .eq. -nx/2 .or. m .eq. -ny/2) then
1330  c_out3(i,j,k)=cmplx(0.d0,0.d0,kind=mytype)
1331  endif
1332 
1333  enddo
1334  enddo
1335  enddo
1336 
1337 
1338  return
1339 
1340  end subroutine
1341 
1342 
1343 
1344  ! Rayleigh damping layer at the top of the boundary layer
1345  subroutine rayleigh_damping_spec(c_varin)
1347  use mpi
1348 
1349  implicit none
1350 
1351  complex(mytype),dimension(:,:,:),intent(inout) :: c_varin
1352 
1353  real(mytype) :: zd,sigz
1354  complex(mytype) :: mean_local,mean_all
1355  integer :: i,j,k
1356 
1357 
1358  ! do the damping for the topmost nzdamp points
1359  do k=kend3-nzdamp,kend3,1
1360 
1361  ! first calculate the mean of the variable on the local processor
1362  mean_local=0.d0
1363  do j=jstr3,jend3,1
1364  do i=istr3,iend3,1
1365  mean_local=mean_local+c_varin(i-ghst,j-ghst,k-ghst)
1366  enddo
1367  enddo
1368  mean_local=mean_local/zsize(1)/zsize(2)
1369 
1370  ! exchange mean_local between processors
1371  ! sum up all the mean_local between all processors
1372  call mpi_allreduce(mean_local,mean_all,1,real_type,mpi_sum, &
1373  mpi_comm_world,ierr)
1374  mean_all=mean_all/nprc ! calculate mean all
1375 
1376  ! do the damping
1377  zd=(k-kend3+nzdamp)*1.d0/nzdamp
1378  sigz=cdamp*0.5d0*( 1.d0-cos(pi*zd) )
1379  do j=jstr3,jend3,1
1380  do i=istr3,iend3,1
1381  c_varin(i-ghst,j-ghst,k-ghst) = c_varin(i-ghst,j-ghst,k-ghst) - &
1382  sigz*( c_varin(i-ghst,j-ghst,k-ghst)-mean_all )
1383  enddo
1384  enddo
1385 
1386  enddo
1387 
1388 
1389  return
1390 
1391  end subroutine
1392 
1393 
1394 end module
integer, save, protected is_ri_var
integer, save, protected myid
real(mytype), save, protected rk_b
real(mytype), dimension(:), allocatable, save, protected l_b
integer, save, protected ochannel
integer, save, protected nprc
real(mytype), save, protected ubar
integer, save, protected ny
subroutine mpi_2d_fft(c_inout3, direction)
real(mytype), save, protected fc
subroutine real_2_cmplx(r_in, c_out)
integer, save, protected issk
complex(mytype), dimension(:,:,:), allocatable, save c_src
real(mytype), parameter pi
real(mytype), save, protected dt
subroutine get_temperature_cnvdiff_spec(u, v, w, t, c_t, c_ss, c_qt)
real(mytype), save, protected rk_a
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ttop
integer, save, protected nz
real(mytype), save, protected ly
real(mytype), save, protected ug
subroutine calculate_mean_fields_spec(u, v, w, t, mean_u, mean_v, mean_w, mean_t, mean_uu, mean_vv, mean_ww, mean_uw, mean_vw, mean_tt, mean_tw)
subroutine spectral_truncation(c_inout3)
real(mytype), dimension(:), allocatable, save, protected dz
integer, save, protected isscalar
real(mytype), save, protected cdamp
real(mytype), save, protected alpha
integer, save, protected istr3
real(mytype), save, protected vg
real(mytype), save, protected t_ref
subroutine get_momentum_cnvdiff_spec(u, v, w, p, t, c_u, c_v, c_w, c_p, c_ff, c_gg, c_hh, c_qu, c_qv, c_qw)
subroutine rayleigh_damping_spec(c_varin)
integer, save, protected jstr3
real(mytype), dimension(:), allocatable, save, protected dz_b
integer, save, protected i_offset
integer, save, protected j_offset
subroutine cmplx_2_real(c_in, r_out)
real(mytype), save, protected u_mrf
complex(mytype), parameter iii
integer, parameter ghst
real(mytype), save, protected lx
real(mytype), save, protected dpx_drive
integer, save, protected nzdamp
real(mytype), dimension(:), allocatable, save, protected dz_t
subroutine poisson_solver_fft_spec(aa, bb, cc, c_src_in, c_var_out, zstagger, delnull, tbnd, ubnd)
integer, save, protected nx
subroutine dev_fft_linear_cmplx(c_in3, c_out3, ord, dir)
subroutine screen_div_error_spec(u, v, w, c_u, c_v, isfinal, ishistout)
real(mytype), save, protected got
subroutine get_div_spec(c_u, c_v, w, c_div)
real(mytype), save div_hist_max
subroutine calculate_new_p_spec(phi, c_phi, c_p)
subroutine calculate_uvw_spec(phi, c_phi, c_u, c_v, c_w)
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), dimension(:), allocatable, save, protected l_t
real(mytype), save, protected tbot
complex(mytype), dimension(:,:,:), allocatable, save c_var
integer, save, protected dp_opt
subroutine dev_fft_nonlinear_cmplx(r_in1, r_in2, c_out3, dir, weqn)
integer, save, protected kstr3
integer, save, protected kend3