diff --git a/MAM.f b/MAM.f
index 70e3184..e6cb482 100644
--- a/MAM.f
+++ b/MAM.f
@@ -896,15 +896,15 @@ subroutine mamposst(di,ve,eve,rso,vso,npg1)
       call sigmaeva(npar,sigma)
       !format(8(f13.5,2x),i2)
       if(nlens.eq.0) then
-      WRITE(iu60,'(9(a13,2x))') 'r200','rc','rs','cbe','A1','A2',
-     & 'cbe0','-log(P)', 'anisotropy'
-      WRITE(iu60,'(9(a13,2x))') '[Mpc]','[Mpc]','[Mpc]','-','-','-',
-     & '-','-', '-' 
+      WRITE(iu60,'(a,a12,2x,8(a13,2x))') '#','r200','rc','rs','cbe',
+     & 'A1','A2','cbe0','-log(P)', 'anisotropy'
+      WRITE(iu60,'(a,a12,2x,8(a13,2x))') '#','[Mpc]','[Mpc]','[Mpc]',
+     & '-','-','-','-','-', '-'
       else
-      WRITE(iu60,'(9(a13,2x))') 'r200','rc','rs','cbe','A1','A2',
-     & 'cbe0','-log(P*Plens)', 'anisotropy'
-      WRITE(iu60,'(9(a13,2x))') '[Mpc]','[Mpc]','[Mpc]','-','-','-',
-     & '-','-', '-'      
+      WRITE(iu60,'(a1,a12,2x,8(a13,2x))') '#','r200','rc','rs','cbe',
+     & 'A1','A2','cbe0','-log(P*Plens)', 'anisotropy'
+      WRITE(iu60,'(a,a12,2x,8(a13,2x))') '#','[Mpc]','[Mpc]','[Mpc]',
+     & '-','-','-','-','-', '-'
       endif
       write(iu60, *) ''
       if (nmcmc.eq.1) then  !Start the MCMC run 
@@ -2582,8 +2582,8 @@ subroutine mamposst(di,ve,eve,rso,vso,npg1)
 
  737  continue
       sigma0=0.
-      write(iu30,'(a4,1x,a13,1x,a9)') 'N(R)','rc','al'
-      write(iu30,'(a4,1x,a11,1x,a9)') '-','[Mpc]','-'  
+      write(iu30,'(a2,a4,1x,a13,1x,a9)') '# ','N(R)','rc','al'
+      write(iu30,'(a2,a4,1x,a13,1x,a9)') '# ','-','[Mpc]','-'
       write(iu30,*) ''    
       write(iu30,330) knfit,rc,al !sigma0,
       
@@ -2596,8 +2596,8 @@ subroutine mamposst(di,ve,eve,rso,vso,npg1)
 
       call vmaxlik(nfv,xfv,fml)
       !header
-      write(iu20,'(3(a13,2x))') 'Ri', 'vz,i','p(Ri,vz,i)' 
-      write(iu20,'(3(a13,2x))') '[Mpc]', '[km/s]','-'
+      write(iu20,'(a1,a12,2x,2(a13,2x))') '#','Ri', 'vz,i','p(Ri,vz,i)'
+      write(iu20,'(a1,a12,2x,2(a13,2x))') '#','[Mpc]', '[km/s]','-'
       write(iu20,*) ''
       do jp=1,npv
          write(iu20,'(3(f13.5,2x))') rpv(jp),vpv(jp),pv(jp)
diff --git a/gomamposstoptS.f b/gomamposstoptS.f
index 7f108f5..69dc3a9 100644
--- a/gomamposstoptS.f
+++ b/gomamposstoptS.f
@@ -851,22 +851,26 @@ program gomamposstoptS
       open(70,file=fsb,status='unknown')
       open(80,file=fsf,status='unknown')
 
-      write(30,*) 'record of N(R) model used'
+      write(30,'(a)') '# Record of N(R) model used'
              
-      write(40,*) 'Binned project N(R) profile'
+      write(40,'(a)') '# Binned project N(R) profile'
       write(40,*)  
       
-      write(50,*) 'Best fit project N(R) profile'
-      write(50,*) 'R[Mpc], N(R)' 
+      write(50,'(a)') '# Best fit project N(R) profile'
+      write(50,'(a1,a25,2x,a19)') '#', 'R', 'N(R)'
+      write(50,'(a1,a25)') '#', '[Mpc]'
        
-      WRITE(60,*) 'Table containing parameters and -log(P)'
-      WRITE(60,*) ''  
+      WRITE(60,'(a)') '# Table containing parameters and -log(P)'
+      WRITE(60,'(a)') '#'
         
-      write(70,*) 'Binned velocity dispersion profile'
-      write(70,*) 'R[Mpc], VDP [km/s], err_up[km/s] err_low [km/s]'
+      write(70,'(a)') '# Binned velocity dispersion profile'
+      write(70,'(a1,a13,a11,a11,a11)') '#','R','VDP','err_up','err_low'
+      write(70,'(a1,a13,a11,a11,a11)') '#','[Mpc]','[km/s]','[km/s]',
+     &'[km/s]'
       
-      write(80,*) 'Fitted velocity dispersion profile'
-      write(80,*) 'R[Mpc], VDP [km/s]'      
+      write(80,'(a)') '# Fitted velocity dispersion profile'
+      write(80,'(a1,a20,a26)') '#', 'R', 'VDP'
+      write(80,'(a1,a20,a26)') '#', '[Mpc]', '[km/s]'
       
 c     read system properties
 
@@ -958,8 +962,11 @@ program gomamposstoptS
       write(*,*) ' '
 
 c     output a binned number density profile N(R)
-      write(40,*) 'Output a binned number density profile N(R)'
-      write(40,*) 'rbin [Mpc], density [1/Mpc^3], err [1/Mpc^3]'
+      write(40,'(a)') '# Output a binned number density profile N(R)'
+      write(40,'(a1,a20,a26,a26)') '#', 'rbin',
+     &'density','err'
+      write(40,'(a1,a20,a26,a26)') '#', '[Mpc]',
+     &'[1/Mpc^3]','[1/Mpc^3]'
       ibwt=0
       nbins=int(dsqrt(dfloat(nga)))
       do j=1,nga