Building Thermal Inertia & τ

program building_thermal_inertia
    implicit none
    integer :: struct_type, iostat_val
    double precision :: Afloor_m2, Vbldg_m3, U_env_Wm2K, Aenv_m2, ACH_rate
    double precision :: Cm_kJm2K, H_tr_WK, H_ve_WK, H_tot_WK, R_th_KW, C_th_kJ_K
    double precision :: tau_hours, omega_rad_h, damping_factor, phase_lag_hours
    double precision :: dTout_degC, dTin_degC, Q_int_gain_W
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) struct_type     ! 1=Lightweight, 2=Medium, 3=Heavyweight
    read(*,*,iostat=iostat_val) Afloor_m2       ! Floor Area [m2] (e.g. 200.0)
    read(*,*,iostat=iostat_val) Vbldg_m3        ! Conditioned Volume [m3] (e.g. 550.0)
    read(*,*,iostat=iostat_val) Aenv_m2         ! Total Envelope Area [m2] (e.g. 450.0)
    read(*,*,iostat=iostat_val) U_env_Wm2K      ! Average Envelope U-value [W/(m2.K)] (e.g. 0.35)
    read(*,*,iostat=iostat_val) ACH_rate        ! Air Change Rate [h-1] (e.g. 0.5)
    read(*,*,iostat=iostat_val) dTout_degC      ! Outdoor Diurnal Temp Amplitude [deg C] (e.g. 14.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for building thermal inertia.'
        stop
    end if

    ! Construction thermal capacitance per m2 floor
    if (struct_type == 1) then
        Cm_kJm2K = 80.0d0   ! Lightweight timber/drywall
    else if (struct_type == 2) then
        Cm_kJm2K = 165.0d0  ! Medium brick/screed
    else
        Cm_kJm2K = 300.0d0  ! Heavyweight reinforced concrete
    end if

    ! Heat loss coefficients
    H_tr_WK = U_env_Wm2K * Aenv_m2
    H_ve_WK = 0.33d0 * ACH_rate * Vbldg_m3
    H_tot_WK = H_tr_WK + H_ve_WK
    R_th_KW = 1.0d0 / max(1.0d-1, H_tot_WK)
    C_th_kJ_K = Cm_kJm2K * Afloor_m2

    ! Thermal time constant [hours]
    tau_hours = (C_th_kJ_K * 1000.0d0) / (3600.0d0 * H_tot_WK)

    ! Diurnal 24h cycle frequency
    omega_rad_h = 2.0d0 * PI / 24.0d0

    ! Damping factor and thermal phase lag
    damping_factor = 1.0d0 / sqrt(1.0d0 + (omega_rad_h * tau_hours)**2)
    phase_lag_hours = atan(omega_rad_h * tau_hours) / omega_rad_h
    dTin_degC = damping_factor * dTout_degC

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — BUILDING THERMAL INERTIA & TIME CONSTANT'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A)')  'Building Floor Area       = ', Afloor_m2, ' m2'
    write(*,'(A,F10.2,A)')  'Envelope Heat Loss H_tr   = ', H_tr_WK, ' W/K'
    write(*,'(A,F10.2,A)')  'Ventilation Loss H_ve     = ', H_ve_WK, ' W/K'
    write(*,'(A,F10.2,A)')  'Total Loss Coeff H_tot    = ', H_tot_WK, ' W/K'
    write(*,'(A,F10.1,A)')  'Effective Heat Capacity   = ', C_th_kJ_K, ' kJ/K'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'THERMAL TIME CONSTANT tau = ', tau_hours, ' hours'
    write(*,'(A,F10.3)')    'Diurnal Damping Factor mu = ', damping_factor
    write(*,'(A,F10.2,A)')  'Thermal Phase Lag Dt_lag  = ', phase_lag_hours, ' hours'
    write(*,'(A,F10.1,A,F6.1,A)') 'Indoor Temp Swing (Out/In)= ', dTout_degC, ' C -> ', dTin_degC, ' C'
    write(*,'(A)') '============================================================'

end program building_thermal_inertia