set memory/size=500
use by4grid  ! original kmt

let ulon=tlong
let ulat=tlat

! Torres Strait 1

let xto1 = if ulon eq ulon[i=163] then 1
let yto1 = if ulat ge ulat[j=148] then 1
let yto2 = if ulat le ulat[j=150] then 1
let mto = missing(xto1*yto1*yto2,0)

! Torres Strait 2

let xto21 = if ulon eq ulon[i=164] then 1
let yto21 = if ulat ge ulat[j=147] then 1
let yto22 = if ulat le ulat[j=151] then 1
let mto2 = missing(xto21*yto21*yto22,0)

! Olom 

let xol1 = if ulon eq ulon[i=131] then 1
let yol1 = if ulat ge ulat[j=165] then 1
let yol2 = if ulat le ulat[j=166] then 1
let mol = missing(xol1*yol1*yol2,0)

! Halmahera1

let xha11 = if ulon eq ulon[i=150] then 1
let yha11 = if ulat ge ulat[j=187] then 1
let yha12 = if ulat le ulat[j=188] then 1
let mha1 = missing(xha11*yha11*yha12,0)

! Halmahera2

let xha21 = if ulon eq ulon[i=152] then 1
let yha21 = if ulat ge ulat[j=185] then 1
let yha22 = if ulat le ulat[j=187] then 1
let mha2 = missing(xha21*yha21*yha22,0)

! Halmahera3

let xha31 = if ulon eq ulon[i=153] then 1
let yha31 = if ulat ge ulat[j=185] then 1
let yha32 = if ulat le ulat[j=185] then 1
let mha3 = missing(xha31*yha31*yha32,0)

! Luzon1

let xlu1 = if ulon eq ulon[i=144] then 1
let ylu1 = if ulat eq ulat[j=251] then 1
let mlu1 = missing(xlu1*ylu1,0)

! Luzon2

let xlu2 = if ulon eq ulon[i=144] then 1
let ylu2 = if ulat eq ulat[j=260] then 1
let mlu2 = missing(xlu2*ylu2,0)

! Lombok1

let xlo11 = if ulon ge ulon[i=146] then 1
let xlo12 = if ulon le ulon[i=147] then 1
let ylo11 = if ulat eq ulat[j=155] then 1
let mlo1  = missing(xlo11*xlo12*ylo11,0)

! Lombok2

let xlo22 = if ulon eq ulon[i=147] then 1
let ylo21 = if ulat eq ulat[j=154] then 1
let mlo2  = missing(xlo22*ylo21,0)

! Lombok3

let xlo31 = if ulon ge ulon[i=140] then 1
let xlo32 = if ulon le ulon[i=142] then 1
let ylo31 = if ulat eq ulat[j=155] then 1
let mlo3  = missing(xlo31*xlo32*ylo31,0)

! Borneo0

let xbo01 = if ulon eq ulon[i=145] then 1
let ybo02 = if ulat eq ulat[j=237] then 1
let mbo0  = missing(xbo01*ybo02,0)

! Borneo1

let xbo11 = if ulon eq ulon[i=144] then 1
let ybo11 = if ulat ge ulat[j=236] then 1
let ybo12 = if ulat le ulat[j=237] then 1
let mbo1  = missing(xbo11*ybo11*ybo12,0)

! Borneo2

let xbo21 = if ulon eq ulon[i=144] then 1
let ybo21 = if ulat ge ulat[j=232] then 1
let ybo22 = if ulat le ulat[j=233] then 1
let mbo2  = missing(xbo21*ybo21*ybo22,0)

! Borneo3

let xbo31 = if ulon eq ulon[i=141] then 1
let ybo31 = if ulat ge ulat[j=209] then 1
let ybo32 = if ulat le ulat[j=216] then 1
let mbo3  = missing(xbo31*ybo31*ybo32,0)

! Borneo4

let xbo41 = if ulon eq ulon[i=142] then 1
let ybo41 = if ulat ge ulat[j=216] then 1
let ybo42 = if ulat le ulat[j=222] then 1
let mbo4  = missing(xbo41*ybo41*ybo42,0)

! Borneo5

let xbo51 = if ulon eq ulon[i=143] then 1
let ybo51 = if ulat ge ulat[j=222] then 1
let ybo52 = if ulat le ulat[j=228] then 1
let mbo5  = missing(xbo51*ybo51*ybo52,0)



let maska=mto+mto2+mol+mha1+mha2+mha3+mlu1+mlu2+mlo1+mlo2+mlo3
let maskb=mbo0+mbo1+mbo2+mbo3+mbo4+mbo5
let mask=maska+maskb


! make new kmt
let km2=3-kmt                        ! difference of kmt to surf.
let kma=km2*mask                      
let nkm0=kmt+kma                      ! kmt identical + k=1 at insets

! raising sill depth in Luzon strait from 2500m(k=29) to 1500m(k=24)

let sil1 = if ulon eq ulon[i=144] then 1
let sil2 = if ulat ge ulat[j=252] then 1
let sil3 = if ulat le ulat[j=259] then 1

let silmask=missing(sil1*sil2*sil3,0) * 5   ! 5 point shallower
let nkm = nkm0 - silmask
! make new elevation
let nele=nkm*100

save/file=ngrid.nc nkm
save/file=ngrid.nc/append nele