// 8-25-2008 

// ET_pbl.c


// Dennis Baldocchi
// Dept Environmental Science, Policy and Mgt
// 151 Hilgard Hall
// University of California, Berkeley
// Berkeley, CA

// Program to compute latent and sensible heat fluxes with
// a quadratic form of the surface energy balance

// It use Isothermal equations to compute LE as a function of Q_In, which allows
// Tsfc and Lout to float and redefine Rnet

// This version is linked to a simple PBL model of McNaughton and Spriggs.
// It considers feedback between fluxes, growth of the PBL and concentrations
// and temperature in the mixed layer.  The current version was for
// a sensitivity test, so it only inputted PAR and estimated the temporal
// change in T and q.  For operational use, T, q and u should be
// inputed to the model.

// Compare Cases with and without PBL feedbacks

// debug

// 8-25-08 Y Ryu found an error in the DES2DT function

// correct Qin for Lin=> ep * Lin

// got energy balance closure and LEiso = LEquad!

// need to re-initialize Tair and pbl for each day loop

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <io.h>
#include <stdlib.h>
#include <process.h>
#include <string.h>

double D2ESDT (T);
void QUADRATIC ();
double ES (T);
double DESDT (T);
double LAMBDA (T);
double ENERB (FTS);
double SKYIR(T);


// declare Structures

struct parameters
{
double ep;     // emissivity
double sigma;  // Stefan Boltzmann constant 
double Cp;     // Specific heat of air, J kg-1 K-1
double Rstar;  // Universal Gas Constant, J mole-1 K-1
double dLdT;   // derivative of latent heat of evaporation with temperature
double Mair;   // Molecular mass of air, g mole-1
double del_t;   // time step, s

int I_start;   // start time
int I_stop;    // stop time

} param;

struct energy_balance
{
double H;     // sensible heat flux density, W m-2
double H_virt; // virtual heat flux density, W m-2
double LE;    // latent heat flux density, W m-2
double LE_PM; // Penman Monteith LE, W m-2
double LE_eq; // equilibrium evaporation, W m-2
double LE_iso; // isothermal PM
double Gsoil; // soil heat flux density, W m-2
double Lout;  // outgoing longwave energy, W m-2
double Lout_Ta;  // outgoing longwave at Ta
double Rnet;  // net radiation, W m-2
double Avail; // available energy, Rg(1-a)+Lin, W m-2
double Rniso; // isothermal radiation
}energy;


struct meteorology
{
double Tair_C;   // air temperature, Celcius
double Tair_K;   // air temperature, Kelvin
double Tsfc_K;   // surface temperature, Kelvin
double eair_kPa; // vapor pressure humidity, kPa
double eair_Pa;  // vapor pressure, Pa
double vpd_Pa;   // vapor pressure deficit, Pa
double RgIn;     // incoming global radiation, W m-2
double RgInMax;  // maximum global radiation
double albedo;   // albedo
double LongIn;   // incoming longwave radiation, W m-2
double Q_In;	 //net incoming short and long wave, Rg(1-a)+ep Lin
double P_kPa;   // pressure, kPa
double P_Pa;    // pressure, Pa
double U;       // wind speed, m s-1
double ustar;   // friction velocity, m s-1
double pbl_z;   // pbl height, m
double del_zi;  // change in pbl height, m
}met;


struct variables
{

double Rs;  // stomatal resistance, s m-1
double Ra;  // boundary layer resistance, s m-1
double Rb;  // quasi laminar boundary resistance, s m-1
double Rheat;  // boundary layer resistance to heat transfer, s m-1
double Rvapor; // boundary layer resistance to vapor transfer, s m-1
double Kc;  // convective heat conductance, m s-1
double air_density;  // air density, kg m-3
double del_T;        // delta time step in temperature, C s-1
double del_h2o;      // delta time step in humidity, Pa s-1
double psych;   // psychrometeric constant
double slope;   //des(T)/dT
double jump;    // temperature difference across pbl
double Gr;     // radiative conductance
double Gw;     // water conductance
double Gh;     // heat conductance
} var;


FILE *fptr1, *fptr2; 

void main()
{
	int I, K;

	double drag;

//	fptr1=fopen("c:\\ET_pbl\\et_savanna_nopbl.dat","w");
	fptr1=fopen("c:\\ET_pbl\\et_grass_nopbl.dat","w");
//	fptr1=fopen("d:\\ET_pbl\\et_nopbl.dat","w");
//	fptr1=fopen("d:\\ET_pbl\\et_nopbl_alb025.dat","w");

		printf("%s, %s, %s, %s, %s, %s, %s \n",
			"I", "LE", "H", "Rnet", "Tsfc_K", "e_air","pbl_z");


	fprintf(fptr1,"%s, %s, %s, %s, %s, %s, %s, %s \n",
			"I", "LE", "H", "Rnet", "Tsfc_K", "e_air","pbl_z", "Rs");

//	fptr2=fopen("c:\\ET_pbl\\etpbl_savanna.dat","w");
	fptr2=fopen("c:\\ET_pbl\\etpbl_grass.dat","w");
//	fptr2=fopen("d:\\ET_pbl\\etpbl.dat","w");
//	fptr2=fopen("d:\\ET_pbl\\etpblalb025.dat","w");

		printf("%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s \n",
			"I", "LE", "LE_PM", "H", "Rnet","Rg", "Tsfc_K", "e_air","Tair", "vpd","pbl_z", "Rc");


	fprintf(fptr2,"%s, %s, %s, %s, %s, %s, %s, %s, %s,  %s, %s, %s, %s, %s \n",
			"I", "LE", "LE_iso", "LE_PM", "LE_eq", "H", "Rnet", "Rg", "Tsfc_K", "e_air", "Tair", 
			"vpd","pbl_z", "Rc");


	// Initialize Parameters

	param.ep = .98;               
	param.sigma = 5.67E-08 ;      
	param.Cp = 1005;              
	param.Rstar = 8.3144;         
	param.Mair = 29.;              
	param.dLdT = -2370.;           
	param.del_t=3600;



	// loop Rs

	K=40;

	while (K <= 3000)
	{

	var.Rs=(double)K;

	K*=(double)2.;

  // Initialize Meteorology

	met.P_kPa = 101.;  
	met.P_Pa=met.P_kPa*1000.;
	met.U = 2.5;
	met.Tair_C=15;
	met.Tair_K=273.15+met.Tair_C;
	met.eair_kPa=.2;
	met.eair_Pa=met.eair_kPa*1000.;
	met.pbl_z = 150.;
	met.RgInMax=950.;
	met.albedo=0.30;  // albedo = 0.25; grass spring, 0.3 summer  // albedo = 0.15; savanna

	fprintf(fptr2," albedo, %f\n",met.albedo);
	
    energy.Gsoil = 50.;

	// Initialize Derived Quantities

	met.LongIn=SKYIR(met.Tair_K);
	met.Q_In=met.RgIn*(1.-met.albedo)+param.ep*met.LongIn-energy.Gsoil;
	met.vpd_Pa=ES(met.Tair_K)-met.eair_Pa;
	drag = 0.1;								// drag = 0.1 grass; drag = 0.2 savanna
	met.ustar = drag* sqrt(met.U*met.U);      // u* savanna = 2 u* grass
 	var.air_density = met.P_kPa * param.Mair / (param.Rstar * met.Tair_K); 
	var.jump=3./1000.;  

	// Resistances and Conductances

	var.Ra = met.U / (met.ustar * met.ustar);
	var.Rb= 2. / (.4 * met.ustar);
	var.Rheat=var.Ra+var.Rb;
	var.Gh=1./var.Rheat;
	var.Gw=1/(var.Rheat+var.Rs);
	
    param.I_start=5;
	param.I_stop=19;


	// Case without PBL Feedback

	for(I=param.I_start;I<=param.I_stop; I++)
	{


	if(I==param.I_start)
	{
	met.Tair_C=15;
	met.Tair_K=273.15+met.Tair_C;
	met.eair_kPa=1.;  //.2
	met.eair_Pa=met.eair_kPa*1000.;
	met.pbl_z = 1000.;
	}
		// assume sinusoidal diurnal pattern in solar radiation

		met.RgIn=met.RgInMax*sin(3.1415*(I-param.I_start)/(param.I_stop-param.I_start));
		
		// for convenience, Gsoil is substracted from Qin
		
		met.Q_In=met.RgIn*(1.-met.albedo)+param.ep*met.LongIn;
		
		QUADRATIC();

		energy.Rnet= met.RgIn*(1-met.albedo)+param.ep*met.LongIn-energy.Lout;
	energy.Avail = energy.Rnet - energy.Gsoil;

	printf("%d, %f, %f, %f, %f, %f, %f, %f, %f\n", 
		I, energy.LE, energy.H,energy.Rnet, met.Tsfc_K, met.eair_Pa,met.pbl_z,var.Rs);
       
	fprintf(fptr1,"%d, %f, %f, %f, %f, %f, %f, %f\n", 
		I, energy.LE, energy.H,energy.Rnet, met.Tsfc_K, met.eair_Pa,met.pbl_z, var.Rs);
	
	} // end loop

	// Case with PBL Feedback

	for(I=param.I_start;I<=param.I_stop; I++)
	{

	

	if(I==param.I_start)
	{
	met.Tair_C=15;
	met.Tair_K=273.15+met.Tair_C;
	met.eair_kPa=1.;  //.2
	met.eair_Pa=met.eair_kPa*1000.;
	met.pbl_z = 150.;
	}
		
		// assume sinusoidal diurnal pattern in solar radiation

	met.RgIn=met.RgInMax*sin(3.1415*(I-param.I_start)/(param.I_stop-param.I_start));

					
	var.psych = met.P_Pa * param.Cp / (.622 * LAMBDA(met.Tair_K));    //Pa K-1
	var.slope = DESDT(met.Tair_K);       
		
	met.Tair_K = met.Tair_C + 273.15;
	met.LongIn=SKYIR(met.Tair_K);
	met.Q_In=met.RgIn*(1.-met.albedo)+param.ep*met.LongIn;
	met.vpd_Pa=ES(met.Tair_K)-met.eair_Pa;
	var.air_density = met.P_kPa * param.Mair / (param.Rstar * met.Tair_K); 
	
	energy.Rniso=met.RgIn*(1.-met.albedo)+param.ep*met.LongIn-param.ep*param.sigma*pow(met.Tair_K,4.);

	var.Gr=4.*param.ep*param.sigma*pow(met.Tair_K,3.)/(var.air_density*param.Cp);

		// energy balance computations using quadradic solutions
		// to the surface energy balance (both for Tsfc and LE)

		QUADRATIC();

//Compute isothermal Penman-Monteith LE and equilibrium LE

		energy.Rnet= met.RgIn*(1-met.albedo)+param.ep*met.LongIn-energy.Lout;
	energy.Avail = energy.Rnet - energy.Gsoil;

energy.LE_iso = (var.slope * (energy.Rniso-energy.Gsoil) + param.Cp * var.air_density * met.vpd_Pa*(var.Gr+var.Gh))/
 (var.slope + var.psych * (var.Gr+var.Gh)/var.Gw);

energy.LE_PM = (var.slope * (energy.Avail) + param.Cp * var.air_density * met.vpd_Pa*var.Gh)/
 (var.slope + var.psych * var.Gh/var.Gw);


energy.LE_eq = var.slope * energy.Avail / (var.slope + var.psych);

	printf("%d, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", 
		I, energy.LE,energy.LE_PM, energy.LE_eq,energy.H,energy.Rnet, met.Tsfc_K, 
		met.eair_Pa,met.Tair_K,met.vpd_Pa,met.pbl_z, var.Rs);
       
    	if(I>=5 && I <= 19)
		{
	    fprintf(fptr2,"%d, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", 
		I, energy.LE,energy.LE_iso, energy.LE_PM,energy.LE_eq,energy.H,energy.Rnet, met.RgIn, 
		met.Tsfc_K,	met.eair_Pa,met.Tair_K,met.vpd_Pa,met.pbl_z, var.Rs);
		}


			// time steps in temperature and vapor pressure

		var.del_T = param.del_t * energy.H / 
			(met.pbl_z * var.air_density * param.Cp);

		var.del_h2o = param.del_t * energy.LE / 
			(met.pbl_z * var.air_density * .622* LAMBDA(met.Tair_K)/met.P_Pa);

		// Compute PBL growth

		energy.H_virt = energy.H + .07 * energy.LE;

		// In future runs scale entrainment as a fraction of surface heat flux

		if(energy.H_virt > 0)
		met.del_zi = energy.H_virt * param.del_t / 
			(var.air_density * param.Cp * met.pbl_z * var.jump);
		met.pbl_z += met.del_zi;


		// update Met variables and energy fluxes

	met.eair_Pa = met.eair_Pa + var.del_h2o;
	met.Tair_C = met.Tair_C + var.del_T;
	met.Tair_K = met.Tair_C + 273.15;
	
       
		} // end time loop


	} // loop Rs

 } // end of main

double D2ESDT (double T)
{
	double fT;


//d(g(x)f(x))/dx = f'(x)g(x) + g'(x)f(x)

// corrected 8/26/08

fT = -2. * ES(T) * LAMBDA(T) * 18. / (param.Rstar * T * T * T * 1000.) +  DESDT(T) * LAMBDA(T) * 18. / (param.Rstar * T * T * 1000.);


return fT;
}

double DESDT(double T)
{
 double fT;

// Slope of the saturation vapor pressure-temperature curve
	
//Thermodynamic relationship from HESS    (Pa K-1)

fT = ES(T) * LAMBDA(T) * 18. / (param.Rstar * T * T * 1000.);
 return fT;


}



double ES(double T) 
{
double fT;

	//SATURATION VAPOR PRESSURE,  Pa

fT = 100. * exp(54.8781919 - 6790.4985 / T - 5.02808 * log(T));
return fT;
}


double LAMBDA (double T)
{
	double fT;

//LATENT HEAT OF VAPORIZATION, J KG-1

fT = 3149000 - 2370 * T;

return fT;
}

void QUADRATIC()
{


	// Quadratic solution of LE at canopy scale

	// aLE^2 + bLE + c = 0


	double Acoef, Bcoef, Ccoef;
	double aT, bT, cT;
	double hcoef,del_Tk;
	double lecoef,le1, le2;
	double dest, d2est;
	double longcoef;


	dest=DESDT(met.Tair_K);
	d2est=D2ESDT(met.Tair_K);

	energy.Lout_Ta = param.ep * param.sigma * pow(met.Tair_K,4.0);
	lecoef=var.air_density*.622*LAMBDA(met.Tair_K)*var.Gw/(met.P_Pa);
	hcoef=var.air_density*param.Cp*var.Gh;
	longcoef=4*param.ep*param.sigma*pow(met.Tair_K,3.);

	Acoef=lecoef*d2est/(2*(hcoef+longcoef));

	Bcoef=-(hcoef+longcoef)-
			lecoef*dest + Acoef*(2.*energy.Lout_Ta +2*energy.Gsoil- 2.*met.Q_In);

	Ccoef=(hcoef+longcoef)*lecoef*met.vpd_Pa+
			lecoef*dest*(met.Q_In - energy.Lout_Ta -energy.Gsoil) +
			Acoef*(met.Q_In*met.Q_In+energy.Lout*energy.Lout_Ta + energy.Gsoil*energy.Gsoil
			-2.*energy.Lout_Ta*met.Q_In-2*met.Q_In*energy.Gsoil+2*energy.Lout_Ta*energy.Gsoil);



// solve for LE

// a LE^2 + bLE + c = 0

le1 = (-Bcoef + pow((Bcoef *Bcoef - 4. * Acoef * Ccoef),.5)) / (2. * Acoef);
le2 = (-Bcoef - pow((Bcoef * Bcoef - 4. * Acoef * Ccoef),.5)) / (2. * Acoef);

energy.LE=le2;

//solve for Ts

// aT^2 + bT + c =0


aT = 6. * param.ep * param.sigma * met.Tair_K * met.Tair_K + d2est * lecoef/2;
bT = longcoef + hcoef + dest * lecoef;
cT =  energy.Lout_Ta + lecoef * met.vpd_Pa - met.Q_In + energy.Gsoil;



del_Tk=(-bT + pow((bT * bT - 4. * aT * cT),.5)) / (2. * aT);

met.Tsfc_K= met.Tair_K + del_Tk;


// H IS SENSIBLE HEAT FLUX

        energy.H = hcoef * (del_Tk);


//LOUT IS LONGWAVE EMITTED RADIATION

        energy.Lout = param.ep * param.sigma * pow(met.Tsfc_K,4.);

return;
}


double RHOS (T)
{ double fT;

//SATURATION VAPOR DENSITY

fT = .2165 * ES(T) / T;
 return fT;
}

double SKYIR(double t)
{
        
        // Infrared radiation from sky, W m-2, using algorithm from Norman
        
        double y, ratrad, product;

		ratrad=0.8;
		product=  ((1. - .261 * exp(-.000777 * pow((273.16 - t), 2.))) * ratrad + 1. - ratrad);
    	y = param.sigma * pow(t,4.); 

        return y;
}