diff --git a/.gitignore b/.gitignore index a93e28036..dff5455ba 100755 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,5 @@ Makefile-* *.backup # site directory from mkdocs site/ +# vscode setting files +.vscode \ No newline at end of file diff --git a/build/Makefile b/build/Makefile index 207e1e9bd..26ab91f95 100644 --- a/build/Makefile +++ b/build/Makefile @@ -216,6 +216,7 @@ SUMMA_SOLVER= \ opSplittin.f90 \ coupled_em.f90 \ run_oneHRU.f90 \ + HDS.f90 \ run_oneGRU.f90 SOLVER = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_SOLVER)) diff --git a/build/source/driver/.gitignore b/build/source/driver/.gitignore deleted file mode 100644 index 4a8d9bbfd..000000000 --- a/build/source/driver/.gitignore +++ /dev/null @@ -1 +0,0 @@ -summaversion.inc diff --git a/build/source/driver/summa_modelRun.f90 b/build/source/driver/summa_modelRun.f90 index ea9eb8df6..84e8a230d 100755 --- a/build/source/driver/summa_modelRun.f90 +++ b/build/source/driver/summa_modelRun.f90 @@ -250,6 +250,7 @@ subroutine summa_runPhysics(modelTimeStep, summa1_struc, err, message) typeStruct%gru(iGRU), & ! intent(in): local classification of soil veg etc. for each HRU idStruct%gru(iGRU), & ! intent(in): local classification of soil veg etc. for each HRU attrStruct%gru(iGRU), & ! intent(in): local attributes for each HRU + bparStruct%gru(iGRU), & ! intent(in): basin-average parameters for HDS ! data structures (input-output) mparStruct%gru(iGRU), & ! intent(inout): local model parameters indxStruct%gru(iGRU), & ! intent(inout): model indices diff --git a/build/source/driver/summa_restart.f90 b/build/source/driver/summa_restart.f90 index 1bef46d47..29049521d 100755 --- a/build/source/driver/summa_restart.f90 +++ b/build/source/driver/summa_restart.f90 @@ -30,6 +30,7 @@ module summa_restart USE var_lookup,only:iLookDIAG ! look-up values for local column model diagnostic variables USE var_lookup,only:iLookFLUX ! look-up values for local column model fluxes USE var_lookup,only:iLookBVAR ! look-up values for basin-average model variables +USE var_lookup,only:iLookBPAR ! look-up values for basin-average model parameters used by HDS USE var_lookup,only:iLookDECISIONS ! look-up values for model decisions ! safety: set private unless specified otherwise @@ -67,7 +68,9 @@ subroutine summa_readRestart(summa1_struc, err, message) ! model decisions USE mDecisions_module,only:& ! look-up values for the choice of method for the spatial representation of groundwater localColumn, & ! separate groundwater representation in each local soil column - singleBasin ! single groundwater store over the entire basin + singleBasin, & ! single groundwater store over the entire basin + HDSmodel ! Hysteretic Depressional Storage model implementation for prairie potholes + USE HDS,only:init_summa_HDS ! initialize HDS variables at the GRU level ! --------------------------------------------------------------------------------------- ! * variables ! --------------------------------------------------------------------------------------- @@ -227,6 +230,12 @@ subroutine summa_readRestart(summa1_struc, err, message) dt_init%gru(iGRU)%hru(iHRU) = progStruct%gru(iGRU)%hru(iHRU)%var(iLookPROG%dt_init)%dat(1) ! seconds end do + ! ***************************************************************************** + ! *** initialize HDS variables at the GRU level + ! ***************************************************************************** + if(model_decisions(iLookDECISIONS%prPotholes)%iDecision == HDSmodel)then + call init_summa_HDS(bvarStruct%gru(iGRU), bparStruct%gru(iGRU)) + endif end do ! end looping through GRUs ! ***************************************************************************** diff --git a/build/source/driver/summa_setup.f90 b/build/source/driver/summa_setup.f90 index e406dfd4b..7b9505c86 100755 --- a/build/source/driver/summa_setup.f90 +++ b/build/source/driver/summa_setup.f90 @@ -29,7 +29,7 @@ module summa_setup USE var_lookup,only:iLookATTR ! look-up values for local attributes USE var_lookup,only:iLookTYPE ! look-up values for classification of veg, soils etc. USE var_lookup,only:iLookPARAM ! look-up values for local column model parameters -USE var_lookup,only:iLookID ! look-up values for local column model parameters +USE var_lookup,only:iLookID ! look-up values for local column model parameters USE var_lookup,only:iLookBVAR ! look-up values for basin-average model variables USE var_lookup,only:iLookDECISIONS ! look-up values for model decisions USE globalData,only:urbanVegCategory ! vegetation category for urban areas diff --git a/build/source/dshare/get_ixname.f90 b/build/source/dshare/get_ixname.f90 index 8243d20a2..38610c4c2 100755 --- a/build/source/dshare/get_ixname.f90 +++ b/build/source/dshare/get_ixname.f90 @@ -95,6 +95,7 @@ function get_ixdecisions(varName) case('subRouting' ); get_ixdecisions=iLookDECISIONS%subRouting ! choice of method for sub-grid routing case('snowDenNew' ); get_ixdecisions=iLookDECISIONS%snowDenNew ! choice of method for new snow density case('snowUnload' ); get_ixdecisions=iLookDECISIONS%snowUnload ! choice of parameterization for snow unloading from canopy + case('prPotholes' ); get_ixdecisions=iLookDECISIONS%prPotholes ! choice of parameterization for prairie potholes ! get to here if cannot find the variable case default get_ixdecisions = integerMissing @@ -403,6 +404,9 @@ function get_ixparam(varName) case('zmaxLayer2_upper' ); get_ixparam = iLookPARAM%zmaxLayer2_upper ! maximum layer depth for the 2nd layer when > 2 layers (m) case('zmaxLayer3_upper' ); get_ixparam = iLookPARAM%zmaxLayer3_upper ! maximum layer depth for the 3rd layer when > 3 layers (m) case('zmaxLayer4_upper' ); get_ixparam = iLookPARAM%zmaxLayer4_upper ! maximum layer depth for the 4th layer when > 4 layers (m) + ! Oudin et al. (2005) PET parameters (for HDS) + case('oudinPETScaleK1' ); get_ixparam = iLookPARAM%oudinPETScaleK1 ! Oudin PET formula scaling factor (deg C) + case('oudinPETTempThrK2' ); get_ixparam = iLookPARAM%oudinPETTempThrK2 ! Oudin PET formula temperature threshold (deg C) ! get to here if cannot find the variable case default get_ixparam = integerMissing @@ -861,6 +865,12 @@ function get_ixbpar(varName) ! sub-grid routing case('routingGammaShape' ); get_ixbpar = iLookBPAR%routingGammaShape ! shape parameter in Gamma distribution used for sub-grid routing (-) case('routingGammaScale' ); get_ixbpar = iLookBPAR%routingGammaScale ! scale parameter in Gamma distribution used for sub-grid routing (s) + ! parameters for HDS pothole storage + case('depressionDepth' ); get_ixbpar = iLookBPAR%depressionDepth ! average depth of depressional storage (depressionVol/depressionArea) (m) + case('depressionAreaFrac' ); get_ixbpar = iLookBPAR%depressionAreaFrac ! fractional depressional area (depressionArea/basinArea) (-) + case('depressionCatchAreaFrac'); get_ixbpar = iLookBPAR%depressionCatchAreaFrac ! fractional area (of the landArea= basinArea - depressionArea) that drains to the depressions (-) + case('depression_p' ); get_ixbpar = iLookBPAR%depression_p ! shape of the slope profile (-) + case('depression_b' ); get_ixbpar = iLookBPAR%depression_b ! shape of contributing fraction curve (-) ! get to here if cannot find the variable case default get_ixbpar = integerMissing @@ -895,6 +905,15 @@ function get_ixbvar(varName) case('routingFractionFuture' ); get_ixbvar = iLookBVAR%routingFractionFuture ! fraction of runoff in future time steps (-) case('averageInstantRunoff' ); get_ixbvar = iLookBVAR%averageInstantRunoff ! instantaneous runoff (m s-1) case('averageRoutedRunoff' ); get_ixbvar = iLookBVAR%averageRoutedRunoff ! routed runoff (m s-1) + ! variables for pothole storage (HDS) + case('vMin' ); get_ixbvar = iLookBVAR%vMin ! volume of water in the meta depression at the start of a fill period (m3) + case('depConAreaFrac' ); get_ixbvar = iLookBVAR%depConAreaFrac ! contributing area fraction of pothole dominated areas [-] + case('basinConAreaFrac' ); get_ixbvar = iLookBVAR%basinConAreaFrac ! contributing area fraction per the entire subbasin from pothole and non-pothole areas [-] + case('pondVolFrac' ); get_ixbvar = iLookBVAR%pondVolFrac ! fractional pond volume = pondVol/depressionVol (-) + case('pondVol' ); get_ixbvar = iLookBVAR%pondVol ! pond volume at the end of time step (m3) + case('pondArea' ); get_ixbvar = iLookBVAR%pondArea ! pond area at the end of the time step (m2) + case('pondOutflow' ); get_ixbvar = iLookBVAR%pondOutflow ! pond outflow (m3) + case('pondEvap' ); get_ixbvar = iLookBVAR%pondEvap ! pond evaporation (kg m-2 s-1) ! get to here if cannot find the variable case default get_ixbvar = integerMissing diff --git a/build/source/dshare/popMetadat.f90 b/build/source/dshare/popMetadat.f90 index 7d741d292..38c8d3bbb 100755 --- a/build/source/dshare/popMetadat.f90 +++ b/build/source/dshare/popMetadat.f90 @@ -284,6 +284,9 @@ subroutine popMetadat(err,message) mpar_meta(iLookPARAM%zmaxLayer2_upper) = var_info('zmaxLayer2_upper' , 'maximum layer depth for the 2nd layer when > 2 layers' , 'm' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) mpar_meta(iLookPARAM%zmaxLayer3_upper) = var_info('zmaxLayer3_upper' , 'maximum layer depth for the 3rd layer when > 3 layers' , 'm' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) mpar_meta(iLookPARAM%zmaxLayer4_upper) = var_info('zmaxLayer4_upper' , 'maximum layer depth for the 4th layer when > 4 layers' , 'm' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + ! + mpar_meta(iLookPARAM%oudinPETScaleK1) = var_info('oudinPETScaleK1' , 'Oudin PET formula scaling factor' , 'deg C' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + mpar_meta(iLookPARAM%oudinPETTempThrK2) = var_info('oudinPETTempThrK2' , 'Oudin PET formula temperature threshold' , 'deg C' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) ! ----- ! * basin parameter data... @@ -293,6 +296,15 @@ subroutine popMetadat(err,message) bpar_meta(iLookBPAR%basin__aquiferBaseflowExp) = var_info('basin__aquiferBaseflowExp', 'baseflow exponent for the big bucket' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) bpar_meta(iLookBPAR%routingGammaShape) = var_info('routingGammaShape' , 'shape parameter in Gamma distribution used for sub-grid routing', '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) bpar_meta(iLookBPAR%routingGammaScale) = var_info('routingGammaScale' , 'scale parameter in Gamma distribution used for sub-grid routing', 's' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + + ! ----- + ! * HDS pothole storage parameters... + ! ----------------------------------- + bpar_meta(iLookBPAR%depressionDepth) = var_info('depressionDepth' , 'average depth of depressional storage (depressionVol/depressionArea)' , 'm' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bpar_meta(iLookBPAR%depressionAreaFrac) = var_info('depressionAreaFrac' , 'fractional depressional area (depressionArea/basinArea)' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bpar_meta(iLookBPAR%depressionCatchAreaFrac) = var_info('depressionCatchAreaFrac' , 'fractional area of the landArea (basinArea - depressionArea) that drains to the depressions', '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bpar_meta(iLookBPAR%depression_p) = var_info('depression_p' , 'shape of the slope profile' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bpar_meta(iLookBPAR%depression_b) = var_info('depression_b' , 'shape of contributing fraction curve' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) ! ----- ! * local model prognostic (state) variables... @@ -596,6 +608,18 @@ subroutine popMetadat(err,message) bvar_meta(iLookBVAR%routingFractionFuture) = var_info('routingFractionFuture' , 'fraction of runoff in future time steps' , '-' , get_ixVarType('routing'), iMissVec, iMissVec, .false.) bvar_meta(iLookBVAR%averageInstantRunoff) = var_info('averageInstantRunoff' , 'instantaneous runoff' , 'm s-1' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) bvar_meta(iLookBVAR%averageRoutedRunoff) = var_info('averageRoutedRunoff' , 'routed runoff' , 'm s-1' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + + ! ----- + ! * basin-wide HDS pothole storage fluxes/variables... + ! ----------------------------------------- + bvar_meta(iLookBVAR%vMin) = var_info('vMin' , 'volume of water in the meta depression at the start of a fill period' , 'm3' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%depConAreaFrac) = var_info('depConAreaFrac' , 'contributing area fraction of pothole dominated areas' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%basinConAreaFrac) = var_info('basinConAreaFrac' , 'contributing area fraction per the entire subbasin from pothole and non-pothole areas' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%pondVolFrac) = var_info('pondVolFrac' , 'fractional pond volume at the end of time step' , '-' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%pondVol) = var_info('pondVol' , 'pond volume at the end of time step' , 'm3' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%pondArea) = var_info('pondArea' , 'pond area at the end of the time step' , 'm2' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%pondOutflow) = var_info('pondOutflow' , 'pond outflow' , 'm3' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) + bvar_meta(iLookBVAR%pondEvap) = var_info('pondEvap' , 'pond evaporation' , 'kg m-2 s-1' , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.) ! ----- ! * model indices... diff --git a/build/source/dshare/var_lookup.f90 b/build/source/dshare/var_lookup.f90 index 73d7d93f7..43f426e82 100755 --- a/build/source/dshare/var_lookup.f90 +++ b/build/source/dshare/var_lookup.f90 @@ -71,6 +71,7 @@ MODULE var_lookup integer(i4b) :: spatial_gw = integerMissing ! choice of method for spatial representation of groundwater integer(i4b) :: subRouting = integerMissing ! choice of method for sub-grid routing integer(i4b) :: snowDenNew = integerMissing ! choice of method for new snow density + integer(i4b) :: prPotholes = integerMissing ! choice of method for prairie potholes representation endtype iLook_decision ! *********************************************************************************************************** @@ -304,6 +305,9 @@ MODULE var_lookup integer(i4b) :: zmaxLayer2_upper = integerMissing ! maximum layer depth for the 2nd layer when > 2 layers (m) integer(i4b) :: zmaxLayer3_upper = integerMissing ! maximum layer depth for the 3rd layer when > 3 layers (m) integer(i4b) :: zmaxLayer4_upper = integerMissing ! maximum layer depth for the 4th layer when > 4 layers (m) + ! Oudin et al. (2005) PET parameters (for HDS) + integer(i4b) :: oudinPETScaleK1 = integerMissing ! Oudin PET formula scaling factor (deg C) + integer(i4b) :: oudinPETTempThrK2 = integerMissing ! Oudin PET formula temperature threshold (deg C) endtype ilook_param ! *********************************************************************************************************** @@ -690,6 +694,12 @@ MODULE var_lookup ! within-grid routing integer(i4b) :: routingGammaShape = integerMissing ! shape parameter in Gamma distribution used for sub-grid routing (-) integer(i4b) :: routingGammaScale = integerMissing ! scale parameter in Gamma distribution used for sub-grid routing (s) + ! define paramters for HDS pothole storage + integer(i4b) :: depressionDepth = integerMissing ! average depth of depressional storage (depressionVol/depressionArea) (m) + integer(i4b) :: depressionAreaFrac = integerMissing ! fractional depressional area (depressionArea/basinArea) (-) + integer(i4b) :: depressionCatchAreaFrac = integerMissing ! fractional area (of the landArea = basinArea - depressionArea) that drains to the depressions (-) + integer(i4b) :: depression_p = integerMissing ! shape of the slope profile (-) + integer(i4b) :: depression_b = integerMissing ! shape of contributing fraction curve (-) endtype iLook_bpar ! *********************************************************************************************************** @@ -712,6 +722,15 @@ MODULE var_lookup integer(i4b) :: routingFractionFuture = integerMissing ! fraction of runoff in future time steps (-) integer(i4b) :: averageInstantRunoff = integerMissing ! instantaneous runoff (m s-1) integer(i4b) :: averageRoutedRunoff = integerMissing ! routed runoff (m s-1) + ! define variables for pothole storage (HDS) + integer(i4b) :: vMin = integerMissing ! volume of water in the meta depression at the start of a fill period (m3) + integer(i4b) :: depConAreaFrac = integerMissing ! contributing area fraction of pothole dominated areas [-] + integer(i4b) :: basinConAreaFrac = integerMissing ! contributing area fraction per the entire subbasin from pothole and non-pothole areas [-] + integer(i4b) :: pondVolFrac = integerMissing ! fractional pond volume at the end of time step (-) + integer(i4b) :: pondVol = integerMissing ! pond volume at the end of time step (m3) + integer(i4b) :: pondArea = integerMissing ! pond area at the end of the time step (m2) + integer(i4b) :: pondOutflow = integerMissing ! pond outflow (m3) + integer(i4b) :: pondEvap = integerMissing ! pond evaporation (kg m-2 s-1) endtype iLook_bvar ! *********************************************************************************************************** @@ -764,7 +783,7 @@ MODULE var_lookup type(iLook_decision),public,parameter :: iLookDECISIONS=iLook_decision( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,& 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,& 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,& - 31, 32, 33, 34, 35, 36, 37, 38) + 31, 32, 33, 34, 35, 36, 37, 38, 39) ! named variables: model time type(iLook_time), public,parameter :: iLookTIME =iLook_time ( 1, 2, 3, 4, 5, 6, 7) @@ -796,7 +815,7 @@ MODULE var_lookup 121,122,123,124,125,126,127,128,129,130,& 131,132,133,134,135,136,137,138,139,140,& 141,142,143,144,145,146,147,148,149,150,& - 151,152,153,154,155,156,157,158,159) + 151,152,153,154,155,156,157,158,159,160,161) ! named variables: model prognostic (state) variables type(iLook_prog), public,parameter :: iLookPROG =iLook_prog ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,& @@ -838,12 +857,13 @@ MODULE var_lookup 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,& 51, 52, 53, 54, 55, 56, 57, 58, 59, 60) - ! named variables: basin-average parameters - type(iLook_bpar), public,parameter :: iLookBPAR =ilook_bpar ( 1, 2, 3, 4, 5) + ! named variables: basin-average parameters (including HDS parameters) + type(iLook_bpar), public,parameter :: iLookBPAR =ilook_bpar ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) - ! named variables: basin-average variables + ! named variables: basin-average variables (including HDS variables) type(iLook_bvar), public,parameter :: iLookBVAR =ilook_bvar ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,& - 11, 12, 13) + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,& + 21) ! named variables in varibale type structure type(iLook_varType), public,parameter :: iLookVarType =ilook_varType ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,& diff --git a/build/source/engine/HDS.f90 b/build/source/engine/HDS.f90 new file mode 100755 index 000000000..f7476672b --- /dev/null +++ b/build/source/engine/HDS.f90 @@ -0,0 +1,366 @@ +module HDS + USE nrtype + ! physical constants + USE multiconst,only:iden_water ! intrinsic density of water (kg m-3) + + ! define data types + USE data_types,only:& + ! no spatial dimension + var_dlength, & ! x%var(:)%dat (dp) + var_d ! var(:) ! for GRU parameters + + ! provide access to the named variables that describe elements of parameter structures + USE var_lookup,only:iLookBVAR ! look-up values for basin-average model variables + USE var_lookup,only:iLookBPAR ! look-up values for basin-average model parameters (for HDS) + + contains + + subroutine init_summa_HDS(bvarData, bparGRU) + ! initialize HDS pothole storage model states for SUMMA + implicit none + !subroutine arguments + type(var_dlength) , intent(inout) :: bvarData ! x%var(:)%dat -- basin-average variables + type(var_d) , intent(in) :: bparGRU ! x%gru(:)%var(:) -- basin-average parameters + + ! local variables + real(rkind) :: depressionArea ! depression area [m2] + real(rkind) :: depressionVol ! depression volume [m3] + !***************************************** + ! initialize HDS pothole storage variables + ! GRU level + !***************************************** + ! start association + associate(pondVolFrac => bvarData%var(iLookBVAR%pondVolFrac)%dat(1) , & ! depression volume fraction [-] + depressionDepth => bparGRU%var(iLookBPAR%depressionDepth) , & ! depression depth [m] + depressionAreaFrac => bparGRU%var(iLookBPAR%depressionAreaFrac) , & ! depression area fraction [-] + totalArea => bvarData%var(iLookBVAR%basin__totalArea)%dat(1) , & ! basin (GRU) total area [m2] + p => bparGRU%var(iLookBPAR%depression_p) , & ! shape of the slope profile [-]. Exponent for calculating the fractional wet area + pondVol => bvarData%var(iLookBVAR%pondVol)%dat(1) , & ! pond volume [m3] + pondArea => bvarData%var(iLookBVAR%pondArea)%dat(1) , & ! pond Area [m2] + depConAreaFrac => bvarData%var(iLookBVAR%depConAreaFrac)%dat(1) , & ! fractional contributing area [-] + vMin => bvarData%var(iLookBVAR%vMin)%dat(1) & ! minimum pond volume [m3] + ) + + pondVolFrac = max(pondVolFrac, zero) + depressionArea = depressionAreaFrac * totalArea + depressionVol = depressionDepth * depressionArea + ! the following checks (<0) are used to set initial values for the variables if there are set to missing value (-999) + ! the following calculations are approximate solutions and will be updated inside the routine + if(pondVol < zero ) pondVol = pondVolFrac * depressionVol ! approximate vol calculation, will be updated later in the subroutine + if(pondArea < zero ) pondArea = depressionArea*((pondVol/depressionVol)**(two/(p + two))) + if(depConAreaFrac < zero) depConAreaFrac = pondVol/depressionVol ! assume that contrib_frac = vol_frac_sml for initialization purposes + if(vMin < zero ) vMin = pondVol + end associate + end subroutine init_summa_HDS + !============================================================= + !============================================================= + subroutine init_pond_Area_Volume(depArea, depVol, totEvap, volFrac, p, pondVol, pondArea) + ! used to initialize pond volume and area (Not used in LSMs) + ! initialize at capacity minus depth of evaporation + implicit none + !subroutine arguments + real(rkind), intent(in) :: depArea ! depression area [m2] + real(rkind), intent(in) :: depVol ! depression volume [m3] + real(rkind), intent(in) :: totEvap ! total evaporation [m] to be subtracted from the ponded water + real(rkind), intent(in) :: volFrac ! volume fraction [-]. Used to fill the depressions using percentage (i.e., 50% full) + real(rkind), intent(in) :: p ! shape of the slope profile [-]. Exponent for calculating the fractional wet area + real(rkind), intent(out) :: pondVol, pondArea ! pond volume [m3], pond area [m2] + ! local variables + real(rkind) :: depHeight ! depression height [m] + real(rkind) :: pondHeight ! height of the water level in depression [m] + + ! estimate the maximum depth of the depression + depHeight = depVol*(one + two/p)/depArea ! depression height [m] + + !subtract ET from ponded water if ET > 0 + if(totEvap > zero) then + + ! estimate the height of the water level after evaporation + pondHeight = max(depHeight - totEvap, zero) + + ! get the volume of water in the depression + pondVol = depVol*((pondHeight/depHeight)**(one + two/p)) + + else !use volume frac to get initial pondVol + pondVol = depVol*volFrac + end if + + ! calculate pondArea + pondArea = depArea*((pondVol/depVol)**(two/(p + two))) + + end subroutine init_pond_Area_Volume + + !============================================================= + !============================================================= + subroutine runDepression(bvarData, & ! basin-average variables + bparGRU, & ! basin-average parameters + basinPrecip, & ! average basin precipitation amount (kg m-2 s-1 = mm s-1) + basinPotentialEvap) ! average basin potential evaporation amount (mm s-1) + + ! global data + USE globalData,only:data_step ! time step of forcing data (s) - used by HDS to accumulate fluxes at each time step (forcing data step) + ! Estimate the volumetric storage at the end of the time interval using the numerical solution + implicit none + ! subroutine arguments + type(var_dlength) , intent(inout) :: bvarData ! x%var(:)%dat -- basin-average variables + type(var_d) , intent(in) :: bparGRU ! x%gru(:)%var(:) -- basin-average parameters + real(rkind) , intent(in) :: basinPrecip ! average basin precipitation amount (kg m-2 s-1 = mm s-1) + real(rkind) , intent(in) :: basinPotentialEvap ! average basin potential evaporation amount (mm s-1) + ! local variables -- subroutine parameters (static or deactivated) + real(rkind) , parameter :: dt=1._rkind ! model time step [length of timestep = 1 -- nosubsteps] + real(rkind) , parameter :: tau=0._rkind ! model parameters = tau time constant linear reservoir [timestep-1] ! currently deactivated + ! local variables -- general + real(rkind) :: qSeas, pRate, etPond ! forcing data = runoff, precipitation, ET [mm/timestep] + real(rkind) :: depArea, depVol, depUpsArea, landArea ! spatial attributes = depression area [m2], depression volume [m3], upstream area of depressions [m2], land area = total area - depression area (m2) + real(rkind) :: rvrUpsArea ! spatial attributes = upland area contributing to the river from non-pothole areas [m2] + ! local variables -- model decisions + integer(i4b), parameter :: implicitEuler=1001 ! named variable for the implicit Euler solution + integer(i4b), parameter :: shortSubsteps=1002 ! named variable for the short substeps solution + integer(i4b) :: solution ! numerical solution method + ! local variables --- implicit Euler solution + integer(i4b) :: iter ! iteration counter + integer(i4b) , parameter :: nIter=100 ! number of iterations (algorithm control parameter) + real(rkind) , parameter :: xConv=1.e-6_rkind ! convergence criteria (algorithm control parameter) + real(rkind) :: xMin, xMax ! min and max storage + real(rkind) :: xRes ! residual + logical(lgt) :: failure ! failure flag + ! local variables -- short substeps solution + integer(i4b) , parameter :: nSub=100 ! number of substeps (algorithm control parameter) + integer(i4b) :: iSub ! counter for shortSubsteps solution + real(rkind) :: dtSub ! dt for shortSubsteps solution + ! local variables -- diagnostic variables (model physics) + real(rkind) :: Q_di, Q_det, Q_dix, Q_do ! fluxes [L3 T-1]: sum of water inputs to the pond, evapotranspiration, infiltration, pond outflow + real(rkind) :: Q_det_adj, Q_dix_adj ! adjusted evapotranspiration & infiltration fluxes [L3 T-1] for mass balance closure (i.e., when losses > pondVol). Zero values mean no adjustment needed. + real(rkind) :: cFrac, g, dgdv ! contributing fraction, net fluxes and derivative + real(rkind) :: xVol ! pond volume at the end of the time step [m3] + + ! start the association + associate(totalArea => bvarData%var(iLookBVAR%basin__totalArea)%dat(1) , & ! basin total area (m2) + basinTotalRunoff => bvarData%var(iLookBVAR%basin__TotalRunoff)%dat(1) , & ! basin total runoff (m s-1) + ! HDS pothole storage variables + vMin => bvarData%var(iLookBVAR%vMin)%dat(1) , & ! volume of water in the meta depression at the start of a fill period (m3) + fArea => bvarData%var(iLookBVAR%depConAreaFrac)%dat(1) , & ! fractional contributing areas of pothole-dominated areas (-) + basinConAreaFrac => bvarData%var(iLookBVAR%basinConAreaFrac)%dat(1) , & ! contributing area fraction per the entire subbasin from pothole and non-pothole areas [-] + fVol => bvarData%var(iLookBVAR%pondVolFrac)%dat(1) , & ! fractional pond volume = pondVol/depressionVol (-) + pondVol => bvarData%var(iLookBVAR%pondVol)%dat(1) , & ! pond volume at the end of time step (m3) + pondArea => bvarData%var(iLookBVAR%pondArea)%dat(1) , & ! pond area at the end of the time step (m2) + pondOutflow => bvarData%var(iLookBVAR%pondOutflow)%dat(1) , & ! pond outflow (m3) + pondEvap => bvarData%var(iLookBVAR%pondEvap)%dat(1) , & ! pond evaporation (kg m-2 s-1) + ! HDS pothole storage parameters + depDepth => bparGRU%var(iLookBPAR%depressionDepth) , & ! depression depth (m) + depAreaFrac => bparGRU%var(iLookBPAR%depressionAreaFrac) , & ! fractional depressional area (depressionArea/basinArea) (-) + depCatchAreaFrac => bparGRU%var(iLookBPAR%depressionCatchAreaFrac) , & ! fractional area (of the landArea= basinArea - depressionArea) that drains to the depressions (-) + p => bparGRU%var(iLookBPAR%depression_p) , & ! shape of the slope profile (-) + b => bparGRU%var(iLookBPAR%depression_p) & ! shape of contributing fraction curve (-) + ) + ! exit the subroutine if depDepth for this GRU is <= 0 (i.e., this GRU does not have depressions) + if(depDepth <= zero) return + + ! convert fluxes to mm/timestep + qSeas = max(basinTotalRunoff, zero) * 1000.0_rkind * data_step ! runoff [m s-1] -> [mm/timestep] + pRate = basinPrecip * data_step ! precipitation [mm s-1] -> [mm/timestep] + etPond = basinPotentialEvap * data_step ! potential evaporation [mm s-1] -> [mm/timestep] + + ! define numerical solution + solution = implicitEuler + !solution = shortSubsteps + + ! initialize pond volume + xVol = pondVol + ! initialize corrected evaporation and infiltration for mass balance closure (not used here) + Q_det_adj = zero + Q_dix_adj = zero + ! initialize pondOutflow + pondOutflow = 0._rkind + ! calculate some spatial attributes + depArea = depAreaFrac * totalArea + depVol = depDepth * depArea + landArea = totalArea - depArea + depUpsArea = max(landArea * depCatchAreaFrac, 0._rkind) + rvrUpsArea = landArea * (1._rkind - depCatchAreaFrac) + + ! ---------- option 1: implicit Euler ---------- + + ! check if implicit Euler is desired + if(solution == implicitEuler)then + + ! initialize brackets + xMin = zero + xMax = depVol + + ! iterate (can start with 1 because the index is not used) + do iter=1,nIter + + ! compute fluxes and derivatives + call computFlux(iter, xVol, qSeas, pRate, etPond, depArea, depVol, depUpsArea, p, tau, b, vMin, Q_di, Q_det, Q_dix, Q_do, cFrac, g, dgdv) + + ! compute residual + xRes = (pondVol + g*dt) - xVol + + ! check convergence + if(iter > 1 .and. abs(xRes) < xConv)then + failure=.false. + exit + endif + + ! update constraints + if(xRes > zero) xMin=xVol + if(xRes < zero) xMax=xVol + + ! special case where xMax is too small + if(xRes > zero .and. xVol > 0.99*xMax) xMax = min(xMax*10.0_rkind, depVol) + + ! update state (pondVol) + xVol = xVol + xRes / (one - dgdv*dt) + + ! use bi-section if violated constraints + if(xVol < xMin .or. xVol > xMax) xVol=(xMin+xMax)/2.0_rkind + + ! assign failure + failure = (iter == nIter) + + enddo ! iterating + + endif ! if implicit Euler + + ! ---------- option 2: short substeps ---------- + + ! check if short substeps are desired + if(solution == shortSubsteps .or. failure)then + + ! define length of the substeps + dtSub = one / real(nSub, kind(rkind)) + + ! loop through substeps + do iSub=1,nSub + call computFlux(iter, xVol, qSeas, pRate, etPond, depArea, depVol, depUpsArea, p, tau, b, vMin, Q_di, Q_det, Q_dix, Q_do, cFrac, g, dgdv) + xVol = xVol + g*dtSub + + ! prevent -ve ponVol values to avoid nans + if(xVol < zero)then + xVol = zero + ! adjust evaporation and infiltration fluxes proprtionally (using weights) to close mass balance (i.e., losses = PondVol) + Q_det_adj = pondVol * (Q_det/(Q_det+Q_dix)) + Q_dix_adj = pondVol * (Q_dix/(Q_det+Q_dix)) + !break the loop as there's no need to continue xvol depVol)then + cFrac = 1._rkind + ! assign excess water (xVol-depVol) as outflow + Q_do = Q_do + (xVol-depVol) + xVol = depVol + !break the loop as there's no need to continue xvol m s-1 + (basinTotalRunoff * rvrUpsArea) / totalArea ! m s-1 -> m3 s-1 -> m s-1 + end associate + end subroutine runDepression + + !============================================================= + !============================================================= + subroutine computFlux(iter, pondVol, & ! input: iteration index, state variable = pond volume [m3] + qSeas, pRate, etPond, & ! input: forcing data = runoff, precipitation, ET [mm/day] + depArea, depVol, upsArea, & ! input: spatial attributes = depression area [m2], depression volume [m3], upstream area [m2] + p, tau, & ! input: model parameters = p [-] shape of the slope profile; tau [day-1] time constant linear reservoir + b, vmin, & ! input: model parameters = b [-] shape of contributing fraction curve; vmin [m3] minimum volume + Q_di, Q_det, Q_dix, Q_do, & ! output: individual model fluxes + cFrac, g, dgdv) ! output: contributing fraction, net fluxes and derivative + implicit none + ! subroutine arguments + integer, intent(in) :: iter ! iteration index + real(rkind), intent(in) :: pondVol ! state variable = pond volume [m3] + real(rkind), intent(in) :: qSeas, pRate, etPond ! input: forcing data = runoff, precipitation, ET [mm/day] + real(rkind), intent(in) :: depArea, depVol, upsArea ! input: spatial attributes = depression area [m2], depression volume [m3], upstream area [m2] + real(rkind), intent(in) :: p, tau ! input: model parameters = p [-] shape of the slope profile; tau [day-1] time constant linear reservoir + real(rkind), intent(in) :: b ! input: model parameters = b [-] shape of contributing fraction curve; vmin [m3] minimum volume + real(rkind), intent(inout) :: vmin ! in/out: vmin [m3] minimum volume + real(rkind), intent(out) :: Q_di, Q_det, Q_dix, Q_do ! output: individual model fluxes + real(rkind), intent(out) :: cFrac, g, dgdv ! output: contributing fraction, net fluxes and derivative, pond area + ! local variables + real(rkind), parameter :: ms=0.0001_rkind ! smoothing parameter (algorithm control) + real(rkind), parameter :: verySmall = 1.0e-12_rkind ! very small value + real(rkind), parameter :: rCoef=zero ! runoff coefficient [-] !HDS_standalone 0.050_rkind + real(rkind) :: pondArea ! calculated pond area + real(rkind) :: pInput ! precipitation (or rain+melt) [m/day] + real(rkind) :: qInput ! surface runoff [m/day] + real(rkind) :: eLosses ! evaporation losses [m/day] + real(rkind) :: vPrime ! smoothed pondVol value for Euler solution + real(rkind) :: vTry ! adjusted pondVol for derivatives calculations + real(rkind) :: dadv, dpdv, dfdv, didv ! derivatives + + ! compute the pond area + pondArea = depArea*((pondVol/depVol)**(two/(p + two))) + + ! get the forcing + pInput = pRate /iden_water ! kg m-2 -> m s-1 + qInput = qSeas /iden_water + rCoef*pRate/iden_water ! surface runoff ! kg m-2 s-1 -> m s-1 + eLosses = etPond/iden_water ! evaporation losses kg m-2 s-1 -> m s-1 + + ! get volume fluxes from the host land model + ! sum of water input to the depression (eq 11) + Q_di = upsArea*qInput + (depArea - pondArea)*qInput + pondArea*pInput + + ! evapotranspiration losses (eq 12) + Q_det = pondArea*eLosses + + ! compute infiltration from the bottom of the pond (eq 13) + Q_dix = tau*pondVol + + ! update the minimum parameter (force hysteresis) + if(iter == 1) then + if(Q_di < Q_det+Q_dix .and. pondVol < depVol-verySmall) vMin = pondVol + endif + + ! compute the outflow from the meta depression + ! smoothing pondVol calculation ! (eq 24) + vPrime = half*(vMin + pondVol + sqrt((pondVol-vMin)**two + ms)) ! (eq 24) + vPrime = min(vPrime, depVol) !MIA -> limit vPrime to be <= depVol + ! calculate contributing fraction !(eq 25) + cFrac = one - ((depVol - vPrime)/(depVol - vMin))**b !(eq 25) + cFrac = max(cFrac, zero ) !prevent -ve values MIA + ! calculate outlfow (eq 16) + Q_do = cFrac*Q_di + + ! compute the net flux (eq 10) + g = Q_di - Q_det - Q_dix - Q_do + + ! compute derivate in pond area w.r.t. pond volume + vTry = max(verySmall, pondVol) ! avoid divide by zero + dadv = ((two *depArea)/((p + two )*depVol)) * ((vTry/depVol)**(-p/(p + two ))) + + ! compute derivatives in volume fluxes w.r.t. pond volume (meta depression) + dpdv = half*((pondVol-vMin)/sqrt((pondVol-vMin)**two + ms) + one ) ! max smoother + dfdv = dpdv*(b*(depVol - pondVol)**(b - one ))/((depVol - vMin)**b) ! contributing fraction (note chain rule) + didv = dadv*(pInput - qInput) + dgdv = didv*(one - cFrac) - dfdv*Q_di - dadv*eLosses - tau + + end subroutine computFlux + +end module HDS diff --git a/build/source/engine/derivforce.f90 b/build/source/engine/derivforce.f90 index caa7e617d..a5920f8f3 100644 --- a/build/source/engine/derivforce.f90 +++ b/build/source/engine/derivforce.f90 @@ -57,7 +57,7 @@ module derivforce_module ! privacy implicit none private -public::derivforce +public::derivforce, calcPotentialEvap_Oudin2005 contains ! ************************************************************************************************ @@ -354,5 +354,33 @@ subroutine derivforce(time_data,forc_data,attr_data,mpar_data,prog_data,diag_dat end subroutine derivforce + ! ************************************************************************************************ + ! function calcPotentialEvap_Oudin2005: compute potential evaporation for HDS pond calculations + ! ************************************************************************************************ + !============================================================= + function calcPotentialEvap_Oudin2005(SWRadAtm, airtemp, K1, K2) result(potentialEvap) + ! calculate potential evaporation using Oudin et al. (2005)'s formula for a specific HRU + USE multiconst,only:LH_vap, & ! latent heat of vaporization (J kg-1) + iden_water ! intrinsic density of water (kg m-3) + implicit none + !function arguments + real(rkind), intent(in) :: SWRadAtm ! downwelling shortwave radiaiton [w/m2] + real(rkind), intent(in) :: airtemp ! air temperature [K] + real(rkind), intent(in) :: K1 ! scaling factor (deg C) + real(rkind), intent(in) :: K2 ! minimum value of air temperature for which PE is not zero (deg C) + + ! local variables + real(rkind) :: potentialEvap ! pontentail evaporation as calculated by Oudin's formula (mm s-1) + + ! Oudin (2005)'s formula + potentialEvap = (1000._rkind * & + (SWRadAtm * 1e-6 / & ! w/m2 to MJ/m2/s + (LH_vap * 1e-6 * iden_water)) * & ! J kg-1 to MJ kg-1 + ((airtemp - 273.15_rkind + K2)/K1)) ! K to deg C + + ! check for negative values + potentialEvap = max(potentialEvap, zero) + +end function calcPotentialEvap_Oudin2005 end module derivforce_module diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90 index b658d7c0d..50b939536 100755 --- a/build/source/engine/mDecisions.f90 +++ b/build/source/engine/mDecisions.f90 @@ -146,6 +146,9 @@ module mDecisions_module ! look-up values for the choice of snow unloading from the canopy integer(i4b),parameter,public :: meltDripUnload = 321 ! Hedstrom and Pomeroy (1998), Storck et al 2002 (snowUnloadingCoeff & ratioDrip2Unloading) integer(i4b),parameter,public :: windUnload = 322 ! Roesch et al 2001, formulate unloading based on wind and temperature +! look-up values for the choice of prairie pothole parameterization +integer(i4b),parameter,public :: noPotholes = 331 ! no explicit pothole parameterization +integer(i4b),parameter,public :: HDSmodel = 332 ! Hysteretic Depressional Storage model implementation ! ----------------------------------------------------------------------------------------------------------- contains @@ -623,11 +626,20 @@ subroutine mDecisions(err,message) ! choice of snow unloading from canopy select case(trim(model_decisions(iLookDECISIONS%snowUnload)%cDecision)) case('meltDripUnload','notPopulatedYet'); model_decisions(iLookDECISIONS%snowUnload)%iDecision = meltDripUnload ! Hedstrom and Pomeroy (1998), Storck et al 2002 (snowUnloadingCoeff & ratioDrip2Unloading) - case('windUnload'); model_decisions(iLookDECISIONS%snowUnload)%iDecision = windUnload ! Roesch et al 2001, formulate unloading based on wind and temperature + case('windUnload'); model_decisions(iLookDECISIONS%snowUnload)%iDecision = windUnload ! Roesch et al 2001, formulate unloading based on wind and temperature case default err=10; message=trim(message)//"unknown option for snow unloading [option="//trim(model_decisions(iLookDECISIONS%snowUnload)%cDecision)//"]"; return end select + ! choice of prairie potholes + ! NOTE: use noPotholes as the default, where prairie pothole method is undefined (not populated yet) + select case(trim(model_decisions(iLookDECISIONS%prPotholes)%cDecision)) + case('noPotholes','notPopulatedYet'); model_decisions(iLookDECISIONS%prPotholes)%iDecision = noPotholes ! No representation of prairie potholes + case('HDSmodel'); model_decisions(iLookDECISIONS%prPotholes)%iDecision = HDSmodel ! Hysteretic Depressional Storage model implementation + case default + err=10; message=trim(message)//"unknown option for prairie potholes [option="//trim(model_decisions(iLookDECISIONS%prPotholes)%cDecision)//"]"; return + end select + ! ----------------------------------------------------------------------------------------------------------------------------------------------- ! check for consistency among options diff --git a/build/source/engine/nrtype.f90 b/build/source/engine/nrtype.f90 index c6ea58b9d..9f7083ace 100755 --- a/build/source/engine/nrtype.f90 +++ b/build/source/engine/nrtype.f90 @@ -26,4 +26,10 @@ MODULE nrtype real(rkind), parameter :: nr_quadMissing=-9999._qp ! missing quadruple precision number real(rkind), parameter :: nr_realMissing=-9999._rkind ! missing double precision number integer(i4b), parameter :: nr_integerMissing=-9999 ! missing integer + ! data types for HDS pothole storage + ! useful shortcuts + real(rkind), parameter :: zero = 0.0_rkind + real(rkind), parameter :: half = 0.5_rkind + real(rkind), parameter :: one = 1.0_rkind + real(rkind), parameter :: two = 2.0_rkind END MODULE nrtype diff --git a/build/source/engine/read_pinit.f90 b/build/source/engine/read_pinit.f90 index d0103d615..82a70b4a7 100755 --- a/build/source/engine/read_pinit.f90 +++ b/build/source/engine/read_pinit.f90 @@ -53,7 +53,7 @@ subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message) integer(i4b),intent(out) :: err ! error code character(*),intent(out) :: message ! error message ! define general variables - logical(lgt),parameter :: backwardsCompatible=.false. ! .true. if skip check that all parameters are populated + logical(lgt),parameter :: backwardsCompatible=.true. ! .true. if skip check that all parameters are populated character(len=256) :: cmessage ! error message for downwind routine character(LEN=256) :: infile ! input filename integer(i4b) :: unt ! file unit (free unit output from file_open) @@ -130,7 +130,7 @@ subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message) end if end do ! (looping through lines in the file) ! check we have populated all variables - ! NOTE: ultimately need a need a parameter dictionary to ensure that the parameters used are populated + ! NOTE: ultimately need a parameter dictionary to ensure that the parameters used are populated if(.not.backwardsCompatible)then ! if we add new variables in future versions of the code, then some may be missing in the input file if(any(parFallback(:)%default_val < 0.99_rkind*realMissing))then do ivar=1,size(parFallback) diff --git a/build/source/engine/run_oneGRU.f90 b/build/source/engine/run_oneGRU.f90 index 17fca209f..42f8928b9 100755 --- a/build/source/engine/run_oneGRU.f90 +++ b/build/source/engine/run_oneGRU.f90 @@ -33,10 +33,12 @@ module run_oneGRU_module ! no spatial dimension var_ilength, & ! x%var(:)%dat (i4b) var_dlength, & ! x%var(:)%dat (dp) + var_d, & ! var(:) ! for GRU parameters ! hru dimension hru_int, & ! x%hru(:)%var(:) (i4b) hru_int8, & ! x%hru(:)%var(:) integer(8) hru_double, & ! x%hru(:)%var(:) (dp) + ! gru_double, & ! x%gru(:)%var(:) (dp) hru_intVec, & ! x%hru(:)%var(:)%dat (i4b) hru_doubleVec ! x%hru(:)%var(:)%dat (dp) @@ -47,6 +49,9 @@ module run_oneGRU_module USE var_lookup,only:iLookINDEX ! look-up values for local column index variables USE var_lookup,only:iLookFLUX ! look-up values for local column model fluxes USE var_lookup,only:iLookBVAR ! look-up values for basin-average model variables +USE var_lookup,only:iLookBPAR ! look-up values for basin-average model parameters (for HDS) +USE var_lookup,only:iLookPARAM ! look-up values for HRU model parameters (for HDS PET calculations) +USE var_lookup, only:iLookFORCE ! look-up values for HRU forcing - used to estimate basin average forcing for HDS ! provide access to model decisions USE globalData,only:model_decisions ! model decision structure @@ -56,7 +61,8 @@ module run_oneGRU_module USE mDecisions_module,only:& ! look-up values for the choice of method for the spatial representation of groundwater localColumn, & ! separate groundwater representation in each local soil column singleBasin, & ! single groundwater store over the entire basin - bigBucket ! a big bucket (lumped aquifer model) + bigBucket, & ! a big bucket (lumped aquifer model) + HDSmodel ! Hysteretic Depressional Storage model implementation for prairie potholes ! ----------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------- @@ -82,6 +88,7 @@ subroutine run_oneGRU(& typeHRU, & ! intent(in): local classification of soil veg etc. for each HRU idHRU, & ! intent(in): local classification of hru and gru IDs attrHRU, & ! intent(in): local attributes for each HRU + bparGRU, & ! intent(in): local attributes for the GRU for HDS calculations ! data structures (input-output) mparHRU, & ! intent(inout): local model parameters indxHRU, & ! intent(inout): model indices @@ -97,6 +104,8 @@ subroutine run_oneGRU(& USE run_oneHRU_module,only:run_oneHRU ! module to run for one HRU USE qTimeDelay_module,only:qOverland ! module to route water through an "unresolved" river network + USE HDS,only:runDepression ! module to run HDS pothole storage dynamics + USE derivforce_module,only:calcPotentialEvap_Oudin2005 ! function to calculate potential evaporation based on Oudin et al. (2005)'s formula ! ----- define dummy variables ------------------------------------------------------------------------------------------ @@ -104,13 +113,14 @@ subroutine run_oneGRU(& ! model control type(gru2hru_map) , intent(inout) :: gruInfo ! HRU information for given GRU (# HRUs, #snow+soil layers) - real(rkind) , intent(inout) :: dt_init(:) ! used to initialize the length of the sub-step for each HRU + real(rkind) , intent(inout) :: dt_init(:) ! used to initialize the length of the sub-step for each HRU integer(i4b) , intent(inout) :: ixComputeVegFlux(:) ! flag to indicate if we are computing fluxes over vegetation (false=no, true=yes) ! data structures (input) integer(i4b) , intent(in) :: timeVec(:) ! integer vector -- model time data type(hru_int) , intent(in) :: typeHRU ! x%hru(:)%var(:) -- local classification of soil veg etc. for each HRU type(hru_int8) , intent(in) :: idHRU ! x%hru(:)%var(:) -- local classification of hru and gru IDs type(hru_double) , intent(in) :: attrHRU ! x%hru(:)%var(:) -- local attributes for each HRU + type(var_d) , intent(in) :: bparGRU ! x%gru(:)%var(:) -- basin-average parameters ! data structures (input-output) type(hru_doubleVec) , intent(inout) :: mparHRU ! x%hru(:)%var(:)%dat -- local (HRU) model parameters type(hru_intVec) , intent(inout) :: indxHRU ! x%hru(:)%var(:)%dat -- model indices @@ -132,8 +142,11 @@ subroutine run_oneGRU(& integer(i4b) :: nSnow ! number of snow layers integer(i4b) :: nSoil ! number of soil layers integer(i4b) :: nLayers ! total number of layers - real(rkind) :: fracHRU ! fractional area of a given HRU (-) + real(rkind) :: fracHRU ! fractional area of a given HRU (-) logical(lgt) :: computeVegFluxFlag ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow) + ! basin average fluxes for HDS -- local variables + real(rkind) :: basinPrecip ! average basin precipitation amount (kg m-2 s-1 = mm s-1) + real(rkind) :: basinPotentialEvap ! average basin potential evaporation amount (mm s-1) ! initialize error control err=0; write(message, '(A24,I0,A2)' ) 'run_oneGRU (gru index = ',gruInfo%gru_nc,')/' @@ -151,6 +164,10 @@ subroutine run_oneGRU(& bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) = 0._rkind ! baseflow from the aquifer (m s-1) bvarData%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = 0._rkind ! transpiration loss from the aquifer (m s-1) + ! initialize basin average forcing + basinPrecip = 0._rkind ! precipitation rate averaged over the basin + basinPotentialEvap = 0._rkind + ! initialize total inflow for each layer in a soil column do iHRU=1,gruInfo%hruCount fluxHRU%hru(iHRU)%var(iLookFLUX%mLayerColumnInflow)%dat(:) = 0._rkind @@ -248,19 +265,39 @@ subroutine run_oneGRU(& bvarData%var(iLookBVAR%basin__AquiferRecharge)%dat(1) = bvarData%var(iLookBVAR%basin__AquiferRecharge)%dat(1) + fluxHRU%hru(iHRU)%var(iLookFLUX%scalarSoilDrainage)%dat(1) * fracHRU bvarData%var(iLookBVAR%basin__AquiferTranspire)%dat(1) = bvarData%var(iLookBVAR%basin__AquiferTranspire)%dat(1) + fluxHRU%hru(iHRU)%var(iLookFLUX%scalarAquiferTranspire)%dat(1) * fracHRU - bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) = bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) & - + fluxHRU%hru(iHRU)%var(iLookFLUX%scalarAquiferBaseflow)%dat(1) * fracHRU + bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) = bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) + fluxHRU%hru(iHRU)%var(iLookFLUX%scalarAquiferBaseflow)%dat(1) * fracHRU end if ! averaging more fluxes (and/or states) can be added to this section as desired + ! averagin fluxes if HDS is active + if(model_decisions(iLookDECISIONS%prPotholes)%iDecision == HDSmodel)then + + associate(pptrate => forcHRU%hru(iHRU)%var(iLookFORCE%pptrate) , & + SWRadAtm => forcHRU%hru(iHRU)%var(iLookFORCE%SWRadAtm), & + airtemp => forcHRU%hru(iHRU)%var(iLookFORCE%airtemp) , & + ! parameters of Oudin PET formula + scaleK1 => mparHRU%hru(iHRU)%var(iLookPARAM%oudinPETScaleK1)%dat(1), & + tempThrK2 => mparHRU%hru(iHRU)%var(iLookPARAM%oudinPETTempThrK2)%dat(1)) + + basinPrecip = basinPrecip + pptrate * fracHRU + ! calculate potential evaporation using Oudin (2005)'s formula + ! check the parameters of Oudin PET formula + if(scaleK1 < 0._rkind .or. tempThrK2 < 0._rkind)then + write(message(len_trim(message) + 1:), '(A7,I0,A75)') 'hruId=',gruInfo%hruInfo(iHRU)%hru_id, '/oudinPETFormula/missing or negative values for ScaleK1 or TempThrK2' + err=20; return + endif + basinPotentialEvap = basinPotentialEvap + calcPotentialEvap_Oudin2005(SWRadAtm, airtemp, scaleK1, tempThrK2) * fracHRU + end associate + end if ! HDS control flag end do ! (looping through HRUs) ! *********************************************************************************************************************** ! ********** END LOOP THROUGH HRUS ************************************************************************************** ! *********************************************************************************************************************** - ! perform the routing - associate(totalArea => bvarData%var(iLookBVAR%basin__totalArea)%dat(1) ) + ! perform the pothole storage and routing + associate(totalArea => bvarData%var(iLookBVAR%basin__totalArea)%dat(1) , & + basinTotalRunoff => bvarData%var(iLookBVAR%basin__TotalRunoff)%dat(1)) ! basin total runoff (m s-1) ! compute water balance for the basin aquifer if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then @@ -271,16 +308,23 @@ subroutine run_oneGRU(& ! calculate total runoff depending on whether aquifer is connected if(model_decisions(iLookDECISIONS%groundwatr)%iDecision == bigBucket) then ! aquifer - bvarData%var(iLookBVAR%basin__TotalRunoff)%dat(1) = bvarData%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) + basinTotalRunoff = bvarData%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + bvarData%var(iLookBVAR%basin__AquiferBaseflow)%dat(1) else ! no aquifer - bvarData%var(iLookBVAR%basin__TotalRunoff)%dat(1) = bvarData%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + bvarData%var(iLookBVAR%basin__SoilDrainage)%dat(1) + basinTotalRunoff = bvarData%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) + bvarData%var(iLookBVAR%basin__ColumnOutflow)%dat(1)/totalArea + bvarData%var(iLookBVAR%basin__SoilDrainage)%dat(1) endif + ! *********************************************************************************************************************** + ! ********** PRAIRIE POTHOLE IMPLEMENTATION (HDS)************************************************************************ + ! *********************************************************************************************************************** + ! run the actual HDS depressional storage model for this GRU + if(model_decisions(iLookDECISIONS%prPotholes)%iDecision == HDSmodel) call runDepression(bvarData, bparGRU, basinPrecip, basinPotentialEvap) + ! *********************************************************************************************************************** + call qOverland(& ! input model_decisions(iLookDECISIONS%subRouting)%iDecision, & ! intent(in): index for routing method - bvarData%var(iLookBVAR%basin__TotalRunoff)%dat(1), & ! intent(in): total runoff to the channel from all active components (m s-1) + basinTotalRunoff, & ! intent(in): total runoff to the channel from all active components (m s-1) bvarData%var(iLookBVAR%routingFractionFuture)%dat, & ! intent(in): fraction of runoff in future time steps (m s-1) bvarData%var(iLookBVAR%routingRunoffFuture)%dat, & ! intent(in): runoff in future time steps (m s-1) ! output diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90 index d66e0df02..06b289adf 100755 --- a/build/source/netcdf/modelwrite.f90 +++ b/build/source/netcdf/modelwrite.f90 @@ -493,6 +493,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file integer(i4b),allocatable :: ncVarID(:) ! netcdf variable id integer(i4b) :: ncSnowID ! index variable id integer(i4b) :: ncSoilID ! index variable id + integer(i4b),dimension(4) :: hdsInitIdx ! intermediate array of loop indices for HDS initial conditions integer(i4b) :: nSoil ! number of soil layers integer(i4b) :: nSnow ! number of snow layers @@ -531,6 +532,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file integer(i4b) :: iHRU ! index of HRUs integer(i4b) :: iGRU ! index of GRUs integer(i4b) :: iVar ! variable index + integer(i4b) :: i ! loop counter logical(lgt) :: okLength ! flag to check if the vector length is OK character(len=256) :: cmessage ! downstream error message ! -------------------------------------------------------------------------------------------------------- @@ -540,7 +542,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! size of prognostic variable vector nProgVars = size(prog_meta) - allocate(ncVarID(nProgVars+1)) ! include 1 additional basin variable in ID array (possibly more later) + allocate(ncVarID(nProgVars+5)) ! include 5 additional basin variable in ID array (routing + HDS variables) ! maximum number of soil layers maxSoil = gru_struc(1)%hruInfo(1)%nSoil @@ -603,6 +605,16 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+1)) err = nf90_put_att(ncid,ncVarID(nProgVars+1),'long_name',trim(bvar_meta(iLookBVAR%routingRunoffFuture)%vardesc)); call netcdf_err(err,message) err = nf90_put_att(ncid,ncVarID(nProgVars+1),'units' ,trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varunit)); call netcdf_err(err,message) + + ! write HDS related states + hdsInitIdx = (/iLookBVAR%pondVolFrac, iLookBVAR%vMin, iLookBVAR%depConAreaFrac, iLookBVAR%pondArea/) ! HDS initial conditions + + do i = 1,size(hdsInitIdx) + iVar = hdsInitIdx(i) + err = nf90_def_var(ncid, trim(bvar_meta(iVar)%varName), nf90_double, (/gruDimID/), ncVarID(nProgVars+i+1)) + err = nf90_put_att(ncid,ncVarID(nProgVars+i+1),'long_name',trim(bvar_meta(iVar)%vardesc)); call netcdf_err(err,message) + err = nf90_put_att(ncid,ncVarID(nProgVars+i+1),'units' ,trim(bvar_meta(iVar)%varunit)); call netcdf_err(err,message) + end do ! writing loop for HDS states ! define index variables - snow err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),nf90_int,(/hruDimID/),ncSnowID); call netcdf_err(err,message) @@ -681,6 +693,12 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! write selected basin variables err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingRunoffFuture)%dat/), start=(/iGRU/),count=(/1,nTimeDelay/)) + ! write pothole storage variables -> HDS + do i = 1,size(hdsInitIdx) + iVar = hdsInitIdx(i) + err=nf90_put_var(ncid,ncVarID(nProgVars+i+1),(/bvar_data%gru(iGRU)%var(iVar)%dat/), start=(/iGRU/),count=(/1/)) + end do ! writing loop for HDS states + end do ! iGRU loop ! close file diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index c66b727d5..531d4216d 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -158,7 +158,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE globalData,only:prog_meta ! metadata for prognostic variables USE globalData,only:bvar_meta ! metadata for basin (GRU) variables USE globalData,only:gru_struc ! gru-hru mapping structures - USE globalData,only:startGRU ! index of first gru for parallel runs + USE globalData,only:startGRU ! index of first gru for parallel runs USE globaldata,only:iname_soil,iname_snow ! named variables to describe the type of layer USE netcdf_util_module,only:nc_file_open ! open netcdf file USE netcdf_util_module,only:nc_file_close ! close netcdf file @@ -170,7 +170,12 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE data_types,only:var_info ! metadata USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages USE updatState_module,only:updateSoil ! update soil states - + ! provide access to model decisions + USE globalData,only:model_decisions ! model decision structure + USE var_lookup,only:iLookDECISIONS ! look-up values for model decisions + ! provide access to the named variables that describe model decisions + USE mDecisions_module,only:HDSmodel ! Hysteretic Depressional Storage model implementation for prairie potholes + implicit none ! -------------------------------------------------------------------------------------------------------- @@ -191,6 +196,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of integer(i4b) :: fileGRU ! number of GRUs in file integer(i4b) :: iVar, i ! loop indices integer(i4b),dimension(1) :: ndx ! intermediate array of loop indices + integer(i4b),dimension(4) :: hdsInitIdx ! intermediate array of loop indices for HDS initial conditions integer(i4b) :: iGRU ! loop index integer(i4b) :: iHRU ! loop index integer(i4b) :: dimID ! varible dimension ids @@ -391,7 +397,54 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! -------------------------------------------------------------------------------------------------------- ! (2) now get the basin variable(s) ! -------------------------------------------------------------------------------------------------------- + + !******************************** + ! get basin variables for HDS + !******************************** + if(model_decisions(iLookDECISIONS%prPotholes)%iDecision == HDSmodel)then + hdsInitIdx = (/iLookBVAR%pondVolFrac, iLookBVAR%vMin, iLookBVAR%depConAreaFrac, iLookBVAR%pondArea/) ! HDS initial conditions + message='updateHDS_States/' + ! get number of GRUs in file + err = nf90_inq_dimid(ncID,"gru",dimID) + if(err/=nf90_noerr)then; message=trim(message)//'must define at least the pondVolFrac initial condition for HDS (per gru)/'//trim(nf90_strerror(err)); return; end if + + if(err==nf90_noerr)then ! proceed if gru dimension exists + err = nf90_inquire_dimension(ncID,dimID,len=fileGRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err)); return; end if + + do i = 1,size(hdsInitIdx) + iVar = hdsInitIdx(i) + + ! get gru-based variable id + err = nf90_inq_varid(ncID,trim(bvar_meta(iVar)%varName),ncVarID) + ! skip any other vairable if missing + if(err/=0)then + write(*,'(A)') ' WARNING: '//trim(bvar_meta(iVar)%varName)//' is not in the initial conditions file ... populating using pondVolFrac' + cycle ! cycle is used to skip to the next hdsInitIdx if not found + endif + + ! initialize the gru variable data + allocate(varData(fileGRU,1),stat=err) + if(err/=0)then; print*, 'err= ',err; message=trim(message)//'problem allocating GRU variable data'; return; endif + + ! get data + err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif + ! store data in basin var (bvar) structure + do iGRU = 1,nGRU + ! put the data into data structures + bvarData%gru(iGRU)%var(iVar)%dat = varData((iGRU+startGRU-1),1) + end do ! end iGRU loop + + ! deallocate temporary data array for next variable + deallocate(varData, stat=err) + if(err/=0)then; message=trim(message)//'problem deallocating GRU variable data'; return; endif + end do ! loop for hdsInitIdx + end if ! end of check that gru dimension exists +endif ! end for model decision HDS + !******************************** + !******************************** + ! get the index in the file: single HRU if(restartFileType/=singleHRU)then