!> FIXME : Add documentation !! !! @note nexp has been changed in a rush. Only known to be correct for nexp=2 module hyper use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN implicit none private public :: init_hyper, finish_hyper, hyper_diff, D_res, D_v, D_eta, nexp public :: read_parameters, wnml_hyper, check_hyper, hypervisc_filter public :: hyper_config_type, set_hyper_config, get_hyper_config ! D_v, D_eta are hyper coefficients, normalized correctly ! i.e., either by unity or by 1/k_perp**2*nexp real :: D_v, D_eta, D_hypervisc, D_hyperres, omega_osc, D_hyper real :: akx4_max, aky4_max, aky_max, akperp4_max integer :: hyper_option_switch, nexp integer, parameter :: hyper_option_none = 1, & hyper_option_visc = 2, & hyper_option_res = 3, & hyper_option_both = 4, & hyper_option_simple3D = 5 character(9) :: hyper_option logical :: const_amp, isotropic_shear, damp_zonal_only, gridnorm real, dimension (:, :), allocatable :: D_res real, dimension (:, :, :), allocatable :: hypervisc_filter ! Coefficients for simple3D hyperviscous model ! the maximum k_perp ^ 2 including geometric coefficients for non-zonal and zonal modes real :: kperp2_max, kperp2_max_zonal ! Hyperviscous coefficient, hyperviscous power, ky cut-off, kx cut-off real :: D_hyper3D, P_hyper3D, ky_cut, kx_cut ! Controls if zonal modes are treated identically to non-zonal modes logical :: isotropic_model logical :: initialized = .false., hyper_on = .false. !> Used to represent the input configuration of hyper type, extends(abstract_config_type) :: hyper_config_type ! namelist : hyper_knobs ! indexed : false !> Determines whether hyperviscosity includes time dependent amplitude !> factor when calculating damping rate. Recommend `true` for !> linear runs and `false` for nolinear runs, since amplutide of !> turbulence grows linearly with time in linear run. logical :: const_amp = .false. !> If `hyper_option = 'both'` is used then this sets both the !> hyperresistivity and hyperviscosity damping coefficients. Can !> override the individual coefficients with !> [[hyper_knobs:D_hyperres]] and [[hyper_knobs:D_hypervisc]]. real :: d_hyper = -10.0 !> Used with the simple3D hyperviscosity model of the form !> D_hyper3D * (|kperp|/ max |kperp|)^P_hyper3D !> and ky_cut, kx_cut set max |kperp| real :: d_hyper3d = -10. !> Sets hyperresistivity parameter multiplying damping term. real :: d_hyperres = -10.0 !> Sets hyperviscosity parameter multiplying damping term. See !> [E. Belli (2006) !> thesis](https://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=50BFD54A8F8D72FC225D025FDEDFFFA5?doi=10.1.1.706.9568&rep=rep1&type=pdf) !> for more information. real :: d_hypervisc = -10.0 !> If `true` then hyperdissipation only applied to the zonal !> mode. logical :: damp_zonal_only = .false. !> If `true` (default) then set wavenumber parameters entering !> the models based on the maximum `ky` and `kx` included in the !> current simulation. If `false` then these values are set to 1. logical :: gridnorm = .true. !> Selects the type of hyper terms included. Should be one of !> !> - 'default' -- no hyper terms included. !> - 'none' -- the same as default. !> - 'visc_only' -- only hyperviscosity included. !> - 'res_only' -- only hyperresistivity included. !> - 'both' -- both viscosity and resistivity included. !> - 'simple3D' -- simple hyperviscous dissipation rate of the !> form D_hyper3D * (|kperp|/ max |kperp|)^P_hyper3D !> described in [“Multiscale turbulence in magnetic !> confinement fusion devices”, M. Hardman, DPhil Thesis, !> appendix !> B.4](https://ora.ox.ac.uk/objects/uuid:233a22cb-3c8b-4fe0-a689-4a37d8fe0314) !> note that the dissipation in this version is applied to g, not h, as !> in the reference. character(len = 9) :: hyper_option = 'default' !> if true damp zonal and drift waves with same dissipation formula logical :: isotropic_model = .true. !> If `true` then use isotropic shear model. logical :: isotropic_shear = .true. !> Used with the simple3D hyperviscosity model of the form !> D_hyper3D * (|kperp|/ max |kperp|)^P_hyper3D !> and ky_cut, kx_cut set max |kperp| real :: kx_cut = -10. !> Used with the simple3D hyperviscosity model of the form !> D_hyper3D * (|kperp|/ max |kperp|)^P_hyper3D !> and ky_cut, kx_cut set max |kperp| real :: ky_cut = -10. !> Sets the power to which \(k_\bot^2\) is raised in the dissipation filter. integer :: nexp = 2 !> Sets a parameter in the anisotropic shearing rate calculation. real :: omega_osc = 0.4 !> Used with the simple3D hyperviscosity model of the form !> D_hyper3D * (|kperp|/ max |kperp|)^P_hyper3D !> and ky_cut, kx_cut set max |kperp| real :: p_hyper3d = 4. contains procedure, public :: read => read_hyper_config procedure, public :: write => write_hyper_config procedure, public :: reset => reset_hyper_config procedure, public :: broadcast => broadcast_hyper_config procedure, public, nopass :: get_default_name => get_default_name_hyper_config procedure, public, nopass :: get_default_requires_index => get_default_requires_index_hyper_config end type hyper_config_type type(hyper_config_type) :: hyper_config contains !> FIXME : Add documentation subroutine check_hyper(report_unit) implicit none integer, intent(in) :: report_unit select case (hyper_option_switch) case (hyper_option_none) if (D_hyperres > 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses no hyperresistivity. & &D_hyperres ignored.')") trim(hyper_option) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) D_hyperres = -10. end if if (D_hypervisc > 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses no hyperviscosity. & &D_hypervisc ignored.')") trim(hyper_option) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) D_hypervisc = -10. endif case (hyper_option_visc) hyper_on = .true. if (D_hyperres > 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses no hyperresistivity. & &D_hyperres ignored.')") trim(hyper_option) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) D_hyperres = -10. end if if (D_hypervisc < 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses hyperviscosity but & &D_hypervisc < 0.')") trim(hyper_option) write (report_unit, fmt="('No hyperviscosity used.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) hyper_on = .false. endif case (hyper_option_res) hyper_on = .true. if (D_hyperres < 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses hyperresistivity but D_hyperres < 0.')") trim(hyper_option) write(report_unit, fmt="('No hyperresistivity used.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) hyper_on = .false. end if if (D_hypervisc > 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses no hyperviscosity. D_hypervisc ignored.')") trim(hyper_option) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) D_hypervisc = -10. endif case (hyper_option_both) hyper_on = .true. if (D_hyperres < 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses hyperresistivity but D_hyperres < 0.')") trim(hyper_option) write (report_unit, fmt="('No hyperresistivity used.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if if (D_hypervisc < 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('hyper_option = ',a,' chooses hyperviscosity but D_hypervisc < 0.')") trim(hyper_option) write (report_unit, fmt="('No hyperviscosity used.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) endif if (D_hypervisc < 0. .and. D_hyperres < 0.) hyper_on = .false. end select if (hyper_on) then write (report_unit, *) write (report_unit, fmt="('------------------------------------------------------------')") write (report_unit, *) select case (hyper_option_switch) case (hyper_option_visc) write (report_unit, *) write (report_unit, fmt="('Hyperviscosity included without hyperresistivity.')") if (const_amp) then write (report_unit, fmt="('Damping rate is ',e11.4,' at highest k_perp.')") D_hypervisc else write (report_unit, fmt="('The damping coefficent is ',e11.4)") D_hypervisc write (report_unit, fmt="('The damping rate is proportional to the RMS amplitude of the turbulence.')") end if if (isotropic_shear) then write (report_unit, fmt="('The hyperviscosity is isotropic in the perpendicular plane.')") write (report_unit, fmt="('This is appropriate for MHD-like calculations.')") else write (report_unit, fmt="('The hyperviscosity is anisotropic in the perpendicular plane.')") write (report_unit, fmt="('This is appropriate for drift-type calculations.')") write (report_unit, fmt="('omega_osc = ',e11.4)") omega_osc end if case (hyper_option_res) write (report_unit, *) write (report_unit, fmt="('Hyperresistivity included without hyperviscosity.')") if (const_amp) then write (report_unit, fmt="('Damping rate is ',e11.4,' at highest k_perp.')") D_hyperres else write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('const_amp = .false. is not implemented for hyperresistivity.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if if (isotropic_shear) then write (report_unit, fmt="('The hyperresistivity is isotropic in the perpendicular plane.')") write (report_unit, fmt="('This is appropriate for MHD-like calculations.')") else write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('isotropic_shear = .false. is not implemented for hyperresistivity.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if case (hyper_option_both) write (report_unit, *) write (report_unit, fmt="('Hyperresistivity and hyperviscosity included.')") if (const_amp) then write (report_unit, fmt="('Damping rate is ',e11.4,' at highest k_perp.')") D_hyperres else write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('const_amp = .false. is not implemented for hyperresistivity.')") write (report_unit, fmt="('THIS IS AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if if (isotropic_shear) then write (report_unit, fmt="('The damping is isotropic in the perpendicular plane.')") write (report_unit, fmt="('This is appropriate for MHD-like calculations.')") else write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('isotropic_shear = .false. is not implemented for hyperresistivity.')") write (report_unit, fmt="('THIS IS AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, *) end if end select end if end subroutine check_hyper !> FIXME : Add documentation subroutine wnml_hyper(unit) implicit none integer, intent(in) :: unit if (.not. hyper_on) return call hyper_config%write(unit) end subroutine wnml_hyper !> FIXME : Add documentation subroutine init_hyper(hyper_config_in) use kt_grids, only: ntheta0, naky, akx, aky, kperp2 use theta_grid, only: gds2, gds21, gds22, shat use gs2_time, only: code_dt use gs2_layouts, only: init_gs2_layouts use warning_helpers, only: is_zero implicit none type(hyper_config_type), intent(in), optional :: hyper_config_in integer :: ik, it if (initialized) return initialized = .true. call init_gs2_layouts call read_parameters(hyper_config_in) call allocate_arrays select case (hyper_option_switch) case(hyper_option_both, hyper_option_res, hyper_option_visc) ! Define variables used in hyperviscosity and hyperresistivity models if (gridnorm) then akx4_max = akx(ntheta0 / 2 + 1)**(2 * nexp) aky_max = aky(naky) aky4_max = aky(naky)**(2 * nexp) akperp4_max = (akx(ntheta0 / 2 + 1)**2 + aky(naky)**2)**nexp else akx4_max = 1. aky_max = 1. aky4_max = 1. akperp4_max = 1. end if ! Get D_res set up if needed if (D_hyperres > 0.) then do ik = 1, naky do it = 1, ntheta0 D_res(it, ik) = D_hyperres*code_dt * (aky(ik)**2 + akx(it)**2)**nexp & / akperp4_max end do end do D_eta = D_hyperres / akperp4_max end if if (D_hypervisc > 0.) D_v = D_hypervisc / akperp4_max case(hyper_option_simple3D) ! Make the default cut-offs the grid scale if (ky_cut <= 0.) ky_cut = aky(naky) if (kx_cut <= 0.) kx_cut = akx((ntheta0+1)/2) ! Here we're evaluating kperp2 for the specified ky, kx value if (is_zero(shat)) then kperp2_max = ky_cut*ky_cut + kx_cut*kx_cut else kperp2_max = maxval(ky_cut * ky_cut * gds2 + & 2.0 * ky_cut * kx_cut * gds21 / shat + & kx_cut * kx_cut * gds22 / (shat * shat)) end if if (isotropic_model .or. is_zero(shat)) then kperp2_max_zonal = kperp2_max else kperp2_max_zonal = maxval(kx_cut * kx_cut * gds22 / (shat * shat)) end if ! Calculate the filter ik =1 ! Zonal modes do it = 1, ntheta0 hypervisc_filter(:, it, ik) = exp(- ( D_hyper3D * code_dt & * ( kperp2(:, it, ik) / kperp2_max_zonal)**(P_hyper3D / 2))) end do do ik = 2, naky ! Non-zonal modes do it = 1, ntheta0 hypervisc_filter(:, it, ik) = exp(- ( D_hyper3D * code_dt & * ( kperp2(:, it, ik) / kperp2_max )**(P_hyper3D / 2)) ) end do end do end select end subroutine init_hyper !> FIXME : Add documentation subroutine read_parameters(hyper_config_in) use file_utils, only: error_unit use text_options, only: text_option, get_option_value use mp, only: proc0 use warning_helpers, only: not_exactly_equal implicit none type(hyper_config_type), intent(in), optional :: hyper_config_in type (text_option), dimension(6), parameter :: hyperopts = & (/ text_option('default', hyper_option_none), & text_option('none', hyper_option_none), & text_option('visc_only', hyper_option_visc), & text_option('res_only', hyper_option_res), & text_option('both', hyper_option_both), & text_option('simple3D', hyper_option_simple3D) /) integer :: ierr if (present(hyper_config_in)) hyper_config = hyper_config_in call hyper_config%init(name = 'hyper_knobs', requires_index = .false.) ! Copy out internal values into module level parameters associate(self => hyper_config) #include "hyper_copy_out_auto_gen.inc" end associate ierr = error_unit() call get_option_value & (hyper_option, hyperopts, hyper_option_switch, & ierr, "hyper_option in hyper_knobs",.true.) if (.not. isotropic_shear .and. nexp /= 2) then if (proc0) write (ierr, *) 'Forcing nexp = 2. Higher values not implemented for anisotropic shear model.' nexp = 2 end if select case (hyper_option_switch) case (hyper_option_none) if (D_hyperres > 0.) then if (proc0) write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses no hyperresistivity. D_hyperres ignored.' D_hyperres = -10. end if if (D_hypervisc > 0.) then if (proc0) write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses no hyperviscosity. D_hypervisc ignored.' D_hypervisc = -10. endif case (hyper_option_visc) hyper_on = .true. if (proc0) write (ierr, *) 'WARNING: It is inconsistent to set D_hypervisc different from ', & 'D_hyperres. Recommend: Set them equal.' if (D_hyperres > 0.) then if (proc0) write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses no hyperresistivity. D_hyperres ignored.' D_hyperres = -10. end if if (D_hypervisc < 0.) then if (proc0) then write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses hyperviscosity but D_hypervisc < 0 is illegal.' write(ierr, *) 'No hyperviscosity used.' end if hyper_on = .false. endif case (hyper_option_res) hyper_on = .true. if (proc0) write (ierr, *) 'WARNING: It is inconsistent to set D_hypervisc different from ', & 'D_hyperres. Recommend: Set them equal.' if (D_hyperres < 0.) then if (proc0) then write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses hyperresistivity but D_hyperres < 0 is illegal.' write(ierr, *) 'No hyperresistivity used.' end if hyper_on = .false. end if if (D_hypervisc > 0.) then if (proc0) write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses no hyperviscosity. D_hypervisc ignored.' D_hypervisc = -10. endif case (hyper_option_both) hyper_on = .true. if (D_hyper < 0.) then if (not_exactly_equal(D_hyperres, D_hypervisc)) then if (proc0) write (ierr, *) 'WARNING: It is inconsistent to set D_hypervisc different from ', & 'D_hyperres. Recommend: Set them equal.' end if else if (proc0) write (ierr, *) 'WARNING: Setting D_hypervisc = D_hyperres, each to value of D_hyper' D_hyperres = D_hyper D_hypervisc = D_hyper end if if (D_hyperres < 0.) then if (proc0) then write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses hyperresistivity but D_hyperres < 0 is illegal.' write(ierr, *) 'No hyperresistivity used.' end if end if if (D_hypervisc < 0.) then if (proc0) then write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses hyperviscosity but D_hypervisc < 0 is illegal.' write(ierr, *) 'No hyperviscosity used.' end if endif if (D_hypervisc < 0. .and. D_hyperres < 0.) hyper_on = .false. case (hyper_option_simple3D) ! MRH hyper_on = .true. if (P_hyper3D < 4.) then hyper_on = .false. if(proc0) write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses a simple hyperviscous filter but P_hyper3D < 4. is illegal.' if(proc0) write(ierr, *) 'No hyperviscosity used.' endif if (D_hyper3D < 0.) then hyper_on = .false. if(proc0) write(ierr, *) 'hyper_option = ',hyper_option, & ' chooses a simple hyperviscous filter but D_hyper3D < 0 is illegal.' if(proc0) write(ierr, *) 'No hyperviscosity used.' endif end select end subroutine read_parameters !> FIXME : Add documentation subroutine allocate_arrays use theta_grid, only: ntgrid use kt_grids, only: ntheta0, naky implicit none if (.not. allocated(D_res)) allocate(D_res(ntheta0, naky)) if (.not. allocated(hypervisc_filter)) & allocate(hypervisc_filter(-ntgrid:ntgrid,ntheta0,naky)) ! Initialise module level variables D_v = 0 ; D_eta = 0 ; D_res = 0 akx4_max = 0 ; aky_max = 0 ; aky4_max = 0 ; akperp4_max = 0 kperp2_max = 0 ; kperp2_max_zonal = 0 ; hypervisc_filter = 1 end subroutine allocate_arrays !> FIXME : Add documentation subroutine hyper_diff (g0, phi) use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx use theta_grid, only: ntgrid use gs2_time, only: code_dt use kt_grids, only: aky, akx, naky, ntheta0 implicit none complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: g0 complex, dimension (-ntgrid:,:,:), intent (in) :: phi real, dimension (-ntgrid:ntgrid) :: shear_rate_nz, shear_rate_z, shear_rate_z_nz integer :: iglo, ik, it !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Includes models by Belli and Hammett ! to calculate the x-y avged shearing rate ! S(theta)^2 = <|grad_perp|^4 |phi|^2> ! = sum_over_k(kperp^4 * |phi|^2) ! (times crazy fac factor due to FFT conventions.) ! ! and to implement this anisotropically in k_perp, taking into ! account properties of zonal flows. ! ! Begun December 2001 ! ! Number conservation added April 2006 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (.not. hyper_on) return ! Calculate hypervisc filter for non-simple3D model select case (hyper_option_switch) case(hyper_option_both, hyper_option_res, hyper_option_visc) if (D_hypervisc < 0.) return if(isotropic_shear) then call iso_shear else call aniso_shear end if end select ! Apply hypervisc filter !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(iglo, ik, it) & !$OMP SHARED(g_lo, damp_zonal_only, g0, hypervisc_filter) & !$OMP SCHEDULE(static) do iglo = g_lo%llim_proc, g_lo%ulim_alloc ik = ik_idx(g_lo, iglo) it = it_idx(g_lo, iglo) if (damp_zonal_only .and. ik /= 1) cycle g0(:, 1, iglo) = g0(:, 1, iglo) * hypervisc_filter(:, it, ik) g0(:, 2, iglo) = g0(:, 2, iglo) * hypervisc_filter(:, it, ik) end do !$OMP END PARALLEL DO contains !> FIXME : Add documentation subroutine aniso_shear real, dimension(naky) :: fac ! Do the Belli-Hammett anisotropic calculation ! which accounts for some zonal/non-zonal differences fac = 0.5 fac(1) = 1.0 ! shearing rate due to non-zonal modes (on nonzonal modes) shear_rate_nz = 0. ! shearing rate due to zonal modes (on nonzonal modes) shear_rate_z = 0 ! shearing rate due to nonzonal modes (on zonal modes) shear_rate_z_nz = 0. ik = 1 !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(it) & !$OMP SHARED(ik, ntheta0, phi, akx, aky, fac) & !$OMP REDUCTION(+: shear_rate_z) & !$OMP SCHEDULE(static) do it = 1, ntheta0 shear_rate_z = shear_rate_z + real(conjg(phi(:, it, ik)) & * phi(:, it, ik)) * (akx(it)**2 + aky(ik)**2)**2 * fac(ik) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(ik, it) & !$OMP SHARED(naky, ntheta0, phi, akx, aky, fac) & !$OMP REDUCTION(+: shear_rate_nz, shear_rate_z_nz) & !$OMP SCHEDULE(static) do ik = 2, naky do it = 1, ntheta0 shear_rate_nz = shear_rate_nz + real(conjg(phi(:, it, ik)) & * phi(:, it, ik)) * (akx(it)**2 + aky(ik)**2)**2 * fac(ik) shear_rate_z_nz = shear_rate_z_nz + real(conjg(phi(:, it, ik)) & * phi(:, it, ik)) * aky(ik)**4 * fac(ik) end do end do !$OMP END PARALLEL DO ! shear_rate = shear_rate ** 0.5 shear_rate_nz = 0.5 * ( -omega_osc + sqrt(omega_osc ** 2 + 2 * shear_rate_nz)) shear_rate_z_nz = 0.5 * ( -omega_osc + sqrt(omega_osc ** 2 + 2 * shear_rate_z_nz)) shear_rate_z = 0.5 * ( -omega_osc + sqrt(omega_osc ** 2 + 2 * shear_rate_z)) ! end of anisotropic shearing rate calculations ik = 1 !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(it) & !$OMP SHARED(ik, ntheta0, hypervisc_filter, D_hypervisc, code_dt, akx, & !$OMP shear_rate_z_nz, akx4_max) & !$OMP SCHEDULE(static) do it = 1, ntheta0 hypervisc_filter(:, it, ik) = exp(- (D_hypervisc * code_dt & * (shear_rate_z_nz * akx(it) ** 4 / akx4_max ))) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO DEFAULT(none) & !$OMP PRIVATE(ik, it) & !$OMP SHARED(naky, ntheta0, hypervisc_filter, D_hypervisc, code_dt, akx, aky, & !$OMP shear_rate_nz, shear_rate_z, nexp, akperp4_max, akx4_max, aky_max) & !$OMP SCHEDULE(static) do ik = 2, naky do it = 1, ntheta0 hypervisc_filter(:, it, ik) = exp(- ( D_hypervisc * code_dt & * ( shear_rate_nz * (aky(ik) ** 2 + akx(it) ** 2 )**nexp / akperp4_max & + shear_rate_z * (akx(it) ** 4) * aky(ik) / (akx4_max * aky_max)))) end do end do !$OMP END PARALLEL DO end subroutine aniso_shear !> FIXME : Add documentation subroutine iso_shear real, dimension(naky) :: fac if (const_amp) then shear_rate_nz = 1. else fac = 0.5 fac(1) = 1.0 shear_rate_nz = 0. do ik = 1, naky if (damp_zonal_only .and. ik > 1) exit do it = 1, ntheta0 shear_rate_nz = shear_rate_nz & + real(conjg(phi(:, it, ik)) * phi(:, it, ik)) & * (akx(it)**2 + aky(ik)**2)**2 * fac(ik) end do end do shear_rate_nz = sqrt(shear_rate_nz) end if do ik = 1, naky !Rely on hypervisc_filter initialised to 1 if (damp_zonal_only .and. ik > 1) exit do it = 1, ntheta0 hypervisc_filter(:, it, ik) = exp(- ( D_hypervisc * code_dt & * ( shear_rate_nz * (aky(ik) ** 2 + akx(it) ** 2 )**nexp / akperp4_max))) end do end do end subroutine iso_shear end subroutine hyper_diff !> FIXME : Add documentation subroutine finish_hyper implicit none hyper_on = .false. initialized = .false. if (allocated(D_res)) deallocate(D_res) if (allocated(hypervisc_filter)) deallocate(hypervisc_filter) call hyper_config%reset() end subroutine finish_hyper #include "hyper_auto_gen.inc" end module hyper