/*
 * SUMMARY:      SensibleHeatFlux.c - Calculate sensible heat flux
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Bart Nijssen
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-96
 * Last Change: Tue Apr  3 09:08:51 2001 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Calculate sensible heat flux
 * DESCRIP-END.
 * FUNCTIONS:    SensibleHeatFlux()
 *               NoSensibleHeatFlux()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: SensibleHeatFlux.c,v 1.3 2001/04/04 16:43:49 waichler Exp $";
#endif /* lint */

#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "massenergy.h"
#include "constants.h"
#include "brent.h"
#include "calckh.h"

/*****************************************************************************
  SensibleHeatFlux()
*****************************************************************************/
void SensibleHeatFlux(int y, int x, int Dt, float Ra, float ZRef, 
                      float Displacement, float Z0, PIXMET *LocalMet, 
                      float NetShort, float LongIn, float ETot, 
                      int NSoilLayers, float *SoilDepth, SOILTABLE *SoilType, 
                      float MeltEnergy, SOILPIX *LocalSoil)
{
  float FluxDepth;		/* Lower boundary for soil heat flux (m) */
  float HeatCapacity;		/* Soil heat capacity */
  float MaxTSurf;		/* Upper bracket for effective surface 
				   temperature (C) */
  float MinTSurf;		/* Lower bracket for effective surface
				   temperature (C) */
  float OldTSurf;		/* Effective surface temperature at the end of
				   the last timestep (C) */
  float KhEff;			/* Thermal conductivity of the soil (W/(m*C)
                                 */ 
  float TMean;			/* Average surface temperature (C) */
  float TSoilLower;             /* Temperature os the soil at FluxDepth (C)
                                 */ 
  float TSoilUpper;             /* Temperature os the soil in top layer (C)
                                 */ 
  double Tmp;                   /* Temporary value */

  OldTSurf = LocalSoil->TSurf;
  MaxTSurf = 0.5 * (LocalSoil->TSurf + LocalMet->Tair) + DELTAT;
  MinTSurf = 0.5 * (LocalSoil->TSurf + LocalMet->Tair) - DELTAT;

  /* WORK IN PROGRESS */
  /* the lower boundary for the soil heat flux is currently fixed at a depth 
     of 1.0 m */

  FluxDepth = 1.0;
  TSoilLower = LocalSoil->Temp[NSoilLayers-1];
  TSoilUpper = LocalSoil->Temp[0];

  /* Calculate the effective thermal conductivity of the soil above between
     FluxDepth and DZ_TOP */

  KhEff = CalcEffectiveKh(NSoilLayers, DZ_TOP, FluxDepth, SoilDepth, 
                          SoilType->KhDry, SoilType->KhSol, LocalSoil->Moist, 
                          SoilType->Porosity, LocalSoil->Temp);

  /*   KhEff = 1; */

  /* Calculate the effective surface temperature that makes sure that the 
     sum of the terms of the energy balance equals 0 */

  LocalSoil->TSurf = 
    RootBrent(y, x, MinTSurf, MaxTSurf, MIN_SOIL_TEMP, MAX_SOIL_TEMP,
              SurfaceEnergyBalance, Dt, Ra, ZRef,
              Displacement, Z0, LocalMet->Wind, NetShort, LongIn,
              LocalMet->AirDens, LocalMet->Lv, ETot, KhEff,
	            SoilType->Ch[0], SoilType->Porosity[0], LocalSoil->Moist[0], 
              FluxDepth, LocalMet->Tair, TSoilUpper, 
              TSoilLower, OldTSurf, MeltEnergy);

  /* Calculate the terms of the energy balance.  This is similar to the
     code in SurfaceEnergyBalance.c */

  TMean = 0.5 * (OldTSurf + LocalSoil->TSurf);

  if (LocalMet->Wind > 0.0)
    Ra /= StabilityCorrection(ZRef, Displacement, TMean, LocalMet->Tair,
                              LocalMet->Wind, Z0);
  else
    Ra = DHSVM_HUGE;

  LocalSoil->Ra = Ra;

  Tmp = TMean+273.15;
  LocalSoil->Qnet = NetShort + LongIn -
    STEFAN * (Tmp*Tmp*Tmp*Tmp);

  LocalSoil->Qs = LocalMet->AirDens * CP * (LocalMet->Tair - TMean)/Ra;

  LocalSoil->Qe = -(LocalMet->Lv * ETot)/Dt * WATER_DENSITY;

  LocalSoil->Qg = KhEff * (TSoilLower - TMean)/FluxDepth;

  HeatCapacity = (1 - SoilType->Porosity[0]) * SoilType->Ch[0];

  if (TSoilUpper >= 0.0)
    HeatCapacity += LocalSoil->Moist[0] * CH_WATER;
  else 
    HeatCapacity += LocalSoil->Moist[0] * CH_ICE;

  LocalSoil->Qst = (HeatCapacity * (OldTSurf - TMean) * DZ_TOP)/Dt;

  LocalSoil->Qrest = LocalSoil->Qnet + LocalSoil->Qs + LocalSoil->Qe +
    LocalSoil->Qg + LocalSoil->Qst + MeltEnergy;
}

/*****************************************************************************
  Function name: NoSensibleHeatFlux()

  Purpose      : Calculate latent heat flux in W/m2

  Required     : 
    int Dt             - Model timestep (seconds)
    PIXMET LocalMet    - Met data for the current pixel 
    float ETot         - Total vapor flux (mm/timestep) 
    SOILPIX *LocalSoil - Structure with soil moisture data and energy data
                         for the current pixel

  Returns      :
    void

  Modifies     :
    members of LocalSoil


  Comments     : This function sets all the energy fluxes at the pixel level 
                 to 0.0, but calculates the evapotranspiration in W/m2
*****************************************************************************/
void NoSensibleHeatFlux(int Dt, PIXMET *LocalMet, float ETot, 
                        SOILPIX *LocalSoil)
{
  LocalSoil->TSurf =  0.0;

  LocalSoil->Ra = 0.0;

  LocalSoil->Qnet = 0.0;

  LocalSoil->Qs = 0.0;
  
  LocalSoil->Qe = -(LocalMet->Lv * ETot)/Dt * WATER_DENSITY;
  
  LocalSoil->Qg = 0.0;
  
  LocalSoil->Qst = 0.0;
  
  LocalSoil->Qrest = 0.0;
}
