Addition of dynamic root scheme (Fan et al 2017; Bieri et al 2025)

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

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

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

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
 !------------------------------------------------------------------------

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

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

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) / &

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

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

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

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