1 /*
2   This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3 
4   Author: Uwe Schulzweida
5 
6 */
7 #ifndef SPECSPACE_H
8 #define SPECSPACE_H
9 
10 #ifdef HAVE_CONFIG_H
11 #include "config.h"
12 #endif
13 
14 #include "cdo_options.h"
15 #include "transform.h"
16 #include "varray.h"
17 
18 enum class PolFlag
19 {
20   UNDEF,
21   SP2FC,
22   FC2SP,
23   UV2DV,
24 };
25 
26 class // FC_Transformation
27 #ifdef WARN_UNUSED
28 [[gnu::warn_unused]]
29 #endif
30 FC_Transformation
31 {
32 public:
33   bool use_fftw = false;
34   long nlon = 0;
35   long nlat = 0;
36   long ntr = 0;
37   long nlev = 0;
38   long ifax[10] = {};
39   Varray<double> vtrig;
40 
FC_Transformation()41   FC_Transformation()
42   {
43     if (Options::Use_FFTW)
44       {
45 #ifdef HAVE_LIBFFTW3
46         if (Options::cdoVerbose) cdo_print("Using fftw3 lib");
47         use_fftw = true;
48 #else
49         if (Options::cdoVerbose) cdo_print("LIBFFTW3 support not compiled in!");
50 #endif
51       }
52   }
53 
54   void
55   init(long _nlon, long _nlat, long _ntr, long _nlev = 0)
56   {
57     if (_nlon <= 0 || _nlat <= 0 || _ntr <= 0)
58       {
59         fprintf(stderr, "SP_Transformation.init(): parameter not initialized\n");
60         return;
61       }
62 
63     nlon = _nlon;
64     nlat = _nlat;
65     ntr = _ntr;
66     nlev = _nlev;
67 
68     //if (nlev > 1) use_fftw = false;
69 
70     if (use_fftw == false)
71       {
72         vtrig.resize(nlon);
73         const auto status = fft_set(vtrig.data(), ifax, nlon);
74         if (status < 0) cdo_abort("FFT error!");
75       }
76   }
77 };
78 
79 class // SP_Transformation
80 #ifdef WARN_UNUSED
81 [[gnu::warn_unused]]
82 #endif
83 SP_Transformation
84 {
85 public:
86   FC_Transformation fcTrans;
87   long nlat = 0;
88   long ntr = 0;
89   Varray<double> poli;
90   Varray<double> pold;
91   Varray<double> pol2;     // only for uv2dv
92   Varray<double> pol3;     // only for uv2dv
93   Varray<double> coslat;   // only for scaluv with uv2dv
94   Varray<double> rcoslat;  // only for scaluv with dv2uv
95 
SP_Transformation()96   SP_Transformation() {}
97 
98   void
99   init(long _nlon, long _nlat, long _ntr, PolFlag polFlag, long _nlev = 0)
100   {
101     if (_nlon <= 0 || _nlat <= 0 || _ntr <= 0)
102       {
103         fprintf(stderr, "SP_Transformation.init(): parameter not initialized\n");
104         return;
105       }
106 
107     fcTrans.init(_nlon, _nlat, _ntr, _nlev);
108 
109     nlat = _nlat;
110     ntr = _ntr;
111 
112     const long nsp = (ntr + 1) * (ntr + 2);
113     const long poldim = (nsp / 2) * nlat;
114 
115     if (polFlag == PolFlag::SP2FC) varrayResize(poli, poldim);
116     if (polFlag == PolFlag::FC2SP) varrayResize(pold, poldim);
117 
118     if (polFlag == PolFlag::UV2DV) varrayResize(pol2, poldim);
119     if (polFlag == PolFlag::UV2DV) varrayResize(pol3, poldim);
120 
121     coslat.resize(nlat);
122     rcoslat.resize(nlat);
123 
124     after_legini_full(ntr, nlat, poli.data(), pold.data(), nullptr, pol2.data(), pol3.data(), coslat.data());
125 
126     for (long jgl = 0; jgl < nlat; ++jgl) rcoslat[jgl] = 1.0 / coslat[jgl];
127   }
128 };
129 
130 class // DV_Transformation
131 #ifdef WARN_UNUSED
132 [[gnu::warn_unused]]
133 #endif
134 DV_Transformation
135 {
136 public:
137   long ntr = 0;
138   long fdim = 0;
139   Varray<double> f1;
140   Varray<double> f2;
141 
DV_Transformation()142   DV_Transformation() {}
143 
144   void
init(long _ntr)145   init(long _ntr)
146   {
147     ntr = _ntr;
148 
149     const long dimsp = (ntr + 1) * (ntr + 2);
150     fdim = dimsp / 2;
151 
152     f1.resize(fdim);
153     f2.resize(fdim);
154 
155     geninx(ntr, f1.data(), f2.data());
156   }
157 };
158 
159 void dv2ps(const double *div, double *pot, long nlev, long ntr);
160 
161 void trans_uv2dv(const SP_Transformation &spTrans, long nlev, int gridID1, double *gu, double *gv, int gridID2, double *sd,
162                  double *svo);
163 
164 void trans_dv2uv(const SP_Transformation &spTrans, const DV_Transformation &dvTrans, long nlev, int gridID1, double *sd,
165                  double *svo, int gridID2, double *gu, double *gv);
166 
167 void grid2spec(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
168 void spec2grid(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
169 void four2spec(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
170 void spec2four(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
171 void four2grid(const FC_Transformation &fcTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
172 void grid2four(const FC_Transformation &fcTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
173 
174 void spec2spec(int gridIDin, double *arrayIn, int gridIDout, double *arrayOut);
175 void speccut(int gridIDin, double *arrayIn, double *arrayOut, int *waves);
176 
177 void spcut(double *arrayIn, double *arrayOut, long ntr, const int *waves);
178 
179 #endif
180