!> This module contains the subroutines which set the initial value of the !> fields and the distribution function. module init_g use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN use constants, only: run_name_size implicit none private public :: ginit, ginitopt_restart_many, restart_file, get_restart_dir public :: init_init_g, finish_init_g, wnml_init_g, check_init_g, tstart public :: init_vnmult, new_field_init, initial_condition_is_nonadiabatic_dfn public :: init_g_config_type, set_init_g_config, get_init_g_config ! knobs integer :: ginitopt_switch integer, parameter :: ginitopt_default = 1, & ginitopt_xi = 3, ginitopt_xi2 = 4, ginitopt_rh = 5, ginitopt_zero = 6, & ginitopt_test3 = 7, ginitopt_convect = 8, & ginitopt_noise = 10, ginitopt_restart_many = 11, ginitopt_continue = 12, & ginitopt_nl = 13, ginitopt_kz0 = 14, ginitopt_restart_small = 15, & ginitopt_nl2 = 16, ginitopt_nl3 = 17, ginitopt_nl4 = 18, & ginitopt_nl5 = 19, ginitopt_alf = 20, ginitopt_kpar = 21, & ginitopt_nl6 = 22, ginitopt_nl7 = 23, ginitopt_gs = 24, ginitopt_recon = 25, & ginitopt_nl3r = 26, ginitopt_smallflat = 27, ginitopt_harris = 28, & ginitopt_recon3 = 29, ginitopt_ot = 30, & ginitopt_zonal_only = 31, ginitopt_single_parallel_mode = 32, & ginitopt_all_modes_equal = 33, & ginitopt_default_odd = 34, ginitopt_restart_eig=35, ginitopt_no_zonal = 36, & ginitopt_stationary_mode = 37, ginitopt_cgyro = 38, ginitopt_random_sine = 39, & ginitopt_constant = 40 real :: width0, dphiinit, phiinit, imfac, refac, zf_init, phifrac real :: init_spec_pow real :: den0, upar0, tpar0, tperp0 real :: den1, upar1, tpar1, tperp1 real :: den2, upar2, tpar2, tperp2 real :: tstart, scale, apar0 logical :: chop_side, clean_init, left, even, new_field_init logical :: initial_condition_is_nonadiabatic_dfn character (len = run_name_size) :: restart_file !WARNING: This is also defined in gs2_save. This should be addressed character (len=150) :: restart_dir integer, dimension(2) :: ikk, itt integer, dimension(3) :: ikkk,ittt complex, dimension(3) :: phiamp, aparamp integer :: restart_eig_id integer :: max_mode ! These are used for the function ginit_single_parallel_mode, and specify the ! kparallel to initialize. In the case of zero magnetic shear, of course, the box ! is periodic in the parallel direction, and so only harmonics of the box size ! (i.e. ikpar_init) are allowed EGH</doc> integer :: ikpar_init real :: kpar_init !> This is used in linear runs with flow shear in order to track the !! evolution of a single Lagrangian mode. integer :: ikx_init ! RN> for recon3 real :: phiinit0 ! amplitude of equilibrium real :: phiinit_rand ! amplitude of random perturbation real :: a0,b0 ! amplitude of equilibrium Apara or By ! if b0 /= 0, u0 is rescaled to give this By amplitude by overriding a0 ! if a0 /= 0, u0 is rescaled to give this Apara amplitude logical :: null_phi, null_bpar, null_apar ! nullify fields ! if use_{Phi,Bpar,Apar} = F, override corresponding flags integer :: adj_spec ! adjust input parameter of this spec ! if = 0, just warn if given conditions are not satisfied character (len=20) :: eq_type = '' ! equilibrium type = 'sinusoidal', 'porcelli', 'doubleharris' real :: prof_width=-.05 ! width of porcelli, doubleharris profile ! if negative, this gives the ratio to the box size integer :: eq_mode_n=0, eq_mode_u=1 ! mode number in x for sinusoidal equilibrium logical :: input_check_recon=.false. complex :: ukxy_pt(3) ! <RN !> Whether to use the constant_random number !! generator, which always gives the same !! random numbers. Used for testing. logical :: constant_random_flag logical, parameter :: debug = .false. logical :: initialized = .false. !> Used to represent the input configuration of init_g type, extends(abstract_config_type) :: init_g_config_type ! namelist : init_g_knobs ! indexed : false !> Amplitude of equilibrium `Apara`: if `a0` is non-zero, \(u_0\) is !> rescaled to give this `Apara` amplitude real :: a0 = 0.0 !> Used by [[init_g_knobs:ginit_option]] `"recon3"`: adjust input parameter !> of this species; if = 0, just warn if given conditions are not satisfied !> !> FIXME: What input parameter, what conditions? integer :: adj_spec = 0 !> Used by [[init_g_knobs:ginit_option]] `"nl3r"`: set amplitude of `apar` !> if [[init_g_knobs:width0]] is negative real :: apar0 = 0.0 !> Used in initialization for the Orszag-Tang 2D vortex problem !> ([[init_g_knobs:ginit_option]] `"ot"`). Controls size of `jpar` term. !> FIXME: Reference equations? complex, dimension(3) :: aparamp = cmplx(0.0,0.0) !> Amplitude of equilibrium `By`: if `b0` is non-zero, \(u_0\) is !> rescaled to give this `By` amplitude by overriding [[init_g_knobs:a0]] real :: b0 = 0.0 !> Rarely needed. Forces asymmetry into initial condition. Warning: This !> does not behave as one may expect in flux tube simulations (see !> clean_init), this can be important as the default is to use chop_side logical :: chop_side = .true. !> Used with [[init_g_knobs:ginit_option]] `"noise"`. Ensures that in flux !> tube simulations the connected boundary points are initialised to the !> same value. Also allows for [[init_g_knobs:chop_side]] to behave !> correctly in flux tube simulations. logical :: clean_init = .true. !> Use the constant_random number generator, which always gives the same !> random numbers. Useful for testing. logical :: constant_random_flag = .false. !> Amplitude of zeroth density perturbation for [[init_g_knobs:ginit_option]] !> options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"ot"` (see [[init_g:ginit_ot]]) !> - `"recon"` (see [[init_g:ginit_recon]]) real :: den0 = 1.0 !> Amplitude of first density perturbation for [[init_g_knobs:ginit_option]] !> options: !> !> - `"gs"` (see [[init_g:ginit_gs]]) !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: den1 = 0.0 !> Amplitude of second density perturbation for [[init_g_knobs:ginit_option]] !> options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: den2 = 0.0 !> Amplitude of \(\phi\) for `"nl7"` (see [[init_g:ginit_nl7]]) real :: dphiinit = 1.0 !> Mode number in \(x\) for the sinusoidal equilbrium option of !> [[init_g_knobs:ginit_option]] `"recon3"` (see [[init_g:ginit_recon3]]), !> applies to density and temperatures integer :: eq_mode_n = 0 !> Mode number in \(x\) for the sinusoidal equilbrium option of !> [[init_g_knobs:ginit_option]] `"recon3"` (see [[init_g:ginit_recon3]]), !> applies to parallel velocity integer :: eq_mode_u = 1 !> Equilbrium choice for [[init_g_knobs:ginit_option]] `"recon3"` (see !> [[init_g:ginit_recon3]]). Choices are: !> !> - `"sinusoidal"` !> - `"porcelli"` !> - `"doubleharris"` character(len = 20) :: eq_type = 'sinusoidal' !> Sometimes initial conditions have definite parity; this picks the parity !> in those cases. Affects the following choices of !> [[init_g_knobs:ginit_option]]: !> !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) logical :: even = .true. !> Method for initialising `g`. There are many possible options, but the !> typical ones are: !> !> - `"default"`: Gaussian in \(\theta\). See [[init_g_knobs:phiinit]], !> [[init_g_knobs:width0]], [[init_g_knobs:chop_side]] !> - `"default_odd"`: Guassian in \(\theta\) multiplied by \(\sin(\theta - !> \theta_0)\). See [[init_g_knobs:phiinit]], [[init_g_knobs:width0]], !> [[init_g_knobs:chop_side]] !> - `"noise"`: Noise along field line. The recommended (but not default) !> selection. See [[init_g_knobs:phiinit]], [[init_g_knobs:zf_init]] !> - `"restart", "many"`: Restart, using `g` from restart files. !> - `"single_parallel_mode"`: Initialise only with a single parallel mode !> specified by either [[init_g_knobs:ikpar_init]] for periodic !> boundary conditions or [[init_g_knobs:kpar_init]] for linked !> boundary conditions. Intended for linear calculations. !> - `"eig_restart"`: Uses the restart files written by the !> eigensolver. Also see [[init_g_knobs:restart_eig_id]]. !> - `"random_sine"`: This option is notable as it attempts to make sure !> that the initial condition satisfies the parallel boundary condition. !> !> All the options are detailed further in the [initialisation options !> documentation](../user_manual/initialisation.html) character(len = 20) :: ginit_option = "default" !> \(k_y\) indices of two modes to initialise for certain options integer, dimension(2) :: ikk = [1, 2] !> \(k_y\) indices of three modes to initialise for certain options integer, dimension(3) :: ikkk = [1, 2, 2] !> Parallel mode number to initialise for [[init_g_knobs:ginit_option]] !> `"single_parallel_mode"` with periodic boundary conditions integer :: ikpar_init = 0 !> If greater than zero, initialise only the given \(k_x\) with !> [[init_g_knobs:ginit_option]] `"noise"`. This is useful for linear runs !> with flow shear, to track the evolution of a single Lagrangian mode. integer :: ikx_init = -1 !> Amplitude of imaginary part of \(\phi\) for various !> [[init_g_knobs:ginit_option]] options real :: imfac = 0.0 !> Initial spectrum power index. with !> [[init_g_knobs:ginit_option]] `"noise"`. real :: init_spec_pow = 0.0 !> Do we want to include the explicit source terms and !> related timesteps when saving/restoring from the restart file. logical :: include_explicit_source_in_restart = .true. !> If true then the initial condition is used to set the non-adiabatic !> distribution function rather than the `g` evolved by GS2. logical :: initial_condition_is_nonadiabatic_dfn = .false. !> Print some debug information for [[init_g_knobs:ginit_option]] `"recon3"` logical :: input_check_recon = .false. !> \(k_x\) indices of two modes to initialise for certain options integer, dimension(2) :: itt = [1, 2] !> \(k_x\) indices of three modes to initialise for certain options integer, dimension(3) :: ittt = [1, 2, 2] !> Parallel mode number to initialise for [[init_g_knobs:ginit_option]] !> `"single_parallel_mode"` with linked boundary conditions real :: kpar_init = 0.0 !> If [[init_g_knobs:chop_side]] is true, chop out left side in theta if !> true, right side if false. Applies to: FIXME: list of initial conditions logical :: left = .true. !> The maximum mode number to initialise for [[init_g_knobs:ginit_option]] !> `"random_sine"`. Actual maximum used will be lower of max_mode and `ntheta/8` integer :: max_mode = 128 !> Change fields implicit initialisation somehow !> !> FIXME: Add documentation logical :: new_field_init = .true. !> "Nullify fields". Used by [[init_g_knobs:ginit_option]] `"recon3"` !> !> FIXME: Add documentation logical :: null_apar = .false. !> "Nullify fields". Used by [[init_g_knobs:ginit_option]] `"recon3"` !> !> FIXME: Add documentation logical :: null_bpar = .false. !> "Nullify fields". Used by [[init_g_knobs:ginit_option]] `"recon3"` !> !> FIXME: Add documentation logical :: null_phi = .false. !> Used in initialization for the Orszag-Tang 2D vortex problem. !> !> FIXME: expand complex, dimension(3) :: phiamp = cmplx(0.0,0.0) !> If doing multiple flux tube calculations, multiply !> [[init_g_knobs:phiinit]] by `(job * phifrac + 1)`. Allows for ensemble !> averaging of multiple flux tube calculations real :: phifrac = 0.1 !> Average amplitude of initial perturbation of each Fourier mode. Used by !> most [[init_g_knobs:ginit_option]] options real :: phiinit = 1.0 !> Amplitude of equilibrium. Used by [[init_g_knobs:ginit_option]] `"recon3"` real :: phiinit0 = 0.0 !> Amplitude of random perturbation. Used by [[init_g_knobs:ginit_option]] `"recon3"` real :: phiinit_rand = 0.0 !> Which processor is responsible for saving/reading potentials to/from restart !> files. If <0 then all procs write to restart files. If >= nproc then set to !> proc_to_save_fields = mod(proc_to_save_fields, nproc) integer :: proc_to_save_fields = 0 !> Width of "porcelli", "doubleharris" profiles in !> [[init_g_knobs:ginit_option]] `"recon3"`. If negative, this gives the !> ratio to the box size real :: prof_width = -0.1 !> Restart from many restart files. If false, use a single restart file: !> this requires GS2 to have been built with parallel netCDF. logical :: read_many = .false. !> Amplitude of real part of \(\phi\) for various !> [[init_g_knobs:ginit_option]] options real :: refac = 1.0 !> Directory to read/write restart files to/from. Make sure this exists !> before running (see [[gs2_diagnostics_knobs:file_safety_check]]) character(len = 150) :: restart_dir = "./" !> Used to select with eigensolver generated restart file to load. Sets !> `<id>` in `restart_file//eig_<id>//.<proc>` string used to set !> filename. Note that this is zero indexed. integer :: restart_eig_id = 0 !> Basename of restart files !> !> @note gets a smart default character(len = run_name_size) :: restart_file = "restart_file_not_set.nc" !> Rescale amplitudes of fields for restart [[init_g_knobs:ginit_option]] !> options such as `"restart"`. Note that if `scale` is negative, the !> fields are scaled by `-scale / max(abs(phi))`! real :: scale = 1.0 !> Amplitude of zeroth parallel temperature perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"ot"` (see [[init_g:ginit_ot]]) !> - `"recon"` (see [[init_g:ginit_recon]]) real :: tpar0 = 0.0 !> Amplitude of first parallel temperature perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"gs"` (see [[init_g:ginit_gs]]) !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: tpar1 = 0.0 !> Amplitude of second parallel temperature perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: tpar2 = 0.0 !> Amplitude of zeroth perpendicular temperature perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"ot"` (see [[init_g:ginit_ot]]) !> - `"recon"` (see [[init_g:ginit_recon]]) real :: tperp0 = 0.0 !> Amplitude of first perpendicular temperature perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"gs"` (see [[init_g:ginit_gs]]) !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: tperp1 = 0.0 !> Amplitude of second perpendicular temperature perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: tperp2 = 0.0 !> Force `t = tstart` at beginning of run real :: tstart = 0.0 !> Perturbation amplitude of parallel velocity? Used by !> [[init_g_knobs:ginit_option]] `"recon3"` !> !> FIXME: Add documentation complex, dimension(3) :: ukxy_pt = cmplx(0.,0.) !> Amplitude of zeroth parallel velocity perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"ot"` (see [[init_g:ginit_ot]]) !> - `"recon"` (see [[init_g:ginit_recon]]) real :: upar0 = 0.0 !> Amplitude of first parallel velocity perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"gs"` (see [[init_g:ginit_gs]]) !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: upar1 = 0.0 !> Amplitude of second parallel velocity perturbation for !> [[init_g_knobs:ginit_option]] options: !> !> - `"kpar"` (see [[init_g:ginit_kpar]]) !> - `"kz0"` (see [[init_g:ginit_kz0]]) !> - `"nl3"` (see [[init_g:ginit_nl3]]) !> - `"nl3r"` (see [[init_g:ginit_nl3r]]) !> - `"nl7"` (see [[init_g:ginit_nl7]]) !> - `"recon"` (see [[init_g:ginit_recon]]) !> - `"recon3"` (see [[init_g:ginit_recon3]]) real :: upar2 = 0.0 !> Width of initial perturbation Gaussian envelope in \(\theta\) real :: width0 = -3.5 !> Amplitude of initial zonal flow perturbations relative to other modes real :: zf_init = 1.0 contains procedure, public :: read => read_init_g_config procedure, public :: write => write_init_g_config procedure, public :: reset => reset_init_g_config procedure, public :: broadcast => broadcast_init_g_config procedure, public, nopass :: get_default_name => get_default_name_init_g_config procedure, public, nopass :: get_default_requires_index => get_default_requires_index_init_g_config procedure :: set_smart_defaults => set_smart_defaults_local end type init_g_config_type type(init_g_config_type) :: init_g_config contains !> FIXME : Add documentation subroutine wnml_init_g(unit) use run_parameters, only: k0 implicit none integer, intent(in) :: unit write (unit, *) write (unit, fmt="(' &',a)") "init_g_knobs" select case (ginitopt_switch) case (ginitopt_default) write (unit, fmt="(' ginit_option = ',a)") '"default"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' width0 = ',e17.10)") width0 write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left case (ginitopt_noise) write (unit, fmt="(' ginit_option = ',a)") '"noise"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' zf_init = ',e17.10)") zf_init write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left write (unit, fmt="(' clean_init = ',L1)") clean_init case (ginitopt_xi) write (unit, fmt="(' ginit_option = ',a)") '"xi"' write (unit, fmt="(' width0 = ',e17.10)") width0 case (ginitopt_xi2) write (unit, fmt="(' ginit_option = ',a)") '"xi2"' write (unit, fmt="(' width0 = ',e17.10)") width0 case (ginitopt_zero) write (unit, fmt="(' ginit_option = ',a)") '"zero"' case (ginitopt_test3) write (unit, fmt="(' ginit_option = ',a)") '"test3"' case (ginitopt_convect) write (unit, fmt="(' ginit_option = ',a)") '"convect"' write (unit, fmt="(' k0 = ',e17.10)") k0 case (ginitopt_rh) write (unit, fmt="(' ginit_option = ',a)") '"rh"' case (ginitopt_restart_many) write (unit, fmt="(' ginit_option = ',a)") '"many"' write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"' write (unit, fmt="(' scale = ',e17.10)") scale case (ginitopt_restart_small) write (unit, fmt="(' ginit_option = ',a)") '"small"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' zf_init = ',e17.10)") zf_init write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"' write (unit, fmt="(' scale = ',e17.10)") scale case (ginitopt_continue) write (unit, fmt="(' ginit_option = ',a)") '"cont"' case (ginitopt_kz0) write (unit, fmt="(' ginit_option = ',a)") '"kz0"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left case (ginitopt_nl) write (unit, fmt="(' ginit_option = ',a)") '"nl"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1) write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2) write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left case (ginitopt_nl2) write (unit, fmt="(' ginit_option = ',a)") '"nl2"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1) write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2) write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left case (ginitopt_nl3) write (unit, fmt="(' ginit_option = ',a)") '"nl3"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' width0 = ',e17.10)") width0 write (unit, fmt="(' refac = ',e17.10)") refac write (unit, fmt="(' imfac = ',e17.10)") imfac write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1) write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2) write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left write (unit, fmt="(' den0 = ',e17.10)") den0 write (unit, fmt="(' den1 = ',e17.10)") den1 write (unit, fmt="(' den2 = ',e17.10)") den2 write (unit, fmt="(' upar0 = ',e17.10)") upar0 write (unit, fmt="(' upar1 = ',e17.10)") upar1 write (unit, fmt="(' upar2 = ',e17.10)") upar2 write (unit, fmt="(' tpar0 = ',e17.10)") tpar0 write (unit, fmt="(' tpar1 = ',e17.10)") tpar1 write (unit, fmt="(' tperp0 = ',e17.10)") tperp0 write (unit, fmt="(' tperp1 = ',e17.10)") tperp1 write (unit, fmt="(' tperp2 = ',e17.10)") tperp2 case (ginitopt_nl4) write (unit, fmt="(' ginit_option = ',a)") '"nl4"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"' write (unit, fmt="(' scale = ',e17.10)") scale write (unit, fmt="(' ikk(1) = ',i3,' itt(1) = ',i3)") ikk(1),itt(1) write (unit, fmt="(' ikk(2) = ',i3,' itt(2) = ',i3)") ikk(2), itt(2) write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left case (ginitopt_nl5) write (unit, fmt="(' ginit_option = ',a)") '"nl5"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"' write (unit, fmt="(' scale = ',e17.10)") scale write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left case (ginitopt_nl6) write (unit, fmt="(' ginit_option = ',a)") '"nl6"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' restart_file = ',a)") '"'//trim(restart_file)//'"' write (unit, fmt="(' scale = ',e17.10)") scale case (ginitopt_alf) write (unit, fmt="(' ginit_option = ',a)") '"alf"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit case (ginitopt_gs) write (unit, fmt="(' ginit_option = ',a)") '"gs"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' refac = ',e17.10)") refac write (unit, fmt="(' imfac = ',e17.10)") imfac write (unit, fmt="(' den1 = ',e17.10)") den1 write (unit, fmt="(' upar1 = ',e17.10)") upar1 write (unit, fmt="(' tpar1 = ',e17.10)") tpar1 write (unit, fmt="(' tperp1 = ',e17.10)") tperp1 case (ginitopt_kpar) write (unit, fmt="(' ginit_option = ',a)") '"kpar"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' width0 = ',e17.10)") width0 write (unit, fmt="(' refac = ',e17.10)") refac write (unit, fmt="(' imfac = ',e17.10)") imfac write (unit, fmt="(' den0 = ',e17.10)") den0 write (unit, fmt="(' den1 = ',e17.10)") den1 write (unit, fmt="(' den2 = ',e17.10)") den2 write (unit, fmt="(' upar0 = ',e17.10)") upar0 write (unit, fmt="(' upar1 = ',e17.10)") upar1 write (unit, fmt="(' upar2 = ',e17.10)") upar2 write (unit, fmt="(' tpar0 = ',e17.10)") tpar0 write (unit, fmt="(' tpar1 = ',e17.10)") tpar1 write (unit, fmt="(' tperp0 = ',e17.10)") tperp0 write (unit, fmt="(' tperp1 = ',e17.10)") tperp1 write (unit, fmt="(' tperp2 = ',e17.10)") tperp2 case (ginitopt_cgyro) write (unit, fmt="(' ginit_option = ',a)") '"cgyro"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' zf_init = ',e17.10)") zf_init write (unit, fmt="(' chop_side = ',L1)") chop_side if (chop_side) write (unit, fmt="(' left = ',L1)") left write (unit, fmt="(' clean_init = ',L1)") clean_init case (ginitopt_random_sine) write (unit, fmt="(' ginit_option = ',a)") '"random_sine"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' zf_init = ',e17.10)") zf_init write (unit, fmt="(' max_mode = ',I0)") max_mode case (ginitopt_constant) write (unit, fmt="(' ginit_option = ',a)") '"constant"' write (unit, fmt="(' phiinit = ',e17.10)") phiinit write (unit, fmt="(' zf_init = ',e17.10)") zf_init end select write (unit, fmt="(' /')") end subroutine wnml_init_g !> FIXME : Add documentation subroutine check_init_g(report_unit) use run_parameters, only : delt_option_switch, delt_option_auto use species, only : spec, has_electron_species implicit none integer, intent(in) :: report_unit select case (ginitopt_switch) case (ginitopt_default) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="(' Width in theta: ',f10.4)") width0 if (chop_side) then write (report_unit, fmt="(' Parity: none')") else write (report_unit, fmt="(' Parity: even')") end if case (ginitopt_kz0) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="(' Constant along field line',f10.4)") width0 if (chop_side) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="(' Parity: none')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('Remedy: set chop_side = .false. in init_g_knobs.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if case (ginitopt_noise) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="(' Noise along field line.')") if (zf_init /= 1.) then write (report_unit, fmt="(' Zonal flows adjusted by factor of zf_init = ',f10.4)") zf_init end if case (ginitopt_kpar) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="(' Real part multiplier: ',f10.4)") refac write (report_unit, fmt="(' Imag part multiplier: ',f10.4)") imfac if (width0 > 0.) then write (report_unit, fmt="(' Gaussian envelope in theta with width: ',f10.4)") width0 end if if (chop_side) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="(' Parity: none')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('Remedy: set chop_side = .false. in init_g_knobs.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if if (den0 > epsilon(0.0) .or. den1 > epsilon(0.0) .or. den2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial density perturbation of the form:')") write (report_unit, fmt="('den0 + den1 * cos(theta) + den2 * cos(2.*theta)')") write (report_unit, *) write (report_unit, fmt="('with den0 =',f7.4,' den1 = ',f7.4,' den2 = ',f7.4)") den0, den1, den2 end if if (upar0 > epsilon(0.0) .or. upar1 > epsilon(0.0) .or. upar2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial parallel velocity perturbation of the form:')") write (report_unit, fmt="('upar0 + upar1 * sin(theta) + upar2 * sin(2.*theta)')") write (report_unit, fmt="('90 degrees out of phase with other perturbations.')") write (report_unit, *) write (report_unit, fmt="('with upar0 =',f7.4,' upar1 = ',f7.4,' upar2 = ',f7.4)") upar0, upar1, upar2 end if if (tpar0 > epsilon(0.0) .or. tpar1 > epsilon(0.0) .or. tpar2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial Tpar perturbation of the form:')") write (report_unit, fmt="('tpar0 + tpar1 * cos(theta) + tpar2 * cos(2.*theta)')") write (report_unit, *) write (report_unit, fmt="('with tpar0 =',f7.4,' tpar1 = ',f7.4,' tpar2 = ',f7.4)") tpar0, tpar1, tpar2 end if if (tperp0 > epsilon(0.0) .or. tperp1 > epsilon(0.0) .or. tperp2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial Tperp perturbation of the form:')") write (report_unit, fmt="('tperp0 + tperp1 * cos(theta) + tperp2 * cos(2.*theta)')") write (report_unit, *) write (report_unit, fmt="('with tperp0 =',f7.4,' tperp1 = ',f7.4,' tperp2 = ',f7.4)") tperp0, tperp1, tperp2 end if if (has_electron_species(spec)) then write (report_unit, *) write (report_unit, fmt="('Field line average of g_electron subtracted off.')") end if case (ginitopt_gs) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Randomly phased kpar=1 sines and cosines')") write (report_unit, fmt="(' in density, upar, tpar, or tperp.')") write (report_unit, fmt="(' Real part amplitude: ',f10.4)") refac*phiinit write (report_unit, fmt="(' Imag part amplitude: ',f10.4)") imfac*phiinit if (abs( den1) > epsilon(0.0)) write (report_unit, fmt="(' Density amplitude: ',f10.4)") den1 if (abs( upar1) > epsilon(0.0)) write (report_unit, fmt="(' Upar amplitude: ',f10.4)") upar1 if (abs( tpar1) > epsilon(0.0)) write (report_unit, fmt="(' Tpar amplitude: ',f10.4)") tpar1 if (abs(tperp1) > epsilon(0.0)) write (report_unit, fmt="(' Tperp amplitude: ',f10.4)") tperp1 case (ginitopt_nl) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="('At most two k_perps excited, with amplitude = ',f10.4)") phiinit write (report_unit, fmt="(' First k_perp has ik = ',i3,' it = ',i3)") ikk(1), itt(1) write (report_unit, fmt="('Second k_perp has ik = ',i3,' it = ',i3)") ikk(2), itt(2) if (chop_side) then write (report_unit, fmt="(' Parity: none')") else write (report_unit, fmt="(' Parity: even')") end if write (report_unit, fmt="('Reality condition is enforced.')") case (ginitopt_nl2) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="('At most two k_perps excited, with amplitude = ',f10.4)") phiinit write (report_unit, fmt="(' First k_perp has ik = ',i3,' it = ',i3)") ikk(1), itt(1) write (report_unit, fmt="('Second k_perp has ik = ',i3,' it = ',i3)") ikk(2), itt(2) if (chop_side) then write (report_unit, fmt="(' Parity: none')") else write (report_unit, fmt="(' Parity: even')") end if write (report_unit, fmt="('Reality condition is enforced.')") write (report_unit, fmt="('g perturbation proportional to (1+v_parallel)*sin(theta)')") case (ginitopt_nl3) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="('At most two k_perps excited, with amplitude = ',f10.4)") phiinit write (report_unit, fmt="(' First k_perp has ik = ',i3,' it = ',i3)") ikk(1), itt(1) write (report_unit, fmt="('Second k_perp has ik = ',i3,' it = ',i3)") ikk(2), itt(2) write (report_unit, fmt="(' Real part multiplied by: ',f10.4)") refac write (report_unit, fmt="(' Imag part multiplied by: ',f10.4)") imfac if (width0 > 0.) then write (report_unit, fmt="(' Gaussian envelope in theta with width: ',f10.4)") width0 end if if (chop_side) then write (report_unit, fmt="(' Parity: none')") else write (report_unit, fmt="(' Parity: even')") end if write (report_unit, fmt="('Reality condition is enforced.')") if (den0 > epsilon(0.0) .or. den1 > epsilon(0.0) .or. den2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial density perturbation of the form:')") write (report_unit, fmt="('den0 + den1 * cos(theta) + den2 * cos(2.*theta)')") write (report_unit, *) write (report_unit, fmt="('with den0 =',f7.4,' den1 = ',f7.4,' den2 = ',f7.4)") den0, den1, den2 end if if (upar0 > epsilon(0.0) .or. upar1 > epsilon(0.0) .or. upar2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial parallel velocity perturbation of the form:')") write (report_unit, fmt="('upar0 + upar1 * cos(theta) + upar2 * cos(2.*theta)')") write (report_unit, fmt="('90 degrees out of phase with other perturbations.')") write (report_unit, *) write (report_unit, fmt="('with upar0 =',f7.4,' upar1 = ',f7.4,' upar2 = ',f7.4)") upar0, upar1, upar2 end if if (tpar0 > epsilon(0.0) .or. tpar1 > epsilon(0.0) .or. tpar2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial Tpar perturbation of the form:')") write (report_unit, fmt="('tpar0 + tpar1 * cos(theta) + tpar2 * cos(2.*theta)')") write (report_unit, *) write (report_unit, fmt="('with tpar0 =',f7.4,' tpar1 = ',f7.4,' tpar2 = ',f7.4)") tpar0, tpar1, tpar2 end if if (tperp0 > epsilon(0.0) .or. tperp1 > epsilon(0.0) .or. tperp2 > epsilon(0.0)) then write (report_unit, *) write (report_unit, fmt="('Initial Tperp perturbation of the form:')") write (report_unit, fmt="('tperp0 + tperp1 * cos(theta) + tperp2 * cos(2.*theta)')") write (report_unit, *) write (report_unit, fmt="('with tperp0 =',f7.4,' tperp1 = ',f7.4,' tperp2 = ',f7.4)") tperp0, tperp1, tperp2 end if case (ginitopt_nl4) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Under development for study of secondary instabilities.')") write (report_unit, fmt="('Scale factor: ',f10.4)") scale write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_nl5) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Under development for study of secondary instabilities.')") write (report_unit, fmt="('Scale factor: ',f10.4)") scale write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_nl6) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Change amplitude of a particular mode.')") write (report_unit, fmt="('Scale factor: ',f10.4)") scale write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_xi) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Perturbation proportional to pitch angle.')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_xi2) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Perturbation proportional to function of pitch angle.')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_rh) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Maxwellian perturbation in ik=1 mode.')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_alf) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Ion dist fn proportional to v_parallel * sin(theta).')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_zero) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('Distribution function = 0.')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_test3) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_convect) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_restart_many) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="('Each PE restarts from its own NetCDF restart file.')") case (ginitopt_restart_small) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="('Each PE restarts from its own NetCDF restart file.')") write (report_unit, fmt="('with amplitudes scaled by factor of scale = ',f10.4)") scale write (report_unit, fmt="('Noise added with amplitude = ',f10.4)") phiinit case (ginitopt_continue) write (report_unit, fmt="('Initial conditions:')") write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) case (ginitopt_cgyro) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="('Constant amplitude for ky!=0')") write (report_unit, fmt="('Zero amplitude for ky=0')") case (ginitopt_random_sine) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="(' Zonal amplitude: ',f10.4)") zf_init write (report_unit, fmt="(' Max mode: ',I0)") max_mode write (report_unit, fmt="(' Random sum of sine waves up to max_mode mode number')") case (ginitopt_constant) write (report_unit, fmt="('Initial conditions:')") write (report_unit, fmt="(' Amplitude: ',f10.4)") phiinit write (report_unit, fmt="(' Zonal amplitude: ',f10.4)") zf_init write (report_unit, fmt="(' Constant value')") end select if (ginitopt_switch == ginitopt_restart_many) then if (delt_option_switch == delt_option_auto) then write (report_unit, *) write (report_unit, fmt="('This run is a continuation of a previous run.')") write (report_unit, fmt="('The time step at the beginning of this run')") write (report_unit, fmt="('will be taken from the end of the previous run.')") else write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('This run is a continuation of a previous run.')") write (report_unit, fmt="('The time step is being set by hand.')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('You probably want to set delt_option to be check_restart in the knobs namelist.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if end if if (delt_option_switch == delt_option_auto) then if (ginitopt_switch /= ginitopt_restart_many) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('This is not a normal continuation run.')") write (report_unit, fmt="('You probably want to set delt_option to be default in the knobs namelist.')") write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if end if end subroutine check_init_g !> FIXME : Add documentation subroutine init_init_g(init_g_config_in) use gs2_save, only: set_restart_file, read_many, proc_to_save_fields use gs2_layouts, only: init_gs2_layouts use mp, only: proc0, broadcast, job, nproc implicit none type(init_g_config_type), intent(in), optional :: init_g_config_in integer :: ind_slash if (initialized) return initialized = .true. call init_gs2_layouts call read_parameters(init_g_config_in) ! prepend restart_dir to restart_file ! append trailing slash if not exists if(restart_dir(len_trim(restart_dir):) /= "/") & restart_dir=trim(restart_dir)//"/" !Determine if restart file contains "/" if so split on this point to give DIR//FILE !so restart files are created in DIR//restart_dir//FILE ind_slash=index(restart_file,"/",.True.) if (ind_slash.EQ.0) then !No slash present restart_file=trim(restart_dir)//trim(restart_file) else !Slash present restart_file=trim(restart_file(1:ind_slash))//trim(restart_dir)//trim(restart_file(ind_slash+1:)) endif ! MAB - allows for ensemble averaging of multiple flux tube calculations ! job=0 if not doing multiple flux tube calculations, so phiinit unaffected !<DD>Only let proc0 do this and then broadcast as other procs don't know phiinit,phifrac etc. if(proc0) phiinit = phiinit * (job*phifrac+1.0) call broadcast (ginitopt_switch) call broadcast (width0) call broadcast (refac) call broadcast (imfac) call broadcast (den0) call broadcast (upar0) call broadcast (tpar0) call broadcast (tperp0) call broadcast (den1) call broadcast (upar1) call broadcast (tpar1) call broadcast (tperp1) call broadcast (den2) call broadcast (upar2) call broadcast (tpar2) call broadcast (tperp2) call broadcast (phiinit) call broadcast (phifrac) call broadcast (dphiinit) call broadcast (zf_init) call broadcast (init_spec_pow) call broadcast (apar0) call broadcast (tstart) call broadcast (chop_side) call broadcast (even) call broadcast (left) call broadcast (max_mode) call broadcast (clean_init) call broadcast (restart_file) call broadcast (proc_to_save_fields) call broadcast (read_many) call broadcast (ikk) call broadcast (itt) call broadcast (ikkk) call broadcast (ittt) call broadcast (phiamp) call broadcast (aparamp) call broadcast (scale) call broadcast (new_field_init) ! RN> call broadcast (phiinit0) call broadcast (phiinit_rand) call broadcast (a0) call broadcast (b0) call broadcast (null_phi) call broadcast (null_bpar) call broadcast (null_apar) call broadcast (adj_spec) call broadcast (eq_type) call broadcast (prof_width) call broadcast (eq_mode_n) call broadcast (eq_mode_u) call broadcast (input_check_recon) call broadcast (ukxy_pt) call broadcast (ikpar_init) call broadcast (ikx_init) call broadcast (kpar_init) call broadcast (restart_eig_id) call broadcast (constant_random_flag) ! <RN call set_restart_file (restart_file) if (proc_to_save_fields >= 0) proc_to_save_fields = mod(proc_to_save_fields, nproc) end subroutine init_init_g !> FIXME : Add documentation subroutine ginit (restarted, override_ginitopt_switch) use gs2_save, only: read_t0_from_restart_file use species, only: has_hybrid_electron_species, spec use mp, only: proc0 use optionals, only: get_option_with_default implicit none logical, intent (out) :: restarted integer, intent(in), optional :: override_ginitopt_switch logical, save :: have_warned_about_hybrid = .false. real :: t0 integer :: istatus integer :: ginitopt_actual ginitopt_actual = get_option_with_default(override_ginitopt_switch, ginitopt_switch) restarted = .false. select case (ginitopt_actual) case (ginitopt_default) call ginit_default case (ginitopt_default_odd) call ginit_default_odd case (ginitopt_kz0) call ginit_kz0 case (ginitopt_noise) call ginit_noise case (ginitopt_single_parallel_mode) call ginit_single_parallel_mode case (ginitopt_all_modes_equal) call ginit_all_modes_equal case (ginitopt_kpar) call ginit_kpar case (ginitopt_gs) call ginit_gs case (ginitopt_nl) call ginit_nl case (ginitopt_nl2) call ginit_nl2 case (ginitopt_nl3) call ginit_nl3 case (ginitopt_nl3r) call ginit_nl3r case (ginitopt_harris) call ginit_harris case (ginitopt_nl4) ! in an old version, this line was commented out. Thus, to recover some old ! results, you might need to comment out this line and... ! t0 = tstart call ginit_nl4 call read_t0_from_restart_file (tstart, istatus) ! this line: ! tstart = t0 restarted = .true. scale = 1. case (ginitopt_nl5) t0 = tstart call read_t0_from_restart_file (tstart, istatus) call ginit_nl5 tstart = t0 restarted = .true. scale = 1. case (ginitopt_recon) t0 = tstart call read_t0_from_restart_file (tstart, istatus) call ginit_recon tstart = t0 restarted = .true. scale = 1. case (ginitopt_nl6) t0 = tstart call read_t0_from_restart_file (tstart, istatus) call ginit_nl6 tstart = t0 restarted = .true. scale = 1. case (ginitopt_nl7) call ginit_nl7 case (ginitopt_xi) call ginit_xi case (ginitopt_xi2) call ginit_xi2 case (ginitopt_rh) call ginit_rh case (ginitopt_alf) call ginit_alf case (ginitopt_zero) call ginit_zero case (ginitopt_test3) call ginit_test3 case (ginitopt_convect) call ginit_convect case (ginitopt_restart_many) call ginit_restart_many call read_t0_from_restart_file (tstart, istatus) restarted = .true. scale = 1. case (ginitopt_restart_eig) call ginit_restart_eig call read_t0_from_restart_file (tstart, istatus) restarted = .true. scale = 1. case (ginitopt_restart_small) call ginit_restart_small call read_t0_from_restart_file (tstart, istatus) restarted = .true. scale = 1. case (ginitopt_zonal_only) call ginit_restart_zonal_only call read_t0_from_restart_file (tstart, istatus) restarted = .true. scale = 1. case (ginitopt_no_zonal) call ginit_restart_no_zonal call read_t0_from_restart_file (tstart, istatus) restarted = .true. scale = 1. case (ginitopt_smallflat) call ginit_restart_smallflat call read_t0_from_restart_file (tstart, istatus) restarted = .true. scale = 1. case (ginitopt_continue) restarted = .true. scale = 1. case (ginitopt_recon3) call ginit_recon3 case (ginitopt_ot) call ginit_ot case (ginitopt_stationary_mode) call ginit_stationary_mode case (ginitopt_cgyro) call ginit_cgyro case (ginitopt_random_sine) call ginit_random_sine case (ginitopt_constant) call ginit_constant end select ! If the initial condition is the nonadiabatic dfn then ! ensure we've actually set `gnew` correctly. ! We avoid doing this if we're restarting as currently we assume ! that the restart files always contain our evolved gnew, rather ! than the nonadiabatic dfn. if (initial_condition_is_nonadiabatic_dfn .and. .not. restarted) then call set_gnew_from_nonadiabatic_dfn else ! If not restarting and have a hybrid electron species then ! issue a warning (once only) as recommend initialising the ! non-adiabatic dfn directly in this situation as this respects ! the hybrid electron h = 0 condition. if ((.not. restarted) .and. has_hybrid_electron_species(spec) & .and. (.not. have_warned_about_hybrid) .and. proc0) then have_warned_about_hybrid = .true. write(*,'("Warning : Running with a hybrid electron species but initial_condition_is_nonadiabatic_dfn is false.")') end if end if end subroutine ginit !> Take the existing gnew value and, assuming it is the nonadiabatic dfn, !> calculate the consistent potentials to set the true gnew value. subroutine set_gnew_from_nonadiabatic_dfn use dist_fn_arrays, only: gnew, g, g_adjust, to_g_gs2 use dist_fn, only: calculate_potentials_from_nonadiabatic_dfn use species, only: spec, has_hybrid_electron_species use le_grids, only: is_passing_hybrid_electron use gs2_layouts, only: g_lo, ik_idx, il_idx, is_idx use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew use array_utils, only: copy implicit none integer :: iglo, ik, il, is ! Respect passing hybrid electrons h=0 if (has_hybrid_electron_species(spec)) then do iglo = g_lo%llim_proc, g_lo%ulim_proc is = is_idx(g_lo, iglo) ik = ik_idx(g_lo, iglo) il = il_idx(g_lo, iglo) if (is_passing_hybrid_electron(is, ik, il)) gnew(:,:,iglo) = 0.0 end do end if ! Calculate the consistent potentials call calculate_potentials_from_nonadiabatic_dfn(gnew, phi, apar, bpar) ! Adjust the nonadiabatic dfn to give the actual gnew we evolve call g_adjust(gnew, phi, bpar, direction = to_g_gs2) ! Ensure both timesteps are set correctly call copy(gnew, g) call copy(phi, phinew) call copy(apar, aparnew) call copy(bpar, bparnew) end subroutine set_gnew_from_nonadiabatic_dfn !> Set the smart defaults for the init_g_config_type instance subroutine set_smart_defaults_local(self) use file_utils, only: run_name implicit none class(init_g_config_type), intent(in out) :: self if (self%is_initialised()) return self%restart_file = trim(run_name)//".nc" end subroutine set_smart_defaults_local !> FIXME : Add documentation subroutine read_parameters(init_g_config_in) use file_utils, only: error_unit use text_options, only: text_option, get_option_value use gs2_save, only: read_many, include_explicit_source_in_restart, proc_to_save_fields implicit none type(init_g_config_type), intent(in), optional :: init_g_config_in type (text_option), dimension (*), parameter :: ginitopts = & (/ text_option('default', ginitopt_default), & text_option('default_odd', ginitopt_default_odd), & text_option('noise', ginitopt_noise), & text_option('xi', ginitopt_xi), & text_option('xi2', ginitopt_xi2), & text_option('zero', ginitopt_zero), & text_option('test3', ginitopt_test3), & text_option('convect', ginitopt_convect), & text_option('rh', ginitopt_rh), & text_option('restart', ginitopt_restart_many), & text_option('many', ginitopt_restart_many), & text_option('small', ginitopt_restart_small), & text_option('file', ginitopt_restart_many), & text_option('cont', ginitopt_continue), & text_option('kz0', ginitopt_kz0), & text_option('nl', ginitopt_nl), & text_option('nl2', ginitopt_nl2), & text_option('nl3', ginitopt_nl3), & text_option('nl3r', ginitopt_nl3r), & text_option('nl4', ginitopt_nl4), & text_option('nl5', ginitopt_nl5), & text_option('nl6', ginitopt_nl6), & text_option('nl7', ginitopt_nl7), & text_option('alf', ginitopt_alf), & text_option('gs', ginitopt_gs), & text_option('kpar', ginitopt_kpar), & text_option('smallflat', ginitopt_smallflat), & text_option('harris', ginitopt_harris), & text_option('recon', ginitopt_recon), & text_option('recon3', ginitopt_recon3), & text_option('ot', ginitopt_ot), & text_option('zonal_only', ginitopt_zonal_only), & text_option('no_zonal', ginitopt_no_zonal), & text_option('single_parallel_mode', ginitopt_single_parallel_mode), & text_option('all_modes_equal', ginitopt_all_modes_equal), & text_option('eig_restart', ginitopt_restart_eig), & text_option('stationary_mode', ginitopt_stationary_mode), & text_option('cgyro', ginitopt_cgyro), & text_option('random_sine', ginitopt_random_sine), & text_option('constant', ginitopt_constant) & /) character(20) :: ginit_option integer :: ierr if (present(init_g_config_in)) init_g_config = init_g_config_in call init_g_config%init(name = 'init_g_knobs', requires_index = .false.) ! Copy out internal values into module level parameters associate(self => init_g_config) #include "init_g_copy_out_auto_gen.inc" end associate ierr = error_unit() call get_option_value & (ginit_option, ginitopts, ginitopt_switch, & ierr, "ginit_option in ginit_knobs",.true.) end subroutine read_parameters !> Set initial condition to gaussian in \(\theta\) !> !> Controlled by: !> - Amplitude: [[init_g_knobs:phiinit]] !> - Width in \(\theta\): [[init_g_knobs:width0]] !> - Parity: even about \(\theta_0\) unless [[init_g_knobs:chop_side]] is true subroutine ginit_default use species, only: spec use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, theta0, aky, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi logical :: right integer :: iglo integer :: ig, ik, it, il, is right = .not. left do ig = -ntgrid, ntgrid phi(ig,:,:) = exp(-((theta(ig)-theta0(:,:))/width0)**2)*cmplx(1.0,1.0) end do if (chop_side .and. left) phi(:-1,:,:) = 0.0 if (chop_side .and. right) phi(1:,:,:) = 0.0 if (reality) then phi(:,1,1) = 0.0 if (naky > 1 .and. aky(1) == 0.0) then phi(:,:,1) = 0.0 end if ! not used: ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) enddo end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phi(:,it,ik)*spec(is)%z*phiinit where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_default !> Set initial condition to gaussian in \(\theta\) !> !> Controlled by: !> - Amplitude: [[init_g_knobs:phiinit]] !> - Width in \(\theta\): [[init_g_knobs:width0]] !> - Parity: odd about \(\theta_0\) unless [[init_g_knobs:chop_side]] is true subroutine ginit_default_odd use species, only: spec use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, theta0, aky, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi logical :: right integer :: iglo integer :: ig, ik, it, il, is right = .not. left do ig = -ntgrid, ntgrid phi(ig,:,:) = sin((theta(ig)-theta0(:,:))/width0)*exp(-((theta(ig)-theta0(:,:))/width0)**2)*cmplx(1.0,1.0) end do if (chop_side .and. left) phi(:-1,:,:) = 0.0 if (chop_side .and. right) phi(1:,:,:) = 0.0 if (reality) then phi(:,1,1) = 0.0 if (naky > 1 .and. aky(1) == 0.0) then phi(:,:,1) = 0.0 end if ! not used: ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) enddo end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phi(:,it,ik)*spec(is)%z*phiinit where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_default_odd !> Initialise with only the kparallel = 0 mode. subroutine ginit_kz0 use species, only: spec use theta_grid, only: ntgrid use kt_grids, only: naky, ntheta0, aky, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi logical :: right integer :: iglo integer :: ik, it, il, is right = .not. left phi = cmplx(refac,imfac) if (chop_side .and. left) phi(:-1,:,:) = 0.0 if (chop_side .and. right) phi(1:,:,:) = 0.0 if (reality) then phi(:,1,1) = 0.0 if (naky > 1 .and. aky(1) == 0.0) then phi(:,:,1) = 0.0 end if ! not used: ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) enddo end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit*(den0 + sqrt(2.)*upar0*vpa(:,1,iglo) + tpar0*(vpa(:,1,iglo)**2-0.5) + tperp0*(vperp2(:,iglo)-1.)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit*(den0 + sqrt(2.)*upar0*vpa(:,2,iglo) + tpar0*(vpa(:,2,iglo)**2-0.5) + tperp0*(vperp2(:,iglo)-1.)) end do call copy(g, gnew) end subroutine ginit_kz0 !> FIXME : Add documentation subroutine single_initial_kx(phi) use theta_grid, only: ntgrid use kt_grids, only: naky, ntheta0 ! use kt_grids, only: akx use mp, only: mp_abort implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky), intent(inout) :: phi integer :: ig, ik, it if (ikx_init < 2 .or. ikx_init > (ntheta0+1)/2) then call mp_abort("The subroutine single_initial_kx should only be called when 1 < ikx_init < (ntheta0+1)/2") end if do it = 1, ntheta0 if (it .ne. ikx_init) then !write (*,*) "zeroing out kx_index: ", it, "at kx: ", akx(it) do ik = 1, naky do ig = -ntgrid, ntgrid phi(ig,it,ik) = 0. end do end do end if end do end subroutine single_initial_kx !> Initialise the distribution function with random noise. This is the !> recommended initialisation option. Each different mode is given a random !> amplitude between zero and one. subroutine ginit_noise use species, only: spec, tracer_species use theta_grid, only: ntgrid use kt_grids, only: naky, ntheta0, aky, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx, proc_id use dist_fn, only: l_links, r_links, pass_right, init_pass_ends, has_linked_boundary use redistribute, only: fill, delete_redist use ran, only: ranf use constant_random, only: init_constant_random use constant_random, only: finish_constant_random use constant_random, only: constant_ranf=>ranf use array_utils, only: copy, zero_array use constants, only: twopi use kt_grids, only: kperp2 implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit real :: a, b, spec_pow integer :: iglo, ig, ik, it, il, is, nn logical :: linked_bc linked_bc = has_linked_boundary() !CMR: need to document tracer phit parameter ;-) call zero_array(phit) do it=2,ntheta0/2+1 nn = it-1 ! extra factor of 4 to reduce the high k wiggles for now phit (:, it, 1) = (-1)**nn*exp(-8.*(real(nn)/ntheta0)**2) end do ! keep old (it, ik) loop order to get old results exactly: !Fill phi with random (complex) numbers between -0.5 and 0.5 if (constant_random_flag) call init_constant_random do it = 1, ntheta0 do ik = 1, naky do ig = -ntgrid, ntgrid if (constant_random_flag) then a = constant_ranf()-0.5 b = constant_ranf()-0.5 else a = ranf()-0.5 b = ranf()-0.5 end if if (init_spec_pow > 0.) then spec_pow = 1./sqrt(1.+kperp2(ig,it,ik)**(0.5*init_spec_pow)) phi(ig,it,ik) = spec_pow*exp(cmplx(0.,twopi*a)) else phi(ig,it,ik) = cmplx(a,b) end if end do !CMR,28/1/13: ! clean_init debrutalises influence of chop_side on phi with linked bc if (chop_side) then if (left) then if (clean_init .and. linked_bc) then !Does this tube have more links to the right than the left? !If so then it's on the left so set to zero if (l_links(it, ik) .lt. r_links(it, ik)) then phi(:,it,ik) = 0.0 !This tube is in the middle so only set the left half to 0 else if (l_links(it, ik) .eq. r_links(it, ik)) then phi(:-1,it,ik) = 0.0 endif else phi(:-1,it,ik) = 0.0 endif else if (clean_init .and. linked_bc) then !Does this tube have more links to the left than the right? !If so then it's on the right so set to zero if (r_links(it, ik) .lt. l_links(it, ik)) then phi(:,it,ik) = 0.0 !This tube is in the middle so only set the right half to 0 else if (l_links(it, ik) .eq. r_links(it, ik)) then phi(1:,it,ik) = 0.0 endif else phi(1:,it,ik) = 0.0 endif endif end if end do end do if (constant_random_flag) call finish_constant_random !Wipe out all but one kx if requested if (ikx_init > 0) call single_initial_kx(phi) !Sort out the zonal/self-periodic modes if (naky .ge. 1 .and. aky(1) == 0.0) then !Apply scaling factor phi(:,:,1) = phi(:,:,1)*zf_init !Set ky=kx=0.0 mode to zero in amplitude phi(:,1,1) = 0.0 !CMR, 25/1/13: clean_init forces periodic zonal flows if (clean_init .and. linked_bc) then phi(ntgrid,:,1)=phi(-ntgrid,:,1) phit(ntgrid,:,1)=phit(-ntgrid,:,1) end if end if !Apply reality condition (i.e. -kx mode is conjugate of +kx mode) if (reality) then do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) phit(:,it+(ntheta0+1)/2,1) = conjg(phit(:,(ntheta0+1)/2+1-it,1)) enddo end if !Now set g using data in phi do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) !Handle tracer_species if (spec(is)%type == tracer_species) then g(:,1,iglo) =-phit(:,it,ik)*spec(is)%z*phiinit else g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit end if !Make sure there's no g in the forbidden regions where (forbid(:,il)) g(:,1,iglo) = 0.0 !Set incoming and outgoing boundary conditions if ( clean_init .and. linked_bc .and. ik > 1) then !CMR: clean_init sets g=0 for incoming particles, ! and outgoing particles ! as initial condition sets g(:,1,:) = g(:,2,:) ) !If no links to the left then we're at the left boundary if ( l_links(it, ik) .eq. 0 ) g(-ntgrid,1,iglo)=0.0 !If no links to the right then we're at the right boundary if ( r_links(it, ik) .eq. 0 ) g(ntgrid,1,iglo)=0.0 endif end do !If clean_init is set and there are any links/connections then make sure repeated points are consistent if ( clean_init .and. linked_bc .and. sum(r_links+l_links) > 0 ) then ! !CMR, 23/1/13: ! clean_init: ensure continuity in || direction with linked bcs ! set up redistribute type pass_right ! this passes g(ntgrid,1,iglo) to g(-ntgrid,1,iglo_r) on right neighbour !Initialise communications object !NOTE This should be moved to dist_fn in future call init_pass_ends(pass_right,'r',1,'c') !Pass right boundaries (ntgrid) of g to linked left boundaries (-ntgrid) of g call fill(pass_right,g,g) !Deallocate comm object as not used anywhere else !NOTE This will need removing if/when initialisation moved to dist_fn call delete_redist(pass_right) !Make leftwards and rightwards particles the same g(:,2,:)=g(:,1,:) !Set gnew to be the "fixed" data call copy(g, gnew) else ! finally initialise isign=2 g(:,2,:) = g(:,1,:) call copy(g, gnew) endif end subroutine ginit_noise !> Initialise only with a single parallel mode \(k_\Vert\) specified by either !> [[init_g_knobs:ikpar_init]] for periodic boundary conditions or !> [[init_g_knobs:kpar_init]] for linked boundary conditions. Only makes sense !> for linear calculations. subroutine ginit_single_parallel_mode use species, only: spec, tracer_species use theta_grid, only: ntgrid, is_effectively_zero_shear, theta use kt_grids, only: naky, ntheta0, aky, reality use le_grids, only: forbid use mp, only: mp_abort use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit real :: a, b integer :: iglo integer :: ig, ik, it, il, is, nn call zero_array(phit) do it=2,ntheta0/2+1 nn = it-1 ! extra factor of 4 to reduce the high k wiggles for now phit (:, it, 1) = (-1)**nn*exp(-8.*(real(nn)/ntheta0)**2) end do if (is_effectively_zero_shear()) then if (ikpar_init+1 > ntgrid) then call mp_abort("Error: this value of k_parallel is too large. Increase ntheta or decrease ikpar_init.") end if kpar_init = ikpar_init end if do it = 1, ntheta0 do ik = 1, naky do ig = -ntgrid, ntgrid !Set the field to cos(kpar_init*theta), where we remember that the gridpoints are not necessarily evenly spaced in the parallel direction, so we use theta(ig). !reality for ky=0 means we must use -kpar for kx < 0 !if (naky == 1 .and. it > (ntheta0+1)/2) then !a = cos(-kpar_init * theta(ig)) !b = sin(-kpar_init * theta(ig)) !else a = cos(kpar_init * theta(ig)) b = sin(kpar_init * theta(ig)) !end if ! a = ranf()-0.5 ! b = ranf()-0.5 ! phi(:,it,ik) = cmplx(a,b) phi(ig,it,ik) = cmplx(a,b) end do if (chop_side) then if (left) then phi(:-1,it,ik) = 0.0 else phi(1:,it,ik) = 0.0 end if end if end do end do if (naky > 1 .and. aky(1) == 0.0) then phi(:,:,1) = phi(:,:,1)*zf_init end if if (ikx_init > 0) call single_initial_kx(phi) !<doc> reality condition for ky = 0 component: </doc> if (reality) then do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) phit(:,it+(ntheta0+1)/2,1) = conjg(phit(:,(ntheta0+1)/2+1-it,1)) enddo end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) if (spec(is)%type == tracer_species) then g(:,1,iglo) =-phit(:,it,ik)*spec(is)%z*phiinit else g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit end if where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_single_parallel_mode !> Initialize with every parallel and perpendicular mode having equal amplitude. !! Only makes sense in a linear calculation. k_parallel is specified with kpar_init !! or with ikpar_init when periodic boundary conditions are used. EGH subroutine ginit_all_modes_equal use species, only: spec, tracer_species use theta_grid, only: ntgrid, theta, ntheta use kt_grids, only: naky, ntheta0, aky, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, phit real :: a, b integer :: iglo integer :: ig, ik, it, il, is, nn, ikpar call zero_array(phit) do it=2,ntheta0/2+1 nn = it-1 ! extra factor of 4 to reduce the high k wiggles for now phit (:, it, 1) = (-1)**nn*exp(-8.*(real(nn)/ntheta0)**2) end do do it = 1, ntheta0 do ik = 1, naky do ig = -ntgrid, ntgrid ! Set the field to cos(kpar*theta) for all kpar, where we remember that the gridpoints are not necessarily evenly spaced in the parallel direction, so we use theta(ig)</doc> a = 0.0 b = 0.0 do ikpar = 0, ntheta - 1 a = a + cos(ikpar * theta(ig)) ! we want to include the positive and negative wave numbers in ! equal measure, which of course means a real sine wave. b = 0.0 !b + cos(ikpar * theta(ig)) end do ! a = ranf()-0.5 ! b = ranf()-0.5 ! phi(:,it,ik) = cmplx(a,b) phi(ig,it,ik) = cmplx(a,b) end do if (chop_side) then if (left) then phi(:-1,it,ik) = 0.0 else phi(1:,it,ik) = 0.0 end if end if end do end do if (naky > 1 .and. aky(1) == 0.0) then phi(:,:,1) = phi(:,:,1)*zf_init end if if (ikx_init > 0) call single_initial_kx(phi) !<doc> reality condition for ky = 0 component: </doc> if (reality) then do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) phit(:,it+(ntheta0+1)/2,1) = conjg(phit(:,(ntheta0+1)/2+1-it,1)) enddo end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) if (spec(is)%type == tracer_species) then g(:,1,iglo) =-phit(:,it,ik)*spec(is)%z*phiinit else g(:,1,iglo) = -phi(:,it,ik)*spec(is)%z*phiinit end if where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_all_modes_equal !> "nl": Two \(k_\perp\) modes set to constant amplitude. Controlled by: !> !> - Overall amplitude: [[init_g_knobs:phiinit]] !> - `ky` indices of modes: [[init_g_knobs:ikk]] !> - `kx` indices of modes: [[init_g_knobs:itt]] !> - Parity: even about \(\theta_0\) unless !> [[init_g_knobs:chop_side]] is true, in which case none subroutine ginit_nl use species, only: spec use theta_grid, only: ntgrid use kt_grids, only: naky, ntheta0 use le_grids, only: forbid, il_is_wfb use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi integer :: iglo integer :: ig, ik, it, il, is, j call zero_array(phi) do j = 1, 2 ik = ikk(j) it = itt(j) do ig = -ntgrid, ntgrid ! phi(ig,it,ik) = cmplx(1.0, 0.0)*sin(theta(ig)) phi(ig,it,ik) = cmplx(1.0, 0.0)!*sin(theta(ig)) end do ! do ig = -ntgrid/2, ntgrid/2 ! phi(ig,it,ik) = cmplx(1.0, 1.0) ! end do if (chop_side) then if (left) then phi(:-1,it,ik) = 0.0 else phi(1:,it,ik) = 0.0 end if end if end do ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) enddo ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = -phi(:,it,ik)*phiinit*spec(is)%z where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) if (il_is_wfb(il)) g(:,:,iglo) = 0.0 end do call copy(g, gnew) end subroutine ginit_nl !> "nl2": Two \(k_\perp\) modes proportional to \(1 + v_\Vert !> \sin(\theta)\). Controlled by: !> !> - Overall amplitude: [[init_g_knobs:phiinit]] !> - `ky` indices of modes: [[init_g_knobs:ikk]] !> - `kx` indices of modes: [[init_g_knobs:itt]] !> - Parity: even about \(\theta_0\) unless !> [[init_g_knobs:chop_side]] is true, in which case none subroutine ginit_nl2 use species, only: spec use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0 use le_grids, only: forbid, il_is_wfb use dist_fn_arrays, only: g, gnew, vpa use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi integer :: iglo integer :: ig, ik, it, il, is, j call zero_array(phi) do j = 1, 2 ik = ikk(j) it = itt(j) do ig = -ntgrid, ntgrid phi(ig,it,ik) = cmplx(1.0, 0.0)!*sin(theta(ig)) end do ! do ig = -ntgrid/2, ntgrid/2 ! phi(ig,it,ik) = cmplx(1.0, 1.0) ! end do if (chop_side) then if (left) then phi(:-1,it,ik) = 0.0 else phi(1:,it,ik) = 0.0 end if end if end do ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) enddo ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = -phi(:,it,ik)*phiinit*(1.+vpa(:,1,iglo)*sin(theta))*spec(is)%z where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = -phi(:,it,ik)*phiinit*(1.+vpa(:,2,iglo)*sin(theta))*spec(is)%z ! g(:,1,iglo)*vpa(:,2,iglo) if (il_is_wfb(il)) g(:,:,iglo) = 0.0 end do call copy(g, gnew) end subroutine ginit_nl2 !> "nl3": Two \(k_\perp\) modes with a Gaussian envelope, and perturbation in !> fields of form \(x_0 + x_1 \cos(\theta) + x_2 \cos(2 \theta)\), except for !> parallel velocity which uses is 90 degrees out of phase and uses !> \(\sin\). Density, parallel velocity, and perpendicular and parallel !> temperatures can all be controlled independently. Controlled by: !> !> - Overall amplitude: [[init_g_knobs:phiinit]] !> - Amplitude of real and imaginary parts of \(\phi\): !> [[init_g_knobs:refac]], [[init_g_knobs:imfac]] !> - Gaussian envelope width in \(\theta\): [[init_g_knobs:width0]] !> - `ky` indices of modes: [[init_g_knobs:ikk]] !> - `kx` indices of modes: [[init_g_knobs:itt]] !> - Parity: even about \(\theta_0\) unless !> [[init_g_knobs:chop_side]] is true, in which case none !> - Density perturbation: [[init_g_knobs:den0]], !> [[init_g_knobs:den1]], [[init_g_knobs:den2]] !> - Parallel velocity perturbation: [[init_g_knobs:upar0]], !> [[init_g_knobs:upar1]], [[init_g_knobs:upar2]] !> - Parallel temperature perturbation: [[init_g_knobs:tpar0]], !> [[init_g_knobs:tpar1]], [[init_g_knobs:tpar2]] !> - Perpendicular temperature perturbation: [[init_g_knobs:tperp0]], !> [[init_g_knobs:tperp1]], [[init_g_knobs:tperp2]] subroutine ginit_nl3 use species, only: spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, theta0, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use constants, only: zi use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t integer :: iglo integer :: ig, ik, it, il, is, j call zero_array(phi) ; call zero_array(odd) do j = 1, 2 ik = ikk(j) it = itt(j) if (width0 > 0.) then do ig = -ntgrid, ntgrid phi(ig,it,ik) = exp(-((theta(ig)-theta0(it,ik))/width0)**2)*cmplx(refac, imfac) end do else do ig = -ntgrid, ntgrid phi(ig,it,ik) = cmplx(refac, imfac) end do end if if (chop_side) then if (left) then phi(:-1,it,ik) = 0.0 else phi(1:,it,ik) = 0.0 end if end if end do odd = zi * phi ! reality condition for k_theta = 0 component: if (reality) then do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1)) enddo end if if (even) then ct = cos(theta) st = sin(theta) c2t = cos(2.*theta) s2t = sin(2.*theta) else ct = sin(theta) st = cos(theta) c2t = sin(2.*theta) s2t = cos(2.*theta) end if dfac = den0 + den1 * ct + den2 * c2t ufac = upar0 + upar1* st + upar2* s2t tparfac = tpar0 + tpar1* ct + tpar2* c2t tperpfac = tperp0 + tperp1*ct + tperp2*c2t ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) ! if (is_electron_species(spec(is))) cycle g(:,1,iglo) = phiinit* &!spec(is)%z* & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,1,iglo) * odd(:,it,ik) & + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = phiinit* &!spec(is)%z* & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,2,iglo) * odd(:,it,ik) & + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,2,iglo) = 0.0 ! if (il == ng2+1) g(:,:,iglo) = 0.0 end do ! if (has_electron_species(spec)) then ! call flae (g, gnew) ! g = g - gnew ! end if call copy(g, gnew) end subroutine ginit_nl3 !> FIXME : Add documentation subroutine ginit_nl3r use species, only: spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, theta0, reality use le_grids, only: forbid use fields_arrays, only: apar use fields_arrays, only: aparnew use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use constants, only: zi use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t integer :: iglo integer :: ig, ik, it, il, is, j call zero_array(phi) ; call zero_array(odd) do j = 1, 2 ik = ikk(j) it = itt(j) if (width0 > 0.) then do ig = -ntgrid, ntgrid phi(ig,it,ik) = exp(-((theta(ig)-theta0(it,ik))/width0)**2)*cmplx(refac, imfac) end do else do ig = -ntgrid, ntgrid phi(ig,it,ik) = cmplx(refac, imfac) apar(ig,it,ik) = apar0*cmplx(refac, imfac) end do end if if (chop_side) then if (left) then phi(:-1,it,ik) = 0.0 else phi(1:,it,ik) = 0.0 end if end if end do odd = zi * phi ! reality condition for k_theta = 0 component: if (reality) then do it = 1, ntheta0/2 apar(:,it+(ntheta0+1)/2,1) = conjg(apar(:,(ntheta0+1)/2+1-it,1)) phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1)) enddo end if aparnew = apar if (even) then ct = cos(theta) st = sin(theta) c2t = cos(2.*theta) s2t = sin(2.*theta) else ct = sin(theta) st = cos(theta) c2t = sin(2.*theta) s2t = cos(2.*theta) end if dfac = den0 + den1 * ct + den2 * c2t ufac = upar0 + upar1* st + upar2* s2t tparfac = tpar0 + tpar1* ct + tpar2* c2t tperpfac = tperp0 + tperp1*ct + tperp2*c2t ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phiinit* &!spec(is)%z* & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,1,iglo)*spec(is)%u0 * odd(:,it,ik) & + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = phiinit* &!spec(is)%z* & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,2,iglo)*spec(is)%u0 * odd(:,it,ik) & + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,2,iglo) = 0.0 ! if (il == ng2+1) g(:,:,iglo) = 0.0 end do ! if (has_electron_species(spec)) then ! call flae (g, gnew) ! g = g - gnew ! end if call copy(g, gnew) end subroutine ginit_nl3r !> FIXME : Add documentation subroutine ginit_harris use species, only: spec, has_electron_species use kt_grids, only: naky, ntheta0, reality, nx use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa, aj0 use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx use constants, only: zi, pi use run_parameters, only: k0 use array_utils, only: copy, zero_array implicit none complex, dimension (ntheta0,naky) :: phi integer :: iglo integer :: ik, it, il, is, j, j1 real, dimension(nx) :: lx_pr, a_pr complex,dimension(nx) :: ff_pr real:: L, dx_pr ! ! Specifying function on x grid, with nx points. Transforming to kx grid, with ! nkx < nx points. Result is function that will be used for initial condition. ! But Fourier transform is the same if you just use nkx points in the first ! place, right? ! ! Can specify x0 along with y0; no need to use this construction, although it is okay. L=2.*pi*k0 ! nx is available from kt_grids: dx_pr=L/real(nx) do j = 1, nx lx_pr(j)=dx_pr*real(j-1) end do do j = 1,nx a_pr(j)=1./cosh((lx_pr(j)-L/4.)/width0)- & 1./cosh((lx_pr(j)-3.*L/4.)/width0) end do a_pr=a_pr/real(nx) do j = 1,nx ff_pr(j) = 0. do j1 = 1,nx ff_pr(j)=ff_pr(j)+a_pr(j1)*exp(-zi*2.*pi*real(j-1)*real(j1-1)/real(nx)) end do end do call zero_array(phi) ! Dealiasing here: do j = 1, ntheta0/2-1 phi(j,1) = ff_pr(j) end do ! ! Note: if the j=1 coefficient is non-zero, it might be a problem, because this ! coefficient is never used. ! do j = ntheta0/2, ntheta0 phi(j,1) = ff_pr(nx-ntheta0+j) end do ! ! presumably this is not needed, but leave it in for now: ! ! reality condition for k_theta = 0 component: if (reality) then do it = 1, ntheta0/2 phi(it+(ntheta0+1)/2,1) = conjg(phi((ntheta0+1)/2+1-it,1)) enddo end if call zero_array(g) do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) ! ions, kx/=0: if ((is==1) .and.(.not.(it==1))) then g(:,1,iglo) = 2.* vpa(:,1,iglo)*spec(is)%u0 * phi(it,ik)/aj0(:,iglo) g(:,2,iglo) = 2.* vpa(:,2,iglo)*spec(is)%u0 * phi(it,ik)/aj0(:,iglo) where (forbid(:,il)) g(:,1,iglo) = 0.0 where (forbid(:,il)) g(:,2,iglo) = 0.0 end if end do call copy(g, gnew) end subroutine ginit_harris !> FIXME : Add documentation subroutine ginit_recon use mp, only: proc0 use species, only: spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, reality use le_grids, only: forbid use gs2_save, only: gs2_restore use dist_fn_arrays, only: g, gnew, vpa, vperp2 use fields_arrays, only: phi, apar, bpar use fields_arrays, only: phinew, aparnew, bparnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use file_utils, only: error_unit use run_parameters, only: has_phi, has_apar, has_bpar use constants, only: zi use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz, odd real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t integer :: iglo, istatus, ierr integer :: ig, ik, it, il, is, j call gs2_restore (g, scale, istatus, has_phi, has_apar, has_bpar) if (istatus /= 0) then ierr = error_unit() if (proc0) write(ierr,*) "Error reading file: ", trim(restart_file) call zero_array(g) end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) if (ik == 1) cycle g (:,1,iglo) = 0. g (:,2,iglo) = 0. end do phinew(:,:,2:naky) = 0. aparnew(:,:,2:naky) = 0. bparnew(:,:,2:naky) = 0. phi(:,:,2:naky) = 0. apar(:,:,2:naky) = 0. bpar(:,:,2:naky) = 0. call zero_array(phiz) ; call zero_array(odd) do j = 1, 2 ik = ikk(j) it = itt(j) do ig = -ntgrid, ntgrid phiz(ig,it,ik) = cmplx(refac, imfac) end do end do odd = zi * phiz ! reality condition for k_theta = 0 component: if (reality) then do it = 1, ntheta0/2 phiz(:,it+(ntheta0+1)/2,1) = conjg(phiz(:,(ntheta0+1)/2+1-it,1)) odd (:,it+(ntheta0+1)/2,1) = conjg(odd (:,(ntheta0+1)/2+1-it,1)) enddo end if if (even) then ct = cos(theta) st = sin(theta) c2t = cos(2.*theta) s2t = sin(2.*theta) else ct = sin(theta) st = cos(theta) c2t = sin(2.*theta) s2t = cos(2.*theta) end if dfac = den0 + den1 * ct + den2 * c2t ufac = upar0 + upar1* st + upar2* s2t tparfac = tpar0 + tpar1* ct + tpar2* c2t tperpfac = tperp0 + tperp1*ct + tperp2*c2t ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phiinit* &!spec(is)%z* & ( dfac*spec(is)%dens0 * phiz(:,it,ik) & + 2.*ufac* vpa(:,1,iglo) * odd (:,it,ik) & + tparfac*(vpa(:,1,iglo)**2-0.5) * phiz(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phiz(:,it,ik)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = phiinit* &!spec(is)%z* & ( dfac*spec(is)%dens0 * phiz(:,it,ik) & + 2.*ufac* vpa(:,2,iglo) * odd (:,it,ik) & + tparfac*(vpa(:,2,iglo)**2-0.5) * phiz(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phiz(:,it,ik)) where (forbid(:,il)) g(:,2,iglo) = 0.0 ! if (il == ng2+1) g(:,:,iglo) = 0.0 end do ! if (has_electron_species(spec)) then ! call flae (g, gnew) ! g = g - gnew ! end if call copy(g, gnew) end subroutine ginit_recon !> FIXME : Add documentation !! !! @note nl4 has been monkeyed with over time. Do not expect to recover old results !! that use this startup routine. subroutine ginit_nl4 use mp, only: proc0 use species, only: spec use gs2_save, only: gs2_restore use theta_grid, only: ntgrid use kt_grids, only: naky, ntheta0 use dist_fn_arrays, only: g, gnew use le_grids, only: forbid use fields_arrays, only: phi, apar, bpar use fields_arrays, only: phinew, aparnew, bparnew use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx use file_utils, only: error_unit use run_parameters, only: has_phi, has_apar, has_bpar use ran, only: ranf use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz integer :: iglo, istatus integer :: ig, ik, it, is, il, ierr call gs2_restore (g, scale, istatus, has_phi, has_apar, has_bpar) if (istatus /= 0) then ierr = error_unit() if (proc0) write(ierr,*) "Error reading file: ", trim(restart_file) call zero_array(g) end if do iglo = g_lo%llim_proc, g_lo%ulim_proc it = it_idx(g_lo,iglo) ik = ik_idx(g_lo,iglo) ! if ((it == 2 .or. it == ntheta0) .and. ik == 1) cycle ! if (ik == 1) cycle if (it == 1 .and. ik == 2) cycle g (:,1,iglo) = 0. g (:,2,iglo) = 0. end do ! do ik = 1, naky ! if (ik /= 2) then ! phinew(:,:,ik) = 0. ! aparnew(:,:,ik) = 0. ! bparnew(:,:,ik) = 0. ! phi(:,:,ik) = 0. ! apar(:,:,ik) = 0. ! bpar(:,:,ik) = 0. ! else ! phinew(:,2:ntheta0,ik) = 0. ! aparnew(:,2:ntheta0,ik) = 0. ! bparnew(:,2:ntheta0,ik) = 0. ! phi(:,2:ntheta0,ik) = 0. ! apar(:,2:ntheta0,ik) = 0. ! bpar(:,2:ntheta0,ik) = 0. ! end if ! end do do ik = 1, naky do it=1,ntheta0 if (it == 1 .and. ik == 2) cycle phinew(:,it,ik) = 0. aparnew(:,it,ik) = 0. bparnew(:,it,ik) = 0. phi(:,it,ik) = 0. apar(:,it,ik) = 0. bpar(:,it,ik) = 0. end do end do do ik = 1, naky do it = 1, ntheta0 do ig = -ntgrid, ntgrid phiz(ig,it,ik) = cmplx(ranf()-0.5,ranf()-0.5) end do if (chop_side) then if (left) then phiz(:-1,it,ik) = 0.0 else phiz(1:,it,ik) = 0.0 end if end if end do end do ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phiz(:,it+(ntheta0+1)/2,1) = conjg(phiz(:,(ntheta0+1)/2+1-it,1)) enddo do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = g(:,1,iglo)-phiz(:,it,ik)*spec(is)%z*phiinit g(:,2,iglo) = g(:,2,iglo)-phiz(:,it,ik)*spec(is)%z*phiinit where (forbid(:,il)) g(:,1,iglo) = 0.0 where (forbid(:,il)) g(:,2,iglo) = 0.0 end do call copy(g, gnew) end subroutine ginit_nl4 !> FIXME : Add documentation subroutine ginit_nl5 use mp, only: proc0 use species, only: spec use gs2_save, only: gs2_restore use theta_grid, only: ntgrid use kt_grids, only: naky, ntheta0 use dist_fn_arrays, only: g, gnew use le_grids, only: forbid use fields_arrays, only: phi, apar, bpar use fields_arrays, only: phinew, aparnew, bparnew use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx use file_utils, only: error_unit use run_parameters, only: has_phi, has_apar, has_bpar use ran, only: ranf use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phiz integer :: iglo, istatus integer :: ig, ik, it, is, il, ierr call gs2_restore (g, scale, istatus, has_phi, has_apar, has_bpar) if (istatus /= 0) then ierr = error_unit() if (proc0) write(ierr,*) "Error reading file: ", trim(restart_file) call zero_array(g) end if do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) if (ik == 1) cycle g (:,1,iglo) = 0. g (:,2,iglo) = 0. end do phinew(:,:,2:naky) = 0. aparnew(:,:,2:naky) = 0. bparnew(:,:,2:naky) = 0. phi(:,:,2:naky) = 0. apar(:,:,2:naky) = 0. bpar(:,:,2:naky) = 0. do ik = 1, naky do it = 1, ntheta0 do ig = -ntgrid, ntgrid phiz(ig,it,ik) = cmplx(ranf()-0.5,ranf()-0.5) end do if (chop_side) then if (left) then phiz(:-1,it,ik) = 0.0 else phiz(1:,it,ik) = 0.0 end if end if end do end do do ik = 1, naky do it = 1, ntheta0 do ig = -ntgrid, ntgrid phiz(ig,it,ik) = cmplx(ranf()-0.5,ranf()-0.5) end do if (chop_side) then if (left) then phiz(:-1,it,ik) = 0.0 else phiz(1:,it,ik) = 0.0 end if end if end do end do ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phiz(:,it+(ntheta0+1)/2,1) = conjg(phiz(:,(ntheta0+1)/2+1-it,1)) enddo do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = g(:,1,iglo)-phiz(:,it,ik)*phiinit*spec(is)%z g(:,2,iglo) = g(:,2,iglo)-phiz(:,it,ik)*phiinit*spec(is)%z where (forbid(:,il)) g(:,1,iglo) = 0.0 where (forbid(:,il)) g(:,2,iglo) = 0.0 end do call copy(g, gnew) end subroutine ginit_nl5 !> FIXME : Add documentation subroutine ginit_nl6 use mp, only: proc0 use gs2_save, only: gs2_restore ! use theta_grid, only: ntgrid use kt_grids, only: ntheta0 ! use kt_grids, only: naky use dist_fn_arrays, only: g, gnew use fields_arrays, only: phi, apar, bpar use fields_arrays, only: phinew, aparnew, bparnew use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx use file_utils, only: error_unit use run_parameters, only: has_phi, has_apar, has_bpar use array_utils, only: copy, zero_array implicit none integer :: iglo, istatus integer :: ik, it, ierr call gs2_restore (g, scale, istatus, has_phi, has_apar, has_bpar) if (istatus /= 0) then ierr = error_unit() if (proc0) write(ierr,*) "Error reading file: ", trim(restart_file) call zero_array(g) end if call zero_array(gnew) do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) if (ik == 1 .and. it == 2) cycle if (ik == 1 .and. it == ntheta0) cycle g (:,:,iglo) = gnew(:,:,iglo) end do phinew(:,2,1) = 0. aparnew(:,2,1) = 0. bparnew(:,2,1) = 0. phi(:,2,1) = 0. apar(:,2,1) = 0. bpar(:,2,1) = 0. phinew(:,ntheta0,1) = 0. aparnew(:,ntheta0,1) = 0. bparnew(:,ntheta0,1) = 0. phi(:,ntheta0,1) = 0. apar(:,ntheta0,1) = 0. bpar(:,ntheta0,1) = 0. call copy(g, gnew) end subroutine ginit_nl6 !> FIXME : Add documentation subroutine ginit_nl7 use species, only: spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, reality use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use constants, only: zi use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac, ct, st, c2t, s2t integer :: iglo integer :: ig, ik, it, il, is, j call zero_array(phi) ; call zero_array(odd) do it=1,ntheta0 do ik=1,naky phi(:,it,ik) = cmplx(refac,imfac)*dphiinit end do end do do j = 1, 2 ik = ikk(j) it = itt(j) do ig = -ntgrid, ntgrid phi(ig,it,ik) = cmplx(refac, imfac) end do end do odd = zi * phi ! reality condition for k_theta = 0 component: if (reality) then do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1)) enddo end if if (even) then ct = cos(theta) st = sin(theta) c2t = cos(2.*theta) s2t = sin(2.*theta) else ct = sin(theta) st = cos(theta) c2t = sin(2.*theta) s2t = cos(2.*theta) end if dfac = den0 + den1 * ct + den2 * c2t ufac = upar0 + upar1* st + upar2* s2t tparfac = tpar0 + tpar1* ct + tpar2* c2t tperpfac = tperp0 + tperp1*ct + tperp2*c2t do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phiinit* & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,1,iglo) * odd(:,it,ik) & + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = phiinit* & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,2,iglo) * odd(:,it,ik) & + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,2,iglo) = 0.0 end do call copy(g, gnew) end subroutine ginit_nl7 !> Orszag-Tang 2D vortex problem subroutine ginit_ot use species, only: spec use theta_grid, only: ntgrid use kt_grids, only: naky, nakx => ntheta0, reality, kperp2 use dist_fn_arrays, only: g, gnew, vpa use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx use fields_arrays, only: phinew, aparnew, bparnew use dist_fn, only: get_init_field use array_utils, only: copy, zero_array implicit none integer :: iglo, ik, it, is, i complex :: fac complex, dimension (-ntgrid:ntgrid,nakx,naky) :: phi, jpar ! local ! real, dimension (-ntgrid:ntgrid) :: dfac, ufac !! phi, jpar are local !! call zero_array(phi) ; call zero_array(jpar) !!$ phi(:,1,2) = phiinit * cmplx(2.0, 0.0) ! 2 cos(y) !!$ phi(:,2,1) = phiinit * cmplx(1.0, 0.0) ! 2 cos(x) !!$ jpar(:,1,2) = apar0 * cmplx(2.0, 0.0) ! 2 cos(y) !!$ jpar(:,3,1) = apar0 * cmplx(2.0, 0.0) ! 4 cos(2x) do i=1, 3 it = ittt(i) ik = ikkk(i) phi(:,it,ik) = phiamp(i) jpar(:,it,ik) = aparamp(i) * kperp2(:,it,ik) end do ! reality condition for ky = 0 component: if (reality) then do it = 1, nakx/2 phi(:,it+(nakx+1)/2,1) = conjg(phi(:,(nakx+1)/2+1-it,1)) jpar(:,it+(nakx+1)/2,1) = conjg(jpar(:,(nakx+1)/2+1-it,1)) end do end if dfac = den0 ufac = upar0 do iglo = g_lo%llim_proc, g_lo%ulim_proc it = it_idx(g_lo,iglo) ik = ik_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,1,iglo)*spec(is)%u0 * jpar(:,it,ik) ) g(:,2,iglo) = & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,2,iglo)*spec(is)%u0 * jpar(:,it,ik) ) end do call copy(g, gnew) ! normalize call get_init_field (phinew, aparnew, bparnew) do i=1, 3 it = ittt(i) ik = ikkk(i) if (abs(phiamp(i)) > epsilon(0.0)) then ! Should this include abs? fac = phiamp(i) / phinew(0,it,ik) phi(:,it,ik) = phi(:,it,ik) * fac end if if (abs(aparamp(i)) > epsilon(0.0)) then fac = aparamp(i) / aparnew(0,it,ik) jpar(:,it,ik) = jpar(:,it,ik) * fac end if end do ! redefine g if (reality) then do it = 1, nakx/2 phi(:,it+(nakx+1)/2,1) = conjg(phi(:,(nakx+1)/2+1-it,1)) jpar(:,it+(nakx+1)/2,1) = conjg(jpar(:,(nakx+1)/2+1-it,1)) end do end if do iglo = g_lo%llim_proc, g_lo%ulim_proc it = it_idx(g_lo,iglo) ik = ik_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,1,iglo)*spec(is)%u0 * jpar(:,it,ik) ) g(:,2,iglo) = & ( dfac*spec(is)%dens0 * phi(:,it,ik) & + 2.*ufac* vpa(:,2,iglo)*spec(is)%u0 * jpar(:,it,ik) ) end do call copy(g, gnew) end subroutine ginit_ot !> FIXME : Add documentation subroutine ginit_kpar use species, only: spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0, theta0 use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use constants, only: zi use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd real, dimension (-ntgrid:ntgrid) :: dfac, ufac, tparfac, tperpfac integer :: iglo integer :: ig, ik, it, il, is call zero_array(phi) ; call zero_array(odd) if (width0 > 0.) then do ig = -ntgrid, ntgrid phi(ig,:,:) = exp(-((theta(ig)-theta0(:,:))/width0)**2)*cmplx(refac, imfac) end do else do ig = -ntgrid, ntgrid phi(ig,:,:) = cmplx(refac, imfac) end do end if if (chop_side) then if (left) then phi(:-1,:,:) = 0.0 else phi(1:,:,:) = 0.0 end if end if odd = zi * phi dfac = den0 + den1 * cos(theta) + den2 * cos(2.*theta) ufac = upar0 + upar1* sin(theta) + upar2* sin(2.*theta) tparfac = tpar0 + tpar1* cos(theta) + tpar2* cos(2.*theta) tperpfac = tperp0 + tperp1*cos(theta) + tperp2*cos(2.*theta) ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phiinit* &!spec(is)%z* & ( dfac * phi(:,it,ik) & + 2.*ufac* vpa(:,1,iglo) * odd(:,it,ik) & + tparfac*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = phiinit* &!spec(is)%z* & ( dfac * phi(:,it,ik) & + 2.*ufac* vpa(:,2,iglo) * odd(:,it,ik) & + tparfac*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) & +tperpfac*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,2,iglo) = 0.0 end do if (has_electron_species(spec)) then call flae (g, gnew) g = g - gnew end if call copy(g, gnew) end subroutine ginit_kpar !> FIXME : Add documentation subroutine ginit_gs use species, only: spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, ntheta0 use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx use constants, only: pi, zi use ran, only: ranf use array_utils, only: copy, zero_array implicit none complex, dimension (-ntgrid:ntgrid,ntheta0,naky) :: phi, odd integer :: iglo integer :: ik, it, il, is real :: phase call zero_array(phi) ; call zero_array(odd) do ik=1,naky do it=1,ntheta0 phase = 2.*pi*ranf() phi(:,it,ik) = cos(theta+phase)*cmplx(refac,imfac) odd(:,it,ik) = sin(theta+phase)*cmplx(refac,imfac) * zi end do end do !<DD>Should this be triggered for kt_grids::reality=.true. only? ! reality condition for k_theta = 0 component: do it = 1, ntheta0/2 phi(:,it+(ntheta0+1)/2,1) = conjg(phi(:,(ntheta0+1)/2+1-it,1)) odd(:,it+(ntheta0+1)/2,1) = conjg(odd(:,(ntheta0+1)/2+1-it,1)) enddo ! charge dependence keeps initial Phi from being too small do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = phiinit* &!spec(is)%z* & ( den1 * phi(:,it,ik) & + 2.*upar1* vpa(:,1,iglo) * odd(:,it,ik) & + tpar1*(vpa(:,1,iglo)**2-0.5) * phi(:,it,ik) & + tperp1*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = phiinit* &!spec(is)%z* & ( den1 * phi(:,it,ik) & + 2.*upar1* vpa(:,2,iglo) * odd(:,it,ik) & + tpar1*(vpa(:,2,iglo)**2-0.5) * phi(:,it,ik) & + tperp1*(vperp2(:,iglo)-1.) * phi(:,it,ik)) where (forbid(:,il)) g(:,2,iglo) = 0.0 end do if (has_electron_species(spec)) then call flae (g, gnew) g = g - gnew end if call copy(g, gnew) end subroutine ginit_gs !> FIXME : Add documentation subroutine ginit_xi use theta_grid, only: ntgrid, theta, bmag use le_grids, only: forbid, al use kt_grids, only: theta0 use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx use array_utils, only: copy implicit none integer :: iglo integer :: ik, it, il real, dimension(-ntgrid:ntgrid) :: xi do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) xi = sqrt(max(1.0-bmag*al(il),0.0)) g(:,1,iglo) = xi*exp(-((theta-theta0(it,ik))/width0)**2) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = -g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_xi !> FIXME : Add documentation subroutine ginit_xi2 use theta_grid, only: ntgrid, theta, bmag use le_grids, only: forbid, al use kt_grids, only: theta0 use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx use array_utils, only: copy implicit none integer :: iglo integer :: ik, it, il real, dimension(-ntgrid:ntgrid) :: xi do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) it = it_idx(g_lo,iglo) il = il_idx(g_lo,iglo) xi = sqrt(max(1.0-bmag*al(il),0.0)) g(:,1,iglo) = (1.0 - 3.0*xi*xi)*exp(-((theta-theta0(it,ik))/width0)**2) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_xi2 !> FIXME : Add documentation subroutine ginit_rh use le_grids, only: forbid, energy use dist_fn_arrays, only: g, gnew use gs2_layouts, only: g_lo, il_idx, ie_idx, is_idx use array_utils, only: copy implicit none integer :: iglo integer :: il, ie, is do iglo = g_lo%llim_proc, g_lo%ulim_proc il = il_idx(g_lo,iglo) ie = ie_idx(g_lo,iglo) is = is_idx(g_lo,iglo) g(:,1,iglo) = exp(-energy(ie,is)) where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_rh !> FIXME : Add documentation subroutine ginit_alf use theta_grid, only: theta use le_grids, only: forbid use dist_fn_arrays, only: g, gnew, vpa use gs2_layouts, only: g_lo, il_idx, is_idx use species, only: spec, is_electron_species use array_utils, only: copy implicit none integer :: iglo integer :: il, is do iglo = g_lo%llim_proc, g_lo%ulim_proc is = is_idx(g_lo,iglo) if (is_electron_species(spec(is))) cycle il = il_idx(g_lo,iglo) g(:,1,iglo) = sin(theta)*vpa(:,1,iglo)*spec(is)%z where (forbid(:,il)) g(:,1,iglo) = 0.0 g(:,2,iglo) = sin(theta)*vpa(:,2,iglo)*spec(is)%z where (forbid(:,il)) g(:,2,iglo) = 0.0 end do g = phiinit * g call copy(g, gnew) end subroutine ginit_alf !> FIXME : Add documentation subroutine ginit_zero use dist_fn_arrays, only: g, gnew use array_utils, only: zero_array implicit none call zero_array(g) call zero_array(gnew) end subroutine ginit_zero !> FIXME : Add documentation subroutine ginit_test3 use dist_fn_arrays, only: g, gnew, vpa use theta_grid, only: ntgrid, delthet, bmag use kt_grids, only: akx use theta_grid_params, only: eps, epsl, pk use gs2_layouts, only: g_lo, ik_idx, il_idx use mp, only: broadcast use constants, only: zi use array_utils, only: copy implicit none integer :: iglo, ik, il real :: c1, c2 call broadcast (epsl) call broadcast (eps) call broadcast (pk) do iglo = g_lo%llim_proc, g_lo%ulim_proc ik = ik_idx(g_lo,iglo) il = il_idx(g_lo,iglo) if (any(vpa(-ntgrid:ntgrid-1,:,iglo) == 0.0)) then c1 = 0.0 c2 = 0.0 else c2 = -akx(ik)*epsl/eps/pk & *sum(delthet(-ntgrid:ntgrid-1)/bmag(-ntgrid:ntgrid-1)) c1 = c2/sum(delthet(-ntgrid:ntgrid-1)/vpa(-ntgrid:ntgrid-1,1,iglo)) c2 = c2/sum(delthet(-ntgrid:ntgrid-1)/vpa(-ntgrid:ntgrid-1,2,iglo)) end if g(:,1,iglo) = -zi*akx(ik)*epsl/eps/pk*vpa(:,1,iglo)/bmag - zi*c1 g(:,2,iglo) = -zi*akx(ik)*epsl/eps/pk*vpa(:,2,iglo)/bmag - zi*c2 end do call copy(g, gnew) end subroutine ginit_test3 !> FIXME : Add documentation subroutine ginit_convect use dist_fn_arrays, only: g, gnew use theta_grid, only: theta use kt_grids, only: theta0 use gs2_layouts, only: g_lo, it_idx, ik_idx use run_parameters, only: k0 use array_utils, only: copy implicit none integer :: it, ik, iglo do iglo = g_lo%llim_proc, g_lo%ulim_proc it = it_idx(g_lo,iglo) ik = ik_idx(g_lo,iglo) g(:,1,iglo) = exp(cmplx(-(theta-theta0(it,ik))**2,k0*theta)) g(:,2,iglo) = g(:,1,iglo) end do call copy(g, gnew) end subroutine ginit_convect !> FIXME : Add documentation subroutine ginit_recon3 use mp, only: proc0, mp_abort use species, only: nspec, spec, has_electron_species use theta_grid, only: ntgrid, theta use kt_grids, only: naky, nakx => ntheta0, akx, aky, reality use kt_grids, only: nx,ny,kperp2 use dist_fn_arrays, only: g, gnew, vpa, vperp2 use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx use gs2_transforms, only: inverse2,transform2 use run_parameters, only: beta use le_grids, only: integrate_moment use diagnostics_moments, only: getmoms_notgc, init_mom_coeff use dist_fn, only: gamtot, gamtot1, gamtot2, get_init_field use fields_arrays, only: phinew, bparnew use constants, only: pi, zi use file_utils, only: error_unit, open_output_file, close_output_file use ran, only: ranf use array_utils, only: copy implicit none real, dimension(:, :, :, :), allocatable :: mom_coeff real, dimension(:, :, :), allocatable :: mom_coeff_npara, mom_coeff_nperp real, dimension(:, :, :), allocatable :: mom_coeff_tpara, mom_coeff_tperp real, dimension(:, :, :), allocatable :: mom_shift_para, mom_shift_perp integer :: iglo integer :: ig, ik, it, is integer :: i,j ! nkxyz and ukxyz determine profiles in kx-ky plane ! nkxyz is for N, T, and ukxyz is for U ! [do not control amplitude by these variables] complex :: nkxyz(-ntgrid:ntgrid,nakx,naky) complex :: ukxyz(-ntgrid:ntgrid,nakx,naky) complex :: nkxyz_eq(-ntgrid:ntgrid,nakx,naky) complex :: ukxyz_eq(-ntgrid:ntgrid,nakx,naky) complex, allocatable :: g_eq(:,:,:) ! equilibrium and perturbation complex :: nkxy_eq(nakx,naky), ukxy_eq(nakx,naky) ! *fac determine z profile ! [do not control amplitude by these variables] real :: dfac(-ntgrid:ntgrid), ufac(-ntgrid:ntgrid) real :: tparfac(-ntgrid:ntgrid), tperpfac(-ntgrid:ntgrid) real :: ct(-ntgrid:ntgrid), st(-ntgrid:ntgrid) real :: c2t(-ntgrid:ntgrid), s2t(-ntgrid:ntgrid) ! normalization real :: fac real, allocatable :: nnrm(:,:,:),unrm(:,:,:) real, allocatable :: tparanrm(:,:,:),tperpnrm(:,:,:) real :: check(3), current character (len=2) :: cfit complex, allocatable :: dens(:,:,:,:),upar(:,:,:,:) complex, allocatable :: tpar(:,:,:,:),tper(:,:,:,:) complex, allocatable :: phi(:,:,:),apar(:,:,:),bpar(:,:,:) complex, allocatable :: jpar(:,:,:) real :: save_dens0, save_tperp0, save_u0, ratio ! real space profile to be Fourier transformed real :: xx,dx,lx,ly integer, parameter :: nfxp=2**10 integer :: nfx real, allocatable :: nxy(:,:),uxy(:,:) real, allocatable :: phixy(:,:,:),aparxy(:,:,:),bparxy(:,:,:) real, allocatable :: jparxy(:,:,:) real, allocatable :: densxy(:,:,:,:),uparxy(:,:,:,:) real, allocatable :: tparxy(:,:,:,:),tperxy(:,:,:,:) complex, allocatable :: by(:,:,:) real, allocatable :: byxy(:,:,:) integer :: nnx,nny integer :: unit real :: xmax, bymax, xa, xl1, xl2 real :: fmin, fmax, a, b if(debug.and.proc0) write(6,*) 'Initialization recon3' cfit = 'A0' check = 0 current = 0 nfx = nfxp if (nfxp < nx) nfx=nx !!! adjust input parameters to kill initial field if wanted if(debug.and.proc0) write(6,*) 'check parameters' do is=1,nspec if(adj_spec == is) cycle check(1)=check(1)+spec(is)%z*spec(is)%dens*spec(is)%dens0 check(2)=check(2)+spec(is)%dens*spec(is)%temp* & & (spec(is)%dens0+spec(is)%tperp0) check(3)=check(3)+spec(is)%z*spec(is)%dens*spec(is)%stm*spec(is)%u0 end do if(adj_spec == 0) then ! just warn if fields don't satisfy given conditions if(null_phi.and.null_bpar) then if(sum(check(1:2)) /= 0.) then if(proc0) write(6,'(a)') & 'WARNING: finite Phi or Bpar in initial condition' endif else if(null_bpar.and..not.null_phi) then ratio=gamtot1(0,eq_mode_n+1,1)/gamtot(0,eq_mode_n+1,1) if(check(1)/check(2) /= ratio) then if(proc0) write(6,'(a)') & 'WARNING: finite Bpar in initial condition' endif else if(null_phi.and..not.null_bpar) then ratio=-(2.*gamtot2(0,eq_mode_n+1,1)+2./beta) & & /gamtot1(0,eq_mode_n+1,1) if(check(1)/check(2) /= ratio) then if(proc0) write(6,'(a)') & 'WARNING: finite Bpar in initial condition' endif endif if(null_apar) then if(check(3) /= 0.) then if(proc0) write(6,'(a)') & 'WARNING: finite Apar in initial condition' endif endif else ! adjust input parameter to satisfy given conditions if(null_phi.and.null_bpar) then save_dens0=spec(adj_spec)%dens0 save_tperp0=spec(adj_spec)%tperp0 spec(adj_spec)%dens0=-check(1)/(spec(adj_spec)%z*spec(adj_spec)%dens) spec(adj_spec)%tperp0=-spec(adj_spec)%dens0 & & -check(2)/(spec(adj_spec)%dens*spec(adj_spec)%temp) if(spec(adj_spec)%dens0 /= save_dens0) then if(proc0) write(6,'(a,i0,a,f10.2,a)') & & 'WARNING: Initial density of spec=', spec(adj_spec)%type, & & ' is adjusted to ', spec(adj_spec)%dens0, & & ' to kill Phi and Bpar' endif if(spec(adj_spec)%tperp0 /= save_tperp0) then if(proc0) write(6,'(a,i0,a,f10.2,a)') & & 'WARNING: Initial Tperp of spec=', spec(adj_spec)%type, & & ' is adjusted to ', spec(adj_spec)%tperp0, & & ' to kill Phi and Bpar' endif else if(null_bpar.and..not.null_phi.and.eq_type.eq.'sinusoidal') then save_tperp0=spec(adj_spec)%tperp0 check(1)=check(1)+ & & spec(adj_spec)%z*spec(adj_spec)%dens*spec(adj_spec)%dens0 ratio=gamtot1(0,eq_mode_n+1,1)/gamtot(0,eq_mode_n+1,1) spec(adj_spec)%tperp0=(-check(1)*ratio-check(2)) & & /(spec(adj_spec)%dens*spec(adj_spec)%temp) & & -spec(adj_spec)%dens0 if(spec(adj_spec)%tperp0 /= save_tperp0) then if(proc0) write(6,'(a,i0,a,f10.2,a)') & & 'WARNING: Initial Tperp of spec=', spec(adj_spec)%type, & & ' is adjusted to ', spec(adj_spec)%tperp0, & & ' to kill Bpar' endif else if(null_phi.and..not.null_bpar.and.eq_type.eq.'sinusoidal') then save_tperp0=spec(adj_spec)%tperp0 check(1)=check(1)+ & & spec(adj_spec)%z*spec(adj_spec)%dens*spec(adj_spec)%dens0 ratio=-(2.*gamtot2(0,eq_mode_n+1,1)+2./beta) & & /gamtot1(0,eq_mode_n+1,1) spec(adj_spec)%tperp0=(-check(1)*ratio-check(2)) & & /(spec(adj_spec)%dens*spec(adj_spec)%temp) & & -spec(adj_spec)%dens0 if(spec(adj_spec)%tperp0 /= save_tperp0) then if(proc0) write(6,'(a,i0,a,f10.2,a)') & & 'WARNING: Initial Tperp of spec=', spec(adj_spec)%type, & & ' is adjusted to ', spec(adj_spec)%tperp0, & & ' to kill Phi' endif endif if (null_apar) then save_u0=spec(adj_spec)%u0 spec(adj_spec)%u0=-check(3)/ & & (spec(adj_spec)%z*spec(adj_spec)%dens*spec(adj_spec)%stm) if(spec(adj_spec)%u0 /= save_u0) then if(proc0) write(6,'(a,i0,a,f10.2,a)') & & 'WARNING: Initial U of spec=', spec(adj_spec)%type, & & ' is adjusted to ', spec(adj_spec)%u0, & & ' to kill Apar' endif endif endif !!! initialize if(debug.and.proc0) write(6,*) 'initialize variable' nkxyz(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.) ukxyz(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.) nkxyz_eq(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.) ukxyz_eq(-ntgrid:ntgrid,1:nakx,1:naky)=cmplx(0.,0.) nkxy_eq(1:nakx,1:naky)=cmplx(0.,0.) ukxy_eq(1:nakx,1:naky)=cmplx(0.,0.) dfac(-ntgrid:ntgrid)=1. ufac(-ntgrid:ntgrid)=1. tparfac(-ntgrid:ntgrid)=1. tperpfac(-ntgrid:ntgrid)=1. !!! equilibrium if(phiinit0 /= 0.) then if(nakx==1 .and. naky==1) then ! grid_option = single case nkxy_eq(1,1)=cmplx(.5,0.) ukxy_eq(1,1)=cmplx(.5,0.) else if(debug.and.proc0) write(6,*) 'set equilibrium profile' allocate(nxy(nfx,ny),uxy(nfx,ny)) ! lx=2.*pi*x0; ly=2.*pi*y0 lx=2.*pi/(akx(2)-akx(1)); ly=2.*pi/(aky(2)-aky(1)) dx=lx/nfx ! if width is negative, it gives the ratio to the box size if(prof_width < 0. ) prof_width=-prof_width*lx select case (eq_type) case ('sinusoidal') ! this gives n,Tpara,Tperp \propto cos^2 (2 pi/Lx) ! nxy(eq_mode1(1),eq_mode1(2))=cmplx(.25, 0.) ! this gives Apara \propto cos(2 pi/Lx), By \propto sin(2 pi x/Lx) ! uxy(eq_mode2(1),eq_mode2(2))=cmplx(.5, 0.) do it=1,nfx xx=dx*(it-1) nxy(it,1:ny)=0. uxy(it,1:ny)=cos(2.*pi/lx*xx*eq_mode_u) end do case ('porcelli') do it=1,nfx xx=dx*(it-1) nxy(it,1:ny)=0. uxy(it,1:ny)=1./cosh((xx-.5*lx)/prof_width)**2 & & * (tanh(2.*pi/lx*xx)**2+tanh(2.*pi/lx*(xx-lx))**2 & & - tanh(2.*pi)**2) / (2.*tanh(pi)**2-tanh(2.*pi)**2) end do case ('doubleharris') do it=1,nfx fmin=tanh(.75*lx/prof_width)-tanh(.25*lx/prof_width) fmax=2.*tanh(.25*lx/prof_width) xx=dx*(it-1) nxy(it,1:ny)=0. uxy(it,1:ny)= 2.*( & & tanh((xx-.25*lx)/prof_width) - & & tanh((xx-.75*lx)/prof_width) - fmin)/(fmax-fmin)-1. end do case default if(proc0) write(error_unit(),'(2a)') & & 'ERROR: Invalid equilibrium type', eq_type call mp_abort('Invalid equilibrium type') end select ! subtract (0,0) mode ! since it (constant part of potential) does nothing do ik=1,ny uxy(1:nfx,ik)=uxy(1:nfx,ik) & - sum(uxy(1:nfx,ik))/nfx end do call inverse2(nxy, nkxy_eq, ny, nfx) call inverse2(uxy, ukxy_eq, ny, nfx) deallocate(nxy,uxy) endif endif do ig = -ntgrid, ntgrid nkxyz(ig,1:nakx,1:naky) = phiinit0*nkxy_eq(1:nakx,1:naky) ukxyz(ig,1:nakx,1:naky) = phiinit0*ukxy_eq(1:nakx,1:naky) end do if(.not.(nakx==1 .and. naky==1</