diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index 40fd682a700b66cc1293b2574903a6bc9a2f7aa4..96aba49eb2d1e320d8b1de9fc4f368563db2f2fc 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -1108,83 +1108,153 @@ 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
-
-    % EC
-    nodenameeff='\results::toray.input:p_gyro';
-    tracetdi=tdi(nodenameeff);
-    if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
-      if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
-      gdat_data.ec.data = [];
-      gdat_data.ec.units = [];
-      gdat_data.ec.dim=[];
-      gdat_data.ec.dimunits=[];
-      gdat_data.ec.t=[];
-      gdat_data.ec.x=[];
-      gdat_data.ec.data_fullpath=[];
-      gdat_data.ec.label='';
+    sources_avail = {,'ohm','ec','nbi'}; % 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 = [];
+      gdat_data.(sources_avail{i}).units = [];
+      gdat_data.(sources_avail{i}).dim=[];
+      gdat_data.(sources_avail{i}).dimunits=[];
+      gdat_data.(sources_avail{i}).t=[];
+      gdat_data.(sources_avail{i}).x=[];
+      gdat_data.(sources_avail{i}).data_fullpath=[];
+      gdat_data.(sources_avail{i}).label=[];
+    end
+    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
+      gdat_data.gdat_params.source = sources_avail;
+    elseif ~iscell(gdat_data.gdat_params.source)
+      if ischar(gdat_data.gdat_params.source)
+        if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail)))
+          if (gdat_params.nverbose>=1)
+            warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
+          end
+          return
+        else
+          gdat_data.gdat_params.source = {lower(gdat_data.gdat_params.source)};
+        end
+      else
+        if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
+        return
+      end
     else
-      gdat_data.ec.data = tracetdi.data*1e3; % at this stage p_gyro is in kW'
-      gdat_data.ec.units = 'W';
-      gdat_data.ec.dim=tracetdi.dim;
-      gdat_data.ec.dimunits=tracetdi.dimunits;
-      gdat_data.ec.t=tracetdi.dim{1};
-      gdat_data.ec.x=tracetdi.dim{2};
-      gdat_data.ec.data_fullpath=[nodenameeff];
-      gdat_data.ec.label='P_{EC}';
-      % set ec time as reference
-      gdat_data.t = gdat_data.ec.t;
-      gdat_data.dim{1} = gdat_data.t;
-      gdat_data.dimunits{1} = 's';
-      gdat_data.units = 'W';
+      for i=1:length(gdat_data.gdat_params.source)
+        gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
+        if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
+          if gdat_data.gdat_params.nverbose>=1
+            warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
+          end
+        end
+      end
     end
-
-    % Ohmic
-    gdat_data.ohm.data = [];
-    gdat_data.ohm.units = '';
-    gdat_data.ohm.dim=gdat_data.dim;
-    gdat_data.ohm.dimunits=gdat_data.dimunits;
-    gdat_data.ohm.t=[];
-    gdat_data.ohm.x=[];
-    gdat_data.ohm.data_fullpath=[];
-    gdat_data.ohm.label='';
-    % get ohmic power simply from vloop*Ip (minus sign for TCV)
-    ip=gdat([],'ip');
-    vloop=gdat([],'vloop');
-    tension = -1e5;
-    if isempty(gdat_data.t)
-      gdat_data.t = vloop.t;
-      gdat_data.dim{1} = gdat_data.t;
-      gdat_data.ohm.data = vloop.data;
+    % always start from ohmic so can use this time as base time since should yield full shot
+    % get ohmic power simply from vloop*Ip (minus sign for TCV), neglect dWp/dt could add it later, see chie_tcv_to_nodes
+    % thus should take it from conf if present
+    
+    mdsopen(shot);
+    ptot_ohm = tdi('\results::conf:ptot_ohm');
+    if ~isempty(ptot_ohm.data) && ~ischar(ptot_ohm.data) && ~isempty(ptot_ohm.dim)
+      gdat_data.ohm.data = ptot_ohm.data;
+      gdat_data.ohm.t = ptot_ohm.dim{1};
+      gdat_data.ohm.dim = ptot_ohm.dim;
+      gdat_data.ohm.dimunits = ptot_ohm.dimunits;
+      gdat_data.ohm.units = ptot_ohm.units;
+      gdat_data.ohm.data_fullpath = '\results::conf:ptot_ohm';
+      gdat_data.ohm.help = ptot_ohm.help;
+    else
+      ip=gdat([],'ip');
+      vloop=gdat([],'vloop');
+      gdat_data.ohm.t = vloop.t;
+      gdat_data.ohm.dim{1} = gdat_data.t;
       gdat_data.dimunits{1} = 's';
       gdat_data.units = 'W';
-    else
-      vloop_smooth=interpos(vloop.t,vloop.data,gdat_data.t,tension);
-      ip_t = interp1(ip.t,ip.data,gdat_data.t);
-      gdat_data.ohm.data = -vloop_smooth.*ip_t;
-    end
-    gdat_data.ohm.units = gdat_data.units;
-    gdat_data.ohm.dim=gdat_data.dim;
-    gdat_data.ohm.dimunits=gdat_data.dimunits;
-    gdat_data.ohm.t=gdat_data.t;
+      tension = -1e5;
+      vloop_smooth=interpos(vloop.t,vloop.data,gdat_data.ohm.t,tension);
+      ip_t = interp1(ip.t,ip.data,gdat_data.ohm.t);
+      gdat_data.ohm.data = -vloop_smooth.*ip_t; % TCV has wrong sign for Vloop
+      gdat_data.ohm.raw_data = -vloop.data.*ip_t;
+      gdat_data.ohm.data_fullpath = 'from vloop*Ip, smoothed vloop in data, unsmoothed in raw_data';
+      gdat_data.ohm.data_fullpath=['from -vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat, in .data, unsmoothed in .raw_data'];
+    end
     gdat_data.ohm.x=[];
-    gdat_data.ohm.data_fullpath=['-vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat'];
-    gdat_data.ohm.label='P_{OHM}';
-    
-    % total power from each and total
-    gdat_data.data(:,1) = gdat_data.ohm.data;
-    if ~isempty(gdat_data.ec.data)
-      gdat_data.data(:,2) = gdat_data.ec.data(:,10);
-      gdat_data.data(:,3) = gdat_data.ec.data(:,10) + gdat_data.ohm.data;
-      gdat_data.dim{2} = [1:3];
-      gdat_data.dimunits{2} = 'Pohm;Pec;Ptot';
-      gdat_data.data_fullpath=['tot power from EC and ohm'];
-      gdat_data.label = 'P_{ohm};P_{EC};P_{tot}';
-    else
-      gdat_data.dim{2} = [1];
-      gdat_data.dimunits{2} = 'Pohm=Ptot';
-      gdat_data.data_fullpath=['tot power from ohm'];
-      gdat_data.label = 'P_{tot}=P_{ohm}';
+    gdat_data.ohm.label='P_{ohm}';
+    %
+    % add each source in main.data, on ohm time array
+    gdat_data.t = linspace(gdat_data.ohm.t(1),gdat_data.ohm.t(end),floor(1e4.*(gdat_data.ohm.t(end)-gdat_data.ohm.t(1))))';
+    gdat_data.data(:,1) = interpos(-21,gdat_data.ohm.t,gdat_data.ohm.data,gdat_data.t);
+    gdat_data.dim{1} = gdat_data.t;
+    gdat_data.x(1) = 1;
+    gdat_data.label={'P_{ohm}'};
+    gdat_data.units = 'W';
+    %
+    if any(strmatch('ec',gdat_data.gdat_params.source))
+      % EC
+      nodenameeff='\results::toray.input:p_gyro';
+      tracetdi=tdi(nodenameeff);
+      if isempty(tracetdi.data) || isempty(tracetdi.dim) || ischar(tracetdi.data)
+        if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+      else
+        gdat_data.ec.data = tracetdi.data*1e3; % at this stage p_gyro is in kW'
+        gdat_data.ec.units = 'W';
+        gdat_data.ec.dim=tracetdi.dim;
+        gdat_data.ec.dimunits=tracetdi.dimunits;
+        gdat_data.ec.t=tracetdi.dim{1};
+        gdat_data.ec.x=tracetdi.dim{2};
+        gdat_data.ec.data_fullpath=[nodenameeff];
+        gdat_data.ec.label='P_{EC}';
+        gdat_data.ec.help = tracetdi.help;
+        % add to main with linear interpolation and 0 for extrapolated values
+        gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t);
+        gdat_data.x(end+1) = size(gdat_data.data,2);
+        gdat_data.label{end+1}=gdat_data.ec.label;
+      end
+    end
+    %
+    if any(strmatch('nb',gdat_data.gdat_params.source))
+      NBH_in_TCV = 0;
+      if shot >= 51641
+        % NBI
+        moderemote = mdsdata('data(\VSYSTEM::TCV_NBHDB_B["NBI:MODEREMOTE_RB"])');
+        lcs_mode   = mdsdata('data(\VSYSTEM::TCV_NBHDB_I["NBI:LCS_MODE"])');
+        % NBH in TCV equiv moderemote='OK' AND lcs_mode = 9
+        NBH_in_TCV = strcmpi(strtrim(moderemote),'OK') && lcs_mode == 9;
+      else
+        % Nodes used in previous block only exist outside of Vista for shots after 51641
+        if any(shot == [51458 51459 51460 51461 51463 51465 51470 51472 ... % 29.JAN.2016
+                        51628 51629 51631 51632 51633 ...                   % 09.FEB.2016
+                        51639 51640 ... 51641                               % 10.FEB.2016
+                       ])
+          NBH_in_TCV = 1;
+        end
+      end
+      if NBH_in_TCV
+        nodenameeff = '\ATLAS::NBH.DATA.MAIN_ADC:DATA';
+        nbh_data_tdi = tdi(nodenameeff);
+        if ~isempty(nbh_data_tdi.data) && ~ischar(nbh_data_tdi.data) && ~isempty(nbh_data_tdi.dim)
+          index_for_tot_pow = 37;
+          nbi_neutral_power_tot = nbh_data_tdi.data(:,index_for_tot_pow).*1e6; % in W
+          nbi_neutral_power_tot = max(nbi_neutral_power_tot,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}=nbh_data_tdi.dim{1};
+          gdat_data.nbi.dimunits{1}=nbh_data_tdi.dimunits{1};
+          gdat_data.nbi.t=gdat_data.nbi.dim{1};
+          gdat_data.nbi.x=[];
+          gdat_data.nbi.data_fullpath=[nodenameeff ', index ' num2str(index_for_tot_pow) ', CLC NEUTRAL POWER'];
+          gdat_data.nbi.label='P_{nbi}';
+          gdat_data.nbi.help = nbh_data_tdi.help;
+          % add to main with linear interpolation and 0 for extrapolated values
+          gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data,gdat_data.t);
+          gdat_data.x(end+1) = size(gdat_data.data,2);
+          gdat_data.label{end+1}=gdat_data.nbi.label;
+        end
+      end
     end
+    % add all to last index of .data(:,i)
+    gdat_data.data(:,end+1) = sum(gdat_data.data,2);
+    gdat_data.x(end+1) = size(gdat_data.data,2);
+    gdat_data.label{end+1}='total';
+    gdat_data.dim{2} = gdat_data.x;
+    gdat_data.dimunits = {'s', 'index for each source + total'};
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'q_rho'}
diff --git a/crpptbx/TCV/tcv_help_parameters.m b/crpptbx/TCV/tcv_help_parameters.m
index ef6e287a370514610812fd851765ed2761a09b03..472eab63ddafe635c0f9b969b929f59b5e709912 100644
--- a/crpptbx/TCV/tcv_help_parameters.m
+++ b/crpptbx/TCV/tcv_help_parameters.m
@@ -38,8 +38,10 @@ help_struct_all.edge = '0 (default), 1 to get edge Thomson values';
 help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type';
 help_struct_all.fit_type = 'provenance of fitted profiles ''conf'' (default) from conf nodes, ''avg'' or ''local'' for resp. proffit: nodes';
 help_struct_all.trialindx = 'value of trialindx desired to get a specific one when relevant, otherwise gets ok_trialindx, that is 1 usually';
-help_struct_all.source = ['cxrs: [1 2 3] (default systems); sxr: main source: ''MPX'' (default) or ''XTOMO'', case insensitive' ...
-                    '; mhd request: ''23'':23 LFS/HFS (default), ''23full'': 23cm sector 3 and 11, ''0'':z=0 LFS/HFS, ''0full'': 0cm sector 3 and 11'];
+help_struct_all.source = sprintf('%s\n%s\n%s\n%s','cxrs: [1 2 3] (default systems);', ...
+          'sxr: main source: ''MPX'' (default) or ''XTOMO'', case insensitive', ...
+          'mhd request: ''23'':23 LFS/HFS (default), ''23full'': 23cm sector 3 and 11, ''0'':z=0 LFS/HFS, ''0full'': 0cm sector 3 and 11', ...
+          'powers: ohmic in any case + ''ec'', ''nbi''');
 help_struct_all.camera = ['sxr: for MPX: ''central'', ''top'' (default), ''bottom'' or ''both'' ; ' ...
                     ' for XTOMO: ''central'' (a central chord only), defaults if empty, [1 3 5] if only camera 1, 3 and 5 are desired'];
 help_struct_all.freq = '''slow'', default, lower sampling; ''fast'' full samples for both mpx and xtomo';