wrfv2_fire/external/fftpack/fftpack5/costf1.F

   1 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   2 !
   3 !   FFTPACK 5.0
   4 !   Copyright (C) 1995-2004, Scientific Computing Division,
   5 !   University Corporation for Atmospheric Research
   6 !   Licensed under the GNU General Public License (GPL)
   7 !
   8 !   Authors:  Paul N. Swarztrauber and Richard A. Valent
   9 !
  10 !   $Id: costf1.f,v 1.2 2004/06/15 21:14:57 rodney Exp $
  11 !
  12 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  13
  14       SUBROUTINE COSTF1(N,INC,X,WSAVE,WORK,IER)
  15       REAL       X(INC,*)       ,WSAVE(*)
  16       DOUBLE PRECISION           DSUM
  17       IER = 0
  18       NM1 = N-1
  19       NP1 = N+1
  20       NS2 = N/2
  21       IF (N-2) 200,101,102
  22   101 X1H = X(1,1)+X(1,2)
  23       X(1,2) = .5*(X(1,1)-X(1,2))
  24       X(1,1) = .5*X1H
  25       GO TO 200
  26   102 IF (N .GT. 3) GO TO 103
  27       X1P3 = X(1,1)+X(1,3)
  28       TX2 = X(1,2)+X(1,2)
  29       X(1,2) = .5*(X(1,1)-X(1,3))
  30       X(1,1) = .25*(X1P3+TX2)
  31       X(1,3) = .25*(X1P3-TX2)
  32       GO TO 200
  33   103 DSUM = X(1,1)-X(1,N)
  34       X(1,1) = X(1,1)+X(1,N)
  35       DO 104 K=2,NS2
  36          KC = NP1-K
  37          T1 = X(1,K)+X(1,KC)
  38          T2 = X(1,K)-X(1,KC)
  39          DSUM = DSUM+WSAVE(KC)*T2
  40          T2 = WSAVE(K)*T2
  41          X(1,K) = T1-T2
  42          X(1,KC) = T1+T2
  43   104 END DO
  44       MODN = MOD(N,2)
  45       IF (MODN .EQ. 0) GO TO 124
  46       X(1,NS2+1) = X(1,NS2+1)+X(1,NS2+1)
  47   124 LENX = INC*(NM1-1)  + 1
  48       LNSV = NM1 + INT(LOG(REAL(NM1))) + 4
  49       LNWK = NM1
  50 !
  51       CALL RFFT1F(NM1,INC,X,LENX,WSAVE(N+1),LNSV,WORK,                  &
  52      &            LNWK,IER1)
  53       IF (IER1 .NE. 0) THEN
  54         IER = 20
  55         CALL XERFFT ('COSTF1',-5)
  56         GO TO 200
  57       ENDIF
  58 !
  59       SNM1 = 1./FLOAT(NM1)
  60       DSUM = SNM1*DSUM
  61       IF(MOD(NM1,2) .NE. 0) GO TO 30
  62       X(1,NM1) = X(1,NM1)+X(1,NM1)
  63    30 DO 105 I=3,N,2
  64          XI = .5*X(1,I)
  65          X(1,I) = .5*X(1,I-1)
  66          X(1,I-1) = DSUM
  67          DSUM = DSUM+XI
  68   105 END DO
  69       IF (MODN .NE. 0) GO TO 117
  70       X(1,N) = DSUM
  71   117 X(1,1) = .5*X(1,1)
  72       X(1,N) = .5*X(1,N)
  73   200 RETURN
  74       END