Line data Source code
1 : module pft_module
2 :
3 : ! This module provides fast-Fourier transforms
4 : !
5 : ! REVISION HISTORY:
6 : ! 01.01.30 Lin Integrated into this module
7 : ! 05.05.25 Sawyer Merged CAM and GEOS5 versions (CAM vectorization)
8 : ! 05.07.26 Worley Revised module using for Cray X1 version
9 :
10 :
11 : use shr_kind_mod, only: r8 => shr_kind_r8
12 :
13 : implicit none
14 : private
15 : save
16 :
17 : #ifdef NO_R16
18 : integer, parameter :: r16= selected_real_kind(12) ! 8 byte real
19 : #else
20 : integer, parameter :: r16= selected_real_kind(24) ! 16 byte real
21 : #endif
22 :
23 : public :: pft2d, pft_cf, fftfax, pftinit, fftrans
24 :
25 : real(r8), parameter :: D0_0 = 0.0_r8
26 : real(r8), parameter :: D1EM20 = 1.0e-20_r8
27 : real(r8), parameter :: D0_5 = 0.5_r8
28 : real(r8), parameter :: D1_0 = 1.0_r8
29 : real(r8), parameter :: D1_01 = 1.01_r8
30 : real(r8), parameter :: D2_0 = 2.0_r8
31 : real(r8), parameter :: D4_0 = 4.0_r8
32 : real(r8), parameter :: D8_0 = 8.0_r8
33 : real(r8), parameter :: D180_0 =180.0_r8
34 :
35 : integer :: ifax(13) !ECMWF fft
36 : real(r8), allocatable :: trigs(:) ! reentrant code??
37 :
38 : integer :: fft_flt ! 0 => FFT/algebraic filter; 1 => FFT filter
39 :
40 : !=========================================================================================
41 : CONTAINS
42 : !=========================================================================================
43 :
44 1536 : subroutine pftinit(im, fft_flt_in)
45 :
46 : ! Two-dimensional FFT initialization
47 :
48 : ! arguments
49 : integer, intent(in) :: im ! Total X dimension
50 : integer, intent(in) :: fft_flt_in
51 :
52 : ! local variables
53 : integer :: icffta
54 : real(r8) :: rcffta
55 : !----------------------------------------------------------------------------
56 :
57 1536 : fft_flt = fft_flt_in
58 :
59 : #if defined( LIBSCI_FFT )
60 : allocate( trigs(2*im+100) )
61 : icffta = 0
62 : rcffta = D0_0
63 : call dzfftm(0, im, icffta, rcffta, rcffta, icffta, &
64 : rcffta, icffta, trigs, rcffta, icffta)
65 : #else
66 4608 : allocate( trigs(3*im/2+1) )
67 1536 : call fftfax(im, ifax, trigs)
68 : #endif
69 :
70 1536 : end subroutine pftinit
71 :
72 :
73 : !-----------------------------------------------------------------------
74 : !BOP
75 : ! !IROUTINE: pft2d --- Two-dimensional fast Fourier transform
76 : !
77 : ! !INTERFACE:
78 3096576 : subroutine pft2d(p, s, damp, im, jp, q1, q2)
79 :
80 : ! !USES:
81 : implicit none
82 :
83 : ! !INPUT PARAMETERS:
84 : integer im ! Total X dimension
85 : integer jp ! Total Y dimension
86 : real(r8) s(jp) ! 3-point algebraic filter
87 : real(r8) damp(im,jp) ! FFT damping coefficients
88 :
89 : ! !INPUT/OUTPUT PARAMETERS:
90 : real(r8) q1( im+2, *) ! Work array
91 : real(r8) q2(*) ! Work array
92 : real(r8) p(im,jp) ! Array to be polar filtered
93 :
94 : ! !DESCRIPTION:
95 : !
96 : ! Perform a two-dimensional fast Fourier transformation.
97 : !
98 : ! !REVISION HISTORY:
99 : ! 01.01.30 Lin Put into this module
100 : ! 01.03.26 Sawyer Added ProTeX documentation
101 : ! 02.04.05 Sawyer Integrated newest FVGCM version
102 : ! 05.05.17 Sawyer Merged CAM and GEOS-5
103 : ! 05.07.26 Worley Removed ifax, trigs from arg list
104 : !
105 : !EOP
106 : !-----------------------------------------------------------------------
107 : !BOC
108 : !
109 : ! !LOCAL VARIABLES:
110 : real(r8) rsc, bt
111 : integer i, j, n, nj
112 :
113 : !Local Auto arrays:
114 6193152 : real(r8) ptmp(0:im+1)
115 : !!! real(r8) q1( im+2, jp)
116 : !!! real(r8) q2( (im+1)*jp )
117 6193152 : integer jf(jp)
118 :
119 3096576 : nj = 0
120 :
121 16353792 : do 200 j=1,jp
122 :
123 13257216 : if(s(j) > D1_01 ) then
124 7133952 : if(fft_flt .eq. 0 .and. s(j) <= D4_0) then
125 :
126 0 : rsc = D1_0/s(j)
127 0 : bt = D0_5*(s(j)-D1_0)
128 :
129 0 : do i=1,im
130 0 : ptmp(i) = p(i,j)
131 : enddo
132 0 : ptmp( 0) = p(im,j)
133 0 : ptmp(im+1) = p(1 ,j)
134 :
135 0 : do i=1,im
136 0 : p(i,j) = rsc * ( ptmp(i) + bt*(ptmp(i-1)+ptmp(i+1)) )
137 : enddo
138 :
139 : else
140 :
141 : ! Packing for FFT
142 7133952 : nj = nj + 1
143 7133952 : jf(nj) = j
144 :
145 2061712128 : do i=1,im
146 2061712128 : q1(i,nj) = p(i,j)
147 : enddo
148 7133952 : q1(im+1,nj) = D0_0
149 7133952 : q1(im+2,nj) = D0_0
150 :
151 : endif
152 : endif
153 3096576 : 200 continue
154 :
155 3096576 : if( nj == 0) return
156 :
157 1763328 : call fftrans(damp, im, jp, nj, jf, q1, q2)
158 :
159 8897280 : do n=1,nj
160 2063475456 : do i=1,im
161 2061712128 : p(i,jf(n)) = q1(i,n)
162 : enddo
163 : enddo
164 :
165 : return
166 : !EOC
167 : end subroutine pft2d
168 : !-----------------------------------------------------------------------
169 :
170 : !-----------------------------------------------------------------------
171 : !BOP
172 : ! !IROUTINE: fftrans --- Two-dimensional fast Fourier transform
173 : !
174 : ! !INTERFACE:
175 1763328 : subroutine fftrans(damp, im, jp, nj, jf, q1, q2)
176 :
177 : ! !USES:
178 : implicit none
179 :
180 : ! !INPUT PARAMETERS:
181 : integer im ! Total X dimension
182 : integer jp ! Total Y dimension
183 : integer nj ! Number of transforms
184 : integer jf(jp) ! J index versus transform number
185 : real(r8) damp(im,jp) ! FFT damping coefficients
186 :
187 : ! !INPUT/OUTPUT PARAMETERS:
188 : real(r8) q1( im+2, *) ! Work array
189 : real(r8) q2(*) ! Work array
190 :
191 : ! !DESCRIPTION:
192 : !
193 : ! Perform a two-dimensional fast Fourier transformation.
194 : !
195 : ! !REVISION HISTORY:
196 : ! 05.05.15 Mirin Initial combined version
197 : !
198 : !EOP
199 : !-----------------------------------------------------------------------
200 : !BOC
201 : !
202 : ! !LOCAL VARIABLES:
203 : integer i, n
204 : real (r8) ooim
205 :
206 : !Local Auto arrays:
207 :
208 : #if defined( LIBSCI_FFT )
209 : real (r8) qwk(2*im+4, jp)
210 : complex(r8) cqf(im/2+1, jp)
211 : integer imo2p
212 : #elif defined( SGI_FFT )
213 : integer*4 im_4, nj_4, imp2_4
214 : #endif
215 :
216 : #if defined( LIBSCI_FFT )
217 : imo2p = im/2 + 1
218 : ooim = D1_0/real(im,r8)
219 :
220 : call dzfftm(-1, im, nj, D1_0, q1, im+2, cqf, imo2p, &
221 : trigs, qwk, 0)
222 :
223 : do n=1,nj
224 : do i=3,imo2p
225 : cqf(i,n) = cqf(i,n) * damp(2*i-2,jf(n))
226 : enddo
227 : enddo
228 :
229 : call zdfftm( 1, im, nj, ooim, cqf, imo2p, q1, im+2, &
230 : trigs, qwk, 0)
231 : #elif defined( SGI_FFT )
232 : im_4 = im
233 : nj_4 = nj
234 : imp2_4 = im+2
235 : call dzfftm1du (-1, im_4, nj_4, q1, 1, imp2_4, trigs)
236 : do n=1,nj
237 : do i=5,im+2
238 : q1(i,n) = q1(i,n) * damp(i-2,jf(n))
239 : enddo
240 : enddo
241 : call dzfftm1du (1, im_4, nj_4, q1, 1, imp2_4, trigs)
242 : ooim = D1_0/real(im,r8)
243 : do n=1,nj
244 : do i=1,im+2
245 : q1(i,n) = ooim*q1(i,n)
246 : enddo
247 : enddo
248 : #else
249 1763328 : call fft991 (q1, q2, trigs, ifax, 1, im+2, im, nj, -1)
250 8897280 : do n=1,nj
251 2049207552 : do i=5,im+2
252 2047444224 : q1(i,n) = q1(i,n) * damp(i-2,jf(n))
253 : enddo
254 : enddo
255 1763328 : call fft991 (q1, q2, trigs, ifax, 1, im+2, im, nj, 1)
256 : #endif
257 :
258 1763328 : return
259 : !EOC
260 : end subroutine fftrans
261 : !-----------------------------------------------------------------------
262 :
263 : !-----------------------------------------------------------------------
264 : !BOP
265 : ! !IROUTINE: pft_cf --- Calculate algebraic and FFT polar filters
266 : !
267 : ! !INTERFACE:
268 3072 : subroutine pft_cf(im, jm, js2g0, jn2g0, jn1g1, sc, se, dc, de, &
269 3072 : cosp, cose, ycrit)
270 :
271 : ! !USES:
272 : implicit none
273 :
274 : ! !INPUT PARAMETERS:
275 : integer im ! Total X dimension
276 : integer jm ! Total Y dimension
277 : integer js2g0 ! j south limit ghosted 0 (SP: from 2)
278 : integer jn2g0 ! j north limit ghosted 0 (NP: from jm-1)
279 : integer jn1g1 ! j north limit ghosted 1 (starts jm)
280 : real (r8) cosp(jm) ! cosine array
281 : real (r8) cose(jm) ! cosine array
282 : real (r8) ycrit ! critical value
283 :
284 : ! !OUTPUT PARAMETERS:
285 : real (r8) sc(js2g0:jn2g0) ! Algebric filter at center
286 : real (r8) se(js2g0:jn1g1) ! Algebric filter at edge
287 : real (r8) dc(im,js2g0:jn2g0) ! FFT filter at center
288 : real (r8) de(im,js2g0:jn1g1) ! FFT filter at edge
289 :
290 : ! !DESCRIPTION:
291 : !
292 : ! Compute coefficients for the 3-point algebraic and the FFT
293 : ! polar filters.
294 : !
295 : ! !REVISION HISTORY:
296 : !
297 : ! 99.01.01 Lin Creation
298 : ! 99.08.20 Sawyer/Lin Changes for SPMD mode
299 : ! 01.01.30 Lin Put into this module
300 : ! 01.03.26 Sawyer Added ProTeX documentation
301 : !
302 : !EOP
303 : !-----------------------------------------------------------------------
304 : !BOC
305 : !
306 : ! !LOCAL VARIABLES:
307 : real (r8), parameter :: pi = 3.14159265358979323846_R8
308 : integer i, j
309 : real (r8) dl, coszc, cutoff, phi, damp
310 :
311 3072 : coszc = cos(ycrit*pi/D180_0)
312 :
313 : ! INIT fft polar coefficients:
314 3072 : dl = pi/real(im,r8)
315 3072 : cutoff = D1EM20
316 :
317 21216 : do j=js2g0,jn2g0
318 5246688 : do i=1,im
319 5243616 : dc(i,j) = D1_0
320 : enddo
321 : enddo
322 :
323 22800 : do j=js2g0,jn1g1
324 5704464 : do i=1,im
325 5701392 : de(i,j) = D1_0
326 : enddo
327 : enddo
328 :
329 : ! write(iulog,*) '3-point polar filter coefficients:'
330 :
331 : !************
332 : ! Cell center
333 : !************
334 21216 : do j=js2g0,jn2g0
335 18144 : sc(j) = (coszc/cosp(j))**2
336 :
337 21216 : if(sc(j) > D1_0 ) then
338 9696 : if(fft_flt .eq. 0 .and. sc(j) <= D2_0) then
339 0 : sc(j) = D1_0 + (sc(j)-D1_0)/(sc(j)+D1_0)
340 9696 : elseif(fft_flt .eq. 0 .and. sc(j) <= D4_0) then
341 0 : sc(j) = D1_0 + sc(j)/(D8_0-sc(j))
342 : else
343 :
344 : ! FFT filter
345 1405920 : do i=1,im/2
346 1396224 : phi = dl * i
347 1396224 : damp = min((cosp(j)/coszc)/sin(phi),D1_0)**2
348 1396224 : if(damp < cutoff) damp = D0_0
349 1396224 : dc(2*i-1,j) = damp
350 1405920 : dc(2*i ,j) = damp
351 : enddo
352 :
353 : endif
354 : endif
355 : enddo
356 :
357 : !************
358 : ! Cell edges
359 : !************
360 22800 : do j=js2g0,jn1g1
361 19728 : se(j) = (coszc/cose(j))**2
362 :
363 22800 : if(se(j) > D1_0 ) then
364 10680 : if(fft_flt .eq. 0 .and. se(j) <= D2_0) then
365 0 : se(j) = D1_0 + (se(j)-D1_0)/(se(j)+D1_0)
366 10680 : elseif(fft_flt .eq. 0 .and. se(j) <= D4_0) then
367 0 : se(j) = D1_0 + se(j)/(D8_0-se(j))
368 : else
369 : ! FFT
370 1548600 : do i=1,im/2
371 1537920 : phi = dl * i
372 1537920 : damp = min((cose(j)/coszc)/sin(phi), D1_0)**2
373 1537920 : if(damp < cutoff) damp = D0_0
374 1537920 : de(2*i-1,j) = damp
375 1548600 : de(2*i ,j) = damp
376 : enddo
377 : endif
378 : endif
379 : enddo
380 3072 : return
381 : !EOC
382 : end subroutine pft_cf
383 : !-----------------------------------------------------------------------
384 :
385 :
386 : !-----------------------------------------------------------------------
387 : !BOP
388 : ! !IROUTINE: fftfax --- Initialize FFT
389 : !
390 : ! !INTERFACE:
391 1536 : subroutine fftfax (n, ifaxx, trigss)
392 :
393 : ! !USES:
394 : implicit none
395 :
396 : ! !DESCRIPTION:
397 : !
398 : ! Initialize the fast Fourier transform. If CPP token SGI_FFT is
399 : ! set, SGI libraries will be used. Otherwise the Fortran code
400 : ! is inlined.
401 : !
402 : ! !REVISION HISTORY:
403 : !
404 : ! 99.11.24 Sawyer Added wrappers for SGI
405 : ! 01.03.26 Sawyer Added ProTeX documentation
406 : ! 05.07.26 Worley Modified version for Cray X1
407 : !
408 : !EOP
409 : !-----------------------------------------------------------------------
410 : !BOC
411 :
412 : integer n
413 :
414 : #if defined( SGI_FFT )
415 : real(r8) trigss(1)
416 : integer ifaxx(*)
417 : ! local
418 : integer*4 nn
419 :
420 : nn=n
421 : call dzfftm1dui (nn,trigss)
422 : #else
423 : integer ifaxx(13)
424 : real(r8) trigss(3*n/2+1)
425 1536 : call set99(trigss,ifaxx,n)
426 : #endif
427 1536 : return
428 : !EOC
429 : end subroutine fftfax
430 : !-----------------------------------------------------------------------
431 :
432 : end module pft_module
|
