minor updates based on comparison with PRMS 5.2.1.1 and 6.0.0

GSFLOW/src/gsflow/gsflow_module.f90

@@ -9,8 +9,8 @@ MODULE PRMS_MODULE
       character(len=*), parameter :: MODDESC = 'PRMS Computation Order'
       character(len=11), parameter :: MODNAME = 'gsflow_prms'
       character(len=*), parameter :: GSFLOW_versn = '2.4.0 12/11/2024'
-      character(len=*), parameter :: PRMS_versn = '2024-12-11'
-      character(len=*), parameter :: PRMS_VERSION = 'Version 6.0.0 12/01/2024'
+      character(len=*), parameter :: PRMS_versn = '2024-12-20'
+      character(len=*), parameter :: PRMS_VERSION = 'Version 6.0.0 12/20/2024'
       character(len=*), parameter :: githash = 'Github Commit Hash a4cffeceecab925507a192e0f4822c89f8f37065'
       character(len=*), parameter :: Version_read_control_file = '2024-08-01'
       character(len=*), parameter :: Version_read_parameter_file = '2024-11-25'

GSFLOW/src/gsflow/gsflow_prms.f90

@@ -51,7 +51,7 @@ SUBROUTINE gsflow_prms(Process_mode, AFR, MS_GSF_converge, Nsegshold, Nlakeshold
       INTEGER, EXTERNAL :: stream_temp, glacr, dynamic_soil_param_read, strmflow_character
       INTEGER, EXTERNAL :: soilzone_ag
       EXTERNAL :: precip_map, temp_map, segment_to_hru, gwflow_inactive_cell
-      EXTERNAL :: water_balance, prms_summary, convert_params
+      EXTERNAL :: water_balance, prms_summary, convert_params, input_error
       EXTERNAL :: gsflow_prms2modsim, gsflow_modsim2prms
       INTEGER, EXTERNAL :: gsflow_prms2mf, gsflow_mf2prms, gsflow_budget, gsflow_sum
 ! Local Variables
@@ -375,6 +375,7 @@ SUBROUTINE gsflow_prms(Process_mode, AFR, MS_GSF_converge, Nsegshold, Nlakeshold
       IF ( Model==FROST ) THEN
         IF ( Process_flag==DECL ) CALL read_parameter_file_params()
         ierr = frost_date()
+        IF ( Inputerror_flag == 1 ) CALL input_error()
         RETURN
       ENDIF
 
@@ -390,6 +391,7 @@ SUBROUTINE gsflow_prms(Process_mode, AFR, MS_GSF_converge, Nsegshold, Nlakeshold
 
       IF ( Model==CLIMATE ) THEN
         IF ( Process_flag==DECL ) CALL read_parameter_file_params()
+        IF ( Inputerror_flag == 1 ) CALL input_error()
         CALL summary_output()
         RETURN
       ENDIF
@@ -410,6 +412,7 @@ SUBROUTINE gsflow_prms(Process_mode, AFR, MS_GSF_converge, Nsegshold, Nlakeshold
 
       IF ( Model==TRANSPIRE ) THEN
         IF ( Process_flag==DECL ) CALL read_parameter_file_params()
+        IF ( Inputerror_flag == 1 ) CALL input_error()
         CALL summary_output()
         RETURN
       ENDIF
@@ -435,11 +438,13 @@ SUBROUTINE gsflow_prms(Process_mode, AFR, MS_GSF_converge, Nsegshold, Nlakeshold
       IF ( Model==WRITE_CLIMATE ) THEN
         ierr = write_climate_hru()
         IF ( Process_flag==DECL ) CALL read_parameter_file_params()
+        IF ( Inputerror_flag == 1 ) CALL input_error()
         RETURN
       ENDIF
 
       IF ( Model==POTET ) THEN
         IF ( Process_flag==DECL ) CALL read_parameter_file_params()
+        IF ( Inputerror_flag == 1 ) CALL input_error()
         CALL summary_output()
         RETURN
       ENDIF
@@ -590,19 +595,8 @@ SUBROUTINE gsflow_prms(Process_mode, AFR, MS_GSF_converge, Nsegshold, Nlakeshold
         IF ( Model==CONVERT ) CALL convert_params()
       ELSEIF ( Process_flag==INIT ) THEN
         CALL check_parameters()
-        IF ( Inputerror_flag==1 ) THEN
-          PRINT '(//,A,//,A,/,A,/,A)', '**Fix input errors in your Parameter File to continue**', &
-     &          '  Set control parameter parameter_check_flag to 0 after', &
-     &          '  all parameter values are valid.'
-          PRINT '(/,A,/,A,/,A,/,A,/,A,/)', &
-     &          'If input errors are related to parameters used for automated', &
-     &          'calibration processes, with CAUTION, set control parameter', &
-     &          'parameter_check_flag to 0. After calibration set the', &
-     &          'parameter_check_flag to 1 to verify that those calibration', &
-     &          'parameters have valid and compatible values.'
-        ENDIF
         IF ( Parameter_check_flag==2 ) STOP
-        IF ( Inputerror_flag==1 ) ERROR STOP ERROR_param
+        IF ( Inputerror_flag==1 ) CALL input_error()
         IF ( Model==CONVERT ) THEN
           CALL convert_params()
           STOP
@@ -1517,6 +1511,25 @@ SUBROUTINE check_module_names()
 
       END SUBROUTINE check_module_names
 
+!***********************************************************************
+      SUBROUTINE input_error()
+!***********************************************************************
+      USE PRMS_CONSTANTS, ONLY: ERROR_param
+      IMPLICIT NONE
+!***********************************************************************
+      PRINT '(//,A,//,A,/,A,/,A)', '**Fix input errors in your Parameter File to continue**', &
+            '  Set control parameter parameter_check_flag to 0 after', &
+            '  all parameter values are valid.'
+      PRINT '(/,A,/,A,/,A,/,A,/,A,/)', &
+            'If input errors are related to parameters used for automated', &
+            'calibration processes, with CAUTION, set control parameter', &
+            'parameter_check_flag to 0. After calibration set the', &
+            'parameter_check_flag to 1 to verify that those calibration', &
+            'parameters have valid and compatible values.'
+      ERROR STOP ERROR_param
+
+      END SUBROUTINE input_error
+
 !***********************************************************************
 !     gsflow_prmsSettings - set MODSIM variables set in PRMS
 !***********************************************************************

GSFLOW/src/gsflow/gwflow_inactive_cell.f90

@@ -17,7 +17,7 @@ MODULE PRMS_GWFLOW_INACTIVE_CELL
 !   Local Variables
       character(len=*), parameter :: MODDESC = 'Groundwater'
       character(len=*), parameter :: MODNAME = 'gwflow_inactive_cell'
-      character(len=*), parameter :: Version_gwflow = '2024-12-02'
+      character(len=*), parameter :: Version_gwflow = '2024-12-01'
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gwstor_minarea(:), Gwin_dprst(:), It0_gwres_stor(:)
       DOUBLE PRECISION, SAVE :: Basin_gw_upslope
       INTEGER, SAVE :: Gwminarea_flag
@@ -206,6 +206,7 @@ SUBROUTINE gwflow_inactivecell_init()
       Gwminarea_flag = OFF
       Gwstor_minarea = 0.0D0
       Gwstor_minarea_wb = 0.0D0
+      Basin_gwstor_minarea_wb = 0.0D0
       IF ( Init_vars_from_file==0 .OR. Init_vars_from_file==2 .OR. Init_vars_from_file==6 ) THEN
         IF ( getparam_real(MODNAME, 'gwstor_init', Ngw, Gwstor_init)/=0 ) CALL read_error(2, 'gwstor_init')
         Gwres_stor = DBLE( Gwstor_init )
@@ -249,6 +250,7 @@ SUBROUTINE gwflow_inactivecell_init()
         Hru_storage(i) = Hru_storage(i) + Gwres_stor(i)
       ENDDO
       IF ( Gwminarea_flag==OFF ) DEALLOCATE ( Gwstor_minarea )
+      DEALLOCATE ( Gwstor_min )
       Basin_gwstor = Basin_gwstor*Basin_area_inv
 
       IF ( Dprst_flag==ACTIVE ) Gwin_dprst = 0.0D0
@@ -311,6 +313,7 @@ SUBROUTINE gwflow_inactivecell_run()
       Gw_upslope_to_MF = 0.0
       Gwres_flow = 0.0
       Gwres_sink = 0.0
+      IF ( Gwminarea_flag==ACTIVE ) Gwstor_minarea_wb = 0.0D0
       DO j = 1, Active_gwrs
         i = Gwr_route_order(j)
         IF ( activeHRU_inactiveCell(i) == OFF ) THEN
@@ -345,7 +348,7 @@ SUBROUTINE gwflow_inactivecell_run()
           IF ( gwstor<Gwstor_minarea(i) ) THEN
             IF ( gwstor<0.0D0 ) THEN
               IF ( Print_debug>DEBUG_less ) PRINT *, 'Warning, groundwater reservoir for HRU:', i, &
-                                                     ' is < 0.0 with gwstor_min active', gwstor
+     &                                               ' is < 0.0 with gwstor_min active', gwstor
 !              ERROR STOP ERROR_var
             ENDIF
             gwstor_last = gwstor
@@ -355,7 +358,7 @@ SUBROUTINE gwflow_inactivecell_run()
             Basin_gwstor_minarea_wb = Basin_gwstor_minarea_wb + Gwstor_minarea_wb(i)
             Gwstor_minarea_wb(i) = Gwstor_minarea_wb(i)/gwarea
             IF ( Print_debug>DEBUG_less ) PRINT *, 'Added to gwres_stor as storage < gwstor_min to GWR:', i, &
-                                                   ' amount:', Gwstor_minarea_wb(i)
+     &                                             ' amount:', Gwstor_minarea_wb(i)
           ELSE
             Gwstor_minarea_wb(i) = 0.0D0
           ENDIF
@@ -365,7 +368,7 @@ SUBROUTINE gwflow_inactivecell_run()
           IF ( Gwr_transfer(i)>0.0 ) THEN
             IF ( SNGL(gwstor*Cfs_conv)<Gwr_transfer(i) ) THEN
               PRINT *, 'ERROR, not enough storage for transfer from groundwater reservoir storage:', &
-                       i, ' Date:', Nowyear, Nowmonth, Nowday
+     &                  i, ' Date:', Nowyear, Nowmonth, Nowday
               PRINT *, '       storage: ', gwstor, '; transfer: ', Gwr_transfer(i)/Cfs_conv
               ERROR STOP ERROR_water_use
             ENDIF
@@ -374,22 +377,19 @@ SUBROUTINE gwflow_inactivecell_run()
         ENDIF
 
         gwsink = 0.0D0
-        IF ( gwstor<0.0D0 ) THEN ! could happen with water use
-          IF ( Print_debug>DEBUG_less ) PRINT *, 'Warning, groundwater reservoir for HRU:', i, ' is < 0.0, set to 0.0', gwstor
-          gwflow = 0.0D0
-          gwstor = 0.0D0
-        ELSEIF ( gwstor>0.0D0 ) THEN
-          ! Compute groundwater discharge
+        IF ( gwstor>0.0D0 ) THEN
+! Compute groundwater discharge
           gwflow = gwstor*DBLE( Gwflow_coef(i) )
-          ! Reduce storage by outflow
+
+! Reduce storage by outflow
           gwstor = gwstor - gwflow
 
           IF ( Gwsink_coef(i)>0.0 ) THEN
             gwsink = MIN( gwstor*DBLE( Gwsink_coef(i) ), gwstor ) ! if gwsink_coef > 1, could have had negative gwstor
             gwstor = gwstor - gwsink
           ENDIF
-          ! if gwr_swale_flag = 1 swale GWR flow goes to sink, 2 included in stream network and cascades
-          ! maybe gwr_swale_flag = 3 abs(hru_segment) so hru_segment could be changed from 0 to allow HRU swales
+! if gwr_swale_flag = 1 swale GWR flow goes to sink, 2 included in stream network and cascades
+! maybe gwr_swale_flag = 3 abs(hru_segment) so hru_segment could be changed from 0 to allow HRU swales
           IF ( Gwr_swale_flag==ACTIVE ) THEN
             IF ( Gwr_type(i)==SWALE ) THEN
               gwsink = gwsink + gwflow
@@ -400,6 +400,9 @@ SUBROUTINE gwflow_inactivecell_run()
           Basin_gwsink = Basin_gwsink + gwsink
           Basin_gwstor = Basin_gwstor + gwstor
           Gwres_flow(i) = SNGL( gwflow/gwarea )
+        ELSEIF ( gwstor<0.0D0 ) THEN ! could happen with water use
+          IF ( Print_debug>DEBUG_less ) PRINT *, 'Warning, groundwater reservoir for HRU:', i, ' is < 0.0, set to 0.0', gwstor
+          gwstor = 0.0D0
         ENDIF
 
         IF ( Cascadegw_flag>CASCADEGW_OFF ) THEN

GSFLOW/src/prms/basin_sum.f90

@@ -7,7 +7,7 @@ MODULE PRMS_BASINSUM
 !   Local Variables
       character(len=*), parameter :: MODDESC = 'Output Summary'
       character(len=9), parameter :: MODNAME = 'basin_sum'
-      character(len=*), parameter :: Version_basin_sum = '2024-09-01'
+      character(len=*), parameter :: Version_basin_sum = '2024-12-01'
 
       INTEGER, SAVE :: BALUNT, Totdays
       INTEGER, SAVE :: Header_prt, Endjday
@@ -311,9 +311,9 @@ END FUNCTION sumbdecl
 !     sumbinit - Initialize basinsum module - get parameter values
 !***********************************************************************
       INTEGER FUNCTION sumbinit()
-      USE PRMS_CONSTANTS, ONLY: OFF, ERROR_param
+      USE PRMS_CONSTANTS, ONLY: OFF, ACTIVE
       use PRMS_READ_PARAM_FILE, only: getparam_int
-      USE PRMS_MODULE, ONLY: Nobs, Init_vars_from_file, Print_debug, Inputerror_flag
+      USE PRMS_MODULE, ONLY: Nobs, Init_vars_from_file, Print_debug, Inputerror_flag, Stream_order_flag
       USE PRMS_BASINSUM
       USE PRMS_FLOWVARS, ONLY: Basin_soil_moist, Basin_ssstor, Basin_lake_stor, Basin_pweqv
       USE PRMS_INTCP, ONLY: Basin_intcp_stor
@@ -451,7 +451,8 @@ INTEGER FUNCTION sumbinit()
       Basin_storage = Basin_soil_moist + Basin_intcp_stor + &
      &                Basin_gwstor + Basin_ssstor + Basin_pweqv + &
      &                Basin_imperv_stor + Basin_lake_stor + &
-     &                Basin_dprst_volop + Basin_dprst_volcl + Basin_segment_storage
+     &                Basin_dprst_volop + Basin_dprst_volcl
+      IF ( Stream_order_flag == ACTIVE ) Basin_storage = Basin_storage + Basin_segment_storage
 !glacier storage not known at start
 
       IF ( Print_freq/=0 ) THEN
@@ -481,7 +482,8 @@ END FUNCTION sumbinit
 !***********************************************************************
       INTEGER FUNCTION sumbrun()
       USE PRMS_CONSTANTS, ONLY: ACTIVE, strmflow_muskingum_lake_module
-      USE PRMS_MODULE, ONLY: Nobs, Print_debug, End_year, Strmflow_flag, Glacier_flag, Nowyear, Nowmonth, Nowday, Nratetbl
+      USE PRMS_MODULE, ONLY: Nobs, Print_debug, End_year, Strmflow_flag, Glacier_flag, Nowyear, Nowmonth, Nowday, Nratetbl, &
+          Stream_order_flag
       USE PRMS_BASINSUM
       USE PRMS_BASIN, ONLY: Active_area, Active_hrus, Hru_route_order
       USE PRMS_FLOWVARS, ONLY: Basin_ssflow, Basin_lakeevap, &
@@ -522,11 +524,12 @@ INTEGER FUNCTION sumbrun()
       Last_basin_stor = Basin_storage
       Basin_storage = Basin_soil_moist + Basin_intcp_stor + &
      &                Basin_gwstor + Basin_ssstor + Basin_pweqv + &
-     &                Basin_imperv_stor + Basin_lake_stor + Basin_dprst_volop + Basin_dprst_volcl + Basin_segment_storage
+     &                Basin_imperv_stor + Basin_lake_stor + Basin_dprst_volop + Basin_dprst_volcl
 ! Basin_storage doesn't include any processes on glacier
 ! In glacier module, Basin_gl_storstart is an estimate for starting glacier volume, but only
 !   includes glaciers that have depth estimates and these are known to be iffy
       IF ( Glacier_flag==ACTIVE ) Basin_storage = Basin_storage + Basin_gl_storage
+      IF ( Stream_order_flag == ACTIVE )  Basin_storage = Basin_storage + Basin_segment_storage
 
 ! volume calculation for storage
       Basin_storvol = Basin_storage*Active_area

GSFLOW/src/prms/glacr_melt.f90

@@ -94,7 +94,7 @@ MODULE PRMS_GLACR
       REAL, SAVE :: Max_gldepth
       REAL, SAVE, ALLOCATABLE :: Glacrva_coef(:), Glacrva_exp(:), Hru_length(:), Hru_width(:)
       REAL, SAVE, ALLOCATABLE :: Stor_ice(:,:), Stor_snow(:,:), Stor_firn(:,:)
-      REAL, SAVE, ALLOCATABLE :: Abl_elev_range(:)
+      REAL, SAVE, ALLOCATABLE :: Hru_slope(:), Abl_elev_range(:)
 
       END MODULE PRMS_GLACR
 
@@ -371,6 +371,13 @@ INTEGER FUNCTION glacrdecl()
      &     ' glacier melt flows, for non-glacier HRUs that do not flow to another HRU enter 0', &
      &     'none')/=0 ) CALL read_error(1, 'tohru')
 
+      ALLOCATE ( Hru_slope(Nhru) )
+      IF ( declparam(MODNAME, 'hru_slope', 'nhru', 'real', &
+     &     '0.0', '0.0', '10.0', &
+     &     'HRU slope', &
+     &     'Slope of each HRU, specified as change in vertical length divided by change in horizontal length', &
+     &     'decimal fraction')/=0 ) CALL read_error(1, 'hru_slope')
+
       IF ( declparam(MODNAME, 'max_gldepth', 'one', 'real', &
            '1.5', '0.1', '3.0', &
            'Upper bound on glacier thickness, thickest glacier measured is Taku at 1.5 km, ice sheet 3 km', &
@@ -446,7 +453,6 @@ INTEGER FUNCTION glacrinit()
       USE PRMS_BASIN, ONLY: Hru_area_dble, Hru_elev_ts, Active_hrus, Hru_route_order, &
      &    Basin_area_inv, Hru_elev_meters
       USE PRMS_FLOWVARS, ONLY: Glacier_frac, Alt_above_ela, Glrette_frac
-      USE PRMS_SOLTAB, ONLY: Hru_slope
       use prms_utils, only: get_ftnunit, read_error
       IMPLICIT NONE
 ! Functions
@@ -479,6 +485,7 @@ INTEGER FUNCTION glacrinit()
       IF ( getparam_real(MODNAME, 'hru_width', Nhru, Hru_width)/=0 ) CALL read_error(2, 'hru_width')
       IF ( getparam_real(MODNAME, 'abl_elev_range', Nhru, Abl_elev_range)/=0 ) CALL read_error(2, 'abl_elev_range')
       IF ( getparam_int(MODNAME, 'tohru', Nhru, Tohru)/=0 ) CALL read_error(2, 'tohru')
+      IF ( getparam_real(MODNAME, 'hru_slope', Nhru, Hru_slope)/=0 ) CALL read_error(2, 'hru_slope')
       IF ( Init_vars_from_file==0 ) THEN
         Alt_above_ela = 0.0
         Prev_out = 0.0
@@ -516,6 +523,7 @@ INTEGER FUNCTION glacrinit()
         Basin_gl_storstart = 0.0D0
         Basin_gl_storvol = 0.0D0
       ENDIF
+      DEALLOCATE ( Hru_slope )
 
       Glac_HRUnum_down = 1 ! 1 is the way Weasel delineation was designed
       ! 1 is terminus is smallest ID and top is largest. IDs are stacked.