Skip to content

Commit

Permalink
Add Xinanjiang Surface Runoff Scheme (#559)
Browse files Browse the repository at this point in the history 
* Xinanjiang runoff scheme added
* spatial soil XAJ params added
* change NWM tests to use new runoff scheme

Co-authored-by: Katelyn FitzGerald <katelynw@ucar.edu>
  • Loading branch information
@prasanthvkrishna @kafitzgerald
prasanthvkrishna and kafitzgerald authored Jun 9, 2021
1 parent df784c8 commit 8a61c8f
Show file tree
Hide file tree
Showing 6 changed files with 234 additions and 55 deletions.
12 changes: 12 additions & 0 deletions trunk/NDHMS/Land_models/NoahMP/IO_code/module_NoahMP_hrldas_driver.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,9 @@ module module_NoahMP_hrldas_driver
REAL, ALLOCATABLE, DIMENSION(:,:) :: hvt_2D ! Canopy Height
REAL, ALLOCATABLE, DIMENSION(:,:) :: mfsno_2D ! Snow cover m parameter
REAL, ALLOCATABLE, DIMENSION(:,:) :: rsurfexp_2D! exponent in the shape parameter for soil resistance option 1
REAL, ALLOCATABLE, DIMENSION(:,:) :: axaj_2D ! Tension water distribution inflection parameter [-]
REAL, ALLOCATABLE, DIMENSION(:,:) :: bxaj_2D ! Tension water distribution shape parameter [-]
REAL, ALLOCATABLE, DIMENSION(:,:) :: xxaj_2D ! Free water distribution shape parameter [-]
#endif

! INOUT (with generic LSM equivalent) (as defined in WRF)
Expand Down Expand Up @@ -741,6 +744,9 @@ subroutine land_driver_ini(NTIME_out, state, wrfits,wrfite,wrfjts,wrfjte)
ALLOCATE ( hvt_2D (XSTART:XEND,YSTART:YEND) ) ! Canopy Height
ALLOCATE ( mfsno_2D (XSTART:XEND,YSTART:YEND) ) ! Snow cover m parameter
ALLOCATE ( rsurfexp_2D(XSTART:XEND,YSTART:YEND) ) ! exponent in the shape parameter for soil resistance option 1
ALLOCATE ( axaj_2D (XSTART:XEND,YSTART:YEND) ) ! Tension water distribution inflection parameter [-]
ALLOCATE ( bxaj_2D (XSTART:XEND,YSTART:YEND) ) ! Tension water distribution shape parameter [-]
ALLOCATE ( xxaj_2D (XSTART:XEND,YSTART:YEND) ) ! Free water distribution shape parameter [-]
#endif
! INOUT (with generic LSM equivalent) (as defined in WRF)
Expand Down Expand Up @@ -1087,6 +1093,11 @@ subroutine land_driver_ini(NTIME_out, state, wrfits,wrfite,wrfjts,wrfjte)
NSOIL,BEXP_3D,SMCDRY_3D,SMCWLT_3D,SMCREF_3D,SMCMAX_3D, &
DKSAT_3D,DWSAT_3D,PSISAT_3D,QUARTZ_3D,REFDK_2D,REFKDT_2D,&
SLOPE_2D,CWPVT_2D,VCMX25_2D,MP_2D,HVT_2D,MFSNO_2D,RSURFEXP_2D)
if (noah_lsm%runoff_option == 7) then
CALL READ_XAJ_RUNOFF(noah_lsm%SPATIAL_FILENAME,XSTART, XEND, YSTART, YEND, &
AXAJ_2D,BXAJ_2D,XXAJ_2D)
end if
#endif
!------------------------------------------------------------------------
Expand Down Expand Up @@ -1658,6 +1669,7 @@ subroutine land_driver_exe(itime, state)
DKSAT_3D,DWSAT_3D,PSISAT_3D,QUARTZ_3D, &
REFDK_2D,REFKDT_2D,SLOPE_2D, &
CWPVT_2D,VCMX25_2D,MP_2D,HVT_2D,MFSNO_2D,RSURFEXP_2D, &
AXAJ_2D,BXAJ_2D,XXAJ_2D, &
#endif
#ifdef WRF_HYDRO
sfcheadrt,INFXSRT,soldrain, & !O
Expand Down
63 changes: 62 additions & 1 deletion trunk/NDHMS/Land_models/NoahMP/IO_code/module_hrldas_netcdf_io.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -771,6 +771,67 @@ subroutine read_mmf_runoff(wrfinput_flnm, &

end subroutine read_mmf_runoff

!-------------------------------------------------------------------------------------------

subroutine read_xaj_runoff(spatial_filename,xstart, xend,ystart, yend, &
axaj_2d,bxaj_2d,xxaj_2d)

implicit none
character(len=*), intent(in) :: spatial_filename
integer, intent(in) :: xstart, xend, ystart, yend
real, dimension(xstart:xend,ystart:yend), intent(out) :: axaj_2d
real, dimension(xstart:xend,ystart:yend), intent(out) :: bxaj_2d
real, dimension(xstart:xend,ystart:yend), intent(out) :: xxaj_2d

character(len=24) :: name
character(len=256) :: units
integer :: ierr,iret, varid
integer :: ncid
real, dimension(xstart:xend,ystart:yend) :: xdum
!------------------------------------------------------------------------------------------------------

ierr = nf90_open(spatial_filename, NF90_NOWRITE, ncid)
if (ierr /= 0) then
write(*,'("FATAL ERROR: In read_xaj_runoff(): Problem opening soil file: ''", A, "''")') trim(spatial_filename)
stop
endif

name = "AXAJ"
iret = nf90_inq_varid(ncid, trim(name), varid)
if (iret /= 0) then
print*, 'ncid = ', ncid
write(*,*) "FATAL ERROR: In read_xaj_runoff(): Problem finding variable '"//trim(name)//"' in NetCDF file: " // trim(spatial_filename)
stop
endif

iret = nf90_get_var(ncid, varid, axaj_2d, start=(/xstart,ystart/), count=(/xend-xstart+1,yend-ystart+1/))
!
name = "BXAJ"
iret = nf90_inq_varid(ncid, trim(name), varid)
if (iret /= 0) then
print*, 'ncid = ', ncid
write(*,*) "FATAL ERROR: In read_xaj_runoff(): Problem finding variable '"//trim(name)//"' in NetCDF file: " // trim(spatial_filename)
stop
endif

iret = nf90_get_var(ncid, varid, bxaj_2d, start=(/xstart,ystart/), count=(/xend-xstart+1,yend-ystart+1/))
!
name = "XXAJ"
iret = nf90_inq_varid(ncid, trim(name), varid)
if (iret /= 0) then
print*, 'ncid = ', ncid
write(*,*) "FATAL ERROR: In read_xaj_runoff(): Problem finding variable '"//trim(name)//"' in NetCDF file: " // trim(spatial_filename)
stop
endif

iret = nf90_get_var(ncid, varid, xxaj_2d, start=(/xstart,ystart/), count=(/xend-xstart+1,yend-ystart+1/))
!
! Close the NetCDF file
ierr = nf90_close(ncid)
if (ierr /= 0) stop "FATAL ERROR: In module_hrldas_netcdf_io.F read_xaj_runoff() - NF90_CLOSE"

end subroutine read_xaj_runoff

!---------------------------------------------------------------------------------------------------------

subroutine read_3d_soil(spatial_filename,xstart, xend,ystart, yend, &
Expand Down Expand Up @@ -1027,7 +1088,7 @@ subroutine read_3d_soil(spatial_filename,xstart, xend,ystart, yend, &
endif

iret = nf90_get_var(ncid, varid, rsurfexp_2d, start=(/xstart,ystart/), count=(/xend-xstart+1,yend-ystart+1/))

!
! Close the NetCDF file
ierr = nf90_close(ncid)
if (ierr /= 0) stop "FATAL ERROR: In module_hrldas_netcdf_io.F read_3d_soil() - NF90_CLOSE"
Expand Down
13 changes: 11 additions & 2 deletions trunk/NDHMS/Land_models/NoahMP/phys/module_sf_noahmpdrv.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN, XLATIN, & ! IN
DKSAT_3D,DWSAT_3D,PSISAT_3D,QUARTZ_3D, &
REFDK_2D,REFKDT_2D,SLOPE_2D, &
CWPVT_2D,VCMX25_2D,MP_2D,HVT_2D,MFSNO_2D,RSURFEXP_2D, &
AXAJ_2D,BXAJ_2D,XXAJ_2D, &
#endif
#ifdef WRF_HYDRO
sfcheadrt,INFXSRT,soldrain, &
Expand Down Expand Up @@ -161,6 +162,9 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN, XLATIN, & ! IN
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: HVT_2D ! Canopy Height
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: MFSNO_2D ! Snow cover m parameter
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: RSURFEXP_2D ! exponent in the shape parameter for soil resistance option 1
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: AXAJ_2D ! Xinanjiang: Tension water distribution inflection parameter [-]
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: BXAJ_2D ! Xinanjiang: Tension water distribution shape parameter [-]
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XXAJ_2D ! Xinanjiang: Free water distribution shape parameter [-]
#endif

! INOUT (with generic LSM equivalent)
Expand Down Expand Up @@ -660,6 +664,9 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN, XLATIN, & ! IN
parameters%hvt = HVT_2D(I,J) ! Canopy Height
parameters%mfsno = MFSNO_2D(I,J) ! Snow cover m parameter
parameters%rsurf_exp = RSURFEXP_2D(I,J) ! exponent in the shape parameter for soil resistance option 1
parameters%axaj = AXAJ_2D(I,J) ! Xinanjiang: Tension water distribution inflection parameter [-]
parameters%bxaj = BXAJ_2D(I,J) ! Xinanjiang: Tension water distribution shape parameter [-]
parameters%xxaj = XXAJ_2D(I,J) ! Xinanjiang: Free water distribution shape parameter [-]
#endif
CALL TRANSFER_MP_PARAMETERS(VEGTYP,SOILTYP,SLOPETYP,SOILCOLOR,CROPTYPE,parameters)

Expand Down Expand Up @@ -1193,9 +1200,11 @@ SUBROUTINE TRANSFER_MP_PARAMETERS(VEGTYPE,SOILTYPE,SLOPETYPE,SOILCOLOR,CROPTYPE,
parameters%MFSNO = MFSNO_TABLE(VEGTYPE) !snowmelt m parameter ()
parameters%HVT = HVT_TABLE(VEGTYPE) !top of canopy (m)
parameters%MP = MP_TABLE(VEGTYPE) !slope of conductance-to-photosynthesis relationship
parameters%RSURF_EXP = RSURF_EXP_TABLE
parameters%RSURF_EXP = RSURF_EXP_TABLE
parameters%AXAJ = AXAJ_TABLE(SOILTYPE)
parameters%BXAJ = BXAJ_TABLE(SOILTYPE)
parameters%XXAJ = XXAJ_TABLE(SOILTYPE)
#endif

! ----------------------------------------------------------------------
! Transfer GENPARM parameters
! ----------------------------------------------------------------------
Expand Down
111 changes: 101 additions & 10 deletions trunk/NDHMS/Land_models/NoahMP/phys/module_sf_noahmplsm.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -361,6 +361,9 @@ MODULE MODULE_SF_NOAHMPLSM
REAL :: DWSAT(NSOIL) !saturated soil hydraulic diffusivity
REAL :: QUARTZ(NSOIL) !soil quartz content
REAL :: F1 !soil thermal diffusivity/conductivity coef (not used MB: 20140718)
REAL :: AXAJ !Xinanjiang: Tension water distribution inflection parameter [-]
REAL :: BXAJ !Xinanjiang: Tension water distribution shape parameter [-]
REAL :: XXAJ !Xinanjiang: Free water distribution shape parameter [-]
!------------------------------------------------------------------------------------------!
! From the GENPARM.TBL file
!------------------------------------------------------------------------------------------!
Expand Down Expand Up @@ -5841,7 +5844,7 @@ SUBROUTINE WATER (parameters,VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , &
RUNSUB = QDIS !mm/s
END IF
IF(OPT_RUN == 3 .or. OPT_RUN == 4) THEN
IF(OPT_RUN == 3 .or. OPT_RUN == 4 .or. OPT_RUN == 7) THEN
RUNSUB = RUNSUB + QDRAIN !mm/s
END IF
Expand Down Expand Up @@ -6995,7 +6998,9 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in
PDDUM = QINSUR - RUNSRF ! m/s
END IF
END IF
IF (OPT_RUN == 7) THEN
CALL COMPUTE_XAJ_SURFRUNOFF(parameters,DT,FCR,NSOIL,SMC,ZSOIL,QINSUR,RUNSRF,PDDUM)
END IF
! determine iteration times and finer time step
NITER = 1
Expand All @@ -7020,6 +7025,11 @@ SUBROUTINE SOILWATER (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in
PDDUM ,RUNSRF ) !out
END IF
IF (QINSUR > 0. .AND. OPT_RUN == 7) THEN
CALL COMPUTE_XAJ_SURFRUNOFF(parameters,DTFINE,FCR,NSOIL,SMC,ZSOIL,QINSUR,& ! in
RUNSRF,PDDUM) ! out
END IF
CALL SRT (parameters,NSOIL ,ZSOIL ,DTFINE ,PDDUM ,ETRANI , & !in
QSEVA ,SH2O ,SMC ,ZWT ,FCR , & !in
SICEMAX,FCRMAX ,ILOC ,JLOC ,SMCWTD , & !in
Expand Down Expand Up @@ -7343,7 +7353,7 @@ SUBROUTINE SRT (parameters,NSOIL ,ZSOIL ,DT ,PDDUM ,ETRANI , & !in
IF(OPT_RUN == 1 .or. OPT_RUN == 2) THEN
QDRAIN = 0.
END IF
IF(OPT_RUN == 3) THEN
IF(OPT_RUN == 3 .OR. OPT_RUN == 7) THEN
QDRAIN = parameters%SLOPE*WCND(K)
END IF
IF(OPT_RUN == 4) THEN
Expand Down Expand Up @@ -7511,6 +7521,82 @@ SUBROUTINE SSTEP (parameters,NSOIL ,NSNOW ,DT ,ZSOIL ,DZSNSO , & !in
END SUBROUTINE SSTEP
!== begin xinanjiag=================================================================================
SUBROUTINE COMPUTE_XAJ_SURFRUNOFF(parameters,DT,FCR,NSOIL,SMC,ZSOIL,QINSUR,RUNSRF,PDDUM)
! ----------------------------------------------------------------------
! Calculate the saturated area and runoff based on Xinanjiag runoff scheme.
! Reference: Knoben, W. J., Freer, J. E., Fowler, K. J., Peel, M. C., & Woods, R. A. (2019).
! Modular Assessment of Rainfall-Runoff Models Toolbox (MARRMoT) v1. 2:
! an open-source, extendable framework providing implementations of 46 conceptual
! hydrologic models as continuous state-space formulations.
! ----------------------------------------------------------------------
! Author: Prasanth Valayamkunnath <prasanth@ucar.edu>
! Date: August 03, 2020
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! Inputs
TYPE (noahmp_parameters), INTENT(IN) :: parameters
INTEGER, INTENT(IN) :: NSOIL
REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC
REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
REAL, DIMENSION(1:NSOIL), INTENT(IN) :: FCR !fraction of imperviousness (-) = IMP
REAL , INTENT(IN) :: QINSUR
REAL , INTENT(IN) :: DT
! Output
REAL , INTENT(INOUT):: RUNSRF
REAL , INTENT(INOUT):: PDDUM
! local
REAL :: WM,WM_MAX,SM,SM_MAX,IRUNOFF,PRUNOFF
INTEGER :: IZ
!------------------------------------------------------------------------
!initialize
WM = 0.0
WM_MAX = 0.0
SM = 0.0
SM_MAX = 0.0
IRUNOFF = 0.0
PRUNOFF = 0.0
RUNSRF = 0.0
DO IZ=1,NSOIL-2
IF ((SMC(IZ)-parameters%SMCREF(IZ)) .GT. 0.) THEN ! soil moisture greater than field capacity
SM = SM + (SMC(IZ) - parameters%SMCREF(IZ) )*-1*ZSOIL(IZ) !m
WM = WM + (parameters%SMCREF(IZ)*-1*ZSOIL(IZ)) !m
ELSE
WM = WM + (SMC(IZ)*-1*ZSOIL(IZ))
END IF
WM_MAX = WM_MAX + (parameters%SMCREF(IZ)*-1*ZSOIL(IZ))
SM_MAX = SM_MAX + (parameters%SMCMAX(IZ) - parameters%SMCREF(IZ))*-1*ZSOIL(IZ)
END DO
WM = MIN(WM,WM_MAX) ! tension water (m)
SM = MIN(SM,SM_MAX) ! free water (m)
! impervious surface runoff R_IMP
IRUNOFF = FCR(1)*QINSUR*DT
! solve pervious surface runoff (m) based on Eq. (310)
IF ((WM/WM_MAX) .LE. (0.5-parameters%AXAJ))THEN
PRUNOFF = (1-FCR(1))*QINSUR*DT*((0.5-parameters%AXAJ)**(1-parameters%BXAJ))*((WM/WM_MAX)**parameters%BXAJ)
ELSE
PRUNOFF = (1-FCR(1))*QINSUR*DT*(1-(((0.5+parameters%AXAJ)**(1-parameters%BXAJ))*((1-(WM/WM_MAX))**parameters%BXAJ)))
END IF
! estimate surface runoff based on Eq. (313)
IF(QINSUR .EQ. 0.0) THEN
RUNSRF = 0.0
ELSE
RUNSRF = PRUNOFF*(1-((1-(SM/SM_MAX))**parameters%XXAJ))+IRUNOFF
END IF
RUNSRF = RUNSRF/DT !m/s
RUNSRF = MAX(0.0, RUNSRF)
RUNSRF = MIN(QINSUR, RUNSRF)
PDDUM = QINSUR - RUNSRF
END SUBROUTINE COMPUTE_XAJ_SURFRUNOFF
!== end xinanjiag ==================================================================================
!== begin wdfcnd1 ==================================================================================
SUBROUTINE WDFCND1 (parameters,WDF,WCND,SMC,FCR,ISOIL)
Expand Down Expand Up @@ -9117,17 +9203,19 @@ MODULE NOAHMP_TABLES
INTEGER :: SLCATS
REAL :: BEXP_TABLE(MAX_SOILTYP) !maximum intercepted h2o per unit lai+sai (mm)
REAL :: BEXP_TABLE(MAX_SOILTYP) !maximum intercepted h2o per unit lai+sai (mm)
REAL :: SMCDRY_TABLE(MAX_SOILTYP) !characteristic leaf dimension (m)
REAL :: F1_TABLE(MAX_SOILTYP) !momentum roughness length (m)
REAL :: F1_TABLE(MAX_SOILTYP) !momentum roughness length (m)
REAL :: SMCMAX_TABLE(MAX_SOILTYP) !top of canopy (m)
REAL :: SMCREF_TABLE(MAX_SOILTYP) !bottom of canopy (m)
REAL :: PSISAT_TABLE(MAX_SOILTYP) !tree density (no. of trunks per m2)
REAL :: DKSAT_TABLE(MAX_SOILTYP) !tree crown radius (m)
REAL :: DWSAT_TABLE(MAX_SOILTYP) !monthly stem area index, one-sided
REAL :: SMCWLT_TABLE(MAX_SOILTYP) !monthly leaf area index, one-sided
REAL :: QUARTZ_TABLE(MAX_SOILTYP) !single-side leaf area per Kg [m2/kg]
REAL :: QUARTZ_TABLE(MAX_SOILTYP) !single-side leaf area per Kg [m2/kg]
REAL :: AXAJ_TABLE(MAX_SOILTYP) !Xinanjiang: Tension water distribution inflection parameter [-]
REAL :: BXAJ_TABLE(MAX_SOILTYP) !Xinanjiang: Tension water distribution shape parameter [-]
REAL :: XXAJ_TABLE(MAX_SOILTYP) !Xinanjiang: Free water distribution shape parameter [-]
! GENPARM.TBL parameters
REAL :: SLOPE_TABLE(9) !slope factor for soil drainage
Expand Down Expand Up @@ -9532,8 +9620,10 @@ subroutine read_mp_soil_parameters()
FRZK_TABLE = -1.E36
ZBOT_TABLE = -1.E36
CZIL_TABLE = -1.E36
!
AXAJ_TABLE = -1.E36
BXAJ_TABLE = -1.E36
XXAJ_TABLE = -1.E36
!
!-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
!
inquire( file='SOILPARM.TBL', exist=file_named )
Expand All @@ -9560,7 +9650,8 @@ subroutine read_mp_soil_parameters()
DO LC=1,SLCATS
READ (21,*) ITMP,BEXP_TABLE(LC),SMCDRY_TABLE(LC),F1_TABLE(LC),SMCMAX_TABLE(LC), &
SMCREF_TABLE(LC),PSISAT_TABLE(LC),DKSAT_TABLE(LC), DWSAT_TABLE(LC), &
SMCWLT_TABLE(LC), QUARTZ_TABLE(LC)
SMCWLT_TABLE(LC), QUARTZ_TABLE(LC), AXAJ_TABLE (LC), BXAJ_TABLE (LC), &
XXAJ_TABLE(LC)
ENDDO
CLOSE (21)
Expand Down
Loading

0 comments on commit 8a61c8f

Please sign in to comment.