module_parameters.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 
24  ! enable C standard for fftw
25  use, intrinsic :: iso_c_binding
26  ! 2d decomposition module
27  use decomp_2d
28 
29  implicit none
30 
31  ! define pi
32  real(mytype),parameter :: pi=3.1415926535897932_mytype
33 
34  !**************************************************************************
35  !********************* User Defined Parameters *********************
36  !**************************************************************************
37 
38  ! p_row/p_col: domain decomposition along the x/y directions
39  ! make sure to use p_row*p_col processors when running the solver
40  integer,protected,save :: p_row
41  integer,protected,save :: p_col
42 
43  ! Lx, Ly, Lz: domain sizes in x, y, and z directions.
44  ! x: streamwise, y: spanwise, z: vertical
45  real(mytype),protected,save :: lx
46  real(mytype),protected,save :: ly
47  real(mytype),protected,save :: lz
48  ! grids in x,y,z direction
49  integer,protected,save :: nx
50  integer,protected,save :: ny
51  integer,protected,save :: nz
52 
53  ! zstretch: this save control the slope in the boundary layer
54  ! the higher this value the finer grids near the walls
55  real(mytype),protected,save :: zstretch
56 
57  ! channel: 0: closed channel, 1: open channel
58  integer,protected,save :: ochannel
59 
60  ! dp_opt: options to drive the flow. dp_opt=3 is for Ekman layer
61  ! dp_opt=1: constant friction.
62  ! dp_opt=2: constant bulk velocity.
63  ! dp_opt=3: constant geostrophic wind.
64  integer,protected,save :: dp_opt
65 
66  ! issk=1: nonlinear term is casted in skew-symmetric form (only for cds=3)
67  integer,protected,save :: issk
68 
69  ! cds: differential scheme for spatial derivatives
70  ! 1: 2nd order CDS
71  ! 2: 4th order CDS
72  ! 3: Spectral method
73  integer,protected,save :: cds
74 
75  ! dts: time advancement scheme for the convective terms. 1: AB2; 2: RK3
76  ! for the viscous terms, Crank-Nicolson is used
77  integer,protected,save :: dts
78 
79  ! dt: time step; imax: how many steps to run
80  real(mytype),protected,save :: dt
81  integer,protected,save :: imax
82 
83  ! nu: kinematic viscosity; alpha: thermal diffusivity; got: g/T
84  ! fc: Coriolis parameter; ug, vg: geostrophic wind speeds
85  real(mytype),protected,save :: nu
86  real(mytype),protected,save :: alpha
87  real(mytype),protected,save :: got
88  real(mytype),protected,save :: fc
89  real(mytype),protected,save :: ug
90  real(mytype),protected,save :: vg
91 
92  ! isscalar: include temperature field? t_ref: a reference temperature
93  integer,protected,save :: isscalar
94  real(mytype),protected,save :: t_ref
95 
96  ! tbot, ttop: temperature boundary conditions at the bottom and top walls
97  real(mytype),protected,save :: tbot
98  real(mytype),protected,save :: ttop
99 
100  ! is_ri_var: varying Richardson number? ri_str, ri_end: ri at the start
101  ! and end of the simulation (only for is_ri_var=1)
102  integer,protected,save :: is_ri_var
103  real(mytype),protected,save :: ri_str
104  real(mytype),protected,save :: ri_end
105 
106  ! isnoise: add noise to the initial fields? if yes, what is the mag? noise_mag
107  integer,protected,save :: isnoise
108  real(mytype),protected,save :: noise_mag
109 
110  ! isdamp: Rayleigh damping? nzdamp: how many points?
111  ! cdamp: damping coefficient
112  integer,protected,save :: isdamp
113  integer,protected,save :: nzdamp
114  real(mytype),protected,save :: cdamp
115 
116  ! u_mrf: moving speed of the reference frame, set it to 0 for stationary frame
117  real(mytype),protected,save :: u_mrf
118 
119  ! ibackup/iinstfl/imeanfl:
120  ! how many steps to save backup/instantaneous/mean fields
121  integer,protected,save :: ibackup
122  integer,protected,save :: iinstfl
123  integer,protected,save :: imeanfl
124 
125  ! restart: read initial fields from files? istmsr: output time-series?
126  integer,protected,save :: restart
127  integer,protected,save :: istmsr
128 
129  ! isxy2d: whether to output 2d x-y slice data.
130  ! xy2d_id: where (k indices) to output
131  integer,protected,save :: isxy2d
132  integer,dimension(50),protected,save :: xy2d_id
133  ! isxz2d: whether to output 2d x-z slice data.
134  ! xz2d_id: where (j indices) to output
135  integer,protected,save :: isxz2d
136  integer,dimension(50),protected,save :: xz2d_id
137  ! isyz2d: whether to output 2d y-z slice data.
138  ! yz2d_id: where (i indices) to output
139  integer,protected,save :: isyz2d
140  integer,dimension(50),protected,save :: yz2d_id
141  ! intv_2d: output frequency (time step)
142  integer,protected,save :: intv_2d
143 
144  ! these are namelist parameters read from "parameters.input"
145  namelist /domain/ &
147 
148  namelist /modeling/ &
150 
151  namelist /constants/ &
152  nu, alpha, got, fc, ug, vg, dt, imax, u_mrf, isscalar, t_ref, &
154 
155  namelist /io/ &
158 
159 
160 
161  !**************************************************************************
162  !************************ Derived Parameters **************************
163  !**************************************************************************
164 
165  ! number of processors, nprc=p_row*p_col
166  integer,protected,save :: nprc
167 
168  ! iii: sqrt(-1)
169  complex(mytype),parameter :: iii=cmplx(0.d0,1.d0,kind=mytype)
170 
171  ! dx, dy: grid spacing in the x and y directions, dx2=dx*dx
172  real(mytype),protected,save :: dx, dy
173  real(mytype),protected,save :: dx2,dy2
174 
175  ! ghst: ghost cell number
176  integer,parameter :: ghst=2
177 
178  ! inner domain start and end indices on z pencil
179  integer,protected,save :: istr3,iend3,jstr3,jend3,kstr3,kend3
180 
181  ! nia, nja, nka: the dimensions of 3D array in i, j, and k indices
182  ! 1, 2, and 3 means x, y, and z pencil
183  integer,protected,save :: nia1, nja1, nka1
184  integer,protected,save :: nia2, nja2, nka2
185  integer,protected,save :: nia3, nja3, nka3
186 
187  ! i and j offset indices, used in mpi
188  integer,protected,save :: i_offset,j_offset
189 
190  ! div_hist_max: historically max divergence errors
191  ! time_sim: simulation time
192  ! cfl_max: maximal CFL number in the domain
193  ! these parameters are non-protected
194  real(mytype),save :: div_hist_max=0.d0
195  real(mytype),save :: time_sim=0.d0
196  real(mytype),save :: cfl_max=0.d0
197 
198  ! dpx_drive: mean pressure difference to drive the flow
199  ! if dp_opt=1, dpx_drive=1; if dp_opt=2 or 3, dpx_drive varies to
200  ! ensure constant flow rate or geostrophic wind.
201  ! ubar: mean velocity in the channel
202  ! dpx_coef: coefficient used for constant mass flow rate case
203  real(mytype),protected,save :: dpx_drive, dpx_coef, ubar
204 
205  ! wx1=e^(i*2pi/nx), wy1=e^(i*2pi/ny), used in fft Poisson solver
206  complex(mytype),protected,save :: wx1
207  complex(mytype),protected,save :: wy1
208 
209  ! coefficients for RK3 scheme, see Orlandi 2000
210  real(mytype),dimension(3),parameter ::rkc1=(/8.d0/15.d0,5.d0/12.d0,3.d0/4.d0/)
211  real(mytype),dimension(3),parameter ::rkc2=(/0.d0,-17.d0/60.d0,-5.d0/12.d0/)
212  real(mytype),dimension(3),parameter ::rkc3=(/8.d0/15.d0,2.d0/15.d0,1.d0/3.d0/)
213  real(mytype),protected,save ::rk_a,rk_b,rk_c
214 
215  ! dz: cell spacing in z direction. these are derived parameters
216  ! dz_t, dz_b: top and bottom spacing between cell vertex
217  ! l_t, l_b : linear interpolation coefficient for top and bottom
218  real(mytype),dimension(:),allocatable,protected,save :: &
219  dz,l_t,l_b,dz_t,dz_b
220 
221  ! myid: the id of the current processor, ierr: error flag
222  ! myid_row/colindx: i/j index of pij for the current processor
223  integer,protected,save :: myid,myid_rowindx,myid_colindx
224  integer,save :: ierr
225 
226  ! pij: processor id for given row and col indices
227  integer,dimension(:,:),allocatable,protected,save :: pij
228 
229  ! fft and ifft plans for fftw. 1 and 2 means x and y pencil
230  type(c_ptr),protected,save :: fft_plan1,ifft_plan1
231  type(c_ptr),protected,save :: fft_plan2,ifft_plan2
232 
233 
234 
235 
236 contains
237 
238 
239 
240  ! this function reads input parameters for the simulation,
241  ! allocates and initiates parameters, calculates derived parameters, etc
242  subroutine initiate_parameters
244  implicit none
245 
246  ! first read parameters from the file "parameters.input"
247  open(10,file='parameters.input')
248  read(10,nml=domain)
249  read(10,nml=modeling)
250  read(10,nml=constants)
251  read(10,nml=io)
252  close(10)
253 
254  ! test if nx and ny are divisible by p_row, and if ny and nz are
255  ! divisible by p_col
256  if (mod(nx,p_row) .ne. 0 .or. mod(ny,p_row) .ne. 0) then
257  print *,'nx and ny are not divisible by p_row! Quit!'
258  stop
259  endif
260 
261  if (mod(ny,p_col) .ne. 0 .or. mod(nz,p_col) .ne. 0) then
262  print *,'ny and nz are not divisible by p_col! Quit!'
263  stop
264  endif
265 
266  ! nx and ny must be even numbers
267  if (mod(nx,2) .ne. 0 .or. mod(ny,2) .ne. 0) then
268  print *,'Please use even numbers for nx and ny! Quit!'
269  stop
270  endif
271 
272  ! calculate some derived parameters
273  ! dx, dy: grid spacing in the x and y directions, dx2=dx*dx
274  dx=lx/nx
275  dy=ly/ny
276  dx2=dx*dx
277  dy2=dy*dy
278 
279  ! dpx_coef: coefficient used for constant mass flow rate case
280  dpx_coef=nu*1.d0
281 
282  ! wx1=e^(i*2pi/nx), wy1=e^(i*2pi/ny), used in fft Poisson solver
283  wx1=exp(iii*2.d0*pi/nx)
284  wy1=exp(iii*2.d0*pi/ny)
285 
286  ! allocate and initiate some arrays
287  ! allocate pij
288  allocate( pij(p_row,p_col) )
289  ! allocate grid parameters
290  allocate( dz(nz+2*ghst) )
291  allocate( l_t(nz+2*ghst) )
292  allocate( l_b(nz+2*ghst) )
293  allocate( dz_t(nz+2*ghst) )
294  allocate( dz_b(nz+2*ghst) )
295 
296  ! initiate them
297  ! grid parameters
298  dz(:)=0.d0
299  l_t(:)=0.d0
300  l_b(:)=0.d0
301  dz_t(:)=0.d0
302  dz_b(:)=0.d0
303  ! pij
304  pij(:,:)=0
305 
306  ! initiate the mean pressure difference for driving the flow,
307  ! if constant friction is used, dpx_drive=1
308  ! if constant mass flow rate is used, dpx_drive varies during the simulation
309  if (dp_opt .eq. 1) then
310  dpx_drive = 1.d0
311  else
312  dpx_drive = 0.d0 ! give inital values for dpx_drive and ubar
313  ubar = 0.d0
314  endif
315 
316  return
317 
318  end subroutine
319 
320 
321 
322  ! this function initiates mpi and tests if nprc=p_row*p_col
323  subroutine initiate_mpi
325  use mpi
326 
327  implicit none
328 
329  ! initialize mpi
330  call mpi_init( ierr )
331  call mpi_comm_rank( mpi_comm_world, myid, ierr )
332  call mpi_comm_size( mpi_comm_world, nprc,ierr )
333 
334  ! before we go, do some tests
335  ! test if nprc=p_row*p_col
336  if (nprc .ne. p_row*p_col) then
337  if (myid .eq. 0) then
338  print *,' '
339  print *,'Error when running mpirun -np nprc hercules.exe'
340  print *,'nprc=',nprc,', while p_row=',p_row,', p_col=',p_col
341  print *,'Please either run mpirun -np ',p_row*p_col,' hercules.exe'
342  print *,'or change p_row and p_col in parameters.info such that'
343  print *,'nprc=p_row*p_col'
344  print *,' '
345  stop
346  endif
347  endif
348 
349  return
350 
351  end subroutine
352 
353 
354 
355  ! this function finalizes mpi
356  subroutine finalize_mpi
358  use mpi
359 
360  implicit none
361 
362  call mpi_finalize( ierr )
363 
364  return
365 
366  end subroutine
367 
368 
369 
370  ! this function initiates i, j, k, indices, e.g., istr3, iend3, i/j_offset, etc.
371  subroutine initiate_indices
373  implicit none
374 
375  integer :: i,j,indx
376 
377  ! initiate indices for z pencil
378  istr3=1+ghst
379  iend3=zsize(1)+ghst
380  jstr3=1+ghst
381  jend3=zsize(2)+ghst
382  kstr3=1+ghst
383  kend3=zsize(3)+ghst
384 
385  ! initiate array size for x pencil
386  nia1=xsize(1)+2*ghst
387  nja1=xsize(2)+2*ghst
388  nka1=xsize(3)+2*ghst
389 
390  ! initiate array size for y pencil
391  nia2=ysize(1)+2*ghst
392  nja2=ysize(2)+2*ghst
393  nka2=ysize(3)+2*ghst
394 
395  ! initiate array size for z pencil
396  nia3=zsize(1)+2*ghst
397  nja3=zsize(2)+2*ghst
398  nka3=zsize(3)+2*ghst
399 
400  !calculate i and j indices offset based on myid
401  i_offset=int(myid/p_col)*nx/p_row
402  j_offset=mod(myid,p_col)*ny/p_col
403 
404  ! calculate the correspondence of processor id and row-column indices
405  ! For example, if p_row=4, p_col=3, and myid=5, then
406  ! myid_rowindx=2,myid_colindx=3, that being said, this processor
407  ! is assigned to the decomposed domain at 2nd row and 3rd column.
408  indx=0
409  do i=1,p_row
410  do j=1,p_col
411 
412  pij(i,j)=indx
413 
414  if ( myid .eq. pij(i,j) ) then
415 
416  myid_rowindx=i
417  myid_colindx=j
418 
419  endif
420 
421  indx=indx+1
422 
423  enddo
424  enddo
425 
426  return
427  end subroutine
428 
429 
430 
431  ! this function calculates mesh parameters, e.g., dz,l_t,dz_t,etc.
432  subroutine get_mesh_param
434  implicit none
435 
436  integer :: k
437  real(mytype) :: zmin,zmax,zmean,stretch_sysm,coeff_a
438  real(mytype),dimension(nka3) :: zz
439 
440  ! calculate z
441  zmin=0.d0
442  zmax=lz
443 
444  ! is it an open-channel or closed-channel?
445  if (ochannel .eq. 0) then
446 
447  ! closed channel
448  zmean=0.5d0*(zmin+zmax)
449  stretch_sysm=0.5d0
450  coeff_a=(2.d0+nz)*stretch_sysm ! this is a parameter
451 
452  elseif (ochannel .eq. 1) then
453 
454  ! open channel
455  zmean=(zmin+zmax)
456  stretch_sysm=1.d0
457  coeff_a=(1.d0+nz)*stretch_sysm ! this is a parameter
458 
459  endif
460 
461  ! calculate z
462  zz(:)=0.d0
463  do k=1,nz+1,1
464  zz(k+ghst)=zmean+(zmin-zmean) / tanh(2.d0*zstretch) * &
465  tanh(2.d0*zstretch*(k-coeff_a)/(1.d0-coeff_a))
466  enddo
467 
468  ! calculate dz
469  dz(:)=0.d0
470  do k=kstr3,kend3,1
471  dz(k)=zz(k+1)-zz(k)
472  enddo
473  dz(kstr3-1)=dz(kstr3)
474  dz(kend3+1)=dz(kend3)
475  dz(kstr3-ghst)=dz(kstr3)
476  dz(kend3+ghst)=dz(kend3)
477 
478  ! calculate interpolation coefficients l_t,l_b
479  l_t(:)=0.d0
480  l_b(:)=0.d0
481  do k=kstr3-1,kend3+1,1 ! we can extend lt and lb
482  l_t(k)=dz(k) / (dz(k)+dz(k+1))
483  l_b(k)=dz(k) / (dz(k)+dz(k-1))
484  enddo
485 
486  ! calculate dz_t and dz_b
487  dz_t(:)=0.d0
488  dz_b(:)=0.d0
489  do k=kstr3-1,kend3+1,1 ! we can extend dz_t and dz_b
490  dz_t(k)=0.5d0 * ( dz(k)+dz(k+1) )
491  dz_b(k)=0.5d0 * ( dz(k)+dz(k-1) )
492  enddo
493 
494  ! output z to files
495  if (myid .eq. 0) then
496 
497  write(*,*) 'Writing Geo Info to Files.'
498  open(11,file='results/geo.info')
499  do k=kstr3,kend3,1
500  write(11,'(4e14.6)') (zz(k+1)+zz(k))/2.d0,dz(k),l_t(k),l_b(k)
501  enddo
502  close(11)
503 
504  endif
505 
506  return
507 
508  end subroutine
509 
510 
511 
512  ! calculate bulk velocity
513  subroutine calculate_ubar(u)
515  use mpi
516 
517  implicit none
518 
519  real(mytype),dimension(:,:,:),intent(in) :: u
520  real(mytype) :: ubar_local
521  integer :: i,j,k
522 
523  ! here we calculate ubar
524  ubar_local=0.d0
525  do k=kstr3,kend3,1
526  do j=jstr3,jend3,1
527  do i=istr3,iend3,1
528 
529  ubar_local=ubar_local+u(i,j,k)*dz(k)
530 
531  enddo
532  enddo
533  enddo
534 
535  ubar_local=ubar_local/lz/zsize(1)/zsize(2)
536 
537  ! exchange ubar_local between processors
538  ! sum up all the ubar_local between all processors
539  ! and get ubar, this will be used in constant mass flow rate control
540  call mpi_allreduce(ubar_local,ubar,1,real_type,mpi_sum,mpi_comm_world,ierr)
541  ubar=ubar/nprc ! we need to do averaging
542 
543 
544  return
545 
546  end subroutine
547 
548 
549 
550  ! calculate fft plan for fftw
551  subroutine get_fft_plan
553  ! for using fftw
554  use, intrinsic :: iso_c_binding
555 
556  implicit none
557 
558  ! for using fftw
559  include 'include/fftw3.f03'
560 
561  complex(mytype),dimension(nx) :: c_tmp1
562  complex(mytype),dimension(ny) :: c_tmp2
563 
564  ! create plan (fft) for x direction
565  fft_plan1 =fftw_plan_dft_1d(nx,c_tmp1,c_tmp1,fftw_forward,fftw_measure)
566  ! create plan (ifft) for x direction
567  ifft_plan1=fftw_plan_dft_1d(nx,c_tmp1,c_tmp1,fftw_backward,fftw_measure)
568 
569  ! create plan (fft) for y direction
570  fft_plan2 =fftw_plan_dft_1d(ny,c_tmp2,c_tmp2,fftw_forward,fftw_measure)
571  ! create plan (ifft) for y direction
572  ifft_plan2=fftw_plan_dft_1d(ny,c_tmp2,c_tmp2,fftw_backward,fftw_measure)
573 
574  return
575 
576  end subroutine
577 
578 
579  ! this function updates ri if is_ri_var=1
580  subroutine update_ri(stps)
582  implicit none
583 
584  integer,intent(in) :: stps
585 
586  got=ri_str+(ri_end-ri_str)*1.d0/imax*stps
587 
588  return
589 
590  end subroutine
591 
592 
593  ! this function updates dpx_drive
594  subroutine update_dpx_drive
596  implicit none
597 
598  ! check which option to drive the flow
599  if (dp_opt .eq. 1) then
600 
601  ! constant friction
602  dpx_drive=1.d0
603 
604  elseif (dp_opt .eq. 2) then
605 
606  ! constant mass flow rate
608 
609  endif
610 
611  return
612 
613  end subroutine
614 
615 
616 
617  ! assign coefficient for the RK3 or AB2 scheme
618  subroutine assign_rk_coeff(a_in,b_in,c_in)
620  implicit none
621 
622  real(mytype), intent(in) :: a_in,b_in,c_in
623 
624  if (dts .eq. 2) then
625 
626  rk_a=a_in
627  rk_b=b_in
628  rk_c=c_in
629 
630  else
631 
632  ! AB2 scheme
633  rk_a= 1.5d0
634  rk_b=-0.5d0
635  rk_c= 1.d0
636 
637  endif
638 
639  return
640 
641  end subroutine
642 
643 
644 
645 end module
646 
647 
648 
real(mytype), save, protected ri_str
integer, save, protected is_ri_var
integer, save, protected imeanfl
integer, save, protected myid
real(mytype), save, protected rk_b
real(mytype), dimension(:), allocatable, save, protected l_b
integer, save, protected isxy2d
integer, save, protected ochannel
integer, save, protected nprc
real(mytype), save, protected ubar
integer, save, protected isnoise
integer, save, protected ny
real(mytype), save, protected dpx_coef
integer, save, protected ibackup
integer, save, protected isyz2d
integer, dimension(50), save, protected yz2d_id
real(mytype), save, protected fc
integer, save, protected issk
real(mytype), parameter pi
real(mytype), save, protected dt
integer, save, protected dts
complex(mytype), save, protected wy1
integer, dimension(:,:), allocatable, save, protected pij
real(mytype), save, protected rk_a
integer, save, protected nja2
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ttop
integer, save, protected nz
integer, save, protected intv_2d
real(mytype), save, protected ly
real(mytype), save, protected ug
integer, save, protected p_row
subroutine calculate_ubar(u)
integer, save, protected iinstfl
integer, save, protected nia3
type(c_ptr), save, protected fft_plan1
real(mytype), dimension(:), allocatable, save, protected dz
integer, save, protected myid_rowindx
integer, dimension(50), save, protected xy2d_id
real(mytype), save, protected dy2
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
integer, save, protected isxz2d
real(mytype), save, protected dx
integer, save, protected jstr3
integer, save, protected istmsr
real(mytype), dimension(:), allocatable, save, protected dz_b
type(c_ptr), save, protected ifft_plan1
integer, save, protected i_offset
subroutine update_ri(stps)
real(mytype), save, protected dx2
integer, save, protected nja3
type(c_ptr), save, protected ifft_plan2
integer, save, protected j_offset
real(mytype), save, protected u_mrf
complex(mytype), parameter iii
integer, save, protected p_col
integer, save, protected nka2
integer, parameter ghst
integer, save, protected cds
integer, save, protected isdamp
integer, save, protected restart
real(mytype), save, protected lz
real(mytype), dimension(3), parameter rkc1
real(mytype), save time_sim
integer, save, protected imax
real(mytype), save, protected lx
real(mytype), save, protected dpx_drive
integer, save, protected nzdamp
real(mytype), dimension(:), allocatable, save, protected dz_t
real(mytype), save, protected noise_mag
real(mytype), save, protected zstretch
real(mytype), save, protected dy
integer, save, protected nx
real(mytype), save, protected ri_end
integer, save, protected nka1
integer, save, protected nia2
integer, save, protected nia1
subroutine initiate_parameters
real(mytype), save, protected got
real(mytype), save div_hist_max
real(mytype), dimension(3), parameter rkc2
type(c_ptr), save, protected fft_plan2
complex(mytype), save, protected wx1
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), dimension(:), allocatable, save, protected l_t
real(mytype), dimension(3), parameter rkc3
real(mytype), save, protected tbot
integer, save, protected nja1
subroutine assign_rk_coeff(a_in, b_in, c_in)
integer, dimension(50), save, protected xz2d_id
integer, save, protected dp_opt
real(mytype), save cfl_max
integer, save, protected kstr3
integer, save, protected myid_colindx
integer, save, protected nka3
integer, save, protected kend3