Changeset 810

Timestamp:

01/07/16 09:18:30 (8 years ago)

Author:

ymipsl

Message:

Improve generation of boundaries for rectilinear grid

YM

File:

• XIOS/trunk/src/node/domain.cpp (modified) (2 diffs)

XIOS/trunk/src/node/domain.cpp

@@ -503 +503 @@
      bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2;
      double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start));
-     if (errorBoundsLon > NumTraits<double>::epsilon()) bounds_lon_end=bounds_lon_start+360;
+
+     // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude
+     if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 )
+     {
+       bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ;
+       bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ;
+     }
+       
      for(j=0;j<nj;++j)
        for(i=0;i<ni;++i)
@@ -516 +523 @@
 
     boundsLat.resize(nvertexValue,nj*ni);
-    bool isNorthPole, isSouthPole ;
+    bool isNorthPole=false ;
+    bool isSouthPole=false ;
     if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
     if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
 
+    // lat boundaries beyond pole the assimilate it to pole
+    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole
     double latStepStart = lat(1)-lat(0);
     if (isNorthPole) bounds_lat_start=lat(0);
-    else bounds_lat_start=lat(0)-latStepStart/2;
-
+    else
+    {
+      bounds_lat_start=lat(0)-latStepStart/2;
+      if (bounds_lat_start >= 90 ) bounds_lat_start=90 ;
+      else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ;
+      else if (bounds_lat_start <= 90 && bounds_lat_start >= lon(0))
+      {
+        if ( 90-bounds_lat_start <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ;
+      }
+      else if (bounds_lat_start >= -90 && bounds_lat_start <= lon(0))
+      {
+        if ( -90 + bounds_lat_start <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ;
+      }
+    }
+    
     double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2);
     if (isSouthPole) bounds_lat_end=lat(nj_glo-1);
-    else bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
-       
-// Work arround for small value close to pole, not too good for remapping    
-    if (bounds_lat_start > 90.-1e-3) bounds_lat_start=90 ;
-    if (bounds_lat_start < -90.+1e-3) bounds_lat_start=-90 ;
-    if (bounds_lat_end > 90.-1e-3) bounds_lat_end=90 ;
-    if (bounds_lat_end < -90.+1e-3) bounds_lat_end=-90 ;
+    else
+    {
+      bounds_lat_end=lat(nj_glo-1)+latStepEnd/2;
+
+      if (bounds_lat_end >= 90 ) bounds_lat_end=90 ;
+      else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ;
+      else if (bounds_lat_end <= 90 && bounds_lat_end >= lon(nj_glo-1))
+      {
+        if ( 90-bounds_lat_end <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ;
+      }
+      else if (bounds_lat_end >= -90 && bounds_lat_end <= lon(nj_glo-1))
+      {
+        if ( -90 + bounds_lat_end <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ;
+      }
+    }      
       
     for(j=0;j<nj;++j)