fixed possibility for divide by zero in snowcomp init for dividing by…

… pk_depth = 0 and den_init value specified = 0

added one_subbasin_flag control parameter to allow simulation of a single subbasin designated by this parameter

GSFLOW/src/gsflow/gsflow_module.f90

@@ -64,7 +64,7 @@ MODULE PRMS_MODULE
       INTEGER, SAVE :: NhruOutON_OFF, Gwr_swale_flag, NsubOutON_OFF, BasinOutON_OFF, NsegmentOutON_OFF
       INTEGER, SAVE :: Stream_temp_flag, Strmtemp_humidity_flag, Stream_temp_shade_flag
       INTEGER, SAVE :: Prms_warmup, statsON_OFF
-      INTEGER, SAVE :: Frozen_flag, Glacier_flag
+      INTEGER, SAVE :: Frozen_flag, Glacier_flag, one_subbasin_flag
       INTEGER, SAVE :: PRMS_land_iteration_flag, Iter_aet_flag, text_restart_flag
       INTEGER, SAVE :: irrigation_apply_flag, Dyn_ag_frac_flag, Dyn_ag_soil_flag, activeHRU_inactiveCELL_flag
       INTEGER, SAVE :: Dprst_add_water_use, Dprst_transfer_water_use

GSFLOW/src/gsflow/gsflow_prms.f90

@@ -1019,6 +1019,7 @@ SUBROUTINE setdims(AFR, Diversions, Idivert, EXCHANGE, DELTAVOL, LAKEVOL, Nsegsh
 ! subbasin dimensions
       IF ( control_integer(Subbasin_flag, 'subbasin_flag')/=0 ) Subbasin_flag = OFF
       IF ( decldim('nsub', 0, MAXDIM, 'Number of internal subbasins')/=0 ) CALL read_error(7, 'nsub')
+      IF ( control_integer(Subbasin_flag, 'one_subbasin_flag')/=0 ) one_subbasin_flag = OFF
 
       IF ( control_integer(Dprst_flag, 'dprst_flag')/=0 ) Dprst_flag = OFF
       IF ( control_integer(Dprst_transfer_water_use, 'dprst_transfer_water_use')/=0 ) Dprst_transfer_water_use = OFF
@@ -1271,6 +1272,8 @@ INTEGER FUNCTION check_dims(Nsegshold, Nlakeshold)
       ! default = 1, turn off if no subbasins
       IF ( Subbasin_flag==1 .AND. Nsub==0 ) Subbasin_flag = 0
 
+      IF ( one_subbasin_flag==1 .AND. Nsub==0 ) one_subbasin_flag = 0
+
       Nsegment = getdim('nsegment')
       IF ( Nsegment==-1 ) CALL read_error(7, 'nsegment')
       Nsegshold = Nsegment

GSFLOW/src/gsflow/gsflow_prms2mf.f90

@@ -153,7 +153,7 @@ END FUNCTION prms2mfdecl
 !***********************************************************************
       INTEGER FUNCTION prms2mfinit()
       USE PRMS_CONSTANTS, ONLY: DEBUG_less, ERROR_param, ACTIVE, OFF
-      USE PRMS_MODULE, ONLY: Hru_type
+      USE PRMS_MODULE, ONLY: Hru_type, one_subbasin_flag
       use PRMS_READ_PARAM_FILE, only: getparam_real
       USE GSFPRMS2MF
       USE GWFUZFMODULE, ONLY: NTRAIL, NWAV, IUZFBND
@@ -165,11 +165,11 @@ INTEGER FUNCTION prms2mfinit()
       USE PRMS_BASIN, ONLY: Active_hrus, Hru_route_order, Basin_area_inv, Hru_area
       USE PRMS_SOILZONE, ONLY: Gvr_hru_id, Gvr_hru_pct_adjusted
       use prms_utils, only: read_error
-      USE GLOBAL, ONLY: NLAY, NROW, NCOL
+      USE GLOBAL, ONLY: NLAY, NROW, NCOL, IBOUND
       IMPLICIT NONE
       INTRINSIC :: ABS, DBLE
 ! Local Variables
-      INTEGER :: is, i, ii, ierr, ihru, icell, irow, icol
+      INTEGER :: is, i, ii, ierr, ihru, icell, irow, icol, jj
 !     INTEGER :: iseg, max_seg, irch
 !     INTEGER, ALLOCATABLE, DIMENSION(:) :: nseg_rch
       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: hru_pct, newpct
@@ -312,8 +312,16 @@ INTEGER FUNCTION prms2mfinit()
         icol = Gwc_col(icell)
         IF ( Print_debug>DEBUG_less ) THEN
           IF ( Hru_type(ihru)==0 ) THEN
-            IF ( IUZFBND(icol, irow)/=0 ) &
-     &           PRINT *, 'WARNING, HRU inactive & UZF cell active, irow:', irow, 'icell:', icell, ' HRU:', ihru
+            IF ( IUZFBND(icol, irow)/=0 ) THEN
+              IF ( one_subbasin_flag==0 ) THEN
+                PRINT *, 'WARNING, HRU inactive & UZF cell active, irow:', irow, 'icell:', icell, ' HRU:', ihru
+              ELSE
+                IUZFBND(icol, irow) = 0
+                DO jj = 1, NLAY
+                  IBOUND(icol, irow, jj) = 0
+                ENDDO
+              ENDIF
+            ENDIF
           ENDIF
           IF ( IUZFBND(icol, irow)==0 ) THEN
             IF ( Hru_type(ihru) > 0 ) then

GSFLOW/src/prms/basin.f90

@@ -32,6 +32,7 @@ MODULE PRMS_BASIN
       REAL, SAVE, ALLOCATABLE :: Hru_area(:), Hru_elev(:), Hru_lat(:) ! , Hru_percent_imperv(:)
       REAL, SAVE, ALLOCATABLE :: Covden_sum(:), Covden_win(:)
       REAL, SAVE, ALLOCATABLE :: Dprst_frac_open(:), Dprst_area(:) !, Dprst_frac(:)
+      INTEGER, SAVE, ALLOCATABLE :: Hru_subbasin(:)
       REAL, SAVE, ALLOCATABLE :: Ag_frac(:)
       END MODULE PRMS_BASIN
 
@@ -65,7 +66,7 @@ END FUNCTION basin
       INTEGER FUNCTION basdecl()
       USE PRMS_CONSTANTS, ONLY: ACTIVE, OFF
       USE PRMS_MODULE, ONLY: Nhru, Nlake, Dprst_flag, Lake_route_flag, &
-     &    PRMS4_flag, gwflow_flag, Glacier_flag, AG_flag, activeHRU_inactiveCELL_flag
+     &    PRMS4_flag, gwflow_flag, Glacier_flag, AG_flag, activeHRU_inactiveCELL_flag, Nsub
       use PRMS_MMFAPI, only: declvar_real, declvar_dble
       use PRMS_READ_PARAM_FILE, only: declparam
       USE PRMS_BASIN
@@ -281,6 +282,15 @@ INTEGER FUNCTION basdecl()
         ALLOCATE ( Lake_area(1) )
       ENDIF
 
+      IF ( Nsub>0 ) THEN
+        ALLOCATE ( Hru_subbasin(Nhru) )
+        IF ( declparam(MODNAME, 'hru_subbasin', 'nhru', 'integer', &
+     &       '0', 'bounded', 'nsub', &
+     &       'Index of subbasin assigned to each HRU', &
+     &       'Index of subbasin assigned to each HRU', &
+     &       'none')/=0 ) CALL read_error(1, 'hru_subbasin')
+      ENDIF
+
       END FUNCTION basdecl
 
 !**********************************************************************
@@ -295,7 +305,7 @@ INTEGER FUNCTION basinit()
       USE PRMS_MODULE, ONLY: Nhru, Nlake, Print_debug, Hru_type, irrigation_apply_flag, &
      &    Dprst_flag, Lake_route_flag, PRMS4_flag, gwflow_flag, PRMS_VERSION, &
      &    Starttime, Endtime, Parameter_check_flag, Inputerror_flag, &
-     &    AG_flag, Ag_package, activeHRU_inactiveCELL_flag !, Frozen_flag
+     &    AG_flag, Ag_package, activeHRU_inactiveCELL_flag, Nsub, one_subbasin_flag !, Frozen_flag
       USE PRMS_BASIN
       use prms_utils, only: checkdim_bounded_limits, error_stop, read_error, write_outfile
       IMPLICIT NONE
@@ -376,6 +386,9 @@ INTEGER FUNCTION basinit()
       ELSE
         Lake_hru_id = 0
       ENDIF
+      IF ( Nsub>0 ) THEN
+        IF ( getparam_int(MODNAME, 'hru_subbasin', Nhru, Hru_subbasin)/=0 ) CALL read_error(2, 'hru_subbasin')
+      ENDIF
 
       Basin_gl_cfs = 0.0D0
       Basin_gl_ice_cfs = 0.0D0
@@ -411,6 +424,11 @@ INTEGER FUNCTION basinit()
         perv_area = harea
 
         IF ( Hru_type(i)==INACTIVE ) CYCLE
+
+        IF ( one_subbasin_flag>0 ) THEN
+          IF ( Hru_subbasin(i) /= one_subbasin_flag ) Hru_type(i) = 0
+        ENDIF
+
 ! ????????? need to fix for lakes with multiple HRUs and PRMS lake routing ????????
         IF ( Hru_type(i)==LAKE ) THEN
           Numlake_hrus = Numlake_hrus + 1

GSFLOW/src/prms/nsub_summary.f90

@@ -20,8 +20,6 @@ MODULE PRMS_NSUB_SUMMARY
       INTEGER, SAVE, ALLOCATABLE :: Monthlyunit(:), Yearlyunit(:)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Nsub_var_monthly(:, :), Nsub_var_yearly(:, :)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Sub_area(:)
-!   Declared Parameters
-      INTEGER, SAVE, ALLOCATABLE :: Hru_subbasin(:)
 ! Control Parameters
       INTEGER, SAVE :: NsubOutVars, NsubOut_freq, NsubOut_format
       CHARACTER(LEN=36), SAVE, ALLOCATABLE :: NsubOutVar_names(:)
@@ -70,7 +68,7 @@ SUBROUTINE nsub_summarydecl()
       USE PRMS_CONSTANTS, ONLY: ERROR_control, DAILY, YEARLY
       use PRMS_CONTROL_FILE, only: control_integer, control_string, control_string_array
       use PRMS_READ_PARAM_FILE, only: declparam
-      USE PRMS_MODULE, ONLY: Nhru, Nsub
+      USE PRMS_MODULE, ONLY: Nsub
       USE PRMS_NSUB_SUMMARY
       use prms_utils, only: error_stop, print_module, read_error
       IMPLICIT NONE
@@ -99,12 +97,7 @@ SUBROUTINE nsub_summarydecl()
         IF ( control_string(NsubOutBaseFileName, 'nsubOutBaseFileName')/=0 ) CALL read_error(5, 'nsubOutBaseFileName')
       ENDIF
 
-      ALLOCATE ( Hru_subbasin(Nhru), Sub_area(Nsub) )
-      IF ( declparam(MODNAME, 'hru_subbasin', 'nhru', 'integer', &
-     &     '0', 'bounded', 'nsub', &
-     &     'Index of subbasin assigned to each HRU', &
-     &     'Index of subbasin assigned to each HRU', &
-     &     'none')/=0 ) CALL read_error(1, 'hru_subbasin')
+      ALLOCATE ( Sub_area(Nsub) )
 
       END SUBROUTINE nsub_summarydecl
 
@@ -118,7 +111,7 @@ SUBROUTINE nsub_summaryinit()
       use PRMS_READ_PARAM_FILE, only: getparam_int
       USE PRMS_MODULE, ONLY: Nhru, Nsub, Inputerror_flag, Start_year, Prms_warmup
       USE PRMS_NSUB_SUMMARY
-      USE PRMS_BASIN, ONLY: Hru_area_dble, Active_hrus, Hru_route_order
+      USE PRMS_BASIN, ONLY: Hru_area_dble, Active_hrus, Hru_route_order, Hru_subbasin
       use prms_utils, only: error_stop, numchars, PRMS_open_output_file, read_error
       IMPLICIT NONE
 ! Local Variables
@@ -255,7 +248,6 @@ SUBROUTINE nsub_summaryinit()
         ENDIF
       ENDDO
 
-      IF ( getparam_int(MODNAME, 'hru_subbasin', Nhru, Hru_subbasin)/=0 ) CALL read_error(2, 'hru_subbasin')
       Sub_area = 0.0D0
       DO i = 1, Active_hrus
         j = Hru_route_order(i)
@@ -291,7 +283,7 @@ SUBROUTINE nsub_summaryrun()
       use PRMS_MMFAPI, only: getvar_dble, getvar_real
       USE PRMS_MODULE, ONLY: Nhru, Nsub, Start_month, Start_day, End_year, End_month, End_day, Nowyear, Nowmonth, Nowday
       USE PRMS_NSUB_SUMMARY
-      USE PRMS_BASIN, ONLY: Active_hrus, Hru_route_order, Hru_area_dble
+      USE PRMS_BASIN, ONLY: Active_hrus, Hru_route_order, Hru_area_dble, Hru_subbasin
       USE PRMS_SET_TIME, ONLY: Modays
       use prms_utils, only: read_error
       IMPLICIT NONE

GSFLOW/src/prms/sm_read_control_file.f90

@@ -369,6 +369,10 @@ module subroutine setup_cont()
     Subbasin_flag = ACTIVE
     Control_parameter_data(i) % values_int(1) = Subbasin_flag
     i = i + 1
+    Control_parameter_data(i) % name = 'one_subbasin_flag'
+    one_subbasin_flag = OFF
+    Control_parameter_data(i) % values_int(1) = one_subbasin_flag
+    i = i + 1
     Control_parameter_data(i) % name = 'cbh_check_flag'
     Cbh_check_flag = OFF
     Control_parameter_data(i) % values_int(1) = Cbh_check_flag

GSFLOW/src/prms/snowcomp.f90

@@ -773,14 +773,18 @@ INTEGER FUNCTION snoinit()
       IF ( Init_vars_from_file==0 .OR. Init_vars_from_file==2 .OR. Init_vars_from_file==3 ) THEN
         IF ( getparam_real(MODNAME, 'snowpack_init', Nhru, Snowpack_init)/=0 ) CALL read_error(2, 'snowpack_init')
         Pkwater_equiv = DBLE( Snowpack_init )
-        Pk_depth = Pkwater_equiv/DBLE( Den_init )
-        Pk_den = SNGL( Pkwater_equiv/Pk_depth )
+        Pk_depth = 0.0
+        Pk_den = 0.0
         Pk_ice = SNGL( Pkwater_equiv )
         Freeh2o = Pk_ice*Freeh2o_cap
         Ai = 0.0D0
         Snowcov_area = 0.0
         DO j = 1, Active_hrus
           i = Hru_route_order(j)
+          IF ( .not.(Den_init(i)>0.0) ) THEN
+            IF ( Print_debug>-1 ) PRINT *, 'WARNING, den_init for HRU:', i, ' specified = 0.0, set to 0.1'
+            Den_init(i) = 0.1 ! to avoid divide by zero potential
+          ENDIF
           IF ( Pkwater_equiv(i)>0.0D0 ) THEN
             Ai(i) = Pkwater_equiv(i) ! [inches]
             IF ( Ai(i)>snarea_thresh_dble(i) ) Ai(i) = snarea_thresh_dble(i) ! [inches]
@@ -789,6 +793,8 @@ INTEGER FUNCTION snoinit()
               Frac_swe(i) = MIN( 1.0, Frac_swe(i) )
             ENDIF
             CALL sca_deplcrv(Snowcov_area(i), Snarea_curve(:,Hru_deplcrv(i)), Frac_swe(i))
+            Pk_depth(i) = Pkwater_equiv(i)/DBLE( Den_init(i) )
+            Pk_den = SNGL( Pkwater_equiv(i)/Pk_depth(i) )
           ENDIF
         ENDDO
         DEALLOCATE ( Snowpack_init )

GSFLOW/src/prms/subbasin.f90

@@ -30,7 +30,7 @@ MODULE PRMS_SUBBASIN
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_rain(:), Subinc_snow(:), Subinc_stor(:)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_recharge(:), Subinc_szstor_frac(:), Subinc_capstor_frac(:)
 !   Declared Parameters
-      INTEGER, SAVE, ALLOCATABLE :: Subbasin_down(:), Hru_subbasin(:)
+      INTEGER, SAVE, ALLOCATABLE :: Subbasin_down(:)
       END MODULE PRMS_SUBBASIN
 
 !***********************************************************************
@@ -65,7 +65,7 @@ INTEGER FUNCTION subdecl()
       USE PRMS_CONSTANTS, ONLY: ACTIVE, ERROR_dim
       use PRMS_MMFAPI, only: declvar_dble
       use PRMS_READ_PARAM_FILE, only: declparam
-      USE PRMS_MODULE, ONLY: Nsub, Nhru, gwflow_flag
+      USE PRMS_MODULE, ONLY: Nsub, gwflow_flag
       USE PRMS_SUBBASIN
       use prms_utils, only: error_stop, print_module, read_error
       IMPLICIT NONE
@@ -220,13 +220,6 @@ INTEGER FUNCTION subdecl()
      &     'Index number for the downstream subbasin whose inflow is outflow from this subbasin', &
      &     'none')/=0 ) CALL read_error(1, 'subbasin_down')
 
-      ALLOCATE ( Hru_subbasin(Nhru) )
-      IF ( declparam(MODNAME, 'hru_subbasin', 'nhru', 'integer', &
-     &     '0', 'bounded', 'nsub', &
-     &     'Index of subbasin assigned to each HRU', &
-     &     'Index of subbasin assigned to each HRU', &
-     &     'none')/=0 ) CALL read_error(1, 'hru_subbasin')
-
 ! Allocate arrays for variables
       ALLOCATE ( Sub_area(Nsub), Laststor(Nsub) )
 
@@ -238,12 +231,12 @@ END FUNCTION subdecl
 !***********************************************************************
       INTEGER FUNCTION subinit()
       USE PRMS_CONSTANTS, ONLY: ACTIVE, CFS2CMS_CONV, LAKE, DNEARZERO
-      USE PRMS_MODULE, ONLY: Nsub, Nhru, Print_debug, &
+      USE PRMS_MODULE, ONLY: Nsub, Print_debug, &
      &    Inputerror_flag, Dprst_flag, Lake_route_flag, Cascade_flag, gwflow_flag, Hru_type
       use PRMS_READ_PARAM_FILE, only: getparam_int
       USE PRMS_SUBBASIN
       USE PRMS_BASIN, ONLY: Hru_area_dble, Active_hrus, Hru_route_order, &
-     &    Hru_frac_perv, Lake_hru_id
+     &    Hru_frac_perv, Lake_hru_id, Hru_subbasin
       USE PRMS_FLOWVARS, ONLY: Ssres_stor, Soil_moist, Pkwater_equiv, Gwres_stor, Sroff, Ssres_flow, Lake_vol, &
      &    Hru_impervstor, Dprst_stor_hru, Hru_intcpstor
       USE PRMS_SET_TIME, ONLY: Cfs_conv, Cfs2inches
@@ -261,7 +254,6 @@ INTEGER FUNCTION subinit()
 !***********************************************************************
       subinit = 0
 
-      IF ( getparam_int(MODNAME, 'hru_subbasin', Nhru, Hru_subbasin)/=0 ) CALL read_error(2, 'hru_subbasin')
       IF ( getparam_int(MODNAME, 'subbasin_down', Nsub, Subbasin_down)/=0 ) CALL read_error(2, 'subbasin_down')
 
 ! Determine the tree structure for the internal nodes
@@ -413,7 +405,7 @@ INTEGER FUNCTION subrun()
       USE PRMS_MODULE, ONLY: Nsub, Dprst_flag, Lake_route_flag, Cascade_flag, Hru_type, gwflow_flag
       USE PRMS_SUBBASIN
       USE PRMS_BASIN, ONLY: Hru_area_dble, Active_hrus, Hru_route_order, &
-     &    Hru_frac_perv, Lake_hru_id
+     &    Hru_frac_perv, Lake_hru_id, Hru_subbasin
       USE PRMS_SET_TIME, ONLY: Cfs_conv, Cfs2inches
       USE PRMS_CLIMATEVARS, ONLY: Hru_ppt, Swrad, Potet, Tminc, Tmaxc, Tavgc, Hru_rain, Hru_snow
       USE PRMS_FLOWVARS, ONLY: Hru_actet, Ssres_flow, Sroff, &