diff --git a/crpptbx/TCV/loadTCVdata.m b/crpptbx/TCV/loadTCVdata.m
index f2ac00b15e9d5542fb8c4f8e02824e8bb0b3f22e..3f37051503d574a5cd1a81439e7400e6fa1afa9d 100644
--- a/crpptbx/TCV/loadTCVdata.m
+++ b/crpptbx/TCV/loadTCVdata.m
@@ -21,6 +21,7 @@
 % 'deltabot', 'triangbot'[_2,_3] =  edge lower (bottom) triangularity vs time
 % 'j_tor'[_2,_3] =  J_TOR vs (R,Z,time)
 % 'neint' =  line-integrated electron density [m/m^3]
+% 'nel' =  line-averaged electron density [1/m^3]
 % 'ne'= ne raw profile on (z,t). ADD error bars in .std
 % 'te'= Te raw profile on (z,t). ADD error bars in .std
 % 'nerho'[_2,_3]= ne profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std
@@ -261,12 +262,12 @@ liuqeFBTE=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\result
 
 % keywords which do not have FBTE equivalence and are returned empty
 noFBTE=[{'ne'} {'te'} {'nerho'} {'terho'} {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
-	{'neft'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'neint'} ...
+	{'neft'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'neint'} {'nel'} ...
         {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'Ti'} {'Tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}];
 
 % all keywords and corresponding case to run below
 TCVkeywrdall=[{'Ip'} {'B0'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
-      {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
+      {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'nel'} ...
       {'ne'} {'te'} {'nerho'} {'terho'}  {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
       {'neft'} {'teft'} {'neftav'} {'teftav'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'}  ...
       {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'pgyro'} {'jtor'} {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'Ti'} {'Tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}];
@@ -285,6 +286,7 @@ TCVsig.idelta=strmatch('delta',TCVkeywrdall,'exact');
 TCVsig.ideltatop=strmatch('deltatop',TCVkeywrdall,'exact');
 TCVsig.ideltabot=strmatch('deltabot',TCVkeywrdall,'exact');
 TCVsig.ineint=strmatch('neint',TCVkeywrdall,'exact');
+TCVsig.inel=strmatch('nel',TCVkeywrdall,'exact');
 TCVsig.ine=strmatch('ne',TCVkeywrdall,'exact');
 TCVsig.ite=strmatch('te',TCVkeywrdall,'exact');
 TCVsig.ine_edge=strmatch('ne_edge',TCVkeywrdall,'exact');
@@ -386,6 +388,7 @@ TCVsiglocation(TCVsig.idelta)={'\results::delta_edge'};
 TCVsiglocation(TCVsig.ideltatop)={'\results::delta_ed_top'};
 TCVsiglocation(TCVsig.ideltabot)={'\results::delta_ed_bot'};
 TCVsiglocation(TCVsig.ineint)={'\results::fir:lin_int_dens'};
+TCVsiglocation(TCVsig.inel)={'\results::fir:n_average'};
 %TCVsiglocation(TCVsig.iprofnerho)={'\results::THOMSON.PROFILES.AUTO:ne'};
 %TCVsiglocation(TCVsig.iprofterho)={'\results::THOMSON.PROFILES.AUTO:te'};
 TCVsiglocation(TCVsig.ineft)={'\results::proffit.local_time:neft_abs'}; TCVsigtimeindx(TCVsig.ineft)=2;
diff --git a/crpptbx/gdat.m b/crpptbx/gdat.m
index 597e04cb169912690330683036cba86a232c228d..b925e8a7ac7ae902435fdf17514fb7224e2ead01 100644
--- a/crpptbx/gdat.m
+++ b/crpptbx/gdat.m
@@ -141,21 +141,22 @@ else
 end
 
 % PLOT DATA (if required)
+if doplot~=0; set_defaults_matlab; end
 if doplot==1 & length(trace.data)>1 & ~ischar(trace.data)
   try
     figure;zoom on
     if length(size(trace.data))<=2
-      plot(trace.t,trace.data);
+      hhh=plot(trace.t,trace.data);
       ylabel(data_type)
     else
       for idim=1:length(trace.dim)
         if length(trace.t)==length(trace.dim{idim}); idim_t=idim; end
       end
       if idim_t<=2
-        plot(trace.t,trace.data(:,:,floor(end/2)));
+        hhh=plot(trace.t,trace.data(:,:,floor(end/2)));
         ylabel([data_type '(:,:,floor(end/2))'])
       elseif idim_t==3;
-        plot(trace.t,reshape(trace.data(:,floor(end/2),:),length(trace.dim{1}),length(trace.t)));
+        hhh=plot(trace.t,reshape(trace.data(:,floor(end/2),:),length(trace.dim{1}),length(trace.t)));
         ylabel([data_type '(:,floor(end/2),:)'])
       end
     end
@@ -168,12 +169,19 @@ if doplot==1 & length(trace.data)>1 & ~ischar(trace.data)
 elseif doplot==-1
   try
     hold on
+    child_h=get(gca,'child');
+    nbplot=length(child_h);
     if length(size(trace.data))<=2
-      plot(trace.t,trace.data,'r');
+      hhh=plotos(trace.t,trace.data,'-',[],[],colos(mod(nbplot,12)+1,:));
     else
-      plot(trace.t,trace.data(:,:,1),'--');
+      hhh=plot(trace.t,trace.data(:,:,1),'--');
     end
   catch
     disp(' error in plotting part, most probably because could not guess time dimension correctly. To check')
   end
 end
+if exist('hhh') && ishandle(hhh); set(hhh(1),'Tag',['gdat: ' num2str(shot)]); end
+h2=findobj(gca,'-regexp','Tag','gdat:*');
+if ~isempty(h2);
+  legend(get(sort(h2),'Tag'));
+end