From 4d99e62d91da76dd0a5c0b1b4247435e500c6c6a Mon Sep 17 00:00:00 2001 From: chan-hoo Date: Fri, 22 Nov 2024 09:44:12 -0600 Subject: [PATCH] add ioda convert capability to prep_obs --- .gitignore | 1 + README.md | 8 +- jobs/JLANDDA_ANALYSIS | 2 - jobs/JLANDDA_FORECAST | 2 - jobs/JLANDDA_PLOT_STATS | 2 - jobs/JLANDDA_POST_ANAL | 2 - jobs/JLANDDA_PREP_OBS | 2 - modulefiles/build_container_intel.lua | 2 +- modulefiles/build_hera_intel.lua | 2 +- modulefiles/build_hercules_intel.lua | 4 +- modulefiles/build_orion_intel.lua | 2 +- modulefiles/build_singularity_intel.lua | 2 +- modulefiles/tasks/hera/task.prep_obs.lua | 11 +- parm/setup_wflow_env.py | 21 +- parm/templates/template.land_analysis.yaml | 27 ++- scripts/exlandda_analysis.sh | 100 ++++----- scripts/exlandda_forecast.sh | 3 + scripts/exlandda_plot_stats.sh | 9 +- scripts/exlandda_post_anal.sh | 6 +- scripts/exlandda_prep_obs.sh | 53 +++-- sorc/test/ci/Dockerfile | 2 +- sorc/test/runtime_vars.sh | 6 +- ush/ghcn_snod2ioda.py | 245 +++++++++++++++++++++ 23 files changed, 376 insertions(+), 138 deletions(-) create mode 100755 ush/ghcn_snod2ioda.py diff --git a/.gitignore b/.gitignore index ef7656fe..cf82bfbc 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ sorc/lib/ sorc/lib64 parm/conda_loc parm/config.yaml +ush/__pycache__/ __pycache__ *.swp diff --git a/README.md b/README.md index b5003f38..78f0f7cb 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,12 @@ -# UFS Offline Land Data Assimilation System +# UFS Land Data Assimilation System The Unified Forecast System (UFS) is a community-based, coupled, comprehensive Earth modeling system. It is designed to be the source system for NOAA's operational numerical weather prediction applications while enabling research, development, and contribution opportunities for the broader Weather Enterprise. For more information about the UFS, visit the UFS Portal at https://ufs.epic.noaa.gov/. -The UFS includes [multiple applications](https://ufs.epic.noaa.gov/applications/) that support different forecast durations and spatial domains. This repository hosts the source code for the UFS Land Data Assimilation (DA) System. Land DA is an offline version of the Noah Multi-Physics (Noah-MP) land surface model (LSM) used in the UFS Weather Model (WM). Its data assimilation framework uses the Joint Effort for Data assimilation Integration (JEDI) software stack, which includes the Object-Oriented Prediction System (OOPS) for the data assimilation algorithm, the Interface for Observation Data Access (IODA) for observation formatting and processing, and the Unified Forward Operator (UFO) for comparing model forecasts and observations. +The UFS includes [multiple applications](https://ufs.epic.noaa.gov/applications/) that support different forecast durations and spatial domains. This repository hosts the source code for the UFS Land Data Assimilation (DA) System. Land DA applies the Noah Multi-Physics (Noah-MP) land surface model (LSM) of the UFS Weather Model (WM) as the key component. Its data assimilation framework uses the Joint Effort for Data assimilation Integration (JEDI) software stack, which includes the Object-Oriented Prediction System (OOPS) for the data assimilation algorithm, the Interface for Observation Data Access (IODA) for observation formatting and processing, and the Unified Forward Operator (UFO) for comparing model forecasts and observations. -The offline Noah-MP LSM is a standalone, uncoupled model used to execute land surface simulations. In this traditional uncoupled mode, near-surface atmospheric forcing data is required as input forcing. This LSM simulates soil moisture (both liquid and frozen), soil temperature, skin temperature, snow depth, snow water equivalent (SWE), snow density, canopy water content, and the energy flux and water flux terms of the surface energy balance and surface water balance. Its data assimilation framework applies the Local Ensemble Transform Kalman Filter-Optimal Interpolation (LETKF-OI) algorithm to combine the state-dependent background error derived from an ensemble forecast with the observations and their corresponding uncertainties to produce an analysis ensemble (Hunt et al., 2007). +The Noah-MP LSM is used to execute land surface simulations. A near-surface atmospheric forcing data is required as input forcing. This LSM simulates soil moisture (both liquid and frozen), soil temperature, skin temperature, snow depth, snow water equivalent (SWE), snow density, canopy water content, and the energy flux and water flux terms of the surface energy balance and surface water balance. Its data assimilation framework applies the Local Ensemble Transform Kalman Filter-Optimal Interpolation (LETKF-OI) algorithm to combine the state-dependent background error derived from an ensemble forecast with the observations and their corresponding uncertainties to produce an analysis ensemble (Hunt et al., 2007). -The Noah-MP LSM has evolved through community efforts to pursue and refine a modern-era LSM suitable for use in the National Centers for Environmental Prediction (NCEP) operational weather and climate prediction models. This collaborative effort continues with participation from entities such as NCAR, NCEP, NASA, and university groups. The development branch of the Land DA System is continually evolving as the system undergoes open development. The latest Land DA release (v1.2.0) represents a snapshot of this continuously evolving system. +The Noah-MP LSM has evolved through community efforts to pursue and refine a modern-era LSM suitable for use in the National Centers for Environmental Prediction (NCEP) operational weather and climate prediction models. This collaborative effort continues with participation from entities such as NCAR, NCEP, NASA, and university groups. The development branch of the Land DA System is continually evolving as the system undergoes open development. The latest Land DA release (v2.0.0) represents a snapshot of this continuously evolving system. The Land DA System User's Guide associated with the development branch is at: https://land-da-workflow.readthedocs.io/en/develop/, while the guide specific to the Land DA v2.0.0 release can be found at: https://land-da-workflow.readthedocs.io/en/release-public-v2.0.0/. Users may download data for use with the most recent release from the [Land DA data bucket](https://registry.opendata.aws/noaa-ufs-land-da/). The [Land DA Docker Hub](https://hub.docker.com/r/noaaepic/ubuntu22.04-intel21.10-landda) hosts Land DA containers. These containers package the Land DA System together with all its software dependencies for an easier experience building and running Land DA. diff --git a/jobs/JLANDDA_ANALYSIS b/jobs/JLANDDA_ANALYSIS index 5dbe8e14..4e7670ad 100755 --- a/jobs/JLANDDA_ANALYSIS +++ b/jobs/JLANDDA_ANALYSIS @@ -80,8 +80,6 @@ mkdir -p ${DATA_SHARE} export DATA_HOFX="${DATA_HOFX:-${DATAROOT}/DATA_SHARE/hofx}" mkdir -p ${DATA_HOFX} -# Set other dates -export PTIME=$($NDATE -24 $PDY$cyc) # #----------------------------------------------------------------------- # diff --git a/jobs/JLANDDA_FORECAST b/jobs/JLANDDA_FORECAST index 0ec431ea..7ea4e8ba 100755 --- a/jobs/JLANDDA_FORECAST +++ b/jobs/JLANDDA_FORECAST @@ -83,8 +83,6 @@ mkdir -p ${DATA_SHARE} export DATA_RESTART="${DATA_RESTART:-${DATAROOT}/DATA_SHARE/RESTART}" mkdir -p ${DATA_RESTART} -# Set other dates -export NTIME=$($NDATE 24 $PDY$cyc) # #----------------------------------------------------------------------- # diff --git a/jobs/JLANDDA_PLOT_STATS b/jobs/JLANDDA_PLOT_STATS index 73fb19b6..f0a83a74 100755 --- a/jobs/JLANDDA_PLOT_STATS +++ b/jobs/JLANDDA_PLOT_STATS @@ -80,8 +80,6 @@ mkdir -p ${DATA_SHARE} export DATA_HOFX="${DATA_HOFX:-${DATAROOT}/DATA_SHARE/hofx}" mkdir -p ${DATA_HOFX} -# Set other dates -export NTIME=$($NDATE 24 $PDY$cyc) # #----------------------------------------------------------------------- # diff --git a/jobs/JLANDDA_POST_ANAL b/jobs/JLANDDA_POST_ANAL index 2bc015bf..aaad6ba3 100755 --- a/jobs/JLANDDA_POST_ANAL +++ b/jobs/JLANDDA_POST_ANAL @@ -76,8 +76,6 @@ mkdir -p ${COMOUT} export DATA_SHARE="${DATA_SHARE:-${DATAROOT}/DATA_SHARE/${PDY}}" mkdir -p ${DATA_SHARE} -# Set other dates -export NTIME=$($NDATE 24 $PDY$cyc) # #----------------------------------------------------------------------- # diff --git a/jobs/JLANDDA_PREP_OBS b/jobs/JLANDDA_PREP_OBS index bcd8a267..3e6bb7ec 100755 --- a/jobs/JLANDDA_PREP_OBS +++ b/jobs/JLANDDA_PREP_OBS @@ -78,8 +78,6 @@ mkdir -p ${COMOUTobs} export DATA_SHARE="${DATA_SHARE:-${DATAROOT}/DATA_SHARE/${PDY}}" mkdir -p ${DATA_SHARE} -# Set other dates -export PTIME=$($NDATE -24 $PDY$cyc) # #----------------------------------------------------------------------- # diff --git a/modulefiles/build_container_intel.lua b/modulefiles/build_container_intel.lua index f7c0036d..d34ea536 100644 --- a/modulefiles/build_container_intel.lua +++ b/modulefiles/build_container_intel.lua @@ -71,5 +71,5 @@ setenv("CC", "mpiicc") setenv("CXX", "mpiicpc") setenv("FC", "mpiifort") -setenv("JEDI_INSTALL", pathJoin(os.getenv("EPICHOME"),"")) +setenv("JEDI_PATH", pathJoin(os.getenv("EPICHOME"),"")) diff --git a/modulefiles/build_hera_intel.lua b/modulefiles/build_hera_intel.lua index 7ce855a2..4b7fa2d0 100644 --- a/modulefiles/build_hera_intel.lua +++ b/modulefiles/build_hera_intel.lua @@ -47,4 +47,4 @@ setenv("FC", "mpiifort") setenv("CMAKE_Platform", "hera.intel") setenv("EPICHOME", "/scratch2/NAGAPE/epic/UFS_Land-DA_Dev") -setenv("JEDI_INSTALL", "/scratch2/NAGAPE/epic/UFS_Land-DA_Dev/jedi_v7") +setenv("JEDI_PATH", "/scratch2/NAGAPE/epic/UFS_Land-DA_Dev/jedi_v7") diff --git a/modulefiles/build_hercules_intel.lua b/modulefiles/build_hercules_intel.lua index 43d2b65e..b6d871d0 100644 --- a/modulefiles/build_hercules_intel.lua +++ b/modulefiles/build_hercules_intel.lua @@ -49,5 +49,5 @@ setenv("CXX", "mpiicpc") setenv("FC", "mpiifort") setenv("CMAKE_Platform", "hercules.intel") -setenv("EPICHOME", "/work/noaa/epic/UFS_Land-DA_Dev") -setenv("JEDI_INSTALL", "/work/noaa/epic/UFS_Land-DA_Dev/jedi_v7_hercules") +setenv("EPICHOME", "/work/noaa/epic/UFS_Land-DA_v2.1") +setenv("JEDI_PATH", "/work/noaa/epic/UFS_Land-DA_v2.1/jedi_v7_ic_hercules") diff --git a/modulefiles/build_orion_intel.lua b/modulefiles/build_orion_intel.lua index dbb94e52..547fca38 100644 --- a/modulefiles/build_orion_intel.lua +++ b/modulefiles/build_orion_intel.lua @@ -50,4 +50,4 @@ setenv("FC", "mpiifort") setenv("CMAKE_Platform", "orion.intel") setenv("EPICHOME", "/work/noaa/epic/UFS_Land-DA_Dev") -setenv("JEDI_INSTALL", "/work/noaa/epic/UFS_Land-DA_Dev/jedi_v7_stack1.6") +setenv("JEDI_PATH", "/work/noaa/epic/UFS_Land-DA_Dev/jedi_v7_stack1.6") diff --git a/modulefiles/build_singularity_intel.lua b/modulefiles/build_singularity_intel.lua index 37d362c5..05de9b1a 100644 --- a/modulefiles/build_singularity_intel.lua +++ b/modulefiles/build_singularity_intel.lua @@ -100,6 +100,6 @@ setenv("CC", "mpiicc") setenv("CXX", "mpiicpc") setenv("FC", "mpiifort") -setenv("JEDI_INSTALL", pathJoin(os.getenv("EPICHOME"),"")) +setenv("JEDI_PATH", pathJoin(os.getenv("EPICHOME"),"")) whatis("Description: UFS build environment") diff --git a/modulefiles/tasks/hera/task.prep_obs.lua b/modulefiles/tasks/hera/task.prep_obs.lua index 4f16978f..e7bcb54c 100644 --- a/modulefiles/tasks/hera/task.prep_obs.lua +++ b/modulefiles/tasks/hera/task.prep_obs.lua @@ -1,7 +1,16 @@ -prepend_path("MODULEPATH", os.getenv("modulepath_spack_stack")) +prepend_path("MODULEPATH", os.getenv("modulepath_spack_stack_unienv")) load(pathJoin("stack-intel", stack_intel_ver)) load(pathJoin("stack-intel-oneapi-mpi", stack_intel_oneapi_mpi_ver)) load(pathJoin("stack-python", stack_python_ver)) load(pathJoin("prod_util", prod_util_ver)) + +load(pathJoin("py-cartopy", py_cartopy_ver)) +load(pathJoin("py-matplotlib", py_matplotlib_ver)) +load(pathJoin("py-netcdf4", py_netcdf4_ver)) +load(pathJoin("py-numpy", py_numpy_ver)) +load(pathJoin("py-pandas", py_pandas_ver)) +load(pathJoin("py-pyyaml", py_pyyaml_ver)) +load(pathJoin("py-xarray", py_xarray_ver)) + diff --git a/parm/setup_wflow_env.py b/parm/setup_wflow_env.py index b6d8fab3..f49ef284 100755 --- a/parm/setup_wflow_env.py +++ b/parm/setup_wflow_env.py @@ -205,7 +205,8 @@ def set_default_parm(): "fhrot": 0, "ic_data_model": "GFS", "imo": 384, - "jedi_install": "/path/to/jedi/install/dir", + "jedi_path": "/path/to/jedi/install/dir", + "jedi_py_ver": "/python/version/used/for/jedi", "jmo": 190, "lnd_calc_snet": ".true.", "lnd_ic_type": "custom", @@ -220,8 +221,7 @@ def set_default_parm(): "nprocs_analysis": 6, "nprocs_fcst_ic": 36, "obsdir": "", - "obsdir_subdir": "", - "obs_types": "GHCN", + "obs_ghcn": "YES", "output_fh": "1 -1", "res": 96, "restart_interval": "12 -1", @@ -242,24 +242,29 @@ def set_machine_parm(machine): lowercase_machine = machine.lower() match lowercase_machine: case "hera": - jedi_install = "/scratch2/NAGAPE/epic/UFS_Land-DA_Dev/jedi_v7" + jedi_path = "/scratch2/NAGAPE/epic/UFS_Land-DA_Dev/jedi_v7" + jedi_py_ver = "python3.11" warmstart_dir = "/scratch2/NAGAPE/epic/UFS_Land-DA_v2.1/inputs/DATA_RESTART" max_cores_per_node = 40 case "orion": - jedi_install = "/work/noaa/epic/UFS_Land-DA_Dev/jedi_v7_stack1.6" + jedi_path = "/work/noaa/epic/UFS_Land-DA_Dev/jedi_v7_stack1.6" + jedi_py_ver = "python3.10" warmstart_dir = "/work/noaa/epic/UFS_Land-DA_v2.1/inputs/DATA_RESTART" max_cores_per_node = 40 case "hercules": - jedi_install = "/work/noaa/epic/UFS_Land-DA_Dev/jedi_v7_hercules" + jedi_path = "/work/noaa/epic/UFS_Land-DA_v2.1/jedi_v7_ic_hercules" + jedi_py_ver = "python3.10" warmstart_dir = "/work/noaa/epic/UFS_Land-DA_v2.1/inputs/DATA_RESTART" max_cores_per_node = 80 case "singularity": - jedi_install = "SINGULARITY_WORKING_DIR" + jedi_path = "SINGULARITY_WORKING_DIR" + jedi_py_ver = "python3.11" warmstart_dir = "SINGULARITY_WORKING_DIR" max_cores_per_node = 40 machine_config = { - "jedi_install": jedi_install, + "jedi_path": jedi_path, + "jedi_py_ver": jedi_py_ver, "warmstart_dir": warmstart_dir, "max_cores_per_node": max_cores_per_node, } diff --git a/parm/templates/template.land_analysis.yaml b/parm/templates/template.land_analysis.yaml index eebc9f6c..ae1b1f9a 100644 --- a/parm/templates/template.land_analysis.yaml +++ b/parm/templates/template.land_analysis.yaml @@ -25,7 +25,7 @@ workflow: CCPP_SUITE: "{{ ccpp_suite }}" COLDSTART: "{{ coldstart }}" COUPLER_CALENDAR: "{{ coupler_calendar }}" - DAtype: "letkfoi_snow" + DATE_CYCLE_FREQ_HR: "{{ date_cycle_freq_hr }}" DATE_FIRST_CYCLE: "{{ date_first_cycle }}" DT_ATMOS: "{{ dt_atmos }}" DT_RUNSEQ: "{{ dt_runseq }}" @@ -35,7 +35,8 @@ workflow: FHROT: "{{ fhrot }}" IC_DATA_MODEL: "{{ ic_data_model }}" IMO: "{{ imo }}" - JEDI_INSTALL: "{{ jedi_install }}" + JEDI_PATH: "{{ jedi_path }}" + JEDI_PY_VER: "{{ jedi_py_ver }}" JMO: "{{ jmo }}" KEEPDATA: "YES" LND_CALC_SNET: "{{ lnd_calc_snet }}" @@ -56,8 +57,7 @@ workflow: NPROCS_FORECAST_LND: "{{ nprocs_forecast_lnd }}" NPROCS_PER_NODE: "{{ nprocs_per_node }}" OBSDIR: "{{ obsdir }}" - OBSDIR_SUBDIR: "{{ obsdir_subdir }}" - OBS_TYPES: "{{ obs_types }}" + OBS_GHCN: "{{ obs_ghcn }}" OUTPUT_FH: "{{ output_fh }}" RES: "{{ res }}" RESTART_INTERVAL: "{{ restart_interval }}" @@ -95,20 +95,22 @@ workflow: COMROOT: "&COMROOT;" cyc: "&cyc;" DATAROOT: "&DATAROOT;" + DATE_CYCLE_FREQ_HR: "&DATE_CYCLE_FREQ_HR;" HOMElandda: "&HOMElandda;" + JEDI_PATH: "&JEDI_PATH;" + JEDI_PY_VER: "&JEDI_PY_VER;" KEEPDATA: "&KEEPDATA;" MACHINE: "&MACHINE;" model_ver: "&model_ver;" OBSDIR: "&OBSDIR;" - OBSDIR_SUBDIR: "&OBSDIR_SUBDIR;" - OBS_TYPES: "&OBS_TYPES;" + OBS_GHCN: "&OBS_GHCN;" PDY: "&PDY;" SCHED: "&SCHED;" account: "&ACCOUNT;" command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "prep_obs" "&HOMElandda;" "&MACHINE;"' jobname: prep_obs cores: 1 - walltime: 00:02:00 + walltime: 00:20:00 queue: batch join: "&LOGDIR;/prep_obs&LOGFN_SUFFIX;" {%- if coldstart == "YES" %} @@ -192,16 +194,16 @@ workflow: COMROOT: "&COMROOT;" COUPLER_CALENDAR: "&COUPLER_CALENDAR;" cyc: "&cyc;" - DAtype: "&DAtype;" DATAROOT: "&DATAROOT;" + DATE_CYCLE_FREQ_HR: "&DATE_CYCLE_FREQ_HR;" KEEPDATA: "&KEEPDATA;" HOMElandda: "&HOMElandda;" - JEDI_INSTALL: "&JEDI_INSTALL;" + JEDI_PATH: "&JEDI_PATH;" LOGDIR: "&LOGDIR;" MACHINE: "&MACHINE;" model_ver: "&model_ver;" NPROCS_ANALYSIS: "&NPROCS_ANALYSIS;" - OBS_TYPES: "&OBS_TYPES;" + OBS_GHCN: "&OBS_GHCN;" PDY: "&PDY;" RES: "&RES;" SCHED: "&SCHED;" @@ -232,6 +234,7 @@ workflow: COMROOT: "&COMROOT;" cyc: "&cyc;" DATAROOT: "&DATAROOT;" + DATE_CYCLE_FREQ_HR: "&DATE_CYCLE_FREQ_HR;" FCSTHR: "&FCSTHR;" HOMElandda: "&HOMElandda;" KEEPDATA: "&KEEPDATA;" @@ -272,9 +275,9 @@ workflow: COLDSTART: "&COLDSTART;" COMROOT: "&COMROOT;" cyc: "&cyc;" - DAtype: "&DAtype;" DATAROOT: "&DATAROOT;" DATE_FIRST_CYCLE: "&DATE_FIRST_CYCLE;" + DATE_CYCLE_FREQ_HR: "&DATE_CYCLE_FREQ_HR;" DT_ATMOS: "&DT_ATMOS;" DT_RUNSEQ: "&DT_RUNSEQ;" FCSTHR: "&FCSTHR;" @@ -299,7 +302,6 @@ workflow: NPROCS_FORECAST_ATM: "&NPROCS_FORECAST_ATM;" NPROCS_FORECAST_LND: "&NPROCS_FORECAST_LND;" NPROCS_PER_NODE: "&NPROCS_PER_NODE;" - OBS_TYPES: "&OBS_TYPES;" OUTPUT_FH: "&OUTPUT_FH;" PDY: "&PDY;" RES: "&RES;" @@ -339,6 +341,7 @@ workflow: COMROOT: "&COMROOT;" cyc: "&cyc;" DATAROOT: "&DATAROOT;" + DATE_CYCLE_FREQ_HR: "&DATE_CYCLE_FREQ_HR;" HOMElandda: "&HOMElandda;" KEEPDATA: "&KEEPDATA;" LOGDIR: "&LOGDIR;" diff --git a/scripts/exlandda_analysis.sh b/scripts/exlandda_analysis.sh index d18003e3..a54876f0 100755 --- a/scripts/exlandda_analysis.sh +++ b/scripts/exlandda_analysis.sh @@ -2,6 +2,9 @@ set -xue +# Set other dates +PTIME=$($NDATE -{DATE_CYCLE_FREQ_HR} $PDY$cyc) + YYYY=${PDY:0:4} MM=${PDY:4:2} DD=${PDY:6:2} @@ -13,8 +16,8 @@ HP=${PTIME:8:2} FILEDATE=${YYYY}${MM}${DD}.${HH}0000 -JEDI_STATICDIR=${JEDI_INSTALL}/jedi-bundle/fv3-jedi/test/Data -JEDI_EXECDIR=${JEDI_INSTALL}/build/bin +JEDI_STATICDIR=${JEDI_PATH}/jedi-bundle/fv3-jedi/test/Data +JEDI_EXECDIR=${JEDI_PATH}/build/bin case $MACHINE in "hera") @@ -31,7 +34,6 @@ case $MACHINE in ;; esac -YAML_DA=construct GFSv17="NO" B=30 # back ground error std for LETKFOI @@ -63,33 +65,28 @@ ${USHlandda}/fill_jinja_template.py -u "${settings}" -t "${fp_template}" -o "${f # CREATE BACKGROUND ENSEMBLE (LETKFOI) ################################################ -if [[ ${DAtype} == "letkfoi_snow" ]]; then - - if [ $GFSv17 == "YES" ]; then - SNOWDEPTHVAR="snodl" - else - SNOWDEPTHVAR="snwdph" - # replace field overwrite file - cp ${PARMlandda}/jedi/gfs-land.yaml ${DATA}/gfs-land.yaml +if [ $GFSv17 == "YES" ]; then + SNOWDEPTHVAR="snodl" +else + SNOWDEPTHVAR="snwdph" + # replace field overwrite file + cp ${PARMlandda}/jedi/gfs-land.yaml ${DATA}/gfs-land.yaml +fi +# FOR LETKFOI, CREATE THE PSEUDO-ENSEMBLE +for ens in pos neg +do + if [ -e $DATA/mem_${ens} ]; then + rm -r $DATA/mem_${ens} fi - # FOR LETKFOI, CREATE THE PSEUDO-ENSEMBLE - for ens in pos neg - do - if [ -e $DATA/mem_${ens} ]; then - rm -r $DATA/mem_${ens} - fi - mkdir -p $DATA/mem_${ens} - cp ${FILEDATE}.sfc_data.tile*.nc ${DATA}/mem_${ens} - cp ${DATA}/${FILEDATE}.coupler.res ${DATA}/mem_${ens}/${FILEDATE}.coupler.res - done - - echo 'do_landDA: calling create ensemble' + mkdir -p $DATA/mem_${ens} + cp ${FILEDATE}.sfc_data.tile*.nc ${DATA}/mem_${ens} + cp ${DATA}/${FILEDATE}.coupler.res ${DATA}/mem_${ens}/${FILEDATE}.coupler.res +done - # using ioda mods to get a python version with netCDF4 - ${USHlandda}/letkf_create_ens.py $FILEDATE $SNOWDEPTHVAR $B - if [[ $? != 0 ]]; then - err_exit "letkf create failed" - fi +# using ioda mods to get a python version with netCDF4 +${USHlandda}/letkf_create_ens.py $FILEDATE $SNOWDEPTHVAR $B +if [[ $? != 0 ]]; then + err_exit "letkf create failed" fi ################################################ @@ -104,15 +101,9 @@ mkdir -p output/DA/hofx # if yaml is specified by user, use that. Otherwise, build the yaml if [[ $do_DA == "YES" ]]; then - if [[ $YAML_DA == "construct" ]];then # construct the yaml - cp ${PARMlandda}/jedi/${DAtype}.yaml ${DATA}/letkf_land.yaml - for obs in "${OBS_TYPES[@]}"; - do - cat ${PARMlandda}/jedi/${obs}.yaml >> letkf_land.yaml - done - else # use specified yaml - echo "Using user specified YAML: ${YAML_DA}" - cp ${PARMlandda}/jedi/${YAML_DA} ${DATA}/letkf_land.yaml + cp "${PARMlandda}/jedi/letkfoi_snow.yaml" "${DATA}/letkf_land.yaml" + if [ "${OBS_GHCN}" = "YES" ]; then + cat ${PARMlandda}/jedi/GHCN.yaml >> letkf_land.yaml fi # update jedi yaml file @@ -141,15 +132,9 @@ fi if [[ $do_HOFX == "YES" ]]; then - if [[ $YAML_HOFX == "construct" ]];then # construct the yaml - cp ${PARMlandda}/jedi/${DAtype}.yaml ${DATA}/hofx_land.yaml - for obs in "${OBS_TYPES[@]}"; - do - cat ${PARMlandda}/jedi/${obs}.yaml >> hofx_land.yaml - done - else # use specified yaml - echo "Using user specified YAML: ${YAML_HOFX}" - cp ${PARMlandda}/jedi/${YAML_HOFX} ${DATA}/hofx_land.yaml + cp "${PARMlandda}/jedi/letkfoi_snow.yaml" "${DATA}/hofx_land.yaml" + if [ "${OBS_GHCN}" = "YES" ]; then + cat ${PARMlandda}/jedi/GHCN.yaml >> hofx_land.yaml fi # update jedi yaml file @@ -179,7 +164,7 @@ fi if [[ "$GFSv17" == "NO" ]]; then cp ${PARMlandda}/jedi/gfs-land.yaml ${DATA}/gfs-land.yaml else - cp ${JEDI_INSTALL}/jedi-bundle/fv3-jedi/test/Data/fieldmetadata/gfs_v17-land.yaml ${DATA}/gfs-land.yaml + cp ${JEDI_PATH}/jedi-bundle/fv3-jedi/test/Data/fieldmetadata/gfs_v17-land.yaml ${DATA}/gfs-land.yaml fi ################################################ @@ -219,8 +204,6 @@ fi if [[ $do_DA == "YES" ]]; then - if [[ $DAtype == "letkfoi_snow" ]]; then - cat << EOF > apply_incr_nml &noahmp_snow date_str=${YYYY}${MM}${DD} @@ -232,17 +215,14 @@ cat << EOF > apply_incr_nml / EOF - echo 'do_landDA: calling apply snow increment' - - export pgm="apply_incr.exe" - . prep_step - # (n=6) -> this is fixed, at one task per tile (with minor code change, could run on a single proc). - ${RUN_CMD} -n 6 ${EXEClandda}/$pgm >>$pgmout 2>errfile - export err=$?; err_chk - cp errfile errfile_apply_incr - if [[ $err != 0 ]]; then - err_exit "apply snow increment failed" - fi + export pgm="apply_incr.exe" + . prep_step + # (n=6) -> this is fixed, at one task per tile (with minor code change, could run on a single proc). + ${RUN_CMD} -n 6 ${EXEClandda}/$pgm >>$pgmout 2>errfile + export err=$?; err_chk + cp errfile errfile_apply_incr + if [[ $err != 0 ]]; then + err_exit "apply snow increment failed" fi for itile in {1..6} diff --git a/scripts/exlandda_forecast.sh b/scripts/exlandda_forecast.sh index 97d6def1..ca852fad 100755 --- a/scripts/exlandda_forecast.sh +++ b/scripts/exlandda_forecast.sh @@ -28,6 +28,9 @@ export ESMF_RUNTIME_PROFILE_OUTPUT="SUMMARY" export PSM_RANKS_PER_CONTEXT=4 export PSM_SHAREDCONTEXTS=1 + +NTIME=$($NDATE ${DATE_CYCLE_FREQ_HR} $PDY$cyc) + YYYY=${PDY:0:4} MM=${PDY:4:2} DD=${PDY:6:2} diff --git a/scripts/exlandda_plot_stats.sh b/scripts/exlandda_plot_stats.sh index 92569f42..67d0bf34 100755 --- a/scripts/exlandda_plot_stats.sh +++ b/scripts/exlandda_plot_stats.sh @@ -2,6 +2,9 @@ set -xue +# Set other dates +NTIME=$($NDATE {DATE_CYCLE_FREQ_HR} $PDY$cyc) + YYYY=${PDY:0:4} MM=${PDY:4:2} DD=${PDY:6:2} @@ -44,7 +47,7 @@ machine: '${MACHINE}' EOF ${USHlandda}/hofx_analysis_stats.py -if [[ $? != 0 ]]; then +if [ $? -ne 0 ]; then err_exit "Scatter/Histogram plots failed" fi @@ -81,7 +84,7 @@ out_fn_time: '${OUT_FN_TIME}' EOF ${USHlandda}/plot_analysis_timehistory.py -if [[ $? != 0 ]]; then +if [ $? -ne 0 ]; then err_exit "Time-history plots failed" fi @@ -112,7 +115,7 @@ machine: '${MACHINE}' EOF ${USHlandda}/plot_forecast_restart.py -if [[ $? != 0 ]]; then +if [ $? -ne 0 ]; then err_exit "Forecast restart plots failed" fi diff --git a/scripts/exlandda_post_anal.sh b/scripts/exlandda_post_anal.sh index 8f1df85e..3b87e144 100755 --- a/scripts/exlandda_post_anal.sh +++ b/scripts/exlandda_post_anal.sh @@ -2,11 +2,11 @@ set -xue -############################ -# copy restarts to workdir, convert to UFS tile for DA (all members) - MACHINE_ID=${MACHINE} +# Set other dates +NTIME=$($NDATE ${DATE_CYCLE_FREQ_HR} $PDY$cyc) + YYYY=${PDY:0:4} MM=${PDY:4:2} DD=${PDY:6:2} diff --git a/scripts/exlandda_prep_obs.sh b/scripts/exlandda_prep_obs.sh index 3b2d39ba..37a6b8ef 100755 --- a/scripts/exlandda_prep_obs.sh +++ b/scripts/exlandda_prep_obs.sh @@ -2,6 +2,8 @@ set -xue +# Set other dates +PTIME=$($NDATE -${DATE_CYCLE_FREQ_HR} $PDY$cyc) YYYY=${PDY:0:4} MM=${PDY:4:2} @@ -13,33 +15,30 @@ DP=${PTIME:6:2} HP=${PTIME:8:2} OBSDIR="${OBSDIR:-${FIXlandda}/DA_obs}" -for obs in "${OBS_TYPES[@]}"; do - # get the obs file name - if [ "${obs}" == "GTS" ]; then - OBSDIR_SUBDIR="${OBSDIR_SUBDIR:-snow_depth/GTS/data_proc}" - obsfile="${OBSDIR}/${OBSDIR_SUBDIR}/${YYYY}${MM}/adpsfc_snow_${YYYY}${MM}${DD}${HH}.nc4" - elif [ "${obs}" == "GHCN" ]; then - # GHCN are time-stamped at 18. If assimilating at 00, need to use previous day's obs, so that - # obs are within DA window. - if [ "${ATMOS_FORC}" == "era5" ]; then - OBSDIR_SUBDIR="${OBSDIR_SUBDIR:-snow_depth/GHCN/data_proc/v3}" - obsfile="${OBSDIR}/${OBSDIR_SUBDIR}/${YYYY}/ghcn_snwd_ioda_${YYYP}${MP}${DP}_jediv7.nc" - elif [ "${ATMOS_FORC}" == "gswp3" ]; then - OBSDIR_SUBDIR="${OBSDIR_SUBDIR:-snow_depth/GHCN/data_proc/v3}" - obsfile="${OBSDIR}/${OBSDIR_SUBDIR}/${YYYY}/ghcn_snwd_ioda_${YYYP}${MP}${DP}.nc" - fi - elif [ ${obs} == "SYNTH" ]; then - OBSDIR_SUBDIR="${OBSDIR_SUBDIR:-synthetic_noahmp}" - obsfile="${OBSDIR}/${OBSDIR_SUBDIR}/IODA.synthetic_gswp_obs.${YYYY}${MM}${DD}${HH}.nc" - else - err_exit "Unknown obs type requested ${obs}, exiting" - fi +DATA_GHCN_RAW="${DATA_GHCN_RAW:-${FIXlandda}/DATA_ghcn}" + +# GHCN snow depth data +if [ "${OBS_GHCN}" = "YES" ]; then + # GHCN are time-stamped at 18. If assimilating at 00, need to use previous day's obs, + # so that obs are within DA window. + obs_fn="ghcn_snwd_ioda_${YYYP}${MP}${DP}${HP}.nc" + obs_fp="${OBSDIR}/GHCN/${YYYY}/${obs_fn}" + out_fn="GHCN_${YYYY}${MM}${DD}${HH}.nc" - # check obs are available - if [[ -e $obsfile ]]; then - echo "$obs observations found: $obsfile" - cp -p $obsfile ${COMOUTobs}/${obs}_${YYYY}${MM}${DD}${HH}.nc + # check obs is available + if [ -e $obs_fp ]; then + echo "GHCN observation file: $obs_fp" + cp -p $obs_fp ${COMOUTobs}/${out_fn} else - err_exit "${obs} observations not found: $obsfile" + input_ghcn_file="${DATA_GHCN_RAW}/${YYYP}.csv" + output_ioda_file="${obs_fn}" + ghcn_station_file="${DATA_GHCN_RAW}/ghcnd-stations.txt" + + ${USHlandda}/ghcn_snod2ioda.py -i ${input_ghcn_file} -o ${output_ioda_file} -f ${ghcn_station_file} -d ${YYYP}${MP}${DP}${HP} -m maskout + if [ $? -ne 0 ]; then + err_exit "Generation of GHCN obs file failed !!!" + fi + cp -p ${output_ioda_file} ${COMOUTobs}/${out_fn} + fi -done +fi diff --git a/sorc/test/ci/Dockerfile b/sorc/test/ci/Dockerfile index 7c0e5ca6..07b4d81d 100644 --- a/sorc/test/ci/Dockerfile +++ b/sorc/test/ci/Dockerfile @@ -11,7 +11,7 @@ RUN wget https://noaa-ufs-land-da-pds.s3.amazonaws.com/develop-20241024/inputs.t # set env vars ENV FIXlandda=$HOME/land-DA_workflow/fix -ENV JEDI_INSTALL=$HOME +ENV JEDI_PATH=$HOME ENV FIXdir=$FIXlandda ENV JEDI_EXECDIR=/opt/jedi-bundle/install/bin diff --git a/sorc/test/runtime_vars.sh b/sorc/test/runtime_vars.sh index df749b08..d0e22f6b 100755 --- a/sorc/test/runtime_vars.sh +++ b/sorc/test/runtime_vars.sh @@ -36,13 +36,13 @@ export EXECDIR=${EXECDIR:-$project_binary_dir/bin} export FIXlandda=${FIXlandda:-"`dirname $project_source_dir`/fix"} # set IODA path -export IODA_BUILD_DIR=${IODA_BUILD_DIR:-"${JEDI_INSTALL}/build"} +export IODA_BUILD_DIR=${IODA_BUILD_DIR:-"${JEDI_PATH}/build"} export PYTHON_VERSION=`python -c 'import sys; version=sys.version_info[:3]; print("{0}.{1}".format(*version))'` export PYTHONPATH=$PYTHONPATH:${IODA_BUILD_DIR}/lib/python${PYTHON_VERSION}/pyioda:${IODA_BUILD_DIR}/lib/pyiodaconv # JEDI directories -export JEDI_EXECDIR=${JEDI_EXECDIR:-"${JEDI_INSTALL}/build/bin"} -export JEDI_STATICDIR=${JEDI_INSTALL}/jedi-bundle/fv3-jedi/test/Data +export JEDI_EXECDIR=${JEDI_EXECDIR:-"${JEDI_PATH}/build/bin"} +export JEDI_STATICDIR=${JEDI_PATH}/jedi-bundle/fv3-jedi/test/Data # set executables export MPIRUN=${MPIRUN:-`which mpiexec`} diff --git a/ush/ghcn_snod2ioda.py b/ush/ghcn_snod2ioda.py new file mode 100755 index 00000000..4c516576 --- /dev/null +++ b/ush/ghcn_snod2ioda.py @@ -0,0 +1,245 @@ +#!/usr/bin/env python3 +# +# (C) Copyright 2021-2022 NOAA/NWS/NCEP/EMC +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# + +import os, sys +import argparse +import numpy as np +import pandas as pd +from datetime import datetime +from dateutil.parser import parse + +jedi_path = os.environ.get('JEDI_PATH') +jedi_py_ver = os.environ.get('JEDI_PY_VER') +sys.path.append(os.path.join(jedi_path, 'build/lib', jedi_py_ver)) +#sys.path.append(os.path.join(jedi_path, 'build/lib', jedi_py_ver, 'pyiodaconv')) + +import pyiodaconv.ioda_conv_engines as iconv +from collections import defaultdict, OrderedDict +from pyiodaconv.orddicts import DefaultOrderedDict +from pyiodaconv.def_jedi_utils import epoch, iso8601_string +from pyiodaconv.def_jedi_utils import record_time + +locationKeyList = [ + ("latitude", "float"), + ("longitude", "float"), + ("stationElevation", "float"), + ("dateTime", "long"), +] + +AttrData = { +} + +DimDict = { +} + +VarDims = { + 'totalSnowDepth': ['Location'], +} + +float_missing_value = iconv.get_default_fill_val(np.float32) +int_missing_value = iconv.get_default_fill_val(np.int32) +long_missing_value = iconv.get_default_fill_val(np.int64) + + +def get_epoch_time(adatetime): + + # take python datetime object and convert to seconds from epoch + time_offset = round((adatetime - epoch).total_seconds()) + + return time_offset + + +def assignValue(colrowValue, df400): + if colrowValue == '' or pd.isnull(colrowValue): + outList = float_missing_value + else: + ml = df400.loc[df400['ID'] == colrowValue, "DATA_VALUE"] + # check if the series is empty + if not ml.empty: + outList = ml.iloc[0] + else: + outList = float_missing_value + return outList + + +class ghcn(object): + + def __init__(self, filename, fixfile, date, mask): + self.filename = filename + self.fixfile = fixfile + self.date = date + self.mask = mask + self.varDict = defaultdict(lambda: defaultdict(dict)) + self.metaDict = defaultdict(lambda: defaultdict(dict)) + self.outdata = defaultdict(lambda: DefaultOrderedDict(OrderedDict)) + self.varAttrs = defaultdict(lambda: DefaultOrderedDict(OrderedDict)) + self._read() + + def _read(self): + + # set up variable names for IODA + iodavar = 'totalSnowDepth' + self.varDict[iodavar]['valKey'] = iodavar, iconv.OvalName() + self.varDict[iodavar]['errKey'] = iodavar, iconv.OerrName() + self.varDict[iodavar]['qcKey'] = iodavar, iconv.OqcName() + self.varAttrs[iodavar, iconv.OvalName()]['coordinates'] = 'longitude latitude' + self.varAttrs[iodavar, iconv.OerrName()]['coordinates'] = 'longitude latitude' + self.varAttrs[iodavar, iconv.OqcName()]['coordinates'] = 'longitude latitude' + self.varAttrs[iodavar, iconv.OvalName()]['units'] = 'mm' + self.varAttrs[iodavar, iconv.OerrName()]['units'] = 'mm' + + # read in the GHCN data csv file + cols = ["ID", "DATETIME", "ELEMENT", "DATA_VALUE", "M_FLAG", "Q_FLAG", "S_FLAG", "OBS_TIME"] + sub_cols = ["ID", "DATETIME", "ELEMENT", "DATA_VALUE"] + df30_list = [] + # Fix dtypeWarning with mixed types via set low_memory=False + df20all = pd.read_csv(self.filename, header=None, names=cols, low_memory=False) + df20 = df20all[sub_cols] + df20all = None + df30_list.append(df20) + df20 = None + + df30 = pd.concat(df30_list, ignore_index=True) + df30 = df30[df30["ELEMENT"] == "SNWD"] + df30["DATETIME"] = df30.apply(lambda row: parse(str(row["DATETIME"])).date(), axis=1) + # select data which matches the Start date + startdate = self.date + valid_date = datetime.strptime(startdate, "%Y%m%d%H") + select_date = valid_date.strftime('%Y%m%d') + new_date = parse(select_date).date() + df30 = df30[df30["DATETIME"] == new_date] + + # Read in the GHCN station files + cols = ["ID", "LATITUDE", "LONGITUDE", "ELEVATION", "STATE", "NAME", "GSN_FLAG", "HCNCRN_FLAG", "WMO_ID"] + df10all = pd.read_csv(self.fixfile, header=None, sep='\r\n', engine='python') + df10all = df10all[0].str.split('\\s+', expand=True) + df10 = df10all.iloc[:, 0:4] + df10all = None + sub_cols = {0: "ID", 1: "LATITUDE", 2: "LONGITUDE", 3: "ELEVATION"} + df10 = df10.rename(columns=sub_cols) + df10 = df10.drop_duplicates(subset=["ID"]) + + # use stations list as the number of obs points + # if no data for a station and a given date, leave float_missing_value + num_obs = len(df10.index) + + # Initialzed data array + vals = np.full((num_obs), float_missing_value) + lats = np.full((num_obs), float_missing_value) + lons = np.full((num_obs), float_missing_value) + alts = np.full((num_obs), float_missing_value) + id_array = np.chararray((num_obs)) + id_array[:] = "UNKNOWN" + + lats = df10["LATITUDE"].values + lons = df10["LONGITUDE"].values + alts = df10["ELEVATION"].values + id_array = df10["ID"].values + + df100 = pd.DataFrame(data=id_array, columns=['ID']) + df100.assign(DATA_VALUE=float_missing_value) + df30Temp = df30.loc[df30["DATETIME"] == new_date] + df100["DATA_VALUE"] = df100.apply(lambda row: assignValue(row['ID'], df30Temp), axis=1) + df30Temp = None + + vals = df100["DATA_VALUE"].values + vals = vals.astype('float32') + lats = lats.astype('float32') + lons = lons.astype('float32') + alts = alts.astype('float32') + qflg = 0*vals.astype('int32') + errs = 0.0*vals + sites = np.empty_like(vals, dtype=object) + times = np.empty_like(vals, dtype='int64') + sites = id_array + + # use maskout options + if self.mask == "maskout": + + with np.errstate(invalid='ignore'): + mask = (vals >= 0.0) & (vals < float_missing_value) + vals = vals[mask] + errs = errs[mask] + qflg = qflg[mask] + lons = lons[mask] + lats = lats[mask] + alts = alts[mask] + sites = sites[mask] + times = times[mask] + + # get datetime from input + my_date = datetime.strptime(startdate, "%Y%m%d%H") + epoch_time = np.int64(get_epoch_time(my_date)) + + # vals[vals >= 0.0] *= 0.001 # mm to meters + # errs[:] = 0.04 # error in meters + errs[:] = 40.0 # error in mm + times[:] = epoch_time + # add metadata variables + self.outdata[('dateTime', 'MetaData')] = times + self.outdata[('stationIdentification', 'MetaData')] = sites + self.outdata[('latitude', 'MetaData')] = lats + self.outdata[('longitude', 'MetaData')] = lons + self.outdata[('stationElevation', 'MetaData')] = alts + self.varAttrs[('stationElevation', 'MetaData')]['units'] = 'm' + self.varAttrs[('dateTime', 'MetaData')]['units'] = iso8601_string + self.varAttrs[('dateTime', 'MetaData')]['_FillValue'] = long_missing_value + + self.outdata[self.varDict[iodavar]['valKey']] = vals + self.outdata[self.varDict[iodavar]['errKey']] = errs + self.outdata[self.varDict[iodavar]['qcKey']] = qflg + + DimDict['Location'] = len(self.outdata[('dateTime', 'MetaData')]) + + +def main(): + + parser = argparse.ArgumentParser( + description=('Read GHCN snow depth file(s) and Converter' + ' of native csv format for observations of snow' + ' depth from GHCN to IODA netCDF format.') + ) + parser.add_argument('-i', '--input', + help="name of ghcn snow depth input file(s)", + type=str, required=True) + parser.add_argument('-f', '--fixfile', + help="name of ghcn station fixed file(s)", + type=str, required=True) + parser.add_argument('-o', '--output', + help="name of ioda output file", + type=str, required=True) + parser.add_argument('-d', '--date', + help="base date (YYYYMMDDHH)", type=str, required=True) + optional = parser.add_argument_group(title='optional arguments') + optional.add_argument( + '-m', '--mask', + help="maskout missing values: maskout/default, default=none", + type=str, required=True) + + args = parser.parse_args() + + # start timer + tic = record_time() + + # Read in the GHCN snow depth data + snod = ghcn(args.input, args.fixfile, args.date, args.mask) + + # report time + toc = record_time(tic=tic) + + # setup the IODA writer + writer = iconv.IodaWriter(args.output, locationKeyList, DimDict) + + # write all data out + writer.BuildIoda(snod.outdata, VarDims, snod.varAttrs, AttrData) + + # report time + toc = record_time(tic=tic) + + +if __name__ == '__main__': + main()