*FAX                                                                         FAX.1
      SUBROUTINE FAX(IFAX,N,MODE)                                            FAX.2
      INCLUDE 'param.h'                                                      FAX.3
      DIMENSION IFAX(10)                                                     FAX.4
      NN=N                                                                   FAX.5
      IF(IABS(MODE).EQ.1)GO TO 10                                            FAX.6
      IF(IABS(MODE).EQ.8)GO TO 10                                            FAX.7
      NN=N/2                                                                 FAX.8
      IF((NN+NN).EQ.N)GO TO 10                                               FAX.9
      IFAX(1)=-99                                                            FAX.10
      RETURN                                                                 FAX.11
   10 K=1                                                                    FAX.12
C     ****                                                                   FAX.13
C     ****    TEST FOR FACTORS OF 4                                          FAX.14
C     ****                                                                   FAX.15
   20 IF(MOD(NN,4).NE.0) GO TO 30                                            FAX.16
      K=K+1                                                                  FAX.17
      IFAX(K)=4                                                              FAX.18
      NN=NN/4                                                                FAX.19
      IF(NN.EQ.1) GO TO 80                                                   FAX.20
      GO TO 20                                                               FAX.21
C     ****                                                                   FAX.22
C     ****    TEST FOR EXTRA FACTOR OF 2                                     FAX.23
C     ****                                                                   FAX.24 
   30 IF(MOD(NN,2).NE.0) GO TO 40                                            FAX.25
      K=K+1                                                                  FAX.26 
      IFAX(K)=2                                                              FAX.27
      NN=NN/2                                                                FAX.28
      IF (NN.EQ.1) GO TO 80                                                  FAX.29
C     ****                                                                   FAX.30
C     ****    TEST FOR FACTORS OF 3                                          FAX.31
C     ****                                                                   FAX.32
   40 IF(MOD(NN,3).NE.0) GO TO 50                                            FAX.33
      K=K+1                                                                  FAX.34
      IFAX(K)=3                                                              FAX.35
      NN=NN/3                                                                FAX.36
      IF(NN.EQ.1) GO TO 80                                                   FAX.37
      GO TO 40                                                               FAX.38
C     ****                                                                   FAX.39
C     ****    NOW FIND REMAINING FACTORS                                     FAX.40
C     ****                                                                   FAX.41
   50 L=5                                                                    FAX.42
      INC=2                                                                  FAX.43
C     ****                                                                   FAX.44
C     ****    INC ALTERNATELY TAKES ON VALUES 2 AND 4                        FAX.45
C     ****                                                                   FAX.46
   60 IF(MOD(NN,L).NE.0) GO TO 70                                            FAX.47
      K=K+1                                                                  FAX.48
      IFAX(K)=L                                                              FAX.49
      NN=NN/L                                                                FAX.50
      IF(NN.EQ.1) GO TO 80                                                   FAX.51
      GO TO 60                                                               FAX.52
   70 L=L+INC                                                                FAX.53
      INC=6-INC                                                              FAX.54
      GO TO 60                                                               FAX.55
   80 IFAX(1)=K-1                                                            FAX.56
C     ****                                                                   FAX.57
C     ****    IFAX(1) CONTAINS NUMBER OF FACTORS                             FAX.58
C     ****                                                                   FAX.59
      NFAX=IFAX(1)                                                           FAX.60
C     ****                                                                   FAX.61
C     ****    SORT FACTORS INTO ASCENDING ORDER                              FAX.62
C     ****                                                                   FAX.63
      IF (NFAX.EQ.1) GO TO 110                                               FAX.64
      DO 100 II=2,NFAX                                                       FAX.65
       ISTOP=NFAX+2-II                                                       FAX.66
       DO 90 I=2,ISTOP                                                       FAX.67
        IF (IFAX(I+1).GE.IFAX(I)) GO TO 90                                   FAX.68
        ITEM=IFAX(I)                                                         FAX.69
        IFAX(I)=IFAX(I+1)                                                    FAX.70
        IFAX(I+1)=ITEM                                                       FAX.71
   90  CONTINUE                                                              FAX.72
  100 CONTINUE                                                               FAX.73
  110 CONTINUE                                                               FAX.74
      RETURN                                                                 FAX.75
      END                                                                    FAX.76
C                                                                            FAX.77
*DK FFT99A                                                                FFT99A.1
      SUBROUTINE FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT)                      FFT99A.2
      IMPLICIT REAL*8 (A-H,O-Z)                                           FFT99A.3
      DIMENSION A(N),WORK(N),TRIGS(N)                                     FFT99A.4
C     ****                                                                FFT99A.5
C     ****    SUBROUTINE FFT99A - PREPROCESSING STEP FOR FFT99, ISIGN=+1  FFT99A.6
C     ****    (SPECTRAL TO GRIDPOINT TRANSFORM)                           FFT99A.7
C     ****                                                                FFT99A.8
      NH=N/2                                                              FFT99A.9
      NX=N+1                                                              FFT99A.10
      INK=INC+INC                                                         FFT99A.11
C     ****                                                                FFT99A.12
C     ****    A(0) AND A(N/2)                                             FFT99A.13
C     ****                                                                FFT99A.14
      IA=1                                                                FFT99A.15
      IB=N*INC+1                                                          FFT99A.16
      JA=1                                                                FFT99A.17
      JB=2                                                                FFT99A.18
      DO 10 L=1,LOT                                                       FFT99A.19
      WORK(JA)=A(IA)+A(IB)                                                FFT99A.20
      WORK(JB)=A(IA)-A(IB)                                                FFT99A.21
      IA=IA+JUMP                                                          FFT99A.22
      IB=IB+JUMP                                                          FFT99A.23
      JA=JA+NX                                                            FFT99A.24
      JB=JB+NX                                                            FFT99A.25
   10 CONTINUE                                                            FFT99A.26
C     ****                                                                FFT99A.27
C     ****    REMAINING WAVENUMBERS                                       FFT99A.28
C     ****                                                                FFT99A.29
      IABASE=2*INC+1                                                      FFT99A.30
      IBBASE=(N-2)*INC+1                                                  FFT99A.31
      JABASE=3                                                            FFT99A.32
      JBBASE=N-1                                                          FFT99A.33
      DO 30 K=3,NH,2                                                      FFT99A.34
      IA=IABASE                                                           FFT99A.35
      IB=IBBASE                                                           FFT99A.36
      JA=JABASE                                                           FFT99A.37
      JB=JBBASE                                                           FFT99A.38
      C=TRIGS(N+K)                                                        FFT99A.39
      S=TRIGS(N+K+1)                                                      FFT99A.40
      DO 20 L=1,LOT                                                       FFT99A.41
      WORK(JA)=(A(IA)+A(IB))-                                             FFT99A.42
     *    (S*(A(IA)-A(IB))+C*(A(IA+INC)+A(IB+INC)))                       FFT99A.43
      WORK(JB)=(A(IA)+A(IB))+                                             FFT99A.44
     *    (S*(A(IA)-A(IB))+C*(A(IA+INC)+A(IB+INC)))                       FFT99A.45
      WORK(JA+1)=(C*(A(IA)-A(IB))-S*(A(IA+INC)+A(IB+INC)))+               FFT99A.46
     *    (A(IA+INC)-A(IB+INC))                                           FFT99A.47
      WORK(JB+1)=(C*(A(IA)-A(IB))-S*(A(IA+INC)+A(IB+INC)))-               FFT99A.48
     *    (A(IA+INC)-A(IB+INC))                                           FFT99A.49
      IA=IA+JUMP                                                          FFT99A.50
      IB=IB+JUMP                                                          FFT99A.51
      JA=JA+NX                                                            FFT99A.52
      JB=JB+NX                                                            FFT99A.53
   20 CONTINUE                                                            FFT99A.54
      IABASE=IABASE+INK                                                   FFT99A.55
      IBBASE=IBBASE-INK                                                   FFT99A.56
      JABASE=JABASE+2                                                     FFT99A.57
      JBBASE=JBBASE-2                                                     FFT99A.58
   30 CONTINUE                                                            FFT99A.59
      IF (IABASE.NE.IBBASE) GO TO 50                                      FFT99A.60
C     ****                                                                FFT99A.61
C     ****    WAVENUMBER N/4 (IF IT EXISTS)                               FFT99A.62
C     ****                                                                FFT99A.63
      IA=IABASE                                                           FFT99A.64
      JA=JABASE                                                           FFT99A.65
      DO 40 L=1,LOT                                                       FFT99A.66
      WORK(JA)=2.0*A(IA)                                                  FFT99A.67
      WORK(JA+1)=-2.0*A(IA+INC)                                           FFT99A.68
      IA=IA+JUMP                                                          FFT99A.69
      JA=JA+NX                                                            FFT99A.70
   40 CONTINUE                                                            FFT99A.71
   50 CONTINUE                                                            FFT99A.72
      RETURN                                                              FFT99A.73
      END                                                                 FFT99A.74
C                                                                         FFT99A.75
*DK  FFT99B                                                               FFT99B.1
      SUBROUTINE FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT)                      FFT99B.2
      IMPLICIT REAL*8 (A-H,O-Z)                                           FFT99B.3
      DIMENSION WORK(N),A(N),TRIGS(N)                                     FFT99B.4
C     ****                                                                FFT99B.5
C     ****    SUBROUTINE FFT99B - POSTPROCESSING STEP FOR FFT99, ISIGN=-1 FFT99B.6
C     ****    (GRIDPOINT TO SPECTRAL TRANSFORM)                           FFT99B.7
C     ****                                                                FFT99B.8
      NH=N/2                                                              FFT99B.9
      NX=N+1                                                              FFT99B.10
      INK=INC+INC                                                         FFT99B.11
C     ****                                                                FFT99B.12
C     ****    A(0) AND A(N/2)                                             FFT99B.13
C     ****                                                                FFT99B.14
      SCALE=1.0/DBLE(N)                                                   FFT99B.15
      IA=1                                                                FFT99B.16
      IB=2                                                                FFT99B.17
      JA=1                                                                FFT99B.18
      JB=N*INC+1                                                          FFT99B.19
      DO 10 L=1,LOT                                                       FFT99B.20
      A(JA)=SCALE*(WORK(IA)+WORK(IB))                                     FFT99B.21
      A(JB)=SCALE*(WORK(IA)-WORK(IB))                                     FFT99B.22
      A(JA+INC)=0.0                                                       FFT99B.23
      A(JB+INC)=0.0                                                       FFT99B.24
      IA=IA+NX                                                            FFT99B.25
      IB=IB+NX                                                            FFT99B.26
      JA=JA+JUMP                                                          FFT99B.27
      JB=JB+JUMP                                                          FFT99B.28
   10 CONTINUE                                                            FFT99B.29
C     ****                                                                FFT99B.30
C     ****     REMAINING WAVENUMBERS                                      FFT99B.31
C     ****                                                                FFT99B.32
      SCALE=0.5*SCALE                                                     FFT99B.33
      IABASE=3                                                            FFT99B.34
      IBBASE=N-1                                                          FFT99B.35
      JABASE=2*INC+1                                                      FFT99B.36
      JBBASE=(N-2)*INC+1                                                  FFT99B.37
      DO 30 K=3,NH,2                                                      FFT99B.38
      IA=IABASE                                                           FFT99B.39
      IB=IBBASE                                                           FFT99B.40
      JA=JABASE                                                           FFT99B.41
      JB=JBBASE                                                           FFT99B.42
      C=TRIGS(N+K)                                                        FFT99B.43
      S=TRIGS(N+K+1)                                                      FFT99B.44
      DO 20 L=1,LOT                                                       FFT99B.45
      A(JA)=SCALE*((WORK(IA)+WORK(IB))                                    FFT99B.46
     *   +(C*(WORK(IA+1)+WORK(IB+1))+S*(WORK(IA)-WORK(IB))))              FFT99B.47
      A(JB)=SCALE*((WORK(IA)+WORK(IB))                                    FFT99B.48
     *   -(C*(WORK(IA+1)+WORK(IB+1))+S*(WORK(IA)-WORK(IB))))              FFT99B.49
      A(JA+INC)=SCALE*((C*(WORK(IA)-WORK(IB))-S*(WORK(IA+1)+WORK(IB+1)))  FFT99B.50
     *    +(WORK(IB+1)-WORK(IA+1)))                                       FFT99B.51
      A(JB+INC)=SCALE*((C*(WORK(IA)-WORK(IB))-S*(WORK(IA+1)+WORK(IB+1)))  FFT99B.52
     *    -(WORK(IB+1)-WORK(IA+1)))                                       FFT99B.53
      IA=IA+NX                                                            FFT99B.54
      IB=IB+NX                                                            FFT99B.55
      JA=JA+JUMP                                                          FFT99B.56
      JB=JB+JUMP                                                          FFT99B.57
   20 CONTINUE                                                            FFT99B.58
      IABASE=IABASE+2                                                     FFT99B.59
      IBBASE=IBBASE-2                                                     FFT99B.60
      JABASE=JABASE+INK                                                   FFT99B.61
      JBBASE=JBBASE-INK                                                   FFT99B.62
   30 CONTINUE                                                            FFT99B.63
      IF (IABASE.NE.IBBASE) GO TO 50                                      FFT99B.64
C     ****                                                                FFT99B.65
C     ****    WAVENUMBER N/4 (IF IT EXISTS)                               FFT99B.66
C     ****                                                                FFT99B.67
      IA=IABASE                                                           FFT99B.68
      JA=JABASE                                                           FFT99B.69
      SCALE=2.0*SCALE                                                     FFT99B.70
      DO 40 L=1,LOT                                                       FFT99B.71
      A(JA)=SCALE*WORK(IA)                                                FFT99B.72
      A(JA+INC)=-SCALE*WORK(IA+1)                                         FFT99B.73
      IA=IA+NX                                                            FFT99B.74
      JA=JA+JUMP                                                          FFT99B.75
   40 CONTINUE                                                            FFT99B.76
   50 CONTINUE                                                            FFT99B.77
      RETURN                                                              FFT99B.78
      END                                                                 FFT99B.79
C                                                                         FFT99B.80
*DK FFT991                                                                FFT991.1
      SUBROUTINE FFT991(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)           FFT991.2
      IMPLICIT REAL*8 (A-H,O-Z)                                           FFT991.3
      DIMENSION A(N),WORK(N),TRIGS(N),IFAX(1)                             FFT991.4
C     ****                                                                FFT991.5
C     ****    SUBROUTINE "FFT991" - MULTIPLE REAL/HALF-COMPLEX PERIODIC   FFT991.6
C     ****    FAST FOURIER TRANSFORM                                      FFT991.7
C     ****                                                                FFT991.8
C     ****    SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO   FFT991.9
C     ****    THAT IN MRFFT2                                              FFT991.10
C     ****                                                                FFT991.11
C     ****    PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM  FFT991.12
C     ****    IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL.12  FFT991.13
C     ****    (1970), 315-337)                                            FFT991.14
C     ****                                                                FFT991.15
C     ****    A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA             FFT991.16
C     ****    WORK IS AN AREA OF SIZE (N+1)*LOT                           FFT991.17
C     ****    TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES FFT991.18
C     ****    IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2        FFT991.19
C     ****    INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'              FFT991.20
C     ****        (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)              FFT991.21
C     ****    JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR FFT991.22
C     ****    N IS THE LENGTH OF THE DATA VECTORS                         FFT991.23
C     ****    LOT IS THE NUMBER OF DATA VECTORS                           FFT991.24
C     ****    ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT         FFT991.25
C     ****          = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL         FFT991.26
C     ****                                                                FFT991.27
C     ****    ORDERING OF COEFFICIENTS:                                   FFT991.28
C     ****       A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)          FFT991.29
C     ****        WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED           FFT991.30
C     ****                                                                FFT991.31
C     ****    ORDERING OF DATA:                                           FFT991.32
C     ****        X(0),X(1),X(2),...,X(N-1)                               FFT991.33
C     ****                                                                FFT991.34
C     *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER                          FFT991.35
C     ****                                                                FFT991.36
C     ****    DEFINITION OF TRANSFORMS:                                   FFT991.37
C     ****                                                                FFT991.38
C     ****    ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))     FFT991.39
C     ****        WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)           FFT991.40
C     ****                                                                FFT991.41
C     ****    ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) FFT991.42
C     ****              B(K)=-(1/N)*SUM(J=0,..,N-1)(X(J)*SIN(2*J*K*PI/N)) FFT991.43
C     ****                                                                FFT991.44
      NFAX=IFAX(1)                                                        FFT991.45
      NX=N+1                                                              FFT991.46
      NH=N/2                                                              FFT991.47
      INK=INC+INC                                                         FFT991.48
      IF (ISIGN.EQ.+1) GO TO 30                                           FFT991.49
C     ****                                                                FFT991.50
C     ****    IF NECESSARY, TRANSFER DATA TO WORK AREA                    FFT991.51
C     ****                                                                FFT991.52
      IGO=50                                                              FFT991.53
      IF (MOD(NFAX,2).EQ.1) GOTO 40                                       FFT991.54
      IBASE=1                                                             FFT991.55
      JBASE=1                                                             FFT991.56
      DO 20 L=1,LOT                                                       FFT991.57
      I=IBASE                                                             FFT991.58
      J=JBASE                                                             FFT991.59
      DO 10 M=1,N                                                         FFT991.60
      WORK(J)=A(I)                                                        FFT991.61
      I=I+INC                                                             FFT991.62
      J=J+1                                                               FFT991.63
   10 CONTINUE                                                            FFT991.64
      IBASE=IBASE+JUMP                                                    FFT991.65
      JBASE=JBASE+NX                                                      FFT991.66
   20 CONTINUE                                                            FFT991.67
      IGO=60                                                              FFT991.68
      GO TO 40                                                            FFT991.69
C     ****                                                                FFT991.70
C     ****    PREPROCESSING (ISIGN=+1)                                    FFT991.71
C     ****                                                                FFT991.72
   30 CONTINUE                                                            FFT991.73
      CALL FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT)                            FFT991.74
      IGO=60                                                              FFT991.75
C     ****                                                                FFT991.76
C     ****    COMPLEX TRANSFORM                                           FFT991.77
C     ****                                                                FFT991.78
   40 CONTINUE                                                            FFT991.79
      IA=1                                                                FFT991.80
      LA=1                                                                FFT991.81
      DO 80 K=1,NFAX                                                      FFT991.82
      IF (IGO.EQ.60) GO TO 60                                             FFT991.83
   50 CONTINUE                                                            FFT991.84
      CALL VPASSM(A(IA),A(IA+INC),WORK(1),WORK(2),TRIGS,                  FFT991.85
     *   INK,2,JUMP,NX,LOT,NH,IFAX(K+1),LA)                               FFT991.86
      IGO=60                                                              FFT991.87
      GO TO 70                                                            FFT991.88
   60 CONTINUE                                                            FFT991.89
      CALL VPASSM(WORK(1),WORK(2),A(IA),A(IA+INC),TRIGS,                  FFT991.90
     *    2,INK,NX,JUMP,LOT,NH,IFAX(K+1),LA)                              FFT991.91
      IGO=50                                                              FFT991.92
   70 CONTINUE                                                            FFT991.93
      LA=LA*IFAX(K+1)                                                     FFT991.94
   80 CONTINUE                                                            FFT991.95
      IF (ISIGN.EQ.-1) GO TO 130                                          FFT991.96
C     ****                                                                FFT991.97
C     ****    IF NECESSARY, TRANSFER DATA FROM WORK AREA                  FFT991.98
C     ****                                                                FFT991.99
      IF (MOD(NFAX,2).EQ.1) GO TO 110                                     FFT991.100
      IBASE=1                                                             FFT991.101
      JBASE=1                                                             FFT991.102
      DO 100 L=1,LOT                                                      FFT991.103
      I=IBASE                                                             FFT991.104
      J=JBASE                                                             FFT991.105
      DO 90 M=1,N                                                         FFT991.106
      A(J)=WORK(I)                                                        FFT991.107
      I=I+1                                                               FFT991.108
      J=J+INC                                                             FFT991.109
   90 CONTINUE                                                            FFT991.110
      IBASE=IBASE+NX                                                      FFT991.111
      JBASE=JBASE+JUMP                                                    FFT991.112
  100 CONTINUE                                                            FFT991.113
C     ****                                                                FFT991.114
C     ****    FILL IN ZEROS AT END                                        FFT991.115
C     ****                                                                FFT991.116
  110 CONTINUE                                                            FFT991.117
      IB=N*INC+1                                                          FFT991.118
      DO 120 L=1,LOT                                                      FFT991.119
      A(IB)=0.0                                                           FFT991.120
      A(IB+INC)=0.0                                                       FFT991.121
      IB=IB+JUMP                                                          FFT991.122
  120 CONTINUE                                                            FFT991.123
      GO TO 140                                                           FFT991.124
C     ****                                                                FFT991.125
C     ****    POSTPROCESSING (ISIGN=-1):                                  FFT991.126
C     ****                                                                FFT991.127
  130 CONTINUE                                                            FFT991.128
      CALL FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT)                            FFT991.129
  140 CONTINUE                                                            FFT991.130
      RETURN                                                              FFT991.131
      END                                                                 FFT991.132
C                                                                         FFT991.133
*DK FFTRIG                                                                FFTRIG.1
      SUBROUTINE FFTRIG(TRIGS,N,MODE)                                     FFTRIG.2
      INCLUDE 'param.h'                                                   FFTRIG.3
      DIMENSION TRIGS(1)                                                  FFTRIG.4
      PI=2.0*DASIN(1.0D0)                                                 FFTRIG.5
      IMODE=IABS(MODE)                                                    FFTRIG.6
      NN=N                                                                FFTRIG.7
      IF(IMODE.GT.1.AND.IMODE.LT.6)NN=N/2                                 FFTRIG.8
      DEL=(PI+PI)/DBLE(NN)                                                FFTRIG.9
      L=NN+NN                                                             FFTRIG.10
      DO I=1,L,2                                                          FFTRIG.11
       ANGLE=0.5*DBLE(I-1)*DEL                                            FFTRIG.12
       TRIGS(I)=DCOS(ANGLE)                                               FFTRIG.13
       TRIGS(I+1)=DSIN(ANGLE)                                             FFTRIG.14
      ENDDO                                                               FFTRIG.15
      IF(IMODE.EQ.1)RETURN                                                FFTRIG.16
      IF(IMODE.EQ.8)RETURN                                                FFTRIG.17
      DEL=0.5*DEL                                                         FFTRIG.18
      NH=(NN+1)/2                                                         FFTRIG.19
      L=NH+NH                                                             FFTRIG.20
      LA=NN+NN                                                            FFTRIG.21
      DO I=1,L,2                                                          FFTRIG.22
       ANGLE=0.5*DBLE(I-1)*DEL                                            FFTRIG.23
       TRIGS(LA+I)=DCOS(ANGLE)                                            FFTRIG.24
       TRIGS(LA+I+1)=DSIN(ANGLE)                                          FFTRIG.25
      ENDDO                                                               FFTRIG.26
      IF(IMODE.LE.3)RETURN                                                FFTRIG.27
      DEL=0.5*DEL                                                         FFTRIG.28
      LA=LA+NN                                                            FFTRIG.29
      IF(MODE.EQ.5)GO TO 40                                               FFTRIG.30
      DO I=2,NN                                                           FFTRIG.31
       ANGLE=DBLE(I-1)*DEL                                                FFTRIG.32
       TRIGS(LA+I)=2.0*DSIN(ANGLE)                                        FFTRIG.33
      ENDDO                                                               FFTRIG.34
      RETURN                                                              FFTRIG.35
   40 CONTINUE                                                            FFTRIG.36
      DEL=0.5*DEL                                                         FFTRIG.37
      DO I=2,N                                                            FFTRIG.38
       ANGLE=DBLE(I-1)*DEL                                                FFTRIG.39
       TRIGS(LA+I)=DSIN(ANGLE)                                            FFTRIG.40
      ENDDO                                                               FFTRIG.41
      RETURN                                                              FFTRIG.42
      END                                                                 FFTRIG.43
C                                                                         FFTRIG.44
*DK FILTER                                                                FILTER.1
      SUBROUTINE FILTER(NX,NLEV,KUTJ)                                     FILTER.2
C     ****                                                                FILTER.3
C     ****         SUBROUTINE TO REMOVE LONGITUDINAL WAVES OF PROGNOSTIC  FILTER.4
C     ****         VARIABLES BY FFTS.                                     FILTER.5
C     ****                                                                FILTER.6
      INCLUDE 'param.h'                                                   FILTER.7
      INCLUDE 'crfft.h'                                                   FILTER.8
      INCLUDE 'const.h'                                                   FILTER.9
      N1=2*KUTJ+3                                                         FILTER.10
      N2=IMAX+2                                                           FILTER.11
      IF(N1.GT.N2)RETURN                                                  FILTER.12
      DO N=1,NLEV                                                         FILTER.13
       NXK=NX+N-1                                                         FILTER.14
C      ****                                                               FILTER.15
C      ****                   SET UP ARRAY OF SEQUENCE TO BE TRANSFORMED  FILTER.16
C      ****                                                               FILTER.17
       DO I=1,IMAX                                                        FILTER.18
        FX(I,N)=F(I+2,NXK)                                                FILTER.19
       ENDDO                                                              FILTER.20
      ENDDO                                                               FILTER.21
      DO I=IMAX+1,IMAXP2                                                  FILTER.22
       DO N=1,NLEV                                                        FILTER.23
        FX(I,N)=0.                                                        FILTER.24
       ENDDO                                                              FILTER.25
      ENDDO                                                               FILTER.26
C     ****                                                                FILTER.27
C     ****                   CALL FFT ROUTINE FOR FORWARD TRANSFORM       FILTER.28
C     ****                                                                FILTER.29
      CALL FFT991(FX,WFFT,TRIGS,IFAX,1,IMAXP2,IMAX,NLEV,-1)               FILTER.30
C     ****                                                                FILTER.31
C     ****                   REMOVE WAVE NUMBERS GREATER THAN KUTJ        FILTER.32
C     ****                                                                FILTER.33
      DO N=1,NLEV                                                         FILTER.34
       DO I=N1,N2                                                         FILTER.35
        FX(I,N)=0.0                                                       FILTER.36
       ENDDO                                                              FILTER.37
c       do i=1,imax
c        write(6,*)FX(I,N),i,n
c       enddo
      ENDDO                                                               FILTER.38
c      pause
      CALL FFT991(FX,WFFT,TRIGS,IFAX,1,IMAXP2,IMAX,NLEV,1)                FILTER.39
C     ****                                                                FILTER.40
C     ****                   STORE SMOOTHED FIELD IN SCM BUFFER           FILTER.41
C     ****                                                                FILTER.42
      DO N=1,NLEV                                                         FILTER.43
       NXK=NX+N-1                                                         FILTER.44
       DO I=1,IMAX                                                        FILTER.45
        F(I+2,NXK)=FX(I,N)                                                FILTER.46
       ENDDO                                                              FILTER.47
      ENDDO                                                               FILTER.48
      RETURN                                                              FILTER.49
      END                                                                 FILTER.50
C                                                                         FILTER.51
*DK FILTER2                                                              FILTER2.1 
      SUBROUTINE FILTER2(NX,NLEV,KUTJ)                                   FILTER2.2
C     ****                                                               FILTER2.3
C     ****         SUBROUTINE TO REMOVE LONGITUDINAL WAVES OF PROGNOSTIC FILTER2.4
C     ****         VARIABLES BY FFTS.                                    FILTER2.5
C     ****                                                               FILTER2.6
      INCLUDE 'param.h'                                                  FILTER2.7
      INCLUDE 'crfft.h'                                                  FILTER2.8
      INCLUDE 'const.h'                                                  FILTER2.9
C     ****                                                               FILTER2.10
C     ****     ARRAY GIVING POWER OF SMOOTHER CORRESPONDING  TO EACH     FILTER2.11
C     ****       CUTOFF FREQUENCY                                        FILTER2.12
C     ****                                                               FILTER2.13
      DIMENSION NN(ZIMX/4)                                               FILTER2.14
C     DATA NN/91,40,23,15,10,8,6,5,4,3,3,3,2,2,2,2,1,1/                  FILTER2.15
      DATA NN/45,20,11,7,5,4,3,2,2,1,1,1,1,1,1,1,1,1/                    FILTER2.16
C     DATA NN/22,10,5,3,2,2,1,1,1,1,1,1,1,1,1,1,1,1/                     FILTER2.17
C     ****                                                               FILTER2.18
C     ****                   SET UP ARRAY OF SEQUENCE TO BE TRANSFORMED  FILTER2.19
C     ****                                                               FILTER2.20
      IF(KUTJ.GT.18)RETURN                                               FILTER2.21
      DO N=1,NLEV                                                        FILTER2.22
       NXK=NX+N-1                                                        FILTER2.23
       DO I=1,IMAX                                                       FILTER2.24
        FX(I,N)=F(I+2,NXK)                                               FILTER2.25
       ENDDO                                                             FILTER2.26
      ENDDO                                                              FILTER2.27
      DO I=IMAX+1,IMAXP2                                                 FILTER2.28
       DO N=1,NLEV                                                       FILTER2.29
        FX(I,N)=0.                                                       FILTER2.30
       ENDDO                                                             FILTER2.31
      ENDDO                                                              FILTER2.32
C     ****                                                               FILTER2.33
C     ****                   CALL FFT ROUTINE FOR FORWARD TRANSFORM      FILTER2.34
C     ****                                                               FILTER2.35
      CALL FFT991(FX,WFFT,TRIGS,IFAX,1,IMAXP2,IMAX,NLEV,-1)              FILTER2.36
C     ****                                                               FILTER2.37
C     ****                   MULTIPLY BY SMOOTHING FUNCTION              FILTER2.38
C     ****                                                               FILTER2.39
      DO I=1,IMAX                                                        FILTER2.40
       DO N=1,NLEV                                                       FILTER2.41
      FX(I+2,N)=FX(I+2,N)*DCOS((DBLE(I+1)/2)*C(1)/2)**(2*NN(KUTJ))       FILTER2.42
       ENDDO                                                             FILTER2.44
      ENDDO                                                              FILTER2.45
C     ****                                                               FILTER2.46
C     ****                   CALL FFT ROUTINE FOR REVERSE TRANSFORM      FILTER2.47
C     ****                                                               FILTER2.48
      CALL FFT991(FX,WFFT,TRIGS,IFAX,1,IMAXP2,IMAX,NLEV,+1)              FILTER2.49
C     ****                                                               FILTER2.50
C     ****                   STORE SMOOTHED FIELD IN SCM BUFFER          FILTER2.51
C     ****                                                               FILTER2.52
      DO N=1,NLEV                                                        FILTER2.53
       NXK=NX+N-1                                                        FILTER2.54
       DO I=1,IMAX                                                       FILTER2.55
        F(I+2,NXK)=FX(I,N)                                               FILTER2.56
       ENDDO                                                             FILTER2.57
      ENDDO                                                              FILTER2.58
      RETURN                                                             FILTER2.59
      END                                                                FILTER2.60
C                                                                        FILTER2.61
*DK SET99                                                                  SET99.1
      SUBROUTINE SET99(TRIGS,IFAX,N)                                       SET99.1
C     ****                                                                 SET99.3
C     ****     INITIALIZE FFT ROUTINE                                      SET99.4
C     ****                                                                 SET99.5
      INCLUDE 'param.h'                                                    SET99.6
      DIMENSION IFAX(13),TRIGS(1)                                          SET99.7
C     ****                                                                 SET99.8
C     ****     MODE 3 IS USED FOR REAL/HALF-COMPLEX TRANSFORMS. IT IS      SET99.9
C     ****     POSSIBLE TO DO COMPLEX/COMPLEX TRANSFORMS WITH OTHER        SET99.10
C     ****     VALUES OF MODE,BUT DOCUMENTATION OF THE DETAILS WERE        SET99.11
C     ****     NOT AVAILABLE WHEN THIS ROUTINE WAS WRITTEN.                SET99.12
C     ****                                                                 SET99.13
      DATA MODE /3/                                                        SET99.14
      CALL FAX (IFAX,N,MODE)                                               SET99.15
      I=IFAX(1)                                                            SET99.16
      IF(IFAX(I+1).GT.5.OR.N.LE.4) IFAX(1)=-99                             SET99.17
      IF(IFAX(1) .LE. 0 ) THEN                                             SET99.18
       WRITE(6,*) ' SET99 -- INVALID N'                                    SET99.19
       STOP'SET99'                                                         SET99.20
      ENDIF                                                                SET99.21
      CALL FFTRIG(TRIGS,N,MODE)                                            SET99.22
      RETURN                                                               SET99.23
      END                                                                  SET99.24
C                                                                          SET99.25
*DK VPASSM                                                                VPASSM.1
      SUBROUTINE VPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA)  VPASSM.2
      IMPLICIT REAL*8 (A-H,O-Z)                                           VPASSM.3
      DIMENSION A(N),B(N),C(N),D(N),TRIGS(N)                              VPASSM.4
C     ****                                                                VPASSM.5
C     ****       SUBROUTINE "VPASSM" - MULTIPLE VERSION OF "VPASSA"       VPASSM.6
C     ****       PERFORMS ONE PASS THROUGH DATA                           VPASSM.7
C     ****       AS PART OF MULTIPLE COMPLEX FFT ROUTINE                  VPASSM.8
C     ****       A IS FIRST REAL INPUT VECTOR                             VPASSM.9
C     ****      B IS FIRST IMAGINARY INPUT VECTOR                         VPASSM.10
C     ****       C IS FIRST REAL OUTPUT VECTOR                            VPASSM.11
C     ****       D IS FIRST IMAGINARY OUTPUT VECTOR                       VPASSM.12
C     ****       TRIGS IS PRECALCULATED TABLE OF SINES " COSINES          VPASSM.13
C     ****       INC1 IS ADDRESSING INCREMENT FOR A AND B                 VPASSM.14
C     ****       INC2 IS ADDRESSING INCREMENT FOR C AND D                 VPASSM.15
C     ****       INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S           VPASSM.16
C     ****       INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S           VPASSM.17
C     ****       LOT IS THE NUMBER OF VECTORS                             VPASSM.18
C     ****       N IS LENGTH OF VECTORS                                   VPASSM.19
C     ****       IFAC IS CURRENT FACTOR OF N                              VPASSM.20
C     ****       LA IS PRODUCT OF PREVIOUS FACTORS                        VPASSM.21
C     ****                                                                VPASSM.22
      DATA SIN36/0.587785252292473/,COS36/0.809016994374947/,             VPASSM.23
     *     SIN72/0.951056516295154/,COS72/0.309016994374947/,             VPASSM.24
     *     SIN60/0.866025403784437/                                       VPASSM.25
      M=N/IFAC                                                            VPASSM.26
      IINK=M*INC1                                                         VPASSM.27
      JINK=LA*INC2                                                        VPASSM.28
      JUMP=(IFAC-1)*JINK                                                  VPASSM.29
      IBASE=0                                                             VPASSM.30
      JBASE=0                                                             VPASSM.31
      IGO=IFAC-1                                                          VPASSM.32
      IF (IGO.GT.4) RETURN                                                VPASSM.33
      GO TO (10,50,90,130),IGO                                            VPASSM.34
C     ****                                                                VPASSM.35
C     ****      CODING FOR FACTOR 2                                       VPASSM.36
C     ****                                                                VPASSM.37
   10 IA=1                                                                VPASSM.38
      JA=1                                                                VPASSM.39
      IB=IA+IINK                                                          VPASSM.40
      JB=JA+JINK                                                          VPASSM.41
      DO 20 L=1,LA                                                        VPASSM.42
      I=IBASE                                                             VPASSM.43
      J=JBASE                                                             VPASSM.44
      DO 15 IJK=1,LOT                                                     VPASSM.45
      C(JA+J)=A(IA+I)+A(IB+I)                                             VPASSM.46
      D(JA+J)=B(IA+I)+B(IB+I)                                             VPASSM.47
      C(JB+J)=A(IA+I)-A(IB+I)                                             VPASSM.48
      D(JB+J)=B(IA+I)-B(IB+I)                                             VPASSM.49
      I=I+INC3                                                            VPASSM.50
      J=J+INC4                                                            VPASSM.51
   15 CONTINUE                                                            VPASSM.52
      IBASE=IBASE+INC1                                                    VPASSM.53
      JBASE=JBASE+INC2                                                    VPASSM.54
   20 CONTINUE                                                            VPASSM.55
      IF (LA.EQ.M) RETURN                                                 VPASSM.56
      LA1=LA+1                                                            VPASSM.57
      JBASE=JBASE+JUMP                                                    VPASSM.58
      DO 40 K=LA1,M,LA                                                    VPASSM.59
      KB=K+K-2                                                            VPASSM.60
      C1=TRIGS(KB+1)                                                      VPASSM.61
      S1=TRIGS(KB+2)                                                      VPASSM.62
      DO 30 L=1,LA                                                        VPASSM.63
      I=IBASE                                                             VPASSM.64
      J=JBASE                                                             VPASSM.65
      DO 25 IJK=1,LOT                                                     VPASSM.66
      C(JA+J)=A(IA+I)+A(IB+I)                                             VPASSM.67
      D(JA+J)=B(IA+I)+B(IB+I)                                             VPASSM.68
      C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))                   VPASSM.69
      D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))                   VPASSM.70
      I=I+INC3                                                            VPASSM.71
      J=J+INC4                                                            VPASSM.72
   25 CONTINUE                                                            VPASSM.73
      IBASE=IBASE+INC1                                                    VPASSM.74
      JBASE=JBASE+INC2                                                    VPASSM.75
   30 CONTINUE                                                            VPASSM.76
      JBASE=JBASE+JUMP                                                    VPASSM.77
   40 CONTINUE                                                            VPASSM.78
      RETURN                                                              VPASSM.79
C     ****                                                                VPASSM.80
C     ****       CODING FOR FACTOR 3                                      VPASSM.81
C     ****                                                                VPASSM.82
   50 IA=1                                                                VPASSM.83
      JA=1                                                                VPASSM.84
      IB=IA+IINK                                                          VPASSM.85
      JB=JA+JINK                                                          VPASSM.86
      IC=IB+IINK                                                          VPASSM.87
      JC=JB+JINK                                                          VPASSM.88
      DO 60 L=1,LA                                                        VPASSM.89
      I=IBASE                                                             VPASSM.90
      J=JBASE                                                             VPASSM.91
      DO 55 IJK=1,LOT                                                     VPASSM.92
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))                                   VPASSM.93
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))                                   VPASSM.94
      C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))   VPASSM.95
      C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))   VPASSM.96
      D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))   VPASSM.97
      D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))   VPASSM.98
      I=I+INC3                                                            VPASSM.99
      J=J+INC4                                                            VPASSM.100
   55 CONTINUE                                                            VPASSM.101
      IBASE=IBASE+INC1                                                    VPASSM.102
      JBASE=JBASE+INC2                                                    VPASSM.103
   60 CONTINUE                                                            VPASSM.104
      IF (LA.EQ.M) RETURN                                                 VPASSM.105
      LA1=LA+1                                                            VPASSM.106
      JBASE=JBASE+JUMP                                                    VPASSM.107
      DO 80 K=LA1,M,LA                                                    VPASSM.108
      KB=K+K-2                                                            VPASSM.109
      KC=KB+KB                                                            VPASSM.110
      C1=TRIGS(KB+1)                                                      VPASSM.111
      S1=TRIGS(KB+2)                                                      VPASSM.112
      C2=TRIGS(KC+1)                                                      VPASSM.113
      S2=TRIGS(KC+2)                                                      VPASSM.114
      DO 70 L=1,LA                                                        VPASSM.115
      I=IBASE                                                             VPASSM.116
      J=JBASE                                                             VPASSM.117
      DO 65 IJK=1,LOT                                                     VPASSM.118
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))                                   VPASSM.119
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))                                   VPASSM.120
      C(JB+J)=                                                            VPASSM.121
     *    C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))  VPASSM.122
     *   -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))  VPASSM.123
      D(JB+J)=                                                            VPASSM.124
     *    S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))  VPASSM.125
     *   +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))  VPASSM.126
      C(JC+J)=                                                            VPASSM.127
     *    C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))  VPASSM.128
     *   -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))  VPASSM.129
      D(JC+J)=                                                            VPASSM.130
     *    S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))  VPASSM.131
     *   +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))  VPASSM.132
      I=I+INC3                                                            VPASSM.133
      J=J+INC4                                                            VPASSM.134
   65 CONTINUE                                                            VPASSM.135
      IBASE=IBASE+INC1                                                    VPASSM.136
      JBASE=JBASE+INC2                                                    VPASSM.137
   70 CONTINUE                                                            VPASSM.138
      JBASE=JBASE+JUMP                                                    VPASSM.139
   80 CONTINUE                                                            VPASSM.140
      RETURN                                                              VPASSM.141
C     ****                                                                VPASSM.142
C     ****       CODING FOR FACTOR 4                                      VPASSM.143
C     ****                                                                VPASSM.144
   90 IA=1                                                                VPASSM.145
      JA=1                                                                VPASSM.146
      IB=IA+IINK                                                          VPASSM.147
      JB=JA+JINK                                                          VPASSM.148
      IC=IB+IINK                                                          VPASSM.149
      JC=JB+JINK                                                          VPASSM.150
      ID=IC+IINK                                                          VPASSM.151
      JD=JC+JINK                                                          VPASSM.152
      DO 100 L=1,LA                                                       VPASSM.153
      I=IBASE                                                             VPASSM.154
      J=JBASE                                                             VPASSM.155
      DO 95 IJK=1,LOT                                                     VPASSM.156
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))                         VPASSM.157
      C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))                         VPASSM.158
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))                         VPASSM.159
      D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))                         VPASSM.160
      C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))                         VPASSM.161
      C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))                         VPASSM.162
      D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))                         VPASSM.163
      D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))                         VPASSM.164
      I=I+INC3                                                            VPASSM.165
      J=J+INC4                                                            VPASSM.166
   95 CONTINUE                                                            VPASSM.167
      IBASE=IBASE+INC1                                                    VPASSM.168
      JBASE=JBASE+INC2                                                    VPASSM.169
  100 CONTINUE                                                            VPASSM.170
      IF (LA.EQ.M) RETURN                                                 VPASSM.171
      LA1=LA+1                                                            VPASSM.172
      JBASE=JBASE+JUMP                                                    VPASSM.173
      DO 120 K=LA1,M,LA                                                   VPASSM.174
      KB=K+K-2                                                            VPASSM.175
      KC=KB+KB                                                            VPASSM.176
      KD=KC+KB                                                            VPASSM.177
      C1=TRIGS(KB+1)                                                      VPASSM.178
      S1=TRIGS(KB+2)                                                      VPASSM.179
      C2=TRIGS(KC+1)                                                      VPASSM.180
      S2=TRIGS(KC+2)                                                      VPASSM.181
      C3=TRIGS(KD+1)                                                      VPASSM.182
      S3=TRIGS(KD+2)                                                      VPASSM.183
      DO 110 L=1,LA                                                       VPASSM.184
      I=IBASE                                                             VPASSM.185
      J=JBASE                                                             VPASSM.186
      DO 105 IJK=1,LOT                                                    VPASSM.187
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))                         VPASSM.188
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))                         VPASSM.189
      C(JC+J)=                                                            VPASSM.190
     *    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))                        VPASSM.191
     *   -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))                        VPASSM.192
      D(JC+J)=                                                            VPASSM.193
     *    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))                        VPASSM.194
     *   +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))                        VPASSM.195
      C(JB+J)=                                                            VPASSM.196
     *    C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))                        VPASSM.197
     *   -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))                        VPASSM.198
      D(JB+J)=                                                            VPASSM.199
     *    S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))                        VPASSM.200
     *   +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))                        VPASSM.201
      C(JD+J)=                                                            VPASSM.202
     *    C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))                        VPASSM.203
     *   -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))                        VPASSM.204
      D(JD+J)=                                                            VPASSM.205
     *    S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))                        VPASSM.206
     *   +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))                        VPASSM.207
      I=I+INC3                                                            VPASSM.208
      J=J+INC4                                                            VPASSM.209
  105 CONTINUE                                                            VPASSM.210
      IBASE=IBASE+INC1                                                    VPASSM.211
      JBASE=JBASE+INC2                                                    VPASSM.212
  110 CONTINUE                                                            VPASSM.213
      JBASE=JBASE+JUMP                                                    VPASSM.214
  120 CONTINUE                                                            VPASSM.215
      RETURN                                                              VPASSM.216
C     ****                                                                VPASSM.217
C     ****       CODING FOR FACTOR 5                                      VPASSM.218
C     ****                                                                VPASSM.219
  130 IA=1                                                                VPASSM.220
      JA=1                                                                VPASSM.221
      IB=IA+IINK                                                          VPASSM.222
      JB=JA+JINK                                                          VPASSM.223
      IC=IB+IINK                                                          VPASSM.224
      JC=JB+JINK                                                          VPASSM.225
      ID=IC+IINK                                                          VPASSM.226
      JD=JC+JINK                                                          VPASSM.227
      IE=ID+IINK                                                          VPASSM.228
      JE=JD+JINK                                                          VPASSM.229
      DO 140 L=1,LA                                                       VPASSM.230
      I=IBASE                                                             VPASSM.231
      J=JBASE                                                             VPASSM.232
      DO 135 IJK=1,LOT                                                    VPASSM.233
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))                 VPASSM.234
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))                 VPASSM.235
      C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))   VPASSM.236
     *  -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))                VPASSM.237
      C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))   VPASSM.238
     *  +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))                VPASSM.239
      D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))   VPASSM.240
     *  +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))                VPASSM.241
      D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))   VPASSM.242
     *  -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))                VPASSM.243
      C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))   VPASSM.244
     *  -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))                VPASSM.245
      C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))   VPASSM.246
     *  +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))                VPASSM.247
      D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))   VPASSM.248
     *  +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))                VPASSM.249
      D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))   VPASSM.250
     *  -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))                VPASSM.251
      I=I+INC3                                                            VPASSM.252
      J=J+INC4                                                            VPASSM.253
  135 CONTINUE                                                            VPASSM.254
      IBASE=IBASE+INC1                                                    VPASSM.255
      JBASE=JBASE+INC2                                                    VPASSM.256
  140 CONTINUE                                                            VPASSM.257
      IF (LA.EQ.M) RETURN                                                 VPASSM.258
      LA1=LA+1                                                            VPASSM.259
      JBASE=JBASE+JUMP                                                    VPASSM.260
      DO 160 K=LA1,M,LA                                                   VPASSM.261
      KB=K+K-2                                                            VPASSM.262
      KC=KB+KB                                                            VPASSM.263
      KD=KC+KB                                                            VPASSM.264
      KE=KD+KB                                                            VPASSM.265
      C1=TRIGS(KB+1)                                                      VPASSM.266
      S1=TRIGS(KB+2)                                                      VPASSM.267
      C2=TRIGS(KC+1)                                                      VPASSM.268
      S2=TRIGS(KC+2)                                                      VPASSM.269
      C3=TRIGS(KD+1)                                                      VPASSM.270
      S3=TRIGS(KD+2)                                                      VPASSM.271
      C4=TRIGS(KE+1)                                                      VPASSM.272
      S4=TRIGS(KE+2)                                                      VPASSM.273
      DO 150 L=1,LA                                                       VPASSM.274
      I=IBASE                                                             VPASSM.275
      J=JBASE                                                             VPASSM.276
      DO 145 IJK=1,LOT                                                    VPASSM.277
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))                 VPASSM.278
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))                 VPASSM.279
      C(JB+J)=                                                            VPASSM.280
     *    C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))   VPASSM.281
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))           VPASSM.282
     *   -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))   VPASSM.283
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))           VPASSM.284
      D(JB+J)=                                                            VPASSM.285
     *    S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))   VPASSM.286
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))           VPASSM.287
     *   +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))   VPASSM.288
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))           VPASSM.289
      C(JE+J)=                                                            VPASSM.290
     *    C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))   VPASSM.291
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))           VPASSM.292
     *   -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))   VPASSM.293
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))           VPASSM.294
      D(JE+J)=                                                            VPASSM.295
     *    S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))   VPASSM.296
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))           VPASSM.297
     *   +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))   VPASSM.298
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))           VPASSM.299
      C(JC+J)=                                                            VPASSM.300
     *    C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))   VPASSM.301
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))           VPASSM.302
     *   -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))   VPASSM.303
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))           VPASSM.304
      D(JC+J)=                                                            VPASSM.305
     *    S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))   VPASSM.306
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))           VPASSM.307
     *   +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))   VPASSM.308
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))           VPASSM.309
      C(JD+J)=                                                            VPASSM.310
     *    C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))   VPASSM.311
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))           VPASSM.312
     *   -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))   VPASSM.313
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))           VPASSM.314
      D(JD+J)=                                                            VPASSM.315
     *    S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))   VPASSM.316
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))           VPASSM.317
     *   +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))   VPASSM.318
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))           VPASSM.319
      I=I+INC3                                                            VPASSM.320
      J=J+INC4                                                            VPASSM.321
  145 CONTINUE                                                            VPASSM.322
      IBASE=IBASE+INC1                                                    VPASSM.323
      JBASE=JBASE+INC2                                                    VPASSM.324
  150 CONTINUE                                                            VPASSM.325
      JBASE=JBASE+JUMP                                                    VPASSM.326
  160 CONTINUE                                                            VPASSM.327
      RETURN                                                              VPASSM.328
      END                                                                 VPASSM.329
C                                                                         VPASSM.330
*DK N1FILTER
      SUBROUTINE N1FILTER(NX,J,NKAX,LEVEL)
      INCLUDE 'n1param.h'
      INCLUDE 'n1array.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n1const.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'index.h'
      DIMENSION XR(IMX),XI(IMX),XFR(IMX),XFI(IMX),WR(IMX),WI(IMX)
      PI=3.14159265358979        
      FDPHI1=PI/36.
      FDPHI=FDPHI1/FLOAT(NESTRI(LEVEL))
      JS=-18
      PID=PI/36.
      FPHI=(NESTJ1(LEVEL)+JS-.5)*FDPHI1+FLOAT(J-1)*FDPHI
      KUT=0
      ARG0=14.0*PID
      IF(DABS(FPHI).LE.ARG0)THEN
        CALL N1SHAPIRO(NX,NKAX)
        RETURN
      ENDIF
      ARG1=14.0*PID
      ARG2=ARG1+PID
      ARG3=ARG2+PID
      ARG4=ARG3+PID
      IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
        KUTF=12
      ELSEIF(ABS(FPHI).GE.ARG2.AND.ABS(FPHI).LT.ARG3)THEN
        KUTF=9
      ELSEIF(ABS(FPHI).GE.ARG3.AND.ABS(FPHI).LT.ARG4)THEN
        KUTF=6
      ELSEIF(ABS(FPHI).GE.ARG4)THEN
        KUTF=1
      ENDIF      
c       KUT2=12
c       KUTF2=14
c       KUT3=8
c       KUTF3=10
c      DO M=1,18
c        ARG1=FLOAT(M-1)*PID
c        ARG2=FLOAT(M)*PID
c        IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
c          KUTF=24-2*M
c          KUTF=7
c        ENDIF
c      ENDDO
      CALL WTABLE(IMAX-2,WR,WI)
      NJNPXK=NJNP+NX
      DO K=1,NKAX
        DO I=1,IMAX-2
          XR(I)=F(I+1,NJNPXK)
          XI(I)=0.
        ENDDO
        CALL DFT1(IMAX-2,XR,XI,XFR,XFI,WR,WI,0)
        DO I=1,IMAX-2
c          IF(I.GE.KUT2.AND.I.LE.KUTF2)THEN
c            XFR(I)=0.5*XFR(I)
c            XFI(I)=0.5*XFI(I)
c          ELSEIF(I.GE.KUT3.AND.I.LE.KUTF3)THEN
c            XFR(I)=0.5*XFR(I)
c            XFI(I)=0.5*XFI(I)
c          ENDIF  
           IF(I.GE.KUTF)THEN
             XFR(I)=0.0*XFR(I)
             XFI(I)=0.0*XFI(I)
           ENDIF
        ENDDO
        CALL DFT1(IMAX-2,XFR,XFI,XR,XI,WR,WI,1)
        DO I=2,IMAX-1
          F(I,NJNPXK)=XR(I-1)
        ENDDO
        NJNPXK=NJNPXK+1
      ENDDO
      RETURN
      END
C       
*DK N1FILTER
      SUBROUTINE N1FILTER2(S,J,NKAX,LEVEL)
      INCLUDE 'n1param.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n1const.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'index.h'
      DIMENSION S(IMX,NKAX),XR(IMX),XI(IMX),XFR(IMX),XFI(IMX),
     1          WR(IMX),WI(IMX)
      PI=3.14159265358979        
      FDPHI1=PI/36.
      FDPHI=FDPHI1/FLOAT(NESTRI(LEVEL))
      JS=-18
      PID=PI/36.
      FPHI=(NESTJ1(LEVEL)+JS-.5)*FDPHI1+FLOAT(J-1)*FDPHI
      KUT=0
      ARG0=14.0*PID
      IF(ABS(FPHI).LE.ARG0)THEN
        RETURN
      ENDIF
      ARG1=14.0*PID
      ARG2=ARG1+PID
      ARG3=ARG2+PID
      ARG4=ARG3+PID            
      IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
        KUTF=12
      ELSEIF(ABS(FPHI).GE.ARG2.AND.ABS(FPHI).LT.ARG3)THEN
        KUTF=9
      ELSEIF(ABS(FPHI).GE.ARG3.AND.ABS(FPHI).LT.ARG4)THEN
        KUTF=6
      ELSEIF(ABS(FPHI).GE.ARG4)THEN
        KUTF=1
      ENDIF      
      CALL WTABLE(IMAX-2,WR,WI)
      DO K=1,NKAX
        DO I=1,IMAX-2
          XR(I)=S(I+1,K)
          XI(I)=0.
        ENDDO
        CALL DFT1(IMAX-2,XR,XI,XFR,XFI,WR,WI,0)
        DO I=1,IMAX-2
c          IF(I.GE.KUT2.AND.I.LE.KUTF2)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          ELSEIF(I.GE.KUT3.AND.I.LE.KUTF3)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          IF(I.GE.KUT3)THEN
c             XFR(I)=0.0*XFR(I)
c             XFI(I)=0.0*XFI(I)
c          ENDIF  
           IF(I.GE.KUTF)THEN
             XFR(I)=0.0*XFR(I)
             XFI(I)=0.0*XFI(I)
           ENDIF
        ENDDO
        CALL DFT1(IMAX-2,XFR,XFI,XR,XI,WR,WI,1)
        DO I=2,IMAX-1
          S(I,K)=XR(I-1)
        ENDDO
      ENDDO
      RETURN
      END
C
*DK N1FILTER
      SUBROUTINE N2FILTER(NX,J,NKAX,LEVEL)
      INCLUDE 'n2param.h'
      INCLUDE 'n2array.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n2const.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'index.h'
      DIMENSION XR(IMX),XI(IMX),XFR(IMX),XFI(IMX),WR(IMX),WI(IMX)
      PI=3.14159265358979        
      FDPHI1=PI/36.
      FDPHI=FDPHI1/FLOAT(NESTRI(LEVEL))
      JS=-18
      PID=PI/36.
      FPHI=(NESTJ1(LEVEL)+JS-.5)*FDPHI1+FLOAT(J-1)*FDPHI
      KUT=0
      ARG0=14.0*PID
      IF(DABS(FPHI).LE.ARG0)THEN
        CALL N1SHAPIRO(NX,NKAX)
        RETURN
      ENDIF
      ARG1=14.0*PID
      ARG2=ARG1+PID
      ARG3=ARG2+PID
      ARG4=ARG3+PID
      IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
        KUTF=12
      ELSEIF(ABS(FPHI).GE.ARG2.AND.ABS(FPHI).LT.ARG3)THEN
        KUTF=9
      ELSEIF(ABS(FPHI).GE.ARG3.AND.ABS(FPHI).LT.ARG4)THEN
        KUTF=6
      ELSEIF(ABS(FPHI).GE.ARG4)THEN
        KUTF=1
      ENDIF      
c       KUT2=12
c       KUTF2=14
c       KUT3=8
c       KUTF3=10
c      DO M=1,18
c        ARG1=FLOAT(M-1)*PID
c        ARG2=FLOAT(M)*PID
c        IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
c          KUTF=24-2*M
c          KUTF=7
c        ENDIF
c      ENDDO
      CALL WTABLE(IMAX-2,WR,WI)
      NJNPXK=NJNP+NX
      DO K=1,NKAX
        DO I=1,IMAX-2
          XR(I)=F(I+1,NJNPXK)
          XI(I)=0.
        ENDDO
        CALL DFT1(IMAX-2,XR,XI,XFR,XFI,WR,WI,0)
        DO I=1,IMAX-2
c          IF(I.GE.KUT2.AND.I.LE.KUTF2)THEN
c            XFR(I)=0.5*XFR(I)
c            XFI(I)=0.5*XFI(I)
c          ELSEIF(I.GE.KUT3.AND.I.LE.KUTF3)THEN
c            XFR(I)=0.5*XFR(I)
c            XFI(I)=0.5*XFI(I)
c          ENDIF  
           IF(I.GE.KUTF)THEN
             XFR(I)=0.0*XFR(I)
             XFI(I)=0.0*XFI(I)
           ENDIF
        ENDDO
        CALL DFT1(IMAX-2,XFR,XFI,XR,XI,WR,WI,1)
        DO I=2,IMAX-1
          F(I,NJNPXK)=XR(I-1)
        ENDDO
        NJNPXK=NJNPXK+1
      ENDDO
      RETURN
      END
C       
*DK N1FILTER
      SUBROUTINE N2FILTER2(S,J,NKAX,LEVEL)
      INCLUDE 'n2param.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n2const.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'index.h'
      DIMENSION S(IMX,NKAX),XR(IMX),XI(IMX),XFR(IMX),XFI(IMX),
     1          WR(IMX),WI(IMX)
      PI=3.14159265358979        
      FDPHI1=PI/36.
      FDPHI=FDPHI1/FLOAT(NESTRI(LEVEL))
      JS=-18
      PID=PI/36.
      FPHI=(NESTJ1(LEVEL)+JS-.5)*FDPHI1+FLOAT(J-1)*FDPHI
      KUT=0
      ARG0=14.0*PID
      IF(ABS(FPHI).LE.ARG0)THEN
        RETURN
      ENDIF
      ARG1=14.0*PID
      ARG2=ARG1+PID
      ARG3=ARG2+PID
      ARG4=ARG3+PID            
      IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
        KUTF=12
      ELSEIF(ABS(FPHI).GE.ARG2.AND.ABS(FPHI).LT.ARG3)THEN
        KUTF=9
      ELSEIF(ABS(FPHI).GE.ARG3.AND.ABS(FPHI).LT.ARG4)THEN
        KUTF=6
      ELSEIF(ABS(FPHI).GE.ARG4)THEN
        KUTF=1
      ENDIF      
      CALL WTABLE(IMAX-2,WR,WI)
      DO K=1,NKAX
        DO I=1,IMAX-2
          XR(I)=S(I+1,K)
          XI(I)=0.
        ENDDO
        CALL DFT1(IMAX-2,XR,XI,XFR,XFI,WR,WI,0)
        DO I=1,IMAX-2
c          IF(I.GE.KUT2.AND.I.LE.KUTF2)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          ELSEIF(I.GE.KUT3.AND.I.LE.KUTF3)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          IF(I.GE.KUT3)THEN
c             XFR(I)=0.0*XFR(I)
c             XFI(I)=0.0*XFI(I)
c          ENDIF  
           IF(I.GE.KUTF)THEN
             XFR(I)=0.0*XFR(I)
             XFI(I)=0.0*XFI(I)
           ENDIF
        ENDDO
        CALL DFT1(IMAX-2,XFR,XFI,XR,XI,WR,WI,1)
        DO I=2,IMAX-1
          S(I,K)=XR(I-1)
        ENDDO
      ENDDO
      RETURN
      END
C
*DK N1FILTER
      SUBROUTINE N3FILTER(S,J,NKAX,LEVEL)
      INCLUDE 'n3param.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n3const.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'index.h'
      DIMENSION S(IMX,NKAX),XR(IMX),XI(IMX),XFR(IMX),XFI(IMX),
     1          WR(IMX),WI(IMX)
      PI=3.14159265358979        
      FDPHI1=PI/36.
      FDPHI=FDPHI1/FLOAT(NESTRI(LEVEL))
      JS=-18
      PID=PI/36.
      FPHI=(NESTJ1(LEVEL)+JS-.5)*FDPHI1+FLOAT(J-1)*FDPHI
      KUT=0
      ARG0=14.0*PID
      IF(ABS(FPHI).LE.ARG0)THEN
        RETURN
      ENDIF
      ARG1=14.0*PID
      ARG2=ARG1+PID
      ARG3=ARG2+PID
      ARG4=ARG3+PID            
      IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
        KUTF=12
      ELSEIF(ABS(FPHI).GE.ARG2.AND.ABS(FPHI).LT.ARG3)THEN
        KUTF=9
      ELSEIF(ABS(FPHI).GE.ARG3.AND.ABS(FPHI).LT.ARG4)THEN
        KUTF=6
      ELSEIF(ABS(FPHI).GE.ARG4)THEN
        KUTF=1
      ENDIF      
      CALL WTABLE(IMAX-2,WR,WI)
      DO K=1,NKAX
        DO I=1,IMAX-2
          XR(I)=S(I+1,K)
          XI(I)=0.
        ENDDO
        CALL DFT1(IMAX-2,XR,XI,XFR,XFI,WR,WI,0)
        DO I=1,IMAX-2
c          IF(I.GE.KUT2.AND.I.LE.KUTF2)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          ELSEIF(I.GE.KUT3.AND.I.LE.KUTF3)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          IF(I.GE.KUT3)THEN
c             XFR(I)=0.0*XFR(I)
c             XFI(I)=0.0*XFI(I)
c          ENDIF  
           IF(I.GE.KUTF)THEN
             XFR(I)=0.0*XFR(I)
             XFI(I)=0.0*XFI(I)
           ENDIF
        ENDDO
        CALL DFT1(IMAX-2,XFR,XFI,XR,XI,WR,WI,1)
        DO I=2,IMAX-1
          S(I,K)=XR(I-1)
        ENDDO
      ENDDO
      RETURN
      END
C
*DK N1FILTER
      SUBROUTINE N4FILTER(S,J,NKAX,LEVEL)
      INCLUDE 'n4param.h'
      INCLUDE 'buffe.h'
      INCLUDE 'n4const.h'
      INCLUDE 'nestcnt.h'
      INCLUDE 'index.h'
      DIMENSION S(IMX,NKAX),XR(IMX),XI(IMX),XFR(IMX),XFI(IMX),
     1          WR(IMX),WI(IMX)
      PI=3.14159265358979        
      FDPHI1=PI/36.
      FDPHI=FDPHI1/FLOAT(NESTRI(LEVEL))
      JS=-18
      PID=PI/36.
      FPHI=(NESTJ1(LEVEL)+JS-.5)*FDPHI1+FLOAT(J-1)*FDPHI
      KUT=0
      ARG0=14.0*PID
      IF(ABS(FPHI).LE.ARG0)THEN
        RETURN
      ENDIF
      ARG1=14.0*PID
      ARG2=ARG1+PID
      ARG3=ARG2+PID
      ARG4=ARG3+PID            
      IF(ABS(FPHI).GE.ARG1.AND.ABS(FPHI).LT.ARG2)THEN
        KUTF=12
      ELSEIF(ABS(FPHI).GE.ARG2.AND.ABS(FPHI).LT.ARG3)THEN
        KUTF=9
      ELSEIF(ABS(FPHI).GE.ARG3.AND.ABS(FPHI).LT.ARG4)THEN
        KUTF=6
      ELSEIF(ABS(FPHI).GE.ARG4)THEN
        KUTF=1
      ENDIF      
      CALL WTABLE(IMAX-2,WR,WI)
      DO K=1,NKAX
        DO I=1,IMAX-2
          XR(I)=S(I+1,K)
          XI(I)=0.
        ENDDO
        CALL DFT1(IMAX-2,XR,XI,XFR,XFI,WR,WI,0)
        DO I=1,IMAX-2
c          IF(I.GE.KUT2.AND.I.LE.KUTF2)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          ELSEIF(I.GE.KUT3.AND.I.LE.KUTF3)THEN
c            XFR(I)=0.0*XFR(I)
c            XFI(I)=0.0*XFI(I)
c          IF(I.GE.KUT3)THEN
c             XFR(I)=0.0*XFR(I)
c             XFI(I)=0.0*XFI(I)
c          ENDIF  
           IF(I.GE.KUTF)THEN
             XFR(I)=0.0*XFR(I)
             XFI(I)=0.0*XFI(I)
           ENDIF
        ENDDO
        CALL DFT1(IMAX-2,XFR,XFI,XR,XI,WR,WI,1)
        DO I=2,IMAX-1
          S(I,K)=XR(I-1)
        ENDDO
      ENDDO
      RETURN
      END
C
*DK DFT
      SUBROUTINE DFT(N,XR,XI,XFR,XFI,ISEL)
      IMPLICIT REAL*8 (A-H,O-Z)
C     ****
C     ****      ISEL=0 : DFT
C     ****      ISEL=1 : INVERSE DFT
C     ****
      DIMENSION XR(N),XI(N),XFR(N),XFI(N)
      PI=3.1415927
      WN=2.*PI/FLOAT(N)
      IF(ISEL.EQ.1)THEN
        WN=-WN
      ENDIF
      DO K=1,N
        XFR(K)=0.
        XFI(K)=0.
        KM1=K-1
        DO I=1,N
          IM1=I-1
          ARG=WN*FLOAT(KM1*IM1)
          C=COS(ARG)
          S=SIN(ARG)
          XFR(K)=XFR(K)+XR(I)*C+XI(I)*S
          XFI(K)=XFI(K)-XR(I)*S+XI(I)*C
        ENDDO
        IF(ISEL.EQ.1)THEN
          XFR(K)=XFR(K)/FLOAT(N)
          XFI(k)=XFI(K)/FLOAT(N)
        ENDIF
      ENDDO
      RETURN
      END
C
*DK DFT
      SUBROUTINE DFT1(N,XR,XI,XFR,XFI,WR,WI,ISEL)
      IMPLICIT REAL*8 (A-H,O-Z)
C     ****
C     ****      ISEL=0 : DFT
C     ****      ISEL=1 : INVERSE DFT
C     ****
      DIMENSION XR(N),XI(N),XFR(N),XFI(N),WR(N),WI(N)
      DO K=1,N
        XFR(K)=XR(1)
        XFI(K)=XI(1)
        L=1
        DO I=2,N
          L=L+K-1
          IF(L.GT.N)THEN
            L=L-N
          ENDIF
          C=WR(L)
          IF(ISEL.EQ.1)THEN
            S=-WI(L)
          ELSE
            S=WI(L)
          ENDIF
          XFR(K)=XFR(K)+XR(I)*C+XI(I)*S
          XFI(K)=XFI(K)-XR(I)*S+XI(I)*C
        ENDDO
        IF(ISEL.EQ.1)THEN
          XFR(K)=XFR(K)/FLOAT(N)
          XFI(k)=XFI(K)/FLOAT(N)
        ENDIF
      ENDDO
      RETURN
      END
C
*DK WTABLE
      SUBROUTINE WTABLE(N,WR,WI)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION WR(N),WI(N)
      PI=3.1415927
      WN=2.*PI/FLOAT(N)
      DO K=1,N
        PHI=WN*FLOAT(K-1)      
        WR(K)=COS(PHI)
        WI(K)=SIN(PHI)
      ENDDO
      RETURN
      END
C
*DK SHAPIRO
      SUBROUTINE N1SHAPIRO(NX,KX)
      include 'n1param.h'                                                 SMOOTH.3
      include 'n1array.h'                                                 SMOOTH.4
      include 'n1const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                 SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMX+1)
      C26=0.125/4
      NXNPK=NJNP+NX
      goto 510
      DO K=1,KX
        DO I=4,IMAX-3
          S(I,K)=F(I,NXNPK)+C26*(F(I+3,NXNPK)+F(I-3,NXNPK)-6.*
     1               (F(I+2,NXNPK)+F(I-2,NXNPK))+15.*(F(I+1,NXNPK)+
     2               F(I-1,NXNPK))-20.*F(I,NXNPK))
        ENDDO
        S(3,K)=F(3,NXNPK)-C26*(F(5,NXNPK)+F(1,NXNPK)-4.*(
     1             F(4,NXNPK)+F(2,NXNPK))+6.*F(3,NXNPK))
        S(LEN1-2,K)=F(LEN1-2,NXNPK)-C26*(F(LEN1,NXNPK)+
     1                  F(LEN1-4,NXNPK)-4.*(F(LEN1-1,NXNPK)+
     2                  F(LEN1-3,NXNPK))+6.*F(LEN1-2,NXNPK))        
        S(2,K)=F(2,NXNPK)-C26*(-4.*(F(1,NXNPK)+F(3,NXNPK))+8.*
     1             F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C26*(-4.*(F(LEN1-2,NXNPK)+
     1                 F(LEN1,NXNPK))+8.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
510   CONTINUE
c      goto 511
      DO K=1,KX
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C26*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C26*(-F(1,NXNPK)-F(3,NXNPK)+2.*
     1             F(2,NXNPK))*4.
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-4.*C26*(-F(LEN1-2,NXNPK)-
     1              F(LEN1,NXNPK)+2.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
511   CONTINUE
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C      
*DK SHAPIRO
      SUBROUTINE N2SHAPIRO(NX,KX)
      include 'n2param.h'                                                 SMOOTH.3
      include 'n2array.h'                                                 SMOOTH.4
      include 'n2const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                 SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMX+1)
      C26=0.125/4
      NXNPK=NJNP+NX
      goto 510
      write(6,*)KX
      DO K=1,KX
        DO I=4,IMAX-3
          S(I,K)=F(I,NXNPK)+C26*(F(I+3,NXNPK)+F(I-3,NXNPK)-6.*
     1               (F(I+2,NXNPK)+F(I-2,NXNPK))+15.*(F(I+1,NXNPK)+
     2               F(I-1,NXNPK))-20.*F(I,NXNPK))
        ENDDO
        S(3,K)=F(3,NXNPK)-C26*(F(5,NXNPK)+F(1,NXNPK)-4.*(
     1             F(4,NXNPK)+F(2,NXNPK))+6.*F(3,NXNPK))
        S(LEN1-2,K)=F(LEN1-2,NXNPK)-C26*(F(LEN1,NXNPK)+
     1                  F(LEN1-4,NXNPK)-4.*(F(LEN1-1,NXNPK)+
     2                  F(LEN1-3,NXNPK))+6.*F(LEN1-2,NXNPK))        
        S(2,K)=F(2,NXNPK)-C26*(-4.*(F(1,NXNPK)+F(3,NXNPK))+8.*
     1             F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C26*(-4.*(F(LEN1-2,NXNPK)+
     1                 F(LEN1,NXNPK))+8.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
510   CONTINUE
c      goto 511
      DO K=1,KX
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C26*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C26*(-F(1,NXNPK)-F(3,NXNPK)+2.*
     1             F(2,NXNPK))*4.
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-4.*C26*(-F(LEN1-2,NXNPK)-
     1              F(LEN1,NXNPK)+2.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
511   CONTINUE
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C      
*DK SHAPIRO
      SUBROUTINE N3SHAPIRO(NX,KX)
      include 'n3param.h'                                                 SMOOTH.3
      include 'n3array.h'                                                 SMOOTH.4
      include 'n3const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                 SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMX+1)
      C26=0.125/4
      NXNPK=NJNP+NX
      goto 510
      write(6,*)KX
      DO K=1,KX
        DO I=4,IMAX-3
          S(I,K)=F(I,NXNPK)+C26*(F(I+3,NXNPK)+F(I-3,NXNPK)-6.*
     1               (F(I+2,NXNPK)+F(I-2,NXNPK))+15.*(F(I+1,NXNPK)+
     2               F(I-1,NXNPK))-20.*F(I,NXNPK))
        ENDDO
        S(3,K)=F(3,NXNPK)-C26*(F(5,NXNPK)+F(1,NXNPK)-4.*(
     1             F(4,NXNPK)+F(2,NXNPK))+6.*F(3,NXNPK))
        S(LEN1-2,K)=F(LEN1-2,NXNPK)-C26*(F(LEN1,NXNPK)+
     1                  F(LEN1-4,NXNPK)-4.*(F(LEN1-1,NXNPK)+
     2                  F(LEN1-3,NXNPK))+6.*F(LEN1-2,NXNPK))        
        S(2,K)=F(2,NXNPK)-C26*(-4.*(F(1,NXNPK)+F(3,NXNPK))+8.*
     1             F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C26*(-4.*(F(LEN1-2,NXNPK)+
     1                 F(LEN1,NXNPK))+8.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
510   CONTINUE
c      goto 511
      DO K=1,KX
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C26*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C26*(-F(1,NXNPK)-F(3,NXNPK)+2.*
     1             F(2,NXNPK))*4.
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-4.*C26*(-F(LEN1-2,NXNPK)-
     1              F(LEN1,NXNPK)+2.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
511   CONTINUE
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C      
*DK SHAPIRO
      SUBROUTINE N4SHAPIRO(NX,KX)
      include 'n4param.h'                                                 SMOOTH.3
      include 'n4array.h'                                                 SMOOTH.4
      include 'n4const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                 SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMX+1)
      C26=0.125/4
      NXNPK=NJNP+NX
      goto 510
      write(6,*)KX
      DO K=1,KX
        DO I=4,IMAX-3
          S(I,K)=F(I,NXNPK)+C26*(F(I+3,NXNPK)+F(I-3,NXNPK)-6.*
     1               (F(I+2,NXNPK)+F(I-2,NXNPK))+15.*(F(I+1,NXNPK)+
     2               F(I-1,NXNPK))-20.*F(I,NXNPK))
        ENDDO
        S(3,K)=F(3,NXNPK)-C26*(F(5,NXNPK)+F(1,NXNPK)-4.*(
     1             F(4,NXNPK)+F(2,NXNPK))+6.*F(3,NXNPK))
        S(LEN1-2,K)=F(LEN1-2,NXNPK)-C26*(F(LEN1,NXNPK)+
     1                  F(LEN1-4,NXNPK)-4.*(F(LEN1-1,NXNPK)+
     2                  F(LEN1-3,NXNPK))+6.*F(LEN1-2,NXNPK))        
        S(2,K)=F(2,NXNPK)-C26*(-4.*(F(1,NXNPK)+F(3,NXNPK))+8.*
     1             F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C26*(-4.*(F(LEN1-2,NXNPK)+
     1                 F(LEN1,NXNPK))+8.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
510   CONTINUE
c      goto 511
      DO K=1,KX
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C26*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C26*(-F(1,NXNPK)-F(3,NXNPK)+2.*
     1             F(2,NXNPK))*4.
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-4.*C26*(-F(LEN1-2,NXNPK)-
     1              F(LEN1,NXNPK)+2.*F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
511   CONTINUE
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C      
*DK SHAPIRO
      SUBROUTINE N1SHAPIRO2(NX,KX)
      include 'n1param.h'                                                 SMOOTH.3
      include 'n1array.h'                                                 SMOOTH.4
      include 'n1const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                    SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMXP)
      C1=1./256.
      C2=1./16.
      C3=4.*C2
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=3,IMAX-2
         S(I,K)=F(I,NXNPK)-C2*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1          (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C2*(2.*F(4,NXNPK)-4.*(F(1,NXNPK)+F(3,NXNPK))
     1         +6.*F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C2*(2.*F(LEN1-3,NXNPK)
     1              -4.*(F(LEN1-2,NXNPK)+F(LEN1,NXNPK))+6.*
     1              F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C
*DK SHAPIRO
      SUBROUTINE N2SHAPIRO2(NX,KX)
      include 'n2param.h'                                                 SMOOTH.3
      include 'n2array.h'                                                 SMOOTH.4
      include 'n2const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMXP)
      C1=1./256.
      C2=1./16.
      C3=4.*C2
      NXNPK=NJNP+NX
      DO K=1,KX
c        DO I=5,IMAX-4
c          S(I,K)=F(I,NXNPK)+C1*(-(F(I+4,NXNPK)+F(I-4,NXNPK))+8.*
c     1               (F(I+3,NXNPK)+F(I-3,NXNPK))-28.*(F(I+2,NXNPK)+
c     2               F(I-2,NXNPK))+56.*(F(I+1,NXNPK)+F(I-1,NXNPK))-
c     3               70.*F(I,NXNPK))
c        ENDDO
c        S(4,K)=F(4,NXNPK)+C2*(-(F(6,NXNPK)+F(2,NXNPK))+4.*(
c     1             F(5,NXNPK)+F(3,NXNPK))-6.*F(4,NXNPK))
c        S(3,K)=F(3,NXNPK)+C2*(-(F(5,NXNPK)+F(1,NXNPK))+4.*(
c     1             F(4,NXNPK)+F(2,NXNPK))-6.*F(3,NXNPK))
c        S(LEN1-3,K)=F(LEN1-3,NXNPK)+C2*(-(F(LEN1-1,NXNPK)+
c     1                  F(LEN1-5,NXNPK))+4.*(F(LEN1-2,NXNPK)+
c     2                  F(LEN1-4,NXNPK))-6.*F(LEN1-1,NXNPK))
c        S(LEN1-2,K)=F(LEN1-2,NXNPK)+C2*(-(F(LEN1,NXNPK)+
c     1                  F(LEN1-4,NXNPK))+4.*(F(LEN1-1,NXNPK)+
c     2                  F(LEN1-3,NXNPK))-6.*F(LEN1-2,NXNPK))
c        
c        S(2,K)=F(2,NXNPK)+C3*((F(1,NXNPK)+F(3,NXNPK))-2.*
c     1             F(2,NXNPK))
c        S(LEN1-1,K)=F(LEN1-1,NXNPK)+C3*((F(LEN1-2,NXNPK)+
c     1                 F(LEN1,NXNPK))-2.*F(LEN1-1,NXNPK))
c        do i=2,len1-1
c          s(i,k)=f(i,nxnpk)+c3*(f(i+1,nxnpk)+f(i-1,nxnpk)-2.*
c     1           f(i,nxnpk))
c        enddo 
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C2*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C2*(2.*F(4,NXNPK)-4.*(F(1,NXNPK)+F(3,NXNPK))
     1         +6.*F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C2*(2.*F(LEN1-3,NXNPK)
     1              -4.*(F(LEN1-2,NXNPK)+F(LEN1,NXNPK))+6.*
     1              F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C
*DK SHAPIRO
      SUBROUTINE N3SHAPIRO2(NX,KX)
      include 'n3param.h'                                                 SMOOTH.3
      include 'n3array.h'                                                 SMOOTH.4
      include 'n1const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMXP)
      C1=1./256.
      C2=1./16.
      C3=4.*C2
      NXNPK=NJNP+NX
      DO K=1,KX
c        DO I=5,IMAX-4
c          S(I,K)=F(I,NXNPK)+C1*(-(F(I+4,NXNPK)+F(I-4,NXNPK))+8.*
c     1               (F(I+3,NXNPK)+F(I-3,NXNPK))-28.*(F(I+2,NXNPK)+
c     2               F(I-2,NXNPK))+56.*(F(I+1,NXNPK)+F(I-1,NXNPK))-
c     3               70.*F(I,NXNPK))
c        ENDDO
c        S(4,K)=F(4,NXNPK)+C2*(-(F(6,NXNPK)+F(2,NXNPK))+4.*(
c     1             F(5,NXNPK)+F(3,NXNPK))-6.*F(4,NXNPK))
c        S(3,K)=F(3,NXNPK)+C2*(-(F(5,NXNPK)+F(1,NXNPK))+4.*(
c     1             F(4,NXNPK)+F(2,NXNPK))-6.*F(3,NXNPK))
c        S(LEN1-3,K)=F(LEN1-3,NXNPK)+C2*(-(F(LEN1-1,NXNPK)+
c     1                  F(LEN1-5,NXNPK))+4.*(F(LEN1-2,NXNPK)+
c     2                  F(LEN1-4,NXNPK))-6.*F(LEN1-1,NXNPK))
c        S(LEN1-2,K)=F(LEN1-2,NXNPK)+C2*(-(F(LEN1,NXNPK)+
c     1                  F(LEN1-4,NXNPK))+4.*(F(LEN1-1,NXNPK)+
c     2                  F(LEN1-3,NXNPK))-6.*F(LEN1-2,NXNPK))
c        
c        S(2,K)=F(2,NXNPK)+C3*((F(1,NXNPK)+F(3,NXNPK))-2.*
c     1             F(2,NXNPK))
c        S(LEN1-1,K)=F(LEN1-1,NXNPK)+C3*((F(LEN1-2,NXNPK)+
c     1                 F(LEN1,NXNPK))-2.*F(LEN1-1,NXNPK))
c        do i=2,len1-1
c          s(i,k)=f(i,nxnpk)+c3*(f(i+1,nxnpk)+f(i-1,nxnpk)-2.*
c     1           f(i,nxnpk))
c        enddo 
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C2*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C2*(2.*F(4,NXNPK)-4.*(F(1,NXNPK)+F(3,NXNPK))
     1         +6.*F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C2*(2.*F(LEN1-3,NXNPK)
     1              -4.*(F(LEN1-2,NXNPK)+F(LEN1,NXNPK))+6.*
     1              F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C
*DK SHAPIRO
      SUBROUTINE N4SHAPIRO2(NX,KX)
      include 'n4param.h'                                                 SMOOTH.3
      include 'n4array.h'                                                 SMOOTH.4
      include 'n4const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      DIMENSION S(IMX,KMXP)
      C1=1./256.
      C2=1./16.
      C3=4.*C2
      NXNPK=NJNP+NX
      DO K=1,KX
c        DO I=5,IMAX-4
c          S(I,K)=F(I,NXNPK)+C1*(-(F(I+4,NXNPK)+F(I-4,NXNPK))+8.*
c     1               (F(I+3,NXNPK)+F(I-3,NXNPK))-28.*(F(I+2,NXNPK)+
c     2               F(I-2,NXNPK))+56.*(F(I+1,NXNPK)+F(I-1,NXNPK))-
c     3               70.*F(I,NXNPK))
c        ENDDO
c        S(4,K)=F(4,NXNPK)+C2*(-(F(6,NXNPK)+F(2,NXNPK))+4.*(
c     1             F(5,NXNPK)+F(3,NXNPK))-6.*F(4,NXNPK))
c        S(3,K)=F(3,NXNPK)+C2*(-(F(5,NXNPK)+F(1,NXNPK))+4.*(
c     1             F(4,NXNPK)+F(2,NXNPK))-6.*F(3,NXNPK))
c        S(LEN1-3,K)=F(LEN1-3,NXNPK)+C2*(-(F(LEN1-1,NXNPK)+
c     1                  F(LEN1-5,NXNPK))+4.*(F(LEN1-2,NXNPK)+
c     2                  F(LEN1-4,NXNPK))-6.*F(LEN1-1,NXNPK))
c        S(LEN1-2,K)=F(LEN1-2,NXNPK)+C2*(-(F(LEN1,NXNPK)+
c     1                  F(LEN1-4,NXNPK))+4.*(F(LEN1-1,NXNPK)+
c     2                  F(LEN1-3,NXNPK))-6.*F(LEN1-2,NXNPK))
c        
c        S(2,K)=F(2,NXNPK)+C3*((F(1,NXNPK)+F(3,NXNPK))-2.*
c     1             F(2,NXNPK))
c        S(LEN1-1,K)=F(LEN1-1,NXNPK)+C3*((F(LEN1-2,NXNPK)+
c     1                 F(LEN1,NXNPK))-2.*F(LEN1-1,NXNPK))
c        do i=2,len1-1
c          s(i,k)=f(i,nxnpk)+c3*(f(i+1,nxnpk)+f(i-1,nxnpk)-2.*
c     1           f(i,nxnpk))
c        enddo 
        DO I=3,IMAX-2
          S(I,K)=F(I,NXNPK)-C2*(F(I+2,NXNPK)+F(I-2,NXNPK)-4.*
     1               (F(I+1,NXNPK)+F(I-1,NXNPK))+6.*F(I,NXNPK))
        ENDDO
        S(2,K)=F(2,NXNPK)-C2*(2.*F(4,NXNPK)-4.*(F(1,NXNPK)+F(3,NXNPK))
     1         +6.*F(2,NXNPK))
        S(LEN1-1,K)=F(LEN1-1,NXNPK)-C2*(2.*F(LEN1-3,NXNPK)
     1              -4.*(F(LEN1-2,NXNPK)+F(LEN1,NXNPK))+6.*
     1              F(LEN1-1,NXNPK))
        NXNPK=NXNPK+1
      ENDDO
      NXNPK=NJNP+NX
      DO K=1,KX
        DO I=2,LEN1-1
          F(I,NXNPK)=S(I,K)
        ENDDO
        NXNPK=NXNPK+1
      ENDDO
      RETURN
      END
C
*DK SMOOTH                                                                SMOOTH.1
      SUBROUTINE N1SMOOTH(NX,NXNM)                                        SMOOTH.2
      include 'n1param.h'                                                 SMOOTH.3
      include 'n1array.h'                                                 SMOOTH.4
      include 'n1const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      include 'n1cvscr.h'                                                 SMOOTH.8
C     ****                                                                SMOOTH.9
C     ****     TIME SMOOTHING PORCESS FOR BOUNDARY VARIABLES              SMOOTH.10
C     ****                                                                SMOOTH.11
      NXNMK=NJ+NXNM                                                       SMOOTH.12
      NXK=NJ+NX                                                           SMOOTH.13
      NXNPK=NJNP+NX                                                       SMOOTH.14
      NXNMNK=NJNP+NXNM                                                    SMOOTH.15
      DO K=1,KMAXP1                                                       SMOOTH.16
        DO I=1,LEN1                                                       SMOOTH.17
           S1(I,K)=(F(I,NXNMK)+F(I,NXNPK))*C(31)                          SMOOTH.18
           F(I,NXNMNK)=F(I,NXK)*C(30)+S1(I,K)                             SMOOTH.19
        ENDDO                                                             SMOOTH.20
        NXNMK=NXNMK+1                                                     SMOOTH.21
        NXK=NXK+1                                                         SMOOTH.22
        NXNPK=NXNPK+1                                                     SMOOTH.23
        NXNMNK=NXNMNK+1                                                   SMOOTH.24
      ENDDO                                                               SMOOTH.25
      RETURN                                                              SMOOTH.26
      END                                                                 SMOOTH.27
C                                                                         SMOOTH.28
*DK SMOOTH                                                                SMOOTH.1
      SUBROUTINE N2SMOOTH(NX,NXNM)                                        SMOOTH.2
      include 'n2param.h'                                                 SMOOTH.3
      include 'n2array.h'                                                 SMOOTH.4
      include 'n2const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      include 'n2cvscr.h'                                                 SMOOTH.8
C     ****                                                                SMOOTH.9
C     ****     TIME SMOOTHING PORCESS FOR BOUNDARY VARIABLES              SMOOTH.10
C     ****                                                                SMOOTH.11
      NXNMK=NJ+NXNM                                                       SMOOTH.12
      NXK=NJ+NX                                                           SMOOTH.13
      NXNPK=NJNP+NX                                                       SMOOTH.14
      NXNMNK=NJNP+NXNM                                                    SMOOTH.15
      DO K=1,KMAXP1                                                       SMOOTH.16
        DO I=1,LEN1                                                       SMOOTH.17
           S1(I,K)=(F(I,NXNMK)+F(I,NXNPK))*C(31)                          SMOOTH.18
           F(I,NXNMNK)=F(I,NXK)*C(30)+S1(I,K)                             SMOOTH.19
        ENDDO                                                             SMOOTH.20
        NXNMK=NXNMK+1                                                     SMOOTH.21
        NXK=NXK+1                                                         SMOOTH.22
        NXNPK=NXNPK+1                                                     SMOOTH.23
        NXNMNK=NXNMNK+1                                                   SMOOTH.24
      ENDDO                                                               SMOOTH.25
      RETURN                                                              SMOOTH.26
      END                                                                 SMOOTH.27
C                                                                         SMOOTH.28
*DK SMOOTH                                                                SMOOTH.1
      SUBROUTINE N3SMOOTH(NX,NXNM)                                        SMOOTH.2
      include 'n3param.h'                                                 SMOOTH.3
      include 'n3array.h'                                                 SMOOTH.4
      include 'n3const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      include 'n3cvscr.h'                                                 SMOOTH.8
C     ****                                                                SMOOTH.9
C     ****     TIME SMOOTHING PORCESS FOR BOUNDARY VARIABLES              SMOOTH.10
C     ****                                                                SMOOTH.11
      NXNMK=NJ+NXNM                                                       SMOOTH.12
      NXK=NJ+NX                                                           SMOOTH.13
      NXNPK=NJNP+NX                                                       SMOOTH.14
      NXNMNK=NJNP+NXNM                                                    SMOOTH.15
      DO K=1,KMAXP1                                                       SMOOTH.16
        DO I=1,LEN1                                                       SMOOTH.17
           S1(I,K)=(F(I,NXNMK)+F(I,NXNPK))*C(31)                          SMOOTH.18
           F(I,NXNMNK)=F(I,NXK)*C(30)+S1(I,K)                             SMOOTH.19
        ENDDO                                                             SMOOTH.20
        NXNMK=NXNMK+1                                                     SMOOTH.21
        NXK=NXK+1                                                         SMOOTH.22
        NXNPK=NXNPK+1                                                     SMOOTH.23
        NXNMNK=NXNMNK+1                                                   SMOOTH.24
      ENDDO                                                               SMOOTH.25
      RETURN                                                              SMOOTH.26
      END                                                                 SMOOTH.27
C                                                                         SMOOTH.28
*DK SMOOTH                                                                SMOOTH.1
      SUBROUTINE N4SMOOTH(NX,NXNM)                                        SMOOTH.2
      include 'n4param.h'                                                 SMOOTH.3
      include 'n4array.h'                                                 SMOOTH.4
      include 'n4const.h'                                                 SMOOTH.5
      include 'buffe.h'                                                   SMOOTH.6
      include 'index.h'                                                   SMOOTH.7
      include 'n4cvscr.h'                                                 SMOOTH.8
C     ****                                                                SMOOTH.9
C     ****     TIME SMOOTHING PORCESS FOR BOUNDARY VARIABLES              SMOOTH.10
C     ****                                                                SMOOTH.11
      NXNMK=NJ+NXNM                                                       SMOOTH.12
      NXK=NJ+NX                                                           SMOOTH.13
      NXNPK=NJNP+NX                                                       SMOOTH.14
      NXNMNK=NJNP+NXNM                                                    SMOOTH.15
      DO K=1,KMAXP1                                                       SMOOTH.16
        DO I=1,LEN1                                                       SMOOTH.17
           S1(I,K)=(F(I,NXNMK)+F(I,NXNPK))*C(31)                          SMOOTH.18
           F(I,NXNMNK)=F(I,NXK)*C(30)+S1(I,K)                             SMOOTH.19
        ENDDO                                                             SMOOTH.20
        NXNMK=NXNMK+1                                                     SMOOTH.21
        NXK=NXK+1                                                         SMOOTH.22
        NXNPK=NXNPK+1                                                     SMOOTH.23
        NXNMNK=NXNMNK+1                                                   SMOOTH.24
      ENDDO                                                               SMOOTH.25
      RETURN                                                              SMOOTH.26
      END                                                                 SMOOTH.27
C                                                                         SMOOTH.28
