FV3 Bundle
fft_f.cc
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2017-2018 UCAR
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  */
7 
8 #include "model/fft_f.h"
9 
10 #include <unsupported/Eigen/FFT>
11 #include <cmath>
12 #include <vector>
13 
14 #include "oops/util/Logger.h"
15 
16 // -----------------------------------------------------------------------------
17 namespace qg {
18 // -----------------------------------------------------------------------------
19 
20 void fft_fwd_f(const int & kk, const double * xx, double * ff) {
21  static Eigen::FFT<double> fft;
22  std::vector<double> grid(kk);
23  const unsigned int size = kk/2+1;
24  std::vector<std::complex<double> > coefs(size);
25 
26  for (int jj = 0; jj < kk; ++jj) grid[jj] = xx[jj];
27 
28  fft.fwd(coefs, grid);
29 
30  for (unsigned int jj = 0; jj < size; ++jj) {
31  ff[2*jj] = coefs[jj].real();
32  ff[2*jj+1] = coefs[jj].imag();
33  }
34 }
35 
36 // -----------------------------------------------------------------------------
37 
38 void fft_inv_f(const int & kk, const double * ff, double * xx) {
39  static Eigen::FFT<double> fft;
40  std::vector<double> grid(kk);
41  std::vector<std::complex<double> > coefs(kk);
42 
43  const std::complex<double> zero(0.0, 0.0);
44  const int size = kk/2+1;
45  for (int jj = 0; jj < size; ++jj) {
46  std::complex<double> zz(ff[2*jj], ff[2*jj+1]);
47  coefs[jj] = zz;
48  }
49  for (int jj = size; jj < kk; ++jj) coefs[jj] = zero;
50 
51  fft.inv(grid, coefs);
52 
53  for (int jj = 0; jj < kk; ++jj) xx[jj] = grid[jj];
54 }
55 
56 // -----------------------------------------------------------------------------
57 
58 } // namespace qg
void fft_inv_f(const int &kk, const double *ff, double *xx)
Definition: fft_f.cc:38
real(double), parameter zero
************************************************************************GNU Lesser General Public License **This file is part of the GFDL Flexible Modeling System(FMS). ! *! *FMS is free software without even the implied warranty of MERCHANTABILITY or *FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License *for more details **You should have received a copy of the GNU Lesser General Public *License along with FMS If see< http:! ***********************************************************************subroutine READ_RECORD_CORE_(unit, field, nwords, data, start, axsiz) integer, intent(in) ::unit type(fieldtype), intent(in) ::field integer, intent(in) ::nwords MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in) ::start(:), axsiz(:) integer(SHORT_KIND) ::i2vals(nwords)!rab used in conjunction with transfer intrinsic to determine size of a variable integer(KIND=1) ::one_byte(8) integer ::word_sz!#ifdef __sgi integer(INT_KIND) ::ivals(nwords) real(FLOAT_KIND) ::rvals(nwords)!#else! integer ::ivals(nwords)! real ::rvals(nwords)!#endif real(DOUBLE_KIND) ::r8vals(nwords) pointer(ptr1, i2vals) pointer(ptr2, ivals) pointer(ptr3, rvals) pointer(ptr4, r8vals) if(mpp_io_stack_size< nwords) call mpp_io_set_stack_size(nwords) call mpp_error(FATAL, 'MPP_READ currently requires use_netCDF option') end subroutine READ_RECORD_CORE_ subroutine READ_RECORD_(unit, field, nwords, data, time_level, domain, position, tile_count, start_in, axsiz_in)!routine that is finally called by all mpp_read routines to perform the read!a non-netCDF record contains:! field ID! a set of 4 coordinates(is:ie, js:je) giving the data subdomain! a timelevel and a timestamp(=NULLTIME if field is static)! 3D real data(stored as 1D)!if you are using direct access I/O, the RECL argument to OPEN must be large enough for the above!in a global direct access file, record position on PE is given by %record.!Treatment of timestamp:! We assume that static fields have been passed without a timestamp.! Here that is converted into a timestamp of NULLTIME.! For non-netCDF fields, field is treated no differently, but is written! with a timestamp of NULLTIME. There is no check in the code to prevent! the user from repeatedly writing a static field. integer, intent(in) ::unit, nwords type(fieldtype), intent(in) ::field MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in), optional ::time_level type(domain2D), intent(in), optional ::domain integer, intent(in), optional ::position, tile_count integer, intent(in), optional ::start_in(:), axsiz_in(:) integer, dimension(size(field%axes(:))) ::start, axsiz integer ::tlevel !, subdomain(4) integer ::i, error, is, ie, js, je, isg, ieg, jsg, jeg type(domain2d), pointer ::io_domain=> tlevel if(PRESENT(start_in) .AND. PRESENT(axsiz_in)) then if(size(start(! the data domain and compute domain must refer to the subdomain being passed ! In this ! since that attempts to gather all data on PE size(field%axes(:)) axsiz(i)
The namespace for the qg model.
void fft_fwd_f(const int &kk, const double *xx, double *ff)
Definition: fft_f.cc:20