wrfv2_fire/external/fftpack/fftpack5/mcsqb1.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: mcsqb1.f,v 1.2 2004/06/15 21:29:19 rodney Exp $
  11 !
  12 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  13
  14       SUBROUTINE MCSQB1 (LOT,JUMP,N,INC,X,WSAVE,WORK,IER)
  15       DIMENSION       X(INC,*)     ,WSAVE(*)     ,WORK(LOT,*)
  16       IER = 0
  17       LJ = (LOT-1)*JUMP+1
  18       NS2 = (N+1)/2
  19       NP2 = N+2
  20       DO 101 I=3,N,2
  21          DO 201 M=1,LJ,JUMP
  22          XIM1 = X(M,I-1)+X(M,I)
  23          X(M,I) = .5*(X(M,I-1)-X(M,I))
  24          X(M,I-1) = .5*XIM1
  25   201    CONTINUE
  26   101 END DO
  27       DO 301 M=1,LJ,JUMP
  28       X(M,1) = .5*X(M,1)
  29   301 END DO
  30       MODN = MOD(N,2)
  31       IF (MODN .NE. 0) GO TO 302
  32       DO 303 M=1,LJ,JUMP
  33       X(M,N) = .5*X(M,N)
  34   303 END DO
  35   302 CONTINUE
  36       LENX = (LOT-1)*JUMP + INC*(N-1)  + 1
  37       LNSV = N + INT(LOG(REAL(N))) + 4
  38       LNWK = LOT*N
  39 !
  40       CALL RFFTMB(LOT,JUMP,N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1)
  41       IF (IER1 .NE. 0) THEN
  42         IER = 20
  43         CALL XERFFT ('MCSQB1',-5)
  44         GO TO 400
  45       ENDIF
  46 !
  47       DO 102 K=2,NS2
  48          KC = NP2-K
  49          M1 = 0
  50          DO 202 M=1,LJ,JUMP
  51          M1 = M1 + 1
  52          WORK(M1,K) = WSAVE(K-1)*X(M,KC)+WSAVE(KC-1)*X(M,K)
  53          WORK(M1,KC) = WSAVE(K-1)*X(M,K)-WSAVE(KC-1)*X(M,KC)
  54   202    CONTINUE
  55   102 END DO
  56       IF (MODN .NE. 0) GO TO 305
  57       DO 304 M=1,LJ,JUMP
  58          X(M,NS2+1) = WSAVE(NS2)*(X(M,NS2+1)+X(M,NS2+1))
  59   304    CONTINUE
  60   305 DO 103 K=2,NS2
  61          KC = NP2-K
  62          M1 = 0
  63          DO 203 M=1,LJ,JUMP
  64             M1 = M1 + 1
  65             X(M,K) = WORK(M1,K)+WORK(M1,KC)
  66             X(M,KC) = WORK(M1,K)-WORK(M1,KC)
  67   203    CONTINUE
  68   103 END DO
  69       DO 104 M=1,LJ,JUMP
  70       X(M,1) = X(M,1)+X(M,1)
  71   104 END DO
  72   400 CONTINUE
  73       RETURN
  74       END