diff --git a/src/ei_colls_mod.f90 b/src/ei_colls_mod.f90
index 9b5e0b2..bfb7c90 100644
--- a/src/ei_colls_mod.f90
+++ b/src/ei_colls_mod.f90
@@ -1,425 +1,425 @@
MODULE ei_colls
! electron-ion collisions module
USE PREC_CONST
USE MODEL
USE BASIC
USE BASIC_MPI
use grid
USE ARRAY, only : CeiTdp,CeiFdp,rCeiT,rCeiF
USE OPERATOR_MODEL
USE OUTPUT
use auxval
USE restart
use gkcoulomb
IMPLICIT NONE
PROCEDURE(eval_CeiT), POINTER :: evaluate_CeiT_collisional_matrix => null()
PROCEDURE(eval_CeiF), POINTER :: evaluate_CeiF_collisional_matrix => null()
ABSTRACT INTERFACE
SUBROUTINE eval_CeiT
! Empty
END SUBROUTINE eval_CeiT
END INTERFACE
ABSTRACT INTERFACE
SUBROUTINE eval_CeiF
! Empty
END SUBROUTINE eval_CeiF
END INTERFACE
CONTAINS
SUBROUTINE ei_colls_control
! Set electron-ions basic/model
IMPLICIT NONE
!
IF(me_e .eq. 0) WRITE(*,*) '===== Electron-Ion Collisions ====='
IF(me_e .eq. 0) WRITE(*,*) 'Set Electron-Ion Collision ...'
CALL set_ei_colls
IF( ETEST .ne. 0) THEN
IF(me_e .eq. 0) WRITE(*,*) 'Construct collisional matrices CeiT ... '
CALL evaluate_auxval_test
! Restart
IF( ifrestart) THEN
IF(me .eq. 0 ) WRITE(*,*) 'Restart ...'
CALL restart_init('eiT')
IF(me .eq. 0 ) WRITE(*,*) '... restart done'
ENDIF
call set_CeiT_collisional_matrix
CALL evaluate_CeiT_collisional_matrix
WRITE(*,*) '... CeiT constructed, me =',me_e
ENDIF
IF( EBACK .ne. 0) THEN
IF(me_e .eq. 0) WRITE(*,*) 'Construct collisional matrices CeiF ... '
call evaluate_auxval_field
IF( ifrestart) THEN
IF(me .eq. 0 ) WRITE(*,*) 'Restart ...'
CALL restart_init('eiF')
IF(me .eq. 0 ) WRITE(*,*) '... restart done'
ENDIF
call set_CeiF_collisional_matrix
CALL evaluate_CeiF_collisional_matrix
WRITE(*,*) '... CeiF constructed, me =',me_e
ENDIF
!
END SUBROUTINE ei_colls_control
!
!------------------------------------------------
!
SUBROUTINE set_ei_colls
! Set electron-ions basic/model
IMPLICIT NONE
!
CALL set_model('ei')
CALL set_basic('ei')
!
END SUBROUTINE set_ei_colls
!
!------------------------------------------------
!
! TEST COMPONENT
!
SUBROUTINE set_CeiT_collisional_matrix
! Evaluate the electron-ion test collisional matrix
IMPLICIT NONE
!
!
IF( ETEST .ne. 0 ) THEN
IF(IFCOULOMB .and. GKE .and. ETEST .eq. 1) THEN
evaluate_CeiT_collisional_matrix => CeiT_gk_coulomb
call set_gkcoulomb
!
IF( .not. GKC_USE_FULL) THEN
call load_gk_coulomb_coeffs('eiT')
ELSE
call set_tabljpmf('t')
ENDIF
!
ELSE
!
evaluate_CeiT_collisional_matrix => CeiT_dflt
!
ENDIF
ELSE
evaluate_CeiT_collisional_matrix => CeiT_zeros_dflt
ENDIF
!
END SUBROUTINE set_CeiT_collisional_matrix
!
!----------------------------------------------------
!
SUBROUTINE CeiT_dflt
! Evaluate the electron-ion collisional matrix
! Computes the test collisional matrix.
!
!
IMPLICIT NONE
!
! Local variables
INTEGER :: irow,irow_tmp
INTEGER, DIMENSION(2) :: lk
INTEGER :: istep, lbar_start
! Allocate restart coll mat and restart indices
IF(restart_from) CALL get_restart(CeiTdp)
lbar_start = lbare_s_l
IF(restart_from) lbar_start = row_rst_l +1
istep=1
!
! Test-part
IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CeiT [x/ ', lbare_e_l ,']:'
!
IF (ETEST .ne. 0) THEN
DO irow=lbar_start,lbare_e_l
! Get Hermite-Laguerre indices
irow_tmp = rows_e(irow)
lk = invnumbe(irow_tmp)
rCeiT(:)= TO_FM(0)
write(nscr,'(i5)', advance = 'no') irow_tmp
flush(nscr)
CALL COLLISIONTE(lk(1),lk(2),rCeiT)
CeiTdp(irow,:) = TO_DP(rCeiT(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ei) THEN
call restart_out(CeiTdp, irow)
istep = 1
ENDIF
ENDDO
ENDIF
!
IF(ifrestart) call restart_out(CeiTdp,irow)
!
END SUBROUTINE CeiT_dflt
!
!--------------------------------------------------------------
!
SUBROUTINE CeiT_gk_coulomb
! Evaluate electron-ion test GK Coulomb
use restart
use gkcoulomb
implicit none
INTEGER :: l,k,irow
INTEGER :: j,n,np,istep
INTEGER :: lbar_start
INTEGER, DIMENSION(4) :: lbarjnnp
INTEGER, DIMENSION(2) :: lk
!
!
! restart
IF(restart_from) CALL get_restart( CeiTdp )
!
! initialize indices
lbar_start = lbare_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbare_s_l
ENDIF
istep=1
!
- IF( me .eq. 0 ) write(*,*) 'CeiT [x / ', lbare_e_l,'] :'
+ IF( me .eq. 0 ) write(*,*) 'CeiT [x / ', numbe(Pmaxe,Jmaxe),'] :'
!
DO irow=lbar_start,lbare_e_l
lbarjnnp = rows_e_4dim(irow,1:4)
! Hermite-Laguerre indices
lk = invnumbe(lbarjnnp(1))
j = lbarjnnp(2)
n = lbarjnnp(3)
np = lbarjnnp(4)
write(nscr,'(i5)', advance = 'no') lbarjnnp(1)
flush(nscr)
rCeiT(:) = ZEROFM
CALL coulombTGK(lk(1),lk(2),j,n,np,rCeiT)
CeiTdp(irow,:) = TO_DP(rCeiT(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ei) THEN
! flush current state every nsave_ei.
call restart_out(CeiTdp, irow)
istep = 1
ENDIF
ENDDO
! flush out the final state
IF(ifrestart) call restart_out(CeiTdp, irow)
!
END SUBROUTINE CeiT_gk_Coulomb
!
!------------------------------------------------
!
SUBROUTINE CeiT_zeros_dflt
! Evaluate the electron-ion collisional matrix
! Computes the test collisional matrix.
!
!
IMPLICIT NONE
!
! Local variables
INTEGER :: irow,irow_tmp
INTEGER, DIMENSION(2) :: lk
!
- IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CeiT [x/ ', lbare_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CeiT [x/ ', numbe(Pmaxe,Jmaxe) ,']:'
! Test-part
IF (ETEST .ne. 0) THEN
DO irow=lbare_s_l,lbare_e_l
! Get Hermite-Laguerre indices
irow_tmp = rows_e(irow)
lk = invnumbe(irow_tmp)
rCeiT(:)= TO_FM(0)
write(nscr,'(i5)', advance = 'no') irow
flush(nscr)
CeiTdp(irow,:) = TO_DP(rCeiT(:))
ENDDO
ENDIF
!
!
END SUBROUTINE CeiT_zeros_dflt
!
!--------------------------------------------------------------
!
! FIELD COMPONENT
!
SUBROUTINE set_CeiF_collisional_matrix
! Evaluate the electron-ion collisional matrix
IMPLICIT NONE
!
!
IF (EBACK .ne. 0 ) THEN
IF(IFCOULOMB .and. GKE .and. EBACK .eq. 1) THEN
evaluate_CeiF_collisional_matrix => CeiF_gk_coulomb
!
call set_gkcoulomb
!
IF (.not. GKC_USE_FULL) THEN
call load_gk_coulomb_coeffs('eiF')
ELSE
call set_tabljpmf('f')
ENDIF
!
ELSE
evaluate_CeiF_collisional_matrix =>CeiF_dflt
!
ENDIF
ELSE
evaluate_CeiF_collisional_matrix => CeiF_zeros_dflt
ENDIF
!
END SUBROUTINE set_CeiF_collisional_matrix
!
!----------------------------------------------------
SUBROUTINE CeiF_dflt
! Evaluate the electron-ion collisional matrix
! Computes the back-reaction collisional matrix.
!
!
IMPLICIT NONE
!
INTEGER :: irow,irow_tmp, lbar_start, istep
INTEGER, DIMENSION(2) :: lk
! Allocate restart coll mat and restart indices
IF(restart_from) CALL get_restart(CeiFdp)
lbar_start = lbare_s_l
IF(restart_from) lbar_start = row_rst_l +1
istep=1
IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CeiF [x/ ', lbare_e_l ,']:'
! Back-Reaction
IF (EBACK .ne. 0) THEN
DO irow=lbar_start,lbare_e_l
! Get Hermite-Laguerre indices
irow_tmp = rows_e(irow)
lk = invnumbe(irow_tmp)
write(nscr,'(i5)', advance = 'no') irow
flush(nscr)
rCeiF(:) = TO_FM(0)
CALL COLLISIONFE(lk(1),lk(2),rCeiF)
CeiFdp(irow,:) = TO_DP(rCeiF(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ei) THEN
call restart_out(CeiFdp, irow)
istep = 1
ENDIF
ENDDO
ENDIF
!
!
IF(ifrestart) call restart_out(CeiFdp, irow)
!
END SUBROUTINE CeiF_dflt
!
!--------------------------------------------
!
SUBROUTINE CeiF_gk_coulomb
! Evaluate electron-ion field GK Coulomb
use restart
use gkcoulomb
implicit none
INTEGER :: l,k,irow
INTEGER :: j,n,np,istep
INTEGER :: lbar_start
INTEGER, DIMENSION(4) :: lbarjnnp
INTEGER, DIMENSION(2) :: lk
!
!
! restart
IF(restart_from) CALL get_restart(CeiFdp)
!
! initialize indices
lbar_start = lbare_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbare_s_l
ENDIF
istep=1
!
- IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CeiF [x/ ', lbare_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CeiF [x/ ', numbe(Pmaxe,Jmaxe) ,']:'
!
DO irow=lbar_start,lbare_e_l
lbarjnnp = rows_e_4dim(irow,1:4)
! Hermite-Laguerre indices
lk = invnumbe(lbarjnnp(1))
j = lbarjnnp(2)
n = lbarjnnp(3)
np = lbarjnnp(4)
write(nscr,'(i5)', advance = 'no') lbarjnnp(1)
flush(nscr)
rCeiF(:) = ZEROFM
CALL coulombFGK(lk(1),lk(2),j,n,np,rCeiF)
CeiFdp(irow,:) = TO_DP(rCeiF(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ei) THEN
! flush current state every nsave_ei.
call restart_out(CeiFdp, irow)
istep = 1
ENDIF
ENDDO
! flush out the final state
IF(ifrestart) call restart_out(CeiFdp, irow)
!
END SUBROUTINE CeiF_gk_Coulomb
!
!--------------------------------------------------
!
SUBROUTINE CeiF_zeros_dflt
! Evaluate the electron-ion collisional matrix
! Computes the back-reaction collisional matrix.
!
!
IMPLICIT NONE
!
INTEGER :: irow,irow_tmp
INTEGER, DIMENSION(2) :: lk
! Back-Reaction
IF (EBACK .ne. 0) THEN
DO irow=lbare_s_l,lbare_e_l
! Get Hermite-Laguerre indices
irow_tmp = rows_e(irow)
lk = invnumbe(irow_tmp)
write(*,*) 'CeiF: ',lk(1),'/',Pmaxe,' ,',lk(2),'/',Jmaxe
rCeiF(:) = TO_FM(0)
CeiFdp(irow,:) = TO_DP(rCeiF(:))
ENDDO
ENDIF
!
!
END SUBROUTINE CeiF_zeros_dflt
!
END MODULE ei_colls
diff --git a/src/ie_colls_mod.f90 b/src/ie_colls_mod.f90
index e464d32..61d085d 100644
--- a/src/ie_colls_mod.f90
+++ b/src/ie_colls_mod.f90
@@ -1,438 +1,438 @@
module ie_colls
! ion-electron collisions module
USE PREC_CONST
USE MODEL
use grid
USE BASIC
USE BASIC_MPI
USE restart
USE ARRAY, only : CieTdp,CieFdp,rCieT,rCieF
USE OPERATOR_MODEL
use auxval
use gkcoulomb
implicit none
PROCEDURE(CieT_sub), POINTER :: evaluate_CieT_collisional_matrix => null()
PROCEDURE(CieF_sub), POINTER :: evaluate_CieF_collisional_matrix => null()
ABSTRACT INTERFACE
SUBROUTINE CieT_sub
! Empty
END SUBROUTINE CieT_sub
END INTERFACE
ABSTRACT INTERFACE
SUBROUTINE CieF_sub
! Empty
END SUBROUTINE CieF_sub
END INTERFACE
contains
!
!------------------------------------------------
!
subroutine ie_colls_control
! ie colls
implicit none
IF(me_i .eq. 0) WRITE(*,*) '====== Ion-Electron Collisions ======='
IF(me_i .eq. 0) WRITE(*,*) 'Set Ion-Electron Collision ...'
CALL set_ie_colls
IF( ITEST .ne. 0) THEN
CALL evaluate_auxval_test
IF(me_i .eq. 0) WRITE(*,*) 'Construct collisional matrices CieT ... '
! Restart
IF( ifrestart) THEN
IF(me .eq. 0 ) WRITE(*,*) 'Restart ...'
CALL restart_init('ieT')
IF(me .eq. 0 ) WRITE(*,*) '... restart done'
ENDIF
call set_CieT_collisional_matrix
CALL evaluate_CieT_collisional_matrix
WRITE(*,*) '... CieT constructed, me =', me_i
ENDIF
IF( IBACK .ne. 0) THEN
IF(me_i .eq. 0) WRITE(*,*) 'Construct collisional matrices CieF ... '
call evaluate_auxval_field
! Restart
IF( ifrestart) THEN
IF(me .eq. 0 ) WRITE(*,*) 'Restart ...'
CALL restart_init('ieF')
IF(me .eq. 0 ) WRITE(*,*) '... restart done'
ENDIF
call set_CieF_collisional_matrix
CALL evaluate_CieF_collisional_matrix
WRITE(*,*) '... CieF constructed, me =', me_i
ENDIF
end subroutine ie_colls_control
!
!------------------------------------------------
!
SUBROUTINE set_ie_colls
! Set ion-electron basic/model
IMPLICIT NONE
!
CALL set_model('ie')
!
CALL set_basic('ie')
!
END SUBROUTINE set_ie_colls
!
!------------------------------------------------
!
! TEST COMPONENT
!
SUBROUTINE set_CieT_collisional_matrix
! Evaluate the electron-ion test collisional matrix
IMPLICIT NONE
!
!
IF( ITEST .ne. 0 ) THEN
IF(IFCOULOMB .and. GKI .and. ITEST .eq. 1) THEN
!
evaluate_CieT_collisional_matrix => CieT_gk_coulomb
call set_gkcoulomb
!
IF ( .not. GKC_USE_FULL) then
call load_gk_coulomb_coeffs('ieT')
else
call set_tabljpmf('t')
endif
!
ELSE
!
evaluate_CieT_collisional_matrix => CieT_dflt
!
ENDIF
ELSE
evaluate_CieT_collisional_matrix => CieT_zeros_dflt
ENDIF
!
END SUBROUTINE set_CieT_collisional_matrix
!
!----------------------------------------------------
!
SUBROUTINE CieT_dflt
! Evaluate the ion-electron collisional matrix
! Computes the test collisional matrix.
!
!
IMPLICIT NONE
!
INTEGER :: irow,irow_tmp, lbar_start, istep
INTEGER, DIMENSION(2) :: lk
! restart
IF(restart_from) CALL get_restart(CieTdp)
!
! initialize indices
lbar_start = lbari_s_l
IF(restart_from) lbar_start = row_rst_l +1
istep=1
!
- IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieT [x/ ', lbari_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieT [x/ ', numbi(Pmaxi,Jmaxi) ,']:'
! Test-part
IF( ITEST .ne. 0) THEN
DO irow=lbar_start,lbari_e_l
! Get Hermite-Laguerre indices
irow_tmp = rows_i(irow)
lk = invnumbi(irow_tmp)
write(nscr,'(i5)', advance = 'no') irow_tmp
flush(nscr)
rCieT(:) = TO_FM('0')
CALL COLLISIONTI(lk(1),lk(2),rCieT)
CieTdp(irow,:) = TO_DP(rCieT(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ie) THEN
! flush current state every nsave_ie.
call restart_out(CieTdp,irow)
istep = 1
ENDIF
ENDDO
ENDIF
!
! flush out the final state
IF(ifrestart) call restart_out(CieTdp,irow)
END SUBROUTINE CieT_dflt
!
!--------------------------------------------------------------
!
SUBROUTINE CieT_gk_coulomb
! Evaluate ion-electron test GK Coulomb
use restart
use gkcoulomb
implicit none
INTEGER :: l,k,irow
INTEGER :: j,n,np,istep
INTEGER :: lbar_start
INTEGER, DIMENSION(4) :: lbarjnnp
INTEGER, DIMENSION(2) :: lk
!
!
! restart
IF(restart_from) CALL get_restart(CieTdp)
!
! initialize indices
lbar_start = lbari_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbari_s_l
ENDIF
istep=1
!
- IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieT [x/ ', lbari_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieT [x/ ', numbi(Pmaxi,Jmaxi) ,']:'
DO irow=lbar_start,lbari_e_l
lbarjnnp = rows_i_4dim(irow,1:4)
! Hermite-Laguerre indices
lk = invnumbi(lbarjnnp(1))
j = lbarjnnp(2)
n = lbarjnnp(3)
np = lbarjnnp(4)
write(nscr,'(i5)', advance = 'no') lbarjnnp(1)
flush(nscr)
rCieT(:) = ZEROFM
CALL coulombTGK(lk(1), lk(2),j, n, np, rCieT)
CieTdp(irow,:) = TO_DP(rCieT(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ie) THEN
! flush current state every nsave_ie.
call restart_out(CieTdp, irow)
istep = 1
ENDIF
ENDDO
! flush out the final state
IF(ifrestart) call restart_out(CieTdp, irow )
!
END SUBROUTINE CieT_gk_Coulomb
!
!------------------------------------------------
!
SUBROUTINE CieT_zeros_dflt
! Evaluate the ion-electron collisional matrix
! Computes the test collisional matrix.
!
!
IMPLICIT NONE
!
INTEGER :: irow,irow_tmp
INTEGER, DIMENSION(2) :: lk
INTEGER :: lbar_start, istep
! restart
IF(restart_from) CALL get_restart(CieTdp)
!
! initialize indices
lbar_start = lbari_s_l
IF(restart_from) lbar_start = row_rst_l +1
istep=1
! Test-part
IF( ITEST .ne. 0) THEN
DO irow=lbari_s_l,lbari_e_l
! Get Hermite-Laguerre indices
irow_tmp = rows_i(irow)
lk = invnumbi(irow_tmp)
write(*,'(a,1x,i3,a,i3,a,i3,a,i3)') 'CieT:', lk(1),'/',Pmaxi, ', ',lk(2),'/',Jmaxi
rCieT(:) = TO_FM('0')
CieTdp(irow,:) = TO_DP(rCieT(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ie) THEN
! flush current state every nsave_ie.
call restart_out(CieTdp, irow)
istep = 1
ENDIF
ENDDO
ENDIF
!
END SUBROUTINE CieT_zeros_dflt
!
!--------------------------------------------------------------
!
! FIELD COMPONENT
!
SUBROUTINE set_CieF_collisional_matrix
! Evaluate the ion-electron collisional matrix
IMPLICIT NONE
!
IF (IBACK .ne. 0 ) THEN
IF(IFCOULOMB .and. GKI .and. IBACK .eq. 1 ) THEN
!
evaluate_CieF_collisional_matrix => CieF_gk_coulomb
call set_gkcoulomb
IF(.not. GKC_USE_FULL) then
call load_gk_coulomb_coeffs('ieF')
else
call set_tabljpmf('f')
endif
!
ELSE
evaluate_CieF_collisional_matrix => CieF_dflt
!
ENDIF
ELSE
evaluate_CieF_collisional_matrix => CieF_zeros_dflt
ENDIF
!
END SUBROUTINE set_CieF_collisional_matrix
!
!--------------------------------------------------------------
!
SUBROUTINE CieF_dflt
! Evaluate the ion-electron collisional matrix
! Computes the test collisional matrix.
!
!
IMPLICIT NONE
!
! Local variables
INTEGER :: irow,irow_tmp
INTEGER :: lbar_start, istep
INTEGER, DIMENSION(2) :: lk
! restart
IF(restart_from) CALL get_restart(CieFdp)
!
! initialize indices
lbar_start = lbari_s_l
IF(restart_from) lbar_start = row_rst_l +1
istep=1
- IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieF [x/ ', lbari_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieF [x/ ', numbi(Pmaxi,Jmaxi) ,']:'
! Back-Reaction
IF (IBACK .ne. 0) THEN
DO irow = lbari_s_l,lbari_e_l
irow_tmp = rows_i(irow)
lk = invnumbi(irow_tmp)
write(nscr,'(i5)', advance = 'no') irow
flush(nscr)
rCieF(:) = TO_FM('0')
CALL COLLISIONFI(lk(1),lk(2),rCieF)
CieFdp(irow,:) = TO_DP(rCieF(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ie) THEN
! flush current state every nsave_ie.
call restart_out(CieFdp, irow )
istep = 1
ENDIF
ENDDO
ENDIF
! flush out the final state
IF(ifrestart) call restart_out(CieFdp, irow)
!
END SUBROUTINE CieF_dflt
!
!--------------------------------------------------------------
!
SUBROUTINE CieF_gk_coulomb
! Evaluate ion-electron field GK Coulomb
use restart
use gkcoulomb
implicit none
INTEGER :: l,k,irow
INTEGER :: j,n,np,istep
INTEGER :: lbar_start
INTEGER, DIMENSION(4) :: lbarjnnp
INTEGER, DIMENSION(2) :: lk
!
!
! restart
IF(restart_from) CALL get_restart(CieFdp)
!
! initialize indices
lbar_start = lbari_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbari_s_l
ENDIF
istep=1
!
IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'CieF [x/ ', lbari_e_l ,']:'
DO irow=lbar_start,lbari_e_l
lbarjnnp = rows_i_4dim(irow,1:4)
! Hermite-Laguerre indices
lk = invnumbi(lbarjnnp(1))
j = lbarjnnp(2)
n = lbarjnnp(3)
np = lbarjnnp(4)
write(nscr,'(i5)', advance = 'no') lbarjnnp(1)
flush(nscr)
rCieF(:) = ZEROFM
CALL coulombFGK(lk(1),lk(2),j,n,np,rCieF)
CieFdp(irow,:) = TO_DP(rCieF(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ie) THEN
! flush current state every nsave_ie.
call restart_out(CieFdp, irow)
istep = 1
ENDIF
ENDDO
! flush out the final state
IF(ifrestart) call restart_out(CieFdp, irow)
!
END SUBROUTINE CieF_gk_Coulomb
!
!---------------------------------------------------------------
!
SUBROUTINE CieF_zeros_dflt
! Evaluate the ion-electron collisional matrix
! Computes the test collisional matrix.
!
!
IMPLICIT NONE
!
! Local variables
INTEGER :: irow,irow_tmp
INTEGER, DIMENSION(2) :: lk
! Back-Reaction
IF (IBACK .ne. 0) THEN
DO irow = lbari_s_l,lbari_e_l
irow_tmp = rows_i(irow)
lk = invnumbi(irow_tmp)
write(*,*) 'CieF:', lk(1),'/',Pmaxi, ', ',lk(2),'/',Jmaxi
rCieF(:) = TO_FM('0')
CieFdp(irow,:) = TO_DP(rCieF(:))
ENDDO
ENDIF
END SUBROUTINE CieF_zeros_dflt
!
end module ie_colls
diff --git a/src/model_mod.f90 b/src/model_mod.f90
index c50bf5b..06e69e2 100644
--- a/src/model_mod.f90
+++ b/src/model_mod.f90
@@ -1,729 +1,730 @@
module model
! phyical parameters
use prec_const
use FMZM
use basic
IMPLICIT NONE
PRIVATE
! physical parameters
real(dp), public, protected :: nu = 1._dp ! electron-ion collision freq.
real(dp), public, protected :: sigmaei =1._dp ! me/ mi
real(dp), public, protected :: tauie = 1._dp ! T_i0/T_e0
real(dp), public, protected :: kperp = 1._dp ! perpendicular wavenumber normalized to cs
! kperp scan
logical, public, protected :: kpscan = .false. ! if .true., run kperp scan
integer, public, protected :: iks, ike, ik_current
integer, public, protected :: ikp_start = 1
-
+
real(dp), dimension(:), allocatable :: kperp_array
- character(len = 64) :: kscanfile
+ character(len = slen) :: kpscanfile = 'kperp.in'
! electron-ion mass/temperature ratio
real(dp), public, protected :: sigmaaa,sigmaie
real(dp), public, protected :: tauei, tauaa
real(dp), public, protected :: sigma_ = 1._dp !
real(dp), public, protected :: tau_ = 1._dp !
real(dp), public, protected :: chi_ = 1._dp !
real(dp), public, protected :: be_ = 0._dp ! perpendicular wavenumber normalized to electron thermal Larmor radius
real(dp), public, protected :: bi_ = 0._dp ! perpendicular wavenumber normalized to electron thermal Larmor radius
real(dp), public, protected :: ba_ = 0._dp
real(dp), public, protected :: bb_ = 0._dp
! Multiple-Precision FM type
TYPE(FM), public, protected :: sigmafm_ ! mass ratio
TYPE(FM), public, protected :: taufm_ ! temperature ratio
TYPE(FM), public, protected :: chifm_ ! SQRT(tau/sigma) = vTa/vTb
TYPE(FM), public, protected :: thetaabfm_ ! Sugama theta_{ab} coefficient [see Eq. (29) in Sugama et al. (2009)]
TYPE(FM), public, protected :: thetabafm_ ! Sugma theta_{ba} coefficient
TYPE(FM), public, protected :: bafm_
TYPE(FM), public, protected :: bbfm_ !
! Test/Back-Reaction Model Collision Operator
! electron collision model
CHARACTER(slen), protected ,public :: etest_coll = 'none'
CHARACTER(slen), protected ,public :: eback_coll = 'none'
CHARACTER(slen), protected ,public :: eself_coll = 'none'
! ion collision model
CHARACTER(slen), protected , public :: itest_coll = 'none'
CHARACTER(slen), protected , public :: iback_coll = 'none'
CHARACTER(slen), protected , public:: iself_coll = 'none'
INTEGER, PUBLIC :: ETEST=0 ! Electron Test Model
INTEGER, PUBLIC :: ITEST=0 ! Ion Test Model
INTEGER, PUBLIC :: EBACK=0 ! Electron Back-Reaction Model
INTEGER, PUBLIC :: IBACK=0 ! Ion Back-Reaction Model
INTEGER, PUBLIC :: ESELF=0 ! Electron self Model
INTEGER, PUBLIC :: ISELF=0 ! Ion Self Model
INTEGER, PUBLIC :: GKE =0 ! Gyrokinetic Electrons
INTEGER, PUBLIC :: GKI =0 ! Gyrokinetic Ions
LOGICAL, PUBLIC :: ADD_IS_TO_OS = .False. ! add improved sugama to originial sugama
! Apply long wavelength limit to GK collision operators
LOGICAL, PUBLIC :: LWL = .FALSE.
! routine to compute gk coulomb
LOGICAL , PUBLIC :: GKC_USE_FULL = .false. ! if true, use explicit routine, false, use precomputed quantities
LOGICAL , PUBLIC :: GKC_USE_OPT = .false. ! if true, use optimized routine, false, use precomputed quantities without optimized version
! Colliding species
CHARACTER(len=2), PUBLIC, PROTECTED :: coll_ab
! Sugama Flag: used in AUXVAL for error functions necessary for back-reaction parts
LOGICAL, PUBLIC, PROTECTED :: IFSUGAMA = .FALSE.
! Improved Sugama Flag: used in AUXVAL and MEMORY for precomputation of CurlI(para/perp)(a/b)pjk in GK mode
LOGICAL, PUBLIC, PROTECTED :: IFIMPROVEDSUGAMA = .FALSE.
! Coulomb flag
LOGICAL, PUBLIC, PROTECTED :: IFCOULOMB = .FALSE.
! Polarization flag
LOGICAL, PUBLIC, PROTECTED :: IFPOL = .FALSE.
! solve Poisson for polarizion term
TYPE(FM), public :: poisson_coeff ! ion polarization and adiabatic electron coeff.
! Zeroth order Hirshman Sigmar Clarke flag
logical, PUBLIC, protected :: IF0HSC = .false.
! GK Abel flag: used in AUXVAL for error functions necessary for back-reaction parts
logical, PUBLIC, protected :: IFABEL = .false.
! DK/GK test or/and field: used only if GKI and/or GKE
! Note: implented only for Abel collision operator
logical, PUBLIC, protected :: DKBACK = .false. ! use DK back reaction
logical, PUBLIC, protected :: DKTEST = .false. ! use DK test part
! Flag turn on/off Test/Field parts
logical, public, protected :: ADDBACK =.true. ! IF 1 add field part, 0 otherwise
logical, public, protected :: ADDTEST=.true. ! IF 1 add test part, 0 otherwise
! Gyrokinetic flag
LOGICAL, PUBLIC, PROTECTED :: IFGK = .FALSE.
public :: model_readinputs, set_model,set_model_fm,operatormodel_readinput,invert_model, set_model_var, set_kpscan, update_kperp
contains
!================================================================================
subroutine model_readinputs
! Read inputs physical input parameters
use basic, only: lu_in,ierr
use prec_const
implicit none
character(len=8) :: fmt = '(f6.4)'
- NAMELIST /MODEL_PAR/ nu, sigmaei ,tauie,kperp, kpscan, ikp_start
+ NAMELIST /MODEL_PAR/ nu, sigmaei ,tauie,kperp, kpscan, ikp_start, kpscanfile
IF( me .ne. 0) GOTO 111
read(lu_in,MODEL_PAR)
! write(*,MODEL_PAR)
111 CONTINUE
! We compute the normalized collisional matrix
nu = 1._dp
!
! Broadcast from PE 0 to
!
CALL MPI_bcast(nu,1, MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(sigmaei, 1, MPI_DOUBLE_PRECISION, 0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(tauie, 1,MPI_DOUBLE_PRECISION, 0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(kperp,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(kpscan,1,MPI_logical, 0, MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ikp_start,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
+ CALL MPI_bcast(kpscanfile, slen , MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)
!
call set_model_var
!
end subroutine model_readinputs
!================================================================================
subroutine set_model_var
implicit none
! set electron-ion mass/temperature ratio
!! defined as me/mi !!!!!!
sigmaie = 1/sigmaei
sigmaaa = 1._dp
tauei = 1/tauie
tauaa = 1._dp
! Set electron/ion perpendicular wavenumber
be_ = kperp*SQRT(2._dp)*SQRT(sigmaei)
bi_ = kperp*SQRT(2._dp)*SQRT(tauei)
if( .not. kpscan) then
iks = 1
ike = 1
endif
end subroutine set_model_var
!
!================================================================================
!
subroutine set_kpscan
! set kperp scan
implicit none
logical :: exist_
integer :: fid, nk, endof
real(dp) :: buff_
- kscanfile = 'kperp.in'
+ !! kscanfile = 'kperp.in'
IF( me .eq. 0) WRITE(*,*) "Run kperp scan ... "
- INQUIRE(FILE=kscanfile, EXIST=exist_)
+ INQUIRE(FILE=kpscanfile, EXIST=exist_)
IF(.not. exist_) THEN
ERROR STOP("kperp.in file not found... ")
ENDIF
! read scanlist
- open(fid,FILE= kscanfile)
+ open(fid,FILE= kpscanfile)
endof = .false.
nk =0
do while( endof .eq. 0)
read( fid, *, iostat = endof) buff_
nk = nk +1
enddo
nk = nk-1
close(fid)
! allocate array
iks = 1
ike = nk
call allocate_array( kperp_array, iks, ike)
- open(fid, FILE=kscanfile)
+ open(fid, FILE=kpscanfile)
endof = .false.
nk =0
do while( endof .eq. 0)
read( fid, *, iostat = endof) buff_
kperp_array( nk+1) = buff_
if( me .eq. 0 .and. endof .eq. 0 ) WRITE(*,*) 'kperp = ', kperp_array( nk+1)
nk = nk +1
enddo
nk = nk-1
close(fid)
!
end subroutine set_kpscan
!
!================================================================================
!
subroutine update_kperp(ik)
!
implicit none
integer, intent(in) :: ik
! update the current ik index to compute
ik_current = ik
kperp = kperp_array( ik)
IF( me .eq. 0) WRITE(*,*) ' update kperp =', kperp
call set_model_var
end subroutine update_kperp
!
!================================================================================
!
subroutine set_model(ab)
! set the species mass/temperature ratios
implicit none
character(len=2), intent(in) :: ab
! Ion-Ion collisions
if(ab .eq. 'ii') then
chi_ = 1._dp
tau_ = 1._dp
sigma_ = 1._dp
ba_ = bi_
bb_ = bi_
coll_ab = 'ii'
end if
! Electron- Electron collisions
if(ab .eq. 'ee') then
chi_ = 1._dp
tau_ = 1._dp
sigma_ = 1._dp
ba_ = be_
bb_ = be_
coll_ab = 'ee'
end if
! Electron-ion collisions
if(ab .eq. 'ei') then
chi_ = SQRT(tauei/sigmaei)
tau_ = tauei
sigma_ = sigmaei
ba_ = be_
bb_ = bi_
coll_ab = 'ei'
end if
! Ion-electron collisions
if(ab .eq. 'ie') then
chi_ = SQRT(tauie/sigmaie)
tau_ = tauie
sigma_ = sigmaie
ba_ = bi_
bb_ = be_
coll_ab = 'ie'
end if
! set model parameter to FM type
CALL set_model_fm
!
end subroutine set_model
!================================================================================
subroutine set_model_fm
! set mass/temperature ratio to FM type
implicit none
!
CALL FM_ENTER_USER_ROUTINE
! Set temperature and mass ratios to FM type
sigmafm_ = TO_FM(sigma_)
taufm_ = TO_FM(tau_)
chifm_ = SQRT(taufm_/sigmafm_)
! Set perpendicular waves numbers to FM type
bafm_ = TO_FM(ba_)
bbfm_ = TO_FM(bb_)
!
!
! Set Sugama parameter
thetaabfm_ = SQRT(TO_FM( taufm_ + chifm_**2)/(TO_FM(1) + chifm_**2))
thetabafm_ = SQRT(TO_FM( TO_FM(1)/taufm_ + TO_FM(1)/chifm_**2)/(TO_FM(1) + TO_FM(1)/chifm_**2))
!
CALL FM_EXIT_USER_ROUTINE
end subroutine set_model_fm
!================================================================================
subroutine invert_model
! inverse ab --> ba
! Used in the Sugama collision operator
IMPLICIT NONE
!
CALL FM_ENTER_USER_ROUTINE
!
sigmafm_ = 1/sigmafm_
taufm_ = 1/taufm_
chifm_ = 1/chifm_
!
! Set Sugama parameter
thetaabfm_ = SQRT( (TO_FM(1) + sigmafm_)/(TO_FM(1) + sigmafm_/taufm_))
thetabafm_ = SQRT( (TO_FM(1) + TO_FM(1)/sigmafm_)/(TO_FM(1) + taufm_/sigmafm_))
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE invert_model
!================================================================================
SUBROUTINE operatormodel_readinput
! Read Model Operator from input file
USE BASIC
IMPLICIT NONE
! Namelist
! NAMELIST /OPERATOR_MODEL/ ETEST,ITEST,EBACK,IBACK,ESELF,ISELF
! NAMELIST /OPERATOR_MODEL/ DKTEST,DKBACK,ADDBACK,ADDTEST
! NAMELIST /OPERATOR_MODEL/ LWL,ADD_IS_TO_OS, GKC_USE_FULL
NAMELIST /OPERATOR_MODEL/ etest_coll, itest_coll, eback_coll, iback_coll, eself_coll, iself_coll
NAMELIST /OPERATOR_MODEL/ DKTEST,DKBACK,ADDBACK,ADDTEST
NAMELIST /OPERATOR_MODEL/ LWL,ADD_IS_TO_OS, GKC_USE_FULL, GKC_USE_OPT
!
IF(me .ne. 0) GOTO 111
READ(lu_in,OPERATOR_MODEL)
WRITE(*,OPERATOR_MODEL)
!
111 CONTINUE
!
! MPI COMM
call mpi_bcast( etest_coll, slen, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
call mpi_bcast( eback_coll, slen, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
call mpi_bcast( eself_coll, slen, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
call mpi_bcast( itest_coll, slen, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
call mpi_bcast( iback_coll, slen, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
call mpi_bcast( iself_coll, slen, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_bcast(LWL,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(DKTEST,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(DKBACK,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ADDBACK,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ADDTEST,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ADD_IS_TO_OS,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(GKC_USE_FULL, 1, MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(GKC_USE_OPT, 1, MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
!
! Set collision operator flags and different other flags
!
call set_eiT_flag
call set_eiF_flag
call set_ieT_flag
call set_ieF_flag
call set_ee_flag
call set_ii_flag
!
call set_coll_flags
!
END SUBROUTINE operatormodel_readinput
!
!------------------------------------------------------------
!
! set collision operator model flags
! 0 = nothing
! 1 = coulomb
! 2 = pitch-angle (with constant coll.freq.)
! 3 = sugama
! 4 = pitch-angle with veloty dependent coll. freq.
! 5 = improved sugama
! 6 = Hirschman-Sigmar Clarke
! 7 = Abel (only for like species)
! 8 = conserving momentun pitch-angle operator (with velocity
! dependent coll. freq.)
! 9 = GK coulomb polarization term
! Note: check operator_model_mod.f90 routines for details.
subroutine set_eiT_flag
implicit none
SELECT CASE( etest_coll )
! electron test
CASE ('dkcoulomb')
ETEST = 1
GKE = 0
CASE ('gkcoulomb')
ETEST =1
GKE =1
CASE ('dksugama')
ETEST =3
GKE =0
CASE ('gksugama')
ETEST =3
GKE =1
CASE ('dkimpsugama')
ETEST =5
GKE =0
CASE ('gkimpsugama')
ETEST =5
GKE =1
CASE ('dklorentz')
ETEST =2
GKE =0
CASE ('dkpitch')
ETEST =4
GKE = 0
CASE ('gkpitch')
ETEST =4
GKE = 1
CASE ('hsclarke')
ETEST = 6
GKE = 0
CASE ('none')
ETEST = 0
GKE = 0
CASE DEFAULT
ERROR STOP('ETEST_coll not defined ... ')
END SELECT
!
end subroutine set_eiT_flag
!
!------------------------------------------------------------
!
subroutine set_eiF_flag
implicit none
SELECT CASE( eback_coll )
! electron field
CASE ('dkcoulomb')
EBACK = 1
GKE = 0
CASE ('gkcoulomb')
EBACK =1
GKE = 1
CASE ('dksugama')
EBACK =3
GKE =0
CASE ('gksugama')
EBACK =3
GKE =1
CASE ('dkimpsugama')
EBACK =5
GKE =0
CASE ('gkimpsugama')
EBACK =5
GKE =1
CASE ('hsclarke')
EBACK = 6
GKE = 0
CASE ('none')
EBACK = 0
GKE = 0
CASE DEFAULT
ERROR STOP('EBACK_coll not defined ... ')
END SELECT
!
end subroutine set_eiF_flag
!
!------------------------------------------------------------
!
subroutine set_ieT_flag
implicit none
SELECT CASE( itest_coll )
! ion test
CASE ('dkcoulomb')
ITEST = 1
GKI = 0
CASE ('gkcoulomb')
ITEST =1
GKI = 1
CASE ('dksugama')
ITEST =3
GKI = 0
CASE ('gksugama')
ITEST =3
GKI = 1
CASE ('dkimpsugama')
ITEST =5
GKI = 0
CASE ('gkimpsugama')
ITEST =5
GKI = 1
CASE ('hsclarke')
ITEST = 6
GKI = 0
CASE ('none')
ITEST = 0
GKI = 0
CASE DEFAULT
ERROR STOP('ITEST_coll not defined ... ')
END SELECT
!
end subroutine set_ieT_flag
!
!------------------------------------------------------------
!
subroutine set_ieF_flag
implicit none
SELECT CASE( iback_coll )
! ion field
CASE ('dkcoulomb')
IBACK = 1
GKI = 0
CASE ('gkcoulomb')
IBACK =1
GKI = 1
CASE ('dksugama')
IBACK =3
GKI = 0
CASE ('gksugama')
IBACK =3
GKI = 1
CASE ('dkimpsugama')
IBACK =5
GKI = 0
CASE ('gkimpsugama')
IBACK =5
GKI = 1
CASE ('hsclarke')
IBACK = 6
GKI = 0
CASE ('none')
IBACK = 0
GKI = 0
CASE DEFAULT
ERROR STOP('IBACK_coll not defined ... ')
END SELECT
!
end subroutine set_ieF_flag
!
!------------------------------------------------------------
!
subroutine set_ii_flag
implicit none
SELECT CASE( iself_coll )
! ion field
CASE ('dkcoulomb')
ISELF = 1
GKI = 0
CASE ('gkcoulomb')
ISELF =1
GKI = 1
CASE ('dksugama')
ISELF =3
GKI = 0
CASE ('gksugama')
ISELF =3
GKI = 1
CASE ('dkimpsugama')
ISELF =5
GKI = 0
CASE ('gkimpsugama')
ISELF =5
GKI = 1
CASE ('dklorentz')
ISELF = 2
GKI = 0
CASE ('dkpitch_wof')
ISELF = 4
GKI = 0
CASE ('gkpitch_wof')
ISELF =4
GKI = 1
CASE ('dkpitch_wf')
ISELF = 8
GKI = 0
CASE ('gkpitch_wf')
ISELF = 8
GKI = 1
CASE ('hsclarke')
ISELF = 6
GKI = 0
CASE('dkdougherty')
ISELF = 10
GKI =0
CASE ('none')
ISELF = 0
GKI = 0
CASE DEFAULT
ERROR STOP('ISELF_coll not defined ... ')
END SELECT
!
end subroutine set_ii_flag
!
!------------------------------------------------------------
!
subroutine set_ee_flag
implicit none
SELECT CASE( eself_coll )
! ion field
CASE ('dkcoulomb')
ESELF = 1
GKE = 0
CASE ('gkcoulomb')
ESELF =1
GKE = 1
CASE ('dksugama')
ESELF =3
GKE = 0
CASE ('gksugama')
ESELF =3
GKE = 1
CASE ('dkimpsugama')
ESELF =5
GKE = 0
CASE ('gkimpsugama')
ESELF =5
GKE = 1
CASE ('dklorentz')
ESELF = 2
GKE = 0
CASE ('dkpitch_wf')
ESELF = 8
GKE = 0
CASE ('gkpitch_wf')
ESELF = 8
GKE = 1
CASE ('hsclarke')
ESELF = 6
GKE = 0
CASE('dkdougherty')
ESELF = 10
GKE =0
CASE ('none')
ESELF = 0
GKE = 0
CASE DEFAULT
ERROR STOP('ESELF_coll not defined ... ')
END SELECT
!
end subroutine set_ee_flag
!
!-----------------------------------------------------------------------------------------
!
subroutine set_coll_flags
! set coll flags
implicit none
! If Sugama COOP
IF(ETEST .eq. 3 .or. EBACK .eq. 3 .or. ITEST .eq. 3 .or. IBACK .eq. 3 .or. ESELF .eq. 3 .or. ISELF .eq. 3) THEN
IFSUGAMA = .true.
ENDIF
! If Improved Sugama COOP
IF(ETEST .eq. 5 .or. EBACK .eq. 5 .or. ITEST .eq. 5 .or. IBACK .eq. 5 .or. ESELF .eq. 5 .or. ISELF .eq. 5) THEN
IFSUGAMA = .true.
IFIMPROVEDSUGAMA = .true.
imaxx = JEmaxx
ENDIF
! IF zeroth order Hirshman-Sigmar-Clarke (0HSC)
IF(ETEST .eq. 6 .or. EBACK .eq. 6 .or. ITEST .eq. 6 .or. IBACK .eq. 6 .or. ESELF .eq. 6 .or. ISELF .eq. 6) THEN
IF0HSC = .true.
imaxx = JEmaxx
ENDIF
! IF GK abel
IF(ISELF .eq. 7 ) THEN
IFABEL = .true.
ENDIF
! IF GK Coulomb operator
IF((ETEST .eq. 1 .or. EBACK .eq. 1 .or. ITEST .eq. 1 .or. IBACK .eq. 1 .or. ESELF .eq. 1 .or. ISELF .eq. 1) &
.and. ((GKI .eq. 1) .or. (GKE .eq. 1)) ) THEN
IFCOULOMB = .true.
ENDIF
! IF gyrokinetic
IF(GKE .eq. 1 .or. GKI .eq. 1) THEN
IFGK = .true.
ENDIF
!
! IF GK Polarization
IF(ISELF .eq. 9) THEN
IFPOL = .true.
ENDIF
!
IF(LWL) THEN
IF( me .eq. 0) WRITE(*,*) 'Apply long wavelength approximation ...'
ENDIF
end subroutine set_coll_flags
END MODULE model
diff --git a/src/self_colls_mod.f90 b/src/self_colls_mod.f90
index 0cf9617..48bae81 100644
--- a/src/self_colls_mod.f90
+++ b/src/self_colls_mod.f90
@@ -1,447 +1,447 @@
MODULE self_colls
! self collisions module
USE PREC_CONST
USE BASIC
USE BASIC_MPI
USE MODEL
USE ARRAY, ONLY: Ceedp,Ciidp,rCee,rCii
USE OPERATOR_MODEL
USE OUTPUT
USE GRID
use restart
use auxval
use gkcoulomb
IMPLICIT NONE
PROCEDURE(Caa_sub), POINTER :: evaluate_Cii_collisional_matrix => null()
PROCEDURE(Caa_sub), POINTER :: evaluate_Cee_collisional_matrix => null()
ABSTRACT INTERFACE
SUBROUTINE Caa_sub
! Empty
END SUBROUTINE Caa_sub
END INTERFACE
CONTAINS
!_______________________________________________________________________
SUBROUTINE set_self_colls(ss)
! set self collisions parameters to species s
IMPLICIT NONE
!
!
CHARACTER(len=2),INTENT(IN):: ss
!
CALL set_model(ss) ! ... set temperature/mass ratio to 1 for both a = [e,i]
!
CALL set_basic(ss)
!
!
END SUBROUTINE set_self_colls
!
!_______________________________________________________________________
! electron-electron collision
SUBROUTINE ee_colls_control
!
IMPLICIT NONE
!
IF(me_self .eq. 0) WRITE(*,*) '====== Electron-Electron Collisions ======='
IF( ESELF .ne. 0 ) THEN
IF(me_self .eq. 0) WRITE(*,*) 'Set Electron-Electron Collisions ...'
CALL set_self_colls('ee') ! ... electrons-slef collision
CALL evaluate_auxval_field
! Restart
IF( ifrestart) THEN
IF(me .eq. 0 ) WRITE(*,*) 'Restart ...'
CALL restart_init('ee')
IF(me .eq. 0 ) WRITE(*,*) '... restart done'
ENDIF
IF(me_self .eq. 0) WRITE(*,*) 'Construct self collisional matrix Cee ... '
call set_Cee_collisional_matrix
call evaluate_Cee_collisional_matrix
WRITE(*,*) '... Cee constructed, me =', me_self
ENDIF
!
!
END SUBROUTINE ee_colls_control
!
!_______________________________________________________________________
!
SUBROUTINE set_Cee_collisional_matrix
! Set evaluate_Cee_collisional_matrix routine
IMPLICIT NONE
!
!
IF (ESELF .ne. 0 ) THEN
IF(IFCOULOMB .and. GKE .and. ESELF .eq. 1) THEN
evaluate_Cee_collisional_matrix => Cee_gk_coulomb
!
call set_gkcoulomb
IF(.not. GKC_USE_FULL ) then
call load_gk_coulomb_coeffs('ee')
else
call set_tabljpmf('s')
ENDIF
!
ELSE
! Evaluate Cee
evaluate_Cee_collisional_matrix => Cee_dflt
!
ENDIF
ELSE
evaluate_Cee_collisional_matrix => Cee_zeros_dflt
ENDIF
!
!
END SUBROUTINE set_Cee_collisional_matrix
!
!_______________________________________________________________________
!
SUBROUTINE Cee_dflt
! evaluate electron-electron collision matrix
use restart
IMPLICIT NONE
! Local variables
INTEGER :: irow_tmp,lbar_start
INTEGER :: irow
INTEGER, DIMENSION(2) :: lk
! current indices for restart
INTEGER :: istep
! Allocate restart coll mat and restart indices
IF(restart_from) CALL get_restart(Ceedp)
lbar_start = lbaree_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbaree_s_l
ENDIF
istep=1
- IF( me .eq. 0 ) write(*, '(a,1x, i3,a)' ) 'Cee [x/ ', lbaree_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i3,a)' ) 'Cee [x/ ', numbe(Pmaxe,Jmaxe) ,']:'
DO irow = lbar_start,lbaree_e_l
irow_tmp = rows_e(irow)
lk = invnumbe(irow_tmp)
write(nscr,'(i3)', advance = 'no') irow
flush(nscr)
rCee(:) = ZEROFM
CALL COLLISIONSELFE(lk(1),lk(2),rCee)
Ceedp(irow,:) = TO_DP(rCee(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ee) THEN
call restart_out(Ceedp, irow)
istep = 1
ENDIF
ENDDO
!
IF(ifrestart) call restart_out(Ceedp, irow)
!
END SUBROUTINE Cee_dflt
!_______________________________________________________________________
!
!
SUBROUTINE Cee_zeros_dflt
! evaluate electron-electron collision matrix
use restart
IMPLICIT NONE
! Local variables
INTEGER :: irow_tmp,lbar_start
INTEGER :: irow
INTEGER, DIMENSION(2) :: lk
! current indices for restart
INTEGER :: istep
lbar_start = lbaree_s_l
DO irow = lbar_start,lbaree_e_l
irow_tmp = rows_e(irow)
lk = invnumbe(irow_tmp)
rCee(:) = ZEROFM
Ceedp(irow,:) = TO_DP(rCee(:))
! save current state
istep = istep +1
ENDDO
!
!
END SUBROUTINE Cee_zeros_dflt
!
!-----------------------------------------------------------------
!
SUBROUTINE Cee_gk_coulomb
! Evaluate GK Coulomb
use restart
use gkcoulomb
implicit none
INTEGER :: l,k,irow
INTEGER :: j,n,np,istep
INTEGER :: lbar_start
INTEGER, DIMENSION(4) :: lbarjnnp
INTEGER, DIMENSION(2) :: lk
!
!
! restart
IF(restart_from) CALL get_restart(Ceedp)
!
lbar_start = lbaree_s_l
! initialize indices
lbar_start = lbaree_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbaree_s_l
ENDIF
istep=1
!
IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'Cee [x/ ', lbaree_e_l ,']:'
!
DO irow=lbar_start,lbaree_e_l
lbarjnnp = rows_e_4dim(irow,1:4)
! Hermite-Laguerre indices
lk = invnumbe(lbarjnnp(1))
j = lbarjnnp(2)
n = lbarjnnp(3)
np = lbarjnnp(4)
write(nscr,'(i5)', advance = 'no') lbarjnnp(1)
flush(nscr)
rCee(:) = ZEROFM
CALL coulombSelfGK(lk(1),lk(2),j,n,np,rCee)
Ceedp(irow,:) = TO_DP(rCee(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ee) THEN
! flush current state every nsave_ee.
call restart_out(Ceedp, irow)
istep = 1
ENDIF
ENDDO
! flush out the final state
IF(ifrestart) call restart_out(Ceedp, irow)
!
END SUBROUTINE Cee_gk_coulomb
!
!_______________________________________________________________________
! ion-ion collision
!
subroutine ii_colls_control
! ii_colls
implicit none
IF(me_self .eq. 0) WRITE(*,*) '====== Ion-Ion Collisions ======='
IF( ISELF .ne. 0 ) THEN
IF(me_self .eq. 0) WRITE(*,*) 'self Collisions ...'
CALL set_self_colls('ii') ! ... self collision
CALL evaluate_auxval_field
! Restart
IF( ifrestart) THEN
IF(me .eq. 0 ) WRITE(*,*) 'Restart ...'
CALL restart_init('ii')
IF(me .eq. 0 ) WRITE(*,*) '... restart done'
ENDIF
IF(me_self .eq. 0) WRITE(*,*) 'Construct ion-ion collisional matrix Cii ... '
call set_Cii_collisional_matrix
CALL evaluate_Cii_collisional_matrix
WRITE(*,*) '... Cii constructed, me =', me_self
ENDIF
!
end subroutine ii_colls_control
!
!__________________________________________________________________
!
SUBROUTINE set_Cii_collisional_matrix
! Set evaluate_Cii_collisional_matrix routine
IMPLICIT NONE
!
!
IF (ISELF .ne. 0 ) THEN
IF(IFCOULOMB .and. GKI .and. ISELF .eq. 1 ) THEN
evaluate_Cii_collisional_matrix => Cii_gk_coulomb
call set_gkcoulomb
IF( .not. GKC_USE_FULL) THEN
call load_gk_coulomb_coeffs('ii')
ELSE
call set_tabljpmf('s')
ENDIF
ELSE
!
evaluate_Cii_collisional_matrix => Cii_dflt
!
ENDIF
ELSE
evaluate_Cii_collisional_matrix => Cii_zeros_dflt
ENDIF
!
END SUBROUTINE set_Cii_collisional_matrix
!
!_______________________________________________________________________
!
SUBROUTINE Cii_dflt
! evaluate ion-ion collision matrix
use restart
IMPLICIT NONE
! Local variables
INTEGER :: irow_tmp,lbar_start
INTEGER :: irow
INTEGER, DIMENSION(2) :: lk
! current indices for restart
INTEGER :: istep
! Allocate restart coll mat and restart indices
IF(restart_from) CALL get_restart(Ciidp)
lbar_start = lbarii_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbarii_s_l
ENDIF
istep=1
- IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'Cii [x/ ', lbarii_e_l ,']:'
+ IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'Cii [x/ ', numbi(Pmaxi,Jmaxi) ,']:'
DO irow = lbar_start,lbarii_e_l
irow_tmp = rows_i(irow)
lk = invnumbi(irow_tmp)
write(nscr,'(i5)', advance = 'no') irow
flush(nscr)
rCii(:) = ZEROFM
CALL COLLISIONSELFI(lk(1),lk(2),rCii)
Ciidp(irow,:) = TO_DP(rCii(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ii) THEN
call restart_out(Ciidp, irow)
istep = 1
ENDIF
ENDDO
!
IF(ifrestart) call restart_out(Ciidp, irow)
!
END SUBROUTINE Cii_dflt
!_______________________________________________________________________
!
!
SUBROUTINE Cii_zeros_dflt
! evaluate ion-ion collision matrix
use restart
IMPLICIT NONE
! Local variables
INTEGER :: irow_tmp,lbar_start
INTEGER :: irow
INTEGER, DIMENSION(2) :: lk
! current indices for restart
INTEGER :: istep
lbar_start = lbarii_s_l
DO irow = lbar_start,lbarii_e_l
irow_tmp = rows_i(irow)
lk = invnumbi(irow_tmp)
write(*,*) 'Cii:', lk(1),'/',Pmaxi, ', ',lk(2),'/',Jmaxi
rCii(:) = ZEROFM
Ciidp(irow,:) = TO_DP(rCii(:))
! save current state
istep = istep +1
ENDDO
!
!
END SUBROUTINE Cii_zeros_dflt
!
!-----------------------------------------------------------------
!
SUBROUTINE Cii_gk_coulomb
! Evaluate GK Coulomb
use restart
use gkcoulomb
implicit none
INTEGER :: l,k,irow
INTEGER :: j,n,np,istep
INTEGER :: lbar_start
INTEGER, DIMENSION(4) :: lbarjnnp
INTEGER, DIMENSION(2) :: lk
!
!
! restart
IF(restart_from) CALL get_restart(Ciidp)
!
! initialize indices
lbar_start = lbarii_s_l
IF(restart_from) THEN
lbar_start = row_rst_l +1
IF(lbar_start .eq. 1) lbar_start = lbarii_s_l
ENDIF
IF( me .eq. 0 ) write(*, '(a,1x, i5,a)' ) 'Cii [x/ ', lbarii_e_l ,']:'
istep=1
!
DO irow=lbar_start,lbarii_e_l
lbarjnnp = rows_i_4dim(irow,1:4)
! Hermite-Laguerre indices
lk = invnumbi(lbarjnnp(1))
j = lbarjnnp(2)
n = lbarjnnp(3)
np = lbarjnnp(4)
write(nscr,'(i5)', advance = 'no') lbarjnnp(1)
flush(nscr)
rCii(:) = ZEROFM
CALL coulombSelfGK(lk(1),lk(2),j,n,np,rCii)
Ciidp(irow,:) = TO_DP(rCii(:))
! save current state
istep = istep +1
IF(ifrestart .and. istep .gt. nsave_ii) THEN
! flush current state every nsave_ii.
call restart_out(Ciidp,irow)
istep = 1
ENDIF
ENDDO
! flush out the final state
IF(ifrestart) call restart_out(Ciidp, irow)
!
END SUBROUTINE Cii_gk_Coulomb
!
!
END MODULE self_colls