diff --git a/build/FUSE_SRC/driver/fuse_evaluate.f90 b/build/FUSE_SRC/driver/fuse_evaluate.f90 index 0638a5e..cd04f94 100644 --- a/build/FUSE_SRC/driver/fuse_evaluate.f90 +++ b/build/FUSE_SRC/driver/fuse_evaluate.f90 @@ -21,6 +21,7 @@ SUBROUTINE fuse_evaluate(XPAR, info, work, domain, OUTPUT_FLAG, METRIC_VAL) ! Modified by Cyril Thébault to allow different metrics as objective function, 2024 ! Modified by Martyn Clark to call differentiable modeling routines, 12/2025 ! Modified by Martyn Clark to simplify/refactor, 02/2026 + ! Modified by Cyril Thébault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -374,7 +375,7 @@ subroutine advance_one_cell(work, sub_idx, iSpat1, iSpat2, dt_sub, dt_full, err, use multistate, only: gState_3d, FSTATE, MSTATE use multiroute, only: MROUTE, AROUTE_3d use multibands - use multi_flux, only: W_FLUX, W_FLUX_3d + use multi_flux, only: W_FLUX, W_FLUX_3d, M_FLUX use set_all_module, only: SET_STATE, SET_FLUXES, SET_ROUTE ! state vector conversions @@ -485,6 +486,34 @@ subroutine advance_one_cell(work, sub_idx, iSpat1, iSpat2, dt_sub, dt_full, err, end select + ! ------------------------- + ! interception + ! ------------------------- + select case(diff_mode) + + case(original) + M_FLUX%PIN0 = M_FLUX%EFF_PPT + call UPDATE_INTERCEPTION(DELTIM, ierr, cmessage) + + if (ierr /= 0) then + err = 1 + message = trim(cmessage) + return + end if + + case(differentiable) + if (SMODL%iINTRC /= iopt_no_intrcep) then + err = 1 + message = 'advance_one_cell: interception not yet implemented for differentiable mode' + return + end if + + case default + err = 1 + message = 'advance_one_cell: cannot identify diff_mode (interception)' + return + + end select ! ------------------------- ! soil physics ! ------------------------- diff --git a/build/FUSE_SRC/physics_orig/evap_lower.f90 b/build/FUSE_SRC/physics_orig/evap_lower.f90 index 2b1e1c4..91288dd 100644 --- a/build/FUSE_SRC/physics_orig/evap_lower.f90 +++ b/build/FUSE_SRC/physics_orig/evap_lower.f90 @@ -20,6 +20,11 @@ SUBROUTINE EVAP_LOWER() USE multistate ! model states USE multi_flux ! model fluxes IMPLICIT NONE +REAL(WP) :: PET_SOIL + +! Potential evaporation remaining after interception evaporation +PET_SOIL = MAX(0._wp, MFORCE%PET - M_FLUX%EVAP_0) + ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH2) ! lower layer architecture CASE(iopt_tens2pll_2,iopt_fixedsiz_2) @@ -32,9 +37,9 @@ SUBROUTINE EVAP_LOWER() ! ----------------------------------------------------- SELECT CASE(SMODL%iESOIL) CASE(iopt_sequential) - M_FLUX%EVAP_2 = (MFORCE%PET-M_FLUX%EVAP_1) * (TSTATE%TENS_2/DPARAM%MAXTENS_2) + M_FLUX%EVAP_2 = MAX(0._wp, PET_SOIL-M_FLUX%EVAP_1) * (TSTATE%TENS_2/DPARAM%MAXTENS_2) CASE(iopt_rootweight) - M_FLUX%EVAP_2 = MFORCE%PET * DPARAM%RTFRAC2 * (TSTATE%TENS_2/DPARAM%MAXTENS_2) + M_FLUX%EVAP_2 = PET_SOIL * DPARAM%RTFRAC2 * (TSTATE%TENS_2/DPARAM%MAXTENS_2) CASE DEFAULT print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" END SELECT ! (evaporation schemes) diff --git a/build/FUSE_SRC/physics_orig/evap_upper.f90 b/build/FUSE_SRC/physics_orig/evap_upper.f90 index d980437..ef5ba75 100644 --- a/build/FUSE_SRC/physics_orig/evap_upper.f90 +++ b/build/FUSE_SRC/physics_orig/evap_upper.f90 @@ -20,6 +20,11 @@ SUBROUTINE EVAP_UPPER() USE multistate ! model states USE multi_flux ! model fluxes IMPLICIT NONE +REAL(WP) :: PET_SOIL + +! Potential evaporation remaining after interception evaporation +PET_SOIL = MAX(0._wp, MFORCE%PET - M_FLUX%EVAP_0) + ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH1) ! upper layer architecture ! -------------------------------------------------------------------------------------- @@ -29,12 +34,12 @@ SUBROUTINE EVAP_UPPER() ! ----------------------------------------------------- SELECT CASE(SMODL%iESOIL) CASE(iopt_sequential) - M_FLUX%EVAP_1A = MFORCE%PET * TSTATE%TENS_1A/DPARAM%MAXTENS_1A - M_FLUX%EVAP_1B = (MFORCE%PET - M_FLUX%EVAP_1A) * TSTATE%TENS_1B/DPARAM%MAXTENS_1B + M_FLUX%EVAP_1A = PET_SOIL * TSTATE%TENS_1A/DPARAM%MAXTENS_1A + M_FLUX%EVAP_1B = MAX(0._wp, PET_SOIL - M_FLUX%EVAP_1A) * TSTATE%TENS_1B/DPARAM%MAXTENS_1B M_FLUX%EVAP_1 = M_FLUX%EVAP_1A + M_FLUX%EVAP_1B CASE(iopt_rootweight) - M_FLUX%EVAP_1A = MFORCE%PET * MPARAM%RTFRAC1 * TSTATE%TENS_1A/DPARAM%MAXTENS_1A - M_FLUX%EVAP_1B = MFORCE%PET * DPARAM%RTFRAC2 * TSTATE%TENS_1B/DPARAM%MAXTENS_1B + M_FLUX%EVAP_1A = PET_SOIL * MPARAM%RTFRAC1 * TSTATE%TENS_1A/DPARAM%MAXTENS_1A + M_FLUX%EVAP_1B = PET_SOIL * DPARAM%RTFRAC2 * TSTATE%TENS_1B/DPARAM%MAXTENS_1B M_FLUX%EVAP_1 = M_FLUX%EVAP_1A + M_FLUX%EVAP_1B CASE DEFAULT print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" @@ -49,11 +54,11 @@ SUBROUTINE EVAP_UPPER() CASE(iopt_sequential) M_FLUX%EVAP_1A = 0._wp M_FLUX%EVAP_1B = 0._wp - M_FLUX%EVAP_1 = MFORCE%PET * TSTATE%TENS_1/DPARAM%MAXTENS_1 + M_FLUX%EVAP_1 = PET_SOIL * TSTATE%TENS_1/DPARAM%MAXTENS_1 CASE(iopt_rootweight) M_FLUX%EVAP_1A = 0._wp M_FLUX%EVAP_1B = 0._wp - M_FLUX%EVAP_1 = MFORCE%PET * MPARAM%RTFRAC1 * TSTATE%TENS_1/DPARAM%MAXTENS_1 + M_FLUX%EVAP_1 = PET_SOIL * MPARAM%RTFRAC1 * TSTATE%TENS_1/DPARAM%MAXTENS_1 CASE DEFAULT print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" END SELECT ! (evaporation schemes) diff --git a/build/FUSE_SRC/physics_orig/update_interception.f90 b/build/FUSE_SRC/physics_orig/update_interception.f90 new file mode 100644 index 0000000..657d62e --- /dev/null +++ b/build/FUSE_SRC/physics_orig/update_interception.f90 @@ -0,0 +1,222 @@ +SUBROUTINE UPDATE_INTERCEPTION(DT, IERR, MESSAGE) +! --------------------------------------------------------------------------------------- +! Creator: +! -------- +! Cyril Thebault, 2026 +! --------------------------------------------------------------------------------------- +! Purpose: +! -------- +! Updates the interception store. +! +! PIN0 : precipitation entering the interception store [depth/time] +! EVAP_0 : evaporation from the interception store [depth/time] +! PTHRU : throughfall leaving the interception store [depth/time] +! SINT_0 : interception storage [depth] +! REFSINT_0 : characteristic interception storage [depth] +! --------------------------------------------------------------------------------------- + +USE nrtype +USE model_defn +USE model_defnames +USE multiparam +USE multiforce, ONLY: MFORCE +USE multistate +USE multi_flux +USE model_numerix, ONLY: ERR_ITER_FUNC, ERR_ITER_DX, NITER_TOTAL + +IMPLICIT NONE + +! Input +REAL(WP), INTENT(IN) :: DT + +! Output +INTEGER(I4B), INTENT(OUT) :: IERR +CHARACTER(*), INTENT(OUT) :: MESSAGE + +! Local variables +INTEGER(I4B) :: ITER + +REAL(WP) :: S_OLD +REAL(WP) :: S_NEW +REAL(WP) :: S_TRIAL +REAL(WP) :: S_LO +REAL(WP) :: S_HI + +REAL(WP) :: PRECIP +REAL(WP) :: PET + +REAL(WP), PARAMETER :: K_SMOOTH = 1.e-5_wp ! It may need to be flexible based on REFSINT_0; e.g. K_SMOOTH=MAX(1.e-6_wp, 0.01_wp * MPARAM%REFSINT_0) +REAL(WP), PARAMETER :: EPS_WET = 1.e-6_wp + +REAL(WP) :: X_LOGISTIC +REAL(WP) :: EXP_X +REAL(WP) :: PHI +REAL(WP) :: WET_FRAC + +REAL(WP) :: DPHI_DS +REAL(WP) :: DWET_DS + +REAL(WP) :: RESIDUAL +REAL(WP) :: DRESIDUAL +REAL(WP) :: NEWTON_STEP + +LOGICAL(LGT) :: CONVERGED + +! --------------------------------------------------------------------------------------- +! Initialise error handling +! --------------------------------------------------------------------------------------- +IERR = 0 +MESSAGE = 'UPDATE_INTERCEPTION/' +CONVERGED = .FALSE. + +IF (DT <= 0._wp) THEN + IERR = 1 + MESSAGE = TRIM(MESSAGE)//'DT must be greater than zero' + RETURN +END IF + +! --------------------------------------------------------------------------------------- +! Interception options +! --------------------------------------------------------------------------------------- +SELECT CASE(SMODL%iINTRC) + + CASE(iopt_no_intrcep) + + FSTATE%SINT_0 = 0._wp + + M_FLUX%EVAP_0 = 0._wp + M_FLUX%PTHRU = M_FLUX%PIN0 + M_FLUX%EFF_PPT = M_FLUX%PTHRU + + RETURN + + CASE(iopt_gr5h_intrc) + + IF (MPARAM%REFSINT_0 < 0._wp) THEN + IERR = 1 + MESSAGE = TRIM(MESSAGE)//'REFSINT_0 must not be negative' + RETURN + END IF + + CASE DEFAULT + + IERR = 1 + MESSAGE = TRIM(MESSAGE)//'SMODL%iINTRC must be iopt_no_intrcep or iopt_gr5h_intrc' + RETURN + +END SELECT + +! --------------------------------------------------------------------------------------- +! Define forcing and smoothing scales +! --------------------------------------------------------------------------------------- + +PRECIP = MAX(0._wp, M_FLUX%PIN0) +PET = MAX(0._wp, MFORCE%PET) + +S_OLD = MAX(0._wp, FSTATE%SINT_0) + +! The storage cannot become negative and cannot exceed the initial storage +! plus all incoming precipitation when evaporation and throughfall are omitted. +S_LO = 0._wp +S_HI = S_OLD + PRECIP * DT + +! Handle the trivial dry case explicitly. +IF (S_HI <= ERR_ITER_DX) THEN + S_NEW = 0._wp + CONVERGED = .TRUE. +ELSE + ! Start Newton from the old storage, constrained to the bracket. + S_NEW = MIN(MAX(S_OLD, S_LO), S_HI) +END IF + +! --------------------------------------------------------------------------------------- +! Solve the backward-Euler residual +! +! R(S) = S - S_OLD - DT * [PRECIP*(1-PHI(S)) - PET*WET_FRAC(S)] +! --------------------------------------------------------------------------------------- +IF (.NOT.CONVERGED) THEN + + DO ITER = 1, NITER_TOTAL + + ! Stable evaluation of the logistic function: + ! PHI = 1 / (1 + exp(-(S-REFSINT_0)/K_SMOOTH)) + X_LOGISTIC = (S_NEW - MPARAM%REFSINT_0) / K_SMOOTH + + IF (X_LOGISTIC >= 0._wp) THEN + PHI = 1._wp / (1._wp + EXP(-X_LOGISTIC)) + ELSE + EXP_X = EXP(X_LOGISTIC) + PHI = EXP_X / (1._wp + EXP_X) + END IF + + WET_FRAC = S_NEW / (S_NEW + EPS_WET) + + DPHI_DS = PHI * (1._wp - PHI) / K_SMOOTH + DWET_DS = EPS_WET / (S_NEW + EPS_WET)**2 + + RESIDUAL = S_NEW - S_OLD - DT * (PRECIP * (1._wp - PHI) - PET * WET_FRAC) + + DRESIDUAL = 1._wp + DT * (PRECIP * DPHI_DS + PET * DWET_DS) + + IF (ABS(RESIDUAL) <= ERR_ITER_FUNC) THEN + CONVERGED = .TRUE. + EXIT + END IF + + ! The residual is monotonic. Update the root bracket first. + IF (RESIDUAL < 0._wp) THEN + S_LO = S_NEW + ELSE + S_HI = S_NEW + END IF + + ! Newton candidate. + NEWTON_STEP = -RESIDUAL / DRESIDUAL + S_TRIAL = S_NEW + NEWTON_STEP + + ! Safeguard Newton with bisection. + IF (S_TRIAL <= S_LO .OR. S_TRIAL >= S_HI) THEN + S_TRIAL = 0.5_wp * (S_LO + S_HI) + END IF + + IF (ABS(S_TRIAL - S_NEW) <= ERR_ITER_DX) THEN + S_NEW = S_TRIAL + CONVERGED = .TRUE. + EXIT + END IF + + S_NEW = S_TRIAL + + END DO + +END IF + +IF (.NOT.CONVERGED) THEN + IERR = 1 + MESSAGE = TRIM(MESSAGE)// & + 'Newton-bisection failed to converge' + RETURN +END IF + +! --------------------------------------------------------------------------------------- +! Evaluate fluxes at the converged end-of-step storage +! --------------------------------------------------------------------------------------- +X_LOGISTIC = (S_NEW - MPARAM%REFSINT_0) / K_SMOOTH + +IF (X_LOGISTIC >= 0._wp) THEN + PHI = 1._wp / (1._wp + EXP(-X_LOGISTIC)) +ELSE + EXP_X = EXP(X_LOGISTIC) + PHI = EXP_X / (1._wp + EXP_X) +END IF + +WET_FRAC = S_NEW / (S_NEW + EPS_WET) + +FSTATE%SINT_0 = MAX(0._wp, S_NEW) + +M_FLUX%PTHRU = PRECIP * PHI +M_FLUX%EVAP_0 = PET * WET_FRAC +M_FLUX%EFF_PPT = M_FLUX%PTHRU + + +END SUBROUTINE UPDATE_INTERCEPTION \ No newline at end of file diff --git a/build/FUSE_SRC/prelim/assign_flx.f90 b/build/FUSE_SRC/prelim/assign_flx.f90 index 5f75894..dba413f 100644 --- a/build/FUSE_SRC/prelim/assign_flx.f90 +++ b/build/FUSE_SRC/prelim/assign_flx.f90 @@ -31,11 +31,29 @@ SUBROUTINE ASSIGN_FLX() N_FLUX=0 C_FLUX(:)%FNAME = ' ' ! --------------------------------------------------------------------------------------- + ! (0) INTERCEPTION + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iINTRC) + + CASE(iopt_no_intrcep) + N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EFF_PPT ' + + CASE(iopt_gr5h_intrc) + N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'PIN0 ' + N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'PTHRU ' + N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EVAP_0 ' + N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EFF_PPT ' + + CASE DEFAULT + print *, "SMODL%iINTRC must be iopt_no_intrcep or iopt_gr5h_intrc" + STOP + + END SELECT + ! --------------------------------------------------------------------------------------- ! (1) DEFINE STATE VARIABLES IN THE UPPER LAYER ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH1) CASE(iopt_tension2_1) - N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EFF_PPT ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EVAP_1A ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EVAP_1B ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'RCHR2EXCS ' @@ -45,7 +63,6 @@ SUBROUTINE ASSIGN_FLX() N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'OFLOW_1 ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'QSURF ' CASE(iopt_tension1_1) - N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EFF_PPT ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EVAP_1 ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'TENS2FREE_1' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'QPERC_12 ' @@ -53,7 +70,6 @@ SUBROUTINE ASSIGN_FLX() N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'OFLOW_1 ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'QSURF ' CASE(iopt_onestate_1) - N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EFF_PPT ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'EVAP_1 ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'QPERC_12 ' N_FLUX=N_FLUX+1; C_FLUX(N_FLUX)%FNAME = 'QINTF_1 ' diff --git a/build/FUSE_SRC/prelim/assign_par.f90 b/build/FUSE_SRC/prelim/assign_par.f90 index 07b880a..99c3a01 100644 --- a/build/FUSE_SRC/prelim/assign_par.f90 +++ b/build/FUSE_SRC/prelim/assign_par.f90 @@ -11,6 +11,7 @@ SUBROUTINE ASSIGN_PAR() ! -------- ! Martyn Clark, 2007 ! Modified by Brian Henn to include snow model, 6/2013 + ! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -80,7 +81,19 @@ SUBROUTINE ASSIGN_PAR() STOP END SELECT ! --------------------------------------------------------------------------------------- - ! (3) UPPER-LAYER ARCHITECTURE + ! (3) INTERCEPTION + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iINTRC) + CASE(iopt_no_intrcep) ! no interception store + CASE(iopt_gr5h_intrc) ! grh5 interception store + MPAR=MPAR+1; LPARAM(MPAR)%PARNAME = 'REFSINT_0' ! characteristic interception storage (50% throughfall) (mm) + CASE DEFAULT + print *, "SMODL%iINTRC must be either iopt_no_intrcep or iopt_gr5h_intrc" + STOP + + END SELECT + ! --------------------------------------------------------------------------------------- + ! (4) UPPER-LAYER ARCHITECTURE ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH1) CASE(iopt_tension2_1) ! tension storage sub-divided into recharge and excess @@ -96,7 +109,7 @@ SUBROUTINE ASSIGN_PAR() STOP END SELECT ! (different upper-layer architechure) ! --------------------------------------------------------------------------------------- - ! (4) LOWER-LAYER ARCHITECTURE / BASEFLOW + ! (5) LOWER-LAYER ARCHITECTURE / BASEFLOW ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH2) CASE(iopt_tens2pll_2) ! tension reservoir plus two parallel tanks @@ -127,7 +140,7 @@ SUBROUTINE ASSIGN_PAR() STOP END SELECT ! different lower-layer architecture / baseflow parameterizations) ! --------------------------------------------------------------------------------------- - ! (5) EVAPORATION + ! (6) EVAPORATION ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iESOIL) CASE(iopt_sequential) @@ -138,7 +151,7 @@ SUBROUTINE ASSIGN_PAR() print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight'" END SELECT ! (different evaporation schemes) ! --------------------------------------------------------------------------------------- - ! (6) PERCOLATION + ! (7) PERCOLATION ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iQPERC) CASE(iopt_perc_f2sat,iopt_perc_w2sat) ! standard equation k(theta)**c @@ -152,7 +165,7 @@ SUBROUTINE ASSIGN_PAR() STOP END SELECT ! (different percolation options) ! --------------------------------------------------------------------------------------- - ! (7) INTERFLOW + ! (8) INTERFLOW ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iQINTF) CASE(iopt_intflwsome) ! interflow @@ -164,7 +177,7 @@ SUBROUTINE ASSIGN_PAR() STOP END SELECT ! (different interflow options) ! --------------------------------------------------------------------------------------- - ! (8) SURFACE RUNOFF + ! (9) SURFACE RUNOFF ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iQSURF) CASE(iopt_arno_x_vic) ! ARNO/Xzang/VIC parameterization (upper zone control) @@ -188,7 +201,7 @@ SUBROUTINE ASSIGN_PAR() STOP END SELECT ! (different surface runoff options) ! --------------------------------------------------------------------------------------- - ! (9) TIME DELAY IN RUNOFF + ! (10) TIME DELAY IN RUNOFF ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iQ_TDH) CASE(iopt_rout_gamma) ! use a Gamma distribution with shape parameter = 2.5 diff --git a/build/FUSE_SRC/prelim/init_state.f90 b/build/FUSE_SRC/prelim/init_state.f90 index 06dfb28..5f7cee8 100644 --- a/build/FUSE_SRC/prelim/init_state.f90 +++ b/build/FUSE_SRC/prelim/init_state.f90 @@ -4,6 +4,7 @@ SUBROUTINE INIT_STATE(FRAC) ! -------- ! Martyn Clark, 2007 ! Modified by Brian Henn to include snow model, 6/2013 +! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -22,6 +23,8 @@ SUBROUTINE INIT_STATE(FRAC) REAL(WP), INTENT(IN) :: FRAC ! fraction of capacity INTEGER(I4B) :: ISNW ! snow band index ! --------------------------------------------------------------------------------------- +! interception layer +FSTATE%SINT_0 = 0._wp ! (upper layer) FSTATE%TENS_1A = DPARAM%MAXTENS_1A * FRAC FSTATE%TENS_1B = DPARAM%MAXTENS_1B * FRAC diff --git a/build/FUSE_SRC/prelim/uniquemodl.f90 b/build/FUSE_SRC/prelim/uniquemodl.f90 index a2ea619..6ce74f0 100644 --- a/build/FUSE_SRC/prelim/uniquemodl.f90 +++ b/build/FUSE_SRC/prelim/uniquemodl.f90 @@ -40,6 +40,7 @@ SUBROUTINE UNIQUEMODL(NMOD) INTEGER(I4B) :: ISW_QINTF ! loop thru interflow INTEGER(I4B) :: ISW_Q_TDH ! loop thru time delay options INTEGER(I4B) :: ISW_SNOWM ! loop thru snow model options + INTEGER(I4B) :: ISW_INTRC ! loop thru interception options ! Start procedure here !err=0; message="UNIQUEMODL/ok" ! --------------------------------------------------------------------------------------- @@ -77,6 +78,9 @@ SUBROUTINE UNIQUEMODL(NMOD) ! snow model switch LIST_SNOWM(1)%MCOMPONENT = 'no_snowmod' ! no snow model LIST_SNOWM(2)%MCOMPONENT = 'temp_index' ! temperature index snow model + ! interception + LIST_INTRC(1)%MCOMPONENT = 'no_intrcep' ! no interception store + LIST_INTRC(2)%MCOMPONENT = 'gr5h_intrc' ! gr5h interception store ! --------------------------------------------------------------------------------------- ! (2) LOOP THROUGH MODEL COMPONENTS AND DEFINE A SET OF UNIQUE MODELS ! --------------------------------------------------------------------------------------- @@ -90,57 +94,61 @@ SUBROUTINE UNIQUEMODL(NMOD) ! g) define interflow method ! h) define time delay in runoff ICOUNT = 0 ! initialize counter - ! loop through snow model options - DO ISW_SNOWM=1,SIZE(LIST_SNOWM) - ! (loop through time delay options) - DO ISW_Q_TDH=1,SIZE(LIST_Q_TDH) - ! (loop through interflow options) - DO ISW_QINTF=1,SIZE(LIST_QINTF) - ! (loop through evaporation options) - DO ISW_ESOIL=1,SIZE(LIST_ESOIL) - ! (loop through percolation options) - DO ISW_QPERC=1,SIZE(LIST_QPERC) - ! (loop through surface runoff options) - DO ISW_QSURF=1,SIZE(LIST_QSURF) - ! (loop through lower-layer architecture options) - DO ISW_ARCH2=1,SIZE(LIST_ARCH2) - ! (loop through upper-layer architecture options) - DO ISW_ARCH1=1,SIZE(LIST_ARCH1) - ! (loop through rainfall error options) - DO ISW_RFERR=1,SIZE(LIST_RFERR) - ! don't allow a lower tension tank when there are two upper ones - IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT(1:10).EQ.'tension2_1'.AND. & - LIST_ARCH2(ISW_ARCH2)%MCOMPONENT(1:10).EQ.'tens2pll_2') CYCLE - ! don't allow percolation below field capacity if there are multiple upper tanks - IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT(1:10).NE.'onestate_1'.AND. & - LIST_QPERC(ISW_QPERC)%MCOMPONENT(1:10).EQ.'perc_w2sat') CYCLE - ICOUNT = ICOUNT + 1 ! (increment counter) - IF (ICOUNT.LE.SIZE(AMODL)) THEN - ! save unique model combinations - AMODL(ICOUNT)%iRFERR = desc_str2int(LIST_RFERR(ISW_RFERR)%MCOMPONENT) - AMODL(ICOUNT)%iARCH1 = desc_str2int(LIST_ARCH1(ISW_ARCH1)%MCOMPONENT) - AMODL(ICOUNT)%iARCH2 = desc_str2int(LIST_ARCH2(ISW_ARCH2)%MCOMPONENT) - AMODL(ICOUNT)%iQSURF = desc_str2int(LIST_QSURF(ISW_QSURF)%MCOMPONENT) - AMODL(ICOUNT)%iQPERC = desc_str2int(LIST_QPERC(ISW_QPERC)%MCOMPONENT) - AMODL(ICOUNT)%iESOIL = desc_str2int(LIST_ESOIL(ISW_ESOIL)%MCOMPONENT) - AMODL(ICOUNT)%iQINTF = desc_str2int(LIST_QINTF(ISW_QINTF)%MCOMPONENT) - AMODL(ICOUNT)%iQ_TDH = desc_str2int(LIST_Q_TDH(ISW_Q_TDH)%MCOMPONENT) - AMODL(ICOUNT)%iSNOWM = desc_str2int(LIST_Q_TDH(ISW_SNOWM)%MCOMPONENT) - !write(*,'(i3,1x,7(a10,1x))') icount, amodl(icount) - ELSE - ! need to allocate more space - print *, 'insufficent space to hold model combinations' - stop - ENDIF - END DO ! RFERR - END DO ! ARCH1 - END DO ! ARCH2 - END DO ! QSURF - END DO ! QPERC - END DO ! ESOIL - END DO ! QINTF - END DO ! Q_TDH - END DO ! SNOWM + ! loop through interception options + DO ISW_INTRC = 1, SIZE(LIST_INTRC) + ! loop through snow model options + DO ISW_SNOWM=1,SIZE(LIST_SNOWM) + ! (loop through time delay options) + DO ISW_Q_TDH=1,SIZE(LIST_Q_TDH) + ! (loop through interflow options) + DO ISW_QINTF=1,SIZE(LIST_QINTF) + ! (loop through evaporation options) + DO ISW_ESOIL=1,SIZE(LIST_ESOIL) + ! (loop through percolation options) + DO ISW_QPERC=1,SIZE(LIST_QPERC) + ! (loop through surface runoff options) + DO ISW_QSURF=1,SIZE(LIST_QSURF) + ! (loop through lower-layer architecture options) + DO ISW_ARCH2=1,SIZE(LIST_ARCH2) + ! (loop through upper-layer architecture options) + DO ISW_ARCH1=1,SIZE(LIST_ARCH1) + ! (loop through rainfall error options) + DO ISW_RFERR=1,SIZE(LIST_RFERR) + ! don't allow a lower tension tank when there are two upper ones + IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT(1:10).EQ.'tension2_1'.AND. & + LIST_ARCH2(ISW_ARCH2)%MCOMPONENT(1:10).EQ.'tens2pll_2') CYCLE + ! don't allow percolation below field capacity if there are multiple upper tanks + IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT(1:10).NE.'onestate_1'.AND. & + LIST_QPERC(ISW_QPERC)%MCOMPONENT(1:10).EQ.'perc_w2sat') CYCLE + ICOUNT = ICOUNT + 1 ! (increment counter) + IF (ICOUNT.LE.SIZE(AMODL)) THEN + ! save unique model combinations + AMODL(ICOUNT)%iRFERR = desc_str2int(LIST_RFERR(ISW_RFERR)%MCOMPONENT) + AMODL(ICOUNT)%iARCH1 = desc_str2int(LIST_ARCH1(ISW_ARCH1)%MCOMPONENT) + AMODL(ICOUNT)%iARCH2 = desc_str2int(LIST_ARCH2(ISW_ARCH2)%MCOMPONENT) + AMODL(ICOUNT)%iQSURF = desc_str2int(LIST_QSURF(ISW_QSURF)%MCOMPONENT) + AMODL(ICOUNT)%iQPERC = desc_str2int(LIST_QPERC(ISW_QPERC)%MCOMPONENT) + AMODL(ICOUNT)%iESOIL = desc_str2int(LIST_ESOIL(ISW_ESOIL)%MCOMPONENT) + AMODL(ICOUNT)%iQINTF = desc_str2int(LIST_QINTF(ISW_QINTF)%MCOMPONENT) + AMODL(ICOUNT)%iQ_TDH = desc_str2int(LIST_Q_TDH(ISW_Q_TDH)%MCOMPONENT) + AMODL(ICOUNT)%iSNOWM = desc_str2int(LIST_SNOWM(ISW_SNOWM)%MCOMPONENT) + AMODL(ICOUNT)%iINTRC = desc_str2int(LIST_INTRC(ISW_INTRC)%MCOMPONENT) + !write(*,'(i3,1x,7(a10,1x))') icount, amodl(icount) + ELSE + ! need to allocate more space + print *, 'insufficent space to hold model combinations' + stop + ENDIF + END DO ! RFERR + END DO ! ARCH1 + END DO ! ARCH2 + END DO ! QSURF + END DO ! QPERC + END DO ! ESOIL + END DO ! QINTF + END DO ! Q_TDH + END DO ! SNOWM + END DO ! INTRC ! --------------------------------------------------------------------------------------- NMOD = ICOUNT !pause diff --git a/build/FUSE_SRC/runtime/initfluxes.f90 b/build/FUSE_SRC/runtime/initfluxes.f90 index 1cd3378..62e02aa 100644 --- a/build/FUSE_SRC/runtime/initfluxes.f90 +++ b/build/FUSE_SRC/runtime/initfluxes.f90 @@ -4,6 +4,7 @@ SUBROUTINE INITFLUXES() ! -------- ! Martyn Clark, 2007 ! Modified by Brian Henn to include snow model, 6/2013 +! Modified by Cyril Thébault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -21,6 +22,9 @@ SUBROUTINE INITFLUXES() IMPLICIT NONE INTEGER(I4B) :: ISNW ! index for looping though SWE ! --------------------------------------------------------------------------------------- +M_FLUX%EVAP_0 = 0._wp; W_FLUX%EVAP_0 = 0._wp +M_FLUX%PIN0 = 0._wp; W_FLUX%PIN0 = 0._wp +M_FLUX%PTHRU = 0._wp; W_FLUX%PTHRU = 0._wp M_FLUX%EFF_PPT = 0._wp; W_FLUX%EFF_PPT = 0._wp M_FLUX%SATAREA = 0._wp; W_FLUX%SATAREA = 0._wp M_FLUX%QSURF = 0._wp; W_FLUX%QSURF = 0._wp diff --git a/build/FUSE_SRC/runtime/set_all.f90 b/build/FUSE_SRC/runtime/set_all.f90 index 8e59064..32c0947 100644 --- a/build/FUSE_SRC/runtime/set_all.f90 +++ b/build/FUSE_SRC/runtime/set_all.f90 @@ -24,6 +24,8 @@ SUBROUTINE SET_STATE(VAL) REAL(WP), INTENT(IN) :: VAL ! value INTEGER(I4B) :: ISNW ! snow band index ! --------------------------------------------------------------------------------------- + ! interception + FSTATE%SINT_0 = VAL ! upper layer FSTATE%TENS_1A = VAL FSTATE%TENS_1B = VAL @@ -67,6 +69,9 @@ SUBROUTINE SET_FLUXES(VAL) REAL(WP), INTENT(IN) :: VAL ! value INTEGER(I4B) :: ISNW ! index for looping though SWE ! --------------------------------------------------------------------------------------- + M_FLUX%EVAP_0 = VAL; W_FLUX%EVAP_0 = VAL + M_FLUX%PIN0 = VAL; W_FLUX%PIN0 = VAL + M_FLUX%PTHRU = VAL; W_FLUX%PTHRU = VAL M_FLUX%EFF_PPT = VAL; W_FLUX%EFF_PPT = VAL M_FLUX%SATAREA = VAL; W_FLUX%SATAREA = VAL M_FLUX%QSURF = VAL; W_FLUX%QSURF = VAL diff --git a/build/FUSE_SRC/share/model_defn_data.f90 b/build/FUSE_SRC/share/model_defn_data.f90 index 80f34a1..8e50bed 100644 --- a/build/FUSE_SRC/share/model_defn_data.f90 +++ b/build/FUSE_SRC/share/model_defn_data.f90 @@ -6,6 +6,7 @@ MODULE model_defn ! Martyn Clark ! Modified by Brian Henn to include snow model, 6/2013 ! Modified by Martyn Clark to separate type definitions from data storage, 01/2026 + ! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- USE nrtype @@ -17,13 +18,13 @@ MODULE model_defn private public :: NDEC, NTDH_MAX, NSTATE, N_FLUX - public :: LIST_RFERR, LIST_ARCH1, LIST_ARCH2, LIST_QSURF, LIST_QPERC, LIST_ESOIL, LIST_QINTF, LIST_Q_TDH, LIST_SNOWM + public :: LIST_RFERR, LIST_ARCH1, LIST_ARCH2, LIST_QSURF, LIST_QPERC, LIST_ESOIL, LIST_QINTF, LIST_Q_TDH, LIST_SNOWM, LIST_INTRC public :: FNAME_PREFIX, FNAME_TEMPRY, FNAME_ASCII public :: FNAME_NETCDF_RUNS, FNAME_NETCDF_PARA, FNAME_NETCDF_PARA_SCE, FNAME_NETCDF_PARA_PRE public :: AMODL, SMODL, CSTATE, C_FLUX ! list of combinations in each model component - INTEGER, PARAMETER :: NDEC = 9 ! number of model decisions + INTEGER, PARAMETER :: NDEC = 10 ! number of model decisions TYPE(DESC), DIMENSION(2) :: LIST_RFERR ! rainfall error TYPE(DESC), DIMENSION(3) :: LIST_ARCH1 ! upper-layer architecture TYPE(DESC), DIMENSION(4) :: LIST_ARCH2 ! lower-layer architecture @@ -33,6 +34,7 @@ MODULE model_defn TYPE(DESC), DIMENSION(2) :: LIST_QINTF ! interflow TYPE(DESC), DIMENSION(2) :: LIST_Q_TDH ! time delay in runoff TYPE(DESC), DIMENSION(2) :: LIST_SNOWM ! snow model + TYPE(DESC), DIMENSION(2) :: LIST_INTRC ! interception ! max steps in routing function INTEGER(I4B),PARAMETER::NTDH_MAX=500 @@ -47,7 +49,7 @@ MODULE model_defn CHARACTER(LEN=256) :: FNAME_ASCII ! ASCII output filename TYPE(UMODEL),DIMENSION(5000) :: AMODL ! (model definition -- all) TYPE(UMODEL) :: SMODL ! (model definition -- single model) - TYPE(SNAMES),DIMENSION(7) :: CSTATE ! (list of model states for SMODL) + TYPE(SNAMES),DIMENSION(8) :: CSTATE ! (list of model states for SMODL) TYPE(FNAMES),DIMENSION(50) :: C_FLUX ! (list of model fluxes for SMODL) INTEGER(I4B) :: NSTATE=0 ! number of model states INTEGER(I4B) :: N_FLUX=0 ! number of model fluxes diff --git a/build/FUSE_SRC/share/model_defnames.f90 b/build/FUSE_SRC/share/model_defnames.f90 index 93322c6..056fd45 100644 --- a/build/FUSE_SRC/share/model_defnames.f90 +++ b/build/FUSE_SRC/share/model_defnames.f90 @@ -1,132 +1,138 @@ -module model_defnames -! Purpose: Contains routines for alternating between char <-> int names -! Programmers: David McInerney and Dmitri Kavetski (University of Adelaide) -USE nrtype -implicit none -! parameterised descriptions -integer(I4B), parameter :: iopt_additive_e = 1001, & - iopt_multiplc_e = 1002, & - iopt_tension1_1 = 2001, & - iopt_tension2_1 = 2002, & - iopt_onestate_1 = 2003, & - iopt_tens2pll_2 = 3001, & - iopt_unlimfrc_2 = 3002, & - iopt_unlimpow_2 = 3003, & - iopt_fixedsiz_2 = 3004, & - iopt_topmdexp_2 = 3005, & - iopt_arno_x_vic = 4001, & - iopt_prms_varnt = 4002, & - iopt_tmdl_param = 4003, & - iopt_perc_f2sat = 5001, & - iopt_perc_w2sat = 5002, & - iopt_perc_lower = 5003, & - iopt_sequential = 6001, & - iopt_rootweight = 6002, & - iopt_intflwnone = 7001, & - iopt_intflwsome = 7002, & - iopt_rout_gamma = 8001, & - iopt_no_routing = 8002, & - iopt_no_snowmod = 8501, & - iopt_temp_index = 8502 -! --- -integer(I4B), parameter :: iopt_TENS1A = 9001, & - iopt_TENS1B = 9002, & - iopt_TENS_1 = 9003, & - iopt_FREE_1 = 9004, & - iopt_WATR_1 = 9005, & - iopt_TENS_2 = 9006, & - iopt_FREE2A = 9007, & - iopt_FREE2B = 9008, & - iopt_WATR_2 = 9009 -! ------------------------------------------ -contains -! ------------------------------------------ -elemental function desc_str2int(name)result(res) -! Purpose: Converts a string description into its corresponding integer value. -implicit none -! dummies -character(*), intent(in) :: name -integer(I4B) :: res -! Start procedure here -selectcase(name) -case("additive_e"); res = iopt_additive_e -case("multiplc_e"); res = iopt_multiplc_e -case("tension1_1"); res = iopt_tension1_1 -case("tension2_1"); res = iopt_tension2_1 -case("onestate_1"); res = iopt_onestate_1 -case("tens2pll_2"); res = iopt_tens2pll_2 -case("unlimfrc_2"); res = iopt_unlimfrc_2 -case("unlimpow_2"); res = iopt_unlimpow_2 -case("fixedsiz_2"); res = iopt_fixedsiz_2 -case("arno_x_vic"); res = iopt_arno_x_vic -case("prms_varnt"); res = iopt_prms_varnt -case("tmdl_param"); res = iopt_tmdl_param -case("perc_f2sat"); res = iopt_perc_f2sat -case("perc_w2sat"); res = iopt_perc_w2sat -case("perc_lower"); res = iopt_perc_lower -case("sequential"); res = iopt_sequential -case("rootweight"); res = iopt_rootweight -case("intflwnone"); res = iopt_intflwnone -case("intflwsome"); res = iopt_intflwsome -case("rout_gamma"); res = iopt_rout_gamma -case("no_routing"); res = iopt_no_routing -case("no_snowmod"); res = iopt_no_snowmod -case("temp_index"); res = iopt_temp_index -case("TENS1B"); res = iopt_TENS1B -case("TENS_1"); res = iopt_TENS_1 -case("FREE_1"); res = iopt_FREE_1 -case("WATR_1"); res = iopt_WATR_1 -case("TENS_2"); res = iopt_TENS_2 -case("FREE2A"); res = iopt_FREE2A -case("FREE2B"); res = iopt_FREE2B -case("WATR_2"); res = iopt_WATR_2 -case default; res = -999 -endselect -! End procedure here -endfunction desc_str2int -! ------------------------------------------ -elemental function desc_int2str(intVal)result(res) -! Purpose: Converts an integer description into corresponding string value -implicit none -! dummies -integer(I4B), intent(in) :: intVal -character(10) :: res -! Start procedure here -selectcase(intVal) -case(iopt_TENS1B); res = "TENS1B" -case(iopt_TENS_1); res = "TENS_1" -case(iopt_FREE_1); res = "FREE_1" -case(iopt_WATR_1); res = "WATR_1" -case(iopt_TENS_2); res = "TENS_2" -case(iopt_FREE2A); res = "FREE2A" -case(iopt_FREE2B); res = "FREE2B" -case(iopt_WATR_2); res = "WATR_2" -case(iopt_additive_e); res = "additive_e" -case(iopt_multiplc_e); res = "multiplc_e" -case(iopt_tension1_1); res = "tension1_1" -case(iopt_tension2_1); res = "tension2_1" -case(iopt_onestate_1); res = "onestate_1" -case(iopt_tens2pll_2); res = "tens2pll_2" -case(iopt_unlimfrc_2); res = "unlimfrc_2" -case(iopt_unlimpow_2); res = "unlimpow_2" -case(iopt_fixedsiz_2); res = "fixedsiz_2" -case(iopt_arno_x_vic); res = "arno_x_vic" -case(iopt_prms_varnt); res = "prms_varnt" -case(iopt_tmdl_param); res = "tmdl_param" -case(iopt_perc_f2sat); res = "perc_f2sat" -case(iopt_perc_w2sat); res = "perc_w2sat" -case(iopt_perc_lower); res = "perc_lower" -case(iopt_sequential); res = "sequential" -case(iopt_rootweight); res = "rootweight" -case(iopt_intflwnone); res = "intflwnone" -case(iopt_intflwsome); res = "intflwsome" -case(iopt_rout_gamma); res = "rout_gamma" -case(iopt_no_routing); res = "no_routing" -case(iopt_no_snowmod); res = "no_snowmod" -case(iopt_temp_index); res = "temp_index" -case default; res = "UNDFND" -endselect -! End procedure here -endfunction desc_int2str -! ------------------------------------------ -endmodule model_defnames +module model_defnames +! Purpose: Contains routines for alternating between char <-> int names +! Programmers: David McInerney and Dmitri Kavetski (University of Adelaide) +USE nrtype +implicit none +! parameterised descriptions +integer(I4B), parameter :: iopt_additive_e = 1001, & + iopt_multiplc_e = 1002, & + iopt_tension1_1 = 2001, & + iopt_tension2_1 = 2002, & + iopt_onestate_1 = 2003, & + iopt_tens2pll_2 = 3001, & + iopt_unlimfrc_2 = 3002, & + iopt_unlimpow_2 = 3003, & + iopt_fixedsiz_2 = 3004, & + iopt_topmdexp_2 = 3005, & + iopt_arno_x_vic = 4001, & + iopt_prms_varnt = 4002, & + iopt_tmdl_param = 4003, & + iopt_perc_f2sat = 5001, & + iopt_perc_w2sat = 5002, & + iopt_perc_lower = 5003, & + iopt_sequential = 6001, & + iopt_rootweight = 6002, & + iopt_intflwnone = 7001, & + iopt_intflwsome = 7002, & + iopt_rout_gamma = 8001, & + iopt_no_routing = 8002, & + iopt_no_snowmod = 8501, & + iopt_temp_index = 8502, & + iopt_no_intrcep = 8601, & + iopt_gr5h_intrc = 8602 +! --- +integer(I4B), parameter :: iopt_TENS1A = 9001, & + iopt_TENS1B = 9002, & + iopt_TENS_1 = 9003, & + iopt_FREE_1 = 9004, & + iopt_WATR_1 = 9005, & + iopt_TENS_2 = 9006, & + iopt_FREE2A = 9007, & + iopt_FREE2B = 9008, & + iopt_WATR_2 = 9009 +! ------------------------------------------ +contains +! ------------------------------------------ +elemental function desc_str2int(name)result(res) +! Purpose: Converts a string description into its corresponding integer value. +implicit none +! dummies +character(*), intent(in) :: name +integer(I4B) :: res +! Start procedure here +selectcase(name) +case("additive_e"); res = iopt_additive_e +case("multiplc_e"); res = iopt_multiplc_e +case("tension1_1"); res = iopt_tension1_1 +case("tension2_1"); res = iopt_tension2_1 +case("onestate_1"); res = iopt_onestate_1 +case("tens2pll_2"); res = iopt_tens2pll_2 +case("unlimfrc_2"); res = iopt_unlimfrc_2 +case("unlimpow_2"); res = iopt_unlimpow_2 +case("fixedsiz_2"); res = iopt_fixedsiz_2 +case("arno_x_vic"); res = iopt_arno_x_vic +case("prms_varnt"); res = iopt_prms_varnt +case("tmdl_param"); res = iopt_tmdl_param +case("perc_f2sat"); res = iopt_perc_f2sat +case("perc_w2sat"); res = iopt_perc_w2sat +case("perc_lower"); res = iopt_perc_lower +case("sequential"); res = iopt_sequential +case("rootweight"); res = iopt_rootweight +case("intflwnone"); res = iopt_intflwnone +case("intflwsome"); res = iopt_intflwsome +case("rout_gamma"); res = iopt_rout_gamma +case("no_routing"); res = iopt_no_routing +case("no_snowmod"); res = iopt_no_snowmod +case("temp_index"); res = iopt_temp_index +case("no_intrcep"); res = iopt_no_intrcep +case("gr5h_intrc"); res = iopt_gr5h_intrc +case("TENS1B"); res = iopt_TENS1B +case("TENS_1"); res = iopt_TENS_1 +case("FREE_1"); res = iopt_FREE_1 +case("WATR_1"); res = iopt_WATR_1 +case("TENS_2"); res = iopt_TENS_2 +case("FREE2A"); res = iopt_FREE2A +case("FREE2B"); res = iopt_FREE2B +case("WATR_2"); res = iopt_WATR_2 +case default; res = -999 +endselect +! End procedure here +endfunction desc_str2int +! ------------------------------------------ +elemental function desc_int2str(intVal)result(res) +! Purpose: Converts an integer description into corresponding string value +implicit none +! dummies +integer(I4B), intent(in) :: intVal +character(10) :: res +! Start procedure here +selectcase(intVal) +case(iopt_TENS1B); res = "TENS1B" +case(iopt_TENS_1); res = "TENS_1" +case(iopt_FREE_1); res = "FREE_1" +case(iopt_WATR_1); res = "WATR_1" +case(iopt_TENS_2); res = "TENS_2" +case(iopt_FREE2A); res = "FREE2A" +case(iopt_FREE2B); res = "FREE2B" +case(iopt_WATR_2); res = "WATR_2" +case(iopt_additive_e); res = "additive_e" +case(iopt_multiplc_e); res = "multiplc_e" +case(iopt_tension1_1); res = "tension1_1" +case(iopt_tension2_1); res = "tension2_1" +case(iopt_onestate_1); res = "onestate_1" +case(iopt_tens2pll_2); res = "tens2pll_2" +case(iopt_unlimfrc_2); res = "unlimfrc_2" +case(iopt_unlimpow_2); res = "unlimpow_2" +case(iopt_fixedsiz_2); res = "fixedsiz_2" +case(iopt_arno_x_vic); res = "arno_x_vic" +case(iopt_prms_varnt); res = "prms_varnt" +case(iopt_tmdl_param); res = "tmdl_param" +case(iopt_perc_f2sat); res = "perc_f2sat" +case(iopt_perc_w2sat); res = "perc_w2sat" +case(iopt_perc_lower); res = "perc_lower" +case(iopt_sequential); res = "sequential" +case(iopt_rootweight); res = "rootweight" +case(iopt_intflwnone); res = "intflwnone" +case(iopt_intflwsome); res = "intflwsome" +case(iopt_rout_gamma); res = "rout_gamma" +case(iopt_no_routing); res = "no_routing" +case(iopt_no_snowmod); res = "no_snowmod" +case(iopt_temp_index); res = "temp_index" +case(iopt_no_intrcep); res = "no_intrcep" +case(iopt_gr5h_intrc); res = "gr5h_intrc" +case default; res = "UNDFND" +endselect +! End procedure here +endfunction desc_int2str +! ------------------------------------------ +endmodule model_defnames diff --git a/build/FUSE_SRC/types/model_defn_types.f90 b/build/FUSE_SRC/types/model_defn_types.f90 index a22acf9..8fb1e7b 100644 --- a/build/FUSE_SRC/types/model_defn_types.f90 +++ b/build/FUSE_SRC/types/model_defn_types.f90 @@ -6,6 +6,7 @@ MODULE model_defn_types ! Martyn Clark ! Modified by Brian Henn to include snow model, 6/2013 ! Modified by Martyn Clark to separate data tyoes from data store, 01/2026 + ! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- USE nrtype @@ -33,6 +34,7 @@ MODULE model_defn_types INTEGER(I4B) :: iQINTF INTEGER(I4B) :: iQ_TDH INTEGER(I4B) :: iSNOWM ! snow + INTEGER(I4B) :: iINTRC ! interception END TYPE UMODEL ! structure to hold model state names diff --git a/build/FUSE_SRC/types/multi_flux_types.f90 b/build/FUSE_SRC/types/multi_flux_types.f90 index 82f3fc0..68176ba 100644 --- a/build/FUSE_SRC/types/multi_flux_types.f90 +++ b/build/FUSE_SRC/types/multi_flux_types.f90 @@ -8,7 +8,10 @@ MODULE multi_flux_types public :: FLUXES TYPE FLUXES - REAL(WP) :: EFF_PPT ! effective precipitation (mm day-1) + REAL(WP) :: EVAP_0 ! evaporation from interception store (mm day-1) + REAL(WP) :: PIN0 ! precipitation entering interception store (mm day-1) + REAL(WP) :: PTHRU ! precipitation through interception store (mm day-1) + REAL(WP) :: EFF_PPT ! effective precipitation available to soil (mm day-1) REAL(WP) :: SATAREA ! saturated area (-) REAL(WP) :: QSURF ! surface runoff (mm day-1) REAL(WP) :: EVAP_1A ! evaporation from soil excess zone (mm day-1) diff --git a/build/FUSE_SRC/types/multiparam_types.f90 b/build/FUSE_SRC/types/multiparam_types.f90 index 23c2746..0efcad4 100644 --- a/build/FUSE_SRC/types/multiparam_types.f90 +++ b/build/FUSE_SRC/types/multiparam_types.f90 @@ -6,6 +6,7 @@ MODULE multiparam_types ! Martyn Clark ! Modified by Brian Henn to include snow model, 6/2013 ! Modified by Martyn Clark to separate type definitions from data storage, 01/2026 + ! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- USE nrtype @@ -50,6 +51,8 @@ MODULE multiparam_types TYPE(PARATT) :: RH1P_SDEV ! prior param2 of hyper param1: prior sdev of hypermean TYPE(PARATT) :: RH2P_MEAN ! prior param1 of hyper param2: lower bound of hypersdev TYPE(PARATT) :: RH2P_SDEV ! prior param2 of hyper param2: upper bound of hypersdev + ! interception (adjustable) + TYPE(PARATT) :: REFSINT_0 ! characteristic interception storage (50% throughfall) (mm) ! bucket sizes (adjustable) TYPE(PARATT) :: MAXWATR_1 ! maximum total storage in layer1 (mm) TYPE(PARATT) :: MAXWATR_2 ! maximum total storage in layer2 (mm) @@ -102,6 +105,8 @@ MODULE multiparam_types REAL(WP) :: RH1P_SDEV ! prior param2 of hyper param1: prior sdev of hypermean REAL(WP) :: RH2P_MEAN ! prior param1 of hyper param2: lower bound of hypersdev REAL(WP) :: RH2P_SDEV ! prior param2 of hyper param2: upper bound of hypersdev + ! interception + REAL(WP) :: REFSINT_0 ! characteristic interception storage (50% throughfall) (mm) ! bucket sizes (adjustable) REAL(WP) :: MAXWATR_1 ! maximum total storage in layer1 (mm) REAL(WP) :: MAXWATR_2 ! maximum total storage in layer2 (mm) diff --git a/build/FUSE_SRC/types/multistate_types.f90 b/build/FUSE_SRC/types/multistate_types.f90 index 480e7a7..776a13b 100644 --- a/build/FUSE_SRC/types/multistate_types.f90 +++ b/build/FUSE_SRC/types/multistate_types.f90 @@ -11,6 +11,8 @@ MODULE multistate_types ! model state structure ! -------------------------------------------------------------------------------------- TYPE STATEV + ! interception layer + REAL(WP) :: SINT_0 ! interception storage (mm) ! snow layer REAL(WP) :: SWE_TOT ! total storage as snow (mm) ! upper layer diff --git a/build/FUSE_SRC/util/getpar_str.f90 b/build/FUSE_SRC/util/getpar_str.f90 index 4b51b3e..8bdb73f 100644 --- a/build/FUSE_SRC/util/getpar_str.f90 +++ b/build/FUSE_SRC/util/getpar_str.f90 @@ -7,6 +7,7 @@ SUBROUTINE GETPAR_STR(PARNAME,METADAT) ! -------- ! Martyn Clark, 2009 ! Modified by Brian Henn to include snow model, 6/2013 +! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -29,6 +30,7 @@ SUBROUTINE GETPAR_STR(PARNAME,METADAT) CASE('RH1P_SDEV'); METADAT = PARMETA%RH1P_SDEV CASE('RH2P_MEAN'); METADAT = PARMETA%RH2P_MEAN CASE('RH2P_SDEV'); METADAT = PARMETA%RH2P_SDEV +CASE('REFSINT_0'); METADAT = PARMETA%REFSINT_0 CASE('MAXWATR_1'); METADAT = PARMETA%MAXWATR_1 CASE('MAXWATR_2'); METADAT = PARMETA%MAXWATR_2 CASE('FRACTEN'); METADAT = PARMETA%FRACTEN diff --git a/build/FUSE_SRC/util/metaparams.f90 b/build/FUSE_SRC/util/metaparams.f90 index 41cc6dd..9267d59 100644 --- a/build/FUSE_SRC/util/metaparams.f90 +++ b/build/FUSE_SRC/util/metaparams.f90 @@ -6,6 +6,7 @@ MODULE metaparams ! Martyn Clark, 2007 ! Modified by Brian Henn to include snow model, 6/2013 ! Modified by Martyn Clark to avoid per-band parameters, 12/2025 + ! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -40,6 +41,7 @@ SUBROUTINE PARDESCRIBE() ! adjustable model parameters I=I+1; PNAME(I)='RFERR_ADD '; PDESC(I)='additive rainfall error '; PUNIT(I)='mm '; isBand(i)=.false. I=I+1; PNAME(I)='RFERR_MLT '; PDESC(I)='multiplicative rainfall error '; PUNIT(I)='- '; isBand(i)=.false. + I=I+1; PNAME(I)='REFSINT_0 '; PDESC(I)='characteristic interception storage '; PUNIT(I)='mm '; isBand(i)=.false. I=I+1; PNAME(I)='MAXWATR_1 '; PDESC(I)='maximum total storage in the upper layer '; PUNIT(I)='mm '; isBand(i)=.false. I=I+1; PNAME(I)='MAXWATR_2 '; PDESC(I)='maximum total storage in the lower layer '; PUNIT(I)='mm '; isBand(i)=.false. I=I+1; PNAME(I)='FRACTEN '; PDESC(I)='fraction total storage as tension storage '; PUNIT(I)='- '; isBand(i)=.false. diff --git a/build/FUSE_SRC/util/par_insert.f90 b/build/FUSE_SRC/util/par_insert.f90 index b5724ab..768b549 100644 --- a/build/FUSE_SRC/util/par_insert.f90 +++ b/build/FUSE_SRC/util/par_insert.f90 @@ -41,6 +41,7 @@ SUBROUTINE PAR_INSERT(XVAR,PARNAME) ! Creator: ! -------- ! Martyn Clark, 2007 +! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -63,6 +64,7 @@ SUBROUTINE PAR_INSERT(XVAR,PARNAME) CASE('RH1P_SDEV'); MPARAM%RH1P_SDEV = XVAR CASE('RH2P_MEAN'); MPARAM%RH2P_MEAN = XVAR CASE('RH2P_SDEV'); MPARAM%RH2P_SDEV = XVAR +CASE('REFSINT_0'); MPARAM%REFSINT_0 = XVAR CASE('MAXWATR_1'); MPARAM%MAXWATR_1 = XVAR CASE('MAXWATR_2'); MPARAM%MAXWATR_2 = XVAR CASE('FRACTEN'); MPARAM%FRACTEN = XVAR diff --git a/build/FUSE_SRC/util/parextract.f90 b/build/FUSE_SRC/util/parextract.f90 index ec15b75..d09f2d9 100644 --- a/build/FUSE_SRC/util/parextract.f90 +++ b/build/FUSE_SRC/util/parextract.f90 @@ -18,6 +18,7 @@ PURE FUNCTION PAREXTRACT(PARNAME) ! -------- ! Martyn Clark, 2007 ! Modified by Martyn Clark to remove elevation band parameters (handled separately) + ! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -46,6 +47,7 @@ PURE FUNCTION PAREXTRACT(PARNAME) CASE ('RH1P_SDEV') ; XVAR = MPARAM%RH1P_SDEV CASE ('RH2P_MEAN') ; XVAR = MPARAM%RH2P_MEAN CASE ('RH2P_SDEV') ; XVAR = MPARAM%RH2P_SDEV + CASE ('REFSINT_0') ; XVAR = MPARAM%REFSINT_0 CASE ('MAXWATR_1') ; XVAR = MPARAM%MAXWATR_1 CASE ('MAXWATR_2') ; XVAR = MPARAM%MAXWATR_2 CASE ('FRACTEN') ; XVAR = MPARAM%FRACTEN diff --git a/build/FUSE_SRC/util/putpar_str.f90 b/build/FUSE_SRC/util/putpar_str.f90 index 139aea3..f3fb69c 100644 --- a/build/FUSE_SRC/util/putpar_str.f90 +++ b/build/FUSE_SRC/util/putpar_str.f90 @@ -7,6 +7,7 @@ SUBROUTINE PUTPAR_STR(METADAT,PARNAME) ! -------- ! Martyn Clark, 2009 ! Modified by Brian Henn to include snow model, 6/2013 +! Modified by Cyril Thebault to include interception, 7/2026 ! --------------------------------------------------------------------------------------- ! Purpose: ! -------- @@ -29,6 +30,7 @@ SUBROUTINE PUTPAR_STR(METADAT,PARNAME) CASE('RH1P_SDEV'); PARMETA%RH1P_SDEV = METADAT CASE('RH2P_MEAN'); PARMETA%RH2P_MEAN = METADAT CASE('RH2P_SDEV'); PARMETA%RH2P_SDEV = METADAT +CASE('REFSINT_0'); PARMETA%REFSINT_0 = METADAT CASE('MAXWATR_1'); PARMETA%MAXWATR_1 = METADAT CASE('MAXWATR_2'); PARMETA%MAXWATR_2 = METADAT CASE('FRACTEN'); PARMETA%FRACTEN = METADAT diff --git a/build/FUSE_SRC/util/selectmodl.f90 b/build/FUSE_SRC/util/selectmodl.f90 index a57518b..d46b1e9 100644 --- a/build/FUSE_SRC/util/selectmodl.f90 +++ b/build/FUSE_SRC/util/selectmodl.f90 @@ -17,7 +17,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) USE fuse_fileManager,only:SETNGS_PATH,M_DECISIONS ! defines data directory USE model_defn,ONLY:NDEC,SMODL,AMODL,& ! defines model decisions LIST_RFERR,LIST_ARCH1,LIST_ARCH2,LIST_QSURF,LIST_QPERC,LIST_ESOIL,& - LIST_QINTF,LIST_Q_TDH,LIST_SNOWM + LIST_QINTF,LIST_Q_TDH,LIST_SNOWM,LIST_INTRC USE model_defnames,ONLY:DESC_STR2INT USE model_numerix,only:solution_method,temporal_error_control IMPLICIT NONE @@ -53,6 +53,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) INTEGER(I4B) :: ISW_QINTF ! loop thru interflow INTEGER(I4B) :: ISW_Q_TDH ! loop thru time delay options INTEGER(I4B) :: ISW_SNOWM ! loop thru snow model options +INTEGER(I4B) :: ISW_INTRC ! loop thru interception options INTEGER(I4B) :: IX_RFERR ! index for rainfall error options INTEGER(I4B) :: IX_ARCH1 ! index for upper layer architecture INTEGER(I4B) :: IX_ARCH2 ! index for lower layer architecture @@ -62,6 +63,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) INTEGER(I4B) :: IX_QINTF ! index for interflow INTEGER(I4B) :: IX_Q_TDH ! index for time delay options INTEGER(I4B) :: IX_SNOWM ! index for snow model options +INTEGER(I4B) :: IX_INTRC ! index for interception options INTEGER(I4B) :: IX_MODEL ! model index ! (3) identify a unique model name CHARACTER(LEN=8) :: CNUM ! model index (converted to text) @@ -73,7 +75,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) ICOUNT =0 IX_MODEL=0 IX_Q_TDH=0; IX_QINTF=0; IX_ESOIL=0; IX_QPERC=0; IX_QSURF=0; IX_ARCH2=0; IX_ARCH1=0; IX_RFERR=0 -IX_SNOWM = 0 +IX_SNOWM = 0; IX_INTRC = 0 ERR =0 MESSAGE ='SELECTMODL/everything is fine' ! --------------------------------------------------------------------------------------- @@ -81,8 +83,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) ! --------------------------------------------------------------------------------------- ! read in control file -!CFILE = TRIM(SETNGS_PATH)//M_DECISIONS ! control file info shared in MODULE ddirectory -CFILE = TRIM(SETNGS_PATH)//'fuse_zDecisions_'//TRIM(FUSE_ID)//'.txt' ! control file info shared in MODULE ddirectory +CFILE = TRIM(SETNGS_PATH)//TRIM(M_DECISIONS) ! control file info shared in MODULE ddirectory INQUIRE(FILE=trim(CFILE),EXIST=LEXIST) ! check that control file exists IF (.not.LEXIST) THEN message="f-SELECTMODL/decisions file '"//trim(CFILE)//"' does not exist" @@ -136,6 +137,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) CASE('QINTF'); SMODL%iQINTF = desc_str2int(M_CHOICE) CASE('Q_TDH'); SMODL%iQ_TDH = desc_str2int(M_CHOICE) CASE('SNOWM'); SMODL%iSNOWM = desc_str2int(M_CHOICE) + CASE('INTRC'); SMODL%iINTRC = desc_str2int(M_CHOICE) CASE DEFAULT message="f-SELECTMODL/UNRECOGNISED[DECISON='"//TRIM(DECISION)//"']& &[M_CHOICE='"//TRIM(M_CHOICE)//"']" @@ -153,56 +155,60 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) ! --------------------------------------------------------------------------------------- ! loop through model options MODEL_OPTIONS: DO ISW_MODEL=1,1 ! (dummy loop to exit) -DO ISW_SNOWM=1,SIZE(LIST_SNOWM) ! (snow model options) - DO ISW_Q_TDH=1,SIZE(LIST_Q_TDH) ! (time delay options) - DO ISW_QINTF=1,SIZE(LIST_QINTF) ! (interflow options) - DO ISW_ESOIL=1,SIZE(LIST_ESOIL) ! (evaporation options) - DO ISW_QPERC=1,SIZE(LIST_QPERC) ! (percolation options) - DO ISW_QSURF=1,SIZE(LIST_QSURF) ! (surface runoff options) - DO ISW_ARCH2=1,SIZE(LIST_ARCH2) ! (lower-layer architecture options) - DO ISW_ARCH1=1,SIZE(LIST_ARCH1) ! (upper-layer architecture options) - DO ISW_RFERR=1,SIZE(LIST_RFERR) ! (rainfall error options) - ! don't allow a lower tension tank when there are two upper ones - IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT.EQ.'tension2_1'.AND. & - LIST_ARCH2(ISW_ARCH2)%MCOMPONENT.EQ.'tens2pll_2') CYCLE - ! don't allow percolation below field capacity if there are multiple upper tanks - IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT.NE.'onestate_1'.AND. & - LIST_QPERC(ISW_QPERC)%MCOMPONENT.EQ.'perc_w2sat') CYCLE - ICOUNT = ICOUNT + 1 ! (increment counter) - ! identify unique model - IF (SMODL%iRFERR.EQ.desc_str2int(LIST_RFERR(ISW_RFERR)%MCOMPONENT) .AND. & - SMODL%iARCH1.EQ.desc_str2int(LIST_ARCH1(ISW_ARCH1)%MCOMPONENT) .AND. & - SMODL%iARCH2.EQ.desc_str2int(LIST_ARCH2(ISW_ARCH2)%MCOMPONENT) .AND. & - SMODL%iQSURF.EQ.desc_str2int(LIST_QSURF(ISW_QSURF)%MCOMPONENT) .AND. & - SMODL%iQPERC.EQ.desc_str2int(LIST_QPERC(ISW_QPERC)%MCOMPONENT) .AND. & - SMODL%iESOIL.EQ.desc_str2int(LIST_ESOIL(ISW_ESOIL)%MCOMPONENT) .AND. & - SMODL%iQINTF.EQ.desc_str2int(LIST_QINTF(ISW_QINTF)%MCOMPONENT) .AND. & - SMODL%iQ_TDH.EQ.desc_str2int(LIST_Q_TDH(ISW_Q_TDH)%MCOMPONENT) .AND. & - SMODL%iSNOWM.EQ.desc_str2int(LIST_SNOWM(ISW_SNOWM)%MCOMPONENT)) THEN - ! identify model components - IX_RFERR = ISW_RFERR - IX_ARCH1 = ISW_ARCH1 - IX_ARCH2 = ISW_ARCH2 - IX_QSURF = ISW_QSURF - IX_QPERC = ISW_QPERC - IX_ESOIL = ISW_ESOIL - IX_QINTF = ISW_QINTF - IX_Q_TDH = ISW_Q_TDH - IX_SNOWM = ISW_SNOWM - ! identify model - IX_MODEL = ICOUNT - ! exit main do loop - EXIT MODEL_OPTIONS - ENDIF - END DO ! RFERR - END DO ! ARCH1 - END DO ! ARCH2 - END DO ! QSURF - END DO ! QPERC - END DO ! ESOIL - END DO ! QINTF - END DO ! Q_TDH -END DO ! SNOWM +DO ISW_INTRC=1,SIZE(LIST_INTRC) ! (interception options) + DO ISW_SNOWM=1,SIZE(LIST_SNOWM) ! (snow model options) + DO ISW_Q_TDH=1,SIZE(LIST_Q_TDH) ! (time delay options) + DO ISW_QINTF=1,SIZE(LIST_QINTF) ! (interflow options) + DO ISW_ESOIL=1,SIZE(LIST_ESOIL) ! (evaporation options) + DO ISW_QPERC=1,SIZE(LIST_QPERC) ! (percolation options) + DO ISW_QSURF=1,SIZE(LIST_QSURF) ! (surface runoff options) + DO ISW_ARCH2=1,SIZE(LIST_ARCH2) ! (lower-layer architecture options) + DO ISW_ARCH1=1,SIZE(LIST_ARCH1) ! (upper-layer architecture options) + DO ISW_RFERR=1,SIZE(LIST_RFERR) ! (rainfall error options) + ! don't allow a lower tension tank when there are two upper ones + IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT.EQ.'tension2_1'.AND. & + LIST_ARCH2(ISW_ARCH2)%MCOMPONENT.EQ.'tens2pll_2') CYCLE + ! don't allow percolation below field capacity if there are multiple upper tanks + IF (LIST_ARCH1(ISW_ARCH1)%MCOMPONENT.NE.'onestate_1'.AND. & + LIST_QPERC(ISW_QPERC)%MCOMPONENT.EQ.'perc_w2sat') CYCLE + ICOUNT = ICOUNT + 1 ! (increment counter) + ! identify unique model + IF (SMODL%iRFERR.EQ.desc_str2int(LIST_RFERR(ISW_RFERR)%MCOMPONENT) .AND. & + SMODL%iARCH1.EQ.desc_str2int(LIST_ARCH1(ISW_ARCH1)%MCOMPONENT) .AND. & + SMODL%iARCH2.EQ.desc_str2int(LIST_ARCH2(ISW_ARCH2)%MCOMPONENT) .AND. & + SMODL%iQSURF.EQ.desc_str2int(LIST_QSURF(ISW_QSURF)%MCOMPONENT) .AND. & + SMODL%iQPERC.EQ.desc_str2int(LIST_QPERC(ISW_QPERC)%MCOMPONENT) .AND. & + SMODL%iESOIL.EQ.desc_str2int(LIST_ESOIL(ISW_ESOIL)%MCOMPONENT) .AND. & + SMODL%iQINTF.EQ.desc_str2int(LIST_QINTF(ISW_QINTF)%MCOMPONENT) .AND. & + SMODL%iQ_TDH.EQ.desc_str2int(LIST_Q_TDH(ISW_Q_TDH)%MCOMPONENT) .AND. & + SMODL%iSNOWM.EQ.desc_str2int(LIST_SNOWM(ISW_SNOWM)%MCOMPONENT) .AND. & + SMODL%iINTRC == desc_str2int(LIST_INTRC(ISW_INTRC)%MCOMPONENT)) THEN + ! identify model components + IX_RFERR = ISW_RFERR + IX_ARCH1 = ISW_ARCH1 + IX_ARCH2 = ISW_ARCH2 + IX_QSURF = ISW_QSURF + IX_QPERC = ISW_QPERC + IX_ESOIL = ISW_ESOIL + IX_QINTF = ISW_QINTF + IX_Q_TDH = ISW_Q_TDH + IX_SNOWM = ISW_SNOWM + IX_INTRC = ISW_INTRC + ! identify model + IX_MODEL = ICOUNT + ! exit main do loop + EXIT MODEL_OPTIONS + ENDIF + END DO ! RFERR + END DO ! ARCH1 + END DO ! ARCH2 + END DO ! QSURF + END DO ! QPERC + END DO ! ESOIL + END DO ! QINTF + END DO ! Q_TDH + END DO ! SNOWM +END DO ! INTRC END DO MODEL_OPTIONS ! check that a model was identified @@ -231,6 +237,7 @@ SUBROUTINE SELECTMODL(FUSE_ID,ERR,MESSAGE) DECISION='QINTF'; CALL MAKE_MODEL(IX_QINTF,DECISION,SMODL%MNAME) DECISION='Q_TDH'; CALL MAKE_MODEL(IX_Q_TDH,DECISION,SMODL%MNAME) DECISION='SNOWM'; CALL MAKE_MODEL(IX_SNOWM,DECISION,SMODL%MNAME) + DECISION='INTRC'; CALL MAKE_MODEL(IX_INTRC,DECISION,SMODL%MNAME) END SELECT ! add the numerix info DECISION='NMETH'; CALL MAKE_MODEL(SOLUTION_METHOD,DECISION,SMODL%MNAME) diff --git a/build/Makefile b/build/Makefile index 37b295e..6fe45df 100755 --- a/build/Makefile +++ b/build/Makefile @@ -285,6 +285,7 @@ PHYSICS = $(patsubst %, $(PHYSICS_DIR)/%, $(FUSE_PHYSICS)) FUSE_MODGUT = FUSE_MODGUT += mod_derivs.f90 FUSE_MODGUT += update_swe.f90 +FUSE_MODGUT += update_interception.f90 FUSE_MODGUT += qrainerror.f90 FUSE_MODGUT += qsatexcess.f90 FUSE_MODGUT += evap_upper.f90 diff --git a/test/CAN_05BB001/lumped/settings/fuse_v1/fuse_zDecisions_2.txt b/test/CAN_05BB001/lumped/settings/fuse_v1/fuse_zDecisions_2.txt index b26f613..2282bb6 100644 --- a/test/CAN_05BB001/lumped/settings/fuse_v1/fuse_zDecisions_2.txt +++ b/test/CAN_05BB001/lumped/settings/fuse_v1/fuse_zDecisions_2.txt @@ -8,6 +8,7 @@ sequential ESOIL ! (6) evaporation intflwnone QINTF ! (7) interflow rout_gamma Q_TDH ! (8) time delay in runoff temp_index SNOWM ! (9) snow model +no_intrcep INTRC ! (10) interception 0 ! naming convention for model (0=full, 1=index + numrx) --------------------------------------------------------------------------- --------------------------------------------------------------------------- @@ -54,3 +55,7 @@ no_routing ! no routing temp_index ! temperature index snow model no_snowmod ! no snow model --------------------------------------------------------------------------- +(10) interception +no_intrcep ! no interception store +gr5h_intrc ! gr5h interception store +--------------------------------------------------------------------------- \ No newline at end of file diff --git a/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zConstraints_snow.txt b/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zConstraints_snow.txt index b255cfd..c511d78 100644 --- a/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zConstraints_snow.txt +++ b/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zConstraints_snow.txt @@ -7,6 +7,7 @@ F 0 0.000 99.000 -99.000 .10 1.0 0 0 0 0 0 0 RH1P_MEAN NO_CHILD1 NO_ F 0 0.400 99.000 -99.000 .10 1.0 0 0 0 0 0 0 RH1P_SDEV NO_CHILD1 NO_CHILD2 ! prior_par2/hyp_par1: prior sdev of hypermean F 0 0.050 99.000 -99.000 .10 1.0 0 0 0 0 0 0 RH2P_MEAN NO_CHILD1 NO_CHILD2 ! prior_par1/hyp_par2: lower bound of hypersdev F 0 0.700 99.000 -99.000 .10 1.0 0 0 0 0 0 0 RH2P_SDEV NO_CHILD1 NO_CHILD2 ! prior_par2/hyp_par2: upper bound of hypersdev +T 0 2.000 0.001 20.000 .10 1.0 0 0 0 0 0 0 REFSINT_0 NO_CHILD1 NO_CHILD2 ! characteristic interception storage (mm) T 0 477.010 25.000 500.000 .10 10.0 0 0 0 0 0 0 MAXWATR_1 NO_CHILD1 NO_CHILD2 ! depth of the upper soil layer (mm) T 0 3314.525 50.000 5000.000 .10 100.0 0 0 0 0 0 0 MAXWATR_2 NO_CHILD1 NO_CHILD2 ! depth of the lower soil layer (mm) T 0 0.856 0.050 0.950 .10 0.1 0 0 0 0 0 0 FRACTEN NO_CHILD1 NO_CHILD2 ! fraction total storage in tension storage (-) diff --git a/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zDecisions_2.txt b/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zDecisions_2.txt index b26f613..2282bb6 100644 --- a/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zDecisions_2.txt +++ b/test/CAN_05BB001/lumped/settings/fuse_v2/fuse_zDecisions_2.txt @@ -8,6 +8,7 @@ sequential ESOIL ! (6) evaporation intflwnone QINTF ! (7) interflow rout_gamma Q_TDH ! (8) time delay in runoff temp_index SNOWM ! (9) snow model +no_intrcep INTRC ! (10) interception 0 ! naming convention for model (0=full, 1=index + numrx) --------------------------------------------------------------------------- --------------------------------------------------------------------------- @@ -54,3 +55,7 @@ no_routing ! no routing temp_index ! temperature index snow model no_snowmod ! no snow model --------------------------------------------------------------------------- +(10) interception +no_intrcep ! no interception store +gr5h_intrc ! gr5h interception store +--------------------------------------------------------------------------- \ No newline at end of file