Ignore:
Timestamp:
12/14/15 09:21:50 (7 years ago)
Author:
ymipsl
Message:

Fix problems for interpolation onto regular domain.

YM

File:
1 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

    r808 r809  
    7878    bool isNorthPole = false; 
    7979    bool isSouthPole = false; 
    80  
    81     if (domainSrc_->latvalue_rectilinear_read_from_file.isEmpty()) 
    82     { 
    83       if (std::abs(poleValue - std::abs(domainSrc_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true; 
    84       if (std::abs(poleValue - std::abs(domainSrc_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true; 
    85     } 
    86     else 
    87     { 
    88       if (std::abs(poleValue - std::abs(domainSrc_->latvalue_rectilinear_read_from_file(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
    89       if (std::abs(poleValue - std::abs(domainSrc_->latvalue_rectilinear_read_from_file(domainSrc_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
    90     } 
     80    CArray<double,1> lon_g ; 
     81    CArray<double,1> lat_g ; 
     82 
     83    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty()) 
     84    { 
     85                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() ) 
     88    { 
     89                lat_g=domainSrc_->latvalue_rectilinear_read_from_file ; 
     90                lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ; 
     91        } 
     92        else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() && 
     93                 !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") ; 
    91101 
    92102    nVertexSrc = constNVertex; 
    93     domainSrc_->fillInRectilinearBoundLonLat(boundsLonSrc, boundsLatSrc, isNorthPole, isSouthPole); 
     103    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc); 
    94104  } 
    95105 
     
    132142    if (std::abs(poleValue - std::abs(domainDest_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true; 
    133143    if (std::abs(poleValue - std::abs(domainDest_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true; 
     144 
     145    CArray<double,1> lon_g ; 
     146    CArray<double,1> lat_g ; 
     147 
     148    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty()) 
     149    { 
     150                domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ; 
     151        } 
     152        else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() ) 
     153    { 
     154                lat_g=domainDest_->latvalue_rectilinear_read_from_file ; 
     155                lon_g=domainDest_->lonvalue_rectilinear_read_from_file ; 
     156        } 
     157        else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() && 
     158                 !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty()) 
     159        { 
     160          double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ; 
     161          for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ; 
     162          step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;  
     163          for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ; 
     164        } 
     165        else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ; 
     166    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
     167    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
     168 
     169     
     170     
     171     
    134172    if (isNorthPole && (0 == domainDest_->jbegin.getValue())) 
    135173    { 
     
    155193    // Ok, fill in boundary values for rectangular domain 
    156194    nVertexDest = constNVertex; 
    157     domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest, isNorthPole, isSouthPole); 
     195    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest); 
    158196  } 
    159197 
Note: See TracChangeset for help on using the changeset viewer.