SetUpDesignDay Subroutine

SUBROUTINE SetUpDesignDay(EnvrnNum)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Linda Lawrie
          !       DATE WRITTEN   February 1977
          !       MODIFIED       June 1997 (RKS); May 2013 (LKL) add temperature profile for drybulb.
          !       RE-ENGINEERED  August 2003;LKL -- to generate timestep weather for design days.

          ! PURPOSE OF THIS SUBROUTINE:
          ! This purpose of this subroutine is to convert the user supplied input
          ! values for the design day parameters into an entire weather day
          ! record.  This now bypasses any file I/O by keeping all of the
          ! weather day record information in the local module level derived type
          ! called DesignDay.

          ! METHODOLOGY EMPLOYED:
          ! Methodology incorporates the design day setup from Tarp as appropriate.

          ! REFERENCES:
          ! ASHRAE Handbook of Fundamentals?

          ! USE STATEMENTS:
  USE General, ONLY: RoundSigDigits, JulianDay
  USE InputProcessor, ONLY: SameString, MakeUPPERcase
  USE ScheduleManager, ONLY:GetSingleDayScheduleValues

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN) :: EnvrnNum  ! Environment number passed into the routine

          ! SUBROUTINE PARAMETER DEFINITIONS:
  REAL(r64), PARAMETER :: GlobalSolarConstant=1367.d0
  REAL(r64), PARAMETER :: ZHGlobalSolarConstant=1355.d0
  CHARACTER(len=*), PARAMETER :: EnvDDHdFormat="('! <Environment:Design Day Data>, Max Dry-Bulb Temp {C}, ',  &
                                                      & 'Temp Range {dC}, Temp Range Ind Type, ',  &
                                                      & 'Hum Ind Value at Max Temp, Hum Ind Type,Pressure {Pa}, ',  &
                                                      & 'Wind Direction {deg CW from N}, ',   &
                                                      & 'Wind Speed {m/s}, Clearness, Rain, Snow')"
  CHARACTER(len=*), PARAMETER :: EnvDDayFormat="('Environment:Design Day Data,')"
  CHARACTER(len=*), PARAMETER :: DDayMiscHdFormat="('! <Environment:Design_Day_Misc>,DayOfYear,ASHRAE A Coeff,',  &
                                                      & 'ASHRAE B Coeff,ASHRAE C Coeff,Solar Constant-Annual Variation,',  &
                                                      & 'Eq of Time {minutes}, Solar Declination Angle {deg}, Solar Model')"
  CHARACTER(len=*), PARAMETER :: DDayMiscFormat="('Environment:Design_Day_Misc,',I3,',')"
  CHARACTER(len=*), PARAMETER :: fmta='(A)'
  CHARACTER(len=*), PARAMETER :: MnDyFmt="(I2.2,'/',I2.2)"
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_C0=.5598d0  !37.6865d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_C1=.4982d0  !13.9263d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_C2=-.6762d0 !-20.2354d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_C3=.02842d0 !0.9695d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_C4=-.00317d0 !-0.2046d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_C5=.014d0    !-0.0980d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_D=-17.853d0  !-10.8568d0
  REAL(r64), PARAMETER :: ZhangHuangModCoeff_K=.843d0  !49.3112d0

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
  TYPE HourlyWeatherData
    REAL(r64), DIMENSION(24) :: BeamSolarRad   = 0.0d0     ! Hourly direct normal solar irradiance
    REAL(r64), DIMENSION(24) :: DifSolarRad    = 0.0d0     ! Hourly sky diffuse horizontal solar irradiance
  END TYPE

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  INTEGER  Hour
  INTEGER  TS
  REAL(r64)     A             ! Apparent solar irradiation at air mass = 0
  REAL(r64)     AVSC          ! Annual variation in the solar constant
  REAL(r64)     B             ! Atmospheric extinction coefficient
  REAL(r64)     C             ! ASHRAE diffuse radiation factor
  REAL(r64)     ETR           ! radiation of an extraterrestrial normal surface, W/m2
  REAL(r64)     HO            ! Radiation on an extraterrestial horizontal surface
  REAL(r64)     KT            ! Radiation ratio
  REAL(r64) SUNCOS(3)         ! Sun direction cosines
  INTEGER CurrentYear
  INTEGER OSky                ! Opaque Sky Cover (tenths)
  REAL(r64) HumidityRatio     ! Humidity Ratio -- when constant for day
  REAL(r64) TDewK             ! Dewpoint in Kelvin
  REAL(r64) ESky              ! Emissivitity of Sky
  REAL(r64) CosZenith         ! Cosine of Zenith Angle of Sun
  REAL(r64) TotHoriz          ! Total Radiation on Horizontal Surface
  REAL(r64) GndReflet         ! Ground Reflectivity
  REAL(r64) CurTime           ! For Solar Calcs
  REAL(r64) WetBulb           ! For calculating
  REAL(r64) DBRange           ! working copy of dry-bulb daily range, C (or 1 if input is difference)
  REAL(r64) WBRange           ! working copy of wet-bulb daily range. C (or 1 if input is difference)

  INTEGER, DIMENSION(8) :: Date0
  LOGICAL, SAVE :: PrintDDHeader
  CHARACTER(len=3) AlpUseRain
  CHARACTER(len=3) AlpUseSnow
  REAL(r64) ::  LastHrBeamSolarRad     ! Direct normal solar irradiance
  REAL(r64) ::  LastHrDifSolarRad      ! Sky diffuse horizontal solar irradiance
  REAL(r64) ::  NextHrBeamSolarRad     ! Direct normal solar irradiance
  REAL(r64) ::  NextHrDifSolarRad      ! Sky diffuse horizontal solar irradiance
  LOGICAL :: ConstantHumidityRatio
  REAL(r64) OutHumRat
  REAL(r64) WgtHourNow
  REAL(r64) WgtPrevHour
  REAL(r64) WgtNextHour
  CHARACTER(len=75) :: StringOut
  TYPE (HourlyWeatherData) :: Wthr
  LOGICAL :: SaveWarmupFlag
  REAL(r64) :: GloHorzRad
  REAL(r64) :: ClearnessIndex_kt
  REAL(r64) :: ClearnessIndex_ktc
  REAL(r64) :: ClearnessIndex_kds
  REAL(r64) :: SinSolarAltitude
  REAL(r64) :: TotSkyCover
  INTEGER :: Hour1Ago,Hour3Ago
  REAL(r64) :: BeamRad, DiffRad     ! working calculated beam and diffuse rad, W/m2
  REAL(r64) :: testval
!     For reporting purposes, set year to current system year
   SaveWarmupFlag=WarmupFlag
   WarmupFlag=.true.

   CALL DATE_AND_TIME(Values=Date0)
   CurrentYear=Date0(1)

   IF (BeginSimFlag) THEN
     PrintDDHeader=.true.
   ENDIF

   DesignDay(EnvrnNum)%Year = CurrentYear ! f90 date_and_time implemented. full 4 digit year !+ 1900
   DesignDay(EnvrnNum)%Month = DesDayInput(EnvrnNum)%Month
   DesignDay(EnvrnNum)%DayOfMonth = DesDayInput(EnvrnNum)%DayOfMonth
   DesignDay(EnvrnNum)%DayOfYear = JulianDay(DesignDay(EnvrnNum)%Month,DesignDay(EnvrnNum)%DayOfMonth,0)
   WRITE(CurMnDy,MnDyFmt) DesDayInput(EnvrnNum)%Month,DesDayInput(EnvrnNum)%DayOfMonth
   EnvironmentName=DesDayInput(EnvrnNum)%Title
   RunPeriodEnvironment=.false.
          ! Following builds Environment start/end for ASHRAE 55 warnings
   EnvironmentStartEnd=TRIM(CurMnDy)//' - '//TRIM(CurMnDy)

   ! Check that barometric pressure is within range
   IF (DesDayInput(EnvrnNum)%PressureEntered) THEN
     IF (ABS((DesDayInput(EnvrnNum)%PressBarom-StdBaroPress)/StdBaroPress) > .1d0) THEN  ! 10% off
       CALL ShowWarningError('SetUpDesignDay: Entered DesignDay Barometric Pressure='//  &
                             TRIM(RoundSigDigits(DesDayInput(EnvrnNum)%PressBarom,0))//   &
                             ' differs by more than 10% from Standard Barometric Pressure='//  &
                             TRIM(RoundSigDigits(StdBaroPress,0))//'.')
       CALL ShowContinueError('...occurs in DesignDay='//TRIM(EnvironmentName)//  &
                              ', Standard Pressure (based on elevation) will be used.')
       DesDayInput(EnvrnNum)%PressBarom=StdBaroPress
     ENDIF
   ELSE
     DesDayInput(EnvrnNum)%PressBarom=StdBaroPress
   ENDIF

   ! verify that design WB or DP <= design DB
   IF (DesDayInput(EnvrnNum)%HumIndType == DDHumIndType_Dewpoint .and. DesDayInput(EnvrnNum)%DewpointNeedsSet) THEN
     ! dew-point
     testval=PsyWFnTdbRhPb(DesDayInput(EnvrnNum)%MaxDryBulb,1.0d0,DesDayInput(EnvrnNum)%PressBarom)
     DesDayInput(EnvrnNum)%HumIndValue=PsyTdpFnWPb(testval,DesDayInput(EnvrnNum)%PressBarom)
   ENDIF

   ! Day of week defaults to Monday, if day type specified, then that is used.
   DesignDay(EnvrnNum)%DayOfWeek = 2
   IF (DesDayInput(EnvrnNum)%DayType <= 7) DesignDay(EnvrnNum)%DayOfWeek = DesDayInput(EnvrnNum)%DayType

   ! set Holiday as indicated by user input
   DesignDay(EnvrnNum)%HolidayIndex = 0
   IF (DesDayInput(EnvrnNum)%DayType > 7) DesignDay(EnvrnNum)%HolidayIndex = DesDayInput(EnvrnNum)%DayType-7

   DesignDay(EnvrnNum)%DaylightSavingIndex = DesDayInput(EnvrnNum)%DSTIndicator

   !  Set up Solar parameters for day
   CALL CalculateDailySolarCoeffs(DesignDay(EnvrnNum)%DayOfYear,A,B,C,AVSC,DesignDay(EnvrnNum)%EquationOfTime, &
                                 DesignDay(EnvrnNum)%SinSolarDeclinAngle,DesignDay(EnvrnNum)%CosSolarDeclinAngle)

   IF (PrintDDHeader .and. DoWeatherInitReporting) THEN
     WRITE(OutputFileInits,EnvDDHdFormat)
     WRITE(OutputFileInits,DDayMiscHdFormat)
     PrintDDHeader=.false.
   ENDIF
   IF (DoWeatherInitReporting) THEN
     IF (DesDayInput(Envrn)%RainInd == 1) THEN
       AlpUseRain='Yes'
     ELSE
       AlpUseRain='No'
     ENDIF
     IF (DesDayInput(Envrn)%SnowInd == 1) THEN
       AlpUseSnow='Yes'
     ELSE
       AlpUseSnow='No'
     ENDIF
     WRITE(OutputFileInits,EnvDDayFormat,advance='No')
     StringOut=RoundSigDigits(DesDayInput(Envrn)%MaxDryBulb,2)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(DesDayInput(Envrn)%DailyDBRange,2)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=','
     IF (DesDayInput(Envrn)%DBTempRangeType == DDDBRangeType_Default) THEN
       StringOut='DefaultMultipliers,'
     ELSEIF (DesDayInput(Envrn)%DBTempRangeType == DDDBRangeType_Multiplier) THEN
       StringOut='MultiplierSchedule,'
     ELSEIF (DesDayInput(Envrn)%DBTempRangeType == DDDBRangeType_Profile) THEN
       StringOut='TemperatureProfile,'
     ELSEIF (DesDayInput(Envrn)%DBTempRangeType == DDDBRangeType_Difference) THEN
       StringOut='DifferenceSchedule,'
     ENDIF
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)
     IF (DesDayInput(Envrn)%HumIndType == DDHumIndType_Wetbulb) THEN
       StringOut='Wetbulb,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,2))//' {C},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_Dewpoint) THEN
       StringOut='Dewpoint,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,2))//' {C},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_Enthalpy) THEN
       StringOut='Enthalpy,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,2))//' {kJ/kg},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_HumRatio) THEN
       StringOut='HumidityRatio,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,4))//' {},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_RelHumSch) THEN
       StringOut='Schedule,<schedule values from 0.0 to 100.0 {percent},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_WBProfDef) THEN
       StringOut='WetBulbProfileDefaultMultipliers,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,2))//' {C},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_WBProfDif) THEN
       StringOut='WetBulbProfileDifferenceSchedule,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,2))//' {C},'
     ELSEIF (DesDayInput(Envrn)%HumIndType == DDHumIndType_WBProfMul) THEN
       StringOut='WetBulbProfileMultiplierSchedule,'//trim(RoundSigDigits(DesDayInput(Envrn)%HumIndValue,2))//' {C},'
     ENDIF
     StringOut=RoundSigDigits(DesDayInput(Envrn)%PressBarom,0)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(DesDayInput(Envrn)%WindDir,0)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(DesDayInput(Envrn)%WindSpeed,1)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(DesDayInput(Envrn)%SkyClear,2)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     WRITE(OutputFileInits,fmta) AlpUseRain//','//AlpUseSnow

     WRITE(OutputFileInits,DDayMiscFormat,advance='No') DesignDay(EnvrnNum)%DayOfYear
     StringOut=RoundSigDigits(A,1)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(B,4)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(C,4)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(AVSC,1)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(DesignDay(EnvrnNum)%EquationOfTime*60.0d0,2)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     StringOut=RoundSigDigits(ASIN(DesignDay(EnvrnNum)%SinSolarDeclinAngle)/DegToRadians,1)
     WRITE(OutputFileInits,fmta,advance='No') TRIM(StringOut)//','
     IF (DesDayInput(EnvrnNum)%SolarModel == ASHRAE_ClearSky) THEN
       StringOut='ASHRAEClearSky'
     ELSEIF (DesDayInput(EnvrnNum)%SolarModel == Zhang_Huang) THEN
       StringOut='ZhangHuang'
     ELSEIF (DesDayInput(EnvrnNum)%SolarModel == SolarModel_Schedule) THEN
       StringOut='User supplied beam/diffuse from schedules'
     ELSEIF (DesDayInput(EnvrnNum)%SolarModel == ASHRAE_Tau) THEN
       StringOut='ASHRAETau'
     ELSE
       StringOut='unknown'
     ENDIF
     WRITE(OutputFileInits,fmta) TRIM(StringOut)
   ENDIF


   ! Must set up weather values for Design Day.  User can specify the "humidity indicator" as
   ! Wetbulb, DewPoint or input the relative humidity schedule.  For both wetbulb and dewpoint indicators, the
   ! humidity for the day will be constant, using the drybulb (max) and humidity indicator temperature to
   ! set the values.  For the scheduled values, these are already set in the DDxxx array.

   CurrentTime=25.0d0

   SELECT CASE(DesDayInput(Envrn)%HumIndType)

   CASE (DDHumIndType_Wetbulb)
     HumidityRatio=PsyWFnTdbTwbPb(DesDayInput(EnvrnNum)%MaxDryBulb,          &
                          DesDayInput(EnvrnNum)%HumIndValue,          &
                          DesDayInput(EnvrnNum)%PressBarom,'SetUpDesignDay:PsyWFnTdbTwbPb')
     ConstantHumidityRatio=.true.

   CASE (DDHumIndType_Dewpoint)
     HumidityRatio=PsyWFnTdpPb(DesDayInput(EnvrnNum)%HumIndValue,          &
                          DesDayInput(EnvrnNum)%PressBarom,'SetUpDesignDay:PsyWFnTdpPb')
     ConstantHumidityRatio=.true.

   CASE (DDHumIndType_HumRatio)
     HumidityRatio=DesDayInput(EnvrnNum)%HumIndValue
     ConstantHumidityRatio=.true.

   CASE (DDHumIndType_Enthalpy)
     HumidityRatio=PsyWFnTdbH(DesDayInput(EnvrnNum)%MaxDryBulb,DesDayInput(EnvrnNum)%HumIndValue*1000.0d0,  &
        'SetUpDesignDay:PsyWFnTdbH')
     ConstantHumidityRatio=.true.

   CASE (DDHumIndType_RelHumSch)
     ! nothing to do -- DDHumIndModifier already contains the scheduled Relative Humidity
     ConstantHumidityRatio=.false.
     TomorrowOutRelHum=DDHumIndModifier( EnvrnNum,:,:)

   CASE (DDHumIndType_WBProfDef, DDHumIndType_WBProfDif, DDHumIndType_WBProfMul)
     ConstantHumidityRatio = .false.

   CASE DEFAULT
     CALL ShowSevereError('SetUpDesignDay: Invalid Humidity Indicator type')
     CALL ShowContinueError('Occurred in Design Day='//TRIM(DesDayInput(Envrn)%Title))

   END SELECT

   IF (DesDayInput(EnvrnNum)%RainInd /= 0) THEN
     TomorrowIsRain(:,:)=.true.
     OSky=10
     TomorrowLiquidPrecip=3.0d0
   ELSE
     TomorrowIsRain(:,:)=.false.
     OSky=0
     TomorrowLiquidPrecip=0.0d0
   ENDIF

   IF (DesDayInput(EnvrnNum)%SnowInd == 0) THEN
     TomorrowIsSnow(:,:)=.false.
     GndReflet=.2d0
   ELSE  ! Snow
     TomorrowIsSnow(:,:)=.true.
     GndReflet = .7d0
   ENDIF

   ! Some values are constant

   TomorrowOutBaroPress(:,:) = DesDayInput(EnvrnNum)%PressBarom
   TomorrowWindSpeed(:,:)    = DesDayInput(EnvrnNum)%WindSpeed
   TomorrowWindDir(:,:)      = DesDayInput(EnvrnNum)%WindDir
   TomorrowAlbedo=0.0d0

   ! resolve daily ranges
   IF (DesDayInput(EnvrnNum)%DBTempRangeType == DDDBRangeType_Difference) THEN
      DBRange = 1.d0      ! use unscaled multiplier values if difference
   ELSEIF (DesDayInput(EnvrnNum)%DBTempRangeType == DDDBRangeType_Profile) THEN
      DBRange = 0.0d0
   ELSE
      DBRange = DesDayInput(EnvrnNum)%DailyDBRange
   ENDIF
   IF (DesDayInput(EnvrnNum)%HumIndType == DDHumIndType_WBProfDif) THEN
      WBRange = 1.d0      ! use unscaled multiplier values if difference
   ELSE
      WBRange = DesDayInput(EnvrnNum)%DailyWBRange
   ENDIF

   DO Hour = 1,24
     DO TS=1,NumOfTimeStepInHour

       IF (DesDayInput(EnvrnNum)%DBTempRangeType /= DDDBRangeType_Profile) THEN
         ! dry-bulb profile
         TomorrowOutDryBulbTemp(Hour,TS)=DesDayInput(EnvrnNum)%MaxDryBulb - DDDBRngModifier(EnvrnNum,Hour,TS)*DBRange
       ELSE ! DesDayInput(EnvrnNum)%DBTempRangeType == DDDBRangeType_Profile
         TomorrowOutDryBulbTemp(Hour,TS)=DDDBRngModifier(EnvrnNum,Hour,TS)
       ENDIF


       ! wet-bulb - generate from profile, humidity ratio, or dew point
       IF (DesDayInput(EnvrnNum)%HumIndType == DDHumIndType_WBProfDef         &
        .or. DesDayInput(EnvrnNum)%HumIndType == DDHumIndType_WBProfDif       &
        .or. DesDayInput(EnvrnNum)%HumIndType == DDHumIndType_WBProfMul) THEN
         WetBulb = DesDayInput(EnvrnNum)%HumIndValue - DDHumIndModifier(EnvrnNum,Hour,TS)*WBRange
         WetBulb = MIN( WetBulb, TomorrowOutDryBulbTemp(Hour,TS))   ! WB must be <= DB
         OutHumRat = PsyWFnTdbTwbPb(TomorrowOutDryBulbTemp(Hour,TS), &
                              WetBulb, DesDayInput(EnvrnNum)%PressBarom)
         TomorrowOutDewPointTemp(Hour,TS)=PsyTdpFnWPb(OutHumRat,DesDayInput(EnvrnNum)%PressBarom)
         TomorrowOutRelHum(Hour,TS)=PsyRhFnTdbWPb(TomorrowOutDryBulbTemp(Hour,TS),OutHumRat, &
                                               DesDayInput(EnvrnNum)%PressBarom,'WeatherManager')*100.0d0
       ELSE IF (ConstantHumidityRatio) THEN
       !  Need Dew Point Temperature.  Use Relative Humidity to get Humidity Ratio, unless Humidity Ratio is constant
        !BG 9-26-07  moved following inside this IF statment; when HumIndType is 'Schedule' HumidityRatio wasn't being initialized
         WetBulb = PsyTwbFnTdbWPb(TomorrowOutDryBulbTemp(Hour,TS),HumidityRatio,DesDayInput(EnvrnNum)%PressBarom, &
         'WeatherManager.f90 subroutine SetUpDesignDay')

         OutHumRat = PsyWFnTdpPb(TomorrowOutDryBulbTemp(Hour,TS),DesDayInput(EnvrnNum)%PressBarom)
         IF (HumidityRatio > OutHumRat) THEN
           WetBulb=TomorrowOutDryBulbTemp(Hour,TS)
         ELSE
           OutHumRat = PsyWFnTdbTwbPb(TomorrowOutDryBulbTemp(Hour,TS), &
                              WetBulb,DesDayInput(EnvrnNum)%PressBarom)
         ENDIF
         TomorrowOutDewPointTemp(Hour,TS)=PsyTdpFnWPb(OutHumRat,DesDayInput(EnvrnNum)%PressBarom)
         TomorrowOutRelHum(Hour,TS)=PsyRhFnTdbWPb(TomorrowOutDryBulbTemp(Hour,TS),OutHumRat, &
                                               DesDayInput(EnvrnNum)%PressBarom,'WeatherManager')*100.0d0
       ELSE
         HumidityRatio=PsyWFnTdbRhPb(TomorrowOutDryBulbTemp(Hour,TS),DDHumIndModifier(EnvrnNum,Hour,TS)/100.0d0,  &
                                               DesDayInput(EnvrnNum)%PressBarom)
         ! TomorrowOutRelHum values set earlier
         TomorrowOutDewPointTemp(Hour,TS)=PsyTdpFnWPb(HumidityRatio,DesDayInput(EnvrnNum)%PressBarom)
       ENDIF

       ! Determine Sky Temp ==>
       ! Function of DryBulb, DewPoint, OpaqueSkyCover
       ! Calculate Sky IR
       !HIR = ESKY * SIGMA * (TOUT**4)
       !
       !where
       !
       !HIR = horizontal IR intensity (W/m2)
       !ESKY = sky emissivity
       !SIGMA = Stefan-Boltzmann constant = 5.6697e-8 W/m2-K4
       !TOUT = drybulb temperature (K)
       !
       !The sky emissivity is given by
       !
       !ESKY = [0.787 + 0.764*ln(TDEW/273)]*[1 + 0.0224*N - 0.0035*(N**2) + 0.00028*(N**3)]
       !
       !where
       !
       !TDEW = dewpoint temperature (K)
       !N = opaque sky cover (tenths)
       !
       !Example: Clear sky (N=0), TOUT = 273+20=293K, TDEW = 273+10=283K:
       !
       !ESKY = 0.787 + 0.764*0.036 = 0.815
       !HIR = 0.815*5.6697e-8*(293**4) = 340.6 W/m2

       !References:
       !
       !George N. Walton, "Thermal Analysis Research Program Reference Manual,"
       !NBSIR 83-2655, March 1983, p. 21.
       !
       !G. Clark and C. Allen, "The Estimation of Atmospheric Radiation for Clear and
       !Cloudy Skies," Proc. 2nd National Passive Solar Conference (AS/ISES), 1978, pp. 675-678.

       IF (Environment(EnvrnNum)%WP_Type1 == 0) THEN
         TDewK=MIN(TomorrowOutDryBulbTemp(Hour,TS),TomorrowOutDewPointTemp(Hour,TS))+TKelvin
         ESky= (.787d0 +.764d0*LOG((TDewK)/TKelvin))*(1.d0 + .0224d0*OSky - 0.0035d0*(OSky**2) + .00028d0*(OSky**3))
         TomorrowHorizIRSky(Hour,TS)=Esky*Sigma*(TomorrowOutDryBulbTemp(Hour,TS)+TKelvin)**4
         TomorrowSkyTemp(Hour,TS)=(TomorrowOutDryBulbTemp(Hour,TS)+TKelvin)*(ESky**.25d0)-TKelvin
       ELSE
         TDewK=MIN(TomorrowOutDryBulbTemp(Hour,TS),TomorrowOutDewPointTemp(Hour,TS))+TKelvin
         ESky= (.787d0 +.764d0*LOG((TDewK)/TKelvin))*(1.d0 + .0224d0*OSky - 0.0035d0*(OSky**2) + .00028d0*(OSky**3))
         TomorrowHorizIRSky(Hour,TS)=Esky*Sigma*(TomorrowOutDryBulbTemp(Hour,TS)+TKelvin)**4
       ENDIF

       ! Generate solar values for timestep
       !    working results = BeamRad and DiffRad
       !    stored to program globals at end of loop
       IF (DesDayInput(EnvrnNum)%SolarModel == SolarModel_Schedule) THEN
         ! scheduled: set value unconditionally (whether sun up or not)
         BeamRad=DDBeamSolarValues(EnvrnNum,Hour,TS)
         DiffRad=DDDiffuseSolarValues(EnvrnNum,Hour,TS)
       ELSE

       ! calc time = fractional hour of day
         IF (NumOfTimeStepInHour /= 1) THEN
           CurTime = REAL(Hour-1,r64) + REAL(TS,r64)*TimeStepFraction
         ELSE
           CurTime = REAL(Hour,r64)+TS1TimeOffset
         ENDIF

         CALL CalculateSunDirectionCosines(CurTime,DesignDay(EnvrnNum)%EquationOfTime,DesignDay(EnvrnNum)%SinSolarDeclinAngle,  &
                                         DesignDay(EnvrnNum)%CosSolarDeclinAngle,SUNCOS)
         CosZenith=SUNCOS(3)
         IF (CosZenith < SunIsUpValue) THEN
           BeamRad = 0.d0
           DiffRad = 0.d0
         ELSE
           SinSolarAltitude=SUNCOS(3)

           SELECT CASE (DesDayInput(EnvrnNum)%SolarModel)

           CASE (Ashrae_ClearSky)
             TotHoriz = DesDayInput(EnvrnNum)%SkyClear * A * (C + CosZenith) * EXP( -B / CosZenith)
             HO=GlobalSolarConstant*AVSC*CosZenith
             KT=TotHoriz/HO
             KT=MIN(KT,.75d0)
             DiffRad = TotHoriz * (1.0045d0 + KT * (.04349d0 + KT * (-3.5227d0 + 2.6313d0 * KT)))
             IF (DesDayInput(EnvrnNum)%SkyClear .GT. 0.70d0) DiffRad = TotHoriz*C/(C+CosZenith)
             BeamRad = (TotHoriz-DiffRad)/CosZenith
             DiffRad = MAX(0.0d0,DiffRad)
             BeamRad = MAX(0.0d0,BeamRad)

           CASE (ASHRAE_Tau)
             ETR = GlobalSolarConstant*AVSC      ! extraterrestrial normal irrad, W/m2
             CALL ASHRAETauModel( ETR, CosZenith, DesDayInput(EnvrnNum)%TauB, DesDayInput(EnvrnNum)%TauD, &
                  BeamRad, DiffRad, GloHorzRad)

           CASE (Zhang_Huang)
             Hour3Ago = MOD( Hour+20, 24)+1       ! hour 3 hours before
             TotSkyCover=MAX( 1.0d0-DesDayInput(EnvrnNum)%SkyClear, 0.0d0)
             GloHorzRad = ( ZHGlobalSolarConstant * SinSolarAltitude * (ZhangHuangModCoeff_C0 &
                 + ZhangHuangModCoeff_C1 * TotSkyCover &
                 + ZhangHuangModCoeff_C2 * (TotSkyCover)**2 &
                 + ZhangHuangModCoeff_C3 * (TomorrowOutDryBulbTemp(Hour,TS) - TomorrowOutDryBulbTemp(Hour3Ago,TS)) &
                 + ZhangHuangModCoeff_C4 * TomorrowOutRelHum(Hour,TS) &
                 + ZhangHuangModCoeff_C5 * TomorrowWindSpeed(Hour,TS)) + ZhangHuangModCoeff_D ) &
                 / ZhangHuangModCoeff_K
             GloHorzRad = MAX( GloHorzRad,0.0d0)
             ClearnessIndex_kt=GloHorzRad/(GlobalSolarConstant * SinSolarAltitude)
  !          ClearnessIndex_kt=DesDayInput(EnvrnNum)%SkyClear
             ClearnessIndex_Ktc = 0.4268d0 + 0.1934d0 * SinSolarAltitude
             IF (ClearnessIndex_Kt < ClearnessIndex_Ktc) THEN
               ClearnessIndex_Kds = (3.996d0 -3.862d0*SinSolarAltitude +1.54d0*SinSolarAltitude**2)* &
                     ClearnessIndex_Kt**3
             ELSE
               ClearnessIndex_Kds = ClearnessIndex_Kt - (1.107d0 + 0.03569d0 * SinSolarAltitude + 1.681d0 * &
                     SinSolarAltitude**2) * (1.d0-ClearnessIndex_Kt)**3
             ENDIF
             ! Calculate direct normal radiation, W/m2
             BeamRad = ZHGlobalSolarConstant * SinSolarAltitude * ClearnessIndex_Kds * &
                  ((1.d0 - ClearnessIndex_Kt) / (1.d0 - ClearnessIndex_Kds))
             ! Calculation diffuse horizontal radiation, W/m2
             DiffRad = ZHGlobalSolarConstant * SinSolarAltitude * &
                  ((ClearnessIndex_Kt - ClearnessIndex_Kds) / (1.d0 - ClearnessIndex_Kds))

           CASE DEFAULT
         END SELECT
       END IF
     END IF

     ! override result to 0 per environment var (for testing)
     IF (IgnoreSolarRadiation .or. IgnoreBeamRadiation) BeamRad = 0.0d0
     IF (IgnoreSolarRadiation .or. IgnoreDiffuseRadiation) DiffRad = 0.0d0;

     TomorrowBeamSolarRad( Hour,TS) = BeamRad
     TomorrowDifSolarRad( Hour,TS)  = DiffRad

     END DO   ! Timestep (TS) Loop
   ENDDO  ! Hour Loop

   ! back-fill hour values from timesteps
   ! hour values = integrated over hour ending at time of hour
   ! insurance: hourly values not known to be needed
   DO Hour = 1,24
     Hour1Ago = MOD( Hour+22, 24)+1
     BeamRad = (TomorrowBeamSolarRad( Hour1Ago, NumOfTimeStepInHour) + TomorrowBeamSolarRad( Hour, NumOfTimeStepInHour)) / 2.0d0
     DiffRad = (TomorrowDifSolarRad(  Hour1Ago, NumOfTimeStepInHour) + TomorrowDifSolarRad(  Hour, NumOfTimeStepInHour)) / 2.0d0
     IF (NumOfTimeStepInHour > 1) THEN
        BeamRad = BeamRad + SUM( TomorrowBeamSolarRad( Hour, 1:NumOfTimeStepInHour-1))
        DiffRad = DiffRad + SUM( TomorrowDifSolarRad(  Hour, 1:NumOfTimeStepInHour-1))
     ENDIF
     Wthr%BeamSolarRad( Hour) = BeamRad / NumOfTimeStepInHour
     Wthr%DifSolarRad( Hour)  = DiffRad / NumOfTimeStepInHour
   END DO

   IF (Environment(EnvrnNum)%WP_Type1 /= 0) THEN

     SELECT CASE(WPSkyTemperature(Environment(EnvrnNum)%WP_Type1)%CalculationType)

       CASE (WP_ScheduleValue)
         CALL GetSingleDayScheduleValues(WPSkyTemperature(Environment(EnvrnNum)%WP_Type1)%SchedulePtr,  &
                             TomorrowSkyTemp)
         DDSkyTempScheduleValues(EnvrnNum,:,:)=TomorrowSkyTemp
       CASE (WP_DryBulbDelta)
         CALL GetSingleDayScheduleValues(WPSkyTemperature(Environment(EnvrnNum)%WP_Type1)%SchedulePtr,  &
                             TomorrowSkyTemp)
         DDSkyTempScheduleValues(EnvrnNum,:,:)=TomorrowSkyTemp
         DO Hour=1,24
           DO TS=1,NumOfTimeStepInHour
             TomorrowSkyTemp(Hour,TS)=TomorrowOutDryBulbTemp(Hour,TS)-TomorrowSkyTemp(Hour,TS)
           ENDDO
         ENDDO

       CASE (WP_DewPointDelta)
         CALL GetSingleDayScheduleValues(WPSkyTemperature(Environment(EnvrnNum)%WP_Type1)%SchedulePtr,  &
                             TomorrowSkyTemp)
         DDSkyTempScheduleValues(EnvrnNum,:,:)=TomorrowSkyTemp
         DO Hour=1,24
           DO TS=1,NumOfTimeStepInHour
             TomorrowSkyTemp(Hour,TS)=TomorrowOutDewPointTemp(Hour,TS)-TomorrowSkyTemp(Hour,TS)
           ENDDO
         ENDDO

       CASE DEFAULT

     END SELECT

   ENDIF


   WarmupFlag=SaveWarmupFlag

   RETURN

END SUBROUTINE SetUpDesignDay