diff --git a/src/inital.F90 b/src/inital.F90 index c50e0bd..8751736 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -1,341 +1,353 @@ !******************************************************************************! !!!!!! initialize the moments and load/build coeff tables !******************************************************************************! SUBROUTINE inital USE basic USE model, ONLY : CO, NON_LIN USE initial_par USE prec_const USE coeff USE time_integration USE array, ONLY : Sepj,Sipj USE collision USE closure USE ghosts USE restarts implicit none CALL set_updatetlevel(1) IF (my_id .EQ. 0) WRITE(*,*) 'Evaluate kernels' CALL evaluate_kernels !!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!! IF ( RESTART ) THEN IF (my_id .EQ. 0) WRITE(*,*) 'Load moments' CALL load_moments ! get N_0 CALL poisson ! compute phi_0=phi(N_0) ELSE IF (INIT_NOISY_PHI) THEN IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi' CALL init_phi ! init noisy phi_0, N_0 = 0 ELSEIF (INIT_ZF .GT. 0) THEN IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi' CALL init_phi ! init ZF ELSE IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments and ghosts' CALL init_moments ! init noisy N_0 IF (my_id .EQ. 0) WRITE(*,*) 'Init phi with Poisson' CALL poisson ! get phi_0 = phi(N_0) ENDIF ENDIF IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure' CALL apply_closure_model IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication' CALL update_ghosts !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!! IF ( NON_LIN ) THEN; IF (my_id .EQ. 0) WRITE(*,*) 'Building Dnjs table' CALL build_dnjs_table IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj' CALL compute_Sapj ! compute S_0 = S(phi_0,N_0) ENDIF !!!!!! Load the COSOlver collision operator coefficients !!!!!! IF (ABS(CO) .GT. 1) THEN CALL load_COSOlver_mat ! Compute collision CALL compute_TColl ! compute C_0 = C(N_0) ENDIF END SUBROUTINE inital !******************************************************************************! !******************************************************************************! !!!!!!! Initialize the moments randomly !******************************************************************************! SUBROUTINE init_moments USE basic USE grid USE fields USE prec_const USE utility, ONLY: checkfield USE initial_par USE model, ONLY : NON_LIN IMPLICIT NONE REAL(dp) :: noise REAL(dp) :: kx, ky, sigma, gain, ky_shift INTEGER, DIMENSION(12) :: iseedarr ! Seed random number generator iseedarr(:)=iseed CALL RANDOM_SEED(PUT=iseedarr+my_id) !**** Broad noise initialization ******************************************* DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e DO ikx=ikxs,ikxe DO iky=ikys,ikye DO iz=izs,ize CALL RANDOM_NUMBER(noise) moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp)) END DO END DO END DO IF ( contains_kx0 ) THEN DO iky=2,Nky/2 !symmetry at kx = 0 for all z moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :) END DO ENDIF END DO END DO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i DO ikx=ikxs,ikxe DO iky=ikys,ikye DO iz=izs,ize CALL RANDOM_NUMBER(noise) moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp)) END DO END DO END DO IF ( contains_kx0 ) THEN DO iky=2,Nky/2 !symmetry at kx = 0 for all z moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:) END DO ENDIF END DO END DO ! Putting to zero modes that are not in the 2/3 Orszag rule IF (NON_LIN) THEN DO ikx=ikxs,ikxe DO iky=ikys,ikye DO iz=izs,ize DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky) ENDDO ENDDO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky) ENDDO ENDDO ENDDO ENDDO ENDDO ENDIF END SUBROUTINE init_moments !******************************************************************************! !******************************************************************************! !!!!!!! Initialize a noisy ES potential and cancel the moments !******************************************************************************! SUBROUTINE init_phi USE basic USE grid USE fields USE prec_const USE initial_par IMPLICIT NONE REAL(dp) :: noise REAL(dp) :: kx, ky, sigma, gain, ky_shift INTEGER, DIMENSION(12) :: iseedarr IF (INIT_NOISY_PHI) THEN ! Seed random number generator iseedarr(:)=iseed CALL RANDOM_SEED(PUT=iseedarr+my_id) !**** noise initialization ******************************************* DO ikx=ikxs,ikxe DO iky=ikys,ikye DO iz=izs,ize CALL RANDOM_NUMBER(noise) phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))*AA_x(ikx)*AA_y(iky) ENDDO END DO END DO !symmetry at kx = 0 to keep real inverse transform IF ( contains_kx0 ) THEN DO iky=2,Nky/2 phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:) END DO phi(ikx_0,Ny/2,:) = REAL(phi(ikx_0,Ny/2,:)) !origin must be real ENDIF !**** Cancel previous moments initialization moments_e = 0._dp; moments_i = 0._dp IF(INIT_ZF .GT. 0) THEN !**** Zonal Flow initialization ******************************************* ! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0 IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN DO iz = izs,ize phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI) moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI) moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp ENDDO ENDIF ENDIF ELSE ! we compute phi from noisy moments and poisson CALL poisson ENDIF END SUBROUTINE init_phi !******************************************************************************! !******************************************************************************! !!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin !******************************************************************************! SUBROUTINE build_dnjs_table USE basic USE array, Only : dnjs USE grid, Only : jmaxe, jmaxi 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 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 IMPLICIT NONE REAL(dp) :: factj, j_dp, j_int REAL(dp) :: be_2, bi_2, alphaD REAL(dp) :: kx, ky, kperp2 !!!!! Electron kernels !!!!! !Precompute species dependant factors factj = 1.0 ! Start of the recursive factorial DO ij = 1, jmaxe+1 j_int = jarray_e(ij) j_dp = REAL(j_int,dp) ! REAL of degree ! Recursive factorial IF (j_dp .GT. 0) THEN factj = factj * j_dp ELSE factj = 1._dp ENDIF DO ikx = ikxs,ikxe kx = kxarray(ikx) DO iky = ikys,ikye ky = kyarray(iky) be_2 = (kx**2 + ky**2) * sigmae2_taue_o2 kernel_e(ij, ikx, iky) = be_2**j_int * exp(-be_2)/factj ENDDO ENDDO ENDDO ! Kernels closure DO ikx = ikxs,ikxe kx = kxarray(ikx) DO iky = ikys,ikye ky = kyarray(iky) be_2 = (kx**2 + ky**2) * sigmae2_taue_o2 - kernel_e(ijeg_e,ikx,iky) = be_2/(real(ijeg_e,dp))*kernel_e(ije_e,ikx,iky) + ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1) + kernel_e(ijeg_e,ikx,iky) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky) + ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0) + IF ( be_2 .NE. 0 ) THEN + kernel_e(ijsg_e,ikx,iky) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky) + ELSE + kernel_e(ijsg_e,ikx,iky) = 0._dp + ENDIF ENDDO ENDDO !!!!! Ion kernels !!!!! factj = 1.0 ! Start of the recursive factorial DO ij = 1, jmaxi+1 j_int = jarray_e(ij) j_dp = REAL(j_int,dp) ! REAL of degree ! Recursive factorial IF (j_dp .GT. 0) THEN factj = factj * j_dp ELSE factj = 1._dp ENDIF DO ikx = ikxs,ikxe kx = kxarray(ikx) DO iky = ikys,ikye ky = kyarray(iky) - bi_2 = (kx**2 + ky**2) * sigmai2_taui_o2 - kernel_i(ij, ikx, iky) = bi_2**j_int * exp(-bi_2)/factj ENDDO ENDDO ENDDO ! Kernels closure DO ikx = ikxs,ikxe kx = kxarray(ikx) DO iky = ikys,ikye ky = kyarray(iky) bi_2 = (kx**2 + ky**2) * sigmai2_taui_o2 - kernel_i(ijeg_i,ikx,iky) = bi_2/(real(ijeg_i,dp))*kernel_e(ije_i,ikx,iky) + ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj + kernel_i(ijeg_i,ikx,iky) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky) + ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0) + IF ( bi_2 .NE. 0 ) THEN + kernel_i(ijsg_i,ikx,iky) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky) + ELSE + kernel_i(ijsg_i,ikx,iky) = 0._dp + ENDIF ENDDO ENDDO END SUBROUTINE evaluate_kernels !******************************************************************************!