setwd("C:/Research and study 2021 spring/Updated R for drydep/for_Bin/ESTAR case update")
#----------------------------------------------------------------------------------------------------------------------------------------
# Resistance functions
#----------------------------------------------------------------------------------------------------------------------------------------
# Aerodynamic resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.ra <- function(ustar, Ubar, VonK) {
# inputs:
# surface friction velocity (ustar) m/s
# mean wind speed (Ubar) m/s
# If ustar is estimated from Ubar without the diabatic corection factor, the correction term will have to mbe added to the 
# expression below. 
   pr0  <- 1.0
   Ra     <- pr0*Ubar/(VonK*(ustar^2))
# P-X
#   pr0  <- 0.95
#   betah  <- 8.21
#   gamah  <- 11.6
#   zntol    <- z0/L_mo
#   zol      <- z1/L_mo
#   psih2    <- ustar*0

#   psih2[L_mo<0] <- 2*log((sqrt(1-gamah*zol[L_mo<0])+1)/(sqrt(1-gamah*zntol[L_mo<0])+1))
#   psih2[L_mo>0 & zol-zntol < 1] <- -betah*(zol-zntol)[L_mo>0 & zol-zntol < 1]   
#   psih2[L_mo>0 & zol-zntol > 1] <- 1-betah-(zol-zntol)[L_mo>0 & zol-zntol > 1]

#   gzozo  <- log((z1-0.7*Ht)/z0)

#   Ra     <- pr0*(gzozo-psih2)/(0.4*ustar)
   return(Ra)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Leaf quasi laminar boundary layer resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.rb_Massad <- function(Name,ustar, lai, Ta, l_leaf){
# Species Name (Name)
# leaf area index (lai) m2/m2
# surface friction velocity (ustar) m/s
# air temperature (Ta) K
# molar volume (MV) l/mol
# molecular weight (MW) g/mol
# characteristic leaf width from Massad et al. 2010 Table 6
   lai[lai<1] <- 1
   k_vis   <- get.k_vis(Ta)
   dif_T   <- get.dif_T(Name,Ta)
# LE error of 84% at Duke forest and 26.6% at Lillington
   karman  <- 0.4
   Rb_leaf <- k_vis / ( dif_T * karman * ustar ) * ( l_leaf * karman * ustar / (k_vis * lai^2 ) )^(1.0/3.0)
#  Corrected Jenson and Hummelshoj 1995/1997
# LE error of 81% at Duke forest and 26.2% at Lillington
#   Rb_leaf <- k_vis / ( dif_T * ustar * lai ) * ( 100 * l_leaf * ustar / (k_vis * lai^2 ) )^(1.0/3.0)
# integrated through the canopy assuming a uniform LAI and Rb(z) = Sc*Re^0.5/u*(z) using Yi 2008 in canopy model
#  LE error of 80% at Duke forest and 26.4% at Lillington
# Improves resutls for all AQMEII 4 O3 flux data.
#   Rb_leaf <- sqrt(l_leaf*k_vis / ustar)/dif_T*4/lai^2*(exp(lai/4)-1)/lai
# Clifton et al. 2020
# LE error 81.8% at Duke forest and of 26.4% at Lillington
#   a <- 3
#   b <- 0.02
#   Rb_leaf <- a/b * sqrt(l_leaf/Ubar)*1/(1-exp(-a/2))*(k_vis/dif_T/0.71)^(2/3)
   return(Rb_leaf)
}
get.rb_JH <- function(Name,ustar, lai, Ta, l_leaf){
# Species Name (Name)
# leaf area index (lai) m2/m2
# surface friction velocity (ustar) m/s
# air temperature (Ta) K
# molar volume (MV) l/mol
# molecular weight (MW) g/mol
# characteristic leaf width from Massad et al. 2010 Table 6
   lai[lai<1] <- 1
   k_vis   <- get.k_vis(Ta)
   dif_T   <- get.dif_T(Name,Ta)

#  Corrected Jenson and Hummelshoj 1995/1997
   Rb_leaf <- k_vis / ( dif_T * ustar * lai ) * ( 100 * l_leaf * ustar / (k_vis * lai^2 ) )^(1.0/3.0)
# integrated through the canopy assuming a uniform LAI and Rb(z) = Sc*Re^0.5/u*(z) using Yi 2008 in canopy model
#  LE error of 80% at Duke forest and 26.4% at Lillington
# Improves resutls for all AQMEII 4 O3 flux data.
#   Rb_leaf <- sqrt(l_leaf*k_vis / ustar)/dif_T*4/lai^2*(exp(lai/4)-1)/lai
# Clifton et al. 2020
# LE error 81.8% at Duke forest and of 26.4% at Lillington
#   a <- 3
#   b <- 0.02
#   Rb_leaf <- a/b * sqrt(l_leaf/Ubar)*1/(1-exp(-a/2))*(k_vis/dif_T/0.71)^(2/3)
   return(Rb_leaf)
}
get.rb_Bash <- function(Name,ustar, lai, Ta, l_leaf){
# Species Name (Name)
# leaf area index (lai) m2/m2
# surface friction velocity (ustar) m/s
# air temperature (Ta) K
# molar volume (MV) l/mol
# molecular weight (MW) g/mol
# characteristic leaf width from Massad et al. 2010 Table 6
   lai[lai<1] <- 1
   k_vis   <- get.k_vis(Ta)
   dif_T   <- get.dif_T(Name,Ta)
# integrated through the canopy assuming a uniform LAI and Rb(z)*lai = Sc*Re^0.5/u*(z) using Yi 2008 in canopy model
# Improves resutls for all AQMEII 4 O3 flux data.
#   Rb_leaf <- sqrt(l_leaf*k_vis / ustar)/dif_T*4/lai^2*(1-exp(-lai/4))#*4/lai^2*(exp(lai/4)-1)
   Ust_mean <- 2*ustar/lai*(1-exp(-lai/2)) # mean canopy ustar
   Re_f     <- Ust_mean*l_leaf/k_vis      # Friction Reynold's number
   nexp     <- 0.65/(1+exp(0.025*(Re_f-45)))+0.33 # Imperical exponent 1 = diffusive, 0.24 = advection influences
   Rb_leaf  <- k_vis/(lai*dif_T*Ust_mean)*Re_f^nexp 
   return(Rb_leaf)
}

#----------------------------------------------------------------------------------------------------------------------------------------
# Soil quasi laminar boundary layer resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.rb_soil_Bash <- function(Name,ustar,Ta, lai){
# Species Name (Name)
# surface friction velocity (ustar) m/s
# air temperature (Ta) K
#   karman  <- 0.4
# Calculate Soil Resistance Nemitz et al 2000 https://doi.org/10.1016/S0168-1923(00)00206-9
   k_vis   <- get.k_vis(Ta)
   dif_T   <- get.dif_T(Name,Ta)
#   scn     <- k_vis / dif_T
# ustar at the soil surface following Yi 2008 https://doi.org/10.1175/2007JAMC1667.1
   ustg    <- ustar * exp( -lai / 2.0 )
#   ustg[ustg<0.001] <- 0.001         
# Bash et al 2010
#   ustg    <- ustar * exp( -lai )
#   ustg[ustg<0.001] <- 0.001         
#   del0    <- dif_T / ( karman * ustg )
# LE error of 81% at Duke forest and 26.2% at Lillington
#   Rb_soil <- ( scn - log( del0 / 0.10 ) ) / ( karman * ustg )

# integrated through the canopy assuming a uniform LAI and Rb(z)*lai = Sc*Re^0.5/u*(z) using Yi 2008 in canopy model
# Improves resutls for all AQMEII 4 O3 flux data.
#   Rb_leaf <- sqrt(l_leaf*k_vis / ustar)/dif_T*4/lai^2*(1-exp(-lai/4))#*4/lai^2*(exp(lai/4)-1)
#   Ust_mean <- 2*ustar/lai*(1-exp(-lai/2)) # mean canopy ustar
   Re_f     <- ustg*0.05/k_vis      # Friction Reynold's number
   nexp     <- 0.65/(1+exp(0.025*(Re_f-45)))+0.33 # Imperical exponent 1 = diffusive, 0.24 = advection influences
   Rb_soil  <- k_vis/(lai*dif_T*ustg)*Re_f^nexp 

# integrating to get z1
# LE error of 76.4% at Duke forest and 24.7% at Lillington
# z1 is the momentum thickness and del0/z1 is the boundary layer shape factor
# from Blasius solution
#   z1 <- 0.664*sqrt(k_vis*Hc*0.07/ustar)
#   z1[lai<1] <- exp(-1)
#   Rb_soil <- ( scn - log( del0 / z1 ) ) / ( karman * ustg )
# Clifton et al. 2020
# LE error of 99.4% at Duke forest and 29.3% at Lillington
#   Rb_soil <- 2 / ( karman * ustg ) * scn^(2/3)
# using z0 soil as the length and Jenson and Hummelshoj 1995 model
#   zsoil <-  (mean(sand)*1.025+mean(silt)*0.0251+mean(clay)*0.001)/100/1e3
#   Rb_soil <- sqrt(0.005*k_vis/ustg)/dif_T
   return(Rb_soil)
}
get.rb_soil <- function(Name,ustar,Ta, lai){
# Species Name (Name)
# surface friction velocity (ustar) m/s
# air temperature (Ta) K
   karman  <- 0.4
# Calculate Soil Resistance Nemitz et al 2000 https://doi.org/10.1016/S0168-1923(00)00206-9
   k_vis   <- get.k_vis(Ta)
   dif_T   <- get.dif_T(Name,Ta)
   scn     <- k_vis / dif_T
# ustar at the soil surface following Yi 2008 https://doi.org/10.1175/2007JAMC1667.1
   ustg    <- karman * ustar * exp( -lai / 2.0 )
   ustg[ustg<0.001] <- 0.001         
# Bash et al 2010
#   ustg    <- ustar * exp( -lai )
#   ustg[ustg<0.001] <- 0.001         
   del0    <- dif_T / ( karman * ustg )
# LE error of 81% at Duke forest and 26.2% at Lillington
   Rb_soil <- ( scn - log( del0 / 0.10 ) ) / ( karman * ustg )

   return(Rb_soil)
}

#----------------------------------------------------------------------------------------------------------------------------------------
# Stomatal resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.r_stom <- function(Name,soilm, wwlt, wfc, Ta, Ts, Rs, Rst_min, lai, ustar, Ubar, RH) {
# inputs:
# Species Name (Name)
# Volumetric soil moisture (soilm) m3/m3
# Soil field capacity at -33 kPa (wfc) m3/m3
# Soil wilting point at -1500 kPa (wwlt) m3/m3
# air temperature (Ta) K
# surface temperature (Ts) K
# Solar radiation (Rs) W/m2
# Minimum stomatal resistance (Rst_min) s/m
# leaf area index (lai) m2/m2
# surface friction velocity (ustar) m/s
# Mean wind speed (Ubar) m/s
# relative humidity (RH) %
# estimate the stomatal resistance based on WRF 4.0 P-X parameterization
   rsmax <- 5000       # m/s
   f3min  <- 0.25      # 
   ftmin <- 0.0000001  # m/s
   dwat  <- 0.2178     # [cm^2/s] at 273.15K
# saturation water vapor 
   svp1 <- 0.611
   svp2 <- 17.67
   svp3 <- 29.65
   s.dat <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   H.Name  <- as.character(s.dat$H)[ns]
   f0    <- s.dat$f0[ns]
   es    <- 10*svp1*exp(svp2*(Ta-273.15)/(Ta-svp3))
   qsm   <- es*0.622/(1000-es)
   qvs   <- 0.01*RH*qsm
   qgs   <- 1000*es*0.622/(1000-es)
   ess   <- 10*svp1*exp(svp2*(Ts-273.15)/(Ts-svp3))
   qss   <- ess*0.622/(1000-ess)

   Ra <- get.ra(ustar, Ubar, VonK)

   w2avail <- soilm*exp( (0.5-zsoil) * 9.81 )^(1.0/bslp)-wwlt
   w2mxav  <- wfc - wwlt
   f2      <- 1.0 / ( 1.0 + exp( -5.0 * ( w2avail / w2mxav - ( w2mxav / 3.0 + wwlt ) ) ) )  
# air temperature function
   f4 <- 1.0 / ( 1.0 + exp( 0.5 * (Ta - 314.0 ) ) )
   f4.low <-which( Ta <= 302.15 ) 
   f4[f4.low] <-  1.0 / ( 1.0 + exp( -0.41 * (Ta[f4.low] - 282.05 ) ) )

# Radiation 
   par <- 0.45*Rs*4.566

   if(length(Rst_min)==1){
      if(Rst_min > 130 ){
         f1max <- 1.0-0.02*lai
      } else {
         f1max <- 1.0-0.07*lai
      }
      f1     <- f1max * ( 1.0 - exp( -0.0017 * par ) )
      f1[f1<Rst_min / rsmax] <- Rst_min / rsmax
   } else {
      f1max <- Rst_min*0
      f1max[Rst_min > 130] <- 1.0-0.02*lai[Rst_min > 130]
      f1max[Rst_min <= 130] <- 1.0-0.07*lai[Rst_min <= 130]
      f1     <- f1max * ( 1.0 - exp( -0.0017 * par ) )
      f1[f1<Rst_min / rsmax] <- Rst_min[f1<Rst_min / rsmax] / rsmax
   }
   ftot   <- lai * f1 * f2 * f4
#   ftot   <- max( ftot,ftmin )
   ftot[ftot<ftmin]   <- ftmin 
   fshelt <- 1.0   # go back to NP89
   gs     <- ftot / ( Rst_min * fshelt )
   Rbw <- get.rb_leaf("H2O",ustar, lai, Ta, l_leaf)
   raw    <- Ra + Rbw
   ga     <- 1.0 / raw
#-- Compute humidity effect according to RH at leaf surf
   f3     <- 0.5 * ( gs - ga + sqrt( ga * ga + ga * gs * ( 4.0 * qvs / qss - 2.0 ) + gs * gs ) ) / gs
   f3[f3<f3min] <- f3min
   f3[f3>1] <- 1
   f3[f3<f3min] <- f3min
   f3[f3>1] <- 1
   Rst    <- 1.0 / ( gs * f3 )   
   Rst[Rst<Rst_min/lai] <- (Rst_min/lai)[Rst<Rst_min/lai]
# consider mesophyll resistance and estimate resistance for the specified pollutant
   dif_T   <- get.dif_T(Name,Ta)
   dwat_T  <- dwat*1.0e-4 * ( Ta/273.15 )^1.81 # Following Massman 1999
#   dwat_T  <- get.dif_T('H2O',Ta)
   if(Name=='H2O'){
      R_stom <- Rst
   } else {
      heff    <- hlconst(H.Name,Ts,6) 
      R_stom <- Rst * dwat_T/dif_T + 1.0/(heff/3000+100*f0)/lai
   }
   return(R_stom)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Cuticular resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.r_cut <- function(Name,Ta,Ts,RH,pH,l_wet, a_cut, snow, lai) {
# Species Name (Name)
# air temperature (Ta) K
# surface temperature (Ts) K
# relative humidity (RH) %
# surface water pH (pH)
# leaf wetness (l_wet) ratio (0-1)
# Land use specific cuticular resistance parameter from Massad et al. 2010
# Snow coverage (snow) ratio (0-1)
# leaf area index (lai) m2/m2
   a0      <- 8
   rcut0   <- 3000.0     # [s/m]
   rwm     <- 31.5       # Minimum NH3 cuticle resistance [s/m] from Massad et al. 2010
   s.dat   <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   rel.rx  <- s.dat$rel.rx[ns]

   if(Name == 'O3'){
# Canopy level wet resistence Rwet to ozone was found to be about 200 s/m on basis of Keysburg exp
# Using LAI(1-sided) of about 6.25 measured at Keysburg gives leaf level rwet about 1250 s/m
# Leaf level rwet estimated from Altimir et al 2006 gives about 1350 s/m                           
      R_wet   <- 1250.0
#      R_wet <- get.r_wet(Name,Ta,Ts,pH)
      rh_func <- ( RH - 70.0 )/30
      rh_func[rh_func<0] <- 0.0
#      R_dry   <- 1.0 / ( ( 1.0 -rh_func) / ( rcut0 * a0 / rel.rx )  + rh_func / R_wet )

#      R_wet <- get.r_wet(Name,Ta,Ts,pH)
#      R_wet   <- 600.0 # same as leaf litter
#      sl <- exp(-0.1*(Ts-278.13))
#      sl[sl<1] <- 1
#      R_dry <- rcut0 * a0 / rel.rx#*sl # v5.3
#      R_dry <- 2108 # optimized constant all
#      R_dry <- 33267/sqrt(Ts) # optimized al
      if(version == 'Opt'){
         R_dry <- lma_fact*lma^ll_exp/sqrt(Ts) # optimized
      } else {
         R_dry <- rcut0 * a0 / rel.rx         
      }
#      R_dry <- ll_fact*ll^ll_exp/sqrt(Ts)
#       R_dry <- log(lma)*ll_fact
#      R_dry <- 289*lma^0.5*sqrt(298/Ts)#*sl
#      R_dry <- 1.0 / ( ( 1.0 -rh_func) / ( R_dry )  + rh_func / R_wet )
   } else {
      if(Name == 'NH3'){
# Massad et al. 2010 Cuticular resistance
         R_dry <- rwm * exp( a_cut * ( 100.0 - RH ) )
      } else {
         R_dry <- rcut0 * a0 / rel.rx
      }
      R_wet <- get.r_wet(Name,Ta,Ts,pH)
   } 
#   rh_func <- ( RH - 70.0 )/30.0 
#   rh_func[rh_func<0] <- 0.0 
#   val_wet <- which(rh_func>l_wet)
#   l_wet[val_wet] <- rh_func[val_wet] 
   R_snow <- get.r_snow( Name, Ta, Ts, pH)
   R_cut  <- snow * R_snow + (1-snow) / ( lai * ( l_wet / R_wet + (1-l_wet)/R_dry ) )
   return(R_cut)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# In canopy + ground resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.r_grnd <- function(Name,ustar,Ubar,Ta,Ts,Tg,lai,s_wet,snow,soilm,wsat,wfc,wwlt,wres,bslp,pH,Ht){
# Species Name (Name)
# surface friction velocity (ustar) m/s
# mean wind speed (Ubar) m/s
# air temperature (Ta) K
# surface temperature (Ts) K
# leaf area index (lai) m2/m2
# surface wetness (s_wet) ratio (0-1)
# Snow coverage (snow) ratio (0-1)
# volumetric soil moisture content (soilm) m3/m3
# volumetric saturation soil moisture content (wsat) m3/m3
# volumetric field capacity (air dry at -33 kPa) soil moisture content (wfc) m3/m3
# volumetric wilting piont (air dry at -1500 kPa) soil moisture content (wwlt) m3/m3
# volumetric residual (air dry at -3.0e5 kPa) soil moisture content (wres) m3/m3
# Exponent on the soil moisture release function using the Campbell 1985 model (bslp) 
# soil water pH (pH)
   rg0  <- 1000.0     # [s/m]
   GRAV <- 9.80622    # [m/s2]
   a0   <- 8
   s.dat <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   rel.rx  <- s.dat$rel.rx[ns]
   dif_T <- get.dif_T(Name, Ta)
# using bidirectional exchange soil diffusive resistance
# The following resistance parameterization is derived from measurements with soil samples of 2 cm thick (Kondo et al 1990)
# https://doi.org/10.1175/1520-0450(1990)029<0385:APOEFB>2.0.CO;2 as discussed in Sakaguchi and Zeng 2009 JGR 
# https://doi.org/10.1029/2008JD010834 According to Swenson and Lawrence 2014 (https://doi.org/10.1002/2014JD022314) and the 
# references therin the dry layer thickness varies from 1 to 3 cm. 
   ldry_max <- 0.02 #zsoil
   w1cm <- soilm #* exp( (0.01-zsoil) * GRAV )^(1.0/bslp)
   s_wet[soilm>=wsat] <- 1
# From Sakaguchi and Zeng 2009 JGR Equation 10
   ldry     <- ldry_max * ( exp( ( 1.0 - w1cm / wsat )^5 ) - 1.0 ) / 1.718         
   if(obs.file == "Grignon_dataset.csv") ldry[which(doy>=73.5 & doy<74)]<-0
   dp       <- dif_T * wsat^2 * ( 1.0 - wres / wsat )^( 2.0 + 3.0 / bslp )         
   R_dif    <- ldry / dp

   if(Name=='NH3' | Name == 'H2O'){
#      dif_T <- get.dif_T(Name, Ta)
# using bidirectional exchange soil diffusive resistance
# The following resistance parameterization is derived from measurements with soil samples of 2 cm thick (Kondo et al 1990)
# https://doi.org/10.1175/1520-0450(1990)029<0385:APOEFB>2.0.CO;2 as discussed in Sakaguchi and Zeng 2009 JGR 
# https://doi.org/10.1029/2008JD010834 According to Swenson and Lawrence 2014 (https://doi.org/10.1002/2014JD022314) and the 
# references therin the dry layer thickness varies from 1 to 3 cm. 
#      ldry_max <- 0.02
# From Sakaguchi and Zeng 2009 JGR Equation 10
#      ldry     <- ldry_max * ( exp( ( 1.0 - soilm / wsat )^5 ) - 1.0 ) / 1.718         
#      if(obs.file == "Grignon_dataset.csv") ldry[which(doy>=73.5 & doy<74)]<-0
#      dp       <- dif_T * wsat^2 * ( 1.0 - wres / wsat )^( 2.0 + 3.0 / bslp )         
      R_dry    <- R_dif
# this is replaced with the source/sink term in bidirectional exchange
      R_wet    <- 0
   } 
   if(Name=='O3'){
# Following based on measurements Fares et al 2014 https://doi.org/10.1016/j.agrformet.2014.08.014 for sandy soil 
# forests at 10cm measured soil moisture and Fumagalli et al. 20016 https://doi.org/10.1016/j.agrformet.2016.07.011 for sandy loam soils
# Here an asymptotic function was applied to set lower and upper bounds in the resistance as repoerted by Fumagalli et al. 2016
      w10cm   <- soilm * exp( (0.1-zsoil) * GRAV )^(1.0/bslp)
      w10cm[w10cm>wsat]   <- wsat
      sm_func <- ( w10cm-wwlt ) / wfc
#      sm_func <- ( soilm-wwlt ) / wfc
      sm_func[sm_func<0] <-  0 
      R_wet <- get.r_wet(Name,Ta,Tg,pH)     # Surface wetness   
      R_dry   <- (250.0 + 2000.0 * atan( sm_func^bslp ) /pi ) # v5.3
#      l_dry  <- (wsat-soilm)^2*((wsat-soilm)/soilm)^(3/bslp)/((wsat-wwlt)^2*((wsat-wwlt)/wwlt)^(3/bslp))
#      l_dry[l_dry>1] <- 1 
      sl <- exp(-0.025*(Ts-278.13))
      sl[sl<1] <- 1
#      R_dry <- rg0 * a0 / rel.rx #* sl
#      R_dry <- Rs_opt #* sl                
#       R_dry  <- 487*sqrt(298.13/Tg)#*sl # optimized Ispra
      if(version == 'Opt' | version == 'v5.2'){
         R_dry  <- soil_fact/Ts^0.25 # optimized
      } else {
         R_dry   <- (250.0 + 2000.0 * atan( sm_func^bslp ) /pi ) # v5.3
      }
#       R_dry   <- 1/(ldry/(ldry_max*R_dry)+1/(R_dif+R_wet))

# new pathway with deposition to soil and soil water
#      R_dry <- 1/(0.5/600+0.5/1000)
#      R_dry <- 1/(1/600+1/(R_dif+R_wet)) # diffusion to soil water
   }
   if(Name != 'NH3' & Name != 'O3' & Name != 'H2O'){
      R_dry <- rg0 * a0 / rel.rx 
      R_wet <- get.r_wet(Name,Ta,Ts,pH)
# Simple scaling from wet to dry. This has the advantage of being a continuous function.
#      R_dry <- 1/(1/(R_dif+R_wet)+ 1/(rg0 * a0 / rel.rx))
   }
   if( Name != 'H2O'){
      R_snow  <- get.r_snow(Name, Ta, Ts, pH)
   } else {
      R_snow  <- 1e6
   }
   R_soil  <- 1/((1-snow) / R_dry + snow / R_snow)
   if(Name != 'NH3' & Name != 'H2O'){
      val.wet <- which(s_wet == 1 )  
      R_soil[val.wet]  <- 1/((1-snow) / (1/(s_wet/R_wet+(1-s_wet)/R_dry)) + snow / R_snow)[val.wet]
   }
   if(version=='Opt'){
      Rb_soil <- get.rb_soil_Bash(Name,ustar,Ta,lai)
   }else{
      Rb_soil <- get.rb_soil(Name,ustar,Ta,lai)
   }
   if(version=='v5.2'){
      R_inc   <- get.r_inc_JWE(ustar,Ht,lai)
   } else {
      R_inc   <- get.r_inc(ustar,Ubar,lai)
   }
# note that this is a big leaf parameterization. If is multilayered model were to be used R_inc would the be the resistance from 
# the lower layer to the soil surface.
   R_grnd  <- R_soil + Rb_soil + R_inc
   return(R_grnd)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Funcitons used in the resistance models
#----------------------------------------------------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------------------------------------------------
# Interum resistances
#---------------------------------------------------------------------------------------------------------------------------------------
# In canopy resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.r_inc <- function(ustar,Ubar,lai){
# inputs:
# surface friction velocity (ustar) m/s
# mean wind speed (Ubar) m/s
# leaf area index (lai) m2/m2
   Ra <- get.ra(ustar, Ubar, VonK)
# Calculate in canopy aerodynamic resistance based on the momentum attenuation coefficient derived 
# by Yi 2008 https://doi.org/10.1175/2007JAMC1667.1
   R_inc <- Ra * ( exp( lai / 2.0 ) - 1.0 )
# Bash et al 2010
#   R_inc <- Ra/2 * ( exp( lai ) - 1.0 )
   return(R_inc)
}
get.r_inc_JWE <- function(ustar,Ht,lai){
# inputs:
# surface friction velocity (ustar) m/s
# canopy height (Ht) m
# leaf area index (lai) m2/m2
   R_inc <- 14 * lai * Ht / ustar
# Bash et al 2010
#   R_inc <- Ra/2 * ( exp( lai ) - 1.0 )
   return(R_inc)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Wet surface resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.r_wet <- function(Name,Ta,Ts,pH) {
# estimates the resitance for gas to parition into the aqueous phase used for soil and cuticular resistances
# Species Name (Name)
# air temperature (Ta) K
# surface temperature (Ts) K
# surface water pH (pH)
   rsnow0  <- 10000.0    # to match the maximum snow deposition (0.01 cm/s) of D Helmig et al. 2007  
   a0      <- 8
   s.dat <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   MW      <- s.dat$MW[ns]
   MV      <- s.dat$MV[ns]
   M_ac    <- s.dat$M.ac[ns]
   H.Name  <- as.character(s.dat$H)[ns]
   rel.rx  <- s.dat$rel.rx[ns]
   heff    <- hlconst(H.Name,Ts,pH) * 0.08205 * Ts
   dif_T   <- get.dif_T(Name,Ta)   
   RGASUNIV <- 8.314510
   rad_wat  <- 1.9e-4

   rmsv     <- sqrt( 3.0 * RGASUNIV * 1.0e3 * Ts / MW ) 
   rawmt    <- rad_wat / dif_T  + 4.0 / ( rmsv  * M_ac )
   R_wet    <- rawmt + rawmt/( heff * rad_wat )
   
# assume a frozen surface. This function will be combined with get.r_snow in the future to remove some inconsistencies, e.g.
# the deposition to partially melted snow.
   T.ice <- which(Ts < 273.15)
   R_wet[T.ice] <- rsnow0 * a0 / rel.rx

   return(R_wet)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Snow resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.r_snow <- function(Name, Ta, Ts, pH) {

# Species Name (Name)
# air temperature (Ta) K
# surface temperature (Ts) K
# surface water pH (pH)

# Liquid snow fraction modeled as a system dominated by van der Walls forces following Dash et al. 1999 S. Rep. Prog. Phys. 
# with a maximum fraction of the disordered interface acting as an aqueous solution as 20% following Conklin et al 1993 with
# the negligible impact of the disordered interface depth of 2 nm following Roth et al 2004. The 2 nm depth was approximated 
# to be around 263 degrees Celsius interpolated from figure 3 in Huthwelker et al 2006 doi:10.1021/cr020506v

   rsnow0  <- 10000.0    # to match the maximum snow deposition (0.01 cm/s) of D Helmig et al. 2007  
   a0      <- 8
   rsndiff <- 10.0       # snow diffusivity fac taken from m3dry
   s.dat   <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   rel.rx  <- s.dat$rel.rx[ns]
   T.slush <- which(273.15-Ts < 0.0)
#   if(273.15-Ts > 0.0 ){
      melt_snow <- 0.025 / (273.15-Ts)**(1.0/3.0)
      melt_snow[melt_snow>0.2] <- 0.2
      melt_snow[melt_snow<0.01] <- 0.01
#      melt_snow <- min(melt_snow, 0.2)
#      melt_snow <- max(melt_snow, 0.01)
#   } else {
      melt_snow[T.slush] <- 0.2
#   }
   R_wet   <- get.r_wet(Name,Ta,Ts,pH)
   R_ice   <- rsnow0 * a0 / rel.rx
   R_slush <- rsndiff + R_wet
   R_snow <- 1/((1-melt_snow)/R_ice+melt_snow/R_slush)
   return(R_snow)
}

#----------------------------------------------------------------------------------------------------------------------------------------
# Physical chemical parameters
#----------------------------------------------------------------------------------------------------------------------------------------
# Kinimatic viscosity of air
#----------------------------------------------------------------------------------------------------------------------------------------
get.k_vis <- function(Ta){
# air temperature (Ta) K
   kvis  <- 0.132      # [cm^2 / s] at 273.15K
   k_vis <- kvis*1.0e-4 * ( Ta/273.15 )^1.81 # Following Massman 1999
   return(k_vis)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Diffusivity of gases
#----------------------------------------------------------------------------------------------------------------------------------------
get.dif_T <- function(Name, Ta){
# Species Name (Name)
# air temperature (Ta) K
# Following Fuller et al 1966. Here we use the LeBas molar volume which is similar to Fuller's diffusive volume
# Molar volume of air at STP [ L/mol ] Non MKS units 
   MOLVOL <- 22.41410
# mean molecular weight for dry air [ g/mol ]
# FSB: 78.06% N2, 21% O2, and 0.943% A on a mole 
# fraction basis ( Source : Hobbs, 1995) pp. 69-70
   MWAIR <- 28.9628
   s.dat <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   MW      <- s.dat$MW[ns]
   MV      <- s.dat$MV[ns]
   dif_T   = 1.0e-7*Ta^1.75 * sqrt( 1.0/MWAIR + 1.0/MW ) /  ( MOLVOL^(1.0/3.0) + MV^(1.0/3.0) )^2
# To match Massman et al. 1998
#   dif_T   = 7.568e-8*Ta^1.81 * sqrt( 1.0/MWAIR + 1.0/MW ) /  ( MOLVOL^(1.0/3.0) + MV^(1.0/3.0) )^2  
   return(dif_T)
}
get.comp <- function(Name,Ts,gamma.g,gamma.st,gamma.cut, pH, soilm, bulk_d){
# Species Name (Name)
# Surface temperature (Ta) K
# gamma.g soil emissions potential (NH4/H+ for NH3)
# gamma.st stomatal emissions potential (NH4/H+ for NH3)
# gamma.cut cuticular emissions potential (NH4/H+ for NH3)
# Soil solution pH (pH)
# volumetric soil moisture (soilm) m3/m3
# soil bulk density (bulk_d) kg/l
   if(Name=='NH3'){
# Estimate the soil NH4 that is in solution and available for evasion following Venteria et al. Sci. Rep. doi:10.1038/srep12153
# Point were soil solution NH4 equals half the maximum sorption capacity Venteria et al. Sci. Rep. doi:10.1038/srep12153. 
# Values estimated from Lillington NC soils.
      half_sol <- 345.0 # measured at Lillington
#      half_sol <- 85.5 # Venteria et al 
# Maximum NH4 soil sorption capacity Venteria et al. Sci. Rep. doi:10.1038/srep12153
      max_sorp <- 550.0 # measured at Lillington
#      max_sorp <- 867  # Venteria et al

      MHp     <- 10^(-pH) # H+ in mol/L
      cg      <- soilm/bulk_d # L water / kg soil 
      mNH4    <- gamma.g*MHp/bulk_d*14.0*1.0e3 # NH4+NH3 in mg N/L
      Ka      <- 5.68853e-10*exp(-6248.151*(1.0/Ts-1.0/298.15))
      coef_a  <- cg*(1.0+Ka/MHp)
      coef_b  <- max_sorp+half_sol*coef_a - mNH4
      coef_c  <- -half_sol*mNH4
      NH4_sol <- (-coef_b+sqrt(coef_b^2.0-4.0*coef_a*coef_c))/(2.0*coef_a)
      NH3_sol <- NH4_sol*Ka/MHp # mg N/l
      cnh3g   <- (NH4_sol+NH3_sol)*1.0e-3/14.0/MHp # gamma of NH3 + NH4 in solution
      a1      <- 161512.0/Ts*10.0^( -4507.11 / Ts )
      Xsoil   <- a1 * cnh3g*17*1e3*1e6 # mol/l -> ug/m3 (17 g/mol * 1e3 l/m3 * 1e6 ug/g)
      Xstom   <- a1 * gamma.st*17*1e3*1e6 # mol/l -> ug/m3 (17 g/mol * 1e3 l/m3 * 1e6 ug/g)
# not used now but could be added to explore NH4 in dew. 
      Xcut    <- a1 * gamma.cut*17*1e3*1e6 # mol/l -> ug/m3 (17 g/mol * 1e3 l/m3 * 1e6 ug/g)
   } else {
# This is a placeholder for other bidirectional species. A parameterization or constants can be placed here to estimate compensation 
# points for vegetation or soil. 
      Xstom = 0
      Xcut  = 0
      Xsoil = 0
   }
   comp <- list('Stom'=Xstom,'Cut'=Xcut,'Soil'=Xsoil)
   return(comp)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Henery's constant
#----------------------------------------------------------------------------------------------------------------------------------------
hlconst <- function(Spc.names,TEMP,pH){
# get henry's law constants and enthalpy 
      h.dat <- read.csv("H.dat.csv")
      Spc   <- as.character(h.dat$Name)
      A_var <- h.dat$A
      E_var <- h.dat$E

# create the dummy output array
      KH <- array(0,c(length(Spc.names),length(TEMP)))

# array indecies
      LSO2       <-  1  # SO2
      LHSO3      <-  2  # HSO3
      LHNO2      <-  3  # HNO3
      LHNO3      <-  4  # HNO3
      LCO2       <-  5  # CO2
      LHCO3      <-  6  # HCO3
      LH2O2      <-  7  # H2O2
      LHCOOH     <-  8  # HCOOH
      LHO2       <-  9  # HO2
      LNH4OH     <- 10  # NH4OH
      LH2O       <- 11  # H2O
      LATRA      <- 12  # Atrazine
      LCL2       <- 13  # CL2
      LHOCL      <- 14  # HOCL
      LHCL       <- 15  # HCL
      LHYDRAZINE <- 16  # Hydrazine
      LHBR       <- 17  # HBR
      LHI        <- 18  # HI
      LACRYACID  <- 19  # CCH3COOH
# dissociation constants 
      diss_name <- c("SO2", "HSO3", "HNO2", "HNO3", "CO2", "H2O2", "FORMIC_ACID", "HO2", "NH3", "HADRAZINE", "ATRA", "CL2", "HCL",  
                     "HOCL", "HBR",  "HI",  "ACRYACID")
      B_var     <- c(1.30E-02, 6.60E-08, 5.10E-04, 1.54E+01, 4.30E-07, 4.68E-11, 2.20E-12, 1.80E-04, 3.50E-05, 1.70E-05, 
                     1.00E-14, 2.09E-02, 5.01E-04, 3.16E-08, 1.74E+06, 1.11E-08, 1.00E+09, 3.20E+09, 5.62E-05)
      D_var     <- c(1.96E+03, 1.50E+03,-1.26E+03, 8.70E+03,-1.00E+03,-1.76E+03,-3.73E+03,-2.00E+01, 0.00E+00,-4.50E+02,
                    -6.71E+03, 0.00E+00, 0.00E+00, 0.00E+00, 6.90E+03, 0.00E+00, 0.00E+00, 0.00E+00, 0.00E+00)
 
# Get vector of species being modeled.

      h.dat.map <- match(Spc.names,Spc)

      TFAC <- ( 298.0 - TEMP) / ( 298.0 * TEMP )
      KH   <- A_var[h.dat.map] * exp( E_var[h.dat.map] * TFAC )
      HLCONST <- KH

      effective  <- which(Spc.names%in%diss_name==T)
      effect_map <- which(diss_name%in%Spc.names==T)
      if(length(effective)>0){
         HPLUS    <- 10^(-pH)
         HPLUSI   <- 1.0 / HPLUS
         HPLUS2I  <- HPLUSI * HPLUSI
         CLMINUS  <- 2.0E-03                # chlorine ion conc [CL-]
         CLMINUSI <- 1.0 / CLMINUS          # 1 / CLMINUS
         if(diss_name[effect_map]=='SO2') {
            AKEQ1   <- B_var[effect_map[LSO2]]  * exp( D_var[ LSO2 ]  * TFAC )
            AKEQ2   <- B_var[ LHSO3 ] * exp( D_var[ LHSO3 ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI + AKEQ1 * AKEQ2 * HPLUS2I )
         }   
         if(diss_name[effect_map]=='HNO2') {
            AKEQ1   <- B_var[ LHNO2 ] * exp( D_var[ LHNO2 ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='HNO3') {
            AKEQ1   <- B_var[ LHNO3 ] * exp( D_var[ LHNO3 ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='CO2') {
          AKEQ1   <- B_var[ LCO2 ]  * exp( D_var[ LCO2 ]  * TFAC )
          AKEQ2   <- B_var[ LHCO3 ] * exp( D_var[ LHCO3 ] * TFAC )
          HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI + AKEQ1 * AKEQ2 * HPLUS2I )
         }
         if(diss_name[effect_map]=='H2O2') {
            AKEQ1   <- B_var[ LH2O2 ] * exp( D_var[ LH2O2 ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='FORMIC_ACID') {
            AKEQ1   <- B_var[ LHCOOH ] * exp( D_var[ LHCOOH ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='HO2') {
            AKEQ1   <- B_var[ LHO2 ] * exp( D_var[ LHO2 ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='NH3') {
            AKEQ1   <- B_var[ LNH4OH ] * exp( D_var[ LNH4OH ] * TFAC )
            AKEQ2   <- B_var[ LH2O ] * exp( D_var[ LH2O ] * TFAC )
            OHION   <- AKEQ2 * HPLUSI
            HLCONST <- KH * ( 1.0 + AKEQ1 / OHION )
         }
         if(diss_name[effect_map]=='HYDRAZINE') {
            AKEQ1   <- B_var[ LHYDRAZINE ] * exp( D_var[ LHYDRAZINE ] * TFAC )
            AKEQ2   <- B_var[ LH2O ] * exp( D_var[ LH2O ] * TFAC )
            OHION   <- AKEQ2 * HPLUSI
            HLCONST <- KH * ( 1.0 + AKEQ1 / OHION )
         }
         if(diss_name[effect_map]=='ATRA' | diss_name[effect_map]=='DATRA' ) {
            AKEQ1   <- B_var[ LATRA ] * exp( D_var[ LATRA ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='CL2') {
            AKEQ1   <- B_var[ LCL2 ]  * exp( D_var[ LCL2 ] * TFAC )
            AKEQ2   <- B_var[ LHOCL ] * exp( D_var[ LHOCL ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI * CLMINUSI + AKEQ1 * AKEQ2 * HPLUS2I * CLMINUSI )
         }
         if(diss_name[effect_map]=='CL2') {
            AKEQ1   <- B_var[ LHCL ] * exp( D_var[ LHCL ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='HI') {
            AKEQ1   <- B_var[ LHI ] * exp( D_var[ LHI ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
         if(diss_name[effect_map]=='ACRYACID') {
            AKEQ1   <- B_var[ LACRYACID ] * exp( D_var[ LACRYACID ] * TFAC )
            HLCONST <- KH * ( 1.0 + AKEQ1 * HPLUSI )
         }
      }
      return(HLCONST)
}
#----------------------------------------------------------------------------------------------------------------------------------------
# Stomatal resistance
#----------------------------------------------------------------------------------------------------------------------------------------
get.le <- function(soilm, wwlt, wfc, Ta, Ts, Tg, Rs, Rst_min, lai, ustar, Ubar, RH, l_leaf,s_wet,snow,wsat,wres,bslp,pH,Ht) {
# inputs:
# Volumetric soil moisture (soilm) m3/m3
# Soil field capacity at -33 kPa (wfc) m3/m3
# Soil wilting point at -1500 kPa (wwlt) m3/m3
# air temperature (Ta) K
# surface temperature (Ts) K
# Solar radiation (Rs) W/m2
# Minimum stomatal resistance (Rst_min) s/m
# leaf area index (lai) m2/m2
# surface friction velocity (ustar) m/s
# Mean wind speed (Ubar) m/s
# relative humidity (RH) %
# estimate the stomatal resistance based on WRF 4.0 P-X parameterization
   rsmax <- 5000       # m/s
   f3min  <- 0.25      # 
   ftmin <- 0.0000001  # m/s
   dwat  <- 0.2178     # [cm^2/s] at 273.15K
# saturation water vapor 
   svp1 <- 0.611
   svp2 <- 17.67
   svp3 <- 29.65
   Name <- 'H2O'
   s.dat <- read.csv("STAGE_DATA.csv")
   Species <- as.character(s.dat$Species)
   ns      <- which(Species%in%Name==T)
   H.Name  <- as.character(s.dat$H)[ns]
   f0    <- s.dat$f0[ns]
   es    <- 10*svp1*exp(svp2*(Ta-273.15)/(Ta-svp3))
   qsm   <- es*0.622/(1000-es)
   qvs   <- 0.01*RH*qsm
   qgs   <- 1000*es*0.622/(1000-es)
   esg   <- 10*svp1*exp(svp2*(Ts-273.15)/(Ts-svp3))
   qss   <- es*0.622/(1000-es)
   qsg   <- esg*0.622/(1000-esg)

   Ra <- get.ra(ustar, Ubar, VonK)

   w2avail <- soilm*exp( (0.5-zsoil) * 9.81 )^(1.0/bslp)-wwlt
   w2mxav  <- wfc - wwlt
   f2      <- 1.0 / ( 1.0 + exp( -5.0 * ( w2avail / w2mxav - ( w2mxav / 3.0 + wwlt ) ) ) )  
#   f2 <- 1
# air temperature function
   f4 <- 1.0 / ( 1.0 + exp( 0.5 * (Ta - 314.0 ) ) )
   f4.low <-which( Ta <= 302.15 ) 
   f4[f4.low] <-  1.0 / ( 1.0 + exp( -0.41 * (Ta[f4.low] - 282.05 ) ) )

# Radiation 
   par <- 0.45*Rs*4.566

   if(length(Rst_min)==1){
      if(Rst_min > 130 ){
         f1max <- 1.0-0.02*lai
      } else {
         f1max <- 1.0-0.07*lai
      }
      f1     <- f1max * ( 1.0 - exp( -0.0017 * par ) )
      f1[f1<Rst_min / rsmax] <- Rst_min / rsmax
   } else {
      f1max <- Rst_min*0
      f1max[Rst_min > 130] <- 1.0-0.02*lai[Rst_min > 130]
      f1max[Rst_min <= 130] <- 1.0-0.07*lai[Rst_min <= 130]
      f1     <- f1max * ( 1.0 - exp( -0.0017 * par ) )
      f1[f1<Rst_min / rsmax] <- Rst_min[f1<Rst_min / rsmax] / rsmax
   }
   ftot   <- lai * f1 * f2 * f4
#   ftot   <- max( ftot,ftmin )
   ftot[ftot<ftmin]   <- ftmin 
   fshelt <- 1.0   # go back to NP89
   gs     <- ftot / ( Rst_min * fshelt )
   Rbw <- get.rb_leaf("H2O",ustar, lai, Ta, l_leaf)
   raw    <- Ra + Rbw
   ga     <- 1.0 / raw
#-- Compute humidity effect according to RH at leaf surf
   f3     <- 0.5 * ( gs - ga + sqrt( ga * ga + ga * gs * ( 4.0 * qvs / qss - 2.0 ) + gs * gs ) ) / gs
   f3[f3<f3min] <- f3min
   f3[f3>1] <- 1
   f3[f3<f3min] <- f3min
   f3[f3>1] <- 1
   Rst    <- 1.0 / ( gs * f3 )   
   Rst[Rst<Rst_min/lai] <- Rst_min/lai[Rst<Rst_min/lai]
# Estimate ET
   btap   <- 0.25*(1-cos(soilm/wfc*pi))
   btap[soilm >= wfc]   <- 1
   sigg  <- array(0,length(l_wet))
   sigg[ qss >= qvs ] <- l_wet[ qss >= qvs ]
   sigg[ qss <  qvs ] <- 1
   Rgrnd   <- get.r_grnd("H2O",ustar,Ubar,Ta,Ts,Tg,lai,s_wet,snow,soilm,wsat,wfc,wwlt,wres,bslp,pH,Ht)

   qss_leaf <- (qvs/(Ra*Rbw)+                                                                # Atmospheric Component
                qss*(1.0/(Ra*Rst)+1.0/(Rbw*Rst)+1.0/(Rgrnd*Rst))+                              # Stomatal Component              
                qsg/(Rbw*Rgrnd))/                                                              # Soil Component
               (1.0/(Ra*Rbw) +1.0/(Ra*Rst) +1.0/(Rbw*Rgrnd)+
                1.0/(Rbw*Rst)+1.0/(Rgrnd*Rst)) # Least common denominator

#Calculate the canopy compensation point follwing Nimitz et al 2001 modified to account for a cuticular compsensation point
   qss_z0     = (qvs/Ra+qss_leaf/Rbw+qss/Rgrnd)/(1.0/Ra+1.0/Rbw+1.0/Rgrnd)           

# calculate the flux
   et     <- -veg  * (qvs-qss_z0)/Ra -                 # air-vegetation flux
           (1-veg) * (qvs-qsg)/( Ra + Rgrnd )     # air-soil flux

# latent heat flux
   lv       <- (2.501-0.00237*(Ta-273.15))*1e6
   lem    <- et*lv
   return(lem)
}

