!script by bencat.xuan@gmail.com


define axis/t=0:`1-1/12`:`1/12` monthax



!!!!! load in data from 3 calibrations
use ../model_and_output_data/noIron_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc	!1 noFe
use ../model_and_output_data/ironMask_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc	!2 FeMAsk
use ../model_and_output_data/ironDyn1_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc !


set mem/size=2048

use "../model_and_output_data/ironDyn1_spinup_fct_drift/data/O_fe_dissolved.nc" !4

let noFe_abio=0.6979944957

let FeMask_kfe=5.020096e-5
let FeMask_kfe_D=0.0004060349
let FeMask_abio=0.6

let FeDyn1_kfe=0.0005845
let FeDyn1_kfe_D=0.0003765222
let FeDyn1_abio=0.6

let FeDyn2_kfe=0.000352
let FeDyn2_kfe_D=0.0005815124
let FeDyn2_abio=0.6

let FeDyn0_kfe=0.0004
let FeDyn0_kfe_D=0.0001
let FeDyn0_abio=0.6

let dummy = compressl(O_no3[d=1,l=1,k=1]*0.)
def sym oper = l=1,k=1,gxyz=dummy
def sym oper2 = l=300,k=1,gxyz=dummy


let k1n=0.7
let redptn=1/16
let k1p=k1n*redptn
let bbio = 1.066  ! b
let cbio = 1.0  ! c [1/deg_C]

let tempNormal=10

let bctNormal=bbio^(cbio*10) ! let 10 degree C as the standard growth temp for phyt
let bctNormal_D=bbio^(cbio*18) ! let 18 degree C as the standard growth temp for diaz

! calculation
!!!!!!!!!!!!!!!!!!!!!temp
let NoFe_bct=compressl(bbio^(cbio*o_temp[d=1,($oper)]))
let FeMask_bct=compressl(bbio^(cbio*o_temp[d=2,($oper)]))
let FeDyn1_bct=compressl(bbio^(cbio*o_temp[d=3,($oper)]))


let NoFe_templim=NoFe_bct/bctNormal
let NoFe_templim_D=NoFe_bct/bctNormal_D
let FeMask_templim=FeMask_bct/bctNormal
let FeMask_templim_D=FeMask_bct/bctNormal_D
let FeDyn1_templim=FeDyn1_bct/bctNormal
let FeDyn1_templim_D=FeDyn1_bct/bctNormal_D
let FeDyn0_templim=FeDyn0_bct/bctNormal
let FeDyn0_templim_D=FeDyn0_bct/bctNormal_D




let xyll=g_latt[d=1]

! parameters that don't change between simulations
let pi = 3.1415
let radian = 360/(2*pi)
let relyr = T[gt=monthax]+0*x[gx=xyll]+0*y[gy=xyll]
let declin = SIN((MOD(relyr,1.) - 0.22)*2.*pi)*0.4 
let rctheta1 = max(-1.5, min(1.5, g_latt[d=1]/radian - declin))
let kw = 0.04*1e-2
let rctheta = kw/(1. - (1. - cos(rctheta1)^2.)/1.33^2.)^0.5
let dayfrac1 = MIN( 1., -1*TAN(g_latt[d=1]/radian)*TAN(declin))
let dayfrac = MAX(1e-12, ACOS(MAX(-1., dayfrac1))/pi)
let dayfracmean=compressl(dayfrac[l=@ave])

!!!!!!!!!!!!!!!!!!!!!!!light
let noFe_lightlim=compressl(o_avej[d=1,($oper)]*3600*24/nofe_abio/NoFe_bct/dayfracmean)*1000
let noFe_lightlim_D=compressl(o_avej_d[d=1,($oper)]*3600*24/nofe_abio/0.5/(NoFe_bct-2.6)/dayfracmean)*1000
let FeMask_lightlim=compressl(o_avej[d=2,($oper)]*3600*24/femask_abio/FeMask_bct/dayfracmean/FeMask_felim)*1000 !/o_avej[d=2,l=@ave,x=@max,y=@max,k=1])!
let FeMask_lightlim_D=compressl(o_avej_d[d=2,($oper)]*3600*24/femask_abio/0.5/(FeMask_bct-2.6)/dayfracmean/FeMask_felim_D)*1000
let FeDyn1_lightlim=compressl(o_avej[d=3,($oper)]*3600*24/fedyn1_abio/FeDyn1_bct/dayfracmean/FeDyn1_felim)*1000
let FeDyn1_lightlim_D=compressl(o_avej_d[d=3,($oper)]*3600*24/fedyn1_abio/0.5/(FeDyn1_bct-2.6)/dayfracmean/FeDyn1_felim_D)*1000



!!!!!!!!!!!!!!!!!!!!!!!Fe
! noFe has no felim

def sym oper3=k=1,l=@ave,gxyz=dummy
let FeMask_felim1= if (o_dissolved_fe[d=4,($oper3)] gt 0) then o_dissolved_fe[d=4,($oper3)]*1e6/(o_dissolved_fe[d=4,($oper3)]*1e6+FeMask_kfe) else dummy
let FeMask_felim=compressl(FeMask_felim1)
let FeMask_felim_D1= if (o_dissolved_fe[d=4,($oper3)] gt 0) then o_dissolved_fe[d=4,($oper3)]*1e6/(o_dissolved_fe[d=4,($oper3)]*1e6+FeMask_kfe_D) else dummy
let FeMask_felim_D=compressl(FeMask_felim_D1)

let FeDyn1_felim= o_deffe[d=3,($oper)]
let FeDyn1_felim_D1=if (o_dfe[d=3,($oper)] gt 0) then o_dfe[d=3,($oper)]*1e3/(o_dfe[d=3,($oper)]*1e3+FeDyn1_kfe_D) else dummy
let FeDyn1_felim_D=compressl(FeDyn1_felim_D1)



!!!!!!!!!!!!!!!!!!!!!!! NO3
let NoFe_nlim=if (o_no3[d=1,($oper)] gt 0) then o_no3[d=1,($oper)]*1e3/(o_no3[d=1,($oper)]*1e3+k1n) else dummy
let FeMask_nlim=if (o_no3[d=2,($oper)] gt 0) then o_no3[d=2,($oper)]*1e3/(o_no3[d=2,($oper)]*1e3+k1n) else dummy
let FeDyn1_nlim=if (o_no3[d=3,($oper)] gt 0) then o_no3[d=3,($oper)]*1e3/(o_no3[d=3,($oper)]*1e3+k1n) else dummy


!!!!!!!!!!!!!!!!!!!!!!!!PO4
let NoFe_plim=if (o_po4[d=1,($oper)] gt 0) then o_po4[d=1,($oper)]*1e3/(o_po4[d=1,($oper)]*1e3+k1p) else dummy
let FeMask_plim=if (o_po4[d=2,($oper)] gt 0) then o_po4[d=2,($oper)]*1e3/(o_po4[d=2,($oper)]*1e3+k1p) else dummy
let FeDyn1_plim=if (o_po4[d=3,($oper)] gt 0) then o_po4[d=3,($oper)]*1e3/(o_po4[d=3,($oper)]*1e3+k1p) else dummy


let NoFe_plim_D=if (o_po4[d=1,($oper)] gt 0) then o_po4[d=1,($oper)]*1e3/(o_po4[d=1,($oper)]*1e3+k1p) else dummy
let FeMask_plim_D=if (o_po4[d=2,($oper)] gt 0) then o_po4[d=2,($oper)]*1e3/(o_po4[d=2,($oper)]*1e3+k1p) else dummy
let FeDyn1_plim_D=if (o_po4[d=3,($oper)] gt 0) then o_po4[d=3,($oper)]*1e3/(o_po4[d=3,($oper)]*1e3+k1p) else dummy


! grazing control
let NoFe_grzlim=o_phytgraz[d=1,($oper)]/o_phytnpp[d=1,($oper)]

! make mask for phyt
!!!! N,P,Fe compare
! noFe
let nofe_phynpfe_n=NoFe_nlim
let nofe_phynpfe_p=NoFe_plim
let nofe_phynpfe_l=NoFe_lightlim
let nofe_phynpfe_t=NoFe_templim
let mask_NoFe_phy_n= if (nofe_phynpfe_n lt nofe_phynpfe_p) and (nofe_phynpfe_n lt nofe_phynpfe_l) and (nofe_phynpfe_n lt nofe_phynpfe_t) then dummy+1 else dummy
let mask_Nofe_phy_p=if (nofe_phynpfe_p lt nofe_phynpfe_n) and (nofe_phynpfe_p lt nofe_phynpfe_l) and (nofe_phynpfe_p lt nofe_phynpfe_t) then dummy+1 else dummy
let mask_Nofe_phy_l=if (nofe_phynpfe_l lt nofe_phynpfe_n) and (nofe_phynpfe_l lt nofe_phynpfe_p) and (nofe_phynpfe_l lt nofe_phynpfe_t) then dummy+1 else dummy
let mask_Nofe_phy_t=if (nofe_phynpfe_t lt nofe_phynpfe_n) and (nofe_phynpfe_t lt nofe_phynpfe_p) and (nofe_phynpfe_t lt nofe_phynpfe_l) then dummy+1 else dummy


let mask_NoFe=mask_NoFe_phy_n+mask_Nofe_phy_p*2+mask_Nofe_phy_l*3+mask_Nofe_phy_t*5



! FeMask

let FeMask_phynpfe_n=FeMask_nlim
let FeMask_phynpfe_p=FeMask_plim
let FeMask_phynpfe_l=femask_lightlim
let FeMask_phynpfe_fe=FeMask_felim
let FeMask_phynpfe_t=FeMask_templim
let FeMask_n_npfe= if (FeMask_phynpfe_n lt FeMask_phynpfe_P) and (FeMask_phynpfe_n lt FeMask_phynpfe_fe) and (FeMask_phynpfe_n lt FeMask_phynpfe_l) and (FeMask_phynpfe_n lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_p_npfe= if (FeMask_phynpfe_p lt FeMask_phynpfe_n) and (FeMask_phynpfe_p lt FeMask_phynpfe_fe) and (FeMask_phynpfe_p lt FeMask_phynpfe_l) and (FeMask_phynpfe_p lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_l_npfe= if (FeMask_phynpfe_l lt FeMask_phynpfe_n) and (FeMask_phynpfe_l lt FeMask_phynpfe_fe) and (FeMask_phynpfe_l lt FeMask_phynpfe_p) and (FeMask_phynpfe_l lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_fe_npfe= if (FeMask_phynpfe_fe lt FeMask_phynpfe_n) and (FeMask_phynpfe_fe lt FeMask_phynpfe_p) and (FeMask_phynpfe_fe lt FeMask_phynpfe_l) and (FeMask_phynpfe_fe lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_t_npfe= if (FeMask_phynpfe_t lt FeMask_phynpfe_n) and (FeMask_phynpfe_t lt FeMask_phynpfe_p) and (FeMask_phynpfe_t lt FeMask_phynpfe_l) and (FeMask_phynpfe_t lt FeMask_phynpfe_fe) then dummy +1 else dummy
let mask_FeMask=FeMask_n_npfe+FeMask_p_npfe*2+FeMask_l_npfe*3+FeMask_fe_npfe*4+FeMask_t_npfe*5



! FeDyn1
let FeDyn1_phynpfe_n=FeDyn1_nlim
let FeDyn1_phynpfe_p=FeDyn1_plim
let FeDyn1_phynpfe_l=FeDyn1_lightlim
let FeDyn1_phynpfe_fe=compressl(FeDyn1_felim)
let FeDyn1_phynpfe_t=compressl(FeDyn1_templim)
let FeDyn1_n_npfe= if (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_P) and (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_fe) and (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_p_npfe= if (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_fe) and (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_l_npfe= if (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_fe) and (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_p) and (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_fe_npfe= if (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_p) and (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_t_npfe= if (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_p) and (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_fe) then dummy +1 else dummy
let mask_FeDyn1=FeDyn1_n_npfe+FeDyn1_p_npfe*2+FeDyn1_l_npfe*3+FeDyn1_fe_npfe*4+FeDyn1_t_npfe*5







go ilandscape4x2_ext.jnl
!!! nlim=1, plim=2, felim =3
go preparefigure.jnl FigA1
go transparency
set win 1
set vie v1
sha/nolab/pal=rainbow/set/lev="(0.5,5.5,1)"/nokey mask_NoFe[X=30E:30E]
ppl ylab "1800"
ppl xlab "NoFe"
ppl axlabp,1,-1
ppl xaxis, 30,390,60 
ppl shade
go land
label -30,110,-1,0.2 (a)
ppl shakey ,,0
set vie v2
sha/nolab/set/lev="(0.5,5.5,1)"/nokey mask_FeMask[X=30E:30E]
ppl xlab "FeMask"
ppl axlabp,1,0
ppl xaxis, 30,390,60 
ppl shade
go land
label -30,110,-1,0.2 (b)
set vie v3
sha/nolab/set/lev="(0.5,5.5,1)" mask_FeDyn1[X=30E:30E]
ppl xlab "FeDyn"
ppl axlabp,1,0
ppl xaxis, 30,390,60 
ppl color 6 100,100,100
ppl shakey 1,1,0.000001,0,0,0,4.6,4.75,1.5,4.8 @P6
ppl shade
go land
label -30,110,-1,0.2 (c)
label/nouser 3.8, 0.4, -1,0,0.15 "NO_3"
label/nouser 3.8, 1.05, -1,0,0.15 "PO_4"
label/nouser 3.8, 1.70, -1,0,0.15 "light"
label/nouser 3.8, 2.35, -1,0,0.15 "Fe"
label/nouser 3.8, 3, -1,0,0.15 "Temp"

!go finalizefigure_psthicken

! mask mask for diaz

!!!! N,P,Fe compare
! noFe

let nofe_phynpfe_p_D=NoFe_plim_D
let nofe_phynpfe_l_D=NoFe_lightlim_D
let nofe_phynpfe_t_D=NoFe_templim_D
let mask_NoFe_p_D= if (nofe_phynpfe_p_D lt nofe_phynpfe_l_D) and (nofe_phynpfe_p_D lt nofe_phynpfe_t_D) then dummy+1 else dummy
let mask_NoFe_l_D= if (nofe_phynpfe_l_D lt nofe_phynpfe_p_D) and (nofe_phynpfe_l_D lt nofe_phynpfe_t_D) then dummy+1 else dummy
let mask_NoFe_t_D= if (nofe_phynpfe_t_D lt nofe_phynpfe_l_D) and (nofe_phynpfe_t_D lt nofe_phynpfe_p_D) then dummy+1 else dummy
let mask_NoFe_D= mask_NoFe_p_D*2+mask_NoFe_l_D*3+mask_NoFe_t_D*5

! FeMask

let FeMask_phynpfe_p_D=FeMask_plim_D
let FeMask_phynpfe_l_D=FeMask_lightlim_D
let FeMask_phynpfe_fe_D=FeMask_felim_D
let FeMask_phynpfe_t_D=FeMask_templim_D
let mask_FeMask_p_D= if (FeMask_phynpfe_p_D lt FeMask_phynpfe_l_D) and (FeMask_phynpfe_p_D lt FeMask_phynpfe_t_D) and (FeMask_phynpfe_p_D lt FeMask_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeMask_l_D= if (FeMask_phynpfe_l_D lt FeMask_phynpfe_p_D) and (FeMask_phynpfe_l_D lt FeMask_phynpfe_t_D) and (FeMask_phynpfe_l_D lt FeMask_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeMask_t_D= if (FeMask_phynpfe_t_D lt FeMask_phynpfe_l_D) and (FeMask_phynpfe_t_D lt FeMask_phynpfe_p_D) and (FeMask_phynpfe_t_D lt FeMask_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeMask_fe_D= if (FeMask_phynpfe_fe_D lt FeMask_phynpfe_l_D) and (FeMask_phynpfe_fe_D lt FeMask_phynpfe_p_D) and (FeMask_phynpfe_fe_D lt FeMask_phynpfe_t_D) then dummy+1 else dummy
let mask_FeMask_D= mask_FeMask_p_D*2+mask_FeMask_l_D*3+mask_FeMask_fe_D*4+mask_FeMask_t_D*5

! FeDyn1
let FeDyn1_phynpfe_p_D=FeDyn1_plim_D
let FeDyn1_phynpfe_l_D=FeDyn1_lightlim_D
let FeDyn1_phynpfe_fe_D=FeDyn1_felim_D
let FeDyn1_phynpfe_t_D=FeDyn1_templim_D
let mask_FeDyn1_p_D= if (FeDyn1_phynpfe_p_D lt FeDyn1_phynpfe_l_D) and (FeDyn1_phynpfe_p_D lt FeDyn1_phynpfe_t_D) and (FeDyn1_phynpfe_p_D lt FeDyn1_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeDyn1_l_D= if (FeDyn1_phynpfe_l_D lt FeDyn1_phynpfe_p_D) and (FeDyn1_phynpfe_l_D lt FeDyn1_phynpfe_t_D) and (FeDyn1_phynpfe_l_D lt FeDyn1_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeDyn1_t_D= if (FeDyn1_phynpfe_t_D lt FeDyn1_phynpfe_l_D) and (FeDyn1_phynpfe_t_D lt FeDyn1_phynpfe_p_D) and (FeDyn1_phynpfe_t_D lt FeDyn1_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeDyn1_fe_D= if (FeDyn1_phynpfe_fe_D lt FeDyn1_phynpfe_l_D) and (FeDyn1_phynpfe_fe_D lt FeDyn1_phynpfe_p_D) and (FeDyn1_phynpfe_fe_D lt FeDyn1_phynpfe_t_D) then dummy+1 else dummy
let mask_FeDyn1_D= mask_FeDyn1_p_D*2+mask_FeDyn1_l_D*3+mask_FeDyn1_fe_D*4+mask_FeDyn1_t_D*5


! FeDyn0
let FeDyn0_phynpfe_p_D=FeDyn0_plim_D
let FeDyn0_phynpfe_l_D=FeDyn0_lightlim_D
let FeDyn0_phynpfe_fe_D=FeDyn0_felim_D
let FeDyn0_phynpfe_t_D=FeDyn0_templim_D
let mask_FeDyn0_p_D= if (FeDyn0_phynpfe_p_D lt FeDyn0_phynpfe_l_D) and (FeDyn0_phynpfe_p_D lt FeDyn0_phynpfe_t_D) and (FeDyn0_phynpfe_p_D lt FeDyn0_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeDyn0_l_D= if (FeDyn0_phynpfe_l_D lt FeDyn0_phynpfe_p_D) and (FeDyn0_phynpfe_l_D lt FeDyn0_phynpfe_t_D) and (FeDyn0_phynpfe_l_D lt FeDyn0_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeDyn0_t_D= if (FeDyn0_phynpfe_t_D lt FeDyn0_phynpfe_l_D) and (FeDyn0_phynpfe_t_D lt FeDyn0_phynpfe_p_D) and (FeDyn0_phynpfe_t_D lt FeDyn0_phynpfe_fe_D) then dummy+1 else dummy
let mask_FeDyn0_fe_D= if (FeDyn0_phynpfe_fe_D lt FeDyn0_phynpfe_l_D) and (FeDyn0_phynpfe_fe_D lt FeDyn0_phynpfe_p_D) and (FeDyn0_phynpfe_fe_D lt FeDyn0_phynpfe_t_D) then dummy+1 else dummy
let mask_FeDyn0_D= mask_FeDyn0_p_D*2+mask_FeDyn0_l_D*3+mask_FeDyn0_fe_D*4+mask_FeDyn0_t_D*5





!!!!!!!!!!!!!!!!!!!!!!!!! 300st year

! calculation
!!!!!!!!!!!!!!!!!!!!!temp
let NoFe_bct=compressl(bbio^(cbio*o_temp[d=1,($oper2)]))
let FeMask_bct=compressl(bbio^(cbio*o_temp[d=2,($oper2)]))
let FeDyn1_bct=compressl(bbio^(cbio*o_temp[d=3,($oper2)]))

let NoFe_templim=NoFe_bct/bctNormal
let NoFe_templim_D=NoFe_bct/bctNormal_D
let FeMask_templim=FeMask_bct/bctNormal
let FeMask_templim_D=FeMask_bct/bctNormal_D
let FeDyn1_templim=FeDyn1_bct/bctNormal
let FeDyn1_templim_D=FeDyn1_bct/bctNormal_D

!!!!!!!!!!!!!!!!!!!!!!!light
let noFe_lightlim=compressl(o_avej[d=1,($oper2)]*3600*24/nofe_abio/NoFe_bct/dayfracmean)*1000
let noFe_lightlim_D=compressl(o_avej_d[d=1,($oper2)]*3600*24/nofe_abio/0.5/(NoFe_bct-2.6)/dayfracmean)*1000
let FeMask_lightlim=compressl(o_avej[d=2,($oper2)]*3600*24/femask_abio/FeMask_bct/dayfracmean/FeMask_felim)*1000 !/o_avej[d=2,l=@ave,x=@max,y=@max,k=1])!
let FeMask_lightlim_D=compressl(o_avej_d[d=2,($oper2)]*3600*24/femask_abio/0.5/(FeMask_bct-2.6)/dayfracmean/FeMask_felim_D)*1000
let FeDyn1_lightlim=compressl(o_avej[d=3,($oper2)]*3600*24/fedyn1_abio/FeDyn1_bct/dayfracmean/FeDyn1_felim)*1000
let FeDyn1_lightlim_D=compressl(o_avej_d[d=3,($oper2)]*3600*24/fedyn1_abio/0.5/(FeDyn1_bct-2.6)/dayfracmean/FeDyn1_felim_D)*1000


!!!!!!!!!!!!!!!!!!!!!!!Fe
! noFe has no felim


let FeMask_felim1= if (o_dissolved_fe[d=4,($oper3)] gt 0) then o_dissolved_fe[d=4,($oper3)]*1e6/(o_dissolved_fe[d=4,($oper3)]*1e6+FeMask_kfe) else dummy
let FeMask_felim=compressl(FeMask_felim1)
let FeMask_felim_D1= if (o_dissolved_fe[d=4,($oper3)] gt 0) then o_dissolved_fe[d=4,($oper3)]*1e6/(o_dissolved_fe[d=4,($oper3)]*1e6+FeMask_kfe_D) else dummy
let FeMask_felim_D=compressl(FeMask_felim_D1)

let FeDyn1_felim= o_deffe[d=3,($oper2)]
let FeDyn1_felim_D1=if (o_dfe[d=3,($oper2)] gt 0) then o_dfe[d=3,($oper2)]*1e3/(o_dfe[d=3,($oper2)]*1e3+FeDyn1_kfe_D) else dummy
let FeDyn1_felim_D=compressl(FeDyn1_felim_D1)




!!!!!!!!!!!!!!!!!!!!!!! NO3
let NoFe_nlim=if (o_no3[d=1,($oper2)] gt 0) then o_no3[d=1,($oper2)]*1e3/(o_no3[d=1,($oper2)]*1e3+k1n) else dummy
let FeMask_nlim=if (o_no3[d=2,($oper2)] gt 0) then o_no3[d=2,($oper2)]*1e3/(o_no3[d=2,($oper2)]*1e3+k1n) else dummy
let FeDyn1_nlim=if (o_no3[d=3,($oper2)] gt 0) then o_no3[d=3,($oper2)]*1e3/(o_no3[d=3,($oper2)]*1e3+k1n) else dummy

!!!!!!!!!!!!!!!!!!!!!!!!PO4
let NoFe_plim=if (o_po4[d=1,($oper2)] gt 0) then o_po4[d=1,($oper2)]*1e3/(o_po4[d=1,($oper2)]*1e3+k1p) else dummy
let FeMask_plim=if (o_po4[d=2,($oper2)] gt 0) then o_po4[d=2,($oper2)]*1e3/(o_po4[d=2,($oper2)]*1e3+k1p) else dummy
let FeDyn1_plim=if (o_po4[d=3,($oper2)] gt 0) then o_po4[d=3,($oper2)]*1e3/(o_po4[d=3,($oper2)]*1e3+k1p) else dummy

let NoFe_plim_D=if (o_po4[d=1,($oper2)] gt 0) then o_po4[d=1,($oper2)]*1e3/(o_po4[d=1,($oper2)]*1e3+k1p) else dummy
let FeMask_plim_D=if (o_po4[d=2,($oper2)] gt 0) then o_po4[d=2,($oper2)]*1e3/(o_po4[d=2,($oper2)]*1e3+k1p) else dummy
let FeDyn1_plim_D=if (o_po4[d=3,($oper2)] gt 0) then o_po4[d=3,($oper2)]*1e3/(o_po4[d=3,($oper2)]*1e3+k1p) else dummy



! make mask for phyt
!!!! N,P,Fe compare
! noFe
let nofe_phynpfe_n=NoFe_nlim
let nofe_phynpfe_p=NoFe_plim
let nofe_phynpfe_l=NoFe_lightlim
let nofe_phynpfe_t=NoFe_templim
let mask_NoFe_phy_n= if (nofe_phynpfe_n lt nofe_phynpfe_p) and (nofe_phynpfe_n lt nofe_phynpfe_l) and (nofe_phynpfe_n lt nofe_phynpfe_t) then dummy+1 else dummy
let mask_Nofe_phy_p=if (nofe_phynpfe_p lt nofe_phynpfe_n) and (nofe_phynpfe_p lt nofe_phynpfe_l) and (nofe_phynpfe_p lt nofe_phynpfe_t) then dummy+1 else dummy
let mask_Nofe_phy_l=if (nofe_phynpfe_l lt nofe_phynpfe_n) and (nofe_phynpfe_l lt nofe_phynpfe_p) and (nofe_phynpfe_l lt nofe_phynpfe_t) then dummy+1 else dummy
let mask_Nofe_phy_t=if (nofe_phynpfe_t lt nofe_phynpfe_n) and (nofe_phynpfe_t lt nofe_phynpfe_p) and (nofe_phynpfe_t lt nofe_phynpfe_l) then dummy+1 else dummy
let mask_NoFe=mask_NoFe_phy_n+mask_Nofe_phy_p*2+mask_Nofe_phy_l*3+mask_Nofe_phy_t*5

! FeMask
let FeMask_phynpfe_n=FeMask_nlim
let FeMask_phynpfe_p=FeMask_plim
let FeMask_phynpfe_l=femask_lightlim
let FeMask_phynpfe_fe=FeMask_felim
let FeMask_phynpfe_t=FeMask_templim
let FeMask_n_npfe= if (FeMask_phynpfe_n lt FeMask_phynpfe_P) and (FeMask_phynpfe_n lt FeMask_phynpfe_fe) and (FeMask_phynpfe_n lt FeMask_phynpfe_l) and (FeMask_phynpfe_n lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_p_npfe= if (FeMask_phynpfe_p lt FeMask_phynpfe_n) and (FeMask_phynpfe_p lt FeMask_phynpfe_fe) and (FeMask_phynpfe_p lt FeMask_phynpfe_l) and (FeMask_phynpfe_p lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_l_npfe= if (FeMask_phynpfe_l lt FeMask_phynpfe_n) and (FeMask_phynpfe_l lt FeMask_phynpfe_fe) and (FeMask_phynpfe_l lt FeMask_phynpfe_p) and (FeMask_phynpfe_l lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_fe_npfe= if (FeMask_phynpfe_fe lt FeMask_phynpfe_n) and (FeMask_phynpfe_fe lt FeMask_phynpfe_p) and (FeMask_phynpfe_fe lt FeMask_phynpfe_l) and (FeMask_phynpfe_fe lt FeMask_phynpfe_t) then dummy +1 else dummy
let FeMask_t_npfe= if (FeMask_phynpfe_t lt FeMask_phynpfe_n) and (FeMask_phynpfe_t lt FeMask_phynpfe_p) and (FeMask_phynpfe_t lt FeMask_phynpfe_l) and (FeMask_phynpfe_t lt FeMask_phynpfe_fe) then dummy +1 else dummy
let mask_FeMask=FeMask_n_npfe+FeMask_p_npfe*2+FeMask_l_npfe*3+FeMask_fe_npfe*4+FeMask_t_npfe*5

! FeDyn1
let FeDyn1_phynpfe_n=FeDyn1_nlim
let FeDyn1_phynpfe_p=FeDyn1_plim
let FeDyn1_phynpfe_l=FeDyn1_lightlim
let FeDyn1_phynpfe_fe=compressl(FeDyn1_felim)
let FeDyn1_phynpfe_t=compressl(FeDyn1_templim)
let FeDyn1_n_npfe= if (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_P) and (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_fe) and (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_n lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_p_npfe= if (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_fe) and (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_p lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_l_npfe= if (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_fe) and (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_p) and (FeDyn1_phynpfe_l lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_fe_npfe= if (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_p) and (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_fe lt FeDyn1_phynpfe_t) then dummy +1 else dummy
let FeDyn1_t_npfe= if (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_n) and (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_p) and (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_l) and (FeDyn1_phynpfe_t lt FeDyn1_phynpfe_fe) then dummy +1 else dummy
let mask_FeDyn1=FeDyn1_n_npfe+FeDyn1_p_npfe*2+FeDyn1_l_npfe*3+FeDyn1_fe_npfe*4+FeDyn1_t_npfe*5




set vie v5
sha/nolab/set/lev="(0.5,5.5,1)"/nokey mask_NoFe[X=30E:30E]
ppl ylab "2100"
ppl xlab "NoFe"
ppl axlabp,-1,-1
ppl xaxis, 30,390,60 
ppl shade
go land
label -30,110,-1,0.2 (d)

set vie v6
sha/nolab/set/lev="(0.5,5.5,1)"/nokey mask_FeMask[X=30E:30E]
ppl axlabp,-1,0
ppl xlab "FeMask"
ppl xaxis, 30,390,60 
ppl shade
go land
label -30,110,-1,0.2 (e)

set vie v7
sha/nolab/set/lev="(0.5,5.5,1)"/nokey mask_FeDyn1[X=30E:30E]
ppl xlab "FeDyn"
ppl axlabp,-1,0
ppl xaxis, 30,390,60 
ppl color 6 100,100,100
ppl shakey 1,1,0.000001,0,0,0,4.6,4.75,1.5,4.8 @P6
ppl shade
go land
label/nouser 3.8, 0.4, -1,0,0.15 "NO_3"
label/nouser 3.8, 1.05, -1,0,0.15 "PO_4"
label/nouser 3.8, 1.70, -1,0,0.15 "light"
label/nouser 3.8, 2.35, -1,0,0.15 "Fe"
label/nouser 3.8, 3, -1,0,0.15 "Temp"
label -30,110,-1,0.2 (f)

go finalizefigure_psthicken