Add wetland scheme of Zhang et al 2022 WRR

drivers/hrldas/ConfigVarInTransferMod.F90

@@ -60,6 +60,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     noahmp%config%nmlist%OptRunoffSubsurface         = NoahmpIO%IOPT_RUNSUB
     noahmp%config%nmlist%OptGlacierTreatment         = NoahmpIO%IOPT_GLA
     noahmp%config%nmlist%OptSnowCompaction           = NoahmpIO%IOPT_COMPACT
+    noahmp%config%nmlist%OptWetlandModel             = NoahmpIO%IOPT_WETLAND
 
     ! config domain variable
     noahmp%config%domain%SurfaceType                 = 1

drivers/hrldas/NoahmpDriverMainMod.F90

@@ -70,6 +70,7 @@ subroutine NoahmpDriverMain(NoahmpIO)
        NoahmpIO%ACC_ECANXY   = 0.0
        NoahmpIO%ACC_ETRANXY  = 0.0
        NoahmpIO%ACC_EDIRXY   = 0.0
+       NoahmpIO%ACC_GLAFLWXY = 0.0
     endif
 
     if ( NoahmpIO%SOIL_UPDATE_STEPS > 1 ) then
@@ -83,6 +84,7 @@ subroutine NoahmpDriverMain(NoahmpIO)
           NoahmpIO%ACC_ECANXY   = 0.0
           NoahmpIO%ACC_ETRANXY  = 0.0
           NoahmpIO%ACC_EDIRXY   = 0.0
+          NoahmpIO%ACC_GLAFLWXY = 0.0
        end if
     endif
 

drivers/hrldas/NoahmpIOVarInitMod.F90

@@ -297,6 +297,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     if ( .not. allocated (NoahmpIO%ACC_QINSURXY)) allocate ( NoahmpIO%ACC_QINSURXY (XSTART:XEND,        YSTART:YEND) )
     if ( .not. allocated (NoahmpIO%ACC_QSEVAXY) ) allocate ( NoahmpIO%ACC_QSEVAXY  (XSTART:XEND,        YSTART:YEND) )
     if ( .not. allocated (NoahmpIO%ACC_ETRANIXY)) allocate ( NoahmpIO%ACC_ETRANIXY (XSTART:XEND,1:NSOIL,YSTART:YEND) )
+    if ( .not. allocated (NoahmpIO%ACC_GLAFLWXY)) allocate ( NoahmpIO%ACC_GLAFLWXY (XSTART:XEND,        YSTART:YEND) )
 
     ! Needed for MMF_RUNOFF (IOPT_RUN = 5); not part of MP driver in WRF
     if ( .not. allocated (NoahmpIO%MSFTX)      ) allocate ( NoahmpIO%MSFTX       (XSTART:XEND,YSTART:YEND) ) 
@@ -325,6 +326,16 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     if ( .not. allocated (NoahmpIO%SEASON_GDD)) allocate ( NoahmpIO%SEASON_GDD (XSTART:XEND,  YSTART:YEND) )
     if ( .not. allocated (NoahmpIO%CROPTYPE)  ) allocate ( NoahmpIO%CROPTYPE   (XSTART:XEND,5,YSTART:YEND) )
 
+    ! Needed for Zhang et al. 2022 wetland model (OPT_WETLAND=1 or 2)
+    if ( NoahmpIO%IOPT_WETLAND > 0 ) then
+       if ( .not. allocated (NoahmpIO%FSATXY) ) allocate ( NoahmpIO%FSATXY     (XSTART:XEND,  YSTART:YEND) ) ! saturated fraction of the grid (-)
+       if ( .not. allocated (NoahmpIO%WSURFXY)) allocate ( NoahmpIO%WSURFXY    (XSTART:XEND,  YSTART:YEND) ) ! wetland water storage [mm]
+    endif
+    if ( NoahmpIO%IOPT_WETLAND == 2 ) then
+       if ( .not. allocated (NoahmpIO%FSATMX) ) allocate ( NoahmpIO%FSATMX     (XSTART:XEND,  YSTART:YEND) ) ! maximum saturated fraction
+       if ( .not. allocated (NoahmpIO%WCAP)   ) allocate ( NoahmpIO%WCAP       (XSTART:XEND,  YSTART:YEND) ) ! maximum wetland capacity [m]
+    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) ) 
@@ -339,7 +350,6 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
        if ( .not. allocated (NoahmpIO%lb_urb2d)   ) allocate ( NoahmpIO%lb_urb2d    (XSTART:XEND,YSTART:YEND) )
        if ( .not. allocated (NoahmpIO%hgt_urb2d)  ) allocate ( NoahmpIO%hgt_urb2d   (XSTART:XEND,YSTART:YEND) )
        if ( .not. allocated (NoahmpIO%ust)        ) allocate ( NoahmpIO%ust         (XSTART:XEND,YSTART:YEND) )
-       !ENDIF
 
        !IF(NoahmpIO%SF_URBAN_PHYSICS == 1 ) THEN  ! single layer urban model  
        if ( .not. allocated (NoahmpIO%cmr_sfcdif)   ) allocate ( NoahmpIO%cmr_sfcdif    (XSTART:XEND,        YSTART:YEND) )
@@ -651,6 +661,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     NoahmpIO%ACC_ECANXY      = 0.0
     NoahmpIO%ACC_ETRANXY     = 0.0
     NoahmpIO%ACC_EDIRXY      = 0.0
+    NoahmpIO%ACC_GLAFLWXY    = 0.0
 
     ! MMF Groundwater
     NoahmpIO%TERRAIN         = undefined_real
@@ -701,6 +712,16 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     NoahmpIO%IRRSPLH         = 0.0
     NoahmpIO%LOCTIM          = undefined_real
 
+    ! wetland model (Zhang et al. 2022)
+    if ( NoahmpIO%IOPT_WETLAND > 0 ) then
+       NoahmpIO%FSATXY       = undefined_real
+       NoahmpIO%WSURFXY      = undefined_real
+    endif
+    if ( NoahmpIO%IOPT_WETLAND == 2 ) then
+       NoahmpIO%FSATMX       = undefined_real
+       NoahmpIO%WCAP         = undefined_real
+    endif
+
     ! spatial varying soil texture
     if ( NoahmpIO%IOPT_SOIL > 1 ) then
        NoahmpIO%SOILCL1      = undefined_real

drivers/hrldas/NoahmpIOVarType.F90

@@ -64,6 +64,7 @@ module NoahmpIOVarType
     integer                                                ::  IOPT_INFDV          ! infiltration options for dynamic VIC (1->Philip; 2-> Green-Ampt;3->Smith-Parlange)
     integer                                                ::  IOPT_TDRN           ! drainage option (0->off; 1->simple scheme; 2->Hooghoudt's scheme)
     integer                                                ::  IOPT_COMPACT        ! snowpack compaction (1->Anderson1976; 2->Abolafia-Rosenzweig2024)
+    integer                                                ::  IOPT_WETLAND        ! wetland model option (0->off; 1->Zhang2022 fixed parameter; 2->Zhang2022 read in 2D parameter)
     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]
@@ -331,6 +332,7 @@ module NoahmpIOVarType
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  ACC_ETRANXY         ! accumulated transpiration per soil timestep [mm]
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  ACC_EDIRXY          ! accumulated net ground (soil/snow) evaporation per soil timestep [mm]
     real(kind=kind_noahmp), allocatable, dimension(:,:,:)  ::  ACC_ETRANIXY        ! accumualted transpiration rate within soil timestep [m/s * dt_soil/dt_main]
+    real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  ACC_GLAFLWXY        ! accumulated glacier excessive flow [mm] per soil timestep
 
 !------------------------------------------------------------------------
 ! Needed for MMF_RUNOFF (IOPT_RUN = 5); not part of MP driver in WRF
@@ -372,6 +374,14 @@ module NoahmpIOVarType
     real(kind=kind_noahmp), allocatable, dimension(:,:)    :: SEASON_GDD           ! seasonal GDD
     real(kind=kind_noahmp), allocatable, dimension(:,:,:)  :: CROPTYPE             ! crop type
 
+!------------------------------------------------------------------------
+! Needed for wetland model (OPT_WETLAND=1 or 2)
+!------------------------------------------------------------------------
+    real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  FSATXY              ! saturated fraction of the grid (-)
+    real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  WSURFXY             ! wetland water storage [mm]
+    real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  FSATMX              ! maximum saturated fraction
+    real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  WCAP                ! maximum wetland capacity [m]
+
 !------------------------------------------------------------------------
 ! Single- and Multi-layer Urban Models
 !------------------------------------------------------------------------
@@ -878,6 +888,9 @@ module NoahmpIOVarType
     real(kind=kind_noahmp), allocatable, dimension(:,:)    :: RTCT_TABLE                ! fraction of carbohydrate translocation from root to grain
     real(kind=kind_noahmp), allocatable, dimension(:)      :: BIO2LAI_TABLE             ! leaf area per living leaf biomass [m2/kg]
 
+    ! wetland parameter (OPT_WETLAND=1)
+    real(kind=kind_noahmp)                                 :: WCAP_TABLE                ! maximum surface wetland capacity
+
     ! soil parameters
     integer                                                :: SLCATS_TABLE              ! number of soil categories
     real(kind=kind_noahmp), allocatable, dimension(:)      :: BEXP_TABLE                ! soil B parameter

drivers/hrldas/NoahmpInitMainMod.F90

@@ -151,6 +151,10 @@ subroutine NoahmpInitMain(NoahmpIO)
              NoahmpIO%QSNOWXY(I,J)  = 0.0
              NoahmpIO%QRAINXY(I,J)  = 0.0
              NoahmpIO%WSLAKEXY(I,J) = 0.0
+             if ( NoahmpIO%IOPT_WETLAND > 0 ) then
+                NoahmpIO%FSATXY(I,J)   = 0.0
+                NoahmpIO%WSURFXY(I,J)  = 0.0
+             endif
              if ( NoahmpIO%IOPT_RUNSUB /= 5 ) then 
                 NoahmpIO%WAXY(I,J)   = 4900.0 
                 NoahmpIO%WTXY(I,J)   = NoahmpIO%WAXY(i,j) 

drivers/hrldas/NoahmpReadNamelistMod.F90

@@ -93,6 +93,7 @@ subroutine NoahmpReadNamelist(NoahmpIO)
     integer                 :: irrigation_method                  = 0
     integer                 :: dvic_infiltration_option           = 1
     integer                 :: tile_drainage_option               = 0
+    integer                 :: wetland_option                     = 0 
     integer                 :: split_output_count                 = 1
     logical                 :: skip_first_output                  = .false.
     integer                 :: khour                              = -9999
@@ -126,7 +127,7 @@ subroutine NoahmpReadNamelist(NoahmpIO)
          btr_option, surface_drag_option, supercooled_water_option,                       &
          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,                      &
+         pcp_partition_option, tbot_option, temp_time_scheme_option, wetland_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,         &
@@ -345,6 +346,7 @@ subroutine NoahmpReadNamelist(NoahmpIO)
     NoahmpIO%IOPT_INFDV                        = dvic_infiltration_option
     NoahmpIO%IOPT_TDRN                         = tile_drainage_option
     NoahmpIO%IOPT_COMPACT                      = snow_compaction_option
+    NoahmpIO%IOPT_WETLAND                      = wetland_option
     ! basic model setup variables
     NoahmpIO%indir                             = indir
     NoahmpIO%forcing_timestep                  = forcing_timestep

drivers/hrldas/NoahmpReadTableMod.F90

@@ -114,7 +114,7 @@ subroutine NoahmpReadTable(NoahmpIO)
                                                 SNOWCOMPACTm_AR24,SNOWCOMPACTb_AR24,SNOWCOMPACT_P1_AR24, SNOWCOMPACT_P2_AR24,&
                                                 SNOWCOMPACT_P3_AR24, SNOWCOMPACT_Up_AR24,                                    &
                                                 WSLMAX, ROUS, CMIC, SNOWDEN_MAX, CLASS_ALB_REF, CLASS_SNO_AGE, CLASS_ALB_NEW,&
-                                                PSIWLT, Z0SOIL, Z0LAKE
+                                                PSIWLT, Z0SOIL, Z0LAKE, WCAP
     namelist / noahmp_global_parameters /       CO2, O2, TIMEAN, FSATMX, Z0SNO, SSI, SNOW_RET_FAC ,SNOW_EMIS, SWEMX, TAU0,   &
                                                 GRAIN_GROWTH, EXTRA_GROWTH, DIRT_SOOT, BATS_COSZ, BATS_VIS_NEW,              &
                                                 BATS_NIR_NEW, BATS_VIS_AGE, BATS_NIR_AGE, BATS_VIS_DIR, BATS_NIR_DIR,        &
@@ -123,7 +123,7 @@ subroutine NoahmpReadTable(NoahmpIO)
                                                 SNOWCOMPACTm_AR24,SNOWCOMPACTb_AR24,SNOWCOMPACT_P1_AR24, SNOWCOMPACT_P2_AR24,&
                                                 SNOWCOMPACT_P3_AR24, SNOWCOMPACT_Up_AR24,                                    &
                                                 WSLMAX, ROUS, CMIC, SNOWDEN_MAX, CLASS_ALB_REF, CLASS_SNO_AGE, CLASS_ALB_NEW,&
-                                                PSIWLT, Z0SOIL, Z0LAKE
+                                                PSIWLT, Z0SOIL, Z0LAKE, WCAP
 
     ! irrigation parameters
     integer                                  :: IRR_HAR
@@ -531,6 +531,7 @@ subroutine NoahmpReadTable(NoahmpIO)
     NoahmpIO%PSIWLT_TABLE              = undefined_real
     NoahmpIO%Z0SOIL_TABLE              = undefined_real
     NoahmpIO%Z0LAKE_TABLE              = undefined_real
+    NoahmpIO%WCAP_TABLE                = undefined_real
 
     ! irrigation parameters
     NoahmpIO%IRR_HAR_TABLE          = undefined_int
@@ -933,6 +934,7 @@ subroutine NoahmpReadTable(NoahmpIO)
     NoahmpIO%PSIWLT_TABLE              = PSIWLT
     NoahmpIO%Z0SOIL_TABLE              = Z0SOIL
     NoahmpIO%Z0LAKE_TABLE              = Z0LAKE
+    NoahmpIO%WCAP_TABLE                = WCAP
 
     !---------------- NoahmpTable.TBL irrigation parameters
     inquire( file='NoahmpTable.TBL', exist=file_named )

drivers/hrldas/WaterVarInTransferMod.F90

@@ -77,6 +77,10 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO)
     noahmp%water%state%SoilMoistureEqui  (1:NumSoilLayer) = NoahmpIO%SMOISEQ    (I,1:NumSoilLayer,J)
     noahmp%water%state%RechargeGwDeepWT                   = 0.0
     noahmp%water%state%RechargeGwShallowWT                = 0.0
+    if ( noahmp%config%nmlist%OptWetlandModel > 0 ) then
+       noahmp%water%state%SoilSaturateFrac                = NoahmpIO%FSATXY     (I,J)
+       noahmp%water%state%WaterStorageWetland             = NoahmpIO%WSURFXY    (I,J)
+    endif
 #ifdef WRF_HYDRO
     noahmp%water%state%WaterTableHydro                    = NoahmpIO%ZWATBLE2D  (I,J)
     noahmp%water%state%WaterHeadSfc                       = NoahmpIO%sfcheadrt  (I,J)
@@ -91,6 +95,7 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO)
     noahmp%water%flux%TranspirationAcc                    = NoahmpIO%ACC_ETRANXY (I,J)
     noahmp%water%flux%EvapGroundNetAcc                    = NoahmpIO%ACC_EDIRXY  (I,J)
     noahmp%water%flux%TranspWatLossSoilAcc(1:NumSoilLayer)= NoahmpIO%ACC_ETRANIXY(I,1:NumSoilLayer,J)
+    noahmp%water%flux%GlacierExcessFlowAcc                = NoahmpIO%ACC_GLAFLWXY(I,J)
 
     ! water parameter variables
     noahmp%water%param%DrainSoilLayerInd                  = NoahmpIO%DRAIN_LAYER_OPT_TABLE
@@ -152,6 +157,7 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO)
     noahmp%water%param%DrainWatDepToImperv                = NoahmpIO%TD_D_TABLE(SoilType(1))
     noahmp%water%param%NumSoilLayerRoot                   = NoahmpIO%NROOT_TABLE(VegType)
     noahmp%water%param%SoilDrainSlope                     = NoahmpIO%SLOPE_TABLE(RunoffSlopeType)
+    noahmp%water%param%WetlandCapMax                      = NoahmpIO%WCAP_TABLE
 
     do IndexSoilLayer = 1, size(SoilType)
        noahmp%water%param%SoilMoistureSat       (IndexSoilLayer) = NoahmpIO%SMCMAX_TABLE(SoilType(IndexSoilLayer))
@@ -205,6 +211,12 @@ subroutine WaterVarInTransfer(noahmp, NoahmpIO)
        noahmp%water%param%DrainTubeDist          = NoahmpIO%TD_SPAC  (I,J)                ! tile spacing
     endif
 
+    ! spatial varying wetland parameters from input
+    if ( noahmp%config%nmlist%OptWetlandModel == 2 ) then
+       noahmp%water%param%SoilSfcSatFracMax      = NoahmpIO%FSATMX(I,J)
+       noahmp%water%param%WetlandCapMax          = NoahmpIO%WCAP(I,J)
+    endif 
+
     ! derived water parameters
     noahmp%water%param%SoilInfilMaxCoeff  = noahmp%water%param%SoilInfilFacRef *           &
                                             noahmp%water%param%SoilWatConductivitySat(1) / &

drivers/hrldas/WaterVarOutTransferMod.F90

@@ -120,6 +120,7 @@ subroutine WaterVarOutTransfer(noahmp, NoahmpIO)
     NoahmpIO%ACC_ECANXY  (I,J) = noahmp%water%flux%EvapCanopyNetAcc
     NoahmpIO%ACC_ETRANXY (I,J) = noahmp%water%flux%TranspirationAcc
     NoahmpIO%ACC_EDIRXY  (I,J) = noahmp%water%flux%EvapGroundNetAcc
+    NoahmpIO%ACC_GLAFLWXY(I,J) = noahmp%water%flux%GlacierExcessFlowAcc
     NoahmpIO%RECHXY      (I,J) = NoahmpIO%RECHXY(I,J) + (noahmp%water%state%RechargeGwShallowWT*1.0e3)
     NoahmpIO%DEEPRECHXY  (I,J) = NoahmpIO%DEEPRECHXY(I,J) + noahmp%water%state%RechargeGwDeepWT
     NoahmpIO%SMCWTDXY    (I,J) = noahmp%water%state%SoilMoistureToWT
@@ -141,6 +142,12 @@ subroutine WaterVarOutTransfer(noahmp, NoahmpIO)
     NoahmpIO%IRFIVOL   (I,J) = NoahmpIO%IRFIVOL(I,J)+(noahmp%water%flux%IrrigationRateFlood*1000.0)
     NoahmpIO%IRELOSS   (I,J) = NoahmpIO%IRELOSS(I,J)+(noahmp%water%flux%EvapIrriSprinkler*NoahmpIO%DTBL)
 
+    ! wetland (Zhang2022)
+    if ( noahmp%config%nmlist%OptWetlandModel > 0 ) then
+       NoahmpIO%WSURFXY(I,J) = noahmp%water%state%WaterStorageWetland
+       NoahmpIO%FSATXY (I,J) = noahmp%water%state%SoilSaturateFrac
+    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

drivers/lis/ConfigVarInTransferMod.F90

@@ -60,6 +60,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     noahmp%config%nmlist%OptRunoffSubsurface         = NoahmpIO%IOPT_RUNSUB
     noahmp%config%nmlist%OptGlacierTreatment         = NoahmpIO%IOPT_GLA
     noahmp%config%nmlist%OptSnowCompaction           = NoahmpIO%IOPT_COMPACT
+    noahmp%config%nmlist%OptWetlandModel             = NoahmpIO%IOPT_WETLAND
 
     ! config domain variable
     noahmp%config%domain%SurfaceType                 = 1

drivers/lis/LisNoahmpParamType.F90

@@ -210,6 +210,9 @@ module LisNoahmpParamType
     real(kind=kind_noahmp)      :: RTCT(NSTAGE)        ! fraction of carbohydrate translocation from root to grain
     real(kind=kind_noahmp)      :: BIO2LAI             ! leaf area per living leaf biomass [m2/kg]
 
+    ! wetland parameters
+    real(kind=kind_noahmp)      :: WCAP                ! maximum wetland water holding capacity [m] (tunable) for opt_wetland=1
+
     ! soil parameters
     real(kind=kind_noahmp)      :: BEXP(NSOIL)         ! soil B parameter
     real(kind=kind_noahmp)      :: SMCDRY(NSOIL)       ! dry soil moisture threshold

drivers/lis/NoahmpDriverMainMod.F90

@@ -68,6 +68,7 @@ subroutine NoahmpDriverMain(NoahmpIO, LISparam)
        NoahmpIO%ACC_ECANXY   = 0.0
        NoahmpIO%ACC_ETRANXY  = 0.0
        NoahmpIO%ACC_EDIRXY   = 0.0
+       NoahmpIO%ACC_GLAFLWXY = 0.0
     endif
 
     if ( NoahmpIO%SOIL_UPDATE_STEPS > 1 ) then
@@ -81,6 +82,7 @@ subroutine NoahmpDriverMain(NoahmpIO, LISparam)
           NoahmpIO%ACC_ECANXY   = 0.0
           NoahmpIO%ACC_ETRANXY  = 0.0
           NoahmpIO%ACC_EDIRXY   = 0.0
+          NoahmpIO%ACC_GLAFLWXY = 0.0
        end if
     endif