Skip to content

Commit

Permalink
Merge pull request #453 from GEOS-ESM/feature/wjiang/ensemble_forcing
Browse files Browse the repository at this point in the history 
Add ensemble forcing functionality
  • Loading branch information
@gmao-rreichle
gmao-rreichle authored Feb 16, 2022
2 parents 63b3464 + 45411d2 commit 1fe076e
Show file tree
Hide file tree
Showing 5 changed files with 120 additions and 31 deletions.
19 changes: 19 additions & 0 deletions src/Applications/LDAS_App/GEOSldas_LDAS.rc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,20 @@ PERTURBATIONS: 0
FIRST_ENS_ID: 0


# ---- Ensemble forcing
#
# NO : Deterministic met forcing (default)
# YES : Ensemble met forcing
# - Typically used in land-atmosphere data assimilation (LADAS) configuration.
# - Must have forcing with matching ensemble ID for each land ensemble member.
# - Forcing files must be stored in member-specific directories MET_PATH[NNN]/,
# where NNN is the 3-digit ensemble ID.
# - User-specified MET_PATH and MET_TAG must not contain ensemble IDs.
# - FIRST_ENS_ID may be used to align ensemble IDs.
#
ENSEMBLE_FORCING: NO


# ---- Path to special namelist input files
#
# Applies only for ensemble simulations. The variable values in the special
Expand Down Expand Up @@ -219,3 +233,8 @@ MAPL_ENABLE_BOOTSTRAP: YES
#----- Write restart or checkpoint by oserver
#
WRITE_RESTART_BY_OSERVER: NO

#
# =================================== EOF ==========================================


28 changes: 21 additions & 7 deletions src/Applications/LDAS_App/ldas_setup
View file
Original file line number Diff line number Diff line change
Expand Up @@ -628,13 +628,25 @@ class LDASsetup:
os.symlink(self.blddir, self.blddirLn)

# met forcing dir
if 'MET_PATH' in self.rqdExeInp:
metpath = self.rqdExeInp['MET_PATH'].rstrip('/')
myMetDir = self.inpdir + '/met_forcing'
myMetPath = myMetDir + '/' + metpath.split('/')[-1]
os.symlink(metpath, myMetPath)
# update 'met_path' to use relative path from outdir
self.rqdExeInp['MET_PATH'] = os.path.relpath(myMetPath, self.rundir)
self.ensemble_forcing = True if self.rqdExeInp.get('ENSEMBLE_FORCING', 'NO').upper() == 'YES' else False

myMetPath =''
for _i in range(self.first_ens_id, _nens + self.first_ens_id) :
str_ens = ''
if ( _nens != 1 and self.ensemble_forcing):
str_ens = '%03d'%(_i)
metpath = self.rqdExeInp['MET_PATH'].rstrip('/')+str_ens
myMetDir = self.inpdir + '/met_forcing'
myMetPath = myMetDir + '/' + metpath.split('/')[-1]
os.symlink(metpath, myMetPath)
# update 'met_path' to use relative path from outdir
if ( not self.ensemble_forcing):
break
if ( _nens !=1 and self.ensemble_forcing) :
# replace last three character with '%s"
self.rqdExeInp['MET_PATH'] = os.path.relpath(myMetPath, self.rundir)[:-3]+'%s'
else:
self.rqdExeInp['MET_PATH'] = os.path.relpath(myMetPath, self.rundir)

# update tile file
tile= self.rqdExeInp['TILING_FILE']
Expand Down Expand Up @@ -1289,6 +1301,8 @@ class LDASsetup:
fout.write(line.replace('MY_FIRST_ENS_ID',str(self.first_ens_id)))
elif 'MY_LADAS_COUPLING' in line :
fout.write(line.replace('MY_LADAS_COUPLING',str(self.ladas_coupling)))
elif 'MY_ENSEMBLE_FORCING' in line :
fout.write(line.replace('MY_ENSEMBLE_FORCING',self.rqdExeInp.get('ENSEMBLE_FORCING', 'NO').upper()))
elif 'MY_ADAS_EXPDIR' in line :
if self.ladas_coupling > 0:
fout.write(line.replace('MY_ADAS_EXPDIR', self.rqdExeInp['ADAS_EXPDIR']))
Expand Down
3 changes: 2 additions & 1 deletion src/Applications/LDAS_App/lenkf.j.template
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ setenv POSTPROC_HIST MY_POSTPROC_HIST
# : 2 -- LDAS coupled to atmospheric ensemble component of ADAS

setenv LADAS_COUPLING MY_LADAS_COUPLING
setenv ENSEMBLE_FORCING MY_ENSEMBLE_FORCING
setenv ADAS_EXPDIR MY_ADAS_EXPDIR

set NENS = `grep NUM_LDAS_ENSEMBLE: $HOMDIR/LDAS.rc | cut -d':' -f2`
Expand All @@ -77,7 +78,7 @@ set NUM_SGMT = `grep NUM_SGMT: $HOMDIR/CAP.rc | cut -d':' -f2`
# if LADAS_COUPLING==2, compute ens avg of atmens forcing
#######################################################################

if ( $LADAS_COUPLING == 2 ) then
if ( $LADAS_COUPLING == 2 && $ENSEMBLE_FORCING == "NO" ) then
cd $HOMDIR
set force_in = $ADAS_EXPDIR
set force_out = `grep MET_PATH: $HOMDIR/LDAS.rc | cut -d ':' -f2`
Expand Down
77 changes: 58 additions & 19 deletions src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,12 +51,14 @@ module GEOS_LdasGridCompMod
! All children
integer,allocatable :: LAND(:)
integer,allocatable :: LANDPERT(:)
integer :: METFORCE, ENSAVG, LANDASSIM
integer,allocatable :: METFORCE(:)
integer :: ENSAVG, LANDASSIM

! other global variables
integer :: NUM_ENSEMBLE
integer :: NUM_ENSEMBLE ! number of land ensemble members
logical :: land_assim
logical :: mwRTM
logical :: ensemble_forcing ! switch between deterministic and ensemble forcing

contains

Expand All @@ -75,7 +77,7 @@ subroutine SetServices(gc, rc)

! ensemble set up:

integer :: i
integer :: i, k
integer,allocatable :: ens_id(:)
type(MAPL_MetaComp), pointer :: MAPL=>null()
type(ESMF_GridComp), pointer :: gcs(:)=>null() ! Children gridcomps
Expand All @@ -85,7 +87,7 @@ subroutine SetServices(gc, rc)
character(len=ESMF_MAXSTR) :: Iam
character(len=ESMF_MAXSTR) :: comp_name
character(len=ESMF_MAXSTR) :: id_string,childname, fmt_str
character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file
character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file, ENS_FORCING_STR
integer :: ens_id_width
! Local variables
type(T_TILECOORD_STATE), pointer :: tcinternal
Expand Down Expand Up @@ -144,6 +146,11 @@ subroutine SetServices(gc, rc)
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)
VERIFY_(STATUS)
ENS_FORCING_STR = ESMF_UtilStringUpperCase(ENS_FORCING_STR, rc=STATUS)
VERIFY_(STATUS)
ensemble_forcing = (trim(ENS_FORCING_STR) == 'YES')

call MAPL_GetResource ( MAPL, LAND_ASSIM_STR, Label="LAND_ASSIM:", DEFAULT="NO", RC=STATUS)
VERIFY_(STATUS)
Expand All @@ -160,12 +167,33 @@ subroutine SetServices(gc, rc)
_ASSERT( .not. (mwRTM .or. land_assim), "CatchCN is Not Ready for assimilation or mwRTM")
endif

METFORCE = MAPL_AddChild(gc, name='METFORCE', ss=MetforceSetServices, rc=status)
VERIFY_(status)
if (ensemble_forcing) then
allocate(METFORCE(NUM_ENSEMBLE))
else
allocate(METFORCE(1))
endif

allocate(ens_id(NUM_ENSEMBLE),LAND(NUM_ENSEMBLE),LANDPERT(NUM_ENSEMBLE))
_ASSERT( ens_id_width < 10, "need 1 billion ensemble members? increase ens_id_width first")
write (fmt_str, "(A2,I1,A1,I1,A1)") "(I", ens_id_width,".",ens_id_width,")"

do i=1,NUM_ENSEMBLE
ens_id(i) = i-1 + FIRST_ENS_ID ! id start form FIRST_ENS_ID
if(NUM_ENSEMBLE == 1 .or. .not. ensemble_forcing) then
id_string=''
else
write(id_string, fmt_str) ens_id(i)
endif

id_string=trim(id_string)

childname='METFORCE'//trim(id_string)
METFORCE(i) = MAPL_AddChild(gc, name=trim(childname), ss=MetforceSetServices, rc=status)
VERIFY_(status)
! exit after i=1 if using deterministic forcing
if (.not. ensemble_forcing ) exit
enddo

do i=1,NUM_ENSEMBLE
ens_id(i) = i-1 + FIRST_ENS_ID ! id start form FIRST_ENS_ID
if(NUM_ENSEMBLE == 1 ) then
Expand Down Expand Up @@ -196,12 +224,14 @@ subroutine SetServices(gc, rc)
! Connections
do i=1,NUM_ENSEMBLE
! -METFORCE-feeds-LANDPERT's-imports-
k = 1
if ( ensemble_forcing ) k = i
call MAPL_AddConnectivity( &
gc, &
SHORT_NAME = ['Tair ', 'Qair ', 'Psurf ', 'Rainf_C', 'Rainf ', &
'Snowf ', 'LWdown ', 'SWdown ', 'PARdrct', 'PARdffs', &
'Wind ', 'RefH '], &
SRC_ID = METFORCE, &
SRC_ID = METFORCE(k), &
DST_ID = LANDPERT(i), &
rc = status &
)
Expand Down Expand Up @@ -230,7 +260,7 @@ subroutine SetServices(gc, rc)
'DUDP ', 'DUSV ', 'DUWT ', 'DUSD ', 'BCDP ', 'BCSV ', &
'BCWT ', 'BCSD ', 'OCDP ', 'OCSV ', 'OCWT ', 'OCSD ', &
'SUDP ', 'SUSV ', 'SUWT ', 'SUSD ', 'SSDP ', 'SSSV ' ], &
SRC_ID = METFORCE, &
SRC_ID = METFORCE(k), &
DST_NAME = ['PS ', 'DZ ', &
'DUDP', 'DUSV', 'DUWT', 'DUSD', 'BCDP', 'BCSV', &
'BCWT', 'BCSD', 'OCDP', 'OCSV', 'OCWT', 'OCSD', &
Expand Down Expand Up @@ -661,17 +691,21 @@ subroutine Initialize(gc, import, export, clock, rc)
tcinternal%grid_f = tile_grid_f
tcinternal%grid_l = tile_grid_l

call MAPL_GetObjectFromGC(gcs(METFORCE), CHILD_MAPL, rc=status)
VERIFY_(status) ! CHILD = METFORCE
call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
VERIFY_(status)
do i = 1, NUM_ENSEMBLE
call MAPL_GetObjectFromGC(gcs(METFORCE(i)), CHILD_MAPL, rc=status)
VERIFY_(status) ! CHILD = METFORCE
call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
VERIFY_(status)
call ESMF_UserCompSetInternalState(gcs(METFORCE(i)), 'TILE_COORD', tcwrap, status)
VERIFY_(status)
! exit after i=1 if using deterministic forcing
if (.not. ensemble_forcing) exit
enddo

call MAPL_GetObjectFromGC(gcs(ENSAVG), CHILD_MAPL, rc=status)
VERIFY_(status) ! CHILD = ens_avg
call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
VERIFY_(status)
call ESMF_UserCompSetInternalState(gcs(METFORCE), 'TILE_COORD', tcwrap, status)
VERIFY_(status)

do i = 1,NUM_ENSEMBLE
call MAPL_GetObjectFromGC(gcs(LAND(i)), CHILD_MAPL, rc=status)
Expand Down Expand Up @@ -835,11 +869,16 @@ subroutine Run(gc, import, export, clock, rc)
enddo


igc = METFORCE
call MAPL_TimerOn(MAPL, gcnames(igc))
call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, userRC=status)
VERIFY_(status)
call MAPL_TimerOff(MAPL, gcnames(igc))
do i = 1, NUM_ENSEMBLE
igc = METFORCE(i)
call MAPL_TimerOn(MAPL, gcnames(igc))
call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, userRC=status)
VERIFY_(status)
call MAPL_TimerOff(MAPL, gcnames(igc))
! exit after i=1 if using deterministic forcing
if (.not. ensemble_forcing) exit
enddo


do i = 1,NUM_ENSEMBLE

Expand Down
24 changes: 20 additions & 4 deletions src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -559,15 +559,15 @@ subroutine Initialize(gc, import, export, clock, rc)
type(tile_coord_type), pointer :: tile_coord(:)=>null()

! Misc variables
integer :: land_nt_local
integer :: land_nt_local, k, NUM_ENSEMBLE, ens_id_width
integer :: ForceDtStep
type(met_force_type) :: mf_nodata
logical :: MERRA_file_specs
logical :: MERRA_file_specs, ensemble_forcing
logical :: backward_looking_fluxes

integer :: AEROSOL_DEPOSITION
type(MAPL_LocStream) :: locstream
character(len=ESMF_MAXSTR) :: grid_type
character(len=ESMF_MAXSTR) :: grid_type, ENS_FORCING_STR, ens_forcing_path, id_string
character(len=ESMF_MAXSTR) :: gridname
type(ESMF_Grid) :: agrid
integer :: dims(ESMF_MAXDIM)
Expand Down Expand Up @@ -669,9 +669,25 @@ subroutine Initialize(gc, import, export, clock, rc)
! -allocate-memory-for-avg-zenith-angle
allocate(mf%zenav(land_nt_local), source=nodata_generic, stat=status)
VERIFY_(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)
call MAPL_GetResource ( MAPL, NUM_ENSEMBLE, Label="NUM_LDAS_ENSEMBLE:", DEFAULT=1, RC=STATUS)
VERIFY_(STATUS)
ensemble_forcing = (trim(ENS_FORCING_STR) == 'YES')
if (ensemble_forcing .and. NUM_ENSEMBLE > 1) then
id_string = ""
call MAPL_GetResource ( MAPL, ens_id_width, Label="ENS_ID_WIDTH:", DEFAULT=0, RC=STATUS)
k = len(trim(comp_name))
id_string = comp_name(k-ens_id_width+1:k)
k = len(trim(id_string))
! hard coded 3 character for forcing
call ESMF_CFIOStrTemplate(ens_forcing_path, trim(adjustl(mf%Path)),'GRADS', xid = trim(id_string(k-2:k)), stat=status)
mf%Path = ens_forcing_path
endif
! Put MetForcing in Ldas' pvt internal state
internal%mf = mf

! Create alarm for MetForcing
! -create-nonsticky-alarm-
MetForcingAlarm = ESMF_AlarmCreate( &
Expand Down

0 comments on commit 1fe076e

Please sign in to comment.