# R BATCH --no-save < /home/grc/scripts/R/compare_CCTM_ensembles.R
  require(ncdf)
  require(maps)
 require(mapdata)
 require(date)
 require(fields)
 source("/home/grc/scripts/R/grid_projection_functions.R")
##################################################################################################
 # Control parameters
 doDataRead<-F
 doTseriesPlot<-T
 doO3ObsRead<-F
 
 convfac<-1000 # If AQ values are in ppm
 convfac<- 1		# If AQ values CO and O3 are in pbb
 
 
 
 COV<-F 
 # Ozone Obs File
 ozone_obs_file<-"/work/METEO/grc/CONUS12_twoway/RD_501_44201_2011-0.txt"
 
 # Arrays for trajectory or/and time series plots
 c_name<-c("Atlanta","Dallas","Baltimore","NYC","Boston","Charlotte","Chicago","Columbus","Houston","Los Angeles","Sacramento")
 t_lat<-c( 33.75,32.77 ,39.3  ,40.78 ,42.35 ,35.23 ,41.83 ,40.00 ,29.75 ,34.05, 38.58)
 t_lon<-c(-84.38,-96.77,-76.63,-73.97,-71.08,-80.83,-87.62,-83.02,-95.35,-118.25,-121.5)
 tindex_traj_start<-113
 tindex_traj_end  <-1
 ens_col<-c(2,2,2,2,2, 3,3,3,3,3, 4,4,4,4,4)
 ens_lwd<-c(2,1,0.75,1,0.75, 2,1,0.75,1,0.75, 2,1,0.75,1,0.75)
 ens_lwd<-c(1,0.75,0.75,0.75,0.75, 1,0.75,0.75,0.75,0.75, 1,0.75,0.75,0.75,0.75)
 # Model Level to analyze
 lev<-1

##################################################################################################
###############################################################
#- - - - - - - - -   START OF FUNCTION  -  - - - - - - - - - ##
###############################################################
# Simple format function, formats numbers from one to two digit (e.g. 7 to 07)
   dform<-function(num) {
     if(num<10){num<-paste("0",num,sep="")}
     return(num)
   }
#####--------------------------		END OF FUNCTION 	--------------------------------------####
##########################################################################################################
##################
# Ensemble arrays
 ens_id<-c("em.ctl","em.p1","em.p2","em.n1","em.n2",
           "nmm.ctl","nmm.p1","nmm.p2","nmm.n1","nmm.n2",
           "eta.ctl1","eta.ctl2","eta.p1","eta.p2","eta.n1","eta.n2")
 ens_id<-c("em.ctl","em.p1","em.p2","em.n1","em.n2",
           "nmm.ctl","nmm.p1","nmm.p2","nmm.n1","nmm.n2",
           "eta.ctl1","eta.ctl2","eta.p1","eta.n1","eta.n2")
 ens_base_dir<-paste("/work/METEO/grc/CONUS12_twoway/output_us_doe_nf_rrtmg_6_5_1_v3350-",ens_id,sep="")
 
 ens_conc_files<-paste(ens_base_dir,"/CCTM_CONC.20110605",sep="")
 ens_wrf_files<-paste(ens_base_dir,"/wrfout_d01_2011-06-05_22:00:00",sep="")

 ens_id<-c("em.ctl","em.p1","em.p2","em.n1","em.n2",
           "nmm.ctl","nmm.p1","nmm.p2","nmm.n1","nmm.n2",
           "eta.ctl1","eta.p1","eta.n1","eta.n2")
 ens_base_dir<-"/work/METEO/grc/CONUS12_twoway/CCTM_WRF_ENS_CTL"
 
 ens_conc_files<-paste(ens_base_dir,"/CCTM_CONC_SMALL.20110605.",ens_id,sep="")
 ens_wrf_files<-paste(ens_base_dir,"/wrfout_d01_2011-06-05_21:00:00.",ens_id,sep="")
##################

 nf<-length(ens_conc_files)


##################################################################
# Start Data read and cosolidation 
if(doDataRead){

 for(f in 1:nf) {
  fileC<-ens_conc_files[f]
  fileW<-ens_wrf_files[f]
  writeLines(paste("CONC file:",fileC))
  writeLines(paste("WRF file:",fileW))

  f1  <-open.ncdf(fileC)
  f2  <-open.ncdf(fileW)

    o3 <- get.var.ncdf(f1,varid="O3", start=c(1,1,lev,1), count=c(-1,-1,1,-1))
    co <- get.var.ncdf(f1,varid="CO", start=c(1,1,lev,1), count=c(-1,-1,1,-1))
    ntr <- get.var.ncdf(f1,varid="NO", start=c(1,1,lev,1), count=c(-1,-1,1,-1))

    t2 <- get.var.ncdf(f2,varid="T2")
    ws <- sqrt( get.var.ncdf(f2,varid="U10")^2 + get.var.ncdf(f2,varid="V10")^2 )
    pbl <- get.var.ncdf(f2,varid="PBLH")

    u <- get.var.ncdf(f2,varid="U", start=c(1,1,lev,1), count=c(-1,-1,1,-1))
    v <- get.var.ncdf(f2,varid="V", start=c(1,1,lev,1), count=c(-1,-1,1,-1))

    nxc<-dim(o3)[1]
    nyc<-dim(o3)[2]
    ntc<-dim(o3)[3]

    nxw<-dim(t2)[1]
    nyw<-dim(t2)[2]
    ntw<-dim(t2)[3]
  
    nxu<-dim(u)[1]
    nyu<-dim(u)[2]
    nxv<-dim(v)[1]
    nyv<-dim(v)[2]

    writeLines(paste("Number of times in file:",ntc))
    if(f ==1){      
      ntimesC    <-array(NA,c(nf))
      o3ens     <-array(NA,c(nxc,nyc,ntc,nf))
      coens     <-array(NA,c(nxc,nyc,ntc,nf))
      ntrens    <-array(NA,c(nxc,nyc,ntc,nf))

      ntimesW    <-array(NA,c(nf))
      t2ens    <-array(NA,c(nxw,nyw,ntw,nf))
      wsens    <-array(NA,c(nxw,nyw,ntw,nf))
      pblens   <-array(NA,c(nxw,nyw,ntw,nf))

      uens     <-array(NA,c(nxu,nyu,ntw,nf))
      vens     <-array(NA,c(nxv,nyv,ntw,nf))

      slat     <- get.var.ncdf(f2, varid="XLAT", start=c(1,1,1), count=c(-1,-1,1))
      slon     <- get.var.ncdf(f2,varid="XLONG", start=c(1,1,1), count=c(-1,-1,1))
      ulat     <- get.var.ncdf(f2, varid="XLAT_U", start=c(1,1,1), count=c(-1,-1,1))
      ulon     <- get.var.ncdf(f2,varid="XLONG_U", start=c(1,1,1), count=c(-1,-1,1))
      vlat     <- get.var.ncdf(f2, varid="XLAT_V", start=c(1,1,1), count=c(-1,-1,1))
      vlon     <- get.var.ncdf(f2,varid="XLONG_V", start=c(1,1,1), count=c(-1,-1,1))
      
    }
    
    ntimesC[f]<-ntc
    ntimesW[f]<-ntw
    o3ens[,,1:ntc,f] <-o3 
    coens[,,1:ntc,f] <-co 
    ntrens[,,1:ntc,f]<-ntr 
    
    t2ens[,,1:ntw,f] <-t2 
    wsens[,,1:ntw,f] <-ws 
    pblens[,,1:ntw,f]<-pbl 
    uens[,,1:ntw,f]  <-u 
    vens[,,1:ntw,f]  <-v 

  close.ncdf(f1) 
  close.ncdf(f2) 
 
 }
}  # End Data read and cosolidation 
#########################################################################

#########################################################################
if(doTseriesPlot){


  for(c in 1:length(c_name)){
   
   id_tmp<-GPLL(slat,slon,t_lat[c],t_lon[c])
     
   i_indexW<-id_tmp$i
   j_indexW<-id_tmp$j
   i_indexC<-i_indexW -6
   j_indexC<-j_indexW -6
   
   o3_ts<-o3ens[i_indexC,j_indexC,,]
   co_ts<-coens[i_indexC,j_indexC,,]
   ntr_ts<-ntrens[i_indexC,j_indexC,,]

   t2_ts<-t2ens[i_indexW,j_indexW,,]
   ws_ts<-wsens[i_indexW,j_indexW,,]
   pbl_ts<-pblens[i_indexW,j_indexW,,]
   #########################################################
   # O3
   figure<-paste("O3.",c_name[c],sep="")
   variable<-o3_ts*convfac
   pdf(file= paste(figure,".pdf",sep=""), width = 10, height = 4)
   par(new=F)
   plot(variable[,1],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[1], lwd=ens_lwd[1])
   for(e in 2:nf){
    par(new=T)
    plot(variable[,e],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[e], lwd=ens_lwd[e])
   }
   dev.off()
   #########################################################
   #########################################################
   # CO
   figure<-paste("CO.",c_name[c],sep="")
   variable<-co_ts*convfac
   pdf(file= paste(figure,".pdf",sep=""), width = 10, height = 3)
   par(new=F)
   plot(variable[,1],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[1], lwd=ens_lwd[1])
   for(e in 2:nf){
    par(new=T)
    plot(variable[,e],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[e], lwd=ens_lwd[e])
   }
   dev.off()
   #########################################################

   #########################################################
   # T2
   figure<-paste("T2.",c_name[c],sep="")
   variable<-t2_ts
   pdf(file= paste(figure,".pdf",sep=""), width = 10, height = 3)
   par(new=F)
   plot(variable[,1],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[1], lwd=ens_lwd[1])
   for(e in 2:nf){
    par(new=T)
    plot(variable[,e],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[e], lwd=ens_lwd[e])
   }
   dev.off()
   #########################################################
   #########################################################
   # WS
   figure<-paste("WS.",c_name[c],sep="")
   variable<-ws_ts
   pdf(file= paste(figure,".pdf",sep=""), width = 10, height = 3)
   par(new=F)
   plot(variable[,1],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[1], lwd=ens_lwd[1])
   for(e in 2:nf){
    par(new=T)
    plot(variable[,e],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[e], lwd=ens_lwd[e])
   }
   dev.off()
   #########################################################
   #########################################################
   # PBL 
   figure<-paste("PBL.",c_name[c],sep="")
   variable<-pbl_ts
   pdf(file= paste(figure,".pdf",sep=""), width = 10, height = 5)
   par(new=F)
   plot(variable[,1],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[1], lwd=ens_lwd[1])
   for(e in 2:nf){
    par(new=T)
    plot(variable[,e],ylim=range(variable,na.rm=T),ylab=figure,xlab="Hour Index",type='l',col=ens_col[e], lwd=ens_lwd[e])
   }
   dev.off()
   #########################################################

  }
}
#########################################################################
##################################################################
# Read Ozone Observations and set up obs arrays
if(doO3ObsRead){
a<-read.delim(file=ozone_obs_file,sep="|", header=T)

#dallas 48-113   69,45,87,55
sget<-48; coget<-113; siget<-c(69,45,87,55)

#baltimore
sget<-24; coget<-510; siget<-c(54)

#atlanta
sget<-13; coget<-121; siget<-c(55)

dayget<-c(20110606,20110607,20110608,20110609)



state_code<-a[,3]
co_code   <-a[,4]
site_id   <-a[,5]
date_obs  <-a[,11]
time_obs  <-a[,12]
o3_obs    <-a[,13]

obs_array<-array(NA,c(24,length(dayget),length(siget)))

for(d in 1:length(dayget)){
 for(s in 1:length(siget)){
  for(h in 1:24){
    hh_str<-paste(dform(h-1),":00",sep="")
    index_site<-which(state_code == sget & co_code == coget & 
                      site_id == siget[s] & date_obs == dayget[d] & time_obs == hh_str)
    writeLines(paste("site=",siget[s],"   Day:",dayget[d],"  Hour=",hh_str,sep=""))                  
    if(length(index_site) != 1) { next }
    obs_array[h,d,s]<- o3_obs[index_site] *convfac        

  }
 }
}

 
 

}