diff --git a/src/model_mod.F90 b/src/model_mod.F90 index 12d16a5..1f1406f 100644 --- a/src/model_mod.F90 +++ b/src/model_mod.F90 @@ -1,141 +1,145 @@ MODULE model ! Module for diagnostic parameters USE prec_const IMPLICIT NONE PRIVATE INTEGER, PUBLIC, PROTECTED :: CLOS = 0 ! linear truncation method INTEGER, PUBLIC, PROTECTED :: NL_CLOS = 0 ! nonlinear truncation method INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model CHARACTER(len=16), & PUBLIC, PROTECTED ::LINEARITY= 'linear' ! To turn on non linear bracket term LOGICAL, PUBLIC, PROTECTED :: KIN_E = .true. ! Simulate kinetic electron (adiabatic otherwise) REAL(dp), PUBLIC, PROTECTED :: mu_x = 0._dp ! spatial x-Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_y = 0._dp ! spatial y-Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_z = 0._dp ! spatial z-Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability) INTEGER, PUBLIC, PROTECTED :: N_HD = 4 ! order of numerical spatial diffusion REAL(dp), PUBLIC, PROTECTED :: nu = 0._dp ! Collision frequency REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: sigma_e = 0.023338_dp! sqrt of electron ion mass ratio REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp ! - REAL(dp), PUBLIC, PROTECTED :: K_n = 1._dp ! Density drive - REAL(dp), PUBLIC, PROTECTED :: K_T = 0._dp ! Temperature drive + REAL(dp), PUBLIC, PROTECTED :: k_N = 1._dp ! Ion density drive + REAL(dp), PUBLIC, PROTECTED :: eta_N = 1._dp ! electron-ion density drive ratio (k_Ne/k_Ni) + REAL(dp), PUBLIC, PROTECTED :: k_T = 0._dp ! Temperature drive + REAL(dp), PUBLIC, PROTECTED :: eta_T = 0._dp ! electron-ion temperature drive ratio (k_Te/k_Ti) REAL(dp), PUBLIC, PROTECTED :: K_E = 0._dp ! Backg. electric field drive REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length REAL(dp), PUBLIC, PROTECTED :: beta = 0._dp ! electron plasma Beta (8piNT_e/B0^2) REAL(dp), PUBLIC, PROTECTED :: taue_qe ! factor of the magnetic moment coupling REAL(dp), PUBLIC, PROTECTED :: taui_qi ! REAL(dp), PUBLIC, PROTECTED :: qi_taui ! REAL(dp), PUBLIC, PROTECTED :: qe_taue ! REAL(dp), PUBLIC, PROTECTED :: sqrtTaue_qe ! factor of parallel moment term REAL(dp), PUBLIC, PROTECTED :: sqrtTaui_qi ! REAL(dp), PUBLIC, PROTECTED :: qe_sigmae_sqrtTaue ! factor of parallel phi term REAL(dp), PUBLIC, PROTECTED :: qi_sigmai_sqrtTaui ! REAL(dp), PUBLIC, PROTECTED :: sigmae2_taue_o2 ! factor of the Kernel argument REAL(dp), PUBLIC, PROTECTED :: sigmai2_taui_o2 ! REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmae2_taue_o2 ! factor of the Kernel argument REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmai2_taui_o2 REAL(dp), PUBLIC, PROTECTED :: nu_e, nu_i ! electron-ion, ion-ion collision frequency REAL(dp), PUBLIC, PROTECTED :: nu_ee, nu_ie ! e-e, i-e coll. frequ. REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum REAL(dp), PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i REAL(dp), PUBLIC, PROTECTED :: sqrt_tau_o_sigma_e, sqrt_tau_o_sigma_i PUBLIC :: model_readinputs, model_outputinputs CONTAINS SUBROUTINE model_readinputs ! Read the input parameters USE basic, ONLY : lu_in, my_id USE prec_const IMPLICIT NONE NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, & mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,& tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,& - K_n, K_T, K_E, GradB, CurvB, lambdaD, beta + k_N, eta_N, k_T, eta_T, K_E, GradB, CurvB, lambdaD, beta READ(lu_in,model_par) IF(.NOT. KIN_E) THEN IF(my_id.EQ.0) print*, 'Adiabatic electron model -> beta = 0' beta = 0._dp ENDIF taue_qe = tau_e/q_e ! factor of the magnetic moment coupling taui_qi = tau_i/q_i ! factor of the magnetic moment coupling qe_taue = q_e/tau_e qi_taui = q_i/tau_i sqrtTaue_qe = sqrt(tau_e)/q_e ! factor of parallel moment term sqrtTaui_qi = sqrt(tau_i)/q_i ! factor of parallel moment term qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i) qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum qi2_taui = (q_i**2)/tau_i sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp sqrt_sigmae2_taue_o2 = SQRT(sigma_e**2 * tau_e/2._dp) ! to avoid multiple SQRT eval sqrt_sigmai2_taui_o2 = SQRT(sigma_i**2 * tau_i/2._dp) q_o_sqrt_tau_sigma_e = q_e/SQRT(tau_e)/sigma_e ! For psi field terms q_o_sqrt_tau_sigma_i = q_i/SQRT(tau_i)/sigma_i ! For psi field terms sqrt_tau_o_sigma_e = SQRT(tau_e)/sigma_e ! For Ampere eq sqrt_tau_o_sigma_i = SQRT(tau_i)/sigma_i !! We use the ion-ion collision as normalization with definition ! nu_ii = 4 sqrt(pi)/3 T_i^(-3/2) m_i^(-1/2) q^4 n_i0 ln(Lambda) ! nu_e = nu/sigma_e * (tau_e)**(3._dp/2._dp) ! electron-ion collision frequency (where already multiplied by 0.532) nu_i = nu ! ion-ion collision frequ. nu_ee = nu_e ! e-e coll. frequ. nu_ie = nu_i ! i-e coll. frequ. ! Old normalization (MOLI Jorge/Frei) ! nu_e = 0.532_dp*nu ! electron-ion collision frequency (where already multiplied by 0.532) ! nu_i = 0.532_dp*nu*sigma_e*tau_e**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ. ! nu_ee = nu_e/SQRT2 ! e-e coll. frequ. ! nu_ie = 0.532_dp*nu*sigma_e**2 ! i-e coll. frequ. END SUBROUTINE model_readinputs SUBROUTINE model_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "CLOS", CLOS) CALL attach(fidres, TRIM(str), "KERN", KERN) CALL attach(fidres, TRIM(str), "LINEARITY", LINEARITY) CALL attach(fidres, TRIM(str), "KIN_E", KIN_E) CALL attach(fidres, TRIM(str), "nu", nu) CALL attach(fidres, TRIM(str), "mu", 0) CALL attach(fidres, TRIM(str), "mu_x", mu_x) CALL attach(fidres, TRIM(str), "mu_y", mu_y) CALL attach(fidres, TRIM(str), "mu_z", mu_z) CALL attach(fidres, TRIM(str), "mu_p", mu_p) CALL attach(fidres, TRIM(str), "mu_j", mu_j) CALL attach(fidres, TRIM(str), "tau_e", tau_e) CALL attach(fidres, TRIM(str), "tau_i", tau_i) CALL attach(fidres, TRIM(str), "sigma_e", sigma_e) CALL attach(fidres, TRIM(str), "sigma_i", sigma_i) CALL attach(fidres, TRIM(str), "q_e", q_e) CALL attach(fidres, TRIM(str), "q_i", q_i) - CALL attach(fidres, TRIM(str), "K_n", K_n) - CALL attach(fidres, TRIM(str), "K_T", K_T) + CALL attach(fidres, TRIM(str), "k_N", k_N) + CALL attach(fidres, TRIM(str), "eta_N", eta_N) + CALL attach(fidres, TRIM(str), "k_T", k_T) + CALL attach(fidres, TRIM(str), "eta_T", eta_T) CALL attach(fidres, TRIM(str), "K_E", K_E) CALL attach(fidres, TRIM(str), "lambdaD", lambdaD) END SUBROUTINE model_outputinputs END MODULE model diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90 index 2fe44a5..092b390 100644 --- a/src/moments_eq_rhs_mod.F90 +++ b/src/moments_eq_rhs_mod.F90 @@ -1,389 +1,389 @@ MODULE moments_eq_rhs IMPLICIT NONE PUBLIC :: compute_moments_eq_rhs CONTAINS SUBROUTINE compute_moments_eq_rhs USE model, only: KIN_E IMPLICIT NONE IF(KIN_E) CALL moments_eq_rhs_e CALL moments_eq_rhs_i !! BETA TESTING !! IF(.false.) CALL add_Maxwellian_background_terms END SUBROUTINE compute_moments_eq_rhs !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Electrons moments RHS !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_e USE basic USE time_integration USE array USE fields USE grid USE model USE prec_const USE collision use geometry USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2, dzlnB_o_J COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: i_kx,i_ky,phikykxz, psikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloope : DO iz = izs,ize kxloope : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector i_kx = imagu * kx ! toroidal derivative kyloope : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value psikykxz = psi(iky,ikx,iz)! tmp psi value ! Kinetic loops jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1 j_int = jarray_e(ij) ploope : DO ip = ips_e, ipe_e ! Hermite loop p_int = parray_e(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_e^{p,j} Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz) ! term propto n_e^{p+2,j} Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz) ! term propto n_e^{p-2,j} Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz) ! term propto n_e^{p,j+1} Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz) ! term propto n_e^{p,j-1} Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz) ! Parallel dynamic ! ddz derivative for Landau damping term Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) & + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz) ! Mirror terms Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz) Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz) Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) ! Trapping terms Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz) Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + xphijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + xphijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz ELSE Tphi = 0._dp ENDIF !! Vector potential term IF ( (p_int .LE. 3) .AND. (p_int .GT. 1) ) THEN ! Kronecker p1 or p3 Tpsi = (xpsij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + xpsijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + xpsijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*psikykxz ELSE Tpsi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)& ! Parallel coupling (Landau Damping) - Tpar*gradz_coeff(iz,eo) & ! Mirror term (parallel magnetic gradient) - gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) & ! Drives (density + temperature gradients) - i_ky * (Tphi - Tpsi) & ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) - mu_x*diff_kx_coeff*kx**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & - mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" see Pueschel 2010 (eq 25) - mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) & ! Collision term + TColl_e(ip,ij,iky,ikx,iz) & ! Nonlinear term - hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz) IF(ip-4 .GT. 0) & ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) & + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel) ELSE moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploope END DO jloope END DO kyloope END DO kxloope END DO zloope ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_e !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_i USE basic USE time_integration, ONLY: updatetlevel USE array USE fields USE grid USE model USE prec_const USE collision USE geometry USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1 COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi COMPLEX(dp) :: i_kx, i_ky, phikykxz, psikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloopi : DO iz = izs,ize kxloopi : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector i_kx = imagu * kx ! toroidal derivative kyloopi : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value psikykxz = psi(iky,ikx,iz)! tmp phi value ! Kinetic loops jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 j_int = jarray_i(ij) ploopi : DO ip = ips_i, ipe_i ! Hermite loop p_int = parray_i(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_i^{p,j} Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz) ! term propto n_i^{p+2,j} Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz) ! term propto n_i^{p-2,j} Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz) ! term propto n_i^{p,j+1} Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz) ! term propto n_i^{p,j-1} Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz) ! Tperp Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 ! Parallel dynamic ! ddz derivative for Landau damping term Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) & + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz) ! Mirror terms Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz) Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz) Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) ! Trapping terms Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz) Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz ELSE Tphi = 0._dp ENDIF !! Vector potential term IF ( (p_int .LE. 3) .AND. (p_int .GT. 1) ) THEN ! Kronecker p1 or p3 Tpsi = (xpsij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + xpsijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + xpsijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*psikykxz ELSE Tpsi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp & ! Parallel coupling (Landau damping) - gradz_coeff(iz,eo) * Tpar & ! Mirror term (parallel magnetic gradient) - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir & ! Drives (density + temperature gradients) - i_ky * (Tphi - Tpsi) & ! Numerical hyperdiffusion (totally artificial, for stability purpose) - mu_x*diff_kx_coeff*kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & - mu_y*diff_ky_coeff*ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" + mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) & ! Collision term + TColl_i(ip,ij,iky,ikx,iz)& ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor - hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz) IF(ip-4 .GT. 0) & ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) & + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel) ELSE moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploopi END DO jloopi END DO kyloopi END DO kxloopi END DO zloopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_i SUBROUTINE add_Maxwellian_background_terms ! This routine is meant to add the terms rising from the magnetic operator, ! i.e. (B x GradB) Grad, applied on the background Maxwellian distribution ! (x_a + spar^2)(b x GradB) GradFaM ! It gives birth to kx=ky=0 sources terms (averages) that hit moments 00, 20, ! 40, 01,02, 21 with background gradient dependences. USE prec_const USE time_integration, ONLY : updatetlevel - USE model, ONLY: taue_qe, taui_qi, K_N, K_T, KIN_E + USE model, ONLY: taue_qe, taui_qi, k_N, eta_N, k_T, eta_T, KIN_E USE array, ONLY: moments_rhs_e, moments_rhs_i USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,& ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,& jarray_e,parray_e,jarray_i,parray_i, zarray, izs,ize,& ip,ij IMPLICIT NONE real(dp), DIMENSION(izs:ize) :: sinz sinz(izs:ize) = SIN(zarray(izs:ize,0)) IF(contains_kx0 .AND. contains_ky0) THEN IF(KIN_E) THEN DO ip = ips_e,ipe_e DO ij = ijs_e,ije_e SELECT CASE(ij-1) CASE(0) ! j = 0 SELECT CASE (ip-1) CASE(0) ! Na00 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taue_qe * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T) + +taue_qe * sinz(izs:ize) * (1.5_dp*eta_N*k_N - 1.125_dp*eta_N*k_T) CASE(2) ! Na20 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T) + +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*eta_N*k_N - 2.75_dp*eta_T*k_T) CASE(4) ! Na40 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*K_T + +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*eta_T*k_T END SELECT CASE(1) ! j = 1 SELECT CASE (ip-1) CASE(0) ! Na01 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - -taue_qe * sinz(izs:ize) * (K_N + 3.5_dp*K_T) + -taue_qe * sinz(izs:ize) * (eta_N*k_N + 3.5_dp*eta_T*k_T) CASE(2) ! Na21 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - -taue_qe * sinz(izs:ize) * SQRT2*K_T + -taue_qe * sinz(izs:ize) * SQRT2*eta_T*k_T END SELECT CASE(2) ! j = 2 SELECT CASE (ip-1) CASE(0) ! Na02 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taue_qe * sinz(izs:ize) * 2._dp*K_T + +taue_qe * sinz(izs:ize) * 2._dp*eta_T*k_T END SELECT END SELECT ENDDO ENDDO ENDIF DO ip = ips_i,ipe_i DO ij = ijs_i,ije_i SELECT CASE(ij-1) CASE(0) ! j = 0 SELECT CASE (ip-1) CASE(0) ! Na00 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taui_qi * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T) + +taui_qi * sinz(izs:ize) * (1.5_dp*k_N - 1.125_dp*k_T) CASE(2) ! Na20 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T) + +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_N - 2.75_dp*k_T) CASE(4) ! Na40 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*K_T + +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_T END SELECT CASE(1) ! j = 1 SELECT CASE (ip-1) CASE(0) ! Na01 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - -taui_qi * sinz(izs:ize) * (K_N + 3.5_dp*K_T) + -taui_qi * sinz(izs:ize) * (k_N + 3.5_dp*k_T) CASE(2) ! Na21 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - -taui_qi * sinz(izs:ize) * SQRT2*K_T + -taui_qi * sinz(izs:ize) * SQRT2*k_T END SELECT CASE(2) ! j = 2 SELECT CASE (ip-1) CASE(0) ! Na02 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& - +taui_qi * sinz(izs:ize) * 2._dp*K_T + +taui_qi * sinz(izs:ize) * 2._dp*k_T END SELECT END SELECT ENDDO ENDDO ENDIF END SUBROUTINE END MODULE moments_eq_rhs diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90 index 1ede5ac..ccfb5f4 100644 --- a/src/numerics_mod.F90 +++ b/src/numerics_mod.F90 @@ -1,475 +1,476 @@ MODULE numerics USE basic USE prec_const USE grid USE utility USE coeff implicit none PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, evaluate_ampere_op PUBLIC :: compute_lin_coeff, play_with_modes, save_EM_ZF_modes CONTAINS !******************************************************************************! !!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin !******************************************************************************! SUBROUTINE build_dnjs_table USE array, Only : dnjs USE coeff IMPLICIT NONE INTEGER :: in, ij, is, J INTEGER :: n_, j_, s_ J = max(jmaxe,jmaxi) DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry n_ = in - 1 DO ij = in,J+1 j_ = ij - 1 DO is = ij,J+1 s_ = is - 1 dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0)) ! By symmetry dnjs(in,is,ij) = dnjs(in,ij,is) dnjs(ij,in,is) = dnjs(in,ij,is) dnjs(ij,is,in) = dnjs(in,ij,is) dnjs(is,ij,in) = dnjs(in,ij,is) dnjs(is,in,ij) = dnjs(in,ij,is) ENDDO ENDDO ENDDO END SUBROUTINE build_dnjs_table !******************************************************************************! !******************************************************************************! !!!!!!! Evaluate the kernels once for all !******************************************************************************! SUBROUTINE evaluate_kernels USE basic USE array, Only : kernel_e, kernel_i, HF_phi_correction_operator USE grid USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, & lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2, KIN_E IMPLICIT NONE INTEGER :: j_int REAL(dp) :: j_dp, y_ DO eo = 0,1 DO ikx = ikxs,ikxe DO iky = ikys,ikye DO iz = izgs,izge !!!!! Electron kernels !!!!! IF(KIN_E) THEN DO ij = ijgs_e, ijge_e j_int = jarray_e(ij) j_dp = REAL(j_int,dp) y_ = sigmae2_taue_o2 * kparray(iky,ikx,iz,eo)**2 kernel_e(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj ENDDO IF (ijs_e .EQ. 1) & kernel_e(ijgs_e,iky,ikx,iz,eo) = 0._dp ENDIF !!!!! Ion kernels !!!!! DO ij = ijgs_i, ijge_i j_int = jarray_i(ij) j_dp = REAL(j_int,dp) y_ = sigmai2_taui_o2 * kparray(iky,ikx,iz,eo)**2 kernel_i(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj ENDDO IF (ijs_i .EQ. 1) & kernel_i(ijgs_i,iky,ikx,iz,eo) = 0._dp ENDDO ENDDO ENDDO ENDDO !! Correction term for the evaluation of the heat flux HF_phi_correction_operator(ikys:ikye,ikxs:ikxe,izs:ize) = & 2._dp * Kernel_i(1,ikys:ikye,ikxs:ikxe,izs:ize,0) & -1._dp * Kernel_i(2,ikys:ikye,ikxs:ikxe,izs:ize,0) DO ij = ijs_i, ije_i j_int = jarray_i(ij) j_dp = REAL(j_int,dp) HF_phi_correction_operator(ikys:ikye,ikxs:ikxe,izs:ize) = HF_phi_correction_operator(ikys:ikye,ikxs:ikxe,izs:ize) & - Kernel_i(ij,ikys:ikye,ikxs:ikxe,izs:ize,0) * (& 2._dp*(j_dp+1.5_dp) * Kernel_i(ij ,ikys:ikye,ikxs:ikxe,izs:ize,0) & - (j_dp+1.0_dp) * Kernel_i(ij+1,ikys:ikye,ikxs:ikxe,izs:ize,0) & - j_dp * Kernel_i(ij-1,ikys:ikye,ikxs:ikxe,izs:ize,0)) ENDDO END SUBROUTINE evaluate_kernels !******************************************************************************! !******************************************************************************! !!!!!!! Evaluate inverse polarisation operator for Poisson equation !******************************************************************************! SUBROUTINE evaluate_poisson_op USE basic USE array, Only : kernel_e, kernel_i, inv_poisson_op, inv_pol_ion USE grid USE model, ONLY : tau_e, tau_i, q_e, q_i, KIN_E IMPLICIT NONE REAL(dp) :: pol_i, pol_e ! (Z_a^2/tau_a (1-sum_n kernel_na^2)) INTEGER :: ini,ine ! This term has no staggered grid dependence. It is evalued for the ! even z grid since poisson uses p=0 moments and phi only. kxloop: DO ikx = ikxs,ikxe kyloop: DO iky = ikys,ikye zloop: DO iz = izs,ize IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN inv_poisson_op(iky, ikx, iz) = 0._dp ELSE !!!!!!!!!!!!!!!!! Ion contribution ! loop over n only if the max polynomial degree pol_i = 0._dp DO ini=1,jmaxi+1 pol_i = pol_i + qi2_taui*kernel_i(ini,iky,ikx,iz,0)**2 ! ... sum recursively ... END DO !!!!!!!!!!!!! Electron contribution pol_e = 0._dp ! Kinetic model IF (KIN_E) THEN ! loop over n only if the max polynomial degree DO ine=1,jmaxe+1 ! ine = n+1 pol_e = pol_e + qe2_taue*kernel_e(ine,iky,ikx,iz,0)**2 ! ... sum recursively ... END DO ! Adiabatic model ELSE pol_e = qe2_taue - 1._dp ENDIF inv_poisson_op(iky, ikx, iz) = 1._dp/(qe2_taue + qi2_taui - pol_i - pol_e) inv_pol_ion (iky, ikx, iz) = 1._dp/(qi2_taui - pol_i) ENDIF END DO zloop END DO kyloop END DO kxloop END SUBROUTINE evaluate_poisson_op !******************************************************************************! !******************************************************************************! !!!!!!! Evaluate inverse polarisation operator for Poisson equation !******************************************************************************! SUBROUTINE evaluate_ampere_op USE basic USE array, Only : kernel_e, kernel_i, inv_ampere_op USE grid USE model, ONLY : tau_e, tau_i, q_e, q_i, KIN_E, beta IMPLICIT NONE REAL(dp) :: pol_i, pol_e, kperp2 ! (Z_a^2/tau_a (1-sum_n kernel_na^2)) INTEGER :: ini,ine ! We do not solve Ampere if beta = 0 to spare waste of ressources IF(SOLVE_AMPERE) THEN ! This term has no staggered grid dependence. It is evalued for the ! even z grid since poisson uses p=0 moments and phi only. kxloop: DO ikx = ikxs,ikxe kyloop: DO iky = ikys,ikye zloop: DO iz = izs,ize kperp2 = kparray(iky,ikx,iz,0)**2 IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN inv_ampere_op(iky, ikx, iz) = 0._dp ELSE !!!!!!!!!!!!!!!!! Ion contribution ! loop over n only if the max polynomial degree pol_i = 0._dp DO ini=1,jmaxi+1 pol_i = pol_i + kernel_i(ini,iky,ikx,iz,0)**2 ! ... sum recursively ... END DO pol_i = q_i**2/(sigma_i**2) * pol_i !!!!!!!!!!!!! Electron contribution pol_e = 0._dp ! loop over n only if the max polynomial degree DO ine=1,jmaxe+1 ! ine = n+1 pol_e = pol_e + kernel_e(ine,iky,ikx,iz,0)**2 ! ... sum recursively ... END DO pol_e = q_e**2/(sigma_e**2) * pol_e inv_ampere_op(iky, ikx, iz) = 1._dp/(2._dp*kperp2 + beta*(pol_i + pol_e)) ENDIF END DO zloop END DO kyloop END DO kxloop ENDIF END SUBROUTINE evaluate_ampere_op !******************************************************************************! SUBROUTINE compute_lin_coeff USE array USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, & - K_T, K_n, CurvB, GradB, KIN_E, tau_e, tau_i, sigma_e, sigma_i + k_T, eta_T, k_N, eta_N, CurvB, GradB, KIN_E,& + tau_e, tau_i, sigma_e, sigma_i USE prec_const USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, & ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i IMPLICIT NONE INTEGER :: p_int, j_int ! polynom. degrees REAL(dp) :: p_dp, j_dp REAL(dp) :: kx, ky, z !! Electrons linear coefficients for moment RHS !!!!!!!!!! IF(KIN_E)THEN DO ip = ips_e, ipe_e p_int= parray_e(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree DO ij = ijs_e, ije_e j_int= jarray_e(ij) ! Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree ! All Napj terms xnepj(ip,ij) = taue_qe*(CurvB*(2._dp*p_dp + 1._dp) & +GradB*(2._dp*j_dp + 1._dp)) ! Mirror force terms ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp) ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp) ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp) ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp) zNepm1j (ip,ij) = +SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp) zNepm1jp1(ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(p_dp) zNepm1jm1(ip,ij) = -SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp) ENDDO ENDDO DO ip = ips_e, ipe_e p_int= parray_e(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree ! Landau damping coefficients (ddz napj term) xnepp1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp + 1_dp) xnepm1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp) ! Magnetic curvature term xnepp2j(ip) = taue_qe * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) xnepm2j(ip) = taue_qe * CurvB * SQRT(p_dp * (p_dp - 1._dp)) ENDDO DO ij = ijs_e, ije_e j_int= jarray_e(ij) ! Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree ! Magnetic gradient term xnepjp1(ij) = -taue_qe * GradB * (j_dp + 1._dp) xnepjm1(ij) = -taue_qe * GradB * j_dp ENDDO ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Ions linear coefficients for moment RHS !!!!!!!!!! DO ip = ips_i, ipe_i p_int= parray_i(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree DO ij = ijs_i, ije_i j_int= jarray_i(ij) ! Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree ! All Napj terms xnipj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) & +GradB*(2._dp*j_dp + 1._dp)) ! Mirror force terms ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp+1._dp) ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp) ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp+1._dp) ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp) ! Trapping terms zNipm1j (ip,ij) = +SQRT(tau_i)/sigma_i* (2._dp*j_dp+1_dp)*SQRT(p_dp) zNipm1jp1(ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1_dp)*SQRT(p_dp) zNipm1jm1(ip,ij) = -SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp) ENDDO ENDDO DO ip = ips_i, ipe_i p_int= parray_i(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree ! Landau damping coefficients (ddz napj term) xnipp1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp + 1._dp) xnipm1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp) ! Magnetic curvature term xnipp2j(ip) = taui_qi * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) xnipm2j(ip) = taui_qi * CurvB * SQRT(p_dp * (p_dp - 1._dp)) ENDDO DO ij = ijs_i, ije_i j_int= jarray_i(ij) ! Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree ! Magnetic gradient term xnipjp1(ij) = -taui_qi * GradB * (j_dp + 1._dp) xnipjm1(ij) = -taui_qi * GradB * j_dp ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! ES linear coefficients for moment RHS !!!!!!!!!! IF (KIN_E) THEN DO ip = ips_e, ipe_e p_int= parray_e(ip) ! Hermite degree DO ij = ijs_e, ije_e j_int= jarray_e(ij) ! REALof Laguerre degree j_dp = REAL(j_int,dp) ! REALof Laguerre degree !! Electrostatic potential pj terms IF (p_int .EQ. 0) THEN ! kronecker p0 - xphij_e(ip,ij) =+K_n + 2.*j_dp*K_T - xphijp1_e(ip,ij) =-K_T*(j_dp+1._dp) - xphijm1_e(ip,ij) =-K_T* j_dp + xphij_e(ip,ij) =+eta_N*k_N + 2.*j_dp*eta_T*k_T + xphijp1_e(ip,ij) =-eta_T*k_T*(j_dp+1._dp) + xphijm1_e(ip,ij) =-eta_T*k_T* j_dp ELSE IF (p_int .EQ. 2) THEN ! kronecker p2 - xphij_e(ip,ij) =+K_T/SQRT2 + xphij_e(ip,ij) =+eta_T*k_T/SQRT2 xphijp1_e(ip,ij) = 0._dp; xphijm1_e(ip,ij) = 0._dp; ELSE xphij_e(ip,ij) = 0._dp; xphijp1_e(ip,ij) = 0._dp xphijm1_e(ip,ij) = 0._dp; ENDIF ENDDO ENDDO ENDIF DO ip = ips_i, ipe_i p_int= parray_i(ip) ! Hermite degree DO ij = ijs_i, ije_i j_int= jarray_i(ij) ! REALof Laguerre degree j_dp = REAL(j_int,dp) ! REALof Laguerre degree !! Electrostatic potential pj terms IF (p_int .EQ. 0) THEN ! kronecker p0 - xphij_i(ip,ij) =+K_n + 2._dp*j_dp*K_T - xphijp1_i(ip,ij) =-K_T*(j_dp+1._dp) - xphijm1_i(ip,ij) =-K_T* j_dp + xphij_i(ip,ij) =+k_N + 2._dp*j_dp*k_T + xphijp1_i(ip,ij) =-k_T*(j_dp+1._dp) + xphijm1_i(ip,ij) =-k_T* j_dp ELSE IF (p_int .EQ. 2) THEN ! kronecker p2 - xphij_i(ip,ij) =+K_T/SQRT2 + xphij_i(ip,ij) =+k_T/SQRT2 xphijp1_i(ip,ij) = 0._dp; xphijm1_i(ip,ij) = 0._dp; ELSE xphij_i(ip,ij) = 0._dp; xphijp1_i(ip,ij) = 0._dp xphijm1_i(ip,ij) = 0._dp; ENDIF ENDDO ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! EM linear coefficients for moment RHS !!!!!!!!!! IF (KIN_E) THEN DO ip = ips_e, ipe_e p_int= parray_e(ip) ! Hermite degree DO ij = ijs_e, ije_e j_int= jarray_e(ij) ! REALof Laguerre degree j_dp = REAL(j_int,dp) ! REALof Laguerre degree !! Electrostatic potential pj terms IF (p_int .EQ. 1) THEN ! kronecker p1 - xpsij_e(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_e)/sigma_e - xpsijp1_e(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e - xpsijm1_e(ip,ij) =- K_T* j_dp * SQRT(tau_e)/sigma_e + xpsij_e(ip,ij) =+(eta_N*k_N + (2._dp*j_dp+1._dp)*eta_T*k_T) * SQRT(tau_e)/sigma_e + xpsijp1_e(ip,ij) =- eta_T*k_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e + xpsijm1_e(ip,ij) =- eta_T*k_T* j_dp * SQRT(tau_e)/sigma_e ELSE IF (p_int .EQ. 3) THEN ! kronecker p3 - xpsij_e(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e + xpsij_e(ip,ij) =+ eta_T*k_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e xpsijp1_e(ip,ij) = 0._dp; xpsijm1_e(ip,ij) = 0._dp; ELSE xpsij_e(ip,ij) = 0._dp; xpsijp1_e(ip,ij) = 0._dp xpsijm1_e(ip,ij) = 0._dp; ENDIF ENDDO ENDDO ENDIF DO ip = ips_i, ipe_i p_int= parray_i(ip) ! Hermite degree DO ij = ijs_i, ije_i j_int= jarray_i(ij) ! REALof Laguerre degree j_dp = REAL(j_int,dp) ! REALof Laguerre degree !! Electrostatic potential pj terms IF (p_int .EQ. 1) THEN ! kronecker p1 - xpsij_i(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_i)/sigma_i - xpsijp1_i(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i - xpsijm1_i(ip,ij) =- K_T* j_dp * SQRT(tau_i)/sigma_i + xpsij_i(ip,ij) =+(k_N + (2._dp*j_dp+1._dp)*k_T) * SQRT(tau_i)/sigma_i + xpsijp1_i(ip,ij) =- k_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i + xpsijm1_i(ip,ij) =- k_T* j_dp * SQRT(tau_i)/sigma_i ELSE IF (p_int .EQ. 3) THEN ! kronecker p3 - xpsij_i(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i + xpsij_i(ip,ij) =+ k_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i xpsijp1_i(ip,ij) = 0._dp; xpsijm1_i(ip,ij) = 0._dp; ELSE xpsij_i(ip,ij) = 0._dp; xpsijp1_i(ip,ij) = 0._dp xpsijm1_i(ip,ij) = 0._dp; ENDIF ENDDO ENDDO END SUBROUTINE compute_lin_coeff !******************************************************************************! !!!!!!! Routine that can artificially increase or wipe modes !******************************************************************************! SUBROUTINE save_EM_ZF_modes USE fields USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM USE grid USE time_integration, ONLY: updatetlevel USE initial_par, ONLY: ACT_ON_MODES USE model, ONLY: KIN_E IMPLICIT NONE ! Store Zonal and entropy modes IF(contains_ky0) THEN IF(KIN_E) & moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) phi_ZF(ikxs:ikxe,izs:ize) = phi(iky_0,ikxs:ikxe,izs:ize) ELSE IF(KIN_E) & moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = 0._dp moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = 0._dp phi_ZF(ikxs:ikxe,izs:ize) = 0._dp ENDIF IF(contains_kx0) THEN IF(KIN_E) & moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) phi_EM(ikys:ikye,izs:ize) = phi(ikys:ikye,ikx_0,izs:ize) ELSE IF(KIN_E) & moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = 0._dp moments_i_EM(ips_e:ipe_e,ijs_i:ije_i,ikys:ikye,izs:ize) = 0._dp phi_EM(ikys:ikye,izs:ize) = 0._dp ENDIF END SUBROUTINE SUBROUTINE play_with_modes USE fields USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM USE grid USE time_integration, ONLY: updatetlevel USE initial_par, ONLY: ACT_ON_MODES USE model, ONLY: KIN_E IMPLICIT NONE REAL(dp) :: AMP = 1.5_dp SELECT CASE(ACT_ON_MODES) CASE('wipe_zonal') ! Errase the zonal flow IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp phi(iky_0,ikxs:ikxe,izs:ize) = 0._dp CASE('wipe_entropymode') IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp phi(ikys:ikye,ikx_0,izs:ize) = 0._dp CASE('wipe_turbulence') DO ikx = ikxs,ikxe DO iky = ikys, ikye IF ( (ikx .NE. ikx_0) .AND. (iky .NE. iky_0) ) THEN IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp phi(iky,ikx,izs:ize) = 0._dp ENDIF ENDDO ENDDO CASE('wipe_nonzonal') DO ikx = ikxs,ikxe DO iky = ikys, ikye IF ( (ikx .NE. ikx_0) ) THEN IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp phi(iky,ikx,izs:ize) = 0._dp ENDIF ENDDO ENDDO CASE('freeze_zonal') IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) phi(iky_0,ikxs:ikxe,izs:ize) = phi_ZF(ikxs:ikxe,izs:ize) CASE('freeze_entropymode') IF(contains_kx0) THEN IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize) phi(ikys:ikye,ikx_0,izs:ize) = phi_EM(ikys:ikye,izs:ize) ENDIF CASE('amplify_zonal') IF(KIN_E) & moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) phi(iky_0,ikxs:ikxe,izs:ize) = AMP*phi_ZF(ikxs:ikxe,izs:ize) END SELECT END SUBROUTINE END MODULE numerics