module_navier_stokes.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 contains
28 
29 
30  ! the main function for HERCULES, it solves the Navier-Stokes for incompressible flows
31  subroutine drive
32 
34  use module_io
35  use module_tools
36  use module_boundary
38  use mpi
39  use module_spectral
40 
41  implicit none
42 
43  ! for using fftw
44  include 'include/fftw3.f03'
45 
46  ! indx_bk/m/ins are indices for output
47  ! backup/mean/instantaneous variables
48  ! indx_slc: indices for output slice data
49  integer :: indx_bk,indx_m,indx_ins,indx_slc
50  integer :: stps
51 
52  ! calculate mesh parameters, e.g., dz
53  call get_mesh_param
54 
55  ! allocate arrays for all variables
56  ! initiate them using default values
57  ! or read them from files (e.g., init_u.dat)
58  call initiate_fields
59 
60  ! assign initial ri if the varying ri is used, i.e., is_ri_var=1
61  if (is_ri_var .eq. 1) then
62  call update_ri(0)
63  endif
64 
65  ! calculate bulk velocity if dp_opt=2
66  if (dp_opt .eq. 2) then
67  call calculate_ubar(u)
68  endif
69 
70  ! create fftw plans, this needs to be done only once
71  call get_fft_plan
72 
73  ! get the coefficient a,b,c for the u, v, w, t, and p equations
74  ! we store them in aa_u,bb_u, etc.,
75  ! and then assign them to aa,bb,cc using subroutine assign_abc
80 
81  ! open new files for outputing time-series
82  if (myid .eq. 0 .and. istmsr .eq. 1) then
83  open(99,file='results/tmsr_u.dat',status='replace')
84  close(99)
85  open(98,file='results/tmsr_v.dat',status='replace')
86  close(98)
87  open(97,file='results/tmsr_w.dat',status='replace')
88  close(97)
89  open(96,file='results/tmsr_t.dat',status='replace')
90  close(96)
91  open(95,file='results/tmsr_p.dat',status='replace')
92  close(95)
93  endif
94 
95  ! for spectral method, we need to calculate the Fourier coefficients
96  ! and do truncation, and update fields in physical space
97  ! note that for spectral method, only the 2/3 truncation dealiasing
98  ! is avaiable right now
99  if (cds .eq. 3) then
100 
101  ! convert variables from physical space to phase space
102  call real_2_cmplx(u,c_u)
103  call real_2_cmplx(v,c_v)
104  call real_2_cmplx(w,c_w)
105  call real_2_cmplx(p,c_p)
106  call real_2_cmplx(t,c_t)
107  call real_2_cmplx(ff,c_ff)
108  call real_2_cmplx(gg,c_gg)
109  call real_2_cmplx(hh,c_hh)
110  call real_2_cmplx(ss,c_ss)
111  !!!!!!!!!!! do spectral truncation !!!!!!!!!!!!
117  !!!!!!!!!! get updated variables after truncation !!!!!!!!
118  call cmplx_2_real(c_u,u)
119  call cmplx_2_real(c_v,v)
120  call cmplx_2_real(c_w,w)
121  call cmplx_2_real(c_p,p)
122  call cmplx_2_real(c_t,t)
123 
124  endif
125 
126  ! update boundary values for u,v,w, t, and p
127  call update_boundary_uvw(u,v,w)
128  call update_boundary_p(p)
129  call update_boundary_t(t)
130 
131 
132  ! indices to save backup/mean/inst fields, or slice data
133  indx_bk=1
134  indx_m=1
135  indx_ins=1
136  indx_slc=1
137 
138  ! main loop, we use fractional time step approach for time advancement
139  do stps=1,imax,1
140 
141  ! print cpu time and CFL number
142  call screen_cpu_time(u,v,w)
143 
144  ! calculate new ri if is_ri_var=1
145  if (is_ri_var .eq. 1) then
146  call update_ri(stps)
147  endif
148 
149  ! do the time advancement
150  if (cds .eq. 3) then
151 
152  ! subroutine for spectral method
154 
155  else
156 
157  ! subroutine for finite-difference method
159 
160  endif
161 
162  ! update bulk velocity if dp_opt=2
163  if (dp_opt .eq. 2) then
164  call calculate_ubar(u)
165  endif
166 
167  !!!!!!!!!!!!!!!!! output data !!!!!!!!!!!!!!!!
168  ! output time-series
169  if (istmsr .eq. 1) then
170  ! output time-series at the channel center
171  if (myid_rowindx .eq. p_row/2+1 .and. &
172  myid_colindx .eq. p_col/2+1 ) then
173  call output_time_series(u,v,w,p,t)
174  endif
175  endif
176 
177  ! write backup files
178  if (indx_bk/ibackup .eq. 1 ) then
179  call output_backup(u,v,w,p,t,ff,gg,hh,ss)
180  indx_bk=0
181  endif
182 
183  ! write mean fields
184  if (indx_m/imeanfl .eq. 1) then
187  mean_tw,stps)
188  indx_m=0
189  endif
190 
191  ! write inst fields
192  if (indx_ins/iinstfl .eq. 1) then
193  call output_inst_fields(u,v,w,p,t,stps)
194  indx_ins=0
195  endif
196 
197  ! output 2d slices
198  if (isxy2d .eq. 1 .or. isxz2d .eq. 1 .or. isyz2d .eq. 1) then
199  if (indx_slc/intv_2d .eq. 1) then
200  call output_2d_slices(u,v,w,p,t,stps)
201  indx_slc=0
202  endif
203  endif
204 
205  indx_bk=indx_bk+1
206  indx_m=indx_m+1
207  indx_ins=indx_ins+1
208  indx_slc=indx_slc+1
209 
210  enddo
211 
212  ! destroy fftw plans
213  call fftw_destroy_plan(fft_plan1)
214  call fftw_destroy_plan(ifft_plan1)
215  call fftw_destroy_plan(fft_plan2)
216  call fftw_destroy_plan(ifft_plan2)
217 
218  end subroutine
219 
220 
221 
222  ! this function does the time-advancement using finite-difference method
223  subroutine time_advancement_fd
225  use module_variables
226  use module_tools
227  use module_boundary
228  use module_cnvdiff
230  use module_io
231  use mpi
232 
233  implicit none
234 
235  integer :: t_its,t_its_max
236 
237  ! if it is RK3 then we need to do 3 intermediate steps
238  ! for fractional time step. While for AB2, 1 step is enough
239  if (dts .eq. 2) then
240  t_its_max=3
241  else
242  t_its_max=1
243  endif
244 
245  ! for RK3 we need to do 3 intermediate sub-steps for each time step
246  ! for AB2, t_its_max=1
247  do t_its=1,t_its_max,1
248 
249  ! if it is RK3 scheme, write some info to screen
250  if (dts .eq. 2 .and. myid .eq. 0) then
251 
252  write (*,218) t_its
253  218 format(' RK3, Step: ',i2)
254 
255  endif
256 
257  ! check if we need to do rayleigh damping at the top
258  if (isdamp .eq. 1) then
259 
260  call rayleigh_damping(u)
261  call rayleigh_damping(v)
262  call rayleigh_damping(w)
263  call rayleigh_damping(p)
264  call rayleigh_damping(t)
265 
266  endif
267 
268  ! before calculating ff, gg, hh, assign them to the previous time step
269  ! ff,gg,hh are used in get_momentum_cnvdiff
271 
272  ! calculate interpolated fields,
273  ! which will be used in cnvdiff calculation
276 
277  ! assign coefficients for RK3 and AB2 scheme
278  ! these coefficients, i.e., rk_a,rk_b,rk_c, will be used
279  ! in get_momentum_cnvdiff and assign_abc
280  call assign_rk_coeff(rkc1(t_its),rkc2(t_its),rkc3(t_its))
281 
282  ! calculate the spatial derivatives
283  ! and get the right hand side for u, v, and w eqns
284  call get_momentum_cnvdiff &
286 
287  ! fractional time step
288  ! first step, calculating u* using p^n
289 
290  ! assign the previously calculated aa_u,bb_u,cc_u, to aa,bb,cc
291  call assign_abc(aa_u,bb_u,cc_u,aa,bb,cc)
292  ! solve the u equation using fft solver
293  call poisson_solver_fft(aa,bb,cc,qu,u,0,0,0,1)
294 
295  ! assign the previously calculated aa_v,bb_v,cc_v, to aa,bb,cc
296  ! coeffs a,b,c for v eqn is same with u eqn,
297  call assign_abc(aa_u,bb_u,cc_u,aa,bb,cc)
298  ! solve the v equation using fft solver
299  call poisson_solver_fft(aa,bb,cc,qv,v,0,0,0,0)
300 
301  ! assign the previously calculated aa_w,bb_w,cc_w, to aa,bb,cc
302  call assign_abc(aa_w,bb_w,cc_w,aa,bb,cc)
303  ! solve the w equation using fft solver
304  call poisson_solver_fft(aa,bb,cc,qw,w,1,0,0,0)
305  ! update the boundary values
306  call update_boundary_uvw(u,v,w)
307 
308  ! based on u*, calculate the pseudo pressure phi
309  ! here we reuse qu for the source term of p eqn to save memory usage
310  call get_p_eqn_src(u,v,w,qu)
311  ! assign the previously calculated aa_p,bb_p,cc_p, to aa,bb,cc
313  ! solve the p equation using fft solver
314  call poisson_solver_fft(aa,bb,cc,qu,phi,0,1,0,0)
315  ! update boundary values
316  call update_boundary_p(phi)
317 
318  ! calculate u^(n+1) using phi
319  call calculate_uvw(u,v,w,phi)
320  call update_boundary_uvw(u,v,w)
321 
322  ! calculate temperature
323  if (isscalar .eq. 1) then
324 
325  ! assign ss to previous steps
326  call assign_s_old(ss,ss0)
327 
328  ! calculate interpolated fields
331 
332  ! here we reuse qu for the source term of t eqn to save memory
334  ! assign the previously calculated aa_t,bb_t,cc_t, to aa,bb,cc
335  call assign_abc(aa_t,bb_t,cc_t,aa,bb,cc)
336  ! solve the t equation using fft solver
337  call poisson_solver_fft(aa,bb,cc,qu,t,0,0,1,0)
338  ! update boundary values
339  call update_boundary_t(t)
340 
341  endif
342 
343  ! calculate the final DIV error
344  if (dts .ne. 2) then
345  call screen_div_error(u,v,w,1,1)
346  elseif (t_its .ne. 3) then
347  call screen_div_error(u,v,w,1,0)
348  else
349  call screen_div_error(u,v,w,1,1)
350  endif
351 
352  ! calculate p^{n+1} using phi
353  call calculate_new_p(phi,p)
354  call update_boundary_p(p)
355 
356  enddo
357 
358  ! calculate mean fields
359  call calculate_mean_fields &
362 
363  return
364 
365  end subroutine
366 
367 
368 
369  ! this function does the time-advancement
370  ! using spectral method in horizontal directions
371  subroutine time_advancement_sp
373  use module_spectral
374  use module_variables
375  use mpi
376  use module_tools
377  use module_boundary
378  use module_io
380 
381  implicit none
382 
383  integer :: t_its,t_its_max
384 
385  ! if it is RK3 then we need to do 3 intermediate steps
386  ! for fractional time step. While for AB2, 1 step is enough
387  if (dts .eq. 2) then
388  t_its_max=3
389  else
390  t_its_max=1
391  endif
392 
393  ! for RK3 we need to do 3 intermediate sub-steps for each time step
394  ! for AB2, t_its_max=1
395  do t_its=1,t_its_max,1
396 
397  ! if it is RK3 scheme, write some info to screen
398  if (dts .eq. 2 .and. myid .eq. 0) then
399 
400  write (*,218) t_its
401  218 format(' RK3, Step: ',i2)
402 
403  endif
404 
405  ! check if we need to do rayleigh damping
406  if (isdamp .eq. 1) then
407 
413 
414  endif
415 
416  ! assign the coefficient for RK3 or AB2
417  call assign_rk_coeff(rkc1(t_its),rkc2(t_its),rkc3(t_its))
418 
419  ! calculate the right-hand side terms
422 
423  ! fractional time step
424  ! first step, calculating u* using p^n
425 
426  ! assign the previously calculated aa_u,bb_u,cc_u, to aa,bb,cc
427  call assign_abc(aa_u,bb_u,cc_u,aa,bb,cc)
428  ! solve the u equation using fft solver
429  call poisson_solver_fft_spec(aa,bb,cc,c_qu,c_u,0,0,0,1)
430 
431  ! assign the previously calculated aa_v,bb_v,cc_v, to aa,bb,cc
432  ! coeffs a,b,c for v eqn is same with u eqn,
433  call assign_abc(aa_u,bb_u,cc_u,aa,bb,cc)
434  ! solve the v equation using fft solver
435  call poisson_solver_fft_spec(aa,bb,cc,c_qv,c_v,0,0,0,0)
436 
437  ! assign the previously calculated aa_w,bb_w,cc_w, to aa,bb,cc
438  call assign_abc(aa_w,bb_w,cc_w,aa,bb,cc)
439  ! solve the w equation using fft solver
440  call poisson_solver_fft_spec(aa,bb,cc,c_qw,c_w,1,0,0,0)
441 
442  ! we only need to update w in physical space, since we need to
443  ! calculate dwdz for the source term of p equation
444  call cmplx_2_real(c_w,w)
445  call update_boundary_uvw(u,v,w)
446 
447  ! based on u*, calculate the pseudo pressure phi
448  ! here we reused qu for the source term of p eqn to save memory usage
449  call get_div_spec(c_u,c_v,w,c_qu)
450  c_qu(:,:,:)=c_qu(:,:,:)/dt
451  ! assign the previously calculated aa_p,bb_p,cc_p, to aa,bb,cc
453  ! solve the p equation using fft solver
454  call poisson_solver_fft_spec(aa,bb,cc,c_qu,c_phi,0,1,0,0)
455 
456  ! we need to update phi in physical space for calculating new p
457  call cmplx_2_real(c_phi,phi)
458  call update_boundary_p(phi)
459 
460  ! calculate p^{n+1} using phi
462 
463  ! calculate u^(n+1) using phi
465 
466  !!!!!!!!!!! do spectral truncation !!!!!!!!!!!!
471 
472  ! update u,v,w,p in physical space
473  call cmplx_2_real(c_u,u)
474  call cmplx_2_real(c_v,v)
475  call cmplx_2_real(c_w,w)
476  call cmplx_2_real(c_p,p)
477  call update_boundary_p(p)
478  call update_boundary_uvw(u,v,w)
479 
480  ! calculate temperature
481  if (isscalar .eq. 1) then
482 
483  ! here we reuse qu for the source term of t eqn to save memory
485  (u,v,w,t,c_t,c_ss,c_qu)
486  ! assign the previously calculated aa_t,bb_t,cc_t, to aa,bb,cc
487  call assign_abc(aa_t,bb_t,cc_t,aa,bb,cc)
488  ! solve the t equation using fft solver
489  call poisson_solver_fft_spec(aa,bb,cc,c_qu,c_t,0,0,1,0)
490 
491  ! update t in physical space
493  call cmplx_2_real(c_t,t)
494  call update_boundary_t(t)
495 
496  endif
497 
498  ! calculate the final DIV error
499  if (dts .ne. 2) then
500  call screen_div_error_spec(u,v,w,c_u,c_v,1,1)
501  elseif (t_its .ne. 3) then
502  call screen_div_error_spec(u,v,w,c_u,c_v,1,0)
503  else
504  call screen_div_error_spec(u,v,w,c_u,c_v,1,1)
505  endif
506 
507 
508  enddo
509 
510  ! calculate mean fields
514 
515  return
516 
517  end subroutine
518 
519 
520 
521  ! this function assigns the previously
522  ! calculated aa_u, bb_u, cc_u, etc. to aa,bb,cc
523  subroutine assign_abc(aa_in,bb_in,cc_in,aa,bb,cc)
525  use decomp_2d
526 
527  implicit none
528 
529  real(mytype), dimension(:), intent(in) :: aa_in,cc_in
530  complex(mytype),dimension(:,:,:),intent(in) :: bb_in
531  real(mytype), dimension(:,:,:),intent(out) :: aa,cc
532  complex(mytype),dimension(:,:,:),intent(out) :: bb
533 
534  integer :: i,j,k
535 
536  ! note that in get_u_eqn_coeff, we did not add 1.0
537  ! to bb, here we add it. Also, if we use RK3, we need to
538  ! multiply rk_c to a,b,c, the rk_c for the current RK sub-step
539  ! has been assign previously by calling assign_rk_coeff
540  if (dts .eq. 2) then
541 
542  do k=1,zsize(3)+2*ghst,1
543  do j=1,zsize(2),1
544  do i=1,zsize(1),1
545 
546  aa(i,j,k)=rk_c*aa_in(k)
547  cc(i,j,k)=rk_c*cc_in(k)
548  bb(i,j,k)=1.d0+rk_c*bb_in(i,j,k)
549 
550  enddo
551  enddo
552  enddo
553 
554  else
555 
556  do k=1,zsize(3)+2*ghst,1
557  do j=1,zsize(2),1
558  do i=1,zsize(1),1
559 
560  aa(i,j,k)=aa_in(k)
561  cc(i,j,k)=cc_in(k)
562  bb(i,j,k)=1.d0+bb_in(i,j,k)
563 
564  enddo
565  enddo
566  enddo
567 
568  endif
569 
570  return
571 
572  end subroutine
573 
574 
575 
576  ! this function assigns the previously calculated
577  ! aa_p, bb_p, cc_p, etc. to aa,bb,cc for p equation
578  subroutine assign_abc_p(aa_in,bb_in,cc_in,aa,bb,cc)
580  use decomp_2d
581 
582  implicit none
583 
584  real(mytype), dimension(:), intent(in) :: aa_in,cc_in
585  complex(mytype),dimension(:,:,:),intent(in) :: bb_in
586  real(mytype), dimension(:,:,:),intent(out) :: aa,cc
587  complex(mytype),dimension(:,:,:),intent(out) :: bb
588 
589  integer :: i,j,k
590 
591  if (dts .eq. 2) then
592 
593  do k=1,zsize(3)+2*ghst,1
594  do j=1,zsize(2),1
595  do i=1,zsize(1),1
596  aa(i,j,k)=rk_c*aa_in(k)
597  cc(i,j,k)=rk_c*cc_in(k)
598  bb(i,j,k)=rk_c*bb_in(i,j,k)
599  enddo
600  enddo
601  enddo
602 
603  else
604 
605  do k=1,zsize(3)+2*ghst,1
606  do j=1,zsize(2),1
607  do i=1,zsize(1),1
608  aa(i,j,k)=aa_in(k)
609  cc(i,j,k)=cc_in(k)
610  bb(i,j,k)=bb_in(i,j,k)
611  enddo
612  enddo
613  enddo
614 
615  endif
616 
617  return
618 
619  end subroutine
620 
621 
622 
623  ! this function assigns ff to the previous time steps, e.g., ff0
624  subroutine assign_fgh_old(ff,gg,hh,ff0,gg0,hh0)
626  implicit none
627 
628  real(mytype),dimension(:,:,:),intent(in) :: ff,gg,hh
629  real(mytype),dimension(:,:,:),intent(inout) :: ff0,gg0,hh0
630 
631  ! assign ff etc to ff0 etc.
632  ff0(:,:,:)=ff(:,:,:)
633  gg0(:,:,:)=gg(:,:,:)
634  hh0(:,:,:)=hh(:,:,:)
635 
636  return
637 
638  end subroutine
639 
640 
641 
642  ! this function assigns ss to the previous time steps, e.g., ss0
643  subroutine assign_s_old(ss,ss0)
645  implicit none
646 
647  real(mytype),dimension(:,:,:),intent(in) :: ss
648  real(mytype),dimension(:,:,:),intent(out) :: ss0
649 
650  ss0(:,:,:)=ss(:,:,:)
651 
652  return
653 
654  end subroutine
655 
656 
657 
658  ! calculate u^(n+1) using phi
659  subroutine calculate_uvw(u,v,w,phi)
660 
661  implicit none
662 
663  real(mytype),dimension(:,:,:),intent(in) :: phi
664  real(mytype),dimension(:,:,:),intent(inout) :: u,v,w
665 
666  integer :: i,j,k
667 
668  ! check to use 2nd or 4th schemes
669  if ( cds .eq. 1) then
670 
671  do k=kstr3,kend3,1
672  do j=jstr3,jend3,1
673  do i=istr3,iend3,1
674 
675  u(i,j,k)= u(i,j,k) - rk_c*dt* &
676  ( phi(i+1,j,k) - phi(i,j,k) ) / dx
677 
678  enddo
679  enddo
680  enddo
681 
682  do k=kstr3,kend3,1
683  do j=jstr3,jend3,1
684  do i=istr3,iend3,1
685 
686  v(i,j,k)= v(i,j,k) - rk_c*dt* &
687  ( phi(i,j+1,k) - phi(i,j,k) ) / dy
688 
689  enddo
690  enddo
691  enddo
692 
693  elseif ( cds .eq. 2) then
694 
695  do k=kstr3,kend3,1
696  do j=jstr3,jend3,1
697  do i=istr3,iend3,1
698 
699  u(i,j,k)= u(i,j,k) - rk_c*dt*( -phi(i+2,j,k) + &
700  27.d0*phi(i+1,j,k) - &
701  27.d0*phi(i,j,k) + &
702  phi(i-1,j,k) ) / 24.d0 / dx
703  enddo
704  enddo
705  enddo
706 
707  do k=kstr3,kend3,1
708  do j=jstr3,jend3,1
709  do i=istr3,iend3,1
710 
711  v(i,j,k)= v(i,j,k) - rk_c*dt*( -phi(i,j+2,k) + &
712  27.d0*phi(i,j+1,k) - &
713  27.d0*phi(i,j,k) + &
714  phi(i,j-1,k) ) / 24.d0 / dy
715  enddo
716  enddo
717  enddo
718 
719  endif
720 
721  do k=kstr3,kend3-1,1 ! note this range is different
722  do j=jstr3,jend3,1
723  do i=istr3,iend3,1
724 
725  w(i,j,k)= w(i,j,k) - rk_c*dt* &
726  ( phi(i,j,k+1) - phi(i,j,k) ) / dz_t(k)
727 
728  enddo
729  enddo
730  enddo
731 
732  return
733 
734  end subroutine
735 
736 
737 
738  ! calculate p^(n+1) using p^n and the pseudo pressure phi
739  subroutine calculate_new_p(phi,p)
741  implicit none
742 
743  real(mytype),dimension(:,:,:),intent(in) :: phi
744  real(mytype),dimension(:,:,:),intent(out) :: p
745 
746  integer :: i,j,k
747  real(mytype) :: phixx,phiyy,phizz
748 
749  ! check to use 2nd or 4th schemes
750  if ( cds .eq. 1) then
751 
752  do k=kstr3,kend3,1
753  do j=jstr3,jend3,1
754  do i=istr3,iend3,1
755 
756  phixx = ( phi(i+1,j,k) - &
757  2.d0*phi(i,j,k) + &
758  phi(i-1,j,k) ) /dx2
759 
760  phiyy = ( phi(i,j+1,k) - &
761  2.d0*phi(i,j,k) + &
762  phi(i,j-1,k) ) /dy2
763 
764  phizz= ( (phi(i,j,k+1)-phi(i,j,k))/dz_t(k) - &
765  (phi(i,j,k)-phi(i,j,k-1))/dz_b(k) ) / dz(k)
766 
767  p(i,j,k)=p(i,j,k)+phi(i,j,k) - &
768  rk_c*0.5d0*dt*nu*(phixx+phiyy+phizz)
769 
770  enddo
771  enddo
772  enddo
773 
774  elseif ( cds .eq. 2) then
775 
776  do k=kstr3,kend3,1
777  do j=jstr3,jend3,1
778  do i=istr3,iend3,1
779 
780  phixx= ( -phi(i+2,j,k) + &
781  16.d0*phi(i+1,j,k) - &
782  30.d0*phi(i,j,k) + &
783  16.d0*phi(i-1,j,k) - &
784  phi(i-2,j,k) ) / 12.d0 / dx2
785 
786  phiyy= ( -phi(i,j+2,k) + &
787  16.d0*phi(i,j+1,k) - &
788  30.d0*phi(i,j,k) + &
789  16.d0*phi(i,j-1,k) - &
790  phi(i,j-2,k) ) / 12.d0 / dy2
791 
792  phizz= ( (phi(i,j,k+1)-phi(i,j,k))/dz_t(k) - &
793  (phi(i,j,k)-phi(i,j,k-1))/dz_b(k) ) / dz(k)
794 
795  p(i,j,k)=p(i,j,k)+phi(i,j,k) - &
796  rk_c*0.5d0*dt*nu*(phixx+phiyy+phizz)
797 
798  enddo
799  enddo
800  enddo
801 
802  endif
803 
804 
805 
806  return
807 
808  end subroutine
809 
810 
811 
812 end module
subroutine output_time_series(u, v, w, p, t)
Definition: module_io.f90:212
real(mytype), dimension(:,:,:), allocatable, save mean_uw
subroutine update_boundary_p(p)
real(mytype), dimension(:), allocatable, save aa_w
integer, save, protected is_ri_var
integer, save, protected imeanfl
subroutine initiate_fields
integer, save, protected myid
real(mytype), dimension(:), allocatable, save aa_t
complex(mytype), dimension(:,:,:), allocatable, save c_qu
integer, save, protected isxy2d
integer, save, protected ibackup
real(mytype), dimension(:,:,:), allocatable, save vhy
complex(mytype), dimension(:,:,:), allocatable, save bb_u
integer, save, protected isyz2d
real(mytype), dimension(:,:,:), allocatable, save thx
subroutine get_u_eqn_coeff(aa, bb, cc)
real(mytype), dimension(:,:,:), allocatable, save ss
complex(mytype), dimension(:,:,:), allocatable, save c_u
complex(mytype), dimension(:,:,:), allocatable, save c_qv
subroutine real_2_cmplx(r_in, c_out)
real(mytype), dimension(:,:,:), allocatable, save uhz
subroutine get_p_eqn_src(u, v, w, qp)
subroutine screen_div_error(u, v, w, isfinal, ishistout)
Definition: module_io.f90:86
real(mytype), dimension(:), allocatable, save aa_p
real(mytype), save, protected dt
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
subroutine get_temperature_cnvdiff_spec(u, v, w, t, c_t, c_ss, c_qt)
real(mytype), dimension(:,:,:), allocatable, save ff0
integer, save, protected jend3
real(mytype), dimension(:,:,:), allocatable, save mean_u
real(mytype), save, protected nu
real(mytype), dimension(:), allocatable, save cc_p
integer, save, protected intv_2d
complex(mytype), dimension(:,:,:), allocatable, save c_p
subroutine get_t_eqn_coeff(aa, bb, cc)
real(mytype), dimension(:), allocatable, save cc_w
real(mytype), dimension(:,:,:), allocatable, save ss0
integer, save, protected p_row
real(mytype), dimension(:,:,:), allocatable, save mean_v
real(mytype), dimension(:,:,:), allocatable, save cc
real(mytype), dimension(:,:,:), allocatable, save p
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)
real(mytype), dimension(:,:,:), allocatable, save thz
subroutine calculate_ubar(u)
subroutine output_2d_slices(u, v, w, p, t, time_step)
Definition: module_io.f90:409
real(mytype), dimension(:,:,:), allocatable, save mean_vw
integer, save, protected iinstfl
type(c_ptr), save, protected fft_plan1
subroutine spectral_truncation(c_inout3)
real(mytype), dimension(:,:,:), allocatable, save aa
real(mytype), dimension(:), allocatable, save, protected dz
complex(mytype), dimension(:,:,:), allocatable, save c_phi
subroutine poisson_solver_fft(aa, bb, cc, r_src, r_var, zstagger, delnull, tbnd, ubnd)
subroutine update_boundary_uvwhxyz(u, v, uhx, uhy, uhz, vhx, vhy, vhz, whx, why)
complex(mytype), dimension(:,:,:), allocatable, save bb_t
integer, save, protected myid_rowindx
real(mytype), dimension(:), allocatable, save cc_u
real(mytype), save, protected dy2
integer, save, protected isscalar
subroutine update_boundary_thxyz(thx, thy, thz)
real(mytype), dimension(:,:,:), allocatable, save mean_tt
complex(mytype), dimension(:,:,:), allocatable, save bb_w
subroutine assign_fgh_old(ff, gg, hh, ff0, gg0, hh0)
complex(mytype), dimension(:,:,:), allocatable, save c_gg
real(mytype), dimension(:,:,:), allocatable, save gg
subroutine rayleigh_damping(varin)
subroutine get_p_eqn_coeff(aa, bb, cc)
real(mytype), dimension(:,:,:), allocatable, save vhx
integer, save, protected istr3
real(mytype), dimension(:,:,:), allocatable, save hh0
subroutine assign_s_old(ss, ss0)
integer, save, protected isxz2d
real(mytype), dimension(:,:,:), allocatable, save mean_ww
real(mytype), save, protected dx
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 mean_t
real(mytype), dimension(:,:,:), allocatable, save w
integer, save, protected istmsr
real(mytype), dimension(:,:,:), allocatable, save u
real(mytype), dimension(:), allocatable, save, protected dz_b
type(c_ptr), save, protected ifft_plan1
subroutine update_ri(stps)
real(mytype), save, protected dx2
subroutine calculate_new_p(phi, p)
type(c_ptr), save, protected ifft_plan2
subroutine get_temperature_cnvdiff(u, v, w, t, thx, thy, thz, ss, ss0, qt)
subroutine cmplx_2_real(c_in, r_out)
subroutine output_backup(u, v, w, p, t, ff, gg, hh, ss)
Definition: module_io.f90:312
integer, save, protected p_col
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 isdamp
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(3), parameter rkc1
real(mytype), dimension(:,:,:), allocatable, save t
integer, save, protected imax
subroutine update_boundary_uvw(u, v, w)
subroutine calculate_uvw(u, v, w, phi)
real(mytype), dimension(:), allocatable, save, protected dz_t
real(mytype), dimension(:,:,:), allocatable, save qw
real(mytype), dimension(:,:,:), allocatable, save whx
real(mytype), dimension(:,:,:), allocatable, save uhx
real(mytype), save, protected dy
subroutine poisson_solver_fft_spec(aa, bb, cc, c_src_in, c_var_out, zstagger, delnull, tbnd, ubnd)
complex(mytype), dimension(:,:,:), allocatable, save c_t
real(mytype), dimension(:,:,:), allocatable, save why
subroutine update_boundary_t(t)
real(mytype), dimension(:,:,:), allocatable, save mean_tw
subroutine assign_abc(aa_in, bb_in, cc_in, aa, bb, cc)
real(mytype), dimension(:), allocatable, save cc_t
subroutine get_w_eqn_coeff(aa, bb, cc)
subroutine screen_div_error_spec(u, v, w, c_u, c_v, isfinal, ishistout)
real(mytype), dimension(:,:,:), allocatable, save ff
subroutine assign_abc_p(aa_in, bb_in, cc_in, aa, bb, cc)
real(mytype), dimension(:,:,:), allocatable, save mean_uu
real(mytype), dimension(:,:,:), allocatable, save v
subroutine output_inst_fields(u, v, w, p, t, time_step)
Definition: module_io.f90:367
subroutine get_div_spec(c_u, c_v, w, c_div)
real(mytype), dimension(3), parameter rkc2
type(c_ptr), save, protected fft_plan2
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), dimension(:,:,:), allocatable, save gg0
real(mytype), save, protected rk_c
real(mytype), dimension(:,:,:), allocatable, save qu
integer, save, protected iend3
real(mytype), dimension(:,:,:), allocatable, save phi
subroutine get_interp_fields_uvwhxyz(u, v, w, uhx, uhy, uhz, vhx, vhy, vhz, whx, why)
subroutine get_momentum_cnvdiff(u, v, w, p, t, uhx, uhy, uhz, vhx, vhy, vhz, whx, why, ff, ff0, gg, gg0, hh, hh0, qu, qv, qw)
subroutine output_mean_fields(mean_u, mean_v, mean_w, mean_t, mean_uu, mean_vv, mean_ww, mean_uw, mean_vw, mean_tt, mean_tw, time_step)
Definition: module_io.f90:548
complex(mytype), dimension(:,:,:), allocatable, save bb
real(mytype), dimension(3), parameter rkc3
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
subroutine screen_cpu_time(u, v, w)
Definition: module_io.f90:33
real(mytype), dimension(:,:,:), allocatable, save qv
subroutine assign_rk_coeff(a_in, b_in, c_in)
integer, save, protected dp_opt
real(mytype), dimension(:,:,:), allocatable, save uhy
complex(mytype), dimension(:,:,:), allocatable, save bb_p
integer, save, protected kstr3
integer, save, protected myid_colindx
real(mytype), dimension(:,:,:), allocatable, save vhz
integer, save, protected kend3