Update and add scripts for analysis stats capabilities

env/HERA.env

@@ -105,6 +105,11 @@ elif [[ "${step}" = "atmanlfv3inc" ]]; then
     export NTHREADS_ATMANLFV3INC=${NTHREADSmax}
     export APRUN_ATMANLFV3INC="${APRUN_default} --cpus-per-task=${NTHREADS_ATMANLFV3INC}"
 
+elif [[ "${step}" = "anlstat" ]]; then
+
+    export NTHREADS_ANLSTAT=${NTHREADSmax}
+    export APRUN_ANLSTAT="${APRUN_default} --cpus-per-task=${NTHREADS_ANLSTAT}"
+
 elif [[ "${step}" = "prepobsaero" ]]; then
 
     export NTHREADS_PREPOBSAERO=${NTHREADS1}

jobs/JGLOBAL_ANALYSIS_STATS

@@ -0,0 +1,38 @@
+#! /usr/bin/env bash
+
+source "${HOMEgfs}/ush/preamble.sh"
+source "${HOMEgfs}/ush/jjob_header.sh" -e "anlstat" -c "base anlstat"
+
+##############################################
+# Set variables used in the script
+##############################################
+
+
+##############################################
+# Begin JOB SPECIFIC work
+##############################################
+
+# Generate COM variables from templates
+YMD=${PDY} HH=${cyc} declare_from_tmpl -rx COM_OBS COM_CHEM_ANALYSIS
+
+
+###############################################################
+# Run relevant script
+
+EXSCRIPT=${ANLSTATSPY:-${SCRgfs}/exglobal_analysis_stats.py}
+${EXSCRIPT}
+status=$?
+[[ ${status} -ne 0 ]] && exit "${status}"
+
+##############################################
+# End JOB SPECIFIC work
+##############################################
+
+##############################################
+# Final processing
+##############################################
+if [[ -e "${pgmout}" ]] ; then
+  cat "${pgmout}"
+fi
+
+exit 0

parm/config/gfs/config.anlstat

@@ -0,0 +1,16 @@
+#!/bin/bash -x
+
+########## config.anlstat ##########
+# Analysis Stat
+
+echo "BEGIN: config.anlstat"
+
+# Get task specific resources
+source "${EXPDIR}/config.resources" anlstat
+
+export JEDI_CONFIG_YAML="${PARMgfs}/gdas/anlstat_jedi_config.yaml.j2"
+export JCB_BASE_YAML="${PARMgfs}/gdas/stat/aero/jcb-base.yaml.j2"
+export JCB_ALGO_YAML="${PARMgfs}/gdas/jcb-algorithms/anlstat.yaml.j2"
+export JEDIEXE=${HOMEgfs}/sorc/gdas.cd/build/bin/ioda-stats.x
+
+echo "END: config.anlstat"

parm/config/gfs/stat/aero/jcb-base.yaml.j2

@@ -0,0 +1,24 @@
+# Search path for model and obs for JCB
+# -------------------------------------
+algorithm_path: "{{PARMgfs}}/gdas/jcb-algorithms"
+app_path_algorithm: "{{PARMgfs}}/gdas/jcb-gdas/algorithm/obstats/aero"
+
+window_begin: "{{ STAT_WINDOW_BEGIN | to_isotime }}"
+window_YMDH: "{{ STAT_WINDOW_BEGIN | to_YMDH }}"
+window_length: "{{ STAT_WINDOW_LENGTH }}"
+
+# Inputted list of ob spaces
+# --------------------------
+obspaces: {{ OBSPACES_LIST }}
+
+# Obspace variable things
+# -----------------------
+aero_obsdatain_path: "{{ DATA }}"
+aero_obsdatatin_simulated_variables: ['aerosolOpticalDepth']
+aero_obsdatain_observed_variables: ['aerosolOpticalDepth']
+
+# Variables
+# ---------
+aero_variables: ['aerosolOpticalDepth']
+aero_file_groups: ['bkgmob', 'bkgmob1']
+aero_file_qc_groups: ['EffectiveQC0', 'EffectiveQC1']

parm/gdas/anlstat_jedi_config.yaml.j2

@@ -0,0 +1,8 @@
+statanl:
+  rundir: '{{ DATA }}'
+  exe_src: '{{ JEDIEXE }}'
+  mpi_cmd: '{{ APRUN_ANLSTAT }}'
+  jedi_args: None
+  jcb_base_yaml: '{{ PARMgfs }}/gdas/stat/aero/jcb-base.yaml.j2'
+  jcb_algo_yaml: '{{ JCB_ALGO_YAML }}'
+  jcb_algo: 'anlstat'

scripts/exglobal_analysis_stats.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+# exglobal_stat_analysis.py
+# This script creates an StatAnalysis class
+# and runs the initialize method
+# which create and stage the runtime directory
+# and create the YAML configuration
+# for a global stat analysis
+import os
+
+from wxflow import Logger, cast_strdict_as_dtypedict
+from pygfs.task.stat_analysis import StatAnalysis
+
+# Initialize root logger
+logger = Logger(level='DEBUG', colored_log=True)
+
+
+if __name__ == '__main__':
+
+    # Take configuration from environment and cast it as python dictionary
+    config = cast_strdict_as_dtypedict(os.environ)
+
+    # Instantiate the atm analysis task
+    StatAnl = StatAnalysis(config)
+
+    # Initialize JEDI variational analysis
+    StatAnl.initialize()

sorc/link_workflow.sh

@@ -356,6 +356,7 @@ if [[ -d "${HOMEgfs}/sorc/gdas.cd/build" ]]; then
                        "gdassoca_obsstats.x" \
                        "gdasapp_land_ensrecenter.x" \
                        "bufr2ioda.x" \
+                       "ioda-stats.x"\
                        "calcfIMS.exe" \
                        "apply_incr.exe" )
   for gdasexe in "${JEDI_EXE[@]}"; do

ush/python/pygfs/task/stat_analysis.py

@@ -0,0 +1,115 @@
+#!/usr/bin/env python3
+
+import os
+import glob
+import gzip
+import tarfile
+import yaml
+from logging import getLogger
+from pprint import pformat
+from typing import Optional, Dict, Any
+
+from wxflow import (AttrDict,
+                    FileHandler,
+                    add_to_datetime, to_fv3time, to_timedelta, to_YMDH,
+                    Task,
+                    parse_j2yaml, save_as_yaml,
+                    logit)
+from pygfs.jedi import Jedi
+
+logger = getLogger(__name__.split('.')[-1])
+
+
+class StatAnalysis(Task):
+    """
+    Class for JEDI-based global stat analysis tasks
+    """
+    @logit(logger, name="StatAnalysis")
+    def __init__(self, config: Dict[str, Any]):
+        """
+        Constructor global stat analysis task
+        This method will construct a global stat analysis task.
+        This includes:
+        - extending the task_config attribute AttrDict to include parameters required for this task
+        - instantiate the Jedi attribute objects
+        Parameters
+        ----------
+        config: Dict
+            dictionary object containing task configuration
+        Returns
+        ----------
+        None
+        """
+        super().__init__(config)
+
+        _res = int(self.task_config.CASE[1:])
+        _window_begin = add_to_datetime(self.task_config.current_cycle, -to_timedelta(f"{self.task_config.assim_freq}H") / 2)
+
+        # Create a local dictionary that is repeatedly used across this class
+        local_dict = AttrDict(
+            {
+                'npx_ges': _res + 1,
+                'npy_ges': _res + 1,
+                'npz_ges': self.task_config.LEVS - 1,
+                'npz': self.task_config.LEVS - 1,
+                'npz_anl': self.task_config.LEVS - 1,
+                'STAT_WINDOW_BEGIN': _window_begin,
+                'STAT_WINDOW_LENGTH': f"PT{self.task_config.assim_freq}H",
+                'OPREFIX': f"{self.task_config.RUN}.t{self.task_config.cyc:02d}z.",
+                'APREFIX': f"{self.task_config.RUN}.t{self.task_config.cyc:02d}z.",
+                'GPREFIX': f"gdas.t{self.task_config.previous_cycle.hour:02d}z."
+            }
+        )
+
+        # Extend task_config with local_dict
+        self.task_config = AttrDict(**self.task_config, **local_dict)
+
+        # Create dictionary of Jedi objects
+        self.jedi_dict = Jedi.get_jedi_dict(self.task_config.JEDI_CONFIG_YAML, self.task_config)
+
+    @logit(logger)
+    def initialize(self) -> None:
+        """
+        Initialize a global stat analysis
+        This method will initialize a global stat analysis.
+        This includes:
+        - initialize JEDI applications
+        - copying stat files
+        Parameters
+        ----------
+        None
+        Returns
+        ----------
+        None
+        """
+        logger.info(f"Copying files to {self.task_config.DATA}/stats")
+
+        # Copy stat files to DATA path
+        aerostat = os.path.join(self.task_config.COM_CHEM_ANALYSIS, f"{self.task_config['APREFIX']}aerostat")
+        dest = os.path.join(self.task_config.DATA, "aerostats")
+        statlist = [[aerostat, dest]]
+        FileHandler({'copy': statlist}).sync()
+
+        # Open tar file
+        logger.info(f"Open tarred stat file in {dest}")
+        with tarfile.open(dest, "r") as tar:
+            # Extract all files to the current directory
+            tar.extractall()
+
+        # Gunzip .nc files
+        logger.info("Gunzip files from tar file")
+        gz_files = glob.glob(os.path.join(self.task_config.DATA, "*gz"))
+
+        for diagfile in gz_files:
+            with gzip.open(diagfile, 'rb') as f_in:
+                with open(diagfile[:-3], 'wb') as f_out:
+                    f_out.write(f_in.read())
+
+        # Get list of .nc4 files
+        obs_space_paths = glob.glob(os.path.join(self.task_config.DATA, "*.nc4"))
+
+        self.task_config.OBSPACES_LIST = ['_'.join(os.path.basename(path).split('_')[1:3]) for path in obs_space_paths]
+
+        # initialize JEDI application
+        logger.info(f"Initializing JEDI variational DA application")
+        self.jedi_dict['statanl'].initialize(self.task_config)

workflow/applications/applications.py

@@ -55,6 +55,7 @@ def __init__(self, conf: Configuration) -> None:
         self.do_verfozn = base.get('DO_VERFOZN', True)
         self.do_verfrad = base.get('DO_VERFRAD', True)
         self.do_vminmon = base.get('DO_VMINMON', True)
+        self.do_anlstat = base.get('DO_ANLSTAT', True)
         self.do_tracker = base.get('DO_TRACKER', True)
         self.do_genesis = base.get('DO_GENESIS', True)
         self.do_genesis_fsu = base.get('DO_GENESIS_FSU', False)

workflow/rocoto/tasks.py

@@ -21,6 +21,7 @@ class Tasks:
                    'atmensanlinit', 'atmensanlobs', 'atmensanlsol', 'atmensanlletkf', 'atmensanlfv3inc', 'atmensanlfinal',
                    'aeroanlinit', 'aeroanlvar', 'aeroanlfinal', 'aeroanlgenb',
                    'snowanl', 'esnowrecen',
+                   'anlstat',
                    'fcst',
                    'atmanlupp', 'atmanlprod', 'atmupp', 'goesupp',
                    'atmos_prod', 'ocean_prod', 'ice_prod',