      SUBROUTINE ELLSOL(IMA, IM, JM, AE, INIT1, INIT2, DIFF
     *,                 C1, C2, C3, C4, C5, C6, FI, FO, WSV, WRK, IX)
C
C*************************************************************************
C                                                                        *
C     (IMA, JM) is the dimension of the forcing FI (input) and the       *
C     solution FO (output).  IM is less than or equal to IMA and         *
C     represents the physical domain (in longitude) over which the       *
C     solution is computed. AE is the radius of the sphere.              *
C                                                                        *
C     INIT1, INIT2 and DIFF are logical variables.                       *
C     When INIT1 is true, the value passed for INIT2 is ignored          *
C                                                                        *
C     C1, C2,...C6 are arrays of dimension JM containing the coefficients*
C     of the elliptic equation.                                          * 
C                                                                        *
C     WSV, WRK and IX are other miscellaneous work arrays.               *
C                                                                        *
C     The dimension of WSV should be at least [4*IM + 9*JM + 7].         *
C     This array stores constants needed in subsequent calls to ELLSOL.  *
c     It should not be overwritten unless INIT1 is true for the          *
c     next call to ELLSOL.                                               *
C                                                                        *
C     The dimension of WRK should be at least [14 * IM*JM].              *
C                                                                        *
C     The dimension of IX should be at least IM and it should not be     *
C     overwritten unless INIT1 is true for the next call to ELLSOL.      *
C                                                                        *
C                                                                        *
C*************************************************************************
C
      DIMENSION    FI(IMA,JM), FO(IMA,JM)
     *,            C1(JM),     C2(JM), C3(JM), C4(JM), C5(JM), C6(JM)
     *,            WSV(1),     WRK(1), IX(IM)
C
      LOGICAL   INIT1, INIT2, DIFF
C
      ITRG = 3*(IM+2) + 1
      IMB2 = IM / 2
      IJM  = (IMB2-1) * JM
      IJMM = IJM*22 + JM
      IJM2 = IJMM + (IM+2)*(JM-2)
C
      CALL DIRSOL(IMA, IM, JM, IMB2, ITRG, AE, INIT1, INIT2, DIFF
     *,           C1, C2, C3, C4, C5, C6, FI, FO, IX 
C
     *,           WSV(1),      WSV(JM+1),   WSV(JM*2+1), WSV(JM*3+1)
     *,           WSV(JM*4+1), WSV(JM*5+1), WSV(JM*6+1), WSV(JM*7+1)
     *,           WSV(JM*9+1), WSV(JM*9+IMB2+1), WSV(JM*9+IM+1) 
C
     *,     WRK(1),       WRK(IJM+1),    WRK(IJM*2+1),  WRK(IJM*3+1)
     *,     WRK(IJM*4+1), WRK(IJM*5+1),  WRK(IJM*6+1),  WRK(IJM*7+1)
     *,     WRK(IJM*8+1), WRK(IJM*9+1),  WRK(IJM*10+1), WRK(IJM*11+1)
     *,     WRK(IJM*12+1),WRK(IJM*13+1), WRK(IJM*14+1), WRK(IJM*15+1)
     *,     WRK(IJM*16+1),WRK(IJM*17+1), WRK(IJM*18+1), WRK(IJM*19+1)
     *,     WRK(IJM*20+1),WRK(IJM*21+1), WRK(IJM*22+1), WRK(IJMM+1)
     *,     WRK(IJM2+1),  WRK(1),        WRK(JM+1),     WRK(JM*2+1)
     *,     WRK(JM*3+1),  WRK(JM*4+1))
C
      RETURN
      END
c
c
      SUBROUTINE DIRSOL(IMA, IM, JM, IMB2, ITRG, AE, INIT1, INIT2, DIFF
     *,                 C1, C2, C3, C4, C5, C6, FI, FO, IX 
C
     *,                 CDS1, CDN1, CDS2, CDN2, CDC1, CDC2, CDC3, CPH
     *,                 WV,   WVSQ, TR
 
     *,                 AN1, BN1, CN1, DN1, FN1, PN1, QN1, RN1
     *,                 AN2, BN2, CN2, DN2, FN2, PN2, QN2, RN2
     *,                 AL1, AL2, GM1, GM2, BT1, BT2
     *,                 FM, A, B, AN,  BN,  CN,  DN, FN)
C
C     ITRG SHOULD AT LEAST BE 3*(IM+2) + 1
C
      INTEGER   FORWARD, BACKWARD
C
      PARAMETER (FORWARD=-1,BACKWARD=1)
      PARAMETER(PI=3.1415926535898, TWOPI=PI+PI, PIO2=0.5*PI)
C
      DIMENSION FI(IMA,JM), FO(IMA,JM)
     *,         C1(JM),     C2(JM), C3(JM), C4(JM), C5(JM), C6(JM)
     *,         TR(ITRG),   IX(IM), WVSQ(IMB2), WV(IMB2)
C
      DIMENSION AN1(IMB2-1,JM), BN1(IMB2-1,JM), CN1(IMB2-1,JM)
     *,         AN2(IMB2-1,JM), BN2(IMB2-1,JM), CN2(IMB2-1,JM)
     *,         QN1(IMB2-1,JM), PN1(IMB2-1,JM), RN1(IMB2-1,JM)
     *,         QN2(IMB2-1,JM), PN2(IMB2-1,JM), RN2(IMB2-1,JM)
     *,         DN1(IMB2-1,JM), DN2(IMB2-1,JM)
     *,         FN1(IMB2-1,JM), FN2(IMB2-1,JM)
C
     *,         AL1(IMB2-1,JM), AL2(IMB2-1,JM)
     *,         GM1(IMB2-1,JM), GM2(IMB2-1,JM)
     *,         BT1(IMB2-1,JM), BT2(IMB2-1,JM)
C
      DIMENSION AN(JM),   BN(JM), CN(JM), DN(JM), FN(JM)
     *,         CDN1(JM), CDS1(JM), CDN2(JM), CDS2(JM)
     *,         CDC1(JM), CDC2(JM), CDC3(JM), FM(JM),   CPH(JM,2)
C
      DIMENSION A(IM+2,JM-2), B(IM*2,JM-2)
C
      LOGICAL   INIT1, INIT2, DIFF, FIRST, SETC0
C
      DATA JMM1, JMM2, IMP1, IMP2, LEN/5*0/, AREA/0.0/
      DATA FIM, DLM, DPH, RDLM, RDPH, DY, AESQ, CDNP, CDSP/9*0.0/
      DATA FIRST/.FALSE./
C
      SAVE 
C
C   INITIALIZATION ON FIRST CALL
C
      IF (INIT1) THEN
         CALL FFTFAX(IM, IX, TR)
C
         JMM1 = JM - 1
         JMM2 = JM - 2
         IMP1 = IM + 1
         IMP2 = IM + 2
         LEN  = IMB2 - 1
         FIM  = 1.0 / FLOAT(IM)
         DLM  = TWOPI * FIM
         DPH  = PI / FLOAT(JMM1)
C
         RDLM  = 1.0 / DLM
         RDPH  = 1.0 / DPH
C
         DY    = AE * DPH
         AESQ  = AE * AE
C
         DO 10 J=1,JM
         TEM      = -PIO2  + (J-1)*DPH
         CPH(J,1) = COS(TEM + 0.5*DPH)
         CPH(J,2) = COS(TEM)
   10    CONTINUE
C
         CPH(1,2)  = 0.0
         CPH(JM,2) = 0.0
         FIRST     = .TRUE.
         INIT2     = .TRUE.
      ENDIF
C
      IF (INIT2) THEN
         DO 20 J=2,JMM1
         TEM1    = 1.0  / (CPH(J,2) * 2.0 * DY)
         TEM2    = 1.0  / (CPH(J,2) * DY*DY)
C
         CDN1(J) = C5(J)   * CPH(J,1)   * TEM1
         CDS1(J) = C5(J-1) * CPH(J-1,1) * TEM1
         CDN2(J) = C3(J)   * CPH(J,1)   * TEM2
         CDS2(J) = C3(J-1) * CPH(J-1,1) * TEM2
C
         CDC1(J) = C1(J) / (CPH(J,2) * CPH(J,2) * AESQ)
         CDC2(J) = C4(J) / (AE * CPH(J,2))
         CDC3(J) = C2(J) / (AE * CPH(J,2) * 2.0 * DY)
C
   20    CONTINUE
         CDC1(1)  = 4.0 * C3(1)    / (DY*DY)
         CDC1(JM) = 4.0 * C3(JMM1) / (DY*DY)
         CDNP     = C5(JMM1) * 2.0 / DY
         CDSP     = C5(1   ) * 2.0 / DY
C
         IF (DIFF) THEN
            DO 30 I=1,IMB2
            WVSQ(I) = (SIN(0.5*I*DLM) * (RDLM * 2.0)) ** 2
            WV  (I) = SIN(I*DLM) * RDLM
   30       CONTINUE
         ELSE
            DO 40 I=1,IMB2
            WVSQ(I) = I * I
            WV  (I) = I
   40       CONTINUE
         ENDIF
      ENDIF
C
      DO 25 J=1,JM
      IF (C6(J) .NE. 0.0) GO TO 26
   25 CONTINUE
      SETC0 = .TRUE.
      CALL REMGLM(IMA, IM, JM, FI, CPH, FIRST, AREA)
      GO TO 27
   26 SETC0 = .FALSE.
   27 CONTINUE
C
      DO 42 J=2,JMM1
      FM(J) = 0.0
      DO 42 I=1,IM
      FM(J) = FM(J) + FI(I,J)
  42  CONTINUE
      DO 43 J=2,JMM1
      FM(J) = FM(J) * FIM
  43  CONTINUE
      DO 45 J=2,JMM1
      DO 45 I=1,IM
      A(I,J-1) = FI(I,J) - FM(J)
C     A(I,J-1) = FI(I,J)
  45  CONTINUE
C
      DO 50 J=1,JMM2
      A(IM+1,J) = A(1,J)
      A(IM+2,J) = A(2,J)
  50  CONTINUE
C
c     CALL RFFTMLT(A,B,TR,IX,1,IMP2,IM,JMM2,FORWARD)
      CALL FFT99(A,B,TR,IX,1,IMP2,IM,JMM2,FORWARD)
C
      DO 60 J=2,JMM1
      DO 60 I=1,LEN
      AN1(I,J) =   CDS2(J) - CDS1(J)
      BN1(I,J) = -(CDS2(J) + CDN2(J) + WVSQ(I) * CDC1(J)
     *           + CDS1(J) - CDN1(J) - C6(J))
      CN1(I,J) =   CDN2(J) + CDN1(J)
      AN2(I,J) =   CDC3(J) * WV(I)
      BN2(I,J) = - CDC2(J) * WV(I)
      CN2(I,J) = - CDC3(J) * WV(I)
C
      DN1(I,J) =   A(I+I+1,J-1)
C
      QN1(I,J) =   CN2(I,J)
      PN1(I,J) = - BN2(I,J)
      RN1(I,J) =   AN2(I,J)
      QN2(I,J) =   AN1(I,J)
      PN2(I,J) =   BN1(I,J)
      RN2(I,J) =   CN1(I,J)
C
      DN2(I,J) = A(I+I+2,J-1)
   60 CONTINUE
C
      CALL SO2ODE(LEN, JM,  AN1, BN1, CN1, DN1, AN2, BN2, CN2, DN2
     *,                 QN1, PN1, RN1, QN2, PN2, RN2, FN1, FN2
     *,                 AL1, AL2, GM1, GM2, BT1, BT2)
C
      DO 70 J=2,JMM1
      DO 70 I=1,LEN
      A(I+I+1,J-1) = FN1(I,J)
      A(I+I+2,J-1) = FN2(I,J)
   70 CONTINUE
C
      DO 80 J=2,JMM1
      AN(J)  = CDS2(J) - CDS1(J)
      BN(J)  = -(CDS2(J) + CDN2(J) + CDS1(J) - CDN1(J)
     *                             + WVSQ(IMB2) * CDC1(J) - C6(J))
      CN(J)  = CDN2(J) + CDN1(J)
      DN(J)  = A(IMP1,J-1)
   80 CONTINUE
C
      CALL SO1ODE(JM,  AN, BN, CN, DN, FN, AL1, BT1)
C
      DO 90 J=2,JMM1
      A(IMP1,J-1) = FN(J)
C     BN(J)       = -(AN(J) + CN(J) - C6(J))
      BN(J)  = -(CDS2(J) + CDN2(J) + CDS1(J) - CDN1(J) - C6(J))
C###  DN(J)       = A(1,J-1)
      DN(J)       = FM(J)
   90 CONTINUE
      BN(1)  = - (-C6(1)  + CDC1(1)  - CDSP)
      BN(JM) = - (-C6(JM) + CDC1(JM) + CDNP)
      AN(JM) = CDC1(JM)  - CDNP
      CN(1)  = CDC1(1)   + CDSP
      DN(1)  = FI(1,1)
      DN(JM) = FI(1,JM)
C
      CALL SO1OD2(JM,  AN, BN, CN, DN, FN, AL1, BT1)
C
      IF (SETC0) THEN
         FLDM = (FN(1) + FN(JM)) * (AREA*CPH(1,1)*0.25)
         DO 92 J=2,JMM1
         FLDM = FLDM + FN(J) * (CPH(J,2)*AREA)
   92    CONTINUE
         DO 94 J=1,JM
         FN(J) = FN(J) - FLDM
   94    CONTINUE
C
      ENDIF
C
      DO 100 J=2,JMM1
C     A(1,   J-1) = FN(J)
      A(2,   J-1) = 0.0
      A(IMP2,J-1) = 0.0
  100 CONTINUE
C
c     CALL RFFTMLT(A,B,TR,IX,1,IMP2,IM,JMM2,BACKWARD)
      CALL FFT99(A,B,TR,IX,1,IMP2,IM,JMM2,BACKWARD)
C
      DO 120 J=1,JMM2
      DO 110 I=1,IM
C     FO(I,J+1) = A(I,J)
      FO(I,J+1) = A(I,J) + FN(J+1)
 110  CONTINUE
 120  CONTINUE
C
      DO 130 I=1,IM
      FO(I,1)  = FN(1)
      FO(I,JM) = FN(JM)
  130 CONTINUE
C
      RETURN
      END
c
      SUBROUTINE SO2ODE(LEN, JM,  AN1, BN1, CN1, DN1, AN2, BN2, CN2, DN2
     *,                 QN1, PN1, RN1, QN2, PN2, RN2, FN1, FN2
     *,                 AL1, AL2, GM1, GM2, BT1, BT2)
C
      DIMENSION AN1(LEN,JM), BN1(LEN,JM), CN1(LEN,JM)
     *,         AN2(LEN,JM), BN2(LEN,JM), CN2(LEN,JM)
     *,         QN1(LEN,JM), PN1(LEN,JM), RN1(LEN,JM)
     *,         QN2(LEN,JM), PN2(LEN,JM), RN2(LEN,JM)
     *,         DN1(LEN,JM), DN2(LEN,JM), FN1(LEN,JM), FN2(LEN,JM)
C
      DIMENSION AL1(LEN,JM), AL2(LEN,JM), GM1(LEN,JM), GM2(LEN,JM)
     *,         BT1(LEN,JM), BT2(LEN,JM)
C
      JMM1 = JM - 1
C
      DO 10 I=1,LEN
      AL1(I,1) = 0.0
      AL2(I,1) = 0.0
      BT1(I,1) = 0.0
      BT2(I,1) = 0.0
      GM1(I,1) = 0.0
      GM2(I,1) = 0.0
   10 CONTINUE
C
      DO 30 J=2,JMM1
      JM1 = J - 1
C
      DO 20 I=1,LEN
      SA1 = AN1(I,J) * AL1(I,JM1) + AN2(I,J) * GM1(I,JM1) + BN1(I,J)
      SA2 = AN1(I,J) * AL2(I,JM1) + AN2(I,J) * GM2(I,JM1) + BN2(I,J)
C
      SB1 = QN1(I,J) * AL1(I,JM1) + QN2(I,J) * GM1(I,JM1) + PN1(I,J)
      SB2 = QN1(I,J) * AL2(I,JM1) + QN2(I,J) * GM2(I,JM1) + PN2(I,J)
C
      SD1 = DN1(I,J) - AN1(I,J) * BT1(I,JM1) - AN2(I,J) * BT2(I,JM1)
      SD2 = DN2(I,J) - QN1(I,J) * BT1(I,JM1) - QN2(I,J) * BT2(I,JM1)
C
      RM  = 1.0 / (SA1*SB2 - SB1*SA2)
C
      AL1(I,J) = - RM * (CN1(I,J) * SB2 - RN1(I,J) * SA2)
      AL2(I,J) = - RM * (CN2(I,J) * SB2 - RN2(I,J) * SA2)
      BT1(I,J) =   RM * (SD1      * SB2 - SD2      * SA2)
C
      GM1(I,J) =   RM * (CN1(I,J) * SB1 - RN1(I,J) * SA1)
      GM2(I,J) =   RM * (CN2(I,J) * SB1 - RN2(I,J) * SA1)
      BT2(I,J) = - RM * (SD1      * SB1 - SD2      * SA1)
C
   20 CONTINUE
   30 CONTINUE
C
      DO 40 I=1,LEN
      FN1(I,JM) = 0.0
      FN2(I,JM) = 0.0
   40 CONTINUE
C
      DO 60 J=JMM1,1,-1
      DO 50 I=1,LEN
      FN1(I,J) = AL1(I,J)*FN1(I,J+1) + AL2(I,J)*FN2(I,J+1) + BT1(I,J)
      FN2(I,J) = GM1(I,J)*FN1(I,J+1) + GM2(I,J)*FN2(I,J+1) + BT2(I,J)
   50 CONTINUE
   60 CONTINUE
C
      RETURN
      END
c
      SUBROUTINE SO1ODE(JM,  AN, BN, CN, DN, FN, AL, BT)
C
      DIMENSION AN(JM), BN(JM), CN(JM), DN(JM), FN(JM)
C
      DIMENSION AL(JM), BT(JM)
C
      JMM1 = JM - 1
C
      AL(1)  = 0.0
      BT(1)  = 0.0
      FN(JM) = 0.0
C
      DO 10 J=2,JMM1
      JM1   = J - 1
      SA1   = 1.0 / (AN(J) * AL(JM1) + BN(J))
      AL(J) = - CN(J) * SA1
      BT(J) = (DN(J) - AN(J) * BT(JM1)) * SA1
   10 CONTINUE
C
      DO 20 J=JMM1,1,-1
      FN(J) = AL(J)*FN(J+1) + BT(J)
   20 CONTINUE
C
      RETURN
      END
c
      SUBROUTINE SO1OD2(JM,  AN, BN, CN, DN, FN, AL, BT)
C
      DIMENSION AN(JM), BN(JM), CN(JM), DN(JM), FN(JM)
C
      DIMENSION AL(JM), BT(JM)
C
      JMM1 = JM - 1
C
      TEM = 1.0 / BN(1)
      AL(1) = - CN(1) * TEM
      BT(1) =   DN(1) * TEM
C
      DO 10 J=2,JMM1
      JM1    = J - 1
      SA1    = 1.0 / (AN(J) * AL(JM1) + BN(J))
      AL(J)  = - CN(J) * SA1
      BT(J)  = (DN(J) - AN(J) * BT(JM1)) * SA1
   10 CONTINUE
C
      FN(JM) = (DN(JM) - AN(JM)*BT(JMM1)) / (BN(JM) + AN(JM)*AL(JMM1))
C
      DO 20 J=JMM1,1,-1
      FN(J) = AL(J)*FN(J+1) + BT(J)
   20 CONTINUE
C
      RETURN
      END
c
      SUBROUTINE REMGLM(IMA, IM, JM, FLD, CPH, FIRST, AREA)
C
      DIMENSION FLD(IMA,JM), CPH(JM,2)
      LOGICAL FIRST
      DATA JMM1/0/, FIM/0.0/
      SAVE
C
      IF (FIRST) THEN
         JMM1   = JM - 1
         FIM    = 1.0 / FLOAT(IM)
         AREA   = 0.0
         DO 10 J=2,JMM1
         AREA   = AREA + CPH(J,2)
   10    CONTINUE
         AREA   = 1.0 / (AREA + CPH(1,1) * 0.5)
         FIRST  = .FALSE.
      ENDIF
C
      FLDM = (FLD(1,1) + FLD(1,JM)) * (AREA*CPH(1,1)*0.25)
      DO 20 J=2,JMM1
      DO 20 I=1,IM
      FLDM = FLDM + FLD(I,J) * (CPH(J,2)*AREA*FIM)
   20 CONTINUE
C
      DO 30 I=1,IMA*JM
      FLD(I,1) = FLD(I,1) - FLDM
   30 CONTINUE
C
      RETURN
      END
