additional cleanup of routing model code:

- moved initializatio of Q_RES (GEOS_RouteGridComp.F90, reservoir.F90)
- renamed "number of" variables for clarity (EASE_pfaf_frac.F90)
- clarified *local* "stream" in comments
- edited comments
- white-space changes for improved vertical alignment

Author: gmao-rreichle

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90

@@ -1030,9 +1030,6 @@ subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC )
        allocate(WTOT_BEFORE(  n_pfaf_local))
        allocate(QINFLOW_LOCAL(n_pfaf_local))
 
-       
-       QRES_OUT    = 0.
-
        WTOT_BEFORE = WSTREAM + WRIVER + WRES
 
        ! Call river_routing_model

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90

@@ -117,6 +117,8 @@ subroutine calc( this, Q_riv_out, Q_res, Wr_res, dt, rc)
     real    :: Qin_res, alp_res, tmpfac   ! Variables for inflow, coefficients, and factors
     integer :: i, status
 
+    Q_res = 0.
+
     do i = 1, size(this%active_res)
 
        if (this%active_res(i) == 1) then                  ! reservoir is active

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/routing_model.F90

@@ -28,6 +28,10 @@ END FUNCTION RESCONST
 ! -------------------------------------------------------------------------------------------------------
 
   RECURSIVE SUBROUTINE SEARCH_DNST (K, NCAT_G, DNST, Pfaf_all, DNST_OUT)
+
+    ! objective:  [best guess by rreichle as of 3 Feb 2026]
+    ! find Pfafstetter index of downstream (DNST) catchment, where "downstream" is really "downriver" in the
+    !   jargon of the routing model with "main rivers" and "local streams"  
 
     implicit none
 
@@ -64,24 +68,24 @@ END SUBROUTINE SEARCH_DNST
   ! Routing Model Input Parameters
   ! ------------------------------
   !**** NCAT       = NUMBER OF CATCHMENTS IN THE STUDY DOMAIN
-  !**** Qrunf0     = RUNOFF PRODUCED BY LAND SURFACE MODEL IN THE CATCHMENT        [m^3/s]
-  !**** llc_ori    = MAIN RIVER LENGTH SCALE                                       [m]
-  !**** lstr       = LOCAL STREAMS LENGTH SCALE                                    [m]
-  !**** qstr_clmt0 = CLIMATOLOGY RUNOFF                                            [m^3/s]  
-  !**** qri_clmt0  = CLIMATOLOGY DISCHAR                                           [m^3/s]
-  !**** qin_clmt0  = CLIMATOLOGY INFLOW                                            [m^3/s]               
-  !**** K          = K PARAMETER FOR MAIN RIVER
-  !**** Kstr0      = K PARAMETER FOR LOCAL STREAM                                  [m^3/s]
-  
-  ! Routing Model Prognostics
-  ! -------------------------
-  !**** Ws0        = AMOUNT OF WATER IN "LOCAL STREAM"                             [m^3]
-  !**** Wr0        = AMOUNT OF WATER IN RIVER                                      [m^3]
+  !**** Qrunf0     = RUNOFF PRODUCED BY LAND SURFACE MODEL IN THE CATCHMENT            [m^3/s]
+  !**** llc_ori    = MAIN RIVER LENGTH SCALE                                           [m]
+  !**** lstr       = LOCAL STREAMS LENGTH SCALE                                        [m]
+  !**** qstr_clmt0 = CLIMATOLOGY RUNOFF                                                [m^3/s]  
+  !**** qri_clmt0  = CLIMATOLOGY DISCHAR                                               [m^3/s]
+  !**** qin_clmt0  = CLIMATOLOGY INFLOW                                                [m^3/s]               
+  !**** K          = K PARAMETER FOR MAIN RIVER                                        
+  !**** Kstr0      = K PARAMETER FOR LOCAL STREAM                                      [m^3/s]
+                                                                                       
+  ! Routing Model Prognostics                                                          
+  ! -------------------------                                                          
+  !**** Ws0        = AMOUNT OF WATER IN "LOCAL STREAM"                                 [m^3]
+  !**** Wr0        = AMOUNT OF WATER IN RIVER                                          [m^3]
 
   ! Routing Model Diagnostics
   ! -------------------------
-  !**** QS         = TRANSFER OF MOISTURE FROM STREAM VARIABLE TO RIVER VARIABLE   [m^3/s]
-  !**** QOUT       = TRANSFER OF RIVER WATER TO THE DOWNSTREAM CATCHMENT           [m^3/s]  
+  !**** QS         = TRANSFER OF MOISTURE FROM STREAM VARIABLE TO RIVER VARIABLE       [m^3/s]
+  !**** QOUT       = TRANSFER OF RIVER WATER TO THE DOWNSTREAM (DOWNRIVER) CATCHMENT   [m^3/s]  
 
   SUBROUTINE RIVER_ROUTING_HYD (          &
        NCAT,ROUTE_DT,                     &
@@ -137,13 +141,13 @@ SUBROUTINE RIVER_ROUTING_HYD (          &
 
     ! Update state variables: ks, Ws, and Qs 
     where(Qrunf<=small)Qrunf=0.                                ! Set runoff to zero if it's too small
-    Qs0=max(0.,alp_s * Ws**(1./(1.-RRM_mm)))                       ! Initial flow from stream storage (kg/s)
+    Qs0=max(0.,alp_s * Ws**(1./(1.-RRM_mm)))                       ! Initial flow from local stream storage (kg/s)
     ks =max(0.,(alp_s/(1.-RRM_mm)) * Ws**(RRM_mm/(1.D0-RRM_mm)))           ! Flow coefficient (s^-1)
     Ws_last=Ws                                                 ! Store the current water storage 
     where(ks>small)  Ws=Ws + (Qrunf-Qs0)/ks*(1.-exp(-ks*dt))   ! Update storage (kg)
     where(ks<=small) Ws=Ws + (Qrunf-Qs0)*dt                    ! Simplified update if ks is small
     Ws=max(0.,Ws)                                              ! Ensure storage is non-negative
-    Qs=max(0.,Qrunf-(Ws-Ws_last)/dt)                           ! Calculate the stream flow (kg/s)
+    Qs=max(0.,Qrunf-(Ws-Ws_last)/dt)                           ! Calculate the local stream flow (kg/s)
 
     ! Calculate variables related to river routing: Qr0, kr
     Wr=Wr+Qs*dt

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_pfaf_frac.F90

@@ -6,7 +6,7 @@ module EASE_pfaf_fracMod
   !               writes information into EASEv[x]_tile2pfaf.nc4
 
   use MAPL
-  use LogRectRasterizeMod, only: nc=>SRTM_maxcat
+  use LogRectRasterizeMod, only: nPfaf=>SRTM_maxcat
 
   implicit none
   private
@@ -68,11 +68,11 @@ subroutine EASE_get_Pfaf_frac(file_out, BCS_PATH, GridName)
     allocate(lon(nlon), lat(nlat))
     allocate(loni(nlon), lati(nlat),loni2(nlon), lati2(nlat))
     allocate(lons(nc_ease),lats(nr_ease))
-    allocate(nsub(nc))
+    allocate(nsub(nPfaf))
 
     call formatter%open(trim(pfafData_file), PFIO_READ)
-    call formatter%get_var("latitude", lat)
-    call formatter%get_var("longitude", lon)
+    call formatter%get_var("latitude",   lat)
+    call formatter%get_var("longitude",  lon)
     call formatter%get_var("CatchIndex", catchind)
     call formatter%close()
 
@@ -92,7 +92,7 @@ subroutine EASE_get_Pfaf_frac(file_out, BCS_PATH, GridName)
     call formatter%close()
     cellarea = cellarea / 1.e6  ! Convert cell area (e.g., from m^2 to km^2)
     ! Initialize aggregation arrays to zero:
-    allocate(xsub0(9999, nc), ysub0(9999, nc), asub0(9999, nc))
+    allocate(xsub0(9999, nPfaf), ysub0(9999, nPfaf), asub0(9999, nPfaf))
     nsub = 0;xsub0 = 0;ysub0 = 0;asub0 = 0.
     ! Loop over all raster grid cells to aggregate cell areas by catchment and sub-area:
     do xi = 1, nlon
@@ -124,52 +124,55 @@ subroutine EASE_get_Pfaf_frac(file_out, BCS_PATH, GridName)
     end do
     nmax = maxval(nsub)
     !print *,nmax
-    allocate(xsub(nmax, nc), ysub(nmax, nc), asub(nmax, nc))
+    allocate(xsub(nmax, nPfaf), ysub(nmax, nPfaf), asub(nmax, nPfaf))
     xsub=xsub0(1:nmax,:)
     ysub=ysub0(1:nmax,:)
     asub=asub0(1:nmax,:)
     deallocate(xsub0,ysub0,asub0)   
     ! Open the catchment definition file for the EASE grid and read the total number of tiles (header)
     open(77, file="clsm/catchment.def");read(77, *) ntot
-    ! Allocate arrays with size nt
+    ! Allocate arrays with size ntot
     allocate(tid(ntot), catid(ntot), lon_left(ntot), lon_right(ntot), lat_bottom(ntot), lat_top(ntot),latc(ntot),lonc(ntot))
-    allocate(lati_tile(ntot),loni_tile(ntot),map_tile(nc_ease,nr_ease))
-    ! Loop over each catchment and read: id, catchment id, left/right longitudes, bottom/top latitudes
+    allocate(lati_tile(ntot),loni_tile(ntot))
+    ! Loop through tiles and read: id, catchment id, left/right longitudes, bottom/top latitudes
     do i = 1, ntot
-      read(77, *) tid(i), catid(i), lon_left(i), lon_right(i), lat_bottom(i), lat_top(i)
+       read(77, *) tid(i), catid(i), lon_left(i), lon_right(i), lat_bottom(i), lat_top(i)
     end do
+    ! get center (not "center-of-mass") lat/lon of each tile based on tile min/max lat/lon
     latc = (lat_bottom + lat_top) / 2.
     lonc = (lon_left + lon_right) / 2.
     call nearest_index_vector(latc,lats,lati_tile)
     call nearest_index_vector(lonc,lons,loni_tile)
+    ! record into map_tile
+    allocate(map_tile(nc_ease,nr_ease))
     map_tile=-9999
     do i =1, ntot
       map_tile(loni_tile(i),lati_tile(i)) = i
     enddo
 
-    allocate(isub(nmax, nc))
+    allocate(isub(nmax, nPfaf))
     isub=0
     ! Loop over each catchment
-    do i = 1, nc
+    do i = 1, nPfaf
       ! Loop over each potential sub-area within the current catchment
       do j = 1, nmax
-        xi = xsub(j, i)  ! Retrieve the x-coordinate for the sub-area
-        yi = ysub(j, i)  ! Retrieve the y-coordinate for the sub-area
-        if (xi /= 0) then  ! Check if a valid sub-area exists (non-zero x-coordinate)
+        xi = xsub(j, i)                  ! Retrieve the x-coordinate for the sub-area
+        yi = ysub(j, i)                  ! Retrieve the y-coordinate for the sub-area
+        if (xi /= 0) then                ! Check if a valid sub-area exists (non-zero x-coordinate)
           isub(j, i) = map_tile(xi, yi)  ! Map the sub-area coordinates to a global tile index using map_tile
         endif
       enddo
     enddo
 
-    allocate(area_cat(nc),frac(nmax,nc),nsub_tile(ntot),flag2(nmax,nc))
+    allocate(area_cat(nPfaf),frac(nmax,nPfaf),nsub_tile(ntot),flag2(nmax,nPfaf))
     where(isub==-9999) asub=0.
     area_cat=0.
-    do i=1,nc
+    do i=1,nPfaf
       area_cat(i)=sum(asub(:,i))
     enddo
     nsub_tile=0
     frac=0.
-    do i=1,nc
+    do i=1,nPfaf
       do j=1,nmax
         if(isub(j,i)>0)then
           it=isub(j,i)
@@ -184,7 +187,7 @@ subroutine EASE_get_Pfaf_frac(file_out, BCS_PATH, GridName)
     csub=0
     nsub_tile=0
     frac_tile=0.
-    do i=1,nc
+    do i=1,nPfaf
       do j=1,nmax
         if(isub(j,i)>0)then
           it=isub(j,i)
@@ -238,19 +241,20 @@ end subroutine EASE_get_Pfaf_frac
 
   subroutine nearest_index_vector(candidates, targets, idx)
     ! For each targets(k), find argmin_j |candidates(j) - targets(k)|
-    real*8, intent(in)  :: targets(:)
-    real*8, intent(in)  :: candidates(:)
+    real*8,  intent(in)  :: targets(:)
+    real*8,  intent(in)  :: candidates(:)
     integer, intent(out) :: idx(size(candidates))   ! 1-based indices
-    integer :: k, j, nT, nC, best_j
-    real*8 :: best_d, d
 
-    nT = size(targets)
-    nC = size(candidates)
+    integer :: k, j, nTarg, nCand, best_j
+    real*8  :: best_d, d
 
-    do k = 1, nC
+    nTarg = size(targets)
+    nCand = size(candidates)
+
+    do k = 1, nCand
        best_d = huge(1.0d0)
        best_j = 1
-       do j = 1, nT
+       do j = 1, nTarg
           d = abs(candidates(k) - targets(j))
           if (d < best_d) then
              best_d = d
@@ -260,5 +264,5 @@ subroutine nearest_index_vector(candidates, targets, idx)
        idx(k) = best_j    ! already 1-based
     end do
   end subroutine nearest_index_vector
-
+  
 end module EASE_pfaf_fracMod

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_latloni_cellarea.py

@@ -52,7 +52,7 @@ def area_global_rectilinear_grid(lat, lon, rearth=6371.):
     lon : numpy.ndarray
         Array of longitude values (degrees).
     rearth : float
-        Earth radius in kilometers. Default is 6371.22 km.
+        Earth radius in kilometers. Default is 6371. km (for consistency with MAPL_Radius as of 3 Feb 2026).
         
     Returns:
     area_grid : numpy.ndarray

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_river_length.f90

@@ -232,13 +232,13 @@ function spherical_distance(lon_dn, lat_dn, lon_up, lat_up) result(distance)
     !   distance       - Great-circle distance between the two points (kilometers)
     !------------------------------------------------------------
     real, intent(in) :: lon_dn, lat_dn   ! Coordinates of downstream point
-    real, intent(in) :: lon_up, lat_up     ! Coordinates of upstream point
+    real, intent(in) :: lon_up, lat_up   ! Coordinates of upstream   point
     real :: distance                     ! Computed distance (km)
     real :: R, dlon, dlat, a, c          ! Intermediate variables
 
-    R = 6371.0                         ! Earth's radius in kilometers
+    R = 6371.0                                       ! Earth's radius in kilometers
     dlon = (lon_up - lon_dn) * (acos(-1.0) / 180.0)  ! Delta longitude (radians)
-    dlat = (lat_up - lat_dn) * (acos(-1.0) / 180.0)  ! Delta latitude (radians)
+    dlat = (lat_up - lat_dn) * (acos(-1.0) / 180.0)  ! Delta latitude  (radians)
 
     a = sin(dlat / 2.0)**2 + cos(lat_dn * (acos(-1.0) / 180.0)) * &
         cos(lat_up * (acos(-1.0) / 180.0)) * sin(dlon / 2.0)**2
@@ -247,4 +247,4 @@ function spherical_distance(lon_dn, lat_dn, lon_up, lat_up) result(distance)
 
 end function spherical_distance
 
-end program main
+end program main

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/k_module_cali.f90

@@ -378,10 +378,10 @@ subroutine get_station_inf(file_pfafmap, ns, nc, nlat, nlon, lati, loni, catid,
     ! Calculate modeled K values for all catchments
     KImodel_all = (Qclmt_all**(exp_clmt)) * (slp_all**(exp_slp))
 
-    ! Calculate stream K values using the scaling factor
+    ! Calculate local stream K values using the scaling factor
     Kstr_all = (Qstr_all**(exp_clmt)) * (slp_all**(exp_slp))
 
-    ! Write stream K values to an output file
+    ! Write local stream K values to an output file
     open(88, file="temp/Pfaf_Kstr_PR_fac1_0p35_0p45_0p2_n0p2.txt")
     do i = 1, nc
       write(88, *) Kstr_all(i)
@@ -1026,4 +1026,4 @@ subroutine sort2(arr)
   end subroutine sort2
 
 !------------------------------------------------------------
-end module k_module
+end module k_module