diff --git a/.gitignore b/.gitignore index 2fa838c3..73c5abf2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ *.o *.mod *.a -prms/prms +prms/prms \ No newline at end of file diff --git a/Makefile b/Makefile old mode 100644 new mode 100755 index dc311de6..687d4ece --- a/Makefile +++ b/Makefile @@ -3,6 +3,9 @@ # # Top-level makefile for the PRMS # +#------------------------------------------------------------------- +# $Id: Makefile +#------------------------------------------------------------------- include ./makelist @@ -10,15 +13,42 @@ include ./makelist # Standard Targets for Users # -all: standard +all: prmsglrip prmsgl -standard: +prmsglrip: +# Create lib directory, if necessary + @if [ ! -d $(MMFDIR) ] ; then \ + mkdir $(MMFDIR) ; \ + echo Created directory $(MMFDIR) ; \ + fi +# Create bin directory, if necessary + @if [ ! -d $(BINDIR) ] ; then \ + mkdir $(BINDIR) ; \ + echo Created directory $(BINDIR) ; \ + fi cd $(MMFDIR); $(MAKE); - cd $(MIZU); $(MAKE); + cd $(MIZUDIR); $(MAKE); + cd $(PRMSRDIR); $(MAKE); + +prmsgl: +# Create lib directory, if necessary + @if [ ! -d $(MMFDIR) ] ; then \ + mkdir $(MMFDIR) ; \ + echo Created directory $(MMFDIR) ; \ + fi +# Create bin directory, if necessary + @if [ ! -d $(BINDIR) ] ; then \ + mkdir $(BINDIR) ; \ + echo Created directory $(BINDIR) ; \ + fi + cd $(MMFDIR); $(MAKE); + cd $(MIZUDIR); $(MAKE); cd $(PRMSDIR); $(MAKE); + clean: cd $(MMFDIR); $(MAKE) clean; - cd $(MIZU); $(MAKE); clean; + cd $(MIZUDIR); $(MAKE) clean; cd $(PRMSDIR); $(MAKE) clean; - + cd $(PRMSRDIR); $(MAKE) clean; + $(RM) $(BINDIR)/prms* diff --git a/README.md b/README.md index aa67641c..61f15279 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,114 @@ # prms Precipitation Runoff Modeling System +This fork makes two executables when doing 'make'. The executable /bin/prms will have capabilities of turning on the glacier_flag, the stream_temp_flag, the frozen_flag (that uses the old CFGI with no frozen depth just binary frozen state), and using the strmflow_module muskingum, muskingum_mann, or mizuroute. The executable /bin/prmsrip will have capabilities of turning on the glacier_flag, the stream_temp_flag, the frozen_flag (that uses the new CFGImod with frozen depth), the ripst_flag, and using the strmflow_module muskingum, muskingum_mann, or mizuroute. + + + +Glacier params, calculate: + +abl_elev_range-- from RGI == Randolph Glacier Inventory and modify, /Roland has code +glacier_frac_init-- from RGI /Roland has code +glrette_frac_init-- from RGI, fraction of HRU that is small glacier and not delineated into special glacier HRU /Roland has code +hru_length--geometric map tools/Roland has code +hru_width--geometric map tools/Roland has code +tohru--from glacier HRU delineation/Roland has code +(Note, some standard PRMS params change on glaciers, like hru_type =4 on all possible to be glacierized HRUs). Roland's code figured out all this + + + +Glacier params, defaults, could calibrate some of these (I did in the Copper paper): + +albedo_coef +albedo_ice +glacr_freeh2o_cap +glacr_layer +glacrva_coef-- will go away if update code with Delta h area change and then might add some +glacrva_exp-- will go away if update code with Delta h area change and then might add some +max_gldepth +stor_firn +stor_ice +stor_snow + + + +Routing params calculate: + +mann_n-- can be from slope if nothing better, that is the method in CONUS and then we calibrate +seg_slope-- DEM +Muskingum: seg_depth-- from regression models or? I got it from a paper in CONUS +Dynamic: seg_width-- from regression models or NARwidth of ? I got it from a paper in CONUS for the streams smaller than 30 m wide, big is NARwidths +seg_length-- DEM? needs to be length including the vertical drop +(does not use K_coef) + + + +Frozen ground params, defaults, can calibrate to permafrost probability maps: + +cfgi_decay -- calibrate this maybe +cfgi_thrshld-- calibrate this definitely +porosity_hru-- soils data porosity +soil_depth-- Unified North American Soil Map? STATSGO? +soil_den --Unified North American Soil Map? STATSGO? + + + +Riparian needs all the new routing params (besides seg_length) and: +calculate: + +porosity_seg-- soils data porosity +ripst_areafr_max-- riparian area data gives fraction of HRU riparian +specyield_seg-- might be data on this somewhere, haven't got it for CONUS yet but default is perhaps close enough +transmiss_seg-- soils data transmissivity + + + +Riparian defaults, most of these have no good reason to calibrate: + +bankfinite_hru +bankst_head_init +ripst_et_coef +ripst_frac_init +tr_ratio +bank_height_fac -- calibrate this or at least think about what want, it's 20 right now + + + +Stream Temp calculate: + +seg_elev-- DEM +seg_lat-- DEM +width_m-- relate to seg_width somehow? + + + +Stream Temp, I'd leave these at defaults ... ask Markstrom: + +albedo +alte +altw +azrh +gw_tau +lat_temp_adj +maxiter_sntemp +melt_temp +ss_tau +vce +vcw +vdemn +vdemx +vdwmn +vdwmx +vdwmx +vhe +vhw +voe +vow +width_alpha + + + +Calving glacier params calculate: +NOTE: currently this isn't going to get put in so you can ignore (there are calving parameters needed for MWBMglacier now) + +ocean_depth -- offshore DEM at end of calving front, current idea is to set up offshore HRUs that contain the glacier tongue area and have properties +(these glaciers will have tohrus that then will be possible to tell that the upstream glacier part feeds a calving front, as long as tongue calculated to float) diff --git a/makelist b/makelist index 225bea55..c1d8f803 100644 --- a/makelist +++ b/makelist @@ -2,55 +2,62 @@ #------------------------------------------------------------------- #------------------------------------------------------------------- -MMFDIR = ./mmf -PRMSDIR = ./prms -BINDIR = ../bin -MMFLIB = .$(MMFDIR)/libmmf.a +F_MASTER = /Users/amedin/Research/NHM_dev +MMFDIR = $(F_MASTER)/mmf +MIZUDIR = $(F_MASTER)/mizu +LIBDIR = $(F_MASTER)/lib +PRMSDIR = $(F_MASTER)/prms +PRMSRDIR = $(F_MASTER)/prmsRip +BINDIR = $(F_MASTER)/bin +MMFLIB = $(LIBDIR)/libmmf.a +MIZULIB = $(LIBDIR)/libmizu.a +INCMIZU = -I$(MIZUDIR)/include ######################################################### # Configure tags for each system ########################################################## -ARC = LINUX +ARC = UNIX +#ARC = LINUX #ARC = WINDOWS #OPTLEVEL = -g OPTLEVEL = -O -Bstatic #for gfortran -#LDFLAGS =$(OPTLEVEL) +LDFLAGS =$(OPTLEVEL) #for ifort -LDFLAGS =$(OPTLEVEL) -nofor_main +#LDFLAGS =$(OPTLEVEL) -nofor_main +MFLAGS = -ffree-line-length-none ########################################################## # Define the Fortran compile flags ########################################################## #for gfortran -#FFLAGS= $(OPTLEVEL) -fbounds-check -fno-second-underscore -Wall +FFLAGS= $(OPTLEVEL) -fbounds-check -Wall -fno-second-underscore $(MFLAGS) #FFLAGS= $(OPTLEVEL) -fno-second-underscore +FC = gfortran #for ifort #FFLAGS= $(OPTLEVEL) -warn all -fltconsistency -FFLAGS= $(OPTLEVEL) -fp-model source -#FC = gfortran -FC = ifort +#FFLAGS= $(OPTLEVEL) -fp-model source +#FC = ifort ########################################################## # Define the C compile flags # -D_UF defines UNIX naming conventions for mixed language compilation. ########################################################## -CFLAGS = $(OPTLEVEL) -D$(ARC) -D_UF -Wall +CFLAGS = $(OPTLEVEL) -D$(ARC) -D_UF -I /usr/include #for gfortran -#CC = gcc +CC = gcc #for ifort -CC = icc +#CC = icc ########################################################## # Define the libraries ########################################################## #for gfortran -#MATHLIB = -lm -#GCLIB = -lgfortran -lgcc $(MATHLIB) +MATHLIB = -lm +GCLIB = -L/opt/local/lib -lgfortran -lgcc_s.1 -L/usr/lib #for ifort -MATHLIB = -GCLIB = +#GCLIB = -lgfortran -lgcc $(MATHLIB) FLIBS = $(GCLIB) ########################################################## diff --git a/prms/call_modules.f90 b/prms/call_modules.f90 index 0f12ebb4..2f4bb880 100644 --- a/prms/call_modules.f90 +++ b/prms/call_modules.f90 @@ -8,7 +8,7 @@ MODULE PRMS_MODULE CHARACTER(LEN=68), PARAMETER :: & & EQULS = '====================================================================' CHARACTER(LEN=12), PARAMETER :: MODNAME = 'call_modules' - CHARACTER(LEN=24), PARAMETER :: PRMS_VERSION = 'Version 5.0.1 06/20/2019' + CHARACTER(LEN=24), PARAMETER :: PRMS_VERSION = 'Version 5.1.0 RC' CHARACTER(LEN=8), SAVE :: Process CHARACTER(LEN=80), SAVE :: PRMS_versn INTEGER, SAVE :: Model, Process_flag, Call_cascade, Ncascade, Ncascdgw diff --git a/prms/climateflow.f90 b/prms/climateflow.f90 index 229b26df..f41b1af8 100644 --- a/prms/climateflow.f90 +++ b/prms/climateflow.f90 @@ -542,17 +542,17 @@ INTEGER FUNCTION climateflow_decl() IF ( Glacier_flag==1 .OR. Model==99 ) THEN ALLOCATE ( Glacier_frac(Nhru) ) IF ( declvar(MODNAME, 'glacier_frac', 'nhru', Nhru, 'real', & - 'Fraction of glaciation (0=none; 1=100%)', & + 'Fraction of glaciation (0=none; 1=100%) in glacier-capable HRU', & 'decimal fraction', Glacier_frac)/=0 ) CALL read_error(3, 'glacier_frac') ALLOCATE ( Glrette_frac(Nhru) ) IF ( declvar(MODNAME, 'glrette_frac', 'nhru', Nhru, 'real', & - 'Fraction of snow field (too small for glacier dynamics)', & + 'Fraction of glacierette (too small for glacier dynamics)', & 'decimal fraction', Glrette_frac)/=0 ) CALL read_error(3, 'glrette_frac') ALLOCATE ( Alt_above_ela(Nhru) ) IF ( declvar(MODNAME, 'alt_above_ela', 'nhru', Nhru, 'real', & - 'Altitude above equilibrium line altitude (ELA)', & + 'Altitude HRU is above equilibrium line altitude (ELA), negative value indicates HRU below ELA', & 'elev_units', Alt_above_ela)/=0 ) CALL read_error(3, 'alt_above_ela') ENDIF diff --git a/prms/dynamic_param_read.f90 b/prms/dynamic_param_read.f90 index 3d32466d..ba254e50 100644 --- a/prms/dynamic_param_read.f90 +++ b/prms/dynamic_param_read.f90 @@ -54,7 +54,7 @@ INTEGER FUNCTION dynamic_param_read() IF ( Process(:3)=='run' ) THEN dynamic_param_read = dynparamrun() ELSEIF ( Process(:4)=='decl' ) THEN - Version_dynamic_param_read = 'dynamic_param_read.f90 2019-05-30 13:50:00Z' + Version_dynamic_param_read = 'dynamic_param_read.f90 2019-09-06 16:05:00Z' CALL print_module(Version_dynamic_param_read, 'Time Series Data ', 90) !MODNAME = 'dynamic_param_read' ELSEIF ( Process(:4)=='init' ) THEN @@ -426,12 +426,14 @@ INTEGER FUNCTION dynparamrun() EXTERNAL write_dynoutput, is_eof, write_dynparam, write_dynparam_int EXTERNAL write_dynparam_potet ! Local Variables - INTEGER :: i, istop, check_dprst_depth_flag + INTEGER :: i, istop, check_dprst_depth_flag, check_sm_max_flag, check_srechr_max_flag REAL :: harea, frac_imperv, tmp, hruperv, dprstfrac, soil_adj CHARACTER(LEN=30), PARAMETER :: fmt1 = '(A, I0, ":", I5, 2("/",I2.2))' !*********************************************************************** dynparamrun = 0 istop = 0 + check_srechr_max_flag = 0 + check_sm_max_flag = 0 IF ( Imperv_frac_flag==1 .OR. Dprst_frac_flag==1 .OR. Dprst_depth_flag==1 ) THEN Check_imperv = 0 @@ -582,6 +584,7 @@ INTEGER FUNCTION dynparamrun() ENDIF Basin_soil_moist = Basin_soil_moist + DBLE( Soil_moist(i)*Hru_perv(i) ) Basin_soil_rechr = Basin_soil_rechr + DBLE( Soil_rechr(i)*Hru_perv(i) ) + Soil_moist_tot(i) = Ssres_stor(i) + Soil_moist(i)*Hru_frac_perv(i) ENDDO Basin_soil_moist = Basin_soil_moist*Basin_area_inv Basin_soil_rechr = Basin_soil_rechr*Basin_area_inv @@ -668,33 +671,34 @@ INTEGER FUNCTION dynparamrun() CALL write_dynparam(Output_unit, Nhru, Updated_hrus, Temp, Potet_coef(1,Nowmonth), 'potet_coef') ENDIF CALL is_eof(Potetcoef_unit, Potetcoef_next_yr, Potetcoef_next_mo, Potetcoef_next_day) - ENDIF - IF ( Et_flag==1 ) THEN ! potet_jh - IF ( Dyn_potet_flag==1 ) THEN - Jh_coef = Potet_coef - ELSE - DO i = 1, Nhru - Jh_coef_hru(i) = Potet_coef(i,Nowmonth) - ENDDO - ENDIF - ELSEIF ( Et_flag==7 ) THEN ! climate_hru - Potet_cbh_adj = Potet_coef - ELSEIF ( Et_flag==11 ) THEN ! potet_pm - IF ( Dyn_potet_flag==1 ) THEN - Pm_n_coef = Potet_coef - ELSE - Pm_d_coef = Potet_coef + + IF ( Et_flag==1 ) THEN ! potet_jh + IF ( Dyn_potet_flag==1 ) THEN + Jh_coef = Potet_coef + ELSE + DO i = 1, Nhru + Jh_coef_hru(i) = Potet_coef(i,Nowmonth) + ENDDO + ENDIF + ELSEIF ( Et_flag==7 ) THEN ! climate_hru + Potet_cbh_adj = Potet_coef + ELSEIF ( Et_flag==11 ) THEN ! potet_pm + IF ( Dyn_potet_flag==1 ) THEN + Pm_n_coef = Potet_coef + ELSE + Pm_d_coef = Potet_coef + ENDIF + ELSEIF ( Et_flag==5 ) THEN ! potet_pt + Pt_alpha = Potet_coef + ELSEIF ( Et_flag==10 ) THEN ! potet_hs + Hs_krs = Potet_coef + ELSEIF ( Et_flag==2 ) THEN ! potet_hamon + Hamon_coef = Potet_coef + !ELSEIF ( Et_flag==6 ) THEN ! potet_jh_hru + !Jh_coef_hru2 = Potet_coef + ELSEIF ( Et_flag==4 ) THEN ! potet_pan + Epan_coef = Potet_coef ENDIF - ELSEIF ( Et_flag==5 ) THEN ! potet_pt - Pt_alpha = Potet_coef - ELSEIF ( Et_flag==10 ) THEN ! potet_hs - Hs_krs = Potet_coef - ELSEIF ( Et_flag==2 ) THEN ! potet_hamon - Hamon_coef = Potet_coef - !ELSEIF ( Et_flag==6 ) THEN ! potet_jh_hru - !Jh_coef_hru2 = Potet_coef - ELSEIF ( Et_flag==4 ) THEN ! potet_pan - Epan_coef = Potet_coef ENDIF ENDIF ENDIF @@ -750,6 +754,7 @@ INTEGER FUNCTION dynparamrun() CALL write_dynparam(Output_unit, Nhru, Updated_hrus, Temp, Soil_rechr_max_frac, 'soil_rechr_max_frac') ENDIF CALL is_eof(Soil_rechr_unit, Soil_rechr_next_yr, Soil_rechr_next_mo, Soil_rechr_next_day) + check_srechr_max_flag = 1 ENDIF ENDIF ENDIF @@ -761,11 +766,12 @@ INTEGER FUNCTION dynparamrun() READ ( Soil_moist_unit, * ) Soil_moist_next_yr, Soil_moist_next_mo, Soil_moist_next_day, Temp CALL write_dynparam(Output_unit, Nhru, Updated_hrus, Temp, Soil_moist_max, 'soil_moist_max') CALL is_eof(Soil_moist_unit, Soil_moist_next_yr, Soil_moist_next_mo, Soil_moist_next_day) + check_sm_max_flag = 1 ENDIF ENDIF ENDIF - IF ( Soilmoist_flag==1 .OR. Soilrechr_flag==1 ) THEN + IF ( check_sm_max_flag==1 .OR. check_srechr_max_flag==1 ) THEN DO i = 1, Nhru IF ( Hru_type(i)==2 .OR. Hru_type(i)==0 ) CYCLE ! skip lake and inactive HRUs @@ -786,7 +792,6 @@ INTEGER FUNCTION dynparamrun() CYCLE ENDIF 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_stor_max(i) = Soil_moist_max(i) - Soil_rechr_max(i) Replenish_frac(i) = Soil_rechr_max(i)/Soil_moist_max(i) ENDDO diff --git a/prms/glacr_melt.f90 b/prms/glacr_melt.f90 index ebc429db..eba55ca8 100644 --- a/prms/glacr_melt.f90 +++ b/prms/glacr_melt.f90 @@ -27,7 +27,7 @@ ! HRUs with glaciers must have parameter glacier_frac(i)=1, unless they ! are at the terminus of the glacier (in which case they can have ! glacier_frac(i)<1). Hru numbering goes from largest HRU ID at top of glacier to -! smallest at ID at bottom (the way Weasel delineation was designed). The parameter +! smallest ID at bottom (the way Weasel delineation was designed). The parameter ! Glac_HRUnum_down = 1 then in the init function. If the opposite direction, ! then set Glac_HRUnum_down = 0. IDs need to be stacked. ! @@ -51,11 +51,10 @@ MODULE PRMS_GLACR ! Local Variables ! Ngl - Number of glaciers counted by termini - ! Ntp - Number of tops of glaciers, so max glaciers that could ever split in two - ! Nhrugl - Number of at least partially glacierized hrus at initiation -!#of cells=Nhrugl,#of streams=Ntp,#of cells/stream<=Ntp, #of glaciers<=Nhru - INTEGER, SAVE :: Nglres, Ngl, Ntp, Nhrugl, MbInit_flag, Output_unit, Fraw_unit, All_unit - INTEGER, SAVE :: Seven, Four, Glac_HRUnum_down + ! Ntp - Number of tops of glaciers, so max glaciers that could ever split in too + ! Nhrugl - Number of at glacier-capable hrus + INTEGER, SAVE :: Nglres, Ngl, Ntp, Nhrugl, Mbinit_flag, Output_unit, Fraw_unit, All_unit + INTEGER, SAVE :: Seven, Four, Glac_hrunum_down DOUBLE PRECISION, SAVE, ALLOCATABLE :: Hru_area_inch2(:) REAL, PARAMETER :: Gravity = 9.8 ! m/s2 REAL, PARAMETER :: Aflow = 1.e-25 ! Pa^-3/s, Farinotti 2009 could be 2.4e-24, could be 1e-26 see Patterson 2010 @@ -70,11 +69,12 @@ MODULE PRMS_GLACR REAL, SAVE, ALLOCATABLE :: Basal_elev(:), Basal_slope(:), Keep_gl(:,:), Prev_outi(:, :), Prev_out(:, :) REAL, SAVE, ALLOCATABLE :: Ode_glacrva_coef(:), Av_basal_slope(:), Av_fgrad(:), Hru_slope_ts(:) REAL, SAVE, ALLOCATABLE :: Hru_mb_yrend(:), Glacr_flow(:), Glacr_slope_init(:), Gl_top_melt(:) + REAL, SAVE, ALLOCATABLE :: Hru_glthick(:) INTEGER, SAVE, ALLOCATABLE :: Top(:), Term(:), Top_tag(:), Ela(:), Order_flowline(:) INTEGER, SAVE, ALLOCATABLE :: Glacr_tag(:), Ikeep_gl(:,:), Tohru(:) DOUBLE PRECISION, SAVE :: Basin_gl_ice_melt, Basin_gl_area, Basin_gl_top_melt DOUBLE PRECISION, SAVE :: Basin_gl_top_gain, Basin_gl_storvol, Basin_gl_storage - DOUBLE PRECISION, SAVE :: Basin_gl_storstart + DOUBLE PRECISION, SAVE :: Basin_gl_storstart, Basin_glthick DOUBLE PRECISION, SAVE, ALLOCATABLE :: Hru_mb_yrcumul(:), Delta_volyr(:), Prev_vol(:) DOUBLE PRECISION, SAVE, ALLOCATABLE :: Prev_area(:), Gl_mb_yrcumul(:), Gl_area(:) DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gl_mb_cumul(:), Glnet_ar_delta(:), Gl_mbc_yrend(:) @@ -85,7 +85,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 :: Hru_slope(:), Abl_elev_range(:) + REAL, SAVE, ALLOCATABLE :: Hru_slope(:), Abl_elev_range(:), Basal_elev_set(:), Basal_slope_set(:) END MODULE PRMS_GLACR @@ -119,7 +119,7 @@ END FUNCTION GLACR ! glacrsetdims - declares glacier module specific dimensions !*********************************************************************** INTEGER FUNCTION glacrsetdims() - USE PRMS_GLACR, ONLY: Nglres, Seven, Four, MbInit_flag + USE PRMS_GLACR, ONLY: Nglres, Seven, Four, Mbinit_flag IMPLICIT NONE ! Functions INTEGER, EXTERNAL :: declfix, control_integer @@ -134,7 +134,8 @@ INTEGER FUNCTION glacrsetdims() IF ( declfix('four',4, 4, 'Need for keeping glacier variable integer array')/=0 ) CALL read_error(7, 'four') Four = 4 - IF ( control_integer(MbInit_flag, 'mbInit_flag')/=0 ) MbInit_flag = 0 + IF ( control_integer(Mbinit_flag, 'mbinit_flag')/=0 ) Mbinit_flag = 0 + print*,Mbinit_flag END FUNCTION glacrsetdims @@ -167,16 +168,25 @@ INTEGER FUNCTION glacrdecl() & 'HRU slope for timestep, which can change for glaciers', & & 'decimal fraction', Hru_slope_ts)/=0 ) CALL read_error(3, 'hru_slope_ts') + ALLOCATE ( Hru_glthick(Nhru) ) + IF ( declvar(MODNAME, 'hru_glthick', 'nhru', Nhru, 'real', & + & 'HRU thickness for glacier, changes through time', & + & 'elev_units', Hru_glthick)/=0 ) CALL read_error(3, 'hru_glthick') + + IF ( declvar(MODNAME, 'basin_glthick', 'one', 1, 'double', & + & 'Basin average thickness of glaciers (glacierettes not included)', & + & 'elev_units', Basin_glthick)/=0 ) CALL read_error(3, 'basin_glthick') + IF ( declvar(MODNAME, 'basin_gl_top_melt', 'one', 1, 'double', & - & 'Basin area-weighted glacier surface melt (snow, ice and rain) coming out of termini of all glaciers and glrettes', & + & 'Basin area-weighted glacier surface melt (snow, ice and rain) coming out of termini of all glaciers and glacierettes', & & 'inches', Basin_gl_top_melt)/=0 ) CALL read_error(3, 'basin_gl_top_melt') IF ( declvar(MODNAME, 'basin_gl_top_gain', 'one', 1, 'double', & - & 'Basin area-weighted glacier surface gain (snow and rain minus evap) for all glaciers and glrettes', & + & 'Basin area-weighted glacier surface gain (snow and rain minus evap) for all glaciers and glacierettes', & & 'inches', Basin_gl_top_gain)/=0 ) CALL read_error(3, 'basin_gl_top_gain') IF ( declvar(MODNAME, 'basin_gl_ice_melt', 'one', 1, 'double', & - & 'Basin area-weighted glacier ice (firn) melt coming out of termini of all glaciers and glrettes', & + & 'Basin area-weighted glacier ice (no snow) melt coming out of termini of all glaciers and glacierettes', & & 'inches', Basin_gl_ice_melt)/=0 ) CALL read_error(3, 'basin_gl_ice_melt') ALLOCATE ( Gl_mb_yrcumul(Nhru) ) @@ -190,7 +200,7 @@ INTEGER FUNCTION glacrdecl() & 'inches', Gl_mb_cumul)/=0 ) CALL read_error(3, 'gl_mb_cumul') IF ( declvar(MODNAME, 'basin_gl_area', 'one', 1, 'double', & - & 'Basin area-weighted average glacier-covered area', & + & 'Basin area-weighted average glacier and glacierette-covered area', & & 'decimal fraction', Basin_gl_area)/=0 ) CALL read_error(3, 'basin_gl_area') ALLOCATE ( Gl_area(Nhru) ) @@ -205,7 +215,7 @@ INTEGER FUNCTION glacrdecl() ALLOCATE ( Glacr_flow(Nhru) ) IF ( declvar(MODNAME, 'glacr_flow', 'nhru', Nhru, 'real', & - & 'Glacier melt and rain from HRU to stream network, only nonzero at termini HRUs and snowfield HRUs', & + & 'Glacier melt and rain from HRU to stream network, only nonzero at termini HRUs and glacierette HRUs', & & 'inches cubed', Glacr_flow)/=0 ) CALL read_error(3, 'glacr_flow') ALLOCATE ( Delta_volyr(Nhru) ) @@ -300,7 +310,7 @@ INTEGER FUNCTION glacrdecl() ALLOCATE ( Gl_ice_melt(Nhru) ) IF ( declvar(MODNAME, 'gl_ice_melt', 'nhru', Nhru, 'real', & - & 'Amount of glacier ice (firn) melt coming out of terminus of glacier, indexed by Glacr_tag', & + & 'Amount of glacier ice (no snow) melt coming out of terminus of glacier, indexed by Glacr_tag', & & 'inches', Gl_ice_melt)/=0 ) CALL read_error(3, 'gl_ice_melt') ALLOCATE ( Basal_elev(Nhru) ) @@ -308,15 +318,7 @@ INTEGER FUNCTION glacrdecl() & 'Glacier basal elevation mean over HRU', & & 'elev_units', Basal_elev)/=0 ) CALL read_error(3, 'basal_elev') - ALLOCATE ( Keep_gl(Nhru,Seven) ) - IF ( declvar(MODNAME, 'keep_gl', 'nhru,seven', Nhru*Seven, 'real', & - & 'Glacier real variables keeping from first year', & - & 'none', Keep_gl)/=0 ) CALL read_error(3, 'keep_gl') - - ALLOCATE ( Ikeep_gl(Nhru,Four) ) - IF ( declvar(MODNAME, 'ikeep_gl', 'nhru,four', Nhru*Four, 'real', & - & 'Glacier integer variables keeping from first year', & - & 'none', Ikeep_gl)/=0 ) CALL read_error(3, 'ikeep_gl') + ALLOCATE ( Keep_gl(Nhru,Seven), Ikeep_gl(Nhru,Four) ) ALLOCATE ( Basal_slope(Nhru) ) IF ( declvar(MODNAME, 'basal_slope', 'nhru', Nhru, 'real', & @@ -329,16 +331,15 @@ INTEGER FUNCTION glacrdecl() & 'decimal fraction', Av_basal_slope)/=0 ) CALL read_error(3, 'av_basal_slope') IF ( declvar(MODNAME, 'basin_gl_storage', 'one', 1, 'double', & - & 'Basin area-weighted average storage change in glacier reservoirs', & + & 'Basin area-weighted average storage change in glacier and glacierette reservoirs', & & 'inches', Basin_gl_storage)/=0 ) CALL read_error(3, 'basin_gl_storage') IF ( declvar(MODNAME, 'basin_gl_storstart', 'one', 1, 'double', & - & 'Basin area-weighted average storage estimated start in glacier reservoirs', & + & 'Basin area-weighted average storage estimated start in glacier and glacierette reservoirs', & & 'inches', Basin_gl_storstart)/=0 ) CALL read_error(3, 'basin_gl_storstart') - IF ( declvar(MODNAME, 'basin_gl_storvol', 'one', 1, 'double', & - & 'Basin storage volume in glacier storage reservoirs', & + & 'Basin storage volume in glacier and glacierette storage reservoirs', & & 'acre-inches', Basin_gl_storvol)/=0 ) CALL read_error(3, 'basin_gl_storvol') IF ( Init_vars_from_file==0 ) THEN @@ -416,25 +417,47 @@ INTEGER FUNCTION glacrdecl() ALLOCATE ( Hru_length(Nhru) ) IF ( declparam(MODNAME, 'hru_length', 'nhru', 'real', & '0.0', '0.0', '10000.0', & - 'Length of segment covering all of glacier-possible HRU', & - 'Length of segment covering all of glacier-possible HRU', & + 'Length of segment covering all of glacier-capable HRU', & + 'Length of segment covering all of glacier-capable HRU', & 'km')/=0 ) CALL read_error(1, 'hru_length') ALLOCATE ( Hru_width(Nhru) ) IF ( declparam(MODNAME, 'hru_width', 'nhru', 'real', & '0.0', '0.0', '10000.0', & - 'Width of glacier-possible HRU', & - 'Width of glacier-possible HRU', & + 'Width of glacier-capable HRU', & + 'Width of glacier-capable HRU', & 'km')/=0 ) CALL read_error(1, 'hru_width') ALLOCATE ( Abl_elev_range(Nhru) ) IF ( declparam(MODNAME, 'abl_elev_range', 'nhru', 'real', & '1000.0', '0.0', '17000.0', & - 'Average HRU snowfield ablation zones elevation range', & - 'Average HRU snowfield ablation zones elevation range or ~ median-min elev', & + 'Average HRU glacierette ablation zones elevation range', & + 'Average HRU glacierette ablation zones elevation range or ~ median-min elev', & 'elev_units')/=0 ) CALL read_error(1, 'abl_elev_range') - END FUNCTION glacrdecl + IF ( declparam(MODNAME, 'glac_hrunum_down', 'one', 'integer', & + '1', '0', '1', & + '1 is Glacier HRU numbering from largest HRU ID at glacier top to smallest ID at terminus, O opposite, IDs stacked', & + '1 is Glacier HRU numbering from largest HRU ID at glacier top to smallest ID at terminus, O opposite, IDs stacked', & + 'none')/=0 ) CALL read_error(1, 'glac_hrunum_down') + + IF (Mbinit_flag==1) THEN + ALLOCATE ( Basal_elev_set(Nhru) ) + IF ( declparam(MODNAME, 'basal_elev_set', 'nhru', 'real', & + '0.0', '-1000.0', '30000.0', & + 'Glacier basal elevation mean over HRU inputted from outside information', & + 'Glacier basal elevation mean over HRU inputted from outside information', & + 'elev_units')/=0 ) CALL read_error(1, 'basal_elev_set') + + ALLOCATE ( Basal_slope_set(Nhru) ) + IF ( declparam(MODNAME, 'basal_slope_set', 'nhru', 'real', & + '0.0', '0.0', '10.0', & + 'Glacier basal slope down flowline mean over HRU inputted from outside information', & + 'Glacier basal slope down flowline mean over HRU inputted from outside information', & + 'decimal fraction')/=0) CALL read_error(1, 'basal_slope_set') + ENDIF + + END FUNCTION glacrdecl !*********************************************************************** ! glacrinit - Initialize glacr module - get parameter values @@ -477,6 +500,7 @@ INTEGER FUNCTION glacrinit() IF ( getparam(MODNAME, 'abl_elev_range', Nhru, 'real', Abl_elev_range)/=0 ) CALL read_error(2, 'abl_elev_range') IF ( getparam(MODNAME, 'tohru', Nhru, 'integer', Tohru)/=0 ) CALL read_error(2, 'tohru') IF ( getparam(MODNAME, 'hru_slope', Nhru, 'real', Hru_slope)/=0 ) CALL read_error(2, 'hru_slope') + IF ( getparam(MODNAME, 'glac_hrunum_down', 1, 'one', Glac_hrunum_down)/=0 ) CALL read_error(2, 'glac_hrunum_down') IF ( Init_vars_from_file==0 ) THEN Alt_above_ela = 0.0 Prev_out = 0.0 @@ -498,10 +522,23 @@ INTEGER FUNCTION glacrinit() Gl_mb_cumul = 0.0D0 Gl_mbc_yrend = 0.0D0 Hru_slope_ts = Hru_slope - Basal_elev = Hru_elev_ts ! Hru_elev_ts always set in basin, need in case of restart - Basal_slope = Hru_slope_ts + IF (Mbinit_flag/=1) THEN + Basal_elev = Hru_elev_ts ! Hru_elev_ts always set in basin, need in case of restart + Basal_slope = Hru_slope_ts + ELSE !get from parameters + IF ( getparam(MODNAME, 'basal_elev_set', Nhru, 'real', Basal_elev_set)/=0 ) CALL read_error(2, 'basal_elev_set') + Basal_elev = Basal_elev_set + IF ( getparam(MODNAME, 'basal_slope_set', Nhru, 'real', Basal_slope_set)/=0 ) CALL read_error(2, 'basal_slope_set') + Basal_slope = Basal_slope_set + ENDIF + Basin_glthick = 0.0D0 + DO jj = 1, Active_hrus + j = Hru_route_order(jj) + Hru_glthick(j) = Hru_elev_ts(j) - Basal_elev(j) + Basin_glthick = Basin_glthick + Hru_glthick(j)*Hru_area(j)*Acre_inch2 + ENDDO Av_basal_slope = 0.0 - Glacr_elev_init = Hru_elev_ts + Glacr_elev_init = Hru_elev_ts !if Mbinit_flag>=2, this will change Glacr_slope_init = Hru_slope_ts Av_fgrad = 0.0 Basin_gl_top_melt = 0.0D0 @@ -515,10 +552,6 @@ INTEGER FUNCTION glacrinit() Basin_gl_storvol = 0.0D0 ENDIF - 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. - ! 0 is terminus is smallest ID and top is largest. IDs are stacked. - hru_flowline = 0 toflowline = 0 str_idm = 1.0E15 @@ -542,9 +575,9 @@ INTEGER FUNCTION glacrinit() count = 1 !has at least one glacier glacier_frac_use(j) = 1.0 !should be end of extensions or branches-- will fail if don't set up with indices stacked - IF ( Glac_HRUnum_down==1) THEN + IF ( Glac_hrunum_down==1) THEN IF (Tohru(j)/=j-1 ) glacier_frac_use(j) = 0.999 - ELSEIF ( Glac_HRUnum_down==0) THEN + ELSEIF ( Glac_hrunum_down==0) THEN IF (Tohru(j)/=j+1 ) glacier_frac_use(j) = 0.999 ENDIF ENDIF @@ -708,9 +741,9 @@ INTEGER FUNCTION glacrinit() glacier_frac_use(j)= Glacier_frac(j) !should be end of extensions or branches-- will fail if don't set up with indices stacked ! making it so has no connected branches because branching bottom calculations don't work - IF ( Glac_HRUnum_down==1) THEN + IF ( Glac_hrunum_down==1) THEN IF (Tohru(j)/=j-1 .AND. glacier_frac_use(j)==1.0 ) glacier_frac_use(j) = 0.999 - ELSEIF ( Glac_HRUnum_down==0) THEN + ELSEIF ( Glac_hrunum_down==0) THEN IF (Tohru(j)/=j+1 .AND. glacier_frac_use(j)==1.0 ) glacier_frac_use(j) = 0.999 ENDIF ENDIF @@ -768,6 +801,8 @@ INTEGER FUNCTION glacrinit() Gl_area(p) = curr_area(Term(o))/Acre_inch2 !print*, 'Glacr_tag', p, ', area acres branches=', Gl_area(p), ', terminus HRU=', Term(o) ENDDO +!average thickness over glaciers only so take area early + Basin_glthick = Basin_glthick/Basin_gl_area!top and bottom /Acre_inch2*Basin_area_inv to give real thing, cancel DO i = 1, Active_hrus j = Hru_route_order(i) IF ( Hru_type(j)==1 ) Basin_gl_area = Basin_gl_area + DBLE(Glrette_frac(j))*Hru_area_inch2(j) @@ -783,9 +818,8 @@ INTEGER FUNCTION glacrinit() ENDDO ENDDO !******Compute basin weighted averages - ! Basin_area_inv is in 1/acres, Basin_gl_area in inches squared + ! Basin_area_inv is in 1/acres, Basin_gl_area in inches squared, put in fraction Basin_gl_area = (Basin_gl_area/Acre_inch2)*Basin_area_inv - !print*, 'Basin area acres=', 1.0/Basin_area_inv ENDIF ! skip all if no glaciers ! @@ -813,7 +847,7 @@ INTEGER FUNCTION glacrrun() i = Hru_route_order(j) IF ( Hru_type(i)==1 ) THEN IF (Glrette_frac(j)>NEARZERO) THEN - count=1 !has at least one snowfield + count=1 !has at least one glacierette EXIT ENDIF ENDIF @@ -833,12 +867,13 @@ INTEGER FUNCTION glacrrun() Gl_mbc_yrend = 0.0D0 Av_basal_slope = 0.0 Av_fgrad = 0.0 + Basin_glthick = 0.0D0 Hru_slope_ts = Basal_slope Alt_above_ela = 0.0 ! doesn't matter if no glaciers dosol = recompute_soltab() ! change soltab tables for Hru_slope_ts IF (count==0) THEN Glacr_flow = 0.0 - Basin_gl_area = 0.D0 !no snowfields either + Basin_gl_area = 0.D0 !no glacierettes either Basin_gl_top_melt = 0.0D0 Basin_gl_top_gain = 0.0D0 Basin_gl_ice_melt = 0.0D0 @@ -903,9 +938,12 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) INTEGER, INTENT(IN) :: glacr_exist, glrette_exist !*********************************************************************** comp_glsurf = 1 - dobot = 1 ! 1 calls bottom calcs, 0 doesn't: Set to 0 for calibrating, then run one extra step with it on - ! Should change so that saves the basal elevations (or reads in as parameter) and then recalibrating does not change + dobot = 1 ! Should change so that saves the basal elevations (or reads in as parameter) and then recalibrating does not change botwrite = 0 ! 1 writes bottom calcs, 0 doesn't: Set to 0 for calibrating + IF (Mbinit_flag<2) THEN ! know the bottom, read as parameter + dobot = 0 + botwrite = 0 + ENDIF ! initialize ela_elevt = 0.0 gt = 0 @@ -950,7 +988,7 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) ENDIF ENDDO ! ELA calculations - IF ( MBinit_flag==2 ) THEN + IF ( Mbinit_flag==3 ) THEN doela = compute_ela_aar() !want steady state ELA estimation for fraw calc DO j = 1, Ntp ela_elevt(j)=Hru_elev(Ela(j)) !will scale inside subroutine, want initial one without _ts @@ -976,9 +1014,9 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) glacier_frac_use(j)= Glacier_frac(j) !should be end of extensions or branches-- will fail if don't set up with indices stacked ! making it so has no connected branches because branching bottom calculations don't work - IF ( Glac_HRUnum_down==1) THEN + IF ( Glac_hrunum_down==1) THEN IF (Tohru(j)/=j-1 .AND. glacier_frac_use(j)==1.0 ) glacier_frac_use(j) = 0.999 - ELSEIF ( Glac_HRUnum_down==0) THEN + ELSEIF ( Glac_hrunum_down==0) THEN IF (Tohru(j)/=j+1 .AND. glacier_frac_use(j)==1.0 ) glacier_frac_use(j) = 0.999 ENDIF ENDIF @@ -1044,8 +1082,9 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) aream(cell_id(i)) = 0.0 !indexed by terminus HRU ENDDO ENDIF +! Need to do this so that Hru_elev_ts is actually the same as Hru_elev before melt in terminus, won't matter if hru_elev = basal_elev +! Do for all dobot cases, all Mbinit_flag choices DO i = 1, Nhrugl -! Need to do this so that Hru_elev_ts is actually the same as Hru_elev before melt in terminus IF ( Glacier_frac(cell_id(i))>NEARZERO) THEN !only effects terminus Glacr_elev_init(cell_id(i)) = (Hru_elev(cell_id(i)) - (1.0-Glacier_frac(cell_id(i))) & & *Basal_elev((cell_id(i))))/Glacier_frac(cell_id(i)) @@ -1301,11 +1340,17 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) ENDIF ENDDO ENDDO + Basin_glthick = 0.0D0 DO jj = 1, Active_hrus j = Hru_route_order(jj) IF ( Hru_type(j)==4 ) THEN curr_area(j) = curr_area(j) + add_area(j) curr_areap(j) = curr_areap(j) + add_areap(j) +! have bottom, compute Hru_elev_ts related stuff + Hru_elev_ts(j) = Glacier_frac(j)*Glacr_elev_init(j) + (1.0-Glacier_frac(j))*Basal_elev(j) + Hru_slope_ts(j) = Glacier_frac(j)*Glacr_slope_init(j) + (1.0-Glacier_frac(j))*Basal_slope(j) + Hru_glthick(j) = Hru_elev_ts(j) - Basal_elev(j) + Basin_glthick = Basin_glthick + Hru_glthick(j)*Hru_area_inch2(j) ENDIF ENDDO Gl_area = 0.D0 @@ -1317,14 +1362,14 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) Basin_gl_area = Basin_gl_area + curr_area(Term(o)) !print*, 'Glacr_tag', p, ', area acres', Nowyear,' =', Gl_area(p), ', terminus HRU=', Term(o) ENDDO -! have bottom, compute Hru_elev_ts related stuff +!average thickness over glaciers only so take area early + Basin_glthick = Basin_glthick/Basin_gl_area!top and bottom /Acre_inch2*Basin_area_inv to give real thing, cancel + DO ii = 1, Ntp DO i = 1, Active_hrus j = Hru_route_order(i) IF ( Hru_type(j)==4 ) THEN IF ( Top_tag(j)==Top_tag(Top(ii)) ) Alt_above_ela(j) = Hru_elev_ts(j)- Hru_elev_ts(Ela(ii)) - Hru_elev_ts(j) = Glacier_frac(j)*Glacr_elev_init(j) + (1.0-Glacier_frac(j))*Basal_elev(j) - Hru_slope_ts(j) = Glacier_frac(j)*Glacr_slope_init(j) + (1.0-Glacier_frac(j))*Basal_slope(j) ENDIF ENDDO ENDDO @@ -1355,9 +1400,9 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) ENDIF ENDDO ENDIF -!Snowfield area change uses Baumann and Winkler 2010 to change area every 10 years; -! technically each snowfield should have own ablation elevation range. - IF (glrette_exist==1) THEN !have snowfields, +!Glacierette area change uses Baumann and Winkler 2010 to change area every 10 years; +! technically each glacierette should have own ablation elevation range. + IF (glrette_exist==1) THEN !have glacierettes, IF ( MOD(Nowyear-Starttime(1),10)==0 ) THEN !change them DO i = 1, Active_hrus j = Hru_route_order(i) @@ -1379,10 +1424,9 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) ! !******Compute basin weighted averages (to units of fraction area) -! Basin_area_inv is in 1/acres, Basin_area_inv is in 1/acres, Basin_gl_area in inches squared +! Basin_area_inv is in 1/acres, Basin_gl_area in inches squared, put in fraction Basin_gl_area = (Basin_gl_area/Acre_inch2)*Basin_area_inv ENDIF !get out of start-of-year computations - ! ! Melt runoff calculations, every day DO i = 1, Active_hrus @@ -1414,6 +1458,7 @@ INTEGER FUNCTION comp_glsurf(glacr_exist, glrette_exist) gl_snow = Glrette_frac(j)*(Net_rain(j)+Net_Snow(j)) !Pk_precip is zero if no snow, so don't use gl_evap = Glrette_frac(j)*(Snow_evap(j) + Glacr_evap(j)/Glrette_frac(j)) gl_gain(j) = DBLE(gl_snow - gl_evap) + Basin_gl_top_gain = Basin_gl_top_gain + gl_gain(j)*Hru_area_inch2(j) ENDIF ENDIF ENDDO @@ -1919,9 +1964,9 @@ END SUBROUTINE tag_count ! subroutine bottom - calculates bottom topo using Salamatin and Mazo ! equations (1985) without optimization for steady state, instead ! needs a proxy for steady state mass balance. Can do this from mass -! balance calculation with climate data first year (MbInit_flag=1) +! balance calculation with climate data first year (Mbinit_flag=2) ! or use max and min balance above and below Ela, respectively and assume -! constant mass balance gradient above and below Ela; e.g. Farinotti (MbInit_flag=2) +! constant mass balance gradient above and below Ela; e.g. Farinotti (Mbinit_flag=3) ! All mass balances are adjust to put glacier in steady state. ! ! The method of Salamatin and Mazo is from conservation of mass, @@ -2061,7 +2106,7 @@ SUBROUTINE bottom(Frawt, Gln, Gt, Topn, Gl_top, Av_elev, Ela_elevt, Flow_slope, ! then add(x)=-int(f(z(x))*area(z(x)))dz/int(area(z(x)))dz ! if add is constant in z, then constant in x because no matter what z is, same add !If only one Hru on glacier, fraw will all be 0 to be in steady state - ! add would be zero if did the full integral with Mbinit_flag==2, but wrong at the moment + ! add would be zero if did the full integral with Mbinit_flag==3, but wrong at the moment hf(1) = hrawe(1)*frawe(1) DO i = 2, len_str_true+1 hf(i) = hraw(i-1)*fraw(i-1) @@ -2392,7 +2437,7 @@ SUBROUTINE getf_fgrad(ela_elev,len_str,sh,urawe,uraw,xraw,hvec,frawe,fraw,fd) ela_i = 0 ela_x = 0.0 fd = 0.0 - IF (Mbinit_flag==1) THEN !use climate data for mass balance + IF (Mbinit_flag==2) THEN !use climate data for mass balance IF (len_str>1) THEN frawe(1) = fraw(1)- (fraw(2)-fraw(1))/hvec(2)*hvec(1) frawe(2) = fraw(len_str) + (fraw(len_str)-fraw(len_str-1))/hvec(len_str)*(1-xraw(len_str)) @@ -2400,7 +2445,7 @@ SUBROUTINE getf_fgrad(ela_elev,len_str,sh,urawe,uraw,xraw,hvec,frawe,fraw,fd) frawe(1) = fraw(1) frawe(2) = fraw(1) ENDIF - ELSEIF (Mbinit_flag==2) THEN !Farinotti method, using ela_x + ELSEIF (Mbinit_flag==3) THEN !Farinotti method, using ela_x DO i = 1, len_str IF (ela_elev-uraw(i)<0.0001) THEN !has rounding errors ela_x=xraw(i) @@ -3038,7 +3083,7 @@ SUBROUTINE glacr_restart(In_out) WRITE ( Restart_outunit ) MODNAME WRITE ( Restart_outunit ) Nhrugl, Basin_gl_top_melt, Gl_mb_yrcumul WRITE ( Restart_outunit ) Gl_mb_cumul, Glnet_ar_delta, Gl_mbc_yrend - WRITE ( Restart_outunit ) Basin_gl_top_gain + WRITE ( Restart_outunit ) Basin_gl_top_gain, Basin_glthick WRITE ( Restart_outunit ) Basin_gl_area, Gl_area, Basin_gl_ice_melt WRITE ( Restart_outunit ) Hru_glres_melt, Basin_gl_storstart WRITE ( Restart_outunit ) Gl_top_melt, Basin_gl_storage, Basin_gl_storvol @@ -3061,13 +3106,13 @@ SUBROUTINE glacr_restart(In_out) WRITE ( Restart_outunit ) Glacr_tag WRITE ( Restart_outunit ) Delta_volyr WRITE ( Restart_outunit ) Hru_mb_yrcumul - WRITE ( Restart_outunit ) Hru_slope_ts + WRITE ( Restart_outunit ) Hru_slope_ts, Hru_glthick ELSE READ ( Restart_inunit ) module_name CALL check_restart(MODNAME, module_name) READ ( Restart_inunit ) Nhrugl, Basin_gl_top_melt, Gl_mb_yrcumul READ ( Restart_inunit ) Gl_mb_cumul, Glnet_ar_delta, Gl_mbc_yrend - READ ( Restart_inunit ) Basin_gl_top_gain + READ ( Restart_inunit ) Basin_gl_top_gain, Basin_glthick READ ( Restart_inunit ) Basin_gl_area, Gl_area, Basin_gl_ice_melt READ ( Restart_inunit ) Hru_glres_melt, Basin_gl_storstart READ ( Restart_inunit ) Gl_top_melt, Basin_gl_storage, Basin_gl_storvol @@ -3090,6 +3135,6 @@ SUBROUTINE glacr_restart(In_out) READ ( Restart_inunit ) Glacr_tag READ ( Restart_inunit ) Delta_volyr READ ( Restart_inunit ) Hru_mb_yrcumul - READ ( Restart_inunit ) Hru_slope_ts + READ ( Restart_inunit ) Hru_slope_ts, Hru_glthick ENDIF END SUBROUTINE glacr_restart diff --git a/prms/routing.f90 b/prms/routing.f90 index 9e700512..baf65715 100644 --- a/prms/routing.f90 +++ b/prms/routing.f90 @@ -148,7 +148,7 @@ INTEGER FUNCTION routingdecl() IF ( declparam( MODNAME, 'seg_length', 'nsegment', 'real', & & '1000.0', '0.001', '200000.0', & & 'Length of each segment', & - & 'Length of each segment, bounds based on CONUS', & + & 'Length of each segment including vertical drop', & & 'meters')/=0 ) CALL read_error(1, 'seg_length') ENDIF @@ -156,8 +156,8 @@ INTEGER FUNCTION routingdecl() ALLOCATE ( Seg_width(Nsegment) ) IF ( declparam(MODNAME, 'seg_width', 'nsegment', 'real', & & '15.0', '0.18', '40000.0', & - & 'Segment river width', & - & 'Segment river width, narrowest observed from Zimmerman 1967, Amazon biggest', & + & 'Segment bankfull river width', & + & 'Segment bankfull river width, narrowest observed from Zimmerman 1967, Amazon biggest', & & 'meters')/=0 ) CALL read_error(1, 'seg_width') ENDIF diff --git a/prms/snowcomp.f90 b/prms/snowcomp.f90 index f3184967..19b230dd 100644 --- a/prms/snowcomp.f90 +++ b/prms/snowcomp.f90 @@ -153,118 +153,118 @@ INTEGER FUNCTION snodecl() ! Glacier declares IF ( Glacier_flag==1 .OR. Model==99 ) THEN + ALLOCATE ( Ann_tempc(Nhru) ) + IF ( declvar(MODNAME, 'ann_tempc', 'nhru', Nhru, 'real', & + & 'Current average year air temperature over HRU', & + & 'degrees Celsius', Ann_tempc)/=0 ) CALL read_error(3, 'ann_tempc') + IF ( declvar(MODNAME, 'yrdays5', 'one', 1, 'integer', & & 'Number of days since last 5 year mark', & & 'none', Yrdays5)/=0 ) CALL read_error(3, 'yrdays5') ALLOCATE ( Glacr_freeh2o_capm(Nhru) ) IF ( declvar(MODNAME, 'glacr_freeh2o_capm', 'nhru', Nhru, 'real', & - & 'Free-water holding capacity of glacier ice, changes to 0 if active layer melts', & + & 'Free-water holding capacity of glacier or glacierette ice, changes to 0 if active layer melts', & & 'decimal fraction', Glacr_freeh2o_capm)/=0 ) CALL read_error(3, 'glacr_freeh2o_capm') ALLOCATE ( Glacrb_melt(Nhru) ) IF ( declvar(MODNAME, 'glacrb_melt', 'nhru', Nhru, 'real', & - 'Glacier basal melt, goes to soil', & + 'Glacier or glacierette basal melt, goes to soil', & 'inches/day', Glacrb_melt)/=0 ) CALL read_error(3, 'glacrb_melt') - ALLOCATE ( Ann_tempc(Nhru) ) - IF ( declvar(MODNAME, 'ann_tempc', 'nhru', Nhru, 'real', & - & 'Current average year air temperature overs HRU', & - & 'degrees Celsius', Ann_tempc)/=0 ) CALL read_error(3, 'ann_tempc') - ALLOCATE ( Glacr_air_5avtemp(Nhru) ) IF ( declvar(MODNAME, 'glacr_air_5avtemp', 'nhru', Nhru, 'real', & - & 'Current 5-yr average summer (June July Aug) air temperature over glacier or glrette HRU', & + & 'Current 5-yr average summer (June July Aug) air temperature over glacier or glacierette HRU', & & 'degrees Celsius', Glacr_air_5avtemp)/=0 ) CALL read_error(3, 'glacr_air_5avtemp') ALLOCATE ( Glacr_air_5avtemp1(Nhru) ) IF ( declvar(MODNAME, 'glacr_air_5avtemp1', 'nhru', Nhru, 'real', & - & 'First 5-yr average summer temperature over glacier or glrette HRU', & + & 'First 5-yr average summer temperature over glacier or glacierette HRU', & & 'degrees Celsius', Glacr_air_5avtemp1)/=0 ) CALL read_error(3, 'glacr_air_5avtemp1') ALLOCATE ( Glacr_air_deltemp(Nhru) ) IF ( declvar(MODNAME, 'glacr_air_deltemp', 'nhru', Nhru, 'real', & - & 'Change in 5-yr average air temperature over glacier or glrette HRU from first', & + & 'Change in 5-yr average air temperature over glacier or glacierette HRU from first', & & 'degrees Celsius', Glacr_air_deltemp)/=0 ) CALL read_error(3, 'glacr_air_deltemp') ALLOCATE ( Glacr_5avsnow(Nhru) ) IF ( declvar(MODNAME, 'glacr_5avsnow', 'nhru', Nhru, 'real', & - & 'Current 5-yr average snow over glacier or glrette HRU', & + & 'Current 5-yr average snow over glacier or glacierette HRU', & & 'inches/yr', Glacr_5avsnow)/=0 ) CALL read_error(3, 'glacr_5avsnow') ALLOCATE ( Glacr_5avsnow1(Nhru) ) IF ( declvar(MODNAME, 'glacr_5avsnow1', 'nhru', Nhru, 'real', & - & 'First 5-yr average snow over glacier or glrette HRU', & + & 'First 5-yr average snow over glacier or glacierette HRU', & & 'inches/yr', Glacr_5avsnow1)/=0 ) CALL read_error(3, 'glacr_5avsnow1') ALLOCATE ( Glacr_delsnow(Nhru) ) IF ( declvar(MODNAME, 'glacr_delsnow', 'nhru', Nhru, 'real', & - & 'Change in 5-yr average snow over glacier or glrette HRU from first', & + & 'Change in 5-yr average snow over glacier or glacierette HRU from first', & & 'inches/yr', Glacr_delsnow)/=0 ) CALL read_error(3, 'glacr_delsnow') ALLOCATE ( Glacr_pk_temp(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_temp', 'nhru', Nhru, 'real', & - & 'Temperature of the glacier on each HRU', & + & 'Temperature of the glacier or glacierette on each HRU', & & 'degrees Celsius', Glacr_pk_temp)/=0 ) CALL read_error(3, 'glacr_pk_temp') ALLOCATE ( Glacr_pk_def(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_def', 'nhru', Nhru, 'real', & - & 'Heat deficit, amount of heat necessary to make the glacier snowpack isothermal at 0 degrees Celsius', & + & 'Heat deficit, amount of heat necessary to make the glacier or or glacierette snowpack isothermal at 0 degrees Celsius', & & 'Langleys', Glacr_pk_def)/=0 ) CALL read_error(3, 'glacr_pk_def') ALLOCATE ( Glacr_pk_den(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_den', 'nhru', Nhru, 'real', & - & 'Density of the icepack on each glacier HRU, hard-coded to equal 0.917', & + & 'Density of the icepack on each glacier or glacierette HRU, hard-coded to equal 0.917', & & 'gm/cm3', Glacr_pk_den)/=0 ) CALL read_error(3, 'glacr_pk_den') ALLOCATE ( Glacr_albedo(Nhru) ) IF ( declvar(MODNAME, 'glacr_albedo', 'nhru', Nhru, 'real', & - & 'Ice surface albedo or the fraction of radiation reflected from the icepack surface for each glacier HRU', & + & 'Ice surface albedo or the fraction of radiation reflected from the icepack surface for each glacier or glacierette HRU', & & 'decimal fraction', Glacr_albedo)/=0 ) CALL read_error(3, 'glacr_albedo') ALLOCATE ( Glacr_evap(Nhru) ) IF ( declvar(MODNAME, 'glacr_evap', 'nhru', Nhru, 'real', & - & 'Evaporation and sublimation from icepack on each glacier HRU', & + & 'Evaporation and sublimation from icepack on each glacier or glacierette HRU', & & 'inches', Glacr_evap)/=0 ) CALL read_error(3, 'glacr_evap') ALLOCATE ( Glacrmelt(Nhru) ) IF ( declvar(MODNAME, 'glacrmelt', 'nhru', Nhru, 'real', & - & 'Melt from icepack on each glacier HRU, includes rain water that does not absorb', & + & 'Melt from icepack on each glacier or glacierette HRU, includes rain water that does not absorb', & & 'inches', Glacrmelt)/=0 ) CALL read_error(3, 'glacrmelt') ALLOCATE ( Glacr_pkwater_equiv(Nhru) ) IF ( declvar(MODNAME, 'glacr_pkwater_equiv', 'nhru', Nhru, 'double', & - & 'Icepack water equivalent on each glacier HRU', & + & 'Icepack water equivalent on each glacier or glacierette HRU', & & 'inches', Glacr_pkwater_equiv)/=0 ) CALL read_error(3, 'glacr_pkwater_equiv') ALLOCATE ( Glacr_pkwater_ante(Nhru) ) IF ( declvar(MODNAME, 'glacr_pkwater_ante', 'nhru', Nhru, 'double', & - & 'Antecedent icepack water equivalent on each glacier HRU', & + & 'Antecedent icepack water equivalent on each glacier or glacierette HRU', & & 'inches', Glacr_pkwater_ante)/=0 ) CALL read_error(3, 'glacr_pkwater_ante') ALLOCATE ( Glacrcov_area(Nhru) ) IF ( declvar(MODNAME, 'glacrcov_area', 'nhru', Nhru, 'real', & - & 'Ice-covered area on each glacier HRU or HRU with glacierette at start of step', & + & 'Ice-covered area (no snowpack) on each glacier HRU or HRU with glacierette at start of step', & & 'decimal fraction', Glacrcov_area)/=0 ) CALL read_error(3, 'glacrcov_area') ALLOCATE ( Glacr_pk_ice(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_ice', 'nhru', Nhru, 'real', & - & 'Storage of frozen water in the icepack on each glacier HRU', & + & 'Storage of frozen water in the icepack on each glacier or glacierette HRU', & & 'inches', Glacr_pk_ice)/=0 ) CALL read_error(3, 'glacr_pk_ice') ALLOCATE ( Glacr_freeh2o(Nhru) ) IF ( declvar(MODNAME, 'glacr_freeh2o', 'nhru', Nhru, 'real', & - & 'Storage of free liquid water in the icepack on each glacier HRU', & + & 'Storage of free liquid water in the icepack on each glacier or glacierette HRU', & & 'inches', Glacr_freeh2o)/=0 ) CALL read_error(3, 'glacr_freeh2o') ALLOCATE ( Glacr_pk_depth(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_depth', 'nhru', Nhru, 'double', & - & 'Depth of icepack on each glacier HRU, make essentially infinite', & + & 'Depth of icepack on each glacier or glacierette HRU, make essentially infinite', & & 'inches', Glacr_pk_depth)/=0 ) CALL read_error(3, 'glacr_pk_depth') ALLOCATE ( Glacr_pss(Nhru) ) IF ( declvar(MODNAME, 'glacr_pss', 'nhru', Nhru, 'double', & - & 'Previous glacier pack water equivalent plus new ice', & + & 'Previous glacier or glacierette pack water equivalent plus new ice', & & 'inches', Glacr_pss)/=0 ) CALL read_error(3, 'glacr_pss') ALLOCATE ( Glacr_pst(Nhru) ) @@ -273,41 +273,9 @@ INTEGER FUNCTION snodecl() & 'inches', Glacr_pst)/=0 ) CALL read_error(3, 'glacr_pst') IF ( declvar(MODNAME, 'basin_snowicecov', 'one', 1, 'double', & - & 'Basin area-weighted average snow and glacier and glrette covered area', & + & 'Basin area-weighted average snow and glacier and glacierette covered area for calibration to satellites', & & 'decimal fraction', Basin_snowicecov)/=0 ) CALL read_error(3, 'basin_snowicecov') - ALLOCATE ( Glacr_freeh2o_cap(Nhru) ) - IF ( declparam(MODNAME, 'glacr_freeh2o_cap', 'nhru', 'real', & - & '0.002', '0.0', '0.01', & - & 'Free-water holding capacity of glacier ice', & - & 'Free-water holding capacity of glacier ice expressed as a' // & - & ' decimal fraction of the frozen water content of the glacier ice (glacr_pk_ice)', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'glacr_freeh2o_cap') - - ALLOCATE ( Glacr_layer(Nhru) ) - IF ( declparam(MODNAME, 'glacr_layer', 'nhru', 'real', & - & '3.94', '0.0', '590.6', & - & 'Active layer on glacier', & - & 'Active layer is 0 to 15 m (590.6 inches) thick at start of year, when' // & - & ' melts will set daily glacr_pk_temp to 0', & - & 'inches')/=0 ) CALL read_error(1, 'glacr_layer') - - IF ( Init_vars_from_file==0 ) THEN - ALLOCATE ( Glacier_frac_init(Nhru) ) - IF ( declparam(MODNAME, 'glacier_frac_init', 'nhru', 'real', & - & '0.0', '0.0', '1.0', & - & 'Inital fraction of glaciation (0=none; 1=100%)', & - & 'Inital fraction of glaciation (0=none; 1=100%)', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'glacier_frac_init') - - ALLOCATE ( Glrette_frac_init(Nhru) ) - IF ( declparam(MODNAME, 'glrette_frac_init', 'nhru', 'real', & - & '0.0', '0.0', '1.0', & - & 'Initial fraction of glacierette (too small for glacier dynamics)', & - & 'Initial fraction of glacierette (too small for glacier dynamics)', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'glrette_frac_init') - - ENDIF ENDIF IF ( declvar(MODNAME, 'basin_snowdepth', 'one', 1, 'double', & @@ -506,6 +474,39 @@ INTEGER FUNCTION snodecl() & 'Ice albedo 300 meters below ELA', & & 'Ice albedo 300 meters below ELA', & & 'decimal fraction')/=0 ) CALL read_error(1, 'albedo_ice') + + ALLOCATE ( Glacr_freeh2o_cap(Nhru) ) + IF ( declparam(MODNAME, 'glacr_freeh2o_cap', 'nhru', 'real', & + & '0.002', '0.0', '0.01', & + & 'Free-water holding capacity of glacier ice', & + & 'Free-water holding capacity of glacier ice expressed as a' // & + & ' decimal fraction of the frozen water content of the glacier ice (glacr_pk_ice)', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'glacr_freeh2o_cap') + + ALLOCATE ( Glacr_layer(Nhru) ) + IF ( declparam(MODNAME, 'glacr_layer', 'nhru', 'real', & + & '3.94', '0.0', '590.6', & + & 'Active layer on glacier', & + & 'Active layer is 0 to 15 m (590.6 inches) thick at start of year, when' // & + & ' melts will set daily glacr_pk_temp to 0', & + & 'inches')/=0 ) CALL read_error(1, 'glacr_layer') + + IF ( Init_vars_from_file==0 ) THEN + ALLOCATE ( Glacier_frac_init(Nhru) ) + IF ( declparam(MODNAME, 'glacier_frac_init', 'nhru', 'real', & + & '0.0', '0.0', '1.0', & + & 'Inital fraction of glaciation (0=none; 1=100%) in glacier-capable HRU', & + & 'Inital fraction of glaciation (0=none; 1=100%) in glacier-capable HRU', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'glacier_frac_init') + + ALLOCATE ( Glrette_frac_init(Nhru) ) + IF ( declparam(MODNAME, 'glrette_frac_init', 'nhru', 'real', & + & '0.0', '0.0', '1.0', & + & 'Initial fraction of glacierette (too small for glacier dynamics)', & + & 'Initial fraction of glacierette (too small for glacier dynamics)', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'glrette_frac_init') + + ENDIF ENDIF IF ( declparam(MODNAME, 'den_init', 'one', 'real', & diff --git a/prms/stream_temp.f90 b/prms/stream_temp.f90 index fc30a24f..a62e6a68 100644 --- a/prms/stream_temp.f90 +++ b/prms/stream_temp.f90 @@ -20,11 +20,13 @@ MODULE PRMS_STRMTEMP REAL, SAVE, ALLOCATABLE :: gw_sum(:), ss_sum(:) REAL, SAVE, ALLOCATABLE :: gw_silo(:,:), ss_silo(:,:) REAL, SAVE, ALLOCATABLE :: hru_area_sum(:) + INTEGER, SAVE, ALLOCATABLE :: upstream_count(:) + INTEGER, SAVE, ALLOCATABLE :: upstream_idx(:,:) INTEGER, SAVE :: gw_index, ss_index ! Declared Variables REAL, SAVE, ALLOCATABLE :: Seg_tave_water(:), seg_tave_upstream(:), Seg_daylight(:) - REAL, SAVE, ALLOCATABLE :: Seg_humid(:), Seg_width_flow(:), Seg_ccov(:), seg_shade(:) + REAL, SAVE, ALLOCATABLE :: Seg_humid(:), Seg_width(:), Seg_ccov(:), seg_shade(:) REAL, SAVE, ALLOCATABLE :: Seg_tave_air(:), Seg_melt(:), Seg_rain(:) DOUBLE PRECISION, ALLOCATABLE :: Seg_potet(:) ! Segment Parameters @@ -41,7 +43,7 @@ MODULE PRMS_STRMTEMP INTEGER, SAVE :: Spring_jday, Summer_jday, Autumn_jday, Winter_jday ! Shade Parameters needed if stream_temp_shade_flag = 2 REAL, SAVE, ALLOCATABLE :: Segshade_sum(:), Segshade_win(:) - REAL, SAVE:: Albedo_str, Melt_temp + REAL, SAVE:: Albedo, Melt_temp ! INTEGER, SAVE :: Shadeflg, now using stream_temp_shade_flag INTEGER, SAVE, ALLOCATABLE :: Ss_tau(:), Gw_tau(:) ! Control parameters @@ -104,10 +106,10 @@ INTEGER FUNCTION stream_temp_decl() IF ( control_integer(Stream_temp_shade_flag, 'stream_temp_shade_flag')/=0 ) Stream_temp_shade_flag = 0 ! Declared Variables - ALLOCATE ( Seg_width_flow(Nsegment) ) - IF ( declvar( MODNAME, 'seg_width_flow', 'nsegment', Nsegment, 'real', & - & 'Width of each segment, flow-dependent', & - & 'meters', Seg_width_flow)/=0 ) CALL read_error(3, 'seg_width_flow') + ALLOCATE ( Seg_width(Nsegment) ) + IF ( declvar( MODNAME, 'seg_width', 'nsegment', Nsegment, 'real', & + & 'Width of each segment', & + & 'meters', Seg_width)/=0 ) CALL read_error(3, 'seg_width') ALLOCATE (Seg_tave_water(Nsegment) ) ! previous ?? IF ( declvar( MODNAME, 'seg_tave_water', 'nsegment', Nsegment, 'real', & @@ -187,11 +189,11 @@ INTEGER FUNCTION stream_temp_decl() ALLOCATE (gw_silo(nsegment,365), ss_silo(nsegment,365)) ALLOCATE (hru_area_sum(nsegment)) - IF ( declparam( MODNAME, 'albedo_str', 'one', 'real', & + IF ( declparam( MODNAME, 'albedo', 'one', 'real', & & '0.10', '0.0', '1.0', & & 'Short-wave solar radiation reflected by streams', & & 'Short-wave solar radiation reflected by streams', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'albedo_str') + & 'decimal fraction')/=0 ) CALL read_error(1, 'albedo') ALLOCATE(lat_temp_adj(Nsegment,12)) IF ( declparam( MODNAME, 'lat_temp_adj', 'nsegment,nmonths', 'real', & @@ -385,13 +387,13 @@ END FUNCTION stream_temp_decl !*********************************************************************** INTEGER FUNCTION stream_temp_init() USE PRMS_STRMTEMP - USE PRMS_MODULE, ONLY: Nsegment, Init_vars_from_file, Strmtemp_humidity_flag + USE PRMS_MODULE, ONLY: Nsegment, Init_vars_from_file, Inputerror_flag, Strmtemp_humidity_flag USE PRMS_BASIN, ONLY: Active_hrus, Hru_route_order, NEARZERO USE PRMS_OBS, ONLY: Nhumid USE PRMS_ROUTING, ONLY: Hru_segment, Tosegment, Segment_order, Segment_up IMPLICIT NONE ! Functions - INTRINSIC :: COS, SIN, ABS, SIGN, ASIN + INTRINSIC :: COS, SIN, ABS, SIGN, ASIN, maxval INTEGER, EXTERNAL :: getparam REAL, EXTERNAL :: solalt EXTERNAL :: read_error, checkdim_param_limits @@ -401,7 +403,7 @@ INTEGER FUNCTION stream_temp_init() !*********************************************************************** stream_temp_init = 0 - IF ( getparam( MODNAME, 'albedo_str', 1, 'real', Albedo_str)/=0 ) CALL read_error(2, 'albedo_str') + IF ( getparam( MODNAME, 'albedo', 1, 'real', Albedo)/=0 ) CALL read_error(2, 'albedo') IF ( getparam( MODNAME, 'lat_temp_adj', Nsegment*12, 'real', lat_temp_adj)/=0 ) CALL read_error(2, 'lat_temp_adj') IF (getparam(MODNAME, 'seg_lat', Nsegment, 'real', Seg_lat)/=0 ) CALL read_error(2, 'seg_lat') @@ -453,7 +455,7 @@ INTEGER FUNCTION stream_temp_init() seg_tave_upstream = 0.0 Seg_potet = 0.0D0 Seg_humid = 0.0 - Seg_width_flow = 0.0 + Seg_width = 0.0 Seg_ccov = 0.0 Seg_tave_air = 0.0 seg_tave_gw = 0.0 @@ -493,6 +495,12 @@ INTEGER FUNCTION stream_temp_init() Seg_hru_count(i) = Seg_hru_count(i) + 1 ENDDO +! exit if there are any segments that are too short + IF ( ierr==1 ) THEN + Inputerror_flag = ierr + RETURN + ENDIF + Seg_close = Segment_up ! assign upstream values DO j = 1, Nsegment ! set values based on routing order for segments without associated HRUs i = Segment_order(j) @@ -627,6 +635,42 @@ INTEGER FUNCTION stream_temp_init() endif enddo enddo + +! For each segment, figure out how many upstream segments. + ALLOCATE(upstream_count(Nsegment)) + upstream_count = 0 + do i = 1, nsegment + do j = 1, nsegment + if (tosegment(j) .eq. i) then + upstream_count(i) = upstream_count(i) + 1 + endif + end do + end do + +! For each segment, figure out the upstream segments. These will be looped over to determine inflows and temps to each segment + ALLOCATE(upstream_idx(Nsegment, maxval(upstream_count))) + upstream_idx = 0 + upstream_count = 0 + do i = 1, nsegment + do j = 1, nsegment + if (tosegment(j) .eq. i) then + upstream_count(i) = upstream_count(i) + 1 + upstream_idx(i,upstream_count(i)) = j + endif + end do + end do + +! do i = 1, nsegment +! write(*, fmt="(1x,a,i0)", advance="no") "segment #", i +! write(*, fmt="(1x,a,i0)", advance="no") " ", upstream_count(i) +! do j = 1, upstream_count(i) +! write(*, fmt="(1x,a,i0)", advance="no") " ", upstream_idx(i,j) +! end do +! write(*, fmt="(1x,a)",advance="yes") " done" +! end do + + + END FUNCTION stream_temp_init @@ -642,7 +686,7 @@ INTEGER FUNCTION stream_temp_run() USE PRMS_CLIMATE_HRU, ONLY: Humidity_hru USE PRMS_FLOWVARS, ONLY: Seg_outflow USE PRMS_SNOW, ONLY: Snowmelt - USE PRMS_ROUTING, ONLY: Hru_segment, Tosegment, Segment_order, Seginc_swrad, Seg_length + USE PRMS_ROUTING, ONLY: Hru_segment, Segment_order, Seginc_swrad, Seg_length USE PRMS_OBS, ONLY: Humidity USE PRMS_SET_TIME, ONLY: Nowyear, Nowmonth, Nowday, Jday USE PRMS_SOLTAB, ONLY: Soltab_potsw, Hru_cossl @@ -654,7 +698,7 @@ INTEGER FUNCTION stream_temp_run() EXTERNAL :: equilb, lat_inflow, shday ! Local Variables REAL :: harea, svi, fs - INTEGER :: i, j, k, iseg + INTEGER :: i, j, k, kk, iseg REAL :: te, ak1, ak2, ccov DOUBLE PRECISION :: qlat REAL :: t_o, up_temp @@ -805,15 +849,14 @@ INTEGER FUNCTION stream_temp_run() ! Find upstream intitial inflow temperature for segment i ! i is the current segment -! k is the upstream segment +! kk is the upstream segment fs = 0.0 up_temp = 0.0 - DO k = 1, Nsegment - IF ( Tosegment(k)==i ) THEN - if (Seg_tave_water(k) > -1.0) then - up_temp = up_temp + (Seg_tave_water(k) * SNGL(Seg_outflow(k))) - fs = fs + SNGL(Seg_outflow(k)) - endif + do k = 1, upstream_count(i) + kk = upstream_idx(i,k) + if (Seg_tave_water(kk) > -1.0) then + up_temp = up_temp + (Seg_tave_water(kk) * SNGL(Seg_outflow(kk))) + fs = fs + SNGL(Seg_outflow(kk)) ENDIF ENDDO @@ -833,9 +876,9 @@ INTEGER FUNCTION stream_temp_run() ! Compute flow-dependent water-in-segment width value if (seg_outflow(i) > NEARZERO) then - Seg_width_flow(i) = width_alpha(i) * sngl(Seg_outflow(i)) ** width_m(i) + Seg_width(i) = width_alpha(i) * sngl(Seg_outflow(i)) ** width_m(i) else - Seg_width_flow(i) = 0.0 + Seg_width(i) = 0.0 if (Seg_tave_water(i) > -99.0) then ! This segment has upstream HRUs somewhere, but the current day's flow is zero Seg_tave_water(i) = -98.9 @@ -921,7 +964,7 @@ INTEGER FUNCTION stream_temp_run() write(*,*) "this is the place: t_o = ", t_o, " ted = ", te, " seg_id = ", i write(*,*) " seg_tave_upstream = ", seg_tave_upstream(i), " fs = ", fs, & & " qlat = ", qlat, " seg_tave_lat = ", seg_tave_lat(i), " lat_temp_adj = ", lat_temp_adj(i,Nowmonth) - write(*,*) " width = ", Seg_width_flow(i), Nowyear, Nowmonth, Nowday + write(*,*) " width = ", Seg_width(i), Nowyear, Nowmonth, Nowday continue exit endif @@ -936,8 +979,8 @@ INTEGER FUNCTION stream_temp_run() CALL equilb(te, ak1, ak2, seg_shade(i), svi, i, t_o) ! Compute the daily mean water temperature - ! In: t_o, qlat, seg_tave_lat(i), te, ak1, ak2, i, seg_width_flow, seg_length/1000 (km) - Seg_tave_water(i) = twavg(fs, t_o, qlat, seg_tave_lat(i), te, ak1, ak2, seg_width_flow(i), seg_length(i)/1000.0) + ! In: t_o, qlat, seg_tave_lat(i), te, ak1, ak2, i, seg_width, seg_length + Seg_tave_water(i) = twavg(fs, t_o, qlat, seg_tave_lat(i), te, ak1, ak2, seg_width(i), seg_length(i)/1000.0) else ! bad t_o value @@ -1022,7 +1065,7 @@ REAL FUNCTION twavg(qup, T0, Qlat, Tl_avg, Te, Ak1, Ak2, width, length) Ql = SNGL( Qlat ) ! This is confused logic coment out here and compute the terms as needed below -! b = (Ql / Seg_length/1000) + ((Ak1 * Seg_width_flow) / 4182.0E03) +! b = (Ql / Seg_length/1000.0) + ((Ak1 * Seg_width) / 4182.0E03) ! IF ( b < NEARZERO ) b = NEARZERO ! rsr, don't know what value this should be to avoid divide by 0 ! r = 1.0 + (Ql / q_init) ! IF ( r < NEARZERO ) r = NEARZERO @@ -1097,11 +1140,11 @@ SUBROUTINE equilb (Ted, Ak1d, Ak2d, Sh, Svi, Seg_id, t_o) ! 1. DETERMINE THE AVERAGE DAILY EQUILIBRIUM WATER TEMPERATURE PARAMETERS ! 2. DETERMINE THE MAXIMUM DAILY EQUILIBRIUM WATER TEMPERATURE PARAMETERS - USE PRMS_STRMTEMP, ONLY: ZERO_C, Seg_width_flow, Seg_humid, Press, MPS_CONVERT, & - & Seg_ccov, Seg_potet, Albedo_str, seg_tave_gw + USE PRMS_STRMTEMP, ONLY: ZERO_C, Seg_width, Seg_humid, Press, MPS_CONVERT, & + & Seg_ccov, Seg_potet, Albedo, seg_tave_gw USE PRMS_BASIN, ONLY: NEARZERO, CFS2CMS_CONV USE PRMS_FLOWVARS, ONLY: Seg_inflow - USE PRMS_ROUTING, ONLY: Seginc_swrad, Seg_slope + USE PRMS_ROUTING, ONLY: Seg_slope, Seginc_swrad IMPLICIT NONE ! Functions INTRINSIC EXP, SQRT, ABS, SNGL, DBLE @@ -1154,8 +1197,8 @@ SUBROUTINE equilb (Ted, Ak1d, Ak2d, Sh, Svi, Seg_id, t_o) & * (1.0 + (0.17*(Seg_ccov(Seg_id)**2)))) ) * (taabs**4) ! hf is heat from stream friction. See eqn. 14. q_init is in CMS - hf = 9805.0 * (q_init/Seg_width_flow(Seg_id)) * Seg_slope(Seg_id) - hs = (1.0 - sh) * sw_power * (1.0 - Albedo_str) + hf = 9805.0 * (q_init/Seg_width(Seg_id)) * Seg_slope(Seg_id) + hs = (1.0 - sh) * sw_power * (1.0 - Albedo) hv = 5.24D-8 * DBLE(Svi) * (taabs**4) ! Stefan-Boltzmann constant = 5.670373D-08; emissivity of water = 0.9526, times each other: 5.4016D-08 @@ -1311,7 +1354,7 @@ SUBROUTINE shday(Seg_id, Shade, Svi) ! Vow = OFFSET, WEST SIDE VEGETATION ! USE PRMS_SET_TIME, ONLY: Jday - USE PRMS_STRMTEMP, ONLY: Azrh, Alte, Altw, Seg_daylight, Seg_width_flow, & + USE PRMS_STRMTEMP, ONLY: Azrh, Alte, Altw, Seg_daylight, Seg_width, & & PI, HALF_PI, Cos_seg_lat, Sin_seg_lat, Cos_lat_decl, Horizontal_hour_angle, & & Level_sunset_azimuth, Max_solar_altitude, Sin_alrs, Sin_declination, Sin_lat_decl, Total_shade USE PRMS_BASIN, ONLY: CFS2CMS_CONV @@ -1382,7 +1425,7 @@ SUBROUTINE shday(Seg_id, Shade, Svi) ! azss = azso hrss = hrso sti = 0.0 - Svi = (rprnvg(hrsr, hrrs, hrss, sino, coso, sin_d, cosod, sinod, Seg_id)) / (Seg_width_flow(Seg_id) * totsh) + Svi = (rprnvg(hrsr, hrrs, hrss, sino, coso, sin_d, cosod, sinod, Seg_id)) / (Seg_width(Seg_id) * totsh) ELSE ! INITIALIZE SHADE VALUES @@ -1458,7 +1501,7 @@ SUBROUTINE shday(Seg_id, Shade, Svi) Seg_daylight(Seg_id) = (hrss - hrsr) * RADTOHOUR sti = 1.0 - ((((hrss - hrsr) * sinod) + ((SIN(hrss) - SIN(hrsr)) * cosod)) / (totsh)) - Svi = ((rprnvg(hrsr, hrrh, hrss, sino, coso, sin_d, cosod, sinod, Seg_id)) / (Seg_width_flow(Seg_id)*totsh)) + Svi = ((rprnvg(hrsr, hrrh, hrss, sino, coso, sin_d, cosod, sinod, Seg_id)) / (Seg_width(Seg_id)*totsh)) ! ! END SUNRISE/SUNSET CALCULATION ENDIF @@ -1647,7 +1690,7 @@ REAL FUNCTION rprnvg (Hrsr, Hrrs, Hrss, Sino, Coso, Sin_d, Cosod, Sinod, Seg_id) ! THIS SUBPROGRAM IS TO COMPUTE THE RIPARIAN VEGETATION SHADE ! SEGMENT BETWEEN THE TWO HOUR ANGLES HRSR & HRSS. ! - USE PRMS_STRMTEMP, ONLY: Azrh, Vce, Vdemx, Vhe, Voe, Vcw, Vdwmx, Vhw, Vow, Seg_width_flow, & + USE PRMS_STRMTEMP, ONLY: Azrh, Vce, Vdemx, Vhe, Voe, Vcw, Vdwmx, Vhw, Vow, Seg_width, & & Vdemn, Vdwmn, HALF_PI USE PRMS_BASIN, ONLY: NEARZERO USE PRMS_SET_TIME, ONLY: Summer_flag @@ -1710,7 +1753,7 @@ REAL FUNCTION rprnvg (Hrsr, Hrrs, Hrss, Sino, Coso, Sin_d, Cosod, Sinod, Seg_id) ! DETERMINE AMOUNT OF STREAM WIDTH SHADED bs = ((Vhe(Seg_id) * (cosals/sinals)) * ABS(SIN(azs-Azrh(Seg_id)))) + vco IF ( bs < 0.0 ) bs = 0.0 - IF ( bs > Seg_width_flow(Seg_id) ) bs = Seg_width_flow(Seg_id) + IF ( bs > Seg_width(Seg_id) ) bs = Seg_width(Seg_id) ! INCREMENT SUNRISE SIDE VEGETATIVE SHADE IF ( Summer_flag == 1 ) THEN ! put back spring and autumn svri = svri + SNGL(Vdemx(Seg_id) * bs * sinals * Weight(n)) @@ -1752,7 +1795,7 @@ REAL FUNCTION rprnvg (Hrsr, Hrrs, Hrss, Sino, Coso, Sin_d, Cosod, Sinod, Seg_id) ! DETERMINE AMOUNT OF STREAM WIDTH SHADED bs = ((Vhw(Seg_id) * (cosals/sinals)) * ABS(SIN(azs-Azrh(Seg_id)))) + vco IF ( bs < 0.0 ) bs = 0.0 - IF ( bs > Seg_width_flow(Seg_id) ) bs = Seg_width_flow(Seg_id) + IF ( bs > Seg_width(Seg_id) ) bs = Seg_width(Seg_id) ! INCREMENT SUNSET SIDE VEGETATIVE SHADE IF ( Summer_flag == 1 ) THEN ! fix for seasons svsi = SNGL(svsi + (Vdwmx(Seg_id) * bs * sinals * Weight(n))) diff --git a/prmsRip/mizurouteRip.f90 b/prmsRip/mizurouteRip.f90 index f49c55d1..539f95a6 100644 --- a/prmsRip/mizurouteRip.f90 +++ b/prmsRip/mizurouteRip.f90 @@ -437,7 +437,7 @@ END FUNCTION mizuroute_init !*********************************************************************** INTEGER FUNCTION mizuroute_run() USE PRMS_MIZUROUTE - USE PRMS_MODULE, ONLY: Nsegment, Ripst_flag, Glacier_flag + USE PRMS_MODULE, ONLY: Nsegment, Ripst_flag, Glacier_flag, Frozen_flag USE PRMS_BASIN, ONLY: CFS2CMS_CONV, Basin_area_inv, METERS2FEET, Active_hrus, Hru_route_order, & & Basin_gl_cfs, Basin_gl_ice_cfs USE PRMS_FLOWVARS, ONLY: Basin_ssflow, Basin_cms, Basin_gwflow_cfs, Basin_ssflow_cfs, & @@ -455,7 +455,7 @@ INTEGER FUNCTION mizuroute_run() & Hru_segment, Seg_length, Basin_ripst_area, Basin_ripst_contrib, Basin_ripst_evap, & & Basin_ripst_vol, Bankst_seep_rate USE PRMS_GLACR, ONLY: Basin_gl_top_melt, Basin_gl_ice_melt - USE PRMS_SRUNOFF, ONLY: Basin_sroff + USE PRMS_SRUNOFF, ONLY: Basin_sroff, Frozen USE PRMS_GWFLOW, ONLY: Basin_gwflow ! mizuroute specific modules USE nrtype ! variable types, etc. @@ -634,10 +634,17 @@ INTEGER FUNCTION mizuroute_run() ENDDO DO j = 1, Active_hrus i = Hru_route_order(j) - IF ( Hru_segment(i)>0 .AND. (Bankfinite_hru(i)==0 .OR. Ripst_areafr_max(i)>0.0)) & - & CALL comp_bank_storage(i) + IF ( Hru_segment(i)>0 .AND. (Bankfinite_hru(i)==0 .OR. Ripst_areafr_max(i)>0.0)) THEN + IF (Frozen_flag==1) THEN + IF ( Frozen(i).ne.1 ) THEN + CALL comp_bank_storage(i) + ENDIF + ELSE + CALL comp_bank_storage(i) + ENDIF ! ******Compute the bank storage component ! transfers water between separate bank storage and stream depending on seepage + ENDIF ENDDO Basin_bankst_seep = Basin_bankst_seep*Basin_area_inv Basin_bankst_head = Basin_bankst_head*Basin_area_inv diff --git a/prmsRip/muskingumRip.f90 b/prmsRip/muskingumRip.f90 index 37c6bfb9..0ac99677 100644 --- a/prmsRip/muskingumRip.f90 +++ b/prmsRip/muskingumRip.f90 @@ -178,7 +178,7 @@ END FUNCTION muskingum_init !*********************************************************************** INTEGER FUNCTION muskingum_run() USE PRMS_MUSKINGUM - USE PRMS_MODULE, ONLY: Nsegment, Ripst_flag, Glacier_flag + USE PRMS_MODULE, ONLY: Nsegment, Ripst_flag, Glacier_flag, Frozen_flag USE PRMS_BASIN, ONLY: CFS2CMS_CONV, Basin_area_inv, Active_hrus, Hru_route_order, & & Basin_gl_cfs, Basin_gl_ice_cfs USE PRMS_FLOWVARS, ONLY: Basin_ssflow, Basin_cms, Basin_gwflow_cfs, Basin_ssflow_cfs, & @@ -196,7 +196,7 @@ INTEGER FUNCTION muskingum_run() & Hru_segment, Seg_length, Basin_ripst_area, Basin_ripst_contrib, Basin_ripst_evap, & & Basin_ripst_vol, Bankst_seep_rate USE PRMS_GLACR, ONLY: Basin_gl_top_melt, Basin_gl_ice_melt - USE PRMS_SRUNOFF, ONLY: Basin_sroff + USE PRMS_SRUNOFF, ONLY: Basin_sroff, Frozen USE PRMS_GWFLOW, ONLY: Basin_gwflow IMPLICIT NONE ! Functions @@ -333,10 +333,17 @@ INTEGER FUNCTION muskingum_run() ENDDO DO j = 1, Active_hrus i = Hru_route_order(j) - IF ( Hru_segment(i)>0 .AND. (Bankfinite_hru(i)==0 .OR. Ripst_areafr_max(i)>0.0)) & - & CALL comp_bank_storage(i) + IF ( Hru_segment(i)>0 .AND. (Bankfinite_hru(i)==0 .OR. Ripst_areafr_max(i)>0.0)) THEN + IF (Frozen_flag==1) THEN + IF ( Frozen(i).ne.1 ) THEN + CALL comp_bank_storage(i) + ENDIF + ELSE + CALL comp_bank_storage(i) + ENDIF ! ******Compute the bank storage component ! transfers water between separate bank storage and stream depending on seepage + ENDIF ENDDO Basin_bankst_seep = Basin_bankst_seep*Basin_area_inv Basin_bankst_head = Basin_bankst_head*Basin_area_inv @@ -344,7 +351,6 @@ INTEGER FUNCTION muskingum_run() DO i = 1, Nsegment Basin_bankst_seep_rate = Basin_bankst_seep_rate + Bankst_seep_rate(i) & & *Seg_length(i)/SUM(Seg_length) !m2/day per stream ft length - !print*, Seg_outflow(i)+Seg_bankflow(i),Seg_outflow(i),Seg_bankflow(i), i Seg_outflow(i) = Seg_outflow(i)+Seg_bankflow(i) IF (Seg_bankflow(i) < 0.0) THEN ! only could go negative because of bankflow if is negative IF (Seg_outflow(i) < 0.0) THEN ! took out more than streamflow, this could also be a water_use problem diff --git a/prmsRip/routingRip.f90 b/prmsRip/routingRip.f90 index 53f10b5f..84d84d42 100644 --- a/prmsRip/routingRip.f90 +++ b/prmsRip/routingRip.f90 @@ -29,10 +29,11 @@ MODULE PRMS_ROUTING REAL, SAVE, ALLOCATABLE :: Seg_depth(:), K_coef(:), X_coef(:), Mann_n(:), Seg_width(:), Segment_flow_init(:) REAL, SAVE, ALLOCATABLE :: Seg_length(:), Seg_slope(:) ! Declared Parameters for Overbank and bank Storage + INTEGER, SAVE :: Two REAL, SAVE, ALLOCATABLE :: Transmiss_seg(:), Ripst_areafr_max(:) REAL, SAVE :: Bank_height_fac ! Declared Parameters for Overbank Storage - REAL, SAVE, ALLOCATABLE :: Tr_ratio(:), Porosity_seg(:), Ripst_et_coef(:), Ripst_frac_init(:) + REAL, SAVE, ALLOCATABLE :: Tr_ratio(:), Porosity_seg(:), Ripst_frac_init(:) ! Declared Variables for Overbank Storage DOUBLE PRECISION, SAVE :: Basin_ripst_evap, Basin_ripst_contrib, Basin_ripst_vol, Basin_ripst_area DOUBLE PRECISION, SAVE, ALLOCATABLE :: Ripst_stor_hru(:), Ripst_vol(:), Seg_ripflow(:) @@ -44,7 +45,7 @@ MODULE PRMS_ROUTING DOUBLE PRECISION, SAVE :: Basin_bankst_head, Basin_bankst_seep_rate DOUBLE PRECISION, SAVE :: Basin_bankst_seep, Basin_bankst_vol, Basin_bankst_area REAL, SAVE, ALLOCATABLE :: Bankst_head(:), Bankst_seep_rate(:), Bankst_seep_hru(:) - REAL, SAVE, ALLOCATABLE :: Bankst_stor_hru(:), Bankst_head_pts(:) + REAL, SAVE, ALLOCATABLE :: Bankst_stor_hru(:), Bankst_head_pts(:,:) DOUBLE PRECISION, SAVE, ALLOCATABLE :: Stage_ante(:), Stage_ts(:), Seg_bankflow(:) END MODULE PRMS_ROUTING @@ -82,7 +83,7 @@ INTEGER FUNCTION routingdecl() & Ripst_flag, Stream_temp_flag, Init_vars_from_file IMPLICIT NONE ! Functions - INTEGER, EXTERNAL :: declparam, declvar + INTEGER, EXTERNAL :: declparam, declvar, declfix EXTERNAL read_error, print_module !*********************************************************************** routingdecl = 0 @@ -148,6 +149,10 @@ INTEGER FUNCTION routingdecl() ! 100 = user normal; 101 - 108 = not used; 109 sink (tosegment used by Lumen) IF ( Ripst_flag==1 .OR. Model==99 ) THEN + + IF ( declfix('two', 2, 2, 'Need for keeping bank storage head points')/=0 ) CALL read_error(7, 'two') + Two = 2 + ! Overbank storage variables IF ( declvar(MODNAME, 'basin_ripst_evap', 'one', 1, 'double', & & 'Basin area-weighted average evaporation from riparian overbank flow storage', & @@ -221,12 +226,12 @@ INTEGER FUNCTION routingdecl() ALLOCATE ( Seg_bankflow(Nsegment) ) IF ( declvar(MODNAME, 'seg_bankflow', 'nsegment', Nsegment, 'double', & - & 'Bank storage area contribution to streamflow can be negative if steam losing water', & + & 'Bank storage area contribution to streamflow, negative if steam losing water', & & 'cfs', Seg_bankflow)/=0 ) CALL read_error(3, 'seg_bankflow') - ALLOCATE ( Bankst_head_pts(Nhru) ) - IF ( declvar(MODNAME, 'bankst_head_pts', 'nhru', Nhru, 'real', & - & 'Head of bank storage above groundwater head: at half width away', & + ALLOCATE ( Bankst_head_pts(Nhru,Two) ) + IF ( declvar(MODNAME, 'bankst_head_pts', 'nhru,two', Nhru*Two, 'real', & + & 'Head of bank storage above groundwater head: at quarter width away and edge of riparian area', & & 'meters', Bankst_head_pts)/=0 ) CALL read_error(3, 'bankst_head_pts') ALLOCATE ( Stage_ante(Nsegment) ) @@ -241,12 +246,12 @@ INTEGER FUNCTION routingdecl() ALLOCATE ( Bankst_seep_hru(Nhru) ) IF ( declvar(MODNAME, 'bankst_seep_hru', 'nhru', Nhru, 'real', & - & 'HRU average seepage from bank storage to associated stream_segment for each HRU', & + & 'HRU average seepage from bank storage to associated stream segment for each HRU', & & 'inches', Bankst_seep_hru)/=0 ) CALL read_error(3, 'bankst_seep_hru') ALLOCATE ( Bankst_stor_hru(Nhru) ) IF ( declvar(MODNAME, 'bankst_stor_hru', 'nhru', Nhru, 'real', & - & 'HRU average bank storage for each HRU', & + & 'Average bank storage for each HRU', & & 'inches', Bankst_stor_hru)/=0 ) CALL read_error(3, 'bankst_stor_hru') ALLOCATE ( Bankst_seep_rate(Nsegment) ) @@ -277,7 +282,7 @@ INTEGER FUNCTION routingdecl() IF ( declparam( MODNAME, 'seg_length', 'nsegment', 'real', & & '1000.0', '0.001', '200000.0', & & 'Length of each segment', & - & 'Length of each segment, bounds based on CONUS', & + & 'Length of each segment including vertical drop', & & 'meters')/=0 ) CALL read_error(1, 'seg_length') ENDIF @@ -285,8 +290,8 @@ INTEGER FUNCTION routingdecl() ALLOCATE ( Seg_width(Nsegment) ) IF ( declparam(MODNAME, 'seg_width', 'nsegment', 'real', & & '15.0', '0.18', '40000.0', & - & 'Segment river width', & - & 'Segment river width, narrowest observed from Zimmerman 1967, Amazon biggest', & + & 'Segment bankfull river width', & + & 'Segment bankfull river width, narrowest observed from Zimmerman 1967, Amazon biggest', & & 'meters')/=0 ) CALL read_error(1, 'seg_width') ENDIF @@ -367,7 +372,7 @@ INTEGER FUNCTION routingdecl() ALLOCATE ( Ripst_frac_init(Nhru) ) IF ( declparam(MODNAME, 'ripst_frac_init', 'nhru', 'real', & & '0.5', '0.0', '1.0', & - & 'Fraction of maximum storage that contains water at the start of a simulation', & + & 'Fraction of maximum storage volume that contains water at the start of a simulation', & & 'Fraction of maximum riparian overbank flow storage that'// & & ' contains water at the start of a simulation', & & 'decimal fraction')/=0 ) CALL read_error(1, 'ripst_frac_init') @@ -385,7 +390,7 @@ INTEGER FUNCTION routingdecl() & 'decimal fraction')/=0 ) CALL read_error(1, 'ripst_areafr_max') IF ( declparam(MODNAME, 'bank_height_fac', 'one', 'real', & - & '20.0', '1.0', '100.0', & + & '20.0', '1.0', '1000.0', & & 'Factor multiplied to Seg_depth to give maximum height of banks', & & 'Factor multiplied to Seg_depth to give maximum height of banks for riparian overbank storage', & & 'none')/=0 ) CALL read_error(1, 'bank_height_fac') @@ -397,13 +402,6 @@ INTEGER FUNCTION routingdecl() & 'Porosity of soil around segment involved in riparian overbank flow storage', & & 'decimal fraction')/=0 ) CALL read_error(1, 'porosity_seg') - ALLOCATE ( Ripst_et_coef(Nhru) ) - IF ( declparam(MODNAME, 'ripst_et_coef', 'nhru', 'real', & - & '1.0', '0.0', '1.0', & - & 'Fraction of unsatisfied potential evapotranspiration to apply to riparian overbank flow storage', & - & 'Fraction of unsatisfied potential evapotranspiration to apply to riparian overbank flow storage', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'ripst_et_coef') - ALLOCATE ( Tr_ratio(Nhru) ) IF ( declparam(MODNAME, 'tr_ratio', 'nhru', 'real', & & '0.5', '0.0', '1.0', & @@ -542,11 +540,12 @@ INTEGER FUNCTION routinginit() USE PRMS_ROUTING USE PRMS_MODULE, ONLY: Nsegment, Nhru, Init_vars_from_file, Strmflow_flag, & & Water_use_flag, Segment_transferON_OFF, Inputerror_flag, Parameter_check_flag , & - & Ripst_flag, Stream_temp_flag !, Print_debug + & Ripst_flag, Stream_temp_flag, PRMS4_flag, Dprst_flag !, Print_debug USE PRMS_SET_TIME, ONLY: Timestep_seconds USE PRMS_BASIN, ONLY: FT2_PER_ACRE, DNEARZERO, Active_hrus, Hru_route_order, Hru_area_dble, NEARZERO, & - & Hru_area, FEET2METERS, CFS2CMS_CONV !, Active_area + & Hru_area, FEET2METERS, CFS2CMS_CONV, Dprst_open_flag !, Active_area USE PRMS_FLOWVARS, ONLY: Seg_outflow + USE PRMS_SRUNOFF, ONLY: Dprst_seep_rate_open, Sro_to_dprst_imperv, Sro_to_dprst_perv, Dprst_et_coef IMPLICIT NONE ! Functions INTRINSIC MOD, DBLE @@ -670,7 +669,7 @@ INTEGER FUNCTION routinginit() Seg_outflow(i) = Segment_flow_init(i) IF ( Ripst_flag==1 ) THEN flow = Seg_outflow(i)*CFS2CMS_CONV - Stage_ts(i) = (Mann_n(i)*flow/( Seg_width(i)*SQRT(Seg_slope(i)) ))**(5./3.) + Stage_ante(i) = (Mann_n(i)*flow/( Seg_width(i)*SQRT(Seg_slope(i)) ))**(3./5.) ENDIF ENDDO DEALLOCATE ( Segment_flow_init ) @@ -705,13 +704,26 @@ INTEGER FUNCTION routinginit() IF ( Ripst_flag==1 ) THEN IF ( getparam(MODNAME, 'ripst_areafr_max', Nhru, 'real', Ripst_areafr_max)/=0 ) CALL read_error(2, 'ripst_areafr_max') IF ( getparam(MODNAME, 'bank_height_fac', 1, 'real', Bank_height_fac)/=0 ) CALL read_error(2, 'bank_height_fac') - IF ( getparam(MODNAME, 'ripst_et_coef', Nhru, 'real', Ripst_et_coef)/=0 ) CALL read_error(2, 'ripst_et_coef') IF ( getparam(MODNAME, 'tr_ratio', Nhru, 'real', Tr_ratio)/=0 ) CALL read_error(2, 'tr_ratio') IF ( getparam(MODNAME, 'bankfinite_hru', Nhru, 'integer', Bankfinite_hru)/=0 ) CALL read_error(2, 'bankfinite_hru') ! might be able to calculate if want bankfinite_hru = 1 or 0 based on ripst_areafr_max and transmiss_seg IF ( getparam(MODNAME, 'transmiss_seg', Nsegment, 'real', Transmiss_seg)/=0 ) CALL read_error(2, 'transmiss_seg') IF ( getparam(MODNAME, 'specyield_seg', Nsegment, 'real', Specyield_seg)/=0 ) CALL read_error(2, 'specyield_seg') IF ( getparam(MODNAME, 'porosity_seg', Nsegment, 'real', Porosity_seg)/=0 ) CALL read_error(2, 'porosity_seg') + + IF (Dprst_flag/=1) THEN !didn't call these already then + IF ( PRMS4_flag==1 ) THEN + IF ( getparam(MODNAME, 'sro_to_dprst', Nhru, 'real', Sro_to_dprst_perv)/=0 ) CALL read_error(2, 'sro_to_dprst') + ELSE + IF ( getparam(MODNAME, 'sro_to_dprst_perv', Nhru, 'real', Sro_to_dprst_perv)/=0 ) CALL read_error(2, 'sro_to_dprst_perv') + ENDIF + IF ( getparam(MODNAME, 'sro_to_dprst_imperv', Nhru, 'real', Sro_to_dprst_imperv)/=0 ) CALL read_error(2, 'sro_to_dprst_imperv') + IF ( getparam(MODNAME, 'dprst_et_coef', Nhru, 'real', Dprst_et_coef)/=0 ) CALL read_error(2, 'dprst_et_coef') + ENDIF + IF (Dprst_open_flag/=1 .OR. Dprst_flag/=1) THEN + IF ( getparam(MODNAME, 'dprst_seep_rate_open', Nhru, 'real', Dprst_seep_rate_open)/=0 ) & + & CALL read_error(2, 'dprst_seep_rate_open') + ENDIF Seg_hru_num = 0 DO i = 1, Active_hrus IF ( Hru_segment(i)>0) THEN @@ -1138,28 +1150,31 @@ END SUBROUTINE init_the_swamp !*********************************************************************** SUBROUTINE drain_the_swamp(Ihru) USE PRMS_ROUTING, ONLY: Seg_width, Seg_depth, Seg_width, Hru_segment, Mann_n, & - & Tr_ratio, Ripst_vol_max, Ripst_et_coef, Ripst_evap_hru, Seg_length, & + & Tr_ratio, Ripst_vol_max, Ripst_evap_hru, Seg_length, & & Basin_ripst_vol, Basin_ripst_evap, Basin_ripst_contrib, Ripst_stor_hru, & & Ripst_frac, Ripst_vol, Ripst_area_max, Ripst_area, Seg_slope, & - & Seg_hru_num, Seg_ripflow, Ripst_depth, Basin_ripst_area, Bank_height_fac !, Transmiss_seg - USE PRMS_MODULE, ONLY: Frozen_flag + & Seg_hru_num, Seg_ripflow, Ripst_depth, Basin_ripst_area, Bank_height_fac, & + & Ripst_areafr_max !, Transmiss_seg + USE PRMS_MODULE, ONLY: Cascade_flag, Frozen_flag USE PRMS_BASIN, ONLY: NEARZERO, DNEARZERO, Hru_area, Hru_area_dble, FEET2METERS, & & FT2_PER_ACRE, CFS2CMS_CONV - USE PRMS_FLOWVARS, ONLY: Seg_outflow + USE PRMS_FLOWVARS, ONLY: Seg_outflow, Pkwater_equiv USE PRMS_CLIMATEVARS, ONLY: Potet USE PRMS_SET_TIME, ONLY: Timestep_seconds USE PRMS_SRUNOFF, ONLY: Hru_impervevap, Dprst_evap_hru, Frozen, Thaw_depth, Soil_depth, & - & Dprst_seep_rate_open - USE PRMS_INTCP, ONLY: Hru_intcpevap - USE PRMS_SNOW, ONLY: Snowcov_area, Snow_evap + & Dprst_seep_rate_open, Upslope_hortonian, Srp, Sri, Imperv_frac, Perv_frac, & + & Sro_to_dprst_imperv, Sro_to_dprst_perv, Dprst_et_coef + USE PRMS_INTCP, ONLY: Hru_intcpevap, Net_rain, Net_snow + USE PRMS_SNOW, ONLY: Snowmelt, Pptmix_nopack, Snowcov_area, Snow_evap IMPLICIT NONE ! Functions INTRINSIC EXP, LOG, MIN, DBLE, SNGL ! Arguments INTEGER, INTENT(IN) :: Ihru ! Local Variables + INTEGER :: seep_surface REAL :: ripst_avail_et, unsatisfied_et, ripst_evap, ripst_wid, thaw_frac - REAL :: inflow, inflow_in, max_depth + REAL :: inflow, inflow_in, tmp, ripst_sri, ripst_srp, max_depth DOUBLE PRECISION :: seep, ripst_grnd, poss, seep_in, ripst_contrib_hru !*********************************************************************** thaw_frac = 1.0 @@ -1170,37 +1185,109 @@ SUBROUTINE drain_the_swamp(Ihru) thaw_frac = Thaw_depth(Ihru)/Soil_depth(Ihru) ENDIF ENDIF + +! add the hortonian flow to the riparian storage volumes: + IF ( Cascade_flag>0 ) THEN + inflow = SNGL( Upslope_hortonian(Ihru) ) + ELSE + inflow = 0.0 + ENDIF + inflow_in = 0.0 + + IF ( Pptmix_nopack(Ihru)==1 ) inflow = inflow + Net_rain(Ihru) + +!******If precipitation on snowpack all water available to the surface is considered to be snowmelt +!******If there is no snowpack and no precip,then check for melt from last of snowpack. +!******If rain/snow mix with no antecedent snowpack, compute snowmelt portion of runoff. + + IF ( Snowmelt(Ihru)>0.0 ) THEN + inflow = inflow + Snowmelt(Ihru) + +!******There was no snowmelt but a snowpack may exist. If there is +!******no snowpack then check for rain on a snowfree HRU. + ELSEIF ( Pkwater_equiv(Ihru)0.0 ) THEN + inflow = inflow + Net_rain(Ihru) + ENDIF + ENDIF + + ! add any pervious surface runoff fraction to riparian areas, use depression storage factor + ripst_srp = 0.0 + ripst_sri = 0.0 + IF ( Srp>0.0 ) THEN + tmp = Srp*Perv_frac*Sro_to_dprst_perv(Ihru)*Hru_area(Ihru) + ripst_srp = tmp*Ripst_areafr_max(Ihru) ! acre-inches + ENDIF + IF ( Sri>0.0 ) THEN + tmp = Sri*Imperv_frac*Sro_to_dprst_imperv(Ihru)*Hru_area(Ihru) + ripst_sri = tmp*Ripst_areafr_max(Ihru) ! acre-inches + ENDIF + !It won't get deeper than this depth, should be close to Seg_depth but not accurate or Seg_width and other terms not accurate max_depth = Seg_depth(Hru_segment(Ihru))*Bank_height_fac ! amount possible in cfs given a river depth poss = Seg_width(Hru_segment(Ihru))*SQRT(Seg_slope(Hru_segment(Ihru)))* & - & max_depth**(3./5.)/ ( CFS2CMS_CONV*Mann_n(Hru_segment(Ihru)) ) + & max_depth**(5./3.)/ ( CFS2CMS_CONV*Mann_n(Hru_segment(Ihru)) ) !inflow is water over bank, remove from Seg_outflow(Hru_segment(Ihru)) and give half to ! each side of bank, in acre inches - inflow = 0.0 - inflow_in = 0.0 + ! in cfs, amount over amount possible, no inflow if everything frozen, and then no outflow either + Ripst_vol(Ihru) = 0.0 + Ripst_frac(Ihru) =0.0 IF (thaw_frac>0.0) THEN - IF ( poss < Seg_outflow(Hru_segment(Ihru))) inflow = SNGL(Seg_outflow(Hru_segment(Ihru)) - poss) + IF ( poss < Seg_outflow(Hru_segment(Ihru))) inflow_in = SNGL(Seg_outflow(Hru_segment(Ihru)) - poss) ! give it equally to each HRU surrounding it - inflow = inflow/REAL(Seg_hru_num(Hru_segment(Ihru))) + inflow_in = inflow_in/REAL(Seg_hru_num(Hru_segment(Ihru))) !negative flow is out of stream into riparian - Seg_ripflow(Hru_segment(Ihru)) = Seg_ripflow(Hru_segment(Ihru)) - inflow - inflow_in = SNGL(inflow*Timestep_seconds/(FT2_PER_ACRE*12.0)) !inch acre - Ripst_vol(Ihru) = Ripst_vol(Ihru) + inflow_in - Ripst_frac(Ihru)= SNGL(Ripst_vol(Ihru)/(Ripst_vol_max(Ihru)*thaw_frac)) - IF (Ripst_frac(Ihru)>1.0) Ripst_frac(Ihru) = 1.0 +! add this in and add Hortonian and Dunnian flow +! + IF (Ripst_frac(Ihru)<1.0) THEN + Seg_ripflow(Hru_segment(Ihru)) = Seg_ripflow(Hru_segment(Ihru)) - inflow_in + inflow_in = SNGL(inflow_in*Timestep_seconds/(FT2_PER_ACRE*12.0)) !inch acre + Ripst_vol(Ihru) = Ripst_vol(Ihru) + inflow*Ripst_areafr_max(Ihru) + inflow_in !inch acre + IF ( Ripst_vol(Ihru) > (Ripst_vol_max(Ihru)*thaw_frac) ) THEN + Ripst_vol(Ihru) = Ripst_vol_max(Ihru)*thaw_frac + ELSE + Ripst_vol(Ihru) = Ripst_vol(Ihru) + ripst_srp + ripst_sri + IF ( Ripst_vol(Ihru) > (Ripst_vol_max(Ihru)*thaw_frac) ) THEN + ripst_srp = SNGL((Ripst_vol(Ihru) - (Ripst_vol_max(Ihru)*thaw_frac))*ripst_srp/(ripst_srp + ripst_sri)) + ripst_sri = SNGL((Ripst_vol(Ihru) - (Ripst_vol_max(Ihru)*thaw_frac))*ripst_sri/(ripst_srp + ripst_sri)) + Ripst_vol(Ihru) = Ripst_vol_max(Ihru)*thaw_frac + ENDIF + IF ( Srp>0.0 ) THEN + Srp = Srp - ripst_srp/Perv_frac/Hru_area(Ihru) + IF ( Srp<0.0 ) THEN + IF ( Srp<-NEARZERO ) PRINT *, 'ripst srp<0.0', Srp, ripst_srp + ! may need to adjust ripst_srp and volumes + Srp = 0.0 + ENDIF + ENDIF + IF ( Sri>0.0 ) THEN + Sri = Sri - ripst_sri/Imperv_frac/Hru_area(Ihru) + IF ( Sri<0.0 ) THEN + IF ( Srp<-NEARZERO ) PRINT *, 'ripst sri<0.0', Sri, ripst_sri + ! may need to adjust ripst_sri and volumes + Sri = 0.0 + ENDIF + ENDIF + ENDIF + Ripst_frac(Ihru)= SNGL(Ripst_vol(Ihru)/(Ripst_vol_max(Ihru)*thaw_frac)) + ENDIF + ! Filled riparian storage surface area for each HRU: ! Fills outward from the river with one edge on river and with same depth and same side shape ! this works out to keeping fraction same for area and volume filled - Ripst_area(Ihru) = Ripst_area_max(Ihru)*Ripst_frac(Ihru) + Ripst_area(Ihru) = Ripst_area_max(Ihru)*Ripst_frac(Ihru) !acres ! evaporate water from riparian area based on snowcov_area - ! ripst_evap_open & ripst_evap_clos = inches-acres on the HRU + ! ripst_evap = inches-acres on the HRU unsatisfied_et = Potet(Ihru) - Snow_evap(Ihru) - Hru_intcpevap(Ihru) & & - Hru_impervevap(Ihru) - Dprst_evap_hru(Ihru) ripst_avail_et = 0.0 - ripst_avail_et = Potet(Ihru)*(1.0-Snowcov_area(Ihru))*Ripst_et_coef(Ihru) + ripst_avail_et = Potet(Ihru)*(1.0-Snowcov_area(Ihru))*Dprst_et_coef(Ihru) Ripst_evap_hru(Ihru) = 0.0 IF ( ripst_avail_et>0.0 ) THEN ripst_evap = 0.0 @@ -1224,29 +1311,33 @@ SUBROUTINE drain_the_swamp(Ihru) ! compute seepage seep = 0.0 seep_in = 0.0 + seep_surface = 0 !0 if just using same as depression storage IF ( Ripst_vol(Ihru)>NEARZERO) THEN - ripst_wid = SNGL(Ripst_area(Ihru)*FT2_PER_ACRE*(FEET2METERS**2.0)/Seg_length(Hru_segment(Ihru))) !meters -!assumed it was a one sided stream, here a headwater with both sides in one HRU - IF ( Seg_hru_num(Hru_segment(Ihru))==1 ) ripst_wid = ripst_wid/2.0 -! Stream ground area is stream side area (flat wall) and other side area (fraction of triangle (1) to rectangle (0)) - ripst_grnd = DBLE( Seg_length(Hru_segment(Ihru))*( ripst_wid*(1.0-Tr_ratio(Ihru)) + & !rectangle + IF (seep_surface==1) THEN !say the seep coefficient has something to do with surface area + ripst_wid = SNGL(Ripst_area(Ihru)*FT2_PER_ACRE*(FEET2METERS**2.0)/Seg_length(Hru_segment(Ihru))) !meters + !assumed it was a one sided stream, here a headwater with both sides in one HRU + IF ( Seg_hru_num(Hru_segment(Ihru))==1 ) ripst_wid = ripst_wid/2.0 + ! Stream ground area is stream side area (flat wall) and other side area (fraction of triangle (1) to rectangle (0)) + ripst_grnd = DBLE( Seg_length(Hru_segment(Ihru))*( ripst_wid*(1.0-Tr_ratio(Ihru)) + & !rectangle & (SQRT( ripst_wid**2.0 + (Ripst_depth(Ihru)*thaw_frac)**2.0 )- Ripst_depth(Ihru)*thaw_frac)*Tr_ratio(Ihru) + & !triangle - & 2.0*Ripst_depth(Ihru)*thaw_frac ) ) !stream and other side -!assumed it was a one sided stream, here a headwater with both sides in one HRU - IF ( Seg_hru_num(Hru_segment(Ihru))==1 ) ripst_grnd = ripst_grnd*2.0 -!seep in a day through ground surface area of riparian, m^3 into ft^3 to acre_in -!Transmissivity way too big - !seep = ripst_grnd*DBLE( Transmiss_seg(Hru_segment(Ihru)) )/(FEET2METERS**3.0) !acre_in -!ground area to total surface area is 5/6, then use depression seep coeff but reduce because surface area smaller - seep = ripst_grnd/(ripst_grnd+Ripst_area(Ihru)*FT2_PER_ACRE/ FEET2METERS**2.0 )/(5.0/6.0) & - & *Ripst_vol(Ihru)*Dprst_seep_rate_open(Ihru)/FT2_PER_ACRE/12.0 + & 2.0*Ripst_depth(Ihru)*thaw_frac ) ) !stream and other side m^2 + !assumed it was a one sided stream, here a headwater with both sides in one HRU + IF ( Seg_hru_num(Hru_segment(Ihru))==1 ) ripst_grnd = ripst_grnd*2.0 + !seep in a day through ground surface area of riparian, m^3 into ft^3 to acre_in + !Transmissivity way too big + !seep = ripst_grnd*DBLE( Transmiss_seg(Hru_segment(Ihru)) )/(FEET2METERS**3.0) !ft^3 + !ground area to total surface area is 5/6, then use depression seep coeff but reduce because surface area smaller + seep = ripst_grnd/(ripst_grnd+Ripst_area(Ihru)*FT2_PER_ACRE/ FEET2METERS**2.0 )/(5.0/6.0) & + & *Ripst_vol(Ihru)*Dprst_seep_rate_open(Ihru)*FT2_PER_ACRE/12.0 !ft^3 + ELSE !assume it is just a volume thing + seep = Ripst_vol(Ihru)*Dprst_seep_rate_open(Ihru)*FT2_PER_ACRE/12.0 !ft^3 + ENDIF !seep = 0.0 !if want to turn off seep - seep_in = seep*FT2_PER_ACRE*12.0 ! inch acres + seep_in = seep*12.0/FT2_PER_ACRE ! inch acres Ripst_vol(Ihru) = Ripst_vol(Ihru) - seep_in IF ( Ripst_vol(Ihru)<0.0D0 ) THEN !IF ( Ripst_vol(Ihru)<-DNEARZERO ) PRINT *, 'issue, ripst_vol<0.0', Ihru, Ripst_vol(Ihru) - seep_in = seep_in + Ripst_vol(Ihru) - seep = seep_in/FT2_PER_ACRE/12.0 !ft^3 + seep_in = seep_in + Ripst_vol(Ihru) !inch acre Ripst_vol(Ihru) = 0.0D0 ENDIF ENDIF @@ -1277,9 +1368,8 @@ SUBROUTINE init_bank_storage() USE PRMS_BASIN, ONLY: NEARZERO, Basin_area_inv, Hru_area_dble, Active_hrus, & & FT2_PER_ACRE, FEET2METERS, CFS2CMS_CONV USE PRMS_ROUTING, ONLY: Basin_bankst_head, Bankst_head_init, Basin_bankst_area, & - & Basin_bankst_vol, Bankst_head, Hru_segment, Seg_width, Seg_length, & - & Bankst_stor_hru, Bankst_head_pts, Ripst_areafr_max, Bankfinite_hru - USE PRMS_FLOWVARS, ONLY: Seg_outflow + & Basin_bankst_vol, Bankst_head, Stage_ante, Bankst_stor_hru, Bankst_head_pts, & + & Ripst_areafr_max, Bankfinite_hru, Hru_segment IMPLICIT NONE ! Functions INTRINSIC SNGL @@ -1289,8 +1379,9 @@ SUBROUTINE init_bank_storage() DO i = 1, Active_hrus IF ( Hru_segment(i)>0) THEN Bankst_head(i) = Bankst_head_init(i) - Bankst_head_pts(i) =SNGL(Seg_outflow(Hru_segment(i))*CFS2CMS_CONV)*60.*60.*24. & - & /Seg_width(Hru_segment(i))/Seg_length(Hru_segment(i)) + Bankst_head_pts(i,1) =( 8.0*Bankst_head(i)-SNGL(Stage_ante(Hru_segment(i))) )/4.0 + IF (Bankst_head_pts(i,1)<0.0) Bankst_head_pts(i,1) = 0.0 + Bankst_head_pts(i,1) = 0.0 IF (Bankfinite_hru(i)==1) THEN Bankst_stor_hru(i) = Ripst_areafr_max(i)*12.0*Bankst_head(i)/FEET2METERS !in inches Basin_bankst_head = Basin_bankst_head + Ripst_areafr_max(i)*Bankst_head(i)*Hru_area_dble(i) ! in meters @@ -1336,13 +1427,15 @@ SUBROUTINE comp_bank_storage(Ihru) INTEGER :: h, t0 INTEGER, PARAMETER :: nbankd = 2 REAL, PARAMETER :: PI = 3.14159 - REAL :: area, str_wid, tot_wid, bank_wid, trans, a, xd, t, td - REAL :: delt, delta_input(nbankd), delta_diff(nbankd), head(nbankd), seep(nbankd) + REAL :: area, str_wid, tot_wid, bank_wid, trans, a, xd, t, td, xd2 + REAL :: delt, delta_input(nbankd), delta_diff(nbankd), head(nbankd,2), seep(nbankd) REAL :: bank(nbankd), bankv(nbankd), ripfrac DOUBLE PRECISION :: input_net(nbankd), diff_net(nbankd), recharge(nbankd), stage(nbankd) DOUBLE PRECISION :: head_step, head_step_grad, seep_sum, head_sum + DOUBLE PRECISION :: head_step2, head_step_grad2, head_sum2 !*********************************************************************** area = Ripst_areafr_max(Ihru)*Hru_area(Ihru) !acres + IF ( Seg_hru_num(Hru_segment(Ihru))==1 ) area = area/2.0 !only take half of area if hru contains all of stream not just one side trans = Transmiss_seg(Hru_segment(Ihru)) !aquifer diffusivity, ratio of the transmissivity/storativity of the aquifer a = trans/Specyield_seg(Hru_segment(Ihru)) @@ -1371,29 +1464,38 @@ SUBROUTINE comp_bank_storage(Ihru) delta_diff(h-1) = SNGL((diff_net(h)-diff_net(h-1))/delt ) ENDDO Bankst_seep_hru(Ihru) = 0.0 - xd = 1.0+ bank_wid/2.0 ! at x = 1.0 is stage which already know, calc at middle of bank storage area - head=Bankst_head_pts(Ihru) !set at last height for initial + xd = 1.0+ bank_wid/4.0 ! at x = 1.0 is stage which already know, calc at 1/4 of bank storage area + DO h = 1, nbankd + head(h,1)=Bankst_head_pts(Ihru,1) !set at last height for initial + head(h,2)=Bankst_head_pts(Ihru,2) !set at last height for initial + ENDDO ! Calculate heads, seepage, and bank storage using convolution ripfrac = Ripst_areafr_max(Ihru) IF (Bankfinite_hru(Ihru)==0) ripfrac = 1.0 DO h = 1, (nbankd-1) head_sum = 0.0 + head_sum2 = 0.0 seep_sum = 0.0 DO t0 = 1,h t = t0*delt td = t*a/(str_wid**2.0) !dimensionless - IF (Bankfinite_hru(Ihru)==1) then !finite solution if transmissivity high, COMPUTATIONALLY EXPENSIVE, might eliminate + IF (Bankfinite_hru(Ihru)==1) then !finite solution if transmissivity high, COMPUTATIONALLY EXPENSIVE + xd2 = bank_wid+1.0 CALL LTST1(td, xd, tot_wid, bank_wid, head_step, head_step_grad) + CALL LTST1(td, xd2, tot_wid, bank_wid, head_step2, head_step_grad2) ELSE IF (Bankfinite_hru(Ihru)==0) then !semi-infinite solution head_step = ERFC( (xd - 1.0)/SQRT((4.0*td)) ) head_step_grad = -( 1.0/SQRT((PI*td)) ) ENDIF !head is a function of xd head_sum = delta_input(h-t0+1)*head_step + head_sum + IF (Bankfinite_hru(Ihru)==1) head_sum2 = delta_input(h-t0+1)*head_step2 + head_sum2 !seep is per unit segment length rate goes out, not a function of xd seep_sum = delta_diff(h-t0+1)*head_step_grad + seep_sum ENDDO - head(h+1)=head(h+1) + SNGL(head_sum*delt) + head(h+1,1)=head(h+1,1) + SNGL(head_sum*delt) + head(h+1,2)=0.0 ! Bankst_head_pts at infinite edge of bank storage area is 0 (xd = 1, so head_step = 0) + IF (Bankfinite_hru(Ihru)==1) head(h+1,2)=head(h+1,2) + SNGL(head_sum2*delt) seep(h+1)=SNGL((trans/str_wid)*seep_sum*delt) bank(h+1)=bank(h) - seep(h+1)*delt bankv(h+1)=bank(h+1)*Seg_length(Hru_segment(Ihru)) @@ -1405,12 +1507,11 @@ SUBROUTINE comp_bank_storage(Ihru) bank = bank*2.0 bankv = bankv*2.0 ENDIF - Bankst_head_pts(Ihru) = head(nbankd) ! meters + Bankst_head_pts(Ihru,1) = head(nbankd,1) ! meters + Bankst_head_pts(Ihru,2) = head(nbankd,2) ! meters !linear interpolation for total average head over bank storage area, meters - Bankst_head(Ihru) = 0.5*(SNGL(stage(nbankd))+Bankst_head_pts(Ihru)) - ! Bankst_head_pts at finite edge of bank storage area is 0 (xd = 1, so head_step = 0) - ! is only saved at the end of the timestep - Bankst_head(Ihru) = Bankst_head(Ihru) + 0.5*Bankst_head_pts(Ihru) + Bankst_head(Ihru) = ( SNGL(stage(nbankd))+4.0*Bankst_head_pts(Ihru,1) & + & + 3.0*Bankst_head_pts(Ihru,2) )/8.0 ! m2 per 24 hr per stream segment for both sides of stream ! seep hru is inch over hru seeping out per day diff --git a/prmsRip/snowcompCfgim.f90 b/prmsRip/snowcompCfgim.f90 index 09a30494..b32519ff 100644 --- a/prmsRip/snowcompCfgim.f90 +++ b/prmsRip/snowcompCfgim.f90 @@ -174,107 +174,107 @@ INTEGER FUNCTION snodecl() ALLOCATE ( Glacr_freeh2o_capm(Nhru) ) IF ( declvar(MODNAME, 'glacr_freeh2o_capm', 'nhru', Nhru, 'real', & - & 'Free-water holding capacity of glacier ice, changes to 0 if active layer melts', & + & 'Free-water holding capacity of glacier or glacierette ice, changes to 0 if active layer melts', & & 'decimal fraction', Glacr_freeh2o_capm)/=0 ) CALL read_error(3, 'glacr_freeh2o_capm') ALLOCATE ( Glacrb_melt(Nhru) ) IF ( declvar(MODNAME, 'glacrb_melt', 'nhru', Nhru, 'real', & - 'Glacier basal melt, goes to soil', & + 'Glacier or glacierette basal melt, goes to soil', & 'inches/day', Glacrb_melt)/=0 ) CALL read_error(3, 'glacrb_melt') ALLOCATE ( Glacr_air_5avtemp(Nhru) ) IF ( declvar(MODNAME, 'glacr_air_5avtemp', 'nhru', Nhru, 'real', & - & 'Current 5-yr average summer (June July Aug) air temperature over glacier or glrette HRU', & + & 'Current 5-yr average summer (June July Aug) air temperature over glacier or glacierette HRU', & & 'degrees Celsius', Glacr_air_5avtemp)/=0 ) CALL read_error(3, 'glacr_air_5avtemp') ALLOCATE ( Glacr_air_5avtemp1(Nhru) ) IF ( declvar(MODNAME, 'glacr_air_5avtemp1', 'nhru', Nhru, 'real', & - & 'First 5-yr average summer temperature over glacier or glrette HRU', & + & 'First 5-yr average summer temperature over glacier or glacierette HRU', & & 'degrees Celsius', Glacr_air_5avtemp1)/=0 ) CALL read_error(3, 'glacr_air_5avtemp1') ALLOCATE ( Glacr_air_deltemp(Nhru) ) IF ( declvar(MODNAME, 'glacr_air_deltemp', 'nhru', Nhru, 'real', & - & 'Change in 5-yr average air temperature over glacier or glrette HRU from first', & + & 'Change in 5-yr average air temperature over glacier or glacierette HRU from first', & & 'degrees Celsius', Glacr_air_deltemp)/=0 ) CALL read_error(3, 'glacr_air_deltemp') ALLOCATE ( Glacr_5avsnow(Nhru) ) IF ( declvar(MODNAME, 'glacr_5avsnow', 'nhru', Nhru, 'real', & - & 'Current 5-yr average snow over glacier or glrette HRU', & + & 'Current 5-yr average snow over glacier or glacierette HRU', & & 'inches/yr', Glacr_5avsnow)/=0 ) CALL read_error(3, 'glacr_5avsnow') ALLOCATE ( Glacr_5avsnow1(Nhru) ) IF ( declvar(MODNAME, 'glacr_5avsnow1', 'nhru', Nhru, 'real', & - & 'First 5-yr average snow over glacier or glrette HRU', & + & 'First 5-yr average snow over glacier or glacierette HRU', & & 'inches/yr', Glacr_5avsnow1)/=0 ) CALL read_error(3, 'glacr_5avsnow1') ALLOCATE ( Glacr_delsnow(Nhru) ) IF ( declvar(MODNAME, 'glacr_delsnow', 'nhru', Nhru, 'real', & - & 'Change in 5-yr average snow over glacier or glrette HRU from first', & + & 'Change in 5-yr average snow over glacier or glacierette HRU from first', & & 'inches/yr', Glacr_delsnow)/=0 ) CALL read_error(3, 'glacr_delsnow') ALLOCATE ( Glacr_pk_temp(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_temp', 'nhru', Nhru, 'real', & - & 'Temperature of the glacier on each HRU', & + & 'Temperature of the glacier or glacierette on each HRU', & & 'degrees Celsius', Glacr_pk_temp)/=0 ) CALL read_error(3, 'glacr_pk_temp') ALLOCATE ( Glacr_pk_def(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_def', 'nhru', Nhru, 'real', & - & 'Heat deficit, amount of heat necessary to make the glacier snowpack isothermal at 0 degrees Celsius', & + & 'Heat deficit, amount of heat necessary to make the glacier or or glacierette snowpack isothermal at 0 degrees Celsius', & & 'Langleys', Glacr_pk_def)/=0 ) CALL read_error(3, 'glacr_pk_def') ALLOCATE ( Glacr_pk_den(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_den', 'nhru', Nhru, 'real', & - & 'Density of the icepack on each glacier HRU, hard-coded to equal 0.917', & + & 'Density of the icepack on each glacier or glacierette HRU, hard-coded to equal 0.917', & & 'gm/cm3', Glacr_pk_den)/=0 ) CALL read_error(3, 'glacr_pk_den') ALLOCATE ( Glacr_albedo(Nhru) ) IF ( declvar(MODNAME, 'glacr_albedo', 'nhru', Nhru, 'real', & - & 'Ice surface albedo or the fraction of radiation reflected from the icepack surface for each glacier HRU', & + & 'Ice surface albedo or the fraction of radiation reflected from the icepack surface for each glacier or glacierette HRU', & & 'decimal fraction', Glacr_albedo)/=0 ) CALL read_error(3, 'glacr_albedo') ALLOCATE ( Glacr_evap(Nhru) ) IF ( declvar(MODNAME, 'glacr_evap', 'nhru', Nhru, 'real', & - & 'Evaporation and sublimation from icepack on each glacier HRU', & + & 'Evaporation and sublimation from icepack on each glacier or glacierette HRU', & & 'inches', Glacr_evap)/=0 ) CALL read_error(3, 'glacr_evap') ALLOCATE ( Glacrmelt(Nhru) ) IF ( declvar(MODNAME, 'glacrmelt', 'nhru', Nhru, 'real', & - & 'Melt from icepack on each glacier HRU, includes rain water that does not absorb', & + & 'Melt from icepack on each glacier or glacierette HRU, includes rain water that does not absorb', & & 'inches', Glacrmelt)/=0 ) CALL read_error(3, 'glacrmelt') ALLOCATE ( Glacr_pkwater_equiv(Nhru) ) IF ( declvar(MODNAME, 'glacr_pkwater_equiv', 'nhru', Nhru, 'double', & - & 'Icepack water equivalent on each glacier HRU', & + & 'Icepack water equivalent on each glacier or glacierette HRU', & & 'inches', Glacr_pkwater_equiv)/=0 ) CALL read_error(3, 'glacr_pkwater_equiv') ALLOCATE ( Glacr_pkwater_ante(Nhru) ) IF ( declvar(MODNAME, 'glacr_pkwater_ante', 'nhru', Nhru, 'double', & - & 'Antecedent icepack water equivalent on each glacier HRU', & + & 'Antecedent icepack water equivalent on each glacier or glacierette HRU', & & 'inches', Glacr_pkwater_ante)/=0 ) CALL read_error(3, 'glacr_pkwater_ante') ALLOCATE ( Glacrcov_area(Nhru) ) IF ( declvar(MODNAME, 'glacrcov_area', 'nhru', Nhru, 'real', & - & 'Ice-covered area on each glacier HRU or HRU with glacierette at start of step', & + & 'Ice-covered area (no snowpack) on each glacier HRU or HRU with glacierette at start of step', & & 'decimal fraction', Glacrcov_area)/=0 ) CALL read_error(3, 'glacrcov_area') ALLOCATE ( Glacr_pk_ice(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_ice', 'nhru', Nhru, 'real', & - & 'Storage of frozen water in the icepack on each glacier HRU', & + & 'Storage of frozen water in the icepack on each glacier or glacierette HRU', & & 'inches', Glacr_pk_ice)/=0 ) CALL read_error(3, 'glacr_pk_ice') ALLOCATE ( Glacr_freeh2o(Nhru) ) IF ( declvar(MODNAME, 'glacr_freeh2o', 'nhru', Nhru, 'real', & - & 'Storage of free liquid water in the icepack on each glacier HRU', & + & 'Storage of free liquid water in the icepack on each glacier or glacierette HRU', & & 'inches', Glacr_freeh2o)/=0 ) CALL read_error(3, 'glacr_freeh2o') ALLOCATE ( Glacr_pk_depth(Nhru) ) IF ( declvar(MODNAME, 'glacr_pk_depth', 'nhru', Nhru, 'double', & - & 'Depth of icepack on each glacier HRU, make essentially infinite', & + & 'Depth of icepack on each glacier or glacierette HRU, make essentially infinite', & & 'inches', Glacr_pk_depth)/=0 ) CALL read_error(3, 'glacr_pk_depth') ALLOCATE ( Glacr_pss(Nhru) ) IF ( declvar(MODNAME, 'glacr_pss', 'nhru', Nhru, 'double', & - & 'Previous glacier pack water equivalent plus new ice', & + & 'Previous glacier or glacierette pack water equivalent plus new ice', & & 'inches', Glacr_pss)/=0 ) CALL read_error(3, 'glacr_pss') ALLOCATE ( Glacr_pst(Nhru) ) @@ -283,41 +283,9 @@ INTEGER FUNCTION snodecl() & 'inches', Glacr_pst)/=0 ) CALL read_error(3, 'glacr_pst') IF ( declvar(MODNAME, 'basin_snowicecov', 'one', 1, 'double', & - & 'Basin area-weighted average snow and glacier and glrette covered area', & + & 'Basin area-weighted average snow and glacier and glacierette covered area for calibration to satellites', & & 'decimal fraction', Basin_snowicecov)/=0 ) CALL read_error(3, 'basin_snowicecov') - ALLOCATE ( Glacr_freeh2o_cap(Nhru) ) - IF ( declparam(MODNAME, 'glacr_freeh2o_cap', 'nhru', 'real', & - & '0.002', '0.0', '0.01', & - & 'Free-water holding capacity of glacier ice', & - & 'Free-water holding capacity of glacier ice expressed as a' // & - & ' decimal fraction of the frozen water content of the glacier ice (glacr_pk_ice)', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'glacr_freeh2o_cap') - - ALLOCATE ( Glacr_layer(Nhru) ) - IF ( declparam(MODNAME, 'glacr_layer', 'nhru', 'real', & - & '3.94', '0.0', '590.6', & - & 'Active layer on glacier', & - & 'Active layer is 0 to 15 m (590.6 inches) thick at start of year, when' // & - & ' melts will set daily glacr_pk_temp to 0', & - & 'inches')/=0 ) CALL read_error(1, 'glacr_layer') - - IF ( Init_vars_from_file==0 ) THEN - ALLOCATE ( Glacier_frac_init(Nhru) ) - IF ( declparam(MODNAME, 'glacier_frac_init', 'nhru', 'real', & - & '0.0', '0.0', '1.0', & - & 'Inital fraction of glaciation (0=none; 1=100%)', & - & 'Inital fraction of glaciation (0=none; 1=100%)', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'glacier_frac_init') - - ALLOCATE ( Glrette_frac_init(Nhru) ) - IF ( declparam(MODNAME, 'glrette_frac_init', 'nhru', 'real', & - & '0.0', '0.0', '1.0', & - & 'Initial fraction of glacierette (too small for glacier dynamics)', & - & 'Initial fraction of glacierette (too small for glacier dynamics)', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'glrette_frac_init') - - ENDIF ENDIF IF ( Frozen_flag==1 .OR. Model==99 ) THEN @@ -526,8 +494,8 @@ INTEGER FUNCTION snodecl() ALLOCATE ( Albedo_coef(Nhru) ) IF ( declparam(MODNAME, 'albedo_coef', 'nhru', 'real', & & '0.137', '0.1', '0.3', & - & 'Coefficient in calculation of ice albedo', & - & 'Coefficient in calculation of ice albedo', & + & 'Coefficient in calculation of ice albedo for glaciers', & + & 'Coefficient in calculation of ice albedo for glaciers', & & 'none')/=0 ) CALL read_error(1, 'albedo_coef') ALLOCATE ( Albedo_ice(Nhru) ) @@ -536,6 +504,39 @@ INTEGER FUNCTION snodecl() & 'Ice albedo 300 meters below ELA', & & 'Ice albedo 300 meters below ELA', & & 'decimal fraction')/=0 ) CALL read_error(1, 'albedo_ice') + + ALLOCATE ( Glacr_freeh2o_cap(Nhru) ) + IF ( declparam(MODNAME, 'glacr_freeh2o_cap', 'nhru', 'real', & + & '0.002', '0.0', '0.01', & + & 'Free-water holding capacity of glacier ice', & + & 'Free-water holding capacity of glacier ice expressed as a' // & + & ' decimal fraction of the frozen water content of the glacr_pk_ice', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'glacr_freeh2o_cap') + + ALLOCATE ( Glacr_layer(Nhru) ) + IF ( declparam(MODNAME, 'glacr_layer', 'nhru', 'real', & + & '3.94', '0.0', '590.6', & + & 'Active layer on glacier', & + & 'Active layer is 0 to 15 m (590.6 inches) thick at start of year, when' // & + & ' melts will set daily glacr_pk_temp to 0', & + & 'inches')/=0 ) CALL read_error(1, 'glacr_layer') + + IF ( Init_vars_from_file==0 ) THEN + ALLOCATE ( Glacier_frac_init(Nhru) ) + IF ( declparam(MODNAME, 'glacier_frac_init', 'nhru', 'real', & + & '0.0', '0.0', '1.0', & + & 'Inital fraction of glaciation (0=none; 1=100%) in glacier-capable HRU', & + & 'Inital fraction of glaciation (0=none; 1=100%) in glacier-capable HRU', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'glacier_frac_init') + + ALLOCATE ( Glrette_frac_init(Nhru) ) + IF ( declparam(MODNAME, 'glrette_frac_init', 'nhru', 'real', & + & '0.0', '0.0', '1.0', & + & 'Initial fraction of glacierette (too small for glacier dynamics)', & + & 'Initial fraction of glacierette (too small for glacier dynamics)', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'glrette_frac_init') + + ENDIF ENDIF IF ( declparam(MODNAME, 'den_init', 'one', 'real', & diff --git a/prmsRip/srunoffCfgim.f90 b/prmsRip/srunoffCfgim.f90 index 24815e85..6a1e2c67 100644 --- a/prmsRip/srunoffCfgim.f90 +++ b/prmsRip/srunoffCfgim.f90 @@ -98,7 +98,7 @@ INTEGER FUNCTION srunoffdecl() USE PRMS_SRUNOFF USE PRMS_MODULE, ONLY: Model, Dprst_flag, Nhru, Nsegment, Print_debug, & & Cascade_flag, Sroff_flag, Nlake, Init_vars_from_file, Call_cascade, PRMS4_flag, & - & Frozen_flag + & Frozen_flag, Ripst_flag IMPLICIT NONE ! Functions INTEGER, EXTERNAL :: declvar, declparam @@ -302,8 +302,8 @@ INTEGER FUNCTION srunoffdecl() IF ( Frozen_flag==1 .OR. Model==99 ) THEN ALLOCATE ( Frozen(Nhru) ) IF ( declvar(MODNAME, 'frozen', 'nhru', Nhru, 'integer', & - & 'Flag for frozen ground (0=no; 1=soil at surface; 2=soil below surf; 3=below soil)', & - & 'dimensionless', Frozen)/=0 ) CALL read_error(3, 'frozen') + & 'Marker for frozen ground layer top (0=none; 1=at surface; 2=in soil below surface; 3=below soil)', & + & 'none', Frozen)/=0 ) CALL read_error(3, 'frozen') ALLOCATE ( Cfgi(Nhru) ) IF ( declvar(MODNAME, 'cfgi', 'nhru', Nhru, 'real', & @@ -422,14 +422,6 @@ INTEGER FUNCTION srunoffdecl() & 'Coefficient in linear flow routing equation for open surface depressions for each HRU', & & 'fraction/day')/=0 ) CALL read_error(1, 'dprst_flow_coef') - ALLOCATE ( Dprst_seep_rate_open(Nhru) ) - IF ( declparam(MODNAME, 'dprst_seep_rate_open', 'nhru', 'real', & - & '0.02', '0.0', '0.2', & - & 'Coefficient used in linear seepage flow equation for open surface depressions', & - & 'Coefficient used in linear seepage flow equation for'// & - & ' open surface depressions for each HRU', & - & 'fraction/day')/=0 ) CALL read_error(1, 'dprst_seep_rate_open') - ALLOCATE ( Dprst_seep_rate_clos(Nhru) ) IF ( declparam(MODNAME, 'dprst_seep_rate_clos', 'nhru', 'real', & & '0.02', '0.0', '0.2', & @@ -447,41 +439,6 @@ INTEGER FUNCTION srunoffdecl() & ' maximum open storage capacity spills as surface runoff', & & 'decimal fraction')/=0 ) CALL read_error(1, 'op_flow_thres') - ALLOCATE ( Sro_to_dprst_perv(Nhru) ) - IF ( PRMS4_flag==1 ) THEN - IF ( declparam(MODNAME, 'sro_to_dprst', 'nhru', 'real', & - & '0.2', '0.0', '1.0', & - & 'Fraction of pervious surface runoff that flows into surface-depression storage', & - & 'Fraction of pervious surface runoff that'// & - & ' flows into surface-depression storage; the remainder'// & - & ' flows to a stream network for each HRU', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'sro_to_dprst') - ELSE - IF ( declparam(MODNAME, 'sro_to_dprst_perv', 'nhru', 'real', & - & '0.2', '0.0', '1.0', & - & 'Fraction of pervious surface runoff that flows into surface-depression storage', & - & 'Fraction of pervious surface runoff that'// & - & ' flows into surface-depression storage; the remainder'// & - & ' flows to a stream network for each HRU', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'sro_to_dprst_perv') - ENDIF - - ALLOCATE ( Sro_to_dprst_imperv(Nhru) ) - IF ( declparam(MODNAME, 'sro_to_dprst_imperv', 'nhru', 'real', & - & '0.2', '0.0', '1.0', & - & 'Fraction of impervious surface runoff that flows into surface-depression storage', & - & 'Fraction of impervious surface runoff that'// & - & ' flows into surface-depression storage; the remainder'// & - & ' flows to a stream network for each HRU', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'sro_to_dprst_imperv') - - ALLOCATE ( Dprst_et_coef(Nhru) ) - IF ( declparam(MODNAME, 'dprst_et_coef', 'nhru', 'real', & - & '1.0', '0.5', '1.5', & - & 'Fraction of unsatisfied potential evapotranspiration to apply to surface-depression storage', & - & 'Fraction of unsatisfied potential evapotranspiration to apply to surface-depression storage', & - & 'decimal fraction')/=0 ) CALL read_error(1, 'dprst_et_coef') - IF ( Init_vars_from_file==0 .OR. Init_vars_from_file==2 .OR. Init_vars_from_file==7 ) THEN ALLOCATE ( Dprst_frac_init(Nhru) ) IF ( declparam(MODNAME, 'dprst_frac_init', 'nhru', 'real', & @@ -515,6 +472,51 @@ INTEGER FUNCTION srunoffdecl() & 'none')/=0 ) CALL read_error(1, 'va_clos_exp') ENDIF + IF ( Dprst_flag==1 .OR. Model==99 .OR. Ripst_flag==1) THEN + ALLOCATE ( Sro_to_dprst_perv(Nhru) ) + IF ( PRMS4_flag==1 ) THEN + IF ( declparam(MODNAME, 'sro_to_dprst', 'nhru', 'real', & + & '0.2', '0.0', '1.0', & + & 'Fraction of pervious surface runoff that flows into surface-depression or riparian storage', & + & 'Fraction of pervious surface runoff that'// & + & ' flows into surface-depression or riparian storage; the remainder'// & + & ' flows to a stream network for each HRU', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'sro_to_dprst') + ELSE + IF ( declparam(MODNAME, 'sro_to_dprst_perv', 'nhru', 'real', & + & '0.2', '0.0', '1.0', & + & 'Fraction of pervious surface runoff that flows into surface-depression or riparian storage', & + & 'Fraction of pervious surface runoff that'// & + & ' flows into surface-depression or riparian storage; the remainder'// & + & ' flows to a stream network for each HRU', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'sro_to_dprst_perv') + ENDIF + + ALLOCATE ( Sro_to_dprst_imperv(Nhru) ) + IF ( declparam(MODNAME, 'sro_to_dprst_imperv', 'nhru', 'real', & + & '0.2', '0.0', '1.0', & + & 'Fraction of impervious surface runoff that flows into surface-depression or riparian storage', & + & 'Fraction of impervious surface runoff that'// & + & ' flows into surface-depression or riparian storage; the remainder'// & + & ' flows to a stream network for each HRU', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'sro_to_dprst_imperv') + + ALLOCATE ( Dprst_seep_rate_open(Nhru) ) + IF ( declparam(MODNAME, 'dprst_seep_rate_open', 'nhru', 'real', & + & '0.02', '0.0', '0.2', & + & 'Coefficient used in linear seepage flow equation for open surface depressions or riparian storage', & + & 'Coefficient used in linear seepage flow equation for'// & + & ' open surface depressions or riparian storage for each HRU', & + & 'fraction/day')/=0 ) CALL read_error(1, 'dprst_seep_rate_open') + + ALLOCATE ( Dprst_et_coef(Nhru) ) + IF ( declparam(MODNAME, 'dprst_et_coef', 'nhru', 'real', & + & '1.0', '0.5', '1.5', & + & 'Fraction of unsatisfied potential evapotranspiration to apply to surface-depression or riparian storage', & + & 'Fraction of unsatisfied potential evapotranspiration to apply to surface-depression or riparian storage', & + & 'decimal fraction')/=0 ) CALL read_error(1, 'dprst_et_coef') + ENDIF + IF ( Print_debug==1 ) THEN ALLOCATE ( Imperv_stor_ante(Nhru) ) IF ( Dprst_flag==1 ) ALLOCATE ( Dprst_stor_ante(Nhru) ) @@ -685,7 +687,7 @@ INTEGER FUNCTION srunoffrun() REAL :: srunoff, avail_et, hperv, sra, availh2o DOUBLE PRECISION :: hru_sroff_down, runoff, apply_sroff, cfgi_sroff REAL :: cfgi_k, depth_cm, nosnow_area, depthg_cm, trad, emiss, emisl ! frozen ground - REAL :: cfgi_kg, soil_cond, latent_soil, nice, ice_cond, beta, thaw_frac ! frozen ground + REAL :: cfgi_kg, soil_cond, latent_soil, nwat, ice_cond, beta, thaw_frac ! frozen ground REAL :: water_cond, sat_cond, mean_cond, lambda, omega, l5, l6, l8 ! frozen ground REAL :: volumetric_soil, thermal_ratio_alp, fusion_param_mu, frz_height ! frozen ground REAL :: glcrmltb, temp, temp2 ! glaciers @@ -814,30 +816,36 @@ INTEGER FUNCTION srunoffrun() depthg_cm = 0.0 !Cov_type =0 bare soil (rock, may be mostly impervious already) IF (Cov_type(i)==1) depthg_cm = 4.0 !grasses (boreal grass, tundra) IF (Cov_type(i)==2) depthg_cm = 3.0 !shrub (tundra) - IF (Cov_type(i)>=3) depthg_cm = 6.0 !trees + IF (Cov_type(i)==3) depthg_cm = 6.0 !trees IF (Cov_type(i)==4) depthg_cm = 2.0 !coniferous ! Continuous frozen ground index Cfgi(i) = Cfgi_decay*Cfgi_prev(i) - trad*( 2.71828**(-0.4*(cfgi_k*depth_cm+cfgi_kg*depthg_cm)) ) IF ( active_glacier==1 ) THEN Cfgi(i) = 0.0 !if glacier over, want ground completely unfrozen, or below threshold, infiltration - IF ( Glacier_frac(i)<1.0 ) Cfgi(i) = Cfgi_thrshld ! glacier with some open fraction + IF ( Glacier_frac(i)<1.0 ) Cfgi(i) = Cfgi_thrshld ! glacier with some open fraction is frozen tongue ENDIF IF ( Cfgi(i)<0.0 ) Cfgi(i) = 0.0 ! If above the threshold to be frozen IF ( Cfgi(i)>=Cfgi_thrshld ) THEN + thaw_frac = 1.0 !use previous frozen state + IF ( Frozen(i)==1 ) THEN + thaw_frac = 0.0 + ELSEIF ( Frozen(i)==2) THEN + thaw_frac = Thaw_depth(i)/Soil_depth(i) + ENDIF ! Use modified Berggren formula to get a depth of frozen ! soil water content % of dry weight is water vol*density / (soil vol*density) omega = Soil_water(i) / (Soil_depth(i)*Soil_den(i)) IF ( omega>1.0 ) omega = 1.0 IF ( omega<0.1 ) omega = 0.1 ! volumetric heat of fusion of the soil - volumetric_soil = Soil_den(i)*(4.187*0.17 + 0.75*omega)*1.e6 ! J/m^3/K, specific heat of rock, water, ice =0.17, 1, 0.5 *4.187 J/g/K , density in g/cm3 + volumetric_soil = Soil_den(i)*(0.71179 + 4.186*(0.5*thaw_frac + 0.5)*omega)*1.e6 ! J/m^3/K, specific heat of rock, water, ice =4.186*0.17,4.186, 0.5 *4.187 J/g/K , density in g/cm3 ! latent heat of fusion of the soil latent_soil = 334.0*Soil_den(i)*omega*1.e6 ! J/m^3, latent heat of fusion of water = 334 J/g , density in g/cm3 thermal_ratio_alp = (Prev_ann_tempc(i) - Freezepoint)/(Cfgi(i) - Cfgi_thrshld) !degree K/ index Ti/Ts IF ( thermal_ratio_alp<0.0 ) thermal_ratio_alp = 0.0 - fusion_param_mu =(Cfgi(i) - Cfgi_thrshld)*volumetric_soil/latent_soil !index/degree K St12 + fusion_param_mu =(Cfgi(i) - Cfgi_thrshld)*volumetric_soil/latent_soil !index/degree K ! lambda corrects the Stefan formula for the effects of volumetric heat which it neglected beta = 1.0 !ranges between 0.95 and 1.3 depending on soil type and soil moisture lambda = 1.0 !Graph in Aldrich 1956, says in Alaska this is usually 1 but if less northern, can be as low as 0.3 @@ -849,13 +857,13 @@ INTEGER FUNCTION srunoffrun() IF ( Cfgi(i)==Cfgi_prev(i) ) lambda = l5 ! dry soil thermal conductivity - soil_cond = ( 486.0*Soil_den(i) + 233.0 )/( 2.7 - 0.947*Soil_den(i) ) !equation Johansen 1975,J/m/hr/K - !from last time step frozen depth - nice = Porosity_hru(i)* Frz_depth(i)/Soil_depth(i) + soil_cond = ( 486.0*Soil_den(i) + 233.0 )/( 2.7 - 0.947*Soil_den(i) ) !Via Fox 1992, equation Johansen 1975,J/m/hr/K, 1000 kg/m3 = 1 g/cm3 + !from last time step thaw_frac + nwat = Porosity_hru(i)* thaw_frac ! soil saturated conductivity is geometric mean of the conductivities of the materials within the soil profile ice_cond = (-0.0176*Tavgc(i) + 2.0526)*3600.0 !Bonales 2017, J/s/m/K to hr water_cond = 0.5918 *3600.0 !J/s/m/K to hr - sat_cond =( soil_cond**(1.0-Porosity_hru(i)) )*( ice_cond**(nice) )*( water_cond**(Porosity_hru(i)-nice) ) + sat_cond =( soil_cond**(1.0-Porosity_hru(i)) )*( water_cond**(nwat) )*( ice_cond**(Porosity_hru(i)-nwat) ) ! mean thermal conductivity of the frozen and unfrozen soil equation of dry and saturated conductivities mean_cond = (sat_cond - soil_cond)*omega + soil_cond !J/m/hr/K ! this is height of frozen soil. Freezes downward from surface. When thaw, also thaw downward from surface so will be thawed area above here @@ -874,7 +882,8 @@ INTEGER FUNCTION srunoffrun() frzen = 2 !soil frozen below top thaw_frac = Thaw_depth(i)/Soil_depth(i) ELSEIF ( Frz_depth(i)>=Soil_depth(i) ) THEN ! Thaw_depth(i)>=Soil_depth(i)) - frzen = 3 !soil not frozen but below is, thaw_frac = 1.0 + frzen = 3 !soil not frozen but below is + thaw_frac = 1.0 ENDIF ENDIF ENDIF @@ -1484,12 +1493,13 @@ SUBROUTINE dprst_comp(Dprst_vol_clos, Dprst_area_clos_max, Dprst_area_clos, & ENDIF Dprst_in = 0.0D0 + IF ( Dprst_area_open_max>0.0 ) THEN - Dprst_in = DBLE( inflow*Dprst_area_open_max*Thaw_frac ) ! inch-acres + Dprst_in = DBLE( inflow*Dprst_area_open_max ) ! inch-acres Dprst_vol_open = Dprst_vol_open + Dprst_in ENDIF IF ( Dprst_area_clos_max>0.0 ) THEN - tmp1 = DBLE( inflow*Dprst_area_clos_max*Thaw_frac ) ! inch-acres + tmp1 = DBLE( inflow*Dprst_area_clos_max ) ! inch-acres Dprst_vol_clos = Dprst_vol_clos + tmp1 Dprst_in = Dprst_in + tmp1 ENDIF @@ -1543,6 +1553,7 @@ SUBROUTINE dprst_comp(Dprst_vol_clos, Dprst_area_clos_max, Dprst_area_clos, & ! Open depression surface area for each HRU: Dprst_area_open = 0.0 IF ( Dprst_vol_open>0.0D0 ) THEN +! Thaw_frac reduces the volume the new water can add to, so the new water will spill open_vol_r = SNGL( Dprst_vol_open/(Dprst_vol_open_max*Thaw_frac) ) IF ( open_vol_rDprst_area_open_max*Thaw_frac ) Dprst_area_open = Dprst_area_open_max*Thaw_frac + Dprst_area_open = Dprst_area_open_max*frac_op_ar + IF ( Dprst_area_open>Dprst_area_open_max ) Dprst_area_open = Dprst_area_open_max ! IF ( Dprst_area_openDprst_area_clos_max*Thaw_frac ) Dprst_area_clos = Dprst_area_clos_max*Thaw_frac + Dprst_area_clos = Dprst_area_clos_max*frac_cl_ar + IF ( Dprst_area_clos>Dprst_area_clos_max) Dprst_area_clos = Dprst_area_clos_max ! IF ( Dprst_area_clos