Add ensemble forcing functionality

src/Applications/LDAS_App/GEOSldas_LDAS.rc

@@ -116,6 +116,18 @@ 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.
+#        - MET_PATH and MET_TAG should not contain the 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 
@@ -219,3 +231,8 @@ MAPL_ENABLE_BOOTSTRAP: YES
 #----- Write restart or checkpoint by oserver
 #
 WRITE_RESTART_BY_OSERVER: NO
+
+#
+# =================================== EOF ==========================================
+
+

src/Applications/LDAS_App/ldas_setup

@@ -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']
@@ -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']))

src/Applications/LDAS_App/lenkf.j.template

@@ -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`
@@ -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`

src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90

@@ -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
 
@@ -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
@@ -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
@@ -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)
@@ -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
@@ -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                                                            &
             )
@@ -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',            &
@@ -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)
@@ -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
 

src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/GEOS_MetforceGridComp.F90

@@ -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)
@@ -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(                                         &