      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
