C                         WD23JP
C                                                         *****************
C                                                         *   OCT.FOR     *
C                                                         *  PURSER 1995  *
C                                                         *****************
C

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
C                   SUBROUTINE INOCT                                           C
C   Initialize the constants of common/CSTOCT/ used to describe grid and to    C
C   perform transformations between map coordinates and the sphere or disk     C
C                                                                              C
C --> N1GRID    number of grid-spaces from center to side of octagon           C
C --> NAGRID    number of grid-spaces from center-line to nearest vertex       C
C --> NTAY      number of Taylor series coefficients retained (NTAY<41)        C
C --> FLON0     longitude of center of octagon-1                               C
C --> FLAT0     latitude of center of octagon-1                                C
C --> SCEN      scale-enhancement parameter = radius of image of octagon-1 in  C
C               a stereographic projection centered on (FLAT0,FLON0) that maps C
C               the concentric hemisphere to the unit disk                     C
C --- RA,QA     work arrays of dimension N+1 (N defined below)                 C
C --> N         number of the form 2**m (eg 128) and > 3*NTAY defining the     C
C               period of FFT used in the construction of Taylor series        C
C --- JUMBLE,WORK   work arrays of dimension N                                 C
C------------------------------------------------------------------------------C
      SUBROUTINE INOCT(N1GRID,NAGRID,NTAY,FLON0,FLAT0,SCEN,RA,QA,N
     *                ,JUMBLE,WORK)
      COMPLEX CI,CIR,CIQ,CIA,CIAP,CIRS2,Z,W,CANG,CDANG1,CDANG2
      COMMON/CSTOCT/A1(40),A2(40),B1(40),B2(40),M1G,MAG,A,B,CX,CY
     *,SSR,N1,ROTM(3,3),SC,SCI,ST,STI,CI,CIR,CIQ
      DIMENSION RA(0:*),QA(0:*),JUMBLE(N),WORK(N)

      IF(3*N1.GT.N)STOP'N TOO SMALL IN INOCT'
      M1G=N1GRID
      MAG=NAGRID
      A=FLOAT(MAG)/M1G
      N1=NTAY
      CALL ININMAP(FLON0,FLAT0,ROTM)
      SC=SCEN
      SCI=1./SC
      ST=4./3.
      STI=3./4.
      CI=CMPLX(0.,1.)
      CIR=(-CI)**.75
      CIQ=CI**.25
      NH=N/2
      R2=SQRT(2.)
      ORC=.5
      SW1=SQRT(1.+A*A)
      SW2=AMIN1(A*2,R2*(1.-A))
      SSR=(SW1/SW2)**2
      RS1=.95*SW1
      RS2=.95*SW2
      RS3=RS1**4
      RS4=RS2**ST

      DO I=1,N1
       A1(I)=0.
      ENDDO
      A1(1)= 1.05060
      A1(2)= -.172629
      A1(3)=  .125716
      A1(4)=  .00015657
      A1(5)= -.0017866
      A1(6)= -.0031168

      DO I=1,N1
       A2(I)=0.
      ENDDO
      A2(1)= .688411
      A2(2)=-.180924
      A2(3)=-.186593
      A2(4)= .024811
      A2(5)=-.044888
      A2(6)=-.049190
      B=-.147179

      JUMBLE(1)=0
      PIO4=ATAN(1.)
      CDANG1=CI*PIO4/NH
      CDANG2=-3.*CDANG1
      CIA=CI*A
      CIAP=CIA+1.
      CIRS2=CI*RS2

      DO KIT=1,50
!          Print *,'JUMBLE(1)=',JUMBLE(1)
!         If(kit.eq.0) JUMBLE(1)=0
ctest
!         Print *, 'kit=',kit
ctest
       DO I=0,NH
        CANG=CDANG1*I
        Z=RS1*CEXP(CANG)
        CALL TOCT(Z,W)
        W=W**4
        U=REAL(W)
        V=AIMAG(W)
        RA(I)=U
        QA(I)=-V
        RA(N-I)=U
        QA(N-I)=V
       ENDDO

       QA(0)=0.
       QA(NH)=0.
ctest
!         Print *, 'Before first CFFT'
!         Print *, 'JUMBLE(1)=',JUMBLE(1)
!         Print *, 'n=',n
ctest
       CALL CFFT(RA,QA,N,1.,WORK,JUMBLE)
         RS3P=1.
       DO I=1,N1
         RS3P=RS3P*RS3
        A1(I)=A1(I) +ORC*(RA(I)/RS3P-A1(I))
       ENDDO

      DO I=0,NH
       CANG=CDANG2*I
       Z=CIAP-CIRS2*CEXP(CANG)
       CALL TOCT(Z,W)
       W=CI*(W-1.)/(W+1.)
       U=REAL(W)
       V=AIMAG(W)
       RA(I)=U
       QA(I)=V
       RA(N-I)=U
       QA(N-I)=-V
      ENDDO
      QA(0)=0.
      QA(NH)=0.
ctest
!         Print *, 'Before second FFT'
!         Print *, 'JUMBLE(1)=',JUMBLE(1)
!         Print *, 'n=',n
ctest
      CALL CFFT(RA,QA,N,1.,WORK,JUMBLE)
      B=B+ORC*(RA(0)-B)
       RS4P=1.
      DO I=1,N1
       RS4P=RS4P*RS4
       A2(I)=A2(I) +ORC*(RA(I)/RS4P-A2(I))
      ENDDO

      ENDDO
      W=B
      W=(W+CI)/(CI-W)
      CX=REAL(W)
      CY=AIMAG(W)
      CALL CINVRT(A1,B1,N1,RA)
      CALL CINVRT(A2,B2,N1,RA)

!      PRINT'('' IMAGE OF OCTAGON-VERTEX ON UNIT-CIRCLE:'')'
!      WRITE(6,62)CX,CY
!      PRINT'('' TAYLOR COEFFICIENTS, SERIES A1,A2,B1,B2:'')'
!      WRITE(6,63)B
!      DO I=1,N1
!       WRITE(6,64)I,A1(I),A2(I),B1(I),B2(I)
!      ENDDO
62    FORMAT(2(1X,E12.6))
63    FORMAT(4X,'0',13X,4(1X,E12.6))
64    FORMAT(1X,I4,4(1X,E12.6))
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
C                   SUBROUTINE OTOC1                                           C
C   Transform to earth-centered cartesians from octagon-1 (OTOC1) or from      C
C   octagon-2 (OTOC2).                                                         C
C                                                                              C
C --> X,Y   x and y map-coordinate of point in octagon-1 (OTOC1) or            C
C           octagon-2 (OTOC2)                                                  C
C <-- XE    3-vector of earth-centered cartesian coordinates corresponding     C
C           to map location (X,Y).                                             C
C------------------------------------------------------------------------------C
      SUBROUTINE OTOC1(X,Y,XE)
      COMPLEX CI,CIR,CIQ,Z,W
      COMMON/CSTOCT/A1(40),A2(40),B1(40),B2(40),M1G,MAG,A,B,CX,CY
     *,SSR,N1,ROTM(3,3),SC,SCI,ST,STI,CI,CIR,CIQ
      DIMENSION XE(3)
      Z=CMPLX(X,Y)
      CALL TOCT(Z,W)
      W=W*SC
      XA=REAL(W)
      YA=AIMAG(W)
      XX=XA**2+YA**2
      S=1./(1.+XX)
      XE(1)=(1.-XX)*S
      XE(2)=2.*XA*S
      XE(3)=2.*YA*S
      CALL NMAPT(ROTM,XE,XE)
      RETURN

      ENTRY OTOC2(X,Y,XE)
      Z=CMPLX(-X,Y)
      CALL TOCT(Z,W)
      W=W*SCI
      XA=REAL(W)
      YA=AIMAG(W)
      XX=XA**2+YA**2
      S=-1./(1.+XX)
      XE(1)=(1.-XX)*S
      XE(2)=2.*XA*S
      XE(3)=2.*YA*S
      CALL NMAPT(ROTM,XE,XE)
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
C                   SUBROUTINE CTOO                                            C
C   Transform from earth-centered cartesians to octagon-1 (KMAP=1) or to       C
C   octagon-2 (KMAP=-1)                                                        C
C                                                                              C
C --> XE    earth-centered coordinates (3-vector) of a point                   C
C <-- X,Y   map-coordinates of this point                                      C
C <-- KMAP  map indicator (=1 for octagon-1, =-1 for octagon-2)                C
C------------------------------------------------------------------------------C
      SUBROUTINE CTOO(XE,X,Y,KMAP)
      LOGICAL KX,KY,KXY,KS
      COMPLEX CI,CIR,CIQ,Z,W,CIAP
      COMMON/CSTOCT/A1(40),A2(40),B1(40),B2(40),M1G,MAG,A,B,CX,CY
     *,SSR,N1,ROTM(3,3),SC,SCI,ST,STI,CI,CIR,CIQ
      DIMENSION XE(3)
      CIAP=CI*A+1.
      CALL NMAP(ROTM,XE,XE)
      SCS=SC*SC
      ZAC=(1.-SCS)/(1.+SCS)
      XA=XE(2)
      YA=XE(3)
      ZA=XE(1)
      IF(ZA.GT.ZAC)THEN
       KMAP=1
       ZAI=1./((1.+ZA)*SC)
       X=XA*ZAI
       Y=YA*ZAI
      ELSE
       KMAP=-1
       ZAI=SC/(1.-ZA)
       X=XA*ZAI
       Y=-YA*ZAI
      ENDIF
      KX=X.LT.0.
      IF(KX)X=-X
      KY=Y.LT.0.
      IF(KY)Y=-Y
      KXY=Y.GT.X
      IF(KXY)THEN
       T=X
       X=Y
       Y=T
      ENDIF
      DD1=X**2+Y**2
      DD2=(CX-X)**2+(CY-Y)**2
      W=CMPLX(X,Y)
      KS=DD1.LT.SSR*DD2
      IF(KS)THEN
       W=W**4
       CALL TAY(W,B1,N1,Z)
       Z=CIQ*(Z/CI)**.25
      ELSE
       W=CI*(W-1.)/(W+1.)
       W=W-B
       CALL TAY(W,B2,N1,Z)
       Z=CIR*(CI*Z)**STI
       Z=CIAP-CI*Z
      ENDIF
      X=REAL(Z)
      Y=AIMAG(Z)
      IF(KXY)THEN
       T=X
       X=Y
       Y=T
      ENDIF
      IF(KY)Y=-Y
      IF(KX)X=-X
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
C                   SUBROUTINE TOCT                                            C
C   Transform from complez-Z in the standard unit-octagon to complex-W in the  C
C   unit-circle                                                                C
C------------------------------------------------------------------------------C
      SUBROUTINE TOCT(Z,W)
      COMPLEX CI,CIR,CIQ,Z,W,ZT
      LOGICAL KX,KY,KXY,KX1,KXY1,KS
      COMMON/CSTOCT/A1(40),A2(40),B1(40),B2(40),M1G,MAG,A,B,CX,CY
     *,SSR,N1,ROTM(3,3),SC,SCI,ST,STI,CI,CIR,CIQ
      X=REAL(Z)
      Y=AIMAG(Z)
      KX=X.LT.0.
      IF(KX)X=-X
      KY=Y.LT.0.
      IF(KY)Y=-Y
      KXY=Y.GT.X
      IF(KXY)THEN
       T=X
       X=Y
       Y=T
      ENDIF
      KX1=X.GT.1.
      KXY1=X+Y.GT.1.+A
      IF(KX1)THEN
       X=2.-X
      ELSEIF(KXY1)THEN
       T=X
       X=1.+A-Y
       Y=1.+A-T
      ENDIF
      DD1=X**2+Y**2
      DD2=(1.-X)**2+(A-Y)**2

      ZT=CMPLX(X,Y)
      KS=DD1.LT.SSR*DD2
      IF(KS)THEN
       ZT=ZT**4
       CALL TAY(ZT,A1,N1,W)
       W=CIQ*(W/CI)**.25
      ELSE
       ZT=A+CI*ZT-CI
       ZT=ZT**ST
       CALL TAY(ZT,A2,N1,W)
       W=W+B
       W=(W+CI)/(CI-W)
      ENDIF
      IF(KX1.OR.KXY1)W=W/CABS(W)**2
      X=REAL(W)
      Y=AIMAG(W)
      IF(KXY)THEN
       T=X
       X=Y
       Y=T
      ENDIF
      IF(KY)Y=-Y
      IF(KX)X=-X
      W=CMPLX(X,Y)
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
C                   SUBROUTINE ININMAP                                         C
C  Initialize the rotation matrix ROT3 needed to transform standard            C
C  earth-centered cartesian components to the alternative cartesian frame      C
C  oriented so as to put geographical point (ALAT0,ALON0) on the projection    C
C  axis.                                                                       C
C------------------------------------------------------------------------------C
      SUBROUTINE ININMAP(ALON0,ALAT0,ROT3)
      DIMENSION ROT3(3,3)
      DIMENSION X(3),T(3),U(3)
      PI=4*ATAN(1.)
      DR=PI/180.
      BLON0=DR*ALON0
      BLAT0=DR*ALAT0
      CLON0=COS(BLON0)
      SLON0=SIN(BLON0)
      CLAT0=COS(BLAT0)
      SLAT0=SIN(BLAT0)
      ROT3(1,1)=CLAT0*CLON0
      ROT3(1,2)=CLAT0*SLON0
      ROT3(1,3)=SLAT0
      ROT3(2,1)=-SLON0
      ROT3(2,2)=CLON0
      ROT3(2,3)=0.
      ROT3(3,1)=-SLAT0*CLON0
      ROT3(3,2)=-SLAT0*SLON0
      ROT3(3,3)=CLAT0
      RETURN

      ENTRY NMAP(ROT3,X,T)
      CALL MULMM(ROT3,X,U,3,3,1,3,3,3)
      DO I=1,3
      T(I)=U(I)
      ENDDO
      RETURN
      ENTRY NMAPT(ROT3,X,T)
      CALL MULTM(ROT3,X,U,3,3,1,3,3,3)
      DO I=1,3
      T(I)=U(I)
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE STOE                                            C
C   Transform latitude and longitude (degrees) to earth-centered cartesian     C
C   coordinates.                                                               C
C  --> DLAT     latitude                                                       C
C  --> DLON     longitude                                                      C
C  <-- XE       three cartesian components.                                    C
C------------------------------------------------------------------------------C
      SUBROUTINE STOE(DLAT,DLON,XE)
      DIMENSION XE(3)
      DATA DTOR/1.7453292E-2/
      RLAT=DTOR*DLAT
      RLON=DTOR*DLON
      SLA=SIN(RLAT)
      CLA=COS(RLAT)
      SLO=SIN(RLON)
      CLO=COS(RLON)
      XE(1)=CLA*CLO
      XE(2)=CLA*SLO
      XE(3)=SLA
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE  TAY                                            C
C  Evaluate the complex function W of Z whose real                             C
C  Taylor series coefficients are RA.                                          C
C                                                                              C
C  --> Z    function argument (complex)                                        C
C  --> RA   Taylor coefficients (real)                                         C
C  --> N    number of coeffients (starting with the linear term)               C
C  <-- W    Taylor-series approximation of the function (complex)              C
C------------------------------------------------------------------------------C
      SUBROUTINE TAY(Z,RA,N,W)
      COMPLEX Z,W
      DIMENSION RA(*)
      W=0.
      DO I=N,1,-1
       W=(W+RA(I))*Z
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE CINVRT                                          C
C  Compute the Taylor series coefficients Z for the functional-inverse of      C
C  the function whose Taylor series coefficients are W, for the case where     C
C  the constant coefficients of both Z and W are 0                             C
C                                                                              C
C  --> W    Taylor coefficients of original function (starting with linear)    C
C  <-- Z    Taylor coefficients of inverse function                            C
C  --> M    number of Taylor series coefficients computed                      C
C  --- WORK workspace array consisting of at least 3*M elements                C
C------------------------------------------------------------------------------C
      SUBROUTINE CINVRT(W,Z,M,WORK)
      DIMENSION W(*),Z(*),WORK(M,3)
      W1=W(1)
      DO I=1,M
       WORK(I,1)=W(I)
       WORK(I,2)=0.
      ENDDO
      WORK(1,2)=1.
      DO J=1,M
       ZJ=WORK(J,2)/WORK(J,1)
       Z(J)=ZJ
       DO I=J,M
        WORK(I,2)=WORK(I,2)-ZJ*WORK(I,1)
       ENDDO
       CALL CONV(WORK(1,1),W,WORK(1,3),M)
       DO I=1,M
        WORK(I,1)=WORK(I,3)
       ENDDO
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE CONV                                            C
C   Convolve double-precision series A with B to form C, up to M terms         C
C   starting with element 1                                                    C
C                                                                              C
C A,B   --> inputs (convolution factors)                                       C
C C     <-- output (convolution product)                                       C
C M     --> number of elements of A, B, C                                      C
C------------------------------------------------------------------------------C
      SUBROUTINE CONV(A,B,C,M)
      DIMENSION A(*),B(*),C(*)
      DO K=1,M
       C(K)=0.
      ENDDO
      DO I=1,M
       DO J=1,M-I
        K=I+J
        C(K)=C(K)+A(I)*B(J)
       ENDDO
      ENDDO
      RETURN
      END

      SUBROUTINE GETROT(ROT1)
      COMPLEX CI,CIR,CIQ
      COMMON/CSTOCT/A1(40),A2(40),B1(40),B2(40),M1G,MAG,A,B,CX,CY
     *,SSR,N1,ROTM(3,3),SC,SCI,ST,STI,CI,CIR,CIQ
      DIMENSION ROT1(3,3)
      DO J=1,3
      DO I=1,3
       ROT1(I,J)=ROTM(I,J)
      ENDDO
      ENDDO
      RETURN
      END
C                                               *****************
C                                               *    MAT1.FOR   *
C                                               *  PURSER 1994  *
C                                               *****************
C   General matrix routines
C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE MULMM                                           C
C  This, and subsequent routines, perform basic algebraic operations on real   C
C  matrices. The task performed by each routine and entry is essentially       C
C  encoded in each routine's name; the first three letters describe the        C
C  operation, the remainder defining the type of operand. This is the "code":  C
C                                                                              C
C  OPERATIONS:                                                                 C
C   ADD     add first two operands, return result as third argument            C
C   CHL     perform cholesky decomposition                                     C
C   CON     copy negative of first operand to second argument                  C
C   COP     copy first operand to second argument                              C
C   DET     evaluate log-determinant                                           C
C   DIV     divide                                                             C
C   DOT     dot product                                                        C
C   INV     invert the matrix, or linear system involving the matrix operand   C
C   MAD     multiply first two operands, but then add result to third          C
C   MUL     multiply first two operands, return result as third argument       C
C   MSB     multiply first two operands, but then subtract result from third   C
C   SUB     subtract first two operands, return result as third argument       C
C   WRT     write out                                                          C
C   ZER     set to zero                                                        C
C           A prefix of "C" on some of the later routines indicates complex    C
C           arguments.                                                         C
C                                                                              C
C  OPERANDS:                                                                   C
C   B       banded matrix                                                      C
C   D       diagonal matrix                                                    C
C   L       lower triangular matrix                                            C
C   M       matrix                                                             C
C   S       scalar                                                             C
C   T       transpose of the matrix                                            C
C   U       upper triangular matrix                                            C
C   V       vector                                                             C
C                                                                              C
C   For matrix routines, the arguments, MI, MJ, etc usually denote the         C
C   sequence of algebraic dimensions of the matrices/vectors used in           C
C   left-to-right order. Thus, for general matrix multiplication, the          C
C   first, MI, is the number of active rows of the first factor, MJ the        C
C   number of active columns of the first factor and the number of active      C
C   rows of the second factor, MK the number of active columns of the          C
C   second factor (and of the product). The arguments, NA, NB, etc usually     C
C   denote the fortran first dimensions of each matrix (but not vector)        C
C   arguments of the routine in the order in which they appear. This allows    C
C   for the fortran dimension of the matrices to exceed the algebraic          C
C   or "logical" dimensions. If it is known in advance that all matrices       C
C   are square with algebraic dimensions identical to the fortran dimensions   C
C   a convenient alternative to these general routine is the set of            C
C   "quick format" routines prefixed by "Q" in which the single dimension N    C
C   replaces the MI,MJ,..., NA,NB,... used in the general routines. Quick      C
C   format routines are grouped in module MAT3.FOR.                            C
C                                                                              C
C   A special collection of various "banded matrix" routine is found in the    C
C   later modules of this series (MAT6 MAT9 MAT10)                             C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE MULMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      LOGICAL OMUL
      DIMENSION A(NA,*),B(NB,*),C(NC,*)
      OMUL=.TRUE.
      GOTO 320
      ENTRY MADMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
320   DO I=1,MI
      DO K=1,MK
      IF(OMUL)C(I,K)=0.
      DO J=1,MJ
        C(I,K)=C(I,K)+A(I,J)*B(J,K)
      end do
      end do
      end do
      RETURN
      ENTRY MULMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.TRUE.
      GOTO 321
      ENTRY MADMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
321   DO I=1,MI
      DO K=1,MK
      IF(OMUL)C(I,K)=0.
      DO J=1,MJ
        C(I,K)=C(I,K)+A(I,J)*B(K,J)
      end do
      end do
      end do
      RETURN
      ENTRY MULTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.TRUE.
      GOTO 322
      ENTRY MADTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
322   DO I=1,MI
      DO K=1,MK
      IF(OMUL)C(I,K)=0.
      DO J=1,MJ
        C(I,K)=C(I,K)+A(J,I)*B(J,K)
      end do
      end do
      end do
      RETURN
      ENTRY MULTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.TRUE.
      GOTO 323
      ENTRY MADTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
323   DO I=1,MI
      DO K=1,MK
      IF(OMUL)C(I,K)=0.
      DO J=1,MJ
         C(I,K)=C(I,K)+A(J,I)*B(K,J)
      end do
      end do
      end do
      RETURN
      ENTRY MSBMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO I=1,MI
      DO K=1,MK
      DO J=1,MJ
        C(I,K)=C(I,K)-A(I,J)*B(J,K)
      end do
      end do
      end do
      RETURN
      ENTRY MSBMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO I=1,MI
      DO K=1,MK
      DO J=1,MJ
        C(I,K)=C(I,K)-A(I,J)*B(K,J)
      end do
      end do
      end do
      RETURN
      ENTRY MSBTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO I=1,MI
      DO K=1,MK
      DO J=1,MJ
        C(I,K)=C(I,K)-A(J,I)*B(J,K)
      end do
      end do
      end do
      RETURN
      ENTRY MSBTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO I=1,MI
      DO K=1,MK
      DO  J=1,MJ
        C(I,K)=C(I,K)-A(J,I)*B(K,J)
      end do
      end do
      end do
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE ADDMM etc                                       C
C   Add and subtract matrices                                                  C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE ADDMM(A,B,C,MI,MJ,NA,NB,NC)
      DIMENSION A(NA,*),B(NB,*),C(NC,*)
      DO I=1,MI
      DO J=1,MJ
        C(I,J)=A(I,J)+B(I,J)
      end do
      end do
      RETURN
      ENTRY ADDMT(A,B,C,MI,MJ,NA,NB,NC)
      DO I=1,MI
      DO  J=1,MJ
        C(I,J)=A(I,J)+B(J,I)
      end do
      end do
      RETURN
      ENTRY ADDTM(A,B,C,MI,MJ,NA,NB,NC)
      DO I=1,MI
      DO J=1,MJ
        C(I,J)=A(J,I)+B(I,J)
      end do
      end do
      RETURN
      ENTRY ADDTT(A,B,C,MI,MJ,NA,NB,NC)
      DO I=1,MI
      DO J=1,MJ
        C(I,J)=A(J,I)+B(J,I)
      end do
      end do
      RETURN
      ENTRY SUBMM(A,B,C,MI,MJ,NA,NB,NC)
      DO I=1,MI
      DO J=1,MJ
        C(I,J)=A(I,J)-B(I,J)
      end do
      end do
      RETURN
      ENTRY SUBMT(A,B,C,MI,MJ,NA,NB,NC)
      DO  I=1,MI
      DO  J=1,MJ
        C(I,J)=A(I,J)-B(J,I)
      end do
      end do
      RETURN
      ENTRY SUBTM(A,B,C,MI,MJ,NA,NB,NC)
      DO I=1,MI
      DO J=1,MJ
        C(I,J)=A(J,I)-B(I,J)
      end do
      end do
      RETURN
      ENTRY SUBTT(A,B,C,MI,MJ,NA,NB,NC)
      DO I=1,MI
      DO J=1,MJ
        C(I,J)=A(J,I)-B(J,I)
      end do
      end do
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE COPM etc                                        C
C  Mainly routines to copy matrices                                            C
C------------------------------------------------------------------------------C
      SUBROUTINE COPM(A,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*)
      ENTRY EQMM(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=A(I,J)
      end do
      end do
      RETURN
      ENTRY NEQMM(A,B,MI,MJ,NA,NB)
      ENTRY CONM(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
         B(I,J)=-A(I,J)
      ENDDO
      ENDDO
      RETURN
      ENTRY MULMS(A,SS,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
         B(I,J)=A(I,J)*SS
      ENDDO
      ENDDO
      RETURN
      ENTRY EQTM(A,B,MI,MJ,NA,NB)
      ENTRY COPT(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
       B(I,J)=A(J,I)
      ENDDO
      ENDDO
      RETURN
      ENTRY CONT(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
       B(I,J)=-A(J,I)
      ENDDO
      ENDDO
      RETURN
      ENTRY ZERM(A,MI,MJ,NA)
      DO I=1,MI
      DO J=1,MJ
       A(I,J)=0.
      ENDDO
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE INVL                                            C
C     Invert lower triangular matrix, possibly in place if A and B are same    C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE INVL(A,B,MI,NA,NB)
      DIMENSION A(NA,*),B(NB,*)
         DO J=MI,1,-1
      JM=J-1
      JP=J+1
      DO I=1,JM
        B(I,J)=0.
      end do
      B(J,J)=1./A(J,J)
      DO I=JP,MI
      IM=I-1
      S=0.
      DO K=J,IM
        S=S+A(I,K)*B(K,J)
      end do
        B(I,J)=-B(I,I)*S
      end do
         end do
      RETURN
      END


C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE INVLV                                           C
C     Solve linear system involving lower triangular (INVLV) or upper          C
C     triangular (INVUV) matrix, right-hand-side vector U and output vector V  C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE INVLV(A,U,V,MI,NA)
      DIMENSION A(NA,*),U(*),V(*)
      DO I=1,MI
       S=U(I)
       DO J=1,I-1
        S=S-A(I,J)*V(J)
       ENDDO
       V(I)=S/A(I,I)
      ENDDO
      RETURN
      ENTRY INVUV(A,U,V,MI,NA)
      DO J=MI,1,-1
       S=U(J)
       DO I=J+1,MI
        S=S-A(I,J)*V(I)
       ENDDO
       V(J)=S/A(J,J)
      ENDDO
      RETURN
      END

      FUNCTION DOT(A,B,M)
      DIMENSION A(M),B(M)
      DOT=0.
      DO I=1,M
       DOT=DOT+A(I)*B(I)
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE MULMD etc                                       C
C  Multiply matrices with diagonal matrices                                    C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE MULMD(A,D,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*),D(*)
      DO I=1,MI
      DO  J=1,MJ
        B(I,J)=A(I,J)*D(J)
      end do
      end do
      RETURN
      ENTRY MULTD(A,D,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO  J=1,MJ
        B(I,J)=A(J,I)*D(J)
      end do
      end do
      RETURN
      ENTRY MULDM(D,A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO  J=1,MJ
        B(I,J)=D(I)*A(I,J)
      end do
      end do
      RETURN
      ENTRY MULDT(D,A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=D(I)*A(J,I)
      end do
      end do
      RETURN
      END

      SUBROUTINE MADMD(A,D,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*),D(*)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=B(I,J)+A(I,J)*D(J)
      end do
      end do
      RETURN
      ENTRY MADTD(A,D,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO  J=1,MJ
        B(I,J)=B(I,J)+A(J,I)*D(J)
      end do
      end do
      RETURN
      ENTRY MADDM(D,A,B,MI,MJ,NA,NB)
      DO  I=1,MI
      DO  J=1,MJ
        B(I,J)=B(I,J)+D(I)*A(I,J)
      end do
      end do
      RETURN
      ENTRY MADDT(D,A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=B(I,J)+D(I)*A(J,I)
      end do
      end do
      RETURN
      ENTRY ADDMD(A,D,B,MI,NA,NB)
      DO I=1,MI
       DO J=1,MI
        B(I,J)=A(I,J)
       ENDDO
       B(I,I)=B(I,I)+D(I)
      ENDDO
      RETURN
      END

      SUBROUTINE MULDD(A,B,C,M)
      DIMENSION A(*),B(*),C(*)
      DO I=1,M
        C(I)=A(I)*B(I)
      end do
      RETURN
      ENTRY ADDDD(A,B,C,M)
      DO I=1,M
        C(I)=A(I)+B(I)
      end do
      RETURN
      ENTRY SUBDD(A,B,C,M)
      DO I=1,M
        C(I)=A(I)-B(I)
      end do
      RETURN
      ENTRY DIVDD(A,B,C,M)
      DO I=1,M
        C(I)=A(I)/B(I)
      end do
      RETURN
      ENTRY EQDD(A,B,M)
      ENTRY COPD(A,B,M)
      DO I=1,M
       B(I)=A(I)
      ENDDO
      RETURN
      END

      SUBROUTINE COPSM(S,A,M,NA)
      DIMENSION A(NA,*),B(NB,*),D(*)
      ENTRY EQSM(S,A,M,NA)
      DO I=1,M
      DO J=1,M
        A(I,J)=0.
      end do
        A(I,I)=S
      end do
      RETURN
      ENTRY EQDM(D,A,M,NA)
      ENTRY COPDM(D,A,M,NA)
      DO I=1,M
      DO J=1,M
        A(I,J)=0.
      end do
        A(I,I)=D(I)
      end do
      RETURN
      ENTRY SUBSM(S,A,B,MI,NA,NB)
C  SUBTRACT MATRIX A FROM SCALAR MATRIX S TO GET B
      DO J=1,MI
      DO I=1,MI
      B(I,J)=-A(I,J)
      ENDDO
      B(J,J)=B(J,J)+S
      ENDDO
      RETURN
      ENTRY ADDSM(S,A,B,M,NA,NB)
C  ADD SCALAR S TO DIAGONALS OF SQUARE MATRIX A, RESULT IN B
      DO I=1,M
      DO J=1,M
        B(I,J)=A(I,J)
      end do
        B(I,I)=B(I,I)+S
      end do
      RETURN
      ENTRY ADDDM(D,A,B,M,NA,NB)
C  ADD DIAGONAL D TO DIAGONALS OF SQUARE MATRIX A, RESULT IN B
      DO I=1,M
      DO J=1,M
         B(I,J)=A(I,J)
      end do
        B(I,I)=B(I,I)+D(I)
      end do
      RETURN
      ENTRY EQLD(A,D,M,NA)
      ENTRY COPLD(A,D,M,NA)
      DO I=1,M
      D(I)=A(I,I)
      ENDDO
      RETURN
      ENTRY DETL(A,DET,MI,NA)
C  GET LOG OF DETERMINANT OF LOWER TRIANGULAR MATRIX
      DET=0.
      DO I=1,MI
         DET=DET+ALOG(A(I,I))
      end do
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE WRTM, WRTIM                                     C
C  Write out contents of a real (WRTM) or integer (WRTIM) matrix               C
C                                                                              C
C  --> A or L   matrix                                                         C
C  --> MI       number of active rows of matrix                                C
C  --> MJ       number of active columns of matrix                             C
C  --> JCOLS    number of columns output at a time                             C
C  --> N        first fortran dimension of matrix                              C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE WRTM(A,MI,MJ,JCOLS,N)
C  WRITE A REAL ARRAY
      DIMENSION A(N,*),L(N,*)
        DO J1=1,MJ,JCOLS
      J2=MIN(MJ,J1+JCOLS-1)
      WRITE(6,600)(J,J=J1,J2)
      WRITE(6,604)('________',J=J1,J2)
      DO I=1,MI
         WRITE(6,64)I,(A(I,J),J=J1,J2)
      end do
      WRITE(6,60)
        end do
      RETURN
      ENTRY WRTIM(L,MI,MJ,JCOLS,N)
C  WRITE AN INTEGER ARRAY
        DO J1=1,MJ,JCOLS
      J2=MIN(MJ,J1+JCOLS-1)
      WRITE(6,601)(J,J=J1,J2)
      WRITE(6,604)('________',J=J1,J2)
      DO I=1,MI
        WRITE(6,61)I,(L(I,J),J=J1,J2)
      end do
      WRITE(6,60)
      end do
60    FORMAT(/)
61    FORMAT(1X,I3,'3',8(1X,I5))
64    FORMAT(1X,I3,10(1X,E11.5))
600   FORMAT(5X,10(I3,9X))
601   FORMAT(9X,20(I3,3X))
604   FORMAT(1X,'________',11A8)
      RETURN
      END
C                                               *****************
C                                               *    MAT2.FOR   *
C                                               *  PURSER 1993  *
C                                               *****************
C    Routines for linear inversions etc.

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE LUFM                                            C
C  perform l-u decomposition of square matrix a in place with                  C
C  partial pivoting                                                            C
C  For DOUBLE PRECISION version see DLUFM                                      C
C                                                                              C
C  --> a    square matrix to be factorized                                     C
C  <-- ipiv array encoding the pivoting sequence                               C
C  <-- d    indicator for possible sign change of determinant                  C
C  --> m    degree of (active part of) a                                       C
C  --> na   first fortran dimension of a                                       C
C                                                                              C
C------------------------------------------------------------------------------C
      subroutine LUFM(a,ipiv,d,m,na)
      dimension a(na,*),ipiv(*)
      d=1.
      ipiv(m)=m
      do j=1,m-1
       jp=j+1
       abig=abs(a(j,j))
       ibig=j
       do i=jp,m
        aa=abs(a(i,j))
        if(aa.gt.abig)then
         ibig=i
         abig=aa
        endif
       enddo
c  swap rows, recording changed sign of determinant
       ipiv(j)=ibig
       if(ibig.ne.j)then
        d=-d
        do k=1,m
         t=a(j,k)
         a(j,k)=a(ibig,k)
         a(ibig,k)=t
        enddo
       endif
       ajj=a(j,j)
       if(ajj.eq.0.)then
                           jm=j-1
                           print'('' failure in lufact:''/
     +                     '' matrix singular, rank='',i3)',jm
                           stop
       endif
       ajji=1./ajj
       do i=jp,m
        aij=ajji*a(i,j)
        a(i,j)=aij
        do k=jp,m
         a(i,k)=a(i,k)-aij*a(j,k)
        enddo
       enddo
      enddo
      return
      end

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE INVMM                                           C
C  invert matrix, possibly in place (a=b), using the l-u decomposition method  C
C  For DOUBLE PRECISION version see DINVMM                                     C
C                                                                              C
C  --> b    square matrix to be inverted                                       C
C  <-- a    inverse of b                                                       C
C  --> m    degree of (active part of) b and a                                 C
C  --> nb   first fortran dimension of b                                       C
C  --> na   first fortran dimension of a                                       C
C                                                                              C
C   LIMITATION:                                                                C
C    ipiv is an index array, internal to this array, encoding the              C
C    pivoting sequence used. It is given a fortran dimension of NN=500         C
C    in the parameter statement below. If the order of the linear system       C
C    exceeds 500, increase this parameter accordingly                          C
C                                                                              C
C------------------------------------------------------------------------------C
      subroutine INVMM(b,a,m,nb,na)
      PARAMETER (NN=500)
      DIMENSION IPIV(NN)
      dimension a(na,*),b(nb,*)
      do j=1,m
      do i=1,m
       a(i,j)=b(i,j)
      enddo
      enddo
      call lufm(a,ipiv,d,m,na)
c  invert u in place:
      do i=1,m
       a(i,i)=1./a(i,i)
      enddo
      do i=1,m-1
       do j=i+1,m
        s=0.
        do k=i,j-1
         s=s-a(i,k)*a(k,j)
        enddo
        a(i,j)=a(j,j)*s
       enddo
      enddo
c  invert l in place assuming implicitly diagonal elements of unity
      do j=1,m-1
       do i=j+1,m
        s=-a(i,j)
        do k=j+1,i-1
         s=s-a(i,k)*a(k,j)
        enddo
        a(i,j)=s
       enddo
      enddo
c  form the product of u**-1 and l**-1 in place
      do j=1,m-1
       do i=1,j
        s=a(i,j)
        do k=j+1,m
         s=s+a(i,k)*a(k,j)
        enddo
        a(i,j)=s
       enddo
       do i=j+1,m
        s=0.
        do k=i,m
         s=s+a(i,k)*a(k,j)
        enddo
        a(i,j)=s
       enddo
      enddo
c  permute columns according to ipiv
      do j=m-1,1,-1
       l=ipiv(j)
       do i=1,m
        t=a(i,j)
        a(i,j)=a(i,l)
        a(i,l)=t
       enddo
      enddo
      return

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE LINVMM                                          C
C  invert linear systems sharing the same square system matrix                 C
C  For DOUBLE PRECISION version see DLINVMM                                    C
C                                                                              C
C  <-> a    square system matrix on entry, its L-D-U factorization on return   C
C  <-> b    right-hand-sides on entry, corresponding matrix of solution        C
C           vectors on return                                                  C
C  --> m    degree of (active part of) b and a                                 C
C  --> mm   number of right-hand-side vectors (active columns of b)            C
C  --> na   first fortran dimension of a                                       C
C  --> nb   first fortran dimension of b                                       C
C                                                                              C
C------------------------------------------------------------------------------C
      entry                            LINVMM(a,b,m,mm,na,nb)
      call lufm(a,ipiv,d,m,na)
      call lubmm(a,b,ipiv,m,mm,na,nb)
      return
      end


C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE LUBMM                                           C
C  use l-u factors in a to back-substitute for mm rhs in b, using ipiv to      C
C  define the pivoting permutation used in the l-u decomposition.              C
C  For DOUBLE PRECISION version see DLUBMM                                     C
C                                                                              C
C  --> a    L-D-U factorization of linear system matrux                        C
C  <-> b    right-hand-sides on entry, corresponding matrix of solution        C
C           vectors on return                                                  C
C  --> ipiv array encoding the pivoting sequence                               C
C  --> m    degree of (active part of) b and a                                 C
C  --> mm   number of right-hand-side vectors (active columns of b)            C
C  --> na   first fortran dimension of a                                       C
C  --> nb   first fortran dimension of b                                       C
C                                                                              C
C------------------------------------------------------------------------------C
      subroutine LUBMM(a,b,ipiv,m,mm,na,nb)
      dimension a(na,*),b(nb,*),ipiv(*)
      do k=1,mm !loop over columns of b
       do i=1,m
        l=ipiv(i)
        s=b(l,k)
        b(l,k)=b(i,k)
        do j=1,i-1
         s=s-a(i,j)*b(j,k)
        enddo
        b(i,k)=s
       enddo
       do i=m,1,-1
        s=b(i,k)
        do j=i+1,m
         s=s-a(i,j)*b(j,k)
        enddo
        b(i,k)=s/a(i,i)
       enddo
      enddo
      return
      end

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE INVHH                                           C
C  Invert, possibly in place, a symmetric matrix                               C
C  For DOUBLE PRECISION version see DINVHH                                     C
C                                                                              C
C  --> B    symmetric square input matrix                                      C
C  <-- A    inverse of B                                                       C
C  --> M    degree of (active part of) A and B                                 C
C  --> NB   first fortran dimension of B                                       C
C  --> NA   first fortran dimension of A                                       C
C                                                                              C
C  LIMITATION                                                                  C
C     This routine incorporates no pivoting - it is intended for matrices      C
C     that are already diagonally dominant                                     C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE INVHH(B,A,M,NB,NA)
      DIMENSION B(NB,*),A(NA,*)
      DO I=1,M
      DO J=1,M
      A(I,J)=B(I,J)
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
C  PERFORM L.D.U DECOMPOSITION OF SYMMETRIC MATRIX ( U = L* )
C-----------------------------------------------------------------------
      DO J=1,M
       JM=J-1
       AJJ=A(J,J)
       DO K=1,JM
        AJJ=AJJ-A(K,K)*A(J,K)**2
       ENDDO
       A(J,J)=AJJ
       DO I=J+1,M
        AIJ=A(I,J)
        DO K=1,JM
         AIJ=AIJ-A(I,K)*A(K,K)*A(J,K)
        ENDDO
        A(I,J)=AIJ/AJJ
        A(J,I)=0.
       ENDDO
      ENDDO
C-----------------------------------------------------------------------
C  INVERT (IN PLACE) THE LOWER TRIANGULAR PART OF A, (ASSUMING UNIT
C  DIAGONAL ELEMENTS), AND INVERT THE DIAGONAL PART OF A (ASSUMING
C  ZERO OFF-DIAGONAL ELEMNTS)
C-----------------------------------------------------------------------
      DO K=1,M
       KP=K+1
       A(K,K)=1./A(K,K)
       IM=K
       DO I=KP,M
        AIK=A(I,K)
        DO J=KP,IM
         AIK=AIK+A(I,J)*A(J,K)
        ENDDO
        A(I,K)=-AIK
        IM=I
       ENDDO
      ENDDO
C-----------------------------------------------------------------------
C  MULTIPLY: THE TRANSPOSE OF THE LOWER PART OF A (ASSUMING UNIT DIAGS),
C  TIMES THE DIAGONAL PART (ASSUMING ZERO OFF-DIAGS), TIMES THE LOWER
C  PART. THIS PRODUCT IS THE SYMMETRIC INVERSE OF THE ORIGINAL B.
C-----------------------------------------------------------------------
      DO I=1,M
       AII=A(I,I)
       DO J=1,I-1
        AIJ=AII*A(I,J)
        DO K=I+1,M
         AIJ=AIJ+A(K,I)*A(K,K)*A(K,J)
        ENDDO
        A(I,J)=AIJ
        A(J,I)=AIJ
       ENDDO
       DO K=I+1,M
        AII=AII+A(K,I)**2*A(K,K)
       ENDDO
       A(I,I)=AII
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE DLUFM                                           C
C  perform l-u decomposition of square matrix a in place with                  C
C  partial pivoting IN DOUBLE PRECISION                                        C
C  For SINGLE PRECISION version see LUFM                                       C
C                                                                              C
C  --> a    square matrix to be factorized                                     C
C  <-- ipiv array encoding the pivoting sequence                               C
C  <-- d    indicator for possible sign change of determinant                  C
C  --> m    degree of (active part of) a                                       C
C  --> na   first fortran dimension of a                                       C
C                                                                              C
C------------------------------------------------------------------------------C
      subroutine DLUFM(a,ipiv,d,m,na)
c      implicit double precision (a-h,o-z)
      dimension a(na,*),ipiv(*)
      d=1.
      ipiv(m)=m
      do j=1,m-1
       jp=j+1
       abig=abs(a(j,j))
       ibig=j
       do i=jp,m
        aa=abs(a(i,j))
        if(aa.gt.abig)then
         ibig=i
         abig=aa
        endif
       enddo
c  swap rows, recording changed sign of determinant
       ipiv(j)=ibig
       if(ibig.ne.j)then
        d=-d
        do k=1,m
         t=a(j,k)
         a(j,k)=a(ibig,k)
         a(ibig,k)=t
        enddo
       endif
       ajj=a(j,j)
c       if(ajj.eq.0.d0)then
       if(ajj.eq.0)then
                           jm=j-1
                           print'('' failure in lufact:''/
     +                     '' matrix singular, rank='',i3)',jm
                           stop
       endif
       ajji=1./ajj
       do i=jp,m
        aij=ajji*a(i,j)
        a(i,j)=aij
        do k=jp,m
         a(i,k)=a(i,k)-aij*a(j,k)
        enddo
       enddo
      enddo
      return
      end

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE DINVMM                                          C
C  invert matrix, possibly in place (a=b), using the l-u decomposition method  C
C  in DOUBLE PRECISION (for single precision version see INVMM                 C
C                                                                              C
C  --> b    square matrix to be inverted                                       C
C  <-- a    inverse of b                                                       C
C  --> m    degree of (active part of) b and a                                 C
C  --> nb   first fortran dimension of b                                       C
C  --> na   first fortran dimension of a                                       C
C                                                                              C
C   LIMITATION:                                                                C
C    ipiv is an index array, internal to this array, encoding the              C
C    pivoting sequence used. It is given a fortran dimension of NN=500         C
C    in the parameter statement below. If the order of the linear system       C
C    exceeds 500, increase this parameter accordingly                          C
C                                                                              C
C------------------------------------------------------------------------------C
      subroutine DINVMM(b,a,m,nb,na)
      PARAMETER (NN=500)
c      implicit double precision (a-h,o-z)
      DIMENSION IPIV(NN)
      dimension a(na,*),b(nb,*)
      do j=1,m
      do i=1,m
       a(i,j)=b(i,j)
      enddo
      enddo
      call DLUFM(a,ipiv,d,m,na)
c  invert u in place:
      do i=1,m
       a(i,i)=1./a(i,i)
      enddo
      do i=1,m-1
       do j=i+1,m
        s=0.
        do k=i,j-1
         s=s-a(i,k)*a(k,j)
        enddo
        a(i,j)=a(j,j)*s
       enddo
      enddo
c  invert l in place assuming implicitly diagonal elements of unity
      do j=1,m-1
       do i=j+1,m
        s=-a(i,j)
        do k=j+1,i-1
         s=s-a(i,k)*a(k,j)
        enddo
        a(i,j)=s
       enddo
      enddo
c  form the product of u**-1 and l**-1 in place
      do j=1,m-1
       do i=1,j
        s=a(i,j)
        do k=j+1,m
         s=s+a(i,k)*a(k,j)
        enddo
        a(i,j)=s
       enddo
       do i=j+1,m
        s=0.
        do k=i,m
         s=s+a(i,k)*a(k,j)
        enddo
        a(i,j)=s
       enddo
      enddo
c  permute columns according to ipiv
      do j=m-1,1,-1
       l=ipiv(j)
       do i=1,m
        t=a(i,j)
        a(i,j)=a(i,l)
        a(i,l)=t
       enddo
      enddo
      return

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE DLINVMM                                         C
C  invert linear systems sharing the same square system matrix in              C
C  DOUBLE PRECISION. (For single precision version see LINVMM                  C
C                                                                              C
C  <-> a    square system matrix on entry, its L-D-U factorization on return   C
C  <-> b    right-hand-sides on entry, corresponding matrix of solution        C
C           vectors on return                                                  C
C  --> m    degree of (active part of) b and a                                 C
C  --> mm   number of right-hand-side vectors (active columns of b)            C
C  --> na   first fortran dimension of a                                       C
C  --> nb   first fortran dimension of b                                       C
C                                                                              C
C------------------------------------------------------------------------------C
      entry                            DLINVMM(a,b,m,mm,na,nb)
      call DLUFM(a,ipiv,d,m,na)
      call DLUBMM(a,b,ipiv,m,mm,na,nb)
      return
      end


C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE DUBMM                                           C
C  use l-u factors in a to back-substitute for mm rhs in b, using ipiv to      C
C  define the pivoting permutation used in the l-u decomposition               C
C  in DOUBLE PRECISION. (For single precision version see LUBMM                C
C                                                                              C
C  --> a    L-D-U factorization of linear system matrux                        C
C  <-> b    right-hand-sides on entry, corresponding matrix of solution        C
C           vectors on return                                                  C
C  --> ipiv array encoding the pivoting sequence                               C
C  --> m    degree of (active part of) b and a                                 C
C  --> mm   number of right-hand-side vectors (active columns of b)            C
C  --> na   first fortran dimension of a                                       C
C  --> nb   first fortran dimension of b                                       C
C                                                                              C
C------------------------------------------------------------------------------C
      subroutine DLUBMM(a,b,ipiv,m,mm,na,nb)
c      implicit double precision (a-h,o-z)
      dimension a(na,*),b(nb,*),ipiv(*)
      do k=1,mm !loop over columns of b
       do i=1,m
        l=ipiv(i)
        s=b(l,k)
        b(l,k)=b(i,k)
        do j=1,i-1
         s=s-a(i,j)*b(j,k)
        enddo
        b(i,k)=s
       enddo
       do i=m,1,-1
        s=b(i,k)
        do j=i+1,m
         s=s-a(i,j)*b(j,k)
        enddo
        b(i,k)=s/a(i,i)
       enddo
      enddo
      return
      end

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1993          C
C                   SUBROUTINE DINVHH                                          C
C  Invert, possibly in place, a symmetric matrix in DOUBLE PRECISION.          C
C  (For single precision version see INVHH                                     C
C                                                                              C
C  --> B    symmetric square input matrix                                      C
C  <-- A    inverse of B                                                       C
C  --> M    degree of (active part of) A and B                                 C
C  --> NB   first fortran dimension of B                                       C
C  --> NA   first fortran dimension of A                                       C
C                                                                              C
C  LIMITATION                                                                  C
C     This routine incorporates no pivoting - it is intended for matrices      C
C     that are already diagonally dominant                                     C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE DINVHH(B,A,M,NB,NA)
c      implicit double precision (a-h,o-z)
      DIMENSION B(NB,*),A(NA,*)
      DO I=1,M
      DO J=1,M
      A(I,J)=B(I,J)
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
C  PERFORM L.D.U DECOMPOSITION OF SYMMETRIC MATRIX ( U = L* )
C-----------------------------------------------------------------------
      DO J=1,M
       JM=J-1
       AJJ=A(J,J)
       DO K=1,JM
        AJJ=AJJ-A(K,K)*A(J,K)**2
       ENDDO
       A(J,J)=AJJ
       DO I=J+1,M
        AIJ=A(I,J)
        DO K=1,JM
         AIJ=AIJ-A(I,K)*A(K,K)*A(J,K)
        ENDDO
        A(I,J)=AIJ/AJJ
        A(J,I)=0.
       ENDDO
      ENDDO
C-----------------------------------------------------------------------
C  INVERT (IN PLACE) THE LOWER TRIANGULAR PART OF A, (ASSUMING UNIT
C  DIAGONAL ELEMENTS), AND INVERT THE DIAGONAL PART OF A (ASSUMING
C  ZERO OFF-DIAGONAL ELEMNTS)
C-----------------------------------------------------------------------
      DO K=1,M
       KP=K+1
       A(K,K)=1./A(K,K)
       IM=K
       DO I=KP,M
        AIK=A(I,K)
        DO J=KP,IM
         AIK=AIK+A(I,J)*A(J,K)
        ENDDO
        A(I,K)=-AIK
        IM=I
       ENDDO
      ENDDO
C-----------------------------------------------------------------------
C  MULTIPLY: THE TRANSPOSE OF THE LOWER PART OF A (ASSUMING UNIT DIAGS),
C  TIMES THE DIAGONAL PART (ASSUMING ZERO OFF-DIAGS), TIMES THE LOWER
C  PART. THIS PRODUCT IS THE SYMMETRIC INVERSE OF THE ORIGINAL B.
C-----------------------------------------------------------------------
      DO I=1,M
       AII=A(I,I)
       DO J=1,I-1
        AIJ=AII*A(I,J)
        DO K=I+1,M
         AIJ=AIJ+A(K,I)*A(K,K)*A(K,J)
        ENDDO
        A(I,J)=AIJ
        A(J,I)=AIJ
       ENDDO
       DO K=I+1,M
        AII=AII+A(K,I)**2*A(K,K)
       ENDDO
       A(I,I)=AII
      ENDDO
      RETURN
      END
C                                               *****************
C                                               *    MAT3.FOR   *
C                                               *  PURSER 1993  *
C                                               *****************
C   "Quick-format" matrix routines with maximal simplification of
C    argument lists (mainly for use with square matrices).

      SUBROUTINE QMULMM(A,B,C,M)
c      LOGICAL ODIAG(*)
      DIMENSION A(M,*),B(M,*),C(M,*),U(*),V(*),D(*),L(*)
      CALL MULMM(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMADMM(A,B,C,M)
      CALL MADMM(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMULMT(A,B,C,M)
      CALL MULMT(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMADMT(A,B,C,M)
      CALL MADMT(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMULTM(A,B,C,M)
      CALL MULTM(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMADTM(A,B,C,M)
      CALL MADTM(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMULTT(A,B,C,M)
      CALL MULTT(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMADTT(A,B,C,M)
      CALL MADTT(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMSBMM(A,B,C,M)
      CALL MSBMM(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMSBMT(A,B,C,M)
      CALL MSBMT(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMSBTM(A,B,C,M)
      CALL MSBTM(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QMSBTT(A,B,C,M)
      CALL MSBTT(A,B,C,M,M,M,M,M,M)
      RETURN
      ENTRY QADDMM(A,B,C,M)
      CALL ADDMM(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QADDMT(A,B,C,M)
      CALL ADDMT(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QADDTM(A,B,C,M)
      CALL ADDTM(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QADDTT(A,B,C,M)
      CALL ADDTT(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QSUBMM(A,B,C,M)
      CALL SUBMM(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QSUBMT(A,B,C,M)
      CALL SUBMT(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QSUBTM(A,B,C,M)
      CALL SUBTM(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QSUBTT(A,B,C,M)
      CALL SUBTT(A,B,C,M,M,M,M,M)
      RETURN
      ENTRY QZERM(A,M)
      CALL ZERM(A,M,M,M)
      RETURN
c      ENTRY QCHLML(A,B,M,ODIAG)
c      CALL CHLML(A,B,M,M,M,ODIAG)
      RETURN
      ENTRY QEQMM(A,B,M)
      CALL EQMM(A,B,M,M,M,M)
      RETURN
      ENTRY QNEQMM(A,B,M)
      CALL NEQMM(A,B,M,M,M,M)
      RETURN
      ENTRY QMULMS(A,SS,B,M)
      CALL MULMS(A,SS,B,M,M,M,M)
      RETURN
      ENTRY QEQTM(A,B,M)
      CALL EQTM(A,B,M,M,M,M)
      RETURN
      ENTRY QINVL(A,B,M)
      CALL INVL(A,B,M,M,M)
      RETURN
      ENTRY QINVLV(A,U,V,M)
      CALL INVLV(A,U,V,M,M)
      RETURN
      ENTRY QINVUV(A,U,V,M)
      CALL INVUV(A,U,V,M,M)
      RETURN
      ENTRY QMULMD(A,D,B,M)
      CALL MULMD(A,D,B,M,M,M,M)
      RETURN
      ENTRY QMULTD(A,D,B,M)
      CALL MULTD(A,D,B,M,M,M,M)
      RETURN
      ENTRY QMULDM(D,A,B,M)
      CALL MULDM(D,A,B,M,M,M,M)
      RETURN
      ENTRY QMULDT(D,A,B,M)
      CALL MULDT(D,A,B,M,M,M,M)
      RETURN
      ENTRY QMADMD(A,D,B,M)
      CALL MADMD(A,D,B,M,M,M,M)
      RETURN
      ENTRY QMADTD(A,D,B,M)
      CALL MADTD(A,D,B,M,M,M,M)
      RETURN
      ENTRY QMADDM(D,A,B,M)
      CALL MADDM(D,A,B,M,M,M,M)
      RETURN
      ENTRY QMADDT(D,A,B,M)
      CALL MADDT(D,A,B,M,M,M,M)
      RETURN
      ENTRY QADDMD(A,D,B,M)
      CALL ADDMD(A,D,B,M,M,M)
      RETURN
      ENTRY QMULDD(A,B,C,M)
      CALL MULDD(A,B,C,M)
      RETURN
      ENTRY QADDDD(A,B,C,M)
      CALL ADDDD(A,B,C,M)
      RETURN
      ENTRY QSUBDD(A,B,C,M)
      CALL SUBDD(A,B,C,M)
      RETURN
      ENTRY QDIVDD(A,B,C,M)
      CALL DIVDD(A,B,C,M)
      RETURN
      ENTRY QEQDD(A,B,M)
      CALL EQDD(A,B,M)
      RETURN
      ENTRY QEQSM(S,A,M)
      CALL EQSM(S,A,M,M)
      RETURN
      ENTRY QEQDM(D,A,M)
      CALL EQDM(D,A,M,M)
      RETURN
      ENTRY QSUBSM(S,A,B,M)
      CALL SUBSM(S,A,B,M,M,M)
      RETURN
      ENTRY QADDSM(S,A,B,M)
      CALL ADDSM(S,A,B,M,M,M)
      RETURN
      ENTRY QADDDM(D,A,B,M)
      CALL ADDDM(D,A,B,M,M,M)
      RETURN
      ENTRY QEQLD(A,D,M)
      CALL EQLD(A,D,M,M)
      RETURN
      ENTRY QDETL(A,DET,M)
      CALL DETL(A,DET,M,M)
      RETURN
      ENTRY QINVM(A,M)
      CALL INVMM(A,A,M,M,M)
      RETURN
      ENTRY QLINVMM(A,B,M,MM)
      CALL LINVMM(A,B,M,MM,M,M)
      RETURN
      ENTRY QLINVMV(A,V,M)
      CALL LINVMM(A,V,M,1,M,M)
      RETURN
      ENTRY QWRTM(A,JCOLS,M)
      CALL WRTM(A,M,M,JCOLS,M)
      RETURN
      ENTRY QWRTIM(L,JCOLS,M)
      CALL WRTIM(L,M,M,JCOLS,M)
      RETURN
      END
      SUBROUTINE QMULMV(A,U,V,M)
      DIMENSION A(M,*),U(*),V(*)
      DO I=1,M
      V(I)=0.
      ENDDO
      ENTRY QMADMV(A,U,V,M)
      DO I=1,M
      DO J=1,M
      V(I)=V(I)+A(I,J)*U(J)
      ENDDO
      ENDDO
      RETURN

      ENTRY QMULTV(A,U,V,M)
      DO I=1,M
      V(I)=0.
      ENDDO
      ENTRY QMADTV(A,U,V,M)
      DO I=1,M
      DO J=1,M
      V(I)=V(I)+A(J,I)*U(J)
      ENDDO
      ENDDO
      RETURN
      END
C                                               *****************
C                                               *   NFFT1.FOR   *
C                                               *  PURSER 1994  *
C                                               *****************
C
      BLOCKDATA DAT235
      COMMON/FFT235/ LN2,LN3,LN5
      DATA LN2/0/,LN3/0/,LN5/0/
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINES CFFT, DFFT                                     C
C                                                                              C
C   Fourier analyze (CFFT) or synthesize (DFFT) a line of complex data         C
C  RB     <-> Real part of data and transform                                  C
C  QB     <-> Imaginary part of data and transform                             C
C  N      --> Number of complex data of series (product of 2's,3's,5's)        C
C  PERIOD --> Period of one cycle                                              C
C  W      <-> Coefficient array of N reals (re)computed only when              C
C             JUMBLE(0)=0                                                      C
C  JUMBLE <-> Coefficient array of N integers (re)computed only when           C
C             JUMBLE(0) set to zero on input to this routine                   C
C                                                                              C
C   For related routines, see subroutine RUMBLE (in which the                  C
C   coefficient arrays are initialized), XCFFT, YCFFT (matrix versions         C
C   of CFFT), XDFFT, YDFFT (matrix versions of DFFT), RFFT (Fourier            C
C   analysis of real data), HFFT (Fourier synthesis of real data), and         C
C   matrix counterparts, XRFFT, YRFFT, XHFFT, YHFFT.                           C
C------------------------------------------------------------------------------C
      SUBROUTINE CFFT(RB,QB,N,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(0:*),QB(0:*),W(0:*),JUMBLE(0:*)
      RFAC=PERIOD/N
      GOTO 300
      ENTRY      DFFT(RB,QB,N,PERIOD,W,JUMBLE)
      RFAC=1./PERIOD

C  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
      DO J=N/2+1,N-1
       I=N-J
       T    =RB(J)
       RB(J)=RB(I)
       RB(I)=T
       T    =QB(J)
       QB(J)=QB(I)
       QB(I)=T
      ENDDO

300   IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

C  SCALE AND PERMUTE THE DATA:
      RB(0)=RB(0)*RFAC
      QB(0)=QB(0)*RFAC
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        T1=RB(I)
        RB(I)=RB(J)*RFAC
        RB(J)=T1
        T1=QB(I)
        QB(I)=QB(J)*RFAC
        QB(J)=T1
       ELSE
        RB(I)=RB(I)*RFAC
        QB(I)=QB(I)*RFAC
       ENDIF
      ENDDO

C  TRANSFORM THE DATA:
      MA=1
      MB=N
C  RADIX 4
      LS=0
      DO L=2,LN2,2
      MB=MB/4
      LS=L
      MA4=MA*4
      MB2=MB*2
      DO J=0,MA-1
      JMB=J*MB
      JMB2=J*MB2
      RF1=W(JMB)
      QF1=W(NH+JMB)
      RF2=W(JMB2)
      QF2=W(NH+JMB2)
      RF3=RF1*RF2-QF1*QF2
      QF3=RF1*QF2+QF1*RF2
      DO I=0,NM,MA4
      K0=I+J
      K1=K0+MA
      K2=K1+MA
      K3=K2+MA
      R0=RB(K0)
      R1=RB(K1)
      R2=RB(K2)
      R3=RB(K3)
      Q0=QB(K0)
      Q1=QB(K1)
      Q2=QB(K2)
      Q3=QB(K3)
      T1=R3*RF3-Q3*QF3  ! Q13
      Q3=R3*QF3+Q3*RF3  ! R13
      R3=R2*RF1-Q2*QF1  ! R12
      R2=R2*QF1+Q2*RF1  ! Q12
      Q2=Q3-R2          ! R23
      R2=Q3+R2          ! Q22
      Q3=R3+T1          ! R22
      T1=R3-T1          ! Q23
      R3=R1*RF2-Q1*QF2  ! R11
      Q1=R1*QF2+Q1*RF2  ! Q11
      R1=R0-R3          ! R21
      R0=R0+R3          ! R20
      RB(K3)=R1-Q2      ! R3
      RB(K1)=R1+Q2      ! R1
      Q2=Q0+Q1          ! Q20
      Q1=Q0-Q1          ! Q21
      QB(K0)=Q2+R2      ! Q0
      QB(K2)=Q2-R2      ! Q2
      RB(K2)=R0-Q3      ! R2
      RB(K0)=R0+Q3      ! R0
      QB(K3)=Q1-T1      ! Q3
      QB(K1)=Q1+T1      ! Q1
      ENDDO
      ENDDO
      MA=MA4
      ENDDO
      IF(LS.NE.LN2)THEN
C  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
      JMB=J*MB
      DO I=0,NM,MA2
      K0=J+I
      K1=K0+MA
      RF1=W(JMB)
      QF1=W(NH+JMB)
      R0=RB(K0)
      Q0=QB(K0)
      R1=RB(K1)
      Q1=QB(K1)
      T1=R1*QF1+Q1*RF1 ! Q11
      Q1=R1*RF1-Q1*QF1 ! R11
      RB(K1)=R0-Q1     ! R1
      RB(K0)=R0+Q1     ! R0
      QB(K1)=Q0-T1     ! Q1
      QB(K0)=Q0+T1     ! Q0
      ENDDO
      ENDDO
      MA=MA2
      ENDIF
C  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
      MB=MB/3
      MA3=MA*3
      DO J=0,MA-1
      JMB=J*MB
      RF1=W(JMB)
      QF1=W(NH+JMB)
      RF2=RF1*RF1-QF1*QF1
      QF2=2*RF1*QF1
      DO I=0,NM,MA3
      K0=I+J
      K1=K0+MA
      K2=K1+MA
      R1=RB(K1)
      Q1=QB(K1)
      R2=RB(K2)
      Q2=QB(K2)
      T1=R2*QF2+Q2*RF2  ! R12
      Q2=R2*RF2-Q2*QF2  ! Q12
      R2=R1*QF1+Q1*RF1  ! Q11
      R1=R1*RF1-Q1*QF1  ! R11
      Q1=R2+T1          ! Q21
      R2=(R2-T1)*QEP    ! R22
      T1=R1+Q2          ! R21
      R1=(R1-Q2)*QEP    ! Q22
      RB(K0)=RB(K0)+T1  ! R0
      QB(K0)=QB(K0)+Q1  ! Q0
      T1=RB(K0)-T1*REC  ! R21
      Q1=QB(K0)-Q1*REC  ! Q21
      QB(K2)=Q1-R1      ! Q2
      QB(K1)=Q1+R1      ! Q1
      RB(K1)=T1-R2      ! R1
      RB(K2)=T1+R2      ! R2
      ENDDO
      ENDDO
      MA=MA3
      ENDDO

      IF(LN5.GT.0)THEN
C  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
      MB=MB/5
      MA5=MA*5
      DO J=0,MA-1
      JMB=J*MB
      JMB2=JMB*2
      RF1=W(JMB)
      QF1=W(NH+JMB)
      RF2=W(JMB2)
      QF2=W(NH+JMB2)
      RF3=RF1*RF2-QF1*QF2
      QF3=RF1*QF2+QF1*RF2
      RF4=RF2*RF2-QF2*QF2
      QF4=2*RF2*QF2
      DO I=0,NM,MA5
      K0=I+J
      K1=K0+MA
      K2=K1+MA
      K3=K2+MA
      K4=K3+MA
      R1=RB(K1)
      R2=RB(K2)
      R3=RB(K3)
      R4=RB(K4)
      Q1=QB(K1)
      Q2=QB(K2)
      Q3=QB(K3)
      Q4=QB(K4)
      T1=R1*QF1+Q1*RF1        ! Q11
      R1=R1*RF1-Q1*QF1        ! R11
      Q1=R4*RF4-Q4*QF4        ! Q14
      R4=R4*QF4+Q4*RF4        ! R14
      Q4=R1-Q1                ! Q24
      R1=R1+Q1                ! R21
      Q1=T1+R4                ! Q21
      R4=T1-R4                ! R24
      T1=R3*RF3-Q3*QF3        ! Q13
      R3=R3*QF3+Q3*RF3        ! R13
      Q3=R2*QF2+Q2*RF2        ! Q12
      R2=R2*RF2-Q2*QF2        ! R12
      Q2=Q3+R3                ! Q22
      R3=Q3-R3                ! R23
      Q3=R2-T1                ! Q23
      R2=R2+T1                ! R22
      RB(K0)=RB(K0)+R1+R2     ! R0
      QB(K0)=QB(K0)+Q1+Q2     ! Q0
      T1=       R4*QZE+R3*QET ! R34
      R3=       R3*QZE-R4*QET ! R33
      R4=RB(K0)-R2*RZC-R1*REC ! R32
      R1=RB(K0)-R1*RZC-R2*REC ! R31
      RB(K2)=R4+R3            ! R2
      RB(K3)=R4-R3            ! R3
      RB(K4)=R1+T1            ! R4
      RB(K1)=R1-T1            ! R1
      T1=QB(K0)-Q1*RZC-Q2*REC ! Q31
      Q2=QB(K0)-Q2*RZC-Q1*REC ! Q32
      Q1=       Q3*QZE-Q4*QET ! Q33
      Q4=       Q4*QZE+Q3*QET ! Q34
      QB(K3)=Q2+Q1            ! Q3
      QB(K2)=Q2-Q1            ! Q2
      QB(K1)=T1+Q4            ! Q1
      QB(K4)=T1-Q4            ! Q4
      ENDDO
      ENDDO
      MA=MA5
      ENDDO
      ENDIF
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE GET235                                          C
C                                                                              C
C   Analyze the facorization of N in terms of 2's, 3's, and 5's and verify     C
C   that these are the only prime factors. Set LN2, LN3, LN5 respectively to   C
C   the number of powers of 2, 3, 5.                                           C
C  N      --> Number of data along the line of Fourier transformation          C
C------------------------------------------------------------------------------C
      SUBROUTINE GET235(N)
      COMMON/FFT235/ LN2,LN3,LN5
      NKT=N
      NM=N-1
      NH=N/2
      LN2=-1
400   NK=NKT
      NKT=NK/2
      LN2=LN2+1
      IF(2*NKT.EQ.NK)GOTO 400
      NKT=NK
      LN3=-1
401   NK=NKT
      NKT=NK/3
      LN3=LN3+1
      IF(3*NKT.EQ.NK)GOTO 401
      NKT=NK
      LN5=-1
402   NK=NKT
      NKT=NK/5
      LN5=LN5+1
      IF(5*NKT.EQ.NK)GOTO 402
ctest      IF(NK.NE.1)STOP'N CANNOT BE FACTORED INTO ONLY 2s, 3s, 5s'
ctest
!      IF(NK.NE.1) Then

!         Stop
!      End If
ctest
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE RUMBLE                                          C
C                                                                              C
C   Initialize coefficient arrays JUMBLE and TUMBLE for use with the           C
C   fast-Fourier-transform routines when the number N of data has the prime-   C
C   factorization 2**LN2*3**LN3*5**LN5                                         C
C                                                                              C
C  JUMBLE <-- Permutation of N data encoded as a sequence of transpositions    C
C  TUMBLE <-- Trigonometric coefficients for use in FFT. The first half are    C
C             the cosines, the second half the sines, of uniformly increasing  C
C             relevant angles                                                  C
C------------------------------------------------------------------------------C
      SUBROUTINE RUMBLE(JUMBLE,TUMBLE)
      PARAMETER(ML=20)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION ND(ML),MD(ML),JUMBLE(0:*),TUMBLE(0:*)
      LN=LN2+LN3+LN5
      N=2**LN2*3**LN3*5**LN5
      NM=N-1
      NH=N/2
      PI2ON=8.*ATAN(1.)/N
      DO I=0,NH-1
      ANG=PI2ON*I
      TUMBLE(I)=COS(ANG)
      TUMBLE(I+NH)=SIN(ANG)
      ENDDO
      ID=1
      IS=0
      DO I=1,LN5
       IS=IS+1
       MD(IS)=ID
       ID=ID*5
      ENDDO
      DO I=1,LN3
       IS=IS+1
       MD(IS)=ID
       ID=ID*3
      ENDDO
      DO I=1,LN2
       IS=IS+1
       MD(IS)=ID
       ID=ID*2
      ENDDO
      ID=1
      DO I=1,LN2
       ND(IS)=ID
       ID=ID*2
       IS=IS-1
      ENDDO
      DO I=1,LN3
       ND(IS)=ID
       ID=ID*3
       IS=IS-1
      ENDDO
      DO I=1,LN5
       ND(IS)=ID
       ID=ID*5
       IS=IS-1
      ENDDO
      JUMBLE(0)=N
      DO I=1,NM
       IR=I
       J=0
       DO L=1,LN
        KD=IR/ND(L)
        IR=IR-KD*ND(L)
        J=J+KD*MD(L)
       ENDDO
       JUMBLE(I)=J
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.LT.I)THEN
400        J=JUMBLE(J)
           IF(J.LT.I)GOTO 400
        JUMBLE(I)=J
       ENDIF
      ENDDO
      RETURN
      END

      SUBROUTINE XCFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(0:NCB-1,*),QB(0:NCB-1,*),W(0:*),JUMBLE(0:*)
      RFAC=PERIOD/N
C      GOTO 300
C
C      ENTRY     XDFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE)
C      RFAC=1./PERIOD
C
CC  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
C      DO J=N/2+1,N-1
C       I=N-J
C       DO IY=1,NY
C       T       =RB(J,IY)
C       RB(J,IY)=RB(I,IY)
C       RB(I,IY)=T
C       T       =QB(J,IY)
C       QB(J,IY)=QB(I,IY)
C       QB(I,IY)=T
C       ENDDO
C      ENDDO
C
      IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

C  SCALE AND PERMUTE THE DATA:
      DO IY=1,NY
       RB(0,IY)=RB(0,IY)*RFAC
       QB(0,IY)=QB(0,IY)*RFAC
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        DO IY=1,NY
         T=RB(I,IY)
         RB(I,IY)=RB(J,IY)*RFAC
         RB(J,IY)=T
         T=QB(I,IY)
         QB(I,IY)=QB(J,IY)*RFAC
         QB(J,IY)=T
        ENDDO
       ELSE
        DO IY=1,NY
         RB(I,IY)=RB(I,IY)*RFAC
         QB(I,IY)=QB(I,IY)*RFAC
        ENDDO
       ENDIF
      ENDDO
C  TRANSFORM THE DATA:
      MA=1
      MB=N
C  RADIX 4
      LS=0
      DO L=2,LN2,2
       MB=MB/4
       LS=L
       MA4=MA*4
       MB2=MB*2
       DO J=0,MA-1
        JMB=J*MB
        JMB2=J*MB2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        DO I=0,NM,MA4
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         DO IY=1,NY
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3    ! Q13
          QB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3    ! R13
          RB(K3,IY)=RB(K2,IY)*RF1-QB(K2,IY)*QF1    ! R12
          RB(K2,IY)=RB(K2,IY)*QF1+QB(K2,IY)*RF1    ! Q12
          QB(K2,IY)=QB(K3,IY)-RB(K2,IY)            ! R23
          RB(K2,IY)=QB(K3,IY)+RB(K2,IY)            ! Q22
          QB(K3,IY)=RB(K3,IY)+T                    ! R22
          T        =RB(K3,IY)-T                    ! Q23
          RB(K3,IY)=RB(K1,IY)*RF2-QB(K1,IY)*QF2    ! R11
          QB(K1,IY)=RB(K1,IY)*QF2+QB(K1,IY)*RF2    ! Q11
          RB(K1,IY)=RB(K0,IY)-RB(K3,IY)            ! R21
          RB(K0,IY)=RB(K0,IY)+RB(K3,IY)            ! R20
          RB(K3,IY)=RB(K1,IY)-QB(K2,IY)            ! R3
          RB(K1,IY)=RB(K1,IY)+QB(K2,IY)            ! R1
          QB(K2,IY)=QB(K0,IY)+QB(K1,IY)            ! Q20
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)            ! Q21
          QB(K0,IY)=QB(K2,IY)+RB(K2,IY)            ! Q0
          QB(K2,IY)=QB(K2,IY)-RB(K2,IY)            ! Q2
          RB(K2,IY)=RB(K0,IY)-QB(K3,IY)            ! R2
          RB(K0,IY)=RB(K0,IY)+QB(K3,IY)            ! R0
          QB(K3,IY)=QB(K1,IY)-T                    ! Q3
          QB(K1,IY)=QB(K1,IY)+T                    ! Q1
         ENDDO
        ENDDO
       ENDDO
       MA=MA4
      ENDDO

      IF(LS.NE.LN2)THEN
C  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
       JMB=J*MB
       DO I=0,NM,MA2
        K0=J+I
        K1=K0+MA
        RF1=W(JMB)
        QF1=W(NH+JMB)
        DO IY=1,NY
         T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11
         QB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11
         RB(K1,IY)=RB(K0,IY)-QB(K1,IY)         ! R1
         RB(K0,IY)=RB(K0,IY)+QB(K1,IY)         ! R0
         QB(K1,IY)=QB(K0,IY)-T                 ! Q1
         QB(K0,IY)=QB(K0,IY)+T                 ! Q0
        ENDDO
       ENDDO
      ENDDO
      MA=MA2
      ENDIF
C  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
       MB=MB/3
       MA3=MA*3
       DO J=0,MA-1
        JMB=J*MB
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=RF1*RF1-QF1*QF1
        QF2=2*RF1*QF1
        DO I=0,NM,MA3
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         DO IY=1,NY
          T        =RB(K2,IY)*QF2+QB(K2,IY)*RF2     ! R12
          QB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2     ! Q12
          RB(K2,IY)=RB(K1,IY)*QF1+QB(K1,IY)*RF1     ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1     ! R11
          QB(K1,IY)=RB(K2,IY)+T                     ! Q21
          RB(K2,IY)=(RB(K2,IY)-T)*QEP               ! R22
          T        =RB(K1,IY)+QB(K2,IY)             ! R21
          RB(K1,IY)=(RB(K1,IY)-QB(K2,IY))*QEP       ! Q22
          RB(K0,IY)=RB(K0,IY)+T                     ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)             ! Q0
          T        =RB(K0,IY)-T*REC                 ! R21
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)*REC         ! Q21
          QB(K2,IY)=QB(K1,IY)-RB(K1,IY)             ! Q2
          QB(K1,IY)=QB(K1,IY)+RB(K1,IY)             ! Q1
          RB(K1,IY)=T        -RB(K2,IY)             ! R1
          RB(K2,IY)=T        +RB(K2,IY)             ! R2
         ENDDO
        ENDDO
       ENDDO
       MA=MA3
      ENDDO

      IF(LN5.GT.0)THEN
C  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
       MB=MB/5
       MA5=MA*5
       DO J=0,MA-1
        JMB=J*MB
        JMB2=JMB*2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        RF4=RF2*RF2-QF2*QF2
        QF4=2*RF2*QF2
        DO I=0,NM,MA5
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         K4=K3+MA
         DO IY=1,NY
          T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1           ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1           ! R11
          QB(K1,IY)=RB(K4,IY)*RF4-QB(K4,IY)*QF4           ! Q14
          RB(K4,IY)=RB(K4,IY)*QF4+QB(K4,IY)*RF4           ! R14
          QB(K4,IY)=RB(K1,IY)-QB(K1,IY)                   ! Q24
          RB(K1,IY)=RB(K1,IY)+QB(K1,IY)                   ! R21
          QB(K1,IY)=T        +RB(K4,IY)                   ! Q21
          RB(K4,IY)=T        -RB(K4,IY)                   ! R24
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3           ! Q13
          RB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3           ! R13
          QB(K3,IY)=RB(K2,IY)*QF2+QB(K2,IY)*RF2           ! Q12
          RB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2           ! R12
          QB(K2,IY)=QB(K3,IY)+RB(K3,IY)                   ! Q22
          RB(K3,IY)=QB(K3,IY)-RB(K3,IY)                   ! R23
          QB(K3,IY)=RB(K2,IY)-T                           ! Q23
          RB(K2,IY)=RB(K2,IY)+T                           ! R22
          RB(K0,IY)=RB(K0,IY)+RB(K1,IY)+RB(K2,IY)         ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)+QB(K2,IY)         ! Q0
          T        =RB(K4,IY)*QZE+RB(K3,IY)*QET           ! R34
          RB(K3,IY)=RB(K3,IY)*QZE-RB(K4,IY)*QET           ! R33
          RB(K4,IY)=RB(K0,IY)-RB(K2,IY)*RZC-RB(K1,IY)*REC ! R32
          RB(K1,IY)=RB(K0,IY)-RB(K1,IY)*RZC-RB(K2,IY)*REC ! R31
          RB(K2,IY)=RB(K4,IY)+RB(K3,IY)                   ! R2
          RB(K3,IY)=RB(K4,IY)-RB(K3,IY)                   ! R3
          RB(K4,IY)=RB(K1,IY)+T                           ! R4
          RB(K1,IY)=RB(K1,IY)-T                           ! R1
          T        =QB(K0,IY)-QB(K1,IY)*RZC-QB(K2,IY)*REC ! Q31
          QB(K2,IY)=QB(K0,IY)-QB(K2,IY)*RZC-QB(K1,IY)*REC ! Q32
          QB(K1,IY)=QB(K3,IY)*QZE-QB(K4,IY)*QET           ! Q33
          QB(K4,IY)=QB(K4,IY)*QZE+QB(K3,IY)*QET           ! Q34
          QB(K3,IY)=QB(K2,IY)+QB(K1,IY)                   ! Q3
          QB(K2,IY)=QB(K2,IY)-QB(K1,IY)                   ! Q2
          QB(K1,IY)=T        +QB(K4,IY)                   ! Q1
          QB(K4,IY)=T        -QB(K4,IY)                   ! Q4
         ENDDO
        ENDDO
       ENDDO
       MA=MA5
      ENDDO
      ENDIF
      RETURN
      END

      SUBROUTINE XDFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(0:NCB-1,*),QB(0:NCB-1,*),W(0:*),JUMBLE(0:*)
      RFAC=1./PERIOD

C  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
      DO J=N/2+1,N-1
       I=N-J
       DO IY=1,NY
        T       =RB(J,IY)
        RB(J,IY)=RB(I,IY)
        RB(I,IY)=T
        T       =QB(J,IY)
        QB(J,IY)=QB(I,IY)
        QB(I,IY)=T
       ENDDO
      ENDDO

      IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

C  SCALE AND PERMUTE THE DATA:
      DO IY=1,NY
       RB(0,IY)=RB(0,IY)*RFAC
       QB(0,IY)=QB(0,IY)*RFAC
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        DO IY=1,NY
         T=RB(I,IY)
         RB(I,IY)=RB(J,IY)*RFAC
         RB(J,IY)=T
         T=QB(I,IY)
         QB(I,IY)=QB(J,IY)*RFAC
         QB(J,IY)=T
        ENDDO
       ELSE
        DO IY=1,NY
         RB(I,IY)=RB(I,IY)*RFAC
         QB(I,IY)=QB(I,IY)*RFAC
        ENDDO
       ENDIF
      ENDDO
C  TRANSFORM THE DATA:
      MA=1
      MB=N
C  RADIX 4
      LS=0
      DO L=2,LN2,2
       MB=MB/4
       LS=L
       MA4=MA*4
       MB2=MB*2
       DO J=0,MA-1
        JMB=J*MB
        JMB2=J*MB2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        DO I=0,NM,MA4
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         DO IY=1,NY
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3    ! Q13
          QB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3    ! R13
          RB(K3,IY)=RB(K2,IY)*RF1-QB(K2,IY)*QF1    ! R12
          RB(K2,IY)=RB(K2,IY)*QF1+QB(K2,IY)*RF1    ! Q12
          QB(K2,IY)=QB(K3,IY)-RB(K2,IY)            ! R23
          RB(K2,IY)=QB(K3,IY)+RB(K2,IY)            ! Q22
          QB(K3,IY)=RB(K3,IY)+T                    ! R22
          T        =RB(K3,IY)-T                    ! Q23
          RB(K3,IY)=RB(K1,IY)*RF2-QB(K1,IY)*QF2    ! R11
          QB(K1,IY)=RB(K1,IY)*QF2+QB(K1,IY)*RF2    ! Q11
          RB(K1,IY)=RB(K0,IY)-RB(K3,IY)            ! R21
          RB(K0,IY)=RB(K0,IY)+RB(K3,IY)            ! R20
          RB(K3,IY)=RB(K1,IY)-QB(K2,IY)            ! R3
          RB(K1,IY)=RB(K1,IY)+QB(K2,IY)            ! R1
          QB(K2,IY)=QB(K0,IY)+QB(K1,IY)            ! Q20
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)            ! Q21
          QB(K0,IY)=QB(K2,IY)+RB(K2,IY)            ! Q0
          QB(K2,IY)=QB(K2,IY)-RB(K2,IY)            ! Q2
          RB(K2,IY)=RB(K0,IY)-QB(K3,IY)            ! R2
          RB(K0,IY)=RB(K0,IY)+QB(K3,IY)            ! R0
          QB(K3,IY)=QB(K1,IY)-T                    ! Q3
          QB(K1,IY)=QB(K1,IY)+T                    ! Q1
         ENDDO
        ENDDO
       ENDDO
       MA=MA4
      ENDDO

      IF(LS.NE.LN2)THEN
C  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
       JMB=J*MB
       DO I=0,NM,MA2
        K0=J+I
        K1=K0+MA
        RF1=W(JMB)
        QF1=W(NH+JMB)
        DO IY=1,NY
         T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11
         QB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11
         RB(K1,IY)=RB(K0,IY)-QB(K1,IY)         ! R1
         RB(K0,IY)=RB(K0,IY)+QB(K1,IY)         ! R0
         QB(K1,IY)=QB(K0,IY)-T                 ! Q1
         QB(K0,IY)=QB(K0,IY)+T                 ! Q0
        ENDDO
       ENDDO
      ENDDO
      MA=MA2
      ENDIF
C  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
       MB=MB/3
       MA3=MA*3
       DO J=0,MA-1
        JMB=J*MB
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=RF1*RF1-QF1*QF1
        QF2=2*RF1*QF1
        DO I=0,NM,MA3
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         DO IY=1,NY
          T        =RB(K2,IY)*QF2+QB(K2,IY)*RF2     ! R12
          QB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2     ! Q12
          RB(K2,IY)=RB(K1,IY)*QF1+QB(K1,IY)*RF1     ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1     ! R11
          QB(K1,IY)=RB(K2,IY)+T                     ! Q21
          RB(K2,IY)=(RB(K2,IY)-T)*QEP               ! R22
          T        =RB(K1,IY)+QB(K2,IY)             ! R21
          RB(K1,IY)=(RB(K1,IY)-QB(K2,IY))*QEP       ! Q22
          RB(K0,IY)=RB(K0,IY)+T                     ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)             ! Q0
          T        =RB(K0,IY)-T*REC                 ! R21
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)*REC         ! Q21
          QB(K2,IY)=QB(K1,IY)-RB(K1,IY)             ! Q2
          QB(K1,IY)=QB(K1,IY)+RB(K1,IY)             ! Q1
          RB(K1,IY)=T        -RB(K2,IY)             ! R1
          RB(K2,IY)=T        +RB(K2,IY)             ! R2
         ENDDO
        ENDDO
       ENDDO
       MA=MA3
      ENDDO

      IF(LN5.GT.0)THEN
C  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
       MB=MB/5
       MA5=MA*5
       DO J=0,MA-1
        JMB=J*MB
        JMB2=JMB*2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        RF4=RF2*RF2-QF2*QF2
        QF4=2*RF2*QF2
        DO I=0,NM,MA5
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         K4=K3+MA
         DO IY=1,NY
          T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1           ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1           ! R11
          QB(K1,IY)=RB(K4,IY)*RF4-QB(K4,IY)*QF4           ! Q14
          RB(K4,IY)=RB(K4,IY)*QF4+QB(K4,IY)*RF4           ! R14
          QB(K4,IY)=RB(K1,IY)-QB(K1,IY)                   ! Q24
          RB(K1,IY)=RB(K1,IY)+QB(K1,IY)                   ! R21
          QB(K1,IY)=T        +RB(K4,IY)                   ! Q21
          RB(K4,IY)=T        -RB(K4,IY)                   ! R24
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3           ! Q13
          RB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3           ! R13
          QB(K3,IY)=RB(K2,IY)*QF2+QB(K2,IY)*RF2           ! Q12
          RB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2           ! R12
          QB(K2,IY)=QB(K3,IY)+RB(K3,IY)                   ! Q22
          RB(K3,IY)=QB(K3,IY)-RB(K3,IY)                   ! R23
          QB(K3,IY)=RB(K2,IY)-T                           ! Q23
          RB(K2,IY)=RB(K2,IY)+T                           ! R22
          RB(K0,IY)=RB(K0,IY)+RB(K1,IY)+RB(K2,IY)         ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)+QB(K2,IY)         ! Q0
          T        =RB(K4,IY)*QZE+RB(K3,IY)*QET           ! R34
          RB(K3,IY)=RB(K3,IY)*QZE-RB(K4,IY)*QET           ! R33
          RB(K4,IY)=RB(K0,IY)-RB(K2,IY)*RZC-RB(K1,IY)*REC ! R32
          RB(K1,IY)=RB(K0,IY)-RB(K1,IY)*RZC-RB(K2,IY)*REC ! R31
          RB(K2,IY)=RB(K4,IY)+RB(K3,IY)                   ! R2
          RB(K3,IY)=RB(K4,IY)-RB(K3,IY)                   ! R3
          RB(K4,IY)=RB(K1,IY)+T                           ! R4
          RB(K1,IY)=RB(K1,IY)-T                           ! R1
          T        =QB(K0,IY)-QB(K1,IY)*RZC-QB(K2,IY)*REC ! Q31
          QB(K2,IY)=QB(K0,IY)-QB(K2,IY)*RZC-QB(K1,IY)*REC ! Q32
          QB(K1,IY)=QB(K3,IY)*QZE-QB(K4,IY)*QET           ! Q33
          QB(K4,IY)=QB(K4,IY)*QZE+QB(K3,IY)*QET           ! Q34
          QB(K3,IY)=QB(K2,IY)+QB(K1,IY)                   ! Q3
          QB(K2,IY)=QB(K2,IY)-QB(K1,IY)                   ! Q2
          QB(K1,IY)=T        +QB(K4,IY)                   ! Q1
          QB(K4,IY)=T        -QB(K4,IY)                   ! Q4
         ENDDO
        ENDDO
       ENDDO
       MA=MA5
      ENDDO
      ENDIF
      RETURN
      END

      SUBROUTINE YCFFT(RB,QB,N,NX,NDX,NCB,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(NCB,0:*),QB(NCB,0:*),W(0:*),JUMBLE(0:*)
      RFAC=PERIOD/N
      GOTO 300

      ENTRY      YDFFT(RB,QB,N,NX,NDX,NCB,PERIOD,W,JUMBLE)
      RFAC=1./PERIOD

C  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
      DO J=N/2+1,N-1
       I=N-J
       DO IX=1,NX,NDX
        T       =RB(IX,J)
        RB(IX,J)=RB(IX,I)
        RB(IX,I)=T
        T       =QB(IX,J)
        QB(IX,J)=QB(IX,I)
        QB(IX,I)=T
       ENDDO
      ENDDO

300   IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

C  SCALE AND PERMUTE THE DATA:
      DO IX=1,NX,NDX
       RB(IX,0)=RB(IX,0)*RFAC
       QB(IX,0)=QB(IX,0)*RFAC
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        DO IX=1,NX,NDX
         T=RB(IX,I)
         RB(IX,I)=RB(IX,J)*RFAC
         RB(IX,J)=T
         T=QB(IX,I)
         QB(IX,I)=QB(IX,J)*RFAC
         QB(IX,J)=T
        ENDDO
       ELSE
        DO IX=1,NX,NDX
         RB(IX,I)=RB(IX,I)*RFAC
         QB(IX,I)=QB(IX,I)*RFAC
        ENDDO
       ENDIF
      ENDDO

C  TRANSFORM THE DATA:
      MA=1
      MB=N
C  RADIX 4
      LS=0
      DO L=2,LN2,2
       MB=MB/4
       LS=L
       MA4=MA*4
       MB2=MB*2
       DO J=0,MA-1
        JMB=J*MB
        JMB2=J*MB2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        DO I=0,NM,MA4
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         DO IX=1,NX,NDX
          T        =RB(IX,K3)*RF3-QB(IX,K3)*QF3     ! Q13
          QB(IX,K3)=RB(IX,K3)*QF3+QB(IX,K3)*RF3     ! R13
          RB(IX,K3)=RB(IX,K2)*RF1-QB(IX,K2)*QF1     ! R12
          RB(IX,K2)=RB(IX,K2)*QF1+QB(IX,K2)*RF1     ! Q12
          QB(IX,K2)=QB(IX,K3)-RB(IX,K2)             ! R23
          RB(IX,K2)=QB(IX,K3)+RB(IX,K2)             ! Q22
          QB(IX,K3)=RB(IX,K3)+T                     ! R22
          T        =RB(IX,K3)-T                     ! Q23
          RB(IX,K3)=RB(IX,K1)*RF2-QB(IX,K1)*QF2     ! R11
          QB(IX,K1)=RB(IX,K1)*QF2+QB(IX,K1)*RF2     ! Q11
          RB(IX,K1)=RB(IX,K0)-RB(IX,K3)             ! R21
          RB(IX,K0)=RB(IX,K0)+RB(IX,K3)             ! R20
          RB(IX,K3)=RB(IX,K1)-QB(IX,K2)             ! R3
          RB(IX,K1)=RB(IX,K1)+QB(IX,K2)             ! R1
          QB(IX,K2)=QB(IX,K0)+QB(IX,K1)             ! Q20
          QB(IX,K1)=QB(IX,K0)-QB(IX,K1)             ! Q21
          QB(IX,K0)=QB(IX,K2)+RB(IX,K2)             ! Q0
          QB(IX,K2)=QB(IX,K2)-RB(IX,K2)             ! Q2
          RB(IX,K2)=RB(IX,K0)-QB(IX,K3)             ! R2
          RB(IX,K0)=RB(IX,K0)+QB(IX,K3)             ! R0
          QB(IX,K3)=QB(IX,K1)-T                     ! Q3
          QB(IX,K1)=QB(IX,K1)+T                     ! Q1
         ENDDO
        ENDDO
       ENDDO
       MA=MA4
      ENDDO
      IF(LS.NE.LN2)THEN

C  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
       JMB=J*MB
       DO I=0,NM,MA2
        K0=J+I
        K1=K0+MA
        RF1=W(JMB)
        QF1=W(NH+JMB)
        DO IX=1,NX,NDX
         T        =RB(IX,K1)*QF1+QB(IX,K1)*RF1 ! Q11
         QB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1 ! R11
         RB(IX,K1)=RB(IX,K0)-QB(IX,K1)         ! R1
         RB(IX,K0)=RB(IX,K0)+QB(IX,K1)         ! R0
         QB(IX,K1)=QB(IX,K0)-T                 ! Q1
         QB(IX,K0)=QB(IX,K0)+T                 ! Q0
        ENDDO
       ENDDO
      ENDDO
      MA=MA2
      ENDIF

C  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
       MB=MB/3
       MA3=MA*3
       DO J=0,MA-1
        JMB=J*MB
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=RF1*RF1-QF1*QF1
        QF2=2*RF1*QF1
        DO I=0,NM,MA3
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         T        =RB(IX,K2)*QF2+QB(IX,K2)*RF2     ! R12
         QB(IX,K2)=RB(IX,K2)*RF2-QB(IX,K2)*QF2     ! Q12
         RB(IX,K2)=RB(IX,K1)*QF1+QB(IX,K1)*RF1     ! Q11
         RB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1     ! R11
         QB(IX,K1)=RB(IX,K2)+T                     ! Q21
         RB(IX,K2)=(RB(IX,K2)-T        )*QEP       ! R22
         T        =RB(IX,K1)+QB(IX,K2)             ! R21
         RB(IX,K1)=(RB(IX,K1)-QB(IX,K2))*QEP       ! Q22
         RB(IX,K0)=RB(IX,K0)+T                     ! R0
         QB(IX,K0)=QB(IX,K0)+QB(IX,K1)             ! Q0
         T        =RB(IX,K0)-T        *REC         ! R21
         QB(IX,K1)=QB(IX,K0)-QB(IX,K1)*REC         ! Q21
         QB(IX,K2)=QB(IX,K1)-RB(IX,K1)             ! Q2
         QB(IX,K1)=QB(IX,K1)+RB(IX,K1)             ! Q1
         RB(IX,K1)=T        -RB(IX,K2)             ! R1
         RB(IX,K2)=T        +RB(IX,K2)             ! R2
        ENDDO
       ENDDO
       MA=MA3
      ENDDO
      IF(LN5.GT.0)THEN

C  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
       MB=MB/5
       MA5=MA*5
       DO J=0,MA-1
        JMB=J*MB
        JMB2=JMB*2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        RF4=RF2*RF2-QF2*QF2
        QF4=2*RF2*QF2
        DO I=0,NM,MA5
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         K4=K3+MA
         T        =RB(IX,K1)*QF1+QB(IX,K1)*RF1           ! Q11
         RB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1           ! R11
         QB(IX,K1)=RB(IX,K4)*RF4-QB(IX,K4)*QF4           ! Q14
         RB(IX,K4)=RB(IX,K4)*QF4+QB(IX,K4)*RF4           ! R14
         QB(IX,K4)=RB(IX,K1)-QB(IX,K1)                   ! Q24
         RB(IX,K1)=RB(IX,K1)+QB(IX,K1)                   ! R21
         QB(IX,K1)=T        +RB(IX,K4)                   ! Q21
         RB(IX,K4)=T        -RB(IX,K4)                   ! R24
         T        =RB(IX,K3)*RF3-QB(IX,K3)*QF3           ! Q13
         RB(IX,K3)=RB(IX,K3)*QF3+QB(IX,K3)*RF3           ! R13
         QB(IX,K3)=RB(IX,K2)*QF2+QB(IX,K2)*RF2           ! Q12
         RB(IX,K2)=RB(IX,K2)*RF2-QB(IX,K2)*QF2           ! R12
         QB(IX,K2)=QB(IX,K3)+RB(IX,K3)                   ! Q22
         RB(IX,K3)=QB(IX,K3)-RB(IX,K3)                   ! R23
         QB(IX,K3)=RB(IX,K2)-T                           ! Q23
         RB(IX,K2)=RB(IX,K2)+T                           ! R22
         RB(IX,K0)=RB(IX,K0)+RB(IX,K1)+RB(IX,K2)         ! R0
         QB(IX,K0)=QB(IX,K0)+QB(IX,K1)+QB(IX,K2)         ! Q0
         T        =RB(IX,K4)*QZE+RB(IX,K3)*QET           ! R34
         RB(IX,K3)=RB(IX,K3)*QZE-RB(IX,K4)*QET           ! R33
         RB(IX,K4)=RB(IX,K0)-RB(IX,K2)*RZC-RB(IX,K1)*REC ! R32
         RB(IX,K1)=RB(IX,K0)-RB(IX,K1)*RZC-RB(IX,K2)*REC ! R31
         RB(IX,K2)=RB(IX,K4)+RB(IX,K3)                   ! R2
         RB(IX,K3)=RB(IX,K4)-RB(IX,K3)                   ! R3
         RB(IX,K4)=RB(IX,K1)+T                           ! R4
         RB(IX,K1)=RB(IX,K1)-T                           ! R1
         T        =QB(IX,K0)-QB(IX,K1)*RZC-QB(IX,K2)*REC ! Q31
         QB(IX,K2)=QB(IX,K0)-QB(IX,K2)*RZC-QB(IX,K1)*REC ! Q32
         QB(IX,K1)=QB(IX,K3)*QZE-QB(IX,K4)*QET           ! Q33
         QB(IX,K4)=QB(IX,K4)*QZE+QB(IX,K3)*QET           ! Q34
         QB(IX,K3)=QB(IX,K2)+QB(IX,K1)                   ! Q3
         QB(IX,K2)=QB(IX,K2)-QB(IX,K1)                   ! Q2
         QB(IX,K1)=T        +QB(IX,K4)                   ! Q1
         QB(IX,K4)=T        -QB(IX,K4)                   ! Q4
        ENDDO
       ENDDO
       MA=MA5
      ENDDO
      ENDIF
      RETURN
      END
