Combine open PRs for v2.0.7 (#92)

* Merge of last small commits that were not already merged in the separate PRs, and cherrypick the viscosity improvements without importing the quadtree infiltration changes

* - Set new default combination that since new advection scheme seems robust: alpha=0.50, theta=1.0, adveciont=1, viscosity=1

* - Fix quadtree hmax output

* - set/check random seed for random_number for wave_maker
- Todo: check whether this is best implementation and results in actual random seed if wmrandom = 1

* - SnapWave doesn't output directional spreading, and in sfincs_snapwave the mean wave direction is actually added, which gives wrong values in hisfile. Turn off for now

* Cherrypick from PR#89/branch 88 determine viscosity per refinement level for quadtree simulation - so quadtree infiltration is not imported too:

- Move whole viscosity value nuvisc determination to sfincs_domain.f90
   - Now determine nuvisc per grid cell, depending on refinement level (min of dx-dy)!
   - crsgeo true or false is still considered
   - removed option for direct 'nuvisc' input, still possible to scale through 'nuviscdim'

* - Fix '-1' mistake for using correct dx/dy when determining viscosity

* - Index was wrong in wavemaker for case subgrid

* - Add again fix in quadtree for weirs, couldn't find back original commit to cherry pick

* Cherrypick from 80-add-option-for-netcdf-infiltration-quadtree-file-input:

- Remove double write cuminf

- add needed 'if (nm>0) ' checks for quadtree max output

* - Cumprcp is already in meters now !? > so then no conversion from mm to m needed anymore
* - Units of 'cumprcp' are changed from mm to m, just like cuminf output, now change it in the netcdf attribute too

* - Fix missing point_prcp because of meteo3d that was turned to false
- fix debug mode issues with extra (nm > 0) checks

* - Give more useful Error where depth SnapWave input data is too shallow by giving the x&y-locations

* - After review Maarten change implementation so that nuvisc is a value per level, not grid cell, to save computation time in sfincs_momentum

* - Change input specification for nuvisc, so we don't put in a default value of 1.0 that is later divided by 100, but make this 0.01 directly
- This is the 'viscosity per meter grid cell length'
- Note, to prevent issues with old input, I use in sfincs.inp the 'nuvisc' keyword instead of 'nuviscdim', so people with nuviscdim=1 in their older sfincs.inp suddenly do not have a lot of visco, but have the same default

* - Update the documentation as well

Author: Leynse

docs/input.rst

@@ -717,23 +717,17 @@ Recommended is to turn the advection term always on.
 
 **viscosity**
 
-'viscosity' turns on the viscosity term in the momentum equation (viscosity = 1, default).
-The recommended value of viscosity 'nuvisc' to add to your model, depends on your grid size.
-For ease, SFINCS internally automatically determines the optimal value for you, which is displayed when running the model:	'Turning on process: Viscosity, with nuvisc=   0.5000000'. In this example corresponding to a grid resolution of 50 meters.
-In case you would want to increase the viscosity term, you can either specify the exact value you want 'nuvisc = XXX', or e.g. multiply it by a factor 2: 	nuviscdim = 2.0 (default = 1.0, dimensionless).
-By default the value of nuvisc is determined like this:
-
-	dx = 50 > nuvisc = 0.5
-	
-        dx = 100 > nuvisc = 1.0
-	
-        dx = 500 > nuvisc = 5.0	
+'viscosity' turns on the viscosity term in the momentum equation (viscosity = 1).
+The recommended value of 'nuvisc', the viscosity coefficient 'per meter of grid cell length' to add to your model is 0.01.
+This coefficient is multiplied internally with the grid cell size (per quadtree level in quadtree mesh mode).
+For ease, this is displayed when running the model:	'Viscosity - nuvisc=   [0.1000000 - 0.5000000]'. 
+Increasing the value of 'nuvisc' increases the viscosity term and effectively internal smoothing of the flow.
 
 .. code-block:: text
 
 	viscosity 	= 1
-	nuviscdim 	= 1.0 (default)
-	nuvisc 		= XXX (automatically determined, or specify a value yourself that overrules this)
+	nuvisc 		= 0.01
+	nuviscdim 	= Depricated after Cauberg release of SFINCS.
 	
 **Drag Coefficients:**
 

docs/parameters.rst

@@ -76,21 +76,18 @@ Parameters for model input
 	  :min:			0.8
 	  :max:			1.0
 	viscosity
-	  :description:		Turns on the viscosity term in the momentum equation (viscosity = 1), advised to combine with theta = 1.0. Value of viscosity term 'nuvisc' automatically determined based on grid size..
+	  :description:		Turns on the viscosity term in the momentum equation (viscosity = 1), advised to combine with theta = 1.0.
 	  :units:		-
 	  :default:		1
 	  :min:			0
 	  :max:			1
 	nuviscdim
-	  :description:		Dimensionless viscosity coefficient, multiplies the automatically determined value for 'nuvisc' with the specified factor for 'nuviscdim'.
-	  :units:		-
-	  :default:		1.0
-	  :min:			0.0
-	  :max:			Inf	  
+	  :description:		Depricated after Cauberg release of SFINCS.
+	  :units:		-  
 	nuvisc
-	  :description:		Viscosity coefficient, by default turned off, but automatically determined if 'viscosity=1'. specifying a value for 'nuvisc' overrule default value.
+	  :description:		Viscosity coefficient 'per meter of grid cell length', used if 'viscosity=1' and multiplied internally with the grid cell size (per quadtree level in quadtree mesh mode).
 	  :units:		-
-	  :default:		-999.0 (=off)
+	  :default:		0.01
 	  :min:			0.0
 	  :max:			Inf	  	  	  
 	zsini

source/src/quadtree.F90

@@ -852,13 +852,19 @@ subroutine find_uv_points_intersected_by_polyline(uv_indices,vertices,nint,xpol,
          do m = m0, m1
             do n = n0, n1            
                !
+               ! First check whether point is on the quadtree grid at all
                nm = find_quadtree_cell_by_index(n, m, iref)
                !
                if (nm==0) then
                   cycle
                endif   
                !
+               ! Second check is whether it is on the active SFINCS mask
                nm = index_sfincs_in_quadtree(nm)
+               !
+               if (nm==0) then
+                  cycle
+               endif        
                ! 
                ! Right (same level or coarser)
                !

source/src/sfincs_data.f90

@@ -72,7 +72,6 @@ module sfincs_data
       real*4 sfacinf
       real*4 dym
       real*4 dym2
-      real*4 nuvisc
       real*4 nuviscdim
       real*4 nuviscinp   
       real*4 spw_merge_frac
@@ -316,6 +315,7 @@ module sfincs_data
       real*4, dimension(:),   allocatable :: dxrinvc
       real*4, dimension(:),   allocatable :: dyrinvc
       real*4, dimension(:),   allocatable :: cell_area
+      real*4, dimension(:),   allocatable :: dxyr      
       !
       ! Cell sizes
       !
@@ -328,6 +328,7 @@ module sfincs_data
       real*4,             dimension(:),   allocatable :: z_xz
       real*4,             dimension(:),   allocatable :: z_yz
       real*4,             dimension(:),   allocatable :: cell_area_m2
+      real*4,             dimension(:),   allocatable :: nuvisc      
       !
       ! UV-points
       !
@@ -844,7 +845,7 @@ subroutine finalize_parameters()
     if(allocated(qinffield)) deallocate(qinffield)
     if(allocated(ksfield)) deallocate(ksfield)
     if(allocated(scs_Se)) deallocate(scs_Se)
-
+    if(allocated(nuvisc)) deallocate(nuvisc)    
     !
     ! Boundary velocity points
     !

source/src/sfincs_domain.f90

@@ -88,7 +88,7 @@ subroutine initialize_processes()
    !
    ! Check for time-spatially varying meteo
    !
-   if (spwfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. amprfile(1:4) /= 'none' .or. ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. netamprfile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then
+   if (prcpfile(1:4) /= 'none' .or. spwfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. amprfile(1:4) /= 'none' .or. ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. netamprfile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then
       meteo3d = .true.
       write(*,*)'Turning on flag: meteo3d'
    endif
@@ -1258,13 +1258,15 @@ subroutine initialize_mesh()
    !
    allocate(dxr(nref))
    allocate(dyr(nref))
+   allocate(dxyr(nref))
    !
    ! Determine grid spacing for different refinement levels
    !
    do iref = 1, nref
       !
       dxr(iref) = dx/(2**(iref - 1))
       dyr(iref) = dy/(2**(iref - 1))
+      dxyr(iref) = min(dxr(iref), dyr(iref))      
       !
    enddo 
    !
@@ -1413,6 +1415,47 @@ subroutine initialize_mesh()
       !
    endif   
    !
+   ! Determine viscosity per refinement level
+   !
+   allocate(nuvisc(nref))  
+   !
+   if (viscosity) then
+        !
+        if (crsgeo) then
+            !
+            do iref = 1, nref
+                !
+                nuvisc(iref) = max(nuviscdim * dxyr(iref) / 0.00001, 0.0) ! take min of dx and dy, don't allow to be negative 
+                !
+                ! with default nuviscdim = 0.01 (dimensionless per meter grid cell length):                
+                ! dx = 1 degree ~ 100km > nuvisc = 1000.0
+                ! dx = 0.001 degree~ 100m > nuvisc = 1.0
+                !
+            enddo           
+            !
+        else
+            !
+            do iref = 1, nref
+                !       
+                nuvisc(iref) = max(nuviscdim * dxyr(iref), 0.0) ! take min of dx and dy, don't allow to be negative    
+                !
+                ! with default nuviscdim = 0.01 (dimensionless per meter grid cell length): 
+                ! dx = 50 > nuvisc = 0.5
+                ! dx = 100 > nuvisc = 1.0
+                ! dx = 500 > nuvisc = 5.0
+                !
+            enddo                    
+            ! 
+        endif
+        !
+        if (use_quadtree) then
+            write(*,*)'Viscosity - nuvisc  = [', minval(nuvisc),'- ',maxval(nuvisc),']'           
+        else
+            write(*,*)'Viscosity - nuvisc  = ', maxval(nuvisc)           
+        endif            
+        !
+   endif   
+   !
    end subroutine
 
    subroutine initialize_bathymetry()

source/src/sfincs_input.f90

@@ -104,9 +104,8 @@ subroutine read_sfincs_input()
    call read_real_input(500,'spwmergefrac',spw_merge_frac,0.5)
    call read_int_input(500,'usespwprecip',ispwprecip,1)   
    call read_int_input(500,'global',iglobal,0)
-   call read_real_input(500,'nuvisc',nuviscinp,-999.0)
-   call read_real_input(500,'nuviscdim',nuviscdim,1.0)   
-   call read_logical_input(500,'viscosity',iviscosity,.true.)
+   call read_real_input(500,'nuvisc',nuviscdim,0.01)      
+   call read_logical_input(500,'viscosity',iviscosity,.false.)
    call read_int_input(500,'spinup_meteo', ispinupmeteo, 0)
    call read_real_input(500,'waveage',waveage,-999.0)
    call read_int_input(500,'snapwave', isnapwave, 0)
@@ -423,41 +422,11 @@ subroutine read_sfincs_input()
       store_tsunami_arrival_time = .true.
    endif      
    !
-   viscosity = .false.   
-   !
-   if (nuviscinp>0.0 .or. iviscosity) then
-      !
-      if (nuviscinp>0.0) then ! if nuvisc given by user, this overrules the dimensionless default use of nuvisc, for backwards compatability
-         !
-         nuvisc = nuviscinp
-         viscosity = .true.         
-         !
-      else ! user defined viscosity = 1
-         !
-         if (qtrfile(1:4) == 'none') then ! this works only for a regular grid model, for quadtree use 'nuvisc' option
-            !
-            viscosity = .true.         
-            !
-            if (crsgeo .eqv. .true.) then ! simplified conversion for spherical grids
-                nuvisc = max(nuviscdim * min(dx,dy) / 0.001, 0.0) ! take min of dx and dy, don't allow to be negative  
-                ! dx = 1 degree ~ 100km
-                ! dx = 0.001 degree~ 100m > nuvisc = 1.0              
-            else                 
-                nuvisc = max(nuviscdim * min(dx,dy) / 100.0, 0.0) ! take min of dx and dy, don't allow to be negative    
-                ! dx = 50 > nuvisc = 0.5
-                ! dx = 100 > nuvisc = 1.0
-                ! dx = 500 > nuvisc = 5.0
-                ! nuviscdim = 1.0
-                ! nuvisc = nuviscdim * dx / 100 
-            endif
-         endif          
-         !
-      endif    
-      !
-      if (viscosity .eqv. .true.) then
-         write(*,*)'Turning on process: Viscosity, with nuvisc = ', nuvisc
-      endif      
-      !
+   !   
+   viscosity = .false. 
+   if (iviscosity) then
+      viscosity = .true. 
+      write(*,*)'Turning on process: Viscosity'               
    endif   
    !
    spinup_meteo = .true. ! Default use data in ampr file as block rather than linear interpolation

source/src/sfincs_lib.f90

@@ -85,7 +85,7 @@ function sfincs_initialize(config_file) result(ierr)
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
-   build_revision = '$Rev: v2.0.6'
+   build_revision = '$Rev: v2.0.7'
    build_date     = '$Date: 2024-05-31'
    !
    write(*,'(a)')''   

source/src/sfincs_momentum.f90

@@ -449,7 +449,7 @@ subroutine compute_fluxes(dt, min_dt, tloop)
             !
             if (viscosity) then
                !
-               frc = frc + nuvisc * hu * ( (uu_nmu - 2*uu_nm + uu_nmd ) * dxuv2inv + (uu_num - 2*uu_nm + uu_ndm ) * dyuv2inv )
+               frc = frc + nuvisc(iref) * hu * ( (uu_nmu - 2*uu_nm + uu_nmd ) * dxuv2inv + (uu_num - 2*uu_nm + uu_ndm ) * dyuv2inv )
                !
             endif
             !