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/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..a90269833
--- /dev/null
+++ b/source/src/sfincs_buildings.f90
@@ -0,0 +1,807 @@
+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)
+   !
+   use sfincs_log
+   !
+   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*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
+   real*8 :: total_detention_stored = 0.0d0
+   !
+   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
+   !
+   call write_log('Info    : Setting up buildings ...', 1)
+   !
+   call initialize_building_masks()
+   !
+   call read_building_file(bldfile, ios)
+   !
+   if (ios /= 0) then
+      call write_log('ERROR: Failed to read building file', 0)
+      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
+         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
+   !
+   call identify_building_cells()
+   call identify_perimeter_cells()
+   !
+   if (use_building_gutters) then
+      call identify_gutter_cells()
+   end if
+   !
+   call create_building_thin_dams()
+   !
+   call write_log('Info    : Building setup complete', 0)
+   !
+   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))
+      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
+   !
+   !
+   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(logstr,'(A,A)') 'ERROR: Could not open building file: ', trim(filename)
+      call write_log(logstr, 0)
+      ios_out = ios_bld
+      return
+   end if
+   !
+   write(logstr,'(A,A)') 'Info    : Reading building file: ', trim(filename)
+   call write_log(logstr, 0)
+   !
+   ! 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
+      call write_log('WARNING: No buildings found in file', 0)
+      close(iunit)
+      return
+   end if
+   !
+   write(logstr,'(A,I0,A)') 'Info    : Found ', nbuildings, ' buildings'
+   call write_log(logstr, 0)
+   !
+   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(logstr,'(A,A)') 'ERROR: Reading coordinates for building ', trim(name)
+            call write_log(logstr, 0)
+            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(logstr,'(A,A,A)') 'WARNING: Building ', trim(name), ' polygon not closed'
+         call write_log(logstr, 0)
+      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(logstr,'(A,I0,A)') 'Info    : Successfully read ', nbuildings, ' buildings'
+   call write_log(logstr, 0)
+   !
+   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(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(logstr,'(A,A)') 'Info    : Reading building properties file: ', trim(filename)
+   call write_log(logstr, 0)
+   !
+   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(logstr,'(A,A)') 'WARNING: Could not parse line: ', trim(line)
+         call write_log(logstr, 0)
+         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(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
+   !
+   !
+   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(:)
+   !
+   call write_log('Info    : Identifying building cells ...', 0)
+   !
+   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
+               !
+               ! 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
+               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
+   !
+   call write_log('Info    : Building cell identification complete', 0)
+   !
+   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
+      ! 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
+   !
+   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
+   !
+   call write_log('Info    : Identifying perimeter cells ...', 0)
+   !
+   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
+   !
+   call write_log('Info    : Perimeter cell identification complete', 0)
+   !
+   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
+   !
+   call write_log('Info    : Identifying gutter outlet cells ...', 0)
+   !
+   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
+   !
+   call write_log('Info    : Gutter cell identification complete', 0)
+   !
+   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(logstr,'(A,I0,A)') 'Info    : Processing ', nbuildings, ' buildings as thin dams'
+   call write_log(logstr, 0)
+   !
+   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)
+      if (allocated(uv_indices)) deallocate(uv_indices)
+      if (allocated(vertices))   deallocate(vertices)
+      !
+   end do
+   !
+   write(logstr,'(A,I0,A)') 'Info    : Processed ', nbuildings, ' buildings as thin dams'
+   call write_log(logstr, 0)
+   !
+   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
+   !
+   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) + dh
+         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) + dh
+         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
+   !
+   ! 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
+   !
+   !
+   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)
+   if (allocated(detention_level_cell))  deallocate(detention_level_cell)
+   !
+   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_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..9455bb1ef 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, finalize_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
@@ -800,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
    !
diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90
index 33f0ba8e3..e4699a9b1 100644
--- a/source/src/sfincs_ncoutput.F90
+++ b/source/src/sfincs_ncoutput.F90
@@ -26,6 +26,8 @@ module sfincs_ncoutput
       integer :: manning_varid
       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
@@ -74,10 +76,11 @@ subroutine ncoutput_regular_map_init()
    !
    use sfincs_date
    use sfincs_data
-   use sfincs_snapwave   
+   use sfincs_snapwave
+   use sfincs_buildings, only: is_building_cell, use_building_detention
+   !
+   implicit none
    !
-   implicit none   
-   !   
    integer                      :: nm, n, m, ntmx
    !
    real*4, dimension(:,:), allocatable :: zsg
@@ -185,15 +188,37 @@ 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
+   !
+   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
@@ -783,7 +808,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,12 +907,13 @@ 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, 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
    !
@@ -1060,13 +1105,33 @@ 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'))    
-   !
-   ! Time variables   
+   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
+   !
+   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)
    !
@@ -1532,7 +1597,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
    !
@@ -2156,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   
    !
@@ -2589,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
@@ -2597,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   
@@ -2990,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