######################################################
#### 
#### R-Code for Indices of Latitudinal Position
#### of tropical edge, S. Broennimann, 23 Aug 2015
#### 
######################################################
#### 
#### INPUT
#### 
#### 1. ascii file of zonally averaged quantity
####    File has tabular format 
####    rows = timesteps [months of full calender years], columns = latitudes 
####    Headerline [cols 3 onward: equally spaced, increasing lat]
####    Col 1 = year, Col 2 = month [1..12], Col 3 onward = values 
####    see sample file "era20c_slp_zonmean.txt"
#### 
#### 2. Requires ranges (latitudes) and definition of seasons given in Table S2
#### 
#### OUTPUT: ascii file with: Year Month Latitude 
#### 
######################################################
#### Modify input below this line 

varName <- "slp"           # used for input/output file name 
dataset <- "era20c"        # used for input/output file name 
indexType <- "edge"        # "itcz", "bulk" or "edge"          
searchfor <- "max"         # "max" or "min" 
minLatWinter <- 24         # for possible lats/spacing see lat/res of input file below
maxLatWinter <- 38
minLatSummer <- 26
maxLatSummer <- 46
summer <- "JJASON"

######################################################
#### Reads file, define latitudes, timesteps, range, season 
######################################################
#setwd("..") 

infile <- paste(dataset,"_",varName,"_zonmean.txt",sep="")
lat <- as.numeric(read.table(infile,header=F,nrows=1)[-c(1,2)]) # latitudes
res <- lat[2]-lat[1]  										  # spacing of latitudes
inp <- as.matrix(read.table(infile,header=T))
noTimesteps <- dim(inp)[1]
noLat <- dim(inp)[2]
ranges <- c(maxLatWinter,minLatWinter,maxLatSummer,minLatSummer)
startMonIndex <- gregexpr(summer,"JFMAMJJASOND")[[1]][1] 
	if(startMonIndex == -1)
  		stop("Error: Defined season is not element of JFMAMJJASOND")
noMonths <- nchar(summer)
summerMonths <- c(rep(0,startMonIndex-1),rep(1,noMonths),rep(0,12-(startMonIndex-1)-noMonths))
m <- rep(summerMonths,noTimesteps/12)
output <- array(NA,dim=c(noTimesteps,3))
output[,1:2] <- inp[,1:2] 

######################################################
#### Main Loop 
######################################################
for (i in 1:noTimesteps){
#### sets the latitude range to search for minima
if (m[i]==0){rangemax <- ranges[1];rangemin <- ranges[2]} else {rangemax <- ranges[3];rangemin <- ranges[4]}
#### converts latitude range to indices of vector lat
idmax <- max(c(1:noLat)[abs(lat-rangemax)==min(abs(lat-rangemax))])
idmin <- min(c(1:noLat)[abs(lat-rangemin)==min(abs(lat-rangemin))])
#### code searches minimum; if maxima are searched for the variable is inverted
if (searchfor=="max"){y <- -inp[i,-c(1:2)]} else {y <- inp[i,-c(1:2)]}
#### searches localminima, ties that are minima, and minima at edge 
localmin <- c(F,y[(idmin+1):(idmax-1)]<y[idmin:(idmax-2)]&y[(idmin+1):(idmax-1)]<y[(idmin+2):idmax],F)
ties <- c(F,(y[(idmin+1):(idmax-2)]<y[idmin:(idmax-3)])&(y[(idmin+1):(idmax-2)]==y[(idmin+2):(idmax-1)])&(y[(idmin+1):(idmax-2)]<y[(idmin+3):(idmax)]),F,F)
southedgemin <- if (y[idmin]==min(y[idmin:idmax])){c(T,rep(F,idmax-idmin))}else{rep(F,(idmax+1)-idmin)}
northedgemin <- if (y[idmax]==min(y[idmin:idmax])){c(rep(F,idmax-idmin),T)}else{rep(F,(idmax+1)-idmin)}
#### decides which of the minima is selected (different for itcz and edge) 
if (indexType=="itcz"){if (sum(localmin|ties)==0){minimum <- (idmin:idmax)[southedgemin|northedgemin]} else 
 {minimum <- min((idmin:idmax)[y[idmin:idmax]==min(y[idmin:idmax][localmin|ties])])}} else {minimum <-  min((idmin:idmax)[southedgemin|localmin|ties|northedgemin])}
#### calculate the latitude of the chosen minimum:  
#### if minimum is a tie, chose average, if its is at the edge, chose edge, else interpolate
#### only double ties are considered 
if ((minimum<idmax)&(y[minimum]==y[minimum+1])){output[i,3] <- lat[minimum]+0.5*res} else {
if (minimum==idmin|minimum==idmax){output[i,3] <- lat[minimum]} else {
output[i,3] <- lat[minimum]-((res*(y[minimum+1]-y[minimum-1]))/(2*abs(y[minimum]-max(y[minimum+1],y[minimum-1]))))}}}

#########################################################
#### write output
#########################################################

write.table(output,paste(indexType,"_",varName,"_",searchfor,"_",dataset,"_",summer,"_win_",minLatWinter,"_",maxLatWinter,"_sum_",minLatSummer,"_",maxLatSummer,".txt",sep=""))