From c0b4be8e4d898a4e913abec33a46243353a05b95 Mon Sep 17 00:00:00 2001
From: Maarten van Ormondt <maarten.vanormondt@deltares.nl>
Date: Fri, 10 Apr 2026 16:15:54 +0200
Subject: [PATCH 1/6] Add building drainage module and integrate

Introduce a new sfincs_buildings module to model building drainage: collects rainfall on building footprints, creates thin dams along edges, optionally stores detention volume, and distributes runoff to perimeter or gutter outlet cells. Added building data and input options (bldfile, bprfile) to sfincs_data and sfincs_input, updated build scripts (Makefile.am and vfproj) to include the new source, and integrated setup and per-timestep update calls (setup_buildings, update_buildings) into sfincs_lib. The module provides routines for reading building and property files, identifying interior/perimeter/gutter cells, creating thin dams via quadtree intersections, accumulating rainfall, and redistributing water each timestep.
---
 source/sfincs_lib/sfincs_lib.vfproj |   1 +
 source/src/Makefile.am              |   1 +
 source/src/sfincs_buildings.f90     | 729 ++++++++++++++++++++++++++++
 source/src/sfincs_data.f90          |   7 +
 source/src/sfincs_input.f90         |   6 +
 source/src/sfincs_lib.f90           |  11 +
 6 files changed, 755 insertions(+)
 create mode 100644 source/src/sfincs_buildings.f90

diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj
index 34bf8f520..83c3738bf 100644
--- a/source/sfincs_lib/sfincs_lib.vfproj
+++ b/source/sfincs_lib/sfincs_lib.vfproj
@@ -60,6 +60,7 @@
 			<File RelativePath="..\src\sfincs_crosssections.f90"/>
 			<File RelativePath="..\src\sfincs_data.f90">
 			</File>
+			<File RelativePath="..\src\sfincs_buildings.f90"/>
 			<File RelativePath="..\src\sfincs_date.f90"/>
 			<File RelativePath="..\src\sfincs_discharges.f90"/>
 			<File RelativePath="..\src\sfincs_domain.f90">
diff --git a/source/src/Makefile.am b/source/src/Makefile.am
index 891652f08..a3f6d13af 100644
--- a/source/src/Makefile.am
+++ b/source/src/Makefile.am
@@ -19,6 +19,7 @@ libsfincs_la_SOURCES = \
 	sfincs_date.f90 \
 	sfincs_spiderweb.f90 \
 	sfincs_data.f90 \
+	sfincs_buildings.f90 \
     ../third_party_open/Delft3D/astro.f90 \
 	../third_party_open/utils/geometry.f90 \
 	sfincs_error.f90 \		
diff --git a/source/src/sfincs_buildings.f90 b/source/src/sfincs_buildings.f90
new file mode 100644
index 000000000..2799b081f
--- /dev/null
+++ b/source/src/sfincs_buildings.f90
@@ -0,0 +1,729 @@
+module sfincs_buildings
+   !
+   ! Building drainage module for SFINCS
+   ! 1. Collect rainfall on building footprint
+   ! 2. Block flow with thin dams on edges
+   ! 3. Redistribute collected water to perimeter cells OR gutter outlets
+   ! 4. Optional detention volume per building
+   ! 5. Optional gutter system with configurable spacing
+   ! IMDC, 2025 (restructured by Deltares)
+   !
+   implicit none
+   !
+   type building_type
+      character(len=256) :: name
+      integer :: npoints
+      real*8, allocatable :: x(:)
+      real*8, allocatable :: y(:)
+      integer, allocatable :: cells_inside(:)
+      integer :: ncells_inside
+      integer, allocatable :: cells_perimeter(:)
+      integer :: ncells_perimeter
+      integer, allocatable :: cells_gutter(:)
+      integer :: ncells_gutter
+      real*8 :: total_area
+      real*8 :: accumulated_rainfall
+      real*8 :: accumulated_runoff
+      real*8 :: detention_volume
+      real*8 :: current_detention
+      real*8 :: gutter_spacing
+      real*8 :: perimeter_length
+   end type building_type
+   !
+   integer :: nbuildings = 0
+   type(building_type), allocatable :: buildings(:)
+   !
+   logical, allocatable :: is_building_cell(:)
+   logical, allocatable :: is_perimeter_cell(:)
+   logical, allocatable :: is_gutter_cell(:)
+   integer, allocatable :: building_id(:)
+   !
+   real*8 :: total_rain_on_buildings = 0.0d0
+   real*8 :: total_runoff_from_buildings = 0.0d0
+   real*8 :: total_detention_stored = 0.0d0
+   !
+   logical :: has_building_properties = .false.
+   logical :: use_building_detention = .false.
+   logical :: use_building_gutters = .false.
+   !
+contains
+   !
+   subroutine setup_buildings()
+   !
+   ! Read files, identify cells, create thin dams.
+   ! Called once from sfincs_domain after mesh initialization.
+   !
+   use sfincs_data, only: has_buildings, bldfile, bprfile
+   !
+   implicit none
+   !
+   integer :: ios
+   !
+   write(*,'(A)') 'Info    : Setting up buildings ...'
+   !
+   call initialize_building_masks()
+   !
+   call read_building_file(bldfile, ios)
+   !
+   if (ios /= 0) then
+      write(*,'(A)') 'ERROR: Failed to read building file'
+      has_buildings = .false.
+      return
+   end if
+   !
+   if (len_trim(bprfile) > 0 .and. trim(bprfile) /= 'none') then
+      call read_building_properties_file(bprfile, ios)
+      if (ios /= 0) then
+         write(*,'(A)') 'WARNING: Failed to read building properties file'
+         write(*,'(A)') '         Using default values (no detention, uniform distribution)'
+      end if
+   end if
+   !
+   call identify_building_cells()
+   call identify_perimeter_cells()
+   !
+   if (use_building_gutters) then
+      call identify_gutter_cells()
+   end if
+   !
+   call create_building_thin_dams()
+   !
+   write(*,'(A)') 'Info    : Building setup complete'
+   !
+   end subroutine setup_buildings
+   !
+   !
+   subroutine update_buildings(precip, zs, dt)
+   !
+   ! Accumulate rainfall and redistribute. Called each timestep.
+   !
+   implicit none
+   !
+   real*4, intent(inout) :: precip(:)
+   real*8, intent(inout) :: zs(:)
+   real*8, intent(in) :: dt
+   !
+   if (nbuildings == 0) return
+   !
+   call accumulate_building_rainfall(precip, dt)
+   call redistribute_building_water(zs, dt)
+   !
+   end subroutine update_buildings
+   !
+   !
+   subroutine initialize_building_masks()
+   !
+   use sfincs_data, only: mmax, nmax
+   !
+   implicit none
+   !
+   integer :: ncells
+   !
+   ncells = mmax * nmax
+   !
+   if (.not. allocated(is_building_cell)) then
+      allocate(is_building_cell(ncells))
+      allocate(is_perimeter_cell(ncells))
+      allocate(is_gutter_cell(ncells))
+      allocate(building_id(ncells))
+   end if
+   !
+   is_building_cell = .false.
+   is_perimeter_cell = .false.
+   is_gutter_cell = .false.
+   building_id = 0
+   !
+   end subroutine initialize_building_masks
+   !
+   !
+   subroutine read_building_file(filename, ios_out)
+   !
+   implicit none
+   !
+   character(len=*), intent(in) :: filename
+   integer, intent(out) :: ios_out
+   !
+   integer :: iunit, ios_bld, i, ibld
+   character(len=1024) :: line
+   character(len=256) :: name
+   integer :: npoints, ncols
+   real*8 :: xtmp, ytmp
+   !
+   ios_out = 0
+   !
+   open(newunit=iunit, file=trim(filename), status='old', action='read', iostat=ios_bld)
+   if (ios_bld /= 0) then
+      write(*,'(A,A)') 'ERROR: Could not open building file: ', trim(filename)
+      ios_out = ios_bld
+      return
+   end if
+   !
+   write(*,'(A,A)') 'Info    : Reading building file: ', trim(filename)
+   !
+   ! Count buildings
+   !
+   nbuildings = 0
+   do
+      read(iunit, '(A)', iostat=ios_bld) line
+      if (ios_bld /= 0) exit
+      if (len_trim(line) == 0) cycle
+      if (line(1:1) == '!') cycle
+      read(line, *, iostat=ios_bld) name, npoints, ncols
+      if (ios_bld == 0 .and. npoints > 0 .and. ncols >= 2) then
+         nbuildings = nbuildings + 1
+         do i = 1, npoints
+            read(iunit, *, iostat=ios_bld)
+            if (ios_bld /= 0) exit
+         end do
+      end if
+   end do
+   !
+   if (nbuildings == 0) then
+      write(*,'(A)') 'WARNING: No buildings found in file'
+      close(iunit)
+      return
+   end if
+   !
+   write(*,'(A,I0,A)') 'Info    : Found ', nbuildings, ' buildings'
+   !
+   allocate(buildings(nbuildings))
+   !
+   rewind(iunit)
+   !
+   ibld = 0
+   do
+      read(iunit, '(A)', iostat=ios_bld) line
+      if (ios_bld /= 0) exit
+      if (len_trim(line) == 0) cycle
+      if (line(1:1) == '!') cycle
+      read(line, *, iostat=ios_bld) name, npoints, ncols
+      if (ios_bld /= 0 .or. npoints <= 0 .or. ncols < 2) cycle
+      !
+      ibld = ibld + 1
+      if (ibld > nbuildings) exit
+      !
+      buildings(ibld)%name = trim(name)
+      buildings(ibld)%npoints = npoints
+      allocate(buildings(ibld)%x(npoints))
+      allocate(buildings(ibld)%y(npoints))
+      !
+      do i = 1, npoints
+         read(iunit, *, iostat=ios_bld) xtmp, ytmp
+         if (ios_bld /= 0) then
+            write(*,'(A,A)') 'ERROR: Reading coordinates for building ', trim(name)
+            ios_out = ios_bld
+            close(iunit)
+            return
+         end if
+         buildings(ibld)%x(i) = xtmp
+         buildings(ibld)%y(i) = ytmp
+      end do
+      !
+      if (buildings(ibld)%x(1) /= buildings(ibld)%x(npoints) .or. &
+          buildings(ibld)%y(1) /= buildings(ibld)%y(npoints)) then
+         write(*,'(A,A,A)') 'WARNING: Building ', trim(name), ' polygon not closed'
+      end if
+      !
+      buildings(ibld)%accumulated_rainfall = 0.0d0
+      buildings(ibld)%accumulated_runoff = 0.0d0
+      buildings(ibld)%ncells_inside = 0
+      buildings(ibld)%ncells_perimeter = 0
+      buildings(ibld)%ncells_gutter = 0
+      buildings(ibld)%total_area = 0.0d0
+      buildings(ibld)%perimeter_length = 0.0d0
+      buildings(ibld)%detention_volume = 0.0d0
+      buildings(ibld)%current_detention = 0.0d0
+      buildings(ibld)%gutter_spacing = 0.0d0
+      !
+   end do
+   !
+   close(iunit)
+   !
+   write(*,'(A,I0,A)') 'Info    : Successfully read ', nbuildings, ' buildings'
+   !
+   end subroutine read_building_file
+   !
+   !
+   subroutine read_building_properties_file(filename, ios_out)
+   !
+   implicit none
+   !
+   character(len=*), intent(in) :: filename
+   integer, intent(out) :: ios_out
+   !
+   integer :: iunit, ios_bpr, ibld, nmatched
+   character(len=1024) :: line
+   character(len=256) :: name
+   real*8 :: det_vol, gut_spacing
+   !
+   ios_out = 0
+   nmatched = 0
+   !
+   open(newunit=iunit, file=trim(filename), status='old', action='read', iostat=ios_bpr)
+   if (ios_bpr /= 0) then
+      write(*,'(A,A)') 'ERROR: Could not open building properties file: ', trim(filename)
+      ios_out = ios_bpr
+      return
+   end if
+   !
+   write(*,'(A,A)') 'Info    : Reading building properties file: ', trim(filename)
+   !
+   do
+      read(iunit, '(A)', iostat=ios_bpr) line
+      if (ios_bpr /= 0) exit
+      if (len_trim(line) == 0) cycle
+      if (line(1:1) == '!' .or. line(1:1) == '#') cycle
+      !
+      read(line, *, iostat=ios_bpr) name, det_vol, gut_spacing
+      if (ios_bpr /= 0) then
+         write(*,'(A,A)') 'WARNING: Could not parse line: ', trim(line)
+         cycle
+      end if
+      !
+      do ibld = 1, nbuildings
+         if (trim(buildings(ibld)%name) == trim(name)) then
+            buildings(ibld)%detention_volume = det_vol
+            buildings(ibld)%gutter_spacing = gut_spacing
+            nmatched = nmatched + 1
+            if (det_vol > 0.0d0) use_building_detention = .true.
+            if (gut_spacing > 0.0d0) use_building_gutters = .true.
+            exit
+         end if
+      end do
+      !
+   end do
+   !
+   close(iunit)
+   !
+   write(*,'(A,I0,A,I0,A)') 'Info    : Matched properties for ', nmatched, ' of ', nbuildings, ' buildings'
+   if (use_building_detention) write(*,'(A)') 'Info    : Building detention volumes enabled'
+   if (use_building_gutters) write(*,'(A)') 'Info    : Building gutter system enabled'
+   !
+   end subroutine read_building_properties_file
+   !
+   !
+   subroutine identify_building_cells()
+   !
+   use sfincs_data, only: mmax, nmax, np, x0, y0, dx, dy, cosrot, sinrot, kcs
+   !
+   implicit none
+   !
+   integer :: m, n, ibld, ip, ncells
+   real*8 :: xc, yc
+   logical :: inside
+   integer, allocatable :: temp_cells(:)
+   !
+   write(*,'(A)') 'Info    : Identifying building cells ...'
+   !
+   ncells = mmax * nmax
+   !
+   do ibld = 1, nbuildings
+      !
+      allocate(temp_cells(ncells))
+      buildings(ibld)%ncells_inside = 0
+      !
+      do n = 1, nmax
+         do m = 1, mmax
+            !
+            xc = x0 + cosrot*(dble(m) - 0.5d0)*dble(dx) - sinrot*(dble(n) - 0.5d0)*dble(dy)
+            yc = y0 + sinrot*(dble(m) - 0.5d0)*dble(dx) + cosrot*(dble(n) - 0.5d0)*dble(dy)
+            !
+            call point_in_polygon(xc, yc, buildings(ibld)%x, buildings(ibld)%y, &
+                                  buildings(ibld)%npoints, inside)
+            !
+            if (inside) then
+               ip = (m - 1)*nmax + n
+               if (ip < 1 .or. ip > np) cycle
+               !
+               is_building_cell(ip) = .true.
+               building_id(ip) = ibld
+               kcs(ip) = 0
+               !
+               buildings(ibld)%ncells_inside = buildings(ibld)%ncells_inside + 1
+               temp_cells(buildings(ibld)%ncells_inside) = ip
+            end if
+            !
+         end do
+      end do
+      !
+      if (buildings(ibld)%ncells_inside > 0) then
+         allocate(buildings(ibld)%cells_inside(buildings(ibld)%ncells_inside))
+         buildings(ibld)%cells_inside = temp_cells(1:buildings(ibld)%ncells_inside)
+      end if
+      !
+      deallocate(temp_cells)
+      !
+      buildings(ibld)%total_area = dble(buildings(ibld)%ncells_inside) * dx * dy
+      call calculate_perimeter_length(ibld)
+      !
+   end do
+   !
+   write(*,'(A)') 'Info    : Building cell identification complete'
+   !
+   end subroutine identify_building_cells
+   !
+   !
+   subroutine calculate_perimeter_length(ibld)
+   !
+   implicit none
+   !
+   integer, intent(in) :: ibld
+   integer :: i
+   real*8 :: length, dx_seg, dy_seg
+   !
+   length = 0.0d0
+   do i = 1, buildings(ibld)%npoints - 1
+      dx_seg = buildings(ibld)%x(i+1) - buildings(ibld)%x(i)
+      dy_seg = buildings(ibld)%y(i+1) - buildings(ibld)%y(i)
+      length = length + sqrt(dx_seg*dx_seg + dy_seg*dy_seg)
+   end do
+   buildings(ibld)%perimeter_length = length
+   !
+   end subroutine calculate_perimeter_length
+   !
+   !
+   subroutine point_in_polygon(xp, yp, xv, yv, nv, inside)
+   !
+   implicit none
+   !
+   real*8, intent(in) :: xp, yp
+   real*8, intent(in) :: xv(:), yv(:)
+   integer, intent(in) :: nv
+   logical, intent(out) :: inside
+   !
+   integer :: i, j
+   logical :: intersect
+   !
+   inside = .false.
+   j = nv
+   do i = 1, nv
+      intersect = ((yv(i) > yp) .neqv. (yv(j) > yp)) .and. &
+                  (xp < (xv(j) - xv(i))*(yp - yv(i))/(yv(j) - yv(i)) + xv(i))
+      if (intersect) inside = .not. inside
+      j = i
+   end do
+   !
+   end subroutine point_in_polygon
+   !
+   !
+   subroutine identify_perimeter_cells()
+   !
+   use sfincs_data, only: mmax, nmax
+   !
+   implicit none
+   !
+   integer :: ibld, ic, ip, m, n, mp, np_local, ipp
+   integer :: dm, dn
+   integer, allocatable :: temp_perimeter(:)
+   integer :: ncells
+   !
+   write(*,'(A)') 'Info    : Identifying perimeter cells ...'
+   !
+   ncells = mmax * nmax
+   !
+   do ibld = 1, nbuildings
+      !
+      allocate(temp_perimeter(ncells))
+      buildings(ibld)%ncells_perimeter = 0
+      !
+      do ic = 1, buildings(ibld)%ncells_inside
+         ip = buildings(ibld)%cells_inside(ic)
+         m = (ip - 1)/nmax + 1
+         n = ip - (m - 1)*nmax
+         !
+         do dn = -1, 1
+            do dm = -1, 1
+               if (abs(dm) + abs(dn) /= 1) cycle
+               mp = m + dm
+               np_local = n + dn
+               if (mp < 1 .or. mp > mmax) cycle
+               if (np_local < 1 .or. np_local > nmax) cycle
+               !
+               ipp = (mp - 1)*nmax + np_local
+               !
+               if (.not. is_building_cell(ipp) .and. .not. is_perimeter_cell(ipp)) then
+                  buildings(ibld)%ncells_perimeter = buildings(ibld)%ncells_perimeter + 1
+                  temp_perimeter(buildings(ibld)%ncells_perimeter) = ipp
+                  is_perimeter_cell(ipp) = .true.
+               end if
+            end do
+         end do
+      end do
+      !
+      if (buildings(ibld)%ncells_perimeter > 0) then
+         allocate(buildings(ibld)%cells_perimeter(buildings(ibld)%ncells_perimeter))
+         buildings(ibld)%cells_perimeter = temp_perimeter(1:buildings(ibld)%ncells_perimeter)
+      end if
+      !
+      deallocate(temp_perimeter)
+      !
+   end do
+   !
+   write(*,'(A)') 'Info    : Perimeter cell identification complete'
+   !
+   end subroutine identify_perimeter_cells
+   !
+   !
+   subroutine identify_gutter_cells()
+   !
+   use sfincs_data, only: mmax, nmax, x0, y0, dx, dy, cosrot, sinrot
+   !
+   implicit none
+   !
+   integer :: ibld, ic, ip, m, n, j, min_idx
+   integer, allocatable :: temp_gutter(:), sorted_perimeter(:)
+   real*8 :: spacing, cum_dist, next_gutter_dist
+   real*8 :: xc, yc, xp_prev, yp_prev
+   real*8 :: x_centroid, y_centroid
+   real*8, allocatable :: angles(:)
+   real*8 :: temp_angle
+   integer :: temp_ip
+   logical :: first_cell
+   !
+   write(*,'(A)') 'Info    : Identifying gutter outlet cells ...'
+   !
+   do ibld = 1, nbuildings
+      !
+      if (buildings(ibld)%gutter_spacing <= 0.0d0) then
+         buildings(ibld)%ncells_gutter = 0
+         cycle
+      end if
+      !
+      if (buildings(ibld)%ncells_perimeter == 0) then
+         buildings(ibld)%ncells_gutter = 0
+         cycle
+      end if
+      !
+      spacing = buildings(ibld)%gutter_spacing
+      !
+      ! Building centroid from polygon vertices (exclude closing point)
+      !
+      if (buildings(ibld)%npoints > 1) then
+         x_centroid = sum(buildings(ibld)%x(1:buildings(ibld)%npoints-1)) / &
+                      dble(buildings(ibld)%npoints - 1)
+         y_centroid = sum(buildings(ibld)%y(1:buildings(ibld)%npoints-1)) / &
+                      dble(buildings(ibld)%npoints - 1)
+      else
+         x_centroid = buildings(ibld)%x(1)
+         y_centroid = buildings(ibld)%y(1)
+      end if
+      !
+      ! Sort perimeter cells by angle from centroid
+      !
+      allocate(sorted_perimeter(buildings(ibld)%ncells_perimeter))
+      allocate(angles(buildings(ibld)%ncells_perimeter))
+      !
+      do ic = 1, buildings(ibld)%ncells_perimeter
+         ip = buildings(ibld)%cells_perimeter(ic)
+         sorted_perimeter(ic) = ip
+         m = (ip - 1)/nmax + 1
+         n = ip - (m - 1)*nmax
+         xc = x0 + cosrot*(dble(m) - 0.5d0)*dble(dx) - sinrot*(dble(n) - 0.5d0)*dble(dy)
+         yc = y0 + sinrot*(dble(m) - 0.5d0)*dble(dx) + cosrot*(dble(n) - 0.5d0)*dble(dy)
+         angles(ic) = atan2(yc - y_centroid, xc - x_centroid)
+      end do
+      !
+      ! Selection sort by angle
+      !
+      do ic = 1, buildings(ibld)%ncells_perimeter - 1
+         min_idx = ic
+         do j = ic + 1, buildings(ibld)%ncells_perimeter
+            if (angles(j) < angles(min_idx)) min_idx = j
+         end do
+         if (min_idx /= ic) then
+            temp_angle = angles(ic)
+            angles(ic) = angles(min_idx)
+            angles(min_idx) = temp_angle
+            temp_ip = sorted_perimeter(ic)
+            sorted_perimeter(ic) = sorted_perimeter(min_idx)
+            sorted_perimeter(min_idx) = temp_ip
+         end if
+      end do
+      !
+      ! Walk sorted perimeter and place gutters at intervals
+      !
+      allocate(temp_gutter(buildings(ibld)%ncells_perimeter))
+      buildings(ibld)%ncells_gutter = 0
+      cum_dist = 0.0d0
+      next_gutter_dist = spacing / 2.0d0
+      first_cell = .true.
+      !
+      do ic = 1, buildings(ibld)%ncells_perimeter
+         ip = sorted_perimeter(ic)
+         m = (ip - 1)/nmax + 1
+         n = ip - (m - 1)*nmax
+         xc = x0 + cosrot*(dble(m) - 0.5d0)*dble(dx) - sinrot*(dble(n) - 0.5d0)*dble(dy)
+         yc = y0 + sinrot*(dble(m) - 0.5d0)*dble(dx) + cosrot*(dble(n) - 0.5d0)*dble(dy)
+         !
+         if (.not. first_cell) then
+            cum_dist = cum_dist + sqrt((xc - xp_prev)**2 + (yc - yp_prev)**2)
+         end if
+         !
+         if (cum_dist >= next_gutter_dist) then
+            buildings(ibld)%ncells_gutter = buildings(ibld)%ncells_gutter + 1
+            temp_gutter(buildings(ibld)%ncells_gutter) = ip
+            is_gutter_cell(ip) = .true.
+            next_gutter_dist = next_gutter_dist + spacing
+         end if
+         !
+         xp_prev = xc
+         yp_prev = yc
+         first_cell = .false.
+      end do
+      !
+      ! Ensure at least one gutter per building
+      !
+      if (buildings(ibld)%ncells_gutter == 0) then
+         buildings(ibld)%ncells_gutter = 1
+         temp_gutter(1) = sorted_perimeter(1)
+         is_gutter_cell(temp_gutter(1)) = .true.
+      end if
+      !
+      if (buildings(ibld)%ncells_gutter > 0) then
+         allocate(buildings(ibld)%cells_gutter(buildings(ibld)%ncells_gutter))
+         buildings(ibld)%cells_gutter = temp_gutter(1:buildings(ibld)%ncells_gutter)
+      end if
+      !
+      deallocate(temp_gutter, sorted_perimeter, angles)
+      !
+   end do
+   !
+   write(*,'(A)') 'Info    : Gutter cell identification complete'
+   !
+   end subroutine identify_gutter_cells
+   !
+   !
+   subroutine create_building_thin_dams()
+   !
+   use sfincs_data, only: kcuv
+   use quadtree, only: find_uv_points_intersected_by_polyline
+   !
+   implicit none
+   !
+   integer :: ibld, iuv, nrows, nr_points, indx
+   real*4, allocatable :: xthd(:), ythd(:)
+   integer, allocatable :: uv_indices(:), vertices(:)
+   !
+   write(*,'(A,I0,A)') 'Info    : Processing ', nbuildings, ' buildings as thin dams'
+   !
+   do ibld = 1, nbuildings
+      !
+      nrows = buildings(ibld)%npoints
+      allocate(xthd(nrows), ythd(nrows))
+      !
+      xthd = real(buildings(ibld)%x, 4)
+      ythd = real(buildings(ibld)%y, 4)
+      !
+      call find_uv_points_intersected_by_polyline(uv_indices, vertices, nr_points, xthd, ythd, nrows)
+      !
+      do iuv = 1, nr_points
+         indx = uv_indices(iuv)
+         kcuv(indx) = 0
+      end do
+      !
+      deallocate(xthd, ythd)
+      !
+   end do
+   !
+   write(*,'(A,I0,A)') 'Info    : Processed ', nbuildings, ' buildings as thin dams'
+   !
+   end subroutine create_building_thin_dams
+   !
+   !
+   subroutine accumulate_building_rainfall(precip, dt)
+   !
+   use sfincs_data, only: dx, dy, netprcp
+   !
+   implicit none
+   !
+   real*4, intent(inout) :: precip(:)
+   real*8, intent(in) :: dt
+   !
+   integer :: ibld, ic, ip
+   real*8 :: rain_volume
+   !
+   do ibld = 1, nbuildings
+      do ic = 1, buildings(ibld)%ncells_inside
+         ip = buildings(ibld)%cells_inside(ic)
+         !
+         if (precip(ip) > 0.0) then
+            rain_volume = dble(precip(ip)) * dt * dx * dy
+            buildings(ibld)%accumulated_rainfall = buildings(ibld)%accumulated_rainfall + rain_volume
+            total_rain_on_buildings = total_rain_on_buildings + rain_volume
+         end if
+         !
+         precip(ip) = 0.0
+         if (allocated(netprcp)) netprcp(ip) = 0.0
+         !
+      end do
+   end do
+   !
+   end subroutine accumulate_building_rainfall
+   !
+   !
+   subroutine redistribute_building_water(zs, dt)
+   !
+   use sfincs_data, only: dx, dy
+   !
+   implicit none
+   !
+   real*8, intent(inout) :: zs(:)
+   real*8, intent(in) :: dt
+   !
+   integer :: ibld, ic, ip
+   real*8 :: volume_to_distribute, volume_per_cell, dh
+   real*8 :: cell_area, available_detention, volume_to_detention
+   !
+   cell_area = dx * dy
+   !
+   do ibld = 1, nbuildings
+      !
+      if (buildings(ibld)%accumulated_rainfall <= 1.0d-100) cycle
+      !
+      volume_to_distribute = buildings(ibld)%accumulated_rainfall
+      !
+      ! Detention storage
+      !
+      if (use_building_detention .and. buildings(ibld)%detention_volume > 0.0d0) then
+         available_detention = buildings(ibld)%detention_volume - buildings(ibld)%current_detention
+         if (available_detention > 0.0d0) then
+            volume_to_detention = min(volume_to_distribute, available_detention)
+            buildings(ibld)%current_detention = buildings(ibld)%current_detention + volume_to_detention
+            volume_to_distribute = volume_to_distribute - volume_to_detention
+            total_detention_stored = total_detention_stored + volume_to_detention
+         end if
+      end if
+      !
+      if (volume_to_distribute <= 1.0d-12) then
+         buildings(ibld)%accumulated_rainfall = 0.0d0
+         cycle
+      end if
+      !
+      ! Redistribute to gutters or perimeter
+      !
+      if (use_building_gutters .and. buildings(ibld)%ncells_gutter > 0) then
+         volume_per_cell = volume_to_distribute / dble(buildings(ibld)%ncells_gutter)
+         dh = volume_per_cell / cell_area
+         do ic = 1, buildings(ibld)%ncells_gutter
+            ip = buildings(ibld)%cells_gutter(ic)
+            zs(ip) = zs(ip) + real(dh, 4)
+         end do
+      else
+         if (buildings(ibld)%ncells_perimeter == 0) cycle
+         volume_per_cell = volume_to_distribute / dble(buildings(ibld)%ncells_perimeter)
+         dh = volume_per_cell / cell_area
+         do ic = 1, buildings(ibld)%ncells_perimeter
+            ip = buildings(ibld)%cells_perimeter(ic)
+            zs(ip) = zs(ip) + real(dh, 4)
+         end do
+      end if
+      !
+      total_runoff_from_buildings = total_runoff_from_buildings + volume_to_distribute
+      buildings(ibld)%accumulated_runoff = buildings(ibld)%accumulated_runoff + volume_to_distribute
+      buildings(ibld)%accumulated_rainfall = 0.0d0
+      !
+   end do
+   !
+   end subroutine redistribute_building_water
+   !
+end module sfincs_buildings
diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90
index 14066306a..a715a4ac4 100644
--- a/source/src/sfincs_data.f90
+++ b/source/src/sfincs_data.f90
@@ -814,6 +814,13 @@ module sfincs_data
       character*256, dimension(:), allocatable :: runup_gauge_name
       integer, dimension(:),     allocatable   :: runup_gauge_nrp
       !
+      !
+      ! Buildings
+      !
+      logical                                  :: has_buildings = .false.
+      character*256                            :: bldfile = ''
+      character*256                            :: bprfile = ''
+      !
       real*4 :: waveage
       !
       real*4 :: bathtub_dt
diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90
index 5aab1f182..d4243d1f5 100644
--- a/source/src/sfincs_input.f90
+++ b/source/src/sfincs_input.f90
@@ -298,6 +298,12 @@ subroutine read_sfincs_input()
    ! Limit to range (0,100)
    percdoneval = max(min(percdoneval,100), 0)
    !
+   ! Buildings
+   !
+   call read_char_input(500,'bldfile',bldfile,'none')
+   if (trim(bldfile) /= 'none') has_buildings = .true.
+   call read_char_input(500,'bprfile',bprfile,'none')
+   !
    ! Coupled SnapWave solver related
    call read_int_input(500,'snapwave_wind',iwind,0)   
    !
diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90
index 8e3707475..b0b27859a 100644
--- a/source/src/sfincs_lib.f90
+++ b/source/src/sfincs_lib.f90
@@ -14,6 +14,7 @@ module sfincs_lib
    use sfincs_meteo
    use sfincs_infiltration
    use sfincs_data
+   use sfincs_buildings, only: setup_buildings, update_buildings
    use sfincs_date
    use sfincs_output
    use sfincs_ncinput
@@ -155,6 +156,10 @@ function sfincs_initialize() result(ierr)
    !
    call read_structures()       ! Reads thd files and sets kcuv to zero where necessary
    !
+   if (has_buildings) then
+      call setup_buildings()    ! Add structures for buildings and set up building parameters
+   endif
+   !
    call read_boundary_data()    ! Reads bnd, bzs, etc files
    !
    call find_boundary_indices()
@@ -529,6 +534,12 @@ function sfincs_update(dtrange) result(ierr)
          !
          call update_meteo_forcing(t, dt, tloopwnd2)
          !
+         ! Buildings: accumulate rainfall and redistribute to perimeter
+         !
+         if (has_buildings) then
+            call update_buildings(prcp, zs, real(dt, 8))
+         endif
+         !
          ! Update infiltration
          !
          if (infiltration) then

From 41bf22ef77ff742495eb9028524cda2ee7c1afbc Mon Sep 17 00:00:00 2001
From: Maarten van Ormondt <maarten.vanormondt@deltares.nl>
Date: Fri, 10 Apr 2026 16:57:19 +0200
Subject: [PATCH 2/6] Use sfincs_log for building module logging

Replace direct write(*) calls in source/src/sfincs_buildings.f90 with calls to write_log by adding "use sfincs_log" and routing informational, warning and error messages through the logging API. Dynamic messages are formatted into a local log string (logstr) before being passed to write_log. This centralizes output handling for the buildings module and preserves existing message content/severity while decoupling from direct stdout writes.
---
 source/src/sfincs_buildings.f90 | 66 ++++++++++++++++++++-------------
 1 file changed, 40 insertions(+), 26 deletions(-)

diff --git a/source/src/sfincs_buildings.f90 b/source/src/sfincs_buildings.f90
index 2799b081f..ab555e068 100644
--- a/source/src/sfincs_buildings.f90
+++ b/source/src/sfincs_buildings.f90
@@ -8,6 +8,8 @@ module sfincs_buildings
    ! 5. Optional gutter system with configurable spacing
    ! IMDC, 2025 (restructured by Deltares)
    !
+   use sfincs_log
+   !
    implicit none
    !
    type building_type
@@ -59,14 +61,14 @@ subroutine setup_buildings()
    !
    integer :: ios
    !
-   write(*,'(A)') 'Info    : Setting up buildings ...'
+   call write_log('Info    : Setting up buildings ...', 1)
    !
    call initialize_building_masks()
    !
    call read_building_file(bldfile, ios)
    !
    if (ios /= 0) then
-      write(*,'(A)') 'ERROR: Failed to read building file'
+      call write_log('ERROR: Failed to read building file', 0)
       has_buildings = .false.
       return
    end if
@@ -74,8 +76,8 @@ subroutine setup_buildings()
    if (len_trim(bprfile) > 0 .and. trim(bprfile) /= 'none') then
       call read_building_properties_file(bprfile, ios)
       if (ios /= 0) then
-         write(*,'(A)') 'WARNING: Failed to read building properties file'
-         write(*,'(A)') '         Using default values (no detention, uniform distribution)'
+         call write_log('WARNING: Failed to read building properties file', 0)
+         call write_log('         Using default values (no detention, uniform distribution)', 0)
       end if
    end if
    !
@@ -88,7 +90,7 @@ subroutine setup_buildings()
    !
    call create_building_thin_dams()
    !
-   write(*,'(A)') 'Info    : Building setup complete'
+   call write_log('Info    : Building setup complete', 0)
    !
    end subroutine setup_buildings
    !
@@ -153,12 +155,14 @@ subroutine read_building_file(filename, ios_out)
    !
    open(newunit=iunit, file=trim(filename), status='old', action='read', iostat=ios_bld)
    if (ios_bld /= 0) then
-      write(*,'(A,A)') 'ERROR: Could not open building file: ', trim(filename)
+      write(logstr,'(A,A)') 'ERROR: Could not open building file: ', trim(filename)
+      call write_log(logstr, 0)
       ios_out = ios_bld
       return
    end if
    !
-   write(*,'(A,A)') 'Info    : Reading building file: ', trim(filename)
+   write(logstr,'(A,A)') 'Info    : Reading building file: ', trim(filename)
+   call write_log(logstr, 0)
    !
    ! Count buildings
    !
@@ -179,12 +183,13 @@ subroutine read_building_file(filename, ios_out)
    end do
    !
    if (nbuildings == 0) then
-      write(*,'(A)') 'WARNING: No buildings found in file'
+      call write_log('WARNING: No buildings found in file', 0)
       close(iunit)
       return
    end if
    !
-   write(*,'(A,I0,A)') 'Info    : Found ', nbuildings, ' buildings'
+   write(logstr,'(A,I0,A)') 'Info    : Found ', nbuildings, ' buildings'
+   call write_log(logstr, 0)
    !
    allocate(buildings(nbuildings))
    !
@@ -210,7 +215,8 @@ subroutine read_building_file(filename, ios_out)
       do i = 1, npoints
          read(iunit, *, iostat=ios_bld) xtmp, ytmp
          if (ios_bld /= 0) then
-            write(*,'(A,A)') 'ERROR: Reading coordinates for building ', trim(name)
+            write(logstr,'(A,A)') 'ERROR: Reading coordinates for building ', trim(name)
+            call write_log(logstr, 0)
             ios_out = ios_bld
             close(iunit)
             return
@@ -221,7 +227,8 @@ subroutine read_building_file(filename, ios_out)
       !
       if (buildings(ibld)%x(1) /= buildings(ibld)%x(npoints) .or. &
           buildings(ibld)%y(1) /= buildings(ibld)%y(npoints)) then
-         write(*,'(A,A,A)') 'WARNING: Building ', trim(name), ' polygon not closed'
+         write(logstr,'(A,A,A)') 'WARNING: Building ', trim(name), ' polygon not closed'
+         call write_log(logstr, 0)
       end if
       !
       buildings(ibld)%accumulated_rainfall = 0.0d0
@@ -239,7 +246,8 @@ subroutine read_building_file(filename, ios_out)
    !
    close(iunit)
    !
-   write(*,'(A,I0,A)') 'Info    : Successfully read ', nbuildings, ' buildings'
+   write(logstr,'(A,I0,A)') 'Info    : Successfully read ', nbuildings, ' buildings'
+   call write_log(logstr, 0)
    !
    end subroutine read_building_file
    !
@@ -261,12 +269,14 @@ subroutine read_building_properties_file(filename, ios_out)
    !
    open(newunit=iunit, file=trim(filename), status='old', action='read', iostat=ios_bpr)
    if (ios_bpr /= 0) then
-      write(*,'(A,A)') 'ERROR: Could not open building properties file: ', trim(filename)
+      write(logstr,'(A,A)') 'ERROR: Could not open building properties file: ', trim(filename)
+      call write_log(logstr, 0)
       ios_out = ios_bpr
       return
    end if
    !
-   write(*,'(A,A)') 'Info    : Reading building properties file: ', trim(filename)
+   write(logstr,'(A,A)') 'Info    : Reading building properties file: ', trim(filename)
+   call write_log(logstr, 0)
    !
    do
       read(iunit, '(A)', iostat=ios_bpr) line
@@ -276,7 +286,8 @@ subroutine read_building_properties_file(filename, ios_out)
       !
       read(line, *, iostat=ios_bpr) name, det_vol, gut_spacing
       if (ios_bpr /= 0) then
-         write(*,'(A,A)') 'WARNING: Could not parse line: ', trim(line)
+         write(logstr,'(A,A)') 'WARNING: Could not parse line: ', trim(line)
+         call write_log(logstr, 0)
          cycle
       end if
       !
@@ -295,9 +306,10 @@ subroutine read_building_properties_file(filename, ios_out)
    !
    close(iunit)
    !
-   write(*,'(A,I0,A,I0,A)') 'Info    : Matched properties for ', nmatched, ' of ', nbuildings, ' buildings'
-   if (use_building_detention) write(*,'(A)') 'Info    : Building detention volumes enabled'
-   if (use_building_gutters) write(*,'(A)') 'Info    : Building gutter system enabled'
+   write(logstr,'(A,I0,A,I0,A)') 'Info    : Matched properties for ', nmatched, ' of ', nbuildings, ' buildings'
+   call write_log(logstr, 0)
+   if (use_building_detention) call write_log('Info    : Building detention volumes enabled', 0)
+   if (use_building_gutters) call write_log('Info    : Building gutter system enabled', 0)
    !
    end subroutine read_building_properties_file
    !
@@ -313,7 +325,7 @@ subroutine identify_building_cells()
    logical :: inside
    integer, allocatable :: temp_cells(:)
    !
-   write(*,'(A)') 'Info    : Identifying building cells ...'
+   call write_log('Info    : Identifying building cells ...', 0)
    !
    ncells = mmax * nmax
    !
@@ -358,7 +370,7 @@ subroutine identify_building_cells()
       !
    end do
    !
-   write(*,'(A)') 'Info    : Building cell identification complete'
+   call write_log('Info    : Building cell identification complete', 0)
    !
    end subroutine identify_building_cells
    !
@@ -417,7 +429,7 @@ subroutine identify_perimeter_cells()
    integer, allocatable :: temp_perimeter(:)
    integer :: ncells
    !
-   write(*,'(A)') 'Info    : Identifying perimeter cells ...'
+   call write_log('Info    : Identifying perimeter cells ...', 0)
    !
    ncells = mmax * nmax
    !
@@ -459,7 +471,7 @@ subroutine identify_perimeter_cells()
       !
    end do
    !
-   write(*,'(A)') 'Info    : Perimeter cell identification complete'
+   call write_log('Info    : Perimeter cell identification complete', 0)
    !
    end subroutine identify_perimeter_cells
    !
@@ -480,7 +492,7 @@ subroutine identify_gutter_cells()
    integer :: temp_ip
    logical :: first_cell
    !
-   write(*,'(A)') 'Info    : Identifying gutter outlet cells ...'
+   call write_log('Info    : Identifying gutter outlet cells ...', 0)
    !
    do ibld = 1, nbuildings
       !
@@ -588,7 +600,7 @@ subroutine identify_gutter_cells()
       !
    end do
    !
-   write(*,'(A)') 'Info    : Gutter cell identification complete'
+   call write_log('Info    : Gutter cell identification complete', 0)
    !
    end subroutine identify_gutter_cells
    !
@@ -604,7 +616,8 @@ subroutine create_building_thin_dams()
    real*4, allocatable :: xthd(:), ythd(:)
    integer, allocatable :: uv_indices(:), vertices(:)
    !
-   write(*,'(A,I0,A)') 'Info    : Processing ', nbuildings, ' buildings as thin dams'
+   write(logstr,'(A,I0,A)') 'Info    : Processing ', nbuildings, ' buildings as thin dams'
+   call write_log(logstr, 0)
    !
    do ibld = 1, nbuildings
       !
@@ -625,7 +638,8 @@ subroutine create_building_thin_dams()
       !
    end do
    !
-   write(*,'(A,I0,A)') 'Info    : Processed ', nbuildings, ' buildings as thin dams'
+   write(logstr,'(A,I0,A)') 'Info    : Processed ', nbuildings, ' buildings as thin dams'
+   call write_log(logstr, 0)
    !
    end subroutine create_building_thin_dams
    !

From ae23dab6e8e5cee9deaba29b248cc4a84d48e19f Mon Sep 17 00:00:00 2001
From: Leynse <Tim.Leijnse@deltares.nl>
Date: Wed, 15 Apr 2026 11:42:31 +0200
Subject: [PATCH 3/6] - Review: added finalize_buildings - Both zs and dh are
 already real*8; the intermediate cast was discarding 7 digits. - Added
 documentation

---
 CLAUDE.md                       | 97 +++++++++++++++++++++++++++++++++
 docs/parameters.rst             | 14 ++++-
 source/src/sfincs_buildings.f90 | 58 ++++++++++++++++++--
 source/src/sfincs_lib.f90       |  6 +-
 4 files changed, 165 insertions(+), 10 deletions(-)
 create mode 100644 CLAUDE.md

diff --git a/CLAUDE.md b/CLAUDE.md
new file mode 100644
index 000000000..d1824b3b2
--- /dev/null
+++ b/CLAUDE.md
@@ -0,0 +1,97 @@
+# SFINCS Fortran Model — Claude Agent Configuration
+
+## Project Overview
+
+SFINCS (Super-Fast INundation of CoastS) is a reduced-complexity hydrodynamic model written in Fortran 90/95.
+It solves the depth-averaged shallow water equations on regular grids, quadtree grids, and subgrid-based meshes
+for coastal/fluvial/pluvial flood simulations.
+
+## Repository Layout
+
+```
+source/
+  src/
+    sfincs.f90                  # Main program entry point
+    sfincs_bmi.f90              # Basic Model Interface (BMI) — public API
+    sfincs_domain.f90           # High-level orchestration (initialize/update/finalize)
+    sfincs_data.f90             # Global variables and arrays
+    sfincs_input.f90            # Input parameter parsing (sfincs.inp)
+    sfincs_read.f90             # File readers (grid, bathymetry, forcings)
+    sfincs_momentum.f90         # Momentum equations (u/v velocity)
+    sfincs_continuity.f90       # Continuity equation (water level h/zs)
+    sfincs_boundaries.f90       # Water level and wave boundary conditions
+    sfincs_discharges.f90       # Discharge source/sink points
+    sfincs_meteo.f90            # Wind, pressure, precipitation forcing
+    sfincs_output.f90           # NetCDF/binary output
+    sfincs_structures.f90       # Weirs, thin dams, drainage structures
+    sfincs_infiltration.f90     # Green-Ampt / curve number infiltration
+    sfincs_vegetation.f90       # Vegetation drag
+    sfincs_snapwave.f90         # SnapWave coupling interface
+    sfincs_advection_diffusion.f90
+    sfincs_nonhydrostatic.f90
+    sfincs_crosssections.f90
+    sfincs_obspoints.f90
+    sfincs_runup_gauges.f90
+    sfincs_wavemaker.f90
+    sfincs_wave_enhanced_roughness.f90
+    sfincs_spiderweb.f90
+    sfincs_openacc.f90          # GPU acceleration
+    sfincs_timestep_analysis.f90
+    sfincs_bathtub.f90
+    sfincs_date.f90
+    sfincs_log.f90
+    sfincs_error.f90
+    snapwave/                   # SnapWave wave propagation submodule (Fortran)
+  third_party_open/             # External libraries (netcdf-fortran, Delft3D utils, bicgstab)
+docs/                           # Sphinx documentation (.rst files)
+```
+
+## Key Conventions
+
+### Code Style
+- Fortran 90/95; modules with `use` statements (never `include` for data sharing)
+- `implicit none` in every subroutine and module procedure
+- All global state lives in `sfincs_data.f90` — no new global variables in other modules
+- Variable naming: grid coordinates `xcor`/`ycor`, bathymetry `zb`/`dep`, water level `zs`/`h`, velocities `u`/`v`
+- Mask arrays: `msk` (active cells), `kfu`/`kfv` (velocity point masks)
+- New input parameters must be added to `sfincs_input.f90` with a default value and read from `sfincs.inp`
+
+### Module Integration Pattern
+When adding a new physics process:
+1. Create `sfincs_<process>.f90` with `initialize_<process>()`, `update_<process>()`, `finalize_<process>()` subroutines
+2. Add `use sfincs_<process>` inside the relevant subroutine in `sfincs_domain.f90` (NOT at module level)
+3. Call `initialize_<process>()` from `initialize_domain()` in `sfincs_domain.f90`
+4. Call `update_<process>()` at the correct point in the time loop
+5. Wire input keywords through `sfincs_input.f90`
+6. **Never break the BMI interface** in `sfincs_bmi.f90` — this is the public API used by coupling frameworks
+
+### BMI Interface
+`sfincs_bmi.f90` is the stable public API. Functions: `initialize`, `update`, `update_until`, `finalize`, `get_var`, `set_var`, `get_time_step`, `get_current_time`, etc. Do not rename or remove existing BMI functions.
+
+### Build System
+CMake-based. Build targets: `sfincs` (executable), `sfincs_dll` (shared library for BMI coupling), `sfincs_lib` (static library).
+
+### NetCDF Output
+All gridded output uses NetCDF via `sfincs_output.f90`. Follow existing patterns for defining dimensions, variables, and attributes. Coordinate reference system metadata must be included.
+
+## Roles and Slash Commands
+
+### `/review-fortran` — Code Reviewer
+Review Fortran changes for:
+- Correctness of numerical schemes (finite differences, time stepping, stability)
+- Array bounds safety and allocation/deallocation hygiene
+- Consistency with `sfincs_data.f90` variable naming and global state
+- BMI interface integrity (no breaking changes to `sfincs_bmi.f90`)
+- Module `use` placement (use inside subroutines, not at module level — except `sfincs_data`)
+- Memory leaks (allocated arrays must be deallocated in `finalize_*` routines)
+- OpenACC GPU directives correctness (in `sfincs_openacc.f90`)
+- Regression risk to existing functionality
+
+### `/develop-fortran` — Code Developer
+When implementing new features:
+- Follow the module integration pattern above
+- Add input keywords with defaults to `sfincs_input.f90`
+- Write the new module as `sfincs_<feature>.f90`
+- Integrate via `sfincs_domain.f90`, not directly in `sfincs.f90`
+- Do not modify third_party_open/ libraries
+- Keep numerical schemes consistent with the existing explicit time-stepping approach
diff --git a/docs/parameters.rst b/docs/parameters.rst
index 5cbf274d5..e1a94af39 100644
--- a/docs/parameters.rst
+++ b/docs/parameters.rst
@@ -676,5 +676,15 @@ Structures
 	  :description:		Drainage pumps, culverts and check valves are both specified using the same format file, put with a different indication of the type (type=1 is drainage pump, type=2 is culvert and type=3 is check valve).
 	  :units:		coordinates: m in projected UTM zone, discharges in m^3/s.
 	  :required:		no
-	  :format:		asc	 
-	  
+	  :format:		asc
+	bldfile = sfincs.bld
+	  :description:		Building footprint polygons. Each building is described by a header line ``name npoints ncols`` followed by ``npoints`` rows of ``x y`` coordinates (m, projected UTM). The polygon must be closed (first and last point identical). Flow through grid edges intersected by a building outline is blocked (thin dam). Rainfall falling on cells inside a building footprint is captured and redistributed to perimeter cells (or gutter outlets if bprfile is provided). Setting bldfile activates the building drainage module.
+	  :units:		coordinates: m in projected UTM zone.
+	  :required:		no
+	  :format:		asc
+	bprfile = sfincs.bpr
+	  :description:		Optional building properties file. Each line contains ``name detention_volume gutter_spacing``, matching building names defined in bldfile. ``detention_volume`` [m³] is the on-roof storage capacity before runoff is released. ``gutter_spacing`` [m] sets the distance between gutter outlet cells along the building perimeter; set to 0 to distribute runoff uniformly over all perimeter cells.
+	  :units:		detention_volume in m³, gutter_spacing in m.
+	  :required:		no
+	  :format:		asc
+
diff --git a/source/src/sfincs_buildings.f90 b/source/src/sfincs_buildings.f90
index ab555e068..15b980b67 100644
--- a/source/src/sfincs_buildings.f90
+++ b/source/src/sfincs_buildings.f90
@@ -44,7 +44,6 @@ module sfincs_buildings
    real*8 :: total_runoff_from_buildings = 0.0d0
    real*8 :: total_detention_stored = 0.0d0
    !
-   logical :: has_building_properties = .false.
    logical :: use_building_detention = .false.
    logical :: use_building_gutters = .false.
    !
@@ -349,6 +348,10 @@ subroutine identify_building_cells()
                !
                is_building_cell(ip) = .true.
                building_id(ip) = ibld
+               !
+               ! Set kcs=0 to make cell inactive: suppresses continuity update
+               ! AND prevents any boundary or inflow assignment on building cells.
+               !
                kcs(ip) = 0
                !
                buildings(ibld)%ncells_inside = buildings(ibld)%ncells_inside + 1
@@ -409,9 +412,12 @@ subroutine point_in_polygon(xp, yp, xv, yv, nv, inside)
    inside = .false.
    j = nv
    do i = 1, nv
-      intersect = ((yv(i) > yp) .neqv. (yv(j) > yp)) .and. &
-                  (xp < (xv(j) - xv(i))*(yp - yv(i))/(yv(j) - yv(i)) + xv(i))
-      if (intersect) inside = .not. inside
+      ! Guard against horizontal edges (yv(j)==yv(i)) before dividing:
+      !
+      if ((yv(i) > yp) .neqv. (yv(j) > yp)) then
+         if (xp < (xv(j) - xv(i))*(yp - yv(i))/(yv(j) - yv(i)) + xv(i)) &
+            inside = .not. inside
+      end if
       j = i
    end do
    !
@@ -635,6 +641,8 @@ subroutine create_building_thin_dams()
       end do
       !
       deallocate(xthd, ythd)
+      if (allocated(uv_indices)) deallocate(uv_indices)
+      if (allocated(vertices))   deallocate(vertices)
       !
    end do
    !
@@ -646,6 +654,9 @@ end subroutine create_building_thin_dams
    !
    subroutine accumulate_building_rainfall(precip, dt)
    !
+   ! precip is rainfall rate in m/s; dt in seconds; dx/dy in m.
+   ! rain_volume = rate [m/s] * dt [s] * area [m^2] => volume in m^3.
+   !
    use sfincs_data, only: dx, dy, netprcp
    !
    implicit none
@@ -720,7 +731,7 @@ subroutine redistribute_building_water(zs, dt)
          dh = volume_per_cell / cell_area
          do ic = 1, buildings(ibld)%ncells_gutter
             ip = buildings(ibld)%cells_gutter(ic)
-            zs(ip) = zs(ip) + real(dh, 4)
+            zs(ip) = zs(ip) + dh
          end do
       else
          if (buildings(ibld)%ncells_perimeter == 0) cycle
@@ -728,7 +739,7 @@ subroutine redistribute_building_water(zs, dt)
          dh = volume_per_cell / cell_area
          do ic = 1, buildings(ibld)%ncells_perimeter
             ip = buildings(ibld)%cells_perimeter(ic)
-            zs(ip) = zs(ip) + real(dh, 4)
+            zs(ip) = zs(ip) + dh
          end do
       end if
       !
@@ -740,4 +751,39 @@ subroutine redistribute_building_water(zs, dt)
    !
    end subroutine redistribute_building_water
    !
+   !
+   subroutine finalize_buildings()
+   !
+   ! Deallocate all module-level and per-building allocatables.
+   ! Called from sfincs_finalize to avoid leaks in DLL/BMI use.
+   !
+   implicit none
+   !
+   integer :: ibld
+   !
+   if (allocated(buildings)) then
+      do ibld = 1, nbuildings
+         if (allocated(buildings(ibld)%x))              deallocate(buildings(ibld)%x)
+         if (allocated(buildings(ibld)%y))              deallocate(buildings(ibld)%y)
+         if (allocated(buildings(ibld)%cells_inside))   deallocate(buildings(ibld)%cells_inside)
+         if (allocated(buildings(ibld)%cells_perimeter))deallocate(buildings(ibld)%cells_perimeter)
+         if (allocated(buildings(ibld)%cells_gutter))   deallocate(buildings(ibld)%cells_gutter)
+      end do
+      deallocate(buildings)
+   end if
+   !
+   if (allocated(is_building_cell))  deallocate(is_building_cell)
+   if (allocated(is_perimeter_cell)) deallocate(is_perimeter_cell)
+   if (allocated(is_gutter_cell))    deallocate(is_gutter_cell)
+   if (allocated(building_id))       deallocate(building_id)
+   !
+   nbuildings = 0
+   use_building_detention = .false.
+   use_building_gutters   = .false.
+   total_rain_on_buildings    = 0.0d0
+   total_runoff_from_buildings = 0.0d0
+   total_detention_stored     = 0.0d0
+   !
+   end subroutine finalize_buildings
+   !
 end module sfincs_buildings
diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90
index b0b27859a..9455bb1ef 100644
--- a/source/src/sfincs_lib.f90
+++ b/source/src/sfincs_lib.f90
@@ -14,7 +14,7 @@ module sfincs_lib
    use sfincs_meteo
    use sfincs_infiltration
    use sfincs_data
-   use sfincs_buildings, only: setup_buildings, update_buildings
+   use sfincs_buildings, only: setup_buildings, update_buildings, finalize_buildings
    use sfincs_date
    use sfincs_output
    use sfincs_ncinput
@@ -811,7 +811,9 @@ function sfincs_finalize() result(ierr)
    !
    call write_log('----------- Closing off SFINCS -----------', 1)
    !
-   ! call finalize_parameters()
+   call finalize_parameters()
+   !
+   if (has_buildings) call finalize_buildings()   
    !
    ierr = 0
    !

From ee987cb6f605205b8df1f22235c9a684cc85a3c4 Mon Sep 17 00:00:00 2001
From: Leynse <Tim.Leijnse@deltares.nl>
Date: Wed, 15 Apr 2026 11:48:35 +0200
Subject: [PATCH 4/6] - Add a building_mask to netcdf output

---
 source/src/sfincs_ncoutput.F90 | 86 +++++++++++++++++++++++++++++-----
 1 file changed, 73 insertions(+), 13 deletions(-)

diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90
index 33f0ba8e3..ff4cb0e8f 100644
--- a/source/src/sfincs_ncoutput.F90
+++ b/source/src/sfincs_ncoutput.F90
@@ -26,6 +26,7 @@ module sfincs_ncoutput
       integer :: manning_varid
       integer :: pnonh_varid
       integer :: subgridslope_varid
+      integer :: building_msk_varid
       !
       integer :: mesh2d_varid
       integer :: mesh2d_node_x_varid, mesh2d_node_y_varid
@@ -74,7 +75,8 @@ subroutine ncoutput_regular_map_init()
    !
    use sfincs_date
    use sfincs_data
-   use sfincs_snapwave   
+   use sfincs_snapwave
+   use sfincs_buildings, only: is_building_cell
    !
    implicit none   
    !   
@@ -185,15 +187,26 @@ subroutine ncoutput_regular_map_init()
    NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'face_coordinates', 'x y'))
    NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'corner_coordinates', 'corner_x corner_y'))   
    !
-   NF90(nf90_def_var(map_file%ncid, 'msk', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%msk_varid)) ! input msk value in cell centre   
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, '_FillValue', FILL_VALUE))      
+   NF90(nf90_def_var(map_file%ncid, 'msk', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%msk_varid)) ! input msk value in cell centre
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, '_FillValue', FILL_VALUE))
    NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'units', '-'))
    NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'standard_name', 'land_binary_mask')) ! land_binary_mask but with added boundary=2
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'long_name', 'msk_active_cells')) 
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3'))    
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'coordinates', 'x y'))   
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'long_name', 'msk_active_cells'))
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3'))
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'coordinates', 'x y'))
    NF90(nf90_def_var_deflate(map_file%ncid, map_file%msk_varid, 1, 1, nc_deflate_level)) ! deflate
    !
+   if (has_buildings) then
+      NF90(nf90_def_var(map_file%ncid, 'building_msk', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%building_msk_varid))
+      NF90(nf90_def_var_deflate(map_file%ncid, map_file%building_msk_varid, 1, 1, nc_deflate_level))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, '_FillValue', FILL_VALUE))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'units', '-'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'standard_name', 'building_binary_mask'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'long_name', 'building_footprint_mask'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'description', 'outside_building=0, inside_building=1'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'coordinates', 'x y'))
+   end if
+   !
    ! Infiltration map
    !
    if (infiltration) then
@@ -783,7 +796,26 @@ subroutine ncoutput_regular_map_init()
       !
    enddo
    !
-   NF90(nf90_put_var(map_file%ncid, map_file%msk_varid, zsg, (/1, 1/))) ! write msk     
+   NF90(nf90_put_var(map_file%ncid, map_file%msk_varid, zsg, (/1, 1/))) ! write msk
+   !
+   ! Write building mask
+   !
+   if (has_buildings) then
+      !
+      zsg = 0
+      !
+      do nm = 1, np
+         !
+         n    = z_index_z_n(nm)
+         m    = z_index_z_m(nm)
+         !
+         if (is_building_cell(nm)) zsg(m, n) = 1.0
+         !
+      enddo
+      !
+      NF90(nf90_put_var(map_file%ncid, map_file%building_msk_varid, zsg, (/1, 1/))) ! write building mask
+      !
+   endif
    !
    ! Write SnapWave msk
    !
@@ -863,9 +895,10 @@ subroutine ncoutput_quadtree_map_init()
    ! 2. write grid/msk/zb to file
    !
    use sfincs_date
-   use sfincs_data   
-   use sfincs_snapwave   
+   use sfincs_data
+   use sfincs_snapwave
    use quadtree
+   use sfincs_buildings, only: is_building_cell
    !
    implicit none   
    !   
@@ -1060,11 +1093,21 @@ subroutine ncoutput_quadtree_map_init()
    !
    NF90(nf90_def_var(map_file%ncid, 'msk', NF90_INT, (/map_file%nmesh2d_face_dimid/), map_file%msk_varid)) ! input msk value in cell centre
    NF90(nf90_def_var_deflate(map_file%ncid, map_file%msk_varid, 1, 1, nc_deflate_level))
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, '_FillValue', -999))      
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, '_FillValue', -999))
    NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'units', '-'))
    NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'standard_name', 'mask'))
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'long_name', 'msk_active_cells')) 
-   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3, wavemaker=4'))    
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'long_name', 'msk_active_cells'))
+   NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3, wavemaker=4'))
+   !
+   if (has_buildings) then
+      NF90(nf90_def_var(map_file%ncid, 'building_msk', NF90_INT, (/map_file%nmesh2d_face_dimid/), map_file%building_msk_varid))
+      NF90(nf90_def_var_deflate(map_file%ncid, map_file%building_msk_varid, 1, 1, nc_deflate_level))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, '_FillValue', -999))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'units', '-'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'standard_name', 'building_binary_mask'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'long_name', 'building_footprint_mask'))
+      NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'description', 'outside_building=0, inside_building=1'))
+   end if
    !
    ! Time variables   
    !
@@ -1532,7 +1575,24 @@ subroutine ncoutput_quadtree_map_init()
       endif
    enddo 
    !
-   NF90(nf90_put_var(map_file%ncid, map_file%msk_varid, vtmpi)) ! write msk 
+   NF90(nf90_put_var(map_file%ncid, map_file%msk_varid, vtmpi)) ! write msk
+   !
+   ! Write building mask
+   !
+   if (has_buildings) then
+      !
+      vtmpi = 0
+      !
+      do nmq = 1, quadtree_nr_points
+         nm = index_sfincs_in_quadtree(nmq)
+         if (nm > 0) then
+            if (is_building_cell(nm)) vtmpi(nmq) = 1
+         endif
+      enddo
+      !
+      NF90(nf90_put_var(map_file%ncid, map_file%building_msk_varid, vtmpi)) ! write building mask
+      !
+   endif
    !
    ! Write SnapWave msk
    !

From 1f47f50edb1500528665a699d3fa03c37e75b3dc Mon Sep 17 00:00:00 2001
From: Leynse <Tim.Leijnse@deltares.nl>
Date: Wed, 15 Apr 2026 11:59:45 +0200
Subject: [PATCH 5/6] - Example of adding current_detention in meters to netcdf
 output - To discuss, is meters or m3 more useful?

---
 source/src/sfincs_buildings.f90 | 26 ++++++++--
 source/src/sfincs_ncoutput.F90  | 84 ++++++++++++++++++++++++++++-----
 2 files changed, 94 insertions(+), 16 deletions(-)

diff --git a/source/src/sfincs_buildings.f90 b/source/src/sfincs_buildings.f90
index 15b980b67..a90269833 100644
--- a/source/src/sfincs_buildings.f90
+++ b/source/src/sfincs_buildings.f90
@@ -39,6 +39,7 @@ module sfincs_buildings
    logical, allocatable :: is_perimeter_cell(:)
    logical, allocatable :: is_gutter_cell(:)
    integer, allocatable :: building_id(:)
+   real*4,  allocatable :: detention_level_cell(:)  ! current detention as equiv. water depth [m] per cell
    !
    real*8 :: total_rain_on_buildings = 0.0d0
    real*8 :: total_runoff_from_buildings = 0.0d0
@@ -127,12 +128,14 @@ subroutine initialize_building_masks()
       allocate(is_perimeter_cell(ncells))
       allocate(is_gutter_cell(ncells))
       allocate(building_id(ncells))
+      allocate(detention_level_cell(ncells))
    end if
    !
    is_building_cell = .false.
    is_perimeter_cell = .false.
    is_gutter_cell = .false.
    building_id = 0
+   detention_level_cell = 0.0
    !
    end subroutine initialize_building_masks
    !
@@ -749,6 +752,20 @@ subroutine redistribute_building_water(zs, dt)
       !
    end do
    !
+   ! Update per-cell detention equivalent depth [m] for output
+   !
+   if (use_building_detention) then
+      do ibld = 1, nbuildings
+         if (buildings(ibld)%total_area > 0.0d0) then
+            do ic = 1, buildings(ibld)%ncells_inside
+               ip = buildings(ibld)%cells_inside(ic)
+               detention_level_cell(ip) = real(buildings(ibld)%current_detention / &
+                                               buildings(ibld)%total_area, 4)
+            end do
+         end if
+      end do
+   end if
+   !
    end subroutine redistribute_building_water
    !
    !
@@ -772,10 +789,11 @@ subroutine finalize_buildings()
       deallocate(buildings)
    end if
    !
-   if (allocated(is_building_cell))  deallocate(is_building_cell)
-   if (allocated(is_perimeter_cell)) deallocate(is_perimeter_cell)
-   if (allocated(is_gutter_cell))    deallocate(is_gutter_cell)
-   if (allocated(building_id))       deallocate(building_id)
+   if (allocated(is_building_cell))      deallocate(is_building_cell)
+   if (allocated(is_perimeter_cell))     deallocate(is_perimeter_cell)
+   if (allocated(is_gutter_cell))        deallocate(is_gutter_cell)
+   if (allocated(building_id))           deallocate(building_id)
+   if (allocated(detention_level_cell))  deallocate(detention_level_cell)
    !
    nbuildings = 0
    use_building_detention = .false.
diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90
index ff4cb0e8f..e4699a9b1 100644
--- a/source/src/sfincs_ncoutput.F90
+++ b/source/src/sfincs_ncoutput.F90
@@ -27,6 +27,7 @@ module sfincs_ncoutput
       integer :: pnonh_varid
       integer :: subgridslope_varid
       integer :: building_msk_varid
+      integer :: detention_level_varid
       !
       integer :: mesh2d_varid
       integer :: mesh2d_node_x_varid, mesh2d_node_y_varid
@@ -76,10 +77,10 @@ subroutine ncoutput_regular_map_init()
    use sfincs_date
    use sfincs_data
    use sfincs_snapwave
-   use sfincs_buildings, only: is_building_cell
+   use sfincs_buildings, only: is_building_cell, use_building_detention
+   !
+   implicit none
    !
-   implicit none   
-   !   
    integer                      :: nm, n, m, ntmx
    !
    real*4, dimension(:,:), allocatable :: zsg
@@ -207,6 +208,17 @@ subroutine ncoutput_regular_map_init()
       NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'coordinates', 'x y'))
    end if
    !
+   if (has_buildings .and. use_building_detention) then
+      NF90(nf90_def_var(map_file%ncid, 'detention_level', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%detention_level_varid))
+      NF90(nf90_def_var_deflate(map_file%ncid, map_file%detention_level_varid, 1, 1, nc_deflate_level))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, '_FillValue', FILL_VALUE))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'units', 'm'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'standard_name', 'building_detention_equivalent_depth'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'long_name', 'building_detention_storage_equivalent_depth'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'description', 'current_detention_volume_divided_by_building_footprint_area'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'coordinates', 'x y'))
+   end if
+   !
    ! Infiltration map
    !
    if (infiltration) then
@@ -898,10 +910,10 @@ subroutine ncoutput_quadtree_map_init()
    use sfincs_data
    use sfincs_snapwave
    use quadtree
-   use sfincs_buildings, only: is_building_cell
+   use sfincs_buildings, only: is_building_cell, use_building_detention
+   !
+   implicit none
    !
-   implicit none   
-   !   
    integer    :: nm, nmq, nmu1, num1, n, m, nn, ntmx, n_nodes, n_faces, iref
    real*4     :: dxx, dyy
    !
@@ -1109,7 +1121,17 @@ subroutine ncoutput_quadtree_map_init()
       NF90(nf90_put_att(map_file%ncid, map_file%building_msk_varid, 'description', 'outside_building=0, inside_building=1'))
    end if
    !
-   ! Time variables   
+   if (has_buildings .and. use_building_detention) then
+      NF90(nf90_def_var(map_file%ncid, 'detention_level', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%detention_level_varid))
+      NF90(nf90_def_var_deflate(map_file%ncid, map_file%detention_level_varid, 1, 1, nc_deflate_level))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, '_FillValue', FILL_VALUE))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'units', 'm'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'standard_name', 'building_detention_equivalent_depth'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'long_name', 'building_detention_storage_equivalent_depth'))
+      NF90(nf90_put_att(map_file%ncid, map_file%detention_level_varid, 'description', 'current_detention_volume_divided_by_building_footprint_area'))
+   end if
+   !
+   ! Time variables
    !
    trefstr_iso8601 = date_to_iso8601(trefstr)
    !
@@ -2216,11 +2238,12 @@ subroutine ncoutput_his_init()
    !  
    subroutine ncoutput_update_regular_map(t,ntmapout)
    !
-   ! Write time, zs, u, v  
+   ! Write time, zs, u, v
    !
    use sfincs_data
    use sfincs_nonhydrostatic
    use sfincs_snapwave
+   use sfincs_buildings, only: detention_level_cell, use_building_detention
    !
    implicit none   
    !
@@ -2649,7 +2672,26 @@ subroutine ncoutput_update_regular_map(t,ntmapout)
       NF90(nf90_put_var(map_file%ncid, map_file%pnonh_varid, zsg, (/1, 1, ntmapout/))) ! write h
       !
    endif
-   !           
+   !
+   ! Write building detention equivalent depth (time-varying)
+   !
+   if (has_buildings .and. use_building_detention) then
+      !
+      zsg = FILL_VALUE
+      !
+      do nm = 1, np
+         !
+         n    = z_index_z_n(nm)
+         m    = z_index_z_m(nm)
+         !
+         zsg(m, n) = detention_level_cell(nm)
+         !
+      enddo
+      !
+      NF90(nf90_put_var(map_file%ncid, map_file%detention_level_varid, zsg, (/1, 1, ntmapout/))) ! write detention_level
+      !
+   endif
+   !
    NF90(nf90_sync(map_file%ncid)) !write away intermediate data ! TL: in first test it seems to be faster to let the file update than keep in memory
    !
    end subroutine
@@ -2657,12 +2699,13 @@ subroutine ncoutput_update_regular_map(t,ntmapout)
    
    subroutine ncoutput_update_quadtree_map(t,ntmapout)
       !
-      ! Write time, zs, u, v  
+      ! Write time, zs, u, v
       !
-      use sfincs_data   
+      use sfincs_data
       use sfincs_snapwave
       use sfincs_nonhydrostatic
       use quadtree
+      use sfincs_buildings, only: detention_level_cell, use_building_detention
       ! use snapwave_data
       !
       implicit none   
@@ -3050,8 +3093,25 @@ subroutine ncoutput_update_quadtree_map(t,ntmapout)
          !
       endif
       !
+      ! Write building detention equivalent depth (time-varying)
+      !
+      if (has_buildings .and. use_building_detention) then
+         !
+         vtmp = FILL_VALUE
+         !
+         do nmq = 1, quadtree_nr_points
+            nm = index_sfincs_in_quadtree(nmq)
+            if (nm > 0) then
+               vtmp(nmq) = detention_level_cell(nm)
+            endif
+         enddo
+         !
+         NF90(nf90_put_var(map_file%ncid, map_file%detention_level_varid, vtmp, (/1, ntmapout/))) ! write detention_level
+         !
+      endif
+      !
       NF90(nf90_sync(map_file%ncid)) !write away intermediate data ! TL: in first test it seems to be faster to let the file update than keep in memory
-      !      
+      !
    end subroutine
    
    

From fe5d50bb00ed64621c60f7cc31e5b2407ab57562 Mon Sep 17 00:00:00 2001
From: Tim Leijnse <timleynse94@gmail.com>
Date: Wed, 8 Jul 2026 13:24:14 +0200
Subject: [PATCH 6/6] Delete CLAUDE.md

---
 CLAUDE.md | 97 -------------------------------------------------------
 1 file changed, 97 deletions(-)
 delete mode 100644 CLAUDE.md

diff --git a/CLAUDE.md b/CLAUDE.md
deleted file mode 100644
index d1824b3b2..000000000
--- a/CLAUDE.md
+++ /dev/null
@@ -1,97 +0,0 @@
-# SFINCS Fortran Model — Claude Agent Configuration
-
-## Project Overview
-
-SFINCS (Super-Fast INundation of CoastS) is a reduced-complexity hydrodynamic model written in Fortran 90/95.
-It solves the depth-averaged shallow water equations on regular grids, quadtree grids, and subgrid-based meshes
-for coastal/fluvial/pluvial flood simulations.
-
-## Repository Layout
-
-```
-source/
-  src/
-    sfincs.f90                  # Main program entry point
-    sfincs_bmi.f90              # Basic Model Interface (BMI) — public API
-    sfincs_domain.f90           # High-level orchestration (initialize/update/finalize)
-    sfincs_data.f90             # Global variables and arrays
-    sfincs_input.f90            # Input parameter parsing (sfincs.inp)
-    sfincs_read.f90             # File readers (grid, bathymetry, forcings)
-    sfincs_momentum.f90         # Momentum equations (u/v velocity)
-    sfincs_continuity.f90       # Continuity equation (water level h/zs)
-    sfincs_boundaries.f90       # Water level and wave boundary conditions
-    sfincs_discharges.f90       # Discharge source/sink points
-    sfincs_meteo.f90            # Wind, pressure, precipitation forcing
-    sfincs_output.f90           # NetCDF/binary output
-    sfincs_structures.f90       # Weirs, thin dams, drainage structures
-    sfincs_infiltration.f90     # Green-Ampt / curve number infiltration
-    sfincs_vegetation.f90       # Vegetation drag
-    sfincs_snapwave.f90         # SnapWave coupling interface
-    sfincs_advection_diffusion.f90
-    sfincs_nonhydrostatic.f90
-    sfincs_crosssections.f90
-    sfincs_obspoints.f90
-    sfincs_runup_gauges.f90
-    sfincs_wavemaker.f90
-    sfincs_wave_enhanced_roughness.f90
-    sfincs_spiderweb.f90
-    sfincs_openacc.f90          # GPU acceleration
-    sfincs_timestep_analysis.f90
-    sfincs_bathtub.f90
-    sfincs_date.f90
-    sfincs_log.f90
-    sfincs_error.f90
-    snapwave/                   # SnapWave wave propagation submodule (Fortran)
-  third_party_open/             # External libraries (netcdf-fortran, Delft3D utils, bicgstab)
-docs/                           # Sphinx documentation (.rst files)
-```
-
-## Key Conventions
-
-### Code Style
-- Fortran 90/95; modules with `use` statements (never `include` for data sharing)
-- `implicit none` in every subroutine and module procedure
-- All global state lives in `sfincs_data.f90` — no new global variables in other modules
-- Variable naming: grid coordinates `xcor`/`ycor`, bathymetry `zb`/`dep`, water level `zs`/`h`, velocities `u`/`v`
-- Mask arrays: `msk` (active cells), `kfu`/`kfv` (velocity point masks)
-- New input parameters must be added to `sfincs_input.f90` with a default value and read from `sfincs.inp`
-
-### Module Integration Pattern
-When adding a new physics process:
-1. Create `sfincs_<process>.f90` with `initialize_<process>()`, `update_<process>()`, `finalize_<process>()` subroutines
-2. Add `use sfincs_<process>` inside the relevant subroutine in `sfincs_domain.f90` (NOT at module level)
-3. Call `initialize_<process>()` from `initialize_domain()` in `sfincs_domain.f90`
-4. Call `update_<process>()` at the correct point in the time loop
-5. Wire input keywords through `sfincs_input.f90`
-6. **Never break the BMI interface** in `sfincs_bmi.f90` — this is the public API used by coupling frameworks
-
-### BMI Interface
-`sfincs_bmi.f90` is the stable public API. Functions: `initialize`, `update`, `update_until`, `finalize`, `get_var`, `set_var`, `get_time_step`, `get_current_time`, etc. Do not rename or remove existing BMI functions.
-
-### Build System
-CMake-based. Build targets: `sfincs` (executable), `sfincs_dll` (shared library for BMI coupling), `sfincs_lib` (static library).
-
-### NetCDF Output
-All gridded output uses NetCDF via `sfincs_output.f90`. Follow existing patterns for defining dimensions, variables, and attributes. Coordinate reference system metadata must be included.
-
-## Roles and Slash Commands
-
-### `/review-fortran` — Code Reviewer
-Review Fortran changes for:
-- Correctness of numerical schemes (finite differences, time stepping, stability)
-- Array bounds safety and allocation/deallocation hygiene
-- Consistency with `sfincs_data.f90` variable naming and global state
-- BMI interface integrity (no breaking changes to `sfincs_bmi.f90`)
-- Module `use` placement (use inside subroutines, not at module level — except `sfincs_data`)
-- Memory leaks (allocated arrays must be deallocated in `finalize_*` routines)
-- OpenACC GPU directives correctness (in `sfincs_openacc.f90`)
-- Regression risk to existing functionality
-
-### `/develop-fortran` — Code Developer
-When implementing new features:
-- Follow the module integration pattern above
-- Add input keywords with defaults to `sfincs_input.f90`
-- Write the new module as `sfincs_<feature>.f90`
-- Integrate via `sfincs_domain.f90`, not directly in `sfincs.f90`
-- Do not modify third_party_open/ libraries
-- Keep numerical schemes consistent with the existing explicit time-stepping approach