diff --git a/drivers/hrldas/ConfigVarInTransferMod.F90 b/drivers/hrldas/ConfigVarInTransferMod.F90 index 5c1faa45..41475f21 100644 --- a/drivers/hrldas/ConfigVarInTransferMod.F90 +++ b/drivers/hrldas/ConfigVarInTransferMod.F90 @@ -62,6 +62,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) noahmp%config%nmlist%OptSnowCompaction = NoahmpIO%IOPT_COMPACT noahmp%config%nmlist%OptWetlandModel = NoahmpIO%IOPT_WETLAND noahmp%config%nmlist%OptSnowCoverGround = NoahmpIO%IOPT_SCF + noahmp%config%nmlist%OptRoot = NoahmpIO%IOPT_ROOT if ( noahmp%config%nmlist%OptSnowAlbedo == 3 ) then ! SNICAR namelist noahmp%config%nmlist%OptSnicarSnowShape = NoahmpIO%SNICAR_SNOWSHAPE_OPT diff --git a/drivers/hrldas/NoahmpDriverMainMod.F90 b/drivers/hrldas/NoahmpDriverMainMod.F90 index 9132b1ac..a1d15df2 100644 --- a/drivers/hrldas/NoahmpDriverMainMod.F90 +++ b/drivers/hrldas/NoahmpDriverMainMod.F90 @@ -110,7 +110,7 @@ subroutine NoahmpDriverMain(NoahmpIO) do K = 2, NoahmpIO%NSOIL NoahmpIO%ZSOIL(K) = -NoahmpIO%DZS(K) + NoahmpIO%ZSOIL(K-1) enddo - + JLOOP : do J = NoahmpIO%JTS, NoahmpIO%JTE NoahmpIO%J = J diff --git a/drivers/hrldas/NoahmpIOVarInitMod.F90 b/drivers/hrldas/NoahmpIOVarInitMod.F90 index 78f6fe59..90dbb1b1 100644 --- a/drivers/hrldas/NoahmpIOVarInitMod.F90 +++ b/drivers/hrldas/NoahmpIOVarInitMod.F90 @@ -424,6 +424,14 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) if ( .not. allocated (NoahmpIO%WCAP) ) allocate ( NoahmpIO%WCAP (XSTART:XEND, YSTART:YEND) ) ! maximum wetland capacity [m] endif + ! Needed for DynaRoot scheme (OPT_ROOT=1) + if ( NoahmpIO%IOPT_ROOT == 1 ) then + if ( .not. allocated (NoahmpIO%EASYXY) ) allocate ( NoahmpIO%EASYXY (XSTART:XEND,1:NSOIL,YSTART:YEND) ) ! ease function (-) + if ( .not. allocated (NoahmpIO%ROOTACTIVITYXY) ) allocate ( NoahmpIO%ROOTACTIVITYXY (XSTART:XEND,1:NSOIL,YSTART:YEND) ) ! root activity function (-) + if ( .not. allocated (NoahmpIO%INACTIVEXY) ) allocate ( NoahmpIO%INACTIVEXY (XSTART:XEND,1:NSOIL,YSTART:YEND) ) ! number of time steps without active roots (-) + if ( .not. allocated (NoahmpIO%MAXUPTAKEXY) ) allocate ( NoahmpIO%MAXUPTAKEXY (XSTART:XEND, YSTART:YEND) ) ! root layer depth (-) + endif + ! Single- and Multi-layer Urban Models if ( NoahmpIO%SF_URBAN_PHYSICS > 0 ) then @@ -891,6 +899,14 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) NoahmpIO%WCAP = undefined_real endif + ! DynaRoot scheme + if ( NoahmpIO%IOPT_ROOT == 1 ) then + NoahmpIO%EASYXY = 0.0 + NoahmpIO%ROOTACTIVITYXY = 0.0 + NoahmpIO%INACTIVEXY = undefined_real + NoahmpIO%MAXUPTAKEXY = undefined_real + endif + ! spatial varying soil texture if ( NoahmpIO%IOPT_SOIL > 1 ) then NoahmpIO%SOILCL1 = undefined_real diff --git a/drivers/hrldas/NoahmpIOVarType.F90 b/drivers/hrldas/NoahmpIOVarType.F90 index 42934326..25f2425c 100644 --- a/drivers/hrldas/NoahmpIOVarType.F90 +++ b/drivers/hrldas/NoahmpIOVarType.F90 @@ -66,6 +66,7 @@ module NoahmpIOVarType integer :: IOPT_COMPACT ! snowpack compaction (1->Anderson1976; 2->Abolafia-Rosenzweig2024) integer :: IOPT_SCF ! snow cover fraction (1->NiuYang07; 2->Abolafia-Rosenzweig2025) integer :: IOPT_WETLAND ! wetland model option (0->off; 1->Zhang2022 fixed parameter; 2->Zhang2022 read in 2D parameter) + integer :: IOPT_ROOT ! root scheme option (0->off; 1->Bieri2025 dynamic root scheme) real(kind=kind_noahmp) :: XICE_THRESHOLD ! fraction of grid determining seaice real(kind=kind_noahmp) :: JULIAN ! Julian day real(kind=kind_noahmp) :: DTBL ! timestep [s] @@ -501,6 +502,14 @@ module NoahmpIOVarType real(kind=kind_noahmp), allocatable, dimension(:,:) :: FSATMX ! maximum saturated fraction real(kind=kind_noahmp), allocatable, dimension(:,:) :: WCAP ! maximum wetland capacity [m] +!------------------------------------------------------------------------ +! Needed for DynaRoot scheme (OPT_ROOT=1) +!------------------------------------------------------------------------ + real(kind=kind_noahmp), allocatable, dimension(:,:,:) :: EASYXY ! ease function (-) + real(kind=kind_noahmp), allocatable, dimension(:,:,:) :: ROOTACTIVITYXY ! root activity function (-) + real(kind=kind_noahmp), allocatable, dimension(:,:,:) :: INACTIVEXY ! number of timesteps without active roots (-) + real(kind=kind_noahmp), allocatable, dimension(:,:) :: MAXUPTAKEXY ! bottom of deepest root water uptake layer (m) + !------------------------------------------------------------------------ ! Single- and Multi-layer Urban Models !------------------------------------------------------------------------ diff --git a/drivers/hrldas/NoahmpInitMainMod.F90 b/drivers/hrldas/NoahmpInitMainMod.F90 index 5ca04da0..d3e8d60f 100644 --- a/drivers/hrldas/NoahmpInitMainMod.F90 +++ b/drivers/hrldas/NoahmpInitMainMod.F90 @@ -269,6 +269,15 @@ subroutine NoahmpInitMain(NoahmpIO) NoahmpIO%STEPWTD = max(NoahmpIO%STEPWTD,1) endif + ! initialize variables for DynaRoot (IOPT_ROOT == 1) + if ( NoahmpIO%IOPT_ROOT == 1 ) then + ! initialize time step counter + NoahmpIO%INACTIVEXY(:,1,:) = 0.0 + NoahmpIO%INACTIVEXY(:,2:NoahmpIO%NSOIL,:) = 366.0*(86400.0 / NoahmpIO%DTBL) + ! initialize root depth + NoahmpIO%NROOT_TABLE = 1 + endif + ! initialize soil albedo NoahmpIO%ALBSOILDIRXY = 0.0 NoahmpIO%ALBSOILDIFXY = 0.0 diff --git a/drivers/hrldas/NoahmpReadNamelistMod.F90 b/drivers/hrldas/NoahmpReadNamelistMod.F90 index 444db52d..3f99aa02 100644 --- a/drivers/hrldas/NoahmpReadNamelistMod.F90 +++ b/drivers/hrldas/NoahmpReadNamelistMod.F90 @@ -96,6 +96,7 @@ subroutine NoahmpReadNamelist(NoahmpIO) integer :: dvic_infiltration_option = 1 integer :: tile_drainage_option = 0 integer :: wetland_option = 0 + integer :: root_option = 0 integer :: split_output_count = 1 logical :: skip_first_output = .false. integer :: khour = -9999 @@ -155,7 +156,7 @@ subroutine NoahmpReadNamelist(NoahmpIO) frozen_soil_option, radiative_transfer_option, snow_albedo_option, & snow_thermal_conductivity, surface_runoff_option, subsurface_runoff_option, & pcp_partition_option, tbot_option, temp_time_scheme_option, wetland_option, & - glacier_option, surface_resistance_option, snow_compaction_option, & + root_option, glacier_option, surface_resistance_option, snow_compaction_option, & irrigation_option, irrigation_method, dvic_infiltration_option, & tile_drainage_option,soil_data_option, pedotransfer_option, crop_option, & sf_urban_physics,use_wudapt_lcz,num_urban_hi,urban_atmosphere_thickness, & @@ -379,6 +380,7 @@ subroutine NoahmpReadNamelist(NoahmpIO) NoahmpIO%IOPT_TDRN = tile_drainage_option NoahmpIO%IOPT_COMPACT = snow_compaction_option NoahmpIO%IOPT_WETLAND = wetland_option + NoahmpIO%IOPT_ROOT = root_option NoahmpIO%IOPT_SCF = snow_cover_option ! basic model setup variables NoahmpIO%indir = indir diff --git a/drivers/hrldas/WaterVarInTransferMod.F90 b/drivers/hrldas/WaterVarInTransferMod.F90 index 3ecdae39..92f7baaa 100644 --- a/drivers/hrldas/WaterVarInTransferMod.F90 +++ b/drivers/hrldas/WaterVarInTransferMod.F90 @@ -285,6 +285,14 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO) noahmp%water%param%SoilMoistureDry = 0.40 endif + ! dynamic root scheme + if ( noahmp%config%nmlist%OptRoot == 1 ) then + noahmp%water%state%EaseFunction = NoahmpIO%EASYXY (I,1:NumSoilLayer,J) + noahmp%water%state%RootActivity = NoahmpIO%ROOTACTIVITYXY(I,1:NumSoilLayer,J) + noahmp%water%state%InactiveTimeSteps = NoahmpIO%INACTIVEXY (I,1:NumSoilLayer,J) + noahmp%water%state%MaxUptakeDepth = NoahmpIO%MAXUPTAKEXY (I,J) + endif + if ( SoilType(1) /= 14 ) then noahmp%water%param%SoilImpervFracCoeff = noahmp%water%param%GroundFrzCoeff * & ((noahmp%water%param%SoilMoistureSat(1) / & diff --git a/drivers/hrldas/WaterVarOutTransferMod.F90 b/drivers/hrldas/WaterVarOutTransferMod.F90 index b404bd31..c9748b4d 100644 --- a/drivers/hrldas/WaterVarOutTransferMod.F90 +++ b/drivers/hrldas/WaterVarOutTransferMod.F90 @@ -170,6 +170,14 @@ subroutine WaterVarOutTransfer(noahmp, NoahmpIO) NoahmpIO%FSATXY (I,J) = noahmp%water%state%SoilSaturateFrac endif + ! DynaRoot scheme + if ( noahmp%config%nmlist%OptRoot == 1 ) then + NoahmpIO%EASYXY (I,1:NumSoilLayer,J) = noahmp%water%state%EaseFunction(1:NumSoilLayer) + NoahmpIO%ROOTACTIVITYXY(I,1:NumSoilLayer,J) = noahmp%water%state%RootActivity(1:NumSoilLayer) + NoahmpIO%INACTIVEXY (I,1:NumSoilLayer,J) = noahmp%water%state%InactiveTimeSteps(1:NumSoilLayer) + NoahmpIO%MAXUPTAKEXY (I,J) = noahmp%water%state%MaxUptakeDepth + endif + #ifdef WRF_HYDRO NoahmpIO%infxsrt (I,J) = max(noahmp%water%flux%RunoffSurface, 0.0) ! mm, surface runoff NoahmpIO%soldrain (I,J) = max(noahmp%water%flux%RunoffSubsurface, 0.0) ! mm, underground runoff diff --git a/src/ConfigVarInitMod.F90 b/src/ConfigVarInitMod.F90 index 55f72965..55597806 100644 --- a/src/ConfigVarInitMod.F90 +++ b/src/ConfigVarInitMod.F90 @@ -49,6 +49,7 @@ subroutine ConfigVarInitDefault(noahmp) noahmp%config%nmlist%OptGlacierTreatment = undefined_int noahmp%config%nmlist%OptSnowCompaction = undefined_int noahmp%config%nmlist%OptWetlandModel = undefined_int + noahmp%config%nmlist%OptRoot = undefined_int noahmp%config%nmlist%OptSnicarSnowShape = undefined_int noahmp%config%nmlist%OptSnicarRTSolver = undefined_int noahmp%config%nmlist%OptSnicarBandNum = undefined_int diff --git a/src/ConfigVarType.F90 b/src/ConfigVarType.F90 index b0a41bf8..e4092ce1 100644 --- a/src/ConfigVarType.F90 +++ b/src/ConfigVarType.F90 @@ -131,6 +131,9 @@ module ConfigVarType ! 0 -> No Wetland model ! 1 -> Single-point/uniform parameter (Zhang, et al. 2022 WRR) ! 2 -> 2-D regional parameter input (Zhang, et al. 2022 WRR) + integer :: OptRoot ! option for root scheme + ! 0 -> No root scheme + ! 1 -> DynaRoot scheme (Bieri et al 2025; Fan et al. 2017) integer :: OptSnicarSnowShape ! options for snow grain shape in SNICAR (He et al. 2017 JC) ! 1 -> sphere ! 2 -> spheroid diff --git a/src/DynaRootGrowDieMod.F90 b/src/DynaRootGrowDieMod.F90 new file mode 100644 index 00000000..ecfbacaf --- /dev/null +++ b/src/DynaRootGrowDieMod.F90 @@ -0,0 +1,165 @@ +module DynaRootGrowDieMod + +!!! determine "growth" or "death" of roots based on root water uptake profile (Bieri et al. 2025) +!!! based on Fan et al. (2017) + + use Machine + use NoahmpVarType + + implicit none + +contains + + subroutine DynaRootGrowDie(noahmp) + +! ------------------------ Code history ----------------------------------- +! Original code: Carolina Bieri +! ------------------------------------------------------------------------- + + implicit none + +! in & out variables + type(noahmp_type), intent(inout) :: noahmp + +! local variables + integer :: IndSoil ! loop index + integer :: WaterTableLayer ! layer that contains water table + integer :: AboveWaterTableLayer ! layer above water table + integer :: DynaRootLayer ! deepest layer with active roots + integer :: MaxUptakeLayer ! deepest layer with active root water uptake + integer :: SoilIceFactor ! multiplication factor corresponding to IndicatorIceSfc + real(kind=kind_noahmp) :: LayerMidpoint ! depth to midpoint of soil layer [m] + real(kind=kind_noahmp) :: MaximumInactiveDays ! max number of days without active roots until layer is inactive [s] + real(kind=kind_noahmp) :: MaximumEaseFunction ! maximum column-wise ease function value + real(kind=kind_noahmp) :: SumEaseFunction ! column-wise sum of ease function + real(kind=kind_noahmp) :: EffectiveCanopyHeight ! canopy height used in ease function calculation [m] + real(kind=kind_noahmp), parameter :: LeafWiltMatPotential = -150.0 ! leaf wilting point matric potential [m] + + real(kind=kind_noahmp), allocatable, dimension(:) :: DynaRootMask ! indicator for presence of active roots (1=roots active; 0=roots not active) + real(kind=kind_noahmp), allocatable, dimension(:) :: ThicknessSoilLayerTmp ! ThicknessSoilLayer accounting for location of water table, only for root scheme [mm] + real(kind=kind_noahmp), allocatable, dimension(:) :: InactiveDays ! number of days without active roots [s] +! -------------------------------------------------------------------- + associate( & + ThicknessSoilLayer => noahmp%config%domain%ThicknessSoilLayer ,& ! in, soil layer thickness [m] + DepthSoilLayer => noahmp%config%domain%DepthSoilLayer ,& ! in, layer-bottom depth from soil surface [m] + NumSoilLayer => noahmp%config%domain%NumSoilLayer ,& ! in, number of soil layers + MainTimeStep => noahmp%config%domain%MainTimeStep ,& ! in, main Noah-MP model time step [s] + IndicatorIceSfc => noahmp%config%domain%IndicatorIceSfc ,& ! in, indicator for sea ice and land ice (1=sea ice; -1=land ice; 0=non-ice) + WaterTableDepth => noahmp%water%state%WaterTableDepth ,& ! in, depth to water table [m] + SoilMatPotential => noahmp%water%state%SoilMatPotential ,& ! in, soil matrix potential [m] + HeightCanopyTop => noahmp%energy%param%HeightCanopyTop ,& ! in, height of canopy top [m] + InactiveTimeSteps => noahmp%water%state%InactiveTimeSteps ,& ! inout, number of time steps without active roots + NumSoilLayerRoot => noahmp%water%param%NumSoilLayerRoot ,& ! inout, number of soil layers with root present + EaseFunction => noahmp%water%state%EaseFunction ,& ! out, ease function + RootActivity => noahmp%water%state%RootActivity ,& ! out, root activity function + MaxUptakeDepth => noahmp%water%state%MaxUptakeDepth & ! out, bottom of deepest root water uptake layer [m] + ) +! ---------------------------------------------------------------------- + + if (.not. allocated(DynaRootMask) ) allocate(DynaRootMask (1:NumSoilLayer)) + if (.not. allocated(ThicknessSoilLayerTmp) ) allocate(ThicknessSoilLayerTmp (1:NumSoilLayer)) + if (.not. allocated(InactiveDays) ) allocate(InactiveDays (1:NumSoilLayer)) + + ! initialize + DynaRootMask = 0.0 + ThicknessSoilLayerTmp = undefined_real + EaseFunction = 0.0 + RootActivity = 0.0 + InactiveDays = 0.0 + + MaximumInactiveDays = 31536000.0 ! seconds in a year + + ThicknessSoilLayerTmp = ThicknessSoilLayer(1:NumSoilLayer) ! soil thickness variable for use in root scheme + + ! determine WaterTableLayer + IF(WaterTableDepth + 1.0E-6 .lt. DepthSoilLayer(NumSoilLayer)) THEN ! set WaterTableLayer to NumSoilLayer if water table below resolved layers + WaterTableLayer = NumSoilLayer + ELSE + DO IndSoil=NumSoilLayer,1,-1 ! find WaterTableLayer if within resolved layers + IF(WaterTableDepth + 1.0E-6 < DepthSoilLayer(IndSoil)) EXIT + ENDDO + AboveWaterTableLayer = IndSoil ! layer above water table + WaterTableLayer = MIN(AboveWaterTableLayer+1,NumSoilLayer) + + ! Adjust layer thickness accordingly + IF (AboveWaterTableLayer .ge. 1) ThicknessSoilLayerTmp(WaterTableLayer) = -1*(WaterTableDepth-DepthSoilLayer(AboveWaterTableLayer)) + END IF + + ! determine DynaRootLayer - bottom layer of root zone + InactiveDays(1:NumSoilLayer) = InactiveTimeSteps(1:NumSoilLayer)*MainTimeStep ! Convert InactiveTimeSteps to InactiveDays (in seconds) + DO IndSoil = NumSoilLayer,1,-1 + IF(InactiveDays(IndSoil) .le. MaximumInactiveDays)EXIT + END DO + DynaRootLayer = MAX(IndSoil,1) + + EffectiveCanopyHeight = HeightCanopyTop*(2./3.) + ! calculation of ease function and root activity + DO IndSoil = 1, MIN(MIN(DynaRootLayer+1,NumSoilLayer), WaterTableLayer) + IF (IndSoil==1) THEN ! get layer midpoints + LayerMidpoint = -0.05 + ELSE + LayerMidpoint = 0.5 * (DepthSoilLayer(IndSoil-1) + DepthSoilLayer(IndSoil)) + ENDIF + + IF(InactiveDays(IndSoil) .le. MaximumInactiveDays) DynaRootMask(IndSoil) = 1. ! set root mask accordingly + + ! account for frozen soil + IF(IndicatorIceSfc.eq.0) THEN + SoilIceFactor = 1. + ELSE + SoilIceFactor = 0. + ENDIF + + ! calculate ease function + EaseFunction(IndSoil) = MAX(-( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ), 0.) + + END DO + + ! to grow roots anew, the layer has to be easiest to get water from than the + ! current active layers with roots + MaximumEaseFunction = MAXVAL(EaseFunction, DynaRootMask == 1.) + + ! eliminate small root activity + WHERE(EaseFunction .lt. 0.001*MaximumEaseFunction) EaseFunction = 0. + + ! allow root water uptake to extend deeper only if needed + DO IndSoil = 1, MIN(DynaRootLayer+1,NumSoilLayer) + IF(InactiveDays(IndSoil) .gt. MaximumInactiveDays .and. EaseFunction(IndSoil) .lt. MaximumEaseFunction) EaseFunction(IndSoil) = 0. + ENDDO + ! update root zone index + IF(DynaRootLayer .lt. NumSoilLayer)THEN + IF(EaseFunction(DynaRootLayer+1) .gt. 0.) DynaRootLayer = DynaRootLayer+1 + ENDIF + + ! calculate root activity + SumEaseFunction = SUM(EaseFunction*ThicknessSoilLayerTmp(1:NumSoilLayer)) + IF(SumEaseFunction .eq. 0.)THEN + RootActivity = 0. + ELSE + RootActivity = MIN(MAX((EaseFunction*ThicknessSoilLayerTmp(1:NumSoilLayer)) / SumEaseFunction , 0. ), 1. ) + ENDIF + + ! Increment or restart inactive timestep count based on root activity + DO IndSoil = 1, NumSoilLayer + IF(EaseFunction(IndSoil) .eq. 0.)THEN + InactiveTimeSteps(IndSoil) = InactiveTimeSteps(IndSoil) + 1 + ELSE + InactiveTimeSteps(IndSoil) = 0 + ENDIF + ENDDO + + ! Set model number of root layers to number of layers with active roots + NumSoilLayerRoot = DynaRootLayer + + ! Find deepest layer with active root water uptake + DO IndSoil = NumSoilLayer,1,-1 + IF(EaseFunction(IndSoil) .gt. 0.)EXIT + END DO + MaxUptakeLayer = MAX(IndSoil,1) + MaxUptakeDepth = DepthSoilLayer(MaxUptakeLayer) + + end associate + + end subroutine DynaRootGrowDie + +end module DynaRootGrowDieMod diff --git a/src/EnergyMainMod.F90 b/src/EnergyMainMod.F90 index 093fecd7..6e7ec73a 100644 --- a/src/EnergyMainMod.F90 +++ b/src/EnergyMainMod.F90 @@ -53,6 +53,8 @@ module EnergyMainMod use SurfaceEnergyFluxBareGroundMod, only : SurfaceEnergyFluxBareGround use SoilSnowTemperatureMainMod, only : SoilSnowTemperatureMain use SoilSnowWaterPhaseChangeMod, only : SoilSnowWaterPhaseChange + use DynaRootGrowDieMod, only : DynaRootGrowDie + use SoilMatricPotentialMod, only : SoilMatricPotential implicit none @@ -80,6 +82,7 @@ subroutine EnergyMain(noahmp) RadSwDownRefHeight => noahmp%forcing%RadSwDownRefHeight ,& ! in, downward shortwave radiation [W/m2] at reference height OptSnowSoilTempTime => noahmp%config%nmlist%OptSnowSoilTempTime ,& ! in, options for snow/soil temperature time scheme OptSnowCoverGround => noahmp%config%nmlist%OptSnowCoverGround ,& ! in, options for ground snow cover + OptRoot => noahmp%config%nmlist%OptRoot ,& ! in, dynamic root option FlagCropland => noahmp%config%domain%FlagCropland ,& ! in, flag to identify croplands FlagSoilProcess => noahmp%config%domain%FlagSoilProcess ,& ! in, flag to determine if calculating soil processes NumSoilTimeStep => noahmp%config%domain%NumSoilTimeStep ,& ! in, number of time step for calculating soil processes @@ -214,6 +217,10 @@ subroutine EnergyMain(noahmp) ! longwave emissivity for vegetation, ground, total net surface call SurfaceEmissivity(noahmp) + call SoilMatricPotential(noahmp) + + if ( OptRoot == 1 ) call DynaRootGrowDie(noahmp) + ! soil water transpiration factor controlling stomatal resistance and evapotranspiration call SoilWaterTranspiration(noahmp) diff --git a/src/Makefile b/src/Makefile index 270ca5bc..d5ff4142 100644 --- a/src/Makefile +++ b/src/Makefile @@ -139,7 +139,9 @@ OBJS = ConstantDefineMod.o \ SurfaceEmissivityGlacierMod.o \ SurfaceEnergyFluxGlacierMod.o \ SurfaceRadiationGlacierMod.o \ - WaterMainGlacierMod.o + WaterMainGlacierMod.o \ + DynaRootGrowDieMod.o \ + SoilMatricPotentialMod.o all: $(OBJS) @@ -278,7 +280,8 @@ EnergyMainMod.o: ../utility/Machine.o NoahmpVarType.o Const SurfaceEnergyFluxBareGroundMod.o SoilSnowTemperatureMainMod.o \ SoilSnowWaterPhaseChangeMod.o SnowCoverGroundNiu07Mod.o SurfaceEmissivityMod.o \ GroundRoughnessPropertyMod.o PsychrometricVariableMod.o ResistanceGroundEvaporationMod.o \ - SoilWaterTranspirationMod.o SurfaceAlbedoMod.o SurfaceRadiationMod.o SnowCoverGroundAR25Mod.o + SoilWaterTranspirationMod.o SurfaceAlbedoMod.o SurfaceRadiationMod.o SnowCoverGroundAR25Mod.o \ + DynaRootGrowDieMod.o SoilMatricPotentialMod.o NoahmpMainMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ IrrigationPrepareMod.o IrrigationSprinklerMod.o CanopyWaterInterceptMod.o \ PrecipitationHeatAdvectMod.o EnergyMainMod.o WaterMainMod.o AtmosForcingMod.o \ @@ -365,4 +368,5 @@ SurfaceEnergyFluxGlacierMod.o: ../utility/Machine.o NoahmpVarType.o Co VaporPressureSaturationMod.o ResistanceBareGroundMostMod.o SurfaceRadiationGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o WaterMainGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowWaterMainGlacierMod.o - +DynaRootGrowDieMod.o: ../utility/Machine.o NoahmpVarType.o +SoilMatricPotentialMod.o: ../utility/Machine.o NoahmpVarType.o diff --git a/src/SoilMatricPotentialMod.F90 b/src/SoilMatricPotentialMod.F90 new file mode 100644 index 00000000..d669864f --- /dev/null +++ b/src/SoilMatricPotentialMod.F90 @@ -0,0 +1,58 @@ +module SoilMatricPotentialMod + +!!! compute soil matric potential that will be used for +!!! soil moisture transpiration factor and dynamic root calculations + + use Machine + use NoahmpVarType + + implicit none + +contains + + subroutine SoilMatricPotential(noahmp) + +! ------------------------ Code history ----------------------------------- +! Original Noah-MP subroutine: None (embedded in ENERGY subroutine) +! Original code: Guo-Yue Niu and Noah-MP team (Niu et al. 2011) +! Refactered code: C. He, P. Valayamkunnath, & refactor team (He et al. 2023) +! ------------------------------------------------------------------------- + + implicit none + +! in & out variables + type(noahmp_type), intent(inout) :: noahmp + +! local variables + integer :: IndSoil ! loop index + +! -------------------------------------------------------------------- + associate( & + SurfaceType => noahmp%config%domain%SurfaceType ,& ! in, surface type 1-soil; 2-lake + NumSoilLayer => noahmp%config%domain%NumSoilLayer ,& ! in, number of soil layers + SoilMatPotentialWilt => noahmp%water%param%SoilMatPotentialWilt ,& ! in, soil metric potential for wilting point [m] + SoilMatPotentialSat => noahmp%water%param%SoilMatPotentialSat ,& ! in, saturated soil matric potential [m] + SoilMoistureSat => noahmp%water%param%SoilMoistureSat ,& ! in, saturated value of soil moisture [m3/m3] + SoilExpCoeffB => noahmp%water%param%SoilExpCoeffB ,& ! in, soil B parameter + SoilLiqWater => noahmp%water%state%SoilLiqWater ,& ! in, soil water content [m3/m3] + SoilMatPotential => noahmp%water%state%SoilMatPotential & ! out, soil matrix potential [m] + ) +! ---------------------------------------------------------------------- + + ! only for soil point + if ( SurfaceType ==1 ) then + + do IndSoil = 1, NumSoilLayer + SoilMatPotential(IndSoil) = max(SoilMatPotentialWilt, -SoilMatPotentialSat(IndSoil) * & + (max(0.01,SoilLiqWater(IndSoil))/SoilMoistureSat(IndSoil)) ** & + (-SoilExpCoeffB(IndSoil))) + enddo + + endif + + end associate + + end subroutine SoilMatricPotential + +end module SoilMatricPotentialMod + diff --git a/src/SoilWaterTranspirationMod.F90 b/src/SoilWaterTranspirationMod.F90 index d5ef583a..bf0fcf14 100644 --- a/src/SoilWaterTranspirationMod.F90 +++ b/src/SoilWaterTranspirationMod.F90 @@ -43,9 +43,9 @@ subroutine SoilWaterTranspiration(noahmp) SoilMoistureSat => noahmp%water%param%SoilMoistureSat ,& ! in, saturated value of soil moisture [m3/m3] SoilExpCoeffB => noahmp%water%param%SoilExpCoeffB ,& ! in, soil B parameter SoilLiqWater => noahmp%water%state%SoilLiqWater ,& ! in, soil water content [m3/m3] + SoilMatPotential => noahmp%water%state%SoilMatPotential ,& ! in, soil matrix potential [m] SoilTranspFac => noahmp%water%state%SoilTranspFac ,& ! out, soil water transpiration factor (0 to 1) - SoilTranspFacAcc => noahmp%water%state%SoilTranspFacAcc ,& ! out, accumulated soil water transpiration factor (0 to 1) - SoilMatPotential => noahmp%water%state%SoilMatPotential & ! out, soil matrix potential [m] + SoilTranspFacAcc => noahmp%water%state%SoilTranspFacAcc & ! out, accumulated soil water transpiration factor (0 to 1) ) ! ---------------------------------------------------------------------- @@ -61,16 +61,16 @@ subroutine SoilWaterTranspiration(noahmp) (SoilMoistureFieldCap(IndSoil) - SoilMoistureWilt(IndSoil)) endif if ( OptSoilWaterTranspiration == 2 ) then ! CLM - SoilMatPotential(IndSoil) = max(SoilMatPotentialWilt, -SoilMatPotentialSat(IndSoil) * & - (max(0.01,SoilLiqWater(IndSoil))/SoilMoistureSat(IndSoil)) ** & - (-SoilExpCoeffB(IndSoil))) + !SoilMatPotential(IndSoil) = max(SoilMatPotentialWilt, -SoilMatPotentialSat(IndSoil) * & + ! (max(0.01,SoilLiqWater(IndSoil))/SoilMoistureSat(IndSoil)) ** & + ! (-SoilExpCoeffB(IndSoil))) SoilWetFac = (1.0 - SoilMatPotential(IndSoil)/SoilMatPotentialWilt) / & (1.0 + SoilMatPotentialSat(IndSoil)/SoilMatPotentialWilt) endif if ( OptSoilWaterTranspiration == 3 ) then ! SSiB - SoilMatPotential(IndSoil) = max(SoilMatPotentialWilt, -SoilMatPotentialSat(IndSoil) * & - (max(0.01,SoilLiqWater(IndSoil))/SoilMoistureSat(IndSoil)) ** & - (-SoilExpCoeffB(IndSoil))) + !SoilMatPotential(IndSoil) = max(SoilMatPotentialWilt, -SoilMatPotentialSat(IndSoil) * & + ! (max(0.01,SoilLiqWater(IndSoil))/SoilMoistureSat(IndSoil)) ** & + ! (-SoilExpCoeffB(IndSoil))) SoilWetFac = 1.0 - exp(-5.8*(log(SoilMatPotentialWilt/SoilMatPotential(IndSoil)))) endif SoilWetFac = min(1.0, max(0.0,SoilWetFac)) diff --git a/src/WaterMainMod.F90 b/src/WaterMainMod.F90 index 988961f5..970f9ce9 100644 --- a/src/WaterMainMod.F90 +++ b/src/WaterMainMod.F90 @@ -43,6 +43,7 @@ subroutine WaterMain(noahmp) FlagSoilProcess => noahmp%config%domain%FlagSoilProcess ,& ! in, flag to calculate soil processes NumSoilTimeStep => noahmp%config%domain%NumSoilTimeStep ,& ! in, number of timesteps for soil process calculation OptWetlandModel => noahmp%config%nmlist%OptWetlandModel ,& ! in, options for wetland model + OptRoot => noahmp%config%nmlist%OptRoot ,& ! in, root option VaporizeGrd => noahmp%water%flux%VaporizeGrd ,& ! in, ground vaporize rate total (evap+sublim) [mm/s] CondenseVapGrd => noahmp%water%flux%CondenseVapGrd ,& ! in, ground vapor condense rate total (dew+frost) [mm/s] RainfallGround => noahmp%water%flux%RainfallGround ,& ! in, ground surface rain rate [mm/s] @@ -55,6 +56,7 @@ subroutine WaterMain(noahmp) ExchCoeffShSfc => noahmp%energy%state%ExchCoeffShSfc ,& ! in, exchange coefficient [m/s] for heat, surface, grid mean SpecHumidityRefHeight => noahmp%forcing%SpecHumidityRefHeight ,& ! in, specific humidity [kg/kg] at reference height HeatLatentGrd => noahmp%energy%flux%HeatLatentGrd ,& ! in, total ground latent heat [W/m2] (+ to atm) + RootActivity => noahmp%water%state%RootActivity ,& ! in, root activity function NumSnowLayerNeg => noahmp%config%domain%NumSnowLayerNeg ,& ! inout, actual number of snow layers (negative) ThicknessSnowSoilLayer => noahmp%config%domain%ThicknessSnowSoilLayer ,& ! inout, thickness of snow/soil layers [m] SnowWaterEquiv => noahmp%water%state%SnowWaterEquiv ,& ! inout, snow water equivalent [mm] @@ -149,9 +151,15 @@ subroutine WaterMain(noahmp) EvapSoilSfcLiq = EvapSoilSfcLiq * 0.001 ! mm/s -> m/s ! transpiration mm/s -> m/s - do LoopInd = 1, NumSoilLayerRoot - TranspWatLossSoil(LoopInd) = Transpiration * SoilTranspFac(LoopInd) * 0.001 - enddo + if ( OptRoot == 1 ) then + do LoopInd = 1, NumSoilLayerRoot + TranspWatLossSoil(LoopInd) = Transpiration * RootActivity(LoopInd) * 0.001 + enddo + else + do LoopInd = 1, NumSoilLayerRoot + TranspWatLossSoil(LoopInd) = Transpiration * SoilTranspFac(LoopInd) * 0.001 + enddo + endif ! total surface input water to soil mm/s -> m/s SoilSfcInflow = (PondSfcThinSnwMelt + PondSfcThinSnwComb + PondSfcThinSnwTrans) / & diff --git a/src/WaterVarType.F90 b/src/WaterVarType.F90 index 1bd80c45..c9a895fa 100644 --- a/src/WaterVarType.F90 +++ b/src/WaterVarType.F90 @@ -139,6 +139,7 @@ module WaterVarType real(kind=kind_noahmp) :: WaterBalanceError ! water balance error [mm] real(kind=kind_noahmp) :: WaterStorageTotEnd ! total water storage [mm] at the end of NoahMP process real(kind=kind_noahmp) :: SnowRadiusFresh ! fresh snow radius [microns] + real(kind=kind_noahmp) :: MaxUptakeDepth ! bottom of deepest root water uptake layer [m] integer , allocatable, dimension(:) :: IndexPhaseChange ! phase change index (0-none;1-melt;2-refreeze) real(kind=kind_noahmp), allocatable, dimension(:) :: SnowIce ! snow layer ice [mm] @@ -160,6 +161,9 @@ module WaterVarType real(kind=kind_noahmp), allocatable, dimension(:) :: SnowLiqWaterVol ! partial volume of snow liquid water [m3/m3] real(kind=kind_noahmp), allocatable, dimension(:) :: SoilSupercoolWater ! supercooled water in soil [kg/m2] real(kind=kind_noahmp), allocatable, dimension(:) :: SoilMatPotential ! soil matric potential [m] + real(kind=kind_noahmp), allocatable, dimension(:) :: EaseFunction ! ease function + real(kind=kind_noahmp), allocatable, dimension(:) :: RootActivity ! root activity function + real(kind=kind_noahmp), allocatable, dimension(:) :: InactiveTimeSteps ! number of time steps without active roots real(kind=kind_noahmp), allocatable, dimension(:) :: SnowRadius ! snow effective grain radius [microns, m-6] real(kind=kind_noahmp), allocatable, dimension(:) :: MassBChydropho ! mass of hydrophobic Black Carbon in snow [kg m-2] real(kind=kind_noahmp), allocatable, dimension(:) :: MassBChydrophi ! mass of hydrophillic Black Carbon in snow [kg m-2]