Support for CatchCN ensemble simulations (#584)

Added CatchCN-specific export variables to EnsGridComp. 
See also GCM GC PR GEOS-ESM/GEOSgcm_GridComp#645

components.yaml

@@ -29,12 +29,12 @@ GMAO_Shared:
 MAPL:
   local: ./src/Shared/@MAPL
   remote: ../MAPL.git
-  tag: v2.25.0
+  tag: v2.27.1
   develop: develop
 
 GEOSgcm_GridComp:
   local: ./src/Components/GEOSldas_GridComp/@GEOSgcm_GridComp
   remote: ../GEOSgcm_GridComp.git
-  branch: develop  
+  branch: develop
   sparse: ./config/GEOSgcm_GridComp_ldas.sparse
   develop: develop

src/Applications/LDAS_App/GEOSldas_HIST.rc

@@ -179,7 +179,7 @@ COLLECTIONS:
 >>>HIST_CATCHCN<<<        'CNTOTC'   , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNVEGC'   , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNROOT'   , 'GridComp'  ,
->>>HIST_CATCHCNCLM45<<<   'CNFROOTC' , 'CATCHCNCLM45'  ,
+>>>HIST_CATCHCNCLM45<<<   'CNFROOTC' , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNNPP'    , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNGPP'    , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNSR'     , 'GridComp'  ,
@@ -281,7 +281,7 @@ COLLECTIONS:
 >>>HIST_CATCHCN<<<        'CNTOTC'   , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNVEGC'   , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNROOT'   , 'GridComp'  ,
->>>HIST_CATCHCNCLM45<<<   'CNFROOTC' , 'CATCHCNCLM45'  ,
+>>>HIST_CATCHCNCLM45<<<   'CNFROOTC' , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNNPP'    , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNGPP'    , 'GridComp'  ,
 >>>HIST_CATCHCN<<<        'CNSR'     , 'GridComp'  ,

src/Applications/LDAS_App/ldas_setup

@@ -113,6 +113,7 @@ class LDASsetup:
         self.ladas_coupling = 0
         self.in_rstfile = None
         self.in_tilefile = 'None' # default string
+        self.ens_id_width = 6     # _eXXXX
         # ------
         # Read exe input file which is required to set up the dir
         # ------
@@ -159,7 +160,10 @@ class LDASsetup:
         if self.nens > 1:
            self.perturb = 1
         self.ensdirs  = ['ens%04d'%iens for iens in range(self.first_ens_id, self.nens  + self.first_ens_id)]
-        self.ensids   = ['%04d'%iens    for iens in range(self.first_ens_id, self.nens  + self.first_ens_id)]
+        # if self.ens_id_width = 4, _width = '_e%04d'
+        _width = '_e%0{}d'.format(self.ens_id_width-2)
+        # self.ensids will be a list of [_e0000, _e0001, ...] 
+        self.ensids   = [ _width%iens for iens in range(self.first_ens_id, self.nens  + self.first_ens_id)]
         if (self.nens == 1) :
             self.ensdirs_avg = self.ensdirs
             self.ensids=['']
@@ -600,12 +604,12 @@ class LDASsetup:
         self._mkdir_p(self.scratchdir)
 
         #-start-shorthand-function-
-        def _getDirName(outtyp, ensid, yyyymm):
+        def _getDirName(outtyp, ensdir, yyyymm):
             return '/'.join([
                     self.outdir,
                     self.rqdExeInp['EXP_DOMAIN'],
                     outtyp, # ana/cat/rs/rc_out
-                    ensid,
+                    ensdir,
                     yyyymm
                     ])
         #-end-shorthand-function-
@@ -628,15 +632,15 @@ class LDASsetup:
                 y4m2_list.append('Y%4d/M%02d' % (newDate.year, newDate.month))
 
             # ExpDomain/ana/, /cat/ directories
-            for ensid in self.ensdirs_avg:
+            for ensdir in self.ensdirs_avg:
                 for y4m2 in y4m2_list:
-                    self._mkdir_p(_getDirName('ana', ensid, y4m2))
-                    self._mkdir_p(_getDirName('cat', ensid, y4m2))
+                    self._mkdir_p(_getDirName('ana', ensdir, y4m2))
+                    self._mkdir_p(_getDirName('cat', ensdir, y4m2))
 
             # ExpDomain/rs/ directories
-            for ensid in self.ensdirs:
+            for ensdir in self.ensdirs:
                 for y4m2 in y4m2_list:
-                    self._mkdir_p(_getDirName('rs', ensid, y4m2))
+                    self._mkdir_p(_getDirName('rs', ensdir, y4m2))
 
             # ExpDomain/rc_out/ - only for _start
             self._mkdir_p(_getDirName('rc_out', '', y4m2_list[0]))
@@ -872,8 +876,8 @@ class LDASsetup:
         catchRstFile0  = ''
         vegdynRstFile0 = ''
         for iens in range(self.nens) :
-            ens   = self.ensdirs[iens]
-            ensid = self.ensids[iens]
+            ensdir   = self.ensdirs[iens]
+            ensid    = self.ensids[iens]
             myCatchRst = myRstDir+'/'+self.catch +ensid +'_internal_rst'
             myVegRst   = myRstDir+'/'+'vegdyn'+ensid +'_internal_rst'
             myPertRst  = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst'
@@ -889,8 +893,8 @@ class LDASsetup:
                     vegdynRstFile = glob.glob(self.rqdExeInp['BCS_PATH']+'vegdyn_*.dat')[0]
                     catchRstFile  = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0]
                 else : # RESTART == 1
-                    catchRstFile = rstpath+ens +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
-                    vegdynRstFile= rstpath+ens +'/'+self.rqdExeInp['RESTART_ID']+ '.vegdyn_internal_rst'
+                    catchRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
+                    vegdynRstFile= rstpath+ensdir +'/'+self.rqdExeInp['RESTART_ID']+ '.vegdyn_internal_rst'
                     if not os.path.isfile(vegdynRstFile): # no vegdyn restart from LDASsa
                        if not os.path.isfile(vegdynRstFile0):
                          vegdynRstFile = glob.glob(self.rqdExeInp['BCS_PATH']+'vegdyn_*.dat')[0]
@@ -900,7 +904,7 @@ class LDASsetup:
 
             # catchment restart file
             if os.path.isfile(catchRstFile) :
-                catchLocal = self.rstdir+ens +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
+                catchLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
                 if self.islocal :
                     print( "Creating local catchment restart file... \n")
                     cmd='./preprocess_ldas.x c_localcatchrst '+ catchRstFile +' ' + catchLocal 
@@ -911,14 +915,14 @@ class LDASsetup:
 
                 catchRstFile = catchLocal
 
-                if '0000' in ens :
+                if '0000' in ensdir :
                     catchRstFile0 = catchRstFile
             else : # re-use 0000 catch file    
                 catchRstFile = catchRstFile0
 
             # vegdyn restart file
             if os.path.isfile(vegdynRstFile) :
-                vegdynLocal = self.rstdir+ens +'/'+self.rqdExeInp['EXP_ID']+'.vegdyn_internal_rst'
+                vegdynLocal = self.rstdir+ensdir +'/'+self.rqdExeInp['EXP_ID']+'.vegdyn_internal_rst'
                 if self.islocal :
                     print ("Creating the local veg restart file... \n")
                     cmd='./preprocess_ldas.x c_localvegrst '+ vegdynRstFile +' ' + vegdynLocal  
@@ -929,14 +933,14 @@ class LDASsetup:
 
                 vegdynRstFile = vegdynLocal
 
-                if '0000' in ens :
+                if '0000' in ensdir :
                     vegdynRstFile0 = vegdynRstFile
             else :    
                 vegdynRstFile = vegdynRstFile0
 
             if (self.has_geos_pert and self.perturb == 1) :
-                pertRstFile = rstpath+ens +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
-                pertLocal   = self.rstdir+ens +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
+                pertRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
+                pertLocal   = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
                 shutil.copy(pertRstFile,pertLocal)
                 pertRstFile = pertLocal
 
@@ -1160,7 +1164,7 @@ class LDASsetup:
 
                 if self.nens > 1 :
                    keyn='ENS_ID_WIDTH'
-                   valn='4'
+                   valn=str(self.ens_id_width)
                    ldasrcInp[keyn]= valn
 
                 if self.has_landassim_seed and self.assim :

src/Applications/LDAS_App/lenkf.j.template

@@ -674,6 +674,7 @@ EOF
           set ENSDIR = `echo ens${inens}`
        endif
        set ENSID = `echo $ENSDIR | cut -c4-7`
+       set ENSID = _e${ENSID}
        if ( $NENS == 1) set ENSID =''
        set THISDIR = $EXPDIR/output/$EXPDOMAIN/rs/$ENSDIR/Y${eYEAR}/M${eMON}/
        if (! -e $THISDIR            ) mkdir -p $THISDIR

src/Applications/LDAS_App/process_hist.csh

@@ -53,7 +53,13 @@ endif
 
 if($NENS > 1) then
    set GridComp = ENSAVG
-   sed -i 's|VEGDYN|'VEGDYN0000'|g' $HISTRC
+   sed -i 's|VEGDYN|'VEGDYN_e0000'|g' $HISTRC
+   sed -i 's|TP1|'TSOIL1TILE'|g' $HISTRC
+   sed -i 's|TP2|'TSOIL2TILE'|g' $HISTRC
+   sed -i 's|TP3|'TSOIL3TILE'|g' $HISTRC
+   sed -i 's|TP4|'TSOIL4TILE'|g' $HISTRC
+   sed -i 's|TP5|'TSOIL5TILE'|g' $HISTRC
+   sed -i 's|TP6|'TSOIL6TILE'|g' $HISTRC
 #   sed -i 's|DATAATM|'DATAATM0000'|g' $HISTRC
 endif
 

src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90

@@ -26,7 +26,7 @@ module GEOS_LdasGridCompMod
   use LDAS_ensdrv_mpi, only: MPI_tile_coord_type, MPI_grid_def_type
   use LDAS_ensdrv_mpi, only: init_MPI_types,mpicomm,numprocs,myid 
   use LDAS_ensdrv_mpi, only: root_proc
-  use LDAS_ensdrv_Globals, only: logunit,logit,root_logit,echo_clsm_ensdrv_glob_param
+  use LDAS_ensdrv_Globals, only: logunit,logit,root_logit,echo_clsm_ensdrv_glob_param, get_ensid_string
   use catch_constants, only: echo_catch_constants  
   use StieglitzSnow, only: StieglitzSnow_echo_constants
   use SurfParams,    only: SurfParams_init
@@ -78,15 +78,15 @@ subroutine SetServices(gc, rc)
     ! ensemble set up:
 
     integer :: i, k
-    integer,allocatable :: ens_id(:)
+    integer :: ens_id
     type(MAPL_MetaComp), pointer :: MAPL=>null()
     type(ESMF_GridComp), pointer :: gcs(:)=>null() ! Children gridcomps
     character(len=ESMF_MAXSTR), pointer :: gcnames(:)=>null() ! Children's names
     ! ErrLog variables
     integer :: status
     character(len=ESMF_MAXSTR) :: Iam
     character(len=ESMF_MAXSTR) :: comp_name
-    character(len=ESMF_MAXSTR) :: id_string,childname, fmt_str
+    character(len=ESMF_MAXSTR) :: ensid_string,childname
     character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file, ENS_FORCING_STR
     integer                    :: ens_id_width
     ! Local variables
@@ -173,42 +173,45 @@ subroutine SetServices(gc, rc)
        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,")" 
+    allocate(LAND(NUM_ENSEMBLE),LANDPERT(NUM_ENSEMBLE))
+
+    ! 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
+    !   hardwired when assembling a format string.
+    ! Assert ens_id_width>=2+3 to avoid user configuration errors when LDAS is coupled into ADAS.
+    !   (Met forcing from the atm ensemble uses hardwired, 3-character ensemble IDs.) 
+
+    if (NUM_ENSEMBLE > 1) then
+       _ASSERT( ens_id_width < 12, "Must use ens_id_width <= 11  (2 for '_e' + 9 digits max)")
+       _ASSERT( ens_id_width >= 5, "Must use ens_id_width >=  5  (2 for '_e' + 3 digits min)")
+    endif
 
     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
+       ens_id = i-1 + FIRST_ENS_ID ! id start form FIRST_ENS_ID
+
+       call get_ensid_string(ensid_string, ens_id, ens_id_width, NUM_ENSEMBLE)   ! "_eXXXX"
 
-       id_string=trim(id_string)
+       ! allow for Catchment ensemble simulation to be forced with single-member met inputs
+       if (.not. ensemble_forcing ) ensid_string = ''    
 
-       childname='METFORCE'//trim(id_string)
+       childname='METFORCE'//trim(ensid_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
-          id_string=''
-       else
-          write(id_string, fmt_str) ens_id(i)
-       endif
+       ens_id = i-1 + FIRST_ENS_ID ! id start form FIRST_ENS_ID
 
-       id_string=trim(id_string)
+       call get_ensid_string(ensid_string, ens_id, ens_id_width, NUM_ENSEMBLE)
 
-       childname='LANDPERT'//trim(id_string)
+       childname='LANDPERT'//trim(ensid_string)
        LANDPERT(i) = MAPL_AddChild(gc, name=childname, ss=LandPertSetServices, rc=status)
        VERIFY_(status)
 
-       childname='LAND'//trim(id_string)
+       childname='LAND'//trim(ensid_string)
        LAND(i) = MAPL_AddChild(gc, name=childname, ss=LandSetServices, rc=status)
        VERIFY_(status)
     enddo