*******************************************************************
      SUBROUTINE GAUSSL(X,W,M)

      DOUBLE PRECISION Z,Z1,P1,P2,P3,PP,EPS
      DIMENSION X(M),W(M)
      PARAMETER (EPS=3.D-14)

c   this is a slightly modified version of a program in Numerical Recipes
c
c   input:
c      m    = number of gaussian latitudes between pole and equator
c   output:
c      x(m) = sin(gaussian latitude) 
c      w(m) = weights in gaussian quadrature (sum should equal 1.0)

      N=2*M

      DO 12 I=1,M
        Z=COS(3.141592654D0*(I-.25D0)/(N+.5D0))

    1   CONTINUE

        P1=1.D0
        P2=0.D0

        DO 11 J=1,N
          P3=P2
          P2=P1
          P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
   11   CONTINUE

        PP=N*(Z*P1-P2)/(Z*Z-1.D0)
        Z1=Z
        Z=Z1-P1/PP

        IF(ABS(Z-Z1).GT.EPS) GO TO 1

        X(I)=Z
        W(I)=2.D0/((1.D0-Z*Z)*PP*PP)

   12 CONTINUE

      RETURN
      END
*******************************************************************
      SUBROUTINE NCWG2(XM1,XM2,X1,X2)     

      INCLUDE 'param.h'

      DIMENSION X1(IX),X2(IX)
      DIMENSION Z(2,IX),XM1(2,NWAVES),XM2(2,NWAVES)
      COMMON/NC/ WORK(1000)

      DO 1 M=1,NWAVES
        Z(1,M)=XM1(1,M)-XM2(2,M)
        Z(2,M)=XM1(2,M)+XM2(1,M)
    1 CONTINUE

      MM1=NWAVES+1
      N2=IX+2

      DO 2 M=MM1,IX
        Z(1,M)=XM1(1,N2-M)+XM2(2,N2-M)
        Z(2,M)=-XM1(2,N2-M)+XM2(1,N2-M)
    2 CONTINUE

      CALL CFFTB(IX,Z,WORK)

      DO 3 I=1,IX
        X1(I)=Z(1,I)
        X2(I)=Z(2,I)
    3 CONTINUE

      RETURN
      END
***************************************************************
      SUBROUTINE NCGW2(X1,X2,XM1,XM2)

      INCLUDE 'param.h'

      DIMENSION X1(IX),X2(IX)
      DIMENSION Z(2,IX),XM1(2,NWAVES),XM2(2,NWAVES)
      COMMON/NC/WORK(1000)

      ANORM=1./FLOAT(IX)
      BNORM=0.5/FLOAT(IX)

      DO 1 I=1,IX
        Z(1,I)=X1(I)
        Z(2,I)=X2(I)  
    1 CONTINUE

      CALL CFFTF(IX,Z,WORK)

      N2=IX+2
      XM1(1,1)=Z(1,1)*ANORM
      XM1(2,1)=0.
      XM2(1,1)=Z(2,1)*ANORM
      XM2(2,1)=0.

      DO 2 M=2,NWAVES
        XM1(1,M)=BNORM*(Z(1,M)+Z(1,N2-M))
        XM1(2,M)=BNORM*(Z(2,M)-Z(2,N2-M))
        XM2(1,M)=BNORM*(Z(2,M)+Z(2,N2-M))
        XM2(2,M)=BNORM*(-Z(1,M)+Z(1,N2-M))
    2 CONTINUE

      RETURN
      END
******************************************************************
      SUBROUTINE PARMTR

      INCLUDE 'param.h'
      INCLUDE 'spec.h'
                                                
      DIMENSION REPSI(MXP,JXP),CONSQ(MXP)
      DIMENSION PALY(MXP,JX),ALP(MXP,JX)
      DIMENSION WT1(IL)
      DIMENSION EMM(MXP),ELL(MXP,JXP),EL2P(MX,JXP),EPSI(MXP,JXP)

      CALL GAUSSL(SIA,WT,IY)

      AM1 = 1./A
      AM2=1./(A*A)

      DO 1 J=1,IY
        COSQR = 1.0-SIA(J)**2
        COA(J)=SQRT(COSQR)
        COAM(J)=1./COA(J)
        WT1(J)=WT(J)/COSQR
    1 CONTINUE

      DO 11 J=1,IY
        JJ=IL+1-J
        COSG(J)=COA(J)
        COSG(JJ)=COA(J)
        COSGR(J)=1./COA(J)
        COSGR(JJ)=1./COA(J)
   11 CONTINUE

      DO 2 M=1,MX
      DO 2 N=1,JX
        MM(M)=ISC*(M-1)
        LL(M,N)=MM(M)+N-1
        L2(M,N)=LL(M,N)*(LL(M,N)+1)
        EL2(M,N)=FLOAT(L2(M,N))*AM2
        EL4(M,N)=EL2(M,N)*EL2(M,N)
    2 CONTINUE

      ELM2(1,1)=0.
      DO 3 M=2,MX
      DO 3 N=1,JX
        ELM2(M,N)=1./EL2(M,N)
    3 CONTINUE
      DO 4 N=2,JX
        ELM2(1,N)=1./EL2(1,N)
    4 CONTINUE

c  for internal use in this subroutine, it is useful to have EL2=l*(l+1) extended 
C    TO JXP=JX+1

      DO 5 M=1,MX
      DO 5 N=1,JXP
        LLP=MM(M)+N-1
        L2P=LLP*(LLP+1)
        EL2P(M,N)=FLOAT(L2P)*AM2
    5 CONTINUE

c    in a sector model retaining m= 0,3,6,9,12,15, then MX = 6, MXP = 16.
c    all m values up to MXP are needed by recursion relation for Legendre polynomials
c    EPSI is a useful quantity when differentiating or generating the polynomials

      DO 10 M=1,MXP
      DO 10 N=1,JXP
        EMM(M)=FLOAT(M-1)
        ELL(M,N)=FLOAT(N+M-2)
        EMM2=EMM(M)**2
        ELL2=ELL(M,N)**2
        IF(N.EQ.JXP) THEN
          EPSI(M,N)=0.0
        ELSE IF(N.EQ.1.AND.M.EQ.1) THEN
          EPSI(M,N)=0.0
        ELSE
          EPSI(M,N)=SQRT((ELL2-EMM2)/(4.*ELL2-1.))
        ENDIF
        REPSI(M,N)=0.0
        IF(EPSI(M,N).GT.0.) REPSI(M,N)=1./EPSI(M,N)
   10 CONTINUE

c  compute quantities required by subroutines GRAD and UVSPEC

      DO 15 M=1,MX
      DO 15 N=1,JX
        M1=MM(M)
        M2=M1+1
        EL1=FLOAT(LL(M,N))
        IF(N.EQ.1) THEN
          GRADX(M)=FLOAT(M1)/A
          UVDX(M,1)=-A/FLOAT(M1+1)
          UVDYM(M,1)=0.0
        ELSE
          UVDX(M,N)=-A*FLOAT(M1)/(EL1*(EL1+1))
          GRADYM(M,N)=(EL1-1.)*EPSI(M2,N)/A
          UVDYM(M,N)=-A*EPSI(M2,N)/EL1
        ENDIF
        GRADYP(M,N)=(EL1+2.)*EPSI(M2,N+1)/A
        UVDYP(M,N)=-A*EPSI(M2,N+1)/(EL1+1.)
   15 CONTINUE

c  generate associated Legendre polynomials

      SQRHLF=SQRT(.5)
      DO 20 M=2,MXP
        CONSQ(M) = SQRT(.5*(2.*EMM(M)+1.)/EMM(M))
   20 CONTINUE

      DO 30 J=1,IY

        Y = COA(J)
        X = SIA(J)

c  define (0,0) polynomial

        ALP(1,1) = SQRHLF

c  use recursive formula for N=1 (m=l) diagonal 

        DO 31 M=2,MXP
          ALP(M,1) = CONSQ(M)*Y*ALP(M-1,1)
   31   CONTINUE

c  derive n=2 polynomials from n=1

        DO 32 M=1,MXP
          ALP(M,2)=(X*ALP(M,1))*REPSI(M,2)
   32   CONTINUE

c  use recursive formula to obtain n=3 to n=jx

        DO 33 N=3,JX
        DO 33 M=1,MXP
          ALP(M,N)=(X*ALP(M,N-1)-EPSI(M,N-1)*ALP(M,N-2))*REPSI(M,N)
   33   CONTINUE

c  zero polynomials with absolute values smaller than 10e-30

        DO 34 N=1,JX
        DO 34 M=1,MXP
          IF(ABS(ALP(M,N)) .GE. 1.0E-30) THEN
            PALY(M,N)=ALP(M,N)
          ELSE
            PALY(M,N)=0.0
          ENDIF
   34   CONTINUE

c   pick off the required plynomials

        DO 35 N=1,JX
        DO 35 M=1,MX
          MM2=MM(M)+1
          POL(M,N,J)=PALY(MM2,N)
   35   CONTINUE

c  differentiate, adding weights

        WGHT=AM1*WT1(J)

        DO 40 M=1,MX
          MM2=MM(M)+1
          POLYDX(M,1,J)= EMM(MM2)*POL(M,1,J)*WGHT
          POLYDY(M,1,J)=-EMM(MM2)*POL(M,2,J)*WGHT*
     *     EPSI(MM2,2)
          POLYDX(M,JX,J)=0.0
          POLYDY(M,JX,J)=0.0
   40   CONTINUE

        DO 41 N=2,JTRUN1                                                     
        DO 41 M=1,MX
          MM2=MM(M)+1
          POLYDX(M,N,J)=EMM(MM2)*POL(M,N,J)*WGHT
          POLYDY(M,N,J)=
     *    ((ELL(MM2,N)+1.)*EPSI(MM2,N  )*POL(M,N-1,J)
     *    - ELL(MM2,N)    *EPSI(MM2,N+1)*POL(M,N+1,J))
     *    *WGHT
   41   CONTINUE
  
   30 CONTINUE

      RETURN
      END
****************************************************************
      SUBROUTINE LAP(STRM,VORM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      COMPLEX VORM(MX,JX),STRM(MX,JX)

c      DO 1 N=1,JX
c      DO 1 M=1,MX
c        VORM(M,N)=-STRM(M,N)*EL2(M,N)
c    1 CONTINUE  

      DO 1 M=1,MXJX
        VORM(M,1)=-STRM(M,1)*EL2(M,1)
    1 CONTINUE

      RETURN
      END
****************************************************************
      SUBROUTINE LAP2(STRM,VORM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      COMPLEX VORM(MX,JX),STRM(MX,JX)

C      DO 1 N=1,JX
C      DO 1 M=1,MX
C        VORM(M,N)=STRM(M,N)*EL4(M,N)
C    1 CONTINUE

      DO 1 M=1,MXJX
        VORM(M,1)=STRM(M,1)*EL4(M,1)
    1 CONTINUE

      RETURN
      END
*******************************************************************
      SUBROUTINE INVLAP(VORM,STRM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      COMPLEX VORM(MX,JX),STRM(MX,JX)

C      DO 1 N=1,JX
C      DO 1 M=1,MX
C        STRM(M,N)=-VORM(M,N)*ELM2(M,N)
C    1 CONTINUE

      DO 1 M=1,MXJX
        STRM(M,1)=-VORM(M,1)*ELM2(M,1)
    1 CONTINUE

      RETURN
      END
*********************************************************************
      SUBROUTINE GRAD(PSI,PSDX,PSDY)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION PSI(2,MX,JX),PSDX(2,MX,JX),PSDY(2,MX,JX)

      DO 1 N=1,JX
      DO 1 M=1,MX
        PSDX(2,M,N)=GRADX(M)*PSI(1,M,N)
        PSDX(1,M,N)=-GRADX(M)*PSI(2,M,N)
    1 CONTINUE

      DO 2 K=1,2
      DO 2 M=1,MX
        PSDY(K,M,1)=GRADYP(M,1)*PSI(K,M,2)
        PSDY(K,M,JX)=-GRADYM(M,JX)*PSI(K,M,JTRUN1)
    2 CONTINUE

      DO 3 K=1,2
      DO 3 N=2,JTRUN1
      DO 3 M=1,MX
        PSDY(K,M,N)=-GRADYM(M,N)*PSI(K,M,N-1)+GRADYP(M,N)*PSI(K,M,N+1)
    3 CONTINUE

      RETURN
      END
******************************************************************
      SUBROUTINE UVSPEC(VORM,DIVM,UCOSM,VCOSM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'
                                                        
      DIMENSION VORM(2,MX,JX),DIVM(2,MX,JX),UCOSM(2,MX,JX),
     * VCOSM(2,MX,JX),ZC(2,MX,JX),ZP(2,MX,JX)

C      DO 1 N=1,JX
C      DO 1 M=1,MX
C        ZP(2,M,N)=UVDX(M,N)*VORM(1,M,N)
C        ZP(1,M,N)=-UVDX(M,N)*VORM(2,M,N)
C        ZC(2,M,N)=UVDX(M,N)*DIVM(1,M,N)
C        ZC(1,M,N)=-UVDX(M,N)*DIVM(2,M,N)
C    1 CONTINUE

      DO 1 M=1,MXJX
        ZP(2,M,1)=UVDX(M,1)*VORM(1,M,1)
        ZP(1,M,1)=-UVDX(M,1)*VORM(2,M,1)
        ZC(2,M,1)=UVDX(M,1)*DIVM(1,M,1)
        ZC(1,M,1)=-UVDX(M,1)*DIVM(2,M,1)
    1 CONTINUE

      DO 3 K=1,2
      DO 3 M=1,MX
        UCOSM(K,M,1)=ZC(K,M,1)-UVDYP(M,1)*VORM(K,M,2)
        UCOSM(K,M,JX)=UVDYM(M,JX)*VORM(K,M,JTRUN1)
        VCOSM(K,M,1)=ZP(K,M,1)+UVDYP(M,1)*DIVM(K,M,2)
        VCOSM(K,M,JX)=-UVDYM(M,JX)*DIVM(K,M,JTRUN1)
    3 CONTINUE

      DO 4 K=1,2
      DO 4 N=2,JTRUN1
      DO 4 M=1,MX
        VCOSM(K,M,N)=-UVDYM(M,N)*DIVM(K,M,N-1)+UVDYP(M,N)*
     *   DIVM(K,M,N+1)+ZP(K,M,N)  
        UCOSM(K,M,N)= UVDYM(M,N)*VORM(K,M,N-1)-UVDYP(M,N)*
     *   VORM(K,M,N+1)+ZC(K,M,N)
    4 CONTINUE

      RETURN     
      END
*******************************************************************
      SUBROUTINE GRID(VORM,VORG,KCOS)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION VORG(IX,IL)
      COMPLEX VORM(MX,JX),VARM(MX,IL)

      CALL GRIDY(VORM,VARM)
      CALL GRIDX(VARM,VORG,KCOS)

      RETURN
      END
*********************************************************************
      SUBROUTINE SPEC(VORG,VORM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION VORG(IX,IL)
      COMPLEX VORM(MX,JX),VARM(MX,IL)

      CALL SPECX(VORG,VARM)
      CALL SPECY(VARM,VORM)

      RETURN
      END
*********************************************************************
      SUBROUTINE VDSPEC(UG,VG,VORM,DIVM,K)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION UG(IX,IL),VG(IX,IL)
      COMPLEX VORM(MX,JX),DIVM(MX,JX),UM(MX,IL),VM(MX,IL)
      DIMENSION UG1(IX,IL),VG1(IX,IL)

      IF (K.EQ.2) THEN
        DO 7 J=1,IL
          DO 8 I=1,IX
            UG1(I,J)=UG(I,J)*COSG(J)
            VG1(I,J)=VG(I,J)*COSG(J)
    8     CONTINUE
    7   CONTINUE
        CALL SPECX(UG1,UM)  
        CALL SPECX(VG1,VM)
      ELSE
        CALL SPECX(UG,UM)
        CALL SPECX(VG,VM)
      ENDIF

      CALL DSPY(UM,VM,VORM,2)
      CALL DSPY(UM,VM,DIVM,1)

      RETURN
      END
*********************************************************************
      SUBROUTINE DSPEC(UG,VG,DIVM,IVD,K)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION UG(IX,IL),VG(IX,IL)
      COMPLEX DIVM(MX,JX),UM(MX,IL),VM(MX,IL)
      DIMENSION UG1(IX,IL),VG1(IX,IL)

*  IVD=1 returns divergence, IVD=2 returns vorticity

      IF (K.EQ.2) THEN
        DO 7 J=1,IL
          DO 8 I=1,IX
            UG1(I,J)=UG(I,J)*COSG(J)
            VG1(I,J)=VG(I,J)*COSG(J)
    8     CONTINUE
    7   CONTINUE
        CALL SPECX(UG1,UM)  
        CALL SPECX(VG1,VM)
      ELSE
        CALL SPECX(UG,UM)
        CALL SPECX(VG,VM)
      ENDIF

      CALL DSPY(UM,VM,DIVM,IVD)

      RETURN
      END
*********************************************************************
      SUBROUTINE GRIDX1(VARM,VORG,KCOS)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION VORG(IX,IL)
      COMPLEX VARM(MX,IL)
      COMPLEX FMVAR(NWAVES,2)
      COMPLEX ZERO

      ZERO=(0.,0.)

      DO 1 M=1,NWAVES
        FMVAR(M,1)=ZERO
        FMVAR(M,2)=ZERO
    1 CONTINUE

      DO 10 J=1,IY

        J1=IL+1-J

        DO 20 M=1,MX
          FMVAR(M,1)=VARM(M,J)
          FMVAR(M,2)=VARM(M,J1)
   20   CONTINUE

        CALL NCWG2(FMVAR(1,1),FMVAR(1,2),vorg(1,j),vorg(1,j1))

        IF(KCOS.EQ.2) THEN
          DO 22 I=1,IX
            VORG(I,J)=Vorg(I,j)*COAM(J)
            VORG(I,J1)=Vorg(I,j1)*COAM(J)
   22     CONTINUE
        ENDIF

   10 CONTINUE

      RETURN
      END
*******************************************************************
      SUBROUTINE GRIDX(VARM,VORG,KCOS)

      INCLUDE 'param.h'
      PARAMETER (NWORK=IL*2*IX,IMAX=IX+3,MX2=MX*2)
      COMMON/FFTCOM/ TRIGS(IX*2),IFAX(19)
      INCLUDE 'spec.h'

      DIMENSION VORG(IX,IL),VARM(MX2,IL),FMVAR(IMAX,IL)
      DIMENSION WORK1(NWORK)

      DO 1 J=1,IL
      DO 1 M=1,IMAX
        FMVAR(M,J)=0.
    1 CONTINUE

      DO 2 J=1,IL
      DO 2 M=1,MX2
        FMVAR(M,J)=VARM(M,J)
    2 CONTINUE

      CALL RFFTMLT(FMVAR,WORK1,TRIGS,IFAX,1,IMAX,IX,IL,+1)

      IF(KCOS.EQ.1) THEN
        DO 4 J=1,IL
        DO 4 I=1,IX
          VORG(I,J)=FMVAR(I,J)
    4   CONTINUE
      ELSE
        DO 3 J=1,IL
        DO 3 I=1,IX
          VORG(I,J)=FMVAR(I,J)*COSGR(J)
    3   CONTINUE
      ENDIF

      RETURN
      END
*******************************************************************
      SUBROUTINE GRIDY(VORM,VARM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      COMPLEX VARM(MX,IL),VORM(MX,JX),POLVM(MX,JX)
      COMPLEX FMVAR1(MX,2),VORM1(MX,JX),ZERO
      
      ZERO=(0.,0.)

      DO 10 J=1,IY

        J1=IL+1-J

C        DO 20 N=1,JX
C        DO 20 M=1,MX
C          POLVM(M,N)=VORM(M,N)*POL(M,N,J)
C  20    CONTINUE

      DO 20 I=1,MXJX
          POLVM(I,1)=VORM(I,1)*POL(I,1,J)
   20 CONTINUE

        DO 21 M=1,MX
          FMVAR1(M,1)=ZERO
          FMVAR1(M,2)=ZERO
   21   CONTINUE

        DO 22 N=1,JX,2
          DO 23 M=1,MX
            FMVAR1(M,1)=FMVAR1(M,1)+POLVM(M,N)
   23     CONTINUE
   22   CONTINUE

        DO 24 N=2,JX,2
          DO 25 M=1,MX
            FMVAR1(M,2)=FMVAR1(M,2)+POLVM(M,N)
   25     CONTINUE
   24   CONTINUE

        DO 26 M=1,MX
          VARM(M,J1)=FMVAR1(M,1)+FMVAR1(M,2)
          VARM(M,J)=FMVAR1(M,1)-FMVAR1(M,2)
   26   CONTINUE

   10 CONTINUE

      RETURN
      END
*********************************************************************
      SUBROUTINE SPECX1(VORG,VARM)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION VORG(IX,IL)
      COMPLEX VARM(MX,IL),VM(NWAVES,IL)

      DO 10 J=1,IY
        J1=IL+1-J
        CALL NCGW2(VORG(1,J),VORG(1,J1),VM(1,J),VM(1,J1))
   10 CONTINUE

      DO 11 J=1,IL 
      DO 11 M=1,MX
        VARM(M,J)=VM(M,J)
   11 CONTINUE

      RETURN
      END
**********************************************************************
      SUBROUTINE SPECX(VORG,VARM)

      INCLUDE 'param.h'
      PARAMETER (NWORK=IL*2*IX,IMAX=IX+3,MX2=MX*2)
      COMMON/FFTCOM/ TRIGS(IX*2),IFAX(19)
      INCLUDE 'spec.h'

      DIMENSION VORG(IX,IL),vorg1(imax,il)
      DIMENSION WORK1(NWORK),VARM(MX2,IL)

      DO 1 J=1,IL
      DO 1 I=1,IX
        VORG1(I,J)=VORG(I,J)
    1 CONTINUE

      CALL RFFTMLT(VORG1,WORK1,TRIGS,IFAX,1,IMAX,IX,IL,-1)

      DO 11 J=1,IL 
      DO 11 M=1,MX2
        VARM(M,J)=VORG1(M,J)
   11 CONTINUE

      RETURN
      END
**********************************************************************
      SUBROUTINE SPECY(VARM,VORM)

      INCLUDE 'param.h'
      PARAMETER (MXJT1=MX*JTRUN1)
      INCLUDE 'spec.h'

      COMPLEX VARM(MX,IL),VORM(MX,JX),SCIFM(MX),DCIFM(MX),
     * DUM(MX,JTRUN1),ZERO

      ZERO=(0.,0.)

      DO 1 N=1,JX
      DO 1 M=1,MX
        VORM(M,N)=ZERO
    1 CONTINUE

      DO 10 J=1,IY

        J1=IL+1-J

        DO 11 M=1,MX
          SCIFM(M)=(VARM(M,J1)+VARM(M,J))*WT(J)
          DCIFM(M)=(VARM(M,J1)-VARM(M,J))*WT(J)
   11   CONTINUE

        DO 22 N = 1,JTRUN1,2
        DO 22 M = 1,MX
   22   DUM(M,N) = POL(M,N,J)*SCIFM(M)

        DO 23 N=2,JTRUN1,2
        DO 23 M=1,MX
   23   DUM(M,N) = POL(M,N,J)*DCIFM(M)

C        DO 12 N = 1,JTRUN1
C        DO 12 M = 1,MX
C          VORM(M,N) = VORM(M,N)+ DUM(M,N)
C   12   CONTINUE

      DO 12 M=1,MXJT1
        VORM(M,1)=VORM(M,1)+DUM(M,1)
   12 CONTINUE

   10 CONTINUE

      RETURN
      END
******************************************************************
      SUBROUTINE DSPY(UM,VM,DIVM,IVD)

      INCLUDE 'param.h'
      PARAMETER (MXJT1=MX*JTRUN1)
      INCLUDE 'spec.h'

      COMPLEX UM(MX,IL),VM(MX,IL),DIVM(MX,JX),
     * ZUM(MX,IL)
      COMPLEX DUM(MX,JX),SZAIFM(MX),
     * DZAIFM(MX),SBIFM(MX),DBIFM(MX),ZERO,II

      II=(0.,1.)
      ZERO=(0.,0.)

      DO 6 N=1,JX
      DO 6 M=1,MX
        DIVM(M,N)=ZERO
    6 CONTINUE

      IF(IVD.EQ.1) THEN
        DO 2 J=1,IL
        DO 2 M=1,MX
          ZUM(M,J)=II*UM(M,J)
    2   CONTINUE
      ELSE
        DO 52 J=1,IL
        DO 52 M=1,MX
          ZUM(M,J)=II*VM(M,J)
   52   CONTINUE
      ENDIF

      DO 1 J=1,IY

        J1=IL+1-J

        IF(IVD.EQ.1) THEN
          DO 5 M=1,MX
            SZAIFM(M)=ZUM(M,J1)+ZUM(M,J)
            DZAIFM(M)=ZUM(M,J1)-ZUM(M,J)
            SBIFM(M)=VM(M,J1)+VM(M,J)
            DBIFM(M)=VM(M,J1)-VM(M,J)
    5     CONTINUE
        ELSE
          DO 55 M=1,MX
            SZAIFM(M)=ZUM(M,J1)+ZUM(M,J)
            DZAIFM(M)=ZUM(M,J1)-ZUM(M,J)
            SBIFM(M)=-UM(M,J1)-UM(M,J)
            DBIFM(M)=-UM(M,J1)+UM(M,J)
   55     CONTINUE
        ENDIF

        DO 22 N = 1,JTRUN1,2
        DO 22 M = 1,MX
          DUM(M,N) = POLYDX(M,N,J)*SZAIFM(M)
     *     - POLYDY(M,N,J)*DBIFM(M)
   22   CONTINUE

        DO 23 N=2,JTRUN1,2
        DO 23 M=1,MX
          DUM(M,N) = POLYDX(M,N,J)*DZAIFM(M)
     *    - POLYDY(M,N,J)*SBIFM(M)
   23 CONTINUE

c        DO 3 N = 1,JTRUN1
c        DO 3 M = 1,MX
c          DIVM(M,N) = DIVM(M,N) + DUM(M,N)
c    3   CONTINUE

        DO 3 M=1,MXJT1
          DIVM(M,1)= DIVM(M,1) + DUM(M,1)
    3   CONTINUE

    1 CONTINUE

      RETURN
      END
****************************************************************
      SUBROUTINE SUB(A,B,C,K)

      DIMENSION A(K),B(K),C(K)
      COMPLEX A

      RADI=180./3.1415926

      DO 1 J=1,K
        C(J)=0.
        B(J)=CABS(A(J))
        AA=AIMAG(A(J))
        BB=REAL(A(J))
        IF(AA.EQ.0..AND.BB.EQ.0.) GO TO 1
        C(J)=ATAN2(AA,BB)*RADI
    1 CONTINUE

      RETURN
      END
*********************************************************************
      SUBROUTINE GLOBAL(VOR,VOR1,AV)

      INCLUDE 'param.h'
      INCLUDE 'spec.h'

      DIMENSION VOR(IX,IL),VOR1(IX,IL)

      AV=0.

      DO 1 I=1,IX
      DO 1 J=1,IY
        JJ=IL+1-J
        AV=AV+(VOR(I,J)+VOR(I,JJ))*WT(J)
    1 CONTINUE

      AV=AV/(2.*FLOAT(IX))

      DO 2 I=1,IX
      DO 2 J=1,IL
        VOR1(I,J)=VOR(I,J)-AV
    2 CONTINUE

      RETURN
      END
**************************************************************
      SUBROUTINE ZON(GG,GZ,GE)

      INCLUDE 'param.h'

      DIMENSION GG(IX,IL),GE(IX,IL),GZ(IL)

      DO 1 J=1,IL
        GZ(J)=0.
        DO 2 I=1,IX
          GZ(J)=GZ(J)+GG(I,J)
    2   CONTINUE
        GZ(J)=GZ(J)/FLOAT(IX)
        DO 3 I=1,IX
          GE(I,J)=GG(I,J)-GZ(J)
    3   CONTINUE
    1 CONTINUE

      RETURN
      END
**********************************************************************
