Add wetland scheme of Zhang et al 2022

hrldas/IO_code/module_NoahMP_hrldas_driver.F

@@ -293,6 +293,16 @@ subroutine land_driver_ini(NTIME_out,wrfits,wrfite,wrfjts,wrfjte)
        call READ_TILE_DRAIN_MAP(NoahmpIO%TDINPUT_FLNM, NoahmpIO%XSTART, NoahmpIO%XEND, NoahmpIO%YSTART, NoahmpIO%YEND, NoahmpIO%TD_FRACTION)
     endif
 
+!------------------------------------------------------------------------
+! For surface wetland scheme, IOPT_WETLAND = 1 or 2, read in necessary extra fields 
+!------------------------------------------------------------------------
+    if (NoahmpIO%IOPT_WETLAND == 2) then ! use wetland_option=2 for reading wetland extra fields
+       NoahmpIO%FSATMX = 0.38
+       NoahmpIO%WCAP   = 0.1
+       call READ_WETLAND_INPUT(NoahmpIO%HRLDAS_SETUP_FILE, NoahmpIO%XSTART, NoahmpIO%XEND, NoahmpIO%YSTART, NoahmpIO%YEND,&
+                               NoahmpIO%FSATMX,NoahmpIO%WCAP)
+    endif
+
 !------------------------------------------------------------------------
 ! For SF_URBAN_PHYSICS > 0 read 2D map of urban parameters
 !------------------------------------------------------------------------
@@ -383,6 +393,7 @@ subroutine land_driver_ini(NTIME_out,wrfits,wrfite,wrfjts,wrfjte)
        call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "ACC_ECAN"  , NoahmpIO%ACC_ECANXY  )
        call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "ACC_ETRAN" , NoahmpIO%ACC_ETRANXY )
        call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "ACC_EDIR"  , NoahmpIO%ACC_EDIRXY  )
+       call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "ACC_GLAFLW", NoahmpIO%ACC_GLAFLWXY)
 
        ! below for irrigation scheme
        if ( NoahmpIO%IOPT_IRR > 0 ) then
@@ -418,6 +429,12 @@ subroutine land_driver_ini(NTIME_out,wrfits,wrfite,wrfjts,wrfjte)
           call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "PEXPXY"     ,NoahmpIO%PEXPXY     )
        endif
 
+       ! for wetland scheme
+       if ( NoahmpIO%IOPT_WETLAND > 0 ) then
+          call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "FSATXY"     , NoahmpIO%FSATXY    )
+          call get_from_restart(xstart, xend, xstart, ixfull, jxfull, "WSURFXY"    , NoahmpIO%WSURFXY   )
+       endif
+
        ! below for urban model
        if ( NoahmpIO%SF_URBAN_PHYSICS > 0 ) then
           call get_from_restart(xstart, xend, xstart, ixfull, jxfull,      "SH_URB2D"  ,     NoahmpIO%SH_URB2D  )
@@ -1228,6 +1245,7 @@ subroutine land_driver_exe(itime)
               call add_to_output(NoahmpIO%ACC_ECANXY  ,"ACC_ECAN"  , "accumulated net canopy evaporation within soil timestep", "mm")
               call add_to_output(NoahmpIO%ACC_ETRANXY ,"ACC_ETRAN" , "accumulated transpiration within soil timestep"         , "mm")
               call add_to_output(NoahmpIO%ACC_EDIRXY  ,"ACC_EDIR"  , "accumulated net ground evaporation within soil timestep", "mm")
+              call add_to_output(NoahmpIO%ACC_GLAFLWXY,"ACC_GLAFLW", "accumuated glacier excessive flow per soil timestep"    , "mm")
               ! additional energy terms
               call add_to_output(NoahmpIO%PAHXY      , "PAH"     , "Precipitation advected heat flux"                         , "W/m2"    )
               call add_to_output(NoahmpIO%PAHGXY     , "PAHG"    , "Precipitation advected heat flux to below-canopy ground"  , "W/m2"    )
@@ -1270,7 +1288,12 @@ subroutine land_driver_exe(itime)
               call add_to_output(NoahmpIO%QSLATXY    , "QSLAT"    , "accumulated lateral flow"                                     , "mm"     )
               call add_to_output(NoahmpIO%QLATXY     , "QLAT"     , "instantaneous lateral flow"                                   , "m"      )
             endif
-
+            ! Wetland model
+            if ( NoahmpIO%IOPT_WETLAND > 0 ) then
+              call add_to_output(NoahmpIO%FSATXY     , "FSAT"    , "saturated fraction of the grid"                           , "-" )
+              call add_to_output(NoahmpIO%WSURFXY    , "WSURF"   , "Wetland Water Storage"                                    , "mm")
+            endif
+
            ! For now, no urban output variables included
 
            enddo DEFINE_MODE_LOOP
@@ -1407,6 +1430,8 @@ subroutine lsm_restart()
   call add_to_restart(NoahmpIO%ACC_ECANXY  ,"ACC_ECAN"  )
   call add_to_restart(NoahmpIO%ACC_ETRANXY ,"ACC_ETRAN" )
   call add_to_restart(NoahmpIO%ACC_EDIRXY  ,"ACC_EDIR"  )
+  call add_to_restart(NoahmpIO%ACC_GLAFLWXY,"ACC_GLAFLW")
+
  ! irrigation scheme
   if ( NoahmpIO%IOPT_IRR > 0 ) then
      call add_to_restart(NoahmpIO%IRNUMSI   , "IRNUMSI")
@@ -1421,6 +1446,7 @@ subroutine lsm_restart()
      call add_to_restart(NoahmpIO%IRELOSS   , "IRELOSS")
      call add_to_restart(NoahmpIO%IRRSPLH   , "IRRSPLH")
   endif
+
   ! below for MMF groundwater scheme
   if ( NoahmpIO%IOPT_RUNSUB == 5 ) then
      call add_to_restart(NoahmpIO%SMOISEQ     , "SMOISEQ"  , layers="SOIL"  )
@@ -1433,12 +1459,19 @@ subroutine lsm_restart()
      call add_to_restart(NoahmpIO%RECHXY      , "RECHXY"     )
      call add_to_restart(NoahmpIO%QRFXY       , "QRFXY"      )
      call add_to_restart(NoahmpIO%QSPRINGXY   , "QSPRINGXY"  )
-     call add_to_restart(NoahmpIO%FDEPTHXY    , "FDEPTHXY"  )
-     call add_to_restart(NoahmpIO%RIVERCONDXY , "RIVERCONDXY"  )
-     call add_to_restart(NoahmpIO%RIVERBEDXY  , "RIVERBEDXY"  )
-     call add_to_restart(NoahmpIO%EQZWT       , "EQZWT"  )
-     call add_to_restart(NoahmpIO%PEXPXY      , "PEXPXY"  )
+     call add_to_restart(NoahmpIO%FDEPTHXY    , "FDEPTHXY"   )
+     call add_to_restart(NoahmpIO%RIVERCONDXY , "RIVERCONDXY")
+     call add_to_restart(NoahmpIO%RIVERBEDXY  , "RIVERBEDXY" )
+     call add_to_restart(NoahmpIO%EQZWT       , "EQZWT"      )
+     call add_to_restart(NoahmpIO%PEXPXY      , "PEXPXY"     )
   endif
+
+  ! for wetland scheme
+  if ( NoahmpIO%IOPT_WETLAND > 0 ) then
+     call add_to_restart(NoahmpIO%FSATXY      , "FSATXY"     )
+     call add_to_restart(NoahmpIO%WSURFXY     , "WSURFXY"    )
+  endif
+
   ! below for urban model
   if ( NoahmpIO%SF_URBAN_PHYSICS > 0 ) then
       call add_to_restart(     NoahmpIO%SH_URB2D  ,      "SH_URB2D" )

hrldas/IO_code/module_hrldas_netcdf_io.F

@@ -666,7 +666,9 @@ subroutine read_crop_input(wrfinput_flnm,                            &
        xstart, xend,                                                   &
        ystart, yend,                                                   &
        croptype,planting,harvest,season_gdd)
+
     implicit none
+
     character(len=*),          intent(in)  :: wrfinput_flnm
     integer,                   intent(in)  :: xstart, xend, ystart, yend
     real,    dimension(xstart:xend,5,ystart:yend), intent(out) :: croptype
@@ -751,6 +753,75 @@ subroutine read_crop_input(wrfinput_flnm,                            &
 
   end subroutine read_crop_input
 
+!---------------------------------------------------------------------------------------------------------
+
+  subroutine read_wetland_input(wrfinput_flnm,                         &
+       xstart, xend,                                                   &
+       ystart, yend,                                                   &
+       fsatmx, wcap)
+
+    implicit none
+
+    character(len=*),          intent(in)  :: wrfinput_flnm
+    integer,                   intent(in)  :: xstart, xend, ystart, yend
+    real,    dimension(xstart:xend,  ystart:yend), intent(out) :: fsatmx
+    real,    dimension(xstart:xend,  ystart:yend), intent(out) :: wcap
+
+    character(len=256) :: units
+    character(len=24)  :: name
+    integer :: ierr,iret
+    integer :: ncid, varid
+    integer :: rank
+
+#ifdef _PARALLEL_  
+    call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
+    if (ierr /= MPI_SUCCESS) stop "MPI_COMM_RANK"
+#else
+    rank = 0
+#endif
+
+
+    ! Open the NetCDF file.
+    if (rank == 0) write(*,'("wrfinput_flnm: ''", A, "''")') trim(wrfinput_flnm)
+#ifdef _PARALLEL_
+    ierr = nf90_open_par(wrfinput_flnm, NF90_NOWRITE, MPI_COMM_WORLD, MPI_INFO_NULL, ncid)
+#else
+    ierr = nf90_open(wrfinput_flnm, NF90_NOWRITE, ncid)
+#endif
+    if (ierr /= 0) then
+       write(*,'("read_wetland_input:  Problem opening wrfinput file: ''", A, "''")') trim(wrfinput_flnm)
+#ifdef _PARALLEL_
+       call mpi_finalize(ierr)
+       if (ierr /= 0) write(*, '("Problem with MPI_finalize.")')
+#endif
+       stop
+    endif
+
+! Get maximum saturated fraction (FSATMX)
+    name = "FSATMX"
+    iret = nf90_inq_varid(ncid,  name,  varid)
+    if (iret == 0) then
+      ierr = nf90_get_var(ncid, varid, fsatmx, start=(/xstart,ystart/), count=(/xend-xstart+1,yend-ystart+1/))
+    else
+      write(*,*) "MODULE_HRLDAS_NETCDF_IO:  Problem finding variable '"//trim(name)//"' in NetCDF file. Using default values."
+    endif
+
+! Get maximum water storeage (WSTORE)
+    name = "WCAP"
+    iret = nf90_inq_varid(ncid,  name,  varid)
+    if (iret == 0) then
+      ierr = nf90_get_var(ncid, varid, wcap, start=(/xstart,ystart/), count=(/xend-xstart+1,yend-ystart+1/))
+    else
+      write(*,*) "MODULE_HRLDAS_NETCDF_IO:  Problem finding variable '"//trim(name)//"' in NetCDF file. Using default values."
+    endif
+
+
+! Close the NetCDF file
+    ierr = nf90_close(ncid)
+    if (ierr /= 0) stop "MODULE_NOAHLSM_HRLDAS_INPUT:  read_wetland_input:  NF90_CLOSE"
+
+  end subroutine read_wetland_input
+
 !---------------------------------------------------------------------------------------------------------
 
   subroutine read_tile_drain_map(tdinput_flnm,                            &

hrldas/run/Makefile

@@ -182,7 +182,9 @@ OBJS = \
         ../../noahmp/src/SurfaceEmissivityGlacierMod.o \
         ../../noahmp/src/SurfaceEnergyFluxGlacierMod.o \
         ../../noahmp/src/SurfaceRadiationGlacierMod.o \
-        ../../noahmp/src/WaterMainGlacierMod.o
+        ../../noahmp/src/RunoffSurfaceWetlandMod.o \
+        ../../noahmp/src/WetlandWaterZhang22Mod.o \
+	../../noahmp/src/WaterMainGlacierMod.o
 
 CMD = hrldas.exe
 all:	$(CMD)

hrldas/run/README.namelist

@@ -261,6 +261,11 @@
                                                              !   1 -> Simple drainage
                                                              !   2 -> Hooghoudt's equation based tile drainage
 
+ WETLAND_OPTION                    = 0                       ! options for wetland model [default = 0]
+                                                             ! wetland water budget model in (Zhang et al. 2022 WRR)
+                                                             !   1 -> Uniform fixed wetland parameters from table 
+                                                             !   2 -> 2-D, read parameters from the input file
+
  SOIL_TIMESTEP                     = 0.0,                    ! Noah-MP soil process timestep (seconds) for solving soil water and temperature
                                                              !   0 -> default, the same as main NoahMP model timestep
                                                              !   N * dt_noahmp -> longer than main NoahMP model timestep (often used for WRF coupled run)

hrldas/run/namelist.hrldas_example

@@ -1,16 +1,16 @@
 &NOAHLSM_OFFLINE
 
- HRLDAS_SETUP_FILE = "/glade/scratch/cenlinhe/NoahMP_NLDAS/HRLDAS_setup_2000010100_d1"
- INDIR = "/glade/scratch/cenlinhe/NoahMP_NLDAS/forcing/"
- OUTDIR = "/glade/scratch/cenlinhe/test/hrldas_output/"
+ HRLDAS_SETUP_FILE = "/init_condition_dir/HRLDAS_setup_2000010100_d1"
+ INDIR = "/atmos_forcing_dir/forcing/"
+ OUTDIR = "/model_output_dir/hrldas_output/"
 
  START_YEAR  = 2000
  START_MONTH = 01
  START_DAY   = 01
  START_HOUR  = 00
  START_MIN   = 00
 
-! RESTART_FILENAME_REQUESTED = "/glade/scratch/cenlinhe/test/hrldas_output/RESTART.2000010100_DOMAIN1"
+! RESTART_FILENAME_REQUESTED = "/model_restart_dir/RESTART.2000010100_DOMAIN1"
 
  KDAY = 365
  SPINUP_LOOPS = 0
@@ -34,20 +34,21 @@
  FROZEN_SOIL_OPTION                = 1
  SUPERCOOLED_WATER_OPTION          = 1
  RADIATIVE_TRANSFER_OPTION         = 3
- SNOW_ALBEDO_OPTION                = 2
+ SNOW_ALBEDO_OPTION                = 1
  SNOW_COMPACTION_OPTION            = 2
  PCP_PARTITION_OPTION              = 1
  SNOW_THERMAL_CONDUCTIVITY         = 1
  TBOT_OPTION                       = 2
- TEMP_TIME_SCHEME_OPTION           = 1
+ TEMP_TIME_SCHEME_OPTION           = 3
  GLACIER_OPTION                    = 1
- SURFACE_RESISTANCE_OPTION         = 1
+ SURFACE_RESISTANCE_OPTION         = 4
  SOIL_DATA_OPTION                  = 1
  PEDOTRANSFER_OPTION               = 1
  CROP_OPTION                       = 0
  IRRIGATION_OPTION                 = 0
  IRRIGATION_METHOD                 = 0
  TILE_DRAINAGE_OPTION              = 0
+ WETLAND_OPTION                    = 0
 
  FORCING_TIMESTEP = 3600
  NOAH_TIMESTEP    = 1800