Added operator grid2fourier.

src/Collgrid.cc

@@ -29,7 +29,7 @@
 
 static int globalGridType = CDI_UNDEFID;
 
-struct ensfileType
+struct EnsfileType
 {
   CdoStreamID streamID;
   int vlistID;
@@ -49,8 +49,8 @@ struct xyinfoType
 static int
 cmpx(const void *a, const void *b)
 {
-  const double x = ((const xyinfoType *) a)->x;
-  const double y = ((const xyinfoType *) b)->x;
+  const auto x = ((const xyinfoType *) a)->x;
+  const auto y = ((const xyinfoType *) b)->x;
   return ((x > y) - (x < y)) * 2 + (x > y) - (x < y);
 }
 
@@ -89,7 +89,7 @@ cmpxy_gt(const void *a, const void *b)
 }
 
 static int
-genGrid(int ngrids, int nfiles, std::vector<ensfileType> &ef, bool ginit, int igrid, long nxblocks)
+genGrid(int ngrids, int nfiles, std::vector<EnsfileType> &ef, bool ginit, int igrid, long nxblocks)
 {
   bool lsouthnorth = true;
   bool lregular = false;
@@ -109,8 +109,7 @@ genGrid(int ngrids, int nfiles, std::vector<ensfileType> &ef, bool ginit, int ig
   const bool lbounds = lunstructured && globalGridType == CDI_UNDEFID && gridInqXbounds(gridID, nullptr) > 0
                        && gridInqYbounds(gridID, nullptr) > 0;
 
-  std::vector<long> xsize(nfiles);
-  std::vector<long> ysize(nfiles);
+  std::vector<long> xsize(nfiles), ysize(nfiles);
   xyinfoType *xyinfo = (xyinfoType *) Malloc(nfiles * sizeof(xyinfoType));
   Varray2D<double> xvals(nfiles), yvals(nfiles);
   Varray2D<double> xbounds(nfiles), ybounds(nfiles);
@@ -264,8 +263,7 @@ genGrid(int ngrids, int nfiles, std::vector<ensfileType> &ef, bool ginit, int ig
         }
     }
 
-  std::vector<long> xoff(nx + 1);
-  std::vector<long> yoff(ny + 1);
+  std::vector<long> xoff(nx + 1), yoff(ny + 1);
 
   xoff[0] = 0;
   for (long i = 0; i < nx; ++i)
@@ -370,7 +368,7 @@ Collgrid(void *process)
   if (!Options::cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename))
     cdoAbort("Outputfile %s already exists!", ofilename);
 
-  std::vector<ensfileType> ef(nfiles);
+  std::vector<EnsfileType> ef(nfiles);
 
   for (int fileID = 0; fileID < nfiles; fileID++)
     {
@@ -453,10 +451,7 @@ Collgrid(void *process)
               cdoWarning("Variable name %s not found!", cdoOperatorArgv(noff + isel).c_str());
             }
         }
-      if (err)
-        {
-          cdoAbort("Could not find all requested variables: (%d/%d)", nsel - err, nsel);
-        }
+      if (err) cdoAbort("Could not find all requested variables: (%d/%d)", nsel - err, nsel);
     }
 
   for (varID = 0; varID < nvars; varID++)

src/FC.cc

@@ -24,6 +24,16 @@
       FC         gp2fc           Gridpoint to fourier
 */
 
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_LIBFFTW3
+#include <fftw3.h>
+#endif
+
+#include <mutex>
+
 #include <cdi.h>
 
 #include "cdo_vlist.h"
@@ -31,6 +41,87 @@
 #include <mpim_grid.h>
 #include "specspace.h"
 #include "griddes.h"
+#include "cimdOmp.h"
+
+static std::mutex fftw_mutex;
+
+
+static int
+vlistGetFirstReg2DGrid(int vlistID)
+{
+  // find first gaussian grid
+  const auto ngrids = vlistNgrids(vlistID);
+  for (int index = 0; index < ngrids; index++)
+    {
+      const auto gridID = vlistGrid(vlistID, index);
+      if (gridInqType(gridID) == GRID_GAUSSIAN ||
+          gridInqType(gridID) == GRID_LONLAT ||
+          gridInqType(gridID) == GRID_CURVILINEAR) return gridID;
+    }
+
+  return -1;
+}
+
+static void
+grid2fourier(int gridID1, const double *array1, int gridID2, double *array2)
+{
+  auto gridsize = gridInqSize(gridID1);
+  auto nlon = gridInqXsize(gridID1);
+  auto nlat = gridInqYsize(gridID1);
+
+  for (size_t i = 0; i < gridsize; ++i) array2[i] = array1[i];
+
+#ifdef HAVE_LIBFFTW3
+  double norm = 1. / nlon;
+  struct FourierMemory
+  {
+    double *in_fft;
+    fftw_complex *out_fft;
+    fftw_plan plan;
+  };
+
+  std::vector<FourierMemory> ompmem(Threading::ompNumThreads);
+
+  for (int i = 0; i < Threading::ompNumThreads; i++)
+    {
+      ompmem[i].in_fft = (double *) fftw_malloc(nlon * sizeof(double));
+      ompmem[i].out_fft = fftw_alloc_complex(nlon);
+      std::unique_lock<std::mutex> locked_mutex(fftw_mutex);
+      ompmem[i].plan = fftw_plan_dft_r2c_1d(nlon, ompmem[i].in_fft, ompmem[i].out_fft, FFTW_ESTIMATE);
+    }
+
+  if (Options::cdoVerbose) fftw_print_plan(ompmem[0].plan);
+
+#ifdef _OPENMP
+#pragma omp parallel for default(shared)
+#endif
+  for (size_t ilat = 0; ilat < nlat; ++ilat)
+    {
+      const auto ompthID = cdo_omp_get_thread_num();
+      auto in_fft = ompmem[ompthID].in_fft;
+      auto out_fft = ompmem[ompthID].out_fft;
+
+      for (size_t ilon = 0; ilon < nlon; ++ilon) in_fft[ilon] = array1[ilat * nlon + ilon];
+
+      fftw_execute(ompmem[ompthID].plan);
+
+      for (size_t ifc = 0; ifc < nlon; ++ifc)
+        {
+          array2[2 * (ilat * nlon + ifc)] = norm * out_fft[ifc][0];
+          array2[2 * (ilat * nlon + ifc) + 1] = norm * out_fft[ifc][1];
+        }
+    }
+
+  for (int i = 0; i < Threading::ompNumThreads; i++)
+    {
+      fftw_free(ompmem[i].in_fft);
+      fftw_free(ompmem[i].out_fft);
+      fftw_destroy_plan(ompmem[i].plan);
+    }
+#else
+  cdoAbort("FFTW support not compiled in!");
+#endif
+}
 
 void *
 FC(void *process)
@@ -53,8 +144,10 @@ FC(void *process)
   const auto SP2FC = cdoOperatorAdd("sp2fc", 0, 0, nullptr);
   const auto FC2GP = cdoOperatorAdd("fc2gp", 0, 0, nullptr);
   const auto GP2FC = cdoOperatorAdd("gp2fc", 0, 0, nullptr);
+  const auto GRID2FOURIER = cdoOperatorAdd("grid2fourier", 1, 0, nullptr);
 
   const auto operatorID = cdoOperatorID();
+  const auto operfunc = cdoOperatorF1(operatorID);
 
   const auto streamID1 = cdoOpenRead(0);
 
@@ -65,9 +158,13 @@ FC(void *process)
   const auto taxisID2 = taxisDuplicate(taxisID1);
   vlistDefTaxis(vlistID2, taxisID2);
 
-  auto gridIDsp = vlistGetFirstSpectralGrid(vlistID1);
-  auto gridIDgp = vlistGetFirstGaussianGrid(vlistID1);
-  auto gridIDfc = vlistGetFirstFourierGrid(vlistID1);
+  int gridIDsp = -1, gridIDgp = -1, gridIDfc = -1;
+  if (operfunc == 0)
+    {
+      gridIDsp = vlistGetFirstSpectralGrid(vlistID1);
+      gridIDgp = vlistGetFirstGaussianGrid(vlistID1);
+      gridIDfc = vlistGetFirstFourierGrid(vlistID1);
+    }
 
   // define output grid
   if (operatorID == FC2SP)
@@ -180,15 +277,28 @@ FC(void *process)
           nlat = gridInqYsize(gridID2);
         }
     }
+  else if (operatorID == GRID2FOURIER)
+    {
+      gridID1 = vlistGetFirstReg2DGrid(vlistID1);
+      if (gridID1 != -1)
+        {
+          nlon = gridInqXsize(gridID1);
+          nlat = gridInqYsize(gridID1);
+          gridID2 = gridID1;
+        }
+    }
 
   FC_Transformation fcTrans;
   SP_Transformation spTrans;
-  if (nlon > 0)
+  if (operfunc == 0)
     {
-      if (operatorID == GP2FC || operatorID == FC2GP)
-        fcTrans.init(nlon, nlat, ntr);
-      else
-        spTrans.init(nlon, nlat, ntr, 0);
+      if (nlon > 0)
+        {
+          if (operatorID == GP2FC || operatorID == FC2GP)
+            fcTrans.init(nlon, nlat, ntr);
+          else
+            spTrans.init(nlon, nlat, ntr, 0);
+        }
     }
 
   // printf("nfc %d, ntr %d, nlat %zu, nlon %zu\n", nfc, ntr, nlat, nlon);
@@ -198,6 +308,12 @@ FC(void *process)
   for (varID = 0; varID < nvars; varID++) vars[varID] = gridID1 == vlistInqVarGrid(vlistID1, varID);
 
   if (gridID1 != -1) vlistChangeGrid(vlistID2, gridID1, gridID2);
+  if (operfunc == 1)
+    {
+      for (varID = 0; varID < nvars; varID++)
+        if (vars[varID])
+          vlistDefVarDatatype(vlistID2, varID, CDI_DATATYPE_CPX32);
+    }
 
   const auto streamID2 = cdoOpenWrite(1);
 
@@ -208,7 +324,7 @@ FC(void *process)
 
   if (gridID2 != -1)
     {
-      const auto gridsize = gridInqSize(gridID2);
+      const auto gridsize = 2 * gridInqSize(gridID2);
       array2.resize(gridsize);
     }
 
@@ -236,6 +352,8 @@ FC(void *process)
                 four2grid(fcTrans, gridID1, array1.data(), gridID2, array2.data());
               else if (operatorID == GP2FC)
                 grid2four(fcTrans, gridID1, array1.data(), gridID2, array2.data());
+              else if (operatorID == GRID2FOURIER)
+                grid2fourier(gridID1, array1.data(), gridID2, array2.data());
 
               cdoDefRecord(streamID2, varID, levelID);
               cdoWriteRecord(streamID2, array2.data(), nmiss);

src/Fourier.cc

@@ -184,8 +184,8 @@ Fourier(void *process)
 
   const auto nvars = vlistNvars(vlistID1);
   FieldVector3D vars;
-  std::vector<int64_t> vdate;
-  std::vector<int> vtime;
+  std::vector<int64_t> vdates;
+  std::vector<int> vtimes;
 
   int tsID = 0;
   while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
@@ -193,13 +193,13 @@ Fourier(void *process)
       if (tsID >= nalloc)
         {
           nalloc += NALLOC_INC;
-          vdate.resize(nalloc);
-          vtime.resize(nalloc);
+          vdates.resize(nalloc);
+          vtimes.resize(nalloc);
           vars.resize(nalloc);
         }
 
-      vdate[tsID] = taxisInqVdate(taxisID1);
-      vtime[tsID] = taxisInqVtime(taxisID1);
+      vdates[tsID] = taxisInqVdate(taxisID1);
+      vtimes[tsID] = taxisInqVtime(taxisID1);
 
       fieldsFromVlist(vlistID1, vars[tsID]);
 
@@ -272,8 +272,8 @@ Fourier(void *process)
 
   for (tsID = 0; tsID < nts; tsID++)
     {
-      taxisDefVdate(taxisID2, vdate[tsID]);
-      taxisDefVtime(taxisID2, vtime[tsID]);
+      taxisDefVdate(taxisID2, vdates[tsID]);
+      taxisDefVtime(taxisID2, vtimes[tsID]);
       cdoDefTimestep(streamID2, tsID);
 
       for (varID = 0; varID < nvars; varID++)

src/module_definitions.h

@@ -156,7 +156,7 @@ void   *Expr(void *argument);
 #define ExprOperators                     {"expr", "exprf", "aexpr", "aexprf"}
 
 void   *FC(void *argument);
-#define FCOperators                       {"fc2sp", "sp2fc", "fc2gp", "gp2fc"}
+#define FCOperators                       {"fc2sp", "sp2fc", "fc2gp", "gp2fc", "grid2fourier"}
 
 void   *Filedes(void *argument);
 #define FiledesOperators                  {"filedes", "griddes", "griddes2", "zaxisdes", "vct", "vct2", "codetab", "vlist", "partab", "partab2", "spartab"}