diff --git a/.gitignore b/.gitignore
index 866cd7b..a43e1e0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,4 +7,7 @@
 *.a
 
 #test output files
-fortran/bfield_library_jdl/tests/*.out
\ No newline at end of file
+fortran/bfield_library_jdl/tests/*.out
+
+# matlab
+*.asv
\ No newline at end of file
diff --git a/matlab/bfield_library_jdl/C2_routines/c2_test_driver.m b/matlab/bfield_library_jdl/C2_routines/c2_test_driver.m
new file mode 100644
index 0000000..a16b222
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/c2_test_driver.m
@@ -0,0 +1,30 @@
+clearvars;
+
+
+%
+%  Setup
+%
+
+run_path = 'C:\Work\C2_example\example_jdl\run01\';
+grid_dir = [run_path,'ingrid\'];
+ndomain = 6;
+time_str = 'final';
+
+
+%
+%  Examples
+% 
+
+%  --> Load c2 output
+grid = read_c2_grid(grid_dir,ndomain);
+data = read_c2_data(run_path,ndomain,time_str);
+
+% --> plot grid
+plot_c2_grid(grid);
+
+% --> plot output 2d
+data_plot = data.Te2d; mytitle = 'T_e [eV]';
+plot_2d_c2(grid,data_plot,mytitle)
+
+% --> Example script to plot Te, Ti, np at OMP
+plot_omp_profiles_c2(grid,data)
\ No newline at end of file
diff --git a/matlab/bfield_library_jdl/C2_routines/plot_2d_c2.m b/matlab/bfield_library_jdl/C2_routines/plot_2d_c2.m
new file mode 100644
index 0000000..0057442
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/plot_2d_c2.m
@@ -0,0 +1,29 @@
+function plot_2d_c2(grid,data_plot,mytitle)
+if nargin < 3
+    mytitle = ' ';
+end
+
+icount = 1;
+for idomain = 1:grid.ndomain
+    for i = 1:grid.ny(idomain)-1
+        for j = 1:grid.nx(idomain)-1
+             xc(:,icount) = [grid.x2d{idomain}(i,j),grid.x2d{idomain}(i,j+1),grid.x2d{idomain}(i+1,j+1),grid.x2d{idomain}(i+1,j)];
+             yc(:,icount) = [grid.y2d{idomain}(i,j),grid.y2d{idomain}(i,j+1),grid.y2d{idomain}(i+1,j+1),grid.y2d{idomain}(i+1,j)];
+             d(icount) = data_plot{idomain}(i,j);
+             icount = icount + 1;  
+        end
+    end
+end
+
+figure; hold on; box on; 
+patch(xc(:,~isnan(d)),yc(:,~isnan(d)),d(~isnan(d)),'edgecolor','none')
+fprintf('Max value in this plot: %f\n',max(d(~isnan(d))))
+
+axis tight;
+colorbar('fontsize',14)
+xlabel('R [m]','fontsize',14)
+ylabel('Z [m]','fontsize',14)
+set(gca,'fontsize',14)
+title(mytitle);
+
+colorbar;
\ No newline at end of file
diff --git a/matlab/bfield_library_jdl/C2_routines/plot_c2_grid.m b/matlab/bfield_library_jdl/C2_routines/plot_c2_grid.m
new file mode 100644
index 0000000..4d50f4a
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/plot_c2_grid.m
@@ -0,0 +1,8 @@
+function plot_c2_grid(grid)
+
+cf = colorflipper(grid.ndomain,'jet');
+figure; hold on; box on;
+for idomain = 1:grid.ndomain
+    plot(grid.x2d{idomain},grid.y2d{idomain},'-','color',cf(idomain,:))
+    plot(grid.x2d{idomain}.',grid.y2d{idomain}.','-','color',cf(idomain,:))
+end
\ No newline at end of file
diff --git a/matlab/bfield_library_jdl/C2_routines/plot_omp_profiles_c2.m b/matlab/bfield_library_jdl/C2_routines/plot_omp_profiles_c2.m
new file mode 100644
index 0000000..cb8f522
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/plot_omp_profiles_c2.m
@@ -0,0 +1,67 @@
+function plot_omp_profiles_c2(grid,data)
+% Currently set up for d3d case and set to give max 1 point per cell
+rax = 2; rend = 2.5; zax = 0; zend = 0; 
+
+nr_prof = 200;
+
+data_plot1 = data.Te2d; fac1 = 1.e-3;
+data_plot2 = data.Ti2d; fac2 = 1.e-3;
+data_plot3 = data.np2d; fac3 = 1;
+
+% intersect line with core inner surface and SOL outer surface
+rmin = 100;
+rmax = 0;
+for iz = 1:grid.ndomain
+    [pint1,ierr1]=int_line_curve([rax,zax],[rend,zend],grid.x2d{iz}(1,:),grid.y2d{iz}(1,:));
+    if ierr1 == 0
+        rmin = min([rmin,pint1(1)]);
+        rmax = max([rmax,pint1(1)]);
+    end
+    [pint1,ierr1]=int_line_curve([rax,zax],[rend,zend],grid.x2d{iz}(grid.ny(iz),:),grid.y2d{iz}(grid.ny(iz),:));
+    if ierr1 == 0
+        rmin = min([rmin,pint1(1)]);
+        rmax = max([rmax,pint1(1)]);
+    end    
+end
+
+r_prof = linspace(rmin-1e-3,rmax+1e-3,nr_prof);
+z_prof = linspace(zax,zax,nr_prof);
+
+data_prof1 = NaN(1,nr_prof);
+data_prof2 = NaN(1,nr_prof);
+data_prof3 = NaN(1,nr_prof);
+
+ig_used = []; icount = 0;
+for irp = 1:nr_prof
+    [izc,irc,jpc,ierrc]=point2cell_c2(grid,r_prof(irp),z_prof(irp),1);
+    if ierrc == 0
+        ig = irc + (jpc-1)*(grid.ny(izc)-1) + grid.Zoffsets(izc);        
+        data_prof1(irp) = data_plot1{izc}(irc,jpc);
+        data_prof2(irp) = data_plot2{izc}(irc,jpc);
+        data_prof3(irp) = data_plot3{izc}(irc,jpc);
+        
+        if ~any(ig_used == ig)
+            icount = icount + 1;
+            r_prof_1pc(icount) = r_prof(irp);
+            z_prof_1pc(icount) = z_prof(irp);
+            data_prof1_1pc(icount) = data_plot1{izc}(irc,jpc);
+            data_prof2_1pc(icount) = data_plot2{izc}(irc,jpc);
+            data_prof3_1pc(icount) = data_plot3{izc}(irc,jpc);
+        end
+        ig_used = [ig_used,ig];
+        
+    end
+end
+
+
+figure; hold on; box on;
+plot(r_prof_1pc,data_prof1_1pc*fac1,'ro-','linewidth',2)
+plot(r_prof_1pc,data_prof2_1pc*fac2,'bo-','linewidth',2)
+plot(r_prof_1pc,data_prof3_1pc*fac3,'ko-','linewidth',2)
+xlabel('R_{OMP} [m]','fontsize',12)
+
+ax_use = axis(gca); ax_use(3) = 0; axis(ax_use);
+
+ylabel('T (keV), n (10^{19} m^{-3})','fontsize',12)
+legend('T_e','T_i','n_e');
+set(gca,'fontsize',12)
diff --git a/matlab/bfield_library_jdl/C2_routines/point2cell_c2.m b/matlab/bfield_library_jdl/C2_routines/point2cell_c2.m
new file mode 100644
index 0000000..c9808fb
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/point2cell_c2.m
@@ -0,0 +1,92 @@
+function [izc,irc,jpc,ierr]=point2cell_c2(grid,r0,z0,quiet)
+% RZ in meters
+
+if nargin < 4
+    quiet = 0;
+end
+
+debug_text = 0;
+debug_plot = 0;
+npts = length(r0);
+ierr = zeros(1,npts);
+
+for ipt = 1:npts
+    
+    % Identify zone by checking in polygon with boundary
+    inpoly0 = zeros(1,grid.ndomain);
+%     figure; hold on;
+    for iz = 1:grid.ndomain
+        nr = grid.ny(iz);
+        np = grid.nx(iz);
+        
+        rpoly0 = [grid.x2d{iz}(1:nr,1);grid.x2d{iz}(nr,2:np).';grid.x2d{iz}(nr-1:-1:1,np);grid.x2d{iz}(1,np-1:-1:1).'];
+        zpoly0 = [grid.y2d{iz}(1:nr,1);grid.y2d{iz}(nr,2:np).';grid.y2d{iz}(nr-1:-1:1,np);grid.y2d{iz}(1,np-1:-1:1).'];        
+        inpoly0(iz) = inpolygon(r0(ipt),z0(ipt),rpoly0,zpoly0);
+        
+%             plot(rpoly0,zpoly0)
+%             plot(r0(ipt),z0(ipt),'rx')
+    end
+    
+    if sum(inpoly0) ~= 1
+        if sum(inpoly0) == 0 && ~quiet
+            fprintf('Point not found in any zone!\n')
+        end
+        if sum(inpoly0) > 1
+            disp(['Point found in more than one zone!? : ',num2str(inpoly0)])
+        end
+        
+        ierr(ipt) = 1;
+        izc(ipt)=NaN;irc(ipt)=NaN;jpc(ipt)=NaN;
+        continue;
+        %     return
+    end
+    
+    %%%%%%%%%%%
+    izc(ipt) = find(inpoly0 == 1);
+    %%%%%%%%%%%
+    
+    % Identify radial and poloidal cell indices
+    nr = grid.ny(izc(ipt));
+    np = grid.nx(izc(ipt));
+    foundit = 0;
+    for ir = 1:nr-1
+        for jp = 1:np-1
+            rpoly0 = [grid.x2d{izc(ipt)}(ir,jp,ipt),grid.x2d{izc(ipt)}(ir+1,jp,ipt),grid.x2d{izc(ipt)}(ir+1,jp+1,ipt),grid.x2d{izc(ipt)}(ir,jp+1,ipt),grid.x2d{izc(ipt)}(ir,jp,ipt)];
+            zpoly0 = [grid.y2d{izc(ipt)}(ir,jp,ipt),grid.y2d{izc(ipt)}(ir+1,jp,ipt),grid.y2d{izc(ipt)}(ir+1,jp+1,ipt),grid.y2d{izc(ipt)}(ir,jp+1,ipt),grid.y2d{izc(ipt)}(ir,jp,ipt)];
+            inpoly0 = inpolygon(r0(ipt),z0(ipt),rpoly0,zpoly0);
+            
+            if inpoly0 == 1
+                foundit = 1;
+                irc(ipt) = ir;
+                jpc(ipt) = jp;
+                break;
+            end
+        end
+        if foundit == 1
+            break
+        end
+    end
+    
+    if foundit == 0
+        ierr(ipt) = 1;
+        fprintf('Point not found in any cell!\n')
+        izc(ipt)=NaN;irc(ipt)=NaN;jpc(ipt)=NaN;
+        continue;
+    end
+    
+    % Debug 
+    if debug_text        
+        fprintf('Found point [R,Z] = [%8.3f, %8.3f] at [iz,ir,jp] = [%d,%d,%d]\n',[r0(ipt),z0(ipt),izc(ipt),irc(ipt),jpc(ipt)])
+    end
+    if debug_plot
+        ZNP = izc(ipt);
+        IR = irc(ipt);
+        JP = jpc(ipt);        
+        
+        rpoly0 = [grid.x2d{ZNP}(IR,JP),grid.x2d{ZNP}(IR+1,JP),grid.x2d{ZNP}(IR+1,JP+1),grid.x2d{ZNP}(IR,JP+1),grid.x2d{ZNP}(IR,JP)];
+        zpoly0 = [grid.y2d{ZNP}(IR,JP),grid.y2d{ZNP}(IR+1,JP),grid.y2d{ZNP}(IR+1,JP+1),grid.y2d{ZNP}(IR,JP+1),grid.y2d{ZNP}(IR,JP)];
+        figure; hold on;        
+        plot(rpoly0,zpoly0,'.k--')
+        plot(r0(ipt),z0(ipt),'rx')
+    end        
+end
\ No newline at end of file
diff --git a/matlab/bfield_library_jdl/C2_routines/read_c2_data.m b/matlab/bfield_library_jdl/C2_routines/read_c2_data.m
new file mode 100644
index 0000000..635b4c7
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/read_c2_data.m
@@ -0,0 +1,41 @@
+function data = read_c2_data(run_path,ndomain,time_str)
+
+for idomain = 1:ndomain
+    istr_app = '';
+    if idomain < 100 && idomain > 9
+        istr_app = '0';
+    elseif idomain < 10
+        istr_app = '00';
+    end
+    
+    data_fname = [run_path,'2d.',time_str,'.',istr_app,num2str(idomain),'.dat'];
+    
+    fid = fopen(data_fname);
+    node_num = fscanf(fid,'%d',2);
+    num_pts = prod(node_num);
+    nx = node_num(1);
+    ny = node_num(2);
+    
+    data_tmp = fscanf(fid,'%d %d %e %e %e %e %e %e %e %e',10*num_pts);
+    data_tmp = reshape(data_tmp,10,num_pts);
+
+    xc1d = data_tmp(3,:);
+    yc1d = data_tmp(4,:);
+    Te1d = data_tmp(5,:);
+    Ti1d = data_tmp(6,:);
+    np1d = data_tmp(7,:);
+    nn1d = data_tmp(8,:);
+    u1d  = data_tmp(9,:);
+    bx1d = data_tmp(10,:);
+    
+    data.xc2d{idomain} = reshape(xc1d,ny,nx);
+    data.yc2d{idomain} = reshape(yc1d,ny,nx);
+    data.Te2d{idomain} = reshape(Te1d,ny,nx);
+    data.Ti2d{idomain} = reshape(Ti1d,ny,nx);
+    data.np2d{idomain} = reshape(np1d,ny,nx);
+    data.nn2d{idomain} = reshape(nn1d,ny,nx);
+    data.u2d{idomain}  = reshape(u1d,ny,nx);
+    data.bx2d{idomain} = reshape(bx1d,ny,nx);
+    
+    fclose(fid);
+end
diff --git a/matlab/bfield_library_jdl/C2_routines/read_c2_grid.m b/matlab/bfield_library_jdl/C2_routines/read_c2_grid.m
new file mode 100644
index 0000000..547a7a9
--- /dev/null
+++ b/matlab/bfield_library_jdl/C2_routines/read_c2_grid.m
@@ -0,0 +1,55 @@
+function grid = read_c2_grid(grid_dir,ndomain)
+
+for idomain = 1:ndomain
+    istr_app = '';
+    if idomain < 100 && idomain > 9
+        istr_app = '0';
+    elseif idomain < 10
+        istr_app = '00';
+    end
+    
+    grid_fname = [grid_dir,'grid.',istr_app,num2str(idomain),'.dat'];
+    
+    fid = fopen(grid_fname);
+    data = fscanf(fid,'%s',1);
+    if isempty(strfind((lower(data)),'domain')); error('output not in assumed order (domain)'); end
+    domain_id = fscanf(fid,'%d',1);
+    
+    data = fscanf(fid,'%s',1);
+    if isempty(strfind((lower(data)),'nodenum')); error('output not in assumed order (nodenum)'); end
+    node_num = fscanf(fid,'%d',2);
+    num_pts = prod(node_num);
+    nx = node_num(1);
+    ny = node_num(2);
+    
+    data = fscanf(fid,'%s',1);
+    if isempty(strfind((lower(data)),'connect')); error('output not in assumed order (connect)'); end
+    connect = fscanf(fid,'%d',4);
+    
+    data = fscanf(fid,'%s',1);
+    if isempty(strfind((lower(data)),'data')); error('output not in assumed order (data)'); end
+    grid_data = fscanf(fid,'%e',num_pts*5);
+    grid_data = reshape(grid_data,5,num_pts);
+    x1d = grid_data(1,:);
+    y1d = grid_data(2,:);
+    x2d = reshape(x1d,ny,nx);
+    y2d = reshape(y1d,ny,nx);
+    
+    grid.x2d{idomain} = x2d;
+    grid.y2d{idomain} = y2d;
+    fclose(fid);
+    
+    grid.nx(idomain) = nx;
+    grid.ny(idomain) = ny;
+    grid.num_pts(idomain) = num_pts;
+    
+end
+grid.ndomain = ndomain;
+
+grid.Zoffsets(1) = 0;
+for iz = 2:grid.ndomain
+    grid.Zoffsets(iz) = grid.Zoffsets(iz-1) + (grid.ny(iz-1) - 1)*(grid.nx(iz-1) - 1);
+end
+% iz = grid.nz+1;
+% Zoff_all = grid.Zoffsets(iz-1) + (grid.nr_array(iz-1) - 1)*(grid.np_array(iz-1) - 1);
+