source: XIOS3/trunk/src/io/nc4_data_output.cpp @ 2384

Last change on this file since 2384 was 2384, checked in by jderouillat, 2 years ago

Replace OutputDomainName? by a CDomain pointer in the hash table which manage element names. Replace unordered_multimap by multimap

  • Property copyright set to
    Software name : XIOS (Xml I/O Server)
    http://forge.ipsl.jussieu.fr/ioserver
    Creation date : January 2009
    Licence : CeCCIL version2
    see license file in root directory : Licence_CeCILL_V2-en.txt
    or http://www.cecill.info/licences/Licence_CeCILL_V2-en.html
    Holder : CEA/LSCE (Laboratoire des Sciences du CLimat et de l'Environnement)
    CNRS/IPSL (Institut Pierre Simon Laplace)
    Project Manager : Yann Meurdesoif
    yann.meurdesoif@cea.fr
  • Property svn:executable set to *
File size: 129.4 KB
RevLine 
[219]1#include "nc4_data_output.hpp"
2
[352]3#include "attribute_template.hpp"
4#include "group_template.hpp"
[219]5
6#include "file.hpp"
7#include "calendar.hpp"
[278]8#include "context.hpp"
[300]9#include "context_server.hpp"
[498]10#include "netCdfException.hpp"
11#include "exception.hpp"
[1158]12#include "timer.hpp"
13#include "uuid.hpp"
[335]14namespace xios
[219]15{
[878]16      /// ////////////////////// Dfinitions ////////////////////// ///
[219]17      CNc4DataOutput::CNc4DataOutput
[1158]18         (CFile* file, const StdString & filename, bool exist)
[219]19            : SuperClass()
20            , SuperClassWriter(filename, exist)
21            , filename(filename)
[2381]22            , file(file),hasTimeInstant(false),hasTimeCentered(false), timeCounterType(none), relElements_()
[219]23      {
[1158]24        SuperClass::type = MULTI_FILE;
[1456]25        compressionLevel= file->compression_level.isEmpty() ? 0 :file->compression_level ;
[219]26      }
27
28      CNc4DataOutput::CNc4DataOutput
[1158]29         (CFile* file, const StdString & filename, bool exist, bool useClassicFormat, bool useCFConvention,
[1639]30          MPI_Comm comm_file, bool multifile, bool isCollective, const StdString& timeCounterName)
[219]31            : SuperClass()
[878]32            , SuperClassWriter(filename, exist, useClassicFormat, useCFConvention, &comm_file, multifile, timeCounterName)
[379]33            , comm_file(comm_file)
[219]34            , filename(filename)
[335]35            , isCollective(isCollective)
[2381]36            , file(file),hasTimeInstant(false),hasTimeCentered(false), timeCounterType(none), relElements_()
[219]37      {
[1158]38        SuperClass::type = (multifile) ? MULTI_FILE : ONE_FILE;
[1459]39        if (file==NULL) compressionLevel = 0 ;
40        else compressionLevel= file->compression_level.isEmpty() ? 0 :file->compression_level ;
[219]41      }
42
43      CNc4DataOutput::~CNc4DataOutput(void)
[879]44    { /* Ne rien faire de plus */ }
[219]45
46      ///--------------------------------------------------------------
47
48      const StdString & CNc4DataOutput::getFileName(void) const
49      {
50         return (this->filename);
51      }
52
53      //---------------------------------------------------------------
54
[347]55      void CNc4DataOutput::writeDomain_(CDomain* domain)
[1622]56      TRY
[219]57      {
[1391]58        StdString lonName,latName ;
[1249]59
[881]60        if (domain->type == CDomain::type_attr::unstructured)
61        {
62          if (SuperClassWriter::useCFConvention)
63            writeUnstructuredDomain(domain) ;
64          else
65            writeUnstructuredDomainUgrid(domain) ;
66          return ;
67        }
[488]68
[347]69         CContext* context = CContext::getCurrent() ;
[219]70         if (domain->IsWritten(this->filename)) return;
71         domain->checkAttributes();
[488]72
73         if (domain->isEmpty())
[881]74           if (SuperClass::type==MULTI_FILE) return;
[219]75
[2381]76
77         // Check that the name associated to the current element is not in conflict with an existing element (due to CGrid::duplicateSentGrid)
78         if (!domain->lonvalue.isEmpty() )
79         {
80           int comm_file_rank(0);
81           MPI_Comm_rank( comm_file, &comm_file_rank );
82           int comm_file_size(1);
83           MPI_Comm_size( comm_file, &comm_file_size );
84         
85           // Get element current FULL view
86           shared_ptr<CLocalView> srcView = domain->getLocalView(CElementView::FULL) ;
87           vector<shared_ptr<CLocalView>> srcViews;
88           srcViews.push_back( srcView );
89         
90           // Compute a without redundancy element FULL view to enable a consistent hash computation
91           vector<shared_ptr<CLocalView>> remoteViews;
92           shared_ptr<CLocalView> remoteView;
93           // define the remote view without redundancy (naive distribution of the remote view)
94           int globalSize = domain->ni_glo.getValue()*domain->nj_glo.getValue();
95           int localSize = globalSize/comm_file_size;
96           if ( (comm_file_rank==comm_file_size-1) && (localSize*comm_file_size != globalSize ) )
97             localSize += globalSize-localSize*comm_file_size;
98           CArray<size_t,1> globalIndex( localSize );
99           CArray<int,1> index( localSize );
100           for (int iloc=0; iloc<localSize ; iloc++ )
101           {
102             globalIndex(iloc) = comm_file_rank*(globalSize/comm_file_size) + iloc;
103             index(iloc) = iloc;
104           }
105           shared_ptr<CLocalElement> localElement = make_shared<CLocalElement>(comm_file_rank, globalSize, globalIndex) ;
106           localElement->addView(CElementView::FULL, index) ;
107           remoteView = localElement->getView(CElementView::FULL) ;
108           remoteViews.push_back( remoteView );
109         
110           // Compute the connector between current and without redundancy FULL views
111           shared_ptr<CGridTransformConnector> gridTransformConnector = make_shared<CGridTransformConnector>(srcViews, remoteViews, comm_file ) ;
112           gridTransformConnector->computeConnector(true) ; // eliminateRedondant = true
113           CArray<double,1> lon_distributedValue, lat_distributedValue ;
114           gridTransformConnector->transfer(domain->lonvalue, lon_distributedValue );
115           gridTransformConnector->transfer(domain->latvalue, lat_distributedValue );
116
117           // Compute the distributed hash (v0) of the element
118           // it will be associated to the default element name (= map key), and to the name really written
119           int localHash = 0;
120           for (int iloc=0; iloc<localSize ; iloc++ ) localHash+=globalIndex(iloc)*lon_distributedValue(iloc)*lat_distributedValue(iloc);
121           int globalHash(0);
122           globalHash = localHash;
123           MPI_Allreduce( &localHash, &globalHash, 1, MPI_INT, MPI_SUM, comm_file  );
124         
125           StdString defaultNameKey = domain->getDomainOutputName();
[2384]126           if ( !relDomains_.count ( defaultNameKey ) )
[2381]127           {
128             // if defaultNameKey not in the map, write the element such as it is defined
[2384]129             relDomains_.insert( make_pair( defaultNameKey, make_pair(globalHash, domain) ) );
[2381]130           }
131           else // look if a hash associated this key is equal
132           {
133             bool elementIsInMap(false);
[2384]134             auto defaultNameKeyElements = relDomains_.equal_range( defaultNameKey );
[2381]135             for (auto it = defaultNameKeyElements.first; it != defaultNameKeyElements.second; it++)
136             {
137               if ( it->second.first == globalHash )
138               {
139                 // if yes, associate the same ids to current element
[2384]140                 domain->name = it->second.second->getDomainOutputName();
141                 domain->lon_name = it->second.second->lon_name;
142                 domain->lat_name = it->second.second->lat_name;
[2381]143                 elementIsInMap = true;
144               }
145             }
146             // if no : inheritance has been excessive, define new names and store it (could be used by another grid)
147             if (!elementIsInMap)  // ! in MAP
148             {
149               domain->name =  domain->getId();
150               domain->lon_name = "lon_"+domain->getId();
151               domain->lat_name = "lat_"+domain->getId();
[2384]152               relDomains_.insert( make_pair( defaultNameKey, make_pair(globalHash, domain) ) ) ;         
[2381]153             }
154           }
155         }
156
157         
[219]158         std::vector<StdString> dim0, dim1;
[772]159         StdString domid = domain->getDomainOutputName();
[318]160         StdString appendDomid  = (singleDomain) ? "" : "_"+domid ;
[706]161         if (isWrittenDomain(domid)) return ;
[774]162         else setWrittenDomain(domid);
[1132]163       
164         int nvertex = (domain->nvertex.isEmpty()) ? 0 : domain->nvertex;
[449]165
[488]166
[1158]167        StdString dimXid, dimYid ;
168
169        nc_type typePrec ;
170        if (domain->prec.isEmpty()) typePrec =  NC_FLOAT ;
171        else if (domain->prec==4)  typePrec =  NC_FLOAT ;
172        else if (domain->prec==8)   typePrec =  NC_DOUBLE ;
173         
[664]174         bool isRegularDomain = (domain->type == CDomain::type_attr::rectilinear);
[449]175         switch (domain->type)
[433]176         {
[449]177           case CDomain::type_attr::curvilinear :
[1391]178
179             if (domain->lon_name.isEmpty()) lonName = "nav_lon";
180             else lonName = domain->lon_name;
181
182             if (domain->lat_name.isEmpty()) latName = "nav_lat";
183             else latName = domain->lat_name;
184
[1430]185             if (domain->dim_i_name.isEmpty()) dimXid=StdString("x").append(appendDomid);
186             else dimXid=domain->dim_i_name.getValue() + appendDomid;
187
188             if (domain->dim_j_name.isEmpty()) dimYid=StdString("y").append(appendDomid);
189             else dimYid=domain->dim_j_name.getValue() + appendDomid;
190
[449]191             break ;
[1391]192
[664]193           case CDomain::type_attr::rectilinear :
[1391]194
[1430]195             if (domain->lon_name.isEmpty())
196             {
197               if (domain->dim_i_name.isEmpty())
198                   lonName = "lon";
199               else
200                 lonName = domain->dim_i_name.getValue();
201             }
[1391]202             else lonName = domain->lon_name;
203
[1430]204             if (domain->lat_name.isEmpty())
205             {
206               if (domain->dim_j_name.isEmpty())
207                 latName = "lat";
208               else
209                 latName = domain->dim_j_name.getValue();
210             }
[1391]211             else latName = domain->lat_name;
212             
[1430]213             if (domain->dim_i_name.isEmpty()) dimXid = lonName+appendDomid;
214             else dimXid = domain->dim_i_name.getValue()+appendDomid;
215
216             if (domain->dim_j_name.isEmpty()) dimYid = latName+appendDomid;
217             else dimYid = domain->dim_j_name.getValue()+appendDomid;
[449]218             break;
[488]219         }
220
[617]221         StdString dimVertId = StdString("nvertex").append(appendDomid);
222
[449]223         string lonid,latid,bounds_lonid,bounds_latid ;
[611]224         string areaId = "area" + appendDomid;
[219]225
[498]226         try
[219]227         {
[498]228           switch (SuperClass::type)
229           {
230              case (MULTI_FILE) :
231              {
232                 switch (domain->type)
233                 {
234                   case CDomain::type_attr::curvilinear :
235                     dim0.push_back(dimYid); dim0.push_back(dimXid);
[1415]236                     lonid = lonName+appendDomid;
237                     latid = latName+appendDomid;
[498]238                     break ;
[664]239                   case CDomain::type_attr::rectilinear :
[1415]240                     lonid = lonName+appendDomid;
241                     latid = latName+appendDomid;
[498]242                     dim0.push_back(dimYid);
243                     dim1.push_back(dimXid);
244                     break;
245                 }
[1430]246                 if (!domain->bounds_lon_name.isEmpty()) bounds_lonid = domain->bounds_lon_name;
247                 else bounds_lonid = "bounds_"+lonName+appendDomid;
248                 if (!domain->bounds_lat_name.isEmpty()) bounds_latid = domain->bounds_lat_name;
249                 else bounds_latid = "bounds_"+latName+appendDomid;
[488]250
[1553]251                 SuperClassWriter::addDimension(dimXid, domain->ni);
252                 SuperClassWriter::addDimension(dimYid, domain->nj);
[488]253
[617]254                 if (domain->hasBounds)
255                   SuperClassWriter::addDimension(dimVertId, domain->nvertex);
256
[1853]257                 if (context->intraCommSize_ > 1)
[498]258                 {
[1553]259                   this->writeLocalAttributes(domain->ibegin,
260                                              domain->ni,
261                                              domain->jbegin,
262                                              domain->nj,
[616]263                                              appendDomid);
[488]264
[628]265                   if (singleDomain)
266                    this->writeLocalAttributes_IOIPSL(dimXid, dimYid,
[1553]267                                                      domain->ibegin,
268                                                      domain->ni,
269                                                      domain->jbegin,
270                                                      domain->nj,
[628]271                                                      domain->ni_glo,domain->nj_glo,
[1853]272                                                      context->intraCommRank_,context->intraCommSize_);
[449]273                 }
[488]274
[665]275                 if (domain->hasLonLat)
[498]276                 {
[665]277                   switch (domain->type)
278                   {
279                     case CDomain::type_attr::curvilinear :
[1456]280                       SuperClassWriter::addVariable(latid, typePrec, dim0, compressionLevel);
281                       SuperClassWriter::addVariable(lonid, typePrec, dim0, compressionLevel);
[665]282                       break ;
283                      case CDomain::type_attr::rectilinear :
[1456]284                        SuperClassWriter::addVariable(latid, typePrec, dim0, compressionLevel);
285                        SuperClassWriter::addVariable(lonid, typePrec, dim1, compressionLevel);
[665]286                        break ;
287                   }
[488]288
[665]289                   this->writeAxisAttributes(lonid, isRegularDomain ? "X" : "", "longitude", "Longitude", "degrees_east", domid);
290                   this->writeAxisAttributes(latid, isRegularDomain ? "Y" : "", "latitude", "Latitude", "degrees_north", domid);
[219]291
[665]292                   if (domain->hasBounds)
293                   {
294                     SuperClassWriter::addAttribute("bounds", bounds_lonid, &lonid);
295                     SuperClassWriter::addAttribute("bounds", bounds_latid, &latid);
[617]296
[665]297                     dim0.clear();
298                     dim0.push_back(dimYid);
299                     dim0.push_back(dimXid);
300                     dim0.push_back(dimVertId);
[1456]301                     SuperClassWriter::addVariable(bounds_lonid, typePrec, dim0, compressionLevel);
302                     SuperClassWriter::addVariable(bounds_latid, typePrec, dim0, compressionLevel);
[665]303                   }
[617]304                 }
305
[498]306                 dim0.clear();
[616]307                 dim0.push_back(dimYid);
[498]308                 dim0.push_back(dimXid);
[219]309
[611]310                 if (domain->hasArea)
311                 {
[1456]312                   SuperClassWriter::addVariable(areaId, typePrec, dim0, compressionLevel);
[614]313                   SuperClassWriter::addAttribute("standard_name", StdString("cell_area"), &areaId);
[611]314                   SuperClassWriter::addAttribute("units", StdString("m2"), &areaId);
315                 }
316
[498]317                 SuperClassWriter::definition_end();
[449]318
[665]319                 if (domain->hasLonLat)
[498]320                 {
[665]321                   switch (domain->type)
322                   {
[1129]323                     case CDomain::type_attr::curvilinear :                       
[1930]324                       SuperClassWriter::writeData(domain->latvalue, latid, isCollective, 0);
325                       SuperClassWriter::writeData(domain->lonvalue, lonid, isCollective, 0);
[665]326                       break;
327                     case CDomain::type_attr::rectilinear :
[1930]328                       CArray<double,1> lat = domain->latvalue(Range(fromStart,toEnd,domain->ni)) ;
[665]329                       SuperClassWriter::writeData(CArray<double,1>(lat.copy()), latid, isCollective, 0);
[1930]330                       CArray<double,1> lon = domain->lonvalue(Range(0,domain->ni-1)) ;
[665]331                       SuperClassWriter::writeData(CArray<double,1>(lon.copy()), lonid, isCollective, 0);
332                       break;
333                   }
[611]334
[665]335                   if (domain->hasBounds)
336                   {
[1930]337                     SuperClassWriter::writeData(domain->bounds_lonvalue, bounds_lonid, isCollective, 0);
338                     SuperClassWriter::writeData(domain->bounds_latvalue, bounds_latid, isCollective, 0);
[665]339                   }
[617]340                 }
341
[611]342                 if (domain->hasArea)
[1129]343                 {
[1930]344                   SuperClassWriter::writeData(domain->areavalue, areaId, isCollective, 0);                   
[1129]345                 }
346
[498]347                 SuperClassWriter::definition_start();
[219]348
[498]349                 break;
350              }
351              case (ONE_FILE) :
352              {
[1553]353                SuperClassWriter::addDimension(dimXid, domain->ni_glo);
354                SuperClassWriter::addDimension(dimYid, domain->nj_glo);
[286]355
[617]356                 if (domain->hasBounds)
357                   SuperClassWriter::addDimension(dimVertId, domain->nvertex);
358
[665]359                 if (domain->hasLonLat)
[498]360                 {
[665]361                   switch (domain->type)
362                   {
363                     case CDomain::type_attr::curvilinear :
364                       dim0.push_back(dimYid); dim0.push_back(dimXid);
[1415]365                       lonid = lonName+appendDomid;
366                       latid = latName+appendDomid;
[1456]367                       SuperClassWriter::addVariable(latid, typePrec, dim0, compressionLevel);
368                       SuperClassWriter::addVariable(lonid, typePrec, dim0, compressionLevel);
[665]369                       break;
[449]370
[665]371                     case CDomain::type_attr::rectilinear :
372                       dim0.push_back(dimYid);
373                       dim1.push_back(dimXid);
[1415]374                       lonid = lonName+appendDomid;
375                       latid = latName+appendDomid;
[1456]376                       SuperClassWriter::addVariable(latid, typePrec, dim0, compressionLevel);
377                       SuperClassWriter::addVariable(lonid, typePrec, dim1, compressionLevel);
[665]378                       break;
379                   }
[1430]380                   if (!domain->bounds_lon_name.isEmpty()) bounds_lonid = domain->bounds_lon_name;
381                   else bounds_lonid = "bounds_"+lonName+appendDomid;
382                   if (!domain->bounds_lat_name.isEmpty()) bounds_latid = domain->bounds_lat_name;
383                   else bounds_latid = "bounds_"+latName+appendDomid;
[665]384
385                   this->writeAxisAttributes
386                      (lonid, isRegularDomain ? "X" : "", "longitude", "Longitude", "degrees_east", domid);
387                   this->writeAxisAttributes
388                      (latid, isRegularDomain ? "Y" : "", "latitude", "Latitude", "degrees_north", domid);
389
390                   if (domain->hasBounds)
391                   {
392                     SuperClassWriter::addAttribute("bounds", bounds_lonid, &lonid);
393                     SuperClassWriter::addAttribute("bounds", bounds_latid, &latid);
394
395                     dim0.clear();
[498]396                     dim0.push_back(dimYid);
[665]397                     dim0.push_back(dimXid);
398                     dim0.push_back(dimVertId);
[1456]399                     SuperClassWriter::addVariable(bounds_lonid, typePrec, dim0, compressionLevel);
400                     SuperClassWriter::addVariable(bounds_latid, typePrec, dim0, compressionLevel);
[665]401                   }
[498]402                 }
[611]403
404                 if (domain->hasArea)
405                 {
406                   dim0.clear();
407                   dim0.push_back(dimYid); dim0.push_back(dimXid);
[1456]408                   SuperClassWriter::addVariable(areaId, typePrec, dim0, compressionLevel);
[614]409                   SuperClassWriter::addAttribute("standard_name", StdString("cell_area"), &areaId);
[611]410                   SuperClassWriter::addAttribute("units", StdString("m2"), &areaId);
411                   dim0.clear();
412                 }
413
414                 SuperClassWriter::definition_end();
[286]415
[498]416                 switch (domain->type)
[384]417                 {
[498]418                   case CDomain::type_attr::curvilinear :
[449]419                   {
[498]420                     std::vector<StdSize> start(2) ;
421                     std::vector<StdSize> count(2) ;
[1959]422                     start[1]=domain->ibegin;
423                     start[0]=domain->jbegin;
424                     count[1]=domain->ni ; count[0]=domain->nj ;
[488]425
[665]426                     if (domain->hasLonLat)
427                     {
[1930]428                       SuperClassWriter::writeData(domain->latvalue, latid, isCollective, 0,&start,&count);
429                       SuperClassWriter::writeData(domain->lonvalue, lonid, isCollective, 0,&start,&count);
[665]430                     }
[498]431                     break;
432                   }
[664]433                   case CDomain::type_attr::rectilinear :
[449]434                   {
[665]435                     if (domain->hasLonLat)
[498]436                     {
[665]437                       std::vector<StdSize> start(1) ;
438                       std::vector<StdSize> count(1) ;
[1930]439                       
440                       start[0]=domain->jbegin;
441                       count[0]=domain->nj;
442                       CArray<double,1> lat = domain->latvalue(Range(fromStart,toEnd,domain->ni));
443                       SuperClassWriter::writeData(CArray<double,1>(lat.copy()), latid, isCollective, 0,&start,&count);
[665]444
[1930]445                       start[0]=domain->ibegin;
446                       count[0]=domain->ni;
447                       CArray<double,1> lon = domain->lonvalue(Range(0,domain->ni-1));
448                       SuperClassWriter::writeData(CArray<double,1>(lon.copy()), lonid, isCollective, 0,&start,&count);
[498]449                     }
450                     break;
[449]451                   }
[384]452                 }
[611]453
[617]454                 if (domain->hasBounds)
455                 {
456                   std::vector<StdSize> start(3);
457                   std::vector<StdSize> count(3);
458                   if (domain->isEmpty())
459                   {
460                     start[2] = start[1] = start[0] = 0;
461                     count[2] = count[1] = count[0] = 0;
462                   }
463                   else
464                   {
465                     start[2] = 0;
[1553]466                     start[1] = domain->ibegin;
467                     start[0] = domain->jbegin;
[1099]468                     count[2] = domain->nvertex;
[1553]469                     count[1] = domain->ni;
470                     count[0] = domain->nj;
[617]471                   }
[1143]472                 
[1930]473                   SuperClassWriter::writeData(domain->bounds_lonvalue, bounds_lonid, isCollective, 0, &start, &count); // will probably not working for rectilinear
474                   SuperClassWriter::writeData(domain->bounds_latvalue, bounds_latid, isCollective, 0, &start, &count);
[617]475                 }
476
[611]477                 if (domain->hasArea)
478                 {
479                   std::vector<StdSize> start(2);
480                   std::vector<StdSize> count(2);
481
[1930]482                   start[1] = domain->ibegin;
483                   start[0] = domain->jbegin;
484                   count[1] = domain->ni;
485                   count[0] = domain->nj;
[1143]486                   
[1930]487                   SuperClassWriter::writeData(domain->areavalue, areaId, isCollective, 0, &start, &count);
[611]488                 }
489
[498]490                 SuperClassWriter::definition_start();
491                 break;
492              }
493              default :
494                 ERROR("CNc4DataOutput::writeDomain(domain)",
495                       << "[ type = " << SuperClass::type << "]"
496                       << " not implemented yet !");
497           }
[449]498         }
[498]499         catch (CNetCdfException& e)
500         {
501           StdString msg("On writing the domain : ");
502           msg.append(domid); msg.append("\n");
503           msg.append("In the context : ");
504           msg.append(context->getId()); msg.append("\n");
505           msg.append(e.what());
506           ERROR("CNc4DataOutput::writeDomain_(CDomain* domain)", << msg);
507         }
508
[449]509         domain->addRelFile(this->filename);
510      }
[1622]511      CATCH
[449]512
[879]513    //--------------------------------------------------------------
[878]514
[879]515    void CNc4DataOutput::writeUnstructuredDomainUgrid(CDomain* domain)
516    {
517      CContext* context = CContext::getCurrent() ;
[878]518
[879]519      if (domain->IsWritten(this->filename)) return;
[1494]520
521      StdString domid = domain->getDomainOutputName();
522
523      // The first domain for the same mesh that will be written is that with the highest value of nvertex.
524      // Thus the entire mesh connectivity will be generated at once.
525      if (isWrittenDomain(domid)) return ;
526      else setWrittenDomain(domid);
527
[879]528      domain->checkAttributes();
529      if (domain->isEmpty())
530        if (SuperClass::type==MULTI_FILE) return ;
[878]531
[1158]532     nc_type typePrec ;
533     if (domain->prec.isEmpty()) typePrec =  NC_FLOAT ;
534     else if (domain->prec==4)  typePrec =  NC_FLOAT ;
535     else if (domain->prec==8)   typePrec =  NC_DOUBLE ;
536
[879]537      std::vector<StdString> dim0;
538      StdString domainName = domain->name;
[924]539      domain->assignMesh(domainName, domain->nvertex);
[1853]540      domain->mesh->createMeshEpsilon(context->intraComm_, domain->lonvalue, domain->latvalue, domain->bounds_lonvalue, domain->bounds_latvalue);
[878]541
[879]542      StdString node_x = domainName + "_node_x";
543      StdString node_y = domainName + "_node_y";
[878]544
[879]545      StdString edge_x = domainName + "_edge_x";
546      StdString edge_y = domainName + "_edge_y";
547      StdString edge_nodes = domainName + "_edge_nodes";
[878]548
[879]549      StdString face_x = domainName + "_face_x";
550      StdString face_y = domainName + "_face_y";
551      StdString face_nodes = domainName + "_face_nodes";
[902]552      StdString face_edges = domainName + "_face_edges";
553      StdString edge_faces = domainName + "_edge_face_links";
554      StdString face_faces = domainName + "_face_links";
[878]555
[879]556      StdString dimNode = "n" + domainName + "_node";
557      StdString dimEdge = "n" + domainName + "_edge";
558      StdString dimFace = "n" + domainName + "_face";
559      StdString dimVertex = "n" + domainName + "_vertex";
560      StdString dimTwo = "Two";
[878]561
[879]562      if (!SuperClassWriter::dimExist(dimTwo)) SuperClassWriter::addDimension(dimTwo, 2);
[1494]563      dim0.clear();
564      SuperClassWriter::addVariable(domainName, NC_INT, dim0, compressionLevel);
565      SuperClassWriter::addAttribute("cf_role", StdString("mesh_topology"), &domainName);
566      SuperClassWriter::addAttribute("long_name", StdString("Topology data of 2D unstructured mesh"), &domainName);
567      SuperClassWriter::addAttribute("topology_dimension", 2, &domainName);
568      SuperClassWriter::addAttribute("node_coordinates", node_x + " " + node_y, &domainName);
[878]569
[879]570      try
571      {
572        switch (SuperClass::type)
573        {
574          case (ONE_FILE) :
575          {
576            // Adding nodes
577            if (domain->nvertex == 1)
578            {
579              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y))
580              {
[902]581                SuperClassWriter::addDimension(dimNode, domain->ni_glo);
[879]582                dim0.clear();
583                dim0.push_back(dimNode);
[1456]584                SuperClassWriter::addVariable(node_x, typePrec, dim0, compressionLevel);
[879]585                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
[924]586                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x);
[879]587                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
[1456]588                SuperClassWriter::addVariable(node_y, typePrec, dim0, compressionLevel);
[879]589                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
[924]590                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y);
[879]591                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
592              }
593            } // domain->nvertex == 1
[878]594
[879]595            // Adding edges and nodes, if nodes have not been defined previously
596            if (domain->nvertex == 2)
597            {
598              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y))
599              {
[924]600                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodesGlo);
[879]601                dim0.clear();
602                dim0.push_back(dimNode);
[1456]603                SuperClassWriter::addVariable(node_x, typePrec, dim0, compressionLevel);
[879]604                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
[924]605                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x);
[879]606                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
[1456]607                SuperClassWriter::addVariable(node_y, typePrec, dim0, compressionLevel);
[879]608                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
[924]609                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y);
[879]610                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
611              }
[924]612              SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName);
613              SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName);
614              SuperClassWriter::addDimension(dimEdge, domain->ni_glo);
615              dim0.clear();
616              dim0.push_back(dimEdge);
[1456]617              SuperClassWriter::addVariable(edge_x, typePrec, dim0, compressionLevel);
[924]618              SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x);
619              SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh edges."), &edge_x);
620              SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x);
[1456]621              SuperClassWriter::addVariable(edge_y, typePrec, dim0, compressionLevel);
[924]622              SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y);
623              SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh edges."), &edge_y);
624              SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y);
625              dim0.clear();
626              dim0.push_back(dimEdge);
627              dim0.push_back(dimTwo);
[1456]628              SuperClassWriter::addVariable(edge_nodes, NC_INT, dim0, compressionLevel);
[924]629              SuperClassWriter::addAttribute("cf_role", StdString("edge_node_connectivity"), &edge_nodes);
630              SuperClassWriter::addAttribute("long_name", StdString("Maps every edge/link to two nodes that it connects."), &edge_nodes);
631              SuperClassWriter::addAttribute("start_index", 0, &edge_nodes);
[879]632            } // domain->nvertex == 2
[878]633
[879]634            // Adding faces, edges, and nodes, if edges and nodes have not been defined previously
635            if (domain->nvertex > 2)
636            {
637              // Nodes
638              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y))
639              {
[924]640                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodesGlo);
[879]641                dim0.clear();
642                dim0.push_back(dimNode);
[1456]643                SuperClassWriter::addVariable(node_x, typePrec, dim0, compressionLevel);
[879]644                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
[924]645                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x);
[879]646                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
[1456]647                SuperClassWriter::addVariable(node_y, typePrec, dim0, compressionLevel);
[879]648                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
[924]649                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y);
[879]650                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
651              }
652              if (!SuperClassWriter::varExist(edge_x) || !SuperClassWriter::varExist(edge_y))
653              {
654                SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName);
655                SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName);
[924]656                SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdgesGlo);
[879]657                dim0.clear();
658                dim0.push_back(dimEdge);
[1456]659                SuperClassWriter::addVariable(edge_x, typePrec, dim0, compressionLevel);
[879]660                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x);
[924]661                SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh edges."), &edge_x);
[879]662                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x);
[1456]663                SuperClassWriter::addVariable(edge_y, typePrec, dim0, compressionLevel);
[879]664                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y);
[924]665                SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh edges."), &edge_y);
[879]666                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y);
667                dim0.clear();
668                dim0.push_back(dimEdge);
669                dim0.push_back(dimTwo);
[1456]670                SuperClassWriter::addVariable(edge_nodes, NC_INT, dim0, compressionLevel);
[879]671                SuperClassWriter::addAttribute("cf_role", StdString("edge_node_connectivity"), &edge_nodes);
672                SuperClassWriter::addAttribute("long_name", StdString("Maps every edge/link to two nodes that it connects."), &edge_nodes);
673                SuperClassWriter::addAttribute("start_index", 0, &edge_nodes);
674              }
[924]675              SuperClassWriter::addAttribute("face_coordinates", face_x + " " + face_y, &domainName);
676              SuperClassWriter::addAttribute("face_node_connectivity", face_nodes, &domainName);
677              SuperClassWriter::addDimension(dimFace, domain->ni_glo);
678              SuperClassWriter::addDimension(dimVertex, domain->nvertex);
679              dim0.clear();
680              dim0.push_back(dimFace);
[1456]681              SuperClassWriter::addVariable(face_x, typePrec, dim0, compressionLevel);
[924]682              SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &face_x);
683              SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh faces."), &face_x);
684              SuperClassWriter::addAttribute("units", StdString("degrees_east"), &face_x);
[1456]685              SuperClassWriter::addVariable(face_y, typePrec, dim0, compressionLevel);
[924]686              SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &face_y);
687              SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh faces."), &face_y);
688              SuperClassWriter::addAttribute("units", StdString("degrees_north"), &face_y);
689              dim0.clear();
690              dim0.push_back(dimFace);
691              dim0.push_back(dimVertex);
[1456]692              SuperClassWriter::addVariable(face_nodes, NC_INT, dim0, compressionLevel);
[924]693              SuperClassWriter::addAttribute("cf_role", StdString("face_node_connectivity"), &face_nodes);
694              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its corner nodes."), &face_nodes);
695              SuperClassWriter::addAttribute("start_index", 0, &face_nodes);
696              dim0.clear();
697              dim0.push_back(dimFace);
698              dim0.push_back(dimVertex);
[1456]699              SuperClassWriter::addVariable(face_edges, NC_INT, dim0, compressionLevel);
[924]700              SuperClassWriter::addAttribute("cf_role", StdString("face_edge_connectivity"), &face_edges);
701              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges);
702              SuperClassWriter::addAttribute("start_index", 0, &face_edges);
[929]703              SuperClassWriter::addAttribute("_FillValue", 999999, &face_edges);
[924]704              dim0.clear();
705              dim0.push_back(dimEdge);
706              dim0.push_back(dimTwo);
[1456]707              SuperClassWriter::addVariable(edge_faces, NC_INT, dim0, compressionLevel);
[924]708              SuperClassWriter::addAttribute("cf_role", StdString("edge_face connectivity"), &edge_faces);
709              SuperClassWriter::addAttribute("long_name", StdString("neighbor faces for edges"), &edge_faces);
710              SuperClassWriter::addAttribute("start_index", 0, &edge_faces);
711              SuperClassWriter::addAttribute("_FillValue", -999, &edge_faces);
712              SuperClassWriter::addAttribute("comment", StdString("missing neighbor faces are indicated using _FillValue"), &edge_faces);
713              dim0.clear();
714              dim0.push_back(dimFace);
715              dim0.push_back(dimVertex);
[1456]716              SuperClassWriter::addVariable(face_faces, NC_INT, dim0, compressionLevel);
[924]717              SuperClassWriter::addAttribute("cf_role", StdString("face_face connectivity"), &face_faces);
718              SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces);
719              SuperClassWriter::addAttribute("start_index", 0, &face_faces);
[929]720              SuperClassWriter::addAttribute("_FillValue", 999999, &face_faces);
[924]721              SuperClassWriter::addAttribute("flag_values", -1, &face_faces);
722              SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces);
[879]723            } // domain->nvertex > 2
[878]724
[879]725            SuperClassWriter::definition_end();
[878]726
[924]727            std::vector<StdSize> startEdges(1) ;
728            std::vector<StdSize> countEdges(1) ;
729            std::vector<StdSize> startNodes(1) ;
730            std::vector<StdSize> countNodes(1) ;
731            std::vector<StdSize> startFaces(1) ;
732            std::vector<StdSize> countFaces(1) ;
733            std::vector<StdSize> startEdgeNodes(2) ;
734            std::vector<StdSize> countEdgeNodes(2) ;
735            std::vector<StdSize> startEdgeFaces(2) ;
736            std::vector<StdSize> countEdgeFaces(2) ;
737            std::vector<StdSize> startFaceConctv(2) ;
738            std::vector<StdSize> countFaceConctv(2) ;
[902]739
[1494]740            if (domain->nvertex == 1)
[879]741            {
[1494]742              if (domain->isEmpty())
743               {
744                 startNodes[0]=0 ;
745                 countNodes[0]=0 ;
746               }
747               else
748               {
[1553]749                 startNodes[0] = domain->ibegin;
750                 countNodes[0] = domain->ni ;
[1494]751               }
[924]752
[1494]753              SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes);
754              SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes);
755            }
756            else if (domain->nvertex == 2)
757            {
758              if (domain->isEmpty())
759               {
760                startEdges[0]=0 ;
761                countEdges[0]=0 ;
762                startNodes[0]=0 ;
763                countNodes[0]=0 ;
764                startEdgeNodes[0]=0;
765                startEdgeNodes[1]=0;
766                countEdgeNodes[0]=0;
767                countEdgeNodes[1]=0;
[924]768
[1494]769               }
770               else
771               {
[1553]772                 startEdges[0] = domain->ibegin;
773                 countEdges[0] = domain->ni;
[1494]774                 startNodes[0] = domain->mesh->node_start;
775                 countNodes[0] = domain->mesh->node_count;
[1553]776                 startEdgeNodes[0] = domain->ibegin;
[1494]777                 startEdgeNodes[1] = 0;
[1553]778                 countEdgeNodes[0] = domain->ni;
[1494]779                 countEdgeNodes[1] = 2;
780               }
781              SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes);
782              SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes);
783              SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges);
784              SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges);
785              SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes);
786            }
[879]787            else
788            {
[1494]789              if (domain->isEmpty())
790               {
791                 startFaces[0] = 0 ;
792                 countFaces[0] = 0 ;
793                 startNodes[0] = 0;
794                 countNodes[0] = 0;
795                 startEdges[0] = 0;
796                 countEdges[0] = 0;
797                 startEdgeFaces[0] = 0;
798                 startEdgeFaces[1] = 0;
799                 countEdgeFaces[0] = 0;
800                 countEdgeFaces[1] = 0;
801                 startFaceConctv[0] = 0;
802                 startFaceConctv[1] = 0;
803                 countFaceConctv[0] = 0;
804                 countFaceConctv[1] = 0;
805               }
806               else
807               {
[1553]808                 startFaces[0] = domain->ibegin;
809                 countFaces[0] = domain->ni ;
[1494]810                 startNodes[0] = domain->mesh->node_start;
811                 countNodes[0] = domain->mesh->node_count;
812                 startEdges[0] = domain->mesh->edge_start;
813                 countEdges[0] = domain->mesh->edge_count;
814                 startEdgeNodes[0] = domain->mesh->edge_start;
815                 startEdgeNodes[1] = 0;
816                 countEdgeNodes[0] = domain->mesh->edge_count;
817                 countEdgeNodes[1]= 2;
818                 startEdgeFaces[0] = domain->mesh->edge_start;
819                 startEdgeFaces[1]= 0;
820                 countEdgeFaces[0] = domain->mesh->edge_count;
821                 countEdgeFaces[1]= 2;
[1553]822                 startFaceConctv[0] = domain->ibegin;
[1494]823                 startFaceConctv[1] = 0;
[1553]824                 countFaceConctv[0] = domain->ni;
[1494]825                 countFaceConctv[1] = domain->nvertex;
826               }
827              SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes);
828              SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes);
829              SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges);
830              SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges);
831              SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes);
832              SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0, &startFaces, &countFaces);
833              SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0, &startFaces, &countFaces);
834              SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes, isCollective, 0, &startFaceConctv, &countFaceConctv);
835              SuperClassWriter::writeData(domain->mesh->face_edges, face_edges, isCollective, 0, &startFaceConctv, &countFaceConctv);
836              SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces, isCollective, 0, &startEdgeFaces, &countEdgeFaces);
837              SuperClassWriter::writeData(domain->mesh->face_faces, face_faces, isCollective, 0, &startFaceConctv, &countFaceConctv);
838            }
[879]839            SuperClassWriter::definition_start();
[878]840
[879]841            break;
842          } // ONE_FILE
[878]843
[879]844          case (MULTI_FILE) :
845          {
[1637]846            ERROR("CNc4DataOutput::writeDomain(domain)",
847            << "[ type = multiple_file ]"
848            << " is not yet implemented for UGRID files !");
[879]849            break;
850          }
[878]851
[879]852          default :
853          ERROR("CNc4DataOutput::writeDomain(domain)",
854          << "[ type = " << SuperClass::type << "]"
855          << " not implemented yet !");
856          } // switch
857        } // try
[878]858
[879]859        catch (CNetCdfException& e)
860        {
861          StdString msg("On writing the domain : ");
862          msg.append(domid); msg.append("\n");
863          msg.append("In the context : ");
864          msg.append(context->getId()); msg.append("\n");
865          msg.append(e.what());
[924]866          ERROR("CNc4DataOutput::writeUnstructuredDomainUgrid(CDomain* domain)", << msg);
[879]867        }
[878]868
[879]869  domain->addRelFile(this->filename);
870  }
[878]871
[879]872    //--------------------------------------------------------------
[878]873
[879]874    void CNc4DataOutput::writeUnstructuredDomain(CDomain* domain)
[449]875      {
876         CContext* context = CContext::getCurrent() ;
[488]877
[449]878         if (domain->IsWritten(this->filename)) return;
879         domain->checkAttributes();
[488]880
881         if (domain->isEmpty())
[449]882           if (SuperClass::type==MULTI_FILE) return ;
883
884         std::vector<StdString> dim0, dim1;
[772]885         StdString domid = domain->getDomainOutputName();
[705]886         if (isWrittenDomain(domid)) return ;
[774]887         else setWrittenDomain(domid);
[705]888
[449]889         StdString appendDomid  = (singleDomain) ? "" : "_"+domid ;
890
[1430]891         StdString lonName,latName, cellName ;
[1391]892         if (domain->lon_name.isEmpty()) lonName = "lon";
893         else lonName = domain->lon_name;
894
895         if (domain->lat_name.isEmpty()) latName = "lat";
896         else latName = domain->lat_name;
897
[1430]898         if (!domain->dim_i_name.isEmpty()) cellName=domain->dim_i_name;
899         else cellName="cell";
900         StdString dimXid = cellName+appendDomid;
[449]901         StdString dimVertId = StdString("nvertex").append(appendDomid);
[488]902
[449]903         string lonid,latid,bounds_lonid,bounds_latid ;
[611]904         string areaId = "area" + appendDomid;
[449]905
[1158]906         nc_type typePrec ;
907         if (domain->prec.isEmpty()) typePrec =  NC_FLOAT ;
908         else if (domain->prec==4)  typePrec =  NC_FLOAT ;
909         else if (domain->prec==8)   typePrec =  NC_DOUBLE ;
910
[1132]911         int nvertex = (domain->nvertex.isEmpty()) ? 0 : domain->nvertex;
[1025]912
[498]913         try
[449]914         {
[498]915           switch (SuperClass::type)
916           {
917              case (MULTI_FILE) :
918              {
919                 dim0.push_back(dimXid);
[1553]920                 SuperClassWriter::addDimension(dimXid, domain->ni);
[488]921
[1415]922                 lonid = lonName+appendDomid;
923                 latid = latName+appendDomid;
[1430]924                 if (!domain->bounds_lon_name.isEmpty()) bounds_lonid = domain->bounds_lon_name;
925                 else bounds_lonid = "bounds_"+lonName+appendDomid;
926                 if (!domain->bounds_lat_name.isEmpty()) bounds_latid = domain->bounds_lat_name;
927                 else bounds_latid = "bounds_"+latName+appendDomid;
928
[665]929                 if (domain->hasLonLat)
930                 {
[1456]931                   SuperClassWriter::addVariable(latid, typePrec, dim0, compressionLevel);
932                   SuperClassWriter::addVariable(lonid, typePrec, dim0, compressionLevel);
[665]933                   this->writeAxisAttributes(lonid, "", "longitude", "Longitude", "degrees_east", domid);
934                   if (domain->hasBounds) SuperClassWriter::addAttribute("bounds",bounds_lonid, &lonid);
935                   this->writeAxisAttributes(latid, "", "latitude", "Latitude", "degrees_north", domid);
936                   if (domain->hasBounds) SuperClassWriter::addAttribute("bounds",bounds_latid, &latid);
937                   if (domain->hasBounds) SuperClassWriter::addDimension(dimVertId, domain->nvertex);
938                 }
[498]939                 dim0.clear();
940                 if (domain->hasBounds)
941                 {
942                   dim0.push_back(dimXid);
943                   dim0.push_back(dimVertId);
[1456]944                   SuperClassWriter::addVariable(bounds_lonid, typePrec, dim0, compressionLevel);
945                   SuperClassWriter::addVariable(bounds_latid, typePrec, dim0, compressionLevel);
[498]946                 }
947
948                 dim0.clear();
[449]949                 dim0.push_back(dimXid);
[611]950                 if (domain->hasArea)
951                 {
[1456]952                   SuperClassWriter::addVariable(areaId, typePrec, dim0, compressionLevel);
[614]953                   SuperClassWriter::addAttribute("standard_name", StdString("cell_area"), &areaId);
[611]954                   SuperClassWriter::addAttribute("units", StdString("m2"), &areaId);
955                 }
956
[498]957                 SuperClassWriter::definition_end();
[449]958
[665]959                 if (domain->hasLonLat)
[498]960                 {
[1930]961                   SuperClassWriter::writeData(domain->latvalue, latid, isCollective, 0);
962                   SuperClassWriter::writeData(domain->lonvalue, lonid, isCollective, 0);
[665]963                   if (domain->hasBounds)
964                   {
[1930]965                     SuperClassWriter::writeData(domain->bounds_lonvalue, bounds_lonid, isCollective, 0);
966                     SuperClassWriter::writeData(domain->bounds_latvalue, bounds_latid, isCollective, 0);
[665]967                   }
[498]968                 }
[611]969
970                 if (domain->hasArea)
[1930]971                   SuperClassWriter::writeData(domain->areavalue, areaId, isCollective, 0);
[611]972
[498]973                 SuperClassWriter::definition_start();
974                 break ;
975              }
[488]976
[498]977              case (ONE_FILE) :
978              {
[1415]979                 lonid = lonName+appendDomid;
980                 latid = latName+appendDomid;
[1430]981                 if (!domain->bounds_lon_name.isEmpty()) bounds_lonid = domain->bounds_lon_name;
982                 else bounds_lonid = "bounds_"+lonName+appendDomid;
983                 if (!domain->bounds_lat_name.isEmpty()) bounds_latid = domain->bounds_lat_name;
984                 else bounds_latid = "bounds_"+latName+appendDomid;
[1391]985
[498]986                 dim0.push_back(dimXid);
[1132]987                 SuperClassWriter::addDimension(dimXid, domain->ni_glo);
[665]988                 if (domain->hasLonLat)
989                 {
[1456]990                   SuperClassWriter::addVariable(latid, typePrec, dim0, compressionLevel);
991                   SuperClassWriter::addVariable(lonid, typePrec, dim0, compressionLevel);
[665]992
993                   this->writeAxisAttributes(lonid, "", "longitude", "Longitude", "degrees_east", domid);
994                   if (domain->hasBounds) SuperClassWriter::addAttribute("bounds",bounds_lonid, &lonid);
995                   this->writeAxisAttributes(latid, "", "latitude", "Latitude", "degrees_north", domid);
996                   if (domain->hasBounds) SuperClassWriter::addAttribute("bounds",bounds_latid, &latid);
[1025]997                   if (domain->hasBounds) SuperClassWriter::addDimension(dimVertId, nvertex);
[665]998                 }
[498]999                 dim0.clear();
[449]1000
[498]1001                 if (domain->hasBounds)
1002                 {
1003                   dim0.push_back(dimXid);
1004                   dim0.push_back(dimVertId);
[1456]1005                   SuperClassWriter::addVariable(bounds_lonid, typePrec, dim0, compressionLevel);
1006                   SuperClassWriter::addVariable(bounds_latid, typePrec, dim0, compressionLevel);
[498]1007                 }
[488]1008
[611]1009                 if (domain->hasArea)
1010                 {
1011                   dim0.clear();
1012                   dim0.push_back(dimXid);
[1456]1013                   SuperClassWriter::addVariable(areaId, typePrec, dim0, compressionLevel);
[614]1014                   SuperClassWriter::addAttribute("standard_name", StdString("cell_area"), &areaId);
[611]1015                   SuperClassWriter::addAttribute("units", StdString("m2"), &areaId);
1016                 }
1017
[498]1018                 SuperClassWriter::definition_end();
[488]1019
[498]1020                 std::vector<StdSize> start(1), startBounds(2) ;
1021                 std::vector<StdSize> count(1), countBounds(2) ;
1022                 if (domain->isEmpty())
1023                 {
1024                   start[0]=0 ;
1025                   count[0]=0 ;
1026                   startBounds[1]=0 ;
[1025]1027                   countBounds[1]=nvertex ;
[498]1028                   startBounds[0]=0 ;
1029                   countBounds[0]=0 ;
1030                 }
1031                 else
1032                 {
[1553]1033                   start[0]=domain->ibegin;
1034                   count[0]=domain->ni;
1035                   startBounds[0]=domain->ibegin;
[498]1036                   startBounds[1]=0 ;
[1553]1037                   countBounds[0]=domain->ni;
[1025]1038                   countBounds[1]=nvertex ;
[498]1039                 }
[665]1040
1041                 if (domain->hasLonLat)
[498]1042                 {
[1930]1043                   SuperClassWriter::writeData(domain->latvalue, latid, isCollective, 0,&start,&count);
1044                   SuperClassWriter::writeData(domain->lonvalue, lonid, isCollective, 0,&start,&count);
[665]1045                   if (domain->hasBounds)
1046                   {
[1930]1047                     SuperClassWriter::writeData(domain->bounds_lonvalue, bounds_lonid, isCollective, 0,&startBounds,&countBounds);
1048                     SuperClassWriter::writeData(domain->bounds_latvalue, bounds_latid, isCollective, 0,&startBounds,&countBounds);
[665]1049                   }
[498]1050                 }
[488]1051
[611]1052                 if (domain->hasArea)
[1930]1053                   SuperClassWriter::writeData(domain->areavalue, areaId, isCollective, 0, &start, &count);
[488]1054
[498]1055                 SuperClassWriter::definition_start();
[488]1056
[498]1057                 break;
1058              }
1059              default :
1060                 ERROR("CNc4DataOutput::writeDomain(domain)",
1061                       << "[ type = " << SuperClass::type << "]"
1062                       << " not implemented yet !");
1063           }
[219]1064         }
[498]1065         catch (CNetCdfException& e)
1066         {
1067           StdString msg("On writing the domain : ");
1068           msg.append(domid); msg.append("\n");
1069           msg.append("In the context : ");
1070           msg.append(context->getId()); msg.append("\n");
1071           msg.append(e.what());
1072           ERROR("CNc4DataOutput::writeUnstructuredDomain(CDomain* domain)", << msg);
1073         }
[219]1074         domain->addRelFile(this->filename);
1075      }
1076      //--------------------------------------------------------------
1077
[347]1078      void CNc4DataOutput::writeAxis_(CAxis* axis)
[219]1079      {
[1030]1080        if (axis->IsWritten(this->filename)) return;
[2381]1081
1082        // Check that the name associated to the current element is not in conflict with an existing element (due to CGrid::duplicateSentGrid)
1083        if (!axis->value.isEmpty() )
1084        {
1085          int comm_file_rank(0);
1086          MPI_Comm_rank( comm_file, &comm_file_rank );
1087          int comm_file_size(1);
1088          MPI_Comm_size( comm_file, &comm_file_size );
1089
1090          // Get element current FULL view
1091          shared_ptr<CLocalView> srcView = axis->getLocalView(CElementView::FULL) ;
1092          vector<shared_ptr<CLocalView>> srcViews;
1093          srcViews.push_back( srcView );
1094         
1095          // Compute a without redundancy element FULL view to enable a consistent hash computation
1096          vector<shared_ptr<CLocalView>> remoteViews;
1097          shared_ptr<CLocalView> remoteView;
1098          // define the remote view without redundancy (naive distribution of the remote view)
1099          int globalSize = axis->n_glo.getValue();
1100          int localSize = globalSize/comm_file_size;
1101          if ( (comm_file_rank==comm_file_size-1) && (localSize*comm_file_size != globalSize ) )
1102            localSize += globalSize-localSize*comm_file_size;
1103          CArray<size_t,1> globalIndex( localSize );
1104          CArray<int,1> index( localSize );
1105          for (int iloc=0; iloc<localSize ; iloc++ )
1106          {
1107            globalIndex(iloc) = comm_file_rank*(globalSize/comm_file_size) + iloc;
1108            index(iloc) = iloc;
1109          }
1110          shared_ptr<CLocalElement> localElement = make_shared<CLocalElement>(comm_file_rank, globalSize, globalIndex) ;
1111          localElement->addView(CElementView::FULL, index) ;
1112          remoteView = localElement->getView(CElementView::FULL) ;
1113          remoteViews.push_back( remoteView );
1114       
1115          // Compute the connector between current and without redundancy FULL views
1116          shared_ptr<CGridTransformConnector> gridTransformConnector = make_shared<CGridTransformConnector>(srcViews, remoteViews, comm_file ) ;
1117          gridTransformConnector->computeConnector(true) ; // eliminateRedondant = true
1118          CArray<double,1> distributedValue ;
1119          gridTransformConnector->transfer(axis->value, distributedValue );
1120       
1121          // Compute the distributed hash (v0) of the element
1122          // it will be associated to the default element name (= map key), and to the name really written
1123          int localHash = 0;
1124          for (int iloc=0; iloc<localSize ; iloc++ ) localHash+=globalIndex(iloc)*distributedValue(iloc);
1125          int globalHash(0);
1126          globalHash = localHash;
1127          MPI_Allreduce( &localHash, &globalHash, 1, MPI_INT, MPI_SUM, comm_file  );
1128
1129          StdString defaultNameKey = axis->getAxisOutputName();
1130          if ( !relElements_.count ( defaultNameKey ) )
1131          {
1132            // if defaultNameKey not in the map, write the element such as it is defined
1133            relElements_.insert( make_pair( defaultNameKey, make_pair(globalHash, defaultNameKey) ) );
1134          }
1135          else // look if a hash associated this key is equal
1136          {
1137            bool elementIsInMap(false);
1138            auto defaultNameKeyElements = relElements_.equal_range( defaultNameKey );
1139            for (auto it = defaultNameKeyElements.first; it != defaultNameKeyElements.second; it++)
1140            {
1141              if ( it->second.first == globalHash )
1142              {
1143                // if yes, associate the same ids to current element
1144                axis->name = it->second.second;
1145                elementIsInMap = true;
1146              }
1147            }
1148             // if no : inheritance has been excessive, define new names and store it (could be used by another grid)
1149            if (!elementIsInMap)  // ! in MAP
1150            {
1151              axis->name =  axis->getId();
1152              relElements_.insert( make_pair( defaultNameKey, make_pair(globalHash, axis->getAxisOutputName()) ) ) ;// = axis->getId()         
1153            }
1154          }
1155        }
1156
[1030]1157        axis->checkAttributes();
[488]1158
[1559]1159        int size  = (MULTI_FILE == SuperClass::type) ? axis->n.getValue()
1160                                                          : axis->n_glo.getValue();
[1099]1161
[1559]1162        if ((0 == axis->n) && (MULTI_FILE == SuperClass::type)) return;
[1235]1163
[1030]1164        std::vector<StdString> dims;
1165        StdString axisid = axis->getAxisOutputName();
[1435]1166        StdString axisDim, axisBoundsId;
[1030]1167        if (isWrittenAxis(axisid)) return ;
1168        else setWrittenAxis(axisid);
[705]1169
[1158]1170        nc_type typePrec ;
1171        if (axis->prec.isEmpty()) typePrec =  NC_FLOAT ;
[1235]1172        else if (axis->prec==4)   typePrec =  NC_FLOAT ;
[1158]1173        else if (axis->prec==8)   typePrec =  NC_DOUBLE ;
1174         
1175        if (!axis->label.isEmpty()) typePrec = NC_CHAR ;
1176        string strId="str_len" ;
[1030]1177        try
1178        {
[1433]1179          if (axis->dim_name.isEmpty()) axisDim = axisid;
1180          else axisDim=axis->dim_name.getValue();
[1559]1181          SuperClassWriter::addDimension(axisDim, size);
[1433]1182          dims.push_back(axisDim);
1183
[1439]1184          if (!axis->label.isEmpty() && !SuperClassWriter::dimExist(strId)) SuperClassWriter::addDimension(strId, stringArrayLen);
[1433]1185
[1435]1186          if (axis->hasValue || !axis->label.isEmpty())
[816]1187          {
[1158]1188            if (!axis->label.isEmpty()) dims.push_back(strId);
[1430]1189
[1456]1190            SuperClassWriter::addVariable(axisid, typePrec, dims, compressionLevel);
[219]1191
[1030]1192            if (!axis->name.isEmpty())
1193              SuperClassWriter::addAttribute("name", axis->name.getValue(), &axisid);
[219]1194
[1030]1195            if (!axis->standard_name.isEmpty())
1196              SuperClassWriter::addAttribute("standard_name", axis->standard_name.getValue(), &axisid);
[540]1197
[1030]1198            if (!axis->long_name.isEmpty())
1199              SuperClassWriter::addAttribute("long_name", axis->long_name.getValue(), &axisid);
[219]1200
[1030]1201            if (!axis->unit.isEmpty())
1202              SuperClassWriter::addAttribute("units", axis->unit.getValue(), &axisid);
[219]1203
[1430]1204            if (!axis->axis_type.isEmpty())
1205            {
1206              switch(axis->axis_type)
1207              {
1208              case CAxis::axis_type_attr::X :
1209                SuperClassWriter::addAttribute("axis", string("X"), &axisid);
1210                break;
1211              case CAxis::axis_type_attr::Y :
1212                SuperClassWriter::addAttribute("axis", string("Y"), &axisid);
1213                break;
1214              case CAxis::axis_type_attr::Z :
1215                SuperClassWriter::addAttribute("axis", string("Z"), &axisid);
1216                break;
1217              case CAxis::axis_type_attr::T :
1218                SuperClassWriter::addAttribute("axis", string("T"), &axisid);
1219                break;
1220              }
1221            }
1222
[1030]1223            if (!axis->positive.isEmpty())
1224            {
1225              SuperClassWriter::addAttribute("positive",
1226                                             (axis->positive == CAxis::positive_attr::up) ? string("up") : string("down"),
1227                                             &axisid);
1228            }
[399]1229
[1266]1230            if (!axis->formula.isEmpty())
1231              SuperClassWriter::addAttribute("formula", axis->formula.getValue(), &axisid);
1232
1233            if (!axis->formula_term.isEmpty())
1234              SuperClassWriter::addAttribute("formula_term", axis->formula_term.getValue(), &axisid);
1235             
[1435]1236            axisBoundsId = (axis->bounds_name.isEmpty()) ? axisid + "_bounds" : axis->bounds_name;
[1249]1237            if (!axis->bounds.isEmpty() && axis->label.isEmpty())
[1030]1238            {
1239              dims.push_back("axis_nbounds");
[1456]1240              SuperClassWriter::addVariable(axisBoundsId, typePrec, dims, compressionLevel);
[1030]1241              SuperClassWriter::addAttribute("bounds", axisBoundsId, &axisid);
[1288]1242
1243              if (!axis->standard_name.isEmpty())
1244                SuperClassWriter::addAttribute("standard_name", axis->standard_name.getValue(), &axisBoundsId);
1245
1246              if (!axis->unit.isEmpty())
1247                SuperClassWriter::addAttribute("units", axis->unit.getValue(), &axisBoundsId);
1248
[1266]1249              if (!axis->formula_bounds.isEmpty())
[1288]1250                SuperClassWriter::addAttribute("formula", axis->formula_bounds.getValue(), &axisBoundsId);
[1266]1251
1252              if (!axis->formula_term_bounds.isEmpty())
[1288]1253                SuperClassWriter::addAttribute("formula_term", axis->formula_term_bounds.getValue(), &axisBoundsId);
[1030]1254            }
[1435]1255          }
[816]1256
[1435]1257          SuperClassWriter::definition_end();
[1950]1258         
[1435]1259          switch (SuperClass::type)
1260          {
1261            case MULTI_FILE:
[1030]1262            {
[1435]1263              if (axis->label.isEmpty())
[1437]1264              {
1265                if (!axis->value.isEmpty())
[1930]1266                  SuperClassWriter::writeData(axis->value, axisid, isCollective, 0);
[633]1267
[1437]1268                if (!axis->bounds.isEmpty())
[1930]1269                  SuperClassWriter::writeData(axis->bounds, axisBoundsId, isCollective, 0);
[1435]1270              }
[1437]1271              else
[1950]1272                SuperClassWriter::writeData(axis->label, axisid, isCollective, 0);
[1435]1273
1274              SuperClassWriter::definition_start();
1275              break;
1276            }
1277            case ONE_FILE:
1278            {
1279              std::vector<StdSize> start(1), startBounds(2) ;
1280              std::vector<StdSize> count(1), countBounds(2) ;
[1559]1281              start[0] = startBounds[0] = axis->begin;
1282              count[0] = countBounds[0] = axis->n;
[1435]1283              startBounds[1] = 0;
1284              countBounds[1] = 2;
1285
1286              if (axis->label.isEmpty())
[1437]1287              {
1288                if (!axis->value.isEmpty())
[1930]1289                  SuperClassWriter::writeData(axis->value, axisid, isCollective, 0, &start, &count);
[1435]1290
[1437]1291                if (!axis->bounds.isEmpty())
[1930]1292                  SuperClassWriter::writeData(axis->bounds, axisBoundsId, isCollective, 0, &startBounds, &countBounds);
[1435]1293              }
[1437]1294              else
[1435]1295              {
1296                std::vector<StdSize> startLabel(2), countLabel(2);
1297                startLabel[0] = start[0]; startLabel[1] = 0;
1298                countLabel[0] = count[0]; countLabel[1] = stringArrayLen;
[1950]1299                SuperClassWriter::writeData(axis->label, axisid, isCollective, 0, &startLabel, &countLabel);
[1435]1300              }
[609]1301
[1435]1302              SuperClassWriter::definition_start();
1303
1304              break;
[609]1305            }
[1435]1306            default :
1307              ERROR("CNc4DataOutput::writeAxis_(CAxis* axis)",
1308                    << "[ type = " << SuperClass::type << "]"
1309                    << " not implemented yet !");
[609]1310          }
[1030]1311        }
1312        catch (CNetCdfException& e)
1313        {
1314          StdString msg("On writing the axis : ");
1315          msg.append(axisid); msg.append("\n");
1316          msg.append("In the context : ");
1317          CContext* context = CContext::getCurrent() ;
1318          msg.append(context->getId()); msg.append("\n");
1319          msg.append(e.what());
1320          ERROR("CNc4DataOutput::writeAxis_(CAxis* axis)", << msg);
1321        }
[1158]1322        axis->addRelFile(this->filename);
[391]1323     }
[488]1324
[887]1325      void CNc4DataOutput::writeScalar_(CScalar* scalar)
1326      {
1327        if (scalar->IsWritten(this->filename)) return;
1328        scalar->checkAttributes();
1329        int scalarSize = 1;
1330
1331        StdString scalaId = scalar->getScalarOutputName();
[1430]1332        StdString boundsId;
[887]1333        if (isWrittenAxis(scalaId)) return ;
1334        else setWrittenAxis(scalaId);
1335
[1158]1336        nc_type typePrec ;
1337        if (scalar->prec.isEmpty()) typePrec =  NC_FLOAT ;
1338        else if (scalar->prec==4)  typePrec =  NC_FLOAT ;
1339        else if (scalar->prec==8)   typePrec =  NC_DOUBLE ;
1340
[1435]1341        if (!scalar->label.isEmpty()) typePrec = NC_CHAR ;
1342        string strId="str_len" ;
1343
[887]1344        try
1345        {
[1439]1346          if (!scalar->label.isEmpty() && !SuperClassWriter::dimExist(strId)) SuperClassWriter::addDimension(strId, stringArrayLen);
[1435]1347
1348          if (!scalar->value.isEmpty() || !scalar->label.isEmpty())
[887]1349          {
[1430]1350            std::vector<StdString> dims;
1351            StdString scalarDim = scalaId;
1352
[1435]1353            if (!scalar->label.isEmpty()) dims.push_back(strId);
1354
[1158]1355            SuperClassWriter::addVariable(scalaId, typePrec, dims);
[887]1356
1357            if (!scalar->name.isEmpty())
1358              SuperClassWriter::addAttribute("name", scalar->name.getValue(), &scalaId);
1359
1360            if (!scalar->standard_name.isEmpty())
1361              SuperClassWriter::addAttribute("standard_name", scalar->standard_name.getValue(), &scalaId);
1362
1363            if (!scalar->long_name.isEmpty())
1364              SuperClassWriter::addAttribute("long_name", scalar->long_name.getValue(), &scalaId);
1365
1366            if (!scalar->unit.isEmpty())
1367              SuperClassWriter::addAttribute("units", scalar->unit.getValue(), &scalaId);
1368
[1430]1369            if (!scalar->axis_type.isEmpty())
1370            {
1371              switch(scalar->axis_type)
1372              {
1373              case CScalar::axis_type_attr::X :
1374                SuperClassWriter::addAttribute("axis", string("X"), &scalaId);
1375                break;
1376              case CScalar::axis_type_attr::Y :
1377                SuperClassWriter::addAttribute("axis", string("Y"), &scalaId);
1378                break;
1379              case CScalar::axis_type_attr::Z :
1380                SuperClassWriter::addAttribute("axis", string("Z"), &scalaId);
1381                break;
1382              case CScalar::axis_type_attr::T :
1383                SuperClassWriter::addAttribute("axis", string("T"), &scalaId);
1384                break;
1385              }
1386            }
1387
[1435]1388            if (!scalar->positive.isEmpty())
[1430]1389            {
[1435]1390              SuperClassWriter::addAttribute("positive",
1391                                             (scalar->positive == CScalar::positive_attr::up) ? string("up") : string("down"),
1392                                             &scalaId);
1393            }
1394
1395            if (!scalar->bounds.isEmpty() && scalar->label.isEmpty())
1396            {
[1436]1397              dims.clear();
1398              dims.push_back("axis_nbounds");
[1435]1399              boundsId = (scalar->bounds_name.isEmpty()) ? (scalaId + "_bounds") : scalar->bounds_name.getValue();
[1430]1400              SuperClassWriter::addVariable(boundsId, typePrec, dims);
[1436]1401              SuperClassWriter::addAttribute("bounds", boundsId, &scalaId);
[1430]1402            }
1403
[887]1404            SuperClassWriter::definition_end();
1405
1406            switch (SuperClass::type)
1407            {
1408              case MULTI_FILE:
1409              {
1410                CArray<double,1> scalarValue(scalarSize);
[1435]1411                CArray<string,1> scalarLabel(scalarSize);
[1436]1412                CArray<double,1> scalarBounds(scalarSize*2);
[1435]1413
1414                if (!scalar->value.isEmpty() && scalar->label.isEmpty())
[1430]1415                {
[1435]1416                  scalarValue(0) = scalar->value;
1417                  SuperClassWriter::writeData(scalarValue, scalaId, isCollective, 0);
1418                }
1419
1420                if (!scalar->bounds.isEmpty() && scalar->label.isEmpty())
1421                {
[1436]1422                  scalarBounds(0) = scalar->bounds(0);
1423                  scalarBounds(1) = scalar->bounds(1);
1424                  SuperClassWriter::writeData(scalarBounds, boundsId, isCollective, 0);
[1430]1425                }
[1435]1426
1427                if (!scalar->label.isEmpty())
1428                {
1429                  scalarLabel(0) = scalar->label;
1430                  SuperClassWriter::writeData(scalarLabel, scalaId, isCollective, 0);
1431                }
1432
[887]1433                SuperClassWriter::definition_start();
1434
1435                break;
1436              }
1437              case ONE_FILE:
1438              {
1439                CArray<double,1> scalarValue(scalarSize);
[1435]1440                CArray<string,1> scalarLabel(scalarSize);
[1436]1441                CArray<double,1> scalarBounds(scalarSize*2);
[887]1442
1443                std::vector<StdSize> start(1);
1444                std::vector<StdSize> count(1);
1445                start[0] = 0;
1446                count[0] = 1;
[1435]1447                if (!scalar->value.isEmpty() && scalar->label.isEmpty())
[1430]1448                {
[1435]1449                  scalarValue(0) = scalar->value;
1450                  SuperClassWriter::writeData(scalarValue, scalaId, isCollective, 0, &start, &count);
1451                }
1452                if (!scalar->bounds.isEmpty() && scalar->label.isEmpty())
1453                {
[1436]1454                  scalarBounds(0) = scalar->bounds(0);
1455                  scalarBounds(1) = scalar->bounds(1);
1456                  count[0] = 2;
1457                  SuperClassWriter::writeData(scalarBounds, boundsId, isCollective, 0, &start, &count);
[1430]1458                }
[1435]1459                if (!scalar->label.isEmpty())
1460                {
1461                  scalarLabel(0) = scalar->label;
1462                  count[0] = stringArrayLen;
1463                  SuperClassWriter::writeData(scalarLabel, scalaId, isCollective, 0, &start, &count);
1464                }
1465
[887]1466                SuperClassWriter::definition_start();
1467
1468                break;
1469              }
1470              default :
1471                ERROR("CNc4DataOutput::writeAxis_(CAxis* scalar)",
1472                      << "[ type = " << SuperClass::type << "]"
1473                      << " not implemented yet !");
1474            }
1475          }
1476        }
1477        catch (CNetCdfException& e)
1478        {
1479          StdString msg("On writing the scalar : ");
1480          msg.append(scalaId); msg.append("\n");
1481          msg.append("In the context : ");
1482          CContext* context = CContext::getCurrent() ;
1483          msg.append(context->getId()); msg.append("\n");
1484          msg.append(e.what());
1485          ERROR("CNc4DataOutput::writeScalar_(CScalar* scalar)", << msg);
1486        }
1487        scalar->addRelFile(this->filename);
1488     }
1489
[676]1490     //--------------------------------------------------------------
1491
1492     void CNc4DataOutput::writeGridCompressed_(CGrid* grid)
1493     {
[1957]1494        if (grid->isScalarGrid() || grid->isWrittenCompressed(this->filename)) return;
1495       
1496        // NOTA : The cuurent algorithm to write compress elements of the grid
1497        //        will work pretting well when on server side you dont't get
1498        //        partial overlap on elements between differents participating process
1499        //        So the element must be totally distributed or non distributed
1500        //        If an element is partially overlaping betwwen process then the
1501        //        total compressed part will apear artificially greater than expected
1502        //        For the current implementation of writer which is decomposed only on
1503        //        one element, it will work as expected, but for future, it must be
1504        //        reconsidered again.
1505        try
1506        {
1507          CArray<int,1> axisDomainOrder = grid->axis_domain_order;
1508          std::vector<StdString> domainList = grid->getDomainList();
1509          std::vector<StdString> axisList   = grid->getAxisList();
1510          std::vector<StdString> scalarList = grid->getScalarList();
1511          int numElement = axisDomainOrder.numElements(), idxDomain = 0, idxAxis = 0, idxScalar = 0;
1512          int commRank ;
1513          MPI_Comm_rank(comm_file,&commRank) ;
[676]1514
[1957]1515          std::vector<StdString> dims;
[676]1516
[1957]1517          for (int i = 0; i < numElement; ++i)
1518          {
1519            StdString varId, compress;
1520            CArray<size_t, 1> indexes;
1521            bool isDistributed;
1522            size_t nbIndexes, totalNbIndexes, offset;
1523            size_t firstGlobalIndex;
1524           
1525            if (2 == axisDomainOrder(i))
1526            {
1527              CDomain* domain = CDomain::get(domainList[idxDomain]);
1528              StdString domId = domain->getDomainOutputName();
[676]1529
[1957]1530              if (!domain->isCompressible()
1531                  || domain->type == CDomain::type_attr::unstructured
1532                  || domain->isWrittenCompressed(this->filename)
1533                  || isWrittenCompressedDomain(domId))
1534                continue;
1535           
1536              // unstructured grid seems not be taken into account why ?
[676]1537
[1957]1538              string lonName,latName ;
[676]1539
[1957]1540              if (domain->lon_name.isEmpty())
1541              { 
1542                if (domain->type==CDomain::type_attr::curvilinear) lonName = "nav_lon";
1543                else lonName = "lon";
1544              }
1545              else lonName = domain->lon_name;
[676]1546
[1957]1547              if (domain->lat_name.isEmpty())
1548              {
1549                if (domain->type==CDomain::type_attr::curvilinear) latName = "nav_lat";
1550                else latName = "lat";
1551              }
1552              else latName = domain->lat_name;
1553             
1554              StdString appendDomId  = singleDomain ? "" : "_" + domId;
[676]1555
[1957]1556              varId = domId + "_points";
1557              compress = latName + appendDomId + " " + lonName + appendDomId;
1558     
[2267]1559              shared_ptr<CLocalView> workflowView = domain->getLocalView(CElementView::WORKFLOW) ;
[1957]1560              workflowView->getGlobalIndexView(indexes) ;
1561              nbIndexes = workflowView->getSize() ;
1562              isDistributed = domain->isDistributed();
1563              if (isDistributed)
1564              {
1565                MPI_Exscan(&nbIndexes, &offset, 1, MPI_SIZE_T, MPI_SUM, comm_file) ;
1566                if (commRank==0) offset=0 ;
1567                MPI_Allreduce(&nbIndexes,&totalNbIndexes,1 , MPI_SIZE_T, MPI_SUM, comm_file) ;
1568              }
1569              else
1570              {
1571                offset=0 ;
1572                totalNbIndexes = nbIndexes ;
1573              }
[676]1574
[1957]1575              firstGlobalIndex = domain->ibegin + domain->jbegin * domain->ni_glo;
[676]1576
[1957]1577              domain->addRelFileCompressed(this->filename);
1578              setWrittenCompressedDomain(domId);
1579              ++idxDomain;
1580            }
1581            else if (1 == axisDomainOrder(i))
1582            {
1583              CAxis* axis = CAxis::get(axisList[idxAxis]);
1584              StdString axisId = axis->getAxisOutputName();
[676]1585
[1957]1586              if (!axis->isCompressible()
1587                  || axis->isWrittenCompressed(this->filename)
1588                  || isWrittenCompressedAxis(axisId))
1589                continue;
[676]1590
[1957]1591              varId = axisId + "_points";
1592              compress = axisId;
[676]1593
[2267]1594              shared_ptr<CLocalView> workflowView = axis->getLocalView(CElementView::WORKFLOW) ;
[1957]1595              workflowView->getGlobalIndexView(indexes) ;
1596              nbIndexes = workflowView->getSize() ;
1597              isDistributed = axis->isDistributed();
1598              if (isDistributed)
1599              {
1600                MPI_Exscan(&nbIndexes, &offset, 1, MPI_SIZE_T, MPI_SUM, comm_file) ;
1601                if (commRank==0) offset=0 ;
1602                MPI_Allreduce(&nbIndexes,&totalNbIndexes,1 , MPI_SIZE_T, MPI_SUM, comm_file) ;
1603              }
1604              else
1605              {
1606                offset=0 ;
1607                totalNbIndexes = nbIndexes ;
1608              }
1609              firstGlobalIndex = axis->begin;
1610             
1611              axis->addRelFileCompressed(this->filename);
1612              setWrittenCompressedAxis(axisId);
1613              ++idxAxis;
1614            }
1615            else
1616            {
1617              //for scalar
1618            }
[676]1619
[1957]1620            if (!varId.empty())
1621            {
1622              SuperClassWriter::addDimension(varId, (SuperClass::type == MULTI_FILE) ? nbIndexes : totalNbIndexes);
[774]1623
[1957]1624              dims.clear();
1625              dims.push_back(varId);
1626              SuperClassWriter::addVariable(varId, NC_UINT64, dims);
[676]1627
[1957]1628              SuperClassWriter::addAttribute("compress", compress, &varId);
[676]1629
[1957]1630              switch (SuperClass::type)
1631              {
1632                case (MULTI_FILE):
1633                {
1634                  indexes -= firstGlobalIndex;
1635                  SuperClassWriter::writeData(indexes, varId, isCollective, 0);
1636                  break;
1637                }
1638                case (ONE_FILE):
1639                {
1640                  std::vector<StdSize> start, count;
1641                  start.push_back(offset);
1642                  count.push_back(nbIndexes);
[676]1643
[1957]1644                  SuperClassWriter::writeData(indexes, varId, isCollective, 0, &start, &count);
1645                  break;
1646                }
1647              }
[1144]1648            }
[1957]1649          }
[676]1650
[1957]1651          grid->addRelFileCompressed(this->filename);
1652        }
1653        catch (CNetCdfException& e)
1654        {
1655          StdString msg("On writing compressed grid : ");
1656          msg.append(grid->getId()); msg.append("\n");
1657          msg.append("In the context : ");
1658          CContext* context = CContext::getCurrent();
1659          msg.append(context->getId()); msg.append("\n");
1660          msg.append(e.what());
1661          ERROR("CNc4DataOutput::writeGridCompressed_(CGrid* grid)", << msg);
1662        }
1663      }
[676]1664
1665     //--------------------------------------------------------------
1666
[391]1667     void CNc4DataOutput::writeTimeDimension_(void)
1668     {
[498]1669       try
1670       {
[802]1671        SuperClassWriter::addDimension(getTimeCounterName());
[498]1672       }
1673       catch (CNetCdfException& e)
1674       {
[614]1675         StdString msg("On writing time dimension : time_couter\n");
[498]1676         msg.append("In the context : ");
1677         CContext* context = CContext::getCurrent() ;
1678         msg.append(context->getId()); msg.append("\n");
1679         msg.append(e.what());
1680         ERROR("CNc4DataOutput::writeTimeDimension_(void)", << msg);
1681       }
[391]1682     }
[676]1683
[219]1684      //--------------------------------------------------------------
1685
[347]1686      void CNc4DataOutput::writeField_(CField* field)
[219]1687      {
[1158]1688        CContext* context = CContext::getCurrent() ;
[300]1689
[1158]1690        std::vector<StdString> dims, coodinates;
[1869]1691        CGrid* grid = field->getGrid();
[1158]1692        if (!grid->doGridHaveDataToWrite())
[567]1693          if (SuperClass::type==MULTI_FILE) return ;
[488]1694
[1158]1695        CArray<int,1> axisDomainOrder = grid->axis_domain_order;
1696        int numElement = axisDomainOrder.numElements(), idxDomain = 0, idxAxis = 0, idxScalar = 0;
1697        std::vector<StdString> domainList = grid->getDomainList();
1698        std::vector<StdString> axisList   = grid->getAxisList();
1699        std::vector<StdString> scalarList = grid->getScalarList();       
[219]1700
[1158]1701        StdString timeid  = getTimeCounterName();
1702        StdString dimXid,dimYid;
1703        std::deque<StdString> dimIdList, dimCoordList;
1704        bool hasArea = false;
1705        StdString cellMeasures = "area:";
1706        bool compressedOutput = !field->indexed_output.isEmpty() && field->indexed_output;
[488]1707
[1158]1708        for (int i = 0; i < numElement; ++i)
1709        {
1710          if (2 == axisDomainOrder(i))
1711          {
1712            CDomain* domain = CDomain::get(domainList[idxDomain]);
1713            StdString domId = domain->getDomainOutputName();
1714            StdString appendDomId  = singleDomain ? "" : "_" + domId ;
[1391]1715            StdString lonName,latName ;
[1430]1716            StdString dimIname,dimJname ;
[676]1717
[1415]1718            if (domain->lon_name.isEmpty())
1719            { 
1720              if (domain->type==CDomain::type_attr::curvilinear) lonName = "nav_lon";
1721              else lonName = "lon";
1722            }
[1391]1723            else lonName = domain->lon_name;
1724
[1415]1725            if (domain->lat_name.isEmpty())
1726            {
1727              if (domain->type==CDomain::type_attr::curvilinear) latName = "nav_lat";
1728              else latName = "lat";
1729            }
[1391]1730            else latName = domain->lat_name;
[1430]1731
1732            if (domain->dim_i_name.isEmpty())
1733            {
1734              if (domain->type==CDomain::type_attr::curvilinear) dimIname = "x";
1735              else if (domain->type==CDomain::type_attr::unstructured) dimIname = "cell";
1736              else dimIname = lonName;
1737            }
1738            else dimIname = domain->dim_i_name;
1739
1740            if (domain->dim_j_name.isEmpty())
1741            {
1742              if (domain->type==CDomain::type_attr::curvilinear) dimJname = "y";
1743              else dimJname = latName;
1744            }
1745            else dimJname = domain->dim_j_name;
[1391]1746       
[1158]1747            if (compressedOutput && domain->isCompressible() && domain->type != CDomain::type_attr::unstructured)
1748            {
1749              dimIdList.push_back(domId + "_points");
1750              field->setUseCompressedOutput();
1751            }
[676]1752
[1158]1753            switch (domain->type)
[879]1754            {
[1158]1755              case CDomain::type_attr::curvilinear:
1756                if (!compressedOutput || !domain->isCompressible())
1757                {
[1430]1758                  dimXid=dimIname+appendDomId;
1759                  dimYid=dimJname+appendDomId;
[1158]1760                  dimIdList.push_back(dimXid);
1761                  dimIdList.push_back(dimYid);
1762                }
[1415]1763                dimCoordList.push_back(lonName+appendDomId);
1764                dimCoordList.push_back(latName+appendDomId);
[1158]1765              break ;
1766              case CDomain::type_attr::rectilinear:
1767                if (!compressedOutput || !domain->isCompressible())
1768                {
[1430]1769                  dimXid     = dimIname+appendDomId;
1770                  dimYid     = dimJname+appendDomId;
[1158]1771                  dimIdList.push_back(dimXid);
1772                  dimIdList.push_back(dimYid);
1773                }
[1430]1774                if (lonName != dimIname)  dimCoordList.push_back(lonName+appendDomId);
1775                if (latName != dimJname)  dimCoordList.push_back(latName+appendDomId);
1776
[1158]1777              break ;
1778              case CDomain::type_attr::unstructured:
1779              {
1780                if (SuperClassWriter::useCFConvention)
1781                {
[1430]1782                  dimXid     = dimIname + appendDomId;
[1158]1783                  dimIdList.push_back(dimXid);
[1415]1784                  dimCoordList.push_back(lonName+appendDomId);
1785                  dimCoordList.push_back(latName+appendDomId);
[1158]1786                }
1787                else
1788                {
1789                  StdString domainName = domain->name;
1790                  if (domain->nvertex == 1)
1791                  {
1792                    dimXid     = "n" + domainName + "_node";
1793                    dimIdList.push_back(dimXid);
1794                    dimCoordList.push_back(StdString(domainName + "_node_x"));
1795                    dimCoordList.push_back(StdString(domainName + "_node_y"));
1796                  }
1797                  else if (domain->nvertex == 2)
1798                  {
1799                    dimXid     = "n" + domainName + "_edge";
1800                    dimIdList.push_back(dimXid);
1801                    dimCoordList.push_back(StdString(domainName + "_edge_x"));
1802                    dimCoordList.push_back(StdString(domainName + "_edge_y"));
1803                  }
1804                  else
1805                  {
1806                    dimXid     = "n" + domainName + "_face";
1807                    dimIdList.push_back(dimXid);
1808                    dimCoordList.push_back(StdString(domainName + "_face_x"));
1809                    dimCoordList.push_back(StdString(domainName + "_face_y"));
1810                  }
1811                }  // ugrid convention
1812              }  // case unstructured domain
[879]1813            }
[1158]1814
1815            if (domain->hasArea)
[879]1816            {
[1158]1817              hasArea = true;
1818              cellMeasures += " area" + appendDomId;
[879]1819            }
[1158]1820            ++idxDomain;
1821          }
1822          else if (1 == axisDomainOrder(i))
1823          {
1824            CAxis* axis = CAxis::get(axisList[idxAxis]);
1825            StdString axisId = axis->getAxisOutputName();
[1430]1826            StdString axisDim;
[1158]1827
[1430]1828            if (axis->dim_name.isEmpty()) axisDim = axisId;
1829            else axisDim=axis->dim_name.getValue();
1830
[1158]1831            if (compressedOutput && axis->isCompressible())
[879]1832            {
[1430]1833              dimIdList.push_back(axisDim + "_points");
[1158]1834              field->setUseCompressedOutput();
[879]1835            }
[1158]1836            else
[1430]1837              dimIdList.push_back(axisDim);
[878]1838
[1430]1839            if (axisDim != axisId) dimCoordList.push_back(axisId);
[1158]1840            ++idxAxis;
1841          }
[1430]1842          else
[1158]1843          {
[1430]1844            CScalar* scalar = CScalar::get(scalarList[idxScalar]);
1845            StdString scalarId = scalar->getScalarOutputName();
[1446]1846            if (!scalar->value.isEmpty() || !scalar->label.isEmpty())
1847              dimCoordList.push_back(scalarId);
[1430]1848            ++idxScalar;
[1158]1849          }
1850        }
[488]1851
[1158]1852        StdString fieldid = field->getFieldOutputName();
[219]1853
[1158]1854        nc_type type ;
1855        if (field->prec.isEmpty()) type =  NC_FLOAT ;
1856        else
1857        {
1858          if (field->prec==2) type = NC_SHORT ;
1859          else if (field->prec==4)  type =  NC_FLOAT ;
1860          else if (field->prec==8)   type =  NC_DOUBLE ;
1861        }
[488]1862
[1158]1863        bool wtime   = !(!field->operation.isEmpty() && field->getOperationTimeType() == func::CFunctor::once);
[488]1864
[1158]1865        if (wtime)
1866        {
1867          if (field->hasTimeInstant && hasTimeInstant) coodinates.push_back(string("time_instant"));
1868          else if (field->hasTimeCentered && hasTimeCentered)  coodinates.push_back(string("time_centered"));
1869          dims.push_back(timeid);
1870        }
[488]1871
[1961]1872        while (!dimIdList.empty())
[1158]1873        {
[1961]1874          dims.push_back(dimIdList.back());
1875          dimIdList.pop_back();
[1158]1876        }
[219]1877
[1158]1878        while (!dimCoordList.empty())
1879        {
1880          coodinates.push_back(dimCoordList.back());
1881          dimCoordList.pop_back();
1882        }
[219]1883
[1158]1884        try
1885        {
[498]1886           SuperClassWriter::addVariable(fieldid, type, dims);
[488]1887
[498]1888           if (!field->standard_name.isEmpty())
1889              SuperClassWriter::addAttribute
1890                 ("standard_name",  field->standard_name.getValue(), &fieldid);
[219]1891
[498]1892           if (!field->long_name.isEmpty())
1893              SuperClassWriter::addAttribute
1894                 ("long_name", field->long_name.getValue(), &fieldid);
[219]1895
[498]1896           if (!field->unit.isEmpty())
1897              SuperClassWriter::addAttribute
1898                 ("units", field->unit.getValue(), &fieldid);
[463]1899
[1158]1900           // Ugrid field attributes "mesh" and "location"
1901           if (!SuperClassWriter::useCFConvention)
1902           {
1903            if (!domainList.empty())
1904            {
1905              CDomain* domain = CDomain::get(domainList[0]); // Suppose that we have only domain
1906              StdString mesh = domain->name;
1907              SuperClassWriter::addAttribute("mesh", mesh, &fieldid);
1908              StdString location;
1909              if (domain->nvertex == 1)
1910                location = "node";
1911              else if (domain->nvertex == 2)
1912                location = "edge";
1913              else if (domain->nvertex > 2)
1914                location = "face";
1915              SuperClassWriter::addAttribute("location", location, &fieldid);
1916            }
1917
1918           }
1919
1920           if (!field->valid_min.isEmpty())
[498]1921              SuperClassWriter::addAttribute
1922                 ("valid_min", field->valid_min.getValue(), &fieldid);
[463]1923
[498]1924           if (!field->valid_max.isEmpty())
1925              SuperClassWriter::addAttribute
1926                 ("valid_max", field->valid_max.getValue(), &fieldid);
[464]1927
[498]1928            if (!field->scale_factor.isEmpty())
1929              SuperClassWriter::addAttribute
1930                 ("scale_factor", field->scale_factor.getValue(), &fieldid);
[464]1931
[498]1932             if (!field->add_offset.isEmpty())
1933              SuperClassWriter::addAttribute
1934                 ("add_offset", field->add_offset.getValue(), &fieldid);
[488]1935
[498]1936           SuperClassWriter::addAttribute
1937                 ("online_operation", field->operation.getValue(), &fieldid);
[472]1938
[498]1939          // write child variables as attributes
[488]1940
1941
[1021]1942           bool alreadyAddCellMethod = false;
1943           StdString cellMethodsPrefix(""), cellMethodsSuffix("");
1944           if (!field->cell_methods.isEmpty())
1945           {
[1158]1946              StdString cellMethodString = field->cell_methods;
1947              if (field->cell_methods_mode.isEmpty() ||
[1021]1948                 (CField::cell_methods_mode_attr::overwrite == field->cell_methods_mode))
[1158]1949              {
1950                SuperClassWriter::addAttribute("cell_methods", cellMethodString, &fieldid);
1951                alreadyAddCellMethod = true;
1952              }
1953              else
1954              {
1955                switch (field->cell_methods_mode)
1956                {
1957                  case (CField::cell_methods_mode_attr::prefix):
1958                    cellMethodsPrefix = cellMethodString;
1959                    cellMethodsPrefix += " ";
1960                    break;
1961                  case (CField::cell_methods_mode_attr::suffix):
1962                    cellMethodsSuffix = " ";
1963                    cellMethodsSuffix += cellMethodString;
1964                    break;
1965                  case (CField::cell_methods_mode_attr::none):
1966                    break;
1967                  default:
1968                    break;
1969                }
1970              }
[1021]1971           }
[488]1972
[1021]1973
[498]1974           if (wtime)
1975           {
[614]1976              CDuration freqOp = field->freq_op.getValue();
1977              freqOp.solveTimeStep(*context->calendar);
1978              StdString freqOpStr = freqOp.toStringUDUnits();
[612]1979              SuperClassWriter::addAttribute("interval_operation", freqOpStr, &fieldid);
[437]1980
[614]1981              CDuration freqOut = field->getRelFile()->output_freq.getValue();
1982              freqOut.solveTimeStep(*context->calendar);
1983              SuperClassWriter::addAttribute("interval_write", freqOut.toStringUDUnits(), &fieldid);
[612]1984
[1021]1985              StdString cellMethods(cellMethodsPrefix + "time: ");
[612]1986              if (field->operation.getValue() == "instant") cellMethods += "point";
1987              else if (field->operation.getValue() == "average") cellMethods += "mean";
1988              else if (field->operation.getValue() == "accumulate") cellMethods += "sum";
1989              else cellMethods += field->operation;
[614]1990              if (freqOp.resolve(*context->calendar) != freqOut.resolve(*context->calendar))
1991                cellMethods += " (interval: " + freqOpStr + ")";
[1021]1992              cellMethods += cellMethodsSuffix;
1993              if (!alreadyAddCellMethod)
1994                SuperClassWriter::addAttribute("cell_methods", cellMethods, &fieldid);
[498]1995           }
[488]1996
[611]1997           if (hasArea)
1998             SuperClassWriter::addAttribute("cell_measures", cellMeasures, &fieldid);
1999
[498]2000           if (!field->default_value.isEmpty())
2001           {
[1424]2002             double default_value = field->default_value.getValue();
2003             if (type == NC_DOUBLE)
2004             {
2005               SuperClassWriter::setDefaultValue(fieldid, &default_value);
2006             }
2007             else if (type == NC_SHORT)
2008             {
2009               short sdefault_value = (short)default_value;
2010               SuperClassWriter::setDefaultValue(fieldid, &sdefault_value);
2011             }
2012             else
2013             {
2014               float fdefault_value = (float)default_value;
2015               SuperClassWriter::setDefaultValue(fieldid, &fdefault_value);
2016             }
[498]2017           }
2018           else
[517]2019              SuperClassWriter::setDefaultValue(fieldid, (double*)NULL);
[219]2020
[606]2021            if (field->compression_level.isEmpty())
[1872]2022              field->compression_level = field->getRelFile()->compression_level.isEmpty() ? 0 : field->getRelFile()->compression_level;
[606]2023            SuperClassWriter::setCompressionLevel(fieldid, field->compression_level);
2024
[878]2025           {  // Ecriture des coordonnes
[488]2026
[498]2027              StdString coordstr; //boost::algorithm::join(coodinates, " ")
2028              std::vector<StdString>::iterator
2029                 itc = coodinates.begin(), endc = coodinates.end();
[488]2030
[498]2031              for (; itc!= endc; itc++)
2032              {
2033                 StdString & coord = *itc;
2034                 if (itc+1 != endc)
2035                       coordstr.append(coord).append(" ");
2036                 else  coordstr.append(coord);
2037              }
[219]2038
[498]2039              SuperClassWriter::addAttribute("coordinates", coordstr, &fieldid);
[219]2040
[498]2041           }
[1222]2042
2043           vector<CVariable*> listVars = field->getAllVariables() ;
2044           for (vector<CVariable*>::iterator it = listVars.begin() ;it != listVars.end(); it++) writeAttribute_(*it, fieldid) ;
2045
[219]2046         }
[498]2047         catch (CNetCdfException& e)
2048         {
2049           StdString msg("On writing field : ");
2050           msg.append(fieldid); msg.append("\n");
2051           msg.append("In the context : ");
2052           msg.append(context->getId()); msg.append("\n");
2053           msg.append(e.what());
2054           ERROR("CNc4DataOutput::writeField_(CField* field)", << msg);
2055         }
[879]2056      } // writeField_()
[219]2057
2058      //--------------------------------------------------------------
2059
[347]2060      void CNc4DataOutput::writeFile_ (CFile* file)
[219]2061      {
[773]2062         StdString filename = file->getFileOutputName();
[219]2063         StdString description = (!file->description.isEmpty())
2064                               ? file->description.getValue()
[335]2065                               : StdString("Created by xios");
[609]2066
2067         singleDomain = (file->nbDomains == 1);
2068
[1158]2069         StdString conv_str ;
2070         if (file->convention_str.isEmpty())
2071         {
2072            if (SuperClassWriter::useCFConvention) conv_str="CF-1.6" ;
2073            else conv_str="UGRID" ;
2074         }
2075         else conv_str=file->convention_str ;
2076           
[498]2077         try
2078         {
[1309]2079           if (!appendMode) this->writeFileAttributes(filename, description,
2080                                                      conv_str,
2081                                                      StdString("An IPSL model"),
2082                                                      this->getTimeStamp());
[609]2083
[701]2084           if (!appendMode)
2085             SuperClassWriter::addDimension("axis_nbounds", 2);
[498]2086         }
2087         catch (CNetCdfException& e)
2088         {
2089           StdString msg("On writing file : ");
2090           msg.append(filename); msg.append("\n");
2091           msg.append("In the context : ");
2092           CContext* context = CContext::getCurrent() ;
2093           msg.append(context->getId()); msg.append("\n");
2094           msg.append(e.what());
2095           ERROR("CNc4DataOutput::writeFile_ (CFile* file)", << msg);
2096         }
[219]2097      }
[488]2098
[472]2099      void CNc4DataOutput::writeAttribute_ (CVariable* var, const string& fieldId)
2100      {
[773]2101        StdString name = var->getVariableOutputName();
[488]2102
[498]2103        try
2104        {
[527]2105          if (var->type.getValue() == CVariable::type_attr::t_int || var->type.getValue() == CVariable::type_attr::t_int32)
2106            addAttribute(name, var->getData<int>(), &fieldId);
2107          else if (var->type.getValue() == CVariable::type_attr::t_int16)
2108            addAttribute(name, var->getData<short int>(), &fieldId);
2109          else if (var->type.getValue() == CVariable::type_attr::t_float)
2110            addAttribute(name, var->getData<float>(), &fieldId);
2111          else if (var->type.getValue() == CVariable::type_attr::t_double)
2112            addAttribute(name, var->getData<double>(), &fieldId);
2113          else if (var->type.getValue() == CVariable::type_attr::t_string)
2114            addAttribute(name, var->getData<string>(), &fieldId);
2115          else
2116            ERROR("CNc4DataOutput::writeAttribute_ (CVariable* var, const string& fieldId)",
2117                  << "Unsupported variable of type " << var->type.getStringValue());
[498]2118        }
2119       catch (CNetCdfException& e)
2120       {
2121         StdString msg("On writing attributes of variable with name : ");
2122         msg.append(name); msg.append("in the field "); msg.append(fieldId); msg.append("\n");
2123         msg.append("In the context : ");
2124         CContext* context = CContext::getCurrent() ;
2125         msg.append(context->getId()); msg.append("\n");
2126         msg.append(e.what());
2127         ERROR("CNc4DataOutput::writeAttribute_ (CVariable* var, const string& fieldId)", << msg);
2128       }
[472]2129     }
[488]2130
[472]2131     void CNc4DataOutput::writeAttribute_ (CVariable* var)
2132     {
[773]2133        StdString name = var->getVariableOutputName();
2134
[498]2135        try
2136        {
[527]2137          if (var->type.getValue() == CVariable::type_attr::t_int || var->type.getValue() == CVariable::type_attr::t_int32)
2138            addAttribute(name, var->getData<int>());
2139          else if (var->type.getValue() == CVariable::type_attr::t_int16)
2140            addAttribute(name, var->getData<short int>());
2141          else if (var->type.getValue() == CVariable::type_attr::t_float)
2142            addAttribute(name, var->getData<float>());
2143          else if (var->type.getValue() == CVariable::type_attr::t_double)
2144            addAttribute(name, var->getData<double>());
2145          else if (var->type.getValue() == CVariable::type_attr::t_string)
2146            addAttribute(name, var->getData<string>());
2147          else
2148            ERROR("CNc4DataOutput::writeAttribute_ (CVariable* var)",
2149                  << "Unsupported variable of type " << var->type.getStringValue());
[498]2150        }
2151       catch (CNetCdfException& e)
2152       {
2153         StdString msg("On writing attributes of variable with name : ");
2154         msg.append(name); msg.append("\n");
2155         msg.append("In the context : ");
2156         CContext* context = CContext::getCurrent() ;
2157         msg.append(context->getId()); msg.append("\n");
2158         msg.append(e.what());
2159         ERROR("CNc4DataOutput::writeAttribute_ (CVariable* var)", << msg);
2160       }
[488]2161     }
2162
[321]2163      void CNc4DataOutput::syncFile_ (void)
2164      {
[498]2165        try
2166        {
2167          SuperClassWriter::sync() ;
2168        }
2169        catch (CNetCdfException& e)
2170        {
2171         StdString msg("On synchronizing the write among processes");
2172         msg.append("In the context : ");
2173         CContext* context = CContext::getCurrent() ;
2174         msg.append(context->getId()); msg.append("\n");
2175         msg.append(e.what());
2176         ERROR("CNc4DataOutput::syncFile_ (void)", << msg);
2177        }
[321]2178      }
[219]2179
[286]2180      void CNc4DataOutput::closeFile_ (void)
2181      {
[498]2182        try
2183        {
2184          SuperClassWriter::close() ;
2185        }
2186        catch (CNetCdfException& e)
2187        {
2188         StdString msg("On closing file");
2189         msg.append("In the context : ");
2190         CContext* context = CContext::getCurrent() ;
2191         msg.append(context->getId()); msg.append("\n");
2192         msg.append(e.what());
2193         ERROR("CNc4DataOutput::syncFile_ (void)", << msg);
2194        }
2195
[286]2196      }
2197
[219]2198      //---------------------------------------------------------------
2199
2200      StdString CNc4DataOutput::getTimeStamp(void) const
2201      {
2202         const int buffer_size = 100;
2203         time_t rawtime;
2204         struct tm * timeinfo = NULL;
2205         char buffer [buffer_size];
[1158]2206         StdString formatStr;
2207         if (file->time_stamp_format.isEmpty()) formatStr="%Y-%b-%d %H:%M:%S %Z" ;
2208         else formatStr=file->time_stamp_format;
[219]2209
[1158]2210//         time ( &rawtime );
2211//         timeinfo = localtime ( &rawtime );
[219]2212         time ( &rawtime );
[1158]2213         timeinfo = gmtime ( &rawtime );
2214         strftime (buffer, buffer_size, formatStr.c_str(), timeinfo);
[219]2215
2216         return (StdString(buffer));
2217      }
[488]2218
[219]2219      //---------------------------------------------------------------
[488]2220
[1961]2221      int CNc4DataOutput::writeFieldData_ (CField*  field, const CArray<double,1>& data, const CDate& lastWrite,
2222                                           const CDate& currentWrite, int nstep)
[219]2223      {
[707]2224        CContext* context = CContext::getCurrent();
[1869]2225        CGrid* grid = field->getGrid();
[1961]2226       
2227        if (nstep<1) 
[1158]2228        {
[1961]2229          return nstep;
[1158]2230        }
[952]2231       
[707]2232        if (!grid->doGridHaveDataToWrite())
[1158]2233          if (SuperClass::type == MULTI_FILE || !isCollective)
2234          {
[1961]2235            return nstep;
[1158]2236          }
[488]2237
[770]2238        StdString fieldid = field->getFieldOutputName();
[286]2239
[707]2240        StdOStringStream oss;
2241        string timeAxisId;
[1158]2242        if (field->hasTimeInstant) timeAxisId = "time_instant";
2243        else if (field->hasTimeCentered) timeAxisId = "time_centered";
[488]2244
[802]2245        StdString timeBoundId = getTimeCounterName() + "_bounds";
[449]2246
[707]2247        StdString timeAxisBoundId;
[1158]2248        if (field->hasTimeInstant) timeAxisBoundId = "time_instant_bounds";
2249        else if (field->hasTimeCentered) timeAxisBoundId = "time_centered_bounds";
[488]2250
[707]2251        if (!field->wasWritten())
2252        {
[1872]2253          if (appendMode && field->getRelFile()->record_offset.isEmpty() && 
[1158]2254              field->getOperationTimeType() != func::CFunctor::once)
[707]2255          {
[1158]2256            double factorUnit;
[1872]2257            if (!field->getRelFile()->time_units.isEmpty() && field->getRelFile()->time_units==CFile::time_units_attr::days)
[1158]2258            factorUnit=context->getCalendar()->getDayLengthInSeconds() ;
2259            else factorUnit=1 ;
[1962]2260            nstep = getRecordFromTime(currentWrite,factorUnit) + 1;
[707]2261          }
[488]2262
[707]2263          field->setWritten();
2264        }
[488]2265
2266
[707]2267        CArray<double,1> time_data(1);
2268        CArray<double,1> time_data_bound(2);
2269        CArray<double,1> time_counter(1);
2270        CArray<double,1> time_counter_bound(2);
2271
2272        bool wtime = (field->getOperationTimeType() != func::CFunctor::once);
[1158]2273        bool wtimeCounter =false ;
2274        bool wtimeData =false ;
2275       
[707]2276
[444]2277        if (wtime)
2278        {
[1158]2279         
2280          if (field->hasTimeInstant)
2281          {
[1961]2282            time_data(0) = time_data_bound(1) = (Time) lastWrite;
2283            time_data_bound(0) = time_data_bound(1) = (Time) currentWrite;
[1158]2284            if (timeCounterType==instant)
2285            {
2286              time_counter(0) = time_data(0);
2287              time_counter_bound(0) = time_data_bound(0);
2288              time_counter_bound(1) = time_data_bound(1);
2289              wtimeCounter=true ;
2290            }
2291            if (hasTimeInstant) wtimeData=true ;
2292          }
2293          else if (field->hasTimeCentered)
[488]2294          {
[1961]2295            time_data(0) = ((Time)currentWrite + (Time)lastWrite) / 2;
2296            time_data_bound(0) = (Time)lastWrite;
2297            time_data_bound(1) = (Time)currentWrite;
[1158]2298            if (timeCounterType==centered)
2299            {
2300              time_counter(0) = time_data(0) ;
2301              time_counter_bound(0) = time_data_bound(0) ;
2302              time_counter_bound(1) = time_data_bound(1) ;
2303              wtimeCounter=true ;
2304            }
2305            if (hasTimeCentered) wtimeData=true ;
[488]2306          }
2307
[1158]2308          if (timeCounterType==record)
2309          {
[1961]2310            time_counter(0) = nstep - 1;
2311            time_counter_bound(0) = time_counter_bound(1) = nstep - 1;
[1158]2312            wtimeCounter=true ;
2313          }
[692]2314
[1872]2315          if (!field->getRelFile()->time_units.isEmpty() && field->getRelFile()->time_units==CFile::time_units_attr::days)
[692]2316          {
[1158]2317            double secByDay=context->getCalendar()->getDayLengthInSeconds() ;
2318            time_data/=secByDay;
2319            time_data_bound/=secByDay;
2320            time_counter/=secByDay;
2321            time_counter_bound/=secByDay;
[692]2322          }
2323        }
2324
[1853]2325         bool isRoot = (context->intraCommRank_ == 0);
[488]2326
[498]2327         try
[219]2328         {
[567]2329           switch (SuperClass::type)
[498]2330           {
[567]2331              case (MULTI_FILE) :
2332              {
[1158]2333                 CTimer::get("Files : writing data").resume();
[1961]2334                 SuperClassWriter::writeData(data, fieldid, isCollective, nstep - 1);
[1158]2335                 CTimer::get("Files : writing data").suspend();
[567]2336                 if (wtime)
2337                 {
[1158]2338                   CTimer::get("Files : writing time axis").resume();
2339                   if ( wtimeData)
[692]2340                   {
[1961]2341                       SuperClassWriter::writeTimeAxisData(time_data, timeAxisId, isCollective, nstep - 1, isRoot);
2342                       SuperClassWriter::writeTimeAxisDataBounds(time_data_bound, timeAxisBoundId, isCollective, nstep - 1, isRoot);
[1158]2343                  }
2344                   if (wtimeCounter)
2345                   {
[1961]2346                     SuperClassWriter::writeTimeAxisData(time_counter, getTimeCounterName(), isCollective, nstep - 1,isRoot);
2347                     if (timeCounterType!=record) SuperClassWriter::writeTimeAxisDataBounds(time_counter_bound, timeBoundId, isCollective, nstep - 1, isRoot);
[692]2348                   }
[1158]2349                   CTimer::get("Files : writing time axis").suspend();
[567]2350                 }
[707]2351                 break;
[567]2352              }
2353              case (ONE_FILE) :
2354              {
[464]2355
[676]2356                std::vector<StdSize> start, count;
[464]2357
[676]2358                if (field->getUseCompressedOutput())
2359                {
[1957]2360                  CArray<int,1> axisDomainOrder = grid->axis_domain_order;
2361                  std::vector<StdString> domainList = grid->getDomainList();
2362                  std::vector<StdString> axisList   = grid->getAxisList();
2363                  int numElement = axisDomainOrder.numElements();
2364                  int idxDomain = domainList.size() - 1, idxAxis = axisList.size() - 1;
2365                  int idx = domainList.size() * 2 + axisList.size() - 1;
2366                  int commRank ;
[464]2367
[1957]2368                  MPI_Comm_rank(comm_file,&commRank) ;
[676]2369
[1957]2370                  start.reserve(idx+1);
2371                  count.reserve(idx+1);
2372
2373                  for (int i = numElement - 1; i >= 0; --i)
2374                  {
2375                    if (2 == axisDomainOrder(i))
[676]2376                    {
[1957]2377                      CDomain* domain = CDomain::get(domainList[idxDomain]);
2378
2379                      if (domain->isCompressible())
[676]2380                      {
[1957]2381                        size_t offset ;
2382                        size_t nbIndexes = domain->getLocalView(CElementView::WORKFLOW)->getSize() ;
2383                        if (domain->isDistributed())
[676]2384                        {
[1957]2385                          MPI_Exscan(&nbIndexes, &offset, 1, MPI_SIZE_T, MPI_SUM, comm_file) ;
2386                          if (commRank==0) offset=0 ;
[676]2387                        }
[1957]2388                        else offset=0 ;
2389
2390                        start.push_back(offset);
2391                        count.push_back(nbIndexes);
2392                        idx -= 2;
2393                      }
2394                      else
2395                      {
2396                        if ((domain->type) != CDomain::type_attr::unstructured)
[676]2397                        {
[1957]2398                          start.push_back(domain->jbegin);
2399                          count.push_back(domain->nj);
[676]2400                        }
[1957]2401                        --idx;
2402                        start.push_back(domain->ibegin);
2403                        count.push_back(domain->ni);
2404                        --idx;
[676]2405                      }
[1957]2406                      --idxDomain;
2407                    }
2408                    else if (1 == axisDomainOrder(i))
2409                    {
2410                      CAxis* axis = CAxis::get(axisList[idxAxis]);
2411
2412                      if (axis->isCompressible())
[676]2413                      {
[1957]2414                        size_t offset ;
2415                        size_t nbIndexes = axis->getLocalView(CElementView::WORKFLOW)->getSize() ;
2416                        if (axis->isDistributed())
[676]2417                        {
[1957]2418                          MPI_Exscan(&nbIndexes, &offset, 1, MPI_SIZE_T, MPI_SUM, comm_file) ;
2419                          if (commRank==0) offset=0 ;
[676]2420                        }
[1957]2421                        else offset=0 ;
2422
2423                        start.push_back(offset);
2424                        count.push_back(nbIndexes);
[676]2425                      }
[1957]2426                      else
2427                      {
2428                        start.push_back(axis->begin);
2429                        count.push_back(axis->n);
2430                      }
2431                      --idxAxis;
2432                      --idx;
[676]2433                    }
[1158]2434                  }
[1961]2435                }
[676]2436                else
[498]2437                {
[887]2438                  CArray<int,1> axisDomainOrder = grid->axis_domain_order;
[705]2439                  std::vector<StdString> domainList = grid->getDomainList();
2440                  std::vector<StdString> axisList   = grid->getAxisList();
[2264]2441                  std::vector<StdString> scalarList  = grid->getScalarList() ;
[705]2442                  int numElement = axisDomainOrder.numElements();
2443                  int idxDomain = domainList.size() - 1, idxAxis = axisList.size() - 1;
[1553]2444                  int idx = domainList.size() * 2 + axisList.size() - 1;
[705]2445
[1553]2446                  start.reserve(idx+1);
2447                  count.reserve(idx+1);
[705]2448
2449                  for (int i = numElement - 1; i >= 0; --i)
2450                  {
[887]2451                    if (2 == axisDomainOrder(i))
[705]2452                    {
2453                      CDomain* domain = CDomain::get(domainList[idxDomain]);
[710]2454                      if ((domain->type) != CDomain::type_attr::unstructured)
[705]2455                      {
[1553]2456                        start.push_back(domain->jbegin);
2457                        count.push_back(domain->nj);
[705]2458                      }
2459                      --idx ;
[1143]2460
[1553]2461                        start.push_back(domain->ibegin);
2462                        count.push_back(domain->ni);
[705]2463                      --idx ;
2464                      --idxDomain;
2465                    }
[887]2466                    else if (1 == axisDomainOrder(i))
[705]2467                    {
[1143]2468                        CAxis* axis = CAxis::get(axisList[idxAxis]);
[1559]2469                        start.push_back(axis->begin);
2470                        count.push_back(axis->n);
[705]2471                      --idx;
[1025]2472                      --idxAxis;
[887]2473                    }
2474                    else
2475                    {
2476                      if (1 == axisDomainOrder.numElements())
2477                      {
[2264]2478                        CScalar* scalar = CScalar::get(scalarList[scalarList.size()-1]);
[887]2479                        start.push_back(0);
[2264]2480                        count.push_back(scalar->n);
[887]2481                      }
2482                      --idx;
2483                    }
[705]2484                  }
[498]2485                }
[488]2486
[1158]2487
2488                CTimer::get("Files : writing data").resume();
[1961]2489                SuperClassWriter::writeData(data, fieldid, isCollective, nstep - 1, &start, &count);
[1158]2490                CTimer::get("Files : writing data").suspend();
2491
2492                 if (wtime)
2493                 {
2494                   CTimer::get("Files : writing time axis").resume();
2495                   if ( wtimeData)
[692]2496                   {
[1961]2497                     SuperClassWriter::writeTimeAxisData(time_data, timeAxisId, isCollective, nstep - 1, isRoot);
2498                     SuperClassWriter::writeTimeAxisDataBounds(time_data_bound, timeAxisBoundId, isCollective, nstep - 1, isRoot);
[692]2499                   }
[1158]2500                   if (wtimeCounter)
2501                   {
[1961]2502                     SuperClassWriter::writeTimeAxisData(time_counter, getTimeCounterName(), isCollective, nstep - 1,isRoot);
2503                     if (timeCounterType!=record) SuperClassWriter::writeTimeAxisDataBounds(time_counter_bound, timeBoundId, isCollective, nstep - 1, isRoot);
[586]2504
[1158]2505                   }
2506                   CTimer::get("Files : writing time axis").suspend(); 
2507                 }
2508
[586]2509                break;
[286]2510              }
[567]2511            }
[1961]2512            return nstep ;
[219]2513         }
[498]2514         catch (CNetCdfException& e)
2515         {
2516           StdString msg("On writing field data: ");
2517           msg.append(fieldid); msg.append("\n");
2518           msg.append("In the context : ");
2519           msg.append(context->getId()); msg.append("\n");
2520           msg.append(e.what());
[1130]2521           ERROR("CNc4DataOutput::writeFieldData_ (CField*  field)", << msg);
[498]2522         }
[219]2523      }
2524
2525      //---------------------------------------------------------------
2526
2527      void CNc4DataOutput::writeTimeAxis_
[347]2528                  (CField*    field,
[1542]2529                   const std::shared_ptr<CCalendar> cal)
[219]2530      {
2531         StdOStringStream oss;
[1158]2532         bool createInstantAxis=false ;
2533         bool createCenteredAxis=false ;
2534         bool createTimeCounterAxis=false ;
2535         
[645]2536         if (field->getOperationTimeType() == func::CFunctor::once) return ;
[488]2537
2538
[1158]2539         StdString axisId ;
2540         StdString axisBoundId;
[802]2541         StdString timeid(getTimeCounterName());
[614]2542         StdString timeBoundId("axis_nbounds");
[488]2543
[1158]2544         StdString strTimeUnits ;
[1872]2545         if (!field->getRelFile()->time_units.isEmpty() && field->getRelFile()->time_units==CFile::time_units_attr::days) strTimeUnits="days since " ;
[1158]2546         else  strTimeUnits="seconds since " ;
2547 
2548         if (field->getOperationTimeType() == func::CFunctor::instant) field->hasTimeInstant = true;
2549         if (field->getOperationTimeType() == func::CFunctor::centered) field->hasTimeCentered = true;
2550
2551
[1872]2552         if (field->getRelFile()->time_counter.isEmpty())
[488]2553         {
[1158]2554           if (timeCounterType==none) createTimeCounterAxis=true ;
2555           if (field->hasTimeCentered)
2556           {
2557             timeCounterType=centered ;
2558             if (!hasTimeCentered) createCenteredAxis=true ;
2559           }
2560           if (field->hasTimeInstant)
2561           {
2562             if (timeCounterType==none) timeCounterType=instant ;
2563             if (!hasTimeInstant) createInstantAxis=true ;
2564           }
[488]2565         }
[1872]2566         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::instant)
[1158]2567         {
2568           if (field->hasTimeCentered)
2569           {
2570             if (!hasTimeCentered) createCenteredAxis=true ;
2571           }
2572           if (field->hasTimeInstant)
2573           {
2574             if (timeCounterType==none) createTimeCounterAxis=true ;
2575             timeCounterType=instant ;
2576             if (!hasTimeInstant) createInstantAxis=true ;
2577           }
2578         }
[1872]2579         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::centered)
[1158]2580         {
2581           if (field->hasTimeCentered)
2582           {
2583             if (timeCounterType==none) createTimeCounterAxis=true ;
2584             timeCounterType=centered ;
2585             if (!hasTimeCentered) createCenteredAxis=true ;
2586           }
2587           if (field->hasTimeInstant)
2588           {
2589             if (!hasTimeInstant) createInstantAxis=true ;
2590           }
2591         }
[1872]2592         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::instant_exclusive)
[1158]2593         {
2594           if (field->hasTimeCentered)
2595           {
2596             if (!hasTimeCentered) createCenteredAxis=true ;
2597           }
2598           if (field->hasTimeInstant)
2599           {
2600             if (timeCounterType==none) createTimeCounterAxis=true ;
2601             timeCounterType=instant ;
2602           }
2603         }
[1872]2604         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::centered_exclusive)
[1158]2605         {
2606           if (field->hasTimeCentered)
2607           {
2608             if (timeCounterType==none) createTimeCounterAxis=true ;
2609             timeCounterType=centered ;
2610           }
2611           if (field->hasTimeInstant)
2612           {
2613             if (!hasTimeInstant) createInstantAxis=true ;
2614           }
2615         }
[1872]2616         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::exclusive)
[1158]2617         {
2618           if (field->hasTimeCentered)
2619           {
2620             if (timeCounterType==none) createTimeCounterAxis=true ;
2621             if (timeCounterType==instant) createInstantAxis=true ;
2622             timeCounterType=centered ;
2623           }
2624           if (field->hasTimeInstant)
2625           {
2626             if (timeCounterType==none)
2627             {
2628               createTimeCounterAxis=true ;
2629               timeCounterType=instant ;
2630             }
2631             if (timeCounterType==centered)
2632             {
2633               if (!hasTimeInstant) createInstantAxis=true ;
2634             }
2635           }
2636         }
[1872]2637         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::none)
[1158]2638         {
2639           if (field->hasTimeCentered)
2640           {
2641             if (!hasTimeCentered) createCenteredAxis=true ;
2642           }
2643           if (field->hasTimeInstant)
2644           {
2645             if (!hasTimeInstant) createInstantAxis=true ;
2646           }
2647         }
[1872]2648         else if (field->getRelFile()->time_counter==CFile::time_counter_attr::record)
[1158]2649         {
2650           if (timeCounterType==none) createTimeCounterAxis=true ;
2651           timeCounterType=record ;
2652           if (field->hasTimeCentered)
2653           {
2654             if (!hasTimeCentered) createCenteredAxis=true ;
2655           }
2656           if (field->hasTimeInstant)
2657           {
2658             if (!hasTimeInstant) createInstantAxis=true ;
2659           }
2660         }
2661         
2662         if (createInstantAxis)
2663         {
2664           axisId="time_instant" ;
2665           axisBoundId="time_instant_bounds";
2666           hasTimeInstant=true ;
2667         }
[488]2668
[1158]2669         if (createCenteredAxis)
2670         {
2671           axisId="time_centered" ;
2672           axisBoundId="time_centered_bounds";
2673           hasTimeCentered=true ;
2674         }
2675
2676         
[498]2677         try
[219]2678         {
[1158]2679            std::vector<StdString> dims;
2680           
2681            if (createInstantAxis || createCenteredAxis)
2682            {
2683              // Adding time_instant or time_centered
2684              dims.push_back(timeid);
2685              if (!SuperClassWriter::varExist(axisId))
2686              {
2687                SuperClassWriter::addVariable(axisId, NC_DOUBLE, dims);
[488]2688
[1158]2689                CDate timeOrigin=cal->getTimeOrigin() ;
2690                StdOStringStream oss2;
2691                StdString strInitdate=oss2.str() ;
2692                StdString strTimeOrigin=timeOrigin.toString() ;
2693                this->writeTimeAxisAttributes(axisId, cal->getType(),strTimeUnits+strTimeOrigin,
2694                                              strTimeOrigin, axisBoundId);
2695             }
[219]2696
[1158]2697             // Adding time_instant_bounds or time_centered_bounds variables
2698             if (!SuperClassWriter::varExist(axisBoundId))
2699             {
2700                dims.clear() ;
2701                dims.push_back(timeid);
2702                dims.push_back(timeBoundId);
2703                SuperClassWriter::addVariable(axisBoundId, NC_DOUBLE, dims);
2704             }
[498]2705           }
[488]2706
[1158]2707           if (createTimeCounterAxis)
[498]2708           {
[692]2709             // Adding time_counter
[1158]2710             axisId = getTimeCounterName();
[802]2711             axisBoundId = getTimeCounterName() + "_bounds";
[692]2712             dims.clear();
2713             dims.push_back(timeid);
[1158]2714             if (!SuperClassWriter::varExist(axisId))
[692]2715             {
[1158]2716                SuperClassWriter::addVariable(axisId, NC_DOUBLE, dims);
2717                SuperClassWriter::addAttribute("axis", string("T"), &axisId);
[488]2718
[1872]2719                if (field->getRelFile()->time_counter.isEmpty() || 
2720                   (field->getRelFile()->time_counter != CFile::time_counter_attr::record))
[692]2721                {
2722                  CDate timeOrigin = cal->getTimeOrigin();
2723                  StdString strTimeOrigin = timeOrigin.toString();
[498]2724
[1158]2725                  this->writeTimeAxisAttributes(axisId, cal->getType(),
2726                                                strTimeUnits+strTimeOrigin,
[692]2727                                                strTimeOrigin, axisBoundId);
2728                }
2729             }
2730
2731             // Adding time_counter_bound dimension
[1872]2732             if (field->getRelFile()->time_counter.isEmpty() || (field->getRelFile()->time_counter != CFile::time_counter_attr::record))
[692]2733             {
2734                if (!SuperClassWriter::varExist(axisBoundId))
2735                {
2736                  dims.clear();
2737                  dims.push_back(timeid);
2738                  dims.push_back(timeBoundId);
2739                  SuperClassWriter::addVariable(axisBoundId, NC_DOUBLE, dims);
2740                }
2741             }
[498]2742           }
[488]2743         }
[498]2744         catch (CNetCdfException& e)
[488]2745         {
[498]2746           StdString msg("On writing time axis data: ");
2747           msg.append("In the context : ");
2748           CContext* context = CContext::getCurrent() ;
2749           msg.append(context->getId()); msg.append("\n");
2750           msg.append(e.what());
2751           ERROR("CNc4DataOutput::writeTimeAxis_ (CField*    field, \
[1542]2752                  const std::shared_ptr<CCalendar> cal)", << msg);
[488]2753         }
[219]2754      }
2755
2756      //---------------------------------------------------------------
[488]2757
[219]2758      void CNc4DataOutput::writeTimeAxisAttributes(const StdString & axis_name,
2759                                                   const StdString & calendar,
2760                                                   const StdString & units,
2761                                                   const StdString & time_origin,
[488]2762                                                   const StdString & time_bounds,
[219]2763                                                   const StdString & standard_name,
[613]2764                                                   const StdString & long_name)
[219]2765      {
[498]2766         try
2767         {
2768           SuperClassWriter::addAttribute("standard_name", standard_name, &axis_name);
2769           SuperClassWriter::addAttribute("long_name",     long_name    , &axis_name);
2770           SuperClassWriter::addAttribute("calendar",      calendar     , &axis_name);
2771           SuperClassWriter::addAttribute("units",         units        , &axis_name);
2772           SuperClassWriter::addAttribute("time_origin",   time_origin  , &axis_name);
2773           SuperClassWriter::addAttribute("bounds",        time_bounds  , &axis_name);
2774         }
2775         catch (CNetCdfException& e)
2776         {
2777           StdString msg("On writing time axis Attribute: ");
2778           msg.append("In the context : ");
2779           CContext* context = CContext::getCurrent() ;
2780           msg.append(context->getId()); msg.append("\n");
2781           msg.append(e.what());
2782           ERROR("CNc4DataOutput::writeTimeAxisAttributes(const StdString & axis_name, \
[613]2783                                                          const StdString & calendar,\
2784                                                          const StdString & units, \
2785                                                          const StdString & time_origin, \
2786                                                          const StdString & time_bounds, \
2787                                                          const StdString & standard_name, \
2788                                                          const StdString & long_name)", << msg);
[498]2789         }
[219]2790      }
[488]2791
[219]2792      //---------------------------------------------------------------
2793
2794      void CNc4DataOutput::writeAxisAttributes(const StdString & axis_name,
2795                                               const StdString & axis,
2796                                               const StdString & standard_name,
2797                                               const StdString & long_name,
2798                                               const StdString & units,
2799                                               const StdString & nav_model)
2800      {
[498]2801         try
2802         {
[613]2803          if (!axis.empty())
2804            SuperClassWriter::addAttribute("axis"       , axis         , &axis_name);
2805
[498]2806          SuperClassWriter::addAttribute("standard_name", standard_name, &axis_name);
2807          SuperClassWriter::addAttribute("long_name"    , long_name    , &axis_name);
2808          SuperClassWriter::addAttribute("units"        , units        , &axis_name);
[973]2809//          SuperClassWriter::addAttribute("nav_model"    , nav_model    , &axis_name);
[498]2810         }
2811         catch (CNetCdfException& e)
2812         {
2813           StdString msg("On writing Axis Attribute: ");
2814           msg.append("In the context : ");
2815           CContext* context = CContext::getCurrent() ;
2816           msg.append(context->getId()); msg.append("\n");
2817           msg.append(e.what());
2818           ERROR("CNc4DataOutput::writeAxisAttributes(const StdString & axis_name, \
[613]2819                                                      const StdString & axis, \
2820                                                      const StdString & standard_name, \
2821                                                      const StdString & long_name, \
2822                                                      const StdString & units, \
2823                                                      const StdString & nav_model)", << msg);
[498]2824         }
[219]2825      }
2826
2827      //---------------------------------------------------------------
[488]2828
[219]2829      void CNc4DataOutput::writeLocalAttributes
2830         (int ibegin, int ni, int jbegin, int nj, StdString domid)
2831      {
[498]2832        try
2833        {
[318]2834         SuperClassWriter::addAttribute(StdString("ibegin").append(domid), ibegin);
2835         SuperClassWriter::addAttribute(StdString("ni"    ).append(domid), ni);
2836         SuperClassWriter::addAttribute(StdString("jbegin").append(domid), jbegin);
2837         SuperClassWriter::addAttribute(StdString("nj"    ).append(domid), nj);
[498]2838        }
2839        catch (CNetCdfException& e)
2840        {
2841           StdString msg("On writing Local Attributes: ");
2842           msg.append("In the context : ");
2843           CContext* context = CContext::getCurrent() ;
2844           msg.append(context->getId()); msg.append("\n");
2845           msg.append(e.what());
2846           ERROR("CNc4DataOutput::writeLocalAttributes \
2847                  (int ibegin, int ni, int jbegin, int nj, StdString domid)", << msg);
2848        }
2849
[219]2850      }
2851
[628]2852      void CNc4DataOutput::writeLocalAttributes_IOIPSL(const StdString& dimXid, const StdString& dimYid,
2853                                                       int ibegin, int ni, int jbegin, int nj, int ni_glo, int nj_glo, int rank, int size)
[391]2854      {
2855         CArray<int,1> array(2) ;
2856
[498]2857         try
2858         {
2859           SuperClassWriter::addAttribute("DOMAIN_number_total",size ) ;
2860           SuperClassWriter::addAttribute("DOMAIN_number", rank) ;
[629]2861           array = SuperClassWriter::getDimension(dimXid) + 1, SuperClassWriter::getDimension(dimYid) + 1;
[498]2862           SuperClassWriter::addAttribute("DOMAIN_dimensions_ids",array) ;
2863           array=ni_glo,nj_glo ;
2864           SuperClassWriter::addAttribute("DOMAIN_size_global", array) ;
2865           array=ni,nj ;
2866           SuperClassWriter::addAttribute("DOMAIN_size_local", array) ;
[819]2867           array=ibegin+1,jbegin+1 ;
[498]2868           SuperClassWriter::addAttribute("DOMAIN_position_first", array) ;
[819]2869           array=ibegin+ni-1+1,jbegin+nj-1+1 ;
[498]2870           SuperClassWriter::addAttribute("DOMAIN_position_last",array) ;
2871           array=0,0 ;
2872           SuperClassWriter::addAttribute("DOMAIN_halo_size_start", array) ;
2873           SuperClassWriter::addAttribute("DOMAIN_halo_size_end", array);
2874           SuperClassWriter::addAttribute("DOMAIN_type",string("box")) ;
2875  /*
2876           SuperClassWriter::addAttribute("DOMAIN_DIM_N001",string("x")) ;
2877           SuperClassWriter::addAttribute("DOMAIN_DIM_N002",string("y")) ;
2878           SuperClassWriter::addAttribute("DOMAIN_DIM_N003",string("axis_A")) ;
2879           SuperClassWriter::addAttribute("DOMAIN_DIM_N004",string("time_counter")) ;
2880  */
2881         }
2882         catch (CNetCdfException& e)
2883         {
[628]2884           StdString msg("On writing Local Attributes IOIPSL \n");
[498]2885           msg.append("In the context : ");
2886           CContext* context = CContext::getCurrent() ;
2887           msg.append(context->getId()); msg.append("\n");
2888           msg.append(e.what());
2889           ERROR("CNc4DataOutput::writeLocalAttributes_IOIPSL \
2890                  (int ibegin, int ni, int jbegin, int nj, int ni_glo, int nj_glo, int rank, int size)", << msg);
2891         }
[391]2892      }
[219]2893      //---------------------------------------------------------------
2894
2895      void CNc4DataOutput:: writeFileAttributes(const StdString & name,
2896                                                const StdString & description,
2897                                                const StdString & conventions,
2898                                                const StdString & production,
2899                                                const StdString & timeStamp)
2900      {
[498]2901         try
2902         {
2903           SuperClassWriter::addAttribute("name"       , name);
2904           SuperClassWriter::addAttribute("description", description);
[613]2905           SuperClassWriter::addAttribute("title"      , description);
2906           SuperClassWriter::addAttribute("Conventions", conventions);
[1158]2907           // SuperClassWriter::addAttribute("production" , production);
2908
2909           StdString timeStampStr ;
2910           if (file->time_stamp_name.isEmpty()) timeStampStr="timeStamp" ;
2911           else timeStampStr=file->time_stamp_name ;
2912           SuperClassWriter::addAttribute(timeStampStr, timeStamp);
2913
2914           StdString uuidName ;
2915           if (file->uuid_name.isEmpty()) uuidName="uuid" ;
2916           else uuidName=file->uuid_name ;
2917
2918           if (file->uuid_format.isEmpty()) SuperClassWriter::addAttribute(uuidName, getUuidStr());
2919           else SuperClassWriter::addAttribute(uuidName, getUuidStr(file->uuid_format));
2920         
[498]2921         }
2922         catch (CNetCdfException& e)
2923         {
2924           StdString msg("On writing File Attributes \n ");
2925           msg.append("In the context : ");
2926           CContext* context = CContext::getCurrent() ;
2927           msg.append(context->getId()); msg.append("\n");
2928           msg.append(e.what());
2929           ERROR("CNc4DataOutput:: writeFileAttributes(const StdString & name, \
2930                                                const StdString & description, \
2931                                                const StdString & conventions, \
2932                                                const StdString & production, \
2933                                                const StdString & timeStamp)", << msg);
2934         }
[219]2935      }
2936
2937      //---------------------------------------------------------------
2938
2939      void CNc4DataOutput::writeMaskAttributes(const StdString & mask_name,
2940                                               int data_dim,
2941                                               int data_ni,
2942                                               int data_nj,
2943                                               int data_ibegin,
2944                                               int data_jbegin)
2945      {
[498]2946         try
2947         {
2948           SuperClassWriter::addAttribute("data_dim"   , data_dim   , &mask_name);
2949           SuperClassWriter::addAttribute("data_ni"    , data_ni    , &mask_name);
2950           SuperClassWriter::addAttribute("data_nj"    , data_nj    , &mask_name);
2951           SuperClassWriter::addAttribute("data_ibegin", data_ibegin, &mask_name);
2952           SuperClassWriter::addAttribute("data_jbegin", data_jbegin, &mask_name);
2953         }
2954         catch (CNetCdfException& e)
2955         {
2956           StdString msg("On writing Mask Attributes \n ");
2957           msg.append("In the context : ");
2958           CContext* context = CContext::getCurrent() ;
2959           msg.append(context->getId()); msg.append("\n");
2960           msg.append(e.what());
2961           ERROR("CNc4DataOutput::writeMaskAttributes(const StdString & mask_name, \
2962                                               int data_dim, \
2963                                               int data_ni, \
2964                                               int data_nj, \
2965                                               int data_ibegin, \
2966                                               int data_jbegin)", << msg);
2967         }
[219]2968      }
2969
2970      ///--------------------------------------------------------------
2971
[1158]2972      StdSize CNc4DataOutput::getRecordFromTime(Time time, double factorUnit)
[707]2973      {
2974        std::map<Time, StdSize>::const_iterator it = timeToRecordCache.find(time);
2975        if (it == timeToRecordCache.end())
2976        {
[802]2977          StdString timeAxisBoundsId(getTimeCounterName() + "_bounds");
[1158]2978          if (!SuperClassWriter::varExist(timeAxisBoundsId)) timeAxisBoundsId = "time_centered_bounds";
2979          if (!SuperClassWriter::varExist(timeAxisBoundsId)) timeAxisBoundsId = "time_instant_bounds";
[707]2980
2981          CArray<double,2> timeAxisBounds;
[1158]2982          std::vector<StdSize> dimSize(SuperClassWriter::getDimensions(timeAxisBoundsId)) ;
2983         
[707]2984          StdSize record = 0;
2985          double dtime(time);
[1158]2986          for (int n = dimSize[0] - 1; n >= 0; n--)
[707]2987          {
[1158]2988            SuperClassWriter::getTimeAxisBounds(timeAxisBounds, timeAxisBoundsId, isCollective, n);
2989            timeAxisBounds*=factorUnit ;
2990            if (timeAxisBounds(1, 0) < dtime)
[707]2991            {
2992              record = n + 1;
2993              break;
2994            }
2995          }
2996          it = timeToRecordCache.insert(std::make_pair(time, record)).first;
2997        }
2998        return it->second;
2999      }
[774]3000
3001      ///--------------------------------------------------------------
3002
3003      bool CNc4DataOutput::isWrittenDomain(const std::string& domainName) const
3004      {
3005        return (this->writtenDomains.find(domainName) != this->writtenDomains.end());
3006      }
3007
3008      bool CNc4DataOutput::isWrittenCompressedDomain(const std::string& domainName) const
3009      {
3010        return (this->writtenCompressedDomains.find(domainName) != this->writtenCompressedDomains.end());
3011      }
3012
3013      bool CNc4DataOutput::isWrittenAxis(const std::string& axisName) const
3014      {
3015        return (this->writtenAxis.find(axisName) != this->writtenAxis.end());
3016      }
3017
3018      bool CNc4DataOutput::isWrittenCompressedAxis(const std::string& axisName) const
3019      {
3020        return (this->writtenCompressedAxis.find(axisName) != this->writtenCompressedAxis.end());
3021      }
3022
[887]3023      bool CNc4DataOutput::isWrittenScalar(const std::string& scalarName) const
3024      {
3025        return (this->writtenScalar.find(scalarName) != this->writtenScalar.end());
3026      }
3027
[774]3028      void CNc4DataOutput::setWrittenDomain(const std::string& domainName)
3029      {
3030        this->writtenDomains.insert(domainName);
3031      }
3032
3033      void CNc4DataOutput::setWrittenCompressedDomain(const std::string& domainName)
3034      {
3035        this->writtenCompressedDomains.insert(domainName);
3036      }
3037
3038      void CNc4DataOutput::setWrittenAxis(const std::string& axisName)
3039      {
3040        this->writtenAxis.insert(axisName);
3041      }
3042
3043      void CNc4DataOutput::setWrittenCompressedAxis(const std::string& axisName)
3044      {
3045        this->writtenCompressedAxis.insert(axisName);
3046      }
[887]3047
3048      void CNc4DataOutput::setWrittenScalar(const std::string& scalarName)
3049      {
3050        this->writtenScalar.insert(scalarName);
3051      }
[335]3052} // namespace xios
Note: See TracBrowser for help on using the repository browser.