Nodes of different colours represent the following:
Solid arrows point from a parent (sub)module to the submodule which is descended from it. Dashed arrows point from a module being used to the module or program unit using it. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=r64), | intent(inout) | :: | Moisture | |||
| real(kind=r64), | intent(inout) | :: | MeanRootMoisture | |||
| real(kind=r64), | intent(in) | :: | MoistureMax | |||
| real(kind=r64), | intent(in) | :: | MoistureResidual | |||
| real(kind=r64), | intent(in) | :: | SoilThickness | |||
| real(kind=r64), | intent(in) | :: | Vfluxf | |||
| real(kind=r64), | intent(in) | :: | Vfluxg | |||
| integer, | intent(inout) | :: | ConstrNum | |||
| real(kind=r64), | intent(inout) | :: | Alphag | |||
| integer, | intent(in) | :: | unit | |||
| real(kind=r64), | intent(in) | :: | Tg | |||
| real(kind=r64), | intent(in) | :: | Tf | |||
| real(kind=r64), | intent(in) | :: | Qsoil |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
SUBROUTINE UpdateSoilProps(Moisture,MeanRootMoisture,MoistureMax,MoistureResidual,SoilThickness,Vfluxf,Vfluxg, &
ConstrNum,Alphag, unit,Tg,Tf,Qsoil)
! SUBROUTINE INFORMATION
! AUTHOR David Sailor
! DATE WRITTEN Jan 2007
! MODIFIED Stephen Forner, Portland State University (SF); 7/15/2010
! RE-ENGINEERED na
! PURPOSE OF THIS MODULE:
! Track moisture input/output to ecoroof soil media (precipitation, irrigation, evapotranspiration, runoff)
! Update soil thermal properties associated with variations in soil moisture and update CTF calculations
! for the ecoroof construction layer.
! METHODOLOGY EMPLOYED:
! Define 2 soil layers (top and root). Moisture redistribution between top and root layers is proportional
! to moisture differences. Input and Output of moisture are partitioned between layers.
! Soil thermal properties vary based on non-dimensionalization of experimental data for 8 typical soils.
! Specifically, THERMAL PROPERTY = Dry Value + (fraction of moisture content)*Wet Value
USE DataGlobals
USE DataInterfaces
USE DataEnvironment
USE DataSurfaces
USE DataWater, ONLY: RainFall
USE General, ONLY: RoundSigDigits
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
REAL(r64), INTENT(INOUT) :: Moisture
REAL(r64), INTENT(INOUT):: MeanRootMoisture
REAL(r64), INTENT(IN) :: MoistureMax
REAL(r64), INTENT(IN) :: MoistureResidual
REAL(r64), INTENT(IN) :: SoilThickness
REAL(r64), INTENT(IN) :: Vfluxf ! Water mass flux from vegetation [m/s]
REAL(r64), INTENT(IN) :: Vfluxg ! Water mass flux from soil surface [m/s]
INTEGER, INTENT(INOUT) :: ConstrNum ! Indicator for contruction index for the current surface
REAL(r64), INTENT(INOUT) :: Alphag
INTEGER, INTENT(IN) :: unit !unused1208
REAL(r64), INTENT(IN) :: Tg !unused1208
REAL(r64), INTENT(IN) :: Tf !unused1208
REAL(r64), INTENT(IN) :: Qsoil !unused1208
! SUBROUTINE PARAMETER DEFINITIONS:
! na
!Soil Parameters from Reference listed in the code:
REAL(r64), PARAMETER :: alpha=23.0d0 !These parameters are empirical constants
REAL(r64), PARAMETER :: n=1.27d0 !These parameters are empirical constants
REAL(r64), PARAMETER :: lambda=0.5d0 !These parameters are empirical constants
!This is another parameter of the soil which describes the soil conductivity at the saturation point (m/s)
REAL(r64), PARAMETER ::SoilConductivitySaturation=5.157d-7
! INTERFACE BLOCK SPECIFICATIONS:
! na
! DERIVED TYPE DEFINITIONS:
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
REAL(r64) :: RatioMax
REAL(r64) :: RatioMin
REAL(r64) :: MoistureDiffusion ! moisture transport down from near-surface to root zone
REAL(r64),SAVE :: TopDepth ! Thickness of "near-surface" soil layer
REAL(r64),SAVE :: RootDepth ! Thickness of "root zone" soil layer
! Note TopDepth+RootDepth = thickness of ecoroof soil layer
REAL(r64),SAVE :: SecondsPerTimeStep ! Seconds per TimeStep
REAL(r64) :: SoilConductivity ! Moisture dependent conductivity to be fed back into CTF Calculator
REAL(r64) :: SoilSpecHeat ! Moisture dependent Spec. Heat to be fed back into CTF Calculator
REAL(r64) :: SoilAbsorpSolar ! Moisture dependent Solar absorptance (1-albedo)
REAL(r64) :: SoilDensity ! Moisture dependent density to be fed back into CTF Calculator
REAL(r64) :: SatRatio
REAL(r64) :: TestRatio ! Ratio to determine if timestep change in properties is too abrupt for CTF
REAL(r64),SAVE :: DryCond ! Dry soil value of conductivity
REAL(r64),SAVE :: DryDens ! Dry soil value of density
REAL(r64),SAVE :: DryAbsorp ! Dry soil value of solar absorptance (1-albedo)
REAL(r64),SAVE :: DrySpecHeat ! Dry soil value of specific heat
REAL(r64) :: AvgMoisture ! Average soil moisture over depth of ecoroof media
LOGICAL,SAVE :: UpdatebeginFlag = .TRUE. ! one time flag
REAL(r64), SAVE ::CapillaryPotentialTop=-3.8997d0 !This variable keeps track of the capillary potential of the soil in both layers and time (m)
REAL(r64), SAVE ::CapillaryPotentialRoot=-3.8997d0 !This variable keeps track of the capillary potential of the soil in both layers and time (m)
REAL(r64), SAVE ::SoilHydroConductivityTop=8.72d-6 !This is the soil water conductivity in the soil (m/s)
REAL(r64), SAVE ::SoilHydroConductivityRoot=8.72d-6 !This is the soil water conductivity in the soil (m/s)
REAL(r64), SAVE ::SoilConductivityAveTop=8.72d-6 !This is the average soil water conductivity (m/s)
REAL(r64), SAVE ::SoilConductivityAveRoot=8.72d-6
REAL(r64), SAVE ::RelativeSoilSaturationTop !Relative Soil Saturation (soil moisture-residual soil moisture)/(saturation soil moisture-residual soil moisture)
REAL(r64), SAVE ::RelativeSoilSaturationRoot
REAL(r64) ::TestMoisture=0.15d0 !This makes sure that the moisture cannot change by too much in each step
INTEGER :: index1
INTEGER :: ErrIndex=0
! NOTE: As Energyplus calls the energy balance manager (and hence CalcEcoroof)
! once for each surface within each zone that has an ecoroof
! --- the CALCECOROOF routine is called multiple times within each time step
! So, it is important that the UpdateSoilProps subroutine is called ONLY ONCE for each time step!!!
! Recall Moisture = near-surface moisture value (m^3/m^3)
! Recall MeanRootMoisture = root zone moisture value (m^3/m^3)
!DJS 2009 set the ratiomax and ratiomin values in the code now (rather than as parameters) so that
!DJS 2009 we can link them to timesteps and make these limits apply to actual RATES...
!DJS 2009 reasonable rates are +/- 10% change in properties per 15 minute period... Otherwise we have
!DJS 2009 stability issues.
!DJS 2011 FEB - Since we no longer use CTF with soil-dependent properties (Do not RECALL INITCONDUCTION...
!DJS 2011 FEB - we may be able to get away with NO limits on rates of change when using CFD routine.
!DJS 2011 FEB - for now we stick with 20% per quarter hour.
RatioMax = 1.0d0 + 0.20d0*minutespertimestep/15.0d0
RatioMin = 1.0d0 - 0.20d0*minutespertimestep/15.0d0
If (UpdatebeginFlag) then
! SET dry values that NEVER CHANGE
DryCond =Material(Construct(ConstrNum)%LayerPoint(1))%Conductivity
DryDens =Material(Construct(ConstrNum)%LayerPoint(1))%Density
DryAbsorp =Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpSolar
DrySpecHeat =Material(Construct(ConstrNum)%LayerPoint(1))%SpecHeat
! DETERMINE RELATIVE THICKNESS OF TWO LAYERS OF SOIL (also unchanging)
If (SoilThickness .gt. 0.12d0) then
TopDepth = 0.06d0 ! For now use 6cm as top depth - everything else is root layer
else
TopDepth = 0.5d0*SoilThickness ! In unusual case of very thin soil make topdepth half of total
endif
!This loop outputs the minimum number of time steps needed to keep the solution stable
!The equation is minimum timestep in seconds=161240*((number of layers)**(-2.3))*(Total thickness of the soil)**2.07
IF (Material(Construct(ConstrNum)%LayerPoint(1))%EcoRoofCalculationMethod == 2) THEN
DO index1=1,20,1
IF ((REAL(MinutesPerTimeStep/index1,r64)) .LE. (161240.d0*2.d0**(-2.3d0)*(TopDepth+RootDepth)**(2.07d0))/60.d0) EXIT
END DO
IF (index1 .GT. 1) THEN
CALL ShowSevereError('CalcEcoRoof: Too few time steps per hour for stability.')
CALL ShowContinueError('...Entered Timesteps per hour=['//trim(RoundSigDigits(NumOfTimeStepInHour))// &
'], Change to some value greater than ['//trim(RoundSigDigits(60*index1/MinutesPerTimeStep))//']'// &
' for assured stability.')
! CALL ShowFatalError('Program terminates due to previous condition.')
END IF
ENDIF
RootDepth = SoilThickness-TopDepth
!Next create a timestep in seconds
SecondsPerTimeStep=MinutesPerTimeStep*60.d0
UpdatebeginFlag = .FALSE.
endif
CurrentRunoff = 0.0d0 ! Initialize current time step runoff as it is used in several spots below...
! FIRST Subtract water evaporated by plants and at soil surface
Moisture = Moisture - (Vfluxg)*MinutesPerTimeStep*60.0d0/TopDepth ! soil surface evaporation
MeanRootMoisture = MeanRootMoisture - (Vfluxf)*MinutesPerTimeStep*60.0d0/RootDepth ! plant extraction from root zone
! NEXT Update evapotranspiration summary variable for print out
CurrentET= (Vfluxg+ Vfluxf)*MinutesPerTimeStep*60.0d0 ! units are meters
IF (.not. WarmupFlag) THEN
CumET = CumET + CurrentET
ENDIF
! NEXT Add Precipitation to surface soil moisture variable (if a schedule exists)
IF (.not. WarmupFlag) THEN
CurrentPrecipitation = 0.0d0 ! first initialize to zero
ENDIF
CurrentPrecipitation = 0.0d0 ! first initialize to zero
If (Rainfall%ModeID ==RainSchedDesign) then
CurrentPrecipitation = Rainfall%CurrentAmount ! units of m
Moisture = Moisture + CurrentPrecipitation/TopDepth ! x (m) evenly put into top layer
IF (.not. WarmupFlag) THEN
CumPrecip = CumPrecip + CurrentPrecipitation
ENDIF
Endif
! NEXT Add Irrigation to surface soil moisture variable (if a schedule exists)
CurrentIrrigation = 0.0d0 ! first initialize to zero
irrigation%actualamount=0.0d0
If (Irrigation%ModeID == IrrSchedDesign) then
CurrentIrrigation = irrigation%ScheduledAmount ! units of m
irrigation%actualamount=CurrentIrrigation
! elseif (Irrigation%ModeID ==IrrSmartSched .and. moisture .lt. 0.4d0*MoistureMax) then
elseif (Irrigation%ModeID ==IrrSmartSched .and. moisture .lt. Irrigation%IrrigationThreshold*MoistureMax) then
! Smart schedule only irrigates when scheduled AND the soil is less than 40% saturated
CurrentIrrigation = irrigation%ScheduledAmount ! units of m
irrigation%actualamount=CurrentIrrigation
endif
Moisture = Moisture + CurrentIrrigation/TopDepth ! irrigation in (m)/timestep put into top layer
IF (.not. WarmupFlag) THEN
CumIrrigation = CumIrrigation + CurrentIrrigation
ENDIF
! Note: If soil top layer gets a massive influx of rain &/or irrigation some of
! the water will simply run right off the top and not penetrate at all!
! At the present time this limit is fairly small due to some minor stability issues
! in EnergyPlus. If the moisture changes too rapidly the code cannot handle the rapid changes in
! surface characteristics and heat fluxes. The result that I've noticed is a non-physical fluctation
! in ground surface temperature that oscillates up to 10 deg C from one hour to the next until the
! code catches up. The temporary solution is to simply limit how much moisture can enter the soil
! in any time step to 0.5"/hour. In the future this might be fixed by running with finer time steps
! or by using a finite difference temperature solution rather than the CTF.
! I suspect that 15 minute intervals may be needed. Another option is to have an internal moisture
! overflow bin that will hold extra moisture and then distribute it in subsequent hours. This way the
! soil still gets the same total moisture... it is just distributed over a longer period.
if (CurrentIrrigation + CurrentPrecipitation .gt. 0.5d0*0.0254d0*MinutesPerTimeStep/60.0d0) then
CurrentRunoff = CurrentIrrigation+CurrentPrecipitation - (0.5d0*0.0254d0*MinutesPerTimeStep/60.0d0)
! If we get here then TOO much moisture has already been added to soil (must now subtract excess)
Moisture = Moisture - CurrentRunoff/TopDepth ! currently any incident moisture in excess of 1/4 " per hour
! simply runs off the top of the soil.
endif
! Now, if top layer is beyond saturation... the excess simply runs off without penetrating into the lower
! layers.
if (Moisture .gt. MoistureMax) then
CurrentRunoff = CurrentRunoff + (Moisture-MoistureMax)*TopDepth
Moisture = MoistureMax
endif
IF (Material(Construct(ConstrNum)%LayerPoint(1))%EcoRoofCalculationMethod == 1) THEN
!THE SECTION BELOW WAS THE INITIAL MOISTURE DISTRIBUTION MODEL.
!Any line with "!-" was code. A line with "!" was just a comment. This is done in case this code needs to be resurected in the future.
!See below this commented out code for the new moisture distribution model.
!*********************************************************************************************************
!*********************************************************************************************************
! NEXT Redistribute moisture based on moisture diffusion.
! The effective diffusivities should be revisted when better moisture transport data in ecoroof soils are
! available.
! Here the diffusion rate is in units of [1/s]
! A value of 0.0001 would be ~ 36% / hour
! A value of 0.00005 would be ~ 18%/hour change in moisture
! A value of 0.000025 gives about a 9%/hour change in moisture
! a value of 0.0000125 gives 5%/hour...
! Note: This formulation allows moisture to have a directional bias - ie, it can sink into soil easier than
! it can be brought up.
If (Moisture .GT. MeanRootMoisture) then
! move moisture from top layer into root zone
MoistureDiffusion = min((MoistureMax-MeanRootMoisture)*RootDepth, (Moisture-MeanRootMoisture)*TopDepth)
MoistureDiffusion = max(0.0d0,MoistureDiffusion) ! Safety net to keep positive (not needed?)
! at this point moistureDiffusion is in units of (m)/timestep
MoistureDiffusion= 0.00005d0*MinutesPerTimeStep*60.0d0*MoistureDiffusion
Moisture = Moisture - MoistureDiffusion/TopDepth
MeanRootMoisture = MeanRootMoisture + MoistureDiffusion/RootDepth
else if (MeanRootMoisture .GT. Moisture) then
! move moisture to top layer from root zone
MoistureDiffusion = min((MoistureMax-Moisture)*TopDepth, (MeanRootMoisture-Moisture)*RootDepth)
MoistureDiffusion = max(0.0d0,MoistureDiffusion) ! Safety net (not needed?)
! at this point moistureDiffusion is in units of (m)/timestep
MoistureDiffusion= 0.00001d0*MinutesPerTimeStep*60.0d0*MoistureDiffusion
Moisture = Moisture + MoistureDiffusion/TopDepth
MeanRootMoisture = MeanRootMoisture - MoistureDiffusion/RootDepth
endif
ELSE
!********************************************************************************************************
!********************************************************************************************************
!THE SECTION ABOVE WAS THE MOISTURE DISTRIBUTION MODEL. REPLACED SF-7/21/2010
!NEXT redistribute the moisture in the soil based on:
!Marcel G Schaap and Martinus Th van Genuchten, 2006, 'A modified Maulem-van
!Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation.
!Written in MATLAB by Vishal Sharma (of Portland State) and modified for FORTRAN by Stephen Forner Summer 2010
!This model is based on curve fit data that describes the capillary motion of the water in combination with the gravitational
!forces on the water.
!This set of equations is unstable if the time step is larger than given by a curve fit equation. This first DO loop is to
!see how many time steps are needed to be stable.
!This method of moisture distribution relies on variables which are experimentally determined: alpha, lambda, n and the
!hydraulic conductivity at saturation.
!Now, solve for the soil parameters for of the top soil layer
RelativeSoilSaturationTop=(Moisture-MoistureResidual)/(MoistureMax-MoistureResidual)
IF (RelativeSoilSaturationTop < 0.0001d0) THEN
IF (ErrIndex == 0) THEN
CALL ShowWarningMessage('EcoRoof: UpdateSoilProps: Relative Soil Saturation Top Moisture <= 0.0001, Value=['// &
trim(RoundSigDigits(RelativeSoilSaturationTop,5))//'].')
CALL ShowContinueError('Value is set to 0.0001 and simulation continues.')
CALL ShowContinueError('You may wish to increase the number of timesteps to attempt to alleviate the problem.')
ENDIF
CALL ShowRecurringWarningErrorAtEnd('EcoRoof: UpdateSoilProps: Relative Soil Saturation Top Moisture < 0. continues', &
ErrIndex,ReportMaxOf=RelativeSoilSaturationTop,ReportMinOf=RelativeSoilSaturationTop)
RelativeSoilSaturationTop=0.0001d0
ENDIF
SoilHydroConductivityTop=SoilConductivitySaturation*(RelativeSoilSaturationTop**lambda)* &
(1.d0-(1.d0-(RelativeSoilSaturationTop)**(n/(n-1.d0)))**((n-1.d0)/n))**2
CapillaryPotentialTop=(-1.d0/alpha)*((((1.d0/RelativeSoilSaturationTop)**(n/(n-1.d0)))-1.d0)**(1.d0/n))
!Then the soil parameters for the root soil layer
RelativeSoilSaturationRoot=(MeanRootMoisture-MoistureResidual)/(MoistureMax-MoistureResidual)
SoilHydroConductivityRoot=SoilConductivitySaturation*(RelativeSoilSaturationRoot**lambda)* &
(1.d0-(1.d0-(RelativeSoilSaturationRoot)**(n/(n-1.d0)))**((n-1.d0)/n))**2
CapillaryPotentialRoot=(-1.d0/alpha)*((((1.d0/RelativeSoilSaturationRoot)**(n/(n-1.d0)))-1.d0)**(1.d0/n))
!Next, using the soil parameters, solve for the soil moisture
SoilConductivityAveTop=(SoilHydroConductivityTop+SoilHydroConductivityRoot)*0.5d0
Moisture=Moisture+(SecondsPerTimeStep/TopDepth)*((SoilConductivityAveTop*(CapillaryPotentialTop&
-CapillaryPotentialRoot)/TopDepth)-SoilConductivityAveTop)
!Now limit the soil from going over the moisture maximum and takes excess to create runoff
IF(Moisture .GE. MoistureMax) THEN !This statement makes sure that the top layer is not over the moisture maximum for the soil.
Moisture=0.9999d0*MoistureMax !then it takes any moisture over the maximum amount and makes it runoff
CurrentRunOff=CurrentRunOff+(Moisture-MoistureMax*0.9999d0)*TopDepth
END IF
!Now make sure that the soil does not go below the moisture minimum
IF(Moisture .LE. (1.01d0*MoistureResidual)) THEN
Moisture=1.01d0*MoistureResidual
END IF
!Next, solve the parameters for the bottom layer
SoilConductivityAveRoot=SoilHydroConductivityRoot
!Now make sure the rate of liquid leaving the soil is more than one drop per hour
IF((SoilConductivityAveRoot*3600.d0) .LE. (2.33d-7)) THEN
SoilConductivityAveRoot=0.0d0
END IF
!Using the parameters above, distribute the Root Layer moisture
TestMoisture=MeanRootMoisture
MeanRootMoisture=MeanRootMoisture+(SecondsPerTimeStep/RootDepth)* &
((SoilConductivityAveTop*(CapillaryPotentialTop-CapillaryPotentialRoot) &
/RootDepth)+SoilConductivityAveTop-SoilConductivityAveRoot)
!Limit the moisture from going over the saturation limit and create runoff:
IF(MeanRootMoisture .GE. MoistureMax) THEN
MeanRootMoisture=0.9999d0*MoistureMax
CurrentRunOff=CurrentRunOff+(Moisture-MoistureMax*0.9999d0)*RootDepth
END IF
!Limit the soil from going below the soil saturation limit:
IF(MeanRootMoisture .LE. (1.01d0*MoistureResidual)) THEN
MeanRootMoisture=1.01d0*MoistureResidual
END IF
!Next, track runoff from the bottom of the soil:
CurrentRunOff=CurrentRunOff+SoilConductivityAveRoot*SecondsPerTimeStep
!~~~END SF EDITS
ENDIF
! NEXT Limit moisture values to saturation (create RUNOFF that we can track)
! CurrentRunoff is sum of "overwatering" in a timestep and excess moisture content
IF (.not. WarmupFlag) THEN
CumRunoff = CumRunoff + CurrentRunoff
ENDIF
if (MeanRootMoisture .LE. MoistureResidual*1.00001d0) then
Moisture = Moisture - (MoistureResidual*1.00001d0 - MeanRootMoisture)*RootDepth/TopDepth
! If the plant has extracted more moisture than is in the root zone soil, then make it come from
! the top layer rather than the root zone... unless top is also dry...
if (Moisture .lt. MoistureResidual*1.00001d0) Moisture = MoistureResidual*1.00001d0
MeanRootMoisture = MoistureResidual*1.00001d0 ! Need to make sure we do not divide by zero later.
endif
! ------------------------------------------------------------------------------------------
! Having updated moisture values now we modify soil thermal properties based on moisture content
! NOTE: Variables SoilAbsorpSolar, SoilDensity, SoilSpecHeat, and SoilConductivity are the values
! that the soil properties OUGHT to attain for the current moisture level. These values are then
! moderated using the TestRatio variable so that from one time step to the next no variable
! changes by more than a certain percentage (typically 5-20%).
! Note wet soil absorptance is generally 25-50% higher than dry soil absorptance (assume linear)
SoilAbsorpSolar = DryAbsorp+ (0.92d0-DryAbsorp)*(Moisture -MoistureResidual)/(MoistureMax - MoistureResidual)
! Limit solar absorptivity to 95% so soil abledo is always above 5%
if (SoilAbsorpSolar .GT. 0.95d0 ) SoilAbsorpSolar = 0.95d0
! Limit solar absorptivity to greater than 20% so that albedo is always less than 80%
if (SoilAbsorpSolar .LT. 0.20d0) SoilAbsorpSolar = 0.20d0
! Need to use for albedo in CalcEcoroof
TestRatio = (1.0d0-SoilAbsorpSolar)/Alphag
If (TestRatio .GT. RatioMax) TestRatio = RatioMax
If (TestRatio .LT. RatioMin) TestRatio = RatioMin
Alphag = TestRatio*Alphag ! included 1.0 - to make it albedo rather than SW absorptivity
!
! Note wet soil density is calculated by simply adding the mass of water...
AvgMoisture = (RootDepth*MeanRootMoisture+TopDepth*Moisture)/SoilThickness
SoilDensity = DryDens + (AvgMoisture-MoistureResidual)*990.0d0
! Note 990 kg/m^3 is water density and the moisture is depth-averaged
! Note wet soil has specific heat that is 40% higher than dry soil (assume linear)
! OLD :: SoilSpecHeat = DrySpecHeat*(1.0d0+ 0.4d0*(AvgMoisture-MoistureResidual)/(MoistureMax-MoistureResidual))
! This is now based on Melos Hagos's results for C (March 2009)
! SoilSpecHeat = DrySpecHeat + 3.09*(AvgMoisture) CLEARLY in ERROR BY FACTOR of 1000
! DJS - Melos report has Spec = Cdry + 1.9 theta (where C is in kJ/kg/K), so...
SoilSpecHeat = DrySpecHeat + 1900.0d0*AvgMoisture
! Note wet soil has thermal conductivity that is up to 3 times that of dry soil ...
! For now simply let it DOUBLE over the range of moisture
! Note wet soil has thermal conductivity that is up to 3 times that of dry soil ...
! OLD :: SoilConductivity = DryCond* (1.0d0 + 1.0d0 * (AvgMoisture-MoistureResidual)/(MoistureMax-MoistureResidual))
! This is now based on Melos Hagos's results for k/kdry (March 2009)
SatRatio = (AvgMoisture - MoistureResidual)/(MoistureMax - MoistureResidual)
SoilConductivity = (DryCond/1.15d0) * (1.45d0 * exp(4.411d0 * SatRatio))/(1.0d0 + 0.45d0 * exp(4.411d0 * SatRatio))
! DJS 2009 - note, this allows the actual conductivity to dip a little below the dry value... corresponding to
! DJS 2009 - "bone dry" if you will, when moisture --> residual value.
! HERE WE RE-RUN THE CONDUCTION TRANSFER FUNCTION (CTF) CALCULATIONS
! NOTE: CTF uses the original Material( )%xxx variables, so must update these for current construction and
! time step...
! TestRatio variable is available just in case there are stability issues. If so, we can limit the amount
! by which soil properties are allowed to vary in one time step (10% in example below).
TestRatio = SoilConductivity/Material(Construct(ConstrNum)%LayerPoint(1))%Conductivity
if (TestRatio .GT. RatioMax) TestRatio=RatioMax
if (TestRatio .LT. RatioMin) TestRatio=RatioMin
Material(Construct(ConstrNum)%LayerPoint(1))%Conductivity = TestRatio*Material(Construct(ConstrNum)%LayerPoint(1))%Conductivity
SoilConductivity =Material(Construct(ConstrNum)%LayerPoint(1))%Conductivity
TestRatio = SoilDensity/Material(Construct(ConstrNum)%LayerPoint(1))%Density
if (TestRatio .GT. RatioMax) TestRatio=RatioMax
if (TestRatio .LT. RatioMin) TestRatio=RatioMin
Material(Construct(ConstrNum)%LayerPoint(1))%Density = Material(Construct(ConstrNum)%LayerPoint(1))%Density*TestRatio
SoilDensity =Material(Construct(ConstrNum)%LayerPoint(1))%Density
TestRatio = SoilSpecHeat/Material(Construct(ConstrNum)%LayerPoint(1))%SpecHeat
if (TestRatio .GT. RatioMax) TestRatio=RatioMax
if (TestRatio .LT. RatioMin) TestRatio=RatioMin
Material(Construct(ConstrNum)%LayerPoint(1))%SpecHeat = Material(Construct(ConstrNum)%LayerPoint(1))%SpecHeat*TestRatio
SoilSpecHeat =Material(Construct(ConstrNum)%LayerPoint(1))%SpecHeat
! Now call InitConductionTransferFunction with the ConstrNum as the argument. As long as the argument is
! non-zero InitConductionTransferFunction will ONLY update this construction. If the argument is 0 it will
! rerun the ENTIRE InitConductionTransferFunction on all constructions (as in initial code start-up mode).
! DJS The following works for most simulations, but has stability issues in some cases.
! DJS - in the present version it seems best to NOT update soil thermal properties.
! Call InitConductionTransferFunctions(ConstrNum)
! DJS In future revision we will implement these modified soil thermal properties in the finite difference
! solution scheme.
! DJS - Note the following write/format statements should be commented out prior to releasing within EnergyPlus
! DJS - they are handy in doing hourly validation/evaluations, though, so leave them in for future development.
! write(unit,799) DayofYear, HourofDay, Qsoil,Tg, Tf, Moisture, MeanRootMoisture,CumPrecip &
! ,CumET,CumRunoff, CumIrrigation, SoilDensity, SoilSpecHeat,SoilConductivity,Alphag
!799 format(' ',I3,' ',I3,' ',' ',f9.3,' ',f6.2,' ',f6.2,' ',f5.3,' ',f5.3,' ',f6.4, ' ' &
! f7.3, ' ', f7.3, ' ',f7.3, ' ',f6.1,' ',f7.1,' ',f6.3,' ',f6.2)
END SUBROUTINE UpdateSoilProps