SimpleChannel

[rburghol]

Overview

• This will be the object for doing the simple routing algorithm. It will have settings:
  • method:
    • 0 - Euler, (code in https://github.com/gopal-bhatt/Dynamic_Model/blob/main/_str_stream/cbpm.f )
    • 1 - 3-d surface, (see Hydrology Lookup/Solver #53)
    • 2 - newtons
    • Do Error comparison amongst methods (see Error Analysis below)
• It needs Local State Variables. Creating these are outlined here: Local Named State Variables (autosetvars) #78
  • Qout
• Part of Local Channel trib network (nhd+ or custom project specific, small pump-store, etc) #16

Error Analysis

• Box-whisker plot of errors for various methods

Development

• Trying to get basic info correct

river_pid = 4431664
el_pid = 4829113

river = model_object_cache['/STATE/RCHRES_R001']
Qout = model_object_cache['/STATE/RCHRES_R001']

lc = model_data['RCHRES_R001']['local_channel']
lco = model_object_cache["/STATE/RCHRES_R001/local_channel"]
cmp = model_object_cache["/STATE/RCHRES_R001/current_monthly_pct"]
lc_Qout = lco.w_var_values['Qout']

lco.r_var_values
lco.inputs
cmp.inputs

# find all variables that have been set in state for local_channel and its children

lc_paths = [key for key, val in state_paths.items() if 'local_channel' in key]
lc_paths

[rburghol]

From fortran to python:

       D_Len = D_Len * 5280 ! miles to feet
       D_W   = ( D_BW + D_BFW ) / 2

       if ( D_Len   .le. 0 ) D_Len = 1.0
       if ( D_slope .le. 0 ) D_slope = 0.0005
       if ( D_H     .le. 0 ) D_H = 0.01
       if ( D_W     .le. 0 ) D_W = 0.01
       if ( D_n     .le. 0 ) D_n = 0.01

       D_tT = 1.0*3600 ! 1 hour in seconds 

       CritCNo = 0.999
       CritVel = CritCNo * D_Len / D_tT
       CritStg = (D_W*(CritVel*D_n)**1.5) / 
     .         ((D_W*(1.49*(D_slope)**0.5)**1.5)-(2*(CritVel*D_n)**1.5))
       if ( CritStg .le. 0 ) CritStg = 9999
       CritVol = CritStg * D_Len * D_W

       !D_storage = 0.0
       D_storage = 0.25 * D_Len * D_W * D_H
       idays =  0
       dmaxh = -9E9
       davgh = 0
       dminh = 9E9
       do it1 = 1,nvals_

          if ( D_storage .gt. dmaxh ) dmaxh = D_storage
          if ( D_storage .lt. dminh ) dminh = D_storage
          davgh = davgh + D_storage/24
          if ( MOD(it1,24) .eq. 0 ) then
             idays = idays + 1
             davgh = davgh / (D_Len * D_W)
             dminh = dminh / (D_Len * D_W)
             dmaxh = dmaxh / (D_Len * D_W)
             print*,idays,davgh,dminh,dmaxh
             dmaxh = -9E9
             davgh = 0
             dminh = 9E9
          end if

          D_tQin = D_RI_WATR_(it1)
          D_tH = D_storage / (D_Len * D_W)

          D_tU = ( 1.49 / D_n ) * 
     .           ( (D_tH*D_W/(D_W + 2*D_tH))**(2.0/3.0) ) * 
     .           D_slope ** 0.5
          D_tQ = D_tU * D_tH * D_W
          D_courant = D_tU * D_tT / D_Len

          D_tvolgrad = ( D_tQin - D_tQ ) * D_tT / D_storage ! factor not %
          if ( D_storage .eq. 0 .and. D_tQin .gt. 0 ) D_tvolgrad =  0
          D_tvol = D_storage + ( D_tQin - D_tQ ) * D_tT

!          if ( D_courant .gt. 0.99 ) then
!             D_tQ = D_RI_WATR_(it1)
!          else
!             if ( D_tH .gt. D_H ) then
!                D_tQ = D_RI_WATR_(it1) + D_tQ
!             end if
!          end if
!          D_RO_WATR_(it1) = D_tQ
!          D_storage = D_storage + (D_RI_WATR_(it1) - D_tQ)*D_tT

!v2
!          if ( D_tvolgrad .lt. -0.99 ) then 
!             Dtvol = 0
!          else
!             if ( D_tvolgrad .gt. 0 .and. Dtvol .gt. CritVol ) then
!                Dtvol = CritVol
!             end if
!          end if

          if ( D_tvolgrad .lt. -0.99 ) D_tvol = 0
          if (D_tvolgrad.gt.0 .and. D_tvol.gt.CritVol) D_tvol = CritVol
          D_tQ = D_tQin - (D_tvol - D_storage)/D_tT

          D_RO_WATR_(it1) = D_tQ
          D_storage = D_tvol

       end do