From 9bea36bdb7f8687bda4185d9fac5cc6e8712d1c6 Mon Sep 17 00:00:00 2001 From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch> Date: Tue, 12 Jan 2021 11:50:39 +0100 Subject: [PATCH] Start writing collision operator in a separated module --- src/collision_mod.F90 | 501 +++++++++++++++++++++++++++--------------- 1 file changed, 318 insertions(+), 183 deletions(-) diff --git a/src/collision_mod.F90 b/src/collision_mod.F90 index 7cbf0aff..35c75e42 100644 --- a/src/collision_mod.F90 +++ b/src/collision_mod.F90 @@ -1,194 +1,329 @@ module collision ! contains the Hermite-Laguerre collision operators. Solved using COSOlver. -use prec_const -use fields -use array -use grid -use basic -use futils -use initial_par -use model - implicit none PUBLIC :: load_FC_mat +PUBLIC :: DoughertyGK_e, DoughertyGK_i CONTAINS -!******************************************************************************! -SUBROUTINE LenardBernsteinDK - -END SUBROUTINE LenardBernsteinDK - -!******************************************************************************! -SUBROUTINE DoughertyDK - -END SUBROUTINE DoughertyDK - -!******************************************************************************! -SUBROUTINE FullCoulombDK - -END SUBROUTINE FullCoulombDK - - -!******************************************************************************! -!!!!!!! Load the Full coulomb coefficient table from COSOlver results -!******************************************************************************! -SUBROUTINE load_FC_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full) - IMPLICIT NONE - - INTEGER :: irow_sub, irow_full, icol_sub, icol_full - INTEGER :: fid1, fid2, fid3, fid4 - - INTEGER :: ip_e, ij_e, il_e, ik_e - INTEGER :: pdime, jdime - REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ceepj_full, CeipjT_full - REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF_full - INTEGER :: ip_i, ij_i, il_i, ik_i - INTEGER :: pdimi, jdimi - REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj_full, CiepjT_full - REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF_full - - !!!!!!!!!!!! Electron matrices !!!!!!!!!!!! - IF (my_id .EQ. 0) WRITE(*,*) 'Load ee FC mat...' - ! get the self electron colision matrix - CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD) - - CALL getatt(fid1,'/Caapj/Ceepj/','Pmaxe',pdime) - CALL getatt(fid1,'/Caapj/Ceepj/','Jmaxe',jdime) - - IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! FC Matrix too small !!' - - CALL allocate_array( Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1)) - CALL getarr(fid1, '/Caapj/Ceepj', Ceepj_full) ! get array (moli format) - - ! Fill sub array with only usefull polynmials degree - DO ip_e = 0,pmaxe ! Loop over rows - DO ij_e = 0,jmaxe - irow_sub = (jmaxe +1)*ip_e + ij_e +1 - irow_full = (jdime +1)*ip_e + ij_e +1 - DO il_e = 0,pmaxe ! Loop over columns - DO ik_e = 0,jmaxe - icol_sub = (jmaxe +1)*il_e + ik_e +1 - icol_full = (jdime +1)*il_e + ik_e +1 - Ceepj (irow_sub,icol_sub) = Ceepj_full (irow_full,icol_full) - ENDDO - ENDDO - ENDDO - ENDDO - - CALL closef(fid1) - DEALLOCATE(Ceepj_full) - - IF (my_id .EQ. 0) WRITE(*,*) 'Load ei FC mat...' - ! get the Test and Back field electron ion collision matrix - CALL openf(eimat_file,fid2, 'r', 'D'); - - CALL getatt(fid2,'/Ceipj/CeipjT/','Pmaxi',pdimi) - CALL getatt(fid2,'/Ceipj/CeipjT/','Jmaxi',jdimi) - IF ( (pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi) ) WRITE(*,*) '!! FC Matrix too small !!' - - CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1)) - CALL allocate_array( CeipjF_full, 1,(pdime+1)*(jdime+1), 1,(pdimi+1)*(jdimi+1)) - - CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT_full) - CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF_full) - - ! Fill sub array with only usefull polynmials degree - DO ip_e = 0,pmaxe ! Loop over rows - DO ij_e = 0,jmaxe - irow_sub = (jmaxe +1)*ip_e + ij_e +1 - irow_full = (jdime +1)*ip_e + ij_e +1 - DO il_e = 0,pmaxe ! Loop over columns - DO ik_e = 0,jmaxe - icol_sub = (jmaxe +1)*il_e + ik_e +1 - icol_full = (jdime +1)*il_e + ik_e +1 - CeipjT(irow_sub,icol_sub) = CeipjT_full(irow_full,icol_full) - ENDDO - ENDDO - DO il_i = 0,pmaxi ! Loop over columns - DO ik_i = 0,jmaxi - icol_sub = (Jmaxi +1)*il_i + ik_i +1 - icol_full = (jdimi +1)*il_i + ik_i +1 - CeipjF(irow_sub,icol_sub) = CeipjF_full(irow_full,icol_full) - ENDDO - ENDDO - ENDDO - ENDDO - - CALL closef(fid2) - DEALLOCATE(CeipjF_full) - DEALLOCATE(CeipjT_full) - - !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!! - IF (my_id .EQ. 0) WRITE(*,*) 'Load ii FC mat...' - ! get the self electron colision matrix - CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD); - - CALL allocate_array( Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1)) - - IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj - CALL getarr(fid3, '/Caapj/Ceepj', Ciipj_full) ! get array (moli format) - ELSE - CALL getarr(fid3, '/Caapj/Ciipj', Ciipj_full) ! get array (moli format) - ENDIF - - ! Fill sub array with only usefull polynmials degree - DO ip_i = 0,Pmaxi ! Loop over rows - DO ij_i = 0,Jmaxi - irow_sub = (Jmaxi +1)*ip_i + ij_i +1 - irow_full = (jdimi +1)*ip_i + ij_i +1 - DO il_i = 0,Pmaxi ! Loop over columns - DO ik_i = 0,Jmaxi - icol_sub = (Jmaxi +1)*il_i + ik_i +1 - icol_full = (jdimi +1)*il_i + ik_i +1 - Ciipj (irow_sub,icol_sub) = Ciipj_full (irow_full,icol_full) - ENDDO - ENDDO - ENDDO - ENDDO - - CALL closef(fid3) - DEALLOCATE(Ciipj_full) - - ! get the Test and Back field electron ion collision matrix - CALL openf(iemat_file,fid4, 'r', 'D'); - - CALL allocate_array( CiepjT_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1)) - CALL allocate_array( CiepjF_full, 1,(pdimi+1)*(jdimi+1), 1,(pdime+1)*(jdime+1)) - - IF (my_id .EQ. 0) WRITE(*,*) 'Load ie FC mat...' - CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT_full) - CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF_full) - - ! Fill sub array with only usefull polynmials degree - DO ip_i = 0,Pmaxi ! Loop over rows - DO ij_i = 0,Jmaxi - irow_sub = (Jmaxi +1)*ip_i + ij_i +1 - irow_full = (jdimi +1)*ip_i + ij_i +1 - DO il_i = 0,Pmaxi ! Loop over columns - DO ik_i = 0,Jmaxi - icol_sub = (Jmaxi +1)*il_i + ik_i +1 - icol_full = (jdimi +1)*il_i + ik_i +1 - CiepjT(irow_sub,icol_sub) = CiepjT_full(irow_full,icol_full) - ENDDO - ENDDO - DO il_e = 0,pmaxe ! Loop over columns - DO ik_e = 0,jmaxe - icol_sub = (jmaxe +1)*il_e + ik_e +1 - icol_full = (jdime +1)*il_e + ik_e +1 - CiepjF(irow_sub,icol_sub) = CiepjF_full(irow_full,icol_full) - ENDDO - ENDDO - ENDDO - ENDDO - - CALL closef(fid4) - DEALLOCATE(CiepjF_full) - DEALLOCATE(CiepjT_full) - -END SUBROUTINE load_FC_mat -!******************************************************************************! + !******************************************************************************! + SUBROUTINE LenardBernsteinDK + + END SUBROUTINE LenardBernsteinDK + + !******************************************************************************! + !! Doughtery gyrokinetic collision operator for electrons + SUBROUTINE DoughertyGK_e(ip,ij,ikr,ikz,TColl) + USE fields + USE array, ONLY : Kernel_e + USE basic + USE grid, ONLY : jmaxe, parray_e, jarray_e, krarray, kzarray + USE prec_const + USE time_integration + USE model + IMPLICIT NONE + INTEGER, INTENT(IN) :: ip,ij,ikr,ikz + COMPLEX(dp), INTENT(INOUT) :: TColl + + COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_ + COMPLEX(dp) :: Dpj, Ppj, T_, ibe_ + INTEGER :: in_ + REAL :: n_dp, j_dp, p_dp, be_2 + + !** Auxiliary variables ** + p_dp = REAL(parray_e(ij),dp) + j_dp = REAL(jarray_e(ij),dp) + be_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmae2_taue_o2 + ibe_ = imagu*SQRT(be_2) + + !** build fluid moments used for collison operator ** + n_ = 0._dp + upar_ = 0._dp + uperp_ = 0._dp + Tpar_ = 0._dp + Tperp_ = 0._dp + DO in_ = 1,jmaxe+1 + n_dp = REAL(in_-1,dp) + n_ = n_ + Kernel_e(in_,ikr,ikz) * moments_e(1,in_,ikr,ikz,updatetlevel) + upar_ = upar_ + Kernel_e(in_,ikr,ikz) * moments_e(2,in_,ikr,ikz,updatetlevel) + uperp_ = uperp_ + ibe_*0.5_dp*Kernel_e(in_,ikr,ikz) * (moments_e(1,in_,ikr,ikz,updatetlevel) -moments_e(1,in_+1,ikr,ikz,updatetlevel)) + Tpar_ = Tpar_ + Kernel_e(in_,ikr,ikz) * (SQRT2*moments_e(3,in_,ikr,ikz,updatetlevel) + moments_e(1,in_,ikr,ikz,updatetlevel)) + Tperp_ = Tperp_ + Kernel_e(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*moments_e(1,in_,ikr,ikz,updatetlevel) - n_dp*moments_e(1,in_-1,ikr,ikz,updatetlevel) - (n_dp+1)*moments_e(1,in_+1,ikr,ikz,updatetlevel)) + ENDDO + T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_ + + ! Velocity-space diffusion + TColl = -(2._dp*p_dp + j_dp + be_2)*moments_e(ip,ij,ikr,ikz,updatetlevel) + + ! Add energy restoring term + IF( ip .eq. 0 ) THEN + TColl = TColl & + + T_*( 2._dp*((2_dp*j_dp + 1._dp)*Kernel_e(ij,ikr,ikz) & + -(j_dp + 1._dp)*Kernel_e(ij+1,ikr,ikz) & + - j_dp*Kernel_e(ij-1,ikr,ikz)) & + - 2._dp*Kernel_e(ij,ikr,ikz)) + ENDIF + + IF( ip .eq. 2 ) THEN + TColl = TColl + T_*SQRT2*Kernel_e(ij,ikr,ikz) + ENDIF + + !Add momentum restoring term + IF( ip .eq. 1) THEN + TColl = TColl + upar_*Kernel_e(ij,ikr,ikz) + ENDIF + + IF(ip .eq. 0 ) THEN + TColl = TColl + uperp_*ibe_*( (j_dp + 1._dp)*Kernel_e(ij,ikr,ikz) - j_dp*Kernel_e(ij-1,ikr,ikz)) + ENDIF + + TColl = TColl*nu_e + + END SUBROUTINE DoughertyGK_e + + !******************************************************************************! + !! Doughtery gyrokinetic collision operator for ions + SUBROUTINE DoughertyGK_i(ip,ij,ikr,ikz,TColl) + USE fields + USE array, ONLY : Kernel_i + USE basic + USE grid, ONLY : jmaxi, parray_i, jarray_i, krarray, kzarray + USE prec_const + USE time_integration + USE model + IMPLICIT NONE + INTEGER, INTENT(IN) :: ip,ij,ikr,ikz + COMPLEX(dp), INTENT(INOUT) :: TColl + + COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_ + COMPLEX(dp) :: Dpj, Ppj, T_, ibi_ + INTEGER :: in_ + REAL :: n_dp, j_dp, p_dp, bi_2 + + !** Auxiliary variables ** + p_dp = REAL(parray_i(ij),dp) + j_dp = REAL(jarray_i(ij),dp) + bi_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmai2_taui_o2 + ibi_ = imagu*SQRT(bi_2) + + !** build fluid moments used for collison operator ** + n_ = 0._dp + upar_ = 0._dp + uperp_ = 0._dp + Tpar_ = 0._dp + Tperp_ = 0._dp + DO in_ = 1,jmaxi+1 + n_dp = REAL(in_-1,dp) + n_ = n_ + Kernel_i(in_,ikr,ikz) * moments_i(1,in_,ikr,ikz,updatetlevel) + upar_ = upar_ + Kernel_i(in_,ikr,ikz) * moments_i(2,in_,ikr,ikz,updatetlevel) + uperp_ = uperp_ + 0.5_dp*Kernel_i(in_,ikr,ikz) * (moments_i(1,in_,ikr,ikz,updatetlevel) -moments_i(1,in_+1,ikr,ikz,updatetlevel)) + Tpar_ = Tpar_ + Kernel_i(in_,ikr,ikz) * (SQRT2*moments_i(3,in_,ikr,ikz,updatetlevel) + moments_i(1,in_,ikr,ikz,updatetlevel)) + Tperp_ = Tperp_ + Kernel_i(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*moments_i(1,in_,ikr,ikz,updatetlevel) - n_dp*moments_i(1,in_-1,ikr,ikz,updatetlevel) - (n_dp+1)*moments_i(1,in_+1,ikr,ikz,updatetlevel)) + ENDDO + T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_ + + ! Velocity-space diffusion + TColl = -(2._dp*p_dp + j_dp + bi_2)*moments_i(ip,ij,ikr,ikz,updatetlevel) + + ! Add energy restoring term + IF( ip .eq. 0 ) THEN + TColl = TColl & + + T_*( 2._dp*((2_dp*j_dp + 1._dp)*Kernel_i(ij,ikr,ikz) & + -(j_dp + 1._dp)*Kernel_i(ij+1,ikr,ikz) & + - j_dp*Kernel_i(ij-1,ikr,ikz)) & + - 2._dp*Kernel_i(ij,ikr,ikz)) + ENDIF + + IF( ip .eq. 2 ) THEN + TColl = TColl + T_*SQRT2*Kernel_i(ij,ikr,ikz) + ENDIF + + !Add momentum restoring term + IF( ip .eq. 1) THEN + TColl = TColl + upar_*Kernel_i(ij,ikr,ikz) + ENDIF + + IF(ip .eq. 0 ) THEN + TColl = TColl + uperp_*imagu*SQRT(bi_2)*( (j_dp + 1._dp)*Kernel_i(ij,ikr,ikz) - j_dp*Kernel_i(ij-1,ikr,ikz)) + ENDIF + + TColl = TColl*nu_i + + END SUBROUTINE DoughertyGK_i + + !******************************************************************************! + SUBROUTINE FullCoulombDK + + END SUBROUTINE FullCoulombDK + + + !******************************************************************************! + !!!!!!! Load the Full coulomb coefficient table from COSOlver results + !******************************************************************************! + SUBROUTINE load_FC_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full) + use prec_const + use fields + use array + use grid + use basic + use futils + use initial_par + use model + IMPLICIT NONE + + INTEGER :: irow_sub, irow_full, icol_sub, icol_full + INTEGER :: fid1, fid2, fid3, fid4 + + INTEGER :: ip_e, ij_e, il_e, ik_e + INTEGER :: pdime, jdime + REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ceepj_full, CeipjT_full + REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF_full + INTEGER :: ip_i, ij_i, il_i, ik_i + INTEGER :: pdimi, jdimi + REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj_full, CiepjT_full + REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF_full + + !!!!!!!!!!!! Electron matrices !!!!!!!!!!!! + IF (my_id .EQ. 0) WRITE(*,*) 'Load ee FC mat...' + ! get the self electron colision matrix + CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD) + + CALL getatt(fid1,'/Caapj/Ceepj/','Pmaxe',pdime) + CALL getatt(fid1,'/Caapj/Ceepj/','Jmaxe',jdime) + + IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! FC Matrix too small !!' + + CALL allocate_array( Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1)) + CALL getarr(fid1, '/Caapj/Ceepj', Ceepj_full) ! get array (moli format) + + ! Fill sub array with only usefull polynmials degree + DO ip_e = 0,pmaxe ! Loop over rows + DO ij_e = 0,jmaxe + irow_sub = (jmaxe +1)*ip_e + ij_e +1 + irow_full = (jdime +1)*ip_e + ij_e +1 + DO il_e = 0,pmaxe ! Loop over columns + DO ik_e = 0,jmaxe + icol_sub = (jmaxe +1)*il_e + ik_e +1 + icol_full = (jdime +1)*il_e + ik_e +1 + Ceepj (irow_sub,icol_sub) = Ceepj_full (irow_full,icol_full) + ENDDO + ENDDO + ENDDO + ENDDO + + CALL closef(fid1) + DEALLOCATE(Ceepj_full) + + IF (my_id .EQ. 0) WRITE(*,*) 'Load ei FC mat...' + ! get the Test and Back field electron ion collision matrix + CALL openf(eimat_file,fid2, 'r', 'D'); + + CALL getatt(fid2,'/Ceipj/CeipjT/','Pmaxi',pdimi) + CALL getatt(fid2,'/Ceipj/CeipjT/','Jmaxi',jdimi) + IF ( (pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi) ) WRITE(*,*) '!! FC Matrix too small !!' + + CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1)) + CALL allocate_array( CeipjF_full, 1,(pdime+1)*(jdime+1), 1,(pdimi+1)*(jdimi+1)) + + CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT_full) + CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF_full) + + ! Fill sub array with only usefull polynmials degree + DO ip_e = 0,pmaxe ! Loop over rows + DO ij_e = 0,jmaxe + irow_sub = (jmaxe +1)*ip_e + ij_e +1 + irow_full = (jdime +1)*ip_e + ij_e +1 + DO il_e = 0,pmaxe ! Loop over columns + DO ik_e = 0,jmaxe + icol_sub = (jmaxe +1)*il_e + ik_e +1 + icol_full = (jdime +1)*il_e + ik_e +1 + CeipjT(irow_sub,icol_sub) = CeipjT_full(irow_full,icol_full) + ENDDO + ENDDO + DO il_i = 0,pmaxi ! Loop over columns + DO ik_i = 0,jmaxi + icol_sub = (Jmaxi +1)*il_i + ik_i +1 + icol_full = (jdimi +1)*il_i + ik_i +1 + CeipjF(irow_sub,icol_sub) = CeipjF_full(irow_full,icol_full) + ENDDO + ENDDO + ENDDO + ENDDO + + CALL closef(fid2) + DEALLOCATE(CeipjF_full) + DEALLOCATE(CeipjT_full) + + !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!! + IF (my_id .EQ. 0) WRITE(*,*) 'Load ii FC mat...' + ! get the self electron colision matrix + CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD); + + CALL allocate_array( Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1)) + + IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj + CALL getarr(fid3, '/Caapj/Ceepj', Ciipj_full) ! get array (moli format) + ELSE + CALL getarr(fid3, '/Caapj/Ciipj', Ciipj_full) ! get array (moli format) + ENDIF + + ! Fill sub array with only usefull polynmials degree + DO ip_i = 0,Pmaxi ! Loop over rows + DO ij_i = 0,Jmaxi + irow_sub = (Jmaxi +1)*ip_i + ij_i +1 + irow_full = (jdimi +1)*ip_i + ij_i +1 + DO il_i = 0,Pmaxi ! Loop over columns + DO ik_i = 0,Jmaxi + icol_sub = (Jmaxi +1)*il_i + ik_i +1 + icol_full = (jdimi +1)*il_i + ik_i +1 + Ciipj (irow_sub,icol_sub) = Ciipj_full (irow_full,icol_full) + ENDDO + ENDDO + ENDDO + ENDDO + + CALL closef(fid3) + DEALLOCATE(Ciipj_full) + + ! get the Test and Back field electron ion collision matrix + CALL openf(iemat_file,fid4, 'r', 'D'); + + CALL allocate_array( CiepjT_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1)) + CALL allocate_array( CiepjF_full, 1,(pdimi+1)*(jdimi+1), 1,(pdime+1)*(jdime+1)) + + IF (my_id .EQ. 0) WRITE(*,*) 'Load ie FC mat...' + CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT_full) + CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF_full) + + ! Fill sub array with only usefull polynmials degree + DO ip_i = 0,Pmaxi ! Loop over rows + DO ij_i = 0,Jmaxi + irow_sub = (Jmaxi +1)*ip_i + ij_i +1 + irow_full = (jdimi +1)*ip_i + ij_i +1 + DO il_i = 0,Pmaxi ! Loop over columns + DO ik_i = 0,Jmaxi + icol_sub = (Jmaxi +1)*il_i + ik_i +1 + icol_full = (jdimi +1)*il_i + ik_i +1 + CiepjT(irow_sub,icol_sub) = CiepjT_full(irow_full,icol_full) + ENDDO + ENDDO + DO il_e = 0,pmaxe ! Loop over columns + DO ik_e = 0,jmaxe + icol_sub = (jmaxe +1)*il_e + ik_e +1 + icol_full = (jdime +1)*il_e + ik_e +1 + CiepjF(irow_sub,icol_sub) = CiepjF_full(irow_full,icol_full) + ENDDO + ENDDO + ENDDO + ENDDO + + CALL closef(fid4) + DEALLOCATE(CiepjF_full) + DEALLOCATE(CiepjT_full) + + END SUBROUTINE load_FC_mat + !******************************************************************************! end module collision -- GitLab