diff --git a/parameters/MPTABLE.TBL b/parameters/MPTABLE.TBL index dfceae24..50f767c1 100644 --- a/parameters/MPTABLE.TBL +++ b/parameters/MPTABLE.TBL @@ -551,6 +551,39 @@ RTCT_S8 = 0.0, 0.0, 0.0, 0.0, 0.0, BIO2LAI = 0.015, 0.030, 0.015, 0.015, 0.015, ! leaf are per living leaf biomass [m^2/kg] +/ + +&noahmp_tiledrain_parameters +!-----------------------------------! +! For simple drainage model ! +!-----------------------------------! +DRAIN_LAYER_OPT = 4 + ! 0 - from one specified layer by TD_DEPTH, + ! 1 - from layers 1 & 2, + ! 2 - from layer layers 1, 2, and 3 + ! 3 - from layer 2 and 3 + ! 4 - from layer layers 3, 4 + ! 5 - from all the four layers +!-------------------------------------------------------------------------------------------------------------------------------------------------------------------- ! +! 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ! +!-------------------------------------------------------------------------------------------------------------------------------------------------------------------- ! +TDSMC_FAC = 0.90, 0.90, 0.90, 0.90, 0.90, 1.25, 0.90, 1.0, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90 ! corresponds to number of soil types SOILPARAM.TBL +TD_DEPTH = 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 ! depth of drain tube from the soil surface +TD_DC = 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20. ! drainage coefficient (mm d^-1) +!-------------------------------------------------------------------------------------------------------------------------------------------------------------------- ! +! +!-------------------------------------! +! For Hooghoudt tile drain model ! +!-------------------------------------! +!-------------------------------------------------------------------------------------------------------------------------------------------------------------------- +TD_DCOEF = 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07 ! m d^-1, drainage coefficent +TD_D = 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00 ! m, depth to impe layer from drain water level (D) +TD_ADEPTH = 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00 ! m, actual depth of imp layer from land surface +TD_RADI = 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07 ! m, effective radius of drains (ro) +TD_SPAC = 60.0, 55.0, 45.0, 20.0, 25.0, 30.0, 40.0, 16.0, 18.0, 50.0, 15.0, 10.0, 35.0, 10.0, 60.0, 60.0, 10.0, 60.0, 60.0 ! m, distance between two drain tubes or tiles (L) +TD_DDRAIN = 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20, 1.20 ! m, Depth of drain +KLAT_FAC = 1.30, 1.80, 2.10, 2.60, 2.90, 2.50, 2.30, 3.00, 2.70, 2.00, 3.10, 3.30, 2.50, 1.00, 1.00, 1.80, 4.00, 1.00, 1.30 ! multiplication factor to lateral hyd.cond +!-------------------------------------------------------------------------------------------------------------------------------------------------------------------- / diff --git a/src/module_sf_noahmplsm.F b/src/module_sf_noahmplsm.F index 0da4e881..1a799d99 100644 --- a/src/module_sf_noahmplsm.F +++ b/src/module_sf_noahmplsm.F @@ -114,6 +114,9 @@ MODULE MODULE_SF_NOAHMPLSM ! 4 -> BATS surface and subsurface runoff (free drainage) ! 5 -> Miguez-Macho&Fan groundwater scheme (Miguez-Macho et al. 2007 JGR; Fan et al. 2007 JGR) ! (needs further testing for public use) + ! 6 -> Variable Infiltration Capacity Model surface runoff scheme (Wood et al., 1992, JGR) + ! 7 -> Xiananjiang Infiltration and surface runoff scheme ((Jayawardena and Zhou, 2000) + ! 8 -> Dynamic VIC surface runoff scheme (Liang and Xie, 2001) INTEGER :: OPT_SFC ! options for surface layer drag coeff (CH & CM) ! **1 -> M-O @@ -185,6 +188,16 @@ MODULE MODULE_SF_NOAHMPLSM ! 2 -> micro/drip irrigation ! 3 -> surface flooding + INTEGER :: OPT_INFDV! options for infiltration in dynamic VIC runoff scheme + ! **1 -> Philip scheme + ! 2 -> Green-Ampt scheme + ! 3 -> Smith-Parlange scheme + + INTEGER :: OPT_TDRN ! options for crop model (currently only tested & calibrated to work with opt_run=3) + ! **0 -> No tile drainage + ! 1 -> on (simple scheme) + ! 2 -> on (Hooghoudt's scheme) + !------------------------------------------------------------------------------------------! ! Physical Constants: ! !------------------------------------------------------------------------------------------! @@ -395,6 +408,14 @@ MODULE MODULE_SF_NOAHMPLSM REAL :: DWSAT(NSOIL) !saturated soil hydraulic diffusivity REAL :: QUARTZ(NSOIL) !soil quartz content REAL :: F1 !soil thermal diffusivity/conductivity coef (not used MB: 20140718) + REAL :: BVIC !VIC model infiltration parameter for opt_run=6 + REAL :: AXAJ !Xinanjiang: Tension water distribution inflection parameter [-] for opt_run=7 + REAL :: BXAJ !Xinanjiang: Tension water distribution shape parameter [-] for opt_run=7 + REAL :: XXAJ !Xinanjiang: Free water distribution shape parameter [-] for opt_run=7 + REAL :: BDVIC !DVIC model infiltration parameter [-] for opt_run=8 + REAL :: GDVIC !Mean Capillary Drive (m) for infiltration models for opt_run=8 + REAL :: BBVIC !DVIC heterogeniety parameter for infiltration for opt_run=8 + !------------------------------------------------------------------------------------------! ! From the GENPARM.TBL file !------------------------------------------------------------------------------------------! @@ -408,6 +429,23 @@ MODULE MODULE_SF_NOAHMPLSM REAL :: KDT !used in compute maximum infiltration rate (in INFIL) REAL :: FRZX !used in compute maximum infiltration rate (in INFIL) +!------------------------------------------------------------------------------------------! +! From the tiledrain section of the MPTABLE.TBL file +!------------------------------------------------------------------------------------------! + REAL :: TDSMC_FAC + REAL :: TD_DC + INTEGER :: TD_DEPTH + INTEGER :: DRAIN_LAYER_OPT + + REAL :: TD_DCOEF ! m d^-1, drainage coefficent + REAL :: TD_D ! m, depth to impervious layer from drain water level (D) + REAL :: TD_ADEPTH ! m, actual depth of impervious layer from land surface + REAL :: TD_RADI ! m, effective radius of drains (ro) + REAL :: TD_SPAC ! m, distance between two drain tubes or tiles (L) + REAL :: TD_DDRAIN ! m, Depth of drain + REAL :: KLAT_FAC ! multiplication factor to determine conk(j1,j) from sol_k(j1,j) for HRU + + END TYPE noahmp_parameters contains @@ -435,6 +473,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & GRAIN , GDD , PGS , & ! IN/OUT SMCWTD ,DEEPRECH , RECH , & ! IN/OUT : GECROS1D, & ! IN/OUT : + QTLDRN , TDFRACMP, & ! IN/OUT : Z0WRF , & IRCNTSI , IRCNTMI , IRCNTFI , IRAMTSI , IRAMTMI , IRAMTFI , & ! IN/OUT : Irrigation: vars IRSIRATE, IRMIRATE, IRFIRATE, FIRR , EIRR , & ! IN/OUT : Irrigation: vars @@ -451,7 +490,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & CHLEAF , CHUC , CHV2 , CHB2 , FPICE , PAHV , & PAHG , PAHB , PAH , LAISUN , LAISHA , RB & ! OUT #ifdef WRF_HYDRO - ,SFCHEADRT & ! IN/OUT : + ,SFCHEADRT, WATBLED & ! IN/OUT : #endif ) @@ -508,7 +547,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & !jref:end #ifdef WRF_HYDRO - REAL , INTENT(INOUT) :: sfcheadrt + REAL , INTENT(INOUT) :: sfcheadrt, WATBLED #endif ! input/output : need arbitary intial values @@ -683,6 +722,10 @@ SUBROUTINE NOAHMP_SFLX (parameters, & REAL , INTENT(INOUT) :: GDD !growing degree days INTEGER , INTENT(INOUT) :: PGS !plant growing stage [-] +! Tile drainage + REAL , INTENT(INOUT) :: QTLDRN !tile drainage [mm/s] + REAL , INTENT(IN) :: TDFRACMP !tile drain fraction map + ! irrigation variables REAL , INTENT(IN) :: IRRFRA ! irrigation fraction REAL , INTENT(IN) :: SIFRA ! sprinkler irrigation fraction @@ -887,7 +930,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & SAV ,SAG ,QMELT ,FSA ,FSR ,TAUX , & !out TAUY ,FIRA ,FSH ,FCEV ,FGEV ,FCTR , & !out TRAD ,PSN ,APAR ,SSOIL ,BTRANI ,BTRAN , & !out - PONDING,TS ,LATHEAV , LATHEAG , frozen_canopy,frozen_ground, & !out + PONDING,TS ,LATHEAV , LATHEAG , frozen_canopy,frozen_ground, & !out TV ,TG ,STC ,SNOWH ,EAH ,TAH , & !inout SNEQVO ,SNEQV ,SH2O ,SMC ,SNICE ,SNLIQ , & !inout ALBOLD ,CM ,CH ,DX ,DZ8W ,Q2 , & !inout @@ -918,16 +961,17 @@ SUBROUTINE NOAHMP_SFLX (parameters, & FICEOLD,PONDING,TG ,IST ,FVEG ,iloc,jloc , SMCEQ , & !in BDFALL ,FP ,RAIN ,SNOW , & !in MB/AN: v3.7 QSNOW ,QRAIN ,SNOWHIN,LATHEAV,LATHEAG,frozen_canopy,frozen_ground, & !in MB + DX ,TDFRACMP, & !in PVK tile drainage ISNOW ,CANLIQ ,CANICE ,TV ,SNOWH ,SNEQV , & !inout SNICE ,SNLIQ ,STC ,ZSNSO ,SH2O ,SMC , & !inout SICE ,ZWT ,WA ,WT ,DZSNSO ,WSLAKE , & !inout SMCWTD ,DEEPRECH,RECH , & !inout IRAMTFI,IRAMTMI ,IRFIRATE ,IRMIRATE, & !inout CMC ,ECAN ,ETRAN ,FWET ,RUNSRF ,RUNSUB , & !out - QIN ,QDIS ,PONDING1 ,PONDING2,& - QSNBOT & + QIN ,QDIS ,PONDING1 ,PONDING2, & + QSNBOT ,QTLDRN & #ifdef WRF_HYDRO - ,sfcheadrt & + ,sfcheadrt,WATBLED & #endif ) !out @@ -974,6 +1018,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & FGEV ,FCTR ,SSOIL ,BEG_WB ,CANLIQ ,CANICE , & !in SNEQV ,WA ,SMC ,DZSNSO ,PRCP ,ECAN , & !in ETRAN ,EDIR ,RUNSRF ,RUNSUB ,DT ,NSOIL , & !in + QTLDRN , & !in NSNOW ,IST ,ERRWAT ,ILOC , JLOC ,FVEG , & SAV ,SAG ,FSRV ,FSRG ,ZWT ,PAH , & PAHV ,PAHG ,PAHB ,FIRR) !in ( Except ERRWAT, which is out ) @@ -1514,6 +1559,7 @@ SUBROUTINE ERROR (parameters,SWDOWN ,FSA ,FSR ,FIRA ,FSH ,FCEV , & FGEV ,FCTR ,SSOIL ,BEG_WB ,CANLIQ ,CANICE , & SNEQV ,WA ,SMC ,DZSNSO ,PRCP ,ECAN , & ETRAN ,EDIR ,RUNSRF ,RUNSUB ,DT ,NSOIL , & + QTLDRN , & NSNOW ,IST ,ERRWAT, ILOC ,JLOC ,FVEG , & SAV ,SAG ,FSRV ,FSRG ,ZWT ,PAH , & PAHV ,PAHG ,PAHB ,FIRR) @@ -1551,6 +1597,7 @@ SUBROUTINE ERROR (parameters,SWDOWN ,FSA ,FSR ,FIRA ,FSH ,FCEV , & REAL , INTENT(IN) :: EDIR !soil surface evaporation rate[mm/s] REAL , INTENT(IN) :: RUNSRF !surface runoff [mm/s] REAL , INTENT(IN) :: RUNSUB !baseflow (saturation excess) [mm/s] + REAL , INTENT(IN) :: QTLDRN !tile drainage [mm/s] REAL , INTENT(IN) :: CANLIQ !intercepted liquid water (mm) REAL , INTENT(IN) :: CANICE !intercepted ice mass (mm) REAL , INTENT(IN) :: SNEQV !snow water eqv. [mm] @@ -1635,7 +1682,7 @@ SUBROUTINE ERROR (parameters,SWDOWN ,FSA ,FSR ,FIRA ,FSH ,FCEV , & DO IZ = 1,NSOIL END_WB = END_WB + SMC(IZ) * DZSNSO(IZ) * 1000. END DO - ERRWAT = END_WB-BEG_WB-(PRCP-ECAN-ETRAN-EDIR-RUNSRF-RUNSUB)*DT + ERRWAT = END_WB-BEG_WB-(PRCP-ECAN-ETRAN-EDIR-RUNSRF-RUNSUB-QTLDRN)*DT #ifndef WRF_HYDRO IF(ABS(ERRWAT) > 0.1) THEN @@ -1647,10 +1694,10 @@ SUBROUTINE ERROR (parameters,SWDOWN ,FSA ,FSR ,FIRA ,FSH ,FCEV , & write(message, *) 'ERRWAT =',ERRWAT, "kg m{-2} timestep{-1}" call wrf_message(trim(message)) WRITE(message, & - '(" I J END_WB BEG_WB PRCP ECAN EDIR ETRAN RUNSRF RUNSUB")') + '(" I J END_WB BEG_WB PRCP ECAN EDIR ETRAN RUNSRF RUNSUB ZWT QTLDRN")') call wrf_message(trim(message)) WRITE(message,'(i6,1x,i6,1x,2f15.3,9f11.5)')ILOC,JLOC,END_WB,BEG_WB,PRCP*DT,ECAN*DT,& - EDIR*DT,ETRAN*DT,RUNSRF*DT,RUNSUB*DT,ZWT + EDIR*DT,ETRAN*DT,RUNSRF*DT,RUNSUB*DT,ZWT,QTLDRN*DT call wrf_message(trim(message)) call wrf_error_fatal("Water budget problem in NOAHMP LSM") END IF @@ -5868,6 +5915,7 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & FICEOLD,PONDING,TG ,IST ,FVEG ,ILOC ,JLOC ,SMCEQ , & !in BDFALL ,FP ,RAIN ,SNOW, & !in MB/AN: v3.7 QSNOW ,QRAIN ,SNOWHIN,LATHEAV,LATHEAG,frozen_canopy,frozen_ground, & !in MB + DX ,TDFRACMP, & ISNOW ,CANLIQ ,CANICE ,TV ,SNOWH ,SNEQV , & !inout SNICE ,SNLIQ ,STC ,ZSNSO ,SH2O ,SMC , & !inout SICE ,ZWT ,WA ,WT ,DZSNSO ,WSLAKE , & !inout @@ -5875,9 +5923,9 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & IRAMTFI,IRAMTMI ,IRFIRATE ,IRMIRATE, & !inout CMC ,ECAN ,ETRAN ,FWET ,RUNSRF ,RUNSUB , & !out QIN ,QDIS ,PONDING1 ,PONDING2, & - QSNBOT & + QSNBOT ,QTLDRN & #ifdef WRF_HYDRO - ,sfcheadrt & + ,sfcheadrt ,WATBLED & #endif ) !out ! ---------------------------------------------------------------------- @@ -5946,6 +5994,9 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & REAL, INTENT(INOUT) :: SMCWTD !soil water content between bottom of the soil and water table [m3/m3] REAL, INTENT(INOUT) :: DEEPRECH !recharge to or from the water table when deep [m] REAL, INTENT(INOUT) :: RECH !recharge to or from the water table when shallow [m] (diagnostic) + REAL, INTENT(INOUT) :: QTLDRN !tile drainage (mm/s) + REAL, INTENT(IN) :: DX + REAL, INTENT(IN) :: TDFRACMP !tile drain fraction map ! output REAL, INTENT(OUT) :: CMC !intercepted water per ground area (mm) @@ -5990,7 +6041,7 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & REAL, PARAMETER :: WSLMAX = 5000. !maximum lake water storage (mm) #ifdef WRF_HYDRO - REAL , INTENT(INOUT) :: sfcheadrt + REAL , INTENT(INOUT) :: sfcheadrt, WATBLED #endif ! ---------------------------------------------------------------------- @@ -6088,10 +6139,15 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & WSLAKE = WSLAKE + (QINSUR-QSEVA)*1000.*DT -RUNSRF*DT !mm ELSE ! soil CALL SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in - QINSUR ,QSEVA ,ETRANI ,SICE ,ILOC , JLOC , & !in - SH2O ,SMC ,ZWT ,VEGTYP , & !inout - SMCWTD, DEEPRECH , & !inout - RUNSRF ,QDRAIN ,RUNSUB ,WCND ,FCRMAX ) !out + QINSUR ,QSEVA ,ETRANI ,SICE ,ILOC , JLOC , & !in + TDFRACMP,DX , & !in + SH2O ,SMC ,ZWT ,VEGTYP , & !inout + SMCWTD, DEEPRECH, & !inout + RUNSRF ,QDRAIN ,RUNSUB ,WCND ,FCRMAX, QTLDRN & !out +#ifdef WRF_HYDRO + ,WATBLED & !in for tile drainage +#endif + ) IF(OPT_RUN == 1) THEN CALL GROUNDWATER (parameters,NSNOW ,NSOIL ,DT ,SICE ,ZSOIL , & !in @@ -6101,7 +6157,8 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & RUNSUB = QDIS !mm/s END IF - IF(OPT_RUN == 3 .or. OPT_RUN == 4) THEN + IF(OPT_RUN == 3 .or. OPT_RUN == 4 .or. & + OPT_RUN == 6 .or. OPT_RUN == 7 .or. OPT_RUN == 8) THEN RUNSUB = RUNSUB + QDRAIN !mm/s END IF @@ -7087,11 +7144,16 @@ END SUBROUTINE SNOWH2O !== begin soilwater ================================================================================ - SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in + SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in QINSUR ,QSEVA ,ETRANI ,SICE ,ILOC , JLOC, & !in - SH2O ,SMC ,ZWT ,VEGTYP ,& !inout - SMCWTD, DEEPRECH ,& !inout - RUNSRF ,QDRAIN ,RUNSUB ,WCND ,FCRMAX ) !out + TDFRACMP, DX , & !in + SH2O ,SMC ,ZWT ,VEGTYP , & !inout + SMCWTD, DEEPRECH, & !inout + RUNSRF ,QDRAIN ,RUNSUB ,WCND ,FCRMAX, QTLDRN & !out +#ifdef WRF_HYDRO + ,WATBLED & !in for tile drainage +#endif + ) ! ---------------------------------------------------------------------- ! calculate surface runoff and soil moisture. @@ -7111,8 +7173,9 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL !depth of soil layer-bottom [m] REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ETRANI !evapotranspiration from soil layers [mm/s] REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO !snow/soil layer depth [m] - REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE !soil ice content [m3/m3] - + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE !soil ice content [m3/m3] + REAL, INTENT(IN) :: DX + REAL, INTENT(IN) :: TDFRACMP! tile drainage map(fraction) INTEGER, INTENT(IN) :: VEGTYP ! input & output @@ -7121,6 +7184,10 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in REAL, INTENT(INOUT) :: ZWT !water table depth [m] REAL, INTENT(INOUT) :: SMCWTD !soil moisture between bottom of the soil and the water table [m3/m3] REAL , INTENT(INOUT) :: DEEPRECH + REAL , INTENT(INOUT) :: QTLDRN ! tile drainage (mm/s) +#ifdef WRF_HYDRO + REAL , INTENT(INOUT) :: WATBLED!in for tile drainage +#endif ! output REAL, INTENT(OUT) :: QDRAIN !soil-bottom free drainage [mm/s] @@ -7159,6 +7226,7 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in INTEGER :: NITER !iteration times soil moisture (-) REAL :: SMCTOT !2-m averaged soil moisture (m3/m3) REAL :: DZTOT !2-m soil depth (m) + REAL :: FACC !accumulated infiltration rate (m/s) REAL, PARAMETER :: A = 4.0 ! ---------------------------------------------------------------------- RUNSRF = 0.0 @@ -7256,6 +7324,19 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in END IF END IF + IF (OPT_RUN == 6) THEN + CALL COMPUTE_VIC_SURFRUNOFF(parameters,DT,NSOIL,SMC,ZSOIL,QINSUR,FSAT,RUNSRF,PDDUM) + END IF + + IF (OPT_RUN == 7) THEN + CALL COMPUTE_XAJ_SURFRUNOFF(parameters,DT,FCR,NSOIL,SMC,ZSOIL,QINSUR,RUNSRF,PDDUM) + END IF + + IF(OPT_RUN == 8)THEN + FACC = 1E-06 + CALL DYNAMIC_VIC(parameters,DT,SMC,SH2O,SICE,SICEMAX,NSOIL,ZSOIL,QINSUR,FACC,PDDUM,RUNSRF) + END IF + ! determine iteration times and finer time step NITER = 1 @@ -7270,9 +7351,10 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in DTFINE = DT / NITER ! solve soil moisture - + FACC = 1E-06 QDRAIN_SAVE = 0.0 RUNSRF_SAVE = 0.0 + DO ITER = 1, NITER IF(QINSUR > 0. .and. OPT_RUN == 3) THEN CALL INFIL (parameters,NSOIL ,DTFINE ,ZSOIL ,SH2O ,SICE , & !in @@ -7280,6 +7362,21 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in PDDUM ,RUNSRF ) !out END IF + IF(QINSUR > 0. .and. OPT_RUN == 6) THEN + CALL COMPUTE_VIC_SURFRUNOFF(parameters,DTFINE,NSOIL,SMC,ZSOIL,QINSUR,& !in + FSAT,RUNSRF,PDDUM) !out + END IF + + IF (QINSUR > 0. .AND. OPT_RUN == 7) THEN + CALL COMPUTE_XAJ_SURFRUNOFF(parameters,DTFINE,FCR,NSOIL,SMC,ZSOIL,QINSUR,& ! in + RUNSRF,PDDUM) ! out + END IF + + IF(QINSUR > 0. .and. OPT_RUN == 8) THEN + CALL DYNAMIC_VIC(parameters,DTFINE,SMC,SH2O,SICE,SICEMAX,NSOIL,& + ZSOIL,QINSUR,FACC,PDDUM,RUNSRF) + END IF + CALL SRT (parameters,NSOIL ,ZSOIL ,DTFINE ,PDDUM ,ETRANI , & !in QSEVA ,SH2O ,SMC ,ZWT ,FCR , & !in SICEMAX,FCRMAX ,ILOC ,JLOC ,SMCWTD , & !in @@ -7302,6 +7399,21 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in RUNSRF = RUNSRF * 1000. + RSAT * 1000./DT ! m/s -> mm/s QDRAIN = QDRAIN * 1000. +! Calling tile drainage ! pvk +! print*, "OPT_TDRN=",OPT_TDRN + IF (OPT_TDRN == 1 .AND. TDFRACMP .GT. 0.3 .AND. OPT_RUN == 3) THEN + print*, "simple tile drain scheme is on" + CALL TILE_DRAIN (parameters,NSOIL,SH2O,SMC,SICE,ZSOIL,QTLDRN,DT) + ELSE IF (OPT_TDRN == 2 .AND. TDFRACMP .GT. 0.1 .AND. OPT_RUN == 3) THEN + print*, "Hooghoudt tile drain scheme is on" + CALL TILE_HOOGHOUDT (parameters,WCND,NSOIL,NSNOW,SH2O,SMC,SICE,& + ZSOIL,DZSNSO,DT,DX,QTLDRN,ZWT & +#ifdef WRF_HYDRO + ,WATBLED & +#endif + ) + END IF + !WRF_HYDRO_DJG... !yw INFXSRT = RUNSRF * DT !mm/s -> mm @@ -7603,7 +7715,7 @@ SUBROUTINE SRT (parameters,NSOIL ,ZSOIL ,DT ,PDDUM ,ETRANI , & !in IF(OPT_RUN == 1 .or. OPT_RUN == 2) THEN QDRAIN = 0. END IF - IF(OPT_RUN == 3) THEN + IF(OPT_RUN == 3 .OR. OPT_RUN == 6 .OR. OPT_RUN == 7 .OR. OPT_RUN == 8) THEN QDRAIN = parameters%SLOPE*WCND(K) END IF IF(OPT_RUN == 4) THEN @@ -7619,7 +7731,8 @@ SUBROUTINE SRT (parameters,NSOIL ,ZSOIL ,DT ,PDDUM ,ETRANI , & !in ENDIF DSMDZ(K) = 2.0 * (SMX(K) - SMXBOT) / TEMP1 QDRAIN = WDF(K ) * DSMDZ(K ) + WCND(K ) - END IF + END IF + WFLUX(K) = -(WDF(K-1)*DSMDZ(K-1))-WCND(K-1)+ETRANI(K) + QDRAIN END IF END DO @@ -7771,6 +7884,1195 @@ SUBROUTINE SSTEP (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in END SUBROUTINE SSTEP +!==begin VIC subroutines ========================================================================== + + SUBROUTINE COMPUTE_VIC_SURFRUNOFF(parameters,DT,NSOIL,SMC,ZSOIL,QINSUR,ASAT,RUNSRF,PDDUM) + +! ---------------------------------------------------------------------- +! Calculate the saturated area and runoff based on VIC runoff scheme. +! This scheme adopted from VIC model +! Author: Prasanth Valayamkunnath +! ---------------------------------------------------------------------- + IMPLICIT NONE +! ---------------------------------------------------------------------- + +! Inputs + TYPE (noahmp_parameters), INTENT(IN) :: parameters + INTEGER, INTENT(IN) :: NSOIL + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL + REAL , INTENT(IN) :: QINSUR + REAL , INTENT(IN) :: DT + REAL , INTENT(OUT) :: ASAT +! Output + REAL , INTENT(INOUT):: RUNSRF + REAL , INTENT(INOUT):: PDDUM +!------------------------------------------------------------------------ +!local + REAL :: EX, I_0, I_MAX, BASIS, TOP_MOIST, TOP_MAX_MOIST + INTEGER :: IZ + +! Initialize Variables + EX = 0.0 + ASAT = 0.0 + I_MAX = 0.0 + I_0 = 0.0 + BASIS = 0.0 + TOP_MOIST = 0.0 + TOP_MAX_MOIST = 0.0 + RUNSRF = 0.0 + + DO IZ=1,NSOIL-2 + TOP_MOIST = TOP_MOIST + (SMC(IZ)*-1*ZSOIL(IZ)) ! m + TOP_MAX_MOIST = TOP_MAX_MOIST + (parameters%SMCMAX(IZ)*-1*ZSOIL(IZ)) ! m + END DO + + ! Saturated area from soil moisture + EX = parameters%BVIC/(1+parameters%BVIC) + ASAT = 1.0 - (( 1.0 - (TOP_MOIST/TOP_MAX_MOIST))**EX) ! + ASAT = MAX(0.0, ASAT) + ASAT = MIN(1.0, ASAT) + + ! Infiltration for the previous time-step soil moisture based on ASAT + I_MAX = (1.0 + parameters%BVIC)*TOP_MAX_MOIST ! m + I_0 = I_MAX*(1.0 - (1.0 - ASAT)**(1.0/parameters%BVIC)) !m + + ! Solve for surface runoff + IF(QINSUR .EQ. 0.0) THEN + RUNSRF = 0.0 + ELSE IF(I_MAX .EQ. 0.0) THEN + RUNSRF = QINSUR*DT + ELSE IF( (I_0 + (QINSUR*DT)) .GT. I_MAX ) THEN + RUNSRF = (QINSUR*DT) - TOP_MAX_MOIST + TOP_MOIST + ELSE + BASIS = 1.0 - ((I_0 + (QINSUR*DT))/I_MAX) + RUNSRF = (QINSUR*DT) - TOP_MAX_MOIST + TOP_MOIST + & + TOP_MAX_MOIST*(basis**(1.0+parameters%BVIC)) + END IF + + RUNSRF = RUNSRF/(DT) ! m/s + IF (RUNSRF .LT. 0.0) RUNSRF = 0.0 + IF (RUNSRF .GT. QINSUR)RUNSRF = QINSUR + + PDDUM = QINSUR - RUNSRF ! m/s + + END SUBROUTINE COMPUTE_VIC_SURFRUNOFF + +! End VIC subroutines + +!== begin xinanjiag================================================================================= + + SUBROUTINE COMPUTE_XAJ_SURFRUNOFF(parameters,DT,FCR,NSOIL,SMC,ZSOIL,QINSUR,RUNSRF,PDDUM) + +! ---------------------------------------------------------------------- +! Calculate the saturated area and runoff based on Xinanjiag runoff scheme. +! Reference: Knoben, W. J., Freer, J. E., Fowler, K. J., Peel, M. C., & Woods, R. A. (2019). +! Modular Assessment of Rainfall-Runoff Models Toolbox (MARRMoT) v1. 2: +! an open-source, extendable framework providing implementations of 46 conceptual +! hydrologic models as continuous state-space formulations. +! ---------------------------------------------------------------------- +! Author: Prasanth Valayamkunnath +! Date: August 03, 2020 +! ---------------------------------------------------------------------- + + IMPLICIT NONE +! ---------------------------------------------------------------------- +! Inputs + TYPE (noahmp_parameters), INTENT(IN) :: parameters + INTEGER, INTENT(IN) :: NSOIL + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: FCR !fraction of imperviousness (-) = IMP + REAL , INTENT(IN) :: QINSUR + REAL , INTENT(IN) :: DT +! Output + REAL , INTENT(INOUT):: RUNSRF + REAL , INTENT(INOUT):: PDDUM +! local + REAL :: WM,WM_MAX,SM,SM_MAX,IRUNOFF,PRUNOFF + INTEGER :: IZ +!------------------------------------------------------------------------ + +!initialize + WM = 0.0 + WM_MAX = 0.0 + SM = 0.0 + SM_MAX = 0.0 + IRUNOFF = 0.0 + PRUNOFF = 0.0 + RUNSRF = 0.0 + + DO IZ=1,NSOIL-2 + IF ((SMC(IZ)-parameters%SMCREF(IZ)) .GT. 0.) THEN ! soil moisture greater than field capacity + SM = SM + (SMC(IZ) - parameters%SMCREF(IZ) )*-1*ZSOIL(IZ) !m + WM = WM + (parameters%SMCREF(IZ)*-1*ZSOIL(IZ)) !m + ELSE + WM = WM + (SMC(IZ)*-1*ZSOIL(IZ)) + END IF + WM_MAX = WM_MAX + (parameters%SMCREF(IZ)*-1*ZSOIL(IZ)) + SM_MAX = SM_MAX + (parameters%SMCMAX(IZ) - parameters%SMCREF(IZ))*-1*ZSOIL(IZ) + END DO + WM = MIN(WM,WM_MAX) ! tension water (m) + SM = MIN(SM,SM_MAX) ! free water (m) + +! impervious surface runoff R_IMP + IRUNOFF = FCR(1)*QINSUR*DT + +! solve pervious surface runoff (m) based on Eq. (310) + IF ((WM/WM_MAX) .LE. (0.5-parameters%AXAJ))THEN + PRUNOFF = (1-FCR(1))*QINSUR*DT*((0.5-parameters%AXAJ)**(1-parameters%BXAJ))*((WM/WM_MAX)**parameters%BXAJ) + ELSE + PRUNOFF = (1-FCR(1))*QINSUR*DT*(1-(((0.5+parameters%AXAJ)**(1-parameters%BXAJ))*((1-(WM/WM_MAX))**parameters%BXAJ))) + END IF + +! estimate surface runoff based on Eq. (313) + IF(QINSUR .EQ. 0.0) THEN + RUNSRF = 0.0 + ELSE + RUNSRF = PRUNOFF*(1-((1-(SM/SM_MAX))**parameters%XXAJ))+IRUNOFF + END IF + RUNSRF = RUNSRF/DT !m/s + RUNSRF = MAX(0.0, RUNSRF) + RUNSRF = MIN(QINSUR, RUNSRF) + PDDUM = QINSUR - RUNSRF + + END SUBROUTINE COMPUTE_XAJ_SURFRUNOFF + +!== end xinanjiag ================================================================================== + +!== begin dynamic VIC ============================================================================== + + SUBROUTINE DYNAMIC_VIC(parameters,DT,SMC,SH2O,SICE,SICEMAX,NSOIL,ZSOIL,QINSUR,FACC,PDDUM,RUNSRF) + +! -------------------------------------------------------------------------------- +! compute inflitration rate at soil surface and estimate surface runoff based on +! Liang, X., & Xie, Z. (2001). A new surface runoff parameterization with subgrid-scale +! soil heterogeneity for land surface models. Advances in Water Resources, 24(9-10), 1173-1193. +! Author: Prasanth Valayamkunnath +! Date : August 3, 2020 +! -------------------------------------------------------------------------------- + + IMPLICIT NONE +! -------------------------------------------------------------------------------- +! inputs + type (noahmp_parameters), intent(in) :: parameters + INTEGER, INTENT(IN) :: NSOIL !no. of soil layers + REAL, INTENT(IN) :: DT !time step (sec) + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL !depth of soil layer-bottom [m] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2O !soil liquid water content [m3/m3] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE !soil ice content [m3/m3] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC !soil moisture content [m3/m3] + REAL, INTENT(IN) :: QINSUR !water input on soil surface [m/s] + REAL, INTENT(IN) :: SICEMAX !maximum soil ice content (m3/m3) +! inouts + REAL, INTENT(INOUT) :: FACC !accumulated infiltration (m) +! outputs + REAL, INTENT(OUT) :: RUNSRF !surface runoff [mm/s] + REAL, INTENT(OUT) :: PDDUM !infiltration rate at surface +! locals + REAL :: BB !B parameter for infiltration scaling curve + REAL :: TOP_MOIST !actual water depth in top layers (m) + REAL :: TOP_MAX_MOIST !water depth in top layers (m) + REAL :: DP !water input on soil surface (m) + REAL :: I_0 !initial water depth (m) + REAL :: I_MAX !maximum water depth (m) + REAL :: FSUR !surface infiltration rate (m/s) + REAL :: FMAX !maximum infiltration rate (m/s) + REAL :: RUNOFFSAT !saturation excess runoff (m/s) + REAL :: RUNOFFINF !infiltration excess runoff (m/s) + REAL :: INFILTRTN !infiltration (m/s) + REAL :: TEMPR1 !temporary saturation excess runoff (m/s) + REAL :: TEMPR2 !temporary infiltration excess runoff (m/s) + REAL :: R1 !saturation excess runoff (m/s) + REAL :: R2 !infiltration excess runoff (m/s) + REAL :: YD !initial depth Y (m) + REAL :: YD_OLD !initial depth Y (m) + REAL :: YD0 !initial depth Y (m) + REAL :: TEMP1, ERROR + INTEGER :: IZ, IZMAX, INFLMAX +!--------------------------------------------------------------------------------- + + TOP_MOIST = 0.0 + TOP_MAX_MOIST = 0.0 + BB = 1.0 + DP = 0.0 + I_MAX = 0.0 + I_0 = 0.0 + RUNOFFSAT = 0.0 + RUNOFFINF = 0.0 + INFILTRTN = 0.0 + RUNSRF = 0.0 + IZMAX = 20 + ERROR = 1.388889E-07*DT ! 0.5 mm per hour time step + BB = parameters%BBVIC + + DO IZ=1,NSOIL-2 + TOP_MOIST = TOP_MOIST + (SMC(IZ)*-1*ZSOIL(IZ)) ! actual moisture in top layers, [m] + TOP_MAX_MOIST = TOP_MAX_MOIST + (parameters%SMCMAX(IZ)*-1*ZSOIL(IZ)) ! maximum moisture in top layers, [m] + END DO + + IF(TOP_MOIST .GT. TOP_MAX_MOIST) TOP_MOIST = TOP_MAX_MOIST + DP = QINSUR * DT ! precipitation depth, [m] + I_MAX = TOP_MAX_MOIST * (parameters%BDVIC+1.0) ! maximum infiltration capacity, im, [m], Eq. 14 + I_0 = I_MAX * (1-(1-(TOP_MOIST/TOP_MAX_MOIST)**(1.0/(1.0+parameters%BDVIC)))) ! infiltration capacity, i [m] in the Eq. 1 + ! I_MAX = CAP_minf ; I_0 = A + INFLMAX = 0 + + IF (OPT_INFDV .EQ. 1) THEN + CALL PHILIP_INFIL(parameters,NSOIL,SMC,SICE,QINSUR,DT,FACC,FSUR,INFLMAX) + ELSE IF (OPT_INFDV .EQ. 2) THEN + CALL GREEN_AMPT_INFIL(parameters,NSOIL,ZSOIL,SMC,SICE,QINSUR,FACC,FSUR,INFLMAX) + ELSE IF (OPT_INFDV .EQ. 3) THEN + CALL SMITH_PARLANGE_INFIL(parameters,NSOIL,ZSOIL,SMC,SICE,QINSUR,FACC,FSUR,INFLMAX) + END IF + + ! I_MM = FSUR; I_M = FMAX + FMAX = (BB+1.0)*FSUR + IF(DP .LE. 0.0) THEN + RUNOFFSAT = 0.0 + RUNOFFINF = 0.0 + INFILTRTN = 0.0 + GOTO 2001 + ELSE + IF((TOP_MOIST .GE. TOP_MAX_MOIST) .AND. (I_0 .GE. I_MAX)) THEN + TOP_MOIST = TOP_MAX_MOIST + I_0 = I_MAX + RUNOFFSAT = DP + RUNOFFINF = 0.0 + INFILTRTN = 0.0 + GOTO 2001 + ELSE + I_0 = I_MAX * (1-(1-(TOP_MOIST/TOP_MAX_MOIST)**(1.0/(1.0+parameters%BDVIC)))) + IF((DP+I_0) .GT. I_MAX)THEN + IF((FMAX*DT) .GE. DP) THEN + YD = I_MAX - I_0 + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + TEMP1 = I_MAX-I_0-TEMPR1-((FSUR*DT) * (1 - (1-((DP-TEMPR1)/(FMAX*DT))**(BB+1.0)))) + IF(TEMP1 .LE. 0.0) THEN + YD = I_MAX - I_0 + INFILTRTN = TOP_MAX_MOIST - TOP_MOIST + RUNOFFSAT = DP - INFILTRTN + RUNOFFINF = 0.0 + TOP_MOIST = TOP_MAX_MOIST + I_0 = I_MAX + GOTO 2001 + ELSE + YD = 0.0 + + DO IZ = 1,IZMAX ! loop : IITERATION1 + YD_OLD = YD + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + YD = TEMPR1 + ((FSUR*DT) * (1 - (1-((DP-TEMPR1)/(FMAX*DT))**(BB+1.0)))) + IF ((ABS(YD-YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + GOTO 1003 + END IF + END DO + END IF + ELSE + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + IF((TEMPR1+(FMAX*DT)) .LE. DP) THEN + IF((I_MAX-I_0-TEMPR1-(FMAX*DT)) .LE. 0.0)THEN + YD = I_MAX - I_0 + INFILTRTN = TOP_MAX_MOIST - TOP_MOIST + RUNOFFSAT = DP - INFILTRTN + RUNOFFINF = 0.0 + TOP_MOIST = TOP_MAX_MOIST + I_0 = I_MAX + GOTO 2001 + ELSE + YD = 0.0 + + DO IZ = 1,IZMAX ! loop : IITERATION2 + YD_OLD = YD + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + YD = TEMPR1 + (FSUR*DT) + IF ((ABS(YD-YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + GOTO 1003 + END IF + END DO + END IF + + ELSE + + YD = DP/2.0 + DO IZ = 1,IZMAX ! loop : IITERATION30 + YD_OLD = YD + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + YD = YD - TEMPR1 - (FSUR*DT) + DP + IF (YD .LE. 0.0) YD = 0.0 + IF (YD .GE. DP) YD = DP + IF ((ABS(YD-YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + YD0 = YD + EXIT + END IF + END DO + DO IZ = 1,IZMAX ! loop : IITERATION3 + YD_OLD = YD + TEMPR1 = 0.0 + TEMPR2 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + CALL RR2(YD,YD0,TEMPR1,FMAX,FSUR,DT,DP,BB,TEMPR2) + YD = DP - TEMPR2 + IF ((ABS(YD-YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + GOTO 1003 + END IF + END DO +1003 IF(YD .LE. 0.0) YD = 0.0 + IF(YD .GE. DP) YD = DP + CALL RR1(parameters,I_0,I_MAX,YD,R1) + RUNOFFSAT = R1 + RUNOFFINF = DP - YD + INFILTRTN = YD - RUNOFFSAT + TOP_MOIST = TOP_MOIST + INFILTRTN + YD = I_0+YD + IF (TOP_MOIST .LE. 0.0) TOP_MOIST=0.0 + IF (TOP_MOIST .GE. TOP_MAX_MOIST) TOP_MOIST = TOP_MAX_MOIST + I_0 = I_MAX * (1-(1-(TOP_MOIST/TOP_MAX_MOIST)**(1.0/(1.0+parameters%BDVIC)))) + GOTO 2001 + END IF + END IF + + ELSE + IF((FMAX*DT) .GE. DP) THEN + YD = DP/2.0 + DO IZ = 1,IZMAX ! ITERATION1 + YD_OLD = YD + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + YD = TEMPR1 + ((FSUR*DT) * (1 - (1-((DP-TEMPR1)/(FMAX*DT))**(BB+1.0)))) + !print*,'YD=',YD,'YD_OLD=',YD_OLD,'ERROR=',ABS(YD - YD_OLD),'IZ=',IZ,'ERROR2=',ERROR,'DT=',DT + IF ((ABS(YD - YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + GOTO 1004 + END IF + END DO + ELSE + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + IF((TEMPR1+(FMAX*DT)) .LE. DP)THEN + YD = DP/2.0 + DO IZ = 1,IZMAX ! ITERATION2 + YD_OLD = YD + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + YD = TEMPR1+(FSUR*DT) + IF((ABS(YD - YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + GOTO 1004 + END IF + END DO + ELSE + YD = 0.0 + DO IZ = 1,IZMAX ! ITERATION30 + YD_OLD = YD + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + YD = (DP - (FMAX*DT)) + YD - TEMPR1 + IF(YD .LE. 0.0) YD = 0.0 + IF(YD .GE. DP) YD = DP + TEMPR1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + IF ((ABS(TEMPR1+(FMAX*DT)-DP) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + YD0 = YD + EXIT + END IF + END DO + DO IZ = 1,IZMAX ! ITERATION3 + YD_OLD = YD + TEMPR1 = 0.0 + TEMPR2 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,TEMPR1) + CALL RR2(YD,YD0,TEMPR1,FMAX,FSUR,DT,DP,BB,TEMPR2) + YD = DP - TEMPR2 + IF ((ABS(YD-YD_OLD) .LE. ERROR) .OR. (IZ .EQ. IZMAX)) THEN + GOTO 1004 + END IF + END DO + END IF + END IF +1004 IF(YD .LE. 0.0) YD = 0.0 + IF(YD .GE. DP) YD = DP + R1 = 0.0 + CALL RR1(parameters,I_0,I_MAX,YD,R1) + RUNOFFSAT = R1 + RUNOFFINF = DP - YD + INFILTRTN = YD - RUNOFFSAT + TOP_MOIST = TOP_MOIST + INFILTRTN + IF (TOP_MOIST .LE. 0.0) TOP_MOIST=0.0 + IF (TOP_MOIST .GE. TOP_MAX_MOIST) TOP_MOIST = TOP_MAX_MOIST + I_0 = I_MAX * (1-(1-(TOP_MOIST/TOP_MAX_MOIST)**(1.0/(1.0+parameters%BDVIC)))) + END IF + END IF + END IF + +2001 RUNSRF = (RUNOFFSAT + RUNOFFINF)/DT + RUNSRF = MIN(RUNSRF,QINSUR) + RUNSRF = MAX(RUNSRF,0.0) + PDDUM = QINSUR - RUNSRF + + END SUBROUTINE DYNAMIC_VIC + +! --------------------------- Runoff subroutines for dynamic VIC ---------------------------- + + SUBROUTINE RR1 (parameters,I_0,I_MAX,YD,R1) +!--------------------------------------------------------------------------------------------- +! This subroutine estimate saturation excess runoff, R1 +! Author: Prasanth Valayamkunnath +!--------------------------------------------------------------------------------------------- + IMPLICIT NONE +! -------------------------------------------------------------------------------------------- + type (noahmp_parameters), intent(in) :: parameters + REAL, INTENT(IN) :: I_0,I_MAX,YD + REAL, INTENT(OUT):: R1 + REAL :: TDEPTH +!------------------------------------------------------ + + TDEPTH = I_0 + YD + IF(TDEPTH .GT. I_MAX) TDEPTH = I_MAX + + !Saturation excess runoff , Eq 5. + R1 = YD - ( (I_MAX/(parameters%BDVIC+1.0)) * ( ((1 - (I_0/I_MAX))**(parameters%BDVIC+1.0)) & + - ((1 - (TDEPTH/I_MAX))**(parameters%BDVIC+1.0)))) + + IF (R1 .LT. 0.0) R1 = 0.0 + + END SUBROUTINE RR1 + +!--------------------------------------------------------------------------------------------- + SUBROUTINE RR2 (YD,Y0,R1,FMAX,FSUR,DT,DP,BB,R2) +!--------------------------------------------------------------------------------------------- +! This subroutine estimate infiltration excess runoff, R1 +! Author: Prasanth Valayamkunnath +!--------------------------------------------------------------------------------------------- + IMPLICIT NONE +! -------------------------------------------------------------------------------------------- + REAL, INTENT(IN) :: YD,Y0,R1,FMAX,FSUR,DT,DP,BB + REAL, INTENT(OUT):: R2 +!------------------------------------------------------ + + IF(YD .GE. Y0)THEN + R2 = DP - R1 - (FMAX*DT* (1 - ((1 - (DP-R1)/(FMAX*DT))**(BB+1.0)))) + ELSE + R2 = DP - R1 - (FMAX*DT) + END IF + + IF(R2 .LT. 0.0) R2 =0.0 + +END SUBROUTINE RR2 + +!== end dynamic VIC ================================================================================ + +!== begin smith-parlange infiltration =============================================================== + + SUBROUTINE SMITH_PARLANGE_INFIL(parameters,NSOIL,ZSOIL,SMC,SICE,QINSUR,FACC,FSUR,INFLMAX) + +!--------------------------------------------------------------------------------------------- +! This function estimate infiltration rate based on Smith-Parlange equation. We use its three +! parameter version of the equation (Eq. 6.25) from Smith, R.E. (2002) Infiltration Theory for +! Hydrologic Applications, Water Resources Monograph 15, AGU. +! Author: Prasanth Valayamkunnath +!--------------------------------------------------------------------------------------------- + IMPLICIT NONE +! -------------------------------------------------------------------------------------------- + type (noahmp_parameters), intent(in) :: parameters + INTEGER, INTENT(IN) :: NSOIL !no of soil layers (4) + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL !depth of soil layer-bottom [m] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC !soil moisture content [m3/m3] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE !soil ice content [m3/m3] + REAL, INTENT(IN) :: QINSUR !water input on soil surface [m/s] + INTEGER, INTENT(IN) :: INFLMAX!check for maximum infiltration at SMCWLT +! in outs + REAL, INTENT(INOUT) :: FACC !accumulated infiltration rate (m/s) +! outputs + REAL, INTENT(OUT) :: FSUR !surface infiltration rate (m/s) +! local variables + REAL :: WDF ! soil water diffusivity + REAL :: WCND ! soil water conductivity[m/s] + REAL :: GAM ! smith-parlang weighing parameter[-] + REAL :: JJ ! dummy variable + INTEGER :: ISOIL +!--------------------------------------------------------------------------------- + + ! smith-parlang weighing parameter, GAMMA + GAM = 0.82 + ISOIL = 1 + ! check whether we are estimating infiltration for current SMC or SMCWLT + IF (INFLMAX .EQ. 1)THEN ! not active for now as the maximum infiltration is estimated based on table values + + ! estimate initial soil hydraulic conductivty (Ki in the equation), WCND (m/s) + CALL WDFCND2 (parameters,WDF,WCND,parameters%SMCWLT(ISOIL),0.0,ISOIL) + + ! Maximum infiltrability based on the Eq. 6.25. (m/s) + JJ = parameters%GDVIC * (parameters%SMCMAX(ISOIL) - parameters%SMCWLT(ISOIL)) * -1 * ZSOIL(ISOIL) + FSUR = parameters%DKSAT(ISOIL) + (GAM * (parameters%DKSAT(ISOIL) - WCND) / (EXP(GAM * 1E-05 / JJ) -1)) + + ! infiltration rate at surface + IF(parameters%DKSAT(ISOIL) .LT. QINSUR)THEN + FSUR = MIN(QINSUR,FSUR) + ELSE + FSUR = QINSUR + END IF + IF(FSUR .LT. 0.0) FSUR = 0.0 + + ELSE + + ! estimate initial soil hydraulic conductivty (Ki in the equation), WCND (m/s) + CALL WDFCND2 (parameters,WDF,WCND,SMC(ISOIL),SICE(ISOIL),ISOIL) + + ! Maximum infiltrability based on the Eq. 6.25. (m/s) + JJ = parameters%GDVIC * (parameters%SMCMAX(ISOIL) - SMC(ISOIL)) * -1 * ZSOIL(ISOIL) + FSUR = parameters%DKSAT(ISOIL) + (GAM * (parameters%DKSAT(ISOIL) - WCND) / (EXP(GAM * FACC / JJ) -1)) + + ! infiltration rate at surface + IF(parameters%DKSAT(ISOIL) .LT. QINSUR)THEN + FSUR = MIN(QINSUR,FSUR) + ELSE + FSUR = QINSUR + END IF + + ! accumulated infiltration function + FACC = FACC + FSUR + + END IF + + END SUBROUTINE SMITH_PARLANGE_INFIL + +!== end smith-parlang infiltration ================================================================= + +!== begin Green_Ampt infiltration ================================================================== + + SUBROUTINE GREEN_AMPT_INFIL(parameters,NSOIL,ZSOIL,SMC,SICE,QINSUR,FACC,FSUR,INFLMAX) + +!------------------------------------------------------------------------------------------------- +! This function estimate infiltration rate based on Green-Ampt equation. We use its three +! parameter version of the smith-parlage equation (Eq. 6.25) from Smith, R.E. (2002) Infiltration Theory for +! Hydrologic Applications, Water Resources Monograph 15, AGU. Where gamma = 0, Eq 6.25 = Green-Ampt. +! Author: Prasanth Valayamkunnath +!------------------------------------------------------------------------------------------------- + IMPLICIT NONE +! ------------------------------------------------------------------------------------------------ + type (noahmp_parameters), intent(in) :: parameters + INTEGER, INTENT(IN) :: NSOIL !no of soil layers (4) + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL !depth of soil layer-bottom [m] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC !soil moisture content [m3/m3] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE !soil ice content [m3/m3] + REAL, INTENT(IN) :: QINSUR !water input on soil surface [m/s] + INTEGER, INTENT(IN) :: INFLMAX!check for maximum infiltration at SMCWLT +! in outs + REAL, INTENT(INOUT) :: FACC !accumulated infiltration rate (m/s) +! outputs + REAL, INTENT(OUT) :: FSUR !surface infiltration rate (m/s) +! local variables + REAL :: WDF ! soil water diffusivity + REAL :: WCND ! soil water conductivity[m/s] + REAL :: JJ ! dummy variable + INTEGER :: ISOIL +!--------------------------------------------------------------------------------- + + ISOIL = 1 + IF(INFLMAX .EQ. 1)THEN + + ! estimate initial soil hydraulic conductivty (Ki in the equation), WCND (m/s) + CALL WDFCND2 (parameters,WDF,WCND,parameters%SMCWLT(ISOIL),0.0,ISOIL) + + ! Maximum infiltrability based on the Eq. 6.25. (m/s) + JJ = parameters%GDVIC * (parameters%SMCMAX(ISOIL) - parameters%SMCWLT(ISOIL)) * -1 * ZSOIL(ISOIL) + FSUR = parameters%DKSAT(ISOIL) + ((JJ/1E-05) * (parameters%DKSAT(ISOIL) - WCND)) + + !maximum infiltration rate at surface + IF(FSUR .LT. 0.0) FSUR = 0.0 + + ELSE + + ! estimate initial soil hydraulic conductivty (Ki in the equation), WCND (m/s) + CALL WDFCND2 (parameters,WDF,WCND,SMC(ISOIL),SICE(ISOIL),ISOIL) + + ! Maximum infiltrability based on the Eq. 6.25. (m/s) + JJ = parameters%GDVIC * (parameters%SMCMAX(ISOIL) - SMC(ISOIL)) * -1 * ZSOIL(ISOIL) + FSUR = parameters%DKSAT(ISOIL) + ((JJ/FACC) * (parameters%DKSAT(ISOIL) - WCND)) + + ! infiltration rate at surface + IF(parameters%DKSAT(ISOIL) .LT. QINSUR)THEN + FSUR = MIN(QINSUR,FSUR) + ELSE + FSUR = QINSUR + END IF + ! accumulated infiltration function + FACC = FACC + FSUR + + END IF + + END SUBROUTINE GREEN_AMPT_INFIL + +!== end Green-Ampt infiltration ==================================================================== + +!== begin Philip's infiltration ==================================================================== + + SUBROUTINE PHILIP_INFIL(parameters,NSOIL,SMC,SICE,QINSUR,DT,FACC,FSUR,INFLMAX) + +!------------------------------------------------------------------------------------------------------- +! This function estimate infiltration rate based on Philip's two parameter equation (Eq. 2) presented in +! Valiantzas (2010). New linearized two-parameter infiltration equation for direct determination +! of conductivity and sorptivity, J. Hydrology. +! Author: Prasanth Valayamkunnath +!--------------------------------------------------------------------------------------------- + IMPLICIT NONE +! -------------------------------------------------------------------------------------------- + type (noahmp_parameters), intent(in) :: parameters + INTEGER, INTENT(IN) :: NSOIL !no of soil layers (4) + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC !soil moisture content [m3/m3] + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE !soil ice content [m3/m3] + REAL, INTENT(IN) :: QINSUR !water input on soil surface [m/s] + REAL, INTENT(IN) :: DT !time-step (sec) + INTEGER, INTENT(IN) :: INFLMAX!check for maximum infiltration at SMCWLT +! in outs + REAL, INTENT(INOUT) :: FACC !accumulated infiltration rate (m/s) +! outputs + REAL, INTENT(OUT) :: FSUR !surface infiltration rate (m/s) +! local variables + REAL :: WDF ! soil water diffusivity (m2/s) + REAL :: WCND ! soil water conductivity[m/s] + REAL :: SP ! sorptivity (LT^-1/2) + REAL :: AP ! intial hydraulic conductivity (m/s,L/T) + REAL :: FMAX ! Maximum infiltration (m/s) + INTEGER :: ISOIL +!--------------------------------------------------------------------------------- + + ISOIL = 1 + IF (INFLMAX .EQ. 1) THEN + + ! estimate initial soil hydraulic conductivty and diffusivity (Ki, D(theta) in the equation) + CALL WDFCND2 (parameters,WDF,WCND,parameters%SMCWLT(ISOIL),0.0,ISOIL) + + ! Sorptivity based on Eq. 10b from Kutílek, Miroslav, and Jana Valentová (1986) + ! Sorptivity approximations. Transport in Porous Media 1.1, 57-62. + SP = SQRT(2 * (parameters%SMCMAX(ISOIL) - parameters%SMCWLT(ISOIL)) * (parameters%DWSAT(ISOIL) - WDF)) + + ! Parameter A in Eq. 9 of Valiantzas (2010) is given by + AP = MIN(WCND, (2/3)*parameters%DKSAT(ISOIL)) + AP = MAX(AP, (1/3)*parameters%DKSAT(ISOIL)) + + ! Maximun infiltration rate, m + FSUR = (1/2)*SP*(DT**(-1/2))+AP ! m/s + IF(FSUR .LT. 0.0) FSUR = 0.0 + + ELSE + + ! estimate initial soil hydraulic conductivty and diffusivity (Ki, D(theta) in the equation) + CALL WDFCND2 (parameters,WDF,WCND,SMC(ISOIL),SICE(ISOIL),ISOIL) + + ! Sorptivity based on Eq. 10b from Kutílek, Miroslav, and Jana Valentová (1986) + ! Sorptivity approximations. Transport in Porous Media 1.1, 57-62. + SP = SQRT(2 * (parameters%SMCMAX(ISOIL) - SMC(ISOIL)) * (parameters%DWSAT(ISOIL) - WDF)) + + ! Parameter A in Eq. 9 of Valiantzas (2010) is given by + AP = MIN(WCND, (2/3)*parameters%DKSAT(ISOIL)) + AP = MAX(AP, (1/3)*parameters%DKSAT(ISOIL)) + + ! Maximun infiltration rate, m + FSUR = (1/2)*SP*(DT**(-1/2))+AP ! m/s + + ! infiltration rate at surface + IF(parameters%DKSAT(ISOIL) .LT. QINSUR)THEN + FSUR = MIN(QINSUR,FSUR) + ELSE + FSUR = QINSUR + END IF + ! accumulated infiltration function + FACC = FACC + FSUR + + END IF + + END SUBROUTINE PHILIP_INFIL + +!== end Phillips infiltration ====================================================================== + +!== begin tile_drain =============================================================================== + + SUBROUTINE TILE_DRAIN (parameters,NSOIL,SH2O,SMC,SICE,ZSOIL,QTLDRN,DT) + +! ---------------------------------------------------------------------- +! Calculate tile drainage discharge (mm) based on simple model ! pvk +! ---------------------------------------------------------------------- + IMPLICIT NONE +! ---------------------------------------------------------------------- +! inout + type (noahmp_parameters), intent(in) :: parameters + INTEGER,INTENT(IN) :: NSOIL + REAL, INTENT(IN) :: DT + REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O + REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC + REAL,INTENT(INOUT) :: QTLDRN +! input + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL +! local + INTEGER :: K + REAL :: TDRVOL ! temp variable for drainage volume (mm) + REAL,DIMENSION(1:NSOIL) :: OVRFC ! temp variable for volume of water above field capacity + REAL,DIMENSION(1:NSOIL) :: AVFC ! Available field capacity = FC - SICE (m3/m3) + REAL,DIMENSION(1:NSOIL) :: ZLAYER ! thickness of soil layer + REAL,DIMENSION(1:NSOIL) :: TDFRAC + REAL :: TDDC + REAL :: TDSUM +! ---------------------------------------------------------------------- + + TDRVOL = 0. + OVRFC = 0. + QTLDRN = 0. + ZLAYER = 0. + AVFC = 0. + TDSUM = 0. + TDFRAC = 0. + TDDC = parameters%TD_DC * DT/(24*3600) + + DO K = 1, NSOIL + IF (K == 1) THEN + ZLAYER(K) = -1 * ZSOIL(K) + ELSE + ZLAYER(K) = (ZSOIL(K-1)-ZSOIL(K)) + END IF + END DO + IF (parameters%DRAIN_LAYER_OPT == 0) THEN ! drainage from one specified layer in MPTABLE.TBL + !print*, "CASE = 1" + K = parameters%TD_DEPTH + AVFC(K) = parameters%SMCREF(K) - SICE (K) + OVRFC(K) = (SH2O(K) - (parameters%TDSMC_FAC*AVFC(K))) * ZLAYER(K) * 1000. ! mm + IF (OVRFC(K) > 0.) THEN + IF (OVRFC(K) > TDDC) OVRFC(K) = TDDC + TDRVOL = TDRVOL + OVRFC(K) + SH2O(K) = SH2O(K) - (OVRFC(K)/(ZLAYER(K) * 1000.)) + SMC(K) = SH2O(K) + SICE (K) + END IF + ELSE IF (parameters%DRAIN_LAYER_OPT == 1) THEN + !print*, "CASE = 2. Draining from layer 1 and 2" + DO K = 1, 2 + AVFC(K) = parameters%SMCREF(K) - SICE (K) + OVRFC(K) = (SH2O(K) - (parameters%TDSMC_FAC*AVFC(K))) * ZLAYER(K) * 1000. ! mm + IF(OVRFC(K) < 0.) OVRFC(K) = 0. + TDSUM = TDSUM + OVRFC(K) + END DO + DO K = 1, 2 + IF(OVRFC(K) .NE. 0.) THEN + TDFRAC(K) = OVRFC(K)/TDSUM + END IF + END DO + IF (TDSUM > 0.) THEN + IF (TDSUM > TDDC) TDSUM = TDDC + TDRVOL = TDRVOL + TDSUM + DO K = 1, 2 + OVRFC(K) = TDFRAC(K)*TDSUM + SH2O(K) = SH2O(K) - (OVRFC(K)/(ZLAYER(K) * 1000.)) + SMC(K) = SH2O(K) + SICE (K) + END DO + END IF + ELSE IF (parameters%DRAIN_LAYER_OPT == 2) THEN + !print*, "CASE = 3. Draining from layer 1 2 and 3" + DO K = 1, 3 + AVFC(K) = parameters%SMCREF(K) - SICE (K) + OVRFC(K) = (SH2O(K) - (parameters%TDSMC_FAC*AVFC(K))) * ZLAYER(K) * 1000. ! mm + IF(OVRFC(K) < 0.) OVRFC(K) = 0. + TDSUM = TDSUM + OVRFC(K) + END DO + DO K = 1, 3 + IF(OVRFC(K) .NE. 0.) THEN + TDFRAC(K) = OVRFC(K)/TDSUM + END IF + END DO + IF (TDSUM > 0.) THEN + IF (TDSUM > TDDC) TDSUM = TDDC + TDRVOL = TDRVOL + TDSUM + DO K = 1, 3 + OVRFC(K) = TDFRAC(K)*TDSUM + SH2O(K) = SH2O(K) - (OVRFC(K)/(ZLAYER(K) * 1000.)) + SMC(K) = SH2O(K) + SICE (K) + END DO + END IF + + ELSE IF (parameters%DRAIN_LAYER_OPT == 3) THEN + !print*, "CASE = 3. Draining from layer 2 and 3" + DO K = 2, 3 + AVFC(K) = parameters%SMCREF(K) - SICE (K) + OVRFC(K) = (SH2O(K) - (parameters%TDSMC_FAC*AVFC(K))) * ZLAYER(K) * 1000. ! mm + IF(OVRFC(K) < 0.) OVRFC(K) = 0. + TDSUM = TDSUM + OVRFC(K) + END DO + DO K = 2, 3 + IF(OVRFC(K) .NE. 0.) THEN + TDFRAC(K) = OVRFC(K)/TDSUM + END IF + END DO + IF (TDSUM > 0.) THEN + IF (TDSUM > TDDC) TDSUM = TDDC + TDRVOL = TDRVOL + TDSUM + DO K = 2, 3 + OVRFC(K) = TDFRAC(K)*TDSUM + SH2O(K) = SH2O(K) - (OVRFC(K)/(ZLAYER(K) * 1000.)) + SMC(K) = SH2O(K) + SICE (K) + END DO + END IF + ELSE IF (parameters%DRAIN_LAYER_OPT == 4) THEN + + !print*, "CASE = 4. Draining from layer 3 and 4" + DO K = 3, 4 + AVFC(K) = parameters%SMCREF(K) - SICE (K) + OVRFC(K) = (SH2O(K) - (parameters%TDSMC_FAC*AVFC(K))) * ZLAYER(K) * 1000. ! mm + IF(OVRFC(K) < 0.) OVRFC(K) = 0. + TDSUM = TDSUM + OVRFC(K) + END DO + DO K = 3, 4 + IF(OVRFC(K) .NE. 0.) THEN + TDFRAC(K) = OVRFC(K)/TDSUM + END IF + END DO + + + IF (TDSUM > 0.) THEN + IF (TDSUM > TDDC) TDSUM = TDDC + TDRVOL = TDRVOL + TDSUM + + DO K = 3, 4 + OVRFC(K) = TDFRAC(K)*TDSUM + SH2O(K) = SH2O(K) - (OVRFC(K)/(ZLAYER(K) * 1000.)) + SMC(K) = SH2O(K) + SICE (K) + END DO + END IF + + ELSE IF (parameters%DRAIN_LAYER_OPT == 5) THEN ! from all the four layers + + !print*, "CASE = 5 Draining from all four layers" + DO K = 1, 4 + AVFC(K) = parameters%SMCREF(K) - SICE (K) + OVRFC(K) = (SH2O(K) - (parameters%TDSMC_FAC*AVFC(K))) * ZLAYER(K) * 1000. ! mm + IF(OVRFC(K) < 0.) OVRFC(K) = 0. + TDSUM = TDSUM + OVRFC(K) + END DO + DO K = 1, 4 + IF(OVRFC(K) .NE. 0.) THEN + TDFRAC(K) = OVRFC(K)/TDSUM + END IF + END DO + + IF (TDSUM > 0.) THEN + IF (TDSUM > TDDC) TDSUM = TDDC + TDRVOL = TDRVOL + TDSUM + DO K = 1, 4 + OVRFC(K) = TDFRAC(K)*TDSUM + SH2O(K) = SH2O(K) - (OVRFC(K)/(ZLAYER(K) * 1000.)) + SMC(K) = SH2O(K) + SICE (K) + END DO + END IF + END IF + + QTLDRN = TDRVOL / DT + !print*,"QTLDRN =",QTLDRN + + END SUBROUTINE TILE_DRAIN + +!=================================================================================================== + + SUBROUTINE TILE_HOOGHOUDT (parameters,WCND,NSOIL,NSNOW,SH2O,SMC,SICE, & + ZSOIL,DZSNSO,DT,DX,QTLDRN,ZWT & +#ifdef WRF_HYDRO + ,WATBLED & +#endif + ) + +!------------------------------------------------------------------------------------------ +! calculate tile drainage discharge (mm) based on Hooghoudt's equation ! pvk +!------------------------------------------------------------------------------------------ + IMPLICIT NONE +!------------------------------------------------------------------------------------------ +! Definitions +!TD_SATZ = Thickness of saturated zone in layer considered ~ W +!TD_SPAC = Tile Drain SPACing ~ SDRAIN +!TD_HAIL = Height of water table in the drain Above Impermeable Layer (de) - +! HDRAIN +!TD_HEMD = Effective Height between water level in the drains to the water +! table MiDpoint ~ EM +!ZLAYER(K) = Thickness of layer K ~DZ(I) +!TD_DTWT = Depth To Water Table (cm) ~ TD_DTWT +!TD_FLUX = Drainge Flux cm/hr~ DFLUX +!TD_DEPTH = Effective Depth to impermeable layer from soil surface~ DEPTH +!TD_TTSZ = Total Thickness of Saturated Zone ~ DEEP +!KLAT = Lateral hydraulic conductivity, CONE +!KLATK = Lateral hydraulic conductivity of a specific layer ~ CONK +!DTOPL = Depth of top of the layer considered~ ABOVE +!TD_DC = Drainage Coefficient ~ DC +!TD_DDRN = Depth of drain ~ DDRAIN +!TD_ADEPTH = Actual depth to impermeable layer form surface. ADEPTH >= TD_DEPTH +!------------------------------------------------------------------------------------------ +! inout + type (noahmp_parameters), intent(in) :: parameters + INTEGER,INTENT(IN) :: NSOIL + INTEGER,INTENT(IN) :: NSNOW + REAL, INTENT(IN) :: DT + REAL, INTENT(IN) :: DX + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SICE + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: DZSNSO + REAL, DIMENSION(1:NSOIL), INTENT(IN) :: WCND + + REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O + REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC + REAL, INTENT(INOUT) :: QTLDRN ! tile drain discharge mm/s + REAL, INTENT(INOUT) :: ZWT ! water table depth + +! local + INTEGER :: K + REAL, DIMENSION(1:NSOIL) :: TD_SATZ ! thickness of saturated zone + REAL, DIMENSION(1:NSOIL) :: KLATK ! lateral hydraulic ocnductivity kth layer + REAL :: TD_TTSZ ! total satutation thickness + REAL :: TD_LQ ! lateral flow + REAL :: DTOPL ! depth to top of the layer + REAL :: XX + REAL :: YY + REAL :: KLAT ! average lateral hydruaic conductivity + REAL :: TD_HAIL + REAL :: TD_DEPTH + REAL :: TD_HEMD + REAL,DIMENSION(1:NSOIL) :: ZLAYER + REAL,DIMENSION(1:NSOIL) :: OVRFC + INTEGER :: NDRAINS + REAL :: TDDC + REAL,DIMENSION(1:NSOIL) :: RMSH2O + REAL :: QTLDRN1 + REAL :: TD_DD + REAL :: OVRFCS +#ifdef WRF_HYDRO + REAL :: WATBLED ! water table depth estimated in WRF-Hydro fine grids +#endif +!---------------------------------------------------------------------------- + + TD_SATZ = 0. + DTOPL = 0. + TD_LQ = 0. + TD_TTSZ = 0. + TDDC = parameters%TD_DCOEF* 1000. * DT/(24. * 3600.) ! m per day to mm per timestep + +! Thickness of soil layers + DO K = 1, NSOIL + IF (K == 1) THEN + ZLAYER(K) = -1 * ZSOIL(K) + ELSE + ZLAYER(K) = (ZSOIL(K-1)-ZSOIL(K)) + END IF + END DO + +#ifdef WRF_HYDRO +! Depth to water table, m + YY = WATBLED +#else + CALL TD_FINDZWAT(parameters,NSOIL,SMC,SH2O,SICE,ZSOIL,ZLAYER,ZWT) + !CALL ZWTEQ (parameters,NSOIL, NSNOW, ZSOIL, DZSNSO, SH2O, ZWT) +! Depth to water table, m + YY = ZWT +#endif + IF (YY .GT. parameters%TD_ADEPTH) YY = parameters%TD_ADEPTH + +! Depth of saturated zone + DO K=1, NSOIL + IF (YY .GT. (-1*ZSOIL(K))) THEN + TD_SATZ(K) = 0. + ELSE + TD_SATZ(K) = (-1*ZSOIL(K)) - YY + XX = (-1*ZSOIL(K)) - DTOPL + IF(TD_SATZ(K) .GT. XX) TD_SATZ(K) = XX + END IF + + !print*,"K = ", K + !print*,"-1*ZSOIL(K)=",-1*ZSOIL(K) + !print*,"TD_SATZ(K)=",TD_SATZ(K) + !print*,"DTOPL=",DTOPL + + DTOPL = -1*ZSOIL(K) + END DO +! amount of water over field capacity + OVRFCS = 0. + DO K=1, NSOIL + OVRFC(K) = (SH2O(K) - (parameters%SMCREF(K)-SICE(K))) * ZLAYER(K) * 1000. !mm + IF(OVRFC(K) .LT. 0.)OVRFC(K) = 0. + OVRFCS = OVRFCS + OVRFC(K) + END DO + +! lateral hydr. conductivity and total lateral flow + DO K=1, NSOIL + KLATK(K)= WCND(K)*parameters%KLAT_FAC * DT ! m/s to m/timestep + TD_LQ = TD_LQ + (TD_SATZ(K) * KLATK(K)) + TD_TTSZ = TD_TTSZ + TD_SATZ(K) + END DO + + IF (TD_TTSZ .LT. 0.001) TD_TTSZ = 0.001 ! unit is m + IF (TD_LQ .LT. 0.001) TD_LQ = 0. ! unit is m + + KLAT = TD_LQ/TD_TTSZ ! lateral hydraulic conductivity per timestep + + TD_DD = parameters%TD_ADEPTH - parameters%TD_DDRAIN + + CALL TD_EQUIVALENT_DEPTH (TD_DD, & + parameters%TD_SPAC, & + parameters%TD_RADI, & + TD_HAIL) + + TD_DEPTH= TD_HAIL + parameters%TD_DDRAIN + TD_HEMD = parameters%TD_DDRAIN - YY + + IF (TD_HEMD .LE. 0.) THEN + QTLDRN = 0. + ELSE + QTLDRN = ((8.*KLAT*TD_HAIL*TD_HEMD) + (4.*KLAT*TD_HEMD*TD_HEMD))& ! m per timestep + /(parameters%TD_SPAC*parameters%TD_SPAC) + END IF + + QTLDRN = QTLDRN * 1000.0 ! m per timestep to mm/timestep /one tile + + IF(QTLDRN .LE. 0.) QTLDRN = 0. + IF(QTLDRN .GT. TDDC) QTLDRN = TDDC + + NDRAINS = INT(DX/parameters%TD_SPAC) + QTLDRN = QTLDRN * NDRAINS + IF(QTLDRN .GT. OVRFCS) QTLDRN = OVRFCS + + !print*,"TD_TTSZ =",TD_TTSZ + !print*,"TD_LQ =",TD_LQ + !print*,"KLAT =",KLAT + !print*,"TD_HAIL =",TD_HAIL + !print*,"TD_DEPTH=",TD_DEPTH + !print*,"YY =",YY + !print*,"TD_HEMD =",TD_HEMD + !print*,"QTLDRN =",QTLDRN + !print*,"NDRAINS =",NDRAINS + !print*,"TD_DD =",TD_DD + !print*,"TDDC =",TDDC + !print*,"OVRFCS =",OVRFCS +! update soil moisture after drainage: moisture drains from top to bottom + + QTLDRN1 = QTLDRN + + DO K=1, NSOIL + IF(QTLDRN1 .GT. 0.) THEN + IF((TD_SATZ(K) .GT. 0.) .AND. (OVRFC(K) .GT. 0.)) THEN + + RMSH2O(K) = OVRFC(K) - QTLDRN1 ! remaining water after tile drain + + IF (RMSH2O(K) .GT. 0.) THEN + SH2O(K) = (parameters%SMCREF(K) - SICE (K)) + RMSH2O(K)/(ZLAYER(K) * 1000.) + SMC(K) = SH2O(K) + SICE (K) + EXIT + ELSE + SH2O(K) = (parameters%SMCREF(K) - SICE (K)) + SMC(K) = SH2O(K) + SICE (K) + QTLDRN1 = QTLDRN1 - OVRFC(K) + END IF + + END IF + END IF + END DO + + QTLDRN = QTLDRN/DT ![mm/s] + + END SUBROUTINE TILE_HOOGHOUDT + +!----------------------------------------------------------------------- + + SUBROUTINE TD_EQUIVALENT_DEPTH (TD_D,TD_L,TD_RD,TD_DE) + +! ---------------------------------------------------------------------- +! calculate tile drainage equivalent depth from d. +! ---------------------------------------------------------------------- + IMPLICIT NONE +! ---------------------------------------------------------------------- + REAL, INTENT(IN) :: TD_D + REAL, INTENT(IN) :: TD_L + REAL, INTENT(IN) :: TD_RD + REAL, INTENT(OUT) :: TD_DE + REAL :: PII = 22./7. + REAL :: TD_X + REAL :: TD_FX, EX,TERM + INTEGER :: I +!------------------------------------- + + TD_FX = 0. + EX = 0. + TERM = 0. + + TD_X = (2.0*PII*TD_D)/TD_L + + IF (TD_X .GT. 0.5) THEN + DO I=1,45,2 + EX = EXP(-2.0*I*TD_X) + TERM = (4.0*EX)/(I*(1-EX)) + TD_FX = TD_FX + TERM + IF(TERM .LT. 1.E-6) THEN + TD_DE = ((PII*TD_L)/8.0)/(LOG(TD_L/(PII*TD_RD))+TD_FX) + EXIT + END IF + END DO + ELSE IF (TD_X .LT. 1.E-8) THEN + TD_DE = TD_D + ELSE + TD_FX = ((PII*PII)/(4*TD_X))+(LOG(TD_X/(2*PII))) + TD_DE = ((PII*TD_L)/8.0)/(LOG(TD_L/(PII*TD_RD))+TD_FX) + END IF + + IF (TD_DE .LT. 0. .AND. I .LE. 2) TD_DE = TD_D + + END SUBROUTINE TD_EQUIVALENT_DEPTH + +!---------------------------------------------------------------------------- + + SUBROUTINE TD_FINDZWAT(parameters,NSOIL,SMC,SH2O,SICE,ZSOIL,SLDPTH,WATBLED) + +!---------------------------------------------------------------------------- +! Calculate watertable depth as on WRF-Hydro/NWM +!---------------------------------------------------------------------------- + IMPLICIT NONE + +!-------- DECLARATIONS ------------------------ + type (noahmp_parameters), intent(in) :: parameters + INTEGER, INTENT(IN) :: NSOIL + REAL, INTENT(IN), DIMENSION(NSOIL) :: SMC,SH2O,SICE + REAL, INTENT(IN), DIMENSION(NSOIL) :: ZSOIL + REAL, INTENT(IN), DIMENSION(NSOIL) :: SLDPTH + REAL, INTENT(OUT) :: WATBLED + REAL :: CWATAVAIL + INTEGER :: SATLYRCHK + +! Local Variables + INTEGER :: K,i,j +!------------------------------------------------------------ + + SATLYRCHK = 0 !set flag for sat. layers + CWATAVAIL = 0. !set wat avail for subsfc rtng = 0. + + DO K=NSOIL,1,-1 + IF ( (SMC(K).GE.parameters%SMCREF(K)).AND.& + (parameters%SMCREF(K) .GT.parameters%SMCWLT(K)) ) THEN + IF((SATLYRCHK .EQ.K+1) .OR. (K.EQ.NSOIL)) SATLYRCHK = K + END IF + END DO + + IF (SATLYRCHK .NE. 0) THEN + IF (SATLYRCHK .NE. 1) then ! soil column is partially sat. + WATBLED = -ZSOIL(SATLYRCHK-1) + ELSE ! soil column is fully saturated to sfc. + WATBLED = 0. + END IF + DO K = SATLYRCHK,NSOIL + CWATAVAIL = CWATAVAIL+(SMC(K)-parameters%SMCREF(K))*SLDPTH(K) + END DO + ELSE ! no saturated layers... + WATBLED = -ZSOIL(NSOIL) + SATLYRCHK = NSOIL + 1 + END IF + + END SUBROUTINE TD_FINDZWAT + +!== end tile drainage subroutines ================================================================== + !== begin wdfcnd1 ================================================================================== SUBROUTINE WDFCND1 (parameters,WDF,WCND,SMC,FCR,ISOIL) @@ -9650,9 +10952,10 @@ END SUBROUTINE EMERG !== begin noahmp_options =========================================================================== - subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & - iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, & - iopt_rsf , iopt_soil, iopt_pedo, iopt_crop, iopt_irr, iopt_irrm) + subroutine noahmp_options( idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & + iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot ,iopt_stc , & + iopt_rsf ,iopt_soil ,iopt_pedo ,iopt_crop ,iopt_irr ,iopt_irrm, & + iopt_infdv, iopt_tdrn) implicit none @@ -9682,6 +10985,8 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ! 1 -> sprinkler ON ! 2 -> micro/drip ON ! 3 -> flood irrigation ON + INTEGER, INTENT(IN) :: iopt_infdv!infiltration options for dynamic VIC (1->Philip; 2-> Green-Ampt;3->Smith-Parlange) + INTEGER, INTENT(IN) :: iopt_tdrn !tile drainage (0->none; 1-> simple 2->Hooghoudt's) ! ------------------------------------------------------------------------------------------------- dveg = idveg @@ -9703,7 +11008,9 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc opt_crop = iopt_crop opt_irr = iopt_irr opt_irrm = iopt_irrm - + opt_infdv= iopt_infdv + opt_tdrn = iopt_tdrn + end subroutine noahmp_options END MODULE MODULE_SF_NOAHMPLSM @@ -9799,16 +11106,23 @@ MODULE NOAHMP_TABLES INTEGER :: SLCATS - REAL :: BEXP_TABLE(MAX_SOILTYP) !maximum intercepted h2o per unit lai+sai (mm) + REAL :: BEXP_TABLE(MAX_SOILTYP) !maximum intercepted h2o per unit lai+sai (mm) REAL :: SMCDRY_TABLE(MAX_SOILTYP) !characteristic leaf dimension (m) - REAL :: F1_TABLE(MAX_SOILTYP) !momentum roughness length (m) + REAL :: F1_TABLE(MAX_SOILTYP) !momentum roughness length (m) REAL :: SMCMAX_TABLE(MAX_SOILTYP) !top of canopy (m) REAL :: SMCREF_TABLE(MAX_SOILTYP) !bottom of canopy (m) REAL :: PSISAT_TABLE(MAX_SOILTYP) !tree density (no. of trunks per m2) REAL :: DKSAT_TABLE(MAX_SOILTYP) !tree crown radius (m) REAL :: DWSAT_TABLE(MAX_SOILTYP) !monthly stem area index, one-sided REAL :: SMCWLT_TABLE(MAX_SOILTYP) !monthly leaf area index, one-sided - REAL :: QUARTZ_TABLE(MAX_SOILTYP) !single-side leaf area per Kg [m2/kg] + REAL :: QUARTZ_TABLE(MAX_SOILTYP) !single-side leaf area per Kg [m2/kg] + REAL :: BVIC_TABLE(MAX_SOILTYP) !VIC model infiltration parameter (-) for opt_run=6 + REAL :: AXAJ_TABLE(MAX_SOILTYP) !Xinanjiang: Tension water distribution inflection parameter [-] for opt_run=7 + REAL :: BXAJ_TABLE(MAX_SOILTYP) !Xinanjiang: Tension water distribution shape parameter [-] for opt_run=7 + REAL :: XXAJ_TABLE(MAX_SOILTYP) !Xinanjiang: Free water distribution shape parameter [-] for opt_run=7 + REAL :: BDVIC_TABLE(MAX_SOILTYP) !VIC model infiltration parameter (-) + REAL :: GDVIC_TABLE(MAX_SOILTYP) !mean capilary drive (m) + REAL :: BBVIC_TABLE(MAX_SOILTYP) !heterogeniety parameter for DVIC infiltration [-] ! GENPARM.TBL parameters @@ -9932,6 +11246,19 @@ MODULE NOAHMP_TABLES REAL :: RTCT_TABLE(NCROP,NSTAGE) ! root to grain REAL :: BIO2LAI_TABLE(NCROP) ! leaf are per living leaf biomass [m^2/kg] +! tile drainage parameters + REAL :: TDSMCFAC_TABLE(MAX_SOILTYP) + REAL :: TD_DC_TABLE(MAX_SOILTYP) + INTEGER :: TD_DEPTH_TABLE(MAX_SOILTYP) + INTEGER :: DRAIN_LAYER_OPT_TABLE + REAL :: TD_DCOEF_TABLE(MAX_SOILTYP) + REAL :: TD_D_TABLE(MAX_SOILTYP) + REAL :: TD_ADEPTH_TABLE(MAX_SOILTYP) + REAL :: TD_RADI_TABLE(MAX_SOILTYP) + REAL :: TD_SPAC_TABLE(MAX_SOILTYP) + REAL :: TD_DDRAIN_TABLE(MAX_SOILTYP) + REAL :: KLAT_FAC_TABLE(MAX_SOILTYP) + ! MPTABLE.TBL optional parameters REAL :: sr2006_theta_1500t_a ! sand coefficient @@ -10268,6 +11595,13 @@ subroutine read_mp_soil_parameters() FRZK_TABLE = -1.E36 ZBOT_TABLE = -1.E36 CZIL_TABLE = -1.E36 + BVIC_TABLE = -1.E36 + AXAJ_TABLE = -1.E36 + BXAJ_TABLE = -1.E36 + XXAJ_TABLE = -1.E36 + BDVIC_TABLE = -1.E36 + GDVIC_TABLE = -1.E36 + BBVIC_TABLE = -1.E36 ! !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL @@ -10292,9 +11626,10 @@ subroutine read_mp_soil_parameters() CALL wrf_message ( message ) DO LC=1,SLCATS - READ (21,*) ITMP,BEXP_TABLE(LC),SMCDRY_TABLE(LC),F1_TABLE(LC),SMCMAX_TABLE(LC), & + READ (21,*) ITMP,BEXP_TABLE(LC),SMCDRY_TABLE(LC),F1_TABLE(LC),SMCMAX_TABLE(LC), & SMCREF_TABLE(LC),PSISAT_TABLE(LC),DKSAT_TABLE(LC), DWSAT_TABLE(LC), & - SMCWLT_TABLE(LC), QUARTZ_TABLE(LC) + SMCWLT_TABLE(LC), QUARTZ_TABLE(LC),BVIC_TABLE(LC), AXAJ_TABLE(LC), & + BXAJ_TABLE(LC),XXAJ_TABLE(LC),BDVIC_TABLE(LC),BBVIC_TABLE(LC),GDVIC_TABLE(LC) ENDDO CLOSE (21) @@ -10832,6 +12167,64 @@ subroutine read_mp_irrigation_parameters() end subroutine read_mp_irrigation_parameters + subroutine read_tiledrain_parameters() + implicit none + integer :: ierr + logical :: file_named + REAL, DIMENSION(MAX_SOILTYP) :: TDSMC_FAC + INTEGER, DIMENSION(MAX_SOILTYP) :: TD_DEPTH + REAL, DIMENSION(MAX_SOILTYP) :: TD_DC + INTEGER :: DRAIN_LAYER_OPT + REAL, DIMENSION(MAX_SOILTYP) :: TD_DCOEF + REAL, DIMENSION(MAX_SOILTYP) :: TD_D + REAL, DIMENSION(MAX_SOILTYP) :: TD_ADEPTH + REAL, DIMENSION(MAX_SOILTYP) :: TD_RADI + REAL, DIMENSION(MAX_SOILTYP) :: TD_SPAC + REAL, DIMENSION(MAX_SOILTYP) :: TD_DDRAIN + REAL, DIMENSION(MAX_SOILTYP) :: KLAT_FAC + NAMELIST / noahmp_tiledrain_parameters /DRAIN_LAYER_OPT,TDSMC_FAC,TD_DEPTH,TD_DC,& + TD_DCOEF,TD_D,TD_ADEPTH,TD_RADI,TD_SPAC,TD_DDRAIN,& + KLAT_FAC + ! Initialize our variables to bad values, so that if the namelist read fails, we come to a screeching halt as soon as we try to use anything. + TDSMCFAC_TABLE = -99999 + TD_DEPTH_TABLE = -99999 + TD_DC_TABLE = -99999 + DRAIN_LAYER_OPT_TABLE = -99999 + TD_DCOEF_TABLE = -99999 + TD_D_TABLE = -99999 + TD_ADEPTH_TABLE = -99999 + TD_RADI_TABLE = -99999 + TD_SPAC_TABLE = -99999 + TD_DDRAIN_TABLE = -99999 + KLAT_FAC_TABLE = -99999 + + inquire( file='MPTABLE.TBL', exist=file_named ) + if ( file_named ) then + open(15, file="MPTABLE.TBL", status='old', form='formatted', action='read', iostat=ierr) + else + open(15, status='old', form='formatted', action='read', iostat=ierr) + end if + if (ierr /= 0) then + write(*,'("WARNING: Cannot find file MPTABLE.TBL")') + call wrf_error_fatal("STOP in Noah-MP read_tiledrain_parameters") + endif + read(15,noahmp_tiledrain_parameters) + close(15) + TDSMCFAC_TABLE = TDSMC_FAC + TD_DEPTH_TABLE = TD_DEPTH + DRAIN_LAYER_OPT_TABLE = DRAIN_LAYER_OPT + TD_DC_TABLE = TD_DC + + TD_DCOEF_TABLE = TD_DCOEF + TD_D_TABLE = TD_D + TD_ADEPTH_TABLE = TD_ADEPTH + TD_RADI_TABLE = TD_RADI + TD_SPAC_TABLE = TD_SPAC + TD_DDRAIN_TABLE = TD_DDRAIN + KLAT_FAC_TABLE = KLAT_FAC + + end subroutine read_tiledrain_parameters + subroutine read_mp_optional_parameters() implicit none integer :: ierr @@ -10871,7 +12264,6 @@ subroutine read_mp_optional_parameters() read(15,noahmp_optional_parameters) close(15) - end subroutine read_mp_optional_parameters END MODULE NOAHMP_TABLES