Changeset 846


Ignore:
Timestamp:
04/27/16 17:28:35 (6 years ago)
Author:
ymipsl
Message:

Management of masked cell for conservative interpolation.

YM

Location:
XIOS/trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/src/config/interpolate_domain_attribute.conf

    r689 r846  
    22DECLARE_ATTRIBUTE(StdString, file) 
    33DECLARE_ATTRIBUTE(int, order) 
     4DECLARE_ATTRIBUTE(bool, renormalize) 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

    r833 r846  
    3838  std::vector<double> srcPole(3,0), dstPole(3,0); 
    3939        int orderInterp = interpDomain_->order.getValue(); 
     40  bool renormalize ; 
     41 
     42  if (interpDomain_->renormalize.isEmpty()) renormalize=true; 
     43  else renormalize = interpDomain_->renormalize; 
    4044 
    4145  const double poleValue = 90.0; 
     
    8488    { 
    8589                domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ; 
    86         } 
    87         else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() ) 
     90          } 
     91          else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() ) 
    8892    { 
    8993                lat_g=domainSrc_->latvalue_rectilinear_read_from_file ; 
    9094                lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ; 
    91         } 
    92         else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() && 
     95          } 
     96          else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() && 
    9397                 !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty()) 
    94         { 
    95           double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ; 
    96           for(int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ; 
    97           step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ; 
    98           for(int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ; 
    99         } 
    100         else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ; 
     98          { 
     99            double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ; 
     100            for(int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ; 
     101            step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ; 
     102            for(int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ; 
     103          } 
     104          else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ; 
    101105 
    102106    nVertexSrc = constNVertex; 
     
    228232  Mapper mapper(client->intraComm); 
    229233  mapper.setVerbosity(PROGRESS) ; 
    230   mapper.setSourceMesh(boundsLonSrc.dataFirst(), boundsLatSrc.dataFirst(), nVertexSrc, nSrcLocal, &srcPole[0], globalSrc); 
    231   mapper.setTargetMesh(boundsLonDest.dataFirst(), boundsLatDest.dataFirst(), nVertexDest, nDstLocal, &dstPole[0], globalDst); 
    232   std::vector<double> timings = mapper.computeWeights(orderInterp); 
     234 
     235   
     236  // supress masked data for the source 
     237  int nSrcLocalUnmasked = 0 ; 
     238  for (int idx=0 ; idx < nSrcLocal; idx++) if (domainSrc_-> mask_1d(idx)) ++nSrcLocalUnmasked ; 
     239 
     240  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
     241  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
     242  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked]; 
     243 
     244  nSrcLocalUnmasked=0 ; 
     245  for (int idx=0 ; idx < nSrcLocal; idx++) 
     246  { 
     247    if (domainSrc_-> mask_1d(idx)) 
     248    { 
     249      for(int n=0;n<nVertexSrc;++n) 
     250      { 
     251        boundsLonSrcUnmasked(n,nSrcLocalUnmasked) = boundsLonSrc(n,idx) ; 
     252        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ; 
     253      } 
     254      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ; 
     255      ++nSrcLocalUnmasked ; 
     256    } 
     257  } 
     258 
     259  int nDstLocalUnmasked = 0 ; 
     260  for (int idx=0 ; idx < nDstLocal; idx++) if (domainDest_-> mask_1d(idx)) ++nDstLocalUnmasked ; 
     261 
     262  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked); 
     263  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked); 
     264  long int * globalDstUnmasked = new long int [nDstLocalUnmasked]; 
     265 
     266  nDstLocalUnmasked=0 ; 
     267  for (int idx=0 ; idx < nDstLocal; idx++) 
     268  { 
     269    if (domainDest_-> mask_1d(idx)) 
     270    { 
     271      for(int n=0;n<nVertexDest;++n) 
     272      { 
     273        boundsLonDestUnmasked(n,nDstLocalUnmasked) = boundsLonDest(n,idx) ; 
     274        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ; 
     275      } 
     276      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ; 
     277      ++nDstLocalUnmasked ; 
     278    } 
     279  } 
     280   
     281//  mapper.setSourceMesh(boundsLonSrc.dataFirst(), boundsLatSrc.dataFirst(), nVertexSrc, nSrcLocal, &srcPole[0], globalSrc); 
     282  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
     283//  mapper.setTargetMesh(boundsLonDest.dataFirst(), boundsLatDest.dataFirst(), nVertexDest, nDstLocal, &dstPole[0], globalDst); 
     284  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
     285 
     286  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize); 
    233287 
    234288  std::map<int,std::vector<std::pair<int,double> > > interpMapValue; 
     
    275329 
    276330  delete [] globalSrc; 
     331  delete [] globalSrcUnmasked; 
    277332  delete [] globalDst; 
     333  delete [] globalDstUnmasked; 
     334 
    278335} 
    279336 
Note: See TracChangeset for help on using the changeset viewer.