Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
@guillaumevernieres
guillaumevernieres committed Jul 31, 2023
1 parent 5176f26 commit f5d06ff
Show file tree
Hide file tree
Showing 4 changed files with 96 additions and 152 deletions.
3 changes: 2 additions & 1 deletion parm/soca/variational/socaincr2mom6.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,9 @@ soca increment:
ocn_filename: 'ocn.3dvarfgat_pseudo.incr.{{ATM_WINDOW_MIDDLE}}.nc'
read_from_file: 1

mom6 iau increment:
output increment:
datadir: ./
date: '{{ATM_WINDOW_BEGIN}}'
exp: mom6_iau
type: incr
output file: inc.nc
27 changes: 23 additions & 4 deletions test/soca/testinput/socaincr2mom6.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ date: 2018-04-15T09:00:00Z

layers variable: [hocn]

increment variables: [tocn, socn, uocn, vocn, ssh]
increment variables: [tocn, socn, uocn, vocn, ssh, hocn]

set increment variables to zero: [uocn, vocn, ssh]

Expand All @@ -16,14 +16,33 @@ vertical geometry:
ocn_filename: MOM.res.nc
read_from_file: 1

soca increment:
date: 2018-04-15T09:00:00Z
soca increments:
- date: 2018-04-15T09:00:00Z
basename: ./Data/
ocn_filename: 'ocn.3dvarfgat_pseudo.incr.2018-04-15T12:00:00Z.nc'
read_from_file: 1
- date: 2018-04-15T09:00:00Z
basename: ./Data/
ocn_filename: 'ocn.3dvarfgat_pseudo.incr.2018-04-15T12:00:00Z.nc'
read_from_file: 1

linear variable change:
linear variable changes:
- linear variable change name: BkgErrFILT
ocean_depth_min: 500 # zero where ocean is shallower than 500m
rescale_bkgerr: 1.0 # rescale perturbation
efold_z: 1500.0 # Apply exponential decay
- linear variable change name: BalanceSOCA
trajectory:
state variables: [tocn, socn, uocn, vocn, ssh, hocn, layer_depth, mld]
date: 2018-04-15T09:00:00Z
basename: ./INPUT/
ocn_filename: MOM.res.nc
read_from_file: 1

mom6 iau increment:
output increment:
datadir: ./
date: 2018-04-15T09:00:00Z
exp: mom6_iau
type: incr
output file: inc.nc
83 changes: 56 additions & 27 deletions utils/gdas_postprocincr.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,14 @@ namespace gdasapp {
class PostProcIncr {
public:
// Constructor
PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry& geom)
PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry& geom,
const eckit::mpi::Comm & comm)
: dt_(getDate(fullConfig)),
layerVar_(getLayerVar(fullConfig)),
geom_(geom),
Layers_(getLayerThickness(fullConfig, geom)){
Layers_(getLayerThickness(fullConfig, geom)),
xTraj_(getTraj(fullConfig, geom)),
comm_(comm) {
oops::Log::info() << "Date: " << std::endl << dt_ << std::endl;

// Increment variables
Expand All @@ -38,11 +41,13 @@ class PostProcIncr {
socaIncrVar_ = socaIncrVar;

// Input increment configuration
eckit::LocalConfiguration inputIncrConfig(fullConfig, "soca increment");
eckit::LocalConfiguration inputIncrConfig = fullConfig.getSubConfiguration("soca increments");
inputIncrConfig_ = inputIncrConfig;

std::cout << "------------------------------------------" << std::endl;
std::cout << inputIncrConfig_ << std::endl;
// Output incrememnt configuration
eckit::LocalConfiguration outputIncrConfig(fullConfig, "mom6 iau increment");
eckit::LocalConfiguration outputIncrConfig(fullConfig, "output increment");
outputIncrConfig_ = outputIncrConfig;

// Variables that should be set to 0
Expand Down Expand Up @@ -136,43 +141,50 @@ class PostProcIncr {
return socaIncr;
}
oops::Log::info() << "====== applying specified change of variables" << std::endl;
const eckit::LocalConfiguration trajConfig(lvcConfig_, "trajectory");
soca::State xTraj(this->geom_, trajConfig);
//const eckit::LocalConfiguration trajConfig(lvcConfig_, "trajectory");
//soca::State xTraj(this->geom_, trajConfig);
soca::LinearVariableChangeParameters params;
params.deserialize(lvcConfig_);
oops::Log::info() << params << std::endl;
soca::LinearVariableChange lvc(this->geom_, params);
lvc.changeVarTraj(xTraj, socaIncrVar_);
oops::Log::info() << "traj: " << xTraj_ << std::endl;
lvc.changeVarTraj(xTraj_, socaIncrVar_);
lvc.changeVarTL(socaIncr, socaIncrVar_);
oops::Log::info() << "incr after var change: " << std::endl << socaIncr << std::endl;

return socaIncr;
}

// Save increment
void save(soca::Increment socaIncr) {
int save(soca::Increment socaIncr) {
oops::Log::info() << "==========================================" << std::endl;
oops::Log::info() << "-------------------- save increment: " << std::endl;
socaIncr.write(outputIncrConfig_);

std::string outputFileName;
outputIncrConfig_.get("output file", outputFileName);

std::string datadir;
outputIncrConfig_.get("datadir", datadir);
std::filesystem::path pathToResolve = datadir;
std::string exp;
outputIncrConfig_.get("exp", exp);
std::string outputType;
outputIncrConfig_.get("type", outputType);
std::string incrFname = std::filesystem::canonical(pathToResolve);
incrFname += "/ocn." + exp + "." + outputType + "." + dt_.toString() + ".nc";
const char* charPtr = incrFname.c_str();
const char* charPtrOut = outputFileName.c_str();
oops::Log::info() << "rename: " << incrFname << " to " << outputFileName << std::endl;
int result = std::rename(charPtr, charPtrOut);

//return;
// wait for everybody to be done
comm_.barrier();

// Change soca standard output name to something specidfied in the config
int result = 0;
if ( comm_.rank() == 0 ) {
std::string outputFileName;
outputIncrConfig_.get("output file", outputFileName);

std::string datadir;
outputIncrConfig_.get("datadir", datadir);
std::filesystem::path pathToResolve = datadir;
std::string exp;
outputIncrConfig_.get("exp", exp);
std::string outputType;
outputIncrConfig_.get("type", outputType);
std::string incrFname = std::filesystem::canonical(pathToResolve);
incrFname += "/ocn." + exp + "." + outputType + "." + dt_.toString() + ".nc";
const char* charPtr = incrFname.c_str();
const char* charPtrOut = outputFileName.c_str();
oops::Log::info() << "rename: " << incrFname << " to " << outputFileName << std::endl;
result = std::rename(charPtr, charPtrOut);
}
return result;
}

// -----------------------------------------------------------------------------
Expand Down Expand Up @@ -200,19 +212,36 @@ class PostProcIncr {
oops::Log::info() << "layerThick: " << std::endl << layerThick << std::endl;
return layerThick;
}
// Initialize the trajectory
soca::State getTraj(const eckit::Configuration& fullConfig,
const soca::Geometry& geom) const {
if ( fullConfig.has("linear variable change") ) {
const eckit::LocalConfiguration trajConfig(fullConfig, "linear variable change.trajectory");
soca::State traj(geom, trajConfig);
oops::Log::info() << "traj:" << traj << std::endl;
return traj;
} else {
oops::Variables tmpVar(fullConfig, "layers variable");
util::DateTime tmpDate(getDate(fullConfig));
soca::State traj(geom, tmpVar, tmpDate);
return traj;
}
}

private:
public:
util::DateTime dt_; // valid date of increment
oops::Variables layerVar_; // layer variable
const soca::Increment Layers_; // layer thicknesses
const soca::Geometry & geom_;
const eckit::mpi::Comm & comm_;
eckit::LocalConfiguration inputIncrConfig_;
eckit::LocalConfiguration outputIncrConfig_;
eckit::LocalConfiguration zeroIncrConfig_;
eckit::LocalConfiguration lvcConfig_;
oops::Variables socaIncrVar_;
bool setToZero_;
bool doLVC_;
const soca::State xTraj_;
oops::Variables socaZeroIncrVar_;
};

Expand Down
135 changes: 15 additions & 120 deletions utils/gdas_socaincr2mom6.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,19 @@ namespace gdasapp {
oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl;
const soca::Geometry geom(geomConfig, this->getComm());


// Initialize the post processing
PostProcIncr postProcIncr(fullConfig, geom, this->getComm());

oops::Log::info() << "soca increments: " << std::endl << postProcIncr.inputIncrConfig_ << std::endl;

std::vector<eckit::LocalConfiguration> increments;
fullConfig.get("soca increments", increments);
oops::Log::info() << "ooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << std::endl;
oops::Log::info() << increments[0] << std::endl;
oops::Log::info() << "ooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << std::endl;

// At the very minimum, we run this script to append the layers state, so do that!
PostProcIncr postProcIncr(fullConfig, geom);
soca::Increment incr = postProcIncr.appendLayer();

// Zero out specified fields
Expand All @@ -46,101 +57,9 @@ namespace gdasapp {
incr = postProcIncr.applyLinVarChange(incr);

// Save final increment
postProcIncr.save(incr);
/*
// get date
std::string strdt;
fullConfig.get("date", strdt);
util::DateTime dt = util::DateTime(strdt);
// Read the vertical geometry from the background (layer thicknesses)
soca::Increment layerThick = readLayerThickness(fullConfig, geom);
// Setup the soca increment
// get the increment variables
oops::Variables socaIncrVar(fullConfig, "increment variables");
ASSERT(socaIncrVar.size() >= 1);
// read the soca increment
soca::Increment socaIncr(geom, socaIncrVar, dt);
const eckit::LocalConfiguration socaIncrConfig(fullConfig, "soca increment");
socaIncr.read(socaIncrConfig);
oops::Log::info() << "socaIncr: " << std::endl << socaIncr << std::endl;
/// Create the MOM6 IAU increment
// concatenate variables
oops::Variables layerVar = getLayerVarName(fullConfig);
oops::Variables mom6IauVar(socaIncrVar);
mom6IauVar += layerVar;
oops::Log::info() << "mom6IauVar: " << std::endl << mom6IauVar << std::endl;
// append layer variable to soca increment
atlas::FieldSet socaIncrFs;
socaIncr.toFieldSet(socaIncrFs);
socaIncr.updateFields(mom6IauVar);
oops::Log::info() << "MOM6 incr: " << std::endl << socaIncr << std::endl;
// pad layer increment with zeros
atlas::FieldSet socaLayerThickFs;
layerThick.toFieldSet(socaLayerThickFs);
layerThick.updateFields(mom6IauVar);
oops::Log::info() << "h: " << std::endl << layerThick << std::endl;
// append layer thinckness to increment
socaIncr += layerThick;
oops::Log::info() << "MOM6 IAU increment: " << std::endl << socaIncr << std::endl;
// set specified increment variables to 0.0
oops::Variables socaZeroIncrVar(fullConfig, "set increment variables to zero");
std::cout << socaZeroIncrVar << std::endl;
for (auto & field : socaIncrFs) {
// only works if rank is 2
ASSERT(field.rank() == 2);
// Set data to zero
if (socaZeroIncrVar.has(field.name())) {
std::cout << "filed.name(): " << field.name() << std::endl;
auto view = atlas::array::make_view<double, 2>(field);
view.assign(0.0);
}
}
socaIncr.fromFieldSet(socaIncrFs);
// Apply change of variable if in the configuration
if ( fullConfig.has("linear variable change") ) {
const eckit::LocalConfiguration trajConfig(fullConfig, "linear variable change.trajectory");
soca::State xTraj(geom, trajConfig);
soca::LinearVariableChangeParameters params;
const eckit::LocalConfiguration lvcConfig(fullConfig, "linear variable change");
params.deserialize(lvcConfig);
oops::Log::info() << params << std::endl;
soca::LinearVariableChange lvc(geom, params);
lvc.changeVarTraj(xTraj, socaIncrVar);
lvc.changeVarTL(socaIncr, socaIncrVar);
oops::Log::info() << "incr after var change: " << std::endl << socaIncr << std::endl;
}
// Save MOM6 IAU Increment
const eckit::LocalConfiguration mom6IauConfig(fullConfig, "mom6 iau increment");
socaIncr.write(mom6IauConfig);
// TODO: the "checkpoint" script expects the ocean increment output to
// be in "inc.nc". Remove what's below, eventually
std::string datadir;
mom6IauConfig.get("datadir", datadir);
std::filesystem::path pathToResolve = datadir;
std::string exp;
mom6IauConfig.get("exp", exp);
std::string outputType;
mom6IauConfig.get("type", outputType);
std::string incrFname = std::filesystem::canonical(pathToResolve);
incrFname += "/ocn." + exp + "." + outputType + "." + dt.toString() + ".nc";
const char* charPtr = incrFname.c_str();
oops::Log::info() << "rename: " << incrFname << " to " << "inc.nc" << std::endl;
int result = std::rename(charPtr, "inc.nc");
*/
return 0; //result;
int result = postProcIncr.save(incr);

return result;
}
// -----------------------------------------------------------------------------
private:
Expand All @@ -150,30 +69,6 @@ namespace gdasapp {
std::string appname() const {
return "gdasapp::SocaIncr2Mom6";
}
// -----------------------------------------------------------------------------
// Read layer thickness from the relevant background
soca::Increment readLayerThickness(const eckit::Configuration& fullConfig,
const soca::Geometry& geom) const {

soca::Increment layerThick(geom, getLayerVarName(fullConfig), getDate(fullConfig));
const eckit::LocalConfiguration vertGeomConfig(fullConfig, "vertical geometry");
layerThick.read(vertGeomConfig);
oops::Log::info() << "layerThick: " << std::endl << layerThick << std::endl;

return layerThick;
}
// -----------------------------------------------------------------------------
util::DateTime getDate(const eckit::Configuration& fullConfig) const {
std::string strdt;
fullConfig.get("date", strdt);
return util::DateTime(strdt);
}
// -----------------------------------------------------------------------------
oops::Variables getLayerVarName(const eckit::Configuration& fullConfig) const {
oops::Variables layerVar(fullConfig, "layers variable");
ASSERT(layerVar.size() == 1);
return layerVar;
}
};

} // namespace gdasapp

0 comments on commit f5d06ff

Please sign in to comment.