c===================================================================
      program STOPOC
c
c     Stokes's potential coefficients of the external gravitational
c     potential of a 3D density layer with non-spherical boundaries
c
c     written by Zdenek Martinec (June 22, 2017)
c===================================================================
      implicit real*8 (a,b,d-h,o-z), complex*16(c)
      parameter (ndimjm=65341)
      dimension rhoa(720,360),rhob(720,360)
      dimension topoa(720,360),topob(720,360),data1(720,360)
      dimension cpotjm(ndimjm),degpwrj(360)
c
      common /meanEarth/ rmean,rhomean
      common /mean_radii/rada_ref,radb_ref
c--------------------------------------------------------------------
      rmean=6371.d03
      rhomean=5.52d03
c--------------------------------------------------------------------
c...Mean depths of the layer boundaries (in km)
ccc      ha=80.d0
ccc      hb=56.d0
ccc      ha=0.d0
ccc      hb=0.d0
      ha=10.d03
      hb=5.d03
c...Mean radii of the layer boundaries
      rada_ref=rmean-ha
      radb_ref=rmean-hb
c--------------------------------------------------------------------
c   INPUTS:
c      jmax  ... cut-off degree of Stokes's potential coefficients
c      h_obs ... height of the observer above the Earth surface (in km)
c
      jmax=359
      h_obs=0.d0                          ! in km
c
      rad_obs=rmean+h_obs*1.d03
      write(6,*) 'jmax=',jmax,' h_obs=',h_obs,' rad_obs=',rad_obs
c--------------------------------------------------------------------
c...Outputs
      open(66,file='stopoc66.lis')         ! Checklist
      open(71,file='stopoc71.out')         ! For plotting of input densities
      open(72,file='stopoc72.out')         ! For plotting of input topographies
      open(73,file='stopoc73.out')         ! alfa and beta arrays
      open(74,file='stopoc74.out')         ! For plotting of the densities
      open(75,file='stopoc75.out')         ! Degree-power spectrum of potential coefficients
      open(76,file='stopoc76.out')         ! Gravitational potential at the
c----------------------------------
c     Densities at the boundaries
c----------------------------------
      write(6,*) 'Reading the densities at the boundaries:'
      open(10,file='rhoBotv05.txt')    ! Input
      open(11,file='rhoTopv05.txt')    ! Input
c
      nth=360
      nph=nth+nth
      step=0.5d0
c
      ithinv=nth+1
      i=0
      do 1 ith=1,nth
      ithinv=ithinv-1
      do 1 iph=1,nph
      read(10,*) qlon1,qlat1,val1
      read(11,*) qlon2,qlat2,val2
      i=i+1
c....
      val1=2.9d03
      val2=val1
c....
      rhoa(iph,ithinv)=val1
      rhob(iph,ithinv)=val2
    1 continue
      write(6,*) 'Number of data=',i
      close(10)
      close(11)
c...Check by plotting
      i=0
      do 2 ith=1,nth
      do 2 iph=1,nph
      i=i+1
      qlon=dfloat(iph-1)*step+step/2.d0
      qlat=-90.d0+dfloat(ith-1)*step+step/2.d0
      write(71,212) qlon,qlat,rhoa(iph,ith),rhob(iph,ith)
    2 continue
  212 format(2f7.2,2x,2f9.1)
      write(6,*) 'Number of data=',i
c----------------------------------
c     Topographies of the boundaries
c----------------------------------
      keytopo=1                    ! keytopo=0 ... spherical boundaries
      write(6,*) 'keytopo=',keytopo
c
c...spherical boundaries
c...
      if(keytopo.eq.0) then
c...
      do 6 ith=1,nth
      do 6 iph=1,nph
      topoa(iph,ith)=0.d0
    6 topob(iph,ith)=0.d0
c...non-spherical boundaries
c...
      else
c...
      write(6,*) 'Reading the topographies of boundaries:'
      open(12,file='botBv05.txt')      ! Input
      open(13,file='topBv05.txt')      ! Input
c
      ithinv=nth+1
      i=0
      do 3 ith=1,nth
      ithinv=ithinv-1
      do 3 iph=1,nph
      i=i+1
      read(12,*) qlon,qlat,pom
      topoa(iph,ithinv)=pom+ha
      read(13,*) qlon,qlat,pom
    3 topob(iph,ithinv)=pom+hb          
      write(6,*) 'Number of data=',i
      close(12)
      close(13)
c...Check by plotting
      i=0
      do 4 ith=1,nth
      do 4 iph=1,nph
      i=i+1
      qlon=dfloat(iph-1)*step+step/2.d0
      qlat=-90.d0+dfloat(ith-1)*step+step/2.d0
    4 write(72,213) qlon,qlat,topoa(iph,ith),topob(iph,ith)
  213 format(2f7.2,2x,2f11.2)
      write(6,*) 'Number of data=',i
c...
      endif
c...
c--------  Computations -------------------------------------
c...Stokes's coefficients
      call STOKES_COEFFS_3D(keytopo,nth,step,topoa,topob,rhoa,rhob,
     : jmax,cpotjm)
c--------  Outputs ------------------------------------------
c...Normalization of potential coeffs wrt the radius of the observer
      p=rmean/rad_obs
      write(6,*) 'Stokes coefficients:'
      write(6,*) 'p=rmean/rad_obs=',p
      do 5 j=0,jmax
      pj=p**(j+1)
      do 5 m=0,j
      jm=j*(j+1)/2+m+1
      cpotjm(jm)=cpotjm(jm)*pj
    5 write(66,201) j,m,cpotjm(jm)
  201 format(2i4,4e13.5)
c...Removing (0,0),(1,0),(1,1) and (2,0) coefficients
      cpotjm(1)=(0.d0,0.d0)
ccc      cpotjm(2)=(0.d0,0.d0)
ccc      cpotjm(3)=(0.d0,0.d0)
ccc      cpotjm(4)=(0.d0,0.d0)
c...Normalization for geoidal heights
      fac=rmean
c----------------------------------------
c     Spatial representation
c----------------------------------------
      jmax=240
c
      call HARMSY(nth,jmax,cpotjm,data1)
      do 7 ith=1,nth
      do 7 iph=1,nph
      qlon=(iph-1)*step+step/2.d0
      qlat=-90.+(ith-1)*step+step/2.d0
      write(76,202) qlon,qlat,data1(iph,ith)*fac
    7 continue
  202 format(2f7.2,f10.3)
c...Degree power in geoidal heights
      call DEGPWR(jmax,fac,cpotjm,degpwrj)
      do 8 j=0,jmax
    8 write(75,*) j,degpwrj(j+1)
      stop
      end
c================================================================
      subroutine STOKES_COEFFS_3D(keytopo,nth,step,topoa,topob,
     : rhoa,rhob,jmax,cpotjm)
c
c     Stokes's coefficients of the gravitational potential
c     of a 3D density layer with non-spherical boundaries
c
c     References
c
c (1) Martinec, Z., Pěč, K. and Burša, M., 1989. The Phobos gravitational 
c     field modeled on the basis of its topography. 
c     Earth, Moon, and Planets, 145, 219–235.
c
c (2) Fullea, J., Lebedev, S., Martinec, Z., Celli, N. L., 2021. 
c     WINTERC-G: mapping the upper mantle thermochemical heterogeneity 
c     from coupled geophysical-petrological inversion of seismic waveforms, 
c     heat flow, surface elevation and gravity satellite data. 
c     Geophys. J. Int., 226, 146–191.
c================================================================
      implicit none
      integer*4 jmax,nth,nph,ith,iph,ndimjm,j,m,jm,keytopo
      real*8 rmean,rhomean
      real*8 step,topoa,topob,rhoa,rhob,rada_ref,radb_ref
      real*8 alfa,beta,rada,radb,qlon,qlat,r,densi,data1
      real*8 qnorm,pom,pomj
      real*8 pa,paj,paj1,paj2,paj3,paj4
      real*8 pb,pbj,pbj1,pbj2,pbj3,pbj4
      real*8 dj,fac1,fac2,fac3,fak1,fak2,fak3
      dimension rhoa(720,360),rhob(720,360)
      dimension topoa(720,360),topob(720,360)
      dimension alfa(720,360),beta(720,360)
      dimension data1(720,360)
      parameter (ndimjm=65341)
      complex*16 cpotjm,cajm,cbjm,cpoma,cpomb,cxjm
      complex*16 cas1jm,cas2jm,cas3jm
      complex*16 cbs1jm,cbs2jm,cbs3jm
      complex*16 cat1jm,cat2jm,cat3jm
      complex*16 cbt1jm,cbt2jm,cbt3jm
      dimension cpotjm(*),cajm(ndimjm),cbjm(ndimjm),cxjm(ndimjm)
      dimension cas1jm(ndimjm),cas2jm(ndimjm),cas3jm(ndimjm)
      dimension cbs1jm(ndimjm),cbs2jm(ndimjm),cbs3jm(ndimjm)
      dimension cat1jm(ndimjm),cat2jm(ndimjm),cat3jm(ndimjm)
      dimension cbt1jm(ndimjm),cbt2jm(ndimjm),cbt3jm(ndimjm)
c
      common /meanEarth/ rmean,rhomean
      common /mean_radii/rada_ref,radb_ref
c--------------------------------------------------------------------
      nph=2*nth
      do 1 ith=1,nth
      do 1 iph=1,nph
      rada=rada_ref+topoa(iph,ith)
      radb=radb_ref+topob(iph,ith)
      alfa(iph,ith)=(rhob(iph,ith)-rhoa(iph,ith))/(radb-rada)
      beta(iph,ith)=rhoa(iph,ith)-alfa(iph,ith)*rada
c...Output for plotting
      qlon=dfloat(iph-1)*step+step/2.d0
      qlat=-90.d0+dfloat(ith-1)*step+step/2.d0
    1 write(73,201) qlon,qlat,alfa(iph,ith),beta(iph,ith)
  201 format(2f7.2,2x,2f15.7)
c...Check the density in  the middle of the layer
      do 2 ith=1,nth
      do 2 iph=1,nph
      r=(rada_ref+topoa(iph,ith)+radb_ref+topob(iph,ith))/2.d0
      densi=alfa(iph,ith)*r+beta(iph,ith)
      qlon=dfloat(iph-1)*step+step/2.d0
      qlat=-90.d0+dfloat(ith-1)*step+step/2.d0
    2 write(74,212) qlon,qlat,densi
  212 format(2f7.2,2x,2f9.1)
c--------------------------------
c...spherical-boundary terms
      call HARMAN(nth,jmax,alfa,cxjm)
      call HARMLS(nth,jmax,cxjm,cajm)
      call HARMAN(nth,jmax,beta,cxjm)
      call HARMLS(nth,jmax,cxjm,cbjm)
c--------------------------
c...Stokes's coefficients
c--------------------------
      pa=rada_ref/rmean
      pb=radb_ref/rmean
      write(6,*) 'pa=',pa,' pb=',pb
c
      do 50 j=0,jmax
      paj3=pa**(j+3)
      paj4=paj3*pa
      pbj3=pb**(j+3)
      pbj4=pbj3*pb
c
      dj=dfloat(j)
      fac1=dj+3.d0
      fak1=dj+4.d0
      qnorm=3.d0/rhomean/(dj+dj+1.d0)
      do 50 m=0,j
      jm=j*(j+1)/2+m+1
      cpoma=cajm(jm)*(pbj4-paj4)
      cpomb=cbjm(jm)*(pbj3-paj3)
   50 cpotjm(jm)=qnorm*(rmean*cpoma/fak1+cpomb/fac1)
c...............................
      if(keytopo.eq.0) return
c...............................
c------------------------------------------------------
c     SH analysis of 12 products for higher-order terms
c------------------------------------------------------
c-------------------------------------------
      write(6,*) 'Compute 1st order terms:'
c...1st order terms
      do 12 ith=1,nth
      do 12 iph=1,nph
      pom=topoa(iph,ith)/rmean
   12 data1(iph,ith)=alfa(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cas1jm)
c
      do 16 ith=1,nth
      do 16 iph=1,nph
      pom=topoa(iph,ith)/rmean
   16 data1(iph,ith)=beta(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cbs1jm)
c
      do 22 ith=1,nth
      do 22 iph=1,nph
      pom=topob(iph,ith)/rmean
   22 data1(iph,ith)=alfa(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cat1jm)
c
      do 26 ith=1,nth
      do 26 iph=1,nph
      pom=topob(iph,ith)/rmean
   26 data1(iph,ith)=beta(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cbt1jm)
c-------------------------------------------
      write(6,*) 'Compute 2nd order terms:'
c...2nd order terms
      do 13 ith=1,nth
      do 13 iph=1,nph
      pom=topoa(iph,ith)/rmean
      pom=pom**2
   13 data1(iph,ith)=alfa(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cas2jm)
c
      do 17 ith=1,nth
      do 17 iph=1,nph
      pom=topoa(iph,ith)/rmean
      pom=pom**2
   17 data1(iph,ith)=beta(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cbs2jm)
c
      do 23 ith=1,nth
      do 23 iph=1,nph
      pom=topob(iph,ith)/rmean
      pom=pom**2
   23 data1(iph,ith)=alfa(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cat2jm)
c
      do 27 ith=1,nth
      do 27 iph=1,nph
      pom=topob(iph,ith)/rmean
      pom=pom**2
   27 data1(iph,ith)=beta(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cbt2jm)
c-------------------------------------------
      write(6,*) 'Compute 3rd order terms:'
c...3rd order terms
      do 14 ith=1,nth
      do 14 iph=1,nph
      pom=topoa(iph,ith)/rmean
      pom=pom**3
   14 data1(iph,ith)=alfa(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cas3jm)
c
      do 18 ith=1,nth
      do 18 iph=1,nph
      pom=topoa(iph,ith)/rmean
      pom=pom**3
   18 data1(iph,ith)=beta(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cbs3jm)
c
      do 24 ith=1,nth
      do 24 iph=1,nph
      pom=topob(iph,ith)/rmean
      pom=pom**3
   24 data1(iph,ith)=alfa(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cat3jm)
c
      do 28 ith=1,nth
      do 28 iph=1,nph
      pom=topob(iph,ith)/rmean
      pom=pom**3
   28 data1(iph,ith)=beta(iph,ith)*pom
      call HARMAN(nth,jmax,data1,cxjm)
      call HARMLS(nth,jmax,cxjm,cbt3jm)
c-------------------------------------------
c--------------------------
c...Stokes's coefficients
c--------------------------
      do 51 j=0,jmax
      paj=pa**j
      paj1=paj*pa
      paj2=paj*pa*pa
      paj3=paj*pa*pa*pa
      pbj=pb**j
      pbj1=pbj*pb
      pbj2=pbj*pb*pb
      pbj3=pbj*pb*pb*pb
c
      dj=dfloat(j)
      fac1=dj+3.d0
      fac2=(dj+3.d0)*(dj+2.d0)/2.d0
      fac3=(dj+3.d0)*(dj+2.d0)*(dj+1.d0)/6.d0
      fak1=dj+4.d0
      fak2=(dj+4.d0)*(dj+3.d0)/2.d0
      fak3=(dj+4.d0)*(dj+3.d0)*(dj+2.d0)/6.d0
c
      qnorm=3.d0/rhomean/(dj+dj+1.d0)
      do 51 m=0,j
      jm=j*(j+1)/2+m+1
      cpoma=fak1*(pbj3*cat1jm(jm)-paj3*cas1jm(jm))
     :     +fak2*(pbj2*cat2jm(jm)-paj2*cas2jm(jm))
     :     +fak3*(pbj1*cat3jm(jm)-paj1*cas3jm(jm))
      cpomb=fac1*(pbj2*cbt1jm(jm)-paj2*cbs1jm(jm))
     :     +fac2*(pbj1*cbt2jm(jm)-paj1*cbs2jm(jm))
     :     +fac3*(pbj*cbt3jm(jm)-paj*cbs3jm(jm))
   51 cpotjm(jm)=cpotjm(jm)+qnorm*(rmean*cpoma/fak1+cpomb/fac1)
      return
      end
c============================================
      subroutine DEGPWR(jmax,fac,cjm,degpwrj)
c
c     Degree powers
c============================================
      implicit none
      integer*4 jmax,j,m,jm
      real*8 fac,degpwrj,pom
      complex*16 cjm,cpom,cpwrj
      dimension cjm(*),degpwrj(*)
c
      do 3 j=0,jmax
      cpwrj=(0.d0,0.d0)
      do 2 m=0,j
      jm=j*(j+1)/2+m+1
      cpom=cjm(jm)*dconjg(cjm(jm))
      if(m.eq.0) then
         cpwrj=cpwrj+cpom
      else
         cpwrj=cpwrj+cpom+cpom
      endif
    2 continue
      pom=dreal(cpwrj)
      pom=dsqrt(pom)
      degpwrj(j+1)=fac*pom
    3 continue
      return
      end
C..................................................................
      SUBROUTINE HARMAN(NTH,JMAX,DATA,CRHS)
C..................................................................
C     THIS SUBROUTINE CREATES THE ARRAY OF THE RIGHT-HAND SIDES
C     OF NORMAL EQUATIONS OCCURING IN THE LEAST SQUARES ADJUSTMENT
C
C     DESCRIPTION OF PARAMETERS:
C     NTH   THE NUMBER OF THE LATITUDE CIRCLES (NYQUIST FREQUENCY)
C           OF AN EQUAL ANGULAR GRID. THE STEP OF THE GRID IS THEN
C           STEP=180/NTH. IN ORDER TO ENSURE THE SYMMETRY OF THE
C           GRID WITH RESPECT TO THE EQUATOR, NTH NUST BE AN EVEN
C           NUMBER.
C     JMAX  THE DESIRED CUT-OFF DEGREE OF THE TRUNCATED SPHERICAL
C           HARMONIC SERIES FOR THE SPECTRAL REPRESENTATION OF A
C           FUNCTION. TO AVOID THE DISTORTION OF THE SPHERICAL
C           HARMONIC COEFFICIENTS BY THE ALIAS EFFECT, IT IS NECE-
C           SSARY TO PUT JMAX.LT.NTH.
C     DATA  TWO DIMENSIONAL ARRAY OF INPUT DATA VALUES. EACH COLUMN
C           OF THIS ARRAY CONTAINS DATA POINTS ALONG ONE PARALLEL.
C           THE PARALLELS ARE STORED FROM SOUTH TO NORTH POLE IN
C           SUCH A WAY THAT THE FIRST COLUMN IS THE CIRCLE WITH THE
C           180-STEP/2 CO-LATITUDE, AND THE LAST COLUMN OF THE
C           DATA ARRAY CONTAINS THE NORTHEST CIRCLE WITH THE
C           CO-LATITUDE STEP/2. THE FIRST ROW OF THE DATA ARRAY
C           CORRESPONDS TO THE ZERO-MERIDIAN. THE LONGITUDE
C           OF THE OTHER MERIDIANS INCREASES EASTWARDS WITH THE
C           STEP-INCREAMENT.
C     CRHS  THE RIGHT-HAND SIDES OF NORMAL EQUATIONS. DUE
C           TO EVEN-ODD SYMMETRY FOR THE NEGATIVE ORDERS M, ONLY
C           THE COEFFICIENTS WITH NON-NEGATIVE ORDERS ARE COMPUTED.
C           THE COUPLE OF INDICES (J,M), J=0,1,...,JMAX; M=0,1,
C           ...,J, IS REDUCED TO THE ONE-DIMENSIONAL INDEX JM AS:
C
C    (1)    JM=J*(J+1)/2+M+1
C
C           THEREFORE, THE COEFFICIENTS ARE ARRANGED FIRST BY THE
C           DEGREE J, AND THEN BY THE ORDER M: 00,10,11,20,21,22,
C           ...,(JMAX,JMAX).
C     DIMENSIONS:
C           DATA(2*NTH,NTH),CRHS((JMAX+1)*(JMAX+2)/2)
C     DIMENS. OF AUXILIARY ARRAYS:
C           PNM(JMAX+1),RE(2*NTH+1),QIM(2*NTH+1)
C     NOTE:
C           THE PRESENT VALUES OF THE DIMENSIONS OF THE AUXILIARY
C           ARRAYS ARE SET-UP FOR NTH.LE.180 AND JMAX.LE.180.
C           IF THESE BOUNDS ARE OVERCOME, THE DIMENSIONS MUST BE
C           INCREASED.
C     THE USED SUBROUTINES:
C           C06ECF - CALCULATES FAST FOURIER TRANSFORM OF MIX-RADIX.
C                    IT WAS TAKEN FROM THE DOUBLE PRECISION NAG
C                    SUBROUTINE LIBRARY.
C           DPMM   - COMPUTES THE FULLY NORMALIZED ASSOCIATED LEGENDRE
C                    FUNCTIONS.
C     REFERENCE
C
C     Martinec, Z., 1991. Program to calculate the least-squares estimates 
C     of the spherical harmonic expansion coefficients of an equally 
C     angular-gridded scalar field. Computer Physics Communications, 64, 140–148.
C........................................................
      IMPLICIT REAL*8(A,B,D-H,O-Z), COMPLEX*16(C)
      DIMENSION DATA(720,360),CRHS(*)
      DIMENSION PNM(361),RE(1441),QIM(1441)
C
C     OUTPUT FORMAT
  200 FORMAT(1X,'IFAIL=',I2)
C
C     SETTING UP OF THE INITIAL VALUES
      STEP=180.D00/DFLOAT(NTH)
      NPH=2*NTH
      NTH2=NTH/2
      FAC=DATAN(1.D00)/45.D00
      TH0=180.D00-STEP/2.D00
      JMAX1=JMAX+1
      DNFR=DFLOAT(NPH)
      SNFR=DSQRT(DNFR)
C
C     SUMMATION OVER MERIDIANS
      DO 5 ITH2=1,NTH2
      THS=TH0-(ITH2-1)*STEP
      TH=THS*FAC
      XTH=DCOS(TH)
      ITHF=ITH2
      ITHS=NTH-ITH2+1
C
C     MIX-RADIX FAST FOURIER TRANSFORM ALONG EACH PARALLEL
         DO 1 IPH=1,NPH
         RE(IPH)=DATA(IPH,ITHF)
    1    QIM(IPH)=DATA(IPH,ITHS)
      IFAIL=0
      CALL C06ECF(RE,QIM,NPH,IFAIL)
      IF(IFAIL.EQ.0) GOTO 2
      WRITE(66,200) IFAIL
      STOP
    2 CONTINUE
      RE(NPH+1)=RE(1)
      QIM(NPH+1)=QIM(1)
C
C     CREATION OF THE RIGHT-HAND SIDES
         DO 3 MS=1,JMAX1
         M=MS-1
         CALL DPMM(XTH,JMAX,M,PNM)
         NFMS=NPH+2-MS
         POM1=(RE(MS)+RE(NFMS))/2.D00
         POM2=(QIM(MS)-QIM(NFMS))/2.D00
         CAF=DCMPLX(POM1,POM2)*SNFR
         POM1=(QIM(MS)+QIM(NFMS))/2.D00
         POM2=(-RE(MS)+RE(NFMS))/2.D00
         CAS=DCMPLX(POM1,POM2)*SNFR
         ZN=-1.D00
            DO 3 JS=MS,JMAX1
            ZN=-ZN
            JM=(JS-1)*JS/2+MS
            JMP=JS-MS+1
            IF(ITH2.EQ.1) CRHS(JM)=(0.D00,0.D00)
    3       CRHS(JM)=CRHS(JM)+PNM(JMP)*(CAF+CAS*ZN)
    5 CONTINUE
      RETURN
      END
C..............................................................
      SUBROUTINE HARMLS(NTH,JMAX,CRHS,COEF)
C..............................................................
C     THE SUBROUTINE COMPUTES THE OPTIMAL LEAST SQUARES ESTIMATES
C     OF THE SPHERICAL HARMONIC COEFFICIENTS. THE ALGORITHM EXPLOITS
C     THE FACT THAT THE NORMAL MATRIX IS SPARSE FOR REGULARLY SPACED
C     GRID. THE NORMAL MATRIX IS REORDERED INTO THE BLOCK-DIAGONAL
C     FORM AND THEN EACH OF THE SUB-SYSTEM IS SOLVED BY THE NAG
C     SUBROUTINE F04ABF PERFORMING CHOLESKY'S FACTORIZATION.
C
C     DESCRIPTION OF PARAMETERS:
C     NTH    THE SAME AS IN HARMAN
C     JMAX   THE SAME AS IN HARMAN
C     CRHS   THE SAME AS IN HARMAN
C     COEF   THE LEAST SQUARES ESTIMATES OF THE SPHERICAL HARMONIC
C            COEFFICIENTS. THE ONE-DIMENSIONAL ARRAY COEF IS ARRANGED
C            ACCORDING TO THE INDEX (1)-VIZ SUBROUTINE HARMAN.
C     DIMENSIONS:
C            CRHS((JMAX+1)*(JMAX+2)/2),COEF((JMAX+1)*(JMAX+2)/2)
C     DIMENS. OF AUXILIARY ARRAYS:
C            PNM(JMAX+1),AMTRX(JMAX/2+1,JMAX/2+1),AUX(JMAX/2+1,2)
C            PARAM(JMAX/2+1,2),BB(JMAX/2+1,2),WORK(JMAX/2+1)
C     NOTE:
C            THE PRESENT VALUES OF THE DIMENSIONS OF THE AUXILIARY
C            ARRAYS ARE SET-UP FOR NTH.LE.180 AND JMAX.LE.180.
C            IF THESE BOUNDS ARE OVERCOME, THE DIMENSIONS MUST BE
C            INCREASED.
C     THE USED SUBROUTINES:
C            FO4ABF - CALCULATES THE ACCURATE SOLUTION OF A SET OF
C                     REAL SYMMETRIC POSITIVE DEFINITE LINEAR
C                     EQUATIONS WITH MULTIPLE RIGHT HAND SIDES BY
C                     CHOLESKY'S DECOMPOSITION METHOD.
C                     IT WAS TAKEN FROM THE DOUBLE PRECISION NAG
C                     SUBROUTINE LIBRARY.
C            DPMM   - COMPUTES THE FULLY NORMALIZED ASSOCIATED LEGENDRE
C                     FUNCTIONS.
C     REFERENCE
C
C     Martinec, Z., 1991. Program to calculate the least-squares estimates 
C     of the spherical harmonic expansion coefficients of an equally 
C     angular-gridded scalar field. Computer Physics Communications, 64, 140–148.
C..............................................................
      IMPLICIT REAL*8(A,B,D-H,O-Z), COMPLEX*16(C)
      DIMENSION CRHS(*),COEF(*)
      DIMENSION PNM(361),AMTRX(181,183),PARAM(181,2)
C
  200 FORMAT(1X,'IFAIL=',I2)
C
      STEP=180.D00/DFLOAT(NTH)
      NTH2=NTH/2
      FAC=DATAN(1.D00)/45.D00
      JMAX1=JMAX+1
      TH0=180.D00-STEP/2.D00
C
      DO 6 M1S=1,JMAX1
      M1=M1S-1
      M1S1=M1S+1
      DO 6 MP=M1S,M1S1
      NPM=(JMAX1-MP+2)/2
      IF(MP.GT.JMAX1) GOTO 6
         DO 2 ITH=1,NTH2
         THS=TH0-(ITH-1)*STEP
         TH=THS*FAC
         XTH=DCOS(TH)
         CALL DPMM(XTH,JMAX,M1,PNM)
            IJ1=0
            DO 1 J1S=MP,JMAX1,2
            IJ1=IJ1+1
            J1M1=J1S-M1S+1
            POM1=4.D00*NTH*PNM(J1M1)
            IF(DABS(POM1).LT.1.D-30) POM1=0.D00
            IJ2=0
            DO 1 J2S=MP,JMAX1,2
            IJ2=IJ2+1
            J2M1=J2S-M1S+1
            POM2=PNM(J2M1)
            IF(DABS(POM2).LT.1.D-30) POM2=0.D00
            IF(ITH.EQ.1) AMTRX(IJ1,IJ2)=0.D00
            AMTRX(IJ1,IJ2)=AMTRX(IJ1,IJ2)+POM1*POM2
    1       CONTINUE
    2    CONTINUE
            IJ1=0
            DO 3 J1S=MP,JMAX1,2
            IJ1=IJ1+1
            J1M1=(J1S-1)*J1S/2+M1S
            AMTRX(IJ1,NPM+1)=DBLE(CRHS(J1M1))
    3       AMTRX(IJ1,NPM+2)=DIMAG(CRHS(J1M1))
      CALL SIMZM2(NPM,AMTRX,2,PARAM)
         IJ1=0
         DO 5 J1S=MP,JMAX1,2
         IJ1=IJ1+1
         J1M1=(J1S-1)*J1S/2+M1S
    5    COEF(J1M1)=DCMPLX(PARAM(IJ1,1),PARAM(IJ1,2))
    6 CONTINUE
      RETURN
      END
C..................................................................
      SUBROUTINE HARMSY(NTH,JMAX,COEF,DATA)
C..................................................................
C     THIS SUBROUTINE EVALUATES THE TRUNCATED SUM OF THE SPHERICAL
C     HARMONICS COMPLETE TO DEGREE AND ORDER JMAX AT EACH ONE OF
C     THE 2*NTH**2 POINTS IN THE GRID.
C
C     DESCRIPTION OF PARAMETERS:
C          NTH   THE SAME AS IN HARMAN
C          JMAX  THE SAME AS IN HARMAN
C          COEF  THE SAME AS IN HARMLS
C          DATA  TWO-DIMENSIONAL ARRAY OF THE DATA VALUES COMPUTED
C                IN THE NODES OF THE REGULAR GRID. THIS ARRAY IS
C                ORGANIZED IN THE SAME WAY AS THE ARRAY DATA IN THE
C                SUBROUTINE HARMAN.
C     DIMENSIONS:
C           DATA(2*NTH,NTH),COEF((JMAX+1)*(JMAX+2)/2)
C     DIMENS. OF AUXILIARY ARRAYS:
C           PNM(JMAX+1),RE1(2*NTH+1),QIM1(2*NTH+1),RE2(2*NTH+1),
C           QIM2(2*NTH+1)
C     NOTE:
C           THE PRESENT VALUES OF THE DIMENSIONS OF THE AUXILIARY
C           ARRAYS ARE SET-UP FOR NTH.LE.180 AND JMAX.LE.180.
C           IF THESE BOUNDS ARE OVERCOME, THE DIMENSIONS MUST BE
C           INCREASED.
C     THE USED SUBROUTINES:
C           C06ECF - CALCULATES FAST FOURIER TRANSFORM OF MIX-RADIX.
C                    IT WAS TAKEN FROM THE DOUBLE PRECISION NAG
C                    SUBROUTINE LIBRARY.
C           DPMM   - COMPUTES THE FULLY NORMALIZED ASSOCIATED LEGENDRE
C                    FUNCTIONS.
C     REFERENCE
C
C     Martinec, Z., 1991. Program to calculate the least-squares estimates 
C     of the spherical harmonic expansion coefficients of an equally 
C     angular-gridded scalar field. Computer Physics Communications, 64, 140–148.
C..................................................................
      IMPLICIT REAL*8(A,B,D-H,O-Z), COMPLEX*16(C)
      DIMENSION DATA(720,360),COEF(*)
      DIMENSION PNM(361),RE1(1441),QIM1(1441),RE2(1441),QIM2(1441)
C
C     OUTPUT FORMAT
  200 FORMAT(1X,'IFAIL=',I2)
C
C     SETTING UP OF THE INITIAL VALUES
      STEP=180.D00/DFLOAT(NTH)
      NPH=2*NTH
      NTH2=NTH/2
      FAC=DATAN(1.D00)/45.D00
      TH0=180.D00-STEP/2.D00
      JMAX1=JMAX+1
      JMAX2=JMAX+2
      DNFR=DFLOAT(NPH)
      SNFR=DSQRT(DNFR)
C
      DO 5 ITH2=1,NTH2
      THS=TH0-(ITH2-1)*STEP
      TH=THS*FAC
      XTH=DCOS(TH)
         DO 2 MS=1,JMAX1
         M=MS-1
         CALL DPMM(XTH,JMAX,M,PNM)
         CSUM1=(0.D00,0.D00)
         CSUM2=(0.D00,0.D00)
         ZN=-1.D00
            DO 1 JS=MS,JMAX1
            ZN=-ZN
            JM=(JS-1)*JS/2+MS
            JMP=JS-MS+1
            CSUM1=CSUM1+COEF(JM)*PNM(JMP)
    1       CSUM2=CSUM2+COEF(JM)*PNM(JMP)*ZN
         RE1(MS)=DBLE(CSUM1)
         QIM1(MS)=-DIMAG(CSUM1)
         RE2(MS)=DBLE(CSUM2)
    2    QIM2(MS)=-DIMAG(CSUM2)
      RE1(1)=RE1(1)/2.D00
      RE2(1)=RE2(1)/2.D00
         DO 3 MS=JMAX2,NPH
         RE1(MS)=0.D00
         QIM1(MS)=0.D00
         RE2(MS)=0.D00
    3    QIM2(MS)=0.D00
      IFAIL=0
      CALL C06ECF(RE1,QIM1,NPH,IFAIL)
      CALL C06ECF(RE2,QIM2,NPH,IFAIL)
      IF(IFAIL.EQ.0) GOTO 4
      WRITE(66,200) IFAIL
      STOP
    4 CONTINUE
         DO 5 IPH=1,NPH
         DATA(IPH,ITH2)=2.D00*RE1(IPH)*SNFR
    5    DATA(IPH,NTH-ITH2+1)=2.D00*RE2(IPH)*SNFR
      RETURN
      END
C..................................................................
      SUBROUTINE DPMM(X,N,M,P)
C..................................................................
C     THIS SUBROUTINE RETURNS THE ARRAY OF THE FULLY NORMALIZED
C     ASSOCIATED LEGENDRE FUNCTIONS FOR THE GIVEN ORDER M
C
C      DESCRIPTION OF PARAMETERS:
C
C         X      THE ARGUMENT OF THE LEGENDRE FUNCTIONS, X=COS(TH)
C         N      THE HIGHEST DEGREE UP TO WHICH THE LEGENDRE FUNCTIONS
C                ARE COMPUTED
C         M      THE DESIRED ORDER FOR WHICH THE LEGENDRE FUNCTIONS
C                ARE COMPUTED
C         P      THE ONE-DIMENSIONAL, REAL*8, ARRAY OF THE ASSOCIATED
C                LEGENDRE FUNCTIONS
C     METHOD
C                RECURSION RELATION:
C         P(J,M)=DSQRT((2*J+1)*(2*J-1)/(J-M)/(J+M))*X*P(J-1,M)-
C         - DSQRT((2*J+1)*(J-M-1)*(J+M-1)/(2*J-3)/(J-M)/(J+M))*P(J-2,M)
C                STARTING VALUES:
C         P(M,M)=(-1)**M*DSQRT((2*M+1)]]/(4*PI)/(2M)]])*(1-X**2)**M/2
C         P(M+1,M)=DSQRT(2*M+3)*X*P(M,M)
C                INDEXING:
C         P(J,M)=P(J-M+1), J=M,M+1,...,N
C     DIMENSION:
C         P(N+1)
C
C     REFERENCE
C
C     Martinec, Z., 1991. Program to calculate the least-squares estimates 
C     of the spherical harmonic expansion coefficients of an equally 
C     angular-gridded scalar field. Computer Physics Communications, 64, 140–148.
C.............................................................
      IMPLICIT REAL*8(D-H,O-Z)
      DIMENSION P(*)
      STH=DSQRT(1.D00-X*X)
      PI4=16.D00*DATAN(1.D00)
      P(1)=1.D00/DSQRT(PI4)
      IF(N.LE.0) RETURN
      F1=0.D00
      IF(M.EQ.0) GOTO 1
         SOU=1.D00
         IFLAG=0
         DO 3 MP=1,M
         F1=DFLOAT(MP+MP)
         SOU=SOU*(F1+1.D00)/F1
         IF(IFLAG-1) 10,20,20
   10    STHM=STH**MP
         IF(STHM.LT.1.D-55) IFLAG=1
         GOTO 25
   20    STHM=0.D00
   25    CONTINUE
    3    CONTINUE
      P(1)=(-1.D00)**M*DSQRT(SOU/PI4)*STHM
      IF(M.EQ.N) GOTO 33
    1 P(2)=DSQRT(F1+3.D00)*X*P(1)
      MP2=M+2
      IF(MP2.GT.N) GOTO 33
         I=1
         DO 2 J=MP2,N
         F1=DFLOAT((J-M)*(J+M))
         F2=DFLOAT((J+J-1)*(J+J+1))
         F3=DFLOAT((J+J+1)*(J-M-1)*(J+M-1))/DFLOAT(J+J-3)
         F2=DSQRT(F2/F1)*X
         F3=DSQRT(F3/F1)
         I=I+1
         P(I+1)=F2*P(I)-F3*P(I-1)
    2    CONTINUE
   33 CONTINUE
      RETURN
      END
c....................................................................
      SUBROUTINE SIMZM2(N,A,NRHS,X)
c     Gauss elimination - no pivoting
C....................................................................
      REAL*8 A,X,PIVOT,AIK
      DIMENSION A(181,183),X(181,2)
      MAX=N+NRHS
C..... BEGIN ELIMINATION PROCEDURE .....
      DO 18 K=1,N
      PIVOT=A(K,K)
c      write(66,900) k,pivot
c  900 format(1x,'k-pivot',i3,2x,d13.5)
C..... NORMALIZE PIVOT ROW ELEMENTS .....
      DO 14 J=1,MAX
   14 A(K,J)=A(K,J)/PIVOT
C..... CARRY OUT ELIMINATION  .......
      A(K,K)=1.D0/PIVOT
      DO 18 I=1,N
      AIK=A(I,K)
      IF(I.EQ.K) GO TO 18
      A(I,K)=-AIK/PIVOT
      DO 17 J=1,MAX
   17 IF(J.NE.K) A(I,J)=A(I,J)-AIK*A(K,J)
   18 CONTINUE
C..... ORDER SOLUTION VALUES .....
      DO 20 IRHS=1,NRHS
      DO 20 I=1,N
   20 X(I,IRHS)=A(I,N+IRHS)
      RETURN
      END
c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE C06ECF(X, Y, PTS, IFAIL)
c
c     FAST FOURIER TRANSFORM SUBROUTINES
c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     COMPLEX FOURIER TRANSFORM
C     .. SCALAR ARGUMENTS ..
      INTEGER IFAIL, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X(PTS), Y(PTS)
C     ..
C     .. LOCAL SCALARS ..
C      DOUBLE PRECISION SRNAME
      REAL*8 SQPTS
      INTEGER IERROR, IPTS, PMAX, PSYM, TWOGRP
C     .. LOCAL ARRAYS ..
      INTEGER FACTOR(21), SYM(21), UNSYM(21)
C     .. FUNCTION REFERENCES ..
      REAL*8 DSQRT
C      INTEGER P01AAE
C     .. SUBROUTINE REFERENCES ..
C     EAZC06, ECWC06, ECYC06
C     ..
C      DATA SRNAME /'  C06ECE'/
      DATA PMAX /19/
      DATA TWOGRP /8/
      IF (PTS.LE.1) GO TO 40
      IERROR = 0
      CALL EAZC06(PTS, PMAX, TWOGRP, FACTOR, SYM, PSYM, UNSYM,
     * IERROR)
      IF (IERROR.NE.0) GO TO 60
      CALL ECWC06(X, Y, PTS, FACTOR)
      CALL ECYC06(X, Y, PTS, SYM, PSYM, UNSYM)
      SQPTS = DSQRT(DFLOAT(PTS))
      DO 20 IPTS=1,PTS
         X(IPTS) = X(IPTS)/SQPTS
         Y(IPTS) = Y(IPTS)/SQPTS
   20 CONTINUE
      IFAIL = 0
      GO TO 80
C
   40 IERROR = 3
C  60 IFAIL = P01AAE(IFAIL,IERROR,SRNAME)
   60 CONTINUE
   80 RETURN
      END
      SUBROUTINE ECWC06(X, Y, PTS, FACTOR)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     COMPLEX FOURIER TRANSFORM KERNEL DRIVER
C     .. SCALAR ARGUMENTS ..
      INTEGER PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X(PTS), Y(PTS)
      INTEGER FACTOR(21)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER F, M1, M2, M3, M4, M5, M6, M7, M, P
C     .. SUBROUTINE REFERENCES ..
C     ECQC06, ECRC06, ECSC06, ECTC06, ECUC06, ECVC06
C     ..
      F = 0
      M = PTS
   20 CONTINUE
      F = F + 1
      P = FACTOR(F)
      IF (P.EQ.0) RETURN
      IF (P.EQ.1) GO TO 20
      M = M/P
      M1 = PTS - M
      M2 = M1 - M
      M3 = M2 - M
      M4 = M3 - M
      M5 = M4 - M
      M6 = M5 - M
      M7 = M6 - M
      IF (P.EQ.2) GO TO 40
      IF (P.EQ.3) GO TO 60
      IF (P.EQ.4) GO TO 80
      IF (P.EQ.5) GO TO 100
      IF (P.EQ.8) GO TO 120
      GO TO 140
C
   40 CONTINUE
      CALL ECVC06(X(1), Y(1), PTS, X(M+1), Y(M+1), M1, M)
      GO TO 20
C
   60 CONTINUE
      CALL ECUC06(X(1), Y(1), PTS, X(M+1), Y(M+1), M1, X(2*M+1),
     * Y(2*M+1), M2, M)
      GO TO 20
C
   80 CONTINUE
      CALL ECTC06(X(1), Y(1), PTS, X(M+1), Y(M+1), M1, X(2*M+1),
     * Y(2*M+1), M2, X(3*M+1), Y(3*M+1), M3, M)
      GO TO 20
C
  100 CONTINUE
      CALL ECSC06(X(1), Y(1), PTS, X(M+1), Y(M+1), M1, X(2*M+1),
     * Y(2*M+1), M2, X(3*M+1), Y(3*M+1), M3, X(4*M+1), Y(4*M+1),
     * M4, M)
      GO TO 20
C
  120 CONTINUE
      CALL ECRC06(X(1), Y(1), PTS, X(M+1), Y(M+1), M1, X(2*M+1),
     * Y(2*M+1), M2, X(3*M+1), Y(3*M+1), M3, X(4*M+1), Y(4*M+1),
     * M4, X(5*M+1), Y(5*M+1), M5, X(6*M+1), Y(6*M+1), M6,
     * X(7*M+1), Y(7*M+1), M7, M)
      GO TO 20
C
  140 CONTINUE
      CALL ECQC06(X, Y, PTS, M, P)
      GO TO 20
C
      END
      SUBROUTINE ECYC06(X, Y, PTS, SYM, PSYM, UNSYM)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     DOUBLE IN PLACE REORDERING PROGRAMME
C     .. SCALAR ARGUMENTS ..
      INTEGER PSYM, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X(PTS), Y(PTS)
      INTEGER SYM(21), UNSYM(21)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 T
      INTEGER DK, I, II, IL, J, JJ, JL, K, KK, KS, LK, MODS, MULT,
     * NEST, PUNSYM, TEST
      LOGICAL ONEMOD
C     .. LOCAL ARRAYS ..
      INTEGER MODULO(20)
C     .. SUBROUTINE REFERENCES ..
C     ECXC06
C     ..
      DATA NEST /20/
      CALL ECXC06(X, Y, PTS, SYM)
C
      IF (UNSYM(1).EQ.0) GO TO 280
      PUNSYM = PTS/PSYM**2
      MULT = PUNSYM/UNSYM(1)
      TEST = (UNSYM(1)*UNSYM(2)-1)*MULT*PSYM
      LK = MULT
      DK = MULT
      DO 20 K=2,NEST
         IF (UNSYM(K).EQ.0) GO TO 40
         LK = LK*UNSYM(K-1)
         DK = DK/UNSYM(K)
         MODULO(K) = (LK-DK)*PSYM
         MODS = K
   20 CONTINUE
   40 CONTINUE
      ONEMOD = MODS.LT.3
      IF (ONEMOD) GO TO 80
      K = (MODS+3)/2
      DO 60 J=3,K
         JJ = MODS + 3 - J
         KK = MODULO(J)
         MODULO(J) = MODULO(JJ)
         MODULO(JJ) = KK
   60 CONTINUE
   80 CONTINUE
      JL = (PUNSYM-3)*PSYM
      KS = PUNSYM*PSYM
C
      DO 260 J=PSYM,JL,PSYM
         JJ = J
C
  100    CONTINUE
         JJ = JJ*MULT
         IF (ONEMOD) GO TO 140
         DO 120 I=3,MODS
            JJ = JJ - (JJ/MODULO(I))*MODULO(I)
  120    CONTINUE
  140    CONTINUE
         IF (JJ.GE.TEST) GO TO 160
         JJ = JJ - (JJ/MODULO(2))*MODULO(2)
         GO TO 180
  160    CONTINUE
         JJ = JJ - (JJ/MODULO(2))*MODULO(2) + MODULO(2)
  180    CONTINUE
         IF (JJ.LT.J) GO TO 100
C
         IF (JJ.EQ.J) GO TO 240
         LK = JJ - J
         II = J + 1
         IL = J + PSYM
         DO 220 I=II,IL
            DO 200 K=I,PTS,KS
               KK = K + LK
               T = X(K)
               X(K) = X(KK)
               X(KK) = T
               T = Y(K)
               Y(K) = Y(KK)
               Y(KK) = T
  200       CONTINUE
  220    CONTINUE
  240    CONTINUE
  260 CONTINUE
C
  280 CONTINUE
      RETURN
      END
      SUBROUTINE ECQC06(X, Y, PTS, M, P)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     RADIX PRIME COMPLEX FOURIER TRANSFORM KERNEL
C     .. SCALAR ARGUMENTS ..
      INTEGER M, P, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X(PTS), Y(PTS)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 ANGLE, IS, IU, RS, RU, T, TWOPI, XT, YT
      INTEGER J, JJ, K0, K, KS1, KS2, MOVER2, MP, PM, PP, U, V
      LOGICAL FOLD, ZERO
C     .. LOCAL ARRAYS ..
      REAL*8 A(18), AA(9,9), B(18), BB(9,9), C(18), IA(9), IB(9),
     * RA(9), RB(9), S(18)
C     .. FUNCTION REFERENCES ..
      REAL*8 DCOS, DSIN
C     ..
      twopi=8.0*datan(1.d00)
      MOVER2 = M/2 + 1
      MP = M*P
      PP = P/2
      PM = P - 1
      DO 20 U=1,PP
         JJ = P - U
         ANGLE = TWOPI*FLOAT(U)/FLOAT(P)
         A(U) = DCOS(ANGLE)
         B(U) = DSIN(ANGLE)
         A(JJ) = A(U)
         B(JJ) = -B(U)
   20 CONTINUE
      DO 60 U=1,PP
         DO 40 V=1,PP
            JJ = U*V - ((U*V)/P)*P
            AA(V,U) = A(JJ)
            BB(V,U) = B(JJ)
   40    CONTINUE
   60 CONTINUE
C
      DO 300 J=1,MOVER2
         FOLD = J.GT.1 .AND. 2*J.LT.M + 2
         K0 = J
         ANGLE = TWOPI*FLOAT(J-1)/FLOAT(MP)
         ZERO = ANGLE.EQ.0.0
         C(1) = DCOS(ANGLE)
         S(1) = DSIN(ANGLE)
         DO 80 U=2,PM
            C(U) = C(U-1)*C(1) - S(U-1)*S(1)
            S(U) = S(U-1)*C(1) + C(U-1)*S(1)
   80    CONTINUE
         GO TO 140
  100    CONTINUE
         FOLD = .FALSE.
         K0 = M + 2 - J
         DO 120 U=1,PM
            T = C(U)*A(U) + S(U)*B(U)
            S(U) = -S(U)*A(U) + C(U)*B(U)
            C(U) = T
  120    CONTINUE
  140    CONTINUE
C
         DO 280 K=K0,PTS,MP
            XT = X(K)
            YT = Y(K)
            KS1 = M + K
            KS2 = (P-1)*M + K
            RS = X(KS1) + X(KS2)
            IS = Y(KS1) + Y(KS2)
            RU = X(KS1) - X(KS2)
            IU = Y(KS1) - Y(KS2)
            DO 160 U=1,PP
               RA(U) = XT + RS*AA(U,1)
               IA(U) = YT + IS*AA(U,1)
               RB(U) = RU*BB(U,1)
               IB(U) = IU*BB(U,1)
  160       CONTINUE
            XT = XT + RS
            YT = YT + IS
            DO 200 U=2,PP
               JJ = P - U
               KS1 = U*M + K
               KS2 = JJ*M + K
               RS = X(KS1) + X(KS2)
               IS = Y(KS1) + Y(KS2)
               RU = X(KS1) - X(KS2)
               IU = Y(KS1) - Y(KS2)
               XT = XT + RS
               YT = YT + IS
               DO 180 V=1,PP
                  RA(V) = RA(V) + RS*AA(V,U)
                  IA(V) = IA(V) + IS*AA(V,U)
                  RB(V) = RB(V) + RU*BB(V,U)
                  IB(V) = IB(V) + IU*BB(V,U)
  180          CONTINUE
  200       CONTINUE
            X(K) = XT
            Y(K) = YT
            DO 260 U=1,PP
               JJ = P - U
               IF (ZERO) GO TO 220
               XT = RA(U) + IB(U)
               YT = IA(U) - RB(U)
               KS1 = U*M + K
               X(KS1) = XT*C(U) + YT*S(U)
               Y(KS1) = YT*C(U) - XT*S(U)
               XT = RA(U) - IB(U)
               YT = IA(U) + RB(U)
               KS1 = JJ*M + K
               X(KS1) = XT*C(JJ) + YT*S(JJ)
               Y(KS1) = YT*C(JJ) - XT*S(JJ)
               GO TO 240
  220          CONTINUE
               KS1 = U*M + K
               X(KS1) = RA(U) + IB(U)
               Y(KS1) = IA(U) - RB(U)
               KS1 = JJ*M + K
               X(KS1) = RA(U) - IB(U)
               Y(KS1) = IA(U) + RB(U)
  240          CONTINUE
  260       CONTINUE
  280    CONTINUE
         IF (FOLD) GO TO 100
  300 CONTINUE
C
      RETURN
      END
      SUBROUTINE ECRC06(X0, Y0, PTS, X1, Y1, M1, X2, Y2, M2, X3,
     * Y3, M3, X4, Y4, M4, X5, Y5, M5, X6, Y6, M6, X7, Y7, M7, M)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     RADIX EIGHT COMPLEX FOURIER TRANSFORM KERNEL
C     .. SCALAR ARGUMENTS ..
      INTEGER M1, M2, M3, M4, M5, M6, M7, M, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X0(PTS), X1(M1), X2(M2), X3(M3), X4(M4), X5(M5), X6(M6),
     * X7(M7), Y0(PTS), Y1(M1), Y2(M2), Y3(M3), Y4(M4), Y5(M5),
     * Y6(M6), Y7(M7)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 ANGLE, C1, C2, C3, C4, C5, C6, C7, E, I1, I2, I3, I4,
     * I5, I6, I7, IS0, IS1, IS2, IS3, ISS0, ISS1, ISU0, ISU1, IU0,
     * IU1, IU2, IU3, IUS0, IUS1, IUU0, IUU1, R1, R2, R3, R4, R5,
     * R6, R7, RS0, RS1, RS2, RS3, RSS0, RSS1, RSU0, RSU1, RU0,
     * RU1, RU2, RU3, RUS0, RUS1, RUU0, RUU1, S1, S2, S3, S4, S5,
     * S6, S7, T, TWOPI
      INTEGER J, K0, K, M8, MOVER2
      LOGICAL FOLD, ZERO
C     .. FUNCTION REFERENCES ..
      REAL*8 DCOS, DSIN
C     ..
      M8 = M*8
      MOVER2 = M/2 + 1
      twopi=8.0*datan(1.d00)
      E = DCOS(TWOPI/8.0D00)
C
      DO 120 J=1,MOVER2
         FOLD = J.GT.1 .AND. 2*J.LT.M + 2
         K0 = J
         ANGLE = TWOPI*FLOAT(J-1)/FLOAT(M8)
         ZERO = ANGLE.EQ.0.0
         C1 = DCOS(ANGLE)
         S1 = DSIN(ANGLE)
         C2 = C1*C1 - S1*S1
         S2 = S1*C1 + C1*S1
         C3 = C2*C1 - S2*S1
         S3 = S2*C1 + C2*S1
         C4 = C2*C2 - S2*S2
         S4 = S2*C2 + C2*S2
         C5 = C4*C1 - S4*S1
         S5 = S4*C1 + C4*S1
         C6 = C4*C2 - S4*S2
         S6 = S4*C2 + C4*S2
         C7 = C4*C3 - S4*S3
         S7 = S4*C3 + C4*S3
         GO TO 40
   20    CONTINUE
         FOLD = .FALSE.
         K0 = M + 2 - J
         T = (C1+S1)*E
         S1 = (C1-S1)*E
         C1 = T
         T = S2
         S2 = C2
         C2 = T
         T = (-C3+S3)*E
         S3 = (C3+S3)*E
         C3 = T
         C4 = -C4
         T = -(C5+S5)*E
         S5 = (-C5+S5)*E
         C5 = T
         T = -S6
         S6 = -C6
         C6 = T
         T = (C7-S7)*E
         S7 = -(C7+S7)*E
         C7 = T
   40    CONTINUE
C
         DO 100 K=K0,PTS,M8
            RS0 = X0(K) + X4(K)
            IS0 = Y0(K) + Y4(K)
            RU0 = X0(K) - X4(K)
            IU0 = Y0(K) - Y4(K)
            RS1 = X1(K) + X5(K)
            IS1 = Y1(K) + Y5(K)
            RU1 = X1(K) - X5(K)
            IU1 = Y1(K) - Y5(K)
            RS2 = X2(K) + X6(K)
            IS2 = Y2(K) + Y6(K)
            RU2 = X2(K) - X6(K)
            IU2 = Y2(K) - Y6(K)
            RS3 = X3(K) + X7(K)
            IS3 = Y3(K) + Y7(K)
            RU3 = X3(K) - X7(K)
            IU3 = Y3(K) - Y7(K)
            RSS0 = RS0 + RS2
            ISS0 = IS0 + IS2
            RSU0 = RS0 - RS2
            ISU0 = IS0 - IS2
            RSS1 = RS1 + RS3
            ISS1 = IS1 + IS3
            RSU1 = RS1 - RS3
            ISU1 = IS1 - IS3
            RUS0 = RU0 - IU2
            IUS0 = IU0 + RU2
            RUU0 = RU0 + IU2
            IUU0 = IU0 - RU2
            RUS1 = RU1 - IU3
            IUS1 = IU1 + RU3
            RUU1 = RU1 + IU3
            IUU1 = IU1 - RU3
            T = (RUS1+IUS1)*E
            IUS1 = (IUS1-RUS1)*E
            RUS1 = T
            T = (RUU1+IUU1)*E
            IUU1 = (IUU1-RUU1)*E
            RUU1 = T
            X0(K) = RSS0 + RSS1
            Y0(K) = ISS0 + ISS1
            IF (ZERO) GO TO 60
            R1 = RUU0 + RUU1
            I1 = IUU0 + IUU1
            R2 = RSU0 + ISU1
            I2 = ISU0 - RSU1
            R3 = RUS0 + IUS1
            I3 = IUS0 - RUS1
            R4 = RSS0 - RSS1
            I4 = ISS0 - ISS1
            R5 = RUU0 - RUU1
            I5 = IUU0 - IUU1
            R6 = RSU0 - ISU1
            I6 = ISU0 + RSU1
            R7 = RUS0 - IUS1
            I7 = IUS0 + RUS1
            X4(K) = R1*C1 + I1*S1
            Y4(K) = I1*C1 - R1*S1
            X2(K) = R2*C2 + I2*S2
            Y2(K) = I2*C2 - R2*S2
            X6(K) = R3*C3 + I3*S3
            Y6(K) = I3*C3 - R3*S3
            X1(K) = R4*C4 + I4*S4
            Y1(K) = I4*C4 - R4*S4
            X5(K) = R5*C5 + I5*S5
            Y5(K) = I5*C5 - R5*S5
            X3(K) = R6*C6 + I6*S6
            Y3(K) = I6*C6 - R6*S6
            X7(K) = R7*C7 + I7*S7
            Y7(K) = I7*C7 - R7*S7
            GO TO 80
   60       CONTINUE
            X4(K) = RUU0 + RUU1
            Y4(K) = IUU0 + IUU1
            X2(K) = RSU0 + ISU1
            Y2(K) = ISU0 - RSU1
            X6(K) = RUS0 + IUS1
            Y6(K) = IUS0 - RUS1
            X1(K) = RSS0 - RSS1
            Y1(K) = ISS0 - ISS1
            X5(K) = RUU0 - RUU1
            Y5(K) = IUU0 - IUU1
            X3(K) = RSU0 - ISU1
            Y3(K) = ISU0 + RSU1
            X7(K) = RUS0 - IUS1
            Y7(K) = IUS0 + RUS1
   80       CONTINUE
  100    CONTINUE
         IF (FOLD) GO TO 20
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE ECSC06(X0, Y0, PTS, X1, Y1, M1, X2, Y2, M2, X3,
     * Y3, M3, X4, Y4, M4, M)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     RADIX FIVE COMPLEX FOURIER TRANSFORM KERNEL
C     .. SCALAR ARGUMENTS ..
      INTEGER M1, M2, M3, M4, M, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X0(PTS), X1(M1), X2(M2), X3(M3), X4(M4), Y0(PTS),
     * Y1(M1), Y2(M2), Y3(M3), Y4(M4)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 A1, A2, ANGLE, AS, AU, B1, B2, C1, C2, C3, C4, I0, I1,
     * I2, I3, I4, IA1, IA2, IAS, IAU, IB1, IB2, IS1, IS2, ISS,
     * IU1, IU2, R0, R1, R2, R3, R4, RA1, RA2, RAS, RAU, RB1, RB2,
     * RS1, RS2, RSS, RU1, RU2, S1, S2, S3, S4, T, TWOPI
      INTEGER J, K0, K, M5, MOVER2
      LOGICAL FOLD, ZERO
C     .. FUNCTION REFERENCES ..
      REAL*8 DCOS, DSIN, DSQRT
C     ..
      M5 = M*5
      MOVER2 = M/2 + 1
      twopi=8.0*datan(1.d00)
      A1 = DCOS(TWOPI/5.0D00)
      B1 = DSIN(TWOPI/5.0D00)
      A2 = DCOS(2.0*TWOPI/5.0D00)
      B2 = DSIN(2.0*TWOPI/5.0D00)
      AS = -1.0/4.0
      AU = DSQRT(5.0D00)/4.0
C
      DO 120 J=1,MOVER2
         FOLD = J.GT.1 .AND. 2*J.LT.M + 2
         K0 = J
         ANGLE = TWOPI*FLOAT(J-1)/FLOAT(M5)
         ZERO = ANGLE.EQ.0.0
         C1 = DCOS(ANGLE)
         S1 = DSIN(ANGLE)
         C2 = C1*C1 - S1*S1
         S2 = S1*C1 + C1*S1
         C3 = C2*C1 - S2*S1
         S3 = S2*C1 + C2*S1
         C4 = C2*C2 - S2*S2
         S4 = S2*C2 + C2*S2
         GO TO 40
   20    CONTINUE
         FOLD = .FALSE.
         K0 = M + 2 - J
         T = C1*A1 + S1*B1
         S1 = C1*B1 - S1*A1
         C1 = T
         T = C2*A2 + S2*B2
         S2 = C2*B2 - S2*A2
         C2 = T
         T = C3*A2 - S3*B2
         S3 = -C3*B2 - S3*A2
         C3 = T
         T = C4*A1 - S4*B1
         S4 = -C4*B1 - S4*A1
         C4 = T
   40    CONTINUE
C
         DO 100 K=K0,PTS,M5
            R0 = X0(K)
            I0 = Y0(K)
            RS1 = X1(K) + X4(K)
            IS1 = Y1(K) + Y4(K)
            RU1 = X1(K) - X4(K)
            IU1 = Y1(K) - Y4(K)
            RS2 = X2(K) + X3(K)
            IS2 = Y2(K) + Y3(K)
            RU2 = X2(K) - X3(K)
            IU2 = Y2(K) - Y3(K)
            RSS = RS1 + RS2
            ISS = IS1 + IS2
            RAS = R0 + RSS*AS
            IAS = I0 + ISS*AS
            RAU = (RS1-RS2)*AU
            IAU = (IS1-IS2)*AU
            RA1 = RAS + RAU
            IA1 = IAS + IAU
            RA2 = RAS - RAU
            IA2 = IAS - IAU
            RB1 = RU1*B1 + RU2*B2
            IB1 = IU1*B1 + IU2*B2
            RB2 = RU1*B2 - RU2*B1
            IB2 = IU1*B2 - IU2*B1
            X0(K) = R0 + RSS
            Y0(K) = I0 + ISS
            IF (ZERO) GO TO 60
            R1 = RA1 + IB1
            I1 = IA1 - RB1
            R2 = RA2 + IB2
            I2 = IA2 - RB2
            R3 = RA2 - IB2
            I3 = IA2 + RB2
            R4 = RA1 - IB1
            I4 = IA1 + RB1
            X1(K) = R1*C1 + I1*S1
            Y1(K) = I1*C1 - R1*S1
            X2(K) = R2*C2 + I2*S2
            Y2(K) = I2*C2 - R2*S2
            X3(K) = R3*C3 + I3*S3
            Y3(K) = I3*C3 - R3*S3
            X4(K) = R4*C4 + I4*S4
            Y4(K) = I4*C4 - R4*S4
            GO TO 80
   60       CONTINUE
            X1(K) = RA1 + IB1
            Y1(K) = IA1 - RB1
            X2(K) = RA2 + IB2
            Y2(K) = IA2 - RB2
            X3(K) = RA2 - IB2
            Y3(K) = IA2 + RB2
            X4(K) = RA1 - IB1
            Y4(K) = IA1 + RB1
   80       CONTINUE
  100    CONTINUE
         IF (FOLD) GO TO 20
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE ECTC06(X0, Y0, PTS, X1, Y1, M1, X2, Y2, M2, X3,
     * Y3, M3, M)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     RADIX FOUR COMPLEX FOURIER TRANSFORM KERNEL
C     .. SCALAR ARGUMENTS ..
      INTEGER M1, M2, M3, M, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X0(PTS), X1(M1), X2(M2), X3(M3), Y0(PTS), Y1(M1),
     * Y2(M2), Y3(M3)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 ANGLE, C1, C2, C3, I1, I2, I3, IS0, IS1, IU0, IU1, R1,
     * R2, R3, RS0, RS1, RU0, RU1, S1, S2, S3, T, TWOPI
      INTEGER J, K0, K, M4, MOVER2
      LOGICAL FOLD, ZERO
C     .. FUNCTION REFERENCES ..
      REAL*8 DCOS, DSIN
C     ..
      M4 = M*4
      MOVER2 = M/2 + 1
      twopi=8.0*datan(1.d00)
C
      DO 120 J=1,MOVER2
         FOLD = J.GT.1 .AND. 2*J.LT.M + 2
         K0 = J
         ANGLE = TWOPI*FLOAT(J-1)/FLOAT(M4)
         ZERO = ANGLE.EQ.0.0
         C1 = DCOS(ANGLE)
         S1 = DSIN(ANGLE)
         C2 = C1*C1 - S1*S1
         S2 = S1*C1 + C1*S1
         C3 = C2*C1 - S2*S1
         S3 = S2*C1 + C2*S1
         GO TO 40
   20    CONTINUE
         FOLD = .FALSE.
         K0 = M + 2 - J
         T = C1
         C1 = S1
         S1 = T
         C2 = -C2
         T = C3
         C3 = -S3
         S3 = -T
   40    CONTINUE
C
         DO 100 K=K0,PTS,M4
            RS0 = X0(K) + X2(K)
            IS0 = Y0(K) + Y2(K)
            RU0 = X0(K) - X2(K)
            IU0 = Y0(K) - Y2(K)
            RS1 = X1(K) + X3(K)
            IS1 = Y1(K) + Y3(K)
            RU1 = X1(K) - X3(K)
            IU1 = Y1(K) - Y3(K)
            X0(K) = RS0 + RS1
            Y0(K) = IS0 + IS1
            IF (ZERO) GO TO 60
            R1 = RU0 + IU1
            I1 = IU0 - RU1
            R2 = RS0 - RS1
            I2 = IS0 - IS1
            R3 = RU0 - IU1
            I3 = IU0 + RU1
            X2(K) = R1*C1 + I1*S1
            Y2(K) = I1*C1 - R1*S1
            X1(K) = R2*C2 + I2*S2
            Y1(K) = I2*C2 - R2*S2
            X3(K) = R3*C3 + I3*S3
            Y3(K) = I3*C3 - R3*S3
            GO TO 80
   60       CONTINUE
            X2(K) = RU0 + IU1
            Y2(K) = IU0 - RU1
            X1(K) = RS0 - RS1
            Y1(K) = IS0 - IS1
            X3(K) = RU0 - IU1
            Y3(K) = IU0 + RU1
   80       CONTINUE
  100    CONTINUE
         IF (FOLD) GO TO 20
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE ECUC06(X0, Y0, PTS, X1, Y1, M1, X2, Y2, M2, M)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     RADIX THREE COMPLEX FOURIER TRANSFORM KERNEL
C     .. SCALAR ARGUMENTS ..
      INTEGER M1, M2, M, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X0(PTS), X1(M1), X2(M2), Y0(PTS), Y1(M1), Y2(M2)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 A, ANGLE, B, C1, C2, I0, I1, I2, IA, IB, IS, R0, R1, R2,
     * RA, RB, RS, S1, S2, T, TWOPI
      INTEGER J, K0, K, M3, MOVER2
      LOGICAL FOLD, ZERO
C     .. FUNCTION REFERENCES ..
      REAL*8 DCOS, DSIN, DSQRT
C     ..
      M3 = M*3
      MOVER2 = M/2 + 1
      twopi=8.0*datan(1.d00)
C     A = DCOS(TWOPI/3.0)
      A = -0.5
C     B = DSIN(TWOPI/3.0)
      B = DSQRT(0.75D00)
C
      DO 120 J=1,MOVER2
         FOLD = J.GT.1 .AND. 2*J.LT.M + 2
         K0 = J
         ANGLE = TWOPI*FLOAT(J-1)/FLOAT(M3)
         ZERO = ANGLE.EQ.0.0
         C1 = DCOS(ANGLE)
         S1 = DSIN(ANGLE)
         C2 = C1*C1 - S1*S1
         S2 = S1*C1 + C1*S1
         GO TO 40
   20    CONTINUE
         FOLD = .FALSE.
         K0 = M + 2 - J
         T = C1*A + S1*B
         S1 = C1*B - S1*A
         C1 = T
         T = C2*A - S2*B
         S2 = -C2*B - S2*A
         C2 = T
   40    CONTINUE
C
         DO 100 K=K0,PTS,M3
            R0 = X0(K)
            I0 = Y0(K)
            RS = X1(K) + X2(K)
            IS = Y1(K) + Y2(K)
            X0(K) = R0 + RS
            Y0(K) = I0 + IS
            RA = R0 + RS*A
            IA = I0 + IS*A
            RB = (X1(K)-X2(K))*B
            IB = (Y1(K)-Y2(K))*B
            IF (ZERO) GO TO 60
            R1 = RA + IB
            I1 = IA - RB
            R2 = RA - IB
            I2 = IA + RB
            X1(K) = R1*C1 + I1*S1
            Y1(K) = I1*C1 - R1*S1
            X2(K) = R2*C2 + I2*S2
            Y2(K) = I2*C2 - R2*S2
            GO TO 80
   60       CONTINUE
            X1(K) = RA + IB
            Y1(K) = IA - RB
            X2(K) = RA - IB
            Y2(K) = IA + RB
   80       CONTINUE
  100    CONTINUE
         IF (FOLD) GO TO 20
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE ECVC06(X0, Y0, PTS, X1, Y1, M1, M)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     RADIX TWO COMPLEX FOURIER TRANSFORM KERNEL
C     .. SCALAR ARGUMENTS ..
      INTEGER M1, M, PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X0(PTS), X1(M1), Y0(PTS), Y1(M1)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 ANGLE, C, IS, IU, RS, RU, S, TWOPI
      INTEGER J, K0, K, M2, MOVER2
      LOGICAL FOLD, ZERO
C     .. FUNCTION REFERENCES ..
      REAL*8 DCOS, DSIN
C     ..
      M2 = M*2
      MOVER2 = M/2 + 1
      twopi=8.0*datan(1.d00)
C
      DO 120 J=1,MOVER2
         FOLD = J.GT.1 .AND. 2*J.LT.M + 2
         K0 = J
         ANGLE = TWOPI*FLOAT(J-1)/FLOAT(M2)
         ZERO = ANGLE.EQ.0.0
         C = DCOS(ANGLE)
         S = DSIN(ANGLE)
         GO TO 40
   20    CONTINUE
         FOLD = .FALSE.
         K0 = M + 2 - J
         C = -C
   40    CONTINUE
C
         DO 100 K=K0,PTS,M2
            RS = X0(K) + X1(K)
            IS = Y0(K) + Y1(K)
            RU = X0(K) - X1(K)
            IU = Y0(K) - Y1(K)
            X0(K) = RS
            Y0(K) = IS
            IF (ZERO) GO TO 60
            X1(K) = RU*C + IU*S
            Y1(K) = IU*C - RU*S
            GO TO 80
   60       CONTINUE
            X1(K) = RU
            Y1(K) = IU
   80       CONTINUE
  100    CONTINUE
         IF (FOLD) GO TO 20
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE ECXC06(X, Y, PTS, SYM)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     DOUBLE SYMMETRIC REORDERING PROGRAMME
C     EQUIVALENCE (I1,I(1)), (K1,K(1)), (L1,L(1))
C     EQUIVALENCE (I2,I(2)), (K2,K(2)), (L2,L(2))
C     EQUIVALENCE (I3,I(3)), (K3,K(3)), (L3,L(3))
C     EQUIVALENCE (I4,I(4)), (K4,K(4)), (L4,L(4))
C     EQUIVALENCE (I5,I(5)), (K5,K(5)), (L5,L(5))
C     EQUIVALENCE (I6,I(6)), (K6,K(6)), (L6,L(6))
C     EQUIVALENCE (I7,I(7)), (K7,K(7)), (L7,L(7))
C     EQUIVALENCE (I8,I(8)), (K8,K(8)), (L8,L(8))
C     EQUIVALENCE (I9,I(9)), (K9,K(9)), (L9,L(9))
C     EQUIVALENCE (I10,I(10)), (K10,K(10)), (L10,L(10))
C     EQUIVALENCE (K11,K(11))
C     .. SCALAR ARGUMENTS ..
      INTEGER PTS
C     .. ARRAY ARGUMENTS ..
      REAL*8 X(PTS), Y(PTS)
      INTEGER SYM(21)
C     ..
C     .. LOCAL SCALARS ..
      REAL*8 T
      INTEGER I10, I1, I2, I3, I4, I5, I6, I7, I8, I9, J, JJ, K10,
     * K11, K1, K2, K3, K4, K5, K6, K7, K8, K9, KK, L10, L1, L2,
     * L3, L4, L5, L6, L7, L8, L9, LEVEL, LOOP, NEST
C     .. LOCAL ARRAYS ..
      INTEGER I(20), K(20), L(20)
C     ..
      DATA NEST /20/
      DATA LOOP /10/
      IF (SYM(1).EQ.0) GO TO 360
      DO 20 J=1,NEST
         L(J) = 1
         I(J) = 1
   20 CONTINUE
      KK = PTS
      DO 40 J=1,NEST
         IF (SYM(J).EQ.0) GO TO 60
         L(J) = KK
         I(J) = KK/SYM(J)
         KK = KK/SYM(J)
   40 CONTINUE
   60 CONTINUE
C
      L1 = L(1)
      L2 = L(2)
      L3 = L(3)
      L4 = L(4)
      L5 = L(5)
      L6 = L(6)
      L7 = L(7)
      L8 = L(8)
      L9 = L(9)
      L10 = L(10)
      I1 = I(1)
      I2 = I(2)
      I3 = I(3)
      I4 = I(4)
      I5 = I(5)
      I6 = I(6)
      I7 = I(7)
      I8 = I(8)
      I9 = I(9)
      I10 = I(10)
C
      KK = 0
      LEVEL = NEST
      K(LEVEL) = 1
      GO TO 100
   80 CONTINUE
      IF (LEVEL.GE.NEST) GO TO 360
      LEVEL = LEVEL + 1
      K(LEVEL) = K(LEVEL) + I(LEVEL)
      IF (K(LEVEL).GT.L(LEVEL)) GO TO 80
  100 CONTINUE
      LEVEL = LEVEL - 1
      DO 120 J=LOOP,LEVEL
         JJ = LEVEL + LOOP - J
         K(JJ) = K(JJ+1)
  120 CONTINUE
      K11 = K(11)
      DO 340 K10=K11,L10,I10
         K(10) = K10
         DO 320 K9=K10,L9,I9
            K(9) = K9
            DO 300 K8=K9,L8,I8
               K(8) = K8
               DO 280 K7=K8,L7,I7
                  K(7) = K7
                  DO 260 K6=K7,L6,I6
                     K(6) = K6
                     DO 240 K5=K6,L5,I5
                        K(5) = K5
                        DO 220 K4=K5,L4,I4
                           K(4) = K4
                           DO 200 K3=K4,L3,I3
                              K(3) = K3
                              DO 180 K2=K3,L2,I2
                                 K(2) = K2
                                 DO 160 K1=K2,L1,I1
                                    K(1) = K1
                                    KK = KK + 1
                                    IF (KK.GE.K1) GO TO 140
                                    T = X(KK)
                                    X(KK) = X(K1)
                                    X(K1) = T
                                    T = Y(KK)
                                    Y(KK) = Y(K1)
                                    Y(K1) = T
  140                               CONTINUE
  160                            CONTINUE
  180                         CONTINUE
  200                      CONTINUE
  220                   CONTINUE
  240                CONTINUE
  260             CONTINUE
  280          CONTINUE
  300       CONTINUE
  320    CONTINUE
  340 CONTINUE
      LEVEL = LOOP
      GO TO 80
C
  360 CONTINUE
      RETURN
      END
      SUBROUTINE EAZC06(PTS, PMAX, TWOGRP, FACTOR, SYM, PSYM,
     * UNSYM, IERROR)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     SYMMETRIZED REORDERING FACTORING PROGRAMME
C     .. SCALAR ARGUMENTS ..
      INTEGER IERROR, PMAX, PSYM, PTS, TWOGRP
C     .. ARRAY ARGUMENTS ..
      INTEGER FACTOR(21), SYM(21), UNSYM(21)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER F, J, JJ, N, NEST, P, PTWO, Q, R
C     .. LOCAL ARRAYS ..
      INTEGER PP(10), QQ(20)
C     ..
      DATA NEST /20/
      N = PTS
      PSYM = 1
      F = 2
      P = 0
      Q = 0
   20 CONTINUE
      IF (N.LE.1) GO TO 100
      DO 40 J=F,PMAX
         IF (N.EQ.(N/J)*J) GO TO 60
   40 CONTINUE
      GO TO 280
   60 CONTINUE
      IF (2*P+Q.GE.NEST) GO TO 300
      F = J
      N = N/F
      IF (N.EQ.(N/F)*F) GO TO 80
      Q = Q + 1
      QQ(Q) = F
      GO TO 20
   80 CONTINUE
      N = N/F
      P = P + 1
      PP(P) = F
      PSYM = PSYM*F
      GO TO 20
C
  100 CONTINUE
      R = 1
      IF (Q.EQ.0) R = 0
      IF (P.LT.1) GO TO 140
      DO 120 J=1,P
         JJ = P + 1 - J
         SYM(J) = PP(JJ)
         FACTOR(J) = PP(JJ)
         JJ = P + Q + J
         FACTOR(JJ) = PP(J)
         JJ = P + R + J
         SYM(JJ) = PP(J)
  120 CONTINUE
  140 CONTINUE
      IF (Q.LT.1) GO TO 180
      DO 160 J=1,Q
         JJ = P + J
         UNSYM(J) = QQ(J)
         FACTOR(JJ) = QQ(J)
  160 CONTINUE
      SYM(P+1) = PTS/PSYM**2
  180 CONTINUE
      JJ = 2*P + Q
      FACTOR(JJ+1) = 0
      PTWO = 1
      J = 0
  200 CONTINUE
      J = J + 1
      IF (FACTOR(J).EQ.0) GO TO 240
      IF (FACTOR(J).NE.2) GO TO 200
      PTWO = PTWO*2
      FACTOR(J) = 1
      IF (PTWO.GE.TWOGRP) GO TO 220
      IF (FACTOR(J+1).EQ.2) GO TO 200
  220 CONTINUE
      FACTOR(J) = PTWO
      PTWO = 1
      GO TO 200
  240 CONTINUE
      IF (P.EQ.0) R = 0
      JJ = 2*P + R
      SYM(JJ+1) = 0
      IF (Q.LE.1) Q = 0
      UNSYM(Q+1) = 0
      IERROR = 0
C
  260 CONTINUE
      RETURN
C
  280 CONTINUE
      IERROR = 1
      GO TO 260
C
  300 CONTINUE
      IERROR = 2
      GO TO 260
C
      END
