diff --git a/src/model_mod.F90 b/src/model_mod.F90 index 12d16a584c117061a45101ae009e2246dca5e291..1f1406fe38b2c3672a22030323d312e1e3cc5bd6 100644 --- a/src/model_mod.F90 +++ b/src/model_mod.F90 @@ -24,8 +24,10 @@ MODULE model 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 @@ -62,7 +64,7 @@ CONTAINS 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) @@ -132,8 +134,10 @@ CONTAINS 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 diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90 index 2fe44a5454cc09fe7eace8128e92f65de54944cd..092b390dfb8cad2170c7880bf3f2cab6f25b2bbd 100644 --- a/src/moments_eq_rhs_mod.F90 +++ b/src/moments_eq_rhs_mod.F90 @@ -300,7 +300,7 @@ SUBROUTINE add_Maxwellian_background_terms ! 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,& @@ -320,28 +320,28 @@ SUBROUTINE add_Maxwellian_background_terms 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 @@ -355,28 +355,28 @@ SUBROUTINE add_Maxwellian_background_terms 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 diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90 index 1ede5ac6f9b392e781d09ea75deb627fb7891e06..ccfb5f4baea6ebe63f035542f239d316afdb5359 100644 --- a/src/numerics_mod.F90 +++ b/src/numerics_mod.F90 @@ -201,7 +201,8 @@ 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 @@ -297,11 +298,11 @@ SUBROUTINE compute_lin_coeff 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 @@ -317,11 +318,11 @@ SUBROUTINE compute_lin_coeff 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 @@ -339,11 +340,11 @@ SUBROUTINE compute_lin_coeff 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 @@ -359,11 +360,11 @@ SUBROUTINE compute_lin_coeff 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