apply_type_bc_crs Subroutine

    subroutine apply_type_bc_crs(self, boundary_target, current_time, A, b, domain, mode)
        class(type_bc), intent(inout) :: self
        character(*), intent(in) :: boundary_target
        real(real64), intent(in) :: current_time
        type(type_crs), intent(inout), optional :: A
        real(real64), intent(inout) :: b(:)
        type(type_domain), intent(in) :: domain
        integer(int32), intent(in), optional :: mode

        integer(int32) :: i

        ! --------------------------------------------------------------------------
        ! 1st Pass: Apply all non-Dirichlet boundary conditions (e.g., Neumann, Adiabatic)
        ! --------------------------------------------------------------------------
        select case (trim(adjustl(boundary_target)))
        case ('thermal')
            do i = 1, self%num_boundaries
                if (allocated(self%bc(i)%t)) then
                    select type (bc => self%bc(i)%t)
                    type is (type_bc_thermal_dirichlet)
                        ! This is a Dirichlet BC, so we skip it in the first pass.
                        cycle
                    class default
                        ! This is any other type of BC, apply it now.
                        if (present(A)) then
                            ! Apply the BC using CRS matrix format.
                            if (present(mode)) then
                                call bc%apply_crs(current_time=current_time, &
                                                  A=A, &
                                                  b=b, &
                                                  domain=domain, &
                                                  mode=mode)
                            else
                                call bc%apply_crs(current_time=current_time, &
                                                  b=b, &
                                                  domain=domain)
                            end if
                        else
                            ! Apply the BC without matrix A.
                            if (present(mode)) then
                                call bc%apply_crs(current_time=current_time, &
                                                  b=b, &
                                                  domain=domain, &
                                                  mode=mode)
                            else
                                call bc%apply_crs(current_time=current_time, &
                                                  b=b, &
                                                  domain=domain)
                            end if
                        end if
                    end select
                end if
            end do

            ! --------------------------------------------------------------------------
            ! 2nd Pass: Apply ONLY the Dirichlet boundary conditions
            ! --------------------------------------------------------------------------
            do i = 1, self%num_boundaries
                if (allocated(self%bc(i)%t)) then
                    select type (bc => self%bc(i)%t)
                    type is (type_bc_thermal_dirichlet)
                        ! This is a Dirichlet BC, apply it in the final pass.
                        if (present(A)) then
                            ! Apply the BC using CRS matrix format.
                            if (present(mode)) then
                                call bc%apply_crs(current_time=current_time, &
                                                  A=A, &
                                                  b=b, &
                                                  domain=domain, &
                                                  mode=mode)
                            else
                                call bc%apply_crs(current_time=current_time, &
                                                  b=b, &
                                                  domain=domain)
                            end if
                        else
                            ! Apply the BC without matrix A.
                            if (present(mode)) then
                                call bc%apply_crs(current_time=current_time, &
                                                  b=b, &
                                                  domain=domain, &
                                                  mode=mode)
                            else
                                call bc%apply_crs(current_time=current_time, &
                                                  b=b, &
                                                  domain=domain)
                            end if
                        end if
                    class default
                        ! All other types were handled in the first pass, so do nothing.
                        cycle
                    end select
                end if
            end do
        end select

    end subroutine apply_type_bc_crs