Add an option to point to the "enkfgdas" backgrounds instead of the p…

…re-generated 30 member ensemble (#771)

#### Work done 
Mainly what the title says. The default ctest is now assuming that
`DOHYBVAR`=`YES`, I don't really think it's worth it to cover the
`DOHYBVAR`=`NO` in the ctest. I may change my mind, but in both cases,
`soca` is configured to use the hybrid envar solver, the only difference
is the origin of the ensemble members.

#### Work left to be done
We need to define an ensemble "COM_OCEAN/ICE_HISTORY_PREV" environment
variable in the prep jjob (in the g-w). Flagged as a `TODO` in the code.

parm/soca/berror/saber_blocks.yaml

@@ -60,7 +60,7 @@ components:
         ice_filename: 'ice.%mem%.nc'
         state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon]
       pattern: '%mem%'
-      nmembers: ${CLIM_ENS_SIZE}
+      nmembers: ${ENS_SIZE}
     localization:
       localization method: SABER
       saber central block:

parm/soca/berror/soca_ensb.yaml

@@ -17,7 +17,7 @@ vertical geometry:
   read_from_file: 3
 
 soca increments:
-  number of increments: ${CLIM_ENS_SIZE}
+  number of increments: ${ENS_SIZE}
   pattern: '%mem%'
   template:
     date: '{{ATM_WINDOW_BEGIN}}'

scripts/exgdas_global_marine_analysis_prep.py

@@ -89,6 +89,10 @@ def cice_hist2fms(input_filename, output_filename):
     # open the CICE history file
     ds = xr.open_dataset(input_filename)
 
+    if 'aicen' in ds.variables and 'hicen' in ds.variables and 'hsnon' in ds.variables:
+        logging.info(f"*** Already reformatted, skipping.")
+        return
+
     # rename the dimensions to xaxis_1 and yaxis_1
     ds = ds.rename({'ni': 'xaxis_1', 'nj': 'yaxis_1'})
 
@@ -224,6 +228,11 @@ def find_clim_ens(input_date):
 comin_obs = os.getenv('COMIN_OBS')  # R2D2 DB for now
 anl_dir = os.getenv('DATA')
 staticsoca_dir = os.getenv('SOCA_INPUT_FIX_DIR')
+if os.getenv('DOHYBVAR') == "YES":
+    dohybvar = True
+    nmem_ens = int(os.getenv('NMEM_ENS'))
+else:
+    dohybvar = False
 
 # create analysis directories
 diags = os.path.join(anl_dir, 'diags')            # output dir for soca DA obs space
@@ -242,6 +251,7 @@ def find_clim_ens(input_date):
 fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H')
 RUN = os.getenv('RUN')
 cyc = os.getenv('cyc')
+gcyc = os.getenv('gcyc')
 PDY = os.getenv('PDY')
 
 ################################################################################
@@ -310,21 +320,43 @@ def find_clim_ens(input_date):
 FileHandler({'copy': bkgerr_list}).sync()
 
 ################################################################################
-# stage static ensemble
-
-logging.info(f"---------------- Stage climatological ensemble")
-clim_ens_member_list = []
-clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None))
-clim_ens_size = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc'))))
-os.environ['CLIM_ENS_SIZE'] = str(clim_ens_size)
-
-for domain in ['ocn', 'ice']:
-    for mem in range(1, clim_ens_size+1):
-        fname = domain+"."+str(mem)+".nc"
-        fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname))
-        fname_out = os.path.abspath(os.path.join(static_ens, fname))
-        clim_ens_member_list.append([fname_in, fname_out])
-FileHandler({'copy': clim_ens_member_list}).sync()
+# stage ensemble members
+if dohybvar:
+    logging.info("---------------- Stage ensemble members")
+    nmem_ens = int(os.getenv('NMEM_ENS'))
+    longname = {'ocn': 'ocean', 'ice': 'ice'}
+    ens_member_list = []
+    for mem in range(1, nmem_ens+1):
+        for domain in ['ocn', 'ice']:
+            # TODO(Guillaume): make use and define ensemble COM in the j-job
+            ensdir = os.path.join(os.getenv('COM_OCEAN_HISTORY_PREV'), '..', '..', '..', '..', '..',
+                                  f'enkf{RUN}.{PDY}', f'{gcyc}', f'mem{str(mem).zfill(3)}',
+                                  'model_data', longname[domain], 'history')
+            ensdir = os.path.normpath(ensdir)
+            f009 = f'enkfgdas.t{gcyc}z.{domain}f009.nc'
+
+            fname_in = os.path.abspath(os.path.join(ensdir, f009))
+            fname_out = os.path.abspath(os.path.join(static_ens, domain+"."+str(mem)+".nc"))
+            ens_member_list.append([fname_in, fname_out])
+    FileHandler({'copy': ens_member_list}).sync()
+
+    # reformat the cice history output
+    for mem in range(1, nmem_ens+1):
+        cice_fname = os.path.abspath(os.path.join(static_ens, "ice."+str(mem)+".nc"))
+        cice_hist2fms(cice_fname, cice_fname)
+else:
+    logging.info("---------------- Stage offline ensemble members")
+    ens_member_list = []
+    clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None))
+    nmem_ens = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc'))))
+    for domain in ['ocn', 'ice']:
+        for mem in range(1, nmem_ens+1):
+            fname = domain+"."+str(mem)+".nc"
+            fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname))
+            fname_out = os.path.abspath(os.path.join(static_ens, fname))
+            ens_member_list.append([fname_in, fname_out])
+    FileHandler({'copy': ens_member_list}).sync()
+os.environ['ENS_SIZE'] = str(nmem_ens)
 
 ################################################################################
 # prepare JEDI yamls
@@ -456,7 +488,6 @@ def find_clim_ens(input_date):
 else:
     logging.info(f"using default SABER blocks yaml")
     os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'saber_blocks.yaml')
-os.environ['CLIM_ENS_SIZE'] = str(clim_ens_size)
 
 # substitute templated variables in the var config
 logging.info(f"{config}")

test/soca/gw/CMakeLists.txt

@@ -39,7 +39,6 @@ file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/test/soca/gw/apps_scratch)
 # Test JGDAS_GLOBAL_OCEAN_ANALYSIS_*
 set(jjob_list "JGDAS_GLOBAL_OCEAN_ANALYSIS_PREP"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_BMAT"
-#              "JGDAS_GLOBAL_OCEAN_ANALYSIS_BMAT_VRFY"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_RUN"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_CHKPT"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_POST"

test/soca/gw/prep.sh

@@ -60,3 +60,21 @@ ice_files=$(ls ${clim_ens_dir}/ice.*.nc)
 for ice_file in ${ice_files}; do
     ncrename -O -d ni,xaxis_1 -d nj,yaxis_1 -v aice_h,aicen -v hi_h,hicen -v hs_h,hsnon ${ice_file}
 done
+
+# Invent ensemble forecast for f009
+fake_ocean_members=`ls`
+fake_ice_members=`ls ${clim_ens_dir/ice.*.nc}`
+COMENS=${project_binary_dir}/test/soca/gw/COM/enkfgdas.20180415
+for mem in  {1..4}
+do
+    echo "member: $mem"
+    # ocean member
+    oceandir=${COMENS}/06/mem00${mem}/model_data/ocean/history
+    mkdir -p $oceandir
+    cp ${clim_ens_dir}/ocn.${mem}.nc $oceandir/enkfgdas.t06z.ocnf009.nc
+    echo ${clim_ens_dir}/ocn.${mem}.nc $oceandir/enkfgdas.t06z.ocnf009.nc
+    # ice member
+    icedir=${COMENS}/06/mem00${mem}/model_data/ice/history
+    mkdir -p $icedir
+    cp ${clim_ens_dir}/ice.${mem}.nc $icedir/enkfgdas.t06z.icef009.nc
+done

test/soca/gw/run_jjobs.yaml.test

@@ -43,6 +43,8 @@ setup_expt config:
     DO_JEDIOCNVAR: "YES"
     DO_JEDILANDDA: "NO"
     DO_MERGENSST: "NO"
+    NMEM_ENS: "4"
+    DOHYBVAR: "YES"
   ocnanal:
     SOCA_INPUT_FIX_DIR: @HOMEgfs@/sorc/gdas.cd/build/soca_static
     CASE_ANL: C48

ush/soca/run_jjobs.py

@@ -190,6 +190,14 @@ def copy_bkgs(self):
         self.f.write(f"cp {com_ice_history_src}/*icef*.nc $COM_ICE_HISTORY_PREV \n")
         self.f.write(f"cp {com_ice_restart_src}/*cice_model*.nc $COM_ICE_RESTART_PREV \n")
 
+        # copy ensemble members
+        ensbkgs = os.path.join(self.com_src, f"enkf{self.RUN}.{self.gPDY}")
+        if os.path.exists(os.path.join(ensbkgs, self.gcyc)):
+            self.f.write(f"cp -r {ensbkgs} $ROTDIR \n")
+        else:
+            print('Aborting, ensemble backgrounds not found')
+            sys.exit()
+
     def fixconfigs(self):
         """
         Replace cone of the env. var. in the configs