Repository: Deltares/Wflow.jl (julia)

Issue:
Function `save_state` should close state file
### Wflow version checks

- [X] I have checked that this issue has not already been reported.

- [X] I have checked that this bug exists on the latest version of Wflow.


### Reproducible Example

Run function `Wflow.save_state(model)`.

### Current behaviour

State file is not closed.

### Desired behaviour

Close state file.

### Additional Context

_No response_

Pull request:
## Issue addressed
Fixes #369 

## Explanation
See issue https://github.com/Deltares/Wflow.jl/issues/369.

## Checklist
- [x] Updated tests or added new tests
- [x] Branch is up to date with `master`
- [x] Tests & pre-commit hooks pass
- [x] Updated documentation if needed
- [x] Updated changelog.md if needed

Patch:
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index cd78cf337..f79c1a8eb 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -21,6 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   routing, the floodplain discharge was not added to the inflow of these locations, 2) the
   `to_river` variable from overland flow and lateral subsurface flow was not added to the
   inflow of these locations.
+- Close netCDF `NCDataset` with state variables in extended BMI function `save_state`.
 
 ### Changed
 - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was
diff --git a/src/bmi.jl b/src/bmi.jl
index 1c7d8f476..17d385ea9 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -76,7 +76,10 @@ end
 "Write state output to netCDF and close files."
 function BMI.finalize(model::Model)
     @unpack config, writer, clock = model
-    write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)
+    # it is possible that the state dataset has been closed by `save_state`  
+    if !isnothing(writer.state_dataset) && isopen(writer.state_dataset)
+        write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)
+    end
     reset_clock!(model.clock, config)
     close_files(model, delete_output = false)
 end
@@ -394,6 +397,7 @@ function save_state(model::Model)
         @info "Write output states to NetCDF file `$(model.writer.state_nc_path)`."
     end
     write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)
+    close(writer.state_dataset)
 end
 
 function get_start_unix_time(model::Model)
diff --git a/src/flow.jl b/src/flow.jl
index a07f15991..a59475a3f 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -524,8 +524,8 @@ end
     bankfull_depth::Vector{T} | "m"                         # bankfull depth
     zb::Vector{T} | "m"                                     # river bed elevation
     froude_limit::Bool | "-" | 0 | "none" | "none"          # if true a check is performed if froude number > 1.0 (algorithm is modified)
-    reservoir_index::Vector{Int} | "-"                      # river cell index with a reservoir (each index of reservoir_index maps to reservoir i in reservoir field)
-    lake_index::Vector{Int} | "-"                           # river cell index with a lake (each index of lake_index maps to lake i in lake field)
+    reservoir_index::Vector{Int} | "-" | 0                  # river cell index with a reservoir (each index of reservoir_index maps to reservoir i in reservoir field)
+    lake_index::Vector{Int} | "-" | 0                       # river cell index with a lake (each index of lake_index maps to lake i in lake field)
     waterbody::Vector{Bool} | "-"                           # water body cells (reservoir or lake)
     reservoir::R | "-" | 0                                  # Reservoir model struct of arrays
     lake::L | "-" | 0                                       # Lake model struct of arrays


Test patch:
diff --git a/test/bmi.jl b/test/bmi.jl
index 772ee352a..3909dfe67 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -201,4 +201,5 @@ end
         :info,
         "Write output states to NetCDF file `$(model.writer.state_nc_path)`.",
     ) Wflow.save_state(model)
+    @test !isopen(model.writer.state_dataset)
 end