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