Skip to content

Commit

Permalink
Merge branch 'develop' into feature/letkf-j2
Browse files Browse the repository at this point in the history
  • Loading branch information
AndrewEichmann-NOAA committed Feb 13, 2025
2 parents 97d4029 + d5083d5 commit 1ee5c50
Show file tree
Hide file tree
Showing 14 changed files with 655 additions and 17 deletions.
63 changes: 59 additions & 4 deletions parm/soca/berror/soca_ensrecenter.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,9 @@ increment variables:
- sea_water_salinity
- eastward_sea_water_velocity
- northward_sea_water_velocity
- sea_ice_area_fraction
- sea_ice_thickness
- sea_ice_snow_thickness
set increment variables to zero:
- eastward_sea_water_velocity
- northward_sea_water_velocity
Expand All @@ -26,6 +29,57 @@ vertical geometry:

add recentering increment: true
recentering around deterministic: true
sea ice recenter: true
sea ice analysis:
pattern: '%mem%'
date: '{{ MARINE_WINDOW_END_ISO }}'
ocn_filename: "ocean.%mem%.nc"
ice_filename: "ice.%mem%.nc"
read_from_file: 1
basename: {{ ENSPERT_RELPATH }}/ens/
state variables:
- sea_water_potential_temperature
- sea_water_salinity
- eastward_sea_water_velocity
- northward_sea_water_velocity
- sea_water_cell_thickness
- sea_water_depth
- ocean_mixed_layer_thickness
- sea_ice_area_fraction
- sea_ice_thickness
- sea_ice_snow_thickness

sea ice variable change:
variable change name: Soca2Cice
arctic:
seaice edge: 0.4
shuffle: true
rescale prior:
rescale: true
min hice: 0.5
min hsno: 0.1
antarctic:
seaice edge: 0.4
shuffle: true
rescale prior:
rescale: true
min hice: 0.5
min hsno: 0.1
cice background state:
pattern: '%mem%'
restart: ens/cice_model.res.%mem%.nc
ncat: 5
ice_lev: 7
sno_lev: 1
cice output:
pattern: '%mem%'
restart: ens/cice_model.res.output.%mem%.nc
output variables:
- sea_water_potential_temperature
- sea_water_salinity
- sea_ice_area_fraction
- sea_ice_thickness
- sea_ice_snow_thickness

soca increments: # Could also be states, but they are read as increments
number of increments: {{ NMEM_ENS }}
Expand All @@ -34,8 +88,7 @@ soca increments: # Could also be states, but they are read as increments
date: '{{ MARINE_WINDOW_END_ISO }}'
basename: ./ens/
ocn_filename: 'ocean.%mem%.nc'
# TODO(AFE) fix ice recentering in cycled da
# ice_filename: 'ice.%mem%.nc'
ice_filename: 'ice.%mem%.nc'
read_from_file: 3

trajectory:
Expand All @@ -48,11 +101,13 @@ trajectory:
- sea_water_cell_thickness
- sea_water_depth
- ocean_mixed_layer_thickness
- sea_ice_area_fraction
- sea_ice_thickness
- sea_ice_snow_thickness
date: '{{ MARINE_WINDOW_END_ISO }}'
basename: ./bkg/
ocn_filename: ocean.bkg.f009.nc
# TODO(AFE) fix ice recentering in cycled da
#ice_filename: cice.res.nc
ice_filename: ice.bkg.f009.nc
read_from_file: 1

output increment:
Expand Down
63 changes: 58 additions & 5 deletions ush/soca/marine_recenter.py
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@ def __init__(self, config: Dict) -> None:

_window_begin = add_to_datetime(self.task_config.current_cycle, -to_timedelta(f"{self.task_config.assim_freq}H") / 2)
_window_end = add_to_datetime(self.task_config.current_cycle, to_timedelta(f"{self.task_config.assim_freq}H") / 2)
_enspert_relpath = os.path.relpath(self.task_config.DATAens, self.task_config.DATA)

local_dict = AttrDict({'window_begin': f"{window_begin.strftime('%Y-%m-%dT%H:%M:%SZ')}",
'PARMsoca': os.path.join(self.task_config.PARMgfs, 'gdas', 'soca'),
Expand All @@ -80,6 +81,8 @@ def __init__(self, config: Dict) -> None:
'stage_dir': DATA,
'soca_input_fix_dir': self.task_config.SOCA_INPUT_FIX_DIR,
'NMEM_ENS': self.task_config.NMEM_ENS,
'GDUMP_ENS': self.task_config.GDUMP_ENS,
'ENSPERT_RELPATH': _enspert_relpath,
'MARINE_WINDOW_LENGTH': f"PT{config['assim_freq']}H",
'recen_yaml_template': os.path.join(berror_yaml_dir, 'soca_ensrecenter.yaml'),
'recen_yaml_file': os.path.join(DATA, 'soca_ensrecenter.yaml'),
Expand Down Expand Up @@ -124,6 +127,9 @@ def initialize(self):
# stage backgrounds
bkg_list = parse_j2yaml(self.task_config.MARINE_DET_STAGE_BKG_YAML_TMPL, self.task_config)
FileHandler(bkg_list).sync()
# stage ensemble backgrounds for soca2cice
ens_bkg_list = parse_j2yaml(self.task_config.MARINE_ENSDA_STAGE_BKG_YAML_TMPL, self.task_config)
FileHandler(ens_bkg_list).sync()

# ################################################################################
# # Copy initial condition
Expand Down Expand Up @@ -155,6 +161,37 @@ def initialize(self):
ens_member_list.append([fname_in, fname_out])

FileHandler({'copy': ens_member_list}).sync()
# stage ensemble ice restarts
# make a copy of the CICE6 restart
logger.info("---------------- Stage ensemble CICE restarts")
# set the restart date, dependent on the cycling type
if self.task_config.DOIAU:
# forecast initialized at the begining of the DA window
fcst_begin = self.task_config.MARINE_WINDOW_BEGIN_ISO
rst_date = self.task_config.MARINE_WINDOW_BEGIN.strftime('%Y%m%d.%H%M%S')
else:
# forecast initialized at the middle of the DA window
fcst_begin = self.task_config.MARINE_WINDOW_MIDDLE_ISO
rst_date = self.task_config.MARINE_WINDOW_MIDDLE.strftime('%Y%m%d.%H%M%S')
ens_cice_list = []
for mem in range(1, nmem_ens+1):
mem_dir = os.path.join(self.task_config.ROTDIR,
f'enkf{RUN}.{gPDYstr}',
f'{gcyc}',
f'mem{str(mem).zfill(3)}',
'model',
'ice',
'restart')
mem_dir_real = os.path.realpath(mem_dir)
f00 = f'{rst_date}.cice_model.res.nc'
fname_in = os.path.abspath(os.path.join(mem_dir_real, f00))
fname_out = os.path.realpath(os.path.join(self.task_config.ens_dir,
"cice_model.res."+str(mem)+".nc"))
ens_cice_list.append([fname_in, fname_out])
fname_out = os.path.realpath(os.path.join(self.task_config.ens_dir,
"cice_model.res.output."+str(mem)+".nc"))
ens_cice_list.append([fname_in, fname_out])
FileHandler({'copy': ens_cice_list}).sync()

################################################################################
# generate the YAML file for recenterer
Expand Down Expand Up @@ -228,10 +265,6 @@ def finalize(self):
mem_dir_list = []
copy_list = []

# Skip the analysis insertion into the CICE restart for now
# TODO (G): Add this back in when we have hardened the soca to cice
# change of variable

# Copy the recentering increment files to the member COMROOT directories
for mem in range(1, nmem_ens+1):
mem_dir = os.path.join(self.task_config.ROTDIR,
Expand All @@ -242,9 +275,29 @@ def finalize(self):
'ocean')
mem_dir_real = os.path.realpath(mem_dir)
mem_dir_list.append(mem_dir_real)

copy_list.append([f'ocn.recenter.incr.{str(mem)}.nc',
os.path.join(mem_dir_real, incr_file)])

FileHandler({'mkdir': mem_dir_list}).sync()
FileHandler({'copy': copy_list}).sync()

# Copy the CICE restart files to the member COMROOT directories
if os.getenv('DOIAU') == "YES":
cice_rst_date = self.task_config.MARINE_WINDOW_BEGIN.strftime('%Y%m%d.%H%M%S')
else:
cice_rst_date = cdate.strftime('%Y%m%d.%H%M%S')
mem_dir_list = []
copy_list = []
for mem in range(1, nmem_ens+1):
mem_dir = os.path.join(self.task_config.ROTDIR,
f'enkf{RUN}.{PDYstr}',
f'{cyc}',
f'mem{str(mem).zfill(3)}',
'analysis',
'ice')
mem_dir_real = os.path.realpath(mem_dir)
mem_dir_list.append(mem_dir_real)
copy_list.append([f'ens/cice_model.res.output.{str(mem)}.nc',
os.path.join(mem_dir_real, f'{cice_rst_date}.cice_model_anl.res.nc')])
FileHandler({'mkdir': mem_dir_list}).sync()
FileHandler({'copy': copy_list}).sync()
162 changes: 162 additions & 0 deletions utils/obsproc/InsituAll2ioda.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
#pragma once

#include <cstdlib>
#include <iostream>
#include <map>
#include <netcdf> // NOLINT (using C API)
#include <string>
#include <vector>

#include "eckit/config/LocalConfiguration.h"

#include <Eigen/Dense> // NOLINT

#include "ioda/Group.h"
#include "ioda/ObsGroup.h"

#include "NetCDFToIodaConverter.h"

namespace gdasapp {

class InsituAll2Ioda : public NetCDFToIodaConverter {
public:
explicit InsituAll2Ioda(const eckit::Configuration &fullConfig, const eckit::mpi::Comm &comm)
: NetCDFToIodaConverter(fullConfig, comm) {
ASSERT(fullConfig_.has("variable"));
fullConfig_.get("variable", variable_);
}

// Read NetCDF file and populate data based on YAML configuration
gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
oops::Log::info() << "Processing files provided from ALL in-situ files" << std::endl;

// Abort the case where the 'error ratio' key is not found
ASSERT(fullConfig_.has("error ratio"));

// Get the obs. error ratio from the configuration (unit per day)
float errRatio;
fullConfig_.get("error ratio", errRatio);
// Convert errRatio from meters per day to its unit per second
errRatio /= 86400.0;

// Open the NetCDF file in read-only mode
netCDF::NcFile ncFile(fileName, netCDF::NcFile::read);
oops::Log::info() << "Reading... " << fileName << std::endl;

// Get the number of obs in the file
int nobs = ncFile.getDim("Location").getSize();

// Set the int metadata names
std::vector<std::string> intMetadataNames = {"oceanBasin"};

// Set the float metadata name
std::vector<std::string> floatMetadataNames = {"depth"};

// Create instance of iodaVars object
gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);

// TODO(Mindo): This is incomplete and needed to leverage ioda for the reading
// Check if the MetaData group is null
netCDF::NcGroup metaDataGroup = ncFile.getGroup("MetaData");
if (metaDataGroup.isNull()) {
oops::Log::error() << "Group 'MetaData' not found! Aborting execution..." << std::endl;
std::abort();
}

// Read non-optional metadata: datetime, longitude, latitude and optional: others
netCDF::NcVar latitudeVar = metaDataGroup.getVar("latitude");
std::vector<float> latitudeData(iodaVars.location_);
latitudeVar.getVar(latitudeData.data());

netCDF::NcVar longitudeVar = metaDataGroup.getVar("longitude");
std::vector<float> longitudeData(iodaVars.location_);
longitudeVar.getVar(longitudeData.data());

netCDF::NcVar datetimeVar = metaDataGroup.getVar("dateTime");
std::vector<int64_t> datetimeData(iodaVars.location_);
datetimeVar.getVar(datetimeData.data());
iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; // Applied to All in-situ obs

netCDF::NcVar depthVar = metaDataGroup.getVar("depth");
std::vector<float> depthData(iodaVars.location_, 0); // Initialize with surface value 0

if (!depthVar.isNull()) { // Checking from surface in-situ obs
oops::Log::info() << "Variable 'depth' found and Reading!" << std::endl;
depthVar.getVar(depthData.data());
} else {
oops::Log::warning()
<< "WARNING: no depth found, assuming the observations are at the surface."
<< std::endl;
}

// Save in optional floatMetadata
for (int i = 0; i < iodaVars.location_; i++) {
iodaVars.floatMetadata_.row(i) << depthData[i];
}

netCDF::NcVar oceanbasinVar = metaDataGroup.getVar("oceanBasin");
std::vector<int> oceanbasinData(iodaVars.location_);
oceanbasinVar.getVar(oceanbasinData.data());

// Define and check obs groups
struct { const char* name; netCDF::NcGroup group; } groups[] = {
{"ObsValue", ncFile.getGroup("ObsValue")},
{"ObsError", ncFile.getGroup("ObsError")},
{"PreQC", ncFile.getGroup("PreQC")}
};

// Validate groups and abort if any is missing
for (const auto& g : groups) {
if (g.group.isNull()) {
oops::Log::error() << "Group '" << g.name
<< "' not found! Aborting execution..." << std::endl;
std::abort();
}
}

// Assign validated groups
netCDF::NcGroup& obsvalGroup = groups[0].group;
netCDF::NcGroup& obserrGroup = groups[1].group;
netCDF::NcGroup& preqcGroup = groups[2].group;

// Get obs values, errors and preqc
netCDF::NcVar obsvalVar = obsvalGroup.getVar(variable_);
std::vector<float> obsvalData(iodaVars.location_);
obsvalVar.getVar(obsvalData.data());

netCDF::NcVar obserrVar = obserrGroup.getVar(variable_);
std::vector<float> obserrData(iodaVars.location_);
obserrVar.getVar(obserrData.data());

netCDF::NcVar preqcVar = preqcGroup.getVar(variable_);
std::vector<int> preqcData(iodaVars.location_);
preqcVar.getVar(preqcData.data());

// Update non-optional Eigen arrays
for (int i = 0; i < iodaVars.location_; i++) {
iodaVars.longitude_(i) = longitudeData[i];
iodaVars.latitude_(i) = latitudeData[i];
iodaVars.datetime_(i) = datetimeData[i];
iodaVars.obsVal_(i) = obsvalData[i];
iodaVars.obsError_(i) = obserrData[i];
iodaVars.preQc_(i) = preqcData[i];
// Save in optional intMetadata
iodaVars.intMetadata_.row(i) << oceanbasinData[i];
}

// Extract EpochTime String Format(1970-01-01T00:00:00Z)
std::string extractedDate = iodaVars.referenceDate_.substr(14);

// Redating and adjusting Errors
if (iodaVars.datetime_.size() == 0) {
oops::Log::info() << "datetime_ is empty" << std::endl;
} else {
// Redating and Adjusting Error
iodaVars.reDate(windowBegin_, windowEnd_, extractedDate, errRatio);
}

return iodaVars;
};
}; // class InsituAll2Ioda
} // namespace gdasapp

5 changes: 4 additions & 1 deletion utils/obsproc/Rads2Ioda.h
Original file line number Diff line number Diff line change
Expand Up @@ -132,12 +132,15 @@ namespace gdasapp {
(iodaVars.obsVal_ > -4.0 && iodaVars.obsVal_ < 4.0);
iodaVars.trim(boundsCheck);

// Extract EpochTime String Format(1858-11-17T00:00:00Z)
std::string extractedDate = iodaVars.referenceDate_.substr(14);

// Redating and adjusting Errors
if (iodaVars.datetime_.size() == 0) {
oops::Log::info() << "datetime_ is empty" << std::endl;
} else {
// Redating and Adjusting Error
iodaVars.reDate(windowBegin_, windowEnd_, errRatio);
iodaVars.reDate(windowBegin_, windowEnd_, extractedDate, errRatio);
}

return iodaVars;
Expand Down
4 changes: 4 additions & 0 deletions utils/obsproc/applications/gdas_obsprovider2ioda.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "../IcecAmsr2Ioda.h"
#include "../IcecJpssrr2Ioda.h"
#include "../IcecMirs2Ioda.h"
#include "../InsituAll2ioda.h"
#include "../Rads2Ioda.h"
#include "../RTOFSSalinity.h"
#include "../RTOFSTemperature.h"
Expand Down Expand Up @@ -65,6 +66,9 @@ namespace gdasapp {
} else if (provider == "VIIRSAOD") {
Viirsaod2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else if (provider == "INSITUOBS") {
InsituAll2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else {
oops::Log::info() << "Provider not implemented" << std::endl;
return 1;
Expand Down
Loading

0 comments on commit 1ee5c50

Please sign in to comment.