Added operator fourier2grid.

src/FC.cc

@@ -63,13 +63,67 @@ vlistGetFirstReg2DGrid(int vlistID)
 }
 
 static void
-grid2fourier(int gridID1, const double *array1, int gridID2, double *array2)
+fourier2grid(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
+  struct FourierMemory
+  {
+    fftw_complex *in_fft;
+    double *out_fft;
+    fftw_plan plan;
+  };
+
+  std::vector<FourierMemory> ompmem(Threading::ompNumThreads);
+
+  for (int i = 0; i < Threading::ompNumThreads; i++)
+    {
+      ompmem[i].in_fft = fftw_alloc_complex(nlon);
+      ompmem[i].out_fft = (double *) fftw_malloc(nlon * sizeof(double));
+      std::unique_lock<std::mutex> locked_mutex(fftw_mutex);
+      ompmem[i].plan = fftw_plan_dft_c2r_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 ifc = 0; ifc < nlon; ++ifc)
+        {
+          in_fft[ifc][0] = array1[2 * (ilat * nlon + ifc)];
+          in_fft[ifc][1] = array1[2 * (ilat * nlon + ifc) + 1];
+        }
+
+      fftw_execute(ompmem[ompthID].plan);
+
+      for (size_t ilon = 0; ilon < nlon; ++ilon) array2[ilat * nlon + ilon] = out_fft[ilon];
+    }
+
+  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
+}
+
+static void
+grid2fourier(int gridID1, const double *array1, int gridID2, double *array2)
+{
+  auto nlon = gridInqXsize(gridID1);
+  auto nlat = gridInqYsize(gridID1);
 
 #ifdef HAVE_LIBFFTW3
   double norm = 1. / nlon;
@@ -145,6 +199,7 @@ FC(void *process)
   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 FOURIER2GRID = cdoOperatorAdd("fourier2grid", 1, 0, nullptr);
 
   const auto operatorID = cdoOperatorID();
   const auto operfunc = cdoOperatorF1(operatorID);
@@ -277,9 +332,22 @@ FC(void *process)
           nlat = gridInqYsize(gridID2);
         }
     }
+  else if (operatorID == FOURIER2GRID)
+    {
+      gridID1 = vlistGetFirstReg2DGrid(vlistID1);
+      if (gridID1 == -1) cdoWarning("No regular 2D data found!");
+      if (gridID1 != -1)
+        {
+          nlon = gridInqXsize(gridID1);
+          nlat = gridInqYsize(gridID1);
+          gridID2 = gridID1;
+        }
+    }
   else if (operatorID == GRID2FOURIER)
     {
       gridID1 = vlistGetFirstReg2DGrid(vlistID1);
+      if (gridID1 == -1) cdoWarning("No regular 2D data found!");
+
       if (gridID1 != -1)
         {
           nlon = gridInqXsize(gridID1);
@@ -308,23 +376,31 @@ 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)
+  if (operatorID == GRID2FOURIER)
     {
       for (varID = 0; varID < nvars; varID++)
         if (vars[varID])
           vlistDefVarDatatype(vlistID2, varID, CDI_DATATYPE_CPX32);
     }
+  else if (operatorID == FOURIER2GRID)
+    {
+      for (varID = 0; varID < nvars; varID++)
+        if (vars[varID])
+          vlistDefVarDatatype(vlistID2, varID, CDI_DATATYPE_FLT32);
+    }
 
   const auto streamID2 = cdoOpenWrite(1);
 
   cdoDefVlist(streamID2, vlistID2);
 
-  const auto gridsizemax = vlistGridsizeMax(vlistID1);
+  auto gridsizemax = vlistGridsizeMax(vlistID1);
+  if (operatorID == FOURIER2GRID) gridsizemax *= 2;
   Varray<double> array1(gridsizemax), array2;
 
   if (gridID2 != -1)
     {
-      const auto gridsize = 2 * gridInqSize(gridID2);
+      auto gridsize = gridInqSize(gridID2);
+      if (operatorID == GRID2FOURIER) gridsize *= 2;
       array2.resize(gridsize);
     }
 
@@ -352,6 +428,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 == FOURIER2GRID)
+                fourier2grid(gridID1, array1.data(), gridID2, array2.data());
               else if (operatorID == GRID2FOURIER)
                 grid2fourier(gridID1, array1.data(), gridID2, array2.data());
 

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", "grid2fourier"}
+#define FCOperators                       {"fc2sp", "sp2fc", "fc2gp", "gp2fc", "fourier2grid", "grid2fourier"}
 
 void   *Filedes(void *argument);
 #define FiledesOperators                  {"filedes", "griddes", "griddes2", "zaxisdes", "vct", "vct2", "codetab", "vlist", "partab", "partab2", "spartab"}

src/module_list.h

@@ -60,9 +60,9 @@ static const module_t module_Eofcoeff       = {Eofcoeff       , EofcoeffHelp
 static const module_t module_Eofcoeff3d     = {Eofcoeff3d     , EofcoeffHelp      , Eofcoeff3dOperators     , EXPOSED  , CDI_REAL , 2  , OBASE , NoRestriction };
 static const module_t module_EOFs           = {EOFs           , EOFsHelp          , EOFsOperators           , EXPOSED  , CDI_REAL , 1  , 2     , OnlyFirst     };
 static const module_t module_EOF3d          = {EOF3d          , EOFsHelp          , EOF3dOperators          , EXPOSED  , CDI_REAL , 1  , 2     , OnlyFirst     };
-static const module_t module_EstFreq        = {EstFreq        , {}                , EstFreqOperators        , INTERNAL  , CDI_REAL , 1  , 1     , NoRestriction };
+static const module_t module_EstFreq        = {EstFreq        , {}                , EstFreqOperators        , INTERNAL , CDI_REAL , 1  , 1     , NoRestriction };
 static const module_t module_Expr           = {Expr           , ExprHelp          , ExprOperators           , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
-static const module_t module_FC             = {FC             , {}                , FCOperators             , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
+static const module_t module_FC             = {FC             , {}                , FCOperators             , EXPOSED  , CDI_BOTH , 1  , 1     , NoRestriction };
 static const module_t module_Filedes        = {Filedes        , FiledesHelp       , FiledesOperators        , EXPOSED  , CDI_BOTH , 1  , 0     , NoRestriction };
 static const module_t module_Fillmiss       = {Fillmiss       , {}                , FillmissOperators       , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };
 static const module_t module_Filter         = {Filter         , FilterHelp        , FilterOperators         , EXPOSED  , CDI_REAL , 1  , 1     , NoRestriction };