diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index d721f86c7ac01fa69257c421dbea7582dbf136d4..d9f2a069b5c6eec90bc3336a903e833a5aa7b2a1 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -2356,7 +2356,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
    case {'powers'}
     % note: same time array for all main, ec, ohm, nbi, ...
     % At this stage fill just ech, later add nbi
-    sources_avail = {'ohm','ec','nbi','rad'}; % note should allow ech, nbh, ohmic in parameter sources
+    sources_avail = {'ohm','ec','nbi','rad', 'nb2'}; % note should allow ech, nbh, ohmic in parameter sources
     % create empty structures for all, so in return one always have same substructres
     for i=1:length(sources_avail)
       gdat_data.(sources_avail{i}).data = [];
@@ -2531,6 +2531,36 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.label{end+1}=gdat_data.rad.label;
       end
     end
+    
+    %
+    if any(strmatch('nb2',gdat_data.gdat_params.source))
+    % NB2  
+    nodenameeff = '\results::NB2:POWR_TCV';
+      nb2_data_tdi = tdi(nodenameeff);
+      if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
+        nbi_neutral_power_tot = nb2_data_tdi.data.*1e6; % in W
+        ij = nbi_neutral_power_tot<100;
+        nbi_neutral_power_tot(ij) = 0.;
+        gdat_data.nbi.data = nbi_neutral_power_tot; % at this stage p_gyro is in kW'
+        gdat_data.nbi.units = 'W';
+        gdat_data.nbi.dim{1}=nb2_data_tdi.dim{1};
+        gdat_data.nbi.dimunits{1}=nb2_data_tdi.dimunits{1};
+        gdat_data.nbi.t=gdat_data.nbi.dim{1};
+        gdat_data.nbi.x=[];
+        gdat_data.nbi.data_fullpath=[nodenameeff];
+        gdat_data.nbi.label='P_{nb2}';
+        gdat_data.nbi.help = nb2_data_tdi.help;
+        nodenameeff2 = '\results::nb2:energy';
+        gdat_data.nbi.data_fullpath=[nodenameeff ' *1e6 for W; and ' nodenameeff2 ' *1e3 for eV'];
+        nbh_energy_data_tdi = tdi(nodenameeff2);
+        gdat_data.nbi.energy = nbh_energy_data_tdi.data*1e3; % for eV
+        % add to main with linear interpolation and 0 for extrapolated values
+        ij=[isfinite(gdat_data.nbi.data)];
+        gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t(ij),gdat_data.nbi.data(ij),gdat_data.t);
+        gdat_data.x(end+1) = size(gdat_data.data,2);
+        gdat_data.label{end+1}=gdat_data.nbi.label;
+      end
+    end
     % add all to last index of .data(:,i)
     ij = setdiff([1:size(gdat_data.data,2)],index_rad);
     gdat_data.data(:,end+1) = sum(gdat_data.data(:,ij),2);