!  set up the initial plotting
set window/aspect=1.5 1
set memory/size=200

!set data "../output/t15_1880_2100_heat.nc" !../output/runR.nc"
set data "./t15_1880_2100_heat.nc" !../output/runR.nc"
!use "/pillar/matear/masks/mark2/model_masks.nc"
use "model_masks.nc"

let gtc=12*1e-15
let gto=gtc/1.1
let gtn=gto*0.2095/.7808
let gta=gtn*500/14.
define axis/t=31:250:1 tpart
define axis/t=1781:2101:1 tyr0
define axis/t=1881:2101:1 tyr1

! model T values are potential temperature
let zref=2000
!let zref=0
let/title=Pot pot=theta_fo(trac12,trac11,z[g=trac12],zref)
!let/title=Pot pot=theta_fo(trac12,trac11,0,zref)
let/title=density densa=rho_un(trac12,pot,zref)-1000
if `zref eq 0` then 
 define axis/z=22:28:.02/unit=density/depth sdepth
elif `zref eq 2000 ` then  
 define axis/z=26:38:.1/unit=density/depth sdepth
endif

! define some useful variables
let sal_den=zaxreplace(trac12,densa,z[gz=sdepth])
let tem_den=ZAXREPLACE(TRAC11,DENSA,Z[GZ=SDEPTH])
let oxy_den= ZAXREPLACE(TRAC13,DENSA,Z[GZ=SDEPTH])
let pho_den=ZAXREPLACE(TRAC16,DENSA,Z[GZ=SDEPTH])