MGBF added by Ting Lei

cmake/compiler_flags_Intel_Fortran.cmake

@@ -0,0 +1,41 @@
+# (C) Copyright 2019 UCAR
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+####################################################################
+# FLAGS COMMON TO ALL BUILD TYPES
+####################################################################
+
+####################################################################
+# RELEASE FLAGS
+####################################################################
+
+#cltorg set( CMAKE_Fortran_FLAGS_RELEASE "-O3 -ip -unroll -inline -no-heap-arrays" )
+set( CMAKE_Fortran_FLAGS_RELEASE "-O3 -ip -unroll -inline -no-heap-arrays -fpp" )
+
+####################################################################
+# DEBUG FLAGS
+####################################################################
+
+#cltorg set( CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g -check bounds -traceback -warn -heap-arrays -fpe-all=0 -fpe:0 -check all" )
+set( CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g -check bounds -traceback -warn -heap-arrays -fpe-all=0 -fpe:0 -check all -fpp" )
+
+####################################################################
+# BIT REPRODUCIBLE FLAGS
+####################################################################
+
+set( CMAKE_Fortran_FLAGS_BIT     "-O2 -ip -ipo -unroll -inline -no-heap-arrays" )
+
+####################################################################
+# LINK FLAGS
+####################################################################
+
+set( CMAKE_Fortran_LINK_FLAGS    "" )
+
+####################################################################
+
+# Meaning of flags
+# ----------------
+# todo
+

src/saber/CMakeLists.txt

@@ -4,7 +4,7 @@
 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 
 # Build list of subdirs with files to add
-set( _subdirs bifourier blocks bump diffusion fastlam generic gsi interpolation oops spectralb torchbalance util vader )
+set( _subdirs bifourier blocks bump diffusion fastlam generic gsi interpolation oops spectralb torchbalance util vader mgbf)
 foreach( _subdir IN LISTS _subdirs )
     add_subdirectory( ${_subdir} )
     list( TRANSFORM ${_subdir}_src_files PREPEND ${_subdir}/ )
@@ -40,10 +40,14 @@ endif()
 if( gsibec_FOUND )
     target_link_libraries( ${PROJECT_NAME} PUBLIC gsibec )
     target_compile_definitions( ${PROJECT_NAME} PUBLIC GSIBEC_FOUND )
-    if ( ip_FOUND )
-        target_link_libraries( ${PROJECT_NAME} PUBLIC ip::ip_d )
+    if ( sp_FOUND )
+        target_link_libraries( ${PROJECT_NAME} PUBLIC sp::sp_d )
     endif()
 endif()
+if( MGBFLIB_FOUND EQUAL 9999 )
+    target_link_libraries( ${PROJECT_NAME} PUBLIC mgbf_lib )
+    target_compile_definitions( ${PROJECT_NAME} PUBLIC MGBF_FOUND)
+endif()
 if ( ENABLE_OFFLINE_CODECOV )
     target_link_libraries( ${PROJECT_NAME} PUBLIC gcov )
 endif()

src/saber/interpolation/Interpolation.cc

@@ -11,6 +11,7 @@
 
 #include "oops/util/FieldSetOperations.h"
 #include "oops/util/Logger.h"
+#include "oops/util/missingValues.h"
 
 namespace saber {
 namespace interpolation {
@@ -19,6 +20,90 @@ namespace interpolation {
 
 static SaberOuterBlockMaker<Interpolation> makerInterpolation_("interpolation");
 
+namespace {
+void fillMissingValuesNearest(const atlas::FieldSet & sourceFieldSet,
+                              atlas::FieldSet & targetFieldSet,
+                              const oops::Variables & vars,
+                              const atlas::FunctionSpace & sourceFs,
+                              const atlas::FunctionSpace & targetFs) {
+  if (vars.size() == 0) {
+    oops::Log::trace()
+        << Interpolation::classname()
+        << "fillMissingValuesNearest: no variables to process"
+        << std::endl;
+    return;
+  }
+
+  const auto src_lonlat = atlas::array::make_view<double, 2>(sourceFs.lonlat());
+  const auto src_ghost = atlas::array::make_view<int, 1>(sourceFs.ghost());
+  std::vector<double> lons;
+  std::vector<double> lats;
+  std::vector<atlas::idx_t> indices;
+  lons.reserve(src_lonlat.shape(0));
+  lats.reserve(src_lonlat.shape(0));
+  indices.reserve(src_lonlat.shape(0));
+  for (atlas::idx_t jj = 0; jj < src_lonlat.shape(0); ++jj) {
+    if (src_ghost(jj) == 0) {
+      lons.push_back(src_lonlat(jj, 0));
+      lats.push_back(src_lonlat(jj, 1));
+      indices.push_back(jj);
+    }
+  }
+  if (indices.empty()) {
+    oops::Log::trace()
+          << Interpolation::classname()
+          << "fillMissingValuesNearest: no owned source points"
+          << std::endl;
+    return;
+  }
+
+  const atlas::Geometry earth(atlas::util::Earth::radius());
+  atlas::util::IndexKDTree2D tree(earth);
+  tree.build(lons, lats, indices);
+
+  const auto tgt_lonlat = atlas::array::make_view<double, 2>(targetFs.lonlat());
+  const auto tgt_ghost = atlas::array::make_view<int, 1>(targetFs.ghost());
+  const double missing = util::missingValue<double>();
+
+
+  for (const auto & var : vars) {
+    if (!targetFieldSet.has(var.name()) || !sourceFieldSet.has(var.name())) {
+      oops::Log::trace()
+            << Interpolation::classname()
+            << "fillMissingValuesNearest: skipping var (missing in fset)" <<var.name()
+            << std::endl;
+      continue;
+    }
+    auto tgt_view = atlas::array::make_view<double, 2>(targetFieldSet[var.name()]);
+    const auto src_view = atlas::array::make_view<double, 2>(
+        sourceFieldSet.field(var.name()));
+
+    for (atlas::idx_t jloc = 0; jloc < tgt_view.shape(0); ++jloc) {
+      if (tgt_ghost(jloc) != 0) {
+        continue;
+      }
+      bool has_missing = false;
+      for (atlas::idx_t jlev = 0; jlev < tgt_view.shape(1); ++jlev) {
+        if (tgt_view(jloc, jlev) == missing) {
+          has_missing = true;
+          break;
+        }
+      }
+      if (!has_missing) {
+        continue;
+      }
+      atlas::PointLonLat pll(tgt_lonlat(jloc, 0), tgt_lonlat(jloc, 1));
+      const auto item = tree.closestPoint(pll);
+      const atlas::idx_t src_index = item.payload();
+      for (atlas::idx_t jlev = 0; jlev < tgt_view.shape(1); ++jlev) {
+        if (tgt_view(jloc, jlev) == missing) {
+          tgt_view(jloc, jlev) = src_view(src_index, jlev);
+       }
+      }
+    }
+  }
+}
+}  // namespace
 // -----------------------------------------------------------------------------
 
 Interpolation::Interpolation(const oops::GeometryData & outerGeometryData,
@@ -183,6 +268,13 @@ void Interpolation::inverseMultiply(oops::FieldSet3D & fieldSet) const {
     }
     inverseRegionalInterp_->execute(sourceFieldSet, targetFieldSet);
   }
+  if (params_.fillMissingValues.value()) {
+    fillMissingValuesNearest(sourceFieldSet, targetFieldSet, invVars,
+                             outerGeometryData_.functionSpace(),
+                             innerGeomData_->functionSpace());
+    // Update halos after filling missing values so boundary points are consistent
+    targetFieldSet.haloExchange();
+  }
 
   // Reset
   fieldSet.fieldSet() = targetFieldSet;

src/saber/interpolation/Interpolation.h

@@ -41,6 +41,7 @@ class InterpolationParameters : public SaberBlockParametersBase {
     eckit::LocalConfiguration(), this};
   oops::Parameter<eckit::LocalConfiguration> inverseInterpConf{"inverse interpolator",
     eckit::LocalConfiguration(), this};
+  oops::Parameter<bool> fillMissingValues{"fill state missing values", false, this};
   oops::Variables mandatoryActiveVars() const override {return oops::Variables();}
 };
 

src/saber/mgbf/CMakeLists.txt

@@ -0,0 +1,58 @@
+# (C) Copyright 2022 United States Government as represented by the Administrator of the National
+#     Aeronautics and Space Administration
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+file(GLOB jbfiles mgbf_lib/*.f90)
+   message(STATUS "thinkdeb-2  "  ${jbfiles}  )
+set (jbfilenames "")
+foreach ( _fname ${jbfiles} ) 
+   get_filename_component( basefilename ${_fname} NAME )
+   list ( APPEND jbfilenames  mgbf_lib/${basefilename} )
+   message(STATUS "thinkdeb-1  " ${basefilename})
+   message(STATUS "thinkdeb0  " ${jbfilenames})
+endforeach ()
+message(STATUS "thinkdeb  " ${jbfilenames})
+#set (jbfilenames "mgbf_lib/jp_pbfil.f90" )
+set (build_saber_mgbf 1)
+if( build_saber_mgbf )
+    list(APPEND mgbf_src_files_list
+
+    # Covariance block
+    covariance/MGBF_Covariance.h
+    covariance/MGBF_Covariance.cc
+    covariance/MGBF_Covariance.interface.F90
+    covariance/MGBF_Covariance.interface.h
+    covariance/mgbf_covariance_mod.f90
+
+#clth    # Grid
+#     covariance/mgbf_Grid.h
+#     covariance/mgbf_Grid.cc
+    # Interpolation block
+#     covariance/mgbf_Interpolation.h
+#     covariance/mgbf_Interpolation.cc
+#    interpolation/MGBF_Interpolation.h
+
+    # Unstructured interpolation code ported from oops (until new interp code can be used)
+#    interpolation/unstructured_interp/saber_unstructured_interpolation_mod.F90
+#    interpolation/unstructured_interp/UnstructuredInterpolation.cc
+#    interpolation/unstructured_interp/UnstructuredInterpolation.h
+#    interpolation/unstructured_interp/UnstructuredInterpolation.interface.F90
+#    interpolation/unstructured_interp/UnstructuredInterpolation.interface.h
+
+    # Utilities
+#    utils/mgbf_utils_mod.f90
+
+    )
+endif()
+#clt find_package(mgbf_lib  REQUIRED )
+message (STATUS "thinkdeb1 " ${mgbf_src_files_list} )
+
+set( mgbf_src_files
+
+${mgbf_src_files_list}
+${jbfilenames}
+
+PARENT_SCOPE
+)
+ message (STATUS "thinkdeb2.4"   ${mgbf_src_files} )

src/saber/mgbf/covariance/MGBF_Covariance.cc

@@ -0,0 +1,36 @@
+/*
+ * (C) Copyright 2022 United States Government as represented by the Administrator of the National
+ *     Aeronautics and Space Administration
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#include "saber/mgbf/covariance/MGBF_Covariance.h"
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "atlas/field.h"
+#include "atlas/functionspace.h"
+#include "atlas/grid.h"
+#include "atlas/library.h"
+#include "atlas/runtime/Log.h"
+
+#include "oops/base/Variables.h"
+#include "oops/util/abor1_cpp.h"
+#include "oops/util/Logger.h"
+#include "oops/util/Timer.h"
+
+#include "saber/mgbf/covariance/MGBF_Covariance.interface.h"
+
+namespace saber {
+namespace mgbf {
+
+// -------------------------------------------------------------------------------------------------
+
+static SaberCentralBlockMaker<MGBF_Covariance> makerCovariance_("MGBF_covariance");
+
+}  // namespace MGBF 
+}  // namespace saber