Add river routing grid comp (#145)

Add  river routing grid comp

Author: gmao-rreichle

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added support for river routing.
 - Added optional NetCDF4 output mode for ObsFcstAna, including NetCDF metadata and runtime context. Changed namelist variable "out_ObsFcstAna" from logical to integer.
 
 ### Changed

CMakeLists.txt

@@ -11,7 +11,7 @@ esma_add_library(${this}
   SRCS GEOS_LdasGridComp.F90
   SUBCOMPONENTS ${alldirs}
   SUBDIRS LDAS_Shared
-  DEPENDENCIES GEOSland_GridComp GEOSlandice_GridComp makebcs MAPL
+  DEPENDENCIES GEOSland_GridComp GEOSlandice_GridComp GEOSroute_GridComp makebcs MAPL
   INCLUDES ${INC_ESMF})
 
 esma_add_subdirectory(GEOSldas_App)

GEOS_LdasGridComp.F90

@@ -15,6 +15,7 @@ module GEOS_LdasGridCompMod
   use GEOS_EnsGridCompMod,       only: EnsSetServices       => SetServices
   use GEOS_LandAssimGridCompMod, only: LandAssimSetServices => SetServices
   use GEOS_LandiceGridCompMod,   only: LandiceSetServices   => SetServices
+  use GEOS_RouteGridCompMod,     only: RouteSetServices     => SetServices
 
   use LDAS_TileCoordType,     only: tile_coord_type , T_TILECOORD_STATE, TILECOORD_WRAP
   use LDAS_TileCoordType,     only: grid_def_type, io_grid_def_type, operator (==)
@@ -51,6 +52,7 @@ module GEOS_LdasGridCompMod
   ! All children
   integer,allocatable :: LAND(:)
   integer,allocatable :: LANDICE(:)
+  integer,allocatable :: ROUTE(:)
   integer,allocatable :: LANDPERT(:)
   integer,allocatable :: METFORCE(:)
   integer             :: ENSAVG, LANDASSIM
@@ -62,6 +64,7 @@ module GEOS_LdasGridCompMod
   logical :: ensemble_forcing   ! switch between deterministic and ensemble forcing
   logical :: with_landice       ! true if landice tiles requested by config
   logical :: with_land          ! true if land    tiles requested by config
+  integer :: RUN_ROUTE          ! 0:, no river routing, 1: routing w/o reservoirs, 2: routing w/ reservoirs
 
 contains
 
@@ -143,13 +146,15 @@ subroutine SetServices(gc, rc)
     !create ensemble children
     call MAPL_GetObjectFromGC(gc, MAPL, rc=status)
     VERIFY_(status)
-    call MAPL_GetResource ( MAPL, NUM_ENSEMBLE, Label="NUM_LDAS_ENSEMBLE:", DEFAULT=1,       RC=STATUS)
+    call MAPL_GetResource ( MAPL, NUM_ENSEMBLE,    Label="NUM_LDAS_ENSEMBLE:", DEFAULT=1,    RC=STATUS)
     VERIFY_(STATUS)
-    call MAPL_GetResource ( MAPL, ens_id_width, Label="ENS_ID_WIDTH:",      DEFAULT=0,       RC=STATUS)
+    call MAPL_GetResource ( MAPL, ens_id_width,    Label="ENS_ID_WIDTH:",      DEFAULT=0,    RC=STATUS)
     VERIFY_(STATUS)
-    call MAPL_GetResource ( MAPL, FIRST_ENS_ID, Label="FIRST_ENS_ID:", DEFAULT=0, RC=STATUS)
+    call MAPL_GetResource ( MAPL, RUN_ROUTE,       Label="RUN_ROUTE:",         DEFAULT=0,    RC=STATUS)
+    VERIFY_(STATUS)
+    call MAPL_GetResource ( MAPL, FIRST_ENS_ID,    Label="FIRST_ENS_ID:",      DEFAULT=0,    RC=STATUS)
     VERIFY_(STATUS)
-    call MAPL_GetResource ( MAPL, ENS_FORCING_STR, Label="ENSEMBLE_FORCING:", DEFAULT="NO", RC=STATUS)
+    call MAPL_GetResource ( MAPL, ENS_FORCING_STR, Label="ENSEMBLE_FORCING:",  DEFAULT="NO", RC=STATUS)
     VERIFY_(STATUS)
     ENS_FORCING_STR =  ESMF_UtilStringUpperCase(ENS_FORCING_STR, rc=STATUS)
     VERIFY_(STATUS)
@@ -174,6 +179,10 @@ subroutine SetServices(gc, rc)
     if (any(tile_types == MAPL_LAND   )) with_land    = .true.
 !   if (any(tile_types == MAPL_LAKE   )) with_lake    = .true.
 
+    if (NUM_ENSEMBLE>1) then
+       _ASSERT( .not. (with_landice .or. RUN_ROUTE>0), "Landice and route not supported in ensemble mode.")
+    endif
+
     call MAPL_GetResource ( MAPL, LAND_ASSIM_STR, Label="LAND_ASSIM:", DEFAULT="NO", RC=STATUS)
     VERIFY_(STATUS)
     LAND_ASSIM_STR =  ESMF_UtilStringUpperCase(LAND_ASSIM_STR, rc=STATUS)
@@ -205,7 +214,10 @@ subroutine SetServices(gc, rc)
 
     if (with_land)    allocate(LAND(   NUM_ENSEMBLE),LANDPERT(NUM_ENSEMBLE))
     if (with_landice) allocate(LANDICE(NUM_ENSEMBLE))
-
+    if (RUN_ROUTE >= 1) then
+       _ASSERT( with_land, "RUNOFF must be from the export of land_gridcomp for now.")
+       allocate(ROUTE(NUM_ENSEMBLE))
+    endif
     ! ens_id_with = 2 + number of digits = total number of chars in ensid_string ("_eXXXX")
     !
     ! Assert ens_id_width<=2+9 so number of digits remains single-digit and "I1" can be
@@ -253,6 +265,12 @@ subroutine SetServices(gc, rc)
           LANDICE(i) = MAPL_AddChild(gc, name=childname, ss=LandiceSetServices, rc=status)
           VERIFY_(status)
        endif
+
+      if (RUN_ROUTE >= 1) then
+          childname='ROUTE'//trim(ensid_string)
+          ROUTE(i) = MAPL_AddChild(gc, name=childname, ss=RouteSetServices, rc=status)
+          VERIFY_(status)
+      endif
     enddo
 
     if (with_land) then
@@ -302,6 +320,17 @@ subroutine SetServices(gc, rc)
                rc = status                                                                            &
                )
           VERIFY_(status)
+
+          IF(RUN_ROUTE >= 1) THEN
+             call MAPL_AddConnectivity (                              &
+                  GC                                                 ,&
+                  SHORT_NAME  = (/'RUNOFF  '/)                       ,&   ! RUNOFF = total runoff = surface runoff + baseflow
+                  SRC_ID =  LAND(I)                                  ,&
+                  DST_ID =  ROUTE(I)                                 ,&
+                  RC=STATUS )
+             VERIFY_(STATUS)
+          ENDIF
+
        enddo
 
        if(land_assim .or. mwRTM) then
@@ -788,6 +817,13 @@ subroutine Initialize(gc, import, export, clock, rc)
           call MAPL_Set(CHILD_MAPL, LocStream=landice_locstream, rc=status)
           VERIFY_(status)
        endif
+
+       if (RUN_ROUTE >= 1) then
+          call MAPL_GetObjectFromGC(gcs(ROUTE(i)), CHILD_MAPL, rc=status)
+          VERIFY_(status) ! CHILD = ens_avg
+          call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
+          VERIFY_(status)
+       endif
 
     enddo
 
@@ -1016,6 +1052,16 @@ subroutine Run(gc, import, export, clock, rc)
           call MAPL_TimerOff(MAPL, gcnames(igc))
        endif ! with_land_ice
 
+       if ( RUN_ROUTE >= 1 ) then
+          igc = ROUTE(i)
+          call MAPL_TimerOn(MAPL, gcnames(igc))
+          call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status)
+          VERIFY_(status)
+          call ESMF_GridCompRun(gcs(ENSAVG), importState=gex(igc), exportState=gex(ENSAVG), clock=clock,phase=4, userRC=status)
+          VERIFY_(status)
+          call MAPL_TimerOff(MAPL, gcnames(igc))
+       endif ! river-routine
+
        if (with_land) then
           ! ApplyPrognPert - moved: now before calculating ensemble average that is picked up by land analysis and HISTORY; reichle 28 May 2020 
           igc = LANDPERT(i)

GEOSens_GridComp/GEOS_EnsGridComp.F90

@@ -2041,9 +2041,6 @@ subroutine Collect_force_ens(gc, import, export, clock, rc)
     character(len=ESMF_MAXSTR) :: Iam
     character(len=ESMF_MAXSTR) :: comp_name
 
-    ! ESMF variables
-    type(ESMF_VM) :: vm
-
     ! MAPL variables
     type(MAPL_MetaComp), pointer :: MAPL=>null() ! MAPL obj
 
@@ -2264,9 +2261,6 @@ subroutine Collect_land_ens(gc, import, export, clock, rc)
     character(len=ESMF_MAXSTR) :: Iam
     character(len=ESMF_MAXSTR) :: comp_name
 
-    ! ESMF variables
-    type(ESMF_VM) :: vm
-
     ! MAPL variables
     type(MAPL_MetaComp), pointer :: MAPL=>null() ! MAPL obj
 

GEOSldas_App/GEOSldas_HIST.rc

@@ -24,6 +24,7 @@ COLLECTIONS:
 #                'const_2d_lnd_Nx'
 #                'tavg24_2d_glc_Nx'
 #                'tavg24_1d_glc_Nt'
+#                'tavg24_1d_route'
            ::
 
 # --------------------------------------------------------------------------------------------------
@@ -580,4 +581,20 @@ EASEv2_M36.LM:           1
                                 'WESNBOT'     , 'LANDICE' ,
                                 'WESNEXT'     , 'LANDICE' ,
                                  ::
+
+  tavg24_1d_route.descr:          'Catchment-space,Daily,Time-Averaged,Single-level,Runoff Routing Diagnostics',
+  tavg24_1d_route.nbits:          12,
+  tavg24_1d_route.template:       '%y4%m2%d2_%h2%n2z.nc4',
+  tavg24_1d_route.archive:        '%c/Y%y4',
+  tavg24_1d_route.format:         'CFIO',
+  tavg24_1d_route.mode:           'time-averaged',
+  tavg24_1d_route.frequency:       240000,
+  tavg24_1d_route.ref_time:        000000,
+  tavg24_1d_route.fields:       'QSFLOW'      , 'ROUTE' ,
+                                'QOUTFLOW'    , 'ROUTE' ,
+                                'QRES'        , 'ROUTE' ,
+                                'QCAT'        , 'ROUTE' ,
+                                'QINFLOW'     , 'ROUTE' ,
+                                'QRUNOFF_CAT' , 'ROUTE' ,
+                                 ::                                 
 # ========================== EOF ==============================================================