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

Timestamp:

11/26/18 10:35:26 (5 years ago)

Author:

ymipsl

Message:

Interpolation enhancement :
Domain area can be used to improved global conservation when there is a discrepency between model area and xios computed area.

New domain attribute :

radius : radius af the spherical domain, used to compute area ond the sphere with a normalized radius of 1 (for remapper).

New domain_interpolate attribute :
use_area : remapper will take model computed area in order to perform a global conservation for flux integrated over the cell (mass for example).
In this case the domain attributes "area" and "radius" must be defined on the source or target domain to be taking into account.

YM

File:

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

XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

@@ -304 +304 @@
   CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
   CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
+  CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked);
+  
   long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked];
 
   nSrcLocalUnmasked=0 ;
+  bool hasSrcArea=!domainSrc_->area.isEmpty() && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true  ;
+  double srcAreaFactor ;
+  if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ;
+  
   for (int idx=0 ; idx < nSrcLocal; idx++)
   {
@@ -316 +322 @@
         boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ;
       }
+      if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = domainSrc_->area(idx)*srcAreaFactor ;
       globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ;
       ++nSrcLocalUnmasked ;
     }
   }
-
+ 
 
   int nDstLocalUnmasked = 0 ;
@@ -327 +334 @@
   CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked);
   CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked);
+  CArray<double,1>   areaDstUnmasked(nSrcLocalUnmasked);
+
   long int * globalDstUnmasked = new long int [nDstLocalUnmasked];
 
   nDstLocalUnmasked=0 ;
+  bool hasDstArea=!domainDest_->area.isEmpty() && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ;
+  double dstAreaFactor ;
+  if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ;
   for (int idx=0 ; idx < nDstLocal; idx++)
   {
@@ -339 +351 @@
         boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ;
       }
+      if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = domainDest_->area(idx)*dstAreaFactor ;
       globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ;
       ++nDstLocalUnmasked ;
@@ -344 +357 @@
   }
 
-  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
-  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
+  double* ptAreaSrcUnmasked = NULL ;
+  if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ;
+
+  double* ptAreaDstUnmasked = NULL ;
+  if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ;
+
+  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
+  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
 
   std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity);
@@ -949 +968 @@
     size_t nbIndex=dataOut.numElements() ; 
 
-    for (int idx = 0; idx < nbIndex; ++idx)
-    {
-      if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan
+    if (allMissing.numElements()>0)
+    {
+      for (int idx = 0; idx < nbIndex; ++idx)
+      {
+        if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan
+      }
     }