### directory containing daily files

outdir=paste("daily/",tile,"/",sep="")  #directory for separate daily files

### directory to hold climatology

outdir2="summary" #directory for combined daily files and summarized files

## Get list of all daily files

if(verbose) print("Checking daily output in preparation for generating climatology")

  fdly$dateid=substr(fdly$file,14,21)

  fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")

  fdly$month=format(fdly$date,"%m")

nrow(fdly)

## check validity (via npar and ntime) of nc files

## print some summaries

if(verbose) print("Summary of available daily files")

## Identify which files failed test

## delete files that fail check?

if(delete) {

  print(paste(sum(fdly$drop),"files will be deleted"))

  file.remove(as.character(fdly$path[fdly$drop]))

}

## remove dropped files from list

#################################################################################

## Combine the year-by-year files into a single daily file in the summary directory (for archiving)

myear=as.integer(max(fdly$year))  #this year will be used in all dates of monthly climatologies (and day will = 15)

## Monthly means

system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmean.nc",sep=""),wait=T)

## Monthly standard deviation

system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))

system(paste(ncopath,"ncrename -v CER,CER_sd -v CLD,CLD_sd -v COT,COT_sd ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))

system(paste(ncopath,"ncatted ",

tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))

## cer > 20

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -ymonmean -gtc,20 -selvar,CER ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))

system(paste(ncopath,"ncrename -v CER,CER20 ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a units,CER20,o,c,\"Proportion\" ",

tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))

## number of observations

system(paste("cdo -O sorttimestamp  -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -eqc,9999 -setmisstoc,9999   -selvar,CER,CLD ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))

system(paste(ncopath,"ncrename -v CER,CER_pmiss -v CLD,CLD_pmiss ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))

system(paste(ncopath,"ncatted ",

system(paste(ncopath,"ncks -O ",tsdir,"/MOD06_",tile,"_ymonmean.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonstd.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncer20.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonmiss.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))

## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes

     }

date="20030102"

tile="h11v08"

platform="pleiades" 

outdir=paste("daily/",tile,"/",sep="")  #directory for separate daily files

  ## location of MOD06 files

  ## path to some executables

  ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"

  swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif"

  projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "

}

  ## specify working directory

  setwd("~/acrobates/projects/interp")

  gisBase="/usr/lib/grass64"

  "Cloud_Effective_Radius_Uncertainty",  "CERU",

  "Cloud_Optical_Thickness",             "COT",

  "Cloud_Optical_Thickness_Uncertainty", "COTU",

  "Cloud_Mask_1km",                      "CM1",

  "Quality_Assurance_1km",               "QA"),

## vector of variables expected to be in final netcdf file.  If these are not present, the file will be deleted at the end.

finalvars=c("CER","COT","CLD")

                CM_clear_",i," =  ((CM1_",i," / 2^0) % 2) == 1  &&  ((CM1_",i," / 2^1) % 2^2) >  2 

                CM_path_",i," =   ((CM1_",i," / 2^6) % 2^2) 

                CM_cloud2_",i," = ((CM1_",i," / 2^1) % 2^2) 

EOF",sep=""))

    ## QA

     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),

         CER2_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"

         CER2_daily=int(CER2_numer/CER2_denom)

         CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),-9999,CM_cloud_",1:nfs,")",sep="",collapse=","),"))) 

EOF",sep=""))

  ### Write the files to a netcdf file

  ## create image group to facilitate export as multiband netcdf

  if(file.exists(ncfile)) file.remove(ncfile)  #if it exists already, delete it

  execGRASS("r.out.gdal",input="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),

system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))

system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))

## add other attributes

  system(paste(ncopath,"ncatted -a scale_factor,CER,o,d,0.01 -a units,CER,o,c,\"micron\" -a missing_value,CER,o,d,-32768 -a long_name,CER,o,c,\"Cloud Particle Effective Radius\" ",ncfile,sep=""))

  system(paste(ncopath,"ncatted -a scale_factor,CER2,o,d,0.01 -a units,CER2,o,c,\"micron\" -a missing_value,CER2,o,d,-32768 -a long_name,CER2,o,c,\"Cloud Particle Effective Radius with clear sky set to zero\" ",ncfile,sep=""))

system(paste(ncopath,"ncatted -a scale_factor,COT,o,d,0.01 -a units,COT,o,c,\"none\" -a missing_value,COT,o,d,-32768 -a long_name,COT,o,c,\"Cloud Optical Thickness\" ",ncfile,sep=""))

system(paste(ncopath,"ncatted -a scale_factor,COT2,o,d,0.01 -a units,COT2,o,c,\"none\" -a missing_value,COT2,o,d,-32768 -a long_name,COT2,o,c,\"Cloud Optical Thickness with clear sky set to zero\" ",ncfile,sep=""))

### delete the temporary files 

#### Script to facilitate processing of MOD06 data

  setwd("/nobackupp1/awilso10/mod06")

library(rgdal)

verbose=T

## get MODLAND tile information

tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")

save(tb,file="modlandTiles.Rdata")

load("modlandTiles.Rdata")

### list of tiles to process

tiles=c("h11v08","h21v09","h08v04","h09v04","h08v05","h09v05","h20v11","h31v11")

tile_bb=tb[tb$tile%in%tiles,]

### get list of files to process

## identify which have been completed

fdone=data.frame(path=list.files(outdir,pattern="nc$",recursive=T))

fdone$date=substr(basename(as.character(fdone$path)),14,21)

## identify which date-tiles have already been run

proclist$done=paste(proclist$tile,proclist$date,sep="_")%in%substr(basename(as.character(fdone$path)),7,21)

### report on what has already been processed

table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done])

script="/u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r"

## write the table processed by mpiexec

write.table(paste("--verbose ",script," --date ",proclist$date[!proclist$done]," --verbose T --tile ",proclist$tile[!proclist$done],sep=""),

CORES=400

HDIR=/u/armichae/pr/

  source /u/awilso10/.bashrc

IDIR=/nobackupp1/awilso10/mod06/

##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt

### Now submit the script to generate the climatologies

tiles

climatescript="/u/awilso10/environmental-layers/climate/procedures/MOD06_Climatology.r"

## write the table processed by mpiexec

write.table(paste("--verbose ",climatescript," --verbose T --tile ",ctiles,sep=""),

file=paste("notdone_climate.txt",sep=""),row.names=F,col.names=F,quote=F)

### qsub script

cat(paste("

#PBS -S /bin/bash

##PBS -l select=2:ncpus=4:mpiprocs=4

#PBS -l walltime=5:00:00

#PBS -j n

#PBS -N mod06_climate

#PBS -q normal

#PBS -V

HDIR=/u/armichae/pr/

  source /u/awilso10/.bashrc

IDIR=/nobackupp1/awilso10/mod06/

##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt