From acf8056423a09fe84f8446a2c6c60f9a721dcaba Mon Sep 17 00:00:00 2001 From: Carolina Bieri Date: Wed, 25 Feb 2026 09:50:16 -0700 Subject: [PATCH 1/4] Addition of dynamic root scheme --- drivers/hrldas/ConfigVarInTransferMod.F90 | 1 + drivers/hrldas/NoahmpDriverMainMod.F90 | 2 +- drivers/hrldas/NoahmpIOVarInitMod.F90 | 16 +++ drivers/hrldas/NoahmpIOVarType.F90 | 9 ++ drivers/hrldas/NoahmpInitMainMod.F90 | 6 + drivers/hrldas/NoahmpReadNamelistMod.F90 | 4 +- drivers/hrldas/WaterVarInTransferMod.F90 | 6 + drivers/hrldas/WaterVarOutTransferMod.F90 | 8 ++ src/ConfigVarInitMod.F90 | 1 + src/ConfigVarType.F90 | 3 + src/DynaRootGrowDieMod.F90 | 156 ++++++++++++++++++++++ src/EnergyMainMod.F90 | 4 + src/Makefile | 8 +- src/WaterMainMod.F90 | 14 +- src/WaterVarType.F90 | 4 + 15 files changed, 234 insertions(+), 8 deletions(-) create mode 100644 src/DynaRootGrowDieMod.F90 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 c3f62825..b5042a15 100644 --- a/drivers/hrldas/NoahmpDriverMainMod.F90 +++ b/drivers/hrldas/NoahmpDriverMainMod.F90 @@ -113,7 +113,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 1f19d86e..3ff2dc6e 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%DRDEPTHXY) ) allocate ( NoahmpIO%DRDEPTHXY (XSTART:XEND, YSTART:YEND) ) ! root layer depth (-) + endif + ! Single- and Multi-layer Urban Models if ( NoahmpIO%SF_URBAN_PHYSICS > 0 ) then if ( .not. allocated (NoahmpIO%sh_urb2d) ) allocate ( NoahmpIO%sh_urb2d (XSTART:XEND,YSTART:YEND) ) @@ -889,6 +897,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%DRDEPTHXY = 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..5bc060bc 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(:,:) :: DRDEPTHXY ! root depth layer (-) + !------------------------------------------------------------------------ ! Single- and Multi-layer Urban Models !------------------------------------------------------------------------ diff --git a/drivers/hrldas/NoahmpInitMainMod.F90 b/drivers/hrldas/NoahmpInitMainMod.F90 index 8c6a9956..ffbe0620 100644 --- a/drivers/hrldas/NoahmpInitMainMod.F90 +++ b/drivers/hrldas/NoahmpInitMainMod.F90 @@ -269,6 +269,12 @@ subroutine NoahmpInitMain(NoahmpIO) NoahmpIO%STEPWTD = max(NoahmpIO%STEPWTD,1) endif + ! initialize time step counter for DynaRoot, if activated + if ( NoahmpIO%IOPT_ROOT == 1 ) then + NoahmpIO%INACTIVEXY(:,1,:) = 0.0 + NoahmpIO%INACTIVEXY(:,2:NoahmpIO%NSOIL,:) = 366.0*(86400.0 / NoahmpIO%DTBL) + endif + endif ! NoahmpIO%restart_flag if ( NoahmpIO%IOPT_ALB == 3 ) then ! initialize SNICAR aerosol content in snow 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..25109089 100644 --- a/drivers/hrldas/WaterVarInTransferMod.F90 +++ b/drivers/hrldas/WaterVarInTransferMod.F90 @@ -81,6 +81,12 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO) noahmp%water%state%SoilSaturateFrac = NoahmpIO%FSATXY (I,J) noahmp%water%state%WaterStorageWetland = NoahmpIO%WSURFXY (I,J) endif + 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%DynamicRootDepth = NoahmpIO%DRDEPTHXY (I,J) + endif #ifdef WRF_HYDRO noahmp%water%state%WaterTableHydro = NoahmpIO%ZWATBLE2D (I,J) noahmp%water%state%WaterHeadSfc = NoahmpIO%sfcheadrt (I,J) diff --git a/drivers/hrldas/WaterVarOutTransferMod.F90 b/drivers/hrldas/WaterVarOutTransferMod.F90 index b404bd31..45a5ca6f 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%DRDEPTHXY (I,J) = noahmp%water%state%DynamicRootDepth + 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..a3117998 --- /dev/null +++ b/src/DynaRootGrowDieMod.F90 @@ -0,0 +1,156 @@ +module DynaRootGrowDieMod + +!!! compute soil water transpiration factor that will be used for +!!! stomata resistance and evapotranspiration calculations + + 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 ! bottom layer of root zone (root depth) + 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] + + integer, 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 + ) +! ---------------------------------------------------------------------- + + 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 + InactiveDays = undefined_real + MaximumInactiveDays = 31536000.0 ! seconds in a year + EaseFunction = 0.0 + RootActivity = 0.0 + InactiveDays = 0.0 + + 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 = MIN(MAX(IndSoil,1),NumSoilLayer) + + 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 + + !INACTIVEDAYS = MIN(INACTIVEDAYS,MaximumInactiveDays+1) + + NumSoilLayerRoot = DynaRootLayer + + end associate + + end subroutine DynaRootGrowDie + +end module DynaRootGrowDieMod diff --git a/src/EnergyMainMod.F90 b/src/EnergyMainMod.F90 index 093fecd7..adb7560c 100644 --- a/src/EnergyMainMod.F90 +++ b/src/EnergyMainMod.F90 @@ -53,6 +53,7 @@ module EnergyMainMod use SurfaceEnergyFluxBareGroundMod, only : SurfaceEnergyFluxBareGround use SoilSnowTemperatureMainMod, only : SoilSnowTemperatureMain use SoilSnowWaterPhaseChangeMod, only : SoilSnowWaterPhaseChange + use DynaRootGrowDieMod, only : DynaRootGrowDie implicit none @@ -80,6 +81,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 +216,8 @@ subroutine EnergyMain(noahmp) ! longwave emissivity for vegetation, ground, total net surface call SurfaceEmissivity(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 880039e8..ae6f83ba 100644 --- a/src/Makefile +++ b/src/Makefile @@ -140,7 +140,8 @@ OBJS = ConstantDefineMod.o \ SurfaceEmissivityGlacierMod.o \ SurfaceEnergyFluxGlacierMod.o \ SurfaceRadiationGlacierMod.o \ - WaterMainGlacierMod.o + WaterMainGlacierMod.o \ + DynaRootGrowDieMod.o all: $(OBJS) @@ -287,7 +288,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 NoahmpMainMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \ IrrigationPrepareMod.o IrrigationSprinklerMod.o CanopyWaterInterceptMod.o \ PrecipitationHeatAdvectMod.o EnergyMainMod.o WaterMainMod.o AtmosForcingMod.o \ @@ -375,4 +377,4 @@ 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 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..73a95bb9 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) :: DynamicRootDepth ! root depth layer 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] From dde04e4043aea38869dc2224ba5cb8c93aa75e02 Mon Sep 17 00:00:00 2001 From: Carolina Bieri Date: Fri, 6 Mar 2026 13:01:49 -0700 Subject: [PATCH 2/4] Bug fixes for root scheme and addition of SoilMatricPotentialMod --- drivers/hrldas/NoahmpInitMainMod.F90 | 5 +- drivers/hrldas/WaterVarInTransferMod.F90 | 13 +++--- src/DynaRootGrowDieMod.F90 | 10 ++-- src/EnergyMainMod.F90 | 3 ++ src/Makefile | 8 ++-- src/SoilMatricPotentialMod.F90 | 58 ++++++++++++++++++++++++ src/SoilWaterTranspirationMod.F90 | 16 +++---- 7 files changed, 92 insertions(+), 21 deletions(-) create mode 100644 src/SoilMatricPotentialMod.F90 diff --git a/drivers/hrldas/NoahmpInitMainMod.F90 b/drivers/hrldas/NoahmpInitMainMod.F90 index ae1e7e31..d3e8d60f 100644 --- a/drivers/hrldas/NoahmpInitMainMod.F90 +++ b/drivers/hrldas/NoahmpInitMainMod.F90 @@ -269,10 +269,13 @@ subroutine NoahmpInitMain(NoahmpIO) NoahmpIO%STEPWTD = max(NoahmpIO%STEPWTD,1) endif - ! initialize time step counter for DynaRoot, if activated + ! 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 diff --git a/drivers/hrldas/WaterVarInTransferMod.F90 b/drivers/hrldas/WaterVarInTransferMod.F90 index 25109089..77ddddaf 100644 --- a/drivers/hrldas/WaterVarInTransferMod.F90 +++ b/drivers/hrldas/WaterVarInTransferMod.F90 @@ -81,12 +81,6 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO) noahmp%water%state%SoilSaturateFrac = NoahmpIO%FSATXY (I,J) noahmp%water%state%WaterStorageWetland = NoahmpIO%WSURFXY (I,J) endif - 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%DynamicRootDepth = NoahmpIO%DRDEPTHXY (I,J) - endif #ifdef WRF_HYDRO noahmp%water%state%WaterTableHydro = NoahmpIO%ZWATBLE2D (I,J) noahmp%water%state%WaterHeadSfc = NoahmpIO%sfcheadrt (I,J) @@ -291,6 +285,13 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO) noahmp%water%param%SoilMoistureDry = 0.40 endif + 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%DynamicRootDepth = NoahmpIO%DRDEPTHXY (I,J) + endif + if ( SoilType(1) /= 14 ) then noahmp%water%param%SoilImpervFracCoeff = noahmp%water%param%GroundFrzCoeff * & ((noahmp%water%param%SoilMoistureSat(1) / & diff --git a/src/DynaRootGrowDieMod.F90 b/src/DynaRootGrowDieMod.F90 index a3117998..a85ac7b7 100644 --- a/src/DynaRootGrowDieMod.F90 +++ b/src/DynaRootGrowDieMod.F90 @@ -34,7 +34,7 @@ subroutine DynaRootGrowDie(noahmp) 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] - integer, allocatable, dimension(:) :: DynaRootMask ! indicator for presence of active roots (1=roots active; 0=roots not active) + 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] ! -------------------------------------------------------------------- @@ -61,7 +61,6 @@ subroutine DynaRootGrowDie(noahmp) ! initialize DynaRootMask = 0.0 ThicknessSoilLayerTmp = undefined_real - InactiveDays = undefined_real MaximumInactiveDays = 31536000.0 ! seconds in a year EaseFunction = 0.0 RootActivity = 0.0 @@ -89,6 +88,7 @@ subroutine DynaRootGrowDie(noahmp) IF(InactiveDays(IndSoil) .le. MaximumInactiveDays)EXIT END DO DynaRootLayer = MIN(MAX(IndSoil,1),NumSoilLayer) + print *, "DynaRootLayer:", DynaRootLayer EffectiveCanopyHeight = HeightCanopyTop*(2./3.) ! calculation of ease function and root activity @@ -109,7 +109,11 @@ subroutine DynaRootGrowDie(noahmp) ENDIF ! calculate ease function - EaseFunction(IndSoil) = MAX(-( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ), 0.) + ! EaseFunction(IndSoil) = MAX(-( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ), 0.) + EaseFunction(IndSoil) = -( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ) + print *, "Ease function:", EaseFunction(IndSoil) + print *, "Matric potential:", SoilMatPotential(IndSoil) + END DO ! to grow roots anew, the layer has to be easiest to get water from than the diff --git a/src/EnergyMainMod.F90 b/src/EnergyMainMod.F90 index adb7560c..6e7ec73a 100644 --- a/src/EnergyMainMod.F90 +++ b/src/EnergyMainMod.F90 @@ -54,6 +54,7 @@ module EnergyMainMod use SoilSnowTemperatureMainMod, only : SoilSnowTemperatureMain use SoilSnowWaterPhaseChangeMod, only : SoilSnowWaterPhaseChange use DynaRootGrowDieMod, only : DynaRootGrowDie + use SoilMatricPotentialMod, only : SoilMatricPotential implicit none @@ -216,6 +217,8 @@ 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 diff --git a/src/Makefile b/src/Makefile index c850a22c..d5ff4142 100644 --- a/src/Makefile +++ b/src/Makefile @@ -140,7 +140,8 @@ OBJS = ConstantDefineMod.o \ SurfaceEnergyFluxGlacierMod.o \ SurfaceRadiationGlacierMod.o \ WaterMainGlacierMod.o \ - DynaRootGrowDieMod.o + DynaRootGrowDieMod.o \ + SoilMatricPotentialMod.o all: $(OBJS) @@ -280,7 +281,7 @@ EnergyMainMod.o: ../utility/Machine.o NoahmpVarType.o Const SoilSnowWaterPhaseChangeMod.o SnowCoverGroundNiu07Mod.o SurfaceEmissivityMod.o \ GroundRoughnessPropertyMod.o PsychrometricVariableMod.o ResistanceGroundEvaporationMod.o \ SoilWaterTranspirationMod.o SurfaceAlbedoMod.o SurfaceRadiationMod.o SnowCoverGroundAR25Mod.o \ - DynaRootGrowDieMod.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 \ @@ -367,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 +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)) From ce461ae9bda7a2c8fac439579c437f4782863bd3 Mon Sep 17 00:00:00 2001 From: Carolina Bieri Date: Mon, 9 Mar 2026 12:02:47 -0600 Subject: [PATCH 3/4] Remove debugging print statements --- src/DynaRootGrowDieMod.F90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/DynaRootGrowDieMod.F90 b/src/DynaRootGrowDieMod.F90 index a85ac7b7..5a1349b5 100644 --- a/src/DynaRootGrowDieMod.F90 +++ b/src/DynaRootGrowDieMod.F90 @@ -1,7 +1,7 @@ module DynaRootGrowDieMod -!!! compute soil water transpiration factor that will be used for -!!! stomata resistance and evapotranspiration calculations +!!! 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 @@ -88,7 +88,6 @@ subroutine DynaRootGrowDie(noahmp) IF(InactiveDays(IndSoil) .le. MaximumInactiveDays)EXIT END DO DynaRootLayer = MIN(MAX(IndSoil,1),NumSoilLayer) - print *, "DynaRootLayer:", DynaRootLayer EffectiveCanopyHeight = HeightCanopyTop*(2./3.) ! calculation of ease function and root activity @@ -109,10 +108,7 @@ subroutine DynaRootGrowDie(noahmp) ENDIF ! calculate ease function - ! EaseFunction(IndSoil) = MAX(-( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ), 0.) - EaseFunction(IndSoil) = -( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ) - print *, "Ease function:", EaseFunction(IndSoil) - print *, "Matric potential:", SoilMatPotential(IndSoil) + EaseFunction(IndSoil) = MAX(-( LeafWiltMatPotential - SoilMatPotential(IndSoil) )*SoilIceFactor / ( EffectiveCanopyHeight-LayerMidpoint ), 0.) END DO From 1bd14e042aa90d6bed70f87b8bac4988c68cd4f0 Mon Sep 17 00:00:00 2001 From: Carolina Bieri Date: Tue, 10 Mar 2026 09:22:58 -0600 Subject: [PATCH 4/4] Added MaxUptakeDepth calculation and updated variable names --- drivers/hrldas/NoahmpIOVarInitMod.F90 | 4 ++-- drivers/hrldas/NoahmpIOVarType.F90 | 2 +- drivers/hrldas/WaterVarInTransferMod.F90 | 3 ++- drivers/hrldas/WaterVarOutTransferMod.F90 | 2 +- src/DynaRootGrowDieMod.F90 | 21 +++++++++++++++------ src/WaterVarType.F90 | 2 +- 6 files changed, 22 insertions(+), 12 deletions(-) diff --git a/drivers/hrldas/NoahmpIOVarInitMod.F90 b/drivers/hrldas/NoahmpIOVarInitMod.F90 index d8d0154f..90dbb1b1 100644 --- a/drivers/hrldas/NoahmpIOVarInitMod.F90 +++ b/drivers/hrldas/NoahmpIOVarInitMod.F90 @@ -429,7 +429,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) 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%DRDEPTHXY) ) allocate ( NoahmpIO%DRDEPTHXY (XSTART:XEND, YSTART:YEND) ) ! root layer depth (-) + if ( .not. allocated (NoahmpIO%MAXUPTAKEXY) ) allocate ( NoahmpIO%MAXUPTAKEXY (XSTART:XEND, YSTART:YEND) ) ! root layer depth (-) endif ! Single- and Multi-layer Urban Models @@ -904,7 +904,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) NoahmpIO%EASYXY = 0.0 NoahmpIO%ROOTACTIVITYXY = 0.0 NoahmpIO%INACTIVEXY = undefined_real - NoahmpIO%DRDEPTHXY = undefined_real + NoahmpIO%MAXUPTAKEXY = undefined_real endif ! spatial varying soil texture diff --git a/drivers/hrldas/NoahmpIOVarType.F90 b/drivers/hrldas/NoahmpIOVarType.F90 index 5bc060bc..25f2425c 100644 --- a/drivers/hrldas/NoahmpIOVarType.F90 +++ b/drivers/hrldas/NoahmpIOVarType.F90 @@ -508,7 +508,7 @@ module NoahmpIOVarType 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(:,:) :: DRDEPTHXY ! root depth layer (-) + 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/WaterVarInTransferMod.F90 b/drivers/hrldas/WaterVarInTransferMod.F90 index 77ddddaf..92f7baaa 100644 --- a/drivers/hrldas/WaterVarInTransferMod.F90 +++ b/drivers/hrldas/WaterVarInTransferMod.F90 @@ -285,11 +285,12 @@ 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%DynamicRootDepth = NoahmpIO%DRDEPTHXY (I,J) + noahmp%water%state%MaxUptakeDepth = NoahmpIO%MAXUPTAKEXY (I,J) endif if ( SoilType(1) /= 14 ) then diff --git a/drivers/hrldas/WaterVarOutTransferMod.F90 b/drivers/hrldas/WaterVarOutTransferMod.F90 index 45a5ca6f..c9748b4d 100644 --- a/drivers/hrldas/WaterVarOutTransferMod.F90 +++ b/drivers/hrldas/WaterVarOutTransferMod.F90 @@ -175,7 +175,7 @@ subroutine WaterVarOutTransfer(noahmp, NoahmpIO) 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%DRDEPTHXY (I,J) = noahmp%water%state%DynamicRootDepth + NoahmpIO%MAXUPTAKEXY (I,J) = noahmp%water%state%MaxUptakeDepth endif #ifdef WRF_HYDRO diff --git a/src/DynaRootGrowDieMod.F90 b/src/DynaRootGrowDieMod.F90 index 5a1349b5..ecfbacaf 100644 --- a/src/DynaRootGrowDieMod.F90 +++ b/src/DynaRootGrowDieMod.F90 @@ -25,7 +25,8 @@ subroutine DynaRootGrowDie(noahmp) integer :: IndSoil ! loop index integer :: WaterTableLayer ! layer that contains water table integer :: AboveWaterTableLayer ! layer above water table - integer :: DynaRootLayer ! bottom layer of root zone (root depth) + 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] @@ -50,7 +51,8 @@ subroutine DynaRootGrowDie(noahmp) 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 + RootActivity => noahmp%water%state%RootActivity ,& ! out, root activity function + MaxUptakeDepth => noahmp%water%state%MaxUptakeDepth & ! out, bottom of deepest root water uptake layer [m] ) ! ---------------------------------------------------------------------- @@ -61,11 +63,12 @@ subroutine DynaRootGrowDie(noahmp) ! initialize DynaRootMask = 0.0 ThicknessSoilLayerTmp = undefined_real - MaximumInactiveDays = 31536000.0 ! seconds in a year 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 @@ -87,7 +90,7 @@ subroutine DynaRootGrowDie(noahmp) DO IndSoil = NumSoilLayer,1,-1 IF(InactiveDays(IndSoil) .le. MaximumInactiveDays)EXIT END DO - DynaRootLayer = MIN(MAX(IndSoil,1),NumSoilLayer) + DynaRootLayer = MAX(IndSoil,1) EffectiveCanopyHeight = HeightCanopyTop*(2./3.) ! calculation of ease function and root activity @@ -145,10 +148,16 @@ subroutine DynaRootGrowDie(noahmp) ENDIF ENDDO - !INACTIVEDAYS = MIN(INACTIVEDAYS,MaximumInactiveDays+1) - + ! 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 diff --git a/src/WaterVarType.F90 b/src/WaterVarType.F90 index 73a95bb9..c9a895fa 100644 --- a/src/WaterVarType.F90 +++ b/src/WaterVarType.F90 @@ -139,7 +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) :: DynamicRootDepth ! root depth layer + 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]