2nd cleanup

src/kernel/diffusive/diffusive.f90

@@ -238,8 +238,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
     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                      
-
+    dsbc_option = para_ar_g(11) ! downstream water depth boundary condition option 1:given water depth data, 2:normal depth                        
   !-----------------------------------------------------------------------------
   ! Some parameters for using natural cross section bathymetry data
     mxnbathy = mxnbathy_g  ! maximum size of bathymetry data points
@@ -324,9 +323,10 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
     z             = z_ar_g    ! node elevation array
     usgs_da_reach = usgs_da_reach_g ! contains indices of reaches where usgs data are avaliable
     usgs_da       = usgs_da_g       ! contains usgs data at a related reach 
+
   !-----------------------------------------------------------------------------
   ! variable initializations
-
+  
     routingNotChanged   = 0
     applyNaturalSection = 1
     x                   = 0.0
@@ -628,30 +628,26 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
 
               ! inflow from upstream mainstem reach
               q_usrch = newQ(frnw_g(usrchj, 1), usrchj)
-              !print *, 'added flow from upstream mainstem'
             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.0
+              tf0 = t +  dtini / 60.
               q_usrch = intp_y(nts_qtrib_g, tarr_qtrib, varr_qtrib, tf0)
-              !print *, 'added flow from tributary'
             end if
 
             ! add upstream flows to reach head
             newQ(1,j)= newQ(1,j) + q_usrch
           end do        
-        !print *, 'sum of upstream reach inflows:', newQ(1,j)
         else
 
           ! There are no links at the upstream of the reach (frnw_g(j,3)==0)
         end if
 
         ! Add lateral inflows to the reach head
         newQ(1, j) = newQ(1, j) + lateralFlow(1, j) * dx(1, j)
-        !print *, 'sum of upstream reach inflows and lateral inflows:', newQ(1,j)
 
         call mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval, j)
 
@@ -688,7 +684,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
             !do n = 1, nts_db_g  <-- already executed when initial values of Q and Y are computed.
             !  varr_db(n) = dbcd(n) + z(ncomp, j) !* when dbcd is water depth [m], channel bottom elev is added.
             !end do          
-            newY(ncomp, j) = intp_y(nts_db_g, tarr_db, varr_db, t+dtini/60.0)
+            newY(ncomp, j) = intp_y(nts_db_g, tarr_db, varr_db, t+dtini/60.)
             xcolID  = 1
             ycolID  = 2
             newArea(ncomp, j) = intp_xsec_tab(ncomp, j, nel, xcolID, ycolID, newY(ncomp,j)) ! area of normal elevation
@@ -752,8 +748,8 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
         call calc_dimensionless_numbers(j)
       end do
 
-      ! test
-      print*, "diffusive min=", t
+      ! diffusive wave simulation time print
+      print*, "diffusive simulatoin time in minute=", t
 
       ! write results to output arrays
       if ( (mod((t - t0 * 60.) * 60., saveInterval) <= TOLERANCE) .or. (t == tfin * 60.)) then
@@ -783,7 +779,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
       end if
 
       ! write initial conditions to output arrays
-      if ( ( t == t0 + dtini / 60 ) ) then
+      if ( ( t == t0 + dtini / 60. ) ) then
         do jm = 1, nmstem_rch  !* mainstem reach only
           j     = mstem_frj(jm)
           ncomp = frnw_g(j, 1)
@@ -1030,7 +1026,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
 
     ! Local variables
       integer :: ncomp
-      integer :: i, irow, flag_da
+      integer :: i, irow, flag_da, n
       double precision :: a1, a2, a3, a4
       double precision :: b1, b2, b3, b4
       double precision :: dd1, dd2, dd3, dd4
@@ -1043,10 +1039,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
       double precision :: currentQ
       double precision :: eei_ghost, ffi_ghost, exi_ghost
       double precision :: fxi_ghost, qp_ghost, qpx_ghost
-
-      ! DA
-      integer :: n 
-
+
     !-----------------------------------------------------------------------------
     !* 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
@@ -1178,14 +1171,15 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
 
       ! 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)
+        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)
+        irow = locate(tarr_da, t + dtini/60.)
         if (irow == nts_da) then
           irow = irow-1
         endif        
@@ -1196,9 +1190,11 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
         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
@@ -1219,8 +1215,9 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
       end do
 
       ! update newQ
-      newQ(1:ncomp, j) = qp(1:ncomp, 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

src/troute-nwm/src/nwm_routing/__main__.py

@@ -127,7 +127,7 @@ def nwm_route(
 
         LOG.debug("MC computation complete in %s seconds." % (time.time() - start_time_mc))
         start_time_diff = time.time()
-        #import pdb; pdb.set_trace()        
+
         # call diffusive wave simulation and append results to MC results
         results.extend(
             compute_diffusive_routing(

src/troute-nwm/src/nwm_routing/preprocess.py

@@ -152,7 +152,7 @@ def nwm_network_preprocess(
     connections, param_df, wbody_conn, gages = nnu.build_connections(
         supernetwork_parameters,
     )
-    
+
     link_gage_df = pd.DataFrame.from_dict(gages)
     link_gage_df.index.name = 'link'
 

src/troute-routing/troute/routing/compute.py

@@ -1109,7 +1109,6 @@ def compute_diffusive_routing(
     waterbodies_df,
     topobathy_data,
     ):
-
     results_diffusive = []
     for tw in diffusive_network_data: # <------- TODO - by-network parallel loop, here.
 
@@ -1162,7 +1161,7 @@ def compute_diffusive_routing(
 
         # run the simulation
         out_q, out_elv = diffusive.compute_diffusive(diffusive_inputs)
-        
+
         # unpack results
         rch_list, dat_all = diff_utils.unpack_output(
             diffusive_inputs['pynw'], 

src/troute-routing/troute/routing/diffusive_utils.py

@@ -558,9 +558,9 @@ def fp_da_map(
                 for seg in range(0, ncomp):
                     segID = seg_list[seg]
                     if segID in usgs_df_complete.index:
-                        usgs_da_g[:,frj] = usgs_df_complete.loc[segID].values
+                        usgs_da_g[:,frj] = usgs_df_complete.loc[segID].values[0:nts_da_g]
                         usgs_da_reach_g[frj] = frj + 1  # Fortran-Python index relationship, that is Python i = Fortran i+1
-                        
+
     return nts_da_g, usgs_da_g, usgs_da_reach_g
 
 def diffusive_input_data_v02(
@@ -604,7 +604,6 @@ def diffusive_input_data_v02(
     -------
     diff_ins -- (dict) formatted inputs for diffusive wave model
     """
-
     # lateral inflow timestep (sec)
     dt_ql_g = dt * qts_subdivisions
     # upstream boundary condition timestep (sec)
@@ -867,8 +866,6 @@ def diffusive_input_data_v02(
 
     # this is a place holder that uses normal depth and the lower boundary.
     # we will need to revisit this 
-    #nts_db_g = 1
-    #dbcd_g = -np.ones(nts_db_g)  
     nts_db_g = int((tfin_g - t0_g) * 3600.0 / dt_db_g)+1 
     dbcd_g = np.zeros(nts_db_g)
 
@@ -910,7 +907,6 @@ def diffusive_input_data_v02(
 
     #       Prepare interpolated USGS streamflow values at every dt_da_g time step [sec]   
     # ---------------------------------------------------------------------------------------------
-    # test: usgs_df.to_csv('./output/usgs_df.txt', header=True, index=True, sep=' ', mode='w')    
     nts_da_g, usgs_da_g, usgs_da_reach_g = fp_da_map(
                                                 mx_jorder,
                                                 ordered_reaches,
@@ -921,8 +917,7 @@ def diffusive_input_data_v02(
                                                 dt_da_g,
                                                 t0_g,
                                                 tfin_g)
-
-
+
     # TODO: Call uniform flow lookup table creation kernel    
     # ---------------------------------------------------------------------------------
     #                              Step 0-12