Merge branch 'main' into GreenAmpt-infiltration

Deltares · Oct 11, 2023 · bcd3cdc · bcd3cdc
2 parents 594218b + 90684f8
commit bcd3cdc
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 5 deletions.
diff --git a/docs/input.rst b/docs/input.rst
@@ -326,6 +326,7 @@ SFINCS allows the specification of the following options for accounting for infi
 Infiltration is specified with either constant in time values in mm/hr (both uniform and spatially varying), or using more detailed parameters for the Curve Number method, The Green-Ampt method or Horton method.
 
 **NOTE - Infiltration in SFINCS is only turned on when any rainfall is forced'** 
+
 **NOTE - Infiltration methods in SFINCS are not designed to be stacked**
 
 
@@ -338,8 +339,9 @@ Specify the keyword:
 
 	qinf = 1.0
 	
-**NOTE - To have some control that no infiltration is added an areas like the sea for this spatially uniform constant in time infiltration method, only infiltration is added to cells above a certain elevation above the bed level reference height **
-**NOTE - By default this is set to 0, qinf_zmin = 0 (default), so below e.g. MSL, no infiltration is added
+**NOTE - To have some control that no infiltration is added an areas like the sea for this spatially uniform constant in time infiltration method, only infiltration is added to cells above a certain elevation above the bed level reference height**
+
+**NOTE - By default this is set to 0, qinf_zmin = 0 (default), so below e.g. MSL, no infiltration is added**
 
 Spatially varying constant in time:
 %%%%%
@@ -383,6 +385,7 @@ where Smax = the soil's maximum moisture storage capacity. Smax typically derive
 
 
 **Without recovery**
+
 For spatially varying infiltration values per cell using the Curve Number method without recovery use the scsfile option, with the same grid based input as the depfile using a binary file. Note here that in pre-processing the wanted CN values should be converted to S values following:
 * scsfile: maximum soil moisture storage capacity in inches
 
@@ -408,7 +411,8 @@ Within SFINCS, the Curve number method with recovery can be used as follows. The
 * ksfile: saturated hydraulic conductivity in mm/hr
 
 Using the saturated hydraulic conductivity, a recovery variables recovery constant (kr) and minimum recovery time before a new rainfall event occurs (hours) are computed similar to SWMM. 
-The Curve Number has been implemented in hydromt-SFINCS. See https://deltares.github.io/hydromt_sfincs/latest/_generated/hydromt_sfincs.SfincsModel.setup_cn_infiltration_with_kr.html for more information.
+The Curve Number has been implemented in hydromt-SFINCS. 
+See https://deltares.github.io/hydromt_sfincs/latest/_generated/hydromt_sfincs.SfincsModel.setup_cn_infiltration_with_kr.html for more information.
 
 This option does support restart functionality. 
 

diff --git a/source/src/sfincs_meteo.f90 b/source/src/sfincs_meteo.f90
@@ -247,6 +247,7 @@ subroutine update_spiderweb_data()
    real*4 meteo_t
    real*4 twfac
    real*4 merge_frac
+   real*4 tspinup_fac
    real*4 wa
    real*4 z0marine
    real*4 z0l
@@ -264,6 +265,7 @@ subroutine update_spiderweb_data()
          meteo_t = meteo_t1
       endif
       !
+      itw0 = 1
       do itspw = 1, spw_nt
          if (spw_times(itspw)<=meteo_t) then
             itw0 = itspw
@@ -275,6 +277,33 @@ subroutine update_spiderweb_data()
       meteo_t = min(meteo_t, spw_times(itw1))
       !
       twfac  = (meteo_t - spw_times(itw0))/max(spw_times(itw1) - spw_times(itw0), 1.0e-6)
+      tspinup_fac = 1.0 
+      !
+      ! Check if spw data is not yet available
+      ! If so, use first time in spw 
+      !
+      if (meteo_t < spw_times(1)) then
+         !
+         itw0 = 1
+         itw1 = 1
+         twfac = 1.0         
+         !
+         tspinup_fac = max(1.0 - (spw_times(1) - meteo_t) / 21600.0, 0.0)
+         !
+      endif
+      !
+      ! Check if spw data is no longer available
+      ! If so, use last time in spw 
+      !
+      if (meteo_t > spw_times(spw_nt)) then
+         !
+         itw0 = spw_nt
+         itw1 = spw_nt
+         twfac = 1.0         
+         !
+         tspinup_fac = max(1.0 - (meteo_t - spw_times(spw_nt)) / 21600.0, 0.0)
+         !
+      endif
       !
       ! Eye
       ! 
@@ -453,6 +482,8 @@ subroutine update_spiderweb_data()
          else
             merge_frac = 1.0
          endif
+         ! 
+         merge_frac = merge_frac * tspinup_fac
          !
          if (wind_reduction) then
             !

diff --git a/source/src/sfincs_wavemaker.f90 b/source/src/sfincs_wavemaker.f90
@@ -1137,7 +1137,7 @@ subroutine read_wavemaker_polylines()
       !
       do ip = 1, npuv
          !
-         if (kcuv(uv_index_u_nmd(ip)) == 4)  uv_flags_adv(ip) = 0
+         if (kcuv(uv_index_u_nmd(ip)) == 4)  uv_flags_adv(ip) = 0 !QUESTION - TL: goes wrong if you don't specify a quadtree grid !? > should make this clear in documentation
          if (kcuv(uv_index_u_nmu(ip)) == 4)  uv_flags_adv(ip) = 0
          if (kcuv(uv_index_u_num(ip)) == 4)  uv_flags_adv(ip) = 0
          if (kcuv(uv_index_u_ndm(ip)) == 4)  uv_flags_adv(ip) = 0
@@ -1176,7 +1176,7 @@ subroutine read_wavemaker_polylines()
       rewind(500)
       allocate(x_wmfp(nwmfp))
       allocate(y_wmfp(nwmfp))
-      do n = 1, nwbnd
+      do n = 1, nwmfp
          read(500,*)x_wmfp(n),y_wmfp(n)
       enddo
       close(500)