diff --git a/compiler.sh b/compiler.sh index 7ffc55a96..9b86d96c2 100755 --- a/compiler.sh +++ b/compiler.sh @@ -61,13 +61,16 @@ if [[ "$build_mc_kernel" == true ]]; then fi if [[ "$build_diffusive_tulane_kernel" == true ]]; then - #building reach and resevoir kernel files .o + #building diffusive kernel files .o cd $REPOROOT/src/kernel/diffusive/ make clean make diffusive.o make pydiffusive.o make chxsec_lookuptable.o make pychxsec_lookuptable.o + make diffusive_lightweight.o + make pydiffusive_lightweight.o + make precis.mod make install || exit fi diff --git a/src/kernel/diffusive/chxsec_lookuptable.f90 b/src/kernel/diffusive/chxsec_lookuptable.f90 index 85b2f3da4..bdbb22b22 100644 --- a/src/kernel/diffusive/chxsec_lookuptable.f90 +++ b/src/kernel/diffusive/chxsec_lookuptable.f90 @@ -1,6 +1,6 @@ !------------------------------------------------------------------------------- module lookuptable - + use precis IMPLICIT NONE !----------------------------------------------------------------------------- @@ -16,20 +16,20 @@ module lookuptable !----------------------------------------------------------------------------- ! Module constants - double precision, parameter :: grav = 9.81 - double precision, parameter :: TOLERANCE = 1e-8 + real(prec), parameter :: grav = 9.81 + real(prec), parameter :: TOLERANCE = 1e-8 ! Module variables integer :: nlinks integer :: mxncomp integer :: maxTableLength integer :: nel - double precision :: so_llm - double precision :: z_g, bo_g, traps_g, tw_g, twcc_g, so_g, mann_g, manncc_g + real(prec) :: so_llm + real(prec) :: z_g, bo_g, traps_g, tw_g, twcc_g, so_g, mann_g, manncc_g integer, dimension(:,:), allocatable :: size_bathy - double precision, dimension(:,:), allocatable :: z - double precision, dimension(:,:,:), allocatable :: x_bathy, z_bathy, mann_bathy - double precision, dimension(:,:,:,:), allocatable :: xsec_tab + real(prec), dimension(:,:), allocatable :: z + real(prec), dimension(:,:,:), allocatable :: x_bathy, z_bathy, mann_bathy + real(prec), dimension(:,:,:,:), allocatable :: xsec_tab contains @@ -77,22 +77,22 @@ subroutine chxsec_lookuptable_calc(mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_ integer, intent(in) :: frnw_col integer, intent(in) :: mxnbathy_g integer, intent(in) :: nrow_chxsec_lookuptable - double precision, intent(in) :: so_lowerlimit_g + real(prec), intent(in) :: so_lowerlimit_g integer, dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g integer, dimension(mxncomp_g, nrch_g), intent(in) :: size_bathy_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: bo_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: traps_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: twcc_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: manncc_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g - double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: x_bathy_g - double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: z_bathy_g - double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: mann_bathy_g - double precision, dimension(11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g), intent(out) :: chxsec_lookuptable - double precision, dimension(mxncomp_g, nrch_g), intent(out) :: z_adj + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: bo_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: traps_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: twcc_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: manncc_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g + real(prec), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: x_bathy_g + real(prec), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: z_bathy_g + real(prec), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: mann_bathy_g + real(prec), dimension(11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g), intent(out) :: chxsec_lookuptable + real(prec), dimension(mxncomp_g, nrch_g), intent(out) :: z_adj ! Local variables integer :: ncomp @@ -109,16 +109,16 @@ subroutine chxsec_lookuptable_calc(mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_ integer, dimension(:), allocatable :: dmy_frj integer, dimension(:), allocatable :: mstem_frj integer, dimension(:,:), allocatable :: frnw_g - double precision :: timesDepth - !double precision :: t - double precision :: slope - double precision :: convey - double precision, dimension(:,:), allocatable :: leftBank - double precision, dimension(:,:), allocatable :: rightBank - double precision, dimension(:,:), allocatable :: skLeft - double precision, dimension(:,:), allocatable :: skMain - double precision, dimension(:,:), allocatable :: skRight - double precision, dimension(:,:), allocatable :: dx + real(prec) :: timesDepth + !real(prec) :: t + real(prec) :: slope + real(prec) :: convey + real(prec), dimension(:,:), allocatable :: leftBank + real(prec), dimension(:,:), allocatable :: rightBank + real(prec), dimension(:,:), allocatable :: skLeft + real(prec), dimension(:,:), allocatable :: skMain + real(prec), dimension(:,:), allocatable :: skRight + real(prec), dimension(:,:), allocatable :: dx !----------------------------------------------------------------------------- ! channel network data @@ -290,24 +290,24 @@ subroutine readXsection_natural_mann_vertices(idx_node, idx_reach, timesDepth) ! subroutine arguments integer, intent(in) :: idx_node, idx_reach - double precision, intent(in) :: timesDepth + real(prec), intent(in) :: timesDepth ! subroutine local variables integer :: i_area, i_find, num integer :: i1, i2 integer :: ic, iel, ii, ii2, iv, iel_start, iel_incr_start, iel_decr_start, ndmy - double precision :: el_min, el_max, el_range, el_incr, el_now, x1, y1, x2, y2, x_start, x_end - double precision :: f2m, cal_area, cal_peri, cal_topW - double precision :: mN_start, mN_end, cal_equiv_mann - double precision :: pos_slope, incr_rate, max_value + real(prec) :: el_min, el_max, el_range, el_incr, el_now, x1, y1, x2, y2, x_start, x_end + real(prec) :: f2m, cal_area, cal_peri, cal_topW + real(prec) :: mN_start, mN_end, cal_equiv_mann + real(prec) :: pos_slope, incr_rate, max_value integer, dimension(:), allocatable :: i_start, i_end - double precision, dimension(:), allocatable :: x_bathy_leftzero - double precision, dimension(:), allocatable :: xcs, ycs, manncs - double precision, dimension(:), allocatable :: el1, a1, peri1, redi1, equiv_mann - double precision, dimension(:), allocatable :: redi1All - double precision, dimension(:), allocatable :: conv1, tpW1 - double precision, dimension(:), allocatable :: newdKdA - double precision, dimension(:), allocatable :: compoundSKK, elev, dmyarr + real(prec), dimension(:), allocatable :: x_bathy_leftzero + real(prec), dimension(:), allocatable :: xcs, ycs, manncs + real(prec), dimension(:), allocatable :: el1, a1, peri1, redi1, equiv_mann + real(prec), dimension(:), allocatable :: redi1All + real(prec), dimension(:), allocatable :: conv1, tpW1 + real(prec), dimension(:), allocatable :: newdKdA + real(prec), dimension(:), allocatable :: compoundSKK, elev, dmyarr allocate(el1(nel), a1(nel), peri1(nel), redi1(nel), redi1All(nel)) allocate(equiv_mann(nel), conv1(nel), tpW1(nel)) @@ -558,7 +558,7 @@ subroutine readXsection_natural_mann_vertices(idx_node, idx_reach, timesDepth) deallocate(x_bathy_leftzero) contains - double precision function cal_dist_x_mann(x1, y1, x2, y2, mN) + real(prec) function cal_dist_x_mann(x1, y1, x2, y2, mN) implicit none @@ -568,16 +568,16 @@ double precision function cal_dist_x_mann(x1, y1, x2, y2, mN) !----------------------------------------------------- ! function arguments - double precision, intent(in) :: x1, y1, x2, y2, mN + real(prec), intent(in) :: x1, y1, x2, y2, mN ! function local variable - double precision :: dist + real(prec) :: dist dist = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + 1.e-32) cal_dist_x_mann = dist * mN**1.50 end function cal_dist_x_mann - double precision function cal_peri_x_mann(xx, yy, mN, n, i1, i2) + real(prec) function cal_peri_x_mann(xx, yy, mN, n, i1, i2) implicit none @@ -588,10 +588,10 @@ double precision function cal_peri_x_mann(xx, yy, mN, n, i1, i2) ! function arguments integer, intent(in) :: n, i1, i2 - double precision, intent(in) :: xx(n), yy(n), mN(n) + real(prec), intent(in) :: xx(n), yy(n), mN(n) ! function local variables integer :: i - double precision :: x1, x2, y1, y2, mN1, pxmN + real(prec) :: x1, x2, y1, y2, mN1, pxmN pxmN = 0.0 @@ -626,29 +626,29 @@ subroutine readXsection(k,lftBnkMann,rmanning_main,rgtBnkMann,leftBnkX_given,rgh ! subroutine arguments integer, intent(in) :: k, num_reach, lmxncomp, lnlinks - double precision, intent(in) :: rmanning_main,lftBnkMann,rgtBnkMann - double precision, intent(in) :: leftBnkX_given,rghtBnkX_given, timesDepth - double precision, dimension(lmxncomp, lnlinks), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g + real(prec), intent(in) :: rmanning_main,lftBnkMann,rgtBnkMann + real(prec), intent(in) :: leftBnkX_given,rghtBnkX_given, timesDepth + real(prec), dimension(lmxncomp, lnlinks), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g ! subroutine local variables integer :: i_area, i_find, i, j, jj, num integer :: i1, i2 integer :: mainChanStrt, mainChanEnd, kkk, startFound, endFound - double precision :: el_min, el_max, el_range, el_incr, el_now - double precision :: x1, y1, x2, y2, x_start, x_end - double precision :: waterElev, leftBnkX,rghtBnkX - double precision :: f2m, cal_area, cal_peri, cal_topW, diffAreaCenter - double precision :: compoundMann, el_min_1 - double precision :: leftBnkY, rghtBnkY,rmanning - double precision :: hbf + real(prec) :: el_min, el_max, el_range, el_incr, el_now + real(prec) :: x1, y1, x2, y2, x_start, x_end + real(prec) :: waterElev, leftBnkX,rghtBnkX + real(prec) :: f2m, cal_area, cal_peri, cal_topW, diffAreaCenter + real(prec) :: compoundMann, el_min_1 + real(prec) :: leftBnkY, rghtBnkY,rmanning + real(prec) :: hbf integer, dimension(:), allocatable :: i_start, i_end, totalNodes - double precision, dimension(:), allocatable :: xcs, ycs - double precision, dimension(:,:), allocatable :: el1, a1, peri1, redi1 - double precision, dimension(:), allocatable :: redi1All - double precision, dimension(:,:), allocatable :: conv1, tpW1, diffArea, newI1, diffPere - double precision, dimension(:), allocatable :: newdPdA, diffAreaAll, diffPereAll, newdKdA - double precision, dimension(:), allocatable :: compoundSKK, elev - double precision, dimension(:,:), allocatable :: allXcs, allYcs + real(prec), dimension(:), allocatable :: xcs, ycs + real(prec), dimension(:,:), allocatable :: el1, a1, peri1, redi1 + real(prec), dimension(:), allocatable :: redi1All + real(prec), dimension(:,:), allocatable :: conv1, tpW1, diffArea, newI1, diffPere + real(prec), dimension(:), allocatable :: newdPdA, diffAreaAll, diffPereAll, newdKdA + real(prec), dimension(:), allocatable :: compoundSKK, elev + real(prec), dimension(:,:), allocatable :: allXcs, allYcs allocate(el1(nel,3), a1(nel,3), peri1(nel,3), redi1(nel,3), redi1All(nel)) allocate(conv1(nel,3), tpW1(nel,3), diffArea(nel,3), newI1(nel,3), diffPere(nel,3)) @@ -948,7 +948,7 @@ subroutine readXsection(k,lftBnkMann,rmanning_main,rgtBnkMann,leftBnkX_given,rgh xsec_tab(9, j, k, num_reach) = newdKdA(j) xsec_tab(11,j, k, num_reach) = compoundSKK(j) end do - + z(k, num_reach) = el_min deallocate(el1, a1, peri1, redi1, redi1All) @@ -962,7 +962,7 @@ subroutine readXsection(k,lftBnkMann,rmanning_main,rgtBnkMann,leftBnkX_given,rgh end subroutine readXsection - double precision function cal_tri_area(el, x0, x1, y1) + real(prec) function cal_tri_area(el, x0, x1, y1) implicit none @@ -972,13 +972,13 @@ double precision function cal_tri_area(el, x0, x1, y1) !---------------------------------------- ! function arguments - doubleprecision, intent(in) :: el, x0, x1, y1 + real(prec), intent(in) :: el, x0, x1, y1 cal_tri_area = abs(0.5 * (x1 - x0) * (el - y1)) end function cal_tri_area - double precision function cal_trap_area(el, x1, y1, x2, y2) + real(prec) function cal_trap_area(el, x1, y1, x2, y2) implicit none @@ -987,13 +987,13 @@ double precision function cal_trap_area(el, x1, y1, x2, y2) ! calculate area of trapezoid !---------------------------------------- - doubleprecision, intent(in) :: el, x1, y1, x2, y2 + real(prec), intent(in) :: el, x1, y1, x2, y2 cal_trap_area = abs(0.5 * (x2 - x1) * (el - y1 + el - y2)) end function cal_trap_area - double precision function cal_multi_area(el, xx, yy, n, i1, i2) + real(prec) function cal_multi_area(el, xx, yy, n, i1, i2) implicit none @@ -1004,11 +1004,11 @@ double precision function cal_multi_area(el, xx, yy, n, i1, i2) ! function arguments integer, intent(in) :: n, i1, i2 - double precision, intent(in) :: el - double precision, intent(in) :: xx(n), yy(n) + real(prec), intent(in) :: el + real(prec), intent(in) :: xx(n), yy(n) ! function local variables integer :: i - double precision :: area, x1, x2, y1, y2 + real(prec) :: area, x1, x2, y1, y2 area = 0.0 @@ -1024,7 +1024,7 @@ double precision function cal_multi_area(el, xx, yy, n, i1, i2) endfunction cal_multi_area - double precision function cal_dist(x1, y1, x2, y2) + real(prec) function cal_dist(x1, y1, x2, y2) implicit none @@ -1034,13 +1034,13 @@ double precision function cal_dist(x1, y1, x2, y2) !------------------------------------------ ! function arguments - doubleprecision, intent(in) :: x1, y1, x2, y2 + real(prec), intent(in) :: x1, y1, x2, y2 cal_dist = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + 1.e-32) end function cal_dist - double precision function cal_perimeter(xx,yy,n,i1,i2) + real(prec) function cal_perimeter(xx,yy,n,i1,i2) implicit none @@ -1051,10 +1051,10 @@ double precision function cal_perimeter(xx,yy,n,i1,i2) ! function arguments integer, intent(in) :: n, i1, i2 - double precision, intent(in) :: xx(n), yy(n) + real(prec), intent(in) :: xx(n), yy(n) ! function local variables integer :: i - double precision :: p, x1, x2, y1, y2 + real(prec) :: p, x1, x2, y1, y2 p = 0. @@ -1069,5 +1069,4 @@ double precision function cal_perimeter(xx,yy,n,i1,i2) cal_perimeter=p end function cal_perimeter - -end module lookuptable \ No newline at end of file +end module lookuptable diff --git a/src/kernel/diffusive/diffusive.f90 b/src/kernel/diffusive/diffusive.f90 index 9fe60b6e9..2b4a779c9 100644 --- a/src/kernel/diffusive/diffusive.f90 +++ b/src/kernel/diffusive/diffusive.f90 @@ -1,6 +1,6 @@ !------------------------------------------------------------------------------- module diffusive - + use precis IMPLICIT NONE !----------------------------------------------------------------------------- @@ -16,8 +16,8 @@ module diffusive !----------------------------------------------------------------------------- ! Module constants - double precision, parameter :: grav = 9.81 - double precision, parameter :: TOLERANCE = 1e-8 + real(prec), parameter :: grav = 9.81 + real(prec), parameter :: TOLERANCE = 1e-8 ! Module variables integer :: nlinks @@ -38,37 +38,37 @@ module diffusive integer, dimension(:,:), allocatable :: routingNotChanged integer, dimension(:,:), allocatable :: frnw_g integer, dimension(:,:), allocatable :: size_bathy - double precision :: dtini, dxini, cfl, minDx, maxCelerity, theta - double precision :: C_llm, D_llm, D_ulm, DD_ulm, DD_llm, q_llm, so_llm - double precision :: frus2, minNotSwitchRouting, minNotSwitchRouting2 - double precision :: dt_qtrib - double precision :: z_g, bo_g, traps_g, tw_g, twcc_g, so_g, mann_g, manncc_g - double precision :: dmyi, dmyj - double precision :: dtini_min - double precision, dimension(:), allocatable :: eei, ffi, exi, fxi, qcx, diffusivity2, celerity2 - double precision, dimension(:), allocatable :: ini_y, ini_q - double precision, dimension(:), allocatable :: lowerLimitCount, higherLimitCount - double precision, dimension(:), allocatable :: elevTable, areaTable, skkkTable - double precision, dimension(:), allocatable :: pereTable, rediTable - double precision, dimension(:), allocatable :: convTable, topwTable - double precision, dimension(:), allocatable :: nwi1Table, dPdATable - double precision, dimension(:), allocatable :: ncompElevTable, ncompAreaTable - double precision, dimension(:), allocatable :: currentSquareDepth - double precision, dimension(:), allocatable :: tarr_qtrib, varr_qtrib - double precision, dimension(:), allocatable :: tarr_da, varr_da - double precision, dimension(:), allocatable :: area, depth, co, froud, courant - double precision, dimension(:,:), allocatable :: bo, dx - double precision, dimension(:,:), allocatable :: areap, qp, z, sk - double precision, dimension(:,:), allocatable :: celerity, diffusivity, qpx - double precision, dimension(:,:), allocatable :: pere, oldQ, newQ, oldArea, newArea, oldY, newY - double precision, dimension(:,:), allocatable :: volRemain - double precision, dimension(:,:), allocatable :: lateralFlow - double precision, dimension(:,:), allocatable :: dimensionless_Cr, dimensionless_Fo, dimensionless_Fi - double precision, dimension(:,:), allocatable :: dimensionless_Di, dimensionless_Fc, dimensionless_D - double precision, dimension(:,:), allocatable :: qtrib - double precision, dimension(:,:), allocatable :: usgs_da - double precision, dimension(:,:,:), allocatable :: x_bathy, z_bathy, mann_bathy - double precision, dimension(:,:,:,:), allocatable :: xsec_tab + real(prec) :: dtini, dxini, cfl, minDx, maxCelerity, theta + real(prec) :: C_llm, D_llm, D_ulm, DD_ulm, DD_llm, q_llm, so_llm + real(prec) :: frus2, minNotSwitchRouting, minNotSwitchRouting2 + real(prec) :: dt_qtrib + real(prec) :: z_g, bo_g, traps_g, tw_g, twcc_g, so_g, mann_g, manncc_g + real(prec) :: dmyi, dmyj + real(prec) :: dtini_min + real(prec), dimension(:), allocatable :: eei, ffi, exi, fxi, qcx, diffusivity2, celerity2 + real(prec), dimension(:), allocatable :: ini_y, ini_q + real(prec), dimension(:), allocatable :: lowerLimitCount, higherLimitCount + real(prec), dimension(:), allocatable :: elevTable, areaTable, skkkTable + real(prec), dimension(:), allocatable :: pereTable, rediTable + real(prec), dimension(:), allocatable :: convTable, topwTable + real(prec), dimension(:), allocatable :: nwi1Table, dPdATable + real(prec), dimension(:), allocatable :: ncompElevTable, ncompAreaTable + real(prec), dimension(:), allocatable :: currentSquareDepth + real(prec), dimension(:), allocatable :: tarr_qtrib, varr_qtrib + real(prec), dimension(:), allocatable :: tarr_da, varr_da + real(prec), dimension(:), allocatable :: area, depth, co, froud, courant + real(prec), dimension(:,:), allocatable :: bo, dx + real(prec), dimension(:,:), allocatable :: areap, qp, z, sk + real(prec), dimension(:,:), allocatable :: celerity, diffusivity, qpx + real(prec), dimension(:,:), allocatable :: pere, oldQ, newQ, oldArea, newArea, oldY, newY + real(prec), dimension(:,:), allocatable :: volRemain + real(prec), dimension(:,:), allocatable :: lateralFlow + real(prec), dimension(:,:), allocatable :: dimensionless_Cr, dimensionless_Fo, dimensionless_Fi + real(prec), dimension(:,:), allocatable :: dimensionless_Di, dimensionless_Fc, dimensionless_D + real(prec), dimension(:,:), allocatable :: qtrib + real(prec), dimension(:,:), allocatable :: usgs_da + real(prec), dimension(:,:,:), allocatable :: x_bathy, z_bathy, mann_bathy + real(prec), dimension(:,:,:,:), allocatable :: xsec_tab contains @@ -132,32 +132,32 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt integer, dimension(nrch_g), intent(in) :: usgs_da_reach_g integer, dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g integer, dimension(mxncomp_g, nrch_g), intent(in) :: size_bathy_g - double precision, dimension(paradim ), intent(in) :: para_ar_g - double precision, dimension(:) , intent(in) :: timestep_ar_g(10) - double precision, dimension(nts_db_g), intent(in) :: dbcd_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: bo_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: traps_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: twcc_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: manncc_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: rdx_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: iniq - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: so_ar_g - double precision, dimension(mxncomp_g, nrch_g), intent(in) :: z_thalweg_g - double precision, dimension(nts_ub_g, nrch_g), intent(in) :: ubcd_g - double precision, dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g - double precision, dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g - double precision, dimension(cwnrow_g, cwncol_g), intent(in) :: crosswalk_g - double precision, dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in ) :: qlat_g - double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: x_bathy_g - double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: z_bathy_g - double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: mann_bathy_g - double precision, dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g - double precision, dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: elv_ev_g - double precision, dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: depth_ev_g + real(prec), dimension(paradim ), intent(in) :: para_ar_g + real(prec), dimension(:) , intent(in) :: timestep_ar_g(10) + real(prec), dimension(nts_db_g), intent(in) :: dbcd_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: bo_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: traps_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: twcc_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: manncc_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: rdx_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: iniq + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: so_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: z_thalweg_g + real(prec), dimension(nts_ub_g, nrch_g), intent(in) :: ubcd_g + real(prec), dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g + real(prec), dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g + real(prec), dimension(cwnrow_g, cwncol_g), intent(in) :: crosswalk_g + real(prec), dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in ) :: qlat_g + real(prec), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: x_bathy_g + real(prec), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: z_bathy_g + real(prec), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: mann_bathy_g + real(prec), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g + real(prec), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: elv_ev_g + real(prec), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: depth_ev_g ! Local variables integer :: ncomp @@ -179,48 +179,49 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt integer :: ri, rj, oi, oj integer :: nlnk, lnk integer :: nusrch + integer :: n_decimal_places integer, dimension(:), allocatable :: dmy_frj integer, dimension(:,:), allocatable :: flag_lfrac - double precision :: x - double precision :: saveInterval, width - !double precision :: maxCourant - double precision :: dtini_given, dtini_divisor - double precision :: timesDepth - double precision :: t - double precision :: tfin - double precision :: t0 - double precision :: q_sk_multi - double precision :: maxCelDx - double precision :: slope - double precision :: y_norm - double precision :: temp - double precision :: dt_ql - double precision :: dt_ub - double precision :: dt_db - double precision :: dt_da - double precision :: wdepth - double precision :: q_usrch - double precision :: tf0 - double precision :: convey - double precision :: equiv_one - double precision :: slopeQ, intcQ, slopeE, intcE, slopeD, intcD - double precision :: lfrac, dst_lnk, dst_top, dst_btm - double precision :: mindepth_nstab - double precision, dimension(:), allocatable :: tarr_ql - double precision, dimension(:), allocatable :: varr_ql - double precision, dimension(:), allocatable :: tarr_ub - double precision, dimension(:), allocatable :: varr_ub - double precision, dimension(:), allocatable :: tarr_db - double precision, dimension(:), allocatable :: varr_db - double precision, dimension(:), allocatable :: dbcd ! temporary - double precision, dimension(:,:), allocatable :: leftBank - double precision, dimension(:,:), allocatable :: rightBank - double precision, dimension(:,:), allocatable :: skLeft - double precision, dimension(:,:), allocatable :: skMain - double precision, dimension(:,:), allocatable :: skRight - double precision, dimension(:,:), allocatable :: used_lfrac - double precision, dimension(:,:,:), allocatable :: temp_q_ev_g - double precision, dimension(:,:,:), allocatable :: temp_elv_ev_g + real(prec) :: x + real(prec) :: saveInterval, width + !real(prec) :: maxCourant + real(prec) :: dtini_given, dtini_divisor + real(prec) :: timesDepth + real(prec) :: t + real(prec) :: tfin + real(prec) :: t0 + real(prec) :: q_sk_multi + real(prec) :: maxCelDx + real(prec) :: slope + real(prec) :: y_norm + real(prec) :: temp + real(prec) :: dt_ql + real(prec) :: dt_ub + real(prec) :: dt_db + real(prec) :: dt_da + real(prec) :: wdepth + real(prec) :: q_usrch + real(prec) :: tf0 + real(prec) :: convey + real(prec) :: equiv_one + real(prec) :: slopeQ, intcQ, slopeE, intcE, slopeD, intcD + real(prec) :: lfrac, dst_lnk, dst_top, dst_btm + real(prec) :: mindepth_nstab + real(prec), dimension(:), allocatable :: tarr_ql + real(prec), dimension(:), allocatable :: varr_ql + real(prec), dimension(:), allocatable :: tarr_ub + real(prec), dimension(:), allocatable :: varr_ub + real(prec), dimension(:), allocatable :: tarr_db + real(prec), dimension(:), allocatable :: varr_db + real(prec), dimension(:), allocatable :: dbcd ! temporary + real(prec), dimension(:,:), allocatable :: leftBank + real(prec), dimension(:,:), allocatable :: rightBank + real(prec), dimension(:,:), allocatable :: skLeft + real(prec), dimension(:,:), allocatable :: skMain + real(prec), dimension(:,:), allocatable :: skRight + real(prec), dimension(:,:), allocatable :: used_lfrac + real(prec), dimension(:,:,:), allocatable :: temp_q_ev_g + real(prec), dimension(:,:,:), allocatable :: temp_elv_ev_g !----------------------------------------------------------------------------- ! Time domain parameters @@ -257,7 +258,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt DD_llm = para_ar_g(5) ! lower limit of dimensionless diffusivity (default -15) DD_ulm = para_ar_g(6) ! upper limit of dimensionless diffusivity (default: -10.0) q_llm = para_ar_g(8) ! lower limit of discharge (default: 0.02831 cms) - so_llm = para_ar_g(9) ! lower limit of channel bed slope (default: 0.0001) + so_llm = para_ar_g(9) ! lower limit of channel bed slope (default: 0.00001) theta = para_ar_g(10) ! weight for computing 2nd derivative: ! 0: explicit, 1: implicit (default: 1.0) dsbc_option = para_ar_g(11) ! downstream water depth boundary condition option 1:given water depth data, 2:normal depth @@ -341,7 +342,6 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt allocate(usgs_da(nts_da, nlinks)) allocate(dbcd(nts_db_g)) allocate(temp_q_ev_g(ntss_ev_g, mxncomp_g, nrch_g), temp_elv_ev_g(ntss_ev_g, mxncomp_g, nrch_g)) - !-------------------------------------------------------------------------------------------- frnw_g = frnw_ar_g ! network mapping matrix @@ -350,8 +350,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt usgs_da = usgs_da_g ! contains usgs data at a related reach !----------------------------------------------------------------------------- - ! variable initializations - + ! variable initializations routingNotChanged = 0 applyNaturalSection = 1 x = 0.0 @@ -359,7 +358,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt newY = -999 t = t0*60.0 ! [min] q_sk_multi = 1.0 - oldQ = iniq + oldQ = iniq newQ = oldQ qp = oldQ dimensionless_Cr = -999 @@ -374,10 +373,10 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt depth_ev_g = 0.0 temp_q_ev_g = 0.0 temp_elv_ev_g = 0.0 + !----------------------------------------------------------------------------- ! Identify mainstem reaches and list their ids in an array - ! Create a dummy array containing mainstem reaches where diffusive wave is applied. nmstem_rch = 0 do j = 1, nlinks !if (frnw_g(j,3) >= 1) then ! mainstem reach identification @@ -394,10 +393,10 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt mstem_frj(jm) = dmy_frj(jm) end do deallocate(dmy_frj) - + !----------------------------------------------------------------------------- ! create dx array from dx_ar_g and determine minimum dx. - + dx = 0. minDx = 1e10 @@ -528,9 +527,17 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt !----------------------------------------------------------------------------- ! Initialize water surface elevation, channel area, and volume + do j = 1, nlinks + ncomp = frnw_g(j, 1) ! number of nodes in reach j + do i=1, ncomp + oldQ(i,j) = max(oldQ(i,j), q_llm) + qp(i,j) = oldQ(i,j) + enddo + enddo + do jm = nmstem_rch, 1, -1 j = mstem_frj(jm) ! reach index - ncomp = frnw_g(j, 1) ! number of nodes in reach j + ncomp = frnw_g(j, 1) ! number of nodes in reach j if (frnw_g(j, 2) < 0) then ! Initial depth at bottom node of tail water reach @@ -563,22 +570,22 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt ! compute newY(i, j) for i=1, ncomp-1 with the given newY(ncomp, j) ! ** At initial time, oldY(i,j) values at i < ncomp, used in subroutine rtsafe, are not defined. ! ** So, let's assume the depth values are all nodes equal to depth at the bottom node. - wdepth = newY(ncomp, j) - z(ncomp, j) + oldY(ncomp, j) = newY(ncomp, j) + wdepth = newY(ncomp, j) - z(ncomp, j) do i = 1, ncomp -1 oldY(i,j) = wdepth + z(i, j) end do call mesh_diffusive_backward(dtini_given, t0, t, tfin, saveInterval, j) - + do i = 1,ncomp ! copy computed initial depths to initial depth array for first timestep oldY(i, j) = newY(i, j) ! Check that node elevation is not lower than bottom node. ! If node elevation is lower than bottom node elevation, correct. - if (oldY(i, j) .lt. oldY(ncomp, nlinks)) oldY(i, j) = oldY(ncomp, nlinks) - end do - + !if (oldY(i, j) .lt. oldY(ncomp, nlinks)) oldY(i, j) = oldY(ncomp, nlinks) + end do end do !----------------------------------------------------------------------------- @@ -621,14 +628,12 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt frus2 = 9999. notSwitchRouting = 0 minNotSwitchRouting = 10000 -! minNotSwitchRouting2 = 000 timestep = 0 ts_ev = 1 t = t0 * 60.0 !----------------------------------------------------------------------------- ! Ordered network routing computations - do while ( t < tfin * 60.) timestep = timestep + 1 ! advance timestep !+------------------------------------------------------------------------- @@ -639,12 +644,12 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt do jm = 1, nmstem_rch ! loop over mainstem reaches [upstream-to-downstream] j = mstem_frj(jm) ! reach index ncomp = frnw_g(j,1) ! number of nodes in reach j - + ! Calculate the duration of this timestep (dtini) ! Timestep duration is selected to maintain numerical stability if (j == mstem_frj(1)) then call calculateDT(t0, t, saveInterval, cfl, tfin, maxCelDx, dtini_given) - end if + end if ! estimate lateral flow at current time t do i = 1, ncomp - 1 @@ -674,14 +679,15 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt do n = 1, nts_qtrib_g varr_qtrib(n) = qtrib_g(n, usrchj) end do - tf0 = t + dtini / 60. + tf0 = t !+ dtini / 60 q_usrch = intp_y(nts_qtrib_g, tarr_qtrib, varr_qtrib, tf0) end if + ! add upstream flows to reach head newQ(1,j)= newQ(1,j) + q_usrch end do else - + ! no stream reaches at the upstream of the reach (frnw_g(j,3)==0) such as a headwater reach among tributaries newQ(1,j)= 0.0 end if @@ -772,7 +778,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt !end if end do - ! Advance model time + ! Advance model time [minute] t = t + dtini/60. ! Calculate dimensionless numbers for each reach @@ -782,13 +788,13 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt !end do ! diffusive wave simulation time print - if (mod(t,30.)==0.) then + !if (mod(t,30.)==0.) then print*, "diffusive simulation time in minute=", t - endif + !endif ! write results to output arrays if ( (mod((t - t0 * 60.) * 60., saveInterval) <= TOLERANCE) .or. (t == tfin * 60.)) then - do jm = 1, nmstem_rch + do jm = 1, nmstem_rch j = mstem_frj(jm) ncomp = frnw_g(j, 1) do i = 1, ncomp @@ -796,7 +802,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt elv_ev_g(ts_ev + 1, i, j) = newY(i, j) depth_ev_g(ts_ev + 1, i, j) = elv_ev_g(ts_ev + 1, i, j) - z(i, j) end do - + !* water elevation for tributaries flowing into the mainstem in the middle or at the upper end do k = 1, frnw_g(j, 3) usrchj = frnw_g(j, 3 + k) @@ -848,7 +854,8 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt !--------------------------------------------------------------------------- ! map routing result from refactored hydrofabric to unrefactored hydrofabric - if (cwnrow_g > 0) then + if (1==2) then + if (cwnrow_g > 0) then allocate(used_lfrac(mxncomp, nlinks)) allocate(flag_lfrac(mxncomp, nlinks)) equiv_one = 0.99 @@ -914,11 +921,12 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt flag_lfrac(oi, oj) = 0 endif end do - end do end do - deallocate(used_lfrac) - deallocate(flag_lfrac) - endif + end do + deallocate(used_lfrac) + deallocate(flag_lfrac) + endif + endif deallocate(frnw_g) deallocate(area, bo, pere, areap, qp, z, depth, sk, co, dx) @@ -961,13 +969,13 @@ subroutine calculateDT(initialTime, time, saveInterval, & !----------------------------------------------------------------------------- ! Subroutine arguments - double precision, intent(in) :: initialTime ! [sec] - double precision, intent(in) :: time ! [hrs] - double precision, intent(in) :: saveInterval ! [sec] - double precision, intent(in) :: tfin ! [hrs] - double precision, intent(in) :: given_dt ! [sec] - double precision, intent(in) :: maxAllowCourantNo ! = cfl - double precision, intent(in) :: max_C_dx ! = maxCelDx + real(prec), intent(in) :: initialTime ! [sec] + real(prec), intent(in) :: time ! [hrs] + real(prec), intent(in) :: saveInterval ! [sec] + real(prec), intent(in) :: tfin ! [hrs] + real(prec), intent(in) :: given_dt ! [sec] + real(prec), intent(in) :: maxAllowCourantNo ! = cfl + real(prec), intent(in) :: max_C_dx ! = maxCelDx ! Local variables integer :: a @@ -1015,29 +1023,29 @@ subroutine calc_dimensionless_numbers(j) ! Local variables integer :: i integer :: ncomp - double precision :: wl_us - double precision :: depth_us - double precision :: q_us - double precision :: v_us - double precision :: pere_us - double precision :: r_us - double precision :: sk_us - double precision :: ch_us - double precision :: wl_ds - double precision :: depth_ds - double precision :: q_ds - double precision :: v_ds - double precision :: pere_ds - double precision :: r_ds - double precision :: sk_ds - double precision :: ch_ds - double precision :: ch_star_avg - double precision :: channel_length - double precision :: avg_celerity - double precision :: avg_velocity - double precision :: avg_depth - double precision :: maxValue - double precision :: dimlessWaveLength + real(prec) :: wl_us + real(prec) :: depth_us + real(prec) :: q_us + real(prec) :: v_us + real(prec) :: pere_us + real(prec) :: r_us + real(prec) :: sk_us + real(prec) :: ch_us + real(prec) :: wl_ds + real(prec) :: depth_ds + real(prec) :: q_ds + real(prec) :: v_ds + real(prec) :: pere_ds + real(prec) :: r_ds + real(prec) :: sk_ds + real(prec) :: ch_ds + real(prec) :: ch_star_avg + real(prec) :: channel_length + real(prec) :: avg_celerity + real(prec) :: avg_velocity + real(prec) :: avg_depth + real(prec) :: maxValue + real(prec) :: dimlessWaveLength !----------------------------------------------------------------------------- maxValue = 1e7 @@ -1125,27 +1133,27 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j) ! Subroutine Arguments integer, intent(in) :: j - double precision, intent(in) :: dtini_given - double precision, intent(in) :: t0 - double precision, intent(in) :: t - double precision, intent(in) :: tfin - double precision, intent(in) :: saveInterval + real(prec), intent(in) :: dtini_given + real(prec), intent(in) :: t0 + real(prec), intent(in) :: t + real(prec), intent(in) :: tfin + real(prec), intent(in) :: saveInterval ! Local variables integer :: ncomp integer :: i, irow, flag_da, n - double precision :: a1, a2, a3, a4 - double precision :: b1, b2, b3, b4 - double precision :: dd1, dd2, dd3, dd4 - double precision :: h1, h2, h3, h4 - double precision :: allqlat - double precision :: qy, qxy, qxxy, qxxxy - double precision :: ppi, qqi, rri, ssi, sxi - double precision :: cour, cour2 - double precision :: alpha - double precision :: currentQ - double precision :: eei_ghost, ffi_ghost, exi_ghost - double precision :: fxi_ghost, qp_ghost, qpx_ghost + real(prec) :: a1, a2, a3, a4 + real(prec) :: b1, b2, b3, b4 + real(prec) :: dd1, dd2, dd3, dd4 + real(prec) :: h1, h2, h3, h4 + real(prec) :: allqlat + real(prec) :: qy, qxy, qxxy, qxxxy + real(prec) :: ppi, qqi, rri, ssi, sxi + real(prec) :: cour, cour2 + real(prec) :: alpha + real(prec) :: currentQ + real(prec) :: eei_ghost, ffi_ghost, exi_ghost + real(prec) :: fxi_ghost, qp_ghost, qpx_ghost !----------------------------------------------------------------------------- !* change 20210228: All qlat to a river reach is applied to the u/s boundary !* Note: lateralFlow(1,j) is already added to the boundary @@ -1224,7 +1232,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j) ffi(i) = ( ssi - ppi * ffi(i-1) ) / ( ppi * eei(i-1) + qqi ) exi(i) = -1.0 * rri / ( ppi * exi(i-1) + qqi ) fxi(i) = ( sxi - ppi * fxi(i-1) ) / ( ppi * exi(i-1) + qqi ) - + end do ! Ghost point calculation @@ -1310,14 +1318,14 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j) do i = ncomp-1, 1, -1 qp(i, j) = eei(i) * qp(i+1, j) + ffi(i) qpx(i, j) = exi(i) * qpx(i+1, j) + fxi(i) - end do + end do ! when a reach hasn't been applied to DA ! if ((usgs_da_reach(j) == 0).or.(flag_da == 0)) then qp(1, j) = newQ(1, j) qp(1, j) = qp(1, j) + allqlat ! endif - + do i = 1, ncomp if (abs(qp(i, j)) < q_llm) then qp(i, j) = q_llm @@ -1328,6 +1336,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j) do i = 1, ncomp newQ(i, j) = qp(i, j) end do + ! ============== DEBUG to find unstable flow calcls ================= ! do i= ncomp, 1, -1 ! if (abs(qp(i,j)) .gt. 2E4) then @@ -1354,7 +1363,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j) end subroutine mesh_diffusive_forward - subroutine mesh_diffusive_backward(dtini_given, t0, t, tfin, saveInterval, j) ! leftBank, rightBank) + subroutine mesh_diffusive_backward(dtini_given, t0, t, tfin, saveInterval, j) ! leftBank, rightBank IMPLICIT NONE @@ -1375,23 +1384,23 @@ subroutine mesh_diffusive_backward(dtini_given, t0, t, tfin, saveInterval, j) ! ! Subroutine arguments integer, intent(in) :: j - double precision, intent(in) :: dtini_given - double precision, intent(in) :: t0 - double precision, intent(in) :: t - double precision, intent(in) :: tfin - double precision, intent(in) :: saveInterval + real(prec), intent(in) :: dtini_given + real(prec), intent(in) :: t0 + real(prec), intent(in) :: t + real(prec), intent(in) :: tfin + real(prec), intent(in) :: saveInterval ! Subroutine variables integer :: depthCalOk(mxncomp) integer :: i, ncomp integer :: jj - double precision :: xt - double precision :: q_sk_multi, sfi - double precision :: vel, currentQ - double precision :: S_ncomp - double precision :: tempDepthi_1 - double precision :: Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds - double precision :: C_ulm + real(prec) :: xt + real(prec) :: q_sk_multi, sfi + real(prec) :: vel, currentQ + real(prec) :: S_ncomp + real(prec) :: tempDepthi_1 + real(prec) :: Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds + real(prec) :: C_ulm !----------------------------------------------------------------------------- ncomp = frnw_g(j, 1) @@ -1547,9 +1556,10 @@ subroutine mesh_diffusive_backward(dtini_given, t0, t, tfin, saveInterval, j) ! diffusivity(1:ncomp, j) = sum(diffusivity2(1:ncomp)) / ncomp do i = 1, ncomp - if (diffusivity(i, j) > D_ulm) diffusivity(i, j) = D_ulm !!! Test - if (diffusivity(i, j) < D_llm) diffusivity(i, j) = D_llm !!! Test + if (diffusivity(i, j) > D_ulm) diffusivity(i, j) = D_ulm + if (diffusivity(i, j) < D_llm) diffusivity(i, j) = D_llm end do + end subroutine mesh_diffusive_backward function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds) @@ -1572,16 +1582,16 @@ function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds) ! Subroutine arguments integer, intent(in) :: i, j - double precision, intent(in) :: Q_cur, Q_ds, z_cur, z_ds, y_ds + real(prec), intent(in) :: Q_cur, Q_ds, z_cur, z_ds, y_ds ! Subroutine local variable integer, parameter :: maxit = 40 - double precision, parameter :: xacc = 1e-4 + real(prec), parameter :: xacc = 1e-4 integer :: iter integer :: xcolID, ycolID - double precision :: x1, x2, df, dxx, dxold, f, fh, fl, temp, xh, xl - double precision :: y_norm, y_ulm_multi, y_llm_multi, elv_norm, y_old - double precision :: rtsafe + real(prec) :: x1, x2, df, dxx, dxold, f, fh, fl, temp, xh, xl + real(prec) :: y_norm, y_ulm_multi, y_llm_multi, elv_norm, y_old + real(prec) :: rtsafe y_ulm_multi = 2.0 y_llm_multi = 0.1 @@ -1675,13 +1685,13 @@ subroutine funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds, f, df) ! subroutine arguments integer, intent(in) :: i, j - double precision, intent(in) :: Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds - double precision, intent(out) :: f, df + real(prec), intent(in) :: Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds + real(prec), intent(out) :: f, df ! subroutine local variables integer :: xcolID, ycolID - double precision :: elv_cur, elv_ds, conv_cur, conv_ds, sf_cur, sf_ds, slope - double precision :: dKdA_cur, topw_cur + real(prec) :: elv_cur, elv_ds, conv_cur, conv_ds, sf_cur, sf_ds, slope + real(prec) :: dKdA_cur, topw_cur xcolID = 1 ! f(y_cur): function of y at the current node @@ -1710,7 +1720,7 @@ subroutine funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds, f, df) end subroutine funcd_diffdepth - double precision function intp_xsec_tab(i, j, nrow, xcolID, ycolID, x) + real(prec) function intp_xsec_tab(i, j, nrow, xcolID, ycolID, x) implicit none @@ -1724,12 +1734,12 @@ double precision function intp_xsec_tab(i, j, nrow, xcolID, ycolID, x) ! subroutine arguments integer , intent(in) :: i, j, nrow, xcolID, ycolID - double precision, intent(in) :: x + real(prec), intent(in) :: x ! subroutine local variables integer :: irow - double precision :: x1, y1, x2, y2, y - double precision, dimension(nrow) :: xarr, yarr + real(prec) :: x1, y1, x2, y2, y + real(prec), dimension(nrow) :: xarr, yarr xarr = xsec_tab(xcolID, 1:nrow, i, j) yarr = xsec_tab(ycolID, 1:nrow, i, j) @@ -1770,24 +1780,24 @@ subroutine readXsection_natural_mann_vertices(idx_node, idx_reach, timesDepth) ! subroutine arguments integer, intent(in) :: idx_node, idx_reach - double precision, intent(in) :: timesDepth + real(prec), intent(in) :: timesDepth ! subroutine local variables integer :: i_area, i_find, num integer :: i1, i2 integer :: ic, iel, ii, ii2, iv, iel_start, iel_incr_start, iel_decr_start, ndmy - double precision :: el_min, el_max, el_range, el_incr, el_now, x1, y1, x2, y2, x_start, x_end - double precision :: f2m, cal_area, cal_peri, cal_topW - double precision :: mN_start, mN_end, cal_equiv_mann - double precision :: pos_slope, incr_rate, max_value + real(prec) :: el_min, el_max, el_range, el_incr, el_now, x1, y1, x2, y2, x_start, x_end + real(prec) :: f2m, cal_area, cal_peri, cal_topW + real(prec) :: mN_start, mN_end, cal_equiv_mann + real(prec) :: pos_slope, incr_rate, max_value integer, dimension(:), allocatable :: i_start, i_end - double precision, dimension(:), allocatable :: x_bathy_leftzero - double precision, dimension(:), allocatable :: xcs, ycs, manncs - double precision, dimension(:), allocatable :: el1, a1, peri1, redi1, equiv_mann - double precision, dimension(:), allocatable :: redi1All - double precision, dimension(:), allocatable :: conv1, tpW1 - double precision, dimension(:), allocatable :: newdKdA - double precision, dimension(:), allocatable :: compoundSKK, elev, dmyarr + real(prec), dimension(:), allocatable :: x_bathy_leftzero + real(prec), dimension(:), allocatable :: xcs, ycs, manncs + real(prec), dimension(:), allocatable :: el1, a1, peri1, redi1, equiv_mann + real(prec), dimension(:), allocatable :: redi1All + real(prec), dimension(:), allocatable :: conv1, tpW1 + real(prec), dimension(:), allocatable :: newdKdA + real(prec), dimension(:), allocatable :: compoundSKK, elev, dmyarr allocate(el1(nel), a1(nel), peri1(nel), redi1(nel), redi1All(nel)) allocate(equiv_mann(nel), conv1(nel), tpW1(nel)) @@ -2038,7 +2048,7 @@ subroutine readXsection_natural_mann_vertices(idx_node, idx_reach, timesDepth) deallocate(x_bathy_leftzero) contains - double precision function cal_dist_x_mann(x1, y1, x2, y2, mN) + real(prec) function cal_dist_x_mann(x1, y1, x2, y2, mN) implicit none @@ -2048,16 +2058,16 @@ double precision function cal_dist_x_mann(x1, y1, x2, y2, mN) !----------------------------------------------------- ! function arguments - double precision, intent(in) :: x1, y1, x2, y2, mN + real(prec), intent(in) :: x1, y1, x2, y2, mN ! function local variable - double precision :: dist + real(prec) :: dist dist = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + 1.e-32) cal_dist_x_mann = dist * mN**1.50 end function cal_dist_x_mann - double precision function cal_peri_x_mann(xx, yy, mN, n, i1, i2) + real(prec) function cal_peri_x_mann(xx, yy, mN, n, i1, i2) implicit none @@ -2068,10 +2078,10 @@ double precision function cal_peri_x_mann(xx, yy, mN, n, i1, i2) ! function arguments integer, intent(in) :: n, i1, i2 - double precision, intent(in) :: xx(n), yy(n), mN(n) + real(prec), intent(in) :: xx(n), yy(n), mN(n) ! function local variables integer :: i - double precision :: x1, x2, y1, y2, mN1, pxmN + real(prec) :: x1, x2, y1, y2, mN1, pxmN pxmN = 0.0 @@ -2106,30 +2116,30 @@ subroutine readXsection(k,lftBnkMann,rmanning_main,rgtBnkMann,leftBnkX_given,rgh ! subroutine arguments integer, intent(in) :: k, num_reach, lmxncomp, lnlinks - double precision, intent(in) :: rmanning_main,lftBnkMann,rgtBnkMann - double precision, intent(in) :: leftBnkX_given,rghtBnkX_given, timesDepth - double precision, dimension(lmxncomp, lnlinks), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g + real(prec), intent(in) :: rmanning_main,lftBnkMann,rgtBnkMann + real(prec), intent(in) :: leftBnkX_given,rghtBnkX_given, timesDepth + real(prec), dimension(lmxncomp, lnlinks), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g ! subroutine local variables integer :: i_area, i_find, i, j, jj, num integer :: i1, i2 integer :: mainChanStrt, mainChanEnd, kkk, startFound, endFound - double precision :: el_min, el_max, el_range, el_incr, el_now - double precision :: x1, y1, x2, y2, x_start, x_end - double precision :: waterElev, leftBnkX,rghtBnkX - double precision :: f2m, cal_area, cal_peri, cal_topW, diffAreaCenter - double precision :: compoundMann, el_min_1 - double precision :: leftBnkY, rghtBnkY,rmanning - double precision :: hbf + real(prec) :: el_min, el_max, el_range, el_incr, el_now + real(prec) :: x1, y1, x2, y2, x_start, x_end + real(prec) :: waterElev, leftBnkX,rghtBnkX + real(prec) :: f2m, cal_area, cal_peri, cal_topW, diffAreaCenter + real(prec) :: compoundMann, el_min_1 + real(prec) :: leftBnkY, rghtBnkY,rmanning + real(prec) :: hbf integer, dimension(:), allocatable :: i_start, i_end, totalNodes - double precision, dimension(:), allocatable :: xcs, ycs - double precision, dimension(:,:), allocatable :: el1, a1, peri1, redi1 - double precision, dimension(:), allocatable :: redi1All - double precision, dimension(:,:), allocatable :: conv1, tpW1, diffArea, newI1, diffPere - double precision, dimension(:), allocatable :: newdPdA, diffAreaAll, diffPereAll, newdKdA - double precision, dimension(:), allocatable :: compoundSKK, elev - double precision, dimension(:,:), allocatable :: allXcs, allYcs - + real(prec), dimension(:), allocatable :: xcs, ycs + real(prec), dimension(:,:), allocatable :: el1, a1, peri1, redi1 + real(prec), dimension(:), allocatable :: redi1All + real(prec), dimension(:,:), allocatable :: conv1, tpW1, diffArea, newI1, diffPere + real(prec), dimension(:), allocatable :: newdPdA, diffAreaAll, diffPereAll, newdKdA + real(prec), dimension(:), allocatable :: compoundSKK, elev + real(prec), dimension(:,:), allocatable :: allXcs, allYcs + allocate(el1(nel,3), a1(nel,3), peri1(nel,3), redi1(nel,3), redi1All(nel)) allocate(conv1(nel,3), tpW1(nel,3), diffArea(nel,3), newI1(nel,3), diffPere(nel,3)) allocate(newdPdA(nel), diffAreaAll(nel), diffPereAll(nel), newdKdA(nel)) ! change Nazmul 20210601 @@ -2428,7 +2438,7 @@ subroutine readXsection(k,lftBnkMann,rmanning_main,rgtBnkMann,leftBnkX_given,rgh xsec_tab(9, j, k, num_reach) = newdKdA(j) xsec_tab(11,j, k, num_reach) = compoundSKK(j) end do - + z(k, num_reach) = el_min deallocate(el1, a1, peri1, redi1, redi1All) @@ -2442,7 +2452,7 @@ subroutine readXsection(k,lftBnkMann,rmanning_main,rgtBnkMann,leftBnkX_given,rgh end subroutine readXsection - double precision function cal_tri_area(el, x0, x1, y1) + real(prec) function cal_tri_area(el, x0, x1, y1) implicit none @@ -2452,13 +2462,13 @@ double precision function cal_tri_area(el, x0, x1, y1) !---------------------------------------- ! function arguments - doubleprecision, intent(in) :: el, x0, x1, y1 + real(prec), intent(in) :: el, x0, x1, y1 cal_tri_area = abs(0.5 * (x1 - x0) * (el - y1)) end function cal_tri_area - double precision function cal_trap_area(el, x1, y1, x2, y2) + real(prec) function cal_trap_area(el, x1, y1, x2, y2) implicit none @@ -2467,13 +2477,13 @@ double precision function cal_trap_area(el, x1, y1, x2, y2) ! calculate area of trapezoid !---------------------------------------- - doubleprecision, intent(in) :: el, x1, y1, x2, y2 + real(prec), intent(in) :: el, x1, y1, x2, y2 cal_trap_area = abs(0.5 * (x2 - x1) * (el - y1 + el - y2)) end function cal_trap_area - double precision function cal_multi_area(el, xx, yy, n, i1, i2) + real(prec) function cal_multi_area(el, xx, yy, n, i1, i2) implicit none @@ -2484,11 +2494,11 @@ double precision function cal_multi_area(el, xx, yy, n, i1, i2) ! function arguments integer, intent(in) :: n, i1, i2 - double precision, intent(in) :: el - double precision, intent(in) :: xx(n), yy(n) + real(prec), intent(in) :: el + real(prec), intent(in) :: xx(n), yy(n) ! function local variables integer :: i - double precision :: area, x1, x2, y1, y2 + real(prec) :: area, x1, x2, y1, y2 area = 0.0 @@ -2504,7 +2514,7 @@ double precision function cal_multi_area(el, xx, yy, n, i1, i2) endfunction cal_multi_area - double precision function cal_dist(x1, y1, x2, y2) + real(prec) function cal_dist(x1, y1, x2, y2) implicit none @@ -2514,13 +2524,13 @@ double precision function cal_dist(x1, y1, x2, y2) !------------------------------------------ ! function arguments - doubleprecision, intent(in) :: x1, y1, x2, y2 + real(prec), intent(in) :: x1, y1, x2, y2 cal_dist = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + 1.e-32) end function cal_dist - double precision function cal_perimeter(xx,yy,n,i1,i2) + real(prec) function cal_perimeter(xx,yy,n,i1,i2) implicit none @@ -2531,10 +2541,10 @@ double precision function cal_perimeter(xx,yy,n,i1,i2) ! function arguments integer, intent(in) :: n, i1, i2 - double precision, intent(in) :: xx(n), yy(n) + real(prec), intent(in) :: xx(n), yy(n) ! function local variables integer :: i - double precision :: p, x1, x2, y1, y2 + real(prec) :: p, x1, x2, y1, y2 p = 0. @@ -2559,8 +2569,8 @@ subroutine r_interpol(x, y, kk, xrt, yt) !--------------------------------------------------------------------------- ! subroutine arguments integer, intent(in) :: kk - double precision, intent(in) :: xrt, x(kk), y(kk) - double precision, intent(out) :: yt + real(prec), intent(in) :: xrt, x(kk), y(kk) + real(prec), intent(out) :: yt ! subroutine local variables integer :: k @@ -2609,11 +2619,11 @@ subroutine normal_crit_y(i, j, q_sk_multi, So, dsc, y_norm, y_crit, area_n, area ! subroutine arguments integer, intent(in) :: i, j - double precision, intent(in) :: q_sk_multi, So, dsc - double precision, intent(out) :: y_norm, y_crit, area_n, area_c + real(prec), intent(in) :: q_sk_multi, So, dsc + real(prec), intent(out) :: y_norm, y_crit, area_n, area_c ! subroutine local variables - double precision :: area_0, width_0, errorY + real(prec) :: area_0, width_0, errorY elevTable = xsec_tab(1, :, i, j) areaTable = xsec_tab(2, :, i, j) @@ -2647,7 +2657,7 @@ subroutine normal_crit_y(i, j, q_sk_multi, So, dsc, y_norm, y_crit, area_n, area end if end subroutine normal_crit_y - double precision function LInterpol(x1, y1, x2, y2, x) + real(prec) function LInterpol(x1, y1, x2, y2, x) implicit none @@ -2657,7 +2667,7 @@ double precision function LInterpol(x1, y1, x2, y2, x) !------------------------------------------------------------------------------------- ! function arguments - double precision, intent(in) :: x1, y1, x2, y2, x + real(prec), intent(in) :: x1, y1, x2, y2, x if (abs(x2-x1).lt.0.0001) then ! to prevent absurdly small value in the denominator @@ -2668,7 +2678,7 @@ double precision function LInterpol(x1, y1, x2, y2, x) end function LInterpol - double precision function intp_y(nrow, xarr, yarr, x) + real(prec) function intp_y(nrow, xarr, yarr, x) implicit none !------------------------------------------------------------------------------------- @@ -2678,11 +2688,11 @@ double precision function intp_y(nrow, xarr, yarr, x) ! function arguments integer, intent(in) :: nrow - double precision, intent(in) :: x - double precision, dimension(nrow), intent(in) :: xarr, yarr + real(prec), intent(in) :: x + real(prec), dimension(nrow), intent(in) :: xarr, yarr ! function local variables integer :: irow - double precision :: x1, y1, x2, y2, y + real(prec) :: x1, y1, x2, y2, y irow = locate(xarr, x) @@ -2718,8 +2728,8 @@ integer function locate(xx, x) !------------------------------------------------------------------------------------ ! function arguments - double precision, intent(in) :: x - double precision, dimension(:), intent(in) :: xx + real(prec), intent(in) :: x + real(prec), dimension(:), intent(in) :: xx ! function local variables integer :: n, jl, jm, ju diff --git a/src/kernel/diffusive/diffusive_lightweight.f90 b/src/kernel/diffusive/diffusive_lightweight.f90 new file mode 100644 index 000000000..4d0afdc7b --- /dev/null +++ b/src/kernel/diffusive/diffusive_lightweight.f90 @@ -0,0 +1,1665 @@ +!------------------------------------------------------------------------------- +module diffusive_lightweight + + use precis + IMPLICIT NONE + +!----------------------------------------------------------------------------- +! Description: +! Numerically solve diffusive wave PDEs using Crank-Nicolson and Hermite +! Interpolation. +! +! Current Code Owner: NOAA-OWP, Inland Hydraulics Team +! +! Code Description: +! Language: Fortran 90. +! This code is written to JULES coding standards v1. +!----------------------------------------------------------------------------- + + ! Module constants + real(prec), parameter :: grav = 9.81 + real(prec), parameter :: TOLERANCE = 1e-8 + + ! Module variables + integer :: nlinks + integer :: mxncomp + integer :: maxTableLength + integer :: nel + integer :: newtonRaphson + integer :: nmstem_rch + integer :: nts_da + integer, dimension(:), allocatable :: mstem_frj + integer, dimension(:), allocatable :: usgs_da_reach + integer, dimension(:,:), allocatable :: frnw_g + real(prec) :: dtini, dxini, cfl, minDx, maxCelerity, theta + real(prec) :: C_llm, D_llm, D_ulm, DD_ulm, DD_llm, q_llm, so_llm + real(prec) :: dt_qtrib + real(prec) :: dmyi, dmyj + real(prec) :: dtini_min + real(prec), dimension(:), allocatable :: eei, ffi, exi, fxi, qcx, diffusivity2, celerity2 + real(prec), dimension(:), allocatable :: elevTable, areaTable, skkkTable + real(prec), dimension(:), allocatable :: pereTable, rediTable + real(prec), dimension(:), allocatable :: convTable, topwTable + real(prec), dimension(:), allocatable :: currentSquareDepth + real(prec), dimension(:), allocatable :: tarr_qtrib, varr_qtrib + real(prec), dimension(:), allocatable :: tarr_da, varr_da + real(prec), dimension(:), allocatable :: area, depth, co, froud, courant + real(prec), dimension(:,:), allocatable :: bo, dx + real(prec), dimension(:,:), allocatable :: areap, qp, z, sk + real(prec), dimension(:,:), allocatable :: celerity, diffusivity, qpx + real(prec), dimension(:,:), allocatable :: pere, oldQ, newQ, oldArea, newArea, oldY, newY + real(prec), dimension(:,:), allocatable :: volRemain + real(prec), dimension(:,:), allocatable :: lateralFlow + real(prec), dimension(:,:), allocatable :: qtrib + real(prec), dimension(:,:), allocatable :: usgs_da + real(prec), dimension(:,:,:,:), allocatable :: xsec_tab + +contains + + subroutine compute_diffusive_couplingtimestep(timestep_ar_g, nts_ql_g, nts_db_g, nts_qtrib_g, nts_da_g, & + mxncomp_g, nrch_g, dx_ar_g, iniq, inidepth, iniqpx, frnw_col, frnw_ar_g, & + qlat_g, dbcd_g, qtrib_g, paradim, para_ar_g, usgs_da_g, usgs_da_reach_g, & + nrow_chxsec_lookuptable, chxsec_lookuptable, z_adj, t_start, t_end, & + q_next_out_time, elv_next_out_time, depth_next_out_time, qpx_next_out_time) + + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Description: + ! Compute diffusive wave routing on the channel networks of either up to version 3.1 + ! of the National Water Model or later versions using the NextGen hydrofabric. This + ! light version skips the computation of hydraulic lookup tables for channel cross sections, but + ! takes already computed ones as input argments instead. Also, instead of computing diffusive + ! wave routing throughout the entire simulation period, this subroutine runs for a coupling + ! time step to interact with other models, including the Muskingum-Cunge model. + ! + ! Method: + ! A Crank-Nicolson solution of the diffusive wave equations is solved with + ! adaptive timestepping to maintain numerical stability. Major operations are + ! as follows: + ! + ! 1. Initialize domain flow and depth. Depth is computed from initial flow + ! 2. For each timestep.... + ! 2.1. Compute domain flow, upstream-to-downstream + ! 2.2. Compute domain depth, downstream-to-upstream + ! 2.3. Determine time step duration needed for stability + ! 2.4. repeat until until final time is reached + ! 3. Record results in output arrays at user-specified time intervals + ! + ! Current Code Owner: NOAA-OWP, Inland Hydraulics Team + ! + ! Development Team: + ! - Dong Ha Kim + ! - Nazmul Azim Beg + ! - Ehab Meselhe + ! - Adam N. Wlostowski + ! - James Halgren + ! - Jacob Hreha + ! + ! Code Description: + ! Language: Fortran 90. + ! This code is written to JULES coding standards v1. + !----------------------------------------------------------------------------- + + ! Subroutine arguments + integer, intent(in) :: mxncomp_g + integer, intent(in) :: nrch_g + integer, intent(in) :: nts_ql_g + integer, intent(in) :: nts_db_g + integer, intent(in) :: nts_qtrib_g + integer, intent(in) :: nts_da_g + integer, intent(in) :: frnw_col + integer, intent(in) :: paradim + integer, intent(in) :: nrow_chxsec_lookuptable + integer, dimension(nrch_g), intent(in) :: usgs_da_reach_g + integer, dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g + real(prec), intent(in) :: t_start + real(prec), intent(in) :: t_end + real(prec), dimension(paradim ), intent(in) :: para_ar_g + real(prec), dimension(:) , intent(in) :: timestep_ar_g(10) + real(prec), dimension(nts_db_g), intent(in) :: dbcd_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: iniq + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: inidepth + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: iniqpx + real(prec), dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g + real(prec), dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g + real(prec), dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in) :: qlat_g + real(prec), dimension(mxncomp_g, nrch_g), intent(in) :: z_adj + real(prec), dimension(11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g), intent(in) :: chxsec_lookuptable + real(prec), dimension(mxncomp_g, nrch_g), intent(out) :: q_next_out_time + real(prec), dimension(mxncomp_g, nrch_g), intent(out) :: elv_next_out_time + real(prec), dimension(mxncomp_g, nrch_g), intent(out) :: depth_next_out_time + real(prec), dimension(mxncomp_g, nrch_g), intent(out) :: qpx_next_out_time + + ! Local variables + integer :: ncomp + integer :: i + integer :: j + integer :: k + integer :: n + integer :: timestep + integer :: kkk + integer :: jm + integer :: usrchj + integer :: downstream_reach_idx + integer :: xcolID + integer :: ycolID + integer :: iel + integer :: dsbc_option + integer :: ts, cwrow + integer :: ri, rj, oi, oj + integer :: nlnk, lnk + integer :: nusrch + integer :: initial_depth_flag + integer :: n_decimal_places + integer, dimension(:), allocatable :: dmy_frj + real(prec) :: x + real(prec) :: saveInterval, width + real(prec) :: maxCourant + real(prec) :: dtini_given, dtini_divisor + real(prec) :: timesDepth + real(prec) :: t + real(prec) :: tfin + real(prec) :: t0 + real(prec) :: q_sk_multi + real(prec) :: maxCelDx + real(prec) :: slope + real(prec) :: y_norm + real(prec) :: temp + real(prec) :: dt_ql + real(prec) :: dt_ub + real(prec) :: dt_db + real(prec) :: dt_da + real(prec) :: wdepth + real(prec) :: q_usrch + real(prec) :: tf0 + real(prec) :: convey + real(prec) :: equiv_one + real(prec) :: mindepth_nstab + real(prec), dimension(:), allocatable :: tarr_ql + real(prec), dimension(:), allocatable :: varr_ql + real(prec), dimension(:), allocatable :: tarr_db + real(prec), dimension(:), allocatable :: varr_db + real(prec), dimension(:), allocatable :: dbcd + real(prec), dimension(:), allocatable :: temp_q + real(prec), dimension(:), allocatable :: temp_elev + + !----------------------------------------------------------------------------- + ! Time domain parameters + dtini = timestep_ar_g(1) ! initial timestep duration [sec] + t0 = timestep_ar_g(2) ! simulation start time [hr] + tfin = timestep_ar_g(3) ! simulation end time [hr] + saveInterval = timestep_ar_g(4) ! output recording interval [sec] + dt_ql = timestep_ar_g(5) ! lateral inflow data time step [sec] + dt_ub = timestep_ar_g(6) ! upstream boundary time step [sec] + dt_db = timestep_ar_g(7) ! downstream boundary time step [sec] + dt_qtrib = timestep_ar_g(8) ! tributary data time step [sec] + dt_da = timestep_ar_g(9) ! data assimilation data time step [sec] + dtini_given = dtini ! preserve user-input timestep duration + dtini_divisor = timestep_ar_g(10) ! numeric value that divide dtini in the next line + dtini_min = dtini/dtini_divisor ! minimum increment in internal time step + + !----------------------------------------------------------------------------- + ! miscellaneous parameters + timesDepth = 4.0 ! water depth multiplier used in readXsection + nel = 501 ! number of sub intervals in look-up tables + nts_da = nts_da_g ! make DA time steps global + + !----------------------------------------------------------------------------- + ! Network mapping (frnw_g) array size parameters + mxncomp = mxncomp_g ! maximum number of nodes in a single reach + nlinks = nrch_g ! number of reaches in the network + + !----------------------------------------------------------------------------- + ! Some essential parameters for Diffusive Wave + cfl = para_ar_g(1) ! maximum Courant number (default: 0.95) + C_llm = para_ar_g(2) ! lower limit of celerity (default: 0.5) + D_llm = para_ar_g(3) ! lower limit of diffusivity (default: 50) + D_ulm = para_ar_g(4) ! upper limit of diffusivity (default: 1000) + DD_llm = para_ar_g(5) ! lower limit of dimensionless diffusivity (default -15) + DD_ulm = para_ar_g(6) ! upper limit of dimensionless diffusivity (default: -10.0) + q_llm = para_ar_g(8) ! lower limit of discharge (default: 0.02831 cms) + so_llm = para_ar_g(9) ! lower limit of channel bed slope (default: 0.0001) + theta = para_ar_g(10) ! weight for computing 2nd derivative: + ! 0: explicit, 1: implicit (default: 1.0) + dsbc_option = para_ar_g(11) ! downstream water depth boundary condition option 1:given water depth data, 2:normal depth + mindepth_nstab = 0.1 ! minimum water depth in meters for numerical stability when estimated/observed downstream boundary depth data is used. + + !----------------------------------------------------------------------------- + ! consider moving variable allocation to a separate module + allocate(area(mxncomp)) + allocate(bo(mxncomp,nlinks)) + allocate(pere(mxncomp,nlinks)) + allocate(areap(mxncomp,nlinks)) + allocate(qp(mxncomp,nlinks)) + allocate(z(mxncomp,nlinks)) + allocate(depth(mxncomp)) + allocate(sk(mxncomp,nlinks)) + allocate(co(mxncomp)) + allocate(dx(mxncomp,nlinks)) + allocate(volRemain(mxncomp-1,nlinks)) + allocate(froud(mxncomp)) + allocate(courant(mxncomp-1)) + allocate(oldQ(mxncomp, nlinks)) + allocate(newQ(mxncomp, nlinks)) + allocate(oldArea(mxncomp, nlinks)) + allocate(newArea(mxncomp, nlinks)) + allocate(oldY(mxncomp, nlinks)) + allocate(newY(mxncomp, nlinks)) + allocate(lateralFlow(mxncomp,nlinks)) + allocate(celerity(mxncomp, nlinks)) + allocate(diffusivity(mxncomp, nlinks)) + allocate(celerity2(mxncomp)) + allocate(diffusivity2(mxncomp)) + allocate(eei(mxncomp)) + allocate(ffi(mxncomp)) + allocate(exi(mxncomp)) + allocate(fxi(mxncomp)) + allocate(qpx(mxncomp, nlinks)) + allocate(qcx(mxncomp)) + allocate(elevTable(nel)) + allocate(areaTable(nel)) + allocate(pereTable(nel)) + allocate(rediTable(nel)) + allocate(convTable(nel)) + allocate(topwTable(nel)) + allocate(skkkTable(nel)) + allocate(xsec_tab(11, nel, mxncomp, nlinks)) + allocate(currentSquareDepth(nel)) + allocate(tarr_ql(nts_ql_g+1), varr_ql(nts_ql_g+1)) + allocate(tarr_qtrib(nts_qtrib_g), varr_qtrib(nts_qtrib_g)) + allocate(tarr_db(nts_db_g), varr_db(nts_db_g)) + allocate(tarr_da(nts_da)) + allocate(dmy_frj(nlinks)) + allocate(frnw_g(nlinks,frnw_col)) + allocate(usgs_da_reach(nlinks)) + allocate(usgs_da(nts_da, nlinks)) + allocate(dbcd(nts_db_g)) + + !-------------------------------------------------------------------------------------------- + frnw_g = frnw_ar_g ! network mapping matrix + usgs_da_reach = usgs_da_reach_g ! contains indices of reaches where usgs data are available + usgs_da = usgs_da_g ! contains usgs data at a related reach + xsec_tab = chxsec_lookuptable ! precomputed hydraulic lookup tables of channel cross sections + z = z_adj ! precomputed channel bottom elevation at each compute node + !----------------------------------------------------------------------------- + ! variable initializations + x = 0.0 + newQ = -999 + newY = -999 + t = t_start ![min] + q_sk_multi = 1.0 + oldQ = iniq + newQ = oldQ + qp = oldQ + dbcd = dbcd_g + q_next_out_time = 0.0 + elv_next_out_time = 0.0 + depth_next_out_time = 0.0 + + !----------------------------------------------------------------------------- + ! Identify mainstem reaches and list their ids in an array + nmstem_rch = 0 + do j = 1, nlinks + nusrch = frnw_g(j,3) ! the number of upstream reaches + if (frnw_g(j, 3+nusrch+1)==555) then ! 555 indicate diffusive reach while -555 already routed reach + nmstem_rch = nmstem_rch + 1 + dmy_frj(nmstem_rch) = j + end if + end do + + ! allocate and populate array for upstream reach ids + allocate (mstem_frj(nmstem_rch)) + do jm = 1, nmstem_rch + mstem_frj(jm) = dmy_frj(jm) + end do + deallocate(dmy_frj) + + !----------------------------------------------------------------------------- + ! create dx array from dx_ar_g and determine minimum dx. + dx = 0. + minDx = 1e10 + + do jm = 1, nmstem_rch !* mainstem reach (=diffusive reach) only + j = mstem_frj(jm) + ncomp = frnw_g(j, 1) + do i = 1, ncomp-1 + dx(i, j) = dx_ar_g(i, j) + end do + minDx = min(minDx, minval(dx(1:ncomp-1, j))) + end do + + !------------------------------------------------------------------------------------------------ + ! Build time arrays for lateral flow, upstream boundary, downstream boundary, and tributary flow + ! time step series for lateral flow. + ! ** ql from t-route starts from the first time step. For example, when t0 and tfin = 0 and 2 hrs with dt=300 sec, + ! ** ql values from t-route starts at 5, 10, 15, ..., 120 min. + do n = 1, nts_ql_g + tarr_ql(n+1) = t0 * 60.0 + dt_ql * & + real(n, KIND(dt_ql)) / 60.0 ! [min] + end do + tarr_ql(1) = t0 * 60 ! initial time step [min] + + ! time step series for tributary flow data + ! ** qtrib from t-route starts from the initial time step . For example, when t0 and tfin = 0 and 2 hrs with dt=300 sec, + ! ** qtrib values from t-route starts at 0, 5, 10, 15, ..., 120 min. + do n=1, nts_qtrib_g + tarr_qtrib(n) = t0 * 60.0 + dt_qtrib * & + real(n-1,KIND(dt_qtrib)) / 60.0 ! [min] + end do + + ! use tailwater downstream boundary observations + ! needed for coastal coupling + ! **** COMING SOON **** + ! time step series for downstream boundary data + do n=1, nts_db_g + tarr_db(n) = t0 * 60.0 + dt_db * & + real(n-1,KIND(dt_db)) / 60.0 ! [min] + end do + + do jm = nmstem_rch, 1, -1 + j = mstem_frj(jm) ! reach index + ncomp = frnw_g(j, 1) ! number of nodes in reach j + if (frnw_g(j, 2) < 0) then + if (dsbc_option == 1) then + ! use tailwater water depth value as downstream boundary condition + do n = 1, nts_db_g + varr_db(n) = dbcd(n) + z(ncomp, j) !* when dbcd is water depth [m], channel bottom elev is added. + end do + endif + endif + enddo + + ! time step series for data assimilation for discharge at bottom node of a related reach + do n=1, nts_da + tarr_da(n) = t0 * 60.0 + dt_da * & + real(n-1,KIND(dt_da)) / 60.0 ! [min] + end do + + !----------------------------------------------------------------------------- + ! Initialize water surface elevation + initial_depth_flag = 0 + !check if initial water depth needs to be computed here or not + do jm = 1, nmstem_rch ! loop over mainstem reaches [upstream-to-downstream] + j = mstem_frj(jm) ! reach index + ncomp = frnw_g(j,1) ! number of nodes in reach j + do i=1, ncomp + if (inidepth(i,j).le.0.0) then + initial_depth_flag = 1 + exit + endif + enddo + enddo + + do j = 1, nlinks + ncomp = frnw_g(j, 1) ! number of nodes in reach j + do i=1, ncomp + oldQ(i,j) = max(oldQ(i,j), q_llm) + qp(i,j) = oldQ(i,j) + enddo + enddo + + if (initial_depth_flag==1) then + ! compute initial water elevation when it is cold start or previously computed values are not provided. + do jm = nmstem_rch, 1, -1 + j = mstem_frj(jm) ! reach index + ncomp = frnw_g(j, 1) ! number of nodes in reach j + + if (frnw_g(j, 2) < 0) then + ! Initial depth at bottom node of tail water reach + if (dsbc_option == 1) then + ! use tailwater downstream boundary observations + oldY(ncomp, j) = intp_y(nts_db_g, tarr_db, varr_db, t) + newY(ncomp, j) = oldY(ncomp, j) + if ((newY(ncomp, j) - z(ncomp, j)).lt.mindepth_nstab) then + newY(ncomp, j) = mindepth_nstab + z(ncomp, j) + end if + else if (dsbc_option == 2) then + ! normal depth as TW boundary condition + xcolID = 10 + ycolID = 1 + oldY(ncomp, j) = intp_xsec_tab(ncomp, j, nel, xcolID, ycolID, oldQ(ncomp,j)) ! normal elevation not depth + newY(ncomp, j) = oldY(ncomp, j) + endif + else + ! Initial depth at bottom node of interior reach is equal to the depth at top node of the downstream reach + downstream_reach_idx = frnw_g(j, 2) + newY(ncomp, j) = newY(1, downstream_reach_idx) + end if + + ! compute newY(i, j) for i=1, ncomp-1 with the given newY(ncomp, j) + ! ** At initial time, oldY(i,j) values at i < ncomp, used in subroutine rtsafe, are not defined. + ! ** So, let's assume the depth values of nodes beyond the bottom node of a reach equal to depth of the bottom node. + oldY(ncomp, j) = newY(ncomp, j) + wdepth = newY(ncomp, j) - z(ncomp, j) + do i = 1, ncomp -1 + oldY(i,j) = wdepth + z(i, j) + end do + + call mesh_diffusive_backward(t, saveInterval, j) + + do i = 1,ncomp + ! copy computed initial depths to initial depth array for first timestep + oldY(i, j) = newY(i, j) + ! Check that node elevation is not lower than bottom node. + ! If node elevation is lower than bottom node elevation, correct. + !if (oldY(i, j) .lt. oldY(ncomp, nlinks)) oldY(i, j) = oldY(ncomp, nlinks) + end do + end do + else + ! Use previous state values when it is warm or hot start. + do jm = 1, nmstem_rch ! loop over mainstem reaches [upstream-to-downstream] + j = mstem_frj(jm) ! reach index + ncomp = frnw_g(j,1) ! number of nodes in reach j + do i=1, ncomp + oldY(i, j) = inidepth(i,j) + z(i,j) + enddo + enddo + endif + + !----------------------------------------------------------------------------- + ! With initial flow and depth, identify the maximum calculated celerity/dx ratio + ! maxCelDx, which is used to determine the duration of the next timestep. Also, + ! compute celerity and diffusivity using initial flow and depth values, which will + ! be used for computing flow at the next time step. + do jm = nmstem_rch, 1, -1 + j = mstem_frj(jm) ! reach index + ncomp = frnw_g(j, 1) ! number of nodes in reach j + allocate(temp_q(ncomp), temp_elev(ncomp)) + do i=1, ncomp + temp_elev(i) = oldY(i,j) + temp_q(i) = oldQ(i,j) + enddo + + call compute_only_celerity_diffusivity(t, ncomp, j, temp_q, temp_elev) + + deallocate(temp_q, temp_elev) + + if (jm == 1) then + maxCelDx = 0. + do i = 1, nmstem_rch + do kkk = 1, frnw_g(mstem_frj(i), 1)-1 + maxCelDx = max(maxCelDx, celerity(kkk, mstem_frj(i)) / dx(kkk, mstem_frj(i))) + end do + end do + endif + end do + + !----------------------------------------------------------------------------- + ! Initializations and re-initializations + qpx = iniqpx + + !----------------------------------------------------------------------------- + ! Ordered network routing computations + do while ( t < t_end) + + !+------------------------------------------------------------------------- + !+ PREDICTOR + !+ + !+ Applies only to mainstem + !+------------------------------------------------------------------------- + do jm = 1, nmstem_rch ! loop over mainstem reaches [upstream-to-downstream] + j = mstem_frj(jm) ! reach index + ncomp = frnw_g(j,1) ! number of nodes in reach j + + ! Calculate the duration of this timestep (dtini) + ! Timestep duration is selected to maintain numerical stability + if (j == mstem_frj(1)) then + call calculateDT(t_start, t, saveInterval, cfl, t_end, maxCelDx, dtini_given) + end if + + ! estimate lateral flow at current time t + do i = 1, ncomp - 1 + do n = 1, nts_ql_g + varr_ql(n+1) = qlat_g(n, i, j) + end do + varr_ql(1) = qlat_g(1, i, j) + lateralFlow(i, j) = intp_y(nts_ql_g+1, tarr_ql, varr_ql, t) + end do + + !+++---------------------------------------------------------------- + !+ Hand over water from upstream to downstream properly according + !+ to network connections, i.e., serial or branching. + !+ Refer to p.52, RM1_MESH + !+++---------------------------------------------------------------- + if (frnw_g(j, 3) > 0) then ! number of upstream reaches, so reach j is not a headwater reach's index + newQ(1, j) = 0.0 + do k = 1, frnw_g(j, 3) ! loop over upstream connected reaches + usrchj = frnw_g(j, 3 + k) + if (any(mstem_frj == usrchj)) then + + ! inflow from upstream mainstem reach's bottom node + q_usrch = newQ(frnw_g(usrchj, 1), usrchj) + else + + ! inflow from upstream tributary reach + do n = 1, nts_qtrib_g + varr_qtrib(n) = qtrib_g(n, usrchj) + end do + tf0 = t !+ dtini / 60. + q_usrch = intp_y(nts_qtrib_g, tarr_qtrib, varr_qtrib, tf0) + end if + + ! add upstream flows to reach head + newQ(1,j)= newQ(1,j) + q_usrch + end do + else + + ! no stream reaches at the upstream of the reach (frnw_g(j,3)==0) such as a headwater reach + newQ(1,j)= 0.0 + end if + + ! Add lateral inflows to the reach head + newQ(1, j) = newQ(1, j) + lateralFlow(1, j) * dx(1, j) + + call mesh_diffusive_forward(t, saveInterval, j) + + end do ! end of j loop for predictor + + !+------------------------------------------------------------------------- + !+ CORRECTOR + !+ + !+ Applies only to mainstem + !+------------------------------------------------------------------------- + do jm = nmstem_rch, 1, -1 ! loop over mainstem reaches [downstream-to-upstream] + j = mstem_frj(jm) + ncomp = frnw_g(j,1) + + !+++------------------------------------------------------------+ + !+ Downstream boundary condition for water elevation either at + !+ a junction or TW. + !+ Refer to p.53-1,RM1_MESH + !+++------------------------------------------------------------+ + if (frnw_g(j, 2) >= 0.0) then + + ! reach of index j has a reach in the downstream + ! set bottom node WSEL equal WSEL in top node of downstream reach + downstream_reach_idx = frnw_g(j,2) + newY(ncomp, j) = newY(1, downstream_reach_idx) + else + + ! reach of index j has NO downstream connection, so it is a tailwater reach + if (dsbc_option == 1) then + ! use downstream boundary data source to set bottom node WSEL value + newY(ncomp, j) = intp_y(nts_db_g, tarr_db, varr_db, t+dtini/60.) + ! to prevent too small water depth at boundary + if ((newY(ncomp, j) - z(ncomp, j)).lt.mindepth_nstab) then + newY(ncomp, j) = mindepth_nstab + z(ncomp, j) + end if + xcolID = 1 + ycolID = 2 + newArea(ncomp, j) = intp_xsec_tab(ncomp, j, nel, xcolID, ycolID, newY(ncomp,j)) ! area of normal elevation + else if (dsbc_option == 2) then + ! Assume normal depth at tailwater downstream boundary + xcolID = 10 + ycolID = 1 + newY(ncomp,j) = intp_xsec_tab(ncomp, j, nel, xcolID, ycolID, abs(newQ(ncomp,j))) ! normal elevation not depth + xcolID = 1 + ycolID = 2 + newArea(ncomp,j) = intp_xsec_tab(ncomp, j, nel, xcolID, ycolID, newY(ncomp,j)) ! area of normal elevation + endif + end if + + ! Calculate WSEL at interior reach nodes + call mesh_diffusive_backward(t, saveInterval, j) + + ! Identify the maximum calculated celerity/dx ratio at this timestep + ! maxCelDx is used to determine the duration of the next timestep + if (jm == 1) then + maxCelDx = 0. + do i = 1, nmstem_rch + do kkk = 1, frnw_g(mstem_frj(i), 1)-1 + maxCelDx = max(maxCelDx, celerity(kkk, mstem_frj(i)) / dx(kkk, mstem_frj(i))) + end do + end do + endif + end do ! end of corrector j loop + + !+------------------------------------------------------------------------- + !+ BOOK KEEPING + !+ + !+------------------------------------------------------------------------- + + ! Calculate Froude and maximum Courant number + do jm = 1, nmstem_rch + j = mstem_frj(jm) + ncomp = frnw_g(j,1) + do i = 1, ncomp + !froud(i) = abs(newQ(i, j)) / sqrt(grav * newArea(i, j) ** 3.0 / bo(i, j)) + if (i < ncomp) then + courant(i) = (newQ(i, j) + newQ(i+1, j)) / (newArea(i, j) + newArea(i+1, j)) * dtini / dx(i, j) + endif + end do + !if (maxCourant < maxval(courant(1:ncomp - 1))) then + ! maxCourant = maxval(courant(1:ncomp-1)) + !end if + end do + + ! Advance model time [minuntes] + t = t + dtini/60. + + ! diffusive wave simulation time print + if (mod(t,5.)==0.) then + print*, "diffusive simulation time in minute=", t + endif + + ! write results to output arrays + if ( (mod((t - t_start ) * 60., saveInterval) <= TOLERANCE) .or. (t == t_end)) then + do jm = 1, nmstem_rch + j = mstem_frj(jm) + ncomp = frnw_g(j, 1) + do i = 1, ncomp + q_next_out_time(i, j) = newQ(i, j) + elv_next_out_time(i, j) = newY(i, j) + depth_next_out_time(i, j) = newY(i, j) - z(i, j) + qpx_next_out_time(i, j) = qpx(i, j) + end do + end do + end if + + ! update of Y, Q and Area vectors + oldY = newY + newY = -999 + oldQ = newQ + newQ = -999 + oldArea = newArea + newArea = -999 + pere = -999 + end do ! end of time loop + + deallocate(frnw_g) + deallocate(area, bo, pere, areap, qp, z, depth, sk, co, dx) + deallocate(volRemain, froud, courant, oldQ, newQ, oldArea, newArea, oldY, newY) + deallocate(lateralFlow, celerity, diffusivity, celerity2, diffusivity2) + deallocate(eei, ffi, exi, fxi, qpx, qcx) + deallocate(elevTable, areaTable, pereTable, rediTable, convTable, topwTable) + deallocate(skkkTable) + deallocate(xsec_tab) + deallocate(currentSquareDepth) + deallocate(tarr_ql, varr_ql, tarr_qtrib, varr_qtrib, tarr_da) + deallocate(mstem_frj) + deallocate(usgs_da_reach, usgs_da) + + end subroutine compute_diffusive_couplingtimestep + + subroutine calculateDT(initialTime, time, saveInterval, maxAllowCourantNo, tfin, max_C_dx, given_dt) + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Description: + ! Calculate time step duration. In order to maintain numerical stability, + ! the model timestep duration must be short enough to keep Courant numbers + ! below a threshold level. Else, numerical instabilities will occur. + ! + ! Method: + ! + ! + ! Current Code Owner: NOAA-OWP, Inland Hydraulics Team + ! + ! Code Description: + ! Language: Fortran 90. + ! This code is written to JULES coding standards v1. + !----------------------------------------------------------------------------- + + ! Subroutine arguments + real(prec), intent(in) :: initialTime ! [minutes] + real(prec), intent(in) :: time ! [minutes] + real(prec), intent(in) :: saveInterval ! [sec] + real(prec), intent(in) :: tfin ! [minutes] + real(prec), intent(in) :: given_dt ! [sec] + real(prec), intent(in) :: maxAllowCourantNo ! = cfl + real(prec), intent(in) :: max_C_dx ! = maxCelDx + + ! Local variables + integer :: a + integer :: b + + !----------------------------------------------------------------------------- + ! Calculate maximum timestep duration for numerical stability + + dtini = maxAllowCourantNo / max_C_dx + + a = floor( (time - initialTime) / ( saveInterval / 60.)) + b = floor(((time - initialTime) + dtini / 60.) / ( saveInterval / 60.)) + if (b > a) then + dtini = (a + 1) * (saveInterval) - (time - initialTime) * 60. + end if + + ! if dtini extends beyond final time, then truncate it + if (time + dtini / 60. > tfin) dtini = (tfin - time) * 60. + + end subroutine + + subroutine mesh_diffusive_forward(t, saveInterval,j) + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Description: + ! Compute discharge using diffusive equation that is numerically solved by + ! Crank-Nicolson + Hermite Interpolation method. + ! + ! Method: + ! + ! Current Code Owner: NOAA-OWP, Inland Hydraulics Team + ! + ! Code Description: + ! Language: Fortran 90. + ! This code is written to JULES coding standards v1. + !----------------------------------------------------------------------------- + + ! Subroutine Arguments + integer, intent(in) :: j + real(prec), intent(in) :: t + real(prec), intent(in) :: saveInterval + + ! Local variables + integer :: ncomp + integer :: i, irow, flag_da, n + real(prec) :: a1, a2, a3, a4 + real(prec) :: b1, b2, b3, b4 + real(prec) :: dd1, dd2, dd3, dd4 + real(prec) :: h1, h2, h3, h4 + real(prec) :: allqlat + real(prec) :: qy, qxy, qxxy, qxxxy + real(prec) :: ppi, qqi, rri, ssi, sxi + real(prec) :: cour, cour2 + real(prec) :: alpha + real(prec) :: currentQ + real(prec) :: eei_ghost, ffi_ghost, exi_ghost + real(prec) :: fxi_ghost, qp_ghost, qpx_ghost + !----------------------------------------------------------------------------- + !* change 20210228: All qlat to a river reach is applied to the u/s boundary + !* Note: lateralFlow(1,j) is already added to the boundary + + eei = -999. + ffi = -999. + exi = -999. + fxi = -999. + eei(1) = 1. + ffi(1) = 0. + exi(1) = 0. + fxi(1) = 0. + + ncomp = frnw_g(j, 1) + + ! sum of lateral inflows along the reach + !allqlat = sum(lateralFlow(2:ncomp - 1, j) * dx(2:ncomp - 1, j)) + allqlat = 0.0 + do i = 2, ncomp - 1 + allqlat = allqlat + lateralFlow(i, j) * dx(i, j) + end do + + !print *, '****** DIFFUSIVE FORWARD *****' + !print *, '---------' + ncomp = frnw_g(j, 1) + + do i = 2, ncomp + + cour = dtini / dx(i - 1, j) + cour2 = abs(celerity(i, j) ) * cour + + a1 = 3.0 * cour2 ** 2.0 - 2.0 * cour2 ** 3.0 + a2 = 1 - a1 + a3 = ( cour2 ** 2.0 - cour2 ** 3.0 ) * dx(i-1,j) + a4 = ( -1.0 * cour2 + 2.0 * cour2 ** 2.0 - cour2 ** 3.0 ) * dx(i-1,j) + + b1 = ( 6.0 * cour2 - 6.0 * cour2 ** 2.0 ) / ( -1.0 * dx(i-1,j) ) + b2 = - b1 + b3 = ( 2.0 * cour2 - 3.0 * cour2 ** 2.0 ) * ( -1.0 ) + b4 = ( -1.0 + 4.0 * cour2 - 3.0 * cour2 ** 2.0 ) * ( -1.0 ) + + dd1 = ( 6.0 - 12.0 * cour2 ) / ( dx(i-1,j) ** 2.0 ) + dd2 = - dd1 + dd3 = ( 2.0 - 6.0 * cour2 ) / dx(i-1,j) + dd4 = ( 4.0 - 6.0 * cour2 ) / dx(i-1,j) + + h1 = 12.0 / ( dx(i-1,j) ** 3.0 ) + h2 = - h1 + h3 = 6.0 / ( dx(i-1,j) ** 2.0 ) + h4 = h3 + + if (i .eq. ncomp) then + alpha = 1.0 + else + alpha = dx(i,j) / dx(i-1,j) + end if + + qy = a1 * oldQ(i-1,j) + a2 * oldQ(i,j) + & + a3 * qpx(i-1,j) + a4 * qpx(i,j) + qxy = b1 * oldQ(i-1,j) + b2 * oldQ(i,j) + & + b3 * qpx(i-1,j) + b4 * qpx(i,j) + qxxy = dd1* oldQ(i-1,j) + dd2* oldQ(i,j) + & + dd3* qpx(i-1,j) + dd4* qpx(i,j) + qxxxy = h1 * oldQ(i-1,j) + h2 * oldQ(i,j) + & + h3 * qpx(i-1,j) + h4 * qpx(i,j) + + ppi = - theta * diffusivity(i,j) * dtini / & + ( dx(i-1,j) ** 2.0 ) * 2.0 / (alpha*(alpha + 1.0)) * alpha + qqi = 1.0 - ppi * (alpha + 1.0) / alpha + rri = ppi / alpha + + ssi = qy + dtini * diffusivity(i,j) * ( 1.0 - theta ) * qxxy + sxi = qxy + dtini * diffusivity(i,j) * ( 1.0 - theta ) * qxxxy + + eei(i) = -1.0 * rri / ( ppi * eei(i-1) + qqi ) + ffi(i) = ( ssi - ppi * ffi(i-1) ) / ( ppi * eei(i-1) + qqi ) + exi(i) = -1.0 * rri / ( ppi * exi(i-1) + qqi ) + fxi(i) = ( sxi - ppi * fxi(i-1) ) / ( ppi * exi(i-1) + qqi ) + end do + + ! Ghost point calculation + cour = dtini / dx(ncomp-1, j) + cour2 = abs(celerity(ncomp-1, j)) * cour + + a1 = 3.0 * cour2 ** 2.0 - 2.0 * cour2 ** 3.0 + a2 = 1 - a1 + a3 = ( cour2 ** 2.0 - cour2 ** 3.0 ) * dx(ncomp-1, j) + a4 = ( -1.0 * cour2 + 2.0 * cour2 ** 2.0 - cour2 ** 3.0 ) * dx(ncomp-1, j) + + b1 = ( 6.0 * cour2 - 6.0 * cour2 ** 2.0 ) / ( -1.0 * dx(ncomp-1, j) ) + b2 = - b1 + b3 = ( 2.0 * cour2 - 3.0 * cour2 ** 2.0 ) * ( -1.0 ) + b4 = ( -1.0 + 4.0 * cour2 - 3.0 * cour2 ** 2.0 ) * ( -1.0 ) + + dd1 = ( 6.0 - 12.0 * cour2 ) / ( dx(ncomp-1, j) ** 2.0 ) + dd2 = - dd1 + dd3 = ( 2.0 - 6.0 * cour2 ) / dx(ncomp-1, j) + dd4 = ( 4.0 - 6.0 * cour2 ) / dx(ncomp-1, j) + + h1 = 12.0 / ( dx(ncomp-1, j) ** 3.0 ) + h2 = - h1 + h3 = 6.0 / ( dx(ncomp-1, j) ** 2.0 ) + h4 = h3 + + alpha = 1.0 + + qy = a1 * oldQ(ncomp, j) + a2 * oldQ(ncomp-1, j) + & + a3 * qpx(ncomp, j) + a4 * qpx(ncomp-1, j) + qxy = b1 * oldQ(ncomp, j) + b2 * oldQ(ncomp-1, j) + & + b3 * qpx(ncomp, j) + b4 * qpx(ncomp-1, j) + qxxy = dd1 * oldQ(ncomp, j) + dd2 * oldQ(ncomp-1, j) + & + dd3 * qpx(ncomp, j) + dd4 * qpx(ncomp-1, j) + qxxxy = h1 * oldQ(ncomp, j) + h2 * oldQ(ncomp-1, j) + & + h3 * qpx(ncomp, j) + h4 * qpx(ncomp-1, j) + + ppi = - theta * diffusivity(ncomp, j) * dtini / & + ( dx(ncomp-1, j) ** 2.0 ) * 2.0 / (alpha*(alpha + 1.0)) * alpha + qqi = 1.0 - ppi * (alpha + 1.0) / alpha + rri = ppi / alpha + + ssi = qy + dtini * diffusivity(ncomp-1, j) * ( 1.0 - theta ) * qxxy + sxi = qxy + dtini * diffusivity(ncomp-1, j) * ( 1.0 - theta ) * qxxxy + + eei_ghost = -1.0 * rri / ( ppi * eei(ncomp) + qqi ) + ffi_ghost = ( ssi - ppi * ffi(ncomp) ) / ( ppi * eei(ncomp) + qqi ) + + exi_ghost = -1.0 * rri / ( ppi * exi(ncomp) + qqi ) + fxi_ghost = ( sxi - ppi * fxi(ncomp) ) / ( ppi * exi(ncomp) + qqi ) + + qp_ghost = oldQ(ncomp-1, j) + qpx_ghost = 0. + + ! when a reach has usgs streamflow data at its location, apply DA + !if (usgs_da_reach(j) /= 0) then + + ! allocate(varr_da(nts_da)) + ! do n = 1, nts_da + ! varr_da(n) = usgs_da(n, j) + ! end do + ! qp(ncomp,j) = intp_y(nts_da, tarr_da, varr_da, t + dtini/60.) + ! flag_da = 1 + ! check usgs_da value is in good quality + ! irow = locate(tarr_da, t + dtini/60.) + ! if (irow == nts_da) then + ! irow = irow-1 + ! endif + ! if ((varr_da(irow)<= -4443.999).or.(varr_da(irow+1)<= -4443.999)) then + ! when usgs data is missing or in poor quality + qp(ncomp,j) = eei(ncomp) * qp_ghost + ffi(ncomp) + flag_da = 0 + ! endif + ! deallocate(varr_da) + ! else + + qp(ncomp,j) = eei(ncomp) * qp_ghost + ffi(ncomp) + flag_da = 0 + ! endif + + qpx(ncomp,j) = exi(ncomp) *qpx_ghost + fxi(ncomp) + + do i = ncomp-1, 1, -1 + qp(i, j) = eei(i) * qp(i+1, j) + ffi(i) + qpx(i, j) = exi(i) * qpx(i+1, j) + fxi(i) + end do + + ! when a reach hasn't been applied to DA + ! if ((usgs_da_reach(j) == 0).or.(flag_da == 0)) then + qp(1, j) = newQ(1, j) + qp(1, j) = qp(1, j) + allqlat + ! endif + + do i = 1, ncomp + if (abs(qp(i, j)) < q_llm) then + qp(i, j) = q_llm + end if + end do + + ! update newQ + do i = 1, ncomp + newQ(i, j) = qp(i, j) + end do + + ! ============== DEBUG to find unstable flow calcls ================= + ! do i= ncomp, 1, -1 + ! if (abs(qp(i,j)) .gt. 2E4) then + ! print *, 'j:', j + ! print *, 'i:', i + ! print *, 'ncomp', ncomp + ! print *, 't:', t + ! print *, 'calculated flow value:', qp(i,j) + ! print *, 'previous upstream flow', oldQ(i-1,j) + ! print *, 'diffusivity(i,j):', diffusivity(i,j) + ! print *, 'eei(i):', eei(i) + ! print *, 'qp(i+1,j):', qp(i+1,j) + ! print *, 'ffi(i):', ffi(i) + ! print *, 'qp_ghost:', qp_ghost + ! print *, 'ffi(ncomp):', ffi(ncomp) + ! print *, 'eei(ncomp):', eei(ncomp) + ! print *, 'cour2:', cour2 + ! print *, 'qp(ncomp,j):', qp(ncomp,j) + ! print *, 'allqlat:', allqlat + ! print *, 'upstream inflow newQ(1,j):', newQ(1,j) + ! stop + ! end if + ! end do + + end subroutine mesh_diffusive_forward + + subroutine mesh_diffusive_backward(t, saveInterval, j) + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Description: + ! Compute water depth using diffusive momentum equation or normal depth + ! with computed Q from the forward step. + ! + ! Method: + ! + ! + ! Current Code Owner: NOAA-OWP, Inland Hydraulics Team + ! + ! Code Description: + ! Language: Fortran 90. + ! This code is written to JULES coding standards v1. + !----------------------------------------------------------------------------- + + ! Subroutine arguments + integer, intent(in) :: j + real(prec), intent(in) :: t + real(prec), intent(in) :: saveInterval + + ! Subroutine variables + integer :: depthCalOk(mxncomp) + integer :: i, ncomp + integer :: jj + real(prec) :: xt + real(prec) :: q_sk_multi, sfi + real(prec) :: vel, currentQ + real(prec) :: S_ncomp + real(prec) :: tempDepthi_1 + real(prec) :: Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds + real(prec) :: C_ulm + + !----------------------------------------------------------------------------- + ncomp = frnw_g(j, 1) + S_ncomp = (-z(ncomp, j) + z(ncomp-1, j)) / dx(ncomp-1, j) + depthCalOk(ncomp) = 1 + elevTable = xsec_tab(1, :,ncomp, j) + areaTable = xsec_tab(2, :,ncomp, j) + topwTable = xsec_tab(6, :,ncomp, j) + + ! Estimate channel area @ downstream boundary + call r_interpol(elevTable, areaTable, nel, & + newY(ncomp, j), newArea(ncomp, j)) + if (newArea(ncomp,j) .eq. -9999) then + print*, 'At j = ',j,'i = ',ncomp, 'time =',t, 'newY(ncomp, j)=', newY(ncomp, j), & + 'newArea(ncomp, j)=', newArea(ncomp,j), 'interpolation of newArea was not possible' +! stop + end if + + ! Estimate channel bottom width @ downstream boundary + call r_interpol(elevTable, topwTable, nel, & + newY(ncomp, j), bo(ncomp, j)) + if (bo(ncomp, j) .eq. -9999) then + print*, 'At j = ',j,', i = ',ncomp, 'time =',t, & + 'interpolation of bo was not possible' +! stop + end if + + do i = ncomp, 1, -1 ! loop through reach nodes [downstream-upstream] + currentQ = qp(i, j) + q_sk_multi=1.0 + + elevTable = xsec_tab(1,:,i,j) + convTable = xsec_tab(5,:,i,j) + areaTable = xsec_tab(2,:,i,j) + pereTable = xsec_tab(3,:,i,j) + topwTable = xsec_tab(6,:,i,j) + skkkTable = xsec_tab(11,:,i,j) + + xt=newY(i, j) + + ! Estimate co(i) (???? what is this?) by interpolation + currentSquareDepth = (elevTable - z(i, j)) ** 2. + call r_interpol(currentSquareDepth, convTable, nel, & + (newY(i, j)-z(i, j)) ** 2.0, co(i)) + + if (co(i) .eq. -9999) then + ! test + print*, t, i, j, newY(i,j), newY(i, j)-z(i, j), co(i) + print*, 'At j = ',j,', i = ',i, 'time =',t, & + 'interpolation of conveyance was not possible, wl', & + newY(i,j), 'z',z(i,j),'previous wl',newY(i+1,j) !, & + !'previous z',z(i+1,j), 'dimensionless_D(i,j)', & + !dimensionless_D(i,j) + ! stop + pause !test + end if + co(i) =q_sk_multi * co(i) + + ! Estimate channel area by interpolation + call r_interpol(elevTable, areaTable, nel, & + xt, newArea(i,j)) + if (newArea(i,j) .eq. -9999) then + print*, 'At j = ',j,', i = ',i, 'time =',t, & + 'interpolation of newArea was not possible' +! stop + end if + + ! Estimate wetted perimeter by interpolation + call r_interpol(elevTable, pereTable, nel, & + xt, pere(i, j)) + + ! Estimate width by interpolation + call r_interpol(elevTable, topwTable, nel, & + xt, bo(i, j)) + + ! Estimate sk(i,j) by interpolation (???? what is sk?) + call r_interpol(elevTable, skkkTable, nel, & + xt, sk(i, j)) + + sfi = qp(i, j) * abs(qp(i, j)) / (co(i) ** 2.0) + + ! ???? What exactly is happening, here ???? + !if (depthCalOk(i) .eq. 1) then + + ! Calculate celerity, diffusivity and velocity + celerity2(i) = 5.0 / 3.0 * abs(sfi) ** 0.3 * & + abs(qp(i, j)) ** 0.4 / bo(i, j) ** 0.4 & + / (1. / (sk(i, j) * q_sk_multi)) ** 0.6 + + ! put upper limit on computed celerity to make sure internal temporal increment not being too small. + if (i.gt.1) then + C_ulm = cfl*dx(i-1,j)/dtini_min + else + C_ulm = cfl*dx(i,j)/dtini_min + endif + + if (celerity2(i).gt.C_ulm) then + celerity2(i) = C_ulm + endif + + diffusivity2(i) = abs(qp(i, j)) / 2.0 / bo(i, j) / abs(sfi) + vel = qp(i, j) / newArea(i, j) + + ! Check celerity value + !if (celerity2(i) .gt. 3.0 * vel) celerity2(i) = vel * 3.0 + !else + ! if (qp(i, j) .lt. 1) then + ! celerity2(i) = C_llm + ! else + ! celerity2(i) = 1.0 + ! end if + ! diffusivity2(i) = diffusivity(i,j) + !end if + + if (i .gt. 1) then + ! ==================================================== + ! Combination of Newton-Raphson and Bisection + ! ==================================================== + Q_cur = qp(i-1, j) + Q_ds = qp(i, j) + z_cur = z(i-1, j) + z_ds = z(i, j) + y_ds = newY(i, j) - z(i, j) + y_ds = max(y_ds, 0.005) + y_cur = rtsafe(i-1, j, Q_cur, Q_ds, z_cur, z_ds, y_ds) + tempDepthi_1 = y_cur + newY(i-1,j) = tempDepthi_1 + z(i-1, j) + + if (newY(i-1, j) > 10.0**5.) newY(i-1, j) = 10.0**5. + + if (newY(i - 1, j) - z(i-1, j) <= 0.) then + print *, ' newY(i-1,j)-z(i-1,j): ', newY(i-1,j)-z(i-1,j) + print *, ' newY(i-1,j): ', newY(i-1,j) + print *, 'z(i-1,j): ', z(i-1,j) + print *, 'qp(i-1,j): ', qp(i-1,j) + print*, 'depth is negative at time=,', t,'j= ', j,'i=',i-1, & + 'newY =', (newY(jj,j),jj=1,ncomp) + !print*, 'dimensionless_D',(dimensionless_D(jj,j),jj=1,ncomp) + print*, 'newQ',(newQ(jj,j),jj=1,ncomp) + print*, 'Bed',(z(jj,j),jj=1,ncomp) + print*, 'dx',(dx(jj,j),jj=1,ncomp-1) + end if + end if + + end do + + celerity(1:ncomp, j) = sum(celerity2(1:ncomp)) / ncomp + + if (celerity(1, j) < C_llm) celerity(1:ncomp,j) = C_llm + + diffusivity(1:ncomp, j) = sum(diffusivity2(1:ncomp)) / ncomp + + do i = 1, ncomp + if (diffusivity(i, j) > D_ulm) diffusivity(i, j) = D_ulm + if (diffusivity(i, j) < D_llm) diffusivity(i, j) = D_llm + end do + + end subroutine mesh_diffusive_backward + + subroutine compute_only_celerity_diffusivity(t, ncomp, j, temp_q, temp_elev) + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Description: + ! Given flow and elevation values at all compute nodes, compute only celerity + ! and diffusivity + ! + ! Method: + ! + ! + ! Current Code Owner: NOAA-OWP, Inland Hydraulics Team + ! + ! Code Description: + ! Language: Fortran 90. + ! This code is written to JULES coding standards v1. + !----------------------------------------------------------------------------- + integer, intent(in) :: j + integer, intent(in) :: ncomp + real(prec), intent(in) :: t + real(prec), dimension(ncomp), intent(in) :: temp_q, temp_elev + + ! Subroutine variables + integer :: i + real(prec) :: q_sk_multi, sfi + real(prec) :: C_ulm + real(prec) :: elev_val, q_val, conv_val, width_val, invmann_val + + q_sk_multi = 1.0 + + do i = 1, ncomp + elevTable = xsec_tab(1,:,i,j) + convTable = xsec_tab(5,:,i,j) + areaTable = xsec_tab(2,:,i,j) + pereTable = xsec_tab(3,:,i,j) + topwTable = xsec_tab(6,:,i,j) + skkkTable = xsec_tab(11,:,i,j) + + q_val = temp_q(i) !qp(i, j) + elev_val = temp_elev(i) !newY(i, j) + + ! Estimat conveyance by interpolation + currentSquareDepth = (elevTable - z(i, j)) ** 2. + + call r_interpol(currentSquareDepth, convTable, nel, & + (elev_val-z(i, j)) ** 2.0, conv_val) + + conv_val =q_sk_multi * conv_val + + ! Estimate width by interpolation + call r_interpol(elevTable, topwTable, nel, & + elev_val, width_val) + + ! Estimate sk(i,j) by interpolation + call r_interpol(elevTable, skkkTable, nel, & + elev_val, invmann_val) + + sfi = q_val* abs(q_val) / (conv_val ** 2.0) + + ! Calculate celerity, diffusivity and velocity + celerity2(i) = 5.0 / 3.0 * abs(sfi) ** 0.3 * & + abs(q_val) ** 0.4 / width_val ** 0.4 & + / (1. / (invmann_val * q_sk_multi)) ** 0.6 + + ! put upper limit on computed celerity to make sure internal temporal increment not being too small. + if (i.gt.1) then + C_ulm = cfl*dx(i-1,j)/dtini_min + else + C_ulm = cfl*dx(i,j)/dtini_min + endif + + if (celerity2(i).gt.C_ulm) then + celerity2(i) = C_ulm + endif + + diffusivity2(i) = abs(q_val) / 2.0 / width_val / abs(sfi) + end do + + celerity(1:ncomp, j) = sum(celerity2(1:ncomp)) / ncomp + + if (celerity(1, j) < C_llm) celerity(1:ncomp,j) = C_llm + + diffusivity(1:ncomp, j) = sum(diffusivity2(1:ncomp)) / ncomp + + do i = 1, ncomp + if (diffusivity(i, j) > D_ulm) diffusivity(i, j) = D_ulm + if (diffusivity(i, j) < D_llm) diffusivity(i, j) = D_llm + end do + + end subroutine compute_only_celerity_diffusivity + + function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds) + + implicit none + + !---------------------------------------------------------------------------------------------------------- + ! Description: + ! Compute water depth using diffusive momentum equation using a combination of Newton-Raphson and Bisection + ! + ! Method: + ! p.1189, Numerical Recipes in F90 + ! Using a combination of newton-raphson and bisection, find the root of a function bracketed + ! between x1 and x2. the root, returned as the function value rtsafe, will be refined until + ! its accuracy is known within ±xacc. + ! - funcd is a user-supplied subroutine that returns both the function value and the first + ! derivative of the function. + ! - parameter: maxit is the maximum allowed number of iterations. + !---------------------------------------------------------------------------------------------------------- + + ! Subroutine arguments + integer, intent(in) :: i, j + real(prec), intent(in) :: Q_cur, Q_ds, z_cur, z_ds, y_ds + + ! Subroutine local variable + integer, parameter :: maxit = 40 + real(prec), parameter :: xacc = 1e-4 + integer :: iter + integer :: xcolID, ycolID + real(prec) :: x1, x2, df, dxx, dxold, f, fh, fl, temp, xh, xl + real(prec) :: y_norm, y_ulm_multi, y_llm_multi, elv_norm, y_old + real(prec) :: rtsafe + + y_ulm_multi = 2.0 + y_llm_multi = 0.1 + + xcolID = 10 + ycolID = 1 + elv_norm = intp_xsec_tab(i, j, nel, xcolID, ycolID, abs(Q_cur)) ! normal elevation not depth + y_norm = elv_norm - z(i, j) + y_old = oldY(i, j) - z(i, j) + + ! option 1 for initial point + !x1 = y_norm * y_llm_multi + !x2 = y_norm * y_ulm_multi + ! option 2 for initial point + x1 = 0.5 * (y_norm + y_old) * y_llm_multi + x2 = 0.5 * (y_norm + y_old) * y_ulm_multi + + call funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, x1, y_ds, fl, df) + + call funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, x2, y_ds, fh, df) + + if ((fl > 0.0 .and. fh > 0.0) .or. (fl < 0.0 .and. fh < 0.0)) then + rtsafe = y_norm + return + endif + + if (fl == 0.0) then + rtsafe = x1 + return + elseif (fh == 0.0) then + rtsafe = x2 + return + elseif (fl < 0.0) then ! orient the search so that f(xl) < 0. + xl = x1 + xh = x2 + else + xh = x1 + xl = x2 + end if + + rtsafe = 0.50 * (x1 + x2) ! initialize the guess for root + dxold = abs(x2 - x1) ! the “stepsize before last,” + dxx = dxold ! and the last step. + + call funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, rtsafe, y_ds, f, df) + + do iter = 1, maxit ! loop over allowed iterations. + if (((rtsafe - xh) * df - f) * ((rtsafe - xl) * df - f) > 0.0 .or. & + abs(2.0 * f) > abs(dxold * df) ) then + ! bisect if newton out of range, or not decreasing fast enough. + dxold = dxx + dxx = 0.50 * (xh - xl) + rtsafe = xl + dxx + if (xl == rtsafe) return + else ! newton step acceptable. take it. + dxold = dxx + dxx = f / df + temp = rtsafe + rtsafe = rtsafe - dxx + if (temp == rtsafe) return + end if + + if (abs(dxx) < xacc) return + + ! one new function evaluation per iteration. + call funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, rtsafe, y_ds, f, df) + + if (f < 0.0) then ! maintain the bracket on the root. + xl = rtsafe + else + xh = rtsafe + end if + end do + + ! when root is not converged: + rtsafe = y_norm + + end function rtsafe + + subroutine funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds, f, df) + + implicit none + + !------------------------------------------------------------------------- + ! Description: + ! Compute diffusive momentum function value and the first derivative + ! + ! Method: + ! Analytical function and its analytical derivative + !------------------------------------------------------------------------- + + ! subroutine arguments + integer, intent(in) :: i, j + real(prec), intent(in) :: Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds + real(prec), intent(out) :: f, df + + ! subroutine local variables + integer :: xcolID, ycolID + real(prec) :: elv_cur, elv_ds, conv_cur, conv_ds, sf_cur, sf_ds, slope + real(prec) :: dKdA_cur, topw_cur + + xcolID = 1 + ! f(y_cur): function of y at the current node + ! - energy slope at downstream node + ycolID = 5 + elv_ds = y_ds + z_ds + conv_ds = intp_xsec_tab(i + 1, j, nel, xcolID, ycolID, elv_ds) + sf_ds = abs(Q_ds) * Q_ds / conv_ds**2.0 + ! - energy slope at current node + elv_cur = y_cur + z_cur + conv_cur = intp_xsec_tab(i, j, nel, xcolID, ycolID, elv_cur) + sf_cur = abs(Q_cur) * Q_cur / conv_cur**2.0 + ! - f(y_cur) + slope = (z(i, j) - z(i+1, j)) / dx(i, j) + slope = max(slope, so_llm) + f = y_cur - y_ds + slope * dx(i, j) - 0.50 * (sf_cur + sf_ds) * dx(i, j) + + ! df/dy at y at current node + ! - dK/dA + ycolID = 9 + dKdA_cur = intp_xsec_tab(i, j, nel, xcolID, ycolID, elv_cur) + ! - top width + ycolID = 6 + topw_cur = intp_xsec_tab(i, j, nel, xcolID, ycolID, elv_cur) + df = 1.0 + (abs(Q_cur) * Q_cur / conv_cur**3.0) * dx(i, j) * topw_cur * dKdA_cur + + end subroutine funcd_diffdepth + + real(prec) function intp_xsec_tab(i, j, nrow, xcolID, ycolID, x) + + implicit none + + !------------------------------------------------------------------- + ! Description: + ! Interpolate given columns of hydraulic lookup table + ! + ! Method: + ! linear interpolation between selected adjacent data points + !------------------------------------------------------------------- + + ! subroutine arguments + integer , intent(in) :: i, j, nrow, xcolID, ycolID + real(prec), intent(in) :: x + + ! subroutine local variables + integer :: irow + real(prec) :: x1, y1, x2, y2, y + real(prec), dimension(nrow) :: xarr, yarr + + xarr = xsec_tab(xcolID, 1:nrow, i, j) + yarr = xsec_tab(ycolID, 1:nrow, i, j) + + irow = locate(xarr, x) + + if (irow == 0) irow = 1 + if (irow == nrow) irow = nrow-1 + + x1 = xarr(irow) + y1 = yarr(irow) + x2 = xarr(irow+1) + y2 = yarr(irow+1) + y = LInterpol(x1, y1, x2, y2, x) + intp_xsec_tab = y + end function intp_xsec_tab + + subroutine r_interpol(x, y, kk, xrt, yt) + + implicit none + !--------------------------------------------------------------------------- + ! Description: + ! Estimate y for a given xrt from x and y arrays using linear interpolation + !--------------------------------------------------------------------------- + ! subroutine arguments + integer, intent(in) :: kk + real(prec), intent(in) :: xrt, x(kk), y(kk) + real(prec), intent(out) :: yt + ! subroutine local variables + integer :: k + + if ((xrt <= maxval(x)) .and. (xrt >= minval(x))) then + + do k = 1, kk-1 + if (((x(k) - xrt) * (x(k+1) - xrt)) <= 0.0) then + yt = (xrt - x(k)) / (x(k+1) - x(k)) * (y(k+1) - y(k)) + y(k) + EXIT + endif + end do + + else if (xrt >= maxval(x)) then +! print*, xrt, ' the given x data point is larger than the upper limit of the set of x data points' +! print*, 'the upper limit: ', maxval(x) + yt = (xrt - x(kk-1)) / (x(kk) - x(kk-1)) * (y(kk) - y(kk-1)) + y(kk-1) ! extrapolation + + else +! print*, xrt, ' the given x data point is less than lower limit of the range of known x data point, ' +! print*, 'so linear interpolation cannot be performed.' +! yt = -9999.0 + yt = minval(y) +! print*, 'The proper range of x is that: ', 'the upper limit: ', maxval(x),& +! ' and lower limit: ', minval(x) +! print*, 'kk', kk +! print*, 't', 'i', dmyi, 'j', dmyj +! print*, 'x', (x(k), k=1, kk) +! print*, 'y', (y(k), k=1, kk) + end if + + end subroutine r_interpol + + subroutine normal_crit_y(i, j, q_sk_multi, So, dsc, y_norm, y_crit, area_n, area_c) + + implicit none + + !------------------------------------------------------------------------------------------------- + ! Description: + ! Estimate normal and critical depth and related areas + ! + ! Method: + ! normal depth by linearly interpolating elevation & conveyance columns w.r.t. given conveyance + ! normal depth area by linearly interpolation elevation & area w.r.t. given normal depth + ! critical depth by iterative method and the computed depth leads to the estimating of the area + !-------------------------------------------------------------------------------------------------- + + ! subroutine arguments + integer, intent(in) :: i, j + real(prec), intent(in) :: q_sk_multi, So, dsc + real(prec), intent(out) :: y_norm, y_crit, area_n, area_c + + ! subroutine local variables + real(prec) :: area_0, width_0, errorY + + elevTable = xsec_tab(1, :, i, j) + areaTable = xsec_tab(2, :, i, j) + pereTable = xsec_tab(3, :, i, j) + convTable = xsec_tab(5, :, i, j) + topwTable = xsec_tab(6, :, i, j) + + call r_interpol(convTable, areaTable, nel, dsc / sqrt(So), area_n) + call r_interpol(convTable, elevTable, nel, dsc / sqrt(So), y_norm) + call r_interpol(elevTable, areaTable, nel, oldY(i, j), area_0) ! initial estimate for critical depth + call r_interpol(elevTable, topwTable, nel, oldY(i, j), width_0) ! initial estimate for critical depth + + area_c = area_0 + errorY = 100. + + do while (errorY > 0.0001) + area_c = (dsc * dsc * width_0 / grav) ** (1./3.) + errorY = abs(area_c - area_0) + + call r_interpol(areaTable, topwTable, nel, area_c, width_0) + + area_0 = area_c + end do + + call r_interpol(areaTable,elevTable,nel,area_c,y_crit) + + if (y_norm .eq. -9999) then + print*, 'At j = ',j,', i = ',i, 'interpolation of y_norm in calculating normal area was not possible, Q', & + dsc,'slope',So +! stop + end if + end subroutine normal_crit_y + + real(prec) function LInterpol(x1, y1, x2, y2, x) + + implicit none + + !------------------------------------------------------------------------------------- + ! Description: + ! Estimate y for a given x applying linear interpolation to two given (x, y) points + !------------------------------------------------------------------------------------- + + ! function arguments + real(prec), intent(in) :: x1, y1, x2, y2, x + + if (abs(x2-x1).lt.0.0001) then + ! to prevent absurdly small value in the denominator + LInterpol = 0.5*(y1 + y2) + else + LInterpol = (y2 - y1) / (x2 - x1) * (x - x1) + y1 + endif + + end function LInterpol + + real(prec) function intp_y(nrow, xarr, yarr, x) + + implicit none + !------------------------------------------------------------------------------------- + ! Description: + ! Estimate y for a given x applying linear interpolation to two x and y arrays + !------------------------------------------------------------------------------------- + + ! function arguments + integer, intent(in) :: nrow + real(prec), intent(in) :: x + real(prec), dimension(nrow), intent(in) :: xarr, yarr + ! function local variables + integer :: irow + real(prec) :: x1, y1, x2, y2, y + + irow = locate(xarr, x) + + if (irow == 0) irow = 1 + if (irow == nrow) irow = nrow - 1 + + x1 = xarr(irow) + y1 = yarr(irow) + x2 = xarr(irow+1) + y2 = yarr(irow+1) + y = LInterpol(x1, y1, x2, y2, x) + intp_y = y + + end function intp_y + + integer function locate(xx, x) + + implicit none + + !------------------------------------------------------------------------------------ + ! Description: + ! Locate function in p.1045,Numerical Recipes in Fortran f90 + ! + ! Method: + ! klo=max(min(locate(xa,x),n-1),1) In the Fortran 77 version of splint, + ! there is in-line code to find the location in the table by bisection. Here + ! we prefer an explicit call to locate, which performs the bisection. On + ! some massively multiprocessor (MMP) machines, one might substitute a different, + ! more parallel algorithm (see next note). + ! Given an array xx(1:N), and given a value x, returns a value j such that x is between + ! xx(j) and xx(j + 1). xx must be monotonic, either increasing or decreasing. + ! j = 0 or j = N is returned to indicate that x is out of range. + !------------------------------------------------------------------------------------ + + ! function arguments + real(prec), intent(in) :: x + real(prec), dimension(:), intent(in) :: xx + + ! function local variables + integer :: n, jl, jm, ju + logical :: ascnd + + n = size(xx) + ascnd = (xx(n) >= xx(1)) ! True if ascending order of table, false otherwise. + jl = 0 ! Initialize lower + ju = n + 1 ! and upper limits. + + do + if ((ju - jl) <= 1) exit ! Repeat until this condition is satisfied. + + jm = (ju + jl) / 2 ! Compute a midpoint, + + if (ascnd .eqv. (x >= xx(jm))) then + jl = jm ! and replace either the lower limit + else + ju = jm ! or the upper limit, as appropriate. + end if + end do + + if (x == xx(1)) then ! Then set the output, being careful with the endpoints. + locate = 1 + else if (x == xx(n)) then + locate = n - 1 + else + locate = jl + end if + + end function locate + +end module diffusive_lightweight diff --git a/src/kernel/diffusive/makefile b/src/kernel/diffusive/makefile index 0844047b9..bf366ada4 100644 --- a/src/kernel/diffusive/makefile +++ b/src/kernel/diffusive/makefile @@ -12,18 +12,26 @@ else FCFLAGS = $(F90FLAGSGFORTRAN) endif -diffusive.o: diffusive.f90 +diffusive.o: diffusive.f90 precis.mod $(FC) $(FCFLAGS) -o $@ $< pydiffusive.o: pydiffusive.f90 $(FC) $(FCFLAGS) -o $@ $< -chxsec_lookuptable.o: chxsec_lookuptable.f90 +chxsec_lookuptable.o: chxsec_lookuptable.f90 precis.mod $(FC) $(FCFLAGS) -o $@ $< pychxsec_lookuptable.o: pychxsec_lookuptable.f90 $(FC) $(FCFLAGS) -o $@ $< +diffusive_lightweight.o: diffusive_lightweight.f90 precis.mod + $(FC) $(FCFLAGS) -o $@ $< + +pydiffusive_lightweight.o: pydiffusive_lightweight.f90 + $(FC) $(FCFLAGS) -o $@ $< + +precis.mod: varPrecision.f90 + $(FC) $(FCFLAGS) $< install: cp *.o ../../../src/troute-routing/troute/routing/fast_reach @@ -37,4 +45,9 @@ clean: rm -f ../../../src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.o rm -f ../../../src/troute-routing/troute/routing/fast_reach/pychxsec_lookuptable.o rm -f ../../../src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.c - rm -f ../../../src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.cpython*.so \ No newline at end of file + rm -f ../../../src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.cpython*.so + rm -f ../../../src/troute-routing/troute/routing/fast_reach/diffusive_lightweight.o + rm -f ../../../src/troute-routing/troute/routing/fast_reach/pydiffusive_lightweight.o + rm -f ../../../src/troute-routing/troute/routing/fast_reach/diffusive_lightweight.c + rm -f ../../../src/troute-routing/troute/routing/fast_reach/diffusive_lightweight.cpython*.so + diff --git a/src/kernel/diffusive/pychxsec_lookuptable.f90 b/src/kernel/diffusive/pychxsec_lookuptable.f90 index 09dcd002b..925d8d1bf 100644 --- a/src/kernel/diffusive/pychxsec_lookuptable.f90 +++ b/src/kernel/diffusive/pychxsec_lookuptable.f90 @@ -1,6 +1,6 @@ module lookuptable_interface -use, intrinsic :: iso_c_binding +use, intrinsic :: iso_c_binding, only: c_float, c_float, c_int use lookuptable, only: chxsec_lookuptable_calc implicit none @@ -16,22 +16,22 @@ subroutine c_chxsec_lookuptable_calc(mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_a integer(c_int), intent(in) :: frnw_col integer(c_int), intent(in) :: mxnbathy_g integer(c_int), intent(in) :: nrow_chxsec_lookuptable - real(c_double), intent(in) :: so_lowerlimit_g + real(c_float), intent(in) :: so_lowerlimit_g integer(c_int), dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g integer(c_int), dimension(mxncomp_g, nrch_g), intent(in) :: size_bathy_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: bo_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: traps_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: twcc_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: manncc_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g - real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: x_bathy_g - real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: z_bathy_g - real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: mann_bathy_g - real(c_double), dimension(11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g), intent(out) :: chxsec_lookuptable - real(c_double), dimension(mxncomp_g, nrch_g), intent(out) :: z_adj + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: bo_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: traps_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: twcc_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: manncc_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g + real(c_float), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: x_bathy_g + real(c_float), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: z_bathy_g + real(c_float), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in) :: mann_bathy_g + real(c_float), dimension(11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g), intent(out) :: chxsec_lookuptable + real(c_float), dimension(mxncomp_g, nrch_g), intent(out) :: z_adj call chxsec_lookuptable_calc(mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g, & mann_ar_g, manncc_ar_g, dx_ar_g, so_lowerlimit_g, frnw_col, frnw_ar_g, & diff --git a/src/kernel/diffusive/pydiffusive.f90 b/src/kernel/diffusive/pydiffusive.f90 index fc854bc42..9bcb8873f 100644 --- a/src/kernel/diffusive/pydiffusive.f90 +++ b/src/kernel/diffusive/pydiffusive.f90 @@ -1,6 +1,6 @@ module diffusive_interface -use, intrinsic :: iso_c_binding +use, intrinsic :: iso_c_binding, only: c_float, c_float, c_int use diffusive, only: diffnw implicit none @@ -24,25 +24,25 @@ subroutine c_diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_ integer(c_int), dimension(nrch_g), intent(in) :: usgs_da_reach_g integer(c_int), dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g integer(c_int), dimension(mxncomp_g, nrch_g), intent(in) :: size_bathy_g - real(c_double), dimension(nts_db_g), intent(in) :: dbcd_g - real(c_double), dimension(paradim), intent(in) :: para_ar_g - real(c_double), dimension(:), intent(in) :: timestep_ar_g(10) - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g, twcc_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g, manncc_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(inout) :: so_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g, rdx_ar_g - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: iniq - real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: z_thalweg_g - real(c_double), dimension(nts_ub_g, nrch_g), intent(in) :: ubcd_g - real(c_double), dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g - real(c_double), dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g - real(c_double), dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in) :: qlat_g - real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: x_bathy_g - real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: z_bathy_g - real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: mann_bathy_g - real(c_double), dimension(cwnrow_g, cwncol_g), intent(in ) :: crosswalk_g - real(c_double), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g, elv_ev_g, depth_ev_g + real(c_float), dimension(nts_db_g), intent(in) :: dbcd_g + real(c_float), dimension(paradim), intent(in) :: para_ar_g + real(c_float), dimension(:), intent(in) :: timestep_ar_g(10) + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g, twcc_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g, manncc_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(inout) :: so_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g, rdx_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: iniq + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: z_thalweg_g + real(c_float), dimension(nts_ub_g, nrch_g), intent(in) :: ubcd_g + real(c_float), dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g + real(c_float), dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g + real(c_float), dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in) :: qlat_g + real(c_float), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: x_bathy_g + real(c_float), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: z_bathy_g + real(c_float), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: mann_bathy_g + real(c_float), dimension(cwnrow_g, cwncol_g), intent(in ) :: crosswalk_g + real(c_float), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g, elv_ev_g, depth_ev_g call diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qtrib_g, nts_da_g, & mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g, mann_ar_g, & diff --git a/src/kernel/diffusive/pydiffusive_lightweight.f90 b/src/kernel/diffusive/pydiffusive_lightweight.f90 new file mode 100644 index 000000000..a138b5238 --- /dev/null +++ b/src/kernel/diffusive/pydiffusive_lightweight.f90 @@ -0,0 +1,51 @@ +module diffusive_lightweight_interface + +use, intrinsic :: iso_c_binding, only: c_double, c_float, c_int +use diffusive_lightweight, only: compute_diffusive_couplingtimestep + +implicit none +contains +subroutine c_compute_diffusive_couplingtimestep(timestep_ar_g, nts_ql_g, nts_db_g, nts_qtrib_g, nts_da_g, & + mxncomp_g, nrch_g, dx_ar_g, iniq, inidepth, iniqpx, frnw_col, frnw_ar_g, & + qlat_g, dbcd_g, qtrib_g, paradim, para_ar_g, usgs_da_g, usgs_da_reach_g, & + nrow_chxsec_lookuptable, chxsec_lookuptable, z_adj, t_start, t_end, & + q_next_out_time, elv_next_out_time, depth_next_out_time, qpx_next_out_time) bind(c) + + integer(c_int), intent(in) :: mxncomp_g + integer(c_int), intent(in) :: nrch_g + integer(c_int), intent(in) :: nts_ql_g + integer(c_int), intent(in) :: nts_db_g + integer(c_int), intent(in) :: nts_qtrib_g + integer(c_int), intent(in) :: nts_da_g + integer(c_int), intent(in) :: frnw_col + integer(c_int), intent(in) :: paradim + integer(c_int), intent(in) :: nrow_chxsec_lookuptable + integer(c_int), dimension(nrch_g), intent(in) :: usgs_da_reach_g + integer(c_int), dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g + real(c_float), intent(in) :: t_start + real(c_float), intent(in) :: t_end + real(c_float), dimension(paradim ), intent(in) :: para_ar_g + real(c_float), dimension(:) , intent(in) :: timestep_ar_g(10) + real(c_float), dimension(nts_db_g), intent(in) :: dbcd_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: iniq + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: inidepth + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: iniqpx + real(c_float), dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g + real(c_float), dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g + real(c_float), dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in) :: qlat_g + real(c_float), dimension(mxncomp_g, nrch_g), intent(in) :: z_adj + real(c_float), dimension(11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g), intent(in) :: chxsec_lookuptable + real(c_float), dimension(mxncomp_g, nrch_g), intent(out) :: q_next_out_time + real(c_float), dimension(mxncomp_g, nrch_g), intent(out) :: elv_next_out_time + real(c_float), dimension(mxncomp_g, nrch_g), intent(out) :: depth_next_out_time + real(c_float), dimension(mxncomp_g, nrch_g), intent(out) :: qpx_next_out_time + + call compute_diffusive_couplingtimestep(timestep_ar_g, nts_ql_g, nts_db_g, nts_qtrib_g, nts_da_g, & + mxncomp_g, nrch_g, dx_ar_g, iniq, inidepth, iniqpx, frnw_col, frnw_ar_g, & + qlat_g, dbcd_g, qtrib_g, paradim, para_ar_g, usgs_da_g, usgs_da_reach_g, & + nrow_chxsec_lookuptable, chxsec_lookuptable, z_adj, t_start, t_end, & + q_next_out_time, elv_next_out_time, depth_next_out_time, qpx_next_out_time) + +end subroutine c_compute_diffusive_couplingtimestep +end module diffusive_lightweight_interface \ No newline at end of file diff --git a/src/kernel/diffusive/varPrecision.f90 b/src/kernel/diffusive/varPrecision.f90 new file mode 100644 index 000000000..0220b33d9 --- /dev/null +++ b/src/kernel/diffusive/varPrecision.f90 @@ -0,0 +1,15 @@ +module precis + implicit none + + !!NOTE: kind(1.0) and selected_real_kind(4) produce equivalent results + integer, parameter :: prec = kind(1.0) + !integer, parameter :: prec = selected_real_kind(4) + + !!NOTE: kind(1.d0) and selected_real_kind(8) produce equivalent results + integer, parameter :: prec_db = kind(1.0d0) + !integer, parameter :: prec = selected_real_kind(8) + + + !integer, parameter :: prec = selected_real_kind(15) + +end module precis \ No newline at end of file diff --git a/src/troute-config/troute/config/compute_parameters.py b/src/troute-config/troute/config/compute_parameters.py index 26b21aa76..04caf2e57 100644 --- a/src/troute-config/troute/config/compute_parameters.py +++ b/src/troute-config/troute/config/compute_parameters.py @@ -17,9 +17,10 @@ "by-subnetwork-jit-clustered", "by-subnetwork-diffusive", "bmi", + "serial-hybrid-routing" ] -ComputeKernel = Literal["V02-structured", "diffusive", "diffusice_cnt"] +ComputeKernel = Literal["V02-structured","V02-structured-hybrid-routing", "diffusive"] class ComputeParameters(BaseModel): diff --git a/src/troute-network/troute/AbstractNetwork.py b/src/troute-network/troute/AbstractNetwork.py index 344bd9919..b5c6622bc 100644 --- a/src/troute-network/troute/AbstractNetwork.py +++ b/src/troute-network/troute/AbstractNetwork.py @@ -327,11 +327,14 @@ def segment_index(self): """ # list of all segments in the domain (MC + diffusive) self._segment_index = self.dataframe.index - if self._routing.diffusive_network_data: - for tw in self._routing.diffusive_network_data: - self._segment_index = self._segment_index.append( - pd.Index(self._routing.diffusive_network_data[tw]['mainstem_segs']) - ) + # Skip this process for enabling the exchange of flow between MC and diffusive during runtime. + if 'hybrid-routing' not in self.compute_parameters['compute_kernel']: + if self._routing.diffusive_network_data: + for tw in self._routing.diffusive_network_data: + self._segment_index = self._segment_index.append( + pd.Index(self._routing.diffusive_network_data[tw]['mainstem_segs']) + ) + return self._segment_index @property @@ -490,7 +493,7 @@ def initialize_routing_scheme(self,): if value==routing_type: routing_scheme = key - routing = routing_scheme(self.hybrid_parameters) + routing = routing_scheme(self.hybrid_parameters, self.compute_parameters) ( self._dataframe, @@ -662,7 +665,6 @@ def initial_warmstate_preprocess(self, from_files, value_dict): #---------------------------------------------------------------------------- start_time = time.time() LOG.info("setting channel initial states ...") - # if lite restart file is provided, the read channel initial states from it if from_files: if restart_parameters.get("lite_channel_restart_file", None): diff --git a/src/troute-network/troute/AbstractRouting.py b/src/troute-network/troute/AbstractRouting.py index 1c46dfcf2..b71725b02 100644 --- a/src/troute-network/troute/AbstractRouting.py +++ b/src/troute-network/troute/AbstractRouting.py @@ -86,7 +86,7 @@ class AbstractRouting(ABC): """ """ - __slots__ = ["hybrid_params", "_diffusive_domain", "_coastal_boundary_depth_df", + __slots__ = ["hybrid_params", "compute_params", "_diffusive_domain", "_coastal_boundary_depth_df", "_diffusive_network_data", "_topobathy_df", "_refactored_diffusive_domain", "_refactored_diffusive_network_data", "_refactored_reaches", "_unrefactored_topobathy_df", "_all_links", "_bad_topobathy_links", "_dataframe"] @@ -154,8 +154,9 @@ def dataframe(self): class MCOnly(AbstractRouting): - def __init__(self, _): + def __init__(self, _, __): self.hybrid_params = None + self.compute_params = None super().__init__() @@ -201,8 +202,9 @@ def dataframe(self): class MCwithDiffusive(AbstractRouting): - def __init__(self, hybrid_params): + def __init__(self, hybrid_params, compute_params): self.hybrid_params = hybrid_params + self.compute_params = compute_params super().__init__() @@ -312,17 +314,19 @@ def update_routing_domain(self, dataframe, connections, waterbody_dataframe): # ==== remove diffusive domain segs from MC domain ==== # drop indices from param_df. Make sure when mainstem_segs accidently includes lake ids, exclude them - # from id list to be dropped from dataframe as dataframe only handles channel parameters. - existing_indicies_in_dataframe = [id for id in mainstem_segs if id in dataframe.index] - dataframe = dataframe.drop(existing_indicies_in_dataframe) - - # remove keys from connections dictionary - for s in mainstem_segs: - connections.pop(s) + # from id list to be dropped from dataframe as dataframe only handles channel parameters. + # Skip this process for enabling the exchange of flow between MC and diffusive during runtime. + if 'hybrid-routing' not in self.compute_params['compute_kernel']: + existing_indicies_in_dataframe = [id for id in mainstem_segs if id in dataframe.index] + dataframe = dataframe.drop(existing_indicies_in_dataframe) + + # remove keys from connections dictionary + for s in mainstem_segs: + connections.pop(s) - # update downstream connections of trib segs - for us in trib_segs: - connections[us] = [] + # update downstream connections of trib segs + for us in trib_segs: + connections[us] = [] return dataframe, connections @@ -381,9 +385,9 @@ def diffusive_domain_by_both_ends_streamid(self, connections, headlink_mainstem, class MCwithDiffusiveNatlXSectionNonRefactored(MCwithDiffusive): - def __init__(self, hybrid_params): + def __init__(self, hybrid_params, compute_params): - super().__init__(hybrid_params = hybrid_params) + super().__init__(hybrid_params = hybrid_params, compute_params=compute_params) @property def topobathy_df(self): @@ -431,9 +435,9 @@ def topobathy_df(self): class MCwithDiffusiveNatlXSectionRefactored(MCwithDiffusive): - def __init__(self, hybrid_params): + def __init__(self, hybrid_params, compute_params): - super().__init__(hybrid_params = hybrid_params) + super().__init__(hybrid_params = hybrid_params, compute_params=compute_params) @property def topobathy_df(self): diff --git a/src/troute-nwm/src/nwm_routing/__main__.py b/src/troute-nwm/src/nwm_routing/__main__.py index 4b2150240..a26c739d0 100644 --- a/src/troute-nwm/src/nwm_routing/__main__.py +++ b/src/troute-nwm/src/nwm_routing/__main__.py @@ -77,8 +77,7 @@ def main_v04(argv): # perform initial warmstate preprocess. network_start_time = time.time() - - #if "ngen_nexus_file" in supernetwork_parameters: + if supernetwork_parameters["network_type"] == 'HYFeaturesNetwork': network = HYFeaturesNetwork(supernetwork_parameters, waterbody_parameters, @@ -106,7 +105,6 @@ def main_v04(argv): ) duplicate_ids_df = pd.DataFrame() - network_end_time = time.time() task_times['network_creation_time'] = network_end_time - network_start_time @@ -1233,6 +1231,9 @@ def nwm_route( waterbody_types_df, waterbody_type_specified, subnetwork_list, + diffusive_network_data, + topobathy_df, + coastal_boundary_depth_df, flowveldepth_interorder, from_files = from_files, ) @@ -1242,71 +1243,72 @@ def nwm_route( results = results[0] # run diffusive side of a hybrid simulation - if diffusive_network_data: - start_time_diff = time.time() - ''' - # retrieve MC-computed streamflow value at upstream boundary of diffusive mainstem - qvd_columns = pd.MultiIndex.from_product( - [range(nts), ["q", "v", "d"]] - ).to_flat_index() - flowveldepth = pd.concat( - [pd.DataFrame(r[1], index=r[0], columns=qvd_columns) for r in results], - copy=False, - ) - ''' - #upstream_boundary_flow={} - #for tw,v in diffusive_network_data.items(): - # upstream_boundary_link = diffusive_network_data[tw]['upstream_boundary_link'] - # flow_ = flowveldepth.loc[upstream_boundary_link][0::3] - # the very first value at time (0,q) is flow value at the first time step after initial time. - # upstream_boundary_flow[tw] = flow_ + if 'hybrid-routing' not in compute_kernel: + if diffusive_network_data: + start_time_diff = time.time() + ''' + # retrieve MC-computed streamflow value at upstream boundary of diffusive mainstem + qvd_columns = pd.MultiIndex.from_product( + [range(nts), ["q", "v", "d"]] + ).to_flat_index() + flowveldepth = pd.concat( + [pd.DataFrame(r[1], index=r[0], columns=qvd_columns) for r in results], + copy=False, + ) + ''' + #upstream_boundary_flow={} + #for tw,v in diffusive_network_data.items(): + # upstream_boundary_link = diffusive_network_data[tw]['upstream_boundary_link'] + # flow_ = flowveldepth.loc[upstream_boundary_link][0::3] + # the very first value at time (0,q) is flow value at the first time step after initial time. + # upstream_boundary_flow[tw] = flow_ - if (firstRun): - compute_log_diff( - logFileName, - diffusive_network_data, - topobathy_df, - refactored_diffusive_domain, - refactored_reaches, - coastal_boundary_depth_df, - unrefactored_topobathy_df, - ) - - # call diffusive wave simulation and append results to MC results - results.extend( - compute_diffusive_routing( - results, - diffusive_network_data, - cpu_pool, - t0, - dt, - nts, - q0, - qlats, - qts_subdivisions, - usgs_df, - lastobs_df, - da_parameter_dict, - waterbodies_df, - topobathy_df, - refactored_diffusive_domain, - refactored_reaches, - coastal_boundary_depth_df, - unrefactored_topobathy_df, + if (firstRun): + compute_log_diff( + logFileName, + diffusive_network_data, + topobathy_df, + refactored_diffusive_domain, + refactored_reaches, + coastal_boundary_depth_df, + unrefactored_topobathy_df, + ) + + # call diffusive wave simulation and append results to MC results + results.extend( + compute_diffusive_routing( + results, + diffusive_network_data, + cpu_pool, + t0, + dt, + nts, + q0, + qlats, + qts_subdivisions, + usgs_df, + lastobs_df, + da_parameter_dict, + waterbodies_df, + topobathy_df, + refactored_diffusive_domain, + refactored_reaches, + coastal_boundary_depth_df, + unrefactored_topobathy_df, + ) ) - ) - LOG.debug("Diffusive computation complete in %s seconds." % (time.time() - start_time_diff)) + LOG.debug("Diffusive computation complete in %s seconds." % (time.time() - start_time_diff)) - else: + else: - if (firstRun): - with open(logFileName, 'a') as preRunLog: - preRunLog.write("**********************\n") - preRunLog.write("No diffusive routing. \n") - preRunLog.write("**********************\n") - preRunLog.close() + if (firstRun): + with open(logFileName, 'a') as preRunLog: + preRunLog.write("**********************\n") + preRunLog.write("No diffusive routing. \n") + preRunLog.write("**********************\n") + preRunLog.close() - LOG.debug("ordered reach computation complete in %s seconds." % (time.time() - start_time)) + LOG.debug("ordered reach computation complete in %s seconds." % (time.time() - start_time)) return results, subnetwork_list diff --git a/src/troute-routing/setup.py b/src/troute-routing/setup.py index a14939f5a..1e4df4cf2 100644 --- a/src/troute-routing/setup.py +++ b/src/troute-routing/setup.py @@ -112,6 +112,16 @@ def build_extensions(self): libraries=[], ) +hybrid_routing_reach = Extension( + "troute.routing.fast_reach.hybrid_routing_reach", + sources=["troute/routing/fast_reach/hybrid_routing_reach.{}".format(ext)], + include_dirs=_include_paths, + libraries=[], + library_dirs=[], + extra_objects=[], + extra_compile_args=["-O2", "-g"], +) + chxsec_lookuptable = Extension( "troute.routing.fast_reach.chxsec_lookuptable", sources=["troute/routing/fast_reach/chxsec_lookuptable.{}".format(ext)], @@ -124,8 +134,20 @@ def build_extensions(self): libraries=[], ) +diffusive_lightweight = Extension( + "troute.routing.fast_reach.diffusive_lightweight", + sources=["troute/routing/fast_reach/diffusive_lightweight.{}".format(ext)], + extra_objects=[ + "troute/routing/fast_reach/diffusive_lightweight.o", + "troute/routing/fast_reach/pydiffusive_lightweight.o", + ], + include_dirs=[np.get_include()], + extra_compile_args=["-O2", "-g"], + libraries=[], +) + package_data = {"troute.fast_reach": ["reach.pxd", "fortran_wrappers.pxd", "utils.pxd"]} -ext_modules = [reach, mc_reach, diffusive, simple_da, chxsec_lookuptable] +ext_modules = [reach, mc_reach, diffusive, simple_da, hybrid_routing_reach, chxsec_lookuptable, diffusive_lightweight] if USE_CYTHON: from Cython.Build import cythonize diff --git a/src/troute-routing/troute/routing/compute.py b/src/troute-routing/troute/routing/compute.py index e96aedf3b..4d8c8eba4 100644 --- a/src/troute-routing/troute/routing/compute.py +++ b/src/troute-routing/troute/routing/compute.py @@ -11,8 +11,11 @@ import troute.nhd_network as nhd_network from troute.routing.fast_reach.mc_reach import compute_network_structured +from troute.routing.fast_reach.hybrid_routing_reach import compute_network_structured_with_hybrid_routing + import troute.routing.diffusive_utils_v02 as diff_utils from troute.routing.fast_reach import diffusive +from troute.routing.fast_reach import chxsec_lookuptable import logging @@ -22,6 +25,7 @@ compute_network_structured, { "V02-structured": compute_network_structured, + "V02-structured-hybrid-routing": compute_network_structured_with_hybrid_routing, }, ) @@ -508,10 +512,13 @@ def compute_nhd_routing_v02( data_assimilation_parameters, waterbody_types_df, waterbody_type_specified, - subnetwork_list, + subnetwork_list, + diffusive_network_data, + topobathy, + coastal_boundary_depth_df, flowveldepth_interorder = {}, from_files = True, -): + ): da_decay_coefficient = da_parameter_dict.get("da_decay_coefficient", 0) param_df["dt"] = dt @@ -519,6 +526,7 @@ def compute_nhd_routing_v02( start_time = time.time() compute_func = _compute_func_map[compute_func_name] + if parallel_compute_method == "by-subnetwork-jit-clustered": # Create subnetwork objects if they have not already been created @@ -757,7 +765,7 @@ def compute_nhd_routing_v02( from_files, offnetwork_upstreams ) - + # results_subn[order].append( # compute_func( jobs.append( @@ -1246,7 +1254,7 @@ def compute_nhd_routing_v02( waterbody_types_df_sub, t0, from_files, - set() + offnetwork_upstreams=set(), # TODO: need to be defined. ) jobs.append( @@ -1414,7 +1422,7 @@ def compute_nhd_routing_v02( t0, from_files, ) - + results.append( compute_func( nts, @@ -1632,6 +1640,282 @@ def compute_nhd_routing_v02( ) ) + elif parallel_compute_method == "serial-hybrid-routing": + + results = [] + for twi, (tw, reach_list) in enumerate(reaches_bytw.items(), 1): + # The X_sub lines use SEGS... + # which becomes invalid with the wbodies included. + # So we define "common_segs" to identify regular routing segments + # and wbodies_segs for the waterbody reaches/segments + segs = list(chain.from_iterable(reach_list)) + common_segs = param_df.index.intersection(segs) + # Assumes everything else is a waterbody... + wbodies_segs = set(segs).symmetric_difference(common_segs) + + #Declare empty dataframe + waterbody_types_df_sub = pd.DataFrame() + + # If waterbody parameters exist + if not waterbodies_df.empty: + + lake_segs = list(waterbodies_df.index.intersection(segs)) + + waterbodies_df_sub = waterbodies_df.loc[ + lake_segs, + [ + "LkArea", + "LkMxE", + "OrificeA", + "OrificeC", + "OrificeE", + "WeirC", + "WeirE", + "WeirL", + "ifd", + "qd0", + "h0", + ], + ] + + #If reservoir types other than Level Pool are active + if not waterbody_types_df.empty: + waterbody_types_df_sub = waterbody_types_df.loc[ + lake_segs, + [ + "reservoir_type", + ], + ] + + else: + lake_segs = [] + waterbodies_df_sub = pd.DataFrame() + + param_df_sub = param_df.loc[ + common_segs, + ["dt", "bw", "tw", "twcc", "dx", "n", "ncc", "cs", "s0", "alt"], + ].sort_index() + + reaches_list_with_type = _build_reach_type_list(reach_list, wbodies_segs) + + # qlat_sub = qlats.loc[common_segs].sort_index() + # q0_sub = q0.loc[common_segs].sort_index() + qlat_sub = qlats.loc[param_df_sub.index] + q0_sub = q0.loc[param_df_sub.index] + + param_df_sub = param_df_sub.reindex( + param_df_sub.index.tolist() + lake_segs + ).sort_index() + + usgs_df_sub, lastobs_df_sub, da_positions_list_byseg = _prep_da_dataframes(usgs_df, lastobs_df, param_df_sub.index) + da_positions_list_byreach, da_positions_list_bygage = _prep_da_positions_byreach(reach_list, lastobs_df_sub.index) + + qlat_sub = qlat_sub.reindex(param_df_sub.index) + q0_sub = q0_sub.reindex(param_df_sub.index) + + # prepare reservoir DA data + (reservoir_usgs_df_sub, + reservoir_usgs_df_time, + reservoir_usgs_update_time, + reservoir_usgs_prev_persisted_flow, + reservoir_usgs_persistence_update_time, + reservoir_usgs_persistence_index, + reservoir_usace_df_sub, + reservoir_usace_df_time, + reservoir_usace_update_time, + reservoir_usace_prev_persisted_flow, + reservoir_usace_persistence_update_time, + reservoir_usace_persistence_index, + reservoir_rfc_df_sub, + reservoir_rfc_totalCounts, + reservoir_rfc_file, + reservoir_rfc_use_forecast, + reservoir_rfc_timeseries_idx, + reservoir_rfc_update_time, + reservoir_rfc_da_timestep, + reservoir_rfc_persist_days, + waterbody_types_df_sub, + ) = _prep_reservoir_da_dataframes( + reservoir_usgs_df, + reservoir_usgs_param_df, + reservoir_usace_df, + reservoir_usace_param_df, + reservoir_rfc_df, + reservoir_rfc_param_df, + waterbody_types_df_sub, + t0, + from_files, + ) + + # prepare diffusive input data for running hybrid routing + for diffusive_tw in diffusive_network_data: # <------- TODO - by-network parallel loop, here. + found = any(tw in sublist for sublist in reach_list) + if found: + + diffusive_segments = diffusive_network_data[diffusive_tw]['mainstem_segs'] + + # create DataFrame of junction inflow data + num_col = nts + trib_segs = diffusive_network_data[diffusive_tw]['tributary_segments'] + junction_inflows = pd.DataFrame(np.zeros((len(trib_segs), num_col)), index=trib_segs, columns=[i for i in range(num_col)]) + + if not topobathy.empty: + # create topobathy data for diffusive mainstem segments of this given tw segment + topobathy_bytw = topobathy.loc[diffusive_segments] + unrefactored_topobathy_bytw = pd.DataFrame() + else: + topobathy_bytw = pd.DataFrame() + unrefactored_topobathy_bytw = pd.DataFrame() + + # diffusive streamflow DA activation switch + if 'diffusive_streamflow_nudging' in da_parameter_dict: + diffusive_usgs_df = usgs_df.loc[usgs_df.index.isin(diffusive_segments)] + else: + diffusive_usgs_df = pd.DataFrame() + + # tw in refactored hydrofabric (NOTE: refactored_hydrofabric not in use any more) + refactored_diffusive_domain_bytw = None + refactored_reaches_byrftw = None + + # coastal boundary depth input data at TW + if tw in coastal_boundary_depth_df.index: + coastal_boundary_depth_bytw_df = coastal_boundary_depth_df.loc[diffusive_tw].to_frame().T + else: + coastal_boundary_depth_bytw_df = pd.DataFrame() + + # temporary: column names of qlats from HYfeature are currently timestamps. To be consistent with qlats from NHD + # the column names need to be changed to intergers from zero incrementing by 1 + diffusive_qlats = qlats.copy() + diffusive_qlats = diffusive_qlats.loc[diffusive_segments] + diffusive_qlats.columns = range(diffusive_qlats.shape[1]) + + # build diffusive inputs + diffusive_inputs = diff_utils.diffusive_input_data_v02( + diffusive_tw, + diffusive_network_data[diffusive_tw]['connections'], + diffusive_network_data[diffusive_tw]['rconn'], + diffusive_network_data[diffusive_tw]['reaches'], + diffusive_network_data[diffusive_tw]['mainstem_segs'], + diffusive_network_data[diffusive_tw]['tributary_segments'], + None, # place holder for diffusive parameters + diffusive_network_data[diffusive_tw]['param_df'], + diffusive_qlats, + q0, + junction_inflows, + qts_subdivisions, + t0, + nts, + dt, + waterbodies_df, + topobathy_bytw, + diffusive_usgs_df, + refactored_diffusive_domain_bytw, + refactored_reaches_byrftw, + coastal_boundary_depth_bytw_df, + unrefactored_topobathy_bytw, + ) + + # diffusive segments and reaches + diffusive_reaches=[] + diffusive_segments = diffusive_network_data[diffusive_tw]['mainstem_segs'] + nondiffusive_segments = diffusive_network_data[diffusive_tw]['tributary_segments'] + + for sublist in reach_list: + if any(item in sublist for item in diffusive_segments): + diffusive_reaches.append(sublist) + + # Compute hydraulic value lookup tables for channel cross sections + out_chxsec_lookuptable, out_z_adj= chxsec_lookuptable.compute_chxsec_lookuptable( + diffusive_inputs) + else: + diffusive_tw = None + diffusive_reaches = [] + nondiffusive_segments = [] + diffusive_inputs = {} + + # Create a dictionary mapping segment ID to a pair of Fotran segment node index and Fortran reach index + diffusive_segment2reach_and_segment_bottom_node_idx={} + total_reaches = diffusive_reaches.copy() + for seg in nondiffusive_segments: + total_reaches.append([seg]) + + for sublist in total_reaches: + first_id = sublist[0] + # find reach order index corresponding head segment id + for key, value in diffusive_inputs['pynw'].items(): + if value == first_id: + first_key = key + break + # Iterate through the sublist and populate the results dictionary + for index, id_value in enumerate(sublist): + diffusive_segment2reach_and_segment_bottom_node_idx[id_value] = [first_key, index+1] + + results.append( + compute_func( + nts, + dt, + qts_subdivisions, + reaches_list_with_type, + independent_networks[tw], + param_df_sub.index.values.astype("int64"), + param_df_sub.columns.values, + param_df_sub.values, + q0_sub.values.astype("float32"), + qlat_sub.values.astype("float32"), + lake_segs, + waterbodies_df_sub.values, + data_assimilation_parameters, + waterbody_types_df_sub.values.astype("int32"), + waterbody_type_specified, + t0.strftime('%Y-%m-%d_%H:%M:%S'), + usgs_df_sub.values.astype("float32"), + np.array(da_positions_list_byseg, dtype="int32"), + np.array(da_positions_list_byreach, dtype="int32"), + np.array(da_positions_list_bygage, dtype="int32"), + lastobs_df_sub.get("lastobs_discharge", pd.Series(index=lastobs_df_sub.index, name="Null", dtype="float32")).values.astype("float32"), + lastobs_df_sub.get("time_since_lastobs", pd.Series(index=lastobs_df_sub.index, name="Null", dtype="float32")).values.astype("float32"), + da_decay_coefficient, + # USGS Hybrid Reservoir DA data + reservoir_usgs_df_sub.values.astype("float32"), + reservoir_usgs_df_sub.index.values.astype("int32"), + reservoir_usgs_df_time.astype('float32'), + reservoir_usgs_update_time.astype('float32'), + reservoir_usgs_prev_persisted_flow.astype('float32'), + reservoir_usgs_persistence_update_time.astype('float32'), + reservoir_usgs_persistence_index.astype('float32'), + # USACE Hybrid Reservoir DA data + reservoir_usace_df_sub.values.astype("float32"), + reservoir_usace_df_sub.index.values.astype("int32"), + reservoir_usace_df_time.astype('float32'), + reservoir_usace_update_time.astype("float32"), + reservoir_usace_prev_persisted_flow.astype("float32"), + reservoir_usace_persistence_update_time.astype("float32"), + reservoir_usace_persistence_index.astype("float32"), + # RFC Reservoir DA data + reservoir_rfc_df_sub.values.astype("float32"), + reservoir_rfc_df_sub.index.values.astype("int32"), + reservoir_rfc_totalCounts.astype("int32"), + reservoir_rfc_file, + reservoir_rfc_use_forecast.astype("int32"), + reservoir_rfc_timeseries_idx.astype("int32"), + reservoir_rfc_update_time.astype("float32"), + reservoir_rfc_da_timestep.astype("int32"), + reservoir_rfc_persist_days.astype("int32"), + # Diffusive wave routing data + diffusive_tw, + diffusive_reaches, + nondiffusive_segments, + diffusive_segment2reach_and_segment_bottom_node_idx, + diffusive_inputs, + out_chxsec_lookuptable, + out_z_adj, + {}, + assume_short_ts, + return_courant, + from_files=from_files, + ) + ) + return results, subnetwork_list def compute_diffusive_routing( diff --git a/src/troute-routing/troute/routing/diffusive_utils_v02.py b/src/troute-routing/troute/routing/diffusive_utils_v02.py index 883515ba3..607e6b08c 100644 --- a/src/troute-routing/troute/routing/diffusive_utils_v02.py +++ b/src/troute-routing/troute/routing/diffusive_utils_v02.py @@ -25,7 +25,6 @@ def adj_alt1( ---------- z_all -- (dict) adjusted altitude dictionary """ - for x in range(mx_jorder, -1, -1): for head_segment, reach in ordered_reaches[x]: seg_list = reach["segments_list"] @@ -130,7 +129,7 @@ def fp_network_map( # frnw_g[frj, 2 + nusrch + r] = usrch_hseg_mainstem_j + 1 # Determine if reach being considered belong to mainstem or tributary reach as - # diffusive wave applies to mainstem reach not tributary reach + # diffusive wave applies only to mainstem reach not tributary reach if head_segment in mainstem_seg_list: frnw_g[frj,2+nusrch+1] = 555 if head_segment in trib_seg_list: @@ -246,6 +245,7 @@ def fp_qlat_map( param_df, qlat, qlat_g, + mainstem_seg_list, ): """ lateral inflow mapping between Python and Fortran @@ -275,17 +275,20 @@ def fp_qlat_map( for seg in range(0, ncomp): segID = seg_list[seg] for tsi in range(0, nts_ql_g): - if seg < ncomp - 1: - - tlf = qlat.loc[segID, tsi] - dx = param_df.loc[segID, 'dx'] - qlat_g[tsi, seg, frj] = tlf / dx # [m^2/sec] - - else: - qlat_g[ - tsi, seg, frj - ] = 0.0 # seg=ncomp is actually for bottom node in Fotran code. - # And, lateral flow enters between adjacent nodes. + if head_segment in mainstem_seg_list: + if seg < ncomp - 1: + + tlf = qlat.loc[segID, tsi] + dx = param_df.loc[segID, 'dx'] + qlat_g[tsi, seg, frj] = tlf / dx # [m^2/sec] + + else: + qlat_g[ + tsi, seg, frj + ] = 0.0 # seg=ncomp is actually for bottom node in Fotran code. + # And, lateral flow enters between adjacent nodes. + else: + qlat_g[tsi, seg, frj] = 0.0 return qlat_g def fp_ubcd_map(frnw_g, pynw, nts_ub_g, nrch_g, ds_seg, upstream_inflows): @@ -442,7 +445,6 @@ def fp_naturalxsec_map( # loop over reach orders. frj = -1 for x in range(mx_jorder, -1, -1): - # loop through all reaches of order x for head_segment, reach in ordered_reaches[x]: frj = frj + 1 @@ -610,7 +612,7 @@ def fp_coastal_boundary_input_map( dt_timeslice = timedelta(minutes=dt_db_g/60.0) tfin = t0 + dt_timeslice*(nts_db_g-1) timestamps = pd.date_range(t0, tfin, freq=dt_timeslice) - + timeslice_dbcd_list = [] tws = coastal_boundary_depth_df.index.values.tolist() for i in range(len(timestamps)): @@ -636,7 +638,7 @@ def fp_coastal_boundary_input_map( dbcd_df.loc[tw, dbcd_df.loc[tw]<= 0] = psmindepth # interpolate missing value in NaN (np.nan) inbetween available values. For extrapolation, the result is the same as - # a last available value either in the left or right. + # a last available value either to the left or to the right. dbcd_df_interpolated = dbcd_df.interpolate(axis='columns', limit_direction='both', limit=6) # if still missing data exists, do not use the coastal depth data as diffusive downstream boundary condition @@ -646,7 +648,7 @@ def fp_coastal_boundary_input_map( else: dsbd_option = 1 # use the data as prepared dbcd_g[:] = dbcd_df_interpolated.loc[tw].values - + else: dt_db_g = 3600.0 # by default in sec nts_db_g = int((tfin_g - t0_g) * 3600.0 / dt_db_g) + 1 # include initial time 0 to the final time @@ -735,7 +737,7 @@ def diffusive_input_data_v02( # timestep_ar_g[6] = dt_db_g <- defined after calling fp_coastal_boundary_input_map timestep_ar_g[7] = dt_qtrib_g timestep_ar_g[8] = dt_da_g - timestep_ar_g[9] = 10.0 # dtini_g/this_number = the min.sim.time interval internally executed within diffusive.f90 + timestep_ar_g[9] = 20.0 # divider for the initial simulation time interval to enforce a min. simulation duration. # CN-mod parameters paradim = 11 @@ -743,17 +745,21 @@ def diffusive_input_data_v02( para_ar_g[0] = 0.95 # Courant number (default: 0.95) para_ar_g[1] = 0.5 # lower limit of celerity (default: 0.5) para_ar_g[2] = 10.0 # lower limit of diffusivity (default: 10) - para_ar_g[3] = 10000.0 # upper limit of diffusivity (default: 10000) + para_ar_g[3] = 10000.0 # upper limit of diffusivity (default: 10000) para_ar_g[4] = -15.0 # lower limit of dimensionless diffusivity, used to determine b/t normal depth and diffusive depth para_ar_g[5] = -10.0 #upper limit of dimensionless diffusivity, used to determine b/t normal depth and diffusive depth para_ar_g[6] = 1.0 # 0:run Bisection to compute water level; 1: Newton Raphson (default: 1.0) - para_ar_g[7] = 0.02831 # lower limit of discharge (default: 0.02831 cms) - para_ar_g[8] = 0.0001 # lower limit of channel bed slope (default: 0.0001) + para_ar_g[7] = 0.02831 # lower limit of discharge (default: 0.02831 cms) + para_ar_g[8] = 0.00001 # lower limit of channel bed slope (default: 0.00001) para_ar_g[9] = 1.0 # weight in numerically computing 2nd derivative: 0: explicit, 1: implicit (default: 1.0) - para_ar_g[10] = 2 # downstream water depth boundary condition: 1: given water depth data, 2: normal depth + para_ar_g[10] = 2 # downstream water depth boundary condition: 1: given water depth data, 2: normal depth + # number of reaches in network nrch_g = len(reach_list) + # the number of rows in the hydraulic value lookup tables for the channel cross sections + nrow_chxsec_lookuptable = 501 + # maximum number of nodes in a reach mxncomp_g = 0 for r in reach_list: @@ -803,7 +809,6 @@ def diffusive_input_data_v02( # Order reaches by junction depth path_func = partial(nhd_network.split_at_waterbodies_and_junctions, set(junction_inflows.index.to_list()),rconn) tr = nhd_network.dfs_decomposition_depth_tuple(rconn, path_func) - jorder_reaches = sorted(tr, key=lambda x: x[0]) mx_jorder = max(jorder_reaches)[0] # maximum junction order of subnetwork of TW @@ -890,6 +895,21 @@ def diffusive_input_data_v02( # covert data type from integer to float for frnw dfrnw_g = frnw_g.astype('float') + # print python-fortran network crosswalk map + df1 = pd.DataFrame(frnw_g) + new_col_names = ['# of compute nodes', 'j of ds.reach', '# of us.reaches'] + ['j of us.reach']*(df1.shape[1]-3) + df1.columns = new_col_names + df1 = df1.rename_axis('reach j') + df1.index += 1 # be compatible with fortran number count starting from 1 + df2 = pd.DataFrame(list(pynw.items()), columns=['reach j', 'segment ID']) + df2.set_index('reach j',inplace=True) + df2.index += 1 # be compatible with fortran number count starting from 1 + df_joined = df1.merge(df2, left_index=True, right_index=True) + last_column = df_joined.columns[-1] + last_column_data = df_joined.pop(last_column) + df_joined.insert(0, last_column, last_column_data) + df_joined.to_csv('python-fortran network crosswalk map.txt', index=True) + # --------------------------------------------------------------------------------- # Step 0-5 # Prepare channel geometry data @@ -917,7 +937,9 @@ def diffusive_input_data_v02( # Step 0-6 # Prepare initial conditions data # --------------------------------------------------------------------------------- - iniq = np.zeros((mxncomp_g, nrch_g)) + iniq = np.zeros((mxncomp_g, nrch_g)) + inidepth = np.zeros((mxncomp_g, nrch_g)) + iniqpx = np.zeros((mxncomp_g, nrch_g)) frj = -1 for x in range(mx_jorder, -1, -1): for head_segment, reach in ordered_reaches[x]: @@ -955,6 +977,7 @@ def diffusive_input_data_v02( param_df, qlat, qlat_g, + mainstem_seg_list, ) # --------------------------------------------------------------------------------- @@ -1042,7 +1065,6 @@ def diffusive_input_data_v02( # Build input dictionary # --------------------------------------------------------------------------------- ntss_ev_g = int((tfin_g - t0_g) * 3600.0 / dt) + 1 - # build a dictionary of diffusive model inputs and helper variables diff_ins = {} if not refactored_diffusive_domain: @@ -1085,8 +1107,10 @@ def diffusive_input_data_v02( diff_ins["z_bathy_g"] = z_bathy_g diff_ins["mann_bathy_g"] = mann_bathy_g diff_ins["size_bathy_g"] = size_bathy_g - # initial flow value + # initial flow/depth value diff_ins["iniq"] = iniq + diff_ins["inidepth"] = inidepth + diff_ins["iniqpx"] = iniqpx # python-fortran crosswalk data diff_ins["pynw"] = pynw diff_ins["ordered_reaches"] = ordered_reaches @@ -1098,6 +1122,7 @@ def diffusive_input_data_v02( diff_ins["cwncol_g"] = crosswalk_ncol diff_ins["crosswalk_g"] = crosswalk_g diff_ins["z_thalweg_g"] = z_thalweg_g + diff_ins["nrow_chxsec_lookuptable"] = nrow_chxsec_lookuptable else: # for refactored hydrofabric # model time steps @@ -1151,6 +1176,7 @@ def diffusive_input_data_v02( diff_ins["cwncol_g"] = crosswalk_ncol diff_ins["crosswalk_g"] = crosswalk_g diff_ins["z_thalweg_g"] = z_thalweg_g + return diff_ins def unpack_output(pynw, ordered_reaches, out_q, out_elv): diff --git a/src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.pyx b/src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.pyx index b3a56e721..ac16bc779 100644 --- a/src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.pyx +++ b/src/troute-routing/troute/routing/fast_reach/chxsec_lookuptable.pyx @@ -8,30 +8,30 @@ from .fortran_wrappers cimport c_chxsec_lookuptable_calc cdef void chxsec_lookuptable( int mxncomp_g, int nrch_g, - double[::1,:] z_ar_g, - double[::1,:] bo_ar_g, - double[::1,:] traps_ar_g, - double[::1,:] tw_ar_g, - double[::1,:] twcc_ar_g, - double[::1,:] mann_ar_g, - double[::1,:] manncc_ar_g, - double[::1,:] dx_ar_g, - double so_lowerlimit_g, + float[::1,:] z_ar_g, + float[::1,:] bo_ar_g, + float[::1,:] traps_ar_g, + float[::1,:] tw_ar_g, + float[::1,:] twcc_ar_g, + float[::1,:] mann_ar_g, + float[::1,:] manncc_ar_g, + float[::1,:] dx_ar_g, + float so_lowerlimit_g, int frnw_col, int[::1,:] frnw_g, int mxnbathy_g, - double[::1,:,:] x_bathy_g, - double[::1,:,:] z_bathy_g, - double[::1,:,:] mann_bathy_g, + float[::1,:,:] x_bathy_g, + float[::1,:,:] z_bathy_g, + float[::1,:,:] mann_bathy_g, int[::1,:] size_bathy_g, int nrow_chxsec_lookuptable, - double[:,:,:,:] out_chxsec_lookuptable, - double[:,:] out_z_adj, + float[:,:,:,:] out_chxsec_lookuptable, + float[:,:] out_z_adj, ): cdef: - double[::1,:,:,:] xsec_tab = np.empty([11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g], dtype = np.double, order = 'F') - double[::1,:] z_adj = np.empty([mxncomp_g, nrch_g], dtype = np.double, order = 'F') + float[::1,:,:,:] xsec_tab = np.empty([11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g], dtype = np.float32, order = 'F') + float[::1,:] z_adj = np.empty([mxncomp_g, nrch_g], dtype = np.float32, order = 'F') c_chxsec_lookuptable_calc( &mxncomp_g, @@ -69,25 +69,25 @@ cpdef object compute_chxsec_lookuptable( cdef: int mxncomp_g = diff_inputs["mxncomp_g"] int nrch_g = diff_inputs["nrch_g"] - double[::1,:] z_ar_g = np.asfortranarray(diff_inputs["z_ar_g"]) - double[::1,:] bo_ar_g = np.asfortranarray(diff_inputs["bo_ar_g"]) - double[::1,:] traps_ar_g = np.asfortranarray(diff_inputs["traps_ar_g"]) - double[::1,:] tw_ar_g = np.asfortranarray(diff_inputs["tw_ar_g"]) - double[::1,:] twcc_ar_g = np.asfortranarray(diff_inputs["twcc_ar_g"]) - double[::1,:] mann_ar_g = np.asfortranarray(diff_inputs["mann_ar_g"]) - double[::1,:] manncc_ar_g = np.asfortranarray(diff_inputs["manncc_ar_g"]) - double[::1,:] dx_ar_g = np.asfortranarray(diff_inputs["dx_ar_g"]) - double so_lowerlimit_g = diff_inputs['para_ar_g'][8] + float[::1,:] z_ar_g = np.asfortranarray(diff_inputs["z_ar_g"].astype(np.float32)) + float[::1,:] bo_ar_g = np.asfortranarray(diff_inputs["bo_ar_g"].astype(np.float32)) + float[::1,:] traps_ar_g = np.asfortranarray(diff_inputs["traps_ar_g"].astype(np.float32)) + float[::1,:] tw_ar_g = np.asfortranarray(diff_inputs["tw_ar_g"].astype(np.float32)) + float[::1,:] twcc_ar_g = np.asfortranarray(diff_inputs["twcc_ar_g"].astype(np.float32)) + float[::1,:] mann_ar_g = np.asfortranarray(diff_inputs["mann_ar_g"].astype(np.float32)) + float[::1,:] manncc_ar_g = np.asfortranarray(diff_inputs["manncc_ar_g"].astype(np.float32)) + float[::1,:] dx_ar_g = np.asfortranarray(diff_inputs["dx_ar_g"].astype(np.float32)) + float so_lowerlimit_g = diff_inputs['para_ar_g'][8].astype(np.float32) int frnw_col = diff_inputs["frnw_col"] int[::1,:] frnw_g = np.asfortranarray(diff_inputs["frnw_g"]) int mxnbathy_g = diff_inputs['mxnbathy_g'] - double[::1,:,:] x_bathy_g = np.asfortranarray(diff_inputs["x_bathy_g"]) - double[::1,:,:] z_bathy_g = np.asfortranarray(diff_inputs["z_bathy_g"]) - double[::1,:,:] mann_bathy_g = np.asfortranarray(diff_inputs["mann_bathy_g"]) + float[::1,:,:] x_bathy_g = np.asfortranarray(diff_inputs["x_bathy_g"].astype(np.float32)) + float[::1,:,:] z_bathy_g = np.asfortranarray(diff_inputs["z_bathy_g"].astype(np.float32)) + float[::1,:,:] mann_bathy_g = np.asfortranarray(diff_inputs["mann_bathy_g"].astype(np.float32)) int[::1,:] size_bathy_g = np.asfortranarray(diff_inputs["size_bathy_g"]) int nrow_chxsec_lookuptable = diff_inputs["nrow_chxsec_lookuptable"] - double[:,:,:,:] out_chxsec_lookuptable = np.empty([11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g], dtype = np.double) - double[:,:] out_z_adj = np.empty([mxncomp_g, nrch_g], dtype = np.double) + float[:,:,:,:] out_chxsec_lookuptable = np.empty([11, nrow_chxsec_lookuptable, mxncomp_g, nrch_g], dtype = np.float32) + float[:,:] out_z_adj = np.empty([mxncomp_g, nrch_g], dtype = np.float32) # call fortran channel cross-section look up table creation subroutine chxsec_lookuptable( diff --git a/src/troute-routing/troute/routing/fast_reach/diffusive.pyx b/src/troute-routing/troute/routing/fast_reach/diffusive.pyx index e7c24075c..40d97cb74 100644 --- a/src/troute-routing/troute/routing/fast_reach/diffusive.pyx +++ b/src/troute-routing/troute/routing/fast_reach/diffusive.pyx @@ -6,7 +6,7 @@ from .fortran_wrappers cimport c_diffnw @cython.boundscheck(False) cdef void diffnw( - double[::1] timestep_ar_g, + float[::1] timestep_ar_g, int nts_ql_g, int nts_ub_g, int nts_db_g, @@ -15,45 +15,45 @@ cdef void diffnw( int nts_da_g, int mxncomp_g, int nrch_g, - double[::1,:] z_ar_g, - double[::1,:] bo_ar_g, - double[::1,:] traps_ar_g, - double[::1,:] tw_ar_g, - double[::1,:] twcc_ar_g, - double[::1,:] mann_ar_g, - double[::1,:] manncc_ar_g, - double[::1,:] so_ar_g, - double[::1,:] dx_ar_g, - double[::1,:] iniq, + float[::1,:] z_ar_g, + float[::1,:] bo_ar_g, + float[::1,:] traps_ar_g, + float[::1,:] tw_ar_g, + float[::1,:] twcc_ar_g, + float[::1,:] mann_ar_g, + float[::1,:] manncc_ar_g, + float[::1,:] so_ar_g, + float[::1,:] dx_ar_g, + float[::1,:] iniq, int frnw_col, int[::1,:] frnw_g, - double[::1,:,:] qlat_g, - double[::1,:] ubcd_g, - double[::1] dbcd_g, - double[::1,:] qtrib_g, + float[::1,:,:] qlat_g, + float[::1,:] ubcd_g, + float[::1] dbcd_g, + float[::1,:] qtrib_g, int paradim, - double[::1] para_ar_g, + float[::1] para_ar_g, int mxnbathy_g, - double[::1,:,:] x_bathy_g, - double[::1,:,:] z_bathy_g, - double[::1,:,:] mann_bathy_g, + float[::1,:,:] x_bathy_g, + float[::1,:,:] z_bathy_g, + float[::1,:,:] mann_bathy_g, int[::1,:] size_bathy_g, - double[::1,:] usgs_da_g, + float[::1,:] usgs_da_g, int[::1] usgs_da_reach_g, - double[::1,:] rdx_ar_g, + float[::1,:] rdx_ar_g, int cwnrow_g, int cwncol_g, - double[::1,:] crosswalk_g, - double[::1,:] z_thalweg_g, - double[:,:,:] out_q, - double[:,:,:] out_elv, - double[:,:,:] out_depth, + float[::1,:] crosswalk_g, + float[::1,:] z_thalweg_g, + float[:,:,:] out_q, + float[:,:,:] out_elv, + float[:,:,:] out_depth, ): cdef: - double[::1,:,:] q_ev_g = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.double, order = 'F') - double[::1,:,:] elv_ev_g = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.double, order = 'F') - double[::1,:,:] depth_ev_g = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.double, order = 'F') + float[::1,:,:] q_ev_g = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.float32, order = 'F') + float[::1,:,:] elv_ev_g = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.float32, order = 'F') + float[::1,:,:] depth_ev_g = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.float32, order = 'F') c_diffnw( ×tep_ar_g[0], @@ -112,7 +112,7 @@ cpdef object compute_diffusive( # unpack/declare diffusive input variables cdef: - double[::1] timestep_ar_g = np.asfortranarray(diff_inputs['timestep_ar_g']) + float[::1] timestep_ar_g = np.asfortranarray(diff_inputs['timestep_ar_g'].astype(np.float32)) int nts_ql_g = diff_inputs["nts_ql_g"] int nts_ub_g = diff_inputs["nts_ub_g"] int nts_db_g = diff_inputs["nts_db_g"] @@ -121,39 +121,39 @@ cpdef object compute_diffusive( int nts_da_g = diff_inputs["nts_da_g"] int mxncomp_g = diff_inputs["mxncomp_g"] int nrch_g = diff_inputs["nrch_g"] - double[::1,:] z_ar_g = np.asfortranarray(diff_inputs["z_ar_g"]) - double[::1,:] bo_ar_g = np.asfortranarray(diff_inputs["bo_ar_g"]) - double[::1,:] traps_ar_g = np.asfortranarray(diff_inputs["traps_ar_g"]) - double[::1,:] tw_ar_g = np.asfortranarray(diff_inputs["tw_ar_g"]) - double[::1,:] twcc_ar_g = np.asfortranarray(diff_inputs["twcc_ar_g"]) - double[::1,:] mann_ar_g = np.asfortranarray(diff_inputs["mann_ar_g"]) - double[::1,:] manncc_ar_g = np.asfortranarray(diff_inputs["manncc_ar_g"]) - double[::1,:] so_ar_g = np.asfortranarray(diff_inputs["so_ar_g"]) - double[::1,:] dx_ar_g = np.asfortranarray(diff_inputs["dx_ar_g"]) - double[::1,:] iniq = np.asfortranarray(diff_inputs["iniq"]) + float[::1,:] z_ar_g = np.asfortranarray(diff_inputs["z_ar_g"].astype(np.float32)) + float[::1,:] bo_ar_g = np.asfortranarray(diff_inputs["bo_ar_g"].astype(np.float32)) + float[::1,:] traps_ar_g = np.asfortranarray(diff_inputs["traps_ar_g"].astype(np.float32)) + float[::1,:] tw_ar_g = np.asfortranarray(diff_inputs["tw_ar_g"].astype(np.float32)) + float[::1,:] twcc_ar_g = np.asfortranarray(diff_inputs["twcc_ar_g"].astype(np.float32)) + float[::1,:] mann_ar_g = np.asfortranarray(diff_inputs["mann_ar_g"].astype(np.float32)) + float[::1,:] manncc_ar_g = np.asfortranarray(diff_inputs["manncc_ar_g"].astype(np.float32)) + float[::1,:] so_ar_g = np.asfortranarray(diff_inputs["so_ar_g"].astype(np.float32)) + float[::1,:] dx_ar_g = np.asfortranarray(diff_inputs["dx_ar_g"].astype(np.float32)) + float[::1,:] iniq = np.asfortranarray(diff_inputs["iniq"].astype(np.float32)) int frnw_col = diff_inputs["frnw_col"] int[::1,:] frnw_g = np.asfortranarray(diff_inputs["frnw_g"]) - double[::1,:,:] qlat_g = np.asfortranarray(diff_inputs["qlat_g"]) - double[::1,:] ubcd_g = np.asfortranarray(diff_inputs["ubcd_g"]) - double[::1] dbcd_g = np.asfortranarray(diff_inputs["dbcd_g"]) - double[::1,:] qtrib_g = np.asfortranarray(diff_inputs["qtrib_g"]) + float[::1,:,:] qlat_g = np.asfortranarray(diff_inputs["qlat_g"].astype(np.float32)) + float[::1,:] ubcd_g = np.asfortranarray(diff_inputs["ubcd_g"].astype(np.float32)) + float[::1] dbcd_g = np.asfortranarray(diff_inputs["dbcd_g"].astype(np.float32)) + float[::1,:] qtrib_g = np.asfortranarray(diff_inputs["qtrib_g"].astype(np.float32)) int paradim = diff_inputs['paradim'] - double[::1] para_ar_g = np.asfortranarray(diff_inputs["para_ar_g"]) + float[::1] para_ar_g = np.asfortranarray(diff_inputs["para_ar_g"].astype(np.float32)) int mxnbathy_g = diff_inputs['mxnbathy_g'] - double[::1,:,:] x_bathy_g = np.asfortranarray(diff_inputs["x_bathy_g"]) - double[::1,:,:] z_bathy_g = np.asfortranarray(diff_inputs["z_bathy_g"]) - double[::1,:,:] mann_bathy_g = np.asfortranarray(diff_inputs["mann_bathy_g"]) + float[::1,:,:] x_bathy_g = np.asfortranarray(diff_inputs["x_bathy_g"].astype(np.float32)) + float[::1,:,:] z_bathy_g = np.asfortranarray(diff_inputs["z_bathy_g"].astype(np.float32)) + float[::1,:,:] mann_bathy_g = np.asfortranarray(diff_inputs["mann_bathy_g"].astype(np.float32)) int[::1,:] size_bathy_g = np.asfortranarray(diff_inputs["size_bathy_g"]) - double[::1,:] usgs_da_g = np.asfortranarray(diff_inputs["usgs_da_g"]) + float[::1,:] usgs_da_g = np.asfortranarray(diff_inputs["usgs_da_g"].astype(np.float32)) int[::1] usgs_da_reach_g = np.asfortranarray(diff_inputs["usgs_da_reach_g"]) - double[::1,:] rdx_ar_g = np.asfortranarray(diff_inputs["rdx_ar_g"]) + float[::1,:] rdx_ar_g = np.asfortranarray(diff_inputs["rdx_ar_g"].astype(np.float32)) int cwnrow_g = diff_inputs["cwnrow_g"] int cwncol_g = diff_inputs["cwncol_g"] - double[::1,:] crosswalk_g = np.asfortranarray(diff_inputs["crosswalk_g"]) - double[::1,:] z_thalweg_g = np.asfortranarray(diff_inputs["z_thalweg_g"]) - double[:,:,:] out_q = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.double) - double[:,:,:] out_elv = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.double) - double[:,:,:] out_depth = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.double) + float[::1,:] crosswalk_g = np.asfortranarray(diff_inputs["crosswalk_g"].astype(np.float32)) + float[::1,:] z_thalweg_g = np.asfortranarray(diff_inputs["z_thalweg_g"].astype(np.float32)) + float[:,:,:] out_q = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.float32) + float[:,:,:] out_elv = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.float32) + float[:,:,:] out_depth = np.empty([ntss_ev_g,mxncomp_g,nrch_g], dtype = np.float32) # call diffusive compute kernel diffnw( diff --git a/src/troute-routing/troute/routing/fast_reach/diffusive_lightweight.pyx b/src/troute-routing/troute/routing/fast_reach/diffusive_lightweight.pyx new file mode 100644 index 000000000..13d36f77a --- /dev/null +++ b/src/troute-routing/troute/routing/fast_reach/diffusive_lightweight.pyx @@ -0,0 +1,161 @@ +import cython +import numpy as np +cimport numpy as np + +from .fortran_wrappers cimport c_compute_diffusive_couplingtimestep + +@cython.boundscheck(False) +cdef void diffusive_couplingtimestep( + float[::1] timestep_ar_g, + int nts_ql_g, + int nts_db_g, + int nts_qtrib_g, + int nts_da_g, + int mxncomp_g, + int nrch_g, + float[::1,:] dx_ar_g, + float[::1,:] iniq, + float[::1,:] inidepth, + float[::1,:] iniqpx, + int frnw_col, + int[::1,:] frnw_ar_g, + float[::1,:,:] qlat_g, + float[::1] dbcd_g, + float[::1,:] qtrib_g, + int paradim, + float[::1] para_ar_g, + float[::1,:] usgs_da_g, + int[::1] usgs_da_reach_g, + int nrow_chxsec_lookuptable, + float[::1,:,:,:] chxsec_lookuptable, + float[::1,:] z_adj, + float t_start, + float t_end, + float[:,:] out_q_next_out_time, + float[:,:] out_elv_next_out_time, + float[:,:] out_depth_next_out_time, + float[:,:] out_qpx_next_out_time, +): + + cdef: + float[::1,:] q_next_out_time = np.empty([mxncomp_g, nrch_g], dtype = np.float32, order = 'F') + float[::1,:] elv_next_out_time = np.empty([mxncomp_g, nrch_g], dtype = np.float32, order = 'F') + float[::1,:] depth_next_out_time = np.empty([mxncomp_g, nrch_g], dtype = np.float32, order = 'F') + float[::1,:] qpx_next_out_time = np.empty([mxncomp_g, nrch_g], dtype = np.float32, order = 'F') + + c_compute_diffusive_couplingtimestep(×tep_ar_g[0], + &nts_ql_g, + &nts_db_g, + &nts_qtrib_g, + &nts_da_g, + &mxncomp_g, + &nrch_g, + &dx_ar_g[0,0], + &iniq[0,0], + &inidepth[0,0], + &iniqpx[0,0], + &frnw_col, + &frnw_ar_g[0,0], + &qlat_g[0,0,0], + &dbcd_g[0], + &qtrib_g[0,0], + ¶dim, + ¶_ar_g[0], + &usgs_da_g[0,0], + &usgs_da_reach_g[0], + &nrow_chxsec_lookuptable, + &chxsec_lookuptable[0,0,0,0], + &z_adj[0,0], + &t_start, + &t_end, + &q_next_out_time[0,0], + &elv_next_out_time[0,0], + &depth_next_out_time[0,0], + &qpx_next_out_time[0,0], + ) + + # copy data from Fortran to Python memory view + out_q_next_out_time[:,:] = q_next_out_time[::1,:] + out_elv_next_out_time[:,:] = elv_next_out_time[::1,:] + out_depth_next_out_time[:,:] = depth_next_out_time[::1,:] + out_qpx_next_out_time[:,:] = qpx_next_out_time[::1,:] + +cpdef object compute_diffusive_couplingtimestep( + dict diff_inputs, + out_chxsec_lookuptable, + out_z_adj, + couplingtime_start, + couplingtime_end, +): + + # unpack/declare diffusive input variables + cdef: + float[::1] timestep_ar_g = np.asfortranarray(diff_inputs['timestep_ar_g'].astype(np.float32)) + int nts_ql_g = diff_inputs["nts_ql_g"] + int nts_db_g = diff_inputs["nts_db_g"] + int nts_qtrib_g = diff_inputs['nts_qtrib_g'] + int nts_da_g = diff_inputs["nts_da_g"] + int mxncomp_g = diff_inputs["mxncomp_g"] + int nrch_g = diff_inputs["nrch_g"] + float[::1,:] dx_ar_g = np.asfortranarray(diff_inputs["dx_ar_g"].astype(np.float32)) + float[::1,:] iniq = np.asfortranarray(diff_inputs["iniq"].astype(np.float32)) + float[::1,:] inidepth = np.asfortranarray(diff_inputs["inidepth"].astype(np.float32)) + float[::1,:] iniqpx = np.asfortranarray(diff_inputs["iniqpx"].astype(np.float32)) + int frnw_col = diff_inputs["frnw_col"] + int[::1,:] frnw_ar_g = np.asfortranarray(diff_inputs["frnw_g"]) + float[::1,:,:] qlat_g = np.asfortranarray(diff_inputs["qlat_g"].astype(np.float32)) + float[::1] dbcd_g = np.asfortranarray(diff_inputs["dbcd_g"].astype(np.float32)) + float[::1,:] qtrib_g = np.asfortranarray(diff_inputs["qtrib_g"].astype(np.float32)) + int paradim = diff_inputs['paradim'] + float[::1] para_ar_g = np.asfortranarray(diff_inputs["para_ar_g"].astype(np.float32)) + float[::1,:] usgs_da_g = np.asfortranarray(diff_inputs["usgs_da_g"].astype(np.float32)) + int[::1] usgs_da_reach_g = np.asfortranarray(diff_inputs["usgs_da_reach_g"]) + int nrow_chxsec_lookuptable = diff_inputs["nrow_chxsec_lookuptable"] + float[::1,:,:,:] chxsec_lookuptable = np.asfortranarray(out_chxsec_lookuptable) #.astype(np.float64)) + float[::1,:] z_adj = np.asfortranarray(out_z_adj) #.astype(np.float64)) + float t_start = float(couplingtime_start) + float t_end = float(couplingtime_end) + float[:,:] out_q_next_out_time = np.empty([mxncomp_g,nrch_g], dtype = np.float32) + float[:,:] out_elv_next_out_time = np.empty([mxncomp_g,nrch_g], dtype = np.float32) + float[:,:] out_depth_next_out_time = np.empty([mxncomp_g,nrch_g], dtype = np.float32) + float[:,:] out_qpx_next_out_time = np.empty([mxncomp_g,nrch_g], dtype = np.float32) + + # call diffusive compute kernel + diffusive_couplingtimestep( + timestep_ar_g, + nts_ql_g, + nts_db_g, + nts_qtrib_g, + nts_da_g, + mxncomp_g, + nrch_g, + dx_ar_g, + iniq, + inidepth, + iniqpx, + frnw_col, + frnw_ar_g, + qlat_g, + dbcd_g, + qtrib_g, + paradim, + para_ar_g, + usgs_da_g, + usgs_da_reach_g, + nrow_chxsec_lookuptable, + chxsec_lookuptable, + z_adj, + t_start, + t_end, + out_q_next_out_time, + out_elv_next_out_time, + out_depth_next_out_time, + out_qpx_next_out_time, + ) + + return ( + np.asarray(out_q_next_out_time), + np.asarray(out_elv_next_out_time), + np.asarray(out_depth_next_out_time), + np.asarray(out_qpx_next_out_time), + ) \ No newline at end of file diff --git a/src/troute-routing/troute/routing/fast_reach/fortran_wrappers.pxd b/src/troute-routing/troute/routing/fast_reach/fortran_wrappers.pxd index 84f6f7a09..24fb7a2af 100644 --- a/src/troute-routing/troute/routing/fast_reach/fortran_wrappers.pxd +++ b/src/troute-routing/troute/routing/fast_reach/fortran_wrappers.pxd @@ -40,7 +40,7 @@ cdef extern from "pyMCsingleSegStime_NoLoop.h": float *X) nogil; cdef extern from "pydiffusive.h": - void c_diffnw(double *timestep_ar_g, + void c_diffnw(float *timestep_ar_g, int *nts_ql_g, int *nts_ub_g, int *nts_db_g, @@ -49,39 +49,39 @@ cdef extern from "pydiffusive.h": int *nts_da_g, int *mxncomp_g, int *nrch_g, - double *z_ar_g, - double *bo_ar_g, - double *traps_ar_g, - double *tw_ar_g, - double *twcc_ar_g, - double *mann_ar_g, - double *manncc_ar_g, - double *so_ar_g, - double *dx_ar_g, - double *iniq, + float *z_ar_g, + float *bo_ar_g, + float *traps_ar_g, + float *tw_ar_g, + float *twcc_ar_g, + float *mann_ar_g, + float *manncc_ar_g, + float *so_ar_g, + float *dx_ar_g, + float *iniq, int *frnw_col, int *frnw_g, - double *qlat_g, - double *ubcd_g, - double *dbcd_g, - double *qtrib_g, + float *qlat_g, + float *ubcd_g, + float *dbcd_g, + float *qtrib_g, int *paradim, - double *para_ar_g, + float *para_ar_g, int *mxnbathy_g, - double *x_bathy_g, - double *z_bathy_g, - double *mann_bathy_g, + float *x_bathy_g, + float *z_bathy_g, + float *mann_bathy_g, int *size_bathy_g, - double *usgs_da_g, + float *usgs_da_g, int *usgs_da_reach_g, - double *rdx_ar_g, + float *rdx_ar_g, int *cwnrow_g, int *cwncol_g, - double *crosswalk_g, - double *z_thalweg_g, - double *q_ev_g, - double *elv_ev_g, - double *depth_ev_g) nogil; + float *crosswalk_g, + float *z_thalweg_g, + float *q_ev_g, + float *elv_ev_g, + float *depth_ev_g) nogil; cdef extern from "pydiffusive_cnt.h": void c_diffnw_cnt(double *dtini_g, @@ -172,22 +172,54 @@ cdef extern from "pydiffusive_cnx.h": cdef extern from "pychxsec_lookuptable.h": void c_chxsec_lookuptable_calc(int *mxncomp_g, int *nrch_g, - double *z_ar_g, - double *bo_ar_g, - double *traps_ar_g, - double *tw_ar_g, - double *twcc_ar_g, - double *mann_ar_g, - double *manncc_ar_g, - double *dx_ar_g, - double *so_lowerlimit_g, + float *z_ar_g, + float *bo_ar_g, + float *traps_ar_g, + float *tw_ar_g, + float *twcc_ar_g, + float *mann_ar_g, + float *manncc_ar_g, + float *dx_ar_g, + float *so_lowerlimit_g, int *frnw_col, int *frnw_g, int *mxnbathy_g, - double *x_bathy_g, - double *z_bathy_g, - double *mann_bathy_g, + float *x_bathy_g, + float *z_bathy_g, + float *mann_bathy_g, int *size_bathy_g, int *nrow_chxsec_lookuptable, - double *chxsec_lookuptable, - double *z_adj) nogil; \ No newline at end of file + float *chxsec_lookuptable, + float *z_adj) nogil; + +cdef extern from "pydiffusive_lightweight.h": + void c_compute_diffusive_couplingtimestep(float *timestep_ar_g, + int *nts_ql_g, + int *nts_db_g, + int *nts_qtrib_g, + int *nts_da_g, + int *mxncomp_g, + int *nrch_g, + float *dx_ar_g, + float *iniq, + float *inidepth, + float *iniqpx, + int *frnw_col, + int *frnw_ar_g, + float *qlat_g, + float *dbcd_g, + float *qtrib_g, + int *paradim, + float *para_ar_g, + float *usgs_da_g, + int *usgs_da_reach_g, + int *nrow_chxsec_lookuptable, + float *chxsec_lookuptable, + float *z_adj, + float *t_start, + float *t_end, + float *q_next_out_time, + float *elv_next_out_time, + float *depth_next_out_time, + float *qpx_next_out_time) nogil; + diff --git a/src/troute-routing/troute/routing/fast_reach/hybrid_routing_reach.pyx b/src/troute-routing/troute/routing/fast_reach/hybrid_routing_reach.pyx new file mode 100644 index 000000000..373f97d2d --- /dev/null +++ b/src/troute-routing/troute/routing/fast_reach/hybrid_routing_reach.pyx @@ -0,0 +1,899 @@ +# cython: language_level=3, boundscheck=True, wraparound=False, profile=True + +import numpy as np +import itertools +from itertools import chain +from operator import itemgetter +from array import array +from numpy cimport ndarray # TODO: Do we need to import numpy and ndarray separately? +from libc.math cimport isnan, NAN +cimport numpy as np # TODO: We are cimporting and importing numpy into the same symbol, 'np'. Problem? +cimport cython +from libc.stdlib cimport malloc, free +from libc.stdio cimport printf +#Note may get slightly better performance using cython mem module (pulls from python's heap) +#from cpython.mem cimport PyMem_Malloc, PyMem_Free +from troute.network.musking.mc_reach cimport MC_Segment, MC_Reach, _MC_Segment, get_mc_segment + +from troute.network.reach cimport Reach, _Reach, compute_type +from troute.network.reservoirs.levelpool.levelpool cimport MC_Levelpool, run_lp_c, update_lp_c +from troute.network.reservoirs.rfc.rfc cimport MC_RFC, run_rfc_c +from troute.routing.fast_reach.reservoir_hybrid_da import reservoir_hybrid_da +from troute.routing.fast_reach.reservoir_RFC_da import reservoir_RFC_da +from cython.parallel import prange + +#import cProfile +#pr = cProfile.Profile() +#NJF For whatever reason, when cimporting muskingcunge from reach, the linker breaks in weird ways +#the mc_reach.so will have an undefined symbol _muskingcunge, and reach.so will have a ____pyx_f_5reach_muskingcunge +#if you cimport reach, then call explicitly reach.muskingcung, then mc_reach.so maps to the correct module symbol +#____pyx_f_5reach_muskingcunge +#from reach cimport muskingcunge, QVD +cimport troute.routing.fast_reach.reach as reach +from troute.routing.fast_reach.simple_da cimport obs_persist_shift, simple_da_with_decay, simple_da +from troute.routing.fast_reach import diffusive +from troute.routing.fast_reach import diffusive_lightweight + +@cython.boundscheck(False) +cpdef object binary_find(object arr, object els): + """ + Find elements in els in arr. + Args: + arr: Array to search. Must be sorted + els: + Returns: + """ + cdef long hi = len(arr) + cdef object idxs = [] + + cdef Py_ssize_t L, R, m + cdef long cand, el + for el in els: + L = 0 + R = hi - 1 + m = 0 + while L <= R: + m = (L + R) // 2 + cand = arr[m] + if cand < el: + L = m + 1 + elif cand > el: + R = m - 1 + else: + break + if arr[m] == el: + idxs.append(m) + else: + raise ValueError(f"element {el} not found in {np.asarray(arr)}") + return idxs + + +@cython.boundscheck(False) +cdef void compute_reach_kernel(float qup, float quc, int nreach, const float[:,:] input_buf, float[:, :] output_buf, bint assume_short_ts, bint return_courant=False) nogil: + """ + Kernel to compute reach. + Input buffer is array matching following description: + axis 0 is reach + axis 1 is inputs in th following order: + qlat, dt, dx, bw, tw, twcc, n, ncc, cs, s0, qdp, velp, depthp + qup and quc are initial conditions. + Output buffer matches the same dimsions as input buffer in axis 0 + Input is nxm (n reaches by m variables) + Ouput is nx3 (n reaches by 3 return values) + 0: current flow, 1: current depth, 2: current velocity + """ + cdef reach.QVD rv + cdef reach.QVD *out = &rv + + cdef: + float dt, qlat, dx, bw, tw, twcc, n, ncc, cs, s0, qdp, velp, depthp + int i + + for i in range(nreach): + qlat = input_buf[i, 0] # n x 1 + dt = input_buf[i, 1] # n x 1 + dx = input_buf[i, 2] # n x 1 + bw = input_buf[i, 3] + tw = input_buf[i, 4] + twcc =input_buf[i, 5] + n = input_buf[i, 6] + ncc = input_buf[i, 7] + cs = input_buf[i, 8] + s0 = input_buf[i, 9] + qdp = input_buf[i, 10] + velp = input_buf[i, 11] + depthp = input_buf[i, 12] + + reach.muskingcunge( + dt, + qup, + quc, + qdp, + qlat, + dx, + bw, + tw, + twcc, + n, + ncc, + cs, + s0, + velp, + depthp, + out) + +# output_buf[i, 0] = quc = out.qdc # this will ignore short TS assumption at seg-to-set scale? + output_buf[i, 0] = out.qdc + output_buf[i, 1] = out.velc + output_buf[i, 2] = out.depthc + + if return_courant: + output_buf[i, 3] = out.cn + output_buf[i, 4] = out.ck + output_buf[i, 5] = out.X + + qup = qdp + + if assume_short_ts: + quc = qup + else: + quc = out.qdc + +cdef void fill_buffer_column(const Py_ssize_t[:] srows, + const Py_ssize_t scol, + const Py_ssize_t[:] drows, + const Py_ssize_t dcol, + const float[:, :] src, float[:, ::1] out) nogil except *: + + cdef Py_ssize_t i + for i in range(srows.shape[0]): + out[drows[i], dcol] = src[srows[i], scol] + +cpdef object column_mapper(object src_cols): + """Map source columns to columns expected by algorithm""" + cdef object index = {} + cdef object i_label + for i_label in enumerate(src_cols): + index[i_label[1]] = i_label[0] + + cdef object rv = [] + cdef object label + #qlat, dt, dx, bw, tw, twcc, n, ncc, cs, s0, qdp, velp, depthp + for label in ['dt', 'dx', 'bw', 'tw', 'twcc', 'n', 'ncc', 'cs', 's0']: + rv.append(index[label]) + return rv + +cpdef object compute_network_structured_with_hybrid_routing( + int nsteps, + float dt, + int qts_subdivisions, + list reaches_wTypes, # a list of tuples + dict upstream_connections, + const long[:] data_idx, + object[:] data_cols, + const float[:,:] data_values, + const float[:,:] initial_conditions, + const float[:,:] qlat_values, + list lake_numbers_col, + const double[:,:] wbody_cols, + dict data_assimilation_parameters, + const int[:,:] reservoir_types, + bint reservoir_type_specified, + str model_start_time, + const float[:,:] usgs_values, + const int[:] usgs_positions, + const int[:] usgs_positions_reach, + const int[:] usgs_positions_gage, + const float[:] lastobs_values_init, + const float[:] time_since_lastobs_init, + const double da_decay_coefficient, + const float[:,:] reservoir_usgs_obs, + const int[:] reservoir_usgs_wbody_idx, + const float[:] reservoir_usgs_time, + const float[:] reservoir_usgs_update_time, + const float[:] reservoir_usgs_prev_persisted_flow, + const float[:] reservoir_usgs_persistence_update_time, + const float[:] reservoir_usgs_persistence_index, + const float[:,:] reservoir_usace_obs, + const int[:] reservoir_usace_wbody_idx, + const float[:] reservoir_usace_time, + const float[:] reservoir_usace_update_time, + const float[:] reservoir_usace_prev_persisted_flow, + const float[:] reservoir_usace_persistence_update_time, + const float[:] reservoir_usace_persistence_index, + const float[:,:] reservoir_rfc_obs, + const int[:] reservoir_rfc_wbody_idx, + const int[:] reservoir_rfc_totalCounts, + list reservoir_rfc_file, + const int[:] reservoir_rfc_use_forecast, + const int[:] reservoir_rfc_timeseries_idx, + const float[:] reservoir_rfc_update_time, + const int[:] reservoir_rfc_da_timestep, + const int[:] reservoir_rfc_persist_days, + long diffusive_tw, + list diffusive_reaches, + list nondiffusive_segments, + dict diffusive_segment2reach_and_segment_bottom_node_idx, + dict diffusive_inputs, + const float[:,:,:,:] out_chxsec_lookuptable, + const float[:,:] out_z_adj, + dict upstream_results={}, + bint assume_short_ts=False, + bint return_courant=False, + int da_check_gage = -1, + bint from_files=True, + ): + + """ + Compute network + Args: + nsteps (int): number of time steps + reaches_wTypes (list): List of tuples: (reach, reach_type), where reach_type is 0 for Muskingum Cunge reach and 1 is a reservoir + upstream_connections (dict): Network + data_idx (ndarray): a 1D sorted index for data_values + data_values (ndarray): a 2D array of data inputs (nodes x variables) + qlats (ndarray): a 2D array of qlat values (nodes x nsteps). The index must be shared with data_values + initial_conditions (ndarray): an n x 3 array of initial conditions. n = nodes, column 1 = qu0, column 2 = qd0, column 3 = h0 + assume_short_ts (bool): Assume short time steps (quc = qup) + Notes: + Array dimensions are checked as a precondition to this method. + This version creates python objects for segments and reaches, + but then uses only the C structures and access for efficiency + """ + # Check shapes + if qlat_values.shape[0] != data_idx.shape[0]: + raise ValueError(f"Number of rows in Qlat is incorrect: expected ({data_idx.shape[0]}), got ({qlat_values.shape[0]})") + + if qlat_values.shape[1] < nsteps/qts_subdivisions: + raise ValueError(f"Number of columns (timesteps) in Qlat is incorrect: expected at most ({data_idx.shape[0]}), got ({qlat_values.shape[1]}). The number of columns in Qlat must be equal to or less than the number of routing timesteps") + + if data_values.shape[0] != data_idx.shape[0] or data_values.shape[1] != data_cols.shape[0]: + raise ValueError(f"data_values shape mismatch") + #define and initialize the final output array, add one extra time step for initial conditions + cdef int qvd_ts_w = 3 # There are 3 values per timestep (corresponding to 3 columns per timestep) + cdef np.ndarray[float, ndim=3] flowveldepth_nd = np.zeros((data_idx.shape[0], nsteps+1, qvd_ts_w), dtype='float32') + #Make ndarrays from the mem views for convience of indexing...may be a better method + cdef np.ndarray[float, ndim=2] data_array = np.asarray(data_values) + cdef np.ndarray[float, ndim=2] init_array = np.asarray(initial_conditions) + cdef np.ndarray[float, ndim=2] qlat_array = np.asarray(qlat_values) + cdef np.ndarray[double, ndim=2] wbody_parameters = np.asarray(wbody_cols) + ###### Declare/type variables ##### + # Source columns + cdef Py_ssize_t[:] scols = np.array(column_mapper(data_cols), dtype=np.intp) + cdef Py_ssize_t max_buff_size = 0 + #lists to hold reach definitions, i.e. list of ids + cdef list reach + cdef list upstream_reach + #lists to hold segment ids + cdef list segment_ids + cdef list upstream_ids + #flow accumulation variables + cdef float upstream_flows, previous_upstream_flows + #starting timestep, shifted by 1 to account for initial conditions + cdef int timestep = 1 + cdef float routing_period = dt # TODO: harmonize the dt with the value from the param_df dt (see line 153) #cdef float routing_period = data_values[0][0] + #buffers to pass to compute_reach_kernel + cdef float[:,:] buf_view + cdef float[:,:] out_buf + cdef float[:] lateral_flows + # list of reach objects to operate on + cdef list reach_objects = [] + cdef list segment_objects + cdef dict segmentid2idx={} + cdef dict segmentidx2id={} + cdef long sid + cdef _MC_Segment segment + + diffusive_segments = list(itertools.chain(*diffusive_reaches)) + + #pr.enable() + #Preprocess the raw reaches, creating MC_Reach/MC_Segments + + for reach_idx, (reach, reach_type) in enumerate(reaches_wTypes): + upstream_reach = upstream_connections.get(reach[0], ()) + upstream_ids = binary_find(data_idx, upstream_reach) + #Check if reach_type is 1 for reservoir + if (reach_type == 1): + my_id = binary_find(data_idx, reach) + wbody_index = binary_find(lake_numbers_col,reach)[0] + #Reservoirs should be singleton list reaches, TODO enforce that here? + + # write initial reservoir flows to flowveldepth array + flowveldepth_nd[my_id, 0, 0] = wbody_parameters[wbody_index, 9] # TODO ref dataframe column label list, rather than hard-coded number + + #Check if reservoir_type is not specified, then initialize default Level Pool reservoir + if (not reservoir_type_specified): + + # Initialize levelpool reservoir object + lp_obj = MC_Levelpool( + my_id[0], # index position of waterbody reach + lake_numbers_col[wbody_index], # lake number + array('l',upstream_ids), # upstream segment IDs + wbody_parameters[wbody_index], # water body parameters + reservoir_types[wbody_index][0], # waterbody type code + ) + reach_objects.append(lp_obj) + + else: + # Check whether to use the fortran or python RFC DA module: + if from_files: + # If reservoir_type is 1, 2, or 3, then initialize Levelpool reservoir + # reservoir_type 1 is a straight levelpool reservoir. + # reservoir_types 2 and 3 are USGS and USACE Hybrid reservoirs, respectively. + if (reservoir_types[wbody_index][0] >= 1 and reservoir_types[wbody_index][0] <= 3): + + # Initialize levelpool reservoir object + lp_obj = MC_Levelpool( + my_id[0], # index position of waterbody reach + lake_numbers_col[wbody_index], # lake number + array('l',upstream_ids), # upstream segment IDs + wbody_parameters[wbody_index], # water body parameters + reservoir_types[wbody_index][0], # waterbody type code + ) + reach_objects.append(lp_obj) + + #If reservoir_type is 4, then initialize RFC forecast reservoir + elif (reservoir_types[wbody_index][0] == 4 or reservoir_types[wbody_index][0] == 5): + + # Initialize rfc reservoir object + rfc_obj = MC_RFC( + my_id[0], + lake_numbers_col[wbody_index], + array('l',upstream_ids), + wbody_parameters[wbody_index], + reservoir_types[wbody_index][0], + data_assimilation_parameters["reservoir_da"]["reservoir_parameter_file"], + model_start_time, + data_assimilation_parameters["reservoir_da"]["reservoir_rfc_da"]["reservoir_rfc_forecasts_time_series_path"], + data_assimilation_parameters["reservoir_da"]["reservoir_rfc_da"]["reservoir_rfc_forecasts_lookback_hours"], + ) + reach_objects.append(rfc_obj) + + else: + + # Initialize levelpool reservoir object + lp_obj = MC_Levelpool( + my_id[0], # index position of waterbody reach + lake_numbers_col[wbody_index], # lake number + array('l',upstream_ids), # upstream segment IDs + wbody_parameters[wbody_index], # water body parameters + reservoir_types[wbody_index][0], # waterbody type code + ) + reach_objects.append(lp_obj) + + else: + segment_ids = binary_find(data_idx, reach) + + # collect network information for applying diffusive wave routing + if diffusive_tw in reach: + diffusive_tw_reach_idx = reach_idx + new_pairs = zip(reach, segment_ids) + # dictionary mapping segment id from hydrofabric to segment index determined here. + segmentid2idx.update(new_pairs) + + #Set the initial condtions before running loop + flowveldepth_nd[segment_ids, 0] = init_array[segment_ids] + segment_objects = [] + #Find the max reach size, used to create buffer for compute_reach_kernel + if len(segment_ids) > max_buff_size: + max_buff_size=len(segment_ids) + + for sid in segment_ids: + #Initialize parameters from the data_array, and set the initial initial_conditions + #These aren't actually used (the initial contions) in the kernel as they are extracted from the + #flowdepthvel array, but they could be used I suppose. Note that velp isn't used anywhere, so + #it is inialized to 0.0 + segment_objects.append( + MC_Segment(sid, *data_array[sid, scols], init_array[sid, 0], 0.0, init_array[sid, 2]) + ) + reach_objects.append( + #tuple of MC_Reach and reach_type + MC_Reach(segment_objects, array('l',upstream_ids)) + ) + + # dictionary mapping segment index determined here to segment id from hydrofabric + segmentidx2id = {value: key for key, value in segmentid2idx.items()} + + # replace initial conditions with gage observations, wherever available + cdef int gages_size = usgs_positions.shape[0] + cdef int gage_maxtimestep = usgs_values.shape[1] + cdef int gage_i, usgs_position_i + cdef float a, da_decay_minutes, da_weighted_shift, replacement_val # , original_val, lastobs_val, + cdef float [:] lastobs_values, lastobs_times + cdef (float, float, float, float) da_buf + cdef int[:] reach_has_gage = np.full(len(reaches_wTypes), np.iinfo(np.int32).min, dtype="int32") + cdef float[:,:] nudge = np.zeros((gages_size, nsteps + 1), dtype="float32") + + lastobs_times = np.full(gages_size, NAN, dtype="float32") + lastobs_values = np.full(gages_size, NAN, dtype="float32") + if gages_size: + for gage_i in range(gages_size): + lastobs_values[gage_i] = lastobs_values_init[gage_i] + lastobs_times[gage_i] = time_since_lastobs_init[gage_i] + reach_has_gage[usgs_positions_reach[gage_i]] = usgs_positions_gage[gage_i] + + if gages_size and gage_maxtimestep > 0: + for gage_i in range(gages_size): + usgs_position_i = usgs_positions[gage_i] + # Handle the instance where there are no values, only gage positions + # TODO: Compare performance with math.isnan (imported for nogil...) + if not np.isnan(usgs_values[gage_i, 0]): + flowveldepth_nd[usgs_position_i, 0, 0] = usgs_values[gage_i, 0] + + + #--------------------------------------------------------------------------------------------- + #--------------------------------------------------------------------------------------------- + + # reservoir id index arrays + cdef np.ndarray[int, ndim=1] usgs_idx = np.asarray(reservoir_usgs_wbody_idx) + cdef np.ndarray[int, ndim=1] usace_idx = np.asarray(reservoir_usace_wbody_idx) + cdef np.ndarray[int, ndim=1] rfc_idx = np.asarray(reservoir_rfc_wbody_idx) + + # reservoir update time arrays + cdef np.ndarray[float, ndim=1] usgs_update_time = np.asarray(reservoir_usgs_update_time) + cdef np.ndarray[float, ndim=1] usace_update_time = np.asarray(reservoir_usace_update_time) + cdef np.ndarray[float, ndim=1] rfc_update_time = np.asarray(reservoir_rfc_update_time) + + # reservoir persisted outflow arrays + cdef np.ndarray[float, ndim=1] usgs_prev_persisted_ouflow = np.asarray(reservoir_usgs_prev_persisted_flow) + cdef np.ndarray[float, ndim=1] usace_prev_persisted_ouflow = np.asarray(reservoir_usace_prev_persisted_flow) + + # reservoir persistence index update time arrays + cdef np.ndarray[float, ndim=1] usgs_persistence_update_time = np.asarray(reservoir_usgs_persistence_update_time) + cdef np.ndarray[float, ndim=1] usace_persistence_update_time = np.asarray(reservoir_usace_persistence_update_time) + + # reservoir persisted outflow period index + cdef np.ndarray[float, ndim=1] usgs_prev_persistence_index = np.asarray(reservoir_usgs_persistence_index) + cdef np.ndarray[float, ndim=1] usace_prev_persistence_index = np.asarray(reservoir_usace_persistence_index) + cdef np.ndarray[int, ndim=1] rfc_timeseries_idx = np.asarray(reservoir_rfc_timeseries_idx) + + #--------------------------------------------------------------------------------------------- + #--------------------------------------------------------------------------------------------- + + cdef np.ndarray fill_index_mask = np.ones_like(data_idx, dtype=bool) + cdef Py_ssize_t fill_index + cdef long upstream_tw_id + cdef dict tmp + cdef int idx + cdef float val + + for upstream_tw_id in upstream_results: + tmp = upstream_results[upstream_tw_id] + fill_index = tmp["position_index"] + fill_index_mask[fill_index] = False + for idx, val in enumerate(tmp["results"]): + flowveldepth_nd[fill_index, (idx//qvd_ts_w) + 1, idx%qvd_ts_w] = val + if data_idx[fill_index] in lake_numbers_col: + res_idx = binary_find(lake_numbers_col, [data_idx[fill_index]]) + flowveldepth_nd[fill_index, 0, 0] = wbody_parameters[res_idx, 9] # TODO ref dataframe column label + else: + flowveldepth_nd[fill_index, 0, 0] = init_array[fill_index, 0] # initial flow condition + flowveldepth_nd[fill_index, 0, 2] = init_array[fill_index, 2] # initial depth condition + + #Init buffers + lateral_flows = np.zeros( max_buff_size, dtype='float32' ) + buf_view = np.zeros( (max_buff_size, 13), dtype='float32') + out_buf = np.full( (max_buff_size, 3), -1, dtype='float32') + + cdef int num_reaches = len(reach_objects) + #Dynamically allocate a C array of reach structs + cdef _Reach* reach_structs = <_Reach*>malloc(sizeof(_Reach)*num_reaches) + + # nstep+1 to add one extra time for initial condition + cdef int num_reaches_diffusive_domain = diffusive_inputs['nrch_g'] + cdef int max_num_node_reach = diffusive_inputs['mxncomp_g'] + #cdef np.ndarray[float, ndim=3] diffusive_reach_head_flowveldepth_nd = np.zeros((diffusive_inputs['nrch_g'], nsteps+1, qvd_ts_w), dtype='float32') + #cdef float[:,:,::1] diffusive_reach_head_flowveldepth = diffusive_reach_head_flowveldepth_nd + cdef float[:,:,:] diffusive_reach_head_flowveldepth = np.zeros((num_reaches_diffusive_domain, nsteps+1, qvd_ts_w), dtype='float32', order='F') + cdef float[:,:] iniqpx_np = np.zeros((max_num_node_reach, num_reaches_diffusive_domain), dtype='float32') + cdef float inflow_to_junction + cdef float min_flow_limit + + #Populate the above array with the structs contained in each reach object + for i in range(num_reaches): + reach_structs[i] = (reach_objects[i])._reach + + #reach iterator + cdef _Reach* r + #create a memory view of the ndarray + cdef float[:,:,::1] flowveldepth = flowveldepth_nd + cdef np.ndarray[float, ndim=3] upstream_array = np.empty((data_idx.shape[0], nsteps+1, 1), dtype='float32') + #cdef np.ndarray[float, ndim=3] upstream_array = 111.0*np.ones((data_idx.shape[0], nsteps+1, 1), dtype='float32') + cdef float reservoir_outflow, reservoir_water_elevation + cdef int id = 0 + + while timestep < nsteps+1: + for i in range(num_reaches): + r = &reach_structs[i] + #Need to get quc and qup + upstream_flows = 0.0 + previous_upstream_flows = 0.0 + + for _i in range(r._num_upstream_ids):#Explicit loop reduces some overhead + id = r._upstream_ids[_i] + upstream_flows += flowveldepth[id, timestep, 0] + previous_upstream_flows += flowveldepth[id, timestep-1, 0] + segmentid = [k for k, v in segmentid2idx.items() if v ==id] + + if assume_short_ts: + upstream_flows = previous_upstream_flows + + if r.type == compute_type.RESERVOIR_LP: + + # water elevation before levelpool calculation + initial_water_elevation = r.reach.lp.water_elevation + + # levelpool reservoir storage/outflow calculation + run_lp_c(r, upstream_flows, 0.0, routing_period, &reservoir_outflow, &reservoir_water_elevation) + + # USGS reservoir hybrid DA inputs + if r.reach.lp.wbody_type_code == 2: + # find index location of waterbody in reservoir_usgs_obs + # and reservoir_usgs_time + res_idx = np.where(usgs_idx == r.reach.lp.lake_number) + wbody_gage_obs = reservoir_usgs_obs[res_idx[0][0],:] + wbody_gage_time = reservoir_usgs_time + prev_persisted_outflow = usgs_prev_persisted_ouflow[res_idx[0][0]] + persistence_update_time = usgs_persistence_update_time[res_idx[0][0]] + persistence_index = usgs_prev_persistence_index[res_idx[0][0]] + update_time = usgs_update_time[res_idx[0][0]] + + # USACE reservoir hybrid DA inputs + if r.reach.lp.wbody_type_code == 3: + # find index location of waterbody in reservoir_usgs_obs + # and reservoir_usgs_time + res_idx = np.where(usace_idx == r.reach.lp.lake_number) + wbody_gage_obs = reservoir_usace_obs[res_idx[0][0],:] + wbody_gage_time = reservoir_usace_time + prev_persisted_outflow = usace_prev_persisted_ouflow[res_idx[0][0]] + persistence_update_time = usace_persistence_update_time[res_idx[0][0]] + persistence_index = usace_prev_persistence_index[res_idx[0][0]] + update_time = usace_update_time[res_idx[0][0]] + + # Execute reservoir DA - both USGS(2) and USACE(3) types + if r.reach.lp.wbody_type_code == 2 or r.reach.lp.wbody_type_code == 3: + + #print('***********************************************************') + #print('calling reservoir DA code for lake_id:', r.reach.lp.lake_number) + #print('before DA, simulated outflow = ', reservoir_outflow) + #print('before DA, simulated water elevation = ', r.reach.lp.water_elevation) + + (new_outflow, + new_persisted_outflow, + new_water_elevation, + new_update_time, + new_persistence_index, + new_persistence_update_time + ) = reservoir_hybrid_da( + r.reach.lp.lake_number, # lake identification number + wbody_gage_obs, # gage observation values (cms) + wbody_gage_time, # gage observation times (sec) + dt * timestep, # model time (sec) + prev_persisted_outflow, # previously persisted outflow (cms) + persistence_update_time, + persistence_index, # number of sequentially persisted update cycles + reservoir_outflow, # levelpool simulated outflow (cms) + upstream_flows, # waterbody inflow (cms) + dt, # model timestep (sec) + r.reach.lp.area, # waterbody surface area (km2) + r.reach.lp.max_depth, # max waterbody depth (m) + r.reach.lp.orifice_elevation, # orifice elevation (m) + initial_water_elevation, # water surface el., previous timestep (m) + 48.0, # gage lookback hours (hrs) + update_time # waterbody update time (sec) + ) + + #print('After DA, outflow = ', new_outflow) + #print('After DA, water elevation =', new_water_elevation) + + # update levelpool water elevation state + update_lp_c(r, new_water_elevation, &reservoir_water_elevation) + + # change reservoir_outflow + reservoir_outflow = new_outflow + + #print('confirming DA elevation replacement:', reservoir_water_elevation) + #print('===========================================================') + + # update USGS DA reservoir state arrays + if r.reach.lp.wbody_type_code == 2: + usgs_update_time[res_idx[0][0]] = new_update_time + usgs_prev_persisted_ouflow[res_idx[0][0]] = new_persisted_outflow + usgs_prev_persistence_index[res_idx[0][0]] = new_persistence_index + usgs_persistence_update_time[res_idx[0][0]] = new_persistence_update_time + + # update USACE DA reservoir state arrays + if r.reach.lp.wbody_type_code == 3: + usace_update_time[res_idx[0][0]] = new_update_time + usace_prev_persisted_ouflow[res_idx[0][0]] = new_persisted_outflow + usace_prev_persistence_index[res_idx[0][0]] = new_persistence_index + usace_persistence_update_time[res_idx[0][0]] = new_persistence_update_time + + # RFC reservoir hybrid DA inputs + if r.reach.lp.wbody_type_code == 4: + # find index location of waterbody in reservoir_rfc_obs + # and reservoir_rfc_time + res_idx = np.where(rfc_idx == r.reach.lp.lake_number) + wbody_gage_obs = reservoir_rfc_obs[res_idx[0][0],:] + totalCounts = reservoir_rfc_totalCounts[res_idx[0][0]] + rfc_file = reservoir_rfc_file[res_idx[0][0]] + use_RFC = reservoir_rfc_use_forecast[res_idx[0][0]] + current_timeseries_idx = rfc_timeseries_idx[res_idx[0][0]] + update_time = rfc_update_time[res_idx[0][0]] + rfc_timestep = reservoir_rfc_da_timestep[res_idx[0][0]] + rfc_persist_days = reservoir_rfc_persist_days[res_idx[0][0]] + + # Execute RFC reservoir DA - both RFC(4) and Glacially Dammed Lake(5) types + if r.reach.lp.wbody_type_code == 4 or r.reach.lp.wbody_type_code == 5: + + #print('***********************************************************') + #print('calling reservoir DA code for lake_id:', r.reach.lp.lake_number) + #print('before DA, simulated outflow = ', reservoir_outflow) + #print('before DA, simulated water elevation = ', r.reach.lp.water_elevation) + + ( + new_outflow, + new_water_elevation, + new_update_time, + new_timeseries_idx, + dynamic_reservoir_type, + assimilated_value, + assimilated_source_file, + ) = reservoir_RFC_da( + use_RFC, # boolean whether to use RFC values or not + wbody_gage_obs, # gage observation values (cms) + current_timeseries_idx, # index of for current time series observation + totalCounts, # total number of observations in RFC timeseries + routing_period, # routing period (sec) + dt * timestep, # model time (sec) + update_time, # time to advance to next time series index + rfc_timestep, # frequency of DA observations (sec) + rfc_persist_days*24*60*60, # max seconds RFC forecasts will be used/persisted (days -> seconds) + r.reach.lp.wbody_type_code, # reservoir type + upstream_flows, # waterbody inflow (cms) + initial_water_elevation, # water surface el., previous timestep (m) + reservoir_outflow, # levelpool simulated outflow (cms) + reservoir_water_elevation, # levelpool simulated water elevation (m) + r.reach.lp.area*1.0e6, # waterbody surface area (km2 -> m2) + r.reach.lp.max_depth, # max waterbody depth (m) + rfc_file, # RFC file name + ) + + #print('After DA, outflow = ', new_outflow) + #print('After DA, water elevation =', new_water_elevation) + + # update levelpool water elevation state + update_lp_c(r, new_water_elevation, &reservoir_water_elevation) + + # change reservoir_outflow + reservoir_outflow = new_outflow + + #print('confirming DA elevation replacement:', reservoir_water_elevation) + #print('===========================================================') + + # update RFC DA reservoir state arrays + rfc_update_time[res_idx[0][0]] = new_update_time + rfc_timeseries_idx[res_idx[0][0]] = new_timeseries_idx + + + # populate flowveldepth array with levelpool or hybrid DA results + flowveldepth[r.id, timestep, 0] = reservoir_outflow + flowveldepth[r.id, timestep, 1] = 0.0 + flowveldepth[r.id, timestep, 2] = reservoir_water_elevation + upstream_array[r.id, timestep, 0] = upstream_flows + print(f"--if res_LP--") + print(f"r.id: {r.id} & timestep: {timestep}") + print(f"upstream_array[r.id, timestep, 0]: {upstream_array[r.id, timestep, 0]}") + + elif r.type == compute_type.RESERVOIR_RFC: + run_rfc_c(r, upstream_flows, 0.0, routing_period, &reservoir_outflow, &reservoir_water_elevation) + flowveldepth[r.id, timestep, 0] = reservoir_outflow + flowveldepth[r.id, timestep, 1] = 0.0 + flowveldepth[r.id, timestep, 2] = reservoir_water_elevation + upstream_array[r.id, timestep, 0] = upstream_flows + print(f"--if res_RFC--") + print(f"r.id: {r.id} & timestep: {timestep}") + print(f"upstream_array[r.id, timestep, 0]: {upstream_array[r.id, timestep, 0]}") + + elif diffusive_segments and i == diffusive_tw_reach_idx: + # run diffusive wave routing when the tailwater segment of diffusive domain is included in the current reach + + # Input forcing 1. Tributary flow: Map MC computed flow to connecting nodes between MC reaches and diffusive reaches + for seg_id in nondiffusive_segments: + seg_idx = segmentid2idx[seg_id] + reach_order_idx = diffusive_segment2reach_and_segment_bottom_node_idx[seg_id][0] + diffusive_inputs["qtrib_g"][timestep-1:, reach_order_idx] = flowveldepth[seg_idx, timestep-1, 0] + + # Input forcing 2. Initial flow and depth. Here, intial time means at each (timestep-1)*(dt or routing_period) + for seg_id in diffusive_segments: + seg_idx = segmentid2idx[seg_id] + reach_order_idx = diffusive_segment2reach_and_segment_bottom_node_idx[seg_id][0] + segment_bottom_node_idx = diffusive_segment2reach_and_segment_bottom_node_idx[seg_id][1] + + # Initial flow for reach head node + if segment_bottom_node_idx==1: + if timestep==1: + inflow_to_junction = 0.0 + for upstream_seg_id in upstream_connections[seg_id]: + upstream_seg_idx = segmentid2idx[upstream_seg_id] + inflow_to_junction += flowveldepth[upstream_seg_idx, timestep-1, 0] + diffusive_inputs["iniq"][0, reach_order_idx] = inflow_to_junction #max(inflow_to_junction, min_flow_limit) + else: + diffusive_inputs["iniq"][0, reach_order_idx] = diffusive_reach_head_flowveldepth[reach_order_idx,timestep-1,0] + + # For initial flow for all the nodes except the head node of a reach, take flow values from the previous time step + diffusive_inputs["iniq"][segment_bottom_node_idx, reach_order_idx] = flowveldepth[seg_idx, timestep-1, 0] #max(flowveldepth[seg_idx, timestep-1, 0], min_flow_limit) + + # For inidepth at head node of a reach, when there are no upstream segments or all upstream segments are MC domain, + # assume the depth of the head node of a reach is equal to that of the next node + if segment_bottom_node_idx==1: + if timestep==1: + if len(upstream_connections[seg_id])==0: + diffusive_inputs["inidepth"][0, reach_order_idx] = flowveldepth[seg_idx,timestep-1,2] + else: + for upstream_seg_id in upstream_connections[seg_id]: + upstream_seg_idx = segmentid2idx[upstream_seg_id] + if upstream_seg_id in diffusive_segments: + diffusive_inputs["inidepth"][0, reach_order_idx] = flowveldepth[upstream_seg_idx, timestep-1, 2] + break + else: + diffusive_inputs["inidepth"][0, reach_order_idx] = flowveldepth[seg_idx,timestep-1,2] + else: + diffusive_inputs["inidepth"][0, reach_order_idx] = diffusive_reach_head_flowveldepth[reach_order_idx,timestep-1,2] + + diffusive_inputs["inidepth"][segment_bottom_node_idx, reach_order_idx] = flowveldepth[seg_idx,timestep-1,2] + + # The remaining forcings, including lateral flow, coastal boundary data, and usgs data, are passed to diffusive fotran kerenl as they are. + + # Internal variables of the diffusive wave: qpx the from previous time step + for fj in range(num_reaches_diffusive_domain): + for fi in range(max_num_node_reach): + diffusive_inputs["iniqpx"][fi,fj] = iniqpx_np[fi,fj] + + (out_q_next_out_time, + out_elv_next_out_time, + out_depth_next_out_time, + out_qpx_next_out_time, + ) = diffusive_lightweight.compute_diffusive_couplingtimestep( + diffusive_inputs, + out_chxsec_lookuptable, + out_z_adj, + dt*(timestep-1)/60.0, # start time of the current coupling time step [minutes] + dt*timestep/60.0, # end time of the current coupling time step [minutes] + ) + + iniqpx_np = out_qpx_next_out_time + + for seg_id in diffusive_segments: + seg_idx = segmentid2idx[seg_id] + reach_order_idx = diffusive_segment2reach_and_segment_bottom_node_idx[seg_id][0] + segment_bottom_node_idx = diffusive_segment2reach_and_segment_bottom_node_idx[seg_id][1] + + flowveldepth[seg_idx, timestep, 0] = out_q_next_out_time[segment_bottom_node_idx, reach_order_idx] + flowveldepth[seg_idx, timestep, 2] = out_depth_next_out_time[segment_bottom_node_idx, reach_order_idx] + + # store computed diffusive flow and depth values only at the head node of a reach + diffusive_reach_head_flowveldepth[reach_order_idx,timestep,0] = out_q_next_out_time[0, reach_order_idx] + diffusive_reach_head_flowveldepth[reach_order_idx,timestep,2] = out_depth_next_out_time[0, reach_order_idx] + + else: + + #Create compute reach kernel input buffer + for _i in range(r.reach.mc_reach.num_segments): + segment = get_mc_segment(r, _i)#r._segments[_i] + buf_view[_i, 0] = qlat_array[ segment.id, ((timestep-1)/qts_subdivisions)] + buf_view[_i, 1] = segment.dt + buf_view[_i, 2] = segment.dx + buf_view[_i, 3] = segment.bw + buf_view[_i, 4] = segment.tw + buf_view[_i, 5] = segment.twcc + buf_view[_i, 6] = segment.n + buf_view[_i, 7] = segment.ncc + buf_view[_i, 8] = segment.cs + buf_view[_i, 9] = segment.s0 + buf_view[_i, 10] = flowveldepth[segment.id, timestep-1, 0] + buf_view[_i, 11] = 0.0 #flowveldepth[segment.id, timestep-1, 1] + buf_view[_i, 12] = flowveldepth[segment.id, timestep-1, 2] + + compute_reach_kernel(previous_upstream_flows, upstream_flows, + r.reach.mc_reach.num_segments, buf_view, + out_buf, + assume_short_ts) + + #Copy the output out + for _i in range(r.reach.mc_reach.num_segments): + segment = get_mc_segment(r, _i) + flowveldepth[segment.id, timestep, 0] = out_buf[_i, 0] # flow [cms] at at the current time and downstream end + if reach_has_gage[i] == da_check_gage: + printf("segment.id: %ld\t", segment.id) + printf("segment.id: %d\t", usgs_positions[reach_has_gage[i]]) + flowveldepth[segment.id, timestep, 1] = out_buf[_i, 1] # mean velocity [m/s] + flowveldepth[segment.id, timestep, 2] = out_buf[_i, 2] # mean water depth [m] + + # For each reach, + # at the end of flow calculation, Check if there is something to assimilate + # by evaluating whether the reach_has_gage array has a value different from + # the initialized value, np.iinfo(np.int32).min (the minimum possible integer). + + # TODO: If it were possible to invert the time and reach loops + # (should be possible for the MC), then this check could be + # performed fewer times -- consider implementing such a change. + + if reach_has_gage[i] > -1: + # We only enter this process for reaches where the + # gage actually exists. + # If assimilation is active for this reach, we touch the + # exactly one gage which is relevant for the reach ... + gage_i = reach_has_gage[i] + usgs_position_i = usgs_positions[gage_i] + da_buf = simple_da( + timestep, + routing_period, + da_decay_coefficient, + gage_maxtimestep, + NAN if timestep >= gage_maxtimestep else usgs_values[gage_i,timestep], + flowveldepth[usgs_position_i, timestep, 0], + lastobs_times[gage_i], + lastobs_values[gage_i], + gage_i == da_check_gage, + ) + if gage_i == da_check_gage: + printf("ts: %d\t", timestep) + printf("gmxt: %d\t", gage_maxtimestep) + printf("gage: %d\t", gage_i) + printf("old: %g\t", flowveldepth[usgs_position_i, timestep, 0]) + printf("exp_gage_val: %g\t", + NAN if timestep >= gage_maxtimestep else usgs_values[gage_i,timestep],) + + flowveldepth[usgs_position_i, timestep, 0] = da_buf[0] + + if gage_i == da_check_gage: + printf("new: %g\t", flowveldepth[usgs_position_i, timestep, 0]) + printf("repl: %g\t", da_buf[0]) + printf("nudg: %g\n", da_buf[1]) + + nudge[gage_i, timestep] = da_buf[1] + lastobs_times[gage_i] = da_buf[2] + lastobs_values[gage_i] = da_buf[3] + + timestep += 1 + + #pr.disable() + #pr.print_stats(sort='time') + #IMPORTANT, free the dynamic array created + free(reach_structs) + #slice off the initial condition timestep and return + output = np.asarray(flowveldepth[:,1:,:], dtype='float32') + #do the same for the upstream_array + output_upstream = np.asarray(upstream_array[:,1:,:], dtype='float32') + #return np.asarray(data_idx, dtype=np.intp), np.asarray(flowveldepth.base.reshape(flowveldepth.shape[0], -1), dtype='float32') + + return ( + np.asarray(data_idx, dtype=np.intp)[fill_index_mask], + output.reshape(output.shape[0], -1)[fill_index_mask], + 0, + ( + np.asarray([data_idx[usgs_position_i] for usgs_position_i in usgs_positions]), + np.asarray(lastobs_times), + np.asarray(lastobs_values) + ), + ( + usgs_idx, + usgs_update_time-((timestep-1)*dt), + usgs_prev_persisted_ouflow, + usgs_prev_persistence_index, + usgs_persistence_update_time-((timestep-1)*dt) + ), + ( + usace_idx, usace_update_time-((timestep-1)*dt), + usace_prev_persisted_ouflow, + usace_prev_persistence_index, + usace_persistence_update_time-((timestep-1)*dt) + ), + output_upstream.reshape(output.shape[0], -1)[fill_index_mask], + ( + rfc_idx, rfc_update_time-((timestep-1)*dt), + rfc_timeseries_idx + ), + np.asarray(nudge) + ) \ No newline at end of file diff --git a/src/troute-routing/troute/routing/fast_reach/pychxsec_lookuptable.h b/src/troute-routing/troute/routing/fast_reach/pychxsec_lookuptable.h index 0bd401d9c..7fcd13a23 100644 --- a/src/troute-routing/troute/routing/fast_reach/pychxsec_lookuptable.h +++ b/src/troute-routing/troute/routing/fast_reach/pychxsec_lookuptable.h @@ -1,21 +1,21 @@ extern void c_chxsec_lookuptable_calc(int *mxncomp_g, int *nrch_g, - double *z_ar_g, - double *bo_ar_g, - double *traps_ar_g, - double *tw_ar_g, - double *twcc_ar_g, - double *mann_ar_g, - double *manncc_ar_g, - double *dx_ar_g, - double *so_lowerlimit_g, + float *z_ar_g, + float *bo_ar_g, + float *traps_ar_g, + float *tw_ar_g, + float *twcc_ar_g, + float *mann_ar_g, + float *manncc_ar_g, + float *dx_ar_g, + float *so_lowerlimit_g, int *frnw_col, int *frnw_g, int *mxnbathy_g, - double *x_bathy_g, - double *z_bathy_g, - double *mann_bathy_g, + float *x_bathy_g, + float *z_bathy_g, + float *mann_bathy_g, int *size_bathy_g, int *nrow_chxsec_lookuptable, - double *chxsec_lookuptable, - double *z_adj); \ No newline at end of file + float *chxsec_lookuptable, + float *z_adj); diff --git a/src/troute-routing/troute/routing/fast_reach/pydiffusive.h b/src/troute-routing/troute/routing/fast_reach/pydiffusive.h index b455e5099..779f21c8b 100644 --- a/src/troute-routing/troute/routing/fast_reach/pydiffusive.h +++ b/src/troute-routing/troute/routing/fast_reach/pydiffusive.h @@ -1,4 +1,4 @@ -extern void c_diffnw(double *timestep_ar_g, +extern void c_diffnw(float *timestep_ar_g, int *nts_ql_g, int *nts_ub_g, int *nts_db_g, @@ -7,36 +7,36 @@ extern void c_diffnw(double *timestep_ar_g, int *nts_da_g, int *mxncomp_g, int *nrch_g, - double *z_ar_g, - double *bo_ar_g, - double *traps_ar_g, - double *tw_ar_g, - double *twcc_ar_g, - double *mann_ar_g, - double *manncc_ar_g, - double *so_ar_g, - double *dx_ar_g, - double *iniq, + float *z_ar_g, + float *bo_ar_g, + float *traps_ar_g, + float *tw_ar_g, + float *twcc_ar_g, + float *mann_ar_g, + float *manncc_ar_g, + float *so_ar_g, + float *dx_ar_g, + float *iniq, int *frnw_col, int *frnw_g, - double *qlat_g, - double *ubcd_g, - double *dbcd_g, - double *qtrib_g, + float *qlat_g, + float *ubcd_g, + float *dbcd_g, + float *qtrib_g, int *paradim, - double *para_ar_g, + float *para_ar_g, int *mxnbathy_g, - double *x_bathy_g, - double *z_bathy_g, - double *mann_bathy_g, + float *x_bathy_g, + float *z_bathy_g, + float *mann_bathy_g, int *size_bathy_g, - double *usgs_da_g, + float *usgs_da_g, int *usgs_da_reach_g, - double *rdx_ar_g, + float *rdx_ar_g, int *cwnrow_g, int *cwncol_g, - double *crosswalk_g, - double *z_thalweg_g, - double *q_ev_g, - double *elv_ev_g, - double *depth_ev_g); + float *crosswalk_g, + float *z_thalweg_g, + float *q_ev_g, + float *elv_ev_g, + float *depth_ev_g); diff --git a/src/troute-routing/troute/routing/fast_reach/pydiffusive_lightweight.h b/src/troute-routing/troute/routing/fast_reach/pydiffusive_lightweight.h new file mode 100644 index 000000000..c7c3fce09 --- /dev/null +++ b/src/troute-routing/troute/routing/fast_reach/pydiffusive_lightweight.h @@ -0,0 +1,29 @@ +extern void c_compute_diffusive_couplingtimestep(float *timestep_ar_g, + int *nts_ql_g, + int *nts_db_g, + int *nts_qtrib_g, + int *nts_da_g, + int *mxncomp_g, + int *nrch_g, + float *dx_ar_g, + float *iniq, + float *inidepth, + float *iniqpx, + int *frnw_col, + int *frnw_ar_g, + float *qlat_g, + float *dbcd_g, + float *qtrib_g, + int *paradim, + float *para_ar_g, + float *usgs_da_g, + int *usgs_da_reach_g, + int *nrow_chxsec_lookuptable, + float *chxsec_lookuptable, + float *z_adj, + float *t_start, + float *t_end, + float *q_next_out_time, + float *elv_next_out_time, + float *depth_next_out_time, + float *qpx_next_out_time); \ No newline at end of file