! ==============================================================================
! 2D Thermal Bridge Calculator (2D Pont Thermique)
! Steady-State Finite Difference Method (FDM) Solver
! References: ISO 10211, EN ISO 14683
! ==============================================================================
program thermal_bridge
    implicit none
    
    ! Grid size
    integer, parameter :: Nx = 50
    integer, parameter :: Ny = 50
    
    ! Inputs
    integer :: bridge_type ! 1=Balcony slab, 2=Wall corner, 3=Window frame, 4=Floor junction, 5=Column penetration
    real(8) :: Ti, Te ! Inside and outside temperatures (C)
    real(8) :: Rsi, Rse ! Surface resistances (m2-K/W)
    real(8) :: bridge_material_k ! Conductivity of penetrating bridge material (W/m-K)
    real(8) :: bridge_thick ! Bridge thickness (mm)
    real(8) :: bridge_length ! Bridge length/extension (mm)
    integer :: num_layers ! Number of wall layers
    real(8) :: layer_thickness(10) ! Layer thickness (mm)
    real(8) :: layer_k(10) ! Layer conductivity (W/m-K)
    real(8) :: bridge_axis_length ! Length of bridge along wall axis (m)
    
    ! Grid and Solver parameters
    real(8) :: H_domain ! Height of domain (m), fixed at 0.5 m
    real(8) :: W_domain ! Width of domain (m), equal to total wall thickness
    real(8) :: dx, dy
    real(8) :: T(0:Nx+1, 0:Ny+1)
    real(8) :: T_new(Nx, Ny)
    real(8) :: k_grid(0:Nx+1, 0:Ny+1)
    integer :: fluid_type(0:Nx+1, 0:Ny+1) ! 0=Solid, 1=Interior fluid, 2=Exterior fluid
    
    ! Derived parameters
    real(8) :: layer_start(11)
    real(8) :: wall_R_1D, wall_U_1D
    real(8) :: q_2D, q_1D, psi, chi, Q_bridge, f_Rsi
    real(8) :: T_si_min
    
    ! Loop & Solver variables
    integer :: i, j, k, iter
    real(8) :: x, y, error, max_diff
    real(8) :: gW, gE, gN, gS, denom
    
    H_domain = 0.5d0
    
    ! Read standard inputs
    read(*, *) bridge_type
    read(*, *) Ti
    read(*, *) Te
    read(*, *) Rsi
    read(*, *) Rse
    read(*, *) bridge_material_k
    read(*, *) bridge_thick
    read(*, *) bridge_length
    read(*, *) num_layers
    
    do i = 1, num_layers
        read(*, *) layer_thickness(i)
        read(*, *) layer_k(i)
    enddo
    read(*, *) bridge_axis_length
    
    ! Compute wall thickness and grid spacing
    W_domain = 0.0d0
    layer_start(1) = 0.0d0
    do i = 1, num_layers
        W_domain = W_domain + layer_thickness(i) / 1000.0d0 ! convert to meters
        layer_start(i+1) = W_domain
    enddo
    
    dx = W_domain / dble(Nx)
    dy = H_domain / dble(Ny)
    
    ! Setup fluid_type and default k_grid
    fluid_type = 0
    k_grid = 1.0d0
    
    ! Boundary default virtual neighbors
    fluid_type(0, :) = 1 ! Left boundary virtual cells behave like Interior fluid (default)
    fluid_type(Nx+1, :) = 2 ! Right boundary virtual cells behave like Exterior fluid (default)
    
    ! Map wall layers to k_grid
    do i = 0, Nx+1
        x = (dble(i) - 0.5d0) * dx
        if (i == 0) x = 0.0d0
        if (i == Nx+1) x = W_domain
        
        ! Find layer
        k = 1
        do j = 1, num_layers
            if (x >= layer_start(j) .and. x <= layer_start(j+1)) then
                k = j
                exit
            endif
        enddo
        k_grid(i, :) = layer_k(k)
    enddo
    
    ! Apply Specific Thermal Bridge Geometries
    if (bridge_type == 1) then
        ! 1. Balcony Slab: slab centered vertically, penetrates all layers
        do j = 1, Ny
            y = (dble(j) - 0.5d0) * dy
            if (abs(y - H_domain/2.0d0) <= (bridge_thick/2000.0d0)) then
                k_grid(1:Nx, j) = bridge_material_k
            endif
        enddo
        
    else if (bridge_type == 2) then
        ! 2. Wall Corner: L-shaped wall
        ! Solid wall is for i <= 20 or j <= 20
        ! Interior room is for i > 20 and j > 20
        do i = 1, Nx
            do j = 1, Ny
                if (i > 20 .and. j > 20) then
                    fluid_type(i, j) = 1 ! Interior fluid
                endif
            enddo
        enddo
        
        ! For Corner, the outside boundaries at i=1 and j=1 are convective with Te
        ! The virtual cells at i=0 and j=0 behave like Exterior fluid
        fluid_type(0, :) = 2
        fluid_type(:, 0) = 2
        
        ! cutoff boundary at i=Nx (for j <= 20) and j=Ny (for i <= 20) are adiabatic
        fluid_type(Nx+1, :) = 0 ! no fluid outside cutoff
        fluid_type(:, Ny+1) = 0
        
    else if (bridge_type == 3) then
        ! 3. Window Frame: wall meets window frame at j >= 25
        ! Frame is made of wood/PVC (low k) or metal (high k) with gas glazing
        do j = 25, Ny
            do i = 1, Nx
                x = (dble(i) - 0.5d0) * dx
                if (x < W_domain / 3.0d0) then
                    k_grid(i, j) = bridge_material_k ! frame outer/inner profiles
                else if (x < 2.0d0 * W_domain / 3.0d0) then
                    k_grid(i, j) = 0.024d0 ! air/glazing cavity
                else
                    k_grid(i, j) = bridge_material_k
                endif
            enddo
        enddo
        
    else if (bridge_type == 4) then
        ! 4. Floor Junction: concrete floor slab joins from interior (left)
        ! but does not fully penetrate the insulation layer.
        do j = 1, Ny
            y = (dble(j) - 0.5d0) * dy
            if (abs(y - H_domain/2.0d0) <= (bridge_thick/2000.0d0)) then
                do i = 1, Nx
                    x = (dble(i) - 0.5d0) * dx
                    if (num_layers >= 2) then
                        ! Penetrates layer 1 (interior structural wall) but stops at layer 2 (insulation)
                        if (x < layer_start(2)) then
                            k_grid(i, j) = bridge_material_k
                        endif
                    else
                        ! Penetrates halfway if single layer
                        if (x < W_domain / 2.0d0) then
                            k_grid(i, j) = bridge_material_k
                        endif
                    endif
                enddo
            endif
        enddo
        
    else if (bridge_type == 5) then
        ! 5. Column Penetration (Point thermal bridge)
        ! Represented in 2D as column cutting through center
        do j = 1, Ny
            y = (dble(j) - 0.5d0) * dy
            if (abs(y - H_domain/2.0d0) <= (bridge_thick/2000.0d0)) then
                k_grid(1:Nx, j) = bridge_material_k
            endif
        enddo
    endif
    
    ! Initialize temperatures based on fluid_type
    T = (Ti + Te) / 2.0d0
    T_new = (Ti + Te) / 2.0d0
    do i = 0, Nx+1
        do j = 0, Ny+1
            if (fluid_type(i, j) == 1) then
                T(i, j) = Ti
            else if (fluid_type(i, j) == 2) then
                T(i, j) = Te
            endif
        enddo
    enddo
    
    ! ==========================================================================
    ! Gauss-Seidel Solver Loop
    ! ==========================================================================
    do iter = 1, 15000
        max_diff = 0.0d0
        
        do i = 1, Nx
            do j = 1, Ny
                ! Skip fluid cells (their temperatures are fixed)
                if (fluid_type(i, j) /= 0) cycle
                
                ! East neighbor link (i+1, j)
                if (i == Nx .and. bridge_type == 2) then
                    ! corner cutoff boundary is adiabatic
                    gE = 0.0d0
                else
                    if (fluid_type(i+1, j) == 0) then
                        gE = (2.0d0 * k_grid(i, j) * k_grid(i+1, j)) / (k_grid(i, j) + k_grid(i+1, j)) / (dx**2)
                    else if (fluid_type(i+1, j) == 1) then
                        gE = 1.0d0 / (Rsi + dx / (2.0d0 * k_grid(i, j))) / dx
                    else
                        gE = 1.0d0 / (Rse + dx / (2.0d0 * k_grid(i, j))) / dx
                    endif
                endif
                
                ! West neighbor link (i-1, j)
                if (fluid_type(i-1, j) == 0) then
                    gW = (2.0d0 * k_grid(i, j) * k_grid(i-1, j)) / (k_grid(i, j) + k_grid(i-1, j)) / (dx**2)
                else if (fluid_type(i-1, j) == 1) then
                    gW = 1.0d0 / (Rsi + dx / (2.0d0 * k_grid(i, j))) / dx
                else
                    gW = 1.0d0 / (Rse + dx / (2.0d0 * k_grid(i, j))) / dx
                endif
                
                ! North neighbor link (i, j+1)
                if (j == Ny) then
                    ! top boundary is adiabatic
                    gN = 0.0d0
                else
                    if (fluid_type(i, j+1) == 0) then
                        gN = (2.0d0 * k_grid(i, j) * k_grid(i, j+1)) / (k_grid(i, j) + k_grid(i, j+1)) / (dy**2)
                    else if (fluid_type(i, j+1) == 1) then
                        gN = 1.0d0 / (Rsi + dy / (2.0d0 * k_grid(i, j))) / dy
                    else
                        gN = 1.0d0 / (Rse + dy / (2.0d0 * k_grid(i, j))) / dy
                    endif
                endif
                
                ! South neighbor link (i, j-1)
                if (j == 1) then
                    gS = 0.0d0
                else
                    if (fluid_type(i, j-1) == 0) then
                        gS = (2.0d0 * k_grid(i, j) * k_grid(i, j-1)) / (k_grid(i, j) + k_grid(i, j-1)) / (dy**2)
                    else if (fluid_type(i, j-1) == 1) then
                        gS = 1.0d0 / (Rsi + dy / (2.0d0 * k_grid(i, j))) / dy
                    else
                        gS = 1.0d0 / (Rse + dy / (2.0d0 * k_grid(i, j))) / dy
                    endif
                endif
                
                denom = gE + gW + gN + gS
                if (denom > 0.0d0) then
                    T_new(i, j) = (gE * T(i+1, j) + gW * T(i-1, j) + gN * T(i, j+1) + gS * T(i, j-1)) / denom
                    
                    error = abs(T_new(i, j) - T(i, j))
                    if (error > max_diff) max_diff = error
                endif
            enddo
        enddo
        
        ! Update T array
        do i = 1, Nx
            do j = 1, Ny
                if (fluid_type(i, j) == 0) then
                    T(i, j) = T_new(i, j)
                endif
            enddo
        enddo
        
        ! Check convergence
        if (max_diff < 1.0d-7) exit
    enddo
    
    ! ==========================================================================
    ! Heat flow integrations & 1D reference calculation
    ! ==========================================================================
    ! 1D Wall R-value
    wall_R_1D = Rsi + Rse
    do i = 1, num_layers
        wall_R_1D = wall_R_1D + (layer_thickness(i) / 1000.0d0) / layer_k(i)
    enddo
    wall_U_1D = 1.0d0 / wall_R_1D
    
    ! Integrate interior heat flow q_2D (W/m)
    q_2D = 0.0d0
    T_si_min = Ti
    
    if (bridge_type == 2) then
        ! L-shaped wall corner heat flow integration on the solid boundary cells
        ! adjacent to the room (i=20 for j > 20 and j=20 for i > 20)
        do j = 21, Ny
            q_2D = q_2D + (Ti - T(20, j)) / (Rsi + dx / (2.0d0 * k_grid(20, j))) * dy
            if (T(20, j) < T_si_min) T_si_min = T(20, j)
        enddo
        do i = 21, Nx
            q_2D = q_2D + (Ti - T(i, 20)) / (Rsi + dy / (2.0d0 * k_grid(i, 20))) * dx
            if (T(i, 20) < T_si_min) T_si_min = T(i, 20)
        enddo
        ! corner reference 1D flow based on internal lengths
        q_1D = wall_U_1D * (dble(Nx - 20) * dx + dble(Ny - 20) * dy) * (Ti - Te)
        
    else
        ! Rectangular domain heat flow integration at i=1
        do j = 1, Ny
            q_2D = q_2D + (Ti - T(1, j)) / (Rsi + dx / (2.0d0 * k_grid(1, j))) * dy
            if (T(1, j) < T_si_min) T_si_min = T(1, j)
        enddo
        q_1D = wall_U_1D * H_domain * (Ti - Te)
    endif
    
    ! Transmittance calculations
    psi = 0.0d0
    chi = 0.0d0
    
    if (bridge_type == 5) then
        ! Point thermal bridge (Column penetration)
        ! Scale to W/K using the 2D plane area slice
        chi = (q_2D * H_domain - q_1D) / (Ti - Te)
        Q_bridge = chi * (Ti - Te)
    else
        ! Linear thermal bridges (Balcony slab, Corner, Window joint, Floor junction)
        psi = (q_2D - q_1D) / (Ti - Te)
        Q_bridge = psi * bridge_axis_length * (Ti - Te)
    endif
    
    ! Condensation factor f_Rsi
    f_Rsi = (T_si_min - Te) / (Ti - Te)
    
    ! ==========================================================================
    ! Output report
    ! ==========================================================================
    print *, "=========================================================================="
    print *, "             THERMAL BRIDGE 2D FDM ANALYSIS REPORT (ISO 10211)            "
    print *, "=========================================================================="
    select case (bridge_type)
    case (1)
        print *, " Bridge Type:             Balcony slab penetration (Linear)"
    case (2)
        print *, " Bridge Type:             Wall Corner junction (Linear)"
    case (3)
        print *, " Bridge Type:             Window frame joint (Linear)"
    case (4)
        print *, " Bridge Type:             Floor slab junction (Linear)"
    case (5)
        print *, " Bridge Type:             Column penetration (Point)"
    end select
    
    print "(A, F8.2, A)", " Inside Temperature (Ti): ", Ti, " C"
    print "(A, F8.2, A)", " Outside Temperature (Te):", Te, " C"
    print "(A, F8.4)", " Inside Resistance (Rsi): ", Rsi
    print "(A, F8.4)", " Outside Resistance (Rse):", Rse
    print *
    print *, "--------------------------------------------------------------------------"
    print *, " UNDISTURBED 1D WALL PERFORMANCE"
    print *, "--------------------------------------------------------------------------"
    print "(A, F10.3, A)", " Total Wall Thickness:    ", W_domain * 1000.0d0, " mm"
    print "(A, F10.4, A)", " 1D Thermal Resistance:   ", wall_R_1D, " m2-K/W"
    print "(A, F10.4, A)", " 1D Heat Transmittance U: ", wall_U_1D, " W/m2-K"
    print *
    print *, "--------------------------------------------------------------------------"
    print *, " 2D NUMERICAL THERMAL BRIDGE RESULTS"
    print *, "--------------------------------------------------------------------------"
    print "(A, F10.3, A)", " Total 2D Heat Flow:      ", q_2D, " W/m"
    print "(A, F10.3, A)", " Undisturbed 1D Heat Flow:", q_1D, " W/m"
    
    if (bridge_type == 5) then
        print "(A, F10.5, A)", " Point Transmittance (chi):", chi, " W/K"
    else
        print "(A, F10.5, A)", " Linear Transmittance (psi):", psi, " W/(m-K)"
        print "(A, F10.2, A)", " Thermal Bridge Length:   ", bridge_axis_length, " m"
    endif
    
    print "(A, F10.2, A)", " Additional Heat Loss:    ", Q_bridge, " W"
    print "(A, F10.2, A)", " Min Inside Surface Temp: ", T_si_min, " C"
    print "(A, F10.4)", " Temp Factor (f_Rsi):      ", f_Rsi
    
    if (f_Rsi < 0.75d0) then
        print *, " WARNING: Condensation & mold growth risk detected (f_Rsi < 0.75)!"
    else
        print *, " Surface temperature factor is safe (f_Rsi >= 0.75)."
    endif
    print *, "=========================================================================="
    
end program thermal_bridge
