R/eb.R

In WaterMap: Mapping evapo-transpiration from satellite images.

# ENERGY BALANCE CODE

deltat_wim <-function(surface_temperature)
{
	# Wim delta T (soil-air) equation
	result <- 0.3225 * surface_temperature - 91.743
	result[result<1]<-1.0
	result[result>13]<-13.0
	return(result)
}

solar_day<-function(latitude, doy, tsw )
{
	#Sun-Earth Distance (ds A.U.)
	# Average Solar Diurnal Radiation after Bastiaanssen (1995) 
	ds <- 1.0 + 0.01672 * sin(2*pi*(doy-93.5)/365.0)
	#Solar declination (delta radians)
	deltarad <- 0.4093*sin((2*pi*doy/365)-1.39)
	#Convert latitude in radians
	latrad <-  latitude * pi / 180.0
	#Convert latitude in radians
	ws <- acos(-tan(latrad)*tan(deltarad))
	cosun <- ws*sin(deltarad)*sin(latrad)+cos(deltarad)*cos(latrad)*sin(ws)
	result <- ( cosun * 1367 * tsw ) / ( pi * ds * ds )
	return(result)
}

rnet_day<-function( bbalb,  solar,  tsw )
{
	#Average Diurnal Net Radiation after Bastiaanssen (1995)
	result <- ((1.0 - bbalb)*solar)-(110.0*tsw)
	return(result)
}

etpot_day<-function( rnetd,  tempk,  roh_w )
{
	#Average Diurnal Potential ET after Bastiaanssen (1995)
	latent<-(2.501-(0.002361*(tempk-273.15)))*1000000.0
	result <- (rnetd*86400*1000.0)/(latent*roh_w)
	return(result)
}

rnet<-function( bbalb,  ndvi,  tempk,  dtair,   e0,  tsw,  doy,  utc,  sunzangle )
{
	#Instantaneous net radiation (W/m2)
	# Tsw <-  atmospheric transmissivity single-way (~0.7 -) */
	# DOY <- Day of Year */
	# utc <- UTC time of sat overpass*/
	# sunzangle <- sun zenith angle at sat. overpass */
	# tair <- air temperature (approximative, or met.station) */
	tsw_for_e_atm<-0.7 #Special tsw, consider it a coefficient
	#Atmospheric emissivity (Bastiaanssen, 1995)
	e_atm	<-  1.08 * (-log10(tsw_for_e_atm))^0.265 
	#Atmospheric emissivity (Pawan, 2004)
#	e_atm	<- 0.85 * (-log10(tsw))^0.09
	ds <- 1.0 + 0.01672 * sin(2*pi*(doy-93.5)/365)
	delta <- 0.4093*sin((2*pi*doy/365)-1.39)
	#Kin is the shortwave incoming radiation
	Kin	<- 1367.0 * (cos(sunzangle*pi/180) * tsw / (ds*ds) )
	#Lin is incoming longwave radiation
	Lin	<- (e_atm) * 5.67 * 10^(-8) * (tempk-dtair^4)
	#Lout is surface grey body emission in Longwave spectrum
	Lout	<- e0 * 5.67 * 10^(-8) * (tempk^4)
	#Lcorr is outgoing longwave radiation "reflected" by the emissivity
	Lcorr	<- (1.0 - e0) * Lin
	result	<- (1.0 - bbalb) * Kin + Lin - Lout - Lcorr  
	return(result)
}

g0<-function(  bbalb,  ndvi,  tempk,  rnet,  time=11.0, roerink=FALSE)
{
	#Soil Heat Flux
	if(time<9.0 || time>15.0){
		r0_coef <- 1.1
	} else if (time>9.0 && time<=11.0){
		r0_coef <- 1.0
	} else if (time>11.0 && time<=13.0){
		r0_coef <- 0.9
	} else if (time>13.0 && time<=15.0){
		r0_coef <- 1.0
	}
	a <- (0.0032 * (bbalb/r0_coef) + 0.0062 * (bbalb/r0_coef) * (bbalb/r0_coef))
	b <- (1 - 0.978 * (ndvi^4))
	#Spain (Bastiaanssen, 1995)
	result <- (rnet * (tempk-273.15) / bbalb) * a * b 
	#HAPEX-Sahel (Roerink, 1995)
	if(roerink==TRUE){
		result <- result * 1.430 - 0.0845
	}
	return(result)
}

dtair0<-function( t0_dem,  tempk_water=298.0,  tempk_desert=333.0)
{
	#DTAIR Initialization
	# Pixel-based input required are: tempk water & desert
	# * additionally, dtair in Desert is vaguely initialized
	ZERO<-273.15	
	if(tempk_desert >= (ZERO+48.0)){
		dtair_desert_0 <- 13.0
	} else if(tempk_desert >= (ZERO+40.0) && tempk_desert < (ZERO+48.0)){
		dtair_desert_0 <- 10.0
	} else if(tempk_desert >= (ZERO+32.0) && tempk_desert < (ZERO+40.0)){
		dtair_desert_0 <- 7.0
	} else if(tempk_desert >= (ZERO+25.0) && tempk_desert < (ZERO+32.0)){
		dtair_desert_0 <- 5.0
	} else if(tempk_desert >= (ZERO+18.0) && tempk_desert < (ZERO+25.0)){
		dtair_desert_0 <- 3.0
	} else if(tempk_desert >= (ZERO+11.0) && tempk_desert < (ZERO+18.0)){
		dtair_desert_0 <- 1.0
	} else {
		dtair_desert_0 <- 0.0
#		printf("WARNING!!! dtair_desert_0 is NOT VALID!\n")
	}
	a <- (dtair_desert_0-0.0)/(tempk_desert-tempk_water)
	b <- 0.0 - a * tempk_water
	result <- t0_dem * a + b
	return(result)
}

dtair<-function( t0_dem,  tempk_water=298.0,  tempk_desert=333.0,  dtair_desert=13.0)
{
	#DTAIR Standard equation
	# Pixel-based input required are: tempk water & desert
	# additionally, dtair in Desert point should be given
	#Only t0_dem (altitude corrected DEM) is Array
	a <- (dtair_desert-0.0)/(tempk_desert-tempk_water)
	b <- 0.0 - a * tempk_water
	result <- t0_dem * a + b
	return(result)
}

dtair_desert<-function( h_desert, roh_air_desert=1.12, rah_desert=10.0)
{
	#DTAIR Dry Pixel
	#Only h_desert input is Array
	result <- (h_desert * rah_desert)/(roh_air_desert * 1004.0)
	return(result)
}

h0<-function( roh_air, rah, dtair)
{
	#Sensible Heat Flux Standard Equation
	#All inputs are Arrays
	result <- roh_air*1004.0*(dtair) / rah
	return (result)
}

psi_h<-function( t0_dem, h, U_0, roh_air)
{
	#Psichrometric parameter for heat momentum
	#All inputs are Arrays
	if(h != 0.0){
		n5_temp <- (-1004* roh_air*(U_0^3)* t0_dem)/(0.41*9.81* h)
	} else {
		n5_temp <- -1000.0
	}
	if(n5_temp < 0.0){
		n12_mem <- ((1-16*(2/n5_temp))^0.25)
		n11_mem <- (2*log10((1+(n12_mem^2))/2))
	} else {
		n12_mem <- 1.0
		n11_mem <- -5*2/n5_temp
	}
	result <- n11_mem
	return (result)
}

rah0<-function( zom_0, u_0)
{
	#Aerodynamic resistance to heat momentum initialization
	result <- log10(2/0.01)/(u_0*0.41)
	return (result)
}

rah1<-function( psih, ustar)
{
	#Aerodynamic resistance to heat momentum standard equation
	result <- (log10(2/0.01)-psih)/(ustar*0.41)
	return (result)
}

rohair0<-function(tempk)
{
	#Air density Initialization
	A <- (18.0*(6.11*exp(17.27*34.8/(34.8+237.3)))/100.0)
	B <- (18.0*(6.11*exp(17.27*34.8/(34.8+237.3)))/100.0)
	result <- (1000.0 - A)/(tempk*2.87)+ B/(tempk*4.61)
	return(result)
}

rohair<-function( dem, tempk, dtair)
{
	#Air density standard paramterization
	a <- tempk - dtair
	b <- (( a - 0.00627*dem)/a)
	result <- 349.467 * ( b ^ 5.26)/ a 
	return(result)
}

U0<-function( zom0, u2m)
{
	#Wind Speed Initialization
	u0<-u2m*0.41*log10(200/(0.15/7))/(log10(2/(0.15/7))*log10(200/zom0)) 
	return(u0)
}

ustar<-function( t0_dem, h, ustar, roh_air, zom, u2m)
{
	#Nominal Wind Speed
# 	 n5_temp # Monin-Obukov Length 		*/
#        n10_mem # psi m 			*/
# 	 n31_mem # x for U200 (that is bh...) 	*/
	 hv<-0.15	# crop height (m) 		*/
	 bh<-200	# blending height (m) 		*/
	
	if(h != 0.0){
		n5_temp <- (-1004* roh_air*(ustar^3)* t0_dem)/(0.41*9.81* h)
	} else {
		n5_temp <- -1000.0
	}
	if(n5_temp < 0.0){
		n31_mem <- ((1-16*(200/n5_temp))^0.25)
		n10_mem <- (2*log10((1+n31_mem)/2)+log10((1+(n31_mem^2))/2)-2*atan(n31_mem)+0.5*pi)
	} else {
#		n31_mem <- 1.0
		n10_mem <- -5*2/n5_temp
	}
	result <- ((u2m*0.41/log10(2/(hv/7)))/0.41*log10(bh /(hv/7)*0.41))/(log10(bh / zom)-n10_mem)
	return(result)
}


zom0<-function( ndvi, ndvi_max)
{
	#Roughness length for heat momentum 
	hv_ndvimax<-1.5 # crop vegetation height (m) */
	hv_desert<-0.002 # desert base vegetation height (m) */
	a <- (log10(hv_desert)-((log10(hv_ndvimax/7)-log10(hv_desert))/(ndvi_max-0.02)*0.02))
	b <- (log10(hv_ndvimax/7)-log10(hv_desert))/(ndvi_max-0.02)* ndvi
	result <- exp(a+b) 
	return(result)
}


sensih<-function( iteration, tempk_water, tempk_desert, t0_dem, tempk, ndvi, ndvi_max, dem, rnet_desert, g0_desert, t0_dem_desert, u2m, dem_desert)
{
	#Sensible Heat flux determination
	# This is the main loop used in SEBAL */
	# Arrays Declarations */
	# define ITER_MAX 10
	ITER_MAX<-10
	# Arrays Declarations */
	# Fat-free junk food */
	if (iteration>ITER_MAX){
		iteration<-ITER_MAX
	}

#	dtair[0] 	<- dt_air_0(t0_dem, tempk_water, tempk_desert)
	dtair[0] <- 5.0
#	printf("*****************************dtair <- %5.3f\n",dtair[0])
	roh_air[0] 	<- rohair0(tempk)
#	printf("*****************************rohair<-%5.3f\n",roh_air[0])
	roh_air_desert 	<- rohair0(tempk_desert)
#	printf("**rohairdesert <- %5.3f\n",roh_air_desert)
	zom0 		<- zom0(ndvi, ndvi_max)
#	printf("*****************************zom <- %5.3f\n",zom0)
	u_0 		<- U0(zom0, u2m)
#	printf("*****************************u0\n")
	rah[0] 		<- rah0(zom0, u_0)
#	printf("*****************************rah <- %5.3f\n",rah[0])
	h[0] 		<- h0(roh_air[0], rah[0], dtair[0])
#	printf("*****************************h\n")

#----------------------------------------------------------------*/
#Main iteration loop of SEBAL*/
	zom[0] <- zom0
	for(ic in 1:(iteration+1)){
		# Where is roh_air[i]? */
		psih <- psi_h(t0_dem,h[ic-1],u_0,roh_air[ic-1])
		ustar[ic] <- ustar(t0_dem,h[ic-1],u_0,roh_air[ic-1],zom[0],u2m)
		rah[ic] <- rah1(psih, ustar[ic])	
		# get desert point values from maps */
		if(ic==1){
			h_desert	<- rnet_desert - g0_desert
			zom_desert	<- 0.002
			psih_desert 	<- psi_h(t0_dem_desert,h_desert,u_0,roh_air_desert)
			ustar_desert	<- ustar(t0_dem_desert,h_desert,u_0,roh_air_desert,zom_desert,u2m)
		} else {
			roh_air_desert	<- rohair(dem_desert,tempk_desert,dtair_desert)
			h_desert	<- h0(roh_air_desert,rah_desert,dtair_desert)
			ustar_desertold <- ustar_desert
			psih_desert 	<- psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert)
			ustar_desert	<- ustar(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,zom_desert,u2m)
		}
		rah_desert	<- rah1(psih_desert,ustar_desert)
		dtair_desert 	<- dtair_desert(h_desert, roh_air_desert, rah_desert)
		# This should find the new dtair from inversed h equation...*/
		dtair[ic] 	<- dtair(t0_dem, tempk_water, tempk_desert, dtair_desert)
		# This produces h[ic] and roh_air[ic+1] */
		roh_air[ic] 	<- rohair(dem, tempk, dtair[ic])
		h[ic] 		<- h0(roh_air[ic], rah[ic], dtair[ic])
	}
	return(h[iteration])
}

.powr <- function(x, y) { return (x^y) }

sensih_SEBAL01<-function(rnet, g0, zom, t0_dem, ustar, ea, row_wet, col_wet, row_dry, col_dry)
{
# rnet, g0, zom, t0_dem are raster data
# Example point data
#     ustar<-0.32407
#     ea<-1.511
    Rn_dry <- xyValues(rnet,cbind(row_dry,col_dry))
    g0_dry <- xyValues(g0,cbind(row_dry,col_dry))
    t0dem_dry <- xyValues(t0_dem,cbind(row_dry,col_dry))
    t0dem_wet <- xyValues(t0_dem,cbind(row_wet,col_wet))
    h_dry <- Rn_dry - g0_dry

    u5 <- (ustar / 0.41) * log(5 / zom)
    rah1<-(1/(u5*.powr(0.41,2)))*log(5/zom)*log(5/(zom*0.1))
    roh1<-((998-ea)/(t0_dem*2.87))+(ea/(t0_dem*4.61))
    roh1[roh1 > 5] <- 1.0
    roh1[roh1 <= 5]<-((1000-4.65)/(t0_dem*2.87))+(4.65/(t0_dem*4.61))
    rah_dry <- xyValues(rah1,cbind(row_dry,col_dry))
    roh_dry <- xyValues(roh1,cbind(row_dry,col_dry))
    Roh <- roh1
    Rah <- rah1

    d_dT_dry <- (h_dry * rah_dry) / (1004 * roh_dry)

    # Calculate coefficients for next dT equation 
    a <- 1.0/ ((d_dT_dry-0.0) / (t0dem_dry-t0dem_wet)) 
    b <- ( a * t0dem_wet ) * (-1.0) 
    sumx <- t0dem_wet + t0dem_dry
    sumy <- d_dT_dry + 0.0
    sumx2 <- .powr(t0dem_wet, 2) + .powr(t0dem_dry, 2)
    sumxy <- (t0dem_wet * 0.0) + (t0dem_dry * d_dT_dry)
    a <- (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (.powr(sumx, 2) / 2.0))
    b <- (sumy - (a * sumx)) / 2.0

    #      ITERATION 1 
    rah1 <- xyValues(Rah,cbind(row,col))
    roh1 <- xyValues(Roh,cbind(row,col))
    d_h1[rah1 < 1.0] <- 0.0
    d_h1[rah1 >= 1.0] <- (1004 * roh1) * (a * t0_dem + b) / rah1
    d_L <- -1004*roh1*.powr(ustar,3)*t0_dem/(d_h1*9.81*0.41)
    d_x <- .powr((1-16*(5/d_L)),0.25)
    d_psim <-2*log((1+d_x)/2)+log((1+.powr(d_x,2))/2)-2*atan(d_x)+0.5*pi
    d_psih <-2*log((1+.powr(d_x,2))/2)
    u5 <-(ustar/0.41)*log(5/zom)
    d_rah2 <- (1/(u5*.powr(0.41,2)))*log((5/zom)-d_psim)*log((5/(zom*0.1))-d_psih)
    rah_dry <- xyValues(d_rah2,cbind(row_dry,col_dry))
    d_h_dry <- xyValues(d_h1,cbind(row_dry,col_dry))
    Rah <- d_rah2

    #     Calculate dT_dry 
    d_dT_dry <- (d_h_dry * rah_dry) / (1004 * roh_dry)
    #     Calculate coefficients for next dT equation 
    a <- (d_dT_dry-0)/(t0dem_dry-t0dem_wet) 
    b <- (-1.0) * ( a * t0dem_wet ) 
    sumx <- t0dem_wet + t0dem_dry
    sumy <- d_dT_dry + 0.0
    sumx2 <- .powr(t0dem_wet, 2) + .powr(t0dem_dry, 2)
    sumxy <- (t0dem_wet * 0.0) + (t0dem_dry * d_dT_dry)
    a <- (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (.powr(sumx, 2) / 2.0))
    b <- (sumy - (a * sumx)) / 2.0

    #      ITERATION 2 
    d_rah2 <- Rah
    roh1 <- Roh
    d_h2[d_rah2 < 1.0] <- 0.0
    d_h2[d_rah2 < 1.0] <-(1004*roh1)*(a*t0_dem+b)/d_rah2
    d_L <- -1004*roh1*.powr(ustar,3)*t0_dem/(d_h2*9.81*0.41)
    d_x <- .powr((1 - 16 * (5 / d_L)), 0.25)
    d_psim <-2*log((1+d_x)/2)+log((1+.powr(d_x,2))/2)-2*atan(d_x)+0.5*pi
    d_psih <-2*log((1+.powr(d_x,2))/2)
    u5 <-(ustar/0.41)*log(5/zom)
    d_rah3<-(1/(u5*.powr(0.41,2)))*log((5/zom)-d_psim)*log((5/(zom*0.1))-d_psih)
    rah_dry <- xyValues(d_rah2,cbind(row_dry,col_dry))
    d_h_dry <- xyValues(d_h2,cbind(row_dry,col_dry))
    Rah <- d_rah2

    #     Calculate dT_dry 
    d_dT_dry <- (d_h_dry * rah_dry) / (1004 * roh_dry)
    #     Calculate coefficients for next dT equation 
    a <- (d_dT_dry-0)/(t0dem_dry-t0dem_wet) 
    b <- (-1.0) * ( a * t0dem_wet ) 
    sumx <- t0dem_wet + t0dem_dry
    sumy <- d_dT_dry + 0.0
    sumx2 <- .powr(t0dem_wet, 2) + .powr(t0dem_dry, 2)
    sumxy <- (t0dem_wet * 0.0) + (t0dem_dry * d_dT_dry)
    a <- (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (.powr(sumx, 2) / 2.0))
    b <- (sumy - (a * sumx)) / 2.0

    #      ITERATION 3 
    d_rah3 <- Rah
    roh1 <- Roh
    d_h3[d_rah3 < 1.0] <- 0.0
    d_h3[d_rah3 >= 1.0] <- (1004 * roh1) * (a * t0_dem + b) / d_rah3
    d_h3[d_h3 < 0 && d_h3 > -50] <- 0.0
    d_h3[d_h3 < -50 || d_h3 > 1000] <- NA
    return (d_h3)
}


evapfr<-function( rnet, g0, h0)
{
	#Evaporative Fraction
	result <- (rnet - g0 - h0) / (rnet - g0)
	return(result)
}


soilmoisture<-function( evapfr )
{
	#soil moisture in the root zone
	#Makin, Molden and Bastiaanssen, 2001
	result <- (exp((evapfr-1.2836)/0.4213))/0.511
	return(result)
}


eta<-function( rnet_day, evapfr, surface_temperature)
{
	#Evapotranspiration from energy balance
	t_celsius <- surface_temperature - 273.15 
	latent 	  <- 86400.0/((2.501-0.002361*t_celsius)*(10^6))
	result 	  <- rnet_day * evapfr * latent 
	return(result)
}