#obs.file <- "/work/MOD3APP/ezv/data/AQMEII/Harvard_Forest/hourly_harvard_forest_aqmeii_filled.csv"; sline <- 2
#obs.dat2 <- read.csv(obs.file)
#obs.file <- "/work/MOD3APP/ezv/data/AQMEII/Harvard_Forest/hourly_harvard_forest_aqmeii.csv"; sline <- 2
#obs.file <- "/work/MOD3APP/ezv/data/AQMEII/Harvard_Forest/hourly_harvard_forest_aqmeii_12022020.csv"; sline <- 2
obs.file <- "C:/Research and study 2021 spring/Updated R for drydep/for_Bin/ESTAR case update/hourly_harvard_forest_aqmeii_05242021.csv"; sline <- 2
obs.dat <- read.csv(obs.file)
header1 <- colnames(obs.dat)                  # Column names for each column
header2 <- as.character(obs.dat[1,])          # Units for each variable
obs.dat <- read.csv(obs.file, skip = 2, header = FALSE) # read observation data without variable names
dim(obs.dat) # dimension of the observation data
names(obs.dat) <- header1 # assign the column names to the observation data
#head(obs.dat)

###read PBLH & PSFC data from input file
PBLH.file <- "C:/Research and study 2021 spring/Updated R for drydep/for_Bin/ESTAR case/PBLHNC01_C01.csv"
PBLH.dat <- read.csv(PBLH.file, header = TRUE)
#head(PBLH.dat)
PSFC.file <- "C:/Research and study 2021 spring/Updated R for drydep/for_Bin/ESTAR case/PSFCNC01_C01.csv"
PSFC.dat <- read.csv(PSFC.file, header = TRUE)
#head(PSFC.dat)

###Create a dummy variable Estar and set it to zero
Estar <- matrix(data = 0, nrow = dim(obs.dat)[1], ncol = 1, byrow = FALSE, dimnames = NULL)

# CMAQ model parameters
Name      <- c("O3")
N.s       <- length(Name)
z1       <- 10 # CMAQ mid layer 1 usually 10
zsoil    <- 0.1 # depth of soil moisture measurement
Rst_min   <- 200 # PX hardwoods
if(version == 'Opt') Rst_min   <- 200 # Optimized
#Rstmin   <- 40 # NOAH

flux.fact <- 1

vd_o3 <- as.numeric(as.character(obs.dat$ozone_deposition_velocity)) # units in meta data are incorrect

# Land use/site specific variables 
# STAGE parameters are listed in run.STAGE.r
#l_leaf   <- 0.05 # Massad et al 2010
l_leaf   <- 0.039*0.8 + 0.001*0.2# Mean of TRY observations
l_life   <- 5.5 # months median of try observations
#lma      <- 65*0.8 + 154*0.2# leaf mass per area g/m2
ll       <- 5.5 # leaf lifetime mon

a_cut    <- 0.176 
snow <- 0 # fraction snow coverage
veg <- 1  # fraction veg coverage
##** Soil data **##
# Sandy Loam (https://esajournals.onlinelibrary.wiley.com/doi/epdf/10.1890/ES13.00183.1) from the WRF v4.0 PX tables
wsat <- 0.435
wfc  <- 0.195
wwlt <- 0.114
bslp <- 4.9 
wres <- 0.041
bulk_d <- 1.6
pH     <- 6.0

# Time and date
obs.date <- paste(obs.dat$month,'/',obs.dat$day,'/',obs.dat$year,sep="")
obs.time <- as.numeric(as.character(obs.dat$hour))

#l_wet    <- as.numeric(as.character(obs.dat2$leaf_wetness_filled))
l_wet    <- as.numeric(as.character(obs.dat$leaf_wetness_filled))

# General Met
rain     <- as.numeric(as.character(obs.dat$precipitation))
#l_wet <- rain
#l_wet[rain>0] <- 1
RH       <- as.numeric(as.character(obs.dat$relative_humidity))*100
Ta       <- as.numeric(as.character(obs.dat$air_temperature))+273.15


# Soil parameters
soilm    <- as.numeric(as.character(obs.dat$soil_moisture))
# ratio of available soil moisture to m3/m3
soilm[soilm>wsat] <- wsat
Tg       <- as.numeric(as.character(obs.dat$soil_temperature))+273.15

# Canopy data
Ht   <- as.numeric(as.character(obs.dat$canopy_height))
z0   <- 0.1*Ht
lai  <- as.numeric(as.character(obs.dat$leaf_area_index))
Ts       <- Ta
#snow <- as.numeric(as.character(obs.dat$X.snow_cover))
#L_mo <- obs.dat$MO

# energy balance
# convert from mmol/(m2 s) to W/m2
L <- as.numeric(as.character(obs.dat$h2oflux))*18*(2.501-0.00237*(Ta-273.15))
H <- as.numeric(as.character(obs.dat$sensible_heat_flux))

# Momentum   
Ubar     <- as.numeric(as.character(obs.dat$wind_speed))
ustar    <- as.numeric(as.character(obs.dat$friction_velocity))

# Radiation budget 
Rs       <- as.numeric(as.character(obs.dat$shortwave_radiation))

# Concentrations
O3    <- as.numeric(as.character(obs.dat$ozone_mixing_ratio))

# Fluxes units of ozone flux are a factor of 1e-6 lower than O3*vd
O3flux  <- as.numeric(as.character(obs.dat$ozone_flux))#*1e6 

PBLH <- as.numeric(as.character(PBLH.dat$PBLH.data.for.Harvard.forest.site.from.wrf.36.km.grids))
PSFC <- as.numeric(as.character(PSFC.dat$PSFC.data.for.Harvard.forest.site.from.wrf.36.km.grids))

# Time and date
date.time <- strptime(paste(obs.date,obs.time),"%m/%d/%Y %H")
summary(date.time)

### remove missing or erroneous values ###
qa_flag       <- as.numeric(as.character(obs.dat$AQMEII_Flag))

# find missing data for vital variabls 
val.miss <- complete.cases(Ubar,ustar,RH,Rs,Ta,O3flux,soilm,Rs,Ts,lai,H,PBLH,PSFC)
#val.miss <- complete.cases(Ubar,ustar,RH,Rs,Ta,O3flux,soilm,Rs,Ts,lai)
# valid data
#valid <- which( RH <= 100 & val.miss == T & soilm > 0 & Ubar > 0 & ustar > 0 & O3 > 0 & l_wet >= 0 &
#                 O3flux >= -2 & O3flux < 0 & Ubar>ustar & Tg > 0 & Tg < 315 & Ta > 0 & Ta < 320 & RH > 0 & lai > 0 )
#valid <- which( RH <= 100 & val.miss == T & soilm > 0 & Ubar > 0 & ustar > 0 & O3 > 0 & l_wet >= 0 &
#                Ubar>ustar & Tg > 0 & Tg < 315 & Ta > 0 & Ta < 320 & RH > 0 & lai > 0  & O3flux > -1)
#valid <- which( val.miss == T & soilm > 0 & Ubar > 0 & ustar > 0 & O3 > 0 & l_wet >= 0 & O3flux < 0.11 & O3flux > -999 & abs(vd_o3) < ustar^2/Ubar &
#                Tg > 0 & Ta > 0 & RH > 0 & lai >= 0)
#valid <- which( val.miss == T & qa_flag ==1 & lai >= 0 & Ta > 0 & Ts > 0 & Tg > 0 & Ubar > 0 & ustar > 0 & RH > 0)
valid <- which( val.miss == T & qa_flag ==1 & soilm > 0 & Ubar > 0 & ustar > 0 & O3 > 0 & l_wet >= 0 & O3flux < 0.11 & O3flux > -999 & abs(vd_o3) < ustar^2/Ubar &
                  Tg > 0 & Ts > 0 & Ta > 0 & RH > 0 & lai >= 0 & H > -999.99 & Rs > -999.99 & !(is.na(date.time)))
PBLH <- PBLH[valid]
PSFC <- PSFC[valid]
# number of time steps
N.t <- length(valid)


# Time and date
date.time <- strptime(paste(obs.date,obs.time),"%m/%d/%Y %H")[valid]

# vegetation parameters
l_wet    <- l_wet[valid]
s_wet    <- l_wet
s_wet[soilm[valid]==wsat] <- 1
# General Met
rain     <- rain[valid]
RH       <- RH[valid]
Ta       <- Ta[valid]

# Soil parameters
soilm <- soilm[valid]
Tg    <- Tg[valid]

# Canopy parameters
Ht   <- Ht[valid]
z0   <- z0[valid]
lai  <- lai[valid]
#lai  <- lai*1.2
#mlai  <- runMean(lai,168)
#val.lai <- which(is.na(mlai)==F)
#lai[val.lai] <- mlai[val.lai]
Ts    <- Ts[valid]
#snow <- snow[valid]


# Energy balance
H <- H[valid]
L <- L[valid]

# Momentum   
Ubar     <- Ubar[valid]
ustar    <- ustar[valid]
Rs       <- Rs[valid]
ustarB   <- ustar
# Concentrations
O3    <- O3[valid]

# Fluxes
O3flux  <- O3flux[valid]
vd_o3   <- vd_o3[valid]

# allocate the arrays in the same order as the model species names listed above
# ambient concentration array 
concs <- array(0,c(N.s,N.t))
concs[1,] <- O3

# Air-surface flux arrays
F.obs <- array(0,c(N.s,N.t))
F.obs[1,] <- O3flux

# assume decidous 65 g/m3 and Eastern Hemlock 154 g/m3
lma <- 0.2/(lai/max(lai))*154+0.8*(lai/max(lai))*65
lma[lma>154] <- 154 

#create new ustar using wind speed data, 
#harvard site parameters: z = 29 meters, z0 = 2.4 meters, zd =16.9, ustar = k*Ubar/ln((z-zd)/z0)
#ustar <- 0.1632*Ubar for just ln(z/zo)
#zinst <- 29.0
#zcanopy <- 24.0
#znot <- 0.1*zcanopy
#zdispht <- 16.9
#0.41/log((29-16.9)/2.4)
#log(2.71828)
ustar <- 0.25344*Ubar
#length(ustar)
ustarN <- ustar

est <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
wst <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
Wf <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
Ff2 <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
Lo  <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
p   <- matrix(data = 0, nrow = length(ustar), ncol = 3, byrow = FALSE, dimnames = NULL)

###Calculate Estar values
for (i in 1:length(ustar)) {
  p1 <- 3.8*ustar[i]*ustar[i]
  Theta <- (1000./PSFC[i])**(287./1004.)
  Theta <- Theta * Ta[i]
  if(H[i] > 0.) {
    wst[i] <- ((9.8/Theta)*PBLH[i]*H[i]/1004.0)**0.33333 
  } else {
    wst[i] = 0.0
  }  
  amol <- (-1.*ustar[i]**3)*Ta[i]
  amol <- amol/(0.4*9.8)
  temp <- H[i]/1004.     # new line
  amol <- amol/temp      # new line
  if(H[i] > 0.) {   
    p2 <- 0.2*wst[i]*wst[i]
    p3 <- 1.9*ustar[i]*ustar[i]
    p4 <- 29.
    p3 <- p3*(-p4/amol)**0.6666
    est[i] <- sqrt(p1+p2+p3)
  } else {
    est[i] <- sqrt(p1)
    p2 <- 0
    p3 <- 0
  }
  p[i,1] <- p1
  p[i,2] <- p2
  p[i,3] <- p3
  if(amol < 0 ) {
  Ff2[i] <- (p2+p3)/p1
  }else {
  Ff2[i] <- 0
  }
  
  if(H[i] > 10) {
  Wf[i] <- 0.24*(wst[i]/ustar[i])^2
  }else {
  Wf[i] <- 0
  }
}

Tf  <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
Tf1 <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
Tf2 <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
Tf3 <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)

###Calculate Tf values
for (i in 1:length(ustar)) {
  Theta <- (1000./PSFC[i])**(287./1004.)
  Theta <- Theta * Ta[i]
  if(H[i] > 0.) {
    wst[i] <- ((9.8/Theta)*PBLH[i]*H[i]/1004.0)**0.33333 
  } else {
    wst[i] = 0.0
  }  
  T1 <- 2.28*ustar[i]*ustar[i]
  amol <- (-1.*ustar[i]**3)*Ta[i]
  amol <- amol/(0.4*9.8)
  temp <- H[i]/1004.     
  amol <- amol/temp      
  if(H[i] > 10) {   
    T2 <- 5.32*ustar[i]*ustar[i] + 0.44*wst[i]*wst[i]
    T3 <- 3.8*ustar[i]*ustar[i]
    T4 <- 29.
    T3 <- T3*(-T4/amol)**0.6666
  } else {
    T2 <- 0
    T3 <- 0
  }
  Tf1[i] <- T1
  Tf2[i] <- T2
  Tf3[i] <- T3
#  if(H[i] > 10 ) {
#    Tf[i] <- (T2+T3)/T1
#  }else {
#    Tf[i] <- 0
#  }
  
  if(H[i] > 10 ) {
    Tf[i] <- (T2+T3)/T1
  }else if(H[i] < 10 && H[i] > -10){
    Tf[i] <- 2.33
  }else {
    Tf[i] <- 2.33
  }
}


VonK <- matrix(data = 0, nrow = length(ustar), ncol = 1, byrow = FALSE, dimnames = NULL)
for (i in 1:length(ustar)) {
  if(H[i] > 0.) {
    VonK[i] <- 0.41^2
  }
  else {
    VonK[i] <- 0.41^2
  }
}
head(VonK)
head(est)
tail(est)
ustar <- est
df = data.frame(date.time, lai)
write.csv(x=df, file="C:/Research and study 2021 spring/Updated R for drydep/for_Bin/LAI analysis/LAI.csv")

###check ustar, estar, and wstar values
#df = data.frame(date.time, ustarB, est, Ubar)


#df = data.frame(date.time, ustarB, est, H)
#write.csv(x=df, file="C:/Research and study 2021 spring/Updated R for drydep/for_Bin/Stratified ESTAR/Stratified ESTAR.csv")

###check Tf, Tf1, Tf2, and Tf3 values
df = data.frame(date.time, Tf, Tf1, Tf2, Tf3, H)
head(df)
write.csv(x=df, file="C:/Research and study 2021 spring/Updated R for drydep/for_Bin/Tf validation/Conversion factor condition.csv")

###check Ff1, Ff2, and Lo values
#df = data.frame(date.time, Ff1, Ff2, Lo, H)
df = data.frame(date.time, Wf)
head(df)
write.csv(x=df, file="C:/Research and study 2021 spring/Updated R for drydep/for_Bin/Tf validation/Wf.csv")
getwd()
#summary(H)
#length(H[H==0])