Changeset 810


Ignore:
Timestamp:
01/07/16 09:18:30 (5 years ago)
Author:
ymipsl
Message:

Improve generation of boundaries for rectilinear grid

YM

File:
1 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/src/node/domain.cpp

    r809 r810  
    503503     bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2; 
    504504     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; 
     505 
     506     // if errorBoundsLon is reasonably small (0.1 x cell size) consider it as closed in longitude 
     507     if (errorBoundsLon < std::abs(lonStepStart)*1e-1 || errorBoundsLon < std::abs(lonStepEnd)*1e-1 ) 
     508     { 
     509       bounds_lon_start= (lon(0) + lon(ni_glo-1)-360)/2 ; 
     510       bounds_lon_end= (lon(0) +360 + lon(ni_glo-1))/2 ; 
     511     } 
     512        
    506513     for(j=0;j<nj;++j) 
    507514       for(i=0;i<ni;++i) 
     
    516523 
    517524    boundsLat.resize(nvertexValue,nj*ni); 
    518     bool isNorthPole, isSouthPole ; 
     525    bool isNorthPole=false ; 
     526    bool isSouthPole=false ; 
    519527    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
    520528    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
    521529 
     530    // lat boundaries beyond pole the assimilate it to pole 
     531    // lat boundarie is relativelly close to pole (0.1 x cell size) assimilate it to pole 
    522532    double latStepStart = lat(1)-lat(0); 
    523533    if (isNorthPole) bounds_lat_start=lat(0); 
    524     else bounds_lat_start=lat(0)-latStepStart/2; 
    525  
     534    else 
     535    { 
     536      bounds_lat_start=lat(0)-latStepStart/2; 
     537      if (bounds_lat_start >= 90 ) bounds_lat_start=90 ; 
     538      else if (bounds_lat_start <= -90 ) bounds_lat_start=-90 ; 
     539      else if (bounds_lat_start <= 90 && bounds_lat_start >= lon(0)) 
     540      { 
     541        if ( 90-bounds_lat_start <= 0.1*std::abs(latStepStart)) bounds_lat_start=90 ; 
     542      } 
     543      else if (bounds_lat_start >= -90 && bounds_lat_start <= lon(0)) 
     544      { 
     545        if ( -90 + bounds_lat_start <= 0.1*std::abs(latStepStart)) bounds_lat_start=-90 ; 
     546      } 
     547    } 
     548     
    526549    double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2); 
    527550    if (isSouthPole) bounds_lat_end=lat(nj_glo-1); 
    528     else bounds_lat_end=lat(nj_glo-1)+latStepEnd/2; 
    529         
    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 ; 
     551    else 
     552    { 
     553      bounds_lat_end=lat(nj_glo-1)+latStepEnd/2; 
     554 
     555      if (bounds_lat_end >= 90 ) bounds_lat_end=90 ; 
     556      else if (bounds_lat_end <= -90 ) bounds_lat_end=-90 ; 
     557      else if (bounds_lat_end <= 90 && bounds_lat_end >= lon(nj_glo-1)) 
     558      { 
     559        if ( 90-bounds_lat_end <= 0.1*std::abs(latStepEnd)) bounds_lat_end=90 ; 
     560      } 
     561      else if (bounds_lat_end >= -90 && bounds_lat_end <= lon(nj_glo-1)) 
     562      { 
     563        if ( -90 + bounds_lat_end <= 0.1*std::abs(latStepEnd)) bounds_lat_end=-90 ; 
     564      } 
     565    }       
    535566       
    536567    for(j=0;j<nj;++j) 
Note: See TracChangeset for help on using the changeset viewer.