diff --git a/CHANGELOG.md b/CHANGELOG.md index a6f7f0d..d49e78b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added support for river routing. + ### Changed ### Fixed diff --git a/CMakeLists.txt b/CMakeLists.txt index 2de54c4..0bd9dee 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/GEOS_LdasGridComp.F90 b/GEOS_LdasGridComp.F90 index ccb6ad8..7250c04 100644 --- a/GEOS_LdasGridComp.F90 +++ b/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 ! 1/2 if run river-routine grid comp without/with 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,14 @@ 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 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) diff --git a/GEOSldas_App/GEOSldas_HIST.rc b/GEOSldas_App/GEOSldas_HIST.rc index 03eebb0..3a66953 100644 --- a/GEOSldas_App/GEOSldas_HIST.rc +++ b/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,17 @@ 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' , + :: # ========================== EOF ============================================================== diff --git a/GEOSldas_App/ldas.py b/GEOSldas_App/ldas.py index bdab4f6..babaa3d 100644 --- a/GEOSldas_App/ldas.py +++ b/GEOSldas_App/ldas.py @@ -32,10 +32,11 @@ def __init__(self, cmdLineArgs): """ """ # These keywords are excluded from LDAS.rc (i.e., only needed in pre- or post-processing) + # or their values are re-assigned in LDAS.rc self.NoneLDASrcKeys=['EXP_ID', 'EXP_DOMAIN', 'BEG_DATE', 'END_DATE','RESTART','RESTART_PATH', 'RESTART_DOMAIN','RESTART_ID','BCS_PATH','TILING_FILE','GRN_FILE','LAI_FILE','LNFM_FILE','NIRDF_FILE', - 'VISDF_FILE','CATCH_DEF_FILE','NDVI_FILE', + 'VISDF_FILE','CATCH_DEF_FILE','NDVI_FILE', 'TILE2PFAF_FILE', 'NML_INPUT_PATH','HISTRC_FILE','RST_FROM_GLOBAL','JOB_SGMT','NUM_SGMT','POSTPROC_HIST', 'MINLON','MAXLON','MINLAT','MAXLAT','EXCLUDE_FILE','INCLUDE_FILE','MWRTM_PATH','GRIDNAME', 'ADAS_EXPDIR', 'BCS_RESOLUTION', 'TILE_FILE_FORMAT' ] @@ -101,6 +102,7 @@ def __init__(self, cmdLineArgs): self.tile_types = '' self.with_land = False self.with_landice = False + self.run_route = 0 self.adas_expdir = '' # assert necessary optional arguments in command line if exeinp file does not exsit @@ -235,6 +237,8 @@ def __init__(self, cmdLineArgs): if "20" in self.tile_types : self.with_landice = True + self.run_route = int(self.ExeInputs.get('RUN_ROUTE',0)) + self.nens = int(self.ExeInputs['NUM_LDAS_ENSEMBLE']) # fails if value of Nens is not an integer self.first_ens_id = int(self.ExeInputs.get('FIRST_ENS_ID',0)) self.perturb = int(self.ExeInputs.get('PERTURBATIONS',0)) @@ -274,6 +278,7 @@ def __init__(self, cmdLineArgs): # assigning BC files self.ExeInputs['LNFM_FILE'] = '' + self.ExeInputs['TILE2PFAF_FILE'] = '' tile_file_format = self.ExeInputs.get('TILE_FILE_FORMAT', 'DEFAULT') domain_ = '' inpdir_ = self.bcs_dir_land @@ -305,28 +310,36 @@ def __init__(self, cmdLineArgs): tmp_ = glob.glob(inpdir_ + 'lnfm_clim_*.data'+domain_) if (len(tmp_) ==1) : self.ExeInputs['LNFM_FILE'] = tmp_[0] + self.ExeInputs['NDVI_FILE'] = glob.glob(inpdir_ + 'ndvi_clim_*.data'+domain_ )[0] self.ExeInputs['NIRDF_FILE'] = glob.glob(inpdir_ + 'nirdf_*.dat' +domain_ )[0] self.ExeInputs['VISDF_FILE'] = glob.glob(inpdir_ + 'visdf_*.dat' +domain_ )[0] - inpdir_ = None - domain_ = None - inpgeom_= None + # assigning Gridname if 'GRIDNAME' not in self.ExeInputs : - tmptile = os.path.realpath(self.ExeInputs['TILING_FILE']) - extension = os.path.splitext(tmptile)[1] - if extension == '.domain': - extension = os.path.splitext(tmptile)[0] - gridname_ ='' - if extension == '.til': - gridname_ = linecache.getline(tmptile, 3).strip() - else: - nc_file = netCDF4.Dataset(tmptile,'r') - gridname_ = nc_file.getncattr('Grid_Name') - # in case it is an old name: SMAP-EASEvx-Mxx - gridname_ = gridname_.replace('SMAP-','').replace('-M','_M') - self.ExeInputs['GRIDNAME'] = gridname_ + tmptile = os.path.realpath(self.ExeInputs['TILING_FILE']) + extension = os.path.splitext(tmptile)[1] + if extension == '.domain': + extension = os.path.splitext(tmptile)[0] + gridname_ ='' + if extension == '.til': + gridname_ = linecache.getline(tmptile, 3).strip() + else: + nc_file = netCDF4.Dataset(tmptile,'r') + gridname_ = nc_file.getncattr('Grid_Name') + # in case it is an old name: SMAP-EASEvx-Mxx + gridname_ = gridname_.replace('SMAP-','').replace('-M','_M') + self.ExeInputs['GRIDNAME'] = gridname_ + if (self.run_route > 0 and 'EASE' in self.ExeInputs['GRIDNAME']): + tmp_= glob.glob(inpgeom_ + '*tile2pfaf*') + if (len(tmp_) > 0) : + self.ExecInput['TILE2PFAF_FILE'] = tmp_[0] + + inpdir_ = None + domain_ = None + inpgeom_= None + if 'POSTPROC_HIST' not in self.ExeInputs: self.ExeInputs['POSTPROC_HIST'] = 0 @@ -395,6 +408,11 @@ def __init__(self, cmdLineArgs): landiceRstFile=self.in_rstdir+'/'+tmpFile assert os.path.isfile(landiceRstFile), 'landice_internal_rst file [%s] does not exist!' %(landiceRstFile) + if self.run_route > 0: + tmpFile=self.ExeInputs['RESTART_ID']+'.route_internal_rst.'+y4m2d2_h2m2 + routeRstFile=self.in_rstdir+'/'+tmpFile + assert os.path.isfile(routeRstFile), 'route_internal_rst file [%s] does not exist!' %(routeRstFile) + if RESTART_str == '0': assert ( self.with_land and not self.with_landice), "RESTART = 0 is only for land" if (self.catch == 'catch'): @@ -760,6 +778,9 @@ def createLnRstBc(self) : bcs += [self.ExeInputs['LNFM_FILE']] if (self.has_vegopacity): bcs += [self.ExeInputs['VEGOPACITY_FILE']] + if (self.ExeInputs['TILE2PFAF_FILE'] != ''): + bcs += [self.ExeInputs['TILE2PFAF_FILE']] + bcstmp=[] for bcf in bcs : shutil.copy(bcf, self.bcsdir+'/') @@ -784,8 +805,12 @@ def createLnRstBc(self) : bcnames += ['lnfm'] if (self.has_vegopacity): bcnames += ['vegopacity'] + if (self.ExeInputs['TILE2PFAF_FILE'] != ''): + bcnames += ['tile2pfaf.nc4'] for bcln,bc in zip(bcnames,bcs) : myBC=self.inpdir+'/'+bcln+'.data' + if '.nc4' in bcln: + myBC=self.inpdir+'/'+bcln os.symlink(bc,myBC) if ("catchcn" in self.catch): @@ -903,10 +928,11 @@ def createLnRstBc(self) : for iens in range(self.nens) : ensdir = self.ensdirs[iens] ensid = self.ensids[iens] - myCatchRst = myRstDir+'/'+self.catch +ensid +'_internal_rst' + myCatchRst = myRstDir+'/'+self.catch +ensid +'_internal_rst' myLandiceRst = myRstDir+'/'+ 'landice' +ensid +'_internal_rst' - myVegRst = myRstDir+'/'+'vegdyn'+ensid +'_internal_rst' - myPertRst = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst' + myVegRst = myRstDir+'/'+ 'vegdyn'+ensid +'_internal_rst' + myPertRst = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst' + myRouteRst = myRstDir+'/'+ 'route' +ensid +'_internal_rst' catchRstFile = '' vegdynRstFile = '' @@ -991,6 +1017,27 @@ def createLnRstBc(self) : landiceRstFile0 = landiceRstFile else : landiceRstFile = landiceRstFile0 + + routeRstFile = '' + if self.run_route > 0 : + if RESTART_str == '0' : + routeRstFile = "/discover/nobackup/yzeng3/data/river_input_weiyuan/route_restart_package/"+"route_internal_rst."+y4m2[-2:]+"01_0000" + if RESTART_str in ['1', '2'] : + routeRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.ExeInputs['RESTART_ID']+'.'+'route_internal_rst.'+y4m2d2_h2m2 + if RESTART_str == 'M': + exit(" RUN_ROUTE does not support MERRA 2 option") + + if os.path.isfile(routeRstFile) : + routeLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.ExeInputs['EXP_ID']+'.route_internal_rst.'+y4m2d2_h2m2 + shutil.copy(routeRstFile,routeLocal) + # WY note: after the copy, depending on in_bc and out_bc version, + # the routeLocal can be changed here + routeRstFile = routeLocal + + if '0000' in ensdir : + routeRstFile0 = routeRstFile + else : + routeRstFile = routeRstFile0 if (self.has_geos_pert and self.perturb == 1) : pertRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.ExeInputs['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2 @@ -1003,9 +1050,15 @@ def createLnRstBc(self) : print ('vegdynRstFile: ' + vegdynRstFile) os.symlink(catchRstFile, myCatchRst) os.symlink(vegdynRstFile, myVegRst) + if self.with_landice : print("link landice restart: " + myLandiceRst) os.symlink(landiceRstFile, myLandiceRst) + + if self.run_route > 0 : + print("link route restart: " + myRouteRst) + os.symlink(routeRstFile, myRouteRst) + if ( self.has_geos_pert and self.perturb == 1 ): os.symlink(pertRstFile, myPertRst) @@ -1219,10 +1272,16 @@ def createRCFiles(self): if self.with_land : bcval=['../input/green','../input/lai','../input/lnfm','../input/ndvi','../input/nirdf','../input/visdf'] bckey=['GREEN','LAI','LNFM','NDVI','NIRDF','VISDF'] + if self.ExeInputs['TILE2PFAF_FILE'] !='': + bcval.append('../input/tile2pfaf.nc4') + bckey.append('TILE2PFAF') for key, val in zip(bckey,bcval): keyn = key+'_FILE' valn = val+'.data' + if '.nc4' in val: + valn = val ldasrcInp[keyn]= valn + if('catchcn' in self.catch): ldasrcInp['CO2_MonthlyMean_DiurnalCycle_FILE']= '../input/CO2_MonthlyMean_DiurnalCycle.nc4' else: @@ -1255,6 +1314,10 @@ def createRCFiles(self): rstkey.append('LANDICE') rstval.append('landice') + if self.run_route > 0: + rstkey.append('ROUTE') + rstval.append('route') + if self.has_mwrtm : keyn='LANDASSIM_INTERNAL_RESTART_FILE' valn='../input/restart/mwrtm_param_rst' @@ -1282,18 +1345,12 @@ def createRCFiles(self): for key,val in zip(rstkey,rstval) : keyn = key+ '_INTERNAL_RESTART_FILE' valn = '../input/restart/'+val+tmpl_+'_internal_rst' - ldasrcInp[keyn]= valn + ldasrcInp[keyn] = valn + if 'VEGDYN' not in key: # not checkpoi for vegdyn + keyn = key + '_INTERNAL_CHECKPOINT_FILE' + valn = val + tmpl_ + '_internal_checkpoint' + ldasrcInp[keyn] = valn - # checkpoint file and its type - if self.with_land : - keyn = catch_ + '_INTERNAL_CHECKPOINT_FILE' - valn = self.catch+tmpl_+'_internal_checkpoint' - ldasrcInp[keyn]= valn - - if self.with_landice : - keyn = 'LANDICE_INTERNAL_CHECKPOINT_FILE' - valn = 'landice'+tmpl_+'_internal_checkpoint' - ldasrcInp[keyn]= valn # specify LANDPERT restart file if (self.perturb == 1): keyn = 'LANDPERT_INTERNAL_RESTART_FILE' diff --git a/GEOSldas_App/lenkf_j_template.py b/GEOSldas_App/lenkf_j_template.py index 60443ff..e9f681e 100644 --- a/GEOSldas_App/lenkf_j_template.py +++ b/GEOSldas_App/lenkf_j_template.py @@ -721,7 +721,7 @@ set THISDIR = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/ if (! -e $THISDIR ) mkdir -p $THISDIR - set rstfs = (${{LANDMODEL}} 'landice') + set rstfs = (${{LANDMODEL}} 'landice' 'route') foreach rstf ( $rstfs ) if (-f ${{rstf}}${{ENSID}}_internal_checkpoint ) then set tmp_file = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${{eYEAR}}/M${{eMON}}/${{EXPID}}.${{rstf}}_internal_rst.${{eYEAR}}${{eMON}}${{eDAY}}_${{eHour}}${{eMin}}