Changeset 809 for XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

Timestamp:

12/14/15 09:21:50 (8 years ago)

Author:

ymipsl

Message:

Fix problems for interpolation onto regular domain.

YM

File:

• XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp (modified) (3 diffs)

XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

@@ -78 +78 @@
     bool isNorthPole = false;
     bool isSouthPole = false;
-
-    if (domainSrc_->latvalue_rectilinear_read_from_file.isEmpty())
-    {
-      if (std::abs(poleValue - std::abs(domainSrc_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true;
-      if (std::abs(poleValue - std::abs(domainSrc_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true;
-    }
-    else
-    {
-      if (std::abs(poleValue - std::abs(domainSrc_->latvalue_rectilinear_read_from_file(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
-      if (std::abs(poleValue - std::abs(domainSrc_->latvalue_rectilinear_read_from_file(domainSrc_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
-    }
+    CArray<double,1> lon_g ;
+    CArray<double,1> lat_g ;
+
+    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty())
+    {
+                domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ;
+        }
+        else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() )
+    {
+                lat_g=domainSrc_->latvalue_rectilinear_read_from_file ;
+                lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ;
+        }
+        else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() &&
+                 !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty())
+        {
+          double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ;
+          for(int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ;
+          step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ; 
+          for(int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ;
+        }
+        else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
 
     nVertexSrc = constNVertex;
-    domainSrc_->fillInRectilinearBoundLonLat(boundsLonSrc, boundsLatSrc, isNorthPole, isSouthPole);
+    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc);
   }
 
@@ -132 +142 @@
     if (std::abs(poleValue - std::abs(domainDest_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true;
     if (std::abs(poleValue - std::abs(domainDest_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true;
+
+    CArray<double,1> lon_g ;
+    CArray<double,1> lat_g ;
+
+    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty())
+    {
+                domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ;
+        }
+        else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() )
+    {
+                lat_g=domainDest_->latvalue_rectilinear_read_from_file ;
+                lon_g=domainDest_->lonvalue_rectilinear_read_from_file ;
+        }
+        else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() &&
+                 !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty())
+        {
+          double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
+          for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
+          step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ; 
+          for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
+        }
+        else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
+    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
+    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
+
+    
+    
+    
     if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
     {
@@ -155 +193 @@
     // Ok, fill in boundary values for rectangular domain
     nVertexDest = constNVertex;
-    domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest, isNorthPole, isSouthPole);
+    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest);
   }