diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 619fb239dacb577a28f7ce5ef201d9975012f708..55aaa1805d942f837b58fe1a7cb53225794ae58e 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1997,6 +1997,20 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         warning('Found negative values of psi_norm at TS positions, replacing them with 0');
         psi_norm(mask) = 0;
       end
+      params_eff = gdat_data.gdat_params;
+      params_eff.doplot=0;
+      params_eff.data_request='z_axis';
+      z_axis = gdat(gdat_data.shot,params_eff);
+      ii=iround_os(psiscatvol.dim{1},z_axis.data);
+      is_psi_plasma = psi_norm < 1;
+      is_in_lcfs = zeros(size(psi_norm));
+      for it=1:size(is_psi_plasma,2)
+        is_in_lcfs(ii(it):-1:1,it)=cumprod(is_psi_plasma(ii(it):-1:1,it));
+        is_in_lcfs(ii(it):end,it)=cumprod(is_psi_plasma(ii(it):end,it));
+      end
+      is_private_flux = is_psi_plasma & ~is_in_lcfs;
+      % at this stage set private flux "rho" to max, could use 1.3^2
+      psi_norm(is_private_flux) = max(max(psi_norm));
       rho = sqrt(psi_norm);
     else
       if gdat_params.nverbose>=1 && gdat_data.gdat_params.edge==0