From 9ed9e3588f9e0d2b1d128d1bd112988d75e6abfd Mon Sep 17 00:00:00 2001 From: Steve Regan Date: Wed, 18 Dec 2024 16:14:12 -0700 Subject: [PATCH] minor updates based on comparison with PRMS 5.2.1.1 and 6.0.0 --- GSFLOW/src/gsflow/gsflow_module.f90 | 4 +- GSFLOW/src/gsflow/gsflow_prms.f90 | 39 ++++++++----- GSFLOW/src/gsflow/gwflow_inactive_cell.f90 | 29 +++++----- GSFLOW/src/prms/basin_sum.f90 | 15 +++-- GSFLOW/src/prms/glacr_melt.f90 | 12 +++- GSFLOW/src/prms/gwflow.f90 | 43 ++++++++------ GSFLOW/src/prms/intcp.f90 | 6 +- GSFLOW/src/prms/muskingum.f90 | 15 ++++- GSFLOW/src/prms/nsegment_summary.f90 | 2 +- GSFLOW/src/prms/routing.f90 | 6 +- GSFLOW/src/prms/snowcomp.f90 | 3 +- GSFLOW/src/prms/soilzone.f90 | 13 +++-- GSFLOW/src/prms/soltab.f90 | 3 +- GSFLOW/src/prms/subbasin.f90 | 67 +++------------------- GSFLOW/src/prms/water_balance.f90 | 10 ++-- 15 files changed, 130 insertions(+), 137 deletions(-) diff --git a/GSFLOW/src/gsflow/gsflow_module.f90 b/GSFLOW/src/gsflow/gsflow_module.f90 index db9d1ea3..b9d35705 100644 --- a/GSFLOW/src/gsflow/gsflow_module.f90 +++ b/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' diff --git a/GSFLOW/src/gsflow/gsflow_prms.f90 b/GSFLOW/src/gsflow/gsflow_prms.f90 index a165bb8b..2577be1e 100644 --- a/GSFLOW/src/gsflow/gsflow_prms.f90 +++ b/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 !*********************************************************************** diff --git a/GSFLOW/src/gsflow/gwflow_inactive_cell.f90 b/GSFLOW/src/gsflow/gwflow_inactive_cell.f90 index c4e2b46d..e0f1a1f6 100644 --- a/GSFLOW/src/gsflow/gwflow_inactive_cell.f90 +++ b/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 ( gwstorDEBUG_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)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 diff --git a/GSFLOW/src/prms/basin_sum.f90 b/GSFLOW/src/prms/basin_sum.f90 index 1a1de859..995f14cb 100644 --- a/GSFLOW/src/prms/basin_sum.f90 +++ b/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 diff --git a/GSFLOW/src/prms/glacr_melt.f90 b/GSFLOW/src/prms/glacr_melt.f90 index d5ff4174..a16427d1 100644 --- a/GSFLOW/src/prms/glacr_melt.f90 +++ b/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. diff --git a/GSFLOW/src/prms/gwflow.f90 b/GSFLOW/src/prms/gwflow.f90 index 97e2efe6..b103c4b5 100644 --- a/GSFLOW/src/prms/gwflow.f90 +++ b/GSFLOW/src/prms/gwflow.f90 @@ -16,8 +16,8 @@ MODULE PRMS_GWFLOW IMPLICIT NONE ! Local Variables character(len=*), parameter :: MODDESC = 'Groundwater' - character(len=*), parameter :: MODNAME = 'gwflow' - character(len=*), parameter :: Version_gwflow = '2024-12-02' + character(len=6), parameter :: MODNAME = 'gwflow' + character(len=*), parameter :: Version_gwflow = '2024-12-01' DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gwstor_minarea(:), Gwin_dprst(:) DOUBLE PRECISION, SAVE :: Basin_gw_upslope INTEGER, SAVE :: Gwminarea_flag @@ -48,13 +48,13 @@ INTEGER FUNCTION gwflow() USE PRMS_MODULE, ONLY: Process_flag, Init_vars_from_file, Save_vars_to_file IMPLICIT NONE ! Functions - INTEGER, EXTERNAL :: gwflowdecl, gwflowinit - EXTERNAL :: gwflow_restart, gwflowrun + INTEGER, EXTERNAL :: gwflowdecl, gwflowinit, gwflowrun + EXTERNAL :: gwflow_restart !*********************************************************************** gwflow = 0 IF ( Process_flag==RUN ) THEN - CALL gwflowrun() + gwflow = gwflowrun() ELSEIF ( Process_flag==DECL ) THEN gwflow = gwflowdecl() ELSEIF ( Process_flag==INIT ) THEN @@ -282,6 +282,7 @@ INTEGER FUNCTION gwflowinit() 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 ) @@ -319,6 +320,7 @@ INTEGER FUNCTION gwflowinit() 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 @@ -363,7 +365,7 @@ END FUNCTION gwflowinit ! gwflowrun - Computes groundwater flow to streamflow and to ! groundwater sink !*********************************************************************** - SUBROUTINE gwflowrun() + INTEGER FUNCTION gwflowrun() USE PRMS_CONSTANTS, ONLY: ACTIVE, LAKE, SWALE, DEBUG_less, CASCADEGW_OFF, ERROR_water_use USE PRMS_MODULE, ONLY: Nlake, Print_debug, Dprst_flag, Cascadegw_flag, Gwr_swale_flag, & & Gwr_add_water_use, Gwr_transfer_water_use, Nowyear, Nowmonth, Nowday @@ -386,6 +388,8 @@ SUBROUTINE gwflowrun() REAL :: dnflow DOUBLE PRECISION :: gwin, gwstor, gwsink, gwflow, gwstor_last, seepage, gwarea, inch2acre_feet !*********************************************************************** + gwflowrun = 0 + IF ( Cascadegw_flag>CASCADEGW_OFF ) THEN Gw_upslope = 0.0D0 Hru_gw_cascadeflow = 0.0D0 @@ -454,6 +458,7 @@ SUBROUTINE gwflowrun() Basin_gwin = 0.0D0 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) gwarea = Hru_area_dble(i) @@ -481,7 +486,7 @@ SUBROUTINE gwflowrun() IF ( gwstorDEBUG_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 @@ -491,7 +496,7 @@ SUBROUTINE gwflowrun() 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 @@ -501,7 +506,7 @@ SUBROUTINE gwflowrun() IF ( Gwr_transfer(i)>0.0 ) THEN IF ( SNGL(gwstor*Cfs_conv)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 @@ -536,6 +538,9 @@ SUBROUTINE gwflowrun() 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 @@ -568,7 +573,7 @@ SUBROUTINE gwflowrun() Basin_gwstor_minarea_wb = Basin_gwstor_minarea_wb*Basin_area_inv Basin_dnflow = Basin_dnflow*Basin_area_inv - END SUBROUTINE gwflowrun + END FUNCTION gwflowrun !*********************************************************************** ! Compute cascading GW flow diff --git a/GSFLOW/src/prms/intcp.f90 b/GSFLOW/src/prms/intcp.f90 index 91f68039..4a621c45 100644 --- a/GSFLOW/src/prms/intcp.f90 +++ b/GSFLOW/src/prms/intcp.f90 @@ -258,13 +258,13 @@ END FUNCTION intinit ! and evaporation for each HRU !*********************************************************************** INTEGER FUNCTION intrun() - USE PRMS_CONSTANTS, ONLY: ACTIVE, OFF, DEBUG_WB, NEARZERO, ZERO_SNOWPACK, & + USE PRMS_CONSTANTS, ONLY: ACTIVE, OFF, DEBUG_WB, NEARZERO, & & DEBUG_less, LAKE, BARESOIL, GRASSES, ERROR_param, CANOPY USE PRMS_MODULE, ONLY: Print_debug, PRMS_land_iteration_flag, Kkiter, Nowyear, Nowmonth, Nowday, & & Ag_package, Hru_ag_irr, irrigation_apply_flag, Hru_type USE PRMS_INTCP USE PRMS_BASIN, ONLY: Basin_area_inv, Active_hrus, Covden_win, Covden_sum, & - & Hru_route_order, Hru_area, Cov_type, Ag_frac, gsflow_ag_area !, Ag_cov_type + & Hru_route_order, Hru_area, Cov_type, Snowpack_threshold, Ag_frac, gsflow_ag_area !, Ag_cov_type USE PRMS_WATER_USE, ONLY: Canopy_gain ! need to add ag apply ??? ! Newsnow and Pptmix can be modfied, WARNING!!! USE PRMS_CLIMATEVARS, ONLY: Newsnow, Pptmix, Hru_rain, Hru_ppt, & @@ -407,7 +407,7 @@ INTEGER FUNCTION intrun() ELSEIF ( Cov_type(i)==GRASSES ) THEN ! cov_type = 1 !rsr, 03/24/2008 intercept rain on snow-free grass, !rsr when not a mixed event - IF ( Pkwater_equiv(i)0.0 ) THEN Slow_stor(i) = MIN( Ssres_stor(i), Pref_flow_thrsh(i) ) Pref_flow_stor(i) = Ssres_stor(i) - Slow_stor(i) + ELSE + Pref_flow_stor(i) = 0.0 ENDIF Ssres_stor(i) = Slow_stor(i) + Pref_flow_stor(i) ENDIF @@ -670,6 +672,7 @@ INTEGER FUNCTION szinit() Soil_zone_max(i) = Sat_threshold(i) + Soil_moist_max(i)*Hru_frac_perv(i) Soil_moist_tot(i) = Ssres_stor(i) + Soil_moist(i)*Hru_frac_perv(i) + Soil_lower(i) = Soil_moist(i) - Soil_rechr(i) Hru_storage(i) = DBLE( Soil_moist_tot(i) + Hru_intcpstor(i) + Hru_impervstor(i) ) + Pkwater_equiv(i) IF ( Dprst_flag==ACTIVE ) Hru_storage(i) = Hru_storage(i) + Dprst_stor_hru(i) ENDDO diff --git a/GSFLOW/src/prms/soltab.f90 b/GSFLOW/src/prms/soltab.f90 index e9b3ee76..c05e7ba1 100644 --- a/GSFLOW/src/prms/soltab.f90 +++ b/GSFLOW/src/prms/soltab.f90 @@ -267,7 +267,8 @@ INTEGER FUNCTION sthinit() ! data jday/356,10,23,38,51,66,80,94,109,123,138,152,173/ ENDIF - IF ( Glacier_flag==OFF ) DEALLOCATE ( Hru_aspect, Hru_slope, Hru_lat ) + DEALLOCATE ( Hru_slope ) + IF ( Glacier_flag==OFF ) DEALLOCATE ( Hru_aspect ) END FUNCTION sthinit diff --git a/GSFLOW/src/prms/subbasin.f90 b/GSFLOW/src/prms/subbasin.f90 index 31d9f3a2..b60903cc 100644 --- a/GSFLOW/src/prms/subbasin.f90 +++ b/GSFLOW/src/prms/subbasin.f90 @@ -15,7 +15,7 @@ MODULE PRMS_SUBBASIN ! Local Variables character(len=*), parameter :: MODDESC = 'Output Summary' character(len=*), parameter :: MODNAME = 'subbasin' - character(len=*), parameter :: Version_subbasin = '2023-11-01' + character(len=*), parameter :: Version_subbasin = '2024-12-01' DOUBLE PRECISION, SAVE, ALLOCATABLE :: Qsub(:), Sub_area(:), Laststor(:) INTEGER, SAVE, ALLOCATABLE :: Tree(:, :) ! Declared Variables @@ -232,25 +232,20 @@ END FUNCTION subdecl INTEGER FUNCTION subinit() USE PRMS_CONSTANTS, ONLY: ACTIVE, CFS2CMS_CONV, LAKE, DNEARZERO USE PRMS_MODULE, ONLY: Nsub, Print_debug, & - & Inputerror_flag, Dprst_flag, Lake_route_flag, Cascade_flag, gwflow_flag, Hru_type + & Inputerror_flag, Dprst_flag, Lake_route_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_subbasin - USE PRMS_FLOWVARS, ONLY: Ssres_stor, Soil_moist, Pkwater_equiv, Gwres_stor, Sroff, Ssres_flow, Lake_vol, & + USE PRMS_FLOWVARS, ONLY: Ssres_stor, Soil_moist, Pkwater_equiv, Gwres_stor, Lake_vol, & & Hru_impervstor, Dprst_stor_hru, Hru_intcpstor - USE PRMS_SET_TIME, ONLY: Cfs_conv, Cfs2inches - USE PRMS_SRUNOFF, ONLY: Hortonian_lakes - USE PRMS_SOILZONE, ONLY: Lakein_sz - USE PRMS_GWFLOW, ONLY: Gwres_flow - USE PRMS_MUSKINGUM_LAKE, ONLY: Lake_outcfs use prms_utils, only: PRMS_open_module_file, read_error IMPLICIT NONE ! Functions INTRINSIC :: DBLE ! Local Variables INTEGER :: i, j, k, kk, TREEUNIT - DOUBLE PRECISION :: harea, gwstor, soilstor, snowstor, landstor, srq, ssq, gwq + DOUBLE PRECISION :: harea, gwstor, soilstor, snowstor, landstor !*********************************************************************** subinit = 0 @@ -304,28 +299,18 @@ INTEGER FUNCTION subinit() CLOSE ( TREEUNIT ) ENDIF -! added some code to allow for restart, but not climate states and fluxes and subinc_deltastor - Subinc_interflow = 0.0D0 - IF ( gwflow_flag==ACTIVE ) Subinc_gwflow = 0.0D0 - Subinc_sroff = 0.0D0 +! added some code to allow for restart Subinc_stor = 0.0D0 Sub_area = 0.0D0 gwstor = 0.0D0 - gwq = 0.0D0 - Qsub = 0.0D0 DO i = 1, Active_hrus j = Hru_route_order(i) ! k indicates which HRU is in which subbasin k = Hru_subbasin(j) IF ( k>0 ) THEN harea = Hru_area_dble(j) - IF ( gwflow_flag==ACTIVE ) THEN - gwstor = Gwres_stor(j)*harea - gwq = DBLE(Gwres_flow(j))*harea - ENDIF + gwstor = Gwres_stor(j)*harea IF ( Hru_type(j)/=LAKE ) THEN - srq = DBLE(Sroff(j))*harea - ssq = DBLE(Ssres_flow(j))*harea soilstor = DBLE(Soil_moist(j)*Hru_frac_perv(j) + Ssres_stor(j))*harea snowstor = Pkwater_equiv(j)*harea landstor = DBLE(Hru_intcpstor(j)+Hru_impervstor(j))*harea @@ -335,22 +320,8 @@ INTEGER FUNCTION subinit() snowstor = 0.0D0 landstor = 0.0D0 ! wrong if multiple HRUs for any lake - IF ( Lake_route_flag==ACTIVE ) THEN - landstor = Lake_vol(Lake_hru_id(j))*12.0D0 - srq = Lake_outcfs(Lake_hru_id(j))*Cfs2inches - ssq = 0.0D0 - ELSEIF ( Cascade_flag>0 ) THEN - srq = Hortonian_lakes(j)*harea - ssq = Lakein_sz(j)*harea - ELSE - srq = 0.0D0 - ssq = 0.0D0 - ENDIF + IF ( Lake_route_flag==ACTIVE ) landstor = Lake_vol(Lake_hru_id(j))*12.0D0 ENDIF - Qsub(k) = Qsub(k) + srq + ssq + gwq - Subinc_interflow(k) = Subinc_interflow(k) + ssq - Subinc_sroff(k) = Subinc_sroff(k) + srq - IF ( gwflow_flag==ACTIVE ) Subinc_gwflow(k) = Subinc_gwflow(k) + gwq Subinc_stor(k) = Subinc_stor(k) + soilstor + gwstor + snowstor + landstor Sub_area(k) = Sub_area(k) + harea ENDIF @@ -363,31 +334,9 @@ INTEGER FUNCTION subinit() ELSE Subinc_stor(i) = Subinc_stor(i)/Sub_area(i) ENDIF - Sub_inq(i) = Qsub(i)*Cfs_conv - Subinc_interflow(i) = Subinc_interflow(i)*Cfs_conv - IF ( gwflow_flag==ACTIVE ) Subinc_gwflow(i) = Subinc_gwflow(i)*Cfs_conv - Subinc_sroff(i) = Subinc_sroff(i)*Cfs_conv ! water balance off if lake or muskingum routing ENDDO -! allow for possible restart - !get cumulative subbasin flows - DO j = 1, Nsub - IF ( gwflow_flag==ACTIVE ) Sub_gwflow(j) = Subinc_gwflow(j) - Sub_sroff(j) = Subinc_sroff(j) - Sub_interflow(j) = Subinc_interflow(j) - Sub_cfs(j) = Sub_inq(j) - DO k = 1, Nsub - IF ( Tree(j,k)/=0 ) THEN - IF ( gwflow_flag==ACTIVE ) Sub_gwflow(j) = Sub_gwflow(j) + Subinc_gwflow(k) - Sub_sroff(j) = Sub_sroff(j) + Subinc_sroff(k) - Sub_interflow(j) = Sub_interflow(j) + Subinc_interflow(k) - Sub_cfs(j) = Sub_cfs(j) + Sub_inq(k) - ENDIF - ENDDO - Sub_cms(j) = Sub_cfs(j)*CFS2CMS_CONV - ENDDO - 9001 FORMAT ('Initial Tree Structure for Internal Subbasins', /, & & ' 1 indicates a node that flows directly into subbasin', /, 5X, 500I3) 9002 FORMAT (I3, 2X, 500I3) @@ -547,10 +496,10 @@ INTEGER FUNCTION subrun() ENDDO !get cumulative subbasin flows - IF ( gwflow_flag==ACTIVE ) Sub_gwflow = Subinc_gwflow DO j = 1, Nsub Sub_sroff(j) = Subinc_sroff(j) Sub_interflow(j) = Subinc_interflow(j) + IF ( gwflow_flag==ACTIVE ) Sub_gwflow(j) = Subinc_gwflow(j) DO k = 1, Nsub IF ( Tree(j,k)/=0 ) THEN IF ( gwflow_flag==ACTIVE ) Sub_gwflow(j) = Sub_gwflow(j) + Subinc_gwflow(k) diff --git a/GSFLOW/src/prms/water_balance.f90 b/GSFLOW/src/prms/water_balance.f90 index 346c4d2b..1dc0768a 100644 --- a/GSFLOW/src/prms/water_balance.f90 +++ b/GSFLOW/src/prms/water_balance.f90 @@ -6,7 +6,7 @@ MODULE PRMS_WATER_BALANCE ! Local Variables character(len=*), parameter :: MODDESC = 'Water Balance Computations' character(len=*), parameter :: MODNAME_WB = 'water_balance' - character(len=*), parameter :: Version_water_balance = '2024-12-02' + character(len=*), parameter :: Version_water_balance = '2024-12-15' INTEGER, SAVE :: BALUNT, SZUNIT, GWUNIT, INTCPUNT, SROUNIT, SNOWUNIT REAL, PARAMETER :: TOOSMALL = 3.1E-05, SMALL = 1.0E-04, BAD = 1.0E-03 DOUBLE PRECISION, PARAMETER :: DSMALL = 1.0D-04, DTOOSMALL = 1.0D-05 @@ -141,13 +141,13 @@ END SUBROUTINE water_balance_init ! water_balance_run - Computes balance for each HRU and model domain !*********************************************************************** SUBROUTINE water_balance_run() - USE PRMS_CONSTANTS, ONLY: ACTIVE, ZERO_SNOWPACK, LAKE, CASCADE_OFF, CASCADEGW_OFF, OFF + USE PRMS_CONSTANTS, ONLY: ACTIVE, LAKE, CASCADE_OFF, CASCADEGW_OFF, OFF USE PRMS_MODULE, ONLY: Cascade_flag, Cascadegw_flag, Dprst_flag, Glacier_flag, Nowyear, Nowmonth, Nowday, & & Hru_type, AG_flag USE PRMS_WATER_BALANCE USE PRMS_BASIN, ONLY: Hru_route_order, Active_hrus, Hru_frac_perv, Hru_area_dble, Hru_perv, & & Basin_area_inv, Dprst_area_max, Hru_frac_imperv, Hru_frac_dprst, Cov_type, Hru_storage, & - & Covden_win, Covden_sum, Hru_area, Ag_frac + & Covden_win, Covden_sum, Hru_area, Snowpack_threshold, Ag_frac USE PRMS_CLIMATEVARS, ONLY: Hru_ppt, Basin_ppt, Hru_rain, Hru_snow, Pptmix, Potet USE PRMS_FLOWVARS, ONLY: Basin_soil_moist, Basin_ssstor, Soil_to_gw, Soil_to_ssr, & & Infil, Soil_moist_max, Ssr_to_gw, Ssres_flow, Basin_soil_to_gw, Soil_moist, Ssres_stor, Pref_flow_stor, & @@ -237,7 +237,7 @@ SUBROUTINE water_balance_run() ENDIF ! Skip the HRU if there is no snowpack and no new snow - IF ( It0_pkwater_equiv(i)>ZERO_SNOWPACK .OR. Net_snow(i)>0.0 ) THEN + IF ( It0_pkwater_equiv(i)>Snowpack_threshold(i) .OR. Net_snow(i)>0.0 ) THEN hrubal = SNGL( It0_pkwater_equiv(i) - Pkwater_equiv(i) ) - Snow_evap(i) - Snowmelt(i) IF ( Pptmix_nopack(i)==ACTIVE ) THEN hrubal = hrubal + Net_snow(i) @@ -268,7 +268,7 @@ SUBROUTINE water_balance_run() IF ( Pptmix_nopack(i) == ACTIVE ) waterin = waterin + Net_rain(i) IF ( Snowmelt(i)>0.0 ) THEN waterin = waterin + Snowmelt(i) - ELSEIF ( Pkwater_equiv(i)0.0) .AND. Net_rain(i)>0.0 ) THEN IF ( .not.(Pk_precip(i)>0.0) .AND. Pptmix_nopack(i) == OFF ) waterin = waterin + Net_rain(i) ENDIF