diff --git a/src/monio/AtlasReader.cc b/src/monio/AtlasReader.cc index 91befc1..0cca25f 100644 --- a/src/monio/AtlasReader.cc +++ b/src/monio/AtlasReader.cc @@ -14,8 +14,7 @@ #include "UtilsAtlas.h" #include "Monio.h" -monio::AtlasReader::AtlasReader(const eckit::mpi::Comm& mpiCommunicator, - const int mpiRankOwner): +monio::AtlasReader::AtlasReader(const eckit::mpi::Comm& mpiCommunicator, const int mpiRankOwner): mpiCommunicator_(mpiCommunicator), mpiRankOwner_(mpiRankOwner) { oops::Log::debug() << "AtlasReader::AtlasReader()" << std::endl; @@ -116,14 +115,16 @@ void monio::AtlasReader::populateField(atlas::Field& field, oops::Log::debug() << "AtlasReader::populateField()" << std::endl; auto fieldView = atlas::array::make_view(field); // Field with noFirstLevel == true should have been adjusted to have 70 levels. - if (noFirstLevel == true && field.levels() == consts::kVerticalFullSize) { + std::vector fieldShape = field.shape(); + atlas::idx_t numLevels = fieldShape[consts::eVertical]; + if (noFirstLevel == true && numLevels == consts::kVerticalFullSize) { Monio::get().closeFiles(); utils::throwException("AtlasReader::populateField()> Field levels misconfiguration..."); // Only valid case for field with noFirstLevel == true. Field is adjusted to have 70 levels but // read data still has enough to fill 71. } else if (isLfricConvention == true && noFirstLevel == true && - field.levels() == consts::kVerticalHalfSize) { + numLevels == consts::kVerticalHalfSize) { for (int j = 1; j < consts::kVerticalFullSize; ++j) { for (std::size_t i = 0; i < lfricToAtlasMap.size(); ++i) { int index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size()); @@ -140,7 +141,7 @@ void monio::AtlasReader::populateField(atlas::Field& field, // Valid case for isLfricConvention == false and fields noFirstLevel == false. Field is filled // with all available data. } else { - for (int j = 0; j < field.levels(); ++j) { + for (atlas::idx_t j = 0; j < numLevels; ++j) { for (std::size_t i = 0; i < lfricToAtlasMap.size(); ++i) { int index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size()); // Bounds checking @@ -177,15 +178,14 @@ void monio::AtlasReader::populateField(atlas::Field& field, const std::vector& dataVec) { oops::Log::debug() << "AtlasReader::populateField()" << std::endl; - std::vector dimVec = field.shape(); + std::vector fieldShape = field.shape(); if (field.metadata().get("global") == false) { - dimVec[consts::eHorizontal] = utilsatlas::getHorizontalSize(field); + fieldShape[consts::eHorizontal] = utilsatlas::getHorizontalSize(field); } auto fieldView = atlas::array::make_view(field); - int numLevels = field.levels(); - for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) { - for (int j = 0; j < numLevels; ++j) { - int index = i + (j * dimVec[consts::eHorizontal]); + for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) { + for (atlas::idx_t j = 0; j < fieldShape[consts::eVertical]; ++j) { + atlas::idx_t index = i + (j * fieldShape[consts::eHorizontal]); if (std::size_t(index) <= dataVec.size()) { fieldView(i, j) = dataVec[index]; } else { @@ -206,8 +206,9 @@ template void monio::AtlasReader::populateField(atlas::Field& field, atlas::Field monio::AtlasReader::getReadField(atlas::Field& field, const bool noFirstLevel) { + std::vector fieldShape = field.shape(); // Check to ensure field has not been initialised with 71 levels - if (noFirstLevel == true && field.levels() == consts::kVerticalFullSize) { + if (noFirstLevel == true && fieldShape[consts::eVertical] == consts::kVerticalFullSize) { atlas::array::DataType atlasType = field.datatype(); if (atlasType != atlasType.KIND_REAL64 && atlasType != atlasType.KIND_REAL32 && diff --git a/src/monio/AtlasWriter.cc b/src/monio/AtlasWriter.cc index e1266ce..47413b5 100644 --- a/src/monio/AtlasWriter.cc +++ b/src/monio/AtlasWriter.cc @@ -22,8 +22,7 @@ #include "UtilsAtlas.h" #include "Writer.h" -monio::AtlasWriter::AtlasWriter(const eckit::mpi::Comm& mpiCommunicator, - const int mpiRankOwner): +monio::AtlasWriter::AtlasWriter(const eckit::mpi::Comm& mpiCommunicator, const int mpiRankOwner): mpiCommunicator_(mpiCommunicator), mpiRankOwner_(mpiRankOwner) { oops::Log::debug() << "AtlasWriter::AtlasWriter()" << std::endl; @@ -60,11 +59,11 @@ void monio::AtlasWriter::populateFileDataWithField(FileData& fileData, Metadata& metadata = fileData.getMetadata(); Data& data = fileData.getData(); // Create dimensions - std::vector dimVec = field.shape(); + std::vector fieldShape = field.shape(); if (field.metadata().get("global") == false) { - dimVec[0] = utilsatlas::getHorizontalSize(field); + fieldShape[0] = utilsatlas::getHorizontalSize(field); } - for (auto& dimSize : dimVec) { + for (auto& dimSize : fieldShape) { std::string dimName = metadata.getDimensionName(dimSize); if (dimName == consts::kNotFoundError) { dimName = "dim" + std::to_string(dimCount_); @@ -91,7 +90,7 @@ void monio::AtlasWriter::populateFileDataWithField(FileData& fileData, metadata.addVariable(consts::kCoordVarNames[consts::eLongitude], lonVar); metadata.addVariable(consts::kCoordVarNames[consts::eLatitude], latVar); - populateDataWithField(data, field, dimVec); + populateDataWithField(data, field, fieldShape); addGlobalAttributes(metadata, false); } } @@ -155,7 +154,7 @@ void monio::AtlasWriter::populateDataWithField(Data& data, void monio::AtlasWriter::populateDataWithField(Data& data, const atlas::Field& field, - const std::vector dimensions) { + const std::vector dimensions) { oops::Log::debug() << "AtlasWriter::populateDataWithField()" << std::endl; std::shared_ptr dataContainer = nullptr; populateDataContainerWithField(dataContainer, field, dimensions); @@ -170,7 +169,7 @@ void monio::AtlasWriter::populateDataContainerWithField( oops::Log::debug() << "AtlasWriter::populateDataContainerWithField()" << std::endl; if (mpiCommunicator_.rank() == mpiRankOwner_) { atlas::array::DataType atlasType = field.datatype(); - int fieldSize = utilsatlas::getGlobalDataSize(field); + atlas::idx_t fieldSize = utilsatlas::getGlobalDataSize(field); switch (atlasType.kind()) { case atlasType.KIND_INT32: { if (dataContainer == nullptr) { @@ -217,13 +216,13 @@ void monio::AtlasWriter::populateDataContainerWithField( void monio::AtlasWriter::populateDataContainerWithField( std::shared_ptr& dataContainer, const atlas::Field& field, - const std::vector& dimensions) { + const std::vector& dimensions) { oops::Log::debug() << "AtlasWriter::populateDataContainerWithField()" << std::endl; oops::Log::info() << "AtlasWriter::populateDataContainerWithField()" << std::endl; if (mpiCommunicator_.rank() == mpiRankOwner_) { std::string fieldName = field.name(); atlas::array::DataType atlasType = field.datatype(); - int fieldSize = utilsatlas::getGlobalDataSize(field); + atlas::idx_t fieldSize = utilsatlas::getGlobalDataSize(field); switch (atlasType.kind()) { case atlasType.KIND_INT32: { if (dataContainer == nullptr) { @@ -272,7 +271,8 @@ void monio::AtlasWriter::populateDataVec(std::vector& dataVec, const atlas::Field& field, const std::vector& lfricToAtlasMap) { oops::Log::debug() << "AtlasWriter::populateDataVec() " << field.name() << std::endl; - int numLevels = field.levels(); + std::vector fieldShape = field.shape(); + atlas::idx_t numLevels = fieldShape[consts::eVertical]; if ((lfricToAtlasMap.size() * numLevels) != dataVec.size()) { Monio::get().closeFiles(); utils::throwException("AtlasWriter::populateDataVec()> " @@ -280,8 +280,8 @@ void monio::AtlasWriter::populateDataVec(std::vector& dataVec, } auto fieldView = atlas::array::make_view(field); for (std::size_t i = 0; i < lfricToAtlasMap.size(); ++i) { - for (int j = 0; j < numLevels; ++j) { - int index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size()); + for (atlas::idx_t j = 0; j < numLevels; ++j) { + atlas::idx_t index = lfricToAtlasMap[i] + (j * lfricToAtlasMap.size()); dataVec[index] = fieldView(i, j); } } @@ -300,12 +300,12 @@ template void monio::AtlasWriter::populateDataVec(std::vector& dataVec template void monio::AtlasWriter::populateDataVec(std::vector& dataVec, const atlas::Field& field, - const std::vector& dimensions) { + const std::vector& dimensions) { oops::Log::debug() << "AtlasWriter::populateDataVec()" << std::endl; auto fieldView = atlas::array::make_view(field); - for (int i = 0; i < dimensions[consts::eHorizontal]; ++i) { - for (int j = 0; j < dimensions[consts::eVertical]; ++j) { - int index = j + (i * dimensions[consts::eVertical]); + for (atlas::idx_t i = 0; i < dimensions[consts::eHorizontal]; ++i) { + for (atlas::idx_t j = 0; j < dimensions[consts::eVertical]; ++j) { + atlas::idx_t index = j + (i * dimensions[consts::eVertical]); dataVec[index] = fieldView(i, j); } } @@ -313,13 +313,13 @@ void monio::AtlasWriter::populateDataVec(std::vector& dataVec, template void monio::AtlasWriter::populateDataVec(std::vector& dataVec, const atlas::Field& field, - const std::vector& dimensions); + const std::vector& dimensions); template void monio::AtlasWriter::populateDataVec(std::vector& dataVec, const atlas::Field& field, - const std::vector& dimensions); + const std::vector& dimensions); template void monio::AtlasWriter::populateDataVec(std::vector& dataVec, const atlas::Field& field, - const std::vector& dimensions); + const std::vector& dimensions); atlas::Field monio::AtlasWriter::getWriteField(atlas::Field& field, const std::string& writeName, @@ -333,14 +333,15 @@ atlas::Field monio::AtlasWriter::getWriteField(atlas::Field& field, Monio::get().closeFiles(); utils::throwException("AtlasWriter::getWriteField())> Data type not coded for..."); } + std::vector fieldShape = field.shape(); // Erroneous case. For noFirstLevel == true field should have 70 levels - if (noFirstLevel == true && field.levels() == consts::kVerticalFullSize) { + if (noFirstLevel == true && fieldShape[consts::eVertical] == consts::kVerticalFullSize) { Monio::get().closeFiles(); utils::throwException("AtlasWriter::getWriteField()> Field levels misconfiguration..."); } // WARNING - This name-check is an LFRic-Lite specific convention... if (utils::findInVector(consts::kMissingVariableNames, writeName) == false) { - if (noFirstLevel == true && field.levels() == consts::kVerticalHalfSize) { + if (noFirstLevel == true && fieldShape[consts::eVertical] == consts::kVerticalHalfSize) { atlas::util::Config atlasOptions = atlas::option::name(writeName) | atlas::option::global(0) | atlas::option::levels(consts::kVerticalFullSize); @@ -373,14 +374,14 @@ atlas::Field monio::AtlasWriter::copySurfaceLevel(const atlas::Field& inputField atlas::Field copiedField = functionSpace.createField(atlasOptions); auto copiedFieldView = atlas::array::make_view(copiedField); auto inputFieldView = atlas::array::make_view(inputField); - std::vector dimVec = inputField.shape(); - for (int j = 0; j < dimVec[consts::eVertical]; ++j) { - for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) { + std::vector fieldShape = inputField.shape(); + for (atlas::idx_t j = 0; j < fieldShape[consts::eVertical]; ++j) { + for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) { copiedFieldView(i, j + 1) = inputFieldView(i, j); } } // Copy surface level of input field - for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) { + for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) { copiedFieldView(i, 0) = inputFieldView(i, 0); } return copiedField; @@ -400,16 +401,16 @@ void monio::AtlasWriter::addVariableDimensions(const atlas::Field& field, const Metadata& metadata, std::shared_ptr var, const std::string& vertConfigName) { - std::vector dimVec = field.shape(); - if (field.metadata().get("global") == false) { - dimVec[0] = utilsatlas::getHorizontalSize(field); // If so, get the 2D size of the Field + std::vector fieldShape = field.shape(); + if (field.metadata().get("global") == false) { // If so, get the 2D size of the Field + fieldShape[consts::eHorizontal] = utilsatlas::getHorizontalSize(field); } // Reversal of dims required for LFRic files. Currently applied to all output files. - std::reverse(dimVec.begin(), dimVec.end()); - for (auto& dimSize : dimVec) { + std::reverse(fieldShape.begin(), fieldShape.end()); + for (auto& dimSize : fieldShape) { std::string dimName; if (vertConfigName != "" && metadata.isDimDefined(vertConfigName) && - dimSize == metadata.getDimension(vertConfigName)) { + dimSize == metadata.getDimension(vertConfigName)) { dimName = vertConfigName; } else { dimName = metadata.getDimensionName(dimSize); @@ -424,7 +425,7 @@ void monio::AtlasWriter::addGlobalAttributes(Metadata& metadata, const bool isLf // Initialise variables std::string variableConvention = isLfricConvention == true ? consts::kNamingConventions[consts::eLfricConvention] : - consts::kNamingConventions[consts::eJediConvention]; + consts::kNamingConventions[consts::eJediConvention]; // Create attribute objects std::shared_ptr namingAttr = std::make_shared(std::string(consts::kVariableConventionName), diff --git a/src/monio/AtlasWriter.h b/src/monio/AtlasWriter.h index 5f6c77c..6a9a18d 100644 --- a/src/monio/AtlasWriter.h +++ b/src/monio/AtlasWriter.h @@ -78,7 +78,7 @@ class AtlasWriter { /// populateFileDataWithField where metadata are created. void populateDataWithField(Data& data, const atlas::Field& field, - const std::vector dimensions); + const std::vector dimensions); /// \brief Derives the container type and makes the call to populate it. Used where metadata are /// provided and data are written in LFRic order. @@ -101,7 +101,7 @@ class AtlasWriter { /// \brief Iterates through field and populates vector with data from field in Atlas order. template void populateDataVec(std::vector& dataVec, const atlas::Field& field, - const std::vector& dimensions); + const std::vector& dimensions); /// \brief Map JEDI fields back into LFRic function space. atlas::Field getWriteField(atlas::Field& inputField, diff --git a/src/monio/UtilsAtlas.cc b/src/monio/UtilsAtlas.cc index ea4eaac..ae316e7 100644 --- a/src/monio/UtilsAtlas.cc +++ b/src/monio/UtilsAtlas.cc @@ -62,7 +62,7 @@ std::vector getAtlasCoords(const atlas::Field& field) { if (field.metadata().get("global") == false) { auto lonLatField = field.functionspace().lonlat(); auto lonLatView = atlas::array::make_view(lonLatField); - for (int i = 0; i < getHorizontalSize(field); ++i) { + for (atlas::idx_t i = 0; i < getHorizontalSize(field); ++i) { atlasCoords.push_back(atlas::PointLonLat(lonLatView(i, consts::eLongitude), lonLatView(i, consts::eLatitude))); } @@ -131,8 +131,9 @@ std::vector createLfricAtlasMap(const std::vector& a atlas::Field getGlobalField(const atlas::Field& field) { if (field.metadata().get("global") == false) { atlas::array::DataType atlasType = field.datatype(); + std::vector fieldShape = field.shape(); atlas::util::Config atlasOptions = atlas::option::name(field.name()) | - atlas::option::levels(field.levels()) | + atlas::option::levels(fieldShape[consts::eVertical]) | atlas::option::datatype(atlasType) | atlas::option::global(0); if (atlasType != atlasType.KIND_REAL64 && @@ -164,11 +165,11 @@ atlas::FieldSet getGlobalFieldSet(const atlas::FieldSet& fieldSet) { } } -int getHorizontalSize(const atlas::Field& field) { +atlas::idx_t getHorizontalSize(const atlas::Field& field) { atlas::Field ghostField = field.functionspace().ghost(); - int size = 0; + atlas::idx_t size = 0; auto ghostView = atlas::array::make_view(ghostField); - for (int i = ghostField.size() - 1; i >= 0; --i) { + for (atlas::idx_t i = ghostField.size() - 1; i >= 0; --i) { if (ghostView(i) == 0) { size = i + 1; break; @@ -177,10 +178,10 @@ int getHorizontalSize(const atlas::Field& field) { return size; } -int getGlobalDataSize(const atlas::Field& field) { - std::vector dimVec = field.shape(); - int size = 1; - for (const auto& dim : dimVec) { +atlas::idx_t getGlobalDataSize(const atlas::Field& field) { + std::vector fieldShape = field.shape(); + atlas::idx_t size = 1; + for (const auto& dim : fieldShape) { size *= dim; } return size; @@ -216,9 +217,9 @@ bool compareFieldSets(const atlas::FieldSet& aSet, const atlas::FieldSet& bSet) bool compareFields(const atlas::Field& a, const atlas::Field& b) { const auto aView = atlas::array::make_view(a); const auto bView = atlas::array::make_view(b); - std::vector dimVec = a.shape(); - for (int j = 0; j < dimVec[consts::eVertical]; ++j) { - for (int i = 0; i < dimVec[consts::eHorizontal]; ++i) { + std::vector fieldShape = a.shape(); + for (atlas::idx_t j = 0; j < fieldShape[consts::eVertical]; ++j) { + for (atlas::idx_t i = 0; i < fieldShape[consts::eHorizontal]; ++i) { if (aView(i, j) != bView(i, j)) { return false; } diff --git a/src/monio/UtilsAtlas.h b/src/monio/UtilsAtlas.h index 0fe1cb4..128b0a1 100644 --- a/src/monio/UtilsAtlas.h +++ b/src/monio/UtilsAtlas.h @@ -43,8 +43,8 @@ namespace utilsatlas { atlas::Field getGlobalField(const atlas::Field& field); - int getHorizontalSize(const atlas::Field& field); // Just 2D size. Any field. - int getGlobalDataSize(const atlas::Field& field); // Full 3D size of data. Global fields only + atlas::idx_t getHorizontalSize(const atlas::Field& field); // Just 2D size. Any field. + atlas::idx_t getGlobalDataSize(const atlas::Field& field); // Full 3D size of global field. int atlasTypeToMonioEnum(atlas::array::DataType atlasType);