a=rtable('stratification.dat',3)
z=reform(a(0,*))
lnrho=reform(a(1,*))
lnTT=reform(a(2,*))
!p.multi=[0,1,2]
z1=min(z)
z2=max(z)
;
;  new mesh
;
zz1=0. & zz2=17. & nzz=576
zz1=0. & zz2=10. & nzz=576
zz1=0. & zz2=20. & nzz=1152  ;(output file)
zz=grange(zz1,zz2,nzz)
;
;  interpolate to new mesh
;
llnrho=interpol(lnrho,z,zz)
llnTT=interpol(lnTT,z,zz)
;
;  compare the two
;
print,'AXEL, zz2=',zz2,z2
print,'AXEL, xr=',!x.range
!x.range=[(z1<zz1),(z2>zz2)]
plot_io,zz,exp(llnrho),ps=-1,yr=[1e-7,1e-3]
oplot,z,exp(lnrho),col=122
plot,zz,exp(llnTT),ps=-1,yr=[0,5000.]
oplot,z,exp(lnTT),col=122
;
;  write result
;
file='stratification.new'
wtable3,file,zz,llnrho,llnTT
END