SUBROUTINE CalcEcoRoof(SurfNum,ZoneNum,ConstrNum,TempExt)
! SUBROUTINE INFORMATION
! AUTHOR David Sailor and Toan Pham
! DATE WRITTEN January 2007
! MODIFIED David Sailor - to fix initialization between DD runs and during warm-up
! RE-ENGINEERED na
! PURPOSE OF THIS MODULE:
! To calculate the heat balance for surfaces with eco roof specified as outside surface
! Note that only ONE ecoroof construction can be employed at present time. If multiple
! surfaces have ecoroof as the outside layer the energy balance is only calculated for
! the first such surface.
! METHODOLOGY EMPLOYED:
! Vikram Madhusudan's Portland State Univ. MS Thesis (Dec 2005) based on FASST model
! of Frankenstein and Koenig (2004) - DRDC/CRREL Technical Report TR-04-25.
! Some data used herein are from: European Centre for Medium-Range Weather Forecasts (ECMWF)
! IFS Documentation, CY25R1 (April 2002), www.ecmwf.int/research/ifsdocs/CY24r1/Physics/
! Physics-08-03.html.
! The Atmospheric Boundary Layer - by J.R. Garratt (Cambridge Atmos. & Space Science Series), 316pp.
!
USE DataGlobals
USE DataEnvironment
USE DataHeatBalFanSys
USE DataHeatBalance
USE DataHeatBalSurface
USE DataSurfaces
! USE DataDaylightingDevices
! USE DaylightingDevices, ONLY: FindTDDPipe
USE Psychrometrics
USE ConvectionCoefficients, ONLY : InitExteriorConvectionCoeff, SetExtConvectionCoeff, SetIntConvectionCoeff
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
!SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: SurfNum ! Indicator of Surface Number for the current surface
INTEGER, INTENT(IN) :: ZoneNum ! Indicator for zone number where the current surface
INTEGER, INTENT(INOUT) :: ConstrNum ! Indicator for contruction index for the current surface
REAL(r64) , INTENT(OUT) :: TempExt ! Exterior temperature boundary condidtion
!SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: Kv = 0.4d0 ! Von Karmen's constant (source FASST)
REAL(r64), PARAMETER :: rch = 0.63d0 ! Turbulent Schimdt Number
REAL(r64), PARAMETER :: Ks = 0.2d0 ! W/m.k. Thermal Conductivity of soil
REAL(r64), PARAMETER :: rche = 0.71d0 ! Turbulent Prandtl Number
REAL(r64), PARAMETER :: Rv = 461.53d0 ! Gas Constant of Water Vapor J/kg K
REAL(r64), PARAMETER :: Rair = 0.286d3 ! Gas Constant of air J/Kg K
REAL(r64), PARAMETER :: g1 = 9.81d0 ! Gravity. In m/sec^2.
REAL(r64), PARAMETER :: Sigma = 5.6697d-08 ! Stefan-Boltzmann constant W/m^2K^4
REAL(r64), PARAMETER :: Cpa = 1005.6d0 ! Specific heat of Water Vapor. (J/Kg.K)
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
INTEGER,SAVE :: FirstEcoSurf=0 ! Indicates lowest numbered surface that is an ecoroof
! used to determine WHEN to updatesoilProps...
INTEGER :: EcoLoop ! an integer loop variable for the simultaneous solution iteration
REAL(r64) :: AbsThermSurf ! Thermal absoptance of the exterior surface
INTEGER :: RoughSurf ! Roughness index of the exterior (ecoroof) surface.
REAL(r64) :: HMovInsul ! "Convection" coefficient of movable insulation
! REAL(r64) :: HSky ! "Convection" coefficient from sky to surface
! REAL(r64) :: HAir ! "Convection" coefficient from air to surface (radiation)
! INTEGER :: OPtr
! INTEGER :: OSCScheduleIndex ! Index number for OSC ConstTempSurfaceName
LOGICAL :: QuickConductionSurf = .FALSE. ! indicator for quick conduction surface
REAL(r64),SAVE :: LAI = 0.2d0 ! Leaf area index
REAL(r64),SAVE :: epsilonf = 0.95d0 ! Leaf Emisivity
REAL(r64),SAVE :: epsilong = 0.95d0 ! Soil Emisivity
REAL(r64),SAVE :: Alphag = 0.3d0 ! Ground Albedo
REAL(r64),SAVE :: Alphaf = 0.2d0 ! Leaf Albedo (reflectivity to solar radiation)
REAL(r64) :: e0 = 2.0d0 ! Windless lower limit of exchange coefficient (from FASST docs)
REAL(r64) :: RH =50.0d0 ! Relative humidity (%)
REAL(r64) :: Pa =101325.d0 ! Atmospheric Pressure (PA)
REAL(r64),SAVE :: Tg =10.0d0 ! Ground Surface temperature C ***** FROM PREVIOUS TIME STEP
REAL(r64),SAVE :: Tf =10.0d0 ! Leaf temperature C ***** FROM PREVIOUS TIME STEP
REAL(r64) :: Tgk ! Ground temperature in Kelvin
REAL(r64),SAVE :: Zf = 0.2d0 ! Height of plants (m)
! DJS Oct 2007 release - note I got rid of the initialization of moisture and meanrootmoisture here as these
! values are now set at beginning of each new DD and each new warm-up loop.
! DJS
REAL(r64),SAVE :: Moisture ! m^3/m^3.The moisture content in the soil is the value provided by a user
REAL(r64),SAVE :: MoistureResidual=0.05d0 ! m^3/m^3. Residual & maximum water contents are unique to each material.
! See Frankenstein et al (2004b) for data.
REAL(r64),SAVE :: MoistureMax =0.5d0 ! Maximum volumetric moisture content (porosity) m^3/m^3
REAL(r64),SAVE :: MeanRootMoisture ! Mean value of root moisture m^3/m^3
REAL(r64),SAVE :: SoilThickness =0.2d0 ! Soil thickness (m)
REAL(r64),SAVE :: StomatalResistanceMin ! s/m . ! Minimum stomatal resistance is unique for each veg. type.
REAL(r64) :: f3 =1.0d0 ! As the value of gd for tall grass is 0, then f3 = 1
! ECMWF 2002 CY25R1 report has gd=0.0 for all veg except trees where gd=0.03.
REAL(r64) :: Ta ! current air temperature
REAL(r64) :: Zog = 0.001d0 ! Ground roughness length scale (m)
REAL(r64) :: Za = 2.0d0 ! Instrument height where atmospheric wind speed is measured (m)
REAL(r64) :: Ws ! Wind Speed (m/s)
REAL(r64) :: Waf ! Windspeed within canopy (m/s)
REAL(r64) :: Latm ! Long Wave Radiation (W/m^2)
REAL(r64) :: qaf ! mixing ratio of air near canopy
REAL(r64) :: qg ! mixing ratio of air at surface.
REAL(r64),SAVE :: Lf ! latent heat flux
REAL(r64),SAVE :: Vfluxf = 0.0d0 ! Water evapotr. rate associated with latent heat from vegetation [m/s]
REAL(r64) :: RS ! shortwave radiation
REAL(r64) :: Qsoil = 0.0d0 ! heat flux from the soil layer
REAL(r64) :: EpsilonOne
!unused1208 REAL(r64) :: e
REAL(r64) :: eair
REAL(r64) :: Rhoa
REAL(r64) :: Tak
REAL(r64) :: qa ! mixing ratio of air
REAL(r64) :: Tafk
REAL(r64) :: Taf
REAL(r64) :: Rhof
REAL(r64) :: Rhoaf ! Average air density
REAL(r64) :: sigmaf
REAL(r64) :: Zd ! zero displacement height (m)
REAL(r64) :: Zo ! foliage roughness length (m)
REAL(r64) :: Cfhn ! transfer coefficient at near-neutral conditions
REAL(r64) :: Cf ! bulk Transfer coefficient, equation 10 page 6 (FASST).
REAL(r64), SAVE :: sheatf ! sensible heat flux coeff for foliage (W/m^2K)
REAL(r64), SAVE :: sensiblef ! sensible heat transfer TO foliage (W/m^2) DJS Jan 2011
REAL(r64) :: ra ! Aerodynamic Resistance
REAL(r64) :: f1inv ! intermediate calculation variable
REAL(r64) :: f2inv ! intermediate calculation variable
REAL(r64) :: f1 ! intermediate calculation variable
REAL(r64) :: f2 ! intermediate calculation variable
REAL(r64) :: r_s ! Minimum Stomatal Resistance, specific to each plant
REAL(r64) :: Mg ! Surface soil moisture content m^3/m^3 (Moisture / MoistureMax)
REAL(r64) :: dOne ! intermediate calculation variable
REAL(r64) :: esf ! the saturation vapor pressure (Pa)
REAL(r64) :: qsf ! Saturation specific humidity at leaf temperature (qfsat)
REAL(r64) :: Lef ! Latent heat of vaporation at leaf surface temperature (J/kg)
REAL(r64) :: Desf ! Derivative of Saturation vapor pressure
REAL(r64) :: dqf ! Derivative of saturation specific humidity
REAL(r64) :: dqg ! this is given by Clausius-Clapeyron equation
REAL(r64) :: esg ! Saturation vapor pressure (Pa)
REAL(r64) :: qsg ! Saturation specific humidity(mixing ratio?) at ground surface temperature
REAL(r64) :: Leg ! Latent heat vaporization at the ground temperature (J/kg)
REAL(r64) :: Desg ! derivative of esg Saturation vapor pressure(?)
REAL(r64) :: F1temp ! intermediate variable in computing flux from the soil layer
REAL(r64) :: P1 ! intermediate variable in the equation for Tf
REAL(r64) :: P2 ! intermediate variable in the equation for Tf and Tg
REAL(r64) :: P3 ! intermediate variable in the equation for Tf and Tg
REAL(r64) :: Rhog ! Density of air at the soil surface temperature
REAL(r64) :: Rhoag ! Average density of air with respect to ground surface and air temperature
REAL(r64) :: Rib ! Richardson Number
REAL(r64) :: Chng ! bulk transfer coefficient near ground
REAL(r64) :: Ce ! bulk transfer coefficient (this is in fact Ceg in equation 28 main report)
REAL(r64) :: Gammah ! latent heat exchange stability correction factor
REAL(r64) :: Chg ! in fact it is the same as Ce (=Ceg) is transfer coefficient (but wot?)
REAL(r64), SAVE :: sheatg ! intermediate calculation variable - sensible flux coef (W/m^2K for ground)
REAL(r64), SAVE :: sensibleg ! sensible heat flux TO ground (w/m^2) DJS Jan 2011
REAL(r64) :: T3G ! intermediate variable in the equation for Tg
REAL(r64) :: T2G ! intermediate variable in the equation for Tg
REAL(r64) :: LeafTK ! the current leaf's temperature (Kelvin)
REAL(r64) :: SoilTK ! the current soil's temperature (Kelvin)
REAL(r64) :: Chne ! is similar to near ground bulk transfer coefficient for latent heat flux (at neutral condition)
REAL(r64) :: Tif ! previous leaf temperature
REAL(r64) :: rn ! rn is the combined effect of both stomatal and aerodynamic resistances
! in fact this is called r'' in the main report
REAL(r64),SAVE :: Lg =0.0d0 ! latent heat flux from ground surface
REAL(r64),SAVE :: Vfluxg = 0.0d0 ! Water evapotr. rate associated with latent heat from ground surface [m/s]
REAL(r64) :: T1G ! intermediate variable in the equation for Tg
REAL(r64) :: Qsoilpart1 ! intermediate variable for evaluating Qsoil (part without the unknown)
REAL(r64) :: Qsoilpart2 ! intermediate variable for evaluating Qsoil (part coeff of the ground temperature)
! INTEGER,EXTERNAL :: getnewunitnumber ! external function to return a new (unique) unit for ecoroof writing
INTEGER :: unit=0
LOGICAL, SAVE :: MyEnvrnFlag=.true.
WS=WindspeedAt(Surface(SurfNum)%Centroid%Z) ! use windspeed at Z of roof
IF (WS < 2.0d0) THEN ! Later we need to adjust for building roof height...
Ws = 2.0d0 ! Set minimum possible wind speed outside vegetation to 2.0 m/s
! consistent with FASST TR-04-25 p. x (W' = 2.0)
END IF
IF(SurfaceWindow(SurfNum)%StormWinFlag==1) ConstrNum = Surface(SurfNum)%StormWinConstruction
RoughSurf = Material(Construct(ConstrNum)%LayerPoint(1))%Roughness
AbsThermSurf = Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpThermal
HMovInsul = 0.0d0
IF (Surface(SurfNum)%ExtWind) THEN
CALL InitExteriorConvectionCoeff(SurfNum,HMovInsul,RoughSurf,AbsThermSurf,TH(SurfNum,1,1), &
HcExtSurf(SurfNum),HSkyExtSurf(SurfNum),HGrdExtSurf(SurfNum),HAirExtSurf(SurfNum))
END IF
RS = BeamSolarRad + AnisoSkyMult(SurfNum)*DifSolarRad
Latm = 1.0d0 * sigma * 1.0d0* Surface(SurfNum)%ViewFactorGround*( (GroundTempKelvin)**4) &
+ 1.0d0 * sigma * 1.0d0*Surface(SurfNum)%ViewFactorSky*(SkyTempKelvin)**4
IF (EcoRoofbeginFlag) THEN
EcoRoofbeginFlag = .FALSE.
IF (Surface(SurfNum)%HeatTransferAlgorithm /= HeatTransferModel_CTF) &
CALL ShowWarningError('CalcEcoRoof: EcoRoof simulation but HeatBalanceAlgorithm is not ConductionTransferFunction(CTF).'// &
' Has not been tested under other solution approaches.')
! ONLY READ ECOROOF PROPERTIES IN THE FIRST TIME
Zf = Material(Construct(ConstrNum)%LayerPoint(1))%HeightofPlants ! Plant height (m)
LAI = Material(Construct(ConstrNum)%LayerPoint(1))%LAI ! Leaf Area Index
Alphag = 1.0d0 - Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpSolar ! albedo rather than absorptivity
Alphaf = Material(Construct(ConstrNum)%LayerPoint(1))%Lreflectivity ! Leaf Reflectivity
epsilonf = Material(Construct(ConstrNum)%LayerPoint(1))%Lemissitivity ! Leaf Emisivity
stomatalresistancemin = Material(Construct(ConstrNum)%LayerPoint(1))%Rstomata ! Leaf min stomatal resistance
epsilong = Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpThermal ! Soil Emisivity
MoistureMax = Material(Construct(ConstrNum)%LayerPoint(1))%porosity ! Max moisture content in soil
MoistureResidual = Material(Construct(ConstrNum)%LayerPoint(1))%MinMoisture ! Min moisture content in soil
Moisture = Material(Construct(ConstrNum)%LayerPoint(1))%InitMoisture ! Initial moisture content in soil
MeanRootMoisture = Moisture ! DJS Oct 2007 Release --> all soil at same initial moisture for Reverse DD fix
SoilThickness = Material(Construct(ConstrNum)%LayerPoint(1))%Thickness ! Total thickness of soil layer (m)
!DJS - This set of statements and the corresponding write statement in the UpdateSoilProps subroutine should
! be removed (or commented out) prior to deployment in a working version of EnergyPlus
!
!Open a unit for writing ecoroof specific data to output file (in EnergyPlus directory)
! unit=getnewunitnumber()
! open(unit,file='ecoroof.txt')
! write(unit,*) " ECOROOF OUTPUT REPORT TRACE - HOURLY "
! write(unit,*) " "
! write(unit,91)
! 91 FORMAT (" Day Hour Flux T_g T_f MoistTop MoistRoot CumRain CumET CumRunoff TotalIrr Dens SpecHeat Cond Albedo")
FirstEcoSurf = Surfnum ! this determines WHEN to updatesoilProps
! DJS NOVEMBER 2010 - Make calls to SetupOutput Variable to allow for reporting of ecoroof variables
CALL SetupOutputVariable('Green Roof Soil Temperature [C]',Tg, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Vegetation Temperature [C]',Tf, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Soil Root Moisture Ratio []',MeanRootMoisture, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Soil Near Surface Moisture Ratio []',Moisture, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Soil Sensible Heat Transfer Rate per Area [W/m2]',sensibleg, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Vegetation Sensible Heat Transfer Rate per Area [W/m2]',sensiblef, &
'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Vegetation Moisture Transfer Rate [m/s]',vfluxf, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Soil Moisture Transfer Rate [m/s]',vfluxg, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Vegetation Latent Heat Transfer Rate per Area [W/m2]',Lf, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Soil Latent Heat Transfer Rate per Area [W/m2]',Lg, 'Zone','State','Environment')
CALL SetupOutputVariable('Green Roof Cumulative Precipitation Depth [m]',CumPrecip, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Cumulative Irrigation Depth [m]',CumIrrigation, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Cumulative Runoff Depth [m]',CumRunoff, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Cumulative Evapotranspiration Depth [m]',CumET, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Current Precipitation Depth [m]',CurrentPrecipitation, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Current Irrigation Depth [m]',CurrentIrrigation, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Current Runoff Depth [m]',CurrentRunoff, 'Zone','Sum','Environment')
CALL SetupOutputVariable('Green Roof Current Evapotranspiration Depth [m]',CurrentET, 'Zone','Sum','Environment')
! DJS NOVEMBER 2010 - end of calls to setup output of ecoroof variables
END IF ! Initialization statements for first entry into ecoroof routines
! DJS July 2007
! Make sure the ecoroof module resets its conditions at start of EVERY warmup day and every new design day
! for Reverse DD testing
IF (BeginEnvrnFlag .or. WarmupFlag) THEN
Moisture = Material(Construct(ConstrNum)%LayerPoint(1))%InitMoisture ! Initial moisture content in soil
MeanRootMoisture= Moisture ! Start the root zone moisture at the same value as the surface.
Alphag = 1.0d0 - Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpSolar ! albedo rather than absorptivity
ENDIF
! DJS July 2007
IF (BeginEnvrnFlag .and. MyEnvrnFlag) THEN
Tgold = OutDryBulbTempAt(Surface(SurfNum)%Centroid%Z) !OutDrybulbTemp ! initial guess
Tfold = OutDryBulbTempAt(Surface(SurfNum)%Centroid%Z) !OutDrybulbTemp ! initial guess
Tg=10.0d0
Tf=10.0d0
VFluxf=0.0d0
VFluxg=0.0d0
CumRunoff = 0.0d0
CumET = 0.0d0
CumPrecip = 0.0d0
CumIrrigation = 0.0d0
CurrentRunoff = 0.0d0
CurrentET = 0.0d0
CurrentPrecipitation = 0.0d0
CurrentIrrigation = 0.0d0
MyEnvrnFlag=.false.
ENDIF
IF (.not. BeginEnvrnFlag) THEN
MyEnvrnFlag=.true.
ENDIF
! If current surface is = FirstEcoSurf then for this time step we need to update the soil moisture
If (Surfnum .eq. FirstEcoSurf) then
Call UpdateSoilProps(Moisture,MeanRootMoisture,MoistureMax,MoistureResidual,SoilThickness,Vfluxf,Vfluxg, &
ConstrNum,Alphag, unit,Tg,Tf,Qsoil)
Ta = OutDryBulbTempAt(Surface(SurfNum)%Centroid%Z) ! temperature outdoor - Surface is dry, use normal correlation
Tg = Tgold
Tf = Tfold
IF (Construct(ConstrNum)%CTFCross(0) > 0.01d0) THEN
QuickConductionSurf = .TRUE.
F1temp = Construct(ConstrNum)%CTFCross(0) &
/( Construct(ConstrNum)%CTFInside(0)+HConvIn(SurfNum) )
Qsoilpart1 = -CTFConstOutPart(SurfNum) &
+F1temp*( CTFConstInPart(SurfNum) &
+QRadSWInAbs(SurfNum) &
+QRadThermInAbs(SurfNum) &
+Construct(ConstrNum)%CTFSourceIn(0)*QsrcHist(SurfNum,1) &
+HConvIn(SurfNum)*MAT(ZoneNum) &
+NetLWRadToSurf(SurfNum) )
ELSE
Qsoilpart1 = -CTFConstOutPart(SurfNum) &
+ Construct(ConstrNum)%CTFCross(0)*TempSurfIn(SurfNum)
F1temp = 0.0d0
END IF
Qsoilpart2 = Construct(ConstrNum)%CTFOutside(0) - F1temp*Construct(ConstrNum)%CTFCross(0)
Pa = StdBaroPress ! standard atmospheric pressure (apparently in Pascals)
Tgk = Tg + KelvinConv
Tak = Ta + KelvinConv
sigmaf = 0.9d0 - 0.7d0*exp(-0.75d0*LAI) ! Fractional veg cover based on (2) from FASST TR-04-25
! Formula for grasses modified to incorporate limits from
! Table 1 for sigmaf_max and min (0.20 to 0.9)
EpsilonOne = epsilonf + epsilong - epsilong * epsilonf ! Checked (eqn. 6 in FASST Veg Models)
RH= OutRelHum ! Get humidity in % from the DataEnvironment.f90
eair = (RH/100.0d0)*611.2d0*EXP( 17.67d0*Ta/(Tak-29.65d0) )
qa = ( 0.622d0*eair)/( Pa-1.000d0*eair ) ! Mixing Ratio of air
Rhoa = Pa/(Rair*Tak) ! Density of air. kg/m^3
Tif = Tf
! Air Temperature within the canopy is given as
! (Deardorff (1987)). Kelvin. based of the previous temperatures
Tafk = (1.d0-sigmaf)*Tak + &
sigmaf*( 0.3d0*Tak+0.6d0*(Tif+KelvinConv)+0.1d0*Tgk)
Taf = Tafk - KelvinConv ! Air Temperature within canopy in Celcius (C).
Rhof= Pa / (Rair*Tafk) ! Density of air at the leaf temperature
Rhoaf = (Rhoa + Rhof )/2.0d0 ! Average of air density
Zd = 0.701d0 * Zf**0.979d0 ! Zero displacement height
Zo = 0.131d0 * Zf**0.997d0 ! Foliage roughness length. (m) Source Ballick (1981)
if (zo .lt. 0.02d0) zo = 0.02d0 ! limit based on p.7 TR-04-25 and Table 2
!transfer coefficient at near-neutral condition Cfhn
Cfhn = (Kv/LOG((Za - Zd)/Zo))**2 !Equation 12, page 7, FASST Model
Waf = 0.83d0*Cfhn**0.5d0*sigmaf*Ws + (1.d0- sigmaf)*Ws ! Wind Speed inside foliage. Equation #6, FASST model
Cf = 0.01d0*(1.0d0+0.3d0/Waf) ! The bulk Transfer coefficient, equation 10 page 6.
sheatf = e0+1.1d0*LAI*Rhoaf*Cpa*Cf*Waf ! Intermediate calculation for Sensible Heat Transfer
sensiblef = sheatf*(Taf - Tf) ! DJS Jan 2011 sensible flux TO foliage into air (Frankenstein 2004, eqn7)
!sourced from Frankenstein et al (2004a). Added e0 windless correction factor.
!an early version had (1.0-0.7)*sigmaf in calc of sensiblef... how did that get there!?! Fixed.
!These parameters were taken from "The Atm Boundary Layer", By J.R. Garratt
!NOTE the Garratt eqn. (A21) gives esf in units of hPA so we have multiplied
!the constant 6.112 by a factor of 100.
esf = 611.2d0 * EXP(17.67d0*Tif/(Tif+KelvinConv-29.65d0) )
! From Garratt - eqn. A21, p284. Note that Tif and Tif+KelvinConv usage is correct.
! Saturation specific humidity at leaf temperature again based on previous temperatures
qsf = 0.622d0 * esf/(Pa-1.000d0*esf) ! "The Atm Boundary Layer", J.R Garrat for Saturation mixing ratio
! Calculate stomatal resistance and atmospheric resistance Part
ra = 1.0d0 / (Cf * Waf) ! Aerodynamic Resistance. Resistance that is caused
! by the boundary layer on a leaf surface to transfer water vapor. It is measured in
! s/m and depends on wind speed, leaf's surface roughness,
! and stability of atsmophere.
f1inv = MIN(1.0d0, (0.004d0*RS + 0.005d0)/(0.81d0 * (0.004d0*RS +1.0d0))) ! SW radiation-related term
f1= 1.d0/f1inv
if (MoistureMax .eq. MoistureResidual) then
f2inv = 1.0d10
else
f2inv = (MeanRootMoisture - MoistureResidual)/(MoistureMax - MoistureResidual) !Equation 19 p. 9 FASST model
endif
!In FASST, Eq 20 is used to compute Moisture.
f2 = 1.d0/f2inv ! In dry areas f2 --> LARGE so that r_s --> LARGE
! and both rn and Latent flux --> 0.0
f3 = 1.d0/( EXP(-0.0d0*(esf-eair)) ) ! Note the 0.0 here is gd which has value of 0.0
! for most plants and 0.03 for trees (see ECMWF
! note given above.
r_s = StomatalResistanceMin*f1*f2*f3/LAI ! Stomatal Resistance (r_s)
rn = ra/(ra+r_s) !rn is foliage surface wetness ... NOT a resistance
! This routine is to calculate ground moisture factor. This factor is from *****
Mg = moisture / MoistureMax ! m^3/m^3.
dOne = 1.0d0- sigmaf *(0.6d0*(1.0d0-rn)+0.1d0*(1.0d0-Mg))
!Latent heat of vaporation at leaf surface temperature. The source of this
!equation is Henderson-Sellers (1984)
Lef = 1.91846d6 * ( (Tif+KelvinConv)/(Tif+KelvinConv-33.91d0) )**2
!Check to see if ice is sublimating or frost is forming.
if (tfold .lt. 0.0d0 ) Lef= 2.838d6 ! per FASST documentation p.15 after eqn. 37.
!Derivative of Saturation vapor pressure, which is used in the calculation of
!derivative of saturation specific humidity.
desf = 611.2d0*exp(17.67d0*(Tf/(Tf+KelvinConv-29.65d0))) * &
( 17.67d0*Tf*(-1.0d0)*(Tf+KelvinConv-29.65d0)**(-2) + 17.67d0/(KelvinConv-29.65d0 + Tf))
dqf = ((0.622d0* Pa)/(Pa - esf)**2)*Desf !Derivative of saturation specific humidity
esg = 611.2d0*EXP( 17.67d0*( Tg/((Tg+KelvinConv)-29.65d0))) !Pa saturation vapor pressure
! From Garratt - eqn. A21, p284.
! Note that Tg and Tg+KelvinConv usage is correct.
qsg = 0.622d0*esg/(Pa-esg) !Saturation mixing ratio at ground surface temperature.
!Latent heat vaporization at the ground temperature
Leg = 1.91846d6*(Tgk/(Tgk-33.91d0))**2
!Check to see if ice is sublimating or frost is forming.
if (tgold .lt. 0.0d0 ) Leg= 2.838d6 ! per FASST documentation p.15 after eqn. 37.
desg = 611.2d0*exp(17.67d0*(Tg/(Tg+KelvinConv-29.65d0))) * &
(17.67d0*Tg*(-1.0d0)*(Tg+KelvinConv-29.65d0)**(-2) + 17.67d0/(KelvinConv-29.65d0 + Tg))
dqg = (0.622d0*Pa/(Pa-esg)**2)*Desg
!Final Ground Atmosphere Energy Balance
!Density of air at the soil surface temperature
Rhog = Pa / (Rair*Tgk)
!Average density of air with respect to ground surface and air temperature
Rhoag = (Rhoa + Rhog)/2.d0
Rib = 2.0d0*g1*Za*(Taf-Tg) / ( (Tafk+Tgk)*Waf**2 ) !Richardson Number
! Compute the stability factor Gammah
IF (Rib < 0.0d0) THEN
Gammah = (1.0d0-16.0d0*Rib)**(-0.5d0)
ELSE
IF (Rib >= 0.19d0) THEN
Rib = .19d0
ENDIF
Gammah = (1.0d0-5.0d0*Rib)**(-0.5d0)
ENDIF
if (roughsurf == VerySmooth) then ! 6= very smooth, 5=smooth, 4= med. sm. ,3= med. rough. , 2= rough, 1= Very rough
zog = 0.0008d0
else if (roughsurf == Smooth) then
zog = 0.0010d0
else if (roughsurf == MediumSmooth) then
zog = 0.0015d0
else if (roughsurf ==MediumRough) then
zog = 0.0020d0
else if (roughsurf ==Rough) then
zog = 0.0030d0
else ! VeryRough
zog = 0.005d0
endif
Chng = (Kv/LOG(Za/Zog))**2 / rch ! bulk transfer coefficient near ground
Chg = Gammah* ( (1.d0-sigmaf)*Chng+sigmaf*Cfhn ) !
sheatg= e0+Rhoag*Cpa*Chg*Waf ! added the e0 windless correction
sensibleg = sheatg*(Taf - Tg) ! sensible flux TO soil (W/m^2) DJS Jan 2011 (eqn. 32 in Frankenstein 2004)
Chne = (Kv/LOG(Za/Zog))**2 / rche
Ce = gammah * ((1.d0-sigmaf)*Chne+sigmaf*Cfhn) ! this is in fact Ceg in eq (28)
!we can approximate Gammae by Gammah (Supported by FASST Veg Models p. 15)
qaf = ( (1.0d0-sigmaf)*qa+sigmaf*(0.3d0*qa+0.6d0*qsf*rn+0.1d0*qsg*Mg) )/&
( 1.0d0-sigmaf*(0.6d0*(1.d0-rn)+0.1d0*(1.0d0-Mg)) )
qg = Mg*qsg + (1.0d0-Mg)*qaf !eq main report (13)
! According to FASST documentation this is correct.
! The following also used Rhoaf and Rhoag respectively... shouldn't these be water densities??!!
Lf= Lef*LAI*Rhoaf*Cf*Waf*rn*(qaf - qsf) ! This had Leg but should be Lef...
Lg = Ce*Leg*Waf*Rhoag*(qaf-qg) *Mg ! In the FASST documentation there is NO Mg. However, in looking
! back at the Deardorff 1978 paper it appears that an alpha = Mg term is
! used to distinguish from POTENTIAL and ACTUAL ground surface evaporation...
! the Lf and Lg calculations are NOT used in this formulation
! rather, the effects are included in terms of dqg and dqf !
! These equations for Lf and Lg are based on Deardorff's paper, but there is
! a sign difference !!! (qsf -qaf) and (qg - qaf) ?!
! These may be useful, however for telling the UpdateSoilProps routine
! how much evaporation has taken place...
Vfluxf= -1.0d0*Lf/Lef/990.0d0 ! water evapotranspire rate [m/s]
Vfluxg= -1.0d0*Lg/Leg/990.0d0 ! water evapotranspire rate [m/s]
if (Vfluxf .LT. 0.0d0) Vfluxf = 0.0d0 ! According to FASST Veg. Models p. 11, eqn 26-27, if Qfsat > qaf the actual
if (Vfluxg .LT. 0.0d0) Vfluxg = 0.0d0 ! evaporative fluxes should be set to zero (delta_c = 1 or 0).
! P1, P2, P3 corespond to first, second and third terms of equation 37 in the main report.
! Note: the FASST model has a term -gamma_p*(1.0-exp...) in first line for P1 (c1_f) where gamma_p is
! a precipitation variable. So, if we assume no precip this term vanishes. We should
! revisit this issue later.
! Implement an iterative solution scheme to solve the simultaneous equations for Leaf and Soil temperature.
! Prior experience suggests that no more than 3 iterations are likely needed
!
LeafTK = Tf+KelvinConv
SoilTK = Tg+KelvinConv
do ecoloop=1,3
P1 = sigmaf *( Rs * (1-Alphaf) + epsilonf * Latm ) &
- 3.d0*sigmaf*epsilonf*epsilong*sigma*(SoilTK)**4/EpsilonOne &
- 3.d0*( -sigmaf*epsilonf*sigma - sigmaf*epsilonf*epsilong*sigma/EpsilonOne )*(LeafTK)**4 &
+ sheatf*(1.d0-0.7d0*sigmaf)*(Ta+KelvinConv) + LAI * &
Rhoaf * Cf * Lef * Waf * rn * ((1.d0 - 0.7d0*sigmaf)/done)*qa + LAI*Rhoaf*Cf*Lef &
* Waf * rn *(((0.6d0 * sigmaf*rn) / done ) - 1.d0) *( qsf-(LeafTK)*dqf)+LAI &
* Rhoaf * Cf * Lef * Waf * rn * ((0.1d0*sigmaf*Mg)/done)*(qsg-(SoilTK)* dqg)
P2 = 4.0d0*(sigmaf*epsilonf*epsilong*sigma)*(SoilTK)**(3)/EpsilonOne + 0.1d0*sigmaf*sheatf+LAI*Rhoaf*Cf*Lef &
*Waf*rn*(0.1d0*sigmaf*Mg)/done *dqg
P3 = 4.0d0* (-sigmaf*epsilonf*sigma - (sigmaf*epsilonf*sigma*epsilong)/EpsilonOne) * (LeafTK)**3 + (0.6d0 * sigmaf - 1.d0) &
* sheatf + LAI * Rhoaf * Cf * Lef * Waf * rn * (((0.6d0 * sigmaf * rn) / done)-1.0d0) * dqf
!T1G, T2G, & T3G corresponds to first, second & third terms of equation 38
!in the main report.
! as with the equations for vegetation the first term in the ground eqn in FASST has a
! term starting with gamma_p --- if no precip this vanishes. Again, revisit this issue later.
T1G = (1.d0- sigmaf) * (Rs * (1.d0- Alphag) + epsilong * Latm)- &
(3.0d0 * (sigmaf * epsilonf * epsilong * sigma) / EpsilonOne) * &
(LeafTK)**4 - 3.0d0*(-(1.0d0-sigmaf) * epsilong * sigma &
- sigmaf*epsilonf*epsilong*sigma/EpsilonOne)*(SoilTK)**4 &
+ sheatg*(1.d0-0.7d0*sigmaf)*(Ta+KelvinConv)+ Rhoag*Ce*Leg*Waf*Mg*((1.d0 - 0.7d0*sigmaf)/dOne)*qa & !finished by T1G
+ Rhoag*Ce*Leg*Waf*Mg*(0.1d0*sigmaf*Mg/dOne - Mg)*(qsg-(SoilTK)*dqg) &
+ Rhoag*Ce*Leg*Waf*Mg*(0.6d0*sigmaf*rn/dOne)*(qsf-(LeafTK)*dqf) + Qsoilpart1 + Qsoilpart2*(KelvinConv)
T2G = 4.0d0*(-(1.0d0-sigmaf)*epsilong*sigma - sigmaf*epsilonf*epsilong*sigma/EpsilonOne)*(SoilTK)**3 &
+ (0.1d0*sigmaf-1.0d0)*sheatg + Rhoag*Ce*Leg*Waf*Mg*( 0.1d0*sigmaf*Mg/dOne - Mg )*dqg - Qsoilpart2
T3G = ( 4.0d0*( sigmaf*epsilong*epsilonf*sigma )/EpsilonOne )*(LeafTK)**3 + 0.6d0*sigmaf*sheatg &
+ Rhoag*Ce*Leg*Waf*Mg*( 0.6d0*sigmaf*rn/dOne ) * dqf
LeafTK = 0.5d0* (LeafTK+( P1*T2G - P2*T1G )/( - P3*T2G + T3G*P2 ) ) ! take avg of old and new each iteration
SoilTK = 0.5d0* (SoilTK+( P1*T3G - P3*T1G )/( - P2*T3G + P3*T2G ) ) ! take avg of old and new each iteration
! in earlier implementations we simply assigned new Leaf and Soil temps once (no iteration -- ecoloop =1,1, but
! with LeafTK = ( P1*T2G - P2*T1G )/( - P3*T2G + T3G*P2 ) and
! SoilTK =(SoilTK+( P1*T3G - P3*T1G )/( - P2*T3G + P3*T2G ). This iterative solution was initially attempted to
! deal with stability issues associated with the CTF. It has virtually no impact on temperatures and it is not
! clear how much it helped with stability problems. Eventually when CTF solution is replaced with a finite
! difference scheme this loop structure should be removed.
end do ! This loop does an iterative solution of the simultaneous equations
Qsoil = -1.0d0*(Qsoilpart1-Qsoilpart2*(SoilTK-KelvinConv)) ! This is heat flux INTO top of the soil
Tfold = LeafTK - KelvinConv
Tgold = SoilTK - KelvinConv
Endif ! if firstecosurface (if not we do NOT need to recalculate ecoroof energybalance as all ecoroof surfaces MUST be the same
! this endif was moved here from the if statement regarding whether we are looking at the first ecoroof surface or not.
TH(SurfNum,1,1) = Tgold ! SoilTemperature
TempExt = Tgold
END SUBROUTINE CalcEcoRoof