diff --git a/parm/soca/variational/socaincr2mom6.yaml b/parm/soca/variational/socaincr2mom6.yaml index 75101244e..9fb1b3679 100644 --- a/parm/soca/variational/socaincr2mom6.yaml +++ b/parm/soca/variational/socaincr2mom6.yaml @@ -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 \ No newline at end of file diff --git a/test/soca/testinput/socaincr2mom6.yaml b/test/soca/testinput/socaincr2mom6.yaml index c66417247..61d3ccd9b 100644 --- a/test/soca/testinput/socaincr2mom6.yaml +++ b/test/soca/testinput/socaincr2mom6.yaml @@ -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] @@ -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 diff --git a/utils/gdas_postprocincr.h b/utils/gdas_postprocincr.h index 560751717..699533d88 100644 --- a/utils/gdas_postprocincr.h +++ b/utils/gdas_postprocincr.h @@ -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 @@ -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 @@ -136,13 +141,14 @@ 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; @@ -150,29 +156,35 @@ class PostProcIncr { } // 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; } // ----------------------------------------------------------------------------- @@ -200,12 +212,28 @@ 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_; @@ -213,6 +241,7 @@ class PostProcIncr { oops::Variables socaIncrVar_; bool setToZero_; bool doLVC_; + const soca::State xTraj_; oops::Variables socaZeroIncrVar_; }; diff --git a/utils/gdas_socaincr2mom6.h b/utils/gdas_socaincr2mom6.h index 49f596bb2..c4978e3b3 100644 --- a/utils/gdas_socaincr2mom6.h +++ b/utils/gdas_socaincr2mom6.h @@ -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 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 @@ -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(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: @@ -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