Title: User guide for library fourpack Author: Leonid Petrov Last modified: 2016.04.28 Library fourpack is for three tasks: a) a convenient user interface to fast Fourier transform kernel either fftw or Intel math kernel library for 1D and 2D cases. b) applying a simple filter in Fourier domain c) direct and inverse spherical harmonic transform on a regular latitude/longitude grid invoking Driscoll and Healy (1994) sampling theorem. It also provides inverse vector spherical harmonics transform. The routines utilizes OpenMP parallelization. The routines do not lose accuracy at high degree/order. They were tested at degree/order 21599. User interface: 1.0 Initialization/quitting. ============================ Internal data structures of fourpack package should be initialized before first usage. There are four procedures for initialization. The first two routines initialize FFT interface. The last two routines initialize both FFT and spherical harmonics transform interface and return the address of the fourpack object that is used for all consecutive calls. 1.1) INIT_FFT_MKL ( NTHREADS, IUER ) -- initializes Intel Math Kernel Library interface for FFT 1.2) INIT_FFTW ( WISDOM_FILE, FFTW_MODE, FFTW_TIMEOUT, NTHREADS, IUER ) -- initializes the FFTW interface for FFT 1.3) SPHE_INIT ( NUM_THR, IUER ) -- initializes the simplified spherical harmonics transform and the FFTW interface for FFT. 1.4) SPHE_INIT_PLAN ( WISDOM_FILE, FFTW_MODE, FFTW_TIMEOUT, NTHREADS, IUER ) -- initializes the full interface for spherical harmonics transform and the FFTW interface for FFT. The allocated data structures are released by call of sphe_quit after the use of Fourpack library. 2.0 Fast Fourier transform. =========================== 2.1) FFT_1D_C16 ( DIM, IEXP, ARR, IUER ) -- performs fast Fourier transform in place, either forward ( IEXP = 1 ), or backward ( IEXP = -1 ) under the complex double precision array ARR of dimension DIM. 2.2) FFT_1D_C8 ( DIM, IEXP, ARR, IUER ) -- the same as FFT_1D_C16, but array ARR is of complex single precision type. 2.3) FFT_2D_C16 ( DIM1, DIM2, IEXP, ARR, IUER ) -- performs fast Fourier transform in place, either forward ( IEXP = 1 ), or backward ( IEXP = -1 ) under the complex double precision 2-dimensional array ARR. 2.4) FFT_2D_C8 ( DIM1, DIM2, IEXP, ARR, IUER ) -- the same as FFT_2D_C16, but array ARR is of complex single precision type. 2.5) FFT_1D_R2C_R8 ( DIM, ARR_R8, ARR_C16, IUER ) -- performs forward fast Fourier transform under the real*8 array ARR_R8 of dimension DIM. Results is complex*16 array ARR_C16 of n/2+1 for non-negative frequencies. 2.6) FFT_1D_R2C_R4 ( DIM, ARR_R4, ARR_C8, IUER ) -- the same as FFT_1D_R2C_R8 but ARR_R4 and ARR_C8 are of real*4 and complex*8 type respectively. 2.7) FFT_1D_C2R_C16 ( DIM, ARR_C16, ARR_R8, IUER ) -- performs backward (inverse) fast Fourier transform under the complex*16 array ARR_C16 of the 1-dimensional Hermitian-conjugated spectrum of dimension DIM. ARR_C16 keeps only dim/2+1 elements for non-negative frequencies. Result, the real*8 array in time order, is put in the output array ARR_R8. 2.8) FFT_1D_C2R_C8 ( DIM, ARR_C8, ARR_R4, IUER ) -- the same as FFT_1D_C2R_C16, but ARR_C8 and ARR_R4 are of complex*8 and real*4 type respectively. 3.0 Power spectrum and filters. =============================== 3.1) POWER_SPECTRUM_R8 ( N, MODE, TIM_ARR, VAL_ARR, FRQ_ARR, POW_ARR, IUER ) -- computes the power spectrum of the real 1D array of data VAL_ARR of dimension of N that is a function of argument TIM_ARR using fast Fourier transform. The result is function POW_ARR of argument FRQ_ARR of dimension N/2+1. Frequency runs from 0 through the Nyquist frequency 2/(TIM_ARR(2)-TIM_ARR(1)) with step 1/(TIM_ARR(N)-TIM_ARR(1)). Depending on MODE, several windows will be applied to the data before Fourier transform. All arguments are of real*8 type. 3.2) POWER_SPECTRUM_R4 ( N, MODE, TIM_ARR, VAL_ARR, & FRQ_ARR, POW_ARR, IUER ) -- the same as POWER_SPECTRUM_R8, but all arguments are of real*8 type. 3.3) ROOT_EDGE_FILTER_R8 ( N, MODE, FRQ_THR, BETA, TIM_ARR, VAL_ARR, IUER ) -- implements either low-pass ( MODE = 1 ), or high=pass ( MODE =2 ) filter. Filtering is performed in three steps: 1) Fourier transform to the frequency domain; 2) multiplication the result by the filter function; 3) inverse Fourier transform of the result. All variables are of real*8 type. 3.4) ROOT_EDGE_FILTER_R4 ( N, MODE, FRQ_THR, BETA, & TIM_ARR, VAL_ARR, IUER ) -- the same as ROOT_EDGE_FILTER_R8, but all arguments are of real*8 type. 3.5) ROOT_BAND_FILTER_R8 ( N, FRQ_LOW, FRQ_HIGH, BETA, TIM_ARR, VAL_ARR, IUER ) -- implements the band-pass digital filter. Filtering is performed in three steps: 1) Fourier transform to the frequency domain; 2) multiplication the result by the filter function; 3) inverse Fourier transform of the result. All variables are of real*8 type. 3.6) ROOT_BAND_FILTER_R4 ( N, FRQ_LOW, FRQ_HIGH, BETA, TIM_ARR, VAL_ARR, IUER ) -- the same as ROOT_BAND_FILTER_R8, but all arguments are of real*8 type. 4.0 Spherical harmonics transform. ================================== 4.1) SPHE_DIR_2NN ( FSH, N, FUN, MD, DEG, NORM, IPHS, SPH, IUER ) expands a grid containing 2*N samples in longitude and N samples in latitude into spherical harmonics. This routine makes use of the sampling theorem presented in Driscoll and Healy (1994) and employs the FFTs when calculating the sin and cos terms. The number of samples, N, must be EVEN for this routine to work, and the spherical harmonic expansion is exact if the function is bandlimited to degree N/2-1. 4.2) SPHE_INV_2NN ( FSH, MD, DEG, NORM, IPHS, SPH, N, FUN, IUER ) -- given the spherical harmonic coefficients SPH of degree/order MD, will evaluate the function on a grid with an equal number of samples N latitude and 2*N in longitude. This is the inverse the routine SPHE_DIR_2NN, both of which are done quickly using FFTs for each degree of each latitude band. The number of samples specified by the spherical harmonic bandwidth DEG may be equal or less than MD. 4.3) SPHE_INV_2NN_VEC ( FSH, MD, DEG, SPH, N, FUN, IUER ) -- evaluates function FUN and its partial derivatives over longitude and latitude on a regular lon/lat grid 2(N-1)*N given its spherical harmonics transform and its tangential derivative. * Four-dimensional array SPH has the first dimension 1:2 that runs over cosine and sine components, second dimension 0:DEG runs over degree, third dimension 0:DEG runs over order, and the fourth dimension 1:2 runs over spherical harmonics and its tangential constituent. SPH contains components nm and mn of degree/order: SPH(a,c,b,d) = SPH(a,b,c,d). The so-called tangential constituent is the Spherical Harmonics Transform of a function that may be different than FUN. 4.4) SPHE_COMP_VAL ( FSH, MD, DEG, LAT_VAL, LON_VAL, NORM, IPHS, SPH, IUER ) -- computes the value at a given latitude and longitude of the function given the set of spherical harmonics of degree/order DEG. Using another language, it computes the value of the the inverse spherical harmonics transform at a given latitude and longitude. 4.5) SPHE_COMP_VEC ( FSH, MD, DEG, LAT_VAL, LON_VAL, NORM, IPHS, SPH, RES, IUER ) -- determines the value and its derivative towards north and east at a given latitude and longitude of the function the given set of spherical harmonics of degree/order DEG. Using another language, it computes the value of the the inverse spherical harmonics transform and its partial derivatives over north and east direction at a given latitude and longitude. 5.0 Auxiliary programs. ====================== 5.1) Program create_fftw_plan creates the so-called wisdom file -- parameter for the tune-up of the FFTW library that performs fast Fourier transform. If hyphen precedes the method, then create_fftw_plan prints timing. Usage: create_fftw_plan where method is MEASURE, PATIENT, or EXHAUSTIVE num_threads -- the number of threads request_file -- the name of the control file with the request Format of the request file: plain ascii file with two 2 or three words separated by one or more blanks 1st word: procedure, one of f -- complex in-place direct/inverse FFT, complex*8 fr2c -- real*4 to complex*8 direct FFT fc2r -- complex*8 to real*4 inverse FFT d -- complex in-place direct/inverse FFT, complex*16 dr2c -- real*8 to complex*16 direct FFT dc2r -- complex*16 to real*8 inverse FFT dim1 -- first dimension; dim2 -- (optional) second dimension. plan_file -- output plan file Using a plan files may increase performance by 1.2-5 times. NB: a plan is thread specific. A plan for K threads is invalid for M threads. References: =========== 1) J.R. Driscoll and D.M. Healy, "Computing Fourier Transforms and Convolutions on the 2-Sphere", (1994), Adv. Applied. Math., 15, 202-250, doi:10.1006/aama.1994.1008. 2) S.A. Holmes, W.E. Featherstone, "A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions, (2002), J. Geodesy, 76, 279-299, 10.1007/s00190-002-0216-2 3) T. Fukushima, "Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers", (2012), J. Geodesy, 86, 271-285, 10.1007/s00190-011-0519-2. 4) M.A. Wieczorek, " SHTOOLS -- Tools for working with spherical harmonics", http://shtools.ipgp.fr/ 6.0 Examples: ============= See source code of two examplse in example/fourpack_example_01.f and example/fourpack_example_02.f See results of timing text in doc/timint 7.0 Interace description: ========================= FUNCTION SPHE_INIT ( NUM_THR, IUER ) Routine SPHE_INIT initializes internal data structure for consecutive computation of direct or inverse spherical transform. It returns the address of the internal data structure and performs some initialization. _________________________ Input parameters: ________________________ NUM_THR ( INTEGER*4 ) -- the number of threads used in subsequent computations. If NUN_THR == -1, then the number of threads is read from the environment variable OMP_NUM_THREADS. _________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. _________________________ Value: ______________________ SPHE_INIT ( INTEGER*8 ) -- address of the created fourpack object --------------------------------------------------------------------------- FUNCTION SPHE_INIT_PLAN ( WISDOM_FILE, FFTW_MODE, FFTW_TIMEOUT, & NTHREADS, IUER ) Routine SPHE_INIT_PLAN initialize internal data structure for consecutive computation of a fast Fourier transform or direct or inverse spherical transform. It returns the address of the internal data structure and performs some initializations. _________________________ Input parameters: ________________________ WISDOM_FILE ( CHARACTER ) -- Name of the so-called wisdom file that contains tune-up parameters of FFTW library. The file name may be blank -- in that case no wisdom is acquired. FFTW_MODE ( INTEGER*4 ) -- Used mode of FFTW. The following modes are supported: 0 -- ESTIMATE -- provides a reasonable speed for execution and very fast plan evaluation. Recommended for using unless, an evidence exists that other modes will be faster. 1 -- MEASURE -- may faster than ESTIMATE, but may be 100 times slower for startup. Recommended only with a wisdom file. 2 -- PATIENT -- may faster than ESTIMATE, nut may be 1000 times slower for startup. Recommended only with a wisdom file. 3 -- EXHAUSTIVE -- may faster than ESTIMATE, but it may take 10000 times slower for a startup. You should be crazy if you try it without a wisdom file. FFTW_TIMEOUT ( REAL*8 ) -- Maximum amount of time in seconds for creation of plan. FFTW may be so slow in deciding which plan to use that you will certainly want to limit this time. 0.01 is good value for dimensions less than 65536. NTHREADS ( INTEGER*4 ) -- The number of threads to be used. Value 0 means no threads. Value -1 means to get the number of threads from the environment variable OMP_NUM_THREADS. NB: FFTW wisdom created for a non-thread environment will be discarded if the number of threads > 1. _________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. _________________________ Value: ______________________ SPHE_INIT_PLAN ( INTEGER*8 ) -- address of the created fourpack object ------------------------------------------------------------------------- SUBROUTINE INIT_FFT_MKL ( NTHREADS, IUER ) Routine INIT_FFT_MKL initializes the interface to fast FFT routines supplied by the Intel proprietary MKL library. Upon calling INIT_FFT_MKL, all subsequent calls to routines of the FAST_FFT library will be served bu MKL subroutines. NB: in order to compile this routine you need to create file mkl_dfti by this command: cd ${path_to_intel}/intel/mkl/include ${fortran_compiler} -c mkl_dfti.f90 Besides, you need to add -I ${path_to_intel}/intel/mkl/include to the compiler options in order to find the mkl_dfti module. _________________________ Input parameters: ________________________ NTHREADS ( INTEGER*4 ) -- The number of threads to be used. Value 0 means no threads. Value -1 means to get the number of threads from the environment variable OMP_NUM_THREADS. NB: FFTW wisdom created for a non-thread environment will be discarded if the number of threads > 1. _________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- FUNCTION INIT_FFTW ( WISDOM_FILE, FFTW_MODE, FFTW_TIMEOUT, & NTHREADS, IUER ) Routine INIT_FFTW reads the so-called WISDOM_FILE with parameters of the tune-up of the FFTW library that implements fast Fourier transform and performs initialization of threads. It reads OMP_NUM_THREADS and determines the number of threads. If OMP_NUM_THREADS is not defined, one thread will be used. The wisdom file should be created by a separate program create_fftw_wisdom. _________________________ Input parameters: ________________________ WISDOM_FILE ( CHARACTER ) -- Name of the so-called wisdom file that contains tune-up parameters of FFTW library. The file name may be blank -- in that case no wisdom is acquired. FFTW_MODE ( INTEGER*4 ) -- Used mode of FFTW. The following modes are supported: 0 -- ESTIMATE -- provides a reasonable speed for execution and very fast plan evaluation. Recommended for using unless, an evidence exists that other modes will be faster. 1 -- MEASURE -- may faster than ESTIMATE, but may be 100 times slower for startup. Recommended only with a wisdom file. 2 -- PATIENT -- may faster than ESTIMATE, nut may be 1000 times slower for startup. Recommended only with a wisdom file. 3 -- EXHAUSTIVE -- may faster than ESTIMATE, but it may take 10000 times slower for a startup. You should be crazy if you try it without a wisdom file. FFTW_TIMEOUT ( REAL*8 ) -- Maximum amount of time in seconds for creation of plan. FFTW may be so slow in deciding which plan to use that you will certainly want to limit this time. 0.01 is good value for dimensions less than 65536. NTHREADS ( INTEGER*4 ) -- The number of threads to be used. Value 0 means no threads. Value -1 means to get the number of threads from the environment variable OMP_NUM_THREADS. NB: FFTW wisdom created for a non-thread environment will be discarded if the number of threads > 1. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. _________________________ Value: ___________________________________ INIT_FFTW ( INTEGER*4 ) -- returns the number of threads initialized. ------------------------------------------------------------------------- SUBROUTINE SPHE_QUIT ( FSH ) Routine SPHE_QUIT releases dynamic memory allocated inside the initialize internal data structure for consecutive computation of direct or inverse spherical transform and then releases memory allocated for this structure itself. Input parameters: FSH ( SPHE_TYPE ) -- Internal data structure of fourpack package that keeps internal arrays with intermediate results and their status for possible re-use. ------------------------------------------------------------------------- SUBROUTINE FFT_1D_C16 ( DIM, IEXP, ARR, IUER ) Routine FFT_1D_C16 performs fast Fourier transform in place, either forward ( IEXP = 1 ), or backward ( IEXP = -1 ) under the complex double precision array ARR of dimension DIM. _________________________ Input parameters: ________________________ DIM ( INTEGER*4 ) -- Dimension of the array. NB: dimension that is the power of 2 yields the best performance. IEXP ( INTEGER*4 ) -- Direction of the transform: 1 -- forward (direct) FFT. -1 -- backward (reverse) FFT. ________________________ Modified parameters: ______________________ ARR ( COMPLEX*16 ) -- Array to be transformed. Input for forward FFT: array ordered in time.* Output of the backward FFT in this order: [1, dim/2-1] -- increase in frequency order starting from the zero frequency [dim/2, dim] -- increase in frequency starting from -n/2 through 1/n IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_1D_C8 ( DIM, IEXP, ARR, IUER ) Routine FFT_1D_C8 performs fast Fourier transform in place, either forward ( IEXP = 1 ), or backward ( IEXP = -1 ) under the complex single precision array ARR of dimension DIM. _________________________ Input parameters: ________________________ DIM ( INTEGER*4 ) -- Dimension of the array. NB: dimension that is the power of 2 yields the best performance. IEXP ( INTEGER*4 ) -- Direction of the transform: 1 -- forward (direct) FFT. -1 -- backward (reverse) FFT. ________________________ Modified parameters: ______________________ ARR ( COMPLEX*8 ) -- Array to be transformed. Input for forward FFT: array ordered in time.* Output of the backward FFT in this order: [1, dim/2-1] -- increase in frequency order starting from the zero frequency [dim/2, dim] -- increase in frequency starting from -n/2 through 1/n IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_2D_C16 ( DIM1, DIM2, IEXP, ARR, IUER ) Routine FFT_2D_C16 performs fast Fourier transform in place, either forward ( IEXP = 1 ), or backward ( IEXP = -1 ) under the complex double precision 2-dimensional array ARR. _________________________ Input parameters: ________________________ DIM1 ( INTEGER*4 ) -- 1st dimension of the array. NB: dimension that is the power of 2 yields the best performance. DIM2 ( INTEGER*4 ) -- 2nd dimension of the array. NB: dimension that is the power of 2 yields the best performance. IEXP ( INTEGER*4 ) -- Direction of the transform: 1 -- forward (direct) FFT. -1 -- backward (reverse) FFT. ________________________ Modified parameters: ______________________ ARR ( COMPLEX*16 ) -- Array to be transformed. Input for forward FFT: array ordered in time. Output of the backward FFT in this order: [1, dim/2-1] -- increase in frequency order starting from the zero frequency [dim/2, dim] -- increase in frequency starting from -n/2 through 1/n IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_2D_C8 ( DIM1, DIM2, IEXP, ARR, IUER ) Routine FFT_2D_C8 performs fast Fourier transform in place, either forward ( IEXP = 1 ), or backward ( IEXP = -1 ) under the complex single precision 2-dimensional array ARR. _________________________ Input parameters: ________________________ DIM1 ( INTEGER*4 ) -- 1st dimension of the array. NB: dimension that is the power of 2 yields the best performance. DIM2 ( INTEGER*4 ) -- 2nd dimension of the array. NB: dimension that is the power of 2 yields the best performance. IEXP ( INTEGER*4 ) -- Direction of the transform: 1 -- forward (direct) FFT. -1 -- backward (reverse) FFT. ________________________ Modified parameters: ______________________ ARR ( COMPLEX*8 ) -- Array to be transformed. Input for forward FFT: array ordered in time.* Output of the backward FFT in this order: [1, dim/2-1] -- increase in frequency order starting from the zero frequency [dim/2, dim] -- increase in frequency starting from -n/2 through 1/n IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_1D_R2C_R8 ( DIM, ARR_R8, ARR_C16, IUER ) Routine FFT_1D_R2C_R8 performs forward fast Fourier transform under the real*8 array ARR_R8 of dimension DIM. Results is complex*16 array ARR_C16 of n/2+1 for non-negative frequencies. _________________________ Input parameters: ________________________ DIM ( INTEGER*4 ) -- Dimension of the array input. NB: dimension that is the power of 2 yields the best performance. ARR_R8 ( REAL*8 ) -- Array to be transformed. Should be time ordered. _________________________ Output parameters: _______________________ ARR_C16 ( COMPLEX*16 ) -- Fourier transform of the input array. Dimension: DIM/2+1. Frequencies start from 0 and run through dim/2 ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_1D_R2C_R4 ( DIM, ARR_R4, ARR_C8, IUER ) Routine FFT_1D_R2C_R4 performs forward fast Fourier transform under the real*4 array ARR_R4 of dimension DIM. Results is complex*16 array ARR_C8 of n/2+1 for non-negative frequencies. _________________________ Input parameters: ________________________ DIM ( INTEGER*4 ) -- Dimension of the array input. NB: dimension that is the power of 2 yields the best performance. ARR_R4 ( REAL*4 ) -- Real array to be transformed. Should be time ordered. _________________________ Output parameters: _______________________ ARR_C8 ( COMPLEX*8 ) -- Fourier transform of the input array. Dimension: DIM/2+1. Frequencies start from 0 and run through dim/2 ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_1D_C2R_C16 ( DIM, ARR_C16, ARR_R8, IUER ) Routine FFT_1D_C2R_C16 performs backward (inverse) fast Fourier transform under the complex*16 array ARR_C16 of the 1-dimensional Hermtian-conjugated spectrum of dimension DIM. ARR_C16 keeps only dim/2+1 elements for non-negative frequencies. Result, the real*8 array in time order, is put in the output array ARR_R8. _________________________ Input parameters: ________________________ DIM ( INTEGER*4 ) -- Dimension of the array input. NB: dimension that is the power of 2 yields the best performance. ARR_C16 ( COMPLEX*16 ) -- Array to be transformed. A portion of the Hermitian-conjugated spectrum for non-negative frequencies. Dimension: DIM/2+1. Frequencies start from 0 and run through dim/2. _________________________ Output parameters: _______________________ ARR_R8 ( REAL*8 ) -- Array to be transformed. Should be time ordered. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE FFT_1D_C2R_C8 ( DIM, ARR_C8, ARR_R4, IUER ) Routine FFT_1D_C2R_C8 performs backward (inverse) fast Fourier transform under the complex*8 array ARR_C8 of the 1-dimensional Hermtian-conjugated spectrum of dimension DIM. ARR_C8 keeps only dim/2+1 elements for non-negative frequencies. Result, the real*8 array in time order, is put in the output array ARR_R8. _________________________ Input parameters: ________________________ DIM ( INTEGER*4 ) -- Dimension of the array input. NB: dimension that is the power of 2 yields the best performance. ARR_C8 ( COMPLEX*8 ) -- Array to be transformed. A portion of the Hermitian-conjugated spectrum for non-negative frequencies. Dimension: DIM/2+1. Frequencies start from 0 and run through dim/2. _________________________ Output parameters: _______________________ ARR_R4 ( REAL*4 ) -- Array to be transformed. Should be time ordered. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE POWER_SPECTRUM_R8 ( N, MODE, TIM_ARR, VAL_ARR, & FRQ_ARR, POW_ARR, IUER ) Routine POWER_SPECTRUM_R8 computes the power spectrum of the real 1D array of data VAL_ARR of dimension of N that is a function of argument TIM_ARR using fast Fourier transform. The result is function POW_ARR of argument FRQ_ARR of dimension N/2+1. Frequency runs from 0 through the Nyquest frequency 2/(TIM_ARR(2)-TIM_ARR(1)) with step 1/(TIM_ARR(N)-TIM_ARR(1)). Depending on MODE, several windows will be applied to the data before Fourier transform. _________________________ Input parameters: ________________________ N ( INTEGER*4 ) -- The number of points of the input array. MODE ( INTEGER*4 ) -- Windowing mode. 0 -- no window applied. 1 -- Hann window (recommended). 2 -- Hamming window. TIM_ARR ( REAL*8 ) -- Argument of the input data. It may not be necessarily time. Dimension: N. NB: the argument MUST be ordered in ascending order and MUST have the equidistant step. Dimension: N. VAL_ARR ( REAL*8 ) -- The function, which power spectrum is to be computed. _________________________ Output parameters: _______________________ FRQ_ARR ( REAL*8 ) -- Argument of the power spectrum. Runs from 0.0 through 2/(TIM_ARR(2)-TIM_ARR(1)) with step 1/(N*(TIM_ARR(N)-TIM_ARR(1))). Dimension: N/2+1. POW_ARR ( REAL*8 ) -- The power spectrum. Dimension: N/2+1. Unit: square of val_arr. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE POWER_SPECTRUM_R4 ( N, MODE, TIM_ARR, VAL_ARR, & FRQ_ARR, POW_ARR, IUER ) Routine POWER_SPECTRUM_R4 computes the power spectrum of the real 1D array of data VAL_ARR of dimension of N that is a function of argument TIM_ARR using fast Fourier transform. The result is function POW_ARR of argument FRQ_ARR of dimension N/2+1. Frequency runs from 0 through the Nyquest frequency 2/(TIM_ARR(2)-TIM_ARR(1)) with step 1/(TIM_ARR(N)-TIM_ARR(1)). Depending on MODE, several windows will be applied to the data before Fourier transform. _________________________ Input parameters: ________________________ N ( INTEGER*4 ) -- The number of points of the input array. MODE ( INTEGER*4 ) -- Windowing mode. 0 -- no window applied. 1 -- Hann window (recommended). 2 -- Hanning window. TIM_ARR ( REAL*4 ) -- Argument of the input data. It may not be necessarily time. Dimension: N. NB: the argument MUST be ordered in ascending order and MUST have the equidistant step. Dimension: N. VAL_ARR ( REAL*4 ) -- The function, which power spectrum is to be computed. _________________________ Output parameters: _______________________ FRQ_ARR ( REAL*4 ) -- Argument of the power spectrum. Runs from 0.0 through 2/(TIM_ARR(2)-TIM_ARR(1)) with step 1/(N*(TIM_ARR(N)-TIM_ARR(1))). Dimension: N/2+1. POW_ARR ( REAL*4 ) -- The power spectrum. Dimension: N/2+1. Unit: square of val_arr. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE ROOT_EDGE_FILTER_R8 ( N, MODE, FRQ_THR, BETA, & TIM_ARR, VAL_ARR, IUER ) Routine ROOT_EDGE_FILTER_R8 implements either low-pass ( MODE = 1 ), or high=pass ( MODE =2 ) filter. Filtering is performed in three steps: 1) Fourier transform to the frequency domain; 2) multiplication the result by the filter function; 3) inverse Fourier transform of the result. The following filter function H(f) for the high-pass filter is used: 1 beta*(f + f_o)/f_o beta*(f - f_o)/f_o - * --------------------------------- - ---------------------------------- 2 sqrt(1 + (beta*(f + f_o)/f_o))**2) sqrt(1 + (beta*(f - f_o)/f_o))**2) The low-pass filter is 1 - H(f) Parameter f_o determines the cut-off frequency, and beta determines the sharpness of the filter: the higher, the sharper. beta=20 is recommended for many applications. _________________________ Input parameters: ________________________ N ( INTEGER*4 ) -- The number of points of the input array. MODE ( INTEGER*4 ) -- Filter mode: 1 -- low-pass filter ( frequencies f < frq_thr are passed, higher frequencies are removed ). 2 -- high pass filter ( frequencies f > frq_thr are passed, lower frequencies are removed ). FRQ_THR ( REAL*8 ) -- The cut-off cyclic frequency of the filter. Units: reciprocal to TIM_ARR units. BETA ( REAL*8 ) -- Parameter that determines the sharpness of the filter. The higher beta, the closer the filter shape in the frequency domain to the capital greek letter "Pi" shape. NB: a very sharp filter results in an extended base in the time domain. TIM_ARR ( REAL*8 ) -- Argument of the input data. It may not be necessarily time. Dimension: N. NB: the argument MUST be ordered in ascending order and MUST have the equidistant step. Dimension: N. ________________________ Modified parameters: ______________________ VAL_ARR ( REAL*8 ) -- The function that is to be transformed. Dimension: N. IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE ROOT_EDGE_FILTER_R4 ( N, MODE, FRQ_THR, BETA, & TIM_ARR, VAL_ARR, IUER ) Routine ROOT_EDGE_FILTER_R4 implements either low-pass ( MODE = 1 ), or high=pass ( MODE =2 ) filter. Filtering is performed in three steps: 1) Fourier transfrom to the frequency domain; 2) multiplication the result by the filter function; 3) inverse Fourier transform of the result. The following filter function H(f) for the high-pass filter is used: 1 beta*(f + f_o)/f_o beta*(f - f_o)/f_o - * --------------------------------- - ---------------------------------- 2 sqrt(1 + (beta*(f + f_o)/f_o))**2) sqrt(1 + (beta*(f - f_o)/f_o))**2) The low-pass filter is 1 - H(f) Parameter f_o determines the cut-off frequency, and beta determines the sharpness of the filter: the higher, the sharper. beta=20 is recommended for many applications. _________________________ Input parameters: ________________________ N ( INTEGER*4 ) -- The number of points of the input array. MODE ( INTEGER*4 ) -- Filter mode: 1 -- low-pass filter ( frequencies f < frq_thr are passed, higher frequencies are removed ). 2 -- high pass filter ( frequencies f > frq_thr are passed, lower frequencies are removed ). FRQ_THR ( REAL*4 ) -- The cut-off cyclic frequency of the filter. Units: reciprocal to TIM_ARR units. BETA ( REAL*4 ) -- Parameter that determines the sharpness of the filter. The higher beta, the closer the filter shape in the frequency domain to the capital greek letter "Pi" shape. NB: a very sharp filter results in an extended base in the time domain. TIM_ARR ( REAL*4 ) -- Argument of the input data. It may not be necessarily time. Dimension: N. NB: the argument MUST be ordered in ascending order and MUST have the equidistant step. Dimension: N. ________________________ Modified parameters: ______________________ VAL_ARR ( REAL*4 ) -- The function that is to be transformed. Dimension: N. IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE ROOT_BAND_FILTER_R8 ( N, FRQ_LOW, FRQ_HIGH, BETA, & TIM_ARR, VAL_ARR, IUER ) Routine ROOT_BAND_FILTER_R8 implements the band-pass digital filter. Filtering is performed in three steps: 1) Fourier transform to the frequency domain; 2) multiplication the result by the filter function; 3) inverse Fourier transform of the result. Filter function K(f,f_l,f_h) = (1 - H(f,f_l)*H(f,f_h), where H(f,f_o): 1 beta*(f + f_o)/f_o beta*(f - f_o)/f_o - * --------------------------------- - ---------------------------------- 2 sqrt(1 + (beta*(f + f_o)/f_o))**2) sqrt(1 + (beta*(f - f_o)/f_o))**2) Parameter f_l, f_h determines the low and high cut-off frequencies, while parameter beta determines the sharpness of the filter: the higher, the sharper. beta=20 is recommended for many applications. _________________________ Input parameters: ________________________ N ( INTEGER*4 ) -- The number of points of the input array. FRQ_LOW ( REAL*8 ) -- The low cut-off cyclic frequency of the filter. Frequencies f < FRE_LOW will not be passed. Units: reciprocal to TIM_ARR units. FRQ_HIGH ( REAL*8 ) -- The high cut-off cyclic frequency of the filter. Frequencies f > FRE_HIGH will not be passed. Units: reciprocal to TIM_ARR units. BETA ( REAL*8 ) -- Parameter that determines the sharpness of the filter. The higher beta, the closer the filter shape in the frequency domain to the capital greek letter "Pi" shape. NB: a very sharp filter results in an extended base in the time domain. TIM_ARR ( REAL*8 ) -- Argument of the input data. It may not be necessarily time. Dimension: N. NB: the argument MUST be ordered in ascending order and MUST have the equidistant step. Dimension: N. ________________________ Modified parameters: ______________________ VAL_ARR ( REAL*8 ) -- The function that is to be transformed. Dimension: N. IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE ROOT_BAND_FILTER_R4 ( N, FRQ_LOW, FRQ_HIGH, BETA, & TIM_ARR, VAL_ARR, IUER ) Routine ROOT_BAND_FILTER_R4 implements the band-pass digital filter. Filtering is performed in three steps: 1) Fourier transform to the frequency domain; 2) multiplication the result by the filter function; 3) inverse Fourier transform of the result. Filter function K(f,f_l,f_h) = (1 - H(f,f_l)*H(f,f_h), where H(f,f_o): 1 beta*(f + f_o)/f_o beta*(f - f_o)/f_o - * --------------------------------- - ---------------------------------- 2 sqrt(1 + (beta*(f + f_o)/f_o))**2) sqrt(1 + (beta*(f - f_o)/f_o))**2) Parameter f_l, f_h determines the low and high cut-off frequencies, while parameter beta determines the sharpness of the filter: the higher, the sharper. beta=20 is recommended for many applications. _________________________ Input parameters: ________________________ N ( INTEGER*4 ) -- The number of points of the input array. FRQ_LOW ( REAL*4 ) -- The low cut-off cyclic frequency of the filter. Frequencies f < FRE_LOW will not be passed. Units: reciprocal to TIM_ARR units. FRQ_HIGH ( REAL*4 ) -- The high cut-off cyclic frequency of the filter. Frequencies f > FRE_HIGH will not be passed. Units: reciprocal to TIM_ARR units. BETA ( REAL*4 ) -- Parameter that determines the sharpness of the filter. The higher beta, the closer the filter shape in the frequency domain to the capital greek letter "Pi" shape. NB: a very sharp filter results in an extended base in the time domain. TIM_ARR ( REAL*4 ) -- Argument of the input data. It may not be necessarily time. Dimension: N. NB: the argument MUST be ordered in ascending order and MUST have the equidistant step. Dimension: N. ________________________ Modified parameters: ______________________ VAL_ARR ( REAL*4 ) -- The function that is to be transformed. Dimension: N. IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE SPHE_DIR_2NN ( FSH, N, FUN, MD, DEG, NORM, IPHS, SPH, IUER ) This routine will expand a grid containing 2*N samples in longitude and N samples in latitude into spherical harmonics. This routine makes use of the sampling theorem presented in Driscoll and Healy (1994) and employs FFTs when calculating the sin and cos terms. The number of samples, N, must be EVEN for this routine to work, and the spherical harmonic expansion is exact if the function is bandlimited to degree N/2-1. In a case if the dimension of transform is less or equal FSH__MAX_SCL, then the Legendre associated functions are computed on the fly using the scaling methodology presented in Holmes and Featherston (2002). When NORM is 1,2 or 4, these are accurate to about degree 2800. When NORM is 3, the routine is only stable to about degree 15. In a case if the dimension of transform is greater than FSH__MAX_SCL=2700, then the Legendre functions are computed on the fly using the X-number algebra introduced by Fukushima (2011). The essence of that approach is that the intermediate variables needed for computation of Legendre associated functions are kept in two numbers: the double precision float number and integer number for additional exponent. The input grid contains N samples in latitude from -90 to 90-interval, and 2*N samples in longitude from 0 to 360-interval, where interval is the latitudinal sampling interval 180/N. The input grid must contain N samples in latitude and 2N samples in longitude. The sampling intervals in latitude and longitude are 180/N and 360/N respectively. When performing the FFTs in longitude, the frequencies greater than N/2-1 are simply discarded to prevent aliasing. Notes: 1. This routine does not use the fast Legendre transforms that are presented in Driscoll and Heally (1994). 2. Use of a N by 2N grid is implemented because many geographic grids are sampled this way. When taking the Fourier transforms in longitude, all of the higher frequencies are ultimately discarded. If, instead, every other column of the grid were discarded to form a NxN grid, higher frequencies could be aliased into lower frequencies. ________________________ Input parameters: _________________________ FSH ( SPHE_TYPE ) -- Internal data structure that keeps internal arrays with intermediate results and their status for possible re-use. N ( INTEGER*4 ) -- Dimension of the input function along latitude. FUN ( REAL*8 ) -- The function that is to be expanded. Dimension 2N*N (longitude,latitude). MD ( INTEGER*4 ) -- Dimension of the array for spherical function. DEG ( INTEGER*4 ) -- Maximum degree of the transform. Should not exceed MD. NORM ( INTEGER*4 ) -- Normalization to be used when calculating Legendre functions 1 -- "geodesy"; 2 -- Schmidt; 3 -- unnormalized; 4 -- orthonormalized; IPHS ( INTEGER*4 ) -- Phase flag. 1: Do not include the Condon-Shortley phase factor of (-1)^m. -1: Apply the Condon-Shortley phase factor of (-1)^m. ________________________ Output parameters: ________________________ SPH ( REAL*8 ) -- Array with spherical transform coefficients. Dimension: (2,0:MD,0:MD). The first dimesion runs over cosine/sine component, the second dimension runs over order degree m, the third dimension runs over order l. NB: only coefficients l =< m are filled! The part of array SPH l > m is filled with zeroes. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- FUNCTION SPHE_INV_2NN ( FSH, MD, DEG, NORM, IPHS, SPH, N, FUN, IUER ) Routine SPHE_INV_2NN given the spherical harmonic coefficients SPH of degree/order MD, will evaluate the function on a grid with an equal number of samples N latitude and 2*N in longitude. This is the inverse the routine SPHE_DIR_2NN, both of which are done quickly using FFTs for each degree of each latitude band. The number of samples specified by the spherical harmonic bandwidth DEG may be equal or less than MD. Note that N is always EVEN for this routine. In a case if the dimension of transform is less or equal FSH__MAX_SCL, then the Legendre associated functions are computed on the fly using the scaling methodology presented in Holmes and Featherston (2002). When NORM is 1,2 or 4, these are accurate to about degree 2800. When NORM is 3, the routine is only stable to about degree 15. In a case if the dimension of transform is greater than FSH__MAX_SCL, then the Legendre functions are computed on the fly using the X-number algebra introduced by Fukushima (2011). The essence of that approach is that the intermediate variables needed for computation of Legendre associated functions are kept in two numbers: the double precision float number and the integer number for additional exponent. The output grid contains N samples in latitude from -90 to +90 deg interval, including south pole, but excluding north pole and in longitude from 0 to 360-1*interval (or 2(N-1) x N), where interval is the sampling interval, and n=2*(deg+1)+1. NB: the grid excludes the northern pole. The output grid contains N samples in latitude from -90 to 90-interval, and in longitude from 0 to 360-2*interval (or 2N x N), where interval is the sampling interval, and n=2*(deg+1). ________________________ Input parameters: _________________________ FSH ( SPHE_TYPE ) -- Internal data structure that keeps internal arrays with intermediate results and their status for possible re-use. MD ( INTEGER*4 ) -- Dimension of the array for spherical function. DEG ( INTEGER*4 ) -- Maximum degree of the transform. Should not exceed MD. NORM ( INTEGER*4 ) -- Normalization to be used when calculating Legendre functions 1 -- "geodesy"; 2 -- Schmidt; 3 -- unnormalized; 4 -- orthonormalized; IPHS ( INTEGER*4 ) -- Phase flag. 1: Do not include the Condon-Shortley phase factor of (-1)^m. -1: Apply the Condon-Shortley phase factor of (-1)^m. SPH ( REAL*8 ) -- Array with spherical transform coefficients. Dimension: (2,0:MD,0:MD). The first dimension runs over cosine/sine component, the second dimension runs degree m, the third dimension runs over order l. NB: only coefficients l =< m are filled! The part of array SPH l > m is filled with zeroes. N ( INTEGER*4 ) -- Dimension of the output function along latitude. ________________________ Output parameters: ________________________ FUN ( REAL*8 ) -- The function that is to be expanded. Dimension 2N*N (longitude,latitude). ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE SPHE_INV_2NN_VEC ( FSH, MD, DEG, SPH, N, FUN, IUER ) Routine SPHE_INV_2NN_VEC evaluates function FUN and its partial derivatives over longitude and latitude on a regular lon/lat grid 2(N-1)*N given its spherical harmonics transform and its tangential derivative. Four-dimensional array SPH has the first dimension 1:2 that runs over cosine and sine components, second dimension 0:DEG runs over degree, third dimension 0:DEG runs over order, and the fourth dimension 1:2 runs over spherical harmonics and its tangential constituent. SPH contains components nm and mn of degree/order: SPH(a,c,b,d) = SPH(a,b,c,d). The so-called tangential constituent is the Spherical Harmonics Transform of a function that may be different than FUN. The output is placed in three-dimensional array FUN defined on a regular longitude/latitude grid at the unit sphere. Its first slide is the value of the function restored from the spherical harmonics defined in the 1st slice of SPH. The second slice of FUN is the partial derivative of function defined by its spherical harmonics stored in the 2nd slice of SPH over the longitudinal direction. The third slice of FUN is the partial derivative of the function defined by its spherical harmonics stored in the 2nd slice of SPH over the latitudinal direction. FUN(i,j,1) = \sum_n \sum_m SPH(b,c,1) * Y_mn FUN(i,j,2) = \sum_n \sum_m SPH(b,c,2) * \d Y_mn \d long / cos\phi FUN(i,j,3) = \sum_n \sum_m SPH(b,c,2) * \d Y_mn \d latitude This is the inverse of the routine SPHE_DIR_2NN, both of which are done quickly using FFTs for each degree of each latitude band. The number of samples is determined by the spherical harmonic bandwidth DEG may be equal or less than MD. Note that N is always ODD for this routine. In a case if the dimension of transform is less or equal FSH__MAX_SCL, then the Legendre associated functions are computed on the fly using the scaling methodology presented in Holmes and Featherston (2002). When NORM is 1,2 or 4, these are accurate to about degree 2800. When NORM is 3, the routine is only stable to about degree 15. In a case if the dimension of transform is greater than FSH__MAX_SCL, then the Legendre functions are computed on the fly using the X-number algebra introduced by Fukushima (2011). The essence of that approach is that the intermediate variables needed for computation of Legendre associated functions are kept in two numbers: the double precision float number and integer number for additional exponent. Geodesy normalization for spherical harmonics us used. The output grid contains N samples in latitude from -90 to +90 deg interval, including both poles, and in longitude from 0 to 360-1*interval (or 2(N-1) x N), where interval is the sampling interval, and n=2*(deg+1)+1. NB: the grid includes both poles. ________________________ Input parameters: _________________________ FSH ( SPHE_TYPE ) -- Internal data structure that keeps internal arrays with intermediate results and their status for possible re-use. MD ( INTEGER*4 ) -- Dimension of the array for spherical function. DEG ( INTEGER*4 ) -- Maximum degree of the transform. Should not exceed MD. SPH ( REAL*8 ) -- Array with spherical transform coefficients. Dimension: (2,0:MD,0:MD,2). The first dimension runs over cosine/sine component, the second dimension runs degree m, the third dimension runs over order l. The fourth dimension runs over 1 and 2. N ( INTEGER*4 ) -- Dimension of the output function along latitude. ________________________ Output parameters: ________________________ FUN ( REAL*8 ) -- The function that is to be expanded. Dimension 2N*N (longitude,latitude). ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- FUNCTION SPHE_COMP_VAL ( FSH, MD, DEG, LAT_VAL, LON_VAL, NORM, & IPHS, SPH, IUER ) Routine SPHE_COMP_VAL determines the value at a given latitude and longitude corresponding to the given set of spherical harmonics. ________________________ Input parameters: _________________________ FSH ( SPHE_TYPE ) -- Internal data structure that keeps internal arrays with intermediate results and their status for possible re-use. MD ( INTEGER*4 ) -- Dimension of the array for spherical function. DEG ( INTEGER*4 ) -- Maximum degree of the transform. Should not exceed MD. LAT_VAL ( REAL*8 ) -- Latitude of the point in radians. LON_VAL ( REAL*8 ) -- Longitude of the point in radians. NORM ( INTEGER*4 ) -- Normalization to be used when calculating Legendre functions 1 -- "geodesy"; 2 -- Schmidt; 3 -- unnormalized; 4 -- orthonormalized; IPHS ( INTEGER*4 ) -- Phase flag. 1: Do not include the Condon-Shortley phase factor of (-1)^m. -1: Apply the Condon-Shortley phase factor of (-1)^m. SPH ( REAL*8 ) -- Array with spherical transform coefficients. Dimension: (2,0:MD,0:MD). The first dimension runs over cosine/sine component, the second dimension runs over order l, the third dimension runs over degree. NB: only coefficients l =< m are filled! The part of array SPH l > m is filled with zeroes. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. ------------------------------------------------------------------------- SUBROUTINE SPHE_COMP_VEC ( FSH, MD, DEG, LAT_VAL, LON_VAL, NORM, IPHS, SPH, RES, IUER ) Routine SPHE_COMP_VEC determines the value and its derivative towards north and east at a given latitude and longitude corresponding to the given set of spherical harmonics. NB: spherical harmonics SPH array is transposed with respect to intuitive expectation: the second dimension runs over order m and the third dimension runs over degree l. ________________________ Input parameters: _________________________ FSH ( SPHE_TYPE ) -- Internal data structure that keeps internal arrays with intermediate results and their status for possible re-use. MD ( INTEGER*4 ) -- Dimension of the array for spherical function. DEG ( INTEGER*4 ) -- Maximum degree of the transform. Should not exceed MD. LAT_VAL ( REAL*8 ) -- Latitude of the point in radians. LON_VAL ( REAL*8 ) -- Longitude of the point in radians. NORM ( INTEGER*4 ) -- Normalization to be used when calculating Legendre functions. Presently, only 1 -- "geodesy" is supported IPHS ( INTEGER*4 ) -- Phase flag. 1: Do not include the Condon-Shortley phase factor of (-1)^m. -1: Apply the Condon-Shortley phase factor of (-1)^m. SPH ( REAL*8 ) -- Array with spherical transform coefficients. Dimension: (2,0:MD,0:MD,2). The first dimension runs over cosine/sine component, the second dimension runs over order m, the third dimension runs over degree l, the fourth dimension runs over radial and transverse components. NB: only coefficients l >= m are filled! The part of array SPH l < m is filled with zeroes. RES ( REAL*8 ) -- Vector of results: value, north derivative, east derivative. Dimension: 3. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error.