| blob | 78a3e01af581da2c9b89b3c9b6b66a9faf4ecb02 |
1 !
2 SUBROUTINE ATMS_Spatial_Average(Num_Obs, NChanl, FOV, Time, BT_InOut, &
3 Scanline, Error_Status)
6 IMPLICIT NONE
8 ! Declare passed variables
9 integer(i_kind) ,intent(in ) :: Num_Obs, NChanl
10 integer(i_kind) ,intent(in ) :: Fov(num_obs)
11 real(r_kind) ,intent(in ) :: Time(Num_Obs)
12 real(r_kind) ,intent(inout) :: BT_InOut(NChanl,Num_Obs)
13 integer(i_kind) ,intent( out) :: Scanline(Num_Obs)
14 integer(i_kind) ,intent( out) :: Error_Status
15 ! Declare local parameters
16 integer(i_kind), parameter :: atms1c_h_wmosatid=224
17 integer(i_kind), parameter :: lninfile=15
18 integer(i_kind), parameter :: max_fov=96
19 integer(i_kind), parameter :: max_obs=1000000
20 real(r_kind), parameter :: scan_interval = 8.0_r_kind/3.0_r_kind
21 ! Maximum number of channels
22 integer(i_kind), parameter :: MaxChans = 22
23 ! Minimum allowed BT as a function of channel number
24 real(r_kind), parameter :: MinBT(MaxChans) = &
25 (/ 120.0_r_kind, 120.0_r_kind, 190.0_r_kind, 190.0_r_kind, &
26 200.0_r_kind, 200.0_r_kind, 200.0_r_kind, 190.0_r_kind, &
27 190.0_r_kind, 180.0_r_kind, 180.0_r_kind, 180.0_r_kind, &
28 190.0_r_kind, 200.0_r_kind, 200.0_r_kind, 120.0_r_kind, &
29 120.0_r_kind, 120.0_r_kind, 150.0_r_kind, 170.0_r_kind, &
30 180.0_r_kind, 190.0_r_kind /)
31 ! Maximum allowed BT as a function of channel number
32 real(r_kind), parameter :: MaxBT(MaxChans) = &
33 (/ 320.0_r_kind, 320.0_r_kind, 300.0_r_kind, 300.0_r_kind, &
34 300.0_r_kind, 270.0_r_kind, 250.0_r_kind, 240.0_r_kind, &
35 240.0_r_kind, 250.0_r_kind, 250.0_r_kind, 270.0_r_kind, &
36 280.0_r_kind, 290.0_r_kind, 300.0_r_kind, 320.0_r_kind, &
37 320.0_r_kind, 300.0_r_kind, 300.0_r_kind, 300.0_r_kind, &
38 300.0_r_kind, 300.0_r_kind /)
40 ! Declare local variables
41 character(30) :: Cline
42 integer(i_kind) :: i, iscan, ifov, ichan, nchannels, wmosatid, version
43 integer(i_kind) :: ios, max_scan, mintime
44 integer(i_kind) :: nxaverage(nchanl), nyaverage(nchanl)
45 integer(i_Kind) :: channelnumber(nchanl),qc_dist(nchanl)
46 integer(i_kind), ALLOCATABLE :: scanline_back(:,:)
47 real(r_kind) :: sampling_dist, beamwidth(nchanl)
48 real(r_kind) :: newwidth(nchanl), cutoff(nchanl)
49 real(r_kind), allocatable, target :: bt_image(:,:,:)
50 real(r_kind), pointer :: bt_image1(:,:)
51 Error_Status=0
52 IF (NChanl > MaxChans) THEN
53 WRITE(0,*) 'Unexpected number of ATMS channels: ',nchanl
54 Error_Status = 1
55 RETURN
56 END IF
57 ! Read the beamwidth requirements
58 OPEN(lninfile,file='radiance_info/atms_beamwidth.txt',form='formatted',status='old', &
59 iostat=ios)
60 IF (ios /= 0) THEN
61 WRITE(*,*) 'Unable to open atms_beamwidth.txt'
62 Error_Status=1
63 RETURN
64 ENDIF
65 wmosatid=999
66 read(lninfile,'(a30)',iostat=ios) Cline
67 DO WHILE (wmosatid /= atms1c_h_wmosatid .AND. ios == 0)
68 DO WHILE (Cline(1:1) == '#')
69 read(lninfile,'(a30)') Cline
70 ENDDO
71 READ(Cline,*) wmosatid
73 read(lninfile,'(a30)') Cline
74 DO WHILE (Cline(1:1) == '#')
75 read(lninfile,'(a30)') Cline
76 ENDDO
77 READ(Cline,*) version
79 read(lninfile,'(a30)') Cline
80 DO WHILE (Cline(1:1) == '#')
81 read(lninfile,'(a30)') Cline
82 ENDDO
83 READ(Cline,*) sampling_dist
85 read(lninfile,'(a30)') Cline
86 DO WHILE (Cline(1:1) == '#')
87 read(lninfile,'(a30)') Cline
88 ENDDO
89 READ(Cline,*) nchannels
91 read(lninfile,'(a30)') Cline
92 if (nchannels > 0) then
93 DO ichan=1,nchannels
94 read(lninfile,'(a30)') Cline
95 DO WHILE (Cline(1:1) == '#')
96 read(lninfile,'(a30)') Cline
97 ENDDO
98 READ(Cline,*) channelnumber(ichan),beamwidth(ichan), &
99 newwidth(ichan),cutoff(ichan),nxaverage(ichan), &
100 nyaverage(ichan), qc_dist(ichan)
101 ENDDO
102 end if
103 read(lninfile,'(a30)',iostat=ios) Cline
104 ENDDO
105 IF (wmosatid /= atms1c_h_wmosatid) THEN
106 WRITE(*,*) 'ATMS_Spatial_Averaging: sat id not matched in atms_beamwidth.dat'
107 Error_Status=1
108 RETURN
109 ENDIF
110 CLOSE(lninfile)
112 ! Determine scanline from time
113 MinTime = MINVAL(Time)
114 Scanline(:) = NINT((Time(1:Num_Obs)-MinTime)/Scan_Interval)+1
115 Max_Scan=MAXVAL(Scanline)
116 write(*,*) 'Max_Scan:',Max_Scan
118 ALLOCATE(BT_Image(Max_FOV,Max_Scan,nchanl))
119 ALLOCATE(Scanline_Back(Max_FOV,Max_Scan))
120 BT_Image(:,:,:) = 1000.0_r_kind
122 ScanLine_Back(:,:) = -1
123 DO I=1,Num_Obs
124 bt_image(FOV(I),Scanline(I),:)=bt_inout(:,I)
125 Scanline_Back(FOV(I),Scanline(I))=I
126 END DO
127 DO IChan=1,nchanl
130 bt_image1 => bt_image(:,:,ichan)
131 ! Set all scan positions to missing in a scanline if one is missing
132 do iscan=1,max_scan
133 if (ANY(bt_image1(:,iscan) > 500.0_r_kind)) &
134 bt_image1(:,iscan)=1000.0_r_kind
135 enddo
136 ! If the channel number is present in the channelnumber array we should process it
137 ! (otherwise bt_inout just keeps the same value):
138 IF (ANY(channelnumber(1:nchannels) == ichan)) THEN
139 CALL MODIFY_BEAMWIDTH ( max_fov, max_scan, bt_image1, &
140 sampling_dist, beamwidth(ichan), newwidth(ichan), &
141 cutoff(ichan), nxaverage(ichan), nyaverage(ichan), &
142 qc_dist(ichan), MinBT(Ichan), MaxBT(IChan), IOS)
144 IF (IOS == 0) THEN
145 do iscan=1,max_scan
146 do ifov=1,max_fov
147 IF (Scanline_Back(IFov, IScan) > 0) &
148 bt_inout(channelnumber(ichan),Scanline_Back(IFov, IScan)) = &
149 BT_Image1(ifov,iscan)
150 end do
151 end do
152 ELSE
153 Error_Status=1
154 RETURN
155 END IF
156 END IF
157 END DO
158 DEALLOCATE(BT_Image, Scanline_Back)
159 NULLIFY(BT_Image1)
161 END Subroutine ATMS_Spatial_Average
163 SUBROUTINE MODIFY_BEAMWIDTH ( nx, ny, image, sampling_dist,&
164 beamwidth, newwidth, mtfcutoff, nxaverage, nyaverage, qc_dist, &
165 Minval, MaxVal, Error)
167 !-----------------------------------------
168 ! Name: $Id: modify_beamwidth.F 222 2010-08-11 14:39:09Z frna $
169 !
170 ! Purpose:
171 ! Manipulate the effective beam width of an image. For example, convert ATMS
172 ! to AMSU-A-like resolution while reducing the noise.
173 !
174 ! Method:
175 ! 1) Pad the image to a power of 2 in each dimension.
176 ! If FFT technique is to be used then:
177 ! 2) Assuming Gaussian beam shapes, calcluate the input and output Modulation
178 ! Transfer Functions (MTF).
179 ! 3) FFT image to frequency domain (2-D).
180 ! 4) Multiply by output MTF divided by input MTF. If a cut-off is specified
181 ! (when attempting to make the beam width narrower), attenuate further
182 ! by an exponential function - factor of 2 at the cutoff.
183 ! 5) FFT back to image domain
184 ! Finally,
185 ! 6) Over-write the input image, with averaging if requested.
186 !
187 ! COPYRIGHT
188 ! This software was developed within the context of the EUMETSAT Satellite
189 ! Application Facility on Numerical Weather Prediction (NWP SAF), under the
190 ! Cooperation Agreement dated 1 December 2006, between EUMETSAT and the
191 ! Met Office, UK, by one or more partners within the NWP SAF. The partners
192 ! in the NWP SAF are the Met Office, ECMWF, KNMI and MeteoFrance.
193 !
194 ! Copyright 2010, EUMETSAT, All Rights Reserved.
195 !
196 ! History:
197 ! Version Date Comment
198 !
199 ! 1.0 22/07/2010 N.C.Atkinson
200 ! 1.1 21/11/2011 Convert to f90. A. Collard
201 !
202 ! Code Description:
203 ! FORTRAN 77, following AAPP standards
204 !
205 ! Declarations:
206 IMPLICIT NONE
207 ! Parameters
208 real(r_kind), parameter :: Missing_Value=-888888.0
209 INTEGER(I_KIND), PARAMETER :: nxmax=128 !Max number of spots per scan line
210 INTEGER(I_KIND), PARAMETER :: nymax=8192 !Max number of lines. Allows 6hrs of ATMS.
211 ! Arguments
212 INTEGER(I_KIND), INTENT(IN) :: nx, ny !Size of image
213 REAL(R_KIND), INTENT(INOUT) :: image(nx,ny) !BT or radiance image
214 REAL(R_KIND), INTENT(IN) :: sampling_dist !typically degrees
215 REAL(R_KIND), INTENT(IN) :: beamwidth !ditto
216 REAL(R_KIND), INTENT(IN) :: newwidth !ditto
217 REAL(R_KIND), INTENT(IN) :: mtfcutoff !0.0 to 1.0
218 INTEGER(I_KIND), INTENT(IN) :: nxaverage !Number of samples to average (or zero)
219 INTEGER(I_KIND), INTENT(IN) :: nyaverage !Number of samples to average (or zero)
220 INTEGER(I_KIND), INTENT(IN) :: qc_dist !Number of samples around missing data to set to
221 REAL(R_KIND), INTENT(IN) :: maxval !BTs above this are considered missing
222 REAL(R_KIND), INTENT(IN) :: minval !BTs below this are considered missing
223 INTEGER(I_KIND), INTENT(OUT) :: Error !Error Status
225 ! Local variables
226 INTEGER(I_KIND) :: nxpad, nypad, dx, dy
227 INTEGER(I_KIND) :: i,j,k,ix,iy, ii, jj
228 INTEGER(I_KIND) :: ifirst
229 INTEGER(I_KIND) :: xpow2, ypow2
230 INTEGER(I_KIND) :: nxav2, nyav2, naverage
231 INTEGER(I_KIND) :: deltax, minii, maxii, minjj, maxjj
232 REAL(R_KIND), ALLOCATABLE :: mtfxin(:),mtfxout(:)
233 REAL(R_KIND), ALLOCATABLE :: mtfyin(:),mtfyout(:)
234 REAL(R_KIND) :: mtfin,mtfout,mtf_constant
235 REAL(R_KIND), ALLOCATABLE :: mtfpad(:,:)
236 REAL(R_KIND), ALLOCATABLE :: imagepad(:,:)
237 REAL(R_KIND), ALLOCATABLE :: work(:)
238 REAL(R_KIND) :: f,df,factor
239 REAL(R_KIND) :: PI, LN2, LNcsquared
240 LOGICAL :: missing
241 LOGICAL, ALLOCATABLE :: gooddata_map(:,:)
242 ! End of declarations
243 !-----------------------------------------
245 PI = 4.0_r_kind*atan(1.0)
246 LN2 = LOG(2.0_r_kind)
247 MTF_Constant=-(PI/(2*sampling_dist))**2/LN2
248 IF (mtfcutoff > 0.0_r_kind) LNcsquared = LOG(mtfcutoff)**2
249 nxav2 = nxaverage/2
250 nyav2 = nyaverage/2
251 naverage = nxaverage*nyaverage
252 Error = 0
253 !1) Pad the image up to the nearest power of 2 in each dimension, by reversing
254 !the points near the edge.
255 xpow2 = INT(LOG(nx*1.0_r_kind)/LN2 + 1.0_r_kind)
256 ypow2 = INT(LOG(ny*1.0_r_kind)/LN2 + 1.0_r_kind)
257 nxpad = 2**xpow2
258 nypad = 2**ypow2
259 dx = (nxpad - nx)/2
260 dy = (nypad - ny)/2
261 IF (nxpad > nxmax) THEN
262 write(*,*) 'ATMS_Spatial_Average: nx too large, maximum allowed value is ',nxmax-1
263 Error = 1
264 RETURN
265 END IF
267 IF (nypad > nymax) THEN
268 write(*,*) 'ATMS_Spatial_Average: ny too large, maximum allowed value is ',nymax-1
269 Error = 1
270 RETURN
271 END IF
272 ALLOCATE(mtfxin(nxpad),mtfxout(nxpad))
273 ALLOCATE(mtfyin(nypad),mtfyout(nypad))
274 ALLOCATE(mtfpad(nxpad,nypad))
275 ALLOCATE(imagepad(nxpad,nypad))
276 ALLOCATE(work(nypad))
277 ALLOCATE(gooddata_map(nxmax,nymax))
278 !Loop over scan positions
279 DO j=dy+1,dy+ny
280 DO i=dx+1,dx+nx
281 if (image(i-dx,j-dy) < minval) &
282 image(i-dx,j-dy) = minval - 1.0_r_kind
283 if (image(i-dx,j-dy) > maxval ) &
284 image(i-dx,j-dy) = maxval + 1.0_r_kind
285 imagepad(i,j) = image(i-dx,j-dy) !Take a copy of the input data
286 gooddata_map(i,j) = .TRUE. ! Initialised for step 6)
287 ENDDO
288 !Interpolate missing points in the along-track direction
289 ifirst = -1
290 missing = .false.
292 DO i=dx+1,dx+nx
293 IF (.not.missing) THEN
294 IF (imagepad(i,j) >= minval .AND. imagepad(i,j) <= maxval) THEN
295 ifirst = i
296 ELSE
297 missing = .true.
298 ENDIF
299 ELSE
300 IF (imagepad(i,j) >= minval .AND. imagepad(i,j) <= maxval) THEN !First good point
301 ! after missing
302 missing = .false.
303 IF (ifirst == -1) THEN
304 DO k=dx+1,i-1
305 imagepad(k,j) = imagepad(i,j) !Constant
306 ENDDO
307 ELSE
308 DO k=ifirst+1,i-1
309 factor = (i-k)*1.0_r_kind/(i-ifirst) !Interpolate
310 imagepad(k,j) = imagepad(ifirst,j)*factor + &
311 imagepad(i,j)*(1.0_r_kind-factor)
312 ENDDO
313 ENDIF
314 ENDIF
315 ENDIF
316 ENDDO
317 IF (missing) THEN !Last scan is missing
318 IF (ifirst >= 1) then
319 DO k=ifirst+1,dx+nx
320 imagepad(k,j) = imagepad(ifirst,j) !Constant
321 ENDDO
322 ENDIF
323 ENDIF
324 !Continue padding the edges
325 DO i=1,dx
326 imagepad(i,j) = imagepad(dx+dx+2-i,j)
327 ENDDO
328 DO i=nx+dx+1,nxpad
329 imagepad(i,j) = imagepad(nx+dx+nx+dx-i,j)
330 ENDDO
331 ENDDO
332 DO j=1,dy
333 DO i=1,nxpad
334 imagepad(i,j) = imagepad(i,dy+dy+2-j)
335 ENDDO
336 ENDDO
337 DO j=ny+dy+1,nypad
338 DO i=1,nxpad
339 imagepad(i,j) = imagepad(i,ny+dy+ny+dy-j)
340 ENDDO
341 ENDDO
342 !2) Compute the MTF modifications. Assume beams are Gaussian.
343 IF (newwidth > 0) THEN
344 df = 1.0_r_kind/nxpad
345 DO i=1,nxpad/2+1
346 f = df*(i-1) !DC to Nyquist
347 mtfxin(i) = exp(MTF_Constant*(f*beamwidth)**2)
348 mtfxout(i) = exp(MTF_Constant*(f*newwidth)**2)
349 IF (i > 1 .AND. i < nxpad/2+1) THEN
350 mtfxin(nxpad-i+2) = mtfxin(i)
351 mtfxout(nxpad-i+2) = mtfxout(i)
352 ENDIF
353 ENDDO
354 df = 1.0_r_kind/nypad
355 DO i=1,nypad/2+1
356 f = df*(i-1) !DC to Nyquist
357 mtfyin(i) = exp(MTF_Constant*(f*beamwidth)**2)
358 mtfyout(i) = exp(MTF_Constant*(f*newwidth)**2)
359 IF (i > 1 .AND. i < nypad/2+1) THEN
360 mtfyin(nypad-i+2) = mtfyin(i)
361 mtfyout(nypad-i+2) = mtfyout(i)
362 ENDIF
363 ENDDO
364 DO i=1,nxpad
365 DO j=1,nypad
366 mtfin = mtfxin(i)*mtfyin(j)
367 mtfout = mtfxout(i)*mtfyout(j)
368 if (mtfcutoff > 0.0_r_kind) THEN
369 mtfpad(i,j) = (mtfout * &
370 exp(-LN2/LNcsquared*(LOG(mtfout))**2))/mtfin
371 else
372 mtfpad(i,j) = mtfout/mtfin
373 endif
374 ENDDO
375 ENDDO
376 !3) Fourier transform, line by line then column by column.
377 !After each FFT, points 1 to nxpad/2+1 contain the real part of the spectrum,
378 !the rest contain the imaginary part in reverse order.
379 DO j=1,nypad
380 CALL SFFTCF(imagepad(:,j),nxpad,xpow2)
381 ENDDO
382 DO i=1,nxpad
383 DO j=1,nypad
384 work(j) = imagepad(i,j)
385 ENDDO
386 CALL SFFTCF(work,nypad,ypow2)
387 DO j=1,nypad
388 imagepad(i,j) = work(j)
389 ENDDO
390 ENDDO
391 !4) Multiply the spectrum by the MTF factor
392 DO j=1,nypad
393 DO i=1,nxpad
394 imagepad(i,j) = imagepad(i,j)*mtfpad(i,j)
395 ENDDO
396 ENDDO
397 !5) Inverse Fourier transform, column by column then line by line
398 DO i=1,nxpad
399 DO j=1,nypad
400 work(j) = imagepad(i,j)
401 ENDDO
402 CALL SFFTCB(work,nypad,ypow2)
403 DO j=1,nypad
404 imagepad(i,j) = work(j)
405 ENDDO
406 ENDDO
407 DO j=1,nypad
408 CALL SFFTCB(imagepad(:,j),nxpad,xpow2)
409 ENDDO
410 ENDIF !New width is specified
411 !6) Reset missing values in gooddata_map, based on qc_dist and the values
412 ! in the input image array
413 ! Set the ends of the image to missing in the along track direction
414 ! (doing the same across track will remove too much data)
415 gooddata_map(:,1:qc_dist)=.FALSE.
416 gooddata_map(:,ny-qc_dist+1:ny)=.FALSE.
418 DO j=1,ny
419 DO i=1,nx
420 IF (image(i,j) <= minval .OR. image(i,j) >= maxval ) THEN
421 minjj=max(j+dy-qc_dist,1)
422 maxjj=min(j+dy+qc_dist,nymax)
423 DO jj=minjj,maxjj
424 deltax=INT(SQRT(REAL(qc_dist**2 - (jj-j-dy)**2 )))
425 minii=max(i+dx-deltax,1)
426 maxii=min(i+dx+deltax,nxmax)
427 DO ii=minii,maxii
428 gooddata_map(ii,jj)=.FALSE.
429 END DO
430 END DO
431 END IF
432 END DO
433 END DO
434 !7) Over-write the input image (points that are not missing)
435 DO j=1,ny
436 DO i=1,nx
437 IF (gooddata_map(i+dx,j+dy)) THEN
438 IF (nxav2 == 0.0_r_kind .AND. nyav2 == 0) THEN
439 image(i,j) = imagepad(i+dx,j+dy)
440 ELSE
441 image(i,j) = 0.0_r_kind !Do averaging
442 DO ix = -nxav2,nxav2
443 DO iy = -nyav2,nyav2
444 image(i,j) = image(i,j) + imagepad(i+dx+ix,j+dy+iy)
445 ENDDO
446 ENDDO
447 image(i,j) = image(i,j)/naverage
448 ENDIF
449 ELSE
450 image(i,j) = missing_value
451 END IF
452 ENDDO
453 ENDDO
454 !8) Deallocate arrays
455 DEALLOCATE(mtfxin,mtfxout)
456 DEALLOCATE(mtfyin,mtfyout)
457 DEALLOCATE(mtfpad)
458 DEALLOCATE(imagepad)
459 DEALLOCATE(work)
460 DEALLOCATE(gooddata_map)
461 RETURN
462 END SUBROUTINE MODIFY_BEAMWIDTH
463 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
464 !
465 ! A real-valued, in place, split-radix FFT program
466 ! Real input and output in data array X
467 ! Length is N = 2 ** M
468 ! Decimation-in-time, cos/sin in second loop
469 ! Output in order:
470 ! [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
471 !
472 ! This FFT computes
473 ! X(k) = sum_{j=0}^{N-1} x(j)*exp(-2ijk*pi/N)
474 !
475 !
476 ! H.V. Sorensen, Rice University, Oct. 1985
477 !
478 ! Reference: H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus;
479 ! Real-Valued Fast Fourier Transform Algorithms; IEEE
480 ! Trans. Acoust., Speech, Signal Process.; Vol ASSP-35,
481 ! June 1987, pp. 849-863.
482 !
483 ! This code was originally named RVFFT.
484 !
485 ! History:
486 ! 21/11/2011 Converted to something resembling f90. A.Collard
487 !
488 !
489 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
490 !
491 SUBROUTINE SFFTCF( X, N, M )
492 IMPLICIT NONE
493 ! ... Parameters ...
494 REAL(R_KIND), PARAMETER :: SQRT2 = 1.4142135623730950488
495 REAL(R_KIND), PARAMETER :: TWOPI = 6.2831853071795864769
496 ! ... Scalar arguments ...
497 INTEGER(I_KIND), INTENT(IN) :: N, M
498 ! ... Array arguments ...
499 REAL(R_KIND), INTENT(INOUT) :: X(N)
500 ! ... Local scalars ...
501 INTEGER(I_KIND) J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8
502 INTEGER(I_KIND) N1, N2, N4, N8
503 REAL(R_KIND) XT, R1, T1, T2, T3, T4, T5, T6
504 REAL(R_KIND) A, A3, E, CC1, SS1, CC3, SS3
505 !
506 ! ... Exe. statements ...
507 !
508 IF ( N == 1 ) RETURN
509 !
510 100 J = 1
511 N1 = N - 1
512 DO 104, I = 1, N1
513 IF ( I >= J ) GOTO 101
514 XT = X(J)
515 X(J) = X(I)
516 X(I) = XT
517 101 K = N / 2
518 102 IF ( K >= J ) GOTO 103
519 J = J - K
520 K = K / 2
521 GOTO 102
522 103 J = J + K
523 104 CONTINUE
524 !
525 IS = 1
526 ID = 4
527 70 DO 60, I0 = IS, N, ID
528 I1 = I0 + 1
529 R1 = X(I0)
530 X(I0) = R1 + X(I1)
531 X(I1) = R1 - X(I1)
532 60 CONTINUE
533 IS = 2 * ID - 1
534 ID = 4 * ID
535 IF ( IS < N ) GOTO 70
536 !
537 N2 = 2
538 DO 10, K = 2, M
539 N2 = N2 * 2
540 N4 = N2 / 4
541 N8 = N2 / 8
542 E = TWOPI / N2
543 IS = 0
544 ID = N2 * 2
545 40 DO 38, I = IS, N-1, ID
546 I1 = I + 1
547 I2 = I1 + N4
548 I3 = I2 + N4
549 I4 = I3 + N4
550 T1 = X(I4) + X(I3)
551 X(I4) = X(I4) - X(I3)
552 X(I3) = X(I1) - T1
553 X(I1) = X(I1) + T1
554 IF ( N4 == 1 ) GOTO 38
555 I1 = I1 + N8
556 I2 = I2 + N8
557 I3 = I3 + N8
558 I4 = I4 + N8
559 T1 = ( X(I3) + X(I4) ) / SQRT2
560 T2 = ( X(I3) - X(I4) ) / SQRT2
561 X(I4) = X(I2) - T1
562 X(I3) = - X(I2) - T1
563 X(I2) = X(I1) - T2
564 X(I1) = X(I1) + T2
565 38 CONTINUE
566 IS = 2 * ID - N2
567 ID = 4 * ID
568 IF ( IS < N ) GOTO 40
569 A = E
570 DO 32, J = 2, N8
571 A3 = 3 * A
572 CC1 = COS(A)
573 SS1 = SIN(A)
574 CC3 = COS(A3)
575 SS3 = SIN(A3)
576 A = J * E
577 IS = 0
578 ID = 2 * N2
579 36 DO 30, I = IS, N-1, ID
580 I1 = I + J
581 I2 = I1 + N4
582 I3 = I2 + N4
583 I4 = I3 + N4
584 I5 = I + N4 - J + 2
585 I6 = I5 + N4
586 I7 = I6 + N4
587 I8 = I7 + N4
588 T1 = X(I3) * CC1 + X(I7) * SS1
589 T2 = X(I7) * CC1 - X(I3) * SS1
590 T3 = X(I4) * CC3 + X(I8) * SS3
591 T4 = X(I8) * CC3 - X(I4) * SS3
592 T5 = T1 + T3
593 T6 = T2 + T4
594 T3 = T1 - T3
595 T4 = T2 - T4
596 T2 = X(I6) + T6
597 X(I3) = T6 - X(I6)
598 X(I8) = T2
599 T2 = X(I2) - T3
600 X(I7) = - X(I2) - T3
601 X(I4) = T2
602 T1 = X(I1) + T5
603 X(I6) = X(I1) - T5
604 X(I1) = T1
605 T1 = X(I5) + T4
606 X(I5) = X(I5) - T4
607 X(I2) = T1
608 30 CONTINUE
609 IS = 2 * ID - N2
610 ID = 4 * ID
611 IF ( IS < N ) GOTO 36
612 32 CONTINUE
613 10 CONTINUE
614 RETURN
615 !
616 ! ... End of subroutine SFFTCF ...
617 !
618 END SUBROUTINE SFFTCF
619 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
620 !
621 ! A real-valued, in place, split-radix IFFT program
622 ! Hermitian symmetric input and real output in array X
623 ! Length is N = 2 ** M
624 ! Decimation-in-frequency, cos/sin in second loop
625 ! Input order:
626 ! [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
627 !
628 ! This FFT computes
629 ! x(j) = (1/N) * sum_{k=0}^{N-1} X(k)*exp(2ijk*pi/N)
630 !
631 !
632 ! H.V. Sorensen, Rice University, Nov. 1985
633 !
634 ! Reference: H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus;
635 ! Real-Valued Fast Fourier Transform Algorithms; IEEE
636 ! Trans. Acoust., Speech, Signal Process.; Vol ASSP-35,
637 ! June 1987, pp. 849-863.
638 !
639 ! This code was originally named IRVFFT.
640 !
641 ! History:
642 ! 21/11/2011 Converted to something resembling f90. A.Collard
643 !
644 !
645 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
646 !
647 SUBROUTINE SFFTCB( X, N, M )
648 IMPLICIT NONE
649 ! ... Parameters ...
650 REAL(R_KIND), PARAMETER :: SQRT2 = 1.4142135623730950488
651 REAL(R_KIND), PARAMETER :: TWOPI = 6.2831853071795864769
652 ! ... Scalar arguments ...
653 INTEGER(I_KIND), INTENT(IN) :: N, M
654 ! ... Array arguments ...
655 REAL(R_KIND), INTENT(INOUT) :: X(N)
656 ! ... Local scalars ...
657 INTEGER(I_KIND) J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8
658 INTEGER(I_KIND) N1, N2, N4, N8
659 REAL(R_KIND) XT, R1, T1, T2, T3, T4, T5
660 REAL(R_KIND) A, A3, E, CC1, SS1, CC3, SS3
661 !
662 ! ... Exe. statements ...
663 !
664 IF ( N == 1 ) RETURN
665 !
666 N2 = 2 * N
667 DO 10, K = 1, M-1
668 IS = 0
669 ID = N2
670 N2 = N2 / 2
671 N4 = N2 / 4
672 N8 = N4 / 2
673 E = TWOPI / N2
674 17 DO 15, I = IS, N-1, ID
675 I1 = I + 1
676 I2 = I1 + N4
677 I3 = I2 + N4
678 I4 = I3 + N4
679 T1 = X(I1) - X(I3)
680 X(I1) = X(I1) + X(I3)
681 X(I2) = 2 * X(I2)
682 X(I3) = T1 - 2 * X(I4)
683 X(I4) = T1 + 2 * X(I4)
684 IF ( N4 == 1 ) GOTO 15
685 I1 = I1 + N8
686 I2 = I2 + N8
687 I3 = I3 + N8
688 I4 = I4 + N8
689 T1 = ( X(I2) - X(I1) ) / SQRT2
690 T2 = ( X(I4) + X(I3) ) / SQRT2
691 X(I1) = X(I1) + X(I2)
692 X(I2) = X(I4) - X(I3)
693 X(I3) = 2 * ( - T2 - T1 )
694 X(I4) = 2 * ( -T2 + T1 )
695 15 CONTINUE
696 IS = 2 * ID - N2
697 ID = 4 * ID
698 IF ( IS < N-1 ) GOTO 17
699 A = E
700 DO 20, J = 2, N8
701 A3 = 3 * A
702 CC1 = COS(A)
703 SS1 = SIN(A)
704 CC3 = COS(A3)
705 SS3 = SIN(A3)
706 A = J * E
707 IS = 0
708 ID = 2 * N2
709 40 DO 30, I = IS, N-1, ID
710 I1 = I + J
711 I2 = I1 + N4
712 I3 = I2 + N4
713 I4 = I3 + N4
714 I5 = I + N4 - J + 2
715 I6 = I5 + N4
716 I7 = I6 + N4
717 I8 = I7 + N4
718 T1 = X(I1) - X(I6)
719 X(I1) = X(I1) + X(I6)
720 T2 = X(I5) - X(I2)
721 X(I5) = X(I2) + X(I5)
722 T3 = X(I8) + X(I3)
723 X(I6) = X(I8) - X(I3)
724 T4 = X(I4) + X(I7)
725 X(I2) = X(I4) - X(I7)
726 T5 = T1 - T4
727 T1 = T1 + T4
728 T4 = T2 - T3
729 T2 = T2 + T3
730 X(I3) = T5 * CC1 + T4 * SS1
731 X(I7) = - T4 * CC1 + T5 * SS1
732 X(I4) = T1 * CC3 - T2 * SS3
733 X(I8) = T2 * CC3 + T1 * SS3
734 30 CONTINUE
735 IS = 2 * ID - N2
736 ID = 4 * ID
737 IF ( IS < N-1 ) GOTO 40
738 20 CONTINUE
739 10 CONTINUE
740 !
741 IS = 1
742 ID = 4
743 70 DO 60, I0 = IS, N, ID
744 I1 = I0 + 1
745 R1 = X(I0)
746 X(I0) = R1 + X(I1)
747 X(I1) = R1 - X(I1)
748 60 CONTINUE
749 IS = 2 * ID - 1
750 ID = 4 * ID
751 IF ( IS < N ) GOTO 70
752 !
753 100 J = 1
754 N1 = N - 1
755 DO 104, I = 1, N1
756 IF ( I >= J ) GOTO 101
757 XT = X(J)
758 X(J) = X(I)
759 X(I) = XT
760 101 K = N / 2
761 102 IF ( K >= J ) GOTO 103
762 J = J - K
763 K = K / 2
764 GOTO 102
765 103 J = J + K
766 104 CONTINUE
767 XT = 1.0_r_kind / FLOAT( N )
768 DO 99, I = 1, N
769 X(I) = XT * X(I)
770 99 CONTINUE
771 RETURN
772 !
773 ! ... End of subroutine SFFTCB ...
774 !
775 END SUBROUTINE SFFTCB