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/node/domain.cpp

    r808 r809  
    457457   } 
    458458 
    459    void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat, 
    460                                               bool isNorthPole, bool isSouthPole) 
     459 
     460 
     461   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g) 
     462   { 
     463          CContext* context = CContext::getCurrent(); 
     464      CContextClient* client = context->client; 
     465          lon_g.resize(ni_glo) ; 
     466          lat_g.resize(nj_glo) ; 
     467           
     468           
     469          int* ibegin_g = new int[client->clientSize] ; 
     470          int* jbegin_g = new int[client->clientSize] ; 
     471          int* ni_g = new int[client->clientSize] ; 
     472          int* nj_g = new int[client->clientSize] ; 
     473          int v ; 
     474          v=ibegin ; 
     475          MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,client->intraComm) ; 
     476          v=jbegin ; 
     477          MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,client->intraComm) ; 
     478          v=ni ; 
     479          MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,client->intraComm) ; 
     480          v=nj ; 
     481          MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,client->intraComm) ; 
     482           
     483          MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,client->intraComm) ; 
     484          MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,client->intraComm) ; 
     485    
     486      delete[] ibegin_g ; 
     487      delete[] jbegin_g ; 
     488      delete[] ni_g ; 
     489      delete[] nj_g ; 
     490   } 
     491 
     492   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat, 
     493                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat) 
    461494   { 
    462495     int i,j,k; 
     496      
    463497     const int nvertexValue = 4; 
    464498     boundsLon.resize(nvertexValue,ni*nj); 
    465499 
    466      if (!lonvalue_rectilinear_read_from_file.isEmpty()) 
    467      { 
    468        double lonStepStart = lonvalue_rectilinear_read_from_file(1)-lonvalue_rectilinear_read_from_file(0); 
    469        bounds_lon_start.setValue(lonvalue_rectilinear_read_from_file(0) - lonStepStart/2); 
    470        double lonStepEnd = (lonvalue_rectilinear_read_from_file(ni_glo-1)-lonvalue_rectilinear_read_from_file(ni_glo-2)); 
    471        bounds_lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1) + lonStepEnd/2); 
    472        double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start)); 
    473        if (errorBoundsLon > NumTraits<double>::epsilon()) bounds_lon_end.setValue(bounds_lon_start+360); 
    474        for(j=0;j<nj;++j) 
    475          for(i=0;i<ni;++i) 
    476          { 
    477            k=j*ni+i; 
    478            boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start 
    479                                                                  : (lonvalue_rectilinear_read_from_file(ibegin + i)+lonvalue_rectilinear_read_from_file(ibegin + i-1))/2; 
    480            boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end 
    481                                                                           : (lonvalue_rectilinear_read_from_file(ibegin + i + 1)+lonvalue_rectilinear_read_from_file(ibegin + i))/2; 
    482          } 
    483      } 
    484      else 
    485      { 
    486        double boundsLonRange = bounds_lon_end - bounds_lon_start; 
    487        double lonStep = boundsLonRange/double(ni_glo.getValue()); 
    488        for(j=0;j<nj;++j) 
    489          for(i=0;i<ni;++i) 
    490          { 
    491            k=j*ni+i; 
    492            boundsLon(0,k) = boundsLon(1,k) = (0 != (ibegin + i)) ? (ibegin + i) * lonStep + bounds_lon_start 
    493                                                                  : bounds_lon_start; 
    494            boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) != ni_glo) ? (ibegin + i +1) * lonStep + bounds_lon_start 
    495                                                                           : bounds_lon_end; 
    496          } 
    497      } 
    498  
    499      boundsLat.resize(nvertexValue,nj*ni); 
    500      if (!latvalue_rectilinear_read_from_file.isEmpty()) 
    501      { 
    502        double latStepStart = latvalue_rectilinear_read_from_file(1)-latvalue_rectilinear_read_from_file(0); 
    503        if (isNorthPole) bounds_lat_start.setValue(latvalue_rectilinear_read_from_file(0) ); 
    504        else bounds_lat_start.setValue(latvalue_rectilinear_read_from_file(0)-latStepStart/2 ); 
    505  
    506  
    507        double latStepEnd = (latvalue_rectilinear_read_from_file(nj_glo-1)-latvalue_rectilinear_read_from_file(nj_glo-2)); 
    508        if (isSouthPole) bounds_lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1)); 
    509        else bounds_lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1)+latStepEnd/2); 
     500     double lonStepStart = lon(1)-lon(0); 
     501     bounds_lon_start=lon(0) - lonStepStart/2; 
     502     double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2); 
     503     bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2; 
     504     double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start)); 
     505     if (errorBoundsLon > NumTraits<double>::epsilon()) bounds_lon_end=bounds_lon_start+360; 
     506     for(j=0;j<nj;++j) 
     507       for(i=0;i<ni;++i) 
     508       { 
     509         k=j*ni+i; 
     510         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start 
     511                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2; 
     512         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end 
     513                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2; 
     514       } 
     515     
     516 
     517    boundsLat.resize(nvertexValue,nj*ni); 
     518    bool isNorthPole, isSouthPole ; 
     519    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
     520    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
     521 
     522    double latStepStart = lat(1)-lat(0); 
     523    if (isNorthPole) bounds_lat_start=lat(0); 
     524    else bounds_lat_start=lat(0)-latStepStart/2; 
     525 
     526    double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2); 
     527    if (isSouthPole) bounds_lat_end=lat(nj_glo-1); 
     528    else bounds_lat_end=lat(nj_glo-1)+latStepEnd/2; 
    510529        
    511        if (bounds_lat_start > 90.-1e-3) bounds_lat_start=90 ; 
    512        if (bounds_lat_start < -90.+1e-3) bounds_lat_start=-90 ; 
    513        if (bounds_lat_end > 90.-1e-3) bounds_lat_end=90 ; 
    514        if (bounds_lat_end < -90.+1e-3) bounds_lat_end=-90 ; 
     530// Work arround for small value close to pole, not too good for remapping     
     531    if (bounds_lat_start > 90.-1e-3) bounds_lat_start=90 ; 
     532    if (bounds_lat_start < -90.+1e-3) bounds_lat_start=-90 ; 
     533    if (bounds_lat_end > 90.-1e-3) bounds_lat_end=90 ; 
     534    if (bounds_lat_end < -90.+1e-3) bounds_lat_end=-90 ; 
    515535       
    516        for(j=0;j<nj;++j) 
    517          for(i=0;i<ni;++i) 
    518          { 
    519            k=j*ni+i; 
    520            boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start 
    521                                                                  : (latvalue_rectilinear_read_from_file(jbegin + j)+latvalue_rectilinear_read_from_file(jbegin + j-1))/2; 
    522            boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end 
    523                                                                  : (latvalue_rectilinear_read_from_file(jbegin + j + 1)+latvalue_rectilinear_read_from_file(jbegin + j))/2; 
    524          } 
    525      } 
    526      else 
    527      { 
    528        double boundsLatRange = bounds_lat_end - bounds_lat_start; 
    529        double latStep = boundsLatRange/double(nj_glo.getValue()); 
    530        double bounds_lat_start_pole = bounds_lat_start; 
    531        double bounds_lat_end_pole   = bounds_lat_end; 
    532        if (isNorthPole) bounds_lat_start_pole = lat_start; 
    533        if (isSouthPole) bounds_lat_end_pole   = lat_end; 
    534  
    535        for(j=0;j<nj;++j) 
    536          for(i=0;i<ni;++i) 
    537          { 
    538            k=j*ni+i; 
    539            boundsLat(1,k) = boundsLat(2,k) = (0 != (jbegin + j)) ? (jbegin + j) * latStep + bounds_lat_start 
    540                                                                  : bounds_lat_start_pole; 
    541            boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) != nj_glo) ? (jbegin + j +1) * latStep + bounds_lat_start 
    542                                                                  : bounds_lat_end_pole; 
    543          } 
    544      } 
    545  
     536    for(j=0;j<nj;++j) 
     537      for(i=0;i<ni;++i) 
     538      { 
     539        k=j*ni+i; 
     540        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start 
     541                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2; 
     542        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end 
     543                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2; 
     544      } 
    546545   } 
    547546 
Note: See TracChangeset for help on using the changeset viewer.