diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 9c1ddef5463607f3390369e159b8c0c50a432fbf..ab6343cdf1c52a5b69cf7a2bd26f4af45094af4a 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -1,44 +1,56 @@
module closure
! Contains the routines to define closures
USE basic
-USE model, ONLY: CLOS
+USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, eta_B, nu
USE grid
USE array, ONLY: kernel_e, kernel_i
USE fields, ONLY: moments_e, moments_i
-
+USE time_integration, ONLY: updatetlevel
IMPLICIT NONE
PUBLIC :: apply_closure_model
CONTAINS
- SUBROUTINE apply_closure_model
- IMPLICIT NONE
-
- ! Negative out of bounds indices are put to zero (analytically correct)
+SUBROUTINE apply_closure_model
+ IMPLICIT NONE
+ complex(dp) :: i_kz
+ real(dp) :: taue_qe_etaB_nu, taui_qi_etaB_nu
+ real(dp) :: sqpp2pp1_e, sqpp2pp1_i, sqpp1p_e, sqpp1p_i
+ real(dp) :: p_dp, j_dp
+ ! Spare some computations
+ taue_qe_etaB_nu = tau_e*eta_B/q_e/nu
+ taui_qi_etaB_nu = tau_i*eta_B/q_i/nu
+ sqpp2pp1_e = SQRT((pmaxe_dp+2)*(pmaxe_dp+1))
+ sqpp2pp1_i = SQRT((pmaxi_dp+2)*(pmaxi_dp+1))
+ sqpp1p_e = SQRT((pmaxe_dp+1)*(pmaxe_dp))
+ sqpp1p_i = SQRT((pmaxi_dp+1)*(pmaxi_dp))
+
+ ! Negative out of bounds indices are put to zero (analytically correct)
DO ikr = ikrs,ikre
DO ikz = ikzs,ikze
DO ip = ipsg_e,ipeg_e
- moments_e(ip,ijs_e-1,ikr,ikz,:) = 0._dp
+ moments_e(ip,ijsg_e,ikr,ikz,:) = 0._dp
ENDDO
DO ij = ijsg_e,ijeg_e
- moments_e(ips_e-1,ij,ikr,ikz,:) = 0._dp
- moments_e(ips_e-2,ij,ikr,ikz,:) = 0._dp
+ moments_e(ipsg_e+1,ij,ikr,ikz,:) = 0._dp
+ moments_e(ipsg_e ,ij,ikr,ikz,:) = 0._dp
ENDDO
- kernel_e(ijs_e-1,ikr,ikz) = 0._dp
+ kernel_e(ijsg_e,ikr,ikz) = 0._dp
DO ip = ipsg_i,ipeg_i
- moments_i(ip,ijs_i-1,ikr,ikz,:) = 0._dp
+ moments_i(ip,ijsg_i,ikr,ikz,:) = 0._dp
ENDDO
DO ij = ijsg_i,ijeg_i
- moments_i(ips_i-1,ij,ikr,ikz,:) = 0._dp
- moments_i(ips_i-2,ij,ikr,ikz,:) = 0._dp
+ moments_i(ipsg_i+1,ij,ikr,ikz,:) = 0._dp
+ moments_i(ipsg_i ,ij,ikr,ikz,:) = 0._dp
ENDDO
- kernel_i(ijs_i-1,ikr,ikz) = 0._dp
+ kernel_i(ijsg_i,ikr,ikz) = 0._dp
ENDDO
ENDDO
+
! Positive Oob indices are approximated with a model
IF (CLOS .EQ. 0) THEN
! zero truncation, An+1=0 for n+1>nmax
@@ -46,52 +58,113 @@ CONTAINS
DO ikz = ikzs,ikze
DO ip = ipsg_e,ipeg_e
- moments_e(ip,ije_e+1,ikr,ikz,:) = 0._dp
+ moments_e(ip,ijeg_e,ikr,ikz,:) = 0._dp
ENDDO
DO ij = ijsg_e,ijeg_e
- moments_e(ipe_e+1,ij,ikr,ikz,:) = 0._dp
- moments_e(ipe_e+2,ij,ikr,ikz,:) = 0._dp
+ moments_e(ipeg_e-1,ij,ikr,ikz,:) = 0._dp
+ moments_e(ipeg_e ,ij,ikr,ikz,:) = 0._dp
ENDDO
- kernel_e(ije_e+1,ikr,ikz) = 0._dp
+ kernel_e(ijeg_e,ikr,ikz) = 0._dp
DO ip = ipsg_i,ipeg_i
- moments_i(ip,ije_i+1,ikr,ikz,:) = 0._dp
+ moments_i(ip,ijeg_i,ikr,ikz,:) = 0._dp
ENDDO
DO ij = ijsg_i,ijeg_i
- moments_i(ipe_i+1,ij,ikr,ikz,:) = 0._dp
- moments_i(ipe_i+2,ij,ikr,ikz,:) = 0._dp
+ moments_i(ipeg_i-1,ij,ikr,ikz,:) = 0._dp
+ moments_i(ipeg_i ,ij,ikr,ikz,:) = 0._dp
ENDDO
- kernel_i(ije_i+1,ikr,ikz) = 0._dp
+ kernel_i(ijeg_i,ikr,ikz) = 0._dp
ENDDO
ENDDO
- ELSEIF (CLOS .EQ. 1) THEN
- ! Copy truncation with n+1 = min(nmax,n+1)
- ! here pmax+1 and pmax_2 are mapped to pmax
- moments_e(ipe_e+1,:,:,:,:) = moments_e(ipe_e,:,:,:,:)
- moments_e(ipe_e+2,:,:,:,:) = moments_e(ipe_e,:,:,:,:)
- moments_e(:,ije_e+1,:,:,:) = moments_e(:,ije_e,:,:,:)
- kernel_e(ije_e+1,:,:) = kernel_e(ije_e,:,:)
-
- moments_i(ipe_i+1,:,:,:,:) = moments_i(ipe_i,:,:,:,:)
- moments_i(ipe_i+2,:,:,:,:) = moments_i(ipe_i,:,:,:,:)
- moments_i(:,ije_i+1,:,:,:) = moments_i(:,ije_i,:,:,:)
- kernel_i(ije_i+1,:,:) = kernel_i(ije_i,:,:)
+ ELSEIF ((CLOS .EQ. 1) .AND. (nu .NE. 0)) THEN
+ !Semi collisional closure, i.e. at high degree, -nu N_M+1 ~ X_lin_M * N_M
+ DO ikz = ikzs,ikze
+ i_kz = imagu*kzarray(ikz)
+ !! ELECTRONS
+ ! Hermite closures
+ DO ij = ijsg_e,ijeg_e
+ j_dp = real(jarray_e(ij),dp)
+ DO ikr = ikrs,ikre
+ ! For p = Pmax + 2
+ moments_e(ipeg_e,ij,ikr,ikz,:) = &
+ -taue_qe_etaB_nu * i_kz * sqpp2pp1_e/(2*(pmaxe_dp+2)+j_dp) &
+ * moments_e(ipe_e,ij,ikr,ikz,:)
+ ! For p = Pmax + 1
+ moments_e(ipeg_e-1,ij,ikr,ikz,:) = &
+ -taue_qe_etaB_nu * i_kz * sqpp1p_e/(2*(pmaxe_dp+1)+j_dp) &
+ * moments_e(ipe_e-1,ij,ikr,ikz,:)
+ ! Kernel closure
+ ENDDO
+ ENDDO
+ ! Laguerre closure
+ DO ip = ipsg_e,ipeg_e
+ p_dp = real(parray_e(ip),dp)
+ DO ikr = ikrs,ikre
+ ! For j = Jmax + 1
+ moments_e(ip,ijeg_e,ikr,ikz,:) = &
+ +taue_qe_etaB_nu * i_kz * (jmaxe_dp+1)/(2*p_dp+jmaxe_dp+1) &
+ * moments_e(ip,ije_e,ikr,ikz,:)
+ ENDDO
+ ENDDO
+ !! IONS
+ ! Hermite closures
+ DO ij = ijsg_i,ijeg_i
+ j_dp = real(jarray_i(ij),dp)
+ DO ikr = ikrs,ikre
+ ! For p = Pmax + 2
+ moments_i(ipeg_i,ij,ikr,ikz,:) = &
+ -taui_qi_etaB_nu * i_kz * sqpp2pp1_i/(2*(pmaxi_dp+2)+j_dp) &
+ * moments_i(ipe_i,ij,ikr,ikz,:)
+ ! For p = Pmax + 1
+ moments_i(ipeg_i-1,ij,ikr,ikz,:) = &
+ -taui_qi_etaB_nu * i_kz * sqpp1p_i/(2*(pmaxi_dp+1)+j_dp) &
+ * moments_i(ipe_i-1,ij,ikr,ikz,:)
+ ENDDO
+ ENDDO
+ ! Laguerre closure
+ DO ip = ipsg_i,ipeg_i
+ p_dp = real(parray_i(ip),dp)
+ DO ikr = ikrs,ikre
+ ! For j = Jmax + 1
+ moments_i(ip,ijeg_i,ikr,ikz,:) = &
+ +taui_qi_etaB_nu * i_kz * (jmaxi_dp+1)/(2*p_dp+jmaxi_dp+1) &
+ * moments_i(ip,ije_i,ikr,ikz,:)
+ ENDDO
+ ENDDO
+ ENDDO
ELSEIF (CLOS .EQ. 2) THEN
- ! Copy truncation with special treatment for Hermite
- ! here pmax+1 is mapped to pmax-1 and pmax+2 to pmax
- moments_e(ipe_e+1,:,:,:,:) = moments_e(ipe_e-1,:,:,:,:)
- moments_e(ipe_e+2,:,:,:,:) = moments_e(ipe_e,:,:,:,:)
- moments_e(:,ije_e+1,:,:,:) = moments_e(:,ije_e,:,:,:)
- kernel_e(ije_e+1,:,:) = kernel_e(ije_e,:,:)
-
- moments_i(ipe_i+1,:,:,:,:) = moments_i(ipe_i-1,:,:,:,:)
- moments_i(ipe_i+2,:,:,:,:) = moments_i(ipe_i,:,:,:,:)
- moments_i(:,ije_i+1,:,:,:) = moments_i(:,ije_i,:,:,:)
- kernel_i(ije_i+1,:,:) = kernel_i(ije_i,:,:)
+ ! Copy closure : P+2 <- P, P+1 <- P-1, J+1 <- J
+ DO ikr = ikrs,ikre
+ DO ikz = ikzs,ikze
+
+ DO ip = ipsg_e,ipeg_e
+ ! J ghost is J+1, so we put moment J to J+1
+ moments_e(ip,ijeg_e,ikr,ikz,:) = moments_e(ip,ije_e,ikr,ikz,:)
+ ENDDO
+ DO ij = ijsg_e,ijeg_e
+ ! P ghosts are P+1 and P+2, P+1 <- P-1 and P+2 <- P
+ moments_e(ipeg_e-1,ij,ikr,ikz,:) = moments_e(ipe_e-1,ij,ikr,ikz,:)
+ moments_e(ipeg_e ,ij,ikr,ikz,:) = moments_e(ipe_e ,ij,ikr,ikz,:)
+ ENDDO
+ ! Same for ions
+ DO ip = ipsg_i,ipeg_i
+ moments_i(ip,ijeg_i,ikr,ikz,:) = moments_i(ip,ije_i,ikr,ikz,:)
+ ENDDO
+ DO ij = ijsg_i,ijeg_i
+ moments_i(ipeg_i-1,ij,ikr,ikz,:) = moments_i(ipe_i-1,ij,ikr,ikz,:)
+ moments_i(ipeg_i ,ij,ikr,ikz,:) = moments_i(ipe_i ,ij,ikr,ikz,:)
+ ENDDO
+
+ ENDDO
+ ENDDO
+
+ ELSE
+ if(my_id .EQ. 0) write(*,*) '! Closure scheme not found !'
+
ENDIF
END SUBROUTINE apply_closure_model