Fix instability in RBM

src/Begin.f90

@@ -13,11 +13,11 @@ Subroutine BEGIN(param_file,spatial_file)
     character (len=200):: param_file,source_file,spatial_file
     !
     integer:: Julian
-    integer:: head_name,trib_cell
+    integer:: head_name,trib_cell
     integer:: jul_start,main_stem,nyear1,nyear2,nc,ncell,nseg
     integer:: ns_max_test,node,ncol,nrow,nr,cum_sgmnt
     !
-    integer:: nreservoir
+    integer:: nreservoir,nseg_temp,nseg_cum
     !
     logical:: first_cell,source
     !
@@ -83,6 +83,7 @@ Subroutine BEGIN(param_file,spatial_file)
     allocate(res_start_node(nres))
     allocate(res_end_node(nres))
     allocate(res_capacity_mcm(nres))
+    allocate(nseg_out(nreach,heat_cells,nseg_out_num))
     !
     !     Start reading the reach date and initialize the reach index, NR
     !     and the cell index, NCELL
@@ -117,12 +118,13 @@ Subroutine BEGIN(param_file,spatial_file)
         !     Initialize NSEG, the total number of segments in this reach
         !
         nseg=0
+        nseg_cum=0
         write(*,*) ' Starting to read reach ',nr
         !
         !     Read the number of cells in this reach, the headwater #,
         !     the number of the cell where it enters the next higher order stream,
         !     the headwater number of the next higher order stream it enters, and
-        !     the river mile of the headwaters.
+        !     the river mile of the headwaters.
         !
         read(90,'(i5,11x,i4,10x,i5,15x,i5,15x,f10.0,i5)') no_cells(nr) &
             ,head_name,trib_cell,main_stem,rmile0
@@ -139,7 +141,7 @@ Subroutine BEGIN(param_file,spatial_file)
         if (trib_cell.gt.0) then
             no_tribs(trib_cell) = no_tribs(trib_cell)+1
             trib(trib_cell,no_tribs(trib_cell)) = nr
-        end if
+        end if
         !
         !     Reading Mohseni parameters for each headwaters (UW_JRY_2011/06/18)
         !
@@ -165,7 +167,7 @@ Subroutine BEGIN(param_file,spatial_file)
             end if
             !
             !     The headwaters index for each cell in this reach is given
-            !     in the order the cells are read
+            !     in the order the cells are read
             !
             !     Card Type 3. Cell indexing #, Node # Row # Column Lat Long RM
             !
@@ -175,14 +177,14 @@ Subroutine BEGIN(param_file,spatial_file)
             if (reservoir) then
                 read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0,i6)')  &
                     node,nrow,ncol,lat,long,rmile1,ndelta(ncell),res_num(ncell)
-                write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell),res_num(ncell)
+                !write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell),res_num(ncell)
                 if(res_num(ncell) .gt. 0) then
                     res_pres(ncell) = .TRUE.
                 end if
             else
                 read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0)')  &
                     node,nrow,ncol,lat,long,rmile1,ndelta(ncell)
-                write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell)
+                !write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell)
             end if
             !
             !    Set the number of segments of the default, if not specified
@@ -197,9 +199,16 @@ Subroutine BEGIN(param_file,spatial_file)
             ! Added variable ndelta (UW_JRY_2011/03/15)
             !
             dx(ncell)=miles_to_ft*(rmile0-rmile1)/ndelta(ncell)
-            rmile0=rmile1
-            nndlta=0
-200     continue
+            rmile0=rmile1
+            !
+            ! Here we define the output segments
+            !
+            do nseg_temp=1,nseg_out_num
+                nseg_out(nr,ncell,nseg_temp)=nseg_cum+ndelta(ncell)*nseg_temp/(nseg_out_num) 
+            end do
+            nseg_cum = nseg_cum+ndelta(ncell)
+            nndlta=0
+200     continue
         nndlta=nndlta+1
         nseg=nseg+1
         segment_cell(nr,nseg)=ncell
@@ -208,20 +217,25 @@ Subroutine BEGIN(param_file,spatial_file)
         !
         !   Write Segment List for mapping to temperature output (UW_JRY_2008/11/19)
         !
-        open(22,file=TRIM(spatial_file),status='unknown') ! (changed by WUR_WF_MvV_2011/01/05)
-        write(22,'(4i6,1x,a8,1x,a10,f5.0)') nr,ncell,nrow,ncol,lat,long,nndlta
+        do nseg_temp=1,nseg_out_num
+            if (nseg_out(nr,ncell,nseg_temp).eq.nseg) then
+                open(22,file=TRIM(spatial_file),status='unknown') ! (changed by WUR_WF_MvV_2011/01/05)
+                write(22,'(4i6,1x,a8,1x,a10,i5)') nr,ncell,nrow,ncol,lat,long,nseg_temp
+            end if
+        end do
         !
         !
         !
         !  Added variable ndelta  (UW_JRY_2011/03/15)
         !
         if(nndlta.lt.ndelta(ncell)) go to 200
-        no_celm(nr)=nseg
+        no_celm(nr)=nseg
         segment_cell(nr,nseg)=ncell
-        x_dist(nr,nseg)=miles_to_ft*rmile1
+        x_dist(nr,nseg)=miles_to_ft*rmile1
+        write(*,*) 'number of segment in reach', nr, nseg
     !
     ! End of cell and segment loop
-    !
+    !
     end do
     !
     ! If this is a reach that is tributary to another, set the confluence cell to the previous
@@ -233,26 +247,26 @@ Subroutine BEGIN(param_file,spatial_file)
         conflnce(trib_cell,no_tribs(trib_cell)) = ncell
     end if
 
-    if(ns_max_test.lt.nseg) ns_max_test=nseg
-!
+    if(ns_max_test.lt.nseg) ns_max_test=nseg
+!
 ! End of reach loop
-!
+!
 end do
 if(ns_max_test.gt.ns_max) then
     write(*,*) 'RBM is terminating because'
     write(*,*) 'NS_MAX exceeded. Change NS_MAX in Block_Network to: ',ns_max_test
     stop
-end if
-!
-nwpd=1
-xwpd=nwpd
-dt_comp=86400./xwpd
+end if
+!
+nwpd=1
+xwpd=nwpd
+dt_comp=86400./xwpd
+!
+!     ******************************************************
+!                         Return to RMAIN
+!     ******************************************************
 !
-!     ******************************************************
-!                         Return to RMAIN
-!     ******************************************************
-!
-900 continue
+900 continue
 !
 !
 end subroutine BEGIN

src/Block_Energy.f90

@@ -32,11 +32,14 @@ module Block_Energy
 !   Some important constants
 !
       real             :: lvp,rb,rho
+      real             :: deriv_2nd,error_EE
+      real             :: deriv_conv,deriv_evap,deriv_ws 
+      real             :: temp_equil,time_equil
       real,parameter   :: evap_coeff=1.5e-9 !Lake Hefner coefficient, 1/meters
       real,parameter   :: pf=0.640,pi=3.14159
       real,parameter   :: rfac=304.8 !rho/Cp kg/meter**3/Kilocalories/kg/Deg K  
       real,parameter   :: sec_day = 86400 !number of seconds in a day
       real,parameter   :: a_z=0.408, b_z=0.392 !Leopold parameter for depth
-      real,parameter   :: a_w=4.346, b_w=0.520 !Leopold parameter for width
+      real,parameter   :: a_w=4.346, b_w=0.520 !Leopold parameter for width
 !
 end module Block_Energy  

src/Block_Hydro.f90

@@ -2,8 +2,8 @@
 ! Module for hydraulic characteristics and water quality constituents of the basin
 !
 module Block_Hydro
-    integer, dimension(2000):: no_dt,nstrt_elm
-    real, dimension(2000)   :: dt_part,x_part
+    integer, dimension(2500):: no_dt,nstrt_elm
+    real, dimension(2500)   :: dt_part,x_part
 !
     real, dimension(:),   allocatable  :: depth
     real, dimension(:),   allocatable  :: width

src/Block_Network.f90

@@ -7,18 +7,29 @@ Module Block_Network
 !
     integer, dimension(:,:), allocatable::conflnce,reach_cell,segment_cell,trib
 !
+    integer, dimension(:,:,:), allocatable::nseg_out
 !
 ! Integer variables 
 !
     integer:: flow_cells,heat_cells
     integer:: ndays,nreach,ntrb,nwpd
-    integer,parameter::ns_max=200
+    integer,parameter::ns_max=3000
+    integer,parameter::nseg_out_num=2
     integer:: start_year,start_month,start_day
     integer:: end_year,end_month,end_day
+    integer:: numsub !number of subdaily timestep
+    integer:: nsub
 !
 ! Real variables
 !
     real   :: delta_n,n_default=2
     real   :: dt_comp
+    real   :: dt_res
     real, dimension(:), allocatable  :: ndelta
+!
+!      Logical variables
+!
+
+
+
 end module Block_Network

src/Block_Reservoir.f90

@@ -12,6 +12,9 @@ module Block_Reservoir
     logical, dimension(:),   allocatable  :: res_trib_calc !
     logical, dimension(:),   allocatable  :: res_stratif_start !
     logical, dimension(:),   allocatable  :: res_turnover !
+    logical :: exceed_error_bound
+    logical :: adjust_timestep
+    logical :: recalculate_volume
     !
     !
     real, parameter :: depth_e_frac=0.4, depth_h_frac=0.6
@@ -23,6 +26,13 @@ module Block_Reservoir
     real :: flow_epi_hyp_x, outflow_x
     real :: res_vol_delta_x, vol_change_hyp_x, vol_change_epi_x
     real :: res_storage_pre, res_storage_post
+    ! 
+    !     Error estimation
+    !
+    real :: m11,m12,m21,m22
+    real :: const1,const2
+    real :: error_e,error_h
+    real, parameter :: error_threshold = 0.1
     real, dimension(:),   allocatable  :: K_z
     real, dimension(:),   allocatable  :: depth_e, depth_h
     real, dimension(:),   allocatable  :: density_epil, density_hypo
@@ -54,6 +64,7 @@ module Block_Reservoir
     real, dimension(:),   allocatable  :: qsurf_tot
     real, dimension(:),   allocatable  :: res_capacity_mcm
     real, dimension(:),   allocatable  :: volume_h_min
+    real, dimension(:),   allocatable  :: volume_e_min
     real, dimension(:,:),   allocatable  :: res_storage
     !
     !