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 *,'nprc error'
339  print *,'please set nprc=p_row*p_col!'
340  stop
341  endif
342  endif
343 
344  return
345 
346  end subroutine
347 
348 
349 
350  ! this function finalizes mpi
351  subroutine finalize_mpi
353  use mpi
354 
355  implicit none
356 
357  call mpi_finalize( ierr )
358 
359  return
360 
361  end subroutine
362 
363 
364 
365  ! this function initiates i, j, k, indices, e.g., istr3, iend3, i/j_offset, etc.
366  subroutine initiate_indices
368  implicit none
369 
370  integer :: i,j,indx
371 
372  ! initiate indices for z pencil
373  istr3=1+ghst
374  iend3=zsize(1)+ghst
375  jstr3=1+ghst
376  jend3=zsize(2)+ghst
377  kstr3=1+ghst
378  kend3=zsize(3)+ghst
379 
380  ! initiate array size for x pencil
381  nia1=xsize(1)+2*ghst
382  nja1=xsize(2)+2*ghst
383  nka1=xsize(3)+2*ghst
384 
385  ! initiate array size for y pencil
386  nia2=ysize(1)+2*ghst
387  nja2=ysize(2)+2*ghst
388  nka2=ysize(3)+2*ghst
389 
390  ! initiate array size for z pencil
391  nia3=zsize(1)+2*ghst
392  nja3=zsize(2)+2*ghst
393  nka3=zsize(3)+2*ghst
394 
395  !calculate i and j indices offset based on myid
396  i_offset=int(myid/p_col)*nx/p_row
397  j_offset=mod(myid,p_col)*ny/p_col
398 
399  ! calculate the correspondence of processor id and row-column indices
400  ! For example, if p_row=4, p_col=3, and myid=5, then
401  ! myid_rowindx=2,myid_colindx=3, that being said, this processor
402  ! is assigned to the decomposed domain at 2nd row and 3rd column.
403  indx=0
404  do i=1,p_row
405  do j=1,p_col
406 
407  pij(i,j)=indx
408 
409  if ( myid .eq. pij(i,j) ) then
410 
411  myid_rowindx=i
412  myid_colindx=j
413 
414  endif
415 
416  indx=indx+1
417 
418  enddo
419  enddo
420 
421  return
422  end subroutine
423 
424 
425 
426  ! this function calculates mesh parameters, e.g., dz,l_t,dz_t,etc.
427  subroutine get_mesh_param
429  implicit none
430 
431  integer :: k
432  real(mytype) :: zmin,zmax,zmean,stretch_sysm,coeff_a
433  real(mytype),dimension(nka3) :: zz
434 
435  ! calculate z
436  zmin=0.d0
437  zmax=lz
438 
439  ! is it an open-channel or closed-channel?
440  if (ochannel .eq. 0) then
441 
442  ! closed channel
443  zmean=0.5d0*(zmin+zmax)
444  stretch_sysm=0.5d0
445  coeff_a=(2.d0+nz)*stretch_sysm ! this is a parameter
446 
447  elseif (ochannel .eq. 1) then
448 
449  ! open channel
450  zmean=(zmin+zmax)
451  stretch_sysm=1.d0
452  coeff_a=(1.d0+nz)*stretch_sysm ! this is a parameter
453 
454  endif
455 
456  ! calculate z
457  zz(:)=0.d0
458  do k=1,nz+1,1
459  zz(k+ghst)=zmean+(zmin-zmean) / tanh(2.d0*zstretch) * &
460  tanh(2.d0*zstretch*(k-coeff_a)/(1.d0-coeff_a))
461  enddo
462 
463  ! calculate dz
464  dz(:)=0.d0
465  do k=kstr3,kend3,1
466  dz(k)=zz(k+1)-zz(k)
467  enddo
468  dz(kstr3-1)=dz(kstr3)
469  dz(kend3+1)=dz(kend3)
470  dz(kstr3-ghst)=dz(kstr3)
471  dz(kend3+ghst)=dz(kend3)
472 
473  ! calculate interpolation coefficients l_t,l_b
474  l_t(:)=0.d0
475  l_b(:)=0.d0
476  do k=kstr3-1,kend3+1,1 ! we can extend lt and lb
477  l_t(k)=dz(k) / (dz(k)+dz(k+1))
478  l_b(k)=dz(k) / (dz(k)+dz(k-1))
479  enddo
480 
481  ! calculate dz_t and dz_b
482  dz_t(:)=0.d0
483  dz_b(:)=0.d0
484  do k=kstr3-1,kend3+1,1 ! we can extend dz_t and dz_b
485  dz_t(k)=0.5d0 * ( dz(k)+dz(k+1) )
486  dz_b(k)=0.5d0 * ( dz(k)+dz(k-1) )
487  enddo
488 
489  ! output z to files
490  if (myid .eq. 0) then
491 
492  write(*,*) 'Writing Geo Info to Files.'
493  open(11,file='results/geo.info')
494  do k=kstr3,kend3,1
495  write(11,'(4e14.6)') (zz(k+1)+zz(k))/2.d0,dz(k),l_t(k),l_b(k)
496  enddo
497  close(11)
498 
499  endif
500 
501  return
502 
503  end subroutine
504 
505 
506 
507  ! calculate bulk velocity
508  subroutine calculate_ubar(u)
510  use mpi
511 
512  implicit none
513 
514  real(mytype),dimension(:,:,:),intent(in) :: u
515  real(mytype) :: ubar_local
516  integer :: i,j,k
517 
518  ! here we calculate ubar
519  ubar_local=0.d0
520  do k=kstr3,kend3,1
521  do j=jstr3,jend3,1
522  do i=istr3,iend3,1
523 
524  ubar_local=ubar_local+u(i,j,k)*dz(k)
525 
526  enddo
527  enddo
528  enddo
529 
530  ubar_local=ubar_local/lz/zsize(1)/zsize(2)
531 
532  ! exchange ubar_local between processors
533  ! sum up all the ubar_local between all processors
534  ! and get ubar, this will be used in constant mass flow rate control
535  call mpi_allreduce(ubar_local,ubar,1,real_type,mpi_sum,mpi_comm_world,ierr)
536  ubar=ubar/nprc ! we need to do averaging
537 
538 
539  return
540 
541  end subroutine
542 
543 
544 
545  ! calculate fft plan for fftw
546  subroutine get_fft_plan
548  ! for using fftw
549  use, intrinsic :: iso_c_binding
550 
551  implicit none
552 
553  ! for using fftw
554  include 'include/fftw3.f03'
555 
556  complex(mytype),dimension(nx) :: c_tmp1
557  complex(mytype),dimension(ny) :: c_tmp2
558 
559  ! create plan (fft) for x direction
560  fft_plan1 =fftw_plan_dft_1d(nx,c_tmp1,c_tmp1,fftw_forward,fftw_measure)
561  ! create plan (ifft) for x direction
562  ifft_plan1=fftw_plan_dft_1d(nx,c_tmp1,c_tmp1,fftw_backward,fftw_measure)
563 
564  ! create plan (fft) for y direction
565  fft_plan2 =fftw_plan_dft_1d(ny,c_tmp2,c_tmp2,fftw_forward,fftw_measure)
566  ! create plan (ifft) for y direction
567  ifft_plan2=fftw_plan_dft_1d(ny,c_tmp2,c_tmp2,fftw_backward,fftw_measure)
568 
569  return
570 
571  end subroutine
572 
573 
574  ! this function updates ri if is_ri_var=1
575  subroutine update_ri(stps)
577  implicit none
578 
579  integer,intent(in) :: stps
580 
581  got=ri_str+(ri_end-ri_str)*1.d0/imax*stps
582 
583  return
584 
585  end subroutine
586 
587 
588  ! this function updates dpx_drive
589  subroutine update_dpx_drive
591  implicit none
592 
593  ! check which option to drive the flow
594  if (dp_opt .eq. 1) then
595 
596  ! constant friction
597  dpx_drive=1.d0
598 
599  elseif (dp_opt .eq. 2) then
600 
601  ! constant mass flow rate
603 
604  endif
605 
606  return
607 
608  end subroutine
609 
610 
611 
612  ! assign coefficient for the RK3 or AB2 scheme
613  subroutine assign_rk_coeff(a_in,b_in,c_in)
615  implicit none
616 
617  real(mytype), intent(in) :: a_in,b_in,c_in
618 
619  if (dts .eq. 2) then
620 
621  rk_a=a_in
622  rk_b=b_in
623  rk_c=c_in
624 
625  else
626 
627  ! AB2 scheme
628  rk_a= 1.5d0
629  rk_b=-0.5d0
630  rk_c= 1.d0
631 
632  endif
633 
634  return
635 
636  end subroutine
637 
638 
639 
640 end module
641 
642 
643 
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