File: graceL1b2Mass.cpp

package info (click to toggle)
groops 0%2Bgit20250907%2Bds-1
  • links: PTS, VCS
  • area: non-free
  • in suites: forky, sid
  • size: 11,140 kB
  • sloc: cpp: 135,607; fortran: 1,603; makefile: 20
file content (135 lines) | stat: -rw-r--r-- 4,400 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/***********************************************/
/**
* @file graceL1b2Mass.cpp
*
* @brief Read GRACE L1B data.
*
* @author Beate Klinger
* @date 2013-11-24
*
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
This program converts mass data (MAS1B or MAS1A) from the GRACE SDS format into \file{instrument file (MASS)}{instrument}.
For further information see \program{GraceL1b2Accelerometer}.
)";

/***********************************************/

#include "programs/program.h"
#include "files/fileInstrument.h"
#include "fileGrace.h"

/***** CLASS ***********************************/

/** @brief Read GRACE L1B data.
* @ingroup programsConversionGroup */
class GraceL1b2Mass
{
public:
  void run(Config &config, Parallel::CommunicatorPtr comm);
};

GROOPS_REGISTER_PROGRAM(GraceL1b2Mass, SINGLEPROCESS, "read GRACE L1B data (MAS1B or MAS1A)", Conversion, Grace, Instrument)

/***********************************************/

void GraceL1b2Mass::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
  try
  {
    FileName              fileNameOut;
    std::vector<FileName> fileNameIn;

    readConfig(config, "outputfileMass", fileNameOut, Config::MUSTSET,  "", "MASS");
    readConfig(config, "inputfile",      fileNameIn,  Config::MUSTSET,  "", "MAS1B or MAS1A");
    if(isCreateSchema(config)) return;

    // =============================================

    logStatus<<"read input files"<<Log::endl;
    Arc arc;
    for(UInt idFile=0; idFile<fileNameIn.size(); idFile++)
    {
      logStatus<<"read file <"<<fileNameIn.at(idFile)<<">"<<Log::endl;
      UInt numberOfRecords;
      FileInGrace file(fileNameIn.at(idFile), numberOfRecords);

      for(UInt idEpoch=0; idEpoch<numberOfRecords; idEpoch++)
      {
        Int32    seconds, time_frac;
        Char     time_ref, GRACE_id;
        Byte     qualflg, prod_flag;
        Double   mass_thr=0, mass_thr_err=0, mass_tnk=0, mass_tnk_err=0;
        Double   gas_mass_thr1=0, gas_mass_thr2=0, gas_mass_tnk1=0, gas_mass_tnk2=0;

        try
        {
          file>>seconds>>time_frac>>time_ref>>GRACE_id>>FileInGrace::flag(qualflg)>>FileInGrace::flag(prod_flag);
          if(prod_flag & (1<<0)) file>>mass_thr;
          if(prod_flag & (1<<1)) file>>mass_thr_err;
          if(prod_flag & (1<<2)) file>>mass_tnk;
          if(prod_flag & (1<<3)) file>>mass_tnk_err;
          if(prod_flag & (1<<4)) file>>gas_mass_thr1;
          if(prod_flag & (1<<5)) file>>gas_mass_thr2;
          if(prod_flag & (1<<6)) file>>gas_mass_tnk1;
          if(prod_flag & (1<<7)) file>>gas_mass_tnk2;
        }
        catch(std::exception &/*e*/)
        {
          // GRACE-FO number of records issue
          logWarning<<arc.back().time.dateTimeStr()<<": file ended at "<<idEpoch<<" of "<<numberOfRecords<<" expected records"<<Log::endl;
          break;
        }

        const Time time = mjd2time(51544.5) + seconds2time(seconds) + seconds2time(1e-6*time_frac);
        if(arc.size() && (time <= arc.back().time))
          logWarning<<"epoch("<<time.dateTimeStr()<<") <= last epoch("<<arc.back().time.dateTimeStr()<<")"<<Log::endl;

        MassEpoch epoch;
        epoch.time     = time;
        if(GRACE_id=='C')
        {
          epoch.massThr = 601.21-(15.62718541166349-gas_mass_tnk1)-(15.63189368778949-gas_mass_tnk2);
        }
        else if(GRACE_id=='D')
        {
          epoch.massThr = 601.21-(15.60994853989272-gas_mass_tnk1)- (15.68091099702227-gas_mass_tnk2);
        }
        else
            epoch.massThr  = mass_thr;
        epoch.massTank = epoch.massThr;
        arc.push_back(epoch);
      } // for(idEpoch)
    } // for(idFile)

    // =============================================

    logStatus<<"sort epochs"<<Log::endl;
    arc.sort();

    logStatus<<"eliminate duplicates"<<Log::endl;
    const UInt oldSize = arc.size();
    arc.removeDuplicateEpochs(TRUE/*keepFirst*/);
    if(arc.size() < oldSize)
      logInfo<<" "<<oldSize-arc.size()<<" duplicates removed!"<<Log::endl;

    Arc::printStatistics(arc);
    if(arc.size() == 0)
      return;

    if(!fileNameOut.empty())
    {
      logInfo<<"write data to <"<<fileNameOut<<">"<<Log::endl;
      InstrumentFile::write(fileNameOut, arc);
    }
  }
  catch(std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}

/***********************************************/