!!! to compare one TMM run output to the observational data for OMZ
!!! two metrics, 1. total volume and volume fraction
!!!			2. area map.
!!!			3. vertical area
!!! by bencat.xuan@gmail.com
!!!   

!———————————————————————————————————————————————————————
!!! process the path, output: name of the run folder, the path to the folder, full path of the run folder.
let foldernameOrfullpath="$1"
let length = strlen(foldernameOrfullpath)
let ferretworkpath=spawn("pwd")
let ferretworkpathslash=strcat(ferretworkpath,"/")
let nameORpath = strrindex(foldernameOrfullpath,"/")
let nameORpathplus= nameORpath+1

let foldername = if nameORpath ge 1 then substring(foldernameOrfullpath,nameORpathplus,length) else foldernameOrfullpath
list foldername
let nameORpathminus = if nameORpath ge 1 then `nameORpath`-1 else 0
let tagetpath= if nameORpath ge 1 then substring(foldernameOrfullpath,1,nameORpathminus) else ferretworkpath

list foldername,tagetpath

let tagetpathslash=strcat(tagetpath, "/")
let fullpath=strcat(tagetpathslash, foldername)

!———————————————————————————————————————————————————————
!!! get the folder name
let slen=strlen(foldername)
let name1=substring(foldername,5,slen)


!———————————————————————————————————————————————————————
!!! the output files form TMM may have no depth_edge, which leads to errors in @DIN and @AVE.so here we define the right boundary for the depth boxes.

let coordin={17.5, 82.5, 177.5, 302.5, 457.5, 642.5, 857.5, 1102.5, 1377.5,1682.5, 2017.5, 2382.5, 2777.5, 3202.5, 3657.5, 4142.5, 4657.5, 5202.5,5777.5}
let bond={ 0, 50, 130, 240, 380, 550, 750, 980, 1240, 1530, 1850, 2200,2580, 2990, 3430, 3900, 4400, 4930, 5490, 6080 }
define axis/z/depth/bounds zax=coordin,bond

!———————————————————————————————————————————————————————
!!! observation data set 
use "/Users/wyao/Desktop/UVicWorking/2.9_ln_iron/dataprocessing/WOA2013/O_no3.nc"
use "/Users/wyao/Desktop/UVicWorking/2.9_ln_iron/dataprocessing/WOA2013/O_o2.nc"
use "/Users/wyao/Desktop/UVicWorking/2.9_ln_iron/dataprocessing/WOA2013/O_po4.nc"
use "/Users/wyao/Desktop/UVicWorking/2.9_ln_iron/dataprocessing/irondatabasecompile/dfe_2017_gauss_mm.nc"
use "/Volumes/BM_2017_36/MODISdata/ave_chl.nc"
use "/Volumes/BM_2017_36/PANGAEA/diaz_obs_rg_luo_2012.nc"


!———————————————————————————————————————————————————————
!!! online run 
use "/Volumes/BM_2017_36/UVicWorking/2.9_ln_iron/quickerTMMphysicOnline/tavg.08100.01.01.nc"

!———————————————————————————————————————————————————————
!!! hand tuned run
use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/NO3mm.nc"
use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/PO4mm.nc"
use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/O2mm.nc"
!use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/DFEmm.nc"
use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/O_chl.nc"
use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/O_chl_diaz.nc"
use "/Volumes/BM_2017_36/UVOKWorking/fifth_refrun1206run3000_chl/DIAZmm.nc"

!———————————————————————————————————————————————————————
!!! load tmm run data
let currentfilepath=strcat(fullpath,"/NO3mm.nc")
use `currentfilepath`
let currentfilepath=strcat(fullpath,"/PO4mm.nc")
use `currentfilepath`
let currentfilepath=strcat(fullpath,"/O2mm.nc")
use `currentfilepath`
!let currentfilepath=strcat(fullpath,"/DFEmm.nc")
!use `currentfilepath`
!let currentfilepath=strcat(fullpath,"/O_chl.nc")
!use `currentfilepath`
!let currentfilepath=strcat(fullpath,"/O_chl_diaz.nc")
!use `currentfilepath`
let currentfilepath=strcat(fullpath,"/DIAZmm.nc")
use `currentfilepath`
!———————————————————————————————————————————————————————
!!! increase memory size
set memory/size=520


let obsO2=o_o2[d=2]
let refo2=o2[d=10,gz=zax,l=@ave]
let opto2=o2[d=16,gz=zax,l=@ave]

let obso2v5=if obso2 le 0.005 then 1 else 0
let refo2v5=if refo2 le 0.005 then 1 else 0
let opto2v5=if opto2 le 0.005 then 1 else 0

let obso2v50=if obso2 le 0.05 then 1 else 0
let refo2v50=if refo2 le 0.05 then 1 else 0
let opto2v50=if opto2 le 0.05 then 1 else 0
let volume= obso2*0+1
!———————————————————————————————————————————————————————
!!! view port 1x4
go isquare2x2.jnl


go preparefigure.jnl `name1`_OMZ_metrics
!———————————————————————————————————————————————————————
!!! 
set win 1
!go set_ocean_basins.jnl O_O2 2 G_MSKHR d=7

!———————————————————————————————————————————————————————
!!! get the area map 
set vie v1
sha/lev="(-inf,0)(0,0.005,0.001)(0.005,0.05,0.005)(0.005,0.4,0.05)" refo2[k=@min]
contour/level="(0.005)(0.050)"/color=black/ov obso2[k=@min]  
go land
set vie v2
sha/lev="(-inf,0)(0,0.005,0.001)(0.005,0.05,0.005)(0.005,0.4,0.05)" opto2[k=@min]
contour/level="(0.005)(0.050)"/color=black/ov obso2[k=@min]  
go land

set vie v3
plot obso2v5[x=@din,y=@din],refo2v5[x=@din,y=@din],opto2v5[x=@din,y=@din]

set vie v4
plot obso2v50[x=@din,y=@din],refo2v50[x=@din,y=@din],opto2v50[x=@din,y=@din]

!plot obso2v5[x=@din,y=@din]/volume[x=@din,y=@din,z=@din],refo2v5[x=@din,y=@din]/volume[x=@din,y=@din,z=@din],opto2v5[x=@din,y=@din]/volume[x=@din,y=@din,z=@din]

go finalizefigure_psthicken


list/file=`name1`.txt  opto2v5[x=@din,y=@din,z=@din],opto2v50[x=@din,y=@din,z=@din],opto2v5[x=@din,y=@din,z=@din]/obso2v5[x=@din,y=@din,z=@din],opto2v50[x=@din,y=@din,z=@din]/obso2v50[x=@din,y=@din,z=@din]

let obso2v5min =if obso2 le 0.0025 then 1 else 0
let obso2v5max =if obso2 le 0.0075 then 1 else 0
let obso2v50min=if obso2 le 0.040 then 1 else 0
let obso2v50max=if obso2 le 0.060 then 1 else 0


list refo2v5[x=@din,y=@din,z=@din],refo2v50[x=@din,y=@din,z=@din],refo2v5[x=@din,y=@din,z=@din]/obso2v5[x=@din,y=@din,z=@din],refo2v50[x=@din,y=@din,z=@din]/obso2v50[x=@din,y=@din,z=@din]
list obso2v5min[x=@din,y=@din,z=@din]/obso2v5[x=@din,y=@din,z=@din],obso2v5max[x=@din,y=@din,z=@din]/obso2v5[x=@din,y=@din,z=@din],obso2v50min[x=@din,y=@din,z=@din]/obso2v50[x=@din,y=@din,z=@din],obso2v50max[x=@din,y=@din,z=@din]/obso2v50[x=@din,y=@din,z=@din]
list obso2v5min[x=@din,y=@din,z=@din],obso2v5max[x=@din,y=@din,z=@din],obso2v50min[x=@din,y=@din,z=@din],obso2v50max[x=@din,y=@din,z=@din]