Create_BGCArgo_data.R

#install.packages("raadfiles",force = T)
#install.packages("raadtools",force = T)
#install.packages(c("oce","seacarb","birk","devtools","segmented","caTools","R.utils))
library(ncdf4)
library(data.table)
library(dplyr)
library(birk)
library(raadfiles)
library(raadtools)
library(lubridate)
#devtools::install_github("KimBaldry/BIOMATE-Rpackage",force = T)
##### Draw from the index file and grey list to subset profiles
## download greylist
url = "ftp://usgodae.org/pub/outgoing/argo/" 
grey = fread(paste(url,"ar_greylist.txt",sep = ""),sep = ",",header =T,sep2 = F)
# load("grey.RData")
# access index file and filter
meta  <- data.table::fread("/rdsi/PUBLIC/raad/data/www.usgodae.org/ftp/outgoing/argo/argo_synthetic-profile_index.txt.gz", skip = 8, data.table = FALSE)
meta <- meta %>%  dplyr::mutate(date2 = substr(date,1,8))
meta$year = as.numeric(substr(meta[,2],1,4))
meta$month = as.numeric(substr(meta[,2],5,6))
meta$day = as.numeric(substr(meta[,2],7,8))
meta <- meta %>% dplyr::filter(latitude <= -30, grepl("CHLA", parameters), year != 2024)
# create some new variables

meta <- meta %>%  dplyr::mutate(PLATFORM_CODE = as.numeric(unlist(strsplit(meta$file,"/"))[c(F,T,F,F)]))


# compare index to grey list, and combine
files = meta
files = files %>% dplyr::mutate(PROFILE_ID = c(1:nrow(files)))
# we only care about some variables, some floats and some quality codes
grey_sub = grey %>% dplyr::filter(PARAMETER_NAME %in% c("CHLA","PSAL","TEMP","PRESS","BBP700"),PLATFORM_CODE %in% files$PLATFORM_CODE, QUALITY_CODE %in% c(3,4))
grey_sub$END_DATE[which(is.na(grey_sub$END_DATE)==T)] = gsub("-", "", Sys.Date())
# combine
joined = dplyr::left_join(files,grey_sub[,c(1:4)], by = "PLATFORM_CODE")
joined = joined %>% dplyr::mutate(on_greylist = as.POSIXct(joined$date2,format = "%Y%m%d", tz = "GMT") < as.POSIXct(joined$END_DATE,format = "%Y%m%d", tz = "GMT") & as.POSIXct(joined$date2,format = "%Y%m%d", tz = "GMT") > as.POSIXct(as.character(joined$START_DATE),format = "%Y%m%d", tz = "GMT") )
# filter based on greylist
good_profiles = joined %>% dplyr::filter(on_greylist == F | is.na(on_greylist))

# filter good profiles based on DAC
dfiles = good_profiles[!duplicated(good_profiles$PROFILE_ID),-c(18:21)] %>% filter(institution %in% c("AO","IF","CS"))
### Sesnor information
# create metadata file path to get sensor information
fl = sapply(strsplit(dfiles$file,"/"),"[[",1)
dfiles$meta_file = paste(fl,"/",dfiles$PLATFORM_CODE,"/",dfiles$PLATFORM_CODE,"_meta.nc",sep = "")
dfiles$sensor = NA
dfiles$sensor_maker = NA
for(fl in unique(dfiles$meta_file)){
  f2 = nc_open(paste("/rdsi/PUBLIC/raad/data/www.usgodae.org/ftp/outgoing/argo/dac/",fl,sep = ""))
  CHL_col = grep("FLUOROMETER_CHLA",ncvar_get(f2,"SENSOR"))
  if(length(CHL_col)){
    dfiles$sensor[which(dfiles$meta_file == fl)] = ncvar_get(f2,"SENSOR_MODEL")[CHL_col]
    dfiles$sensor_maker[which(dfiles$meta_file == fl)] = ncvar_get(f2,"SENSOR_MAKER")[CHL_col]}
  nc_close(f2)
}

# filter by sensor and manufacturer
dfiles = dfiles %>% filter(grepl("ECO_FLBB",sensor) | institution == "AO"| grepl("MCOMS_FLBB",sensor)  )
ids <- dfiles$file
A_files = argo_files()$fullname[basename(argo_files()$fullname) %in% basename(ids)]
A_files = A_files[grepl("usgodae",A_files)]
idx_list = 1:length(A_files)
dfiles = dfiles %>% filter(basename(file) %in% basename(A_files))
save(dfiles, joined,A_files,idx_list,file = "./Final_Data/BGCArgoIndex.RData")


# function to run in paralell to compile BGCArgo data
Argo_calcs = function(idx){
  tryCatch({
    # I couldnt get the functions to load without re-calling them in clusters
    library(ncdf4)
    library(data.table)
    library(dplyr)
    library(oce)
    library(seacarb)
    library(birk)
    library(raadfiles)
    library(raadtools)
    library(maptools)
    library(gsw)
    library(BIOMATE)
    #library(peaks)
    path = "/perm_storage/home/kbaldry/SRC/OSM"
    source(file.path(path,"front_class.R"))
    source(file.path("/perm_storage/home/kbaldry/Ocean_code/ocean_calcs","MLD_calcs.R"))
    source(file.path("/perm_storage/home/kbaldry/Ocean_code/cleaning/BSM","fit_bp_segmented.R"))
    source(file.path("/perm_storage/home/kbaldry/Ocean_code/cleaning/BSM","apply_bsm.R"))
    files.sources = list.files(c(file.path("/perm_storage/home/kbaldry/Ocean_code/cleaning"), file.path("/perm_storage/home/kbaldry/Ocean_code/NPQ_correction")),pattern = "*.R",full.names = T)
    sapply(files.sources, source)
    source(file.path("/perm_storage/home/kbaldry/SRC","abstract_profile.R"))
    
    # open the BGC-Argo profile
    f1 = nc_open(A_files[idx])
    Param = ncvar_get(f1, "STATION_PARAMETERS")
    BGC_col = col(Param)[grep("CHLA",Param)]
    if(length(BGC_col > 1)){BGC_col = BGC_col[1]}
    #### META DATA ####
    prof_data = data.frame("DIRECTION"= substr(ncvar_get(f1,"DIRECTION"),1,1))
    prof_data$idx = idx
    prof_data$float = dfiles$PLATFORM_CODE[idx]
    prof_data$inst = dfiles$institution[idx]
    # Date, location and time
    if(BGC_col == 1){ 
      prof_data$JULD = ncvar_get(f1, "JULD")
      prof_data$LAT = ncvar_get(f1,"LATITUDE")
      prof_data$LON = ncvar_get(f1,"LONGITUDE")}else{
        prof_data$JULD = ncvar_get(f1, "JULD")[1]
        prof_data$LAT = ncvar_get(f1,"LATITUDE")[1]
        prof_data$LON = ncvar_get(f1,"LONGITUDE")[1]
      }
    # convert Julien day to Rdate and year, month, day, hour
    DATE = lubridate::ymd_hms("1950-01-01 00:00:00") + prof_data$JULD*3600*24
    prof_data$YYYY = year(DATE)
    prof_data$MM = month(DATE)
    prof_data$DD = mday(DATE)
    prof_data$HH = hour(DATE)
    prof_data$DATE = DATE
    prof_data$YD = yday(DATE)
    prof_data$isNight = NA
    prof_data$QAnote = NA
    prof_data$HighSF_flag = F
    prof_data$HighSF_flag2 = F
    prof_data$HighSF_flag_changed = F
    prof_data$EMLD_flag_changed  = F
    # prof_data$isNight_5 = NA
    ## insert latitude/longitude interpolation here if needed
    # check if interpolated #
    # if(length(grep(8,ncvar_get(f1,POSITION_QC)))>0){f1.2D$Position_interpolated = "Y"}
    # else{f1.2D$Position_interpolated = "N"}
    
    # sun angle
    prof_data$SUN_ANGLE = sunAngle(DATE,longitude = prof_data$LON,latitude = prof_data$LAT)$altitude
    
    ### day night information
    coord <- SpatialPoints(data.frame(prof_data$LON,prof_data$LAT),proj4string=CRS("+proj=longlat +datum=WGS84"))
    rise = sunriset(coord,DATE, direction = "sunrise",POSIXct.out = T)$time
    set = sunriset(coord, DATE, direction = "sunset",POSIXct.out = T)$time
    
    rm(coord)
    if(!is.na(rise) & !is.na(set)){
    # if(DATE  > (set + 2*60*60) | DATE < (rise)){
    #   prof_data$isNight = "Y"}else{prof_data$isNight = "N"}}else{prof_data$isNight = NA}
    if(prof_data$SUN_ANGLE > -5 | ((DATE - set > 0) & (DATE - set <= 2.5*60*60))){prof_data$isNight = "N"}else{prof_data$isNight = "Y"}}else{
      if(prof_data$SUN_ANGLE > -5){prof_data$isNight = "N"}else{prof_data$isNight = "Y"}
    }
        
    prof_data$DATE = as.Date(DATE)
    prof_data$TIME = DATE
    prof_data$RISE = rise
    prof_data$SET = set
    
    ### ACCESS OCEAN VARIALES ### this could probably be simplified now using s-files
    if(!all(is.na(ncvar_get(f1,"PRES_ADJUSTED")))){
      
      if(BGC_col == 1){
        CTDSAL = as.vector(ncvar_get(f1,"PSAL_ADJUSTED"))
        CTDPRS = as.vector(ncvar_get(f1,"PRES_ADJUSTED"))
        CTDTMP = as.vector(ncvar_get(f1,"TEMP_ADJUSTED"))
        CHL = as.vector(ncvar_get(f1,"CHLA"))/2
        if(is.na(ncol(ncvar_get(f1,"BBP700")))){BBP = as.vector(ncvar_get(f1,"BBP700"))}
        if(any(grepl("CP660",names(f1$var)))){if(is.na(ncol(ncvar_get(f1,"CP660")))){CP = as.vector(ncvar_get(f1,"CP660"))}}
        CHL_P = CTDPRS
      }else{
        if(is.na(ncol(ncvar_get(f1,"CHLA")))){nc_close(f1)
          prof_data$QAnote = "no CHLA" 
          return(list(prof_data$QAnote,idx))}
        CTDSAL = as.vector(ncvar_get(f1,"PSAL_ADJUSTED")[,1])
        CTDPRS = as.vector(ncvar_get(f1,"PRES_ADJUSTED")[,1])
        CTDTMP = as.vector(ncvar_get(f1,"TEMP_ADJUSTED")[,1])
        CHL = as.vector(ncvar_get(f1,"CHLA")[,BGC_col])/2
        if(is.na(ncol(ncvar_get(f1,"BBP700")))){BBP = as.vector(ncvar_get(f1,"BBP700")[,BGC_col])}
        if(any(grepl("CP660",names(f1$var)))){if(is.na(ncol(ncvar_get(f1,"CP660")))){CP = as.vector(ncvar_get(f1,"CP660")[,BGC_col])}}
        CHL_P = as.vector(ncvar_get(f1,"PRES_ADJUSTED")[,BGC_col])}}else{
          if(BGC_col == 1){
            CTDSAL = as.vector(ncvar_get(f1,"PSAL"))
            CTDPRS = as.vector(ncvar_get(f1,"PRES"))
            CTDTMP = as.vector(ncvar_get(f1,"TEMP"))
            CHL = as.vector(ncvar_get(f1,"CHLA"))/2
            if(is.na(ncol(ncvar_get(f1,"BBP700")))){BBP = as.vector(ncvar_get(f1,"BBP700"))}
            if(any(grepl("CP660",names(f1$var)))){if(is.na(ncol(ncvar_get(f1,"CP660")))){CP = as.vector(ncvar_get(f1,"CP660"))}}
            CHL_P = CTDPRS
          }else{
            if(is.na(ncol(ncvar_get(f1,"CHLA")))){nc_close(f1)
              prof_data$QAnote = "no CHLA" 
              return(list(prof_data$QAnote,idx))}
            CTDSAL = as.vector(ncvar_get(f1,"PSAL")[,1])
            CTDPRS = as.vector(ncvar_get(f1,"PRES")[,1])
            CTDTMP = as.vector(ncvar_get(f1,"TEMP")[,1])
            CHL = as.vector(ncvar_get(f1,"CHLA")[,BGC_col])/2
            if(is.na(ncol(ncvar_get(f1,"BBP700")))){BBP = as.vector(ncvar_get(f1,"BBP700")[,BGC_col])}
            if(any(grepl("CP660",names(f1$var)))){if(is.na(ncol(ncvar_get(f1,"CP660")))){CP = as.vector(ncvar_get(f1,"CP660")[,BGC_col])}}
            CHL_P = as.vector(ncvar_get(f1,"PRES")[,BGC_col])}
        }
    
    if(length(which(!is.na(CHL))) == 0){nc_close(f1)
      prof_data$QAnote = "no CHL"
      return(list(prof_data$QAnote,idx))}
    nc_close(f1)
    
    
    # characterise front zone
    prof_data$FZ_insitu = front_class_insitu(CTDPRS, CTDTMP, CTDSAL,prof_data$LAT,prof_data$LON )
    
    # calculate density
    DENS = as.numeric(rho(S = CTDSAL, T = CTDTMP, P = 0))
    if(length(which(is.nan(DENS) == T)) > 0){DENS[which(is.nan(DENS))] = NA }
    if(length(which(is.finite(DENS) == T)) == 0){
      prof_data$QAnote = "no DENS" 
      return(list(prof_data$QAnote,idx))}
    
    #calculate depth
    DEPTH = swDepth(CTDPRS, latitude = prof_data$LAT)
    DENS_adj = DENS
    DENS_adj =  c(fit_bp_segmented(DEPTH[which(DEPTH < 500)], DENS[which(DEPTH < 500)], CD_thresh = 0)$fluor.out,DENS[DEPTH>=500])
    
    # reference temperature and density for calculations
    if(min(DEPTH,na.rm = T) > 10){
      DENS_10 = DENS_adj[!is.na(DENS_adj)][1]
      T_10 = CTDTMP[!is.na(CTDTMP)][1]}else{DENS_10 = approx(x = DEPTH, y = DENS_adj, xout = 10)$y
      T_10 = approx(x = DEPTH, y = CTDTMP, xout = 10)$y}
    
    if(length(which(is.finite(DENS) == T)) != 0){
    sigmatheta = swSigmaTheta(
      CTDSAL,
      CTDTMP,
      CTDPRS,
      referencePressure = 0,
      longitude = prof_data$LON,
      latitude = prof_data$LAT,
      eos = getOption("oceEOS", default = "gsw")
    )
    N2 = swN2(CTDPRS, sigmatheta)
      #as.numeric(gsw_Nsquared(CTDSAL[complete.cases(CTDSAL,CTDTMP)], CTDTMP[complete.cases(CTDSAL,CTDTMP)],CTDPRS[complete.cases(CTDSAL,CTDTMP)],prof_data$LAT)$N2)
    prof_data$MLD_N2 = CTDPRS[complete.cases(CTDSAL,CTDTMP)][which.max(N2)]
    rm(N2,sigmatheta)}else{prof_data$MLD_N2 = NA}
   Spice = swSpice(salinity = CTDSAL, temperature =  CTDTMP, pressure = CTDPRS, longitude = prof_data$LON,latitude = prof_data$LAT)
    # physical structure calculations 
    prof_data$MLD = approx(x = DENS_adj[DEPTH>=10], y = DEPTH[DEPTH>=10], xout = DENS_10 + 0.03,rule = 2)$y
    if(is.na(prof_data$MLD)){
      prof_data$QAnote = "no DENS in top 10 m" 
      return(list(prof_data$QAnote,idx))}
    
    prof_data$TML_diff = CTDTMP[which(is.finite(CTDTMP))[1]] - min(CTDTMP, na.rm = T)
    prof_data$TML_depth = DEPTH[which.min(CTDTMP)]
    prof_data$TML_depth2 = DEPTH[which.min(CTDSAL)]
    prof_data$FWL = prof_data$TML_depth < prof_data$MLD & prof_data$TML_depth2 < prof_data$MLD
    prof_data$TML =  prof_data$TML_diff > 0.2 & prof_data$TML_depth < 150 & prof_data$TML_depth > prof_data$MLD
    
    
    prof_data$TLD = TLD(prof_data$MLD,DENS_adj, DEPTH)
    # relic from merged profiles
    if(BGC_col == 1){CHL_D = DEPTH}else{CHL_D =  swDepth(CHL_P, latitude = prof_data$LAT)}
    
    # CHL = CHL[CHL_D > 5]
    # BBP = BBP[CHL_D > 5]
    # CHL_D = CHL_D[CHL_D > 5]
    # CHL_QC
    # no NPQ as only night time prof (2 hours after sunset + 2 hours before sunrise)
    # res considered (only <5m)
    
    # quality control
    if(sd(CHL, na.rm = T) == 0){print("skipped bad prof")
      prof_data$QAnote = "no CHL var"
      return(list(prof_data$QAnote,idx))} # no variance
    CHL_old = CHL[CHL_D < 500 & CHL_D > 0]
    # smooth CHL and BBP
    BBP_old = BBP[CHL_D < 500 & CHL_D > 0]
    BBP_min = runmin(BBP[CHL_D < 500 & CHL_D > 0],5,endrule = "constant")
    BBP_max = runmax(BBP[CHL_D < 500 & CHL_D > 0],5,endrule = "constant")
    
    CHL = abstract_fluor(CHL_D,CHL)
    
    if(!exists("BBP")){BBP = NA}else{
      BBP = abstract_bbp(CHL_D,BBP, CHL)
      BBP = BBP[CHL_D < 500 & CHL_D > 0]
      }
    if(!exists("CP")){CP = NA}else{
      CP_old = CP[CHL_D < 500 & CHL_D > 0]
      CP = abstract_fluor(CHL_D,CP)
      CP = CP[CHL_D < 500 & CHL_D > 0]
      }
    
    # quality control
    CHL_D = CHL_D[CHL_D < 500]
    CHL = CHL[CHL_D > 0]
    CHL_D = CHL_D[CHL_D > 0]
    
    emld = Eco_MLD(CHL_D, CHL)
    prof_data$EMLD = emld$EMLD
    prof_data$EMLD_QI = emld$QI
    prof_data$CHL50 = CHL_50(CHL_D, CHL)
    prof_data$EMLD2 = emld$EMLD
    prof_data$max_d = max(CHL_D[!is.na(CHL_old)], na.rm = T)
    if(is.na(emld$EMLD)){
      prof_data$EMLD = prof_data$CHL50 
      prof_data$EMLD_flag_changed = T
      }
     if((max(CHL[CHL_D < 0.25*prof_data$EMLD], na.rm = T) > (1.5*min(CHL[CHL_D > 0.25*prof_data$EMLD & CHL_D < prof_data$EMLD], na.rm = T)) & CHL_D[which.max(CHL)[1]] < 0.25*prof_data$EMLD)){
       if(prof_data$EMLD/prof_data$CHL50 > 1.5  | (prof_data$EMLD - prof_data$CHL50) > 50){prof_data$EMLD = prof_data$CHL50 
       prof_data$HighSF_flag_changed = T}
       prof_data$HighSF_flag = T
     }
    if(max(CHL[CHL_D < 0.025*prof_data$EMLD], na.rm = T) > (1.5*min(CHL[CHL_D > 0.025*prof_data$EMLD & CHL_D < 0.1*prof_data$EMLD], na.rm = T)) & (max(CHL[CHL_D < 0.025*prof_data$EMLD], na.rm = T)/max(CHL[CHL_D < prof_data$EMLD], na.rm = T) > 0.5)){
      if(prof_data$EMLD/prof_data$CHL50 > 1.5  | (prof_data$EMLD - prof_data$CHL50) > 50){prof_data$EMLD = prof_data$CHL50 
      prof_data$HighSF_flag_changed = T}
      prof_data$HighSF_flag2 = T
    }
    
    
    CHL_D2 = CHL_D[is.finite(CHL)]
    prof_data$min_d1 = min(CHL_D2, na.rm = T)

    prof_data$min_d2 = min(CHL_D2, na.rm = T)
    
    prof_data$res_50 = mean(CHL_D2[2:(which.closest(CHL_D2,50)[1])] - CHL_D2[1:((which.closest(CHL_D2,50))[1]-1)], na.rm = T)
    prof_data$res_100 = mean(CHL_D2[(which.closest(CHL_D2,75)[1] + 1):(which.closest(CHL_D2,100)[1])] - CHL_D2[(which.closest(CHL_D2,75)[1]):(which.closest(CHL_D2,100)[1] -1)], na.rm = T)
    
    prof_data$Fvar = sd(CHL - CHL_old, na.rm = T)/diff(range(CHL, na.rm = T))
    if( prof_data$Fvar > 0.2){print("skipped bad prof")
      prof_data$QAnote = paste("high CHL var. Fvar =",prof_data$Fvar)
      return(list(prof_data$QAnote,idx))}
    
    if(prof_data$Fvar == 0 ){print("skipped bad prof")
      prof_data$QAnote = "Abstracted profile no different to measured"
      return(list(prof_data$QAnote,idx))}
    
    if(CHL[which(!is.na(CHL))[1]] > 5*median(CHL[which(!is.na(CHL))[3:7]], na.rm = T)){print("skipped bad prof")
      prof_data$QAnote = "questionably high surface chl"
      return(list(prof_data$QAnote,idx))}
    prof_data$LowCHL = ifelse(max(CHL , na.rm = T) < 0.05,T,F)
    
    # if(any(CHL[CHL_D < prof_data$EMLD] < min(CHL[CHL_D > prof_data$EMLD], na.rm = T), na.rm = T)){print("skipped bad prof")
    #   prof_data$QAnote = "negative CHL above EMLD"
    #   return(list(prof_data$QAnote,idx))}
    
    # if(any(CHL[CHL_D < prof_data$EMLD] < min(CHL[CHL_D > prof_data$EMLD], na.rm = T), na.rm = T)){print("skipped bad prof")
    #   prof_data$QAnote = "negative CHL above EMLD"
    #   return(list(prof_data$QAnote,idx))}
    
    ### Optical calucilations
    prof_data$Zeu2 = Zeu(CHL_D, CHL)
    prof_data$Zeu0.15 = Zeu(CHL_D, CHL,0.15)
    prof_data$Zeu0.2 = Zeu(CHL_D, CHL,0.2)
    prof_data$Zeu0.1 = Zeu(CHL_D, CHL,0.1)
    prof_data$Zeu0.25 = Zeu(CHL_D, CHL,0.25)
    prof_data$Zeu0.05 = Zeu(CHL_D, CHL,0.05)
    prof_data$SFM_d = CHL_D[which.max(CHL)]
    prof_data$SBM_d = CHL_D[which.max(BBP)]
    prof_data$DENS_SFM = approx(x = DEPTH, y = DENS_adj, xout = prof_data$SFM_d)$y
    prof_data$TEMP_SFM = approx(x = DEPTH, y = CTDTMP, xout = prof_data$SFM_d)$y
    prof_data$SAL_SFM = approx(x = DEPTH, y = CTDSAL, xout = prof_data$SFM_d)$y
    prof_data$DENS_surf = approx(x = DEPTH, y = DENS_adj, xout = 0,rule = 2)$y
    prof_data$TEMP_surf = approx(x = DEPTH, y = CTDTMP, xout = 0,rule = 2)$y
    prof_data$SAL_surf = approx(x = DEPTH, y = CTDSAL, xout = 0,rule = 2)$y
    prof_data$SPI_surf = approx(x = DEPTH, y = Spice, xout = 0,rule = 2)$y
    prof_data$SPI_SFM = approx(x = DEPTH, y = Spice, xout = prof_data$SFM_d)$y
    prof_data$SPI_MLD = approx(x = DEPTH, y = Spice, xout = prof_data$MLD)$y
    prof_data$SPI_10 = approx(x = DEPTH, y = Spice, xout = 10)$y
    if(prof_data$isNight == "N"){
    NPQ = NPQ_SCC(CHL_D, CHL, BBP,prof_data$MLD, prof_data$FWL, prof_data$SFM_d)
    CHL_cor = NPQ$corr_fluor
    prof_data$NPQ_d = NPQ$NPQdepth}else{
      CHL_cor = CHL
      prof_data$NPQ_d = NA
    }
    prof_data$CHL_surf = ifelse(min(CHL_D,na.rm = T)<20,mean(CHL_cor[CHL_D<20],na.rm = T),CHL_cor[which(!is.na(CHL_cor))[1]])
    prof_data$CHL_surf1 = CHL_cor[which(!is.na(CHL_cor))[1]]
    if(min(CHL_D,na.rm = T)<20 & is.nan(prof_data$CHL_surf)){
      print("skipped bad prof")
      prof_data$QAnote = "No CHL measurement in top 20 m"
      return(list(prof_data$QAnote,idx))
    }
    
    prof_data$CHL_15m = ifelse(min(CHL_D,na.rm = T)<15,median(CHL_cor[CHL_D<15],na.rm = T),NA)
    prof_data$CHL_surf30 = ifelse(min(CHL_D,na.rm = T)<30,mean(CHL_cor[CHL_D<30],na.rm = T),CHL_cor[which(!is.na(CHL_cor))[1]])
    prof_data$BBP_surf = ifelse(min(CHL_D,na.rm = T)<20,mean(BBP[CHL_D<20],na.rm = T),BBP[which(!is.na(BBP))[1]])
    
    if(!all(is.na(CP))){prof_data$SPM_d = CHL_D[which.max(CP)]
    prof_data$CP_surf = ifelse(min(CHL_D,na.rm = T)<20,mean(CP[CHL_D<20],na.rm = T),CP[which(!is.na(CP))[1]])
    prof_data$SPM_val = CP[which.max(CP)]}else{
      prof_data$SPM_d = NA
      prof_data$CP_surf = NA
      prof_data$SPM_val = NA
    }
    ucCHL_surf = ifelse(min(CHL_D,na.rm = T)<20,mean(CHL[CHL_D<20],na.rm = T),CHL[1])  
    prof_data$NPQsignal = prof_data$CHL_surf - ucCHL_surf
    
    prof_data$NPQ = ifelse(prof_data$CHL_surf > 1, ( prof_data$NPQsignal/prof_data$CHL_surf ) > 0.1, prof_data$NPQsignal > 0.1)
    prof_data$Zeu = Zeu(CHL_D, CHL_cor)
    prof_data$SFM_val = CHL_cor[which.max(CHL_cor)]
    prof_data$SBM_val = BBP[which.max(BBP)]
    
    prof_data$SFM_bbp = BBP[which.max(CHL_cor)]
    prof_data$BBP_R_surf = ifelse(min(CHL_D,na.rm = T)<20,mean((BBP/CHL_cor)[CHL_D<20],na.rm = T),BBP[which(!is.na(CHL_cor))[1]]/CHL_cor[which(!is.na(CHL_cor))[1]])
    prof_data$BBP_R_scm = (BBP/CHL_cor)[which.max(CHL_cor)]
    prof_data$SFM_thick = NA
    prof_data$SFM = NA
    if(prof_data$SFM_d <= 20){prof_data$SFM = FALSE}
    if(prof_data$SFM_d  > 20){
      if((prof_data$CHL_surf + 0.2) < prof_data$SFM_val & prof_data$CHL_surf < 1){
        prof_data$SFM = TRUE
      }else{if(1.2*prof_data$CHL_surf  < prof_data$SFM_val &  prof_data$CHL_surf > 1){
        prof_data$SFM = TRUE
      }else{prof_data$SFM  = FALSE}}
      if(prof_data$SFM){
        
        t_depth2 = CHL_D[CHL_D > prof_data$SFM_d][which.closest(CHL_cor[CHL_D > prof_data$SFM_d], prof_data$CHL_surf + ((prof_data$SFM_val - prof_data$CHL_surf)/2 ))[1]]
        
        t_depth1 = CHL_D[CHL_D < prof_data$SFM_d][which.closest(CHL_cor[CHL_D < prof_data$SFM_d], prof_data$CHL_surf + ((prof_data$SFM_val - prof_data$CHL_surf)/2 ))[1]]
        
        prof_data$SFM_thick = (prof_data$SFM_d - t_depth1[which.closest(t_depth1,prof_data$SFM_d)]) + (t_depth2[which.closest(t_depth2,prof_data$SFM_d)] - prof_data$SFM_d)
        rm(t_depth1, t_depth2)}else{prof_data$SFM_thick  = NA}}else{prof_data$SFM  = FALSE}
    
    
    
    a1 = approxfun(CHL_D,CHL_cor, rule = 2) 
    prof_data$int_CHL = integrate(a1, min(CHL_D,na.rm = T), 300, subdivisions = 500)$value
    
    prof_data$SFM_mag = prof_data$SFM_val/prof_data$CHL_surf
    
    # comparison
    prof_data$SFM_Cornec = ifelse(is.na(prof_data$CHL_15m), NA, ifelse(prof_data$SFM_val > 2*prof_data$CHL_15m, "SCM","no SCM"))

    # Initial parameters for model fits - Carranza and Mignot
    Fsurf = prof_data$CHL_surf30
    Z0.5 = max(prof_data$CHL50, prof_data$MLD)
    slope = approx(x = CHL_D, y = CHL_cor, xout = Z0.5 + 0.5, rule = 2)$y - approx(x = CHL_D, y = CHL_cor, xout = Z0.5 - 0.5, rule = 2)$y
    peak_ids = localMaxima(CHL_cor[is.finite(CHL_cor)])
    if(length(peak_ids) == 0){peak_ids = 20}
    peak_vals = CHL_cor[is.finite(CHL_cor)][peak_ids]
    peak_depths = CHL_D[is.finite(CHL_cor)][peak_ids]
    Zmax = peak_depths[which.max(peak_vals)]
    Fmax = peak_vals[which.max(peak_vals)]
    if(Fmax < 0.02){
      Fmax = prof_data$SFM_val
      Zmax = prof_data$SFM_d
    }
    dz = 5


    Ma_fit = optim(par = c(Fsurf,Z0.5,slope) , fn = function(x){sum((CHL_cor - (x[1]/(1 + exp((x[2]-CHL_D)*x[3]))))^2,na.rm = T)},method =  "L-BFGS-B", lower = c(0,1,-500), upper = 500)
    Mb_fit = optim(par = c(Fsurf,Z0.5) , fn = function(x){sum((CHL_cor - (x[1]*(exp(-CHL_D*(log(CHL_D)/x[2])))))^2,na.rm = T)},method =  "L-BFGS-B", lower = c(0,1), upper = 500)
    Mc_fit = optim(par = c(Fmax,Zmax,dz) , fn = function(x){sum((CHL_cor - (x[1]*(exp(-((CHL_D - x[2])^2)/(x[3]^2)))))^2,na.rm = T)},method =  "L-BFGS-B", lower = c(0,0.0001), upper = 500)
    Md_fit = optim(par = c(Fsurf,Z0.5,Fmax,Zmax,dz) , fn = function(x){sum((CHL_cor - (x[1]*(exp(-CHL_D*(log(CHL_D)/x[2])))) - x[3]*(exp(-((CHL_D - x[4])^2)/(x[5]^2))))^2,na.rm = T)},method =  "L-BFGS-B", lower = c(0,1,0.00001,0.00001,0.00001), upper = 500)
    Me_fit = optim(par = c(Fsurf,Z0.5,slope,Fmax,Zmax,dz) , fn = function(x){sum((CHL_cor - (x[1]/(1 + exp((x[2]-CHL_D)*x[3]))) -  x[4]*(exp(-((CHL_D - x[5])^2)/(x[6]^2))))^2,na.rm = T)},method =  "L-BFGS-B", lower = c(0,1,-500,0.00001,0.00001,0.000001), upper = 500)

    res_ma = (Ma_fit$par[1]/(1 + exp((Ma_fit$par[2]-CHL_D)*Ma_fit$par[3])))
    res_mb = (Mb_fit$par[1]*(exp(-CHL_D*(log(CHL_D)/Mb_fit$par[2]))))
    res_mc = (Mc_fit$par[1]*(exp(-((CHL_D - Mc_fit$par[2])^2)/(Mc_fit$par[3]^2))))
    res_md = (Md_fit$par[1]*(exp(-CHL_D*(log(CHL_D)/Md_fit$par[2])))) + Md_fit$par[3]*(exp(-((CHL_D - Md_fit$par[4])^2)/(Md_fit$par[5]^2)))
    res_me = (Me_fit$par[1]/(1 + exp((Me_fit$par[2]-CHL_D)*Me_fit$par[3]))) +  Me_fit$par[4]*(exp(-((CHL_D - Me_fit$par[5])^2)/(Me_fit$par[6]^2)))

    chi_a = sum((CHL_cor - res_ma)^2/(0.0016), na.rm = T)
    crit_a = qchisq(0.95,length(which(!is.na(CHL_cor))) - 3)
    chi_b = sum((CHL_cor - res_mb)^2/(0.0016), na.rm = T)
    crit_b = qchisq(0.95,length(which(!is.na(CHL_cor))) - 2)
    chi_c = sum((CHL_cor - res_mc)^2/(0.0016), na.rm = T)
    crit_c = qchisq(0.95,length(which(!is.na(CHL_cor))) - 3)
    chi_d = sum((CHL_cor - res_md)^2/(0.0016), na.rm = T)
    crit_d = qchisq(0.95,length(which(!is.na(CHL_cor))) - 5)
    chi_e = sum((CHL_cor - res_me)^2/(0.0016), na.rm = T)
    crit_e = qchisq(0.95,length(which(!is.na(CHL_cor))) - 6)

    chi = c(chi_a, chi_b, chi_c, chi_d, chi_e)
    crit = c(crit_a, crit_b, crit_c, crit_d, crit_d)
    names = c("sigmoid","exponential","gaussian","gaussian+exponential","gaussian+sigmoid")
    prof_data$Gauss_fit_r2 = max(cor(CHL_cor , res_mc,use = "complete.obs")^2,cor(CHL_cor , res_md,use = "complete.obs")^2,cor(CHL_cor , res_me,use = "complete.obs")^2)

    ### chi square to choose
    prof_data$shape = ifelse(chi[which.min(chi)] < crit[which.min(chi)],names[which.min(chi)],NA)
    prof_data$isGaussian = ifelse(prof_data$shape %in% c("gaussian","gaussian+exponential","gaussian+sigmoid") & prof_data$Gauss_fit_r2 > 0.8, T, F)

    prof_data$SFM_Carranza = ifelse(prof_data$isGaussian & (prof_data$SFM_val - prof_data$CHL_surf30 > 0.04) & (Zmax > 20), "SCM", "no SCM")

    prof_data$SFM_Mignot = ifelse(prof_data$isGaussian & prof_data$Gauss_fit_r2 > 0.8, "SCM", "no SCM")
    rm(crit,crit_a, chi, crit_b, crit_c, crit_d, crit_e, chi_a, chi_b, chi_c, chi_d, chi_e, res_ma, res_mb, res_mc, res_md, res_me, Ma_fit, Mb_fit, Mc_fit, Md_fit, Me_fit, Fsurf,Z0.5,slope,Fmax,Zmax,dz)

    # fPCA matrix
    depths = seq(0,2,0.0005)
    depths2 = seq(0,500,0.125)
    if(exists("BBP")){
      if(!all(is.na(CP))){
      data = cbind(approx(CHL_D/prof_data$EMLD,CHL_cor,depths,rule = 2)$y,approx(DEPTH,CTDSAL,depths2,rule = 2)$y,approx(DEPTH,CTDTMP,depths2,rule = 2)$y,approx(CHL_D/prof_data$EMLD,BBP,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,CHL,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,CHL_old,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,BBP_old,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,BBP_min,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,CP,depths,rule = 2)$y)
      }else{
        data = cbind(approx(CHL_D/prof_data$EMLD,CHL_cor,depths,rule = 2)$y,approx(DEPTH,CTDSAL,depths2,rule = 2)$y,approx(DEPTH,CTDTMP,depths2,rule = 2)$y,approx(CHL_D/prof_data$EMLD,BBP,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,CHL,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,CHL_old,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,BBP_old,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,BBP_min,depths,rule = 2)$y,NA)
      }
      }else{
        data = cbind(approx(CHL_D/prof_data$EMLD,CHL_cor,depths,rule = 2)$y,approx(DEPTH,CTDSAL,depths2,rule = 2)$y,approx(DEPTH,CTDTMP,depths2,rule = 2)$y,NA,approx(CHL_D/prof_data$EMLD,CHL,depths,rule = 2)$y,approx(CHL_D/prof_data$EMLD,CHL_old,depths,rule = 2)$y,NA,NA,NA)
      }
    rm(CP)
    system(paste('echo', idx))
    return(list(prof_data,data))
    
  },error = function(e){
    return( list(paste("Note",idx,"returned error:",e),idx))
    
    #return(NULL)
  })
}


### Run calculations in parallel
# set up run parameters
# remove objects to save memory
rm(list=setdiff(ls(), c("idx_list","Argo_calcs","A_files","dfiles","ids")))

#idx_list = which(BGCArgoResults$HighSF_flag2 == T)
numJobs = length(idx_list)
numCores <- min(numJobs, 8)
#theCluster <- parallel::makeForkCluster(getOption("cl.cores", numCores))
#discard_output <- parallel::clusterApplyLB(theCluster
#                                               ,x=idx_list
#                                               ,fun=Argo_calcs)                 
#parallel::stopCluster(theCluster)       


for(sp in 1:84){
  idx_list2 = idx_list[which(idx_list <= sp*500 & idx_list > (sp-1)*500)]
  theCluster <- parallel::makeForkCluster(getOption("cl.cores", numCores))
  discard_output_int <- parallel::clusterApplyLB(theCluster
                                           ,x=idx_list2
                                           ,fun=Argo_calcs)
  if(sp == 1){discard_output = discard_output_int}else{
    discard_output = append(discard_output_int, discard_output)
  }
  save(discard_output, file = "discard_output_final.RData")
  parallel::stopCluster(theCluster)
}


 ### extract outputs from paralell calculations
library(abind)

discard_output1 =  sapply(discard_output, "[[", 1)
QA_idx = unlist(lapply(discard_output1,FUN = function(x){!is.data.frame(x)}))
QA_notes =  unlist(discard_output1[QA_idx])
#discard_output1 = discard_output1[!grepl("Note",discard_output1)]
discard_output1 = discard_output1[unlist(lapply(discard_output1,FUN = function(x){is.data.frame(x)}))]

BGCArgoResults = do.call('rbind', discard_output1) # BGC-Argo calculation results
#BGCArgoResults$LowCHL = ifelse( BGCArgoResults$SFM_val < 0.1,T,F)



discard_output2 = sapply(discard_output, "[[", 2)
QA_id = discard_output2[QA_idx]
QA_id[sapply(QA_id, is.null)] <- NA
QA_id = unlist(QA_id)
discard_output2 = discard_output2[unlist(lapply(discard_output2,FUN = function(x){!is.function(x) & !is.null(x) & length(x) > 1}))]
PCAData = abind(discard_output2,along = 3) # Data to use in fPC analysis
PCAData = aperm(PCAData,c(1,3,2))

QA_notes = data.frame(id = QA_id, note = QA_notes)


# tab = sapply(str_split(QA_notes$note,":",2),"[",2)

#### QA to remove profiles
# remove empty profiles, profiles with less than 50 filed depth levels, 
# profiles in the STZ, profles with surface chl >10 , profiles with surface chl < 0
p = seq(0,2,0.0005)
rel_file = read.table("releaseFile")
colnames(rel_file) = c("ID", "LON", "LAT", "n1","n2","YYYY","MM","DD","HH")
rel_file$rel_match = unlist(lapply(1:nrow(rel_file),function(x){paste(rel_file[x,c("LAT","LON","YYYY","MM","DD")], collapse = " ")}))
BGCArgoResults$rel_match = unlist(lapply(1:nrow(BGCArgoResults),function(x){paste(c(round(BGCArgoResults$LAT[x],digits = 3),round(BGCArgoResults$LON[x],digits = 3),BGCArgoResults[x,c("YYYY","MM","DD")]), collapse = " ")}))

remove_profiles = NULL
# count_chl = NULL
# count_tmp = NULL
# count_sal = NULL
# for(i in 1:dim(PCAData)[2]){
#   count_chl = c(count_chl,length(which(is.finite(PCAData[,i,1]))))
#   count_sal = c(count_sal,length(which(is.finite(PCAData[,i,2]))))
#   count_tmp =c(count_tmp,length(which(is.finite(PCAData[,i,3]))))
#   if(all(is.na(PCAData[,i,2])) | all(is.na(PCAData[,i,3])) | length(which(!is.na(PCAData[,i,3]))) < 50 | length(which(!is.na(PCAData[,i,3]))) < 50){
#     #remove_profiles = c(remove_profiles,i)
#   }else{
#     PCAData[,i,2] = approx(p,PCAData[,i,2],p,rule = 2)$y
#     PCAData[,i,1] = approx(p,PCAData[,i,1],p,rule = 2)$y
#     PCAData[,i,3] = approx(p,PCAData[,i,3],p,rule = 2)$y}
#   #if(PCAData[,i,1][which(!is.na(PCAData[,i,1]))[1]] > 5*min(PCAData[,i,1][which(!is.na(PCAData[,i,1]))[3:7]], na.rm = T)){remove_profiles = c(remove_profiles,i)}
# }
# 332 rm
rm_bgc_prof = NULL
surf_spike = NULL
PCAData_n = PCAData[,,1]
#PCAData_1 = PCAData_1[,!BGCArgoResults_PCA$LowCHL] # remove low chlorophyll profiles
for(i in 1:dim(PCAData)[2]){
  
  PCAData_n[,i] = PCAData[,i,1]/max(PCAData[,i,1],na.rm = T) # scale to maximum chl value of 1
  
  #if(any(PCAData_n[1:1600,i] < 0)){
    #plot(PCAData_1[,i], main = i) # test for negative values in top 80% of EMLD
    #rm_bgc_prof = c(rm_bgc_prof,i)
    #}
  #if(any(PCAData_n[2000:3000,i] < -0.2)){
    #plot(PCAData_1[,i], main = i) # test for negative values in middle of EMLD
    #rm_bgc_prof = c(rm_bgc_prof,i)
   # }
  #if(length(which(PCAData_n[2000:4000,i] < -0.1)) > 500){
    #plot(PCAData_1[,i], main = i) # test for negative values below EMLD
    #rm_bgc_prof = c(rm_bgc_prof,i)
  #}
  if(any(PCAData_n[3000:4000,i] > 0.8)){
    #plot(PCAData_1[,i], main = i) # test for negative values below EMLD
    rm_bgc_prof = c(rm_bgc_prof,i)
  }
  if(max(PCAData_n[0:50,i]) > (2*min(PCAData_n[50:200,i])) & max(PCAData_n[0:50,i]) > 0.5){
    #plot(PCAData_n[,i], main = i) # test for surface spikes
    surf_spike = c(surf_spike,i)
  }
  
}
BGCArgoResults$spike_surf = FALSE
BGCArgoResults$spike_surf[surf_spike] = TRUE
# the following QA was tagged using the fPCA analysis from the SCM ID. We saved the data to remove the spiky cluster, preserving thin SCMs
#load("QAbbp.RData")
#load("QAbbp2.RData")
#load("QAbbp4.RData")
#load("QAbbp_temp.RData")
#spiky = setdiff(bad_pf_remove,thin_SCM_info)
#spiky_rm = which(BGCArgoResults$DIRECTION %in% spiky$DIRECTION & BGCArgoResults$float %in% spiky$float & BGCArgoResults$JULD %in% spiky$JULD & round(BGCArgoResults$LAT,3) %in% round(spiky$LAT,3) & round(BGCArgoResults$LON,3) %in% round(spiky$LON,3))
#spiky_rm2 = which(BGCArgoResults$DIRECTION %in% bad_pf_remove2$DIRECTION & BGCArgoResults$float %in% bad_pf_remove2$float & BGCArgoResults$JULD %in% bad_pf_remove2$JULD & round(BGCArgoResults$LAT,3) %in% round(bad_pf_remove2$LAT,3) & round(BGCArgoResults$LON,3) %in% round(bad_pf_remove2$LON,3))
#spiky_rm4 = which(BGCArgoResults$DIRECTION %in% bad_pf_remove4$DIRECTION & BGCArgoResults$float %in% bad_pf_remove4$float & BGCArgoResults$JULD %in% bad_pf_remove4$JULD & round(BGCArgoResults$LAT,3) %in% round(bad_pf_remove4$LAT,3) & round(BGCArgoResults$LON,3) %in% round(bad_pf_remove4$LON,3))
#spiky_rm_temp = which(BGCArgoResults$DIRECTION %in% bad_pf_remove3$DIRECTION & BGCArgoResults$float %in% bad_pf_remove3$float & BGCArgoResults$JULD %in% bad_pf_remove3$JULD & round(BGCArgoResults$LAT,3) %in% round(bad_pf_remove3$LAT,3) & round(BGCArgoResults$LON,3) %in% round(bad_pf_remove3$LON,3))
#rerun_idx = which( dfiles$PLATFORM_CODE %in% bad_pf_remove3$float & round(dfiles$latitude,3) %in% round(bad_pf_remove3$LAT,3) & round(dfiles$longitude,3) %in% round(bad_pf_remove3$LON,3))
#save(rerun_idx,file = "Rerun_idx.RData")

### check surface threshold again
remove_profiles = unique(c(surf_spike,rm_bgc_prof, which(BGCArgoResults$FZ_insitu == "STZ"), which(BGCArgoResults$TEMP_surf >= 11),which(BGCArgoResults$CHL_surf > 20),which(BGCArgoResults$CHL_surf < 0),which(PCAData[1,,1] < 0), which(PCAData[1,,1] > 20))) # spiky_rm,spiky_rm2,spiky_rm4,which(BGCArgoResults$HighSF_flag),
PCAData_final = PCAData[,-remove_profiles,]
BGCArgoResults_PCA = BGCArgoResults[-remove_profiles,]
## need some way of keeping really low chl.

# save data as back up
save(BGCArgoResults,PCAData,remove_profiles, file = "./BGC-ArgoDatasets.RData")
save(PCAData_final,BGCArgoResults_PCA, file = "./PCA.RData")
save(discard_output, file = "./BGCArgo_calcs_output.RData")