diff --git a/python_scripts/gkcoulomb.py b/python_scripts/gkcoulomb.py
index 2789096..a735894 100644
--- a/python_scripts/gkcoulomb.py
+++ b/python_scripts/gkcoulomb.py
@@ -1,229 +1,243 @@
 #!/usr/bin/env python
+# This scripts compute the GK Coulomb collision operator from the COSOlver output files
+# Performs the summation over the last indices, i.e. (p,j,m,n,np) -> (p,j)
+# Writes the results matrices 
+
 import numpy as np
 import os
 import sys
 import h5py as h5
 
+
 # Define directories
 cwd = os.getcwd()
+cwd = '/misc/bjfrei/cosolver_marconi/gk.coulomb.TEM.kperp.0.7080'
+
+print('Create GK Coulomb from ' + cwd)
 
 # Define h5 files 
 eicolls_file = os.path.join(cwd, 'ei.h5')
 iecolls_file = os.path.join(cwd,'ie.h5')
 selfcolls_file =  os.path.join( cwd,  'self.h5')
 
 # Defines new h5 files
-eicolls_file_out = os.path.join(cwd, 'ei.0.h5')
-iecolls_file_out = os.path.join(cwd,'ie.0.h5')
-selfcolls_file_out =  os.path.join( cwd,  'self.0.h5')
+eicolls_file_out = os.path.join(cwd, 'ei.pj.h5')
+iecolls_file_out = os.path.join(cwd,'ie.pj.h5')
+selfcolls_file_out =  os.path.join( cwd,  'self.pj.h5')
 
 
 ############################################################
 # Self collisions
+print('Self collisions....')
 
 # Open file
 selfh5 = h5.File(selfcolls_file, 'r')
+#selfeeh5=  h5.File('/home/bjfrei/Documents/cosolver/dk.Coulomb.Cee.P16J8/self.h5', 'r')
 selfh5_out = h5.File(selfcolls_file_out, 'w')
 
 # Ion-ion
 
 # get dimension
 Pmaxi = selfh5['/Caapj/Ciipjpjnn'].attrs['Pmaxi']
 Jmaxi = selfh5['/Caapj/Ciipjpjnn'].attrs['Jmaxi']
 
 def numb(pp,jj):
     lbari = (Jmaxi +1)*pp + jj +1
     return lbari
 
 # Get 4 dim rows
 rows_4dim= selfh5['/Caapj/rows_i']
 # Get Cippj
 Cpj_4dim=selfh5['/Caapj/Ciipjpjnn']
 #Cpj_4dim
 
 Cpj = np.zeros((numb(Pmaxi,Jmaxi), numb(Pmaxi,Jmaxi)))
 
 for irow in np.arange(1, numb(Pmaxi,Jmaxi)+1,1):
     idx_of_row, = np.where(rows_4dim[0,:] == irow)
     Cpj_ = Cpj_4dim[:, idx_of_row]
     Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
 
 # Create output    
 Cpj_dst = selfh5_out.create_dataset("/Caapj/Ciipj", data = Cpj)
 Cpj_dst.attrs['Pmaxi'] =  Pmaxi
 Cpj_dst.attrs['Jmaxi'] = Jmaxi
 Cpj_dst.attrs['kperp'] = selfh5['/Caapj/Ceepjpjnn'].attrs['kperp']
 
-# Electron-Electron
+#Electron-Electron
 
-# get dimension
+#get dimension
 Pmaxe = selfh5['/Caapj/Ceepjpjnn'].attrs['Pmaxe']
 Jmaxe = selfh5['/Caapj/Ceepjpjnn'].attrs['Jmaxe']
 
 def numb(pp,jj):
     lbari = (Jmaxe +1)*pp + jj +1
     return lbari
 
 # Get 4 dim rows
 rows_4dim= selfh5['/Caapj/rows_e']
 # Get Cippj
 Cpj_4dim=selfh5['/Caapj/Ceepjpjnn']
 #Cpj_4dim
 
 Cpj = np.zeros((numb(Pmaxe,Jmaxe), numb(Pmaxe,Jmaxe)))
 
 for irow in np.arange(1, numb(Pmaxe,Jmaxe)+1,1):
     idx_of_row, = np.where(rows_4dim[0,:] == irow)
     Cpj_ = Cpj_4dim[:, idx_of_row]
     Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
 
-# Create output    
+#Create output    
 Cpj_dst = selfh5_out.create_dataset("/Caapj/Ceepj", data = Cpj)
+#Cpj_dst = selfh5_out.create_dataset("/Caapj/Ceepj", data = selfeeh5['/Caapj/Ceepj'])
+
+
 Cpj_dst.attrs['Pmaxe'] =  Pmaxe
 Cpj_dst.attrs['Jmaxe'] = Jmaxe
 Cpj_dst.attrs['kperp'] = selfh5['/Caapj/Ceepjpjnn'].attrs['kperp']
 
-selfh5_out.create_dataset("/files/fort.90", data = selfh5['/files/fort.f90'])
+selfh5_out.create_dataset("/files/fort.90", data = selfh5['/files/fort.90'])
 
 # close file
 selfh5_out.close()
 selfh5.close()
 
 ############################################################
 # Ion-elctron collision
+print('Self ion-electron collisions....')
 
 ieh5 = h5.File(iecolls_file, 'r')
 ieh5_out = h5.File(iecolls_file_out, 'w')
 
 ### Test component #####
 
 # get dimension
 Pmaxi = ieh5['/Ciepj/CiepjTpjnn'].attrs['Pmaxi']
 Jmaxi = ieh5['/Ciepj/CiepjTpjnn'].attrs['Jmaxi']
 
 def numb(pp,jj):
     lbari = (Jmaxi +1)*pp + jj +1
     return lbari
 
 # Get 4 dim rows
 rows_4dim= ieh5['/Ciepj/rows_iT']
 # Get Cippj
 Cpj_4dim= ieh5['/Ciepj/CiepjTpjnn']
 
 Cpj = np.zeros((numb(Pmaxi,Jmaxi), numb(Pmaxi,Jmaxi)))
 
 for irow in np.arange(1, numb(Pmaxi,Jmaxi)+1,1):
     idx_of_row, = np.where(rows_4dim[0,:] == irow)
     Cpj_ = Cpj_4dim[:, idx_of_row]
     Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
 
 Cpj_dst = ieh5_out.create_dataset("/Ciepj/CiepjT", data = Cpj)
 Cpj_dst.attrs['Pmaxi'] =  Pmaxi
 Cpj_dst.attrs['Jmaxi'] = Jmaxi
 Cpj_dst.attrs['kperp'] = ieh5['/Ciepj/CiepjTpjnn'].attrs['kperp']
 
 
 ### Field component #####
 
 # get dimension
 Pmaxi = ieh5['/Ciepj/CiepjFpjnn'].attrs['Pmaxi']
 Jmaxi = ieh5['/Ciepj/CiepjFpjnn'].attrs['Jmaxi']
 
 def numb(pp,jj):
     lbari = (Jmaxi +1)*pp + jj +1
     return lbari
 
 # Get 4 dim rows
 rows_4dim= ieh5['/Ciepj/rows_iF']
 # Get Cippj
 Cpj_4dim= ieh5['/Ciepj/CiepjFpjnn']
 
 Cpj = np.zeros((numb(Pmaxi,Jmaxi), numb(Pmaxi,Jmaxi)))
 
 for irow in np.arange(1, numb(Pmaxi,Jmaxi)+1,1):
     idx_of_row, = np.where(rows_4dim[0,:] == irow)
     Cpj_ = Cpj_4dim[:, idx_of_row]
     Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
 
 Cpj_dst = ieh5_out.create_dataset("/Ciepj/CiepjF", data = Cpj)
 Cpj_dst.attrs['Pmaxi'] =  Pmaxi
 Cpj_dst.attrs['Jmaxi'] = Jmaxi
 Cpj_dst.attrs['kperp'] = ieh5['/Ciepj/CiepjFpjnn'].attrs['kperp']
 
-ieh5_out.create_dataset("/files/fort.90", data = ieh5['/files/fort.f90'])
+ieh5_out.create_dataset("/files/fort.90", data = ieh5['/files/fort.90'])
 
 # close file
 ieh5.close()
 ieh5_out.close()
 
 ############################################################
 # Electron-ion collision
-
+print('Self electron-ion collisions....')
 ### Test component #####
 
 # Open file
 eih5 = h5.File(eicolls_file, 'r')
 eih5_out = h5.File(eicolls_file_out, 'w')
 
 # get dimension
 Pmaxe = eih5['/Ceipj/CeipjTpjnn'].attrs['Pmaxi']
 Jmaxe = eih5['/Ceipj/CeipjTpjnn'].attrs['Jmaxi']
 
 def numb(pp,jj):
     lbar = (Jmaxe +1)*pp + jj +1
     return lbar
 
 # Get 4 dim rows
 rows_4dim= eih5['/Ceipj/rows_eT']
 # Get Cippj
 Cpj_4dim= eih5['/Ceipj/CeipjTpjnn']
 
 Cpj = np.zeros((numb(Pmaxe,Jmaxe), numb(Pmaxe,Jmaxe)))
 
 for irow in np.arange(1, numb(Pmaxe,Jmaxe)+1,1):
     idx_of_row, = np.where(rows_4dim[0,:] == irow)
     Cpj_ = Cpj_4dim[:, idx_of_row]
     Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
 
 Cpj_dst = eih5_out.create_dataset("/Ceipj/CeipjT", data = Cpj)
 Cpj_dst.attrs['Pmaxe'] =  Pmaxe
 Cpj_dst.attrs['Jmaxe'] = Jmaxe
 Cpj_dst.attrs['kperp'] = eih5['/Ceipj/CeipjTpjnn'].attrs['kperp']
 
 
 ### Field component #####
 
 # get dimension
 Pmaxe = eih5['/Ceipj/CeipjFpjnn'].attrs['Pmaxe']
 Jmaxe = eih5['/Ceipj/CeipjFpjnn'].attrs['Jmaxe']
 
 def numb(pp,jj):
     lbar = (Jmaxe +1)*pp + jj +1
     return lbar
 
 # Get 4 dim rows
 rows_4dim= eih5['/Ceipj/rows_eF']
 
 # Get Cippj
 Cpj_4dim= eih5['/Ceipj/CeipjFpjnn']
 
 Cpj = np.zeros((numb(Pmaxe, Jmaxe), numb(Pmaxe,Jmaxe)))
 
 for irow in np.arange(1, numb(Pmaxe,Jmaxe)+1,1):
     idx_of_row, = np.where(rows_4dim[0,:] == irow)
     Cpj_ = Cpj_4dim[:, idx_of_row]
     Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
 
 Cpj_dst = eih5_out.create_dataset("/Ceipj/CeipjF", data = Cpj)
 Cpj_dst.attrs['Pmaxe'] =  Pmaxe
 Cpj_dst.attrs['Jmaxe'] = Jmaxe
 Cpj_dst.attrs['kperp'] = eih5['/Ceipj/CeipjFpjnn'].attrs['kperp']
 
-eih5_out.create_dataset("/files/fort.90", data = eih5['/files/fort.f90'])
+eih5_out.create_dataset("/files/fort.90", data = eih5['/files/fort.90'])
 
 # close file
 eih5.close()
 eih5_out.close()
 
 
 
diff --git a/python_scripts/gkcoulomb_ee.py b/python_scripts/gkcoulomb_ee.py
new file mode 100644
index 0000000..ed1e6ce
--- /dev/null
+++ b/python_scripts/gkcoulomb_ee.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+# This scripts compute the GK Coulomb collision operator from the COSOlver output files
+# Performs the summation over the last indices, i.e. (p,j,m,n,np) -> (p,j)
+# Writes the results matrices 
+
+import numpy as np
+import os
+import sys
+import h5py as h5
+
+
+# Define directories
+cwd = os.getcwd()
+cwd = '/misc/bjfrei/cosolver_marconi/gk.coulomb.TEM.kperp.1.18.m.4.NFLR.8.ee'
+
+print('Create GK Coulomb from ' + cwd)
+
+# Define h5 files 
+selfcolls_file =  os.path.join( cwd,  'self.h5')
+
+# Defines new h5 files
+selfcolls_file_out =  os.path.join( cwd,  'self.pj.h5')
+
+
+############################################################
+# Self collisions
+print('Self collisions....')
+
+# Open file
+selfh5 = h5.File(selfcolls_file, 'r')
+selfeeh5=  h5.File('/home/bjfrei/Documents/cosolver/dk.Coulomb.Cee.P16J8/self.h5', 'r')
+selfh5_out = h5.File(selfcolls_file_out, 'w')
+
+
+#Electron-Electron
+
+#get dimension
+Pmaxe = selfh5['/Caapj/Ceepjpjnn'].attrs['Pmaxe']
+Jmaxe = selfh5['/Caapj/Ceepjpjnn'].attrs['Jmaxe']
+
+def numb(pp,jj):
+    lbari = (Jmaxe +1)*pp + jj +1
+    return lbari
+
+# Get 4 dim rows
+rows_4dim= selfh5['/Caapj/rows_e']
+# Get Cippj
+Cpj_4dim=selfh5['/Caapj/Ceepjpjnn']
+#Cpj_4dim
+
+Cpj = np.zeros((numb(Pmaxe,Jmaxe), numb(Pmaxe,Jmaxe)))
+
+for irow in np.arange(1, numb(Pmaxe,Jmaxe)+1,1):
+    idx_of_row, = np.where(rows_4dim[0,:] == irow)
+    Cpj_ = Cpj_4dim[:, idx_of_row]
+    Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
+
+#Create output    
+Cpj_dst = selfh5_out.create_dataset("/Caapj/Ceepj", data = Cpj)
+#Cpj_dst = selfh5_out.create_dataset("/Caapj/Ceepj", data = selfeeh5['/Caapj/Ceepj'])
+
+
+Cpj_dst.attrs['Pmaxe'] =  Pmaxe
+Cpj_dst.attrs['Jmaxe'] = Jmaxe
+Cpj_dst.attrs['kperp'] = selfh5['/Caapj/Ceepjpjnn'].attrs['kperp']
+
+selfh5_out.create_dataset("/files/fort.90", data = selfh5['/files/fort.90'])
+
+# close file
+selfh5_out.close()
+selfh5.close()
diff --git a/python_scripts/gkcoulomb_ei.py b/python_scripts/gkcoulomb_ei.py
new file mode 100644
index 0000000..d12e4d6
--- /dev/null
+++ b/python_scripts/gkcoulomb_ei.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+# This scripts compute the GK Coulomb collision operator from the COSOlver output files
+# Performs the summation over the last indices, i.e. (p,j,m,n,np) -> (p,j)
+# Writes the results matrices 
+
+import numpy as np
+import os
+import sys
+import h5py as h5
+
+
+# Define directories
+cwd = os.getcwd()
+
+print('Create GK Coulomb from ' + cwd)
+
+# Define h5 files input
+eicolls_file = os.path.join(cwd, 'ei.h5')
+
+# Defines new h5 files
+eicolls_file_out = os.path.join(cwd, 'ei.pj.h5')
+
+
+############################################################
+# Electron-ion collision
+print('Self electron-ion collisions....')
+
+### Test component #####
+
+# Open file
+eih5 = h5.File(eicolls_file, 'r')
+eih5_out = h5.File(eicolls_file_out, 'w')
+
+# get dimension
+Pmaxe = eih5['/Ceipj/CeipjTpjnn'].attrs['Pmaxi']
+Jmaxe = eih5['/Ceipj/CeipjTpjnn'].attrs['Jmaxi']
+
+def numb(pp,jj):
+    lbar = (Jmaxe +1)*pp + jj +1
+    return lbar
+
+# Get 4 dim rows
+rows_4dim= eih5['/Ceipj/rows_eT']
+# Get Cippj
+Cpj_4dim= eih5['/Ceipj/CeipjTpjnn']
+
+Cpj = np.zeros((numb(Pmaxe,Jmaxe), numb(Pmaxe,Jmaxe)))
+
+for irow in np.arange(1, numb(Pmaxe,Jmaxe)+1,1):
+    idx_of_row, = np.where(rows_4dim[0,:] == irow)
+    Cpj_ = Cpj_4dim[:, idx_of_row]
+    Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
+
+Cpj_dst = eih5_out.create_dataset("/Ceipj/CeipjT", data = Cpj)
+Cpj_dst.attrs['Pmaxe'] =  Pmaxe
+Cpj_dst.attrs['Jmaxe'] = Jmaxe
+Cpj_dst.attrs['kperp'] = eih5['/Ceipj/CeipjTpjnn'].attrs['kperp']
+
+
+### Field component #####
+
+# get dimension
+Pmaxe = eih5['/Ceipj/CeipjFpjnn'].attrs['Pmaxe']
+Jmaxe = eih5['/Ceipj/CeipjFpjnn'].attrs['Jmaxe']
+
+def numb(pp,jj):
+    lbar = (Jmaxe +1)*pp + jj +1
+    return lbar
+
+# Get 4 dim rows
+rows_4dim= eih5['/Ceipj/rows_eF']
+
+# Get Cippj
+Cpj_4dim= eih5['/Ceipj/CeipjFpjnn']
+
+Cpj = np.zeros((numb(Pmaxe, Jmaxe), numb(Pmaxe,Jmaxe)))
+
+for irow in np.arange(1, numb(Pmaxe,Jmaxe)+1,1):
+    idx_of_row, = np.where(rows_4dim[0,:] == irow)
+    Cpj_ = Cpj_4dim[:, idx_of_row]
+    Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
+
+Cpj_dst = eih5_out.create_dataset("/Ceipj/CeipjF", data = Cpj)
+Cpj_dst.attrs['Pmaxe'] =  Pmaxe
+Cpj_dst.attrs['Jmaxe'] = Jmaxe
+Cpj_dst.attrs['kperp'] = eih5['/Ceipj/CeipjFpjnn'].attrs['kperp']
+
+eih5_out.create_dataset("/files/fort.90", data = eih5['/files/fort.90'])
+
+# close file
+eih5.close()
+eih5_out.close()
+
+
+
diff --git a/python_scripts/gkcoulomb_ie.py b/python_scripts/gkcoulomb_ie.py
new file mode 100644
index 0000000..8b2a2a6
--- /dev/null
+++ b/python_scripts/gkcoulomb_ie.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+# This scripts compute the GK Coulomb collision operator from the COSOlver output files
+# Performs the summation over the last indices, i.e. (p,j,m,n,np) -> (p,j)
+# Writes the results matrices 
+
+import numpy as np
+import os
+import sys
+import h5py as h5
+
+
+# Define directories
+cwd = os.getcwd()
+
+print('Create GK Coulomb from ' + cwd)
+
+# Define h5 files 
+iecolls_file = os.path.join(cwd,'ie.h5')
+
+# Defines new h5 files
+iecolls_file_out = os.path.join(cwd,'ie.pj.h5')
+
+############################################################
+# Ion-electron collision
+print('Self ion-electron collisions....')
+
+ieh5 = h5.File(iecolls_file, 'r')
+ieh5_out = h5.File(iecolls_file_out, 'w')
+
+### Test component #####
+
+# get dimension
+Pmaxi = ieh5['/Ciepj/CiepjTpjnn'].attrs['Pmaxi']
+Jmaxi = ieh5['/Ciepj/CiepjTpjnn'].attrs['Jmaxi']
+
+def numb(pp,jj):
+    lbari = (Jmaxi +1)*pp + jj +1
+    return lbari
+
+# Get 4 dim rows
+rows_4dim= ieh5['/Ciepj/rows_iT']
+# Get Cippj
+Cpj_4dim= ieh5['/Ciepj/CiepjTpjnn']
+
+Cpj = np.zeros((numb(Pmaxi,Jmaxi), numb(Pmaxi,Jmaxi)))
+
+for irow in np.arange(1, numb(Pmaxi,Jmaxi)+1,1):
+    idx_of_row, = np.where(rows_4dim[0,:] == irow)
+    Cpj_ = Cpj_4dim[:, idx_of_row]
+    Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
+
+Cpj_dst = ieh5_out.create_dataset("/Ciepj/CiepjT", data = Cpj)
+Cpj_dst.attrs['Pmaxi'] =  Pmaxi
+Cpj_dst.attrs['Jmaxi'] = Jmaxi
+Cpj_dst.attrs['kperp'] = ieh5['/Ciepj/CiepjTpjnn'].attrs['kperp']
+
+
+### Field component #####
+
+# get dimension
+Pmaxi = ieh5['/Ciepj/CiepjFpjnn'].attrs['Pmaxi']
+Jmaxi = ieh5['/Ciepj/CiepjFpjnn'].attrs['Jmaxi']
+
+def numb(pp,jj):
+    lbari = (Jmaxi +1)*pp + jj +1
+    return lbari
+
+# Get 4 dim rows
+rows_4dim= ieh5['/Ciepj/rows_iF']
+# Get Cippj
+Cpj_4dim= ieh5['/Ciepj/CiepjFpjnn']
+
+Cpj = np.zeros((numb(Pmaxi,Jmaxi), numb(Pmaxi,Jmaxi)))
+
+for irow in np.arange(1, numb(Pmaxi,Jmaxi)+1,1):
+    idx_of_row, = np.where(rows_4dim[0,:] == irow)
+    Cpj_ = Cpj_4dim[:, idx_of_row]
+    Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
+
+Cpj_dst = ieh5_out.create_dataset("/Ciepj/CiepjF", data = Cpj)
+Cpj_dst.attrs['Pmaxi'] =  Pmaxi
+Cpj_dst.attrs['Jmaxi'] = Jmaxi
+Cpj_dst.attrs['kperp'] = ieh5['/Ciepj/CiepjFpjnn'].attrs['kperp']
+
+ieh5_out.create_dataset("/files/fort.90", data = ieh5['/files/fort.90'])
+
+# close file
+ieh5.close()
+ieh5_out.close()
diff --git a/python_scripts/gkcoulomb_ii.py b/python_scripts/gkcoulomb_ii.py
new file mode 100644
index 0000000..b3c8b38
--- /dev/null
+++ b/python_scripts/gkcoulomb_ii.py
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+# This scripts compute the GK Coulomb collision operator from the COSOlver output files
+# Performs the summation over the last indices, i.e. (p,j,m,n,np) -> (p,j)
+# Writes the results matrices 
+
+import numpy as np
+import os
+import sys
+import h5py as h5
+
+
+# Define directories
+cwd = os.getcwd()
+cwd = '/misc/bjfrei/cosolver_marconi/gk.coulomb.TEM.kperp.1.18.m.4.NFLR.8.ee'
+
+print('Create GK Coulomb from ' + cwd)
+
+# Define h5 files 
+selfcolls_file =  os.path.join( cwd,  'self.h5')
+
+# Defines new h5 files
+selfcolls_file_out =  os.path.join( cwd,  'self.pj.h5')
+
+
+############################################################
+# Self collisions
+print('Self collisions....')
+
+# Open file
+selfh5 = h5.File(selfcolls_file, 'r')
+selfh5_out = h5.File(selfcolls_file_out, 'w')
+
+# Ion-ion
+
+# get dimension
+Pmaxi = selfh5['/Caapj/Ciipjpjnn'].attrs['Pmaxi']
+Jmaxi = selfh5['/Caapj/Ciipjpjnn'].attrs['Jmaxi']
+
+def numb(pp,jj):
+    lbari = (Jmaxi +1)*pp + jj +1
+    return lbari
+
+# Get 4 dim rows
+rows_4dim= selfh5['/Caapj/rows_i']
+# Get Cippj
+Cpj_4dim=selfh5['/Caapj/Ciipjpjnn']
+#Cpj_4dim
+
+Cpj = np.zeros((numb(Pmaxi,Jmaxi), numb(Pmaxi,Jmaxi)))
+
+for irow in np.arange(1, numb(Pmaxi,Jmaxi)+1,1):
+    idx_of_row, = np.where(rows_4dim[0,:] == irow)
+    Cpj_ = Cpj_4dim[:, idx_of_row]
+    Cpj[irow-1,:] = np.sum(Cpj_,axis = 1, initial = 0)
+
+# Create output    
+Cpj_dst = selfh5_out.create_dataset("/Caapj/Ciipj", data = Cpj)
+Cpj_dst.attrs['Pmaxi'] =  Pmaxi
+Cpj_dst.attrs['Jmaxi'] = Jmaxi
+Cpj_dst.attrs['kperp'] = selfh5['/Caapj/Ceepjpjnn'].attrs['kperp']
+
+selfh5_out.create_dataset("/files/fort.90", data = selfh5['/files/fort.90'])
+
+# close file
+selfh5_out.close()
+selfh5.close()