      SUBROUTINE ALPHA(CALPHA,CXALPH,CXE0R,CXEPS,ICOMP,AKR,NAT0,NAT3,
     &                 MXCOMP,MXN3)
C*** Arguments:
      CHARACTER*6 CALPHA
      INTEGER MXN3,MXCOMP,NAT0,NAT3
      COMPLEX CXALPH(MXN3),CXE0R(3),CXEPS(MXCOMP)
      INTEGER*2 ICOMP(MXN3)
      REAL AKR(3)
C
C*** Local Variables:
      CHARACTER CMSGNM*70
      INTEGER IC,L
      COMPLEX CXI,CXRR
      REAL AK1,AK2,B1,B2,EMKD,PI,SUM
      DATA CXI/(0.,1.)/
C***********************************************************************
C Given:
C       CALPHA = 'LATTDR' or 'DRAI88' or 'GOBR88' or 'LDRISO'
C                to specify method ('LATTDR' is recommended)
C       CXE0R(1-3) = polarization vector in lattice coordinates
C                    (assumed to be normalized)
C       CXEPS(1-NCOMP)=distinct values of dielectric constant
C       ICOMP(1-NAT3)=composition identifier for each lattice site and
C                     direction (ICOMP=0 if lattice site is unoccupied)
C       MXCOMP = dimensioning information
C       MXN3 = dimensioning information
C       NAT0 = number of sites in original target (i.e., nonvacuum sites)
C       NAT3 = 3*number of sites in extended target (including vacuum
C              sites)
C Returns:
C       CXALPHA(1-NAT3)=(polarizability/d^3) for each of 3 directions at
C                       each lattice site
C***
C If CALPHA = LATTDR:
C    Compute dipole polarizability using "Lattice Dispersion Relation"
C    of Draine & Goodman (1993,ApJ,March 10).  It is required that
C    polarizability be such that an infinite lattice of such dipoles
C    reproduce the continuum dispersion relation for radiation
C    propagating with direction and polarization of radiation incident
C    on the DDA target.  This is the recommended option.
C
C If CALPHA = DRAI88:
C    Compute dipole polarizability using
C           (1) eq. (2.03) of Draine (1988) (Clausius-Mossotti plus
C               radiative reaction correction)
C
C If CALPHA = GOBR88:
C    Compute dipole polarizability using
C           (1) radiative reaction correction from Draine (1988)
C           (2) O[(kd)^2] correction term used by
C               Goedecke & OBrien (1988) and Hage & Greenberg (1990)
C
C If CALPHA = LDRISO:
C    Compute dipole polarizability using "Isotropized Lattice Dispersion
C    Relation" of Draine & Goodman (1993,ApJ,March 10).  This replaces
C    direction- and polarization-dependent term in Lattice Dispersion
C    Relation with the average value for random direction and polarization.
C
C***
C    Note: CXALPH = polarizability/d^3
C    In the event that ICOMP=0, then we set CXALPH=1.
C***********************************************************************
C B.T.Draine, Princeton Univ. Obs.
C History:
C 90. 9.13 (BTD): Corrected error in DO loop limit.
C 90.11. 1 (BTD): Special treatment for "vacuum" sites in order
C                 to allow FFT treatment.
C 90.11. 2 (BTD): Set CXALPH=(1.,0.) at vacuum sites.
C 91. 4.30 (BTD): Modified to include both O[(kd)^2] correction term
C                 from Goedecke & Obrien (1988) in addition to
C                 radiative reaction correction.
C 91. 5. 7 (BTD): Modified to allow easy choice among three methods
C                 for computing polarizability.
C 91. 5. 7 (BTD): Added call to WRIMSG to record which method in use
C 91.05.08 (BTD): Added CALPHA to argument list.
C 91.05.09 (BTD): Added printing of kd and magnitude of correction
C                 terms.
C 91.05.13 (BTD): Experiment with use of approximate directional
C                 average for coefficient of CXEPS in Lattice
C                 Disperision Relation theory
C 91.05.13 (BTD): Corrected error in radiative reaction correction
C                 (error presumably introduced during last week)
C 91.05.14 (BTD): Added separate options LDRXYZ and LDRAVG
C 91.05.15 (BTD): Added option LDR000 to omit CXEPS-dependent
C                 correction term
C 91.05.30 (BTD): Changed LDRXYZ -> LDR100
C                 Added option LDR111
C 91.09.12 (BTD): Corrected error in numerical coefficient for
C                 LDRAVG case
C 91.09.17 (BTD): Eliminate options LDR000,LDRXYZ,LDRAVG
C                 Introduce option LATTDR (LATTice Dispersion Relation)
C                 which includes explicit dependence
C                 on direction and polarization state
C                 remove AK1 from argument list
C                 add AKR, CXE0R to argument list
C 92.05.14 (BTD): Introduce option LDRISO to experiment with
C                 using average directional correction term
C                 rather than using correction for incident direction
C                 and polarization
C 92.05.16 (BTD): Correct error in LDRISO option: average SUM should be
C                 3/15 rather than 4/15
C
C Copyright (C) 1993,1996 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      PI=4.*ATAN(1.)
      AK2=AKR(1)*AKR(1)+AKR(2)*AKR(2)+AKR(3)*AKR(3)
      AK1=SQRT(AK2)
      CXRR=-(AK1*AK2/1.5)*CXI
C*** EMKD = |m|*k_0*d , where |m|=refractive index
C                             k_0 = wave vector in vacuo
C                             d = lattice spacing
      EMKD=SQRT(CABS(CXEPS(1)))*AK1
C
C***********************************************************************
C*** Draine 1988:
C
      IF(CALPHA.EQ.'DRAI88')THEN
         WRITE(CMSGNM,FMT='(A,F6.4)')
     &         'Clausius-Mossotti plus radiative reaction for |m|kd=',
     &         EMKD
         CALL WRIMSG('ALPHA ',CMSGNM)
         DO 1000 L=1,NAT3
            IC=ICOMP(L)
            IF(IC.GT.0)THEN
C*** First compute Clausius-Mossotti polarizability:
               CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.)
C*** Now correct this for
C    (1) radiative-reaction
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR)
            ELSEIF(IC.EQ.0)THEN
C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites.
               CXALPH(L)=1.
            ENDIF
 1000    CONTINUE
         RETURN
C
C***********************************************************************
C*** Goedecke & OBrien 1988:
C
      ELSEIF(CALPHA.EQ.'GOBR88')THEN
         B1=-(4.*PI/3.)**(1./3.)*AK2
         WRITE(CMSGNM,FMT='(A,F6.4)')
     &         'Goedecke & OBrien (1988) correction for |m|kd=',EMKD
         CALL WRIMSG('ALPHA ',CMSGNM)
         DO 2000 L=1,NAT3
            IC=ICOMP(L)
            IF(IC.GT.0)THEN
C*** First compute Clausius-Mossotti polarizability:
               CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.)
C*** Now correct this for
C    (1) O[(kd)**2] correction from Goedecke and Obrien (1988)
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*B1)
C    (2) radiative-reaction
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR)
            ELSEIF(IC.EQ.0)THEN
C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites.
               CXALPH(L)=1.
            ENDIF
 2000    CONTINUE
         RETURN
C
C***********************************************************************
C*** Lattice dispersion relation
C
      ELSEIF(CALPHA.EQ.'LATTDR')THEN
         WRITE(CMSGNM,FMT='(A,F6.4)')
     &         'Lattice dispersion relation for |m|k_0d=',EMKD
         CALL WRIMSG('ALPHA ',CMSGNM)
C*** Compute sum (a_j*e_j)^2 , where a_j=unit propagation vector
C                                    e_j=unit polarization vector
         SUM=0.
         DO 2700 L=1,3
            SUM=SUM+(AKR(L)*CABS(CXE0R(L)))**2
 2700    CONTINUE
         SUM=SUM/AK2
         B1=-1.8915316*AK2
         B2=(0.1648469-1.7700004*SUM)*AK2
         DO 3000 L=1,NAT3
            IC=ICOMP(L)
            IF(IC.GT.0)THEN
C*** First compute Clausius-Mossotti polarizability:
               CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.)
C*** Determine polarizability by requiring that infinite lattice of
C    dipoles have dipersion relation of continuum.
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*(B1+CXEPS(IC)*B2))
C*** Radiative-reaction correction:
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR)
            ELSEIF(IC.EQ.0)THEN
C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites.
               CXALPH(L)=1.
            ENDIF
 3000    CONTINUE
         RETURN
      ELSEIF(CALPHA.EQ.'LDRISO')THEN
         WRITE(CMSGNM,FMT='(A,F6.4)')
     &         'Isotropic Lattice disp. rel. for |m|k_0d=',EMKD
         CALL WRIMSG('ALPHA ',CMSGNM)
C*** Compute sum (a_j*e_j)^2 , where a_j=unit propagation vector
C                                    e_j=unit polarization vector
C    Here experiment by taking average value of sum (a_j*e_j)^2 = 1/5
         SUM=0.2D0
         B1=-1.8915316*AK2
         B2=(0.1648469-1.7700004*SUM)*AK2
         DO 4000 L=1,NAT3
            IC=ICOMP(L)
            IF(IC.GT.0)THEN
C*** First compute Clausius-Mossotti polarizability:
               CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.)
C*** Determine polarizability by requiring that infinite lattice of
C    dipoles have dipersion relation of continuum.
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*(B1+CXEPS(IC)*B2))
C*** Radiative-reaction correction:
               CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR)
            ELSEIF(IC.EQ.0)THEN
C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites.
               CXALPH(L)=1.
            ENDIF
 4000    CONTINUE
         RETURN
      ENDIF
      END

      SUBROUTINE PETR(X,B,WRK,LDA,IPAR,SPAR,MATVEC,CMATVEC)
C solves Ax=b
C
C Input:
C B(LOCLEN) - RHS (doesn't change on output)
C WRK(LDA, LOCLEN) - scratch array 
C IPAR - defines integer paramteres (see documentation of PIM)
C SPAR - define real parameters
C LDA - leading dimension of arrays 
C MATVEC -- name of procedure for computing y=Ax
C CMATVEC -- name of procedure for computing y= conj(A') x
C
C Output:
C X(LOCLEN)
C 
C Reference:
C    Petravic,M and Kuo-Petravic,G. 1979, JCP, 32,263
C History: 
C    Late 80'  (Coded by Draine)
C    Changes to include "cprod" (PJF+BTD)
C    95.06.01 (PJF) re-written to conform to PIM
C    95.08.11 (BTD) minor editing (comments, etc.)
C
C Parameters:
      INTEGER IACE,IAXI,IGI,IPI,IQI,IR
      PARAMETER (IACE=1,IGI=2,IPI=3,IQI=4,IAXI=5,IR=6)
C Arguments:
      INTEGER LDA
      INTEGER IPAR(*)
      REAL SPAR(*)
      COMPLEX B(*), WRK(LDA,*), X(*)
      EXTERNAL CMATVEC,MATVEC
C
C Local variables
      COMPLEX CDUMMY
      INTEGER IDUMMY,ITNO,L,LOCLEN,MAXIT,STATUS
      REAL
     &   ALPHAI,BETAI1,BNORM,DUMMY,EPSILON,EXITNORM,
     &   GIGI,GI1GI1,QIQI,RHSSTOP
      REAL SCSETRHSSTOP
      EXTERNAL PSCNRM2
C
C***********************************************************************
      LOCLEN = IPAR(4)
      MAXIT = IPAR(10)
      EPSILON=SPAR(1)
C
C Compute RHSSTOP = EPSILON*|B|
C
      RHSSTOP = SCSETRHSSTOP(B,CDUMMY,EPSILON,IPAR,CDUMMY,PSCNRM2)
C
C Compute conjg(A')*B
C
      CALL CMATVEC(B,WRK(1,IACE),IDUMMY)
      BNORM = 0.
      DO 170 L = 1,LOCLEN
         BNORM = BNORM + B(L)*CONJG(B(L))
         WRK(L,IGI) = WRK(L,IACE)
         WRK(L,IPI) = WRK(L,IGI)
  170 CONTINUE
C
C Compute |QI>=A|PI>
C
      CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY)
      QIQI = 0.
      GIGI = 0.
      DO 190 L = 1,LOCLEN
         GIGI = GIGI + CONJG(WRK(L,IGI))*WRK(l,IGI)
         QIQI = QIQI + CONJG(WRK(L,IQI))*WRK(L,IQI)
  190 CONTINUE
      ALPHAI = GIGI/QIQI
      DO 200 L = 1,LOCLEN
         X(L) = X(L) + ALPHAI*WRK(L,IPI)
  200 CONTINUE
C
C compute |AX1>:
C
      CALL MATVEC(X,WRK(1,IAXI),IDUMMY)
C
      DO 340 ITNO = 2,MAXIT
         IPAR(11) = ITNO
C
C Transfer <GI|GI> -> <GI-1|GI-1>:
C
         GI1GI1 = GIGI
C
C compute |GI>=AC|E>-AC|AXI>:
C
         CALL CMATVEC(WRK(1,IAXI), WRK(1,IGI), IDUMMY)
C
C Compute GIGI=<GI|GI>:
         GIGI = 0.
         DO 230 L = 1,LOCLEN
            WRK(L,IGI) = WRK(L,IACE) - WRK(L,IGI)
  230    CONTINUE
         DO 235 L = 1,LOCLEN
            GIGI = GIGI + CONJG(WRK(L,IGI))*WRK(L,IGI)
  235    CONTINUE
C
C Compute BETAI-1=<GI|GI>/<GI-1|GI-1>:
C
         BETAI1 = GIGI/GI1GI1
C
C Compute |PI>=|GI>+BETAI-1|PI-1>:
C
         DO 240 L = 1,LOCLEN
            WRK(L,IPI) = WRK(L,IGI) + BETAI1*WRK(L,IPI)
  240    CONTINUE
C
C Compute |QI>=A|PI>:
C
         CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY)
C
C Compute <QI|QI>:
C
         QIQI = 0.
         DO 250 L = 1,LOCLEN
            QIQI = QIQI + CONJG(WRK(L,IQI))*WRK(L,IQI)
  250    CONTINUE
C
C Compute ALPHAI = <GI|GI>/<QI|QI>:
C
         ALPHAI = GIGI/QIQI
C
C Compute |XI+1>=|XI>+ALPHAI*|PI>:
C
         DO 260 L = 1,LOCLEN
            X(L) = X(L) + ALPHAI*WRK(L,IPI)
  260    CONTINUE
C
C Except every 10TH iteration, compute |AXI+1>=|AXI>+ALPHAI*|QI>
C Draine: warning this part perhaps not checked 
C
         IF (ITNO.NE.10* (ITNO/10)) THEN
         DO 270 L = 1,LOCLEN
            WRK(L,IAXI) = WRK(L,IAXI) + ALPHAI*WRK(L,IQI)
  270    CONTINUE
         ELSE 
            CALL MATVEC(X,WRK(1,IAXI),IDUMMY)
         ENDIF
C
C Compute residual vector |RI>=A|XI>-|B>
C
         DO 300 L=1, LOCLEN
            WRK(L,IR) = WRK(L,IAXI)-B(L)
 300     CONTINUE
C
C Call STOPCRIT to check whether to terminate iteration.
C Criterion: terminate if <RI|RI> < TOL**2 * <B|B>
C Return with STATUS=0 if this condition is satisfied.
C
         CALL STOPCRIT(B,WRK(1,IR),CDUMMY,CDUMMY,CDUMMY,WRK,
     +                 RHSSTOP,IDUMMY,EXITNORM,
     +                 STATUS,IPAR,DUMMY,DUMMY,DUMMY,DUMMY,PSCNRM2)
C
C Call PROGRESS to report "error" = sqrt(<RI|RI>)/sqrt(<B|B>)
C
         CALL PROGRESS(LOCLEN,ITNO,EXITNORM,CDUMMY,CDUMMY,CDUMMY)
C      
      IF(STATUS.EQ.0)GOTO 500
 340  CONTINUE
 500  CONTINUE
      IPAR(12) = 0
      RETURN
      END



      SUBROUTINE COPYIT(CXA,CXB,NPTS)
C*** Arguments:
      INTEGER NPTS
      COMPLEX CXA(*), CXB(*)
C
C*** Local variables:
      INTEGER IPTS
C***********************************************************************
C Purpose: To copy a complex vector from one memory location to another.
C***********************************************************************
      DO 10 IPTS=1,NPTS
         CXB(IPTS)=CXA(IPTS)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE DIAGL(U,V,IPAR)
C cg package
C Left preconditioning subroutine - division by diagonal
C Arguments:
      COMPLEX U(*), V(*)
      INTEGER IPAR(*)
C
C Common variables:
      CHARACTER*6 CMETHD
      INTEGER
     &   IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &   NAT,NAT0,NAT3,NX,NY,NZ
      COMPLEX CXADIA(1)
C-----------------------------------------------------------------------
      COMMON /M2/ CXADIA
      COMMON /M6/ MXNATF, MXNXF, MXNYF, MXNZF, NAT, NAT3, NAT0, 
     $  NX, NY, NZ, MXMEMF, MXN3F, IDVOUT, CMETHD
C-----------------------------------------------------------------------
      CALL CDIV(U,V,CXADIA,NAT3)
      RETURN
      END
      SUBROUTINE CDIV(U,V,CXA,N)
C Arguments:
      INTEGER N
      COMPLEX U(*), V(*), CXA(*)
C
C Local variables:
      INTEGER I
C
C cg package 
C needed by left preconditioning subroutine
      DO 10 I=1,N
        V(I)=U(I)/CXA(I)
 10   CONTINUE
      RETURN
      END
      SUBROUTINE MATVEC(CXX, CXY, IPAR)
C Arguments:
      COMPLEX CXX(*), CXY(*), IPAR(*)
C Common:
      INTEGER*2 iocc(1)
      INTEGER
     &   IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &   NAT,NAT0,NAT3,NX,NY,NZ
      CHARACTER*6 CMETHD
      REAL AKR(3)
      COMPLEX CXADIA(1),CXZC(1), CXZW(1)
C-----------------------------------------------------------------------
      COMMON /M1/ AKR, /M2/ CXADIA, /M3/ CXZC, /M4/CXZW, /M5/ IOCC
      COMMON /M6/ MXNATF, MXNXF, MXNYF, MXNZF, NAT, NAT3, NAT0, 
     &  NX, NY, NZ, MXMEMF, MXN3F, IDVOUT, CMETHD
C-----------------------------------------------------------------------
C
C***********************************************************************
C matvec  --> cxy = A cxx         'N'
C tmatvec --> cxy = A' cxx        'T'
C cmatvec --> cxy = conjg(A') cxx  'C'
C 
C Notice that DDSCAT has cases 'N' and 'C' implemented in CPROD
C Notice, however, that in DDSCAT 'T'='N'.
C Some CG methods use cases N, C (like Petravic).
C Some CG methods use cases N, T (like PIM)
C Some CG use all three products N, C, T  (parts of CCGPAK)
C
C Information about matrix A is transfered via commons /m1/-/m6/ with
C the main program. Arrays cxadia, cxzc, cxzwm, iocc are properly
C dimensioned inside "cprod". You may get warning during the compilation
C about "misaligned common". Ignore these. Also notice mxnat-f,
C mxn3-f,mxnx-f, mxny-f, mxnz-f. These are defined by "parameter
C statement" in main code and this is the reason why they have "f" suffix.
C
C ipar - information array, not used, defined for compatibility
C***********************************************************************
      CALL CPROD(AKR,CMETHD,'N',CXADIA,
     &           CXX,CXY,CXZC,CXZW,
     &           IDVOUT,IOCC,MXN3f,MXNATf,MXNXf,MXNYf,MXNZf,
     &           NAT,NAT0,NAT3,NX,NY,NZ)
      RETURN
      END

      SUBROUTINE CMATVEC(CXX, CXY, IPAR)
C Arguments:
      COMPLEX CXX(*), CXY(*), IPAR(*)
C
C Common:
      INTEGER*2 iocc(1)
      INTEGER
     &   IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &   NAT,NAT0,NAT3,NX,NY,NZ
      CHARACTER*6 CMETHD
      REAL akr(3)
      COMPLEX CXADIA(1),CXZC(1), CXZW(1)
C-----------------------------------------------------------------------
      COMMON /M1/ AKR, /M2/ CXADIA, /M3/ CXZC, /M4/CXZW, /M5/ IOCC
      COMMON /M6/ MXNATF, MXNXF, MXNYF, MXNZF, NAT, NAT3, NAT0, 
     &  NX, NY, NZ, MXMEMF, MXN3F, IDVOUT, CMETHD
C-----------------------------------------------------------------------
C Local variables:
C***********************************************************************
C matvec  --> cxy = A cxx   
C tmatvec --> cxy = A' cxx   
C cmatvec --> cxy = conjg(A') cxx   
C***********************************************************************
      CALL CPROD(AKR,CMETHD,'C',CXADIA,
     &           CXX,CXY,CXZC,CXZW,
     &           IDVOUT,IOCC,MXN3f,MXNATf,MXNXf,MXNYf,MXNZf,
     &           NAT,NAT0,NAT3,NX,NY,NZ)
      RETURN
      END

      SUBROUTINE CPROD(AK,CMETHD,CWHAT,CXADIA,
     &                 CXX,CXY,CXZC,CXZW,
     &                 IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ,
     &                 NAT,NAT0,NAT3,NX,NY,NZ)
C
C If CWHAT = 'N', this returns CXY(j)=A(j,k)*CXX(k) (sum over k)
C          = 'C', this returns CXY(j)=conjg(A(k,j))*CXX(k) (sum over k)
C
C It is assumed that matrix A(j,k) is symmetric
C Diagonal elements of A(j,k) are in vector ADIAG(j)
C Off-diagonal elements of A(j,k) are produced internally by
C   subroutine ESELF, which uses 3d-FFT to accomplish convolution
C   of CXX with A-matrix.
C
C Original version of CPROD created by P.J.Flatau, Colorado State Univ.
C Subsequently modified by B.T.Draine, Princeton Univ. Obs.
C History:
C 90/11/ 4 (BTD): Major modification to incorporate FFT method (FFT
C                 computations are handled by ESELF).
C 90/11/ 8 (BTD): Remove DO...ENDDO structures; renumber DO...CONTINUEs
C 90/11/21 (BTD): Remove "CCGMVM" option: only FFT stuff retained.
C 90/12/03 (BTD): No longer necessary to change ordering of CXX and CXY
C 90/12/05 (BTD): Changed call to ESELF to eliminate MXNX,MXNY,MXNZ
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C*** Arguments:
      INTEGER IDVOUT,MXN3,MXNAT,MXNX,MXNY,MXNZ,
     &        NAT,NAT0,NAT3,NX,NY,NZ
      CHARACTER CMETHD*6,CWHAT*1
      INTEGER*2 IOCC(MXNAT)
      REAL AK(3)
      COMPLEX CXADIA(MXN3),
     &        CXX(MXN3),CXY(MXN3),CXZC(NX+1,NY+1,NZ+1,6),
     &        CXZW(MXNX*MXNY*MXNZ,*)
C
C*** Local variables:
      INTEGER J
      REAL AKD
C
C Intrinsic functions:
      INTRINSIC CONJG
C
C***********************************************************************
      AKD=SQRT(AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3))
C We assume that input vector CXX has zeroes for elements corresponding
C to vacuum sites.
      IF(CWHAT.EQ.'N')THEN
C**********************************************************************
C
          CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,CXZC,CXZW,CXY)
C
C*************************************************************************
C ESELF omitted contribution from diagonal terms: add these.
C Also remember that ESELF computes -A_jk*x_k
          DO 160 J=1,NAT3
              CXY(J)=CXADIA(J)*CXX(J)-CXY(J)
  160     CONTINUE
C
C***
C
      ELSEIF(CWHAT.EQ.'C')THEN
C Need to temporarily replace CXX by conjg(CXX)
          DO 170 J=1,NAT3
              CXX(J)=CONJG(CXX(J))
  170     CONTINUE
C**********************************************************************
C
          CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,CXZC,CXZW,CXY)
C
C************************************************************************
C
C ESELF omitted contribution from diagonal terms: add these.
C Also remember that ESELF computed -A_jk*conjg(x_k), while
C we seek to compute conjg(A_kj)* x_k.  Since A_kj=A_jk, we have
C conjg(A_kj)*x_k=conjg(A_jk*conjg(x_k))
          DO 230 J=1,NAT3
              CXY(J)=CXADIA(J)*CXX(J)-CXY(J)
              CXY(J)=CONJG(CXY(J))
  230     CONTINUE
C Restore CXX:
          DO 240 J=1,NAT3
              CXX(J)=CONJG(CXX(J))
  240     CONTINUE
      ELSE
          WRITE(IDVOUT,FMT=*)' Error in CPROD: CWHAT= ',CWHAT
      ENDIF
C
C*** Now must "clean" product vector CXY:
C    (zero out elements corresponding to vacuum sites)
C
      IF(NAT0.LT.NAT)CALL NULLER(CXY,IOCC,MXNAT,MXN3,NAT)
      RETURN
      END
      SUBROUTINE CFFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C
C PURPOSE      PERFORMS MULTIPLE FAST FOURIER TRANSFORMS.  THIS PACKAGE
C              WILL PERFORM A NUMBER OF SIMULTANEOUS COMPLEX PERIODIC
C              FOURIER TRANSFORMS OR CORRESPONDING INVERSE TRANSFORMS.
C              THAT IS, GIVEN A SET OF COMPLEX GRIDPOINT VECTORS, THE
C              PACKAGE RETURNS A SET OF COMPLEX FOURIER
C              COEFFICIENT VECTORS, OR VICE VERSA.  THE LENGTH OF THE
C              TRANSFORMS MUST BE A NUMBER GREATER THAN 1 THAT HAS
C              NO PRIME FACTORS OTHER THAN 2, 3, AND 5.
C
C              THE PACKAGE CFFT99 CONTAINS SEVERAL USER-LEVEL ROUTINES:
C
C            SUBROUTINE CFTFAX
C                AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE
C                BEFORE A SEQUENCE OF CALLS TO CFFT99
C                (PROVIDED THAT N IS NOT CHANGED).
C
C            SUBROUTINE CFFT99
C                THE ACTUAL TRANSFORM ROUTINE ROUTINE, CABABLE OF
C                PERFORMING BOTH THE TRANSFORM AND ITS INVERSE.
C                HOWEVER, AS THE TRANSFORMS ARE NOT NORMALIZED,
C                THE APPLICATION OF A TRANSFORM FOLLOWED BY ITS
C                INVERSE WILL YIELD THE ORIGINAL VALUES MULTIPLIED
C                BY N.
C
C
C ACCESS       *FORTRAN,P=XLIB,SN=CFFT99
C
C
C USAGE        LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 0,
C              Q .GE. 0, AND R .GE. 0.  THEN A TYPICAL SEQUENCE OF
C              CALLS TO TRANSFORM A GIVEN SET OF COMPLEX VECTORS OF
C              LENGTH N TO A SET OF (UNSCALED) COMPLEX FOURIER
C              COEFFICIENT VECTORS OF LENGTH N IS
C
C                   DIMENSION IFAX(13),TRIGS(2*N)
C                   COMPLEX A(...), WORK(...)
C
C                   CALL CFTFAX (N, IFAX, TRIGS)
C                   CALL CFFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C
C              THE OUTPUT VECTORS OVERWRITE THE INPUT VECTORS, AND
C              THESE ARE STORED IN A.  WITH APPROPRIATE CHOICES FOR
C              THE OTHER ARGUMENTS, THESE VECTORS MAY BE CONSIDERED
C              EITHER THE ROWS OR THE COLUMNS OF THE ARRAY A.
C              SEE THE INDIVIDUAL WRITE-UPS FOR CFTFAX AND
C              CFFT99 BELOW, FOR A DETAILED DESCRIPTION OF THE
C              ARGUMENTS.
C
C HISTORY      THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN
C              NOVEMBER, 1978.  IT WAS MODIFIED, DOCUMENTED, AND TESTED
C              FOR NCAR BY RUSS REW IN SEPTEMBER, 1980.  IT WAS
C              FURTHER MODIFIED FOR THE FULLY COMPLEX CASE BY DAVE
C              FULKER IN NOVEMBER, 1980.
C
C-----------------------------------------------------------------------
C
C SUBROUTINE CFTFAX (N,IFAX,TRIGS)
C
C PURPOSE      A SET-UP ROUTINE FOR CFFT99.  IT NEED ONLY BE
C              CALLED ONCE BEFORE A SEQUENCE OF CALLS TO CFFT99,
C              PROVIDED THAT N IS NOT CHANGED.
C
C ARGUMENT     IFAX(13),TRIGS(2*N)
C DIMENSIONS
C
C ARGUMENTS
C
C ON INPUT     N
C               AN EVEN NUMBER GREATER THAN 1 THAT HAS NO PRIME FACTOR
C               GREATER THAN 5.  N IS THE LENGTH OF THE TRANSFORMS (SEE
C               THE DOCUMENTATION FOR CFFT99 FOR THE DEFINITION OF
C               THE TRANSFORMS).
C
C              IFAX
C               AN INTEGER ARRAY.  THE NUMBER OF ELEMENTS ACTUALLY USED
C               WILL DEPEND ON THE FACTORIZATION OF N.  DIMENSIONING
C               IFAX FOR 13 SUFFICES FOR ALL N LESS THAN 1 MILLION.
C
C              TRIGS
C               A REAL ARRAY OF DIMENSION 2*N
C
C ON OUTPUT    IFAX
C               CONTAINS THE FACTORIZATION OF N.  IFAX(1) IS THE
C               NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED
C               IN IFAX(2),IFAX(3),...  IF N HAS ANY PRIME FACTORS
C               GREATER THAN 5, IFAX(1) IS SET TO -99.
C
C              TRIGS
C               AN ARRAY OF TRIGONOMETRIC FUNCTION VALUES SUBSEQUENTLY
C               USED BY THE CFT ROUTINES.
C
C-----------------------------------------------------------------------
C
C SUBROUTINE CFFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C
C PURPOSE      PERFORM A NUMBER OF SIMULTANEOUS (UNNORMALIZED) COMPLEX
C              PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
C              TRANSFORMS.  GIVEN A SET OF COMPLEX GRIDPOINT
C              VECTORS, THE PACKAGE RETURNS A SET OF
C              COMPLEX FOURIER COEFFICIENT VECTORS, OR VICE
C              VERSA.  THE LENGTH OF THE TRANSFORMS MUST BE A
C              NUMBER HAVING NO PRIME FACTORS OTHER THAN
C              2, 3, AND 5.  THIS ROUTINE IS
C              OPTIMIZED FOR USE ON THE CRAY-1.
C
C ARGUMENT     COMPLEX A(N*INC+(LOT-1)*JUMP), WORK(N*LOT)
C DIMENSIONS   REAL TRIGS(2*N), INTEGER IFAX(13)
C
C ARGUMENTS
C
C ON INPUT     A
C               A COMPLEX ARRAY OF LENGTH N*INC+(LOT-1)*JUMP CONTAINING
C               THE INPUT GRIDPOINT OR COEFFICIENT VECTORS.  THIS ARRAY
C               OVERWRITTEN BY THE RESULTS.
C
C               N.B. ALTHOUGH THE ARRAY A IS USUALLY CONSIDERED TO BE OF
C               TYPE COMPLEX IN THE CALLING PROGRAM, IT IS TREATED AS
C               REAL WITHIN THE TRANSFORM PACKAGE.  THIS REQUIRES THAT
C               SUCH TYPE CONFLICTS ARE PERMITTED IN THE USER"S
C               ENVIRONMENT, AND THAT THE STORAGE OF COMPLEX NUMBERS
C               MATCHES THE ASSUMPTIONS OF THIS ROUTINE.  THIS ROUTINE
C               ASSUMES THAT THE REAL AND IMAGINARY PORTIONS OF A
C               COMPLEX NUMBER OCCUPY ADJACENT ELEMENTS OF MEMORY.  IF
C               THESE CONDITIONS ARE NOT MET, THE USER MUST TREAT THE
C               ARRAY A AS REAL (AND OF TWICE THE ABOVE LENGTH), AND
C               WRITE THE CALLING PROGRAM TO TREAT THE REAL AND
C               IMAGINARY PORTIONS EXPLICITLY.
C
C              WORK
C               A COMPLEX WORK ARRAY OF LENGTH N*LOT OR A REAL ARRAY
C               OF LENGTH 2*N*LOT.  SEE N.B. ABOVE.
C
C              TRIGS
C               AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST.
C
C              IFAX
C               AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST.
C
C
C               N.B. IN THE FOLLOWING ARGUMENTS, INCREMENTS ARE MEASURED
C               IN WORD PAIRS, BECAUSE EACH COMPLEX ELEMENT IS ASSUMED
C               TO OCCUPY AN ADJACENT PAIR OF WORDS IN MEMORY.
C
C              INC
C               THE INCREMENT (IN WORD PAIRS) BETWEEN SUCCESSIVE ELEMENT
C               OF EACH (COMPLEX) GRIDPOINT OR COEFFICIENT VECTOR
C               (E.G.  INC=1 FOR CONSECUTIVELY STORED DATA).
C
C              JUMP
C               THE INCREMENT (IN WORD PAIRS) BETWEEN THE FIRST ELEMENTS
C               OF SUCCESSIVE DATA OR COEFFICIENT VECTORS.  ON THE CRAY-
C               TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8
C               (TO AVOID MEMORY BANK CONFLICTS).  FOR CLARIFICATION OF
C               INC AND JUMP, SEE THE EXAMPLES BELOW.
C
C              N
C               THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF
C               TRANSFORMS, BELOW).
C
C              LOT
C               THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY.
C
C              ISIGN
C               = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER
C                    COEFFICIENTS.
C               = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO
C                    GRIDPOINT VALUES.
C
C ON OUTPUT    A
C               IF ISIGN = -1, AND LOT GRIDPOINT VECTORS ARE SUPPLIED,
C               EACH CONTAINING THE COMPLEX SEQUENCE:
C
C               G(0),G(1), ... ,G(N-1)  (N COMPLEX VALUES)
C
C               THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH
C               CONTAINING THE CORRESPONDING N COEFFICIENT VALUES:
C
C               C(0),C(1), ... ,C(N-1)  (N COMPLEX VALUES)
C
C               DEFINED BY:
C                 C(K) = SUM(J=0,...,N-1)( G(J)*EXP(-2*I*J*K*PI/N) )
C                 WHERE I = SQRT(-1)
C
C
C               IF ISIGN = +1, AND LOT COEFFICIENT VECTORS ARE SUPPLIED,
C               EACH CONTAINING THE COMPLEX SEQUENCE:
C
C               C(0),C(1), ... ,C(N-1)  (N COMPLEX VALUES)
C
C               THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH
C               CONTAINING THE CORRESPONDING N GRIDPOINT VALUES:
C
C               G(0),G(1), ... ,G(N-1)  (N COMPLEX VALUES)
C
C               DEFINED BY:
C                 G(J) = SUM(K=0,...,N-1)( G(K)*EXP(+2*I*J*K*PI/N) )
C                 WHERE I = SQRT(-1)
C
C
C               A CALL WITH ISIGN=-1 FOLLOWED BY A CALL WITH ISIGN=+1
C               (OR VICE VERSA) RETURNS THE ORIGINAL DATA, MULTIPLIED
C               BY THE FACTOR N.
C
C
C EXAMPLE       GIVEN A 64 BY 9 GRID OF COMPLEX VALUES, STORED IN
C               A 66 BY 9 COMPLEX ARRAY, A, COMPUTE THE TWO DIMENSIONAL
C               FOURIER TRANSFORM OF THE GRID.  FROM TRANSFORM THEORY,
C               IT IS KNOWN THAT A TWO DIMENSIONAL TRANSFORM CAN BE
C               OBTAINED BY FIRST TRANSFORMING THE GRID ALONG ONE
C               DIRECTION, THEN TRANSFORMING THESE RESULTS ALONG THE
C               ORTHOGONAL DIRECTION.
C
C               COMPLEX A(66,9), WORK(64,9)
C               REAL TRIGS1(128), TRIGS2(18)
C               INTEGER IFAX1(13), IFAX2(13)
C
C               SET UP THE IFAX AND TRIGS ARRAYS FOR EACH DIRECTION:
C
C               CALL CFTFAX(64, IFAX1, TRIGS1)
C               CALL CFTFAX( 9, IFAX2, TRIGS2)
C
C               IN THIS CASE, THE COMPLEX VALUES OF THE GRID ARE
C               STORED IN MEMORY AS FOLLOWS (USING U AND V TO
C               DENOTE THE REAL AND IMAGINARY COMPONENTS, AND
C               ASSUMING CONVENTIONAL FORTRAN STORAGE):
C
C   U(1,1), V(1,1), U(2,1), V(2,1),  ...  U(64,1), V(64,1), 4 NULLS,
C
C   U(1,2), V(1,2), U(2,2), V(2,2),  ...  U(64,2), V(64,2), 4 NULLS,
C
C   .       .       .       .         .   .        .        .
C   .       .       .       .         .   .        .        .
C   .       .       .       .         .   .        .        .
C
C   U(1,9), V(1,9), U(2,9), V(2,9),  ...  U(64,9), V(64,9), 4 NULLS.
C
C               WE CHOOSE (ARBITRARILY) TO TRANSORM FIRST ALONG THE
C               DIRECTION OF THE FIRST SUBSCRIPT.  THUS WE DEFINE
C               THE LENGTH OF THE TRANSFORMS, N, TO BE 64, THE
C               NUMBER OF TRANSFORMS, LOT, TO BE 9, THE INCREMENT
C               BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 1,
C               AND THE INCREMENT BETWEEN THE STARTING POINTS
C               FOR EACH TRANSFORM, JUMP, TO BE 66 (THE FIRST
C               DIMENSION OF A).
C
C               CALL CFFT99( A, WORK, TRIGS1, IFAX1, 1, 66, 64, 9, -1)
C
C               TO TRANSFORM ALONG THE DIRECTION OF THE SECOND SUBSCRIPT
C               THE ROLES OF THE INCREMENTS ARE REVERSED.  THUS WE DEFIN
C               THE LENGTH OF THE TRANSFORMS, N, TO BE 9, THE
C               NUMBER OF TRANSFORMS, LOT, TO BE 64, THE INCREMENT
C               BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 66,
C               AND THE INCREMENT BETWEEN THE STARTING POINTS
C               FOR EACH TRANSFORM, JUMP, TO BE 1
C
C               CALL CFFT99( A, WORK, TRIGS2, IFAX2, 66, 1, 9, 64, -1)
C
C               THESE TWO SEQUENTIAL STEPS RESULTS IN THE TWO-DIMENSIONA
C               FOURIER COEFFICIENT ARRAY OVERWRITING THE INPUT
C               GRIDPOINT ARRAY, A.  THE SAME TWO STEPS APPLIED AGAIN
C               WITH ISIGN = +1 WOULD RESULT IN THE RECONSTRUCTION OF
C               THE GRIDPOINT ARRAY (MULTIPLIED BY A FACTOR OF 64*9).
C
C
C-----------------------------------------------------------------------
C*** Modification 90/11/29 (BTD)
C      DIMENSION A(1),WORK(1),TRIGS(1),IFAX(1)
      DIMENSION A(*),WORK(*),TRIGS(*),IFAX(*)
C
C     SUBROUTINE "CFFT99" - MULTIPLE FAST COMPLEX FOURIER TRANSFORM
C
C     A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
C     WORK IS AN AREA OF SIZE N*LOT
C     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
C     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N
C     INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
C         (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
C     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
C     N IS THE LENGTH OF THE DATA VECTORS
C     LOT IS THE NUMBER OF DATA VECTORS
C     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
C           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
C
C
C     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
C     PARALLEL.
C
C
C
C
C
      NN = N+N
      INK=INC+INC
      JUM = JUMP+JUMP
      NFAX=IFAX(1)
      JNK = 2
      JST = 2
      IF (ISIGN.GE.0) GO TO 30
C
C     THE INNERMOST TEMPERTON ROUTINES HAVE NO FACILITY FOR THE
C     FORWARD (ISIGN = -1) TRANSFORM.  THEREFORE, THE INPUT MUST BE
C     REARRANGED AS FOLLOWS:
C
C     THE ORDER OF EACH INPUT VECTOR,
C
C     G(0), G(1), G(2), ... , G(N-2), G(N-1)
C
C     IS REVERSED (EXCLUDING G(0)) TO YIELD
C
C     G(0), G(N-1), G(N-2), ... , G(2), G(1).
C
C     WITHIN THE TRANSFORM, THE CORRESPONDING EXPONENTIAL MULTIPLIER
C     IS THEN PRECISELY THE CONJUGATE OF THAT FOR THE NORMAL
C     ORDERING.  THUS THE FORWARD (ISIGN = -1) TRANSFORM IS
C     ACCOMPLISHED
C
C     FOR NFAX ODD, THE INPUT MUST BE TRANSFERRED TO THE WORK ARRAY,
C     AND THE REARRANGEMENT CAN BE DONE DURING THE MOVE.
C
      JNK = -2
      JST = NN-2
      IF (MOD(NFAX,2).EQ.1) GOTO 40
C
C     FOR NFAX EVEN, THE REARRANGEMENT MUST BE APPLIED DIRECTLY TO
C     THE INPUT ARRAY.  THIS CAN BE DONE BY SWAPPING ELEMENTS.
C
      IBASE = 1
      ILAST = (N-1)*INK
      NH = N/2
      DO 20 L=1,LOT
      I1 = IBASE+INK
      I2 = IBASE+ILAST
CDIR$ IVDEP
      DO 10 M=1,NH
C     SWAP REAL AND IMAGINARY PORTIONS
      HREAL = A(I1)
      HIMAG = A(I1+1)
      A(I1) = A(I2)
      A(I1+1) = A(I2+1)
      A(I2) = HREAL
      A(I2+1) = HIMAG
      I1 = I1+INK
      I2 = I2-INK
   10 CONTINUE
      IBASE = IBASE+JUM
   20 CONTINUE
      GOTO 100
C
   30 CONTINUE
      IF (MOD(NFAX,2).EQ.0) GOTO 100
C
   40 CONTINUE
C
C     DURING THE TRANSFORM PROCESS, NFAX STEPS ARE TAKEN, AND THE
C     RESULTS ARE STORED ALTERNATELY IN WORK AND IN A.  IF NFAX IS
C     ODD, THE INPUT DATA ARE FIRST MOVED TO WORK SO THAT THE FINAL
C     RESULT (AFTER NFAX STEPS) IS STORED IN ARRAY A.
C
      IBASE=1
      JBASE=1
      DO 60 L=1,LOT
C     MOVE REAL AND IMAGINARY PORTIONS OF ELEMENT ZERO
      WORK(JBASE) = A(IBASE)
      WORK(JBASE+1) = A(IBASE+1)
      I=IBASE+INK
      J=JBASE+JST
CDIR$ IVDEP
      DO 50 M=2,N
C     MOVE REAL AND IMAGINARY PORTIONS OF OTHER ELEMENTS (POSSIBLY IN
C     REVERSE ORDER, DEPENDING ON JST AND JNK)
      WORK(J) = A(I)
      WORK(J+1) = A(I+1)
      I=I+INK
      J=J+JNK
   50 CONTINUE
      IBASE=IBASE+JUM
      JBASE=JBASE+NN
   60 CONTINUE
C
  100 CONTINUE
C
C     PERFORM THE TRANSFORM PASSES, ONE PASS FOR EACH FACTOR.  DURING
C     EACH PASS THE DATA ARE MOVED FROM A TO WORK OR FROM WORK TO A.
C
C     FOR NFAX EVEN, THE FIRST PASS MOVES FROM A TO WORK
      IGO = 110
C     FOR NFAX ODD, THE FIRST PASS MOVES FROM WORK TO A
      IF (MOD(NFAX,2).EQ.1) IGO = 120
      LA=1
      DO 140 K=1,NFAX
      IF (IGO.EQ.120) GO TO 120
  110 CONTINUE
      CALL VPASSM(A(1),A(2),WORK(1),WORK(2),TRIGS,
     *   INK,2,JUM,NN,LOT,N,IFAX(K+1),LA)
      IGO=120
      GO TO 130
  120 CONTINUE
      CALL VPASSM(WORK(1),WORK(2),A(1),A(2),TRIGS,
     *    2,INK,NN,JUM,LOT,N,IFAX(K+1),LA)
      IGO=110
  130 CONTINUE
      LA=LA*IFAX(K+1)
  140 CONTINUE
C
C     AT THIS POINT THE FINAL TRANSFORM RESULT IS STORED IN A.
C
      RETURN
      END
      SUBROUTINE CFTFAX(N,IFAX,TRIGS)
C Modification 90/11/29 (BTD):
C      DIMENSION IFAX(13),TRIGS(1)
      DIMENSION IFAX(13),TRIGS(*)
C
C     THIS ROUTINE WAS MODIFIED FROM TEMPERTON"S ORIGINAL
C     BY DAVE FULKER.  IT NO LONGER PRODUCES FACTORS IN ASCENDING
C     ORDER, AND THERE ARE NONE OF THE ORIGINAL 'MODE' OPTIONS.
C
C ON INPUT     N
C               THE LENGTH OF EACH COMPLEX TRANSFORM TO BE PERFORMED
C
C               N MUST BE GREATER THAN 1 AND CONTAIN NO PRIME
C               FACTORS GREATER THAN 5.
C
C ON OUTPUT    IFAX
C               IFAX(1)
C                 THE NUMBER OF FACTORS CHOSEN OR -99 IN CASE OF ERROR
C               IFAX(2) THRU IFAX( IFAX(1)+1 )
C                 THE FACTORS OF N IN THE FOLLOWIN ORDER:  APPEARING
C                 FIRST ARE AS MANY FACTORS OF 4 AS CAN BE OBTAINED.
C                 SUBSEQUENT FACTORS ARE PRIMES, AND APPEAR IN
C                 ASCENDING ORDER, EXCEPT FOR MULTIPLE FACTORS.
C
C              TRIGS
C               2N SIN AND COS VALUES FOR USE BY THE TRANSFORM ROUTINE
C
      CALL FACT(N,IFAX)
      K = IFAX(1)
      IF (K .LT. 1 .OR. IFAX(K+1) .GT. 5) IFAX(1) = -99
      IF (IFAX(1) .LE. 0 )  CALL ERRMSG('FATAL','CFFT99F',
     $ '  FFTFAX - INVALID N')
      CALL CFTRIG (N, TRIGS)
      RETURN
      END
      SUBROUTINE FACT(N,IFAX)
C     FACTORIZATION ROUTINE THAT FIRST EXTRACTS ALL FACTORS OF 4
      DIMENSION IFAX(13)
      IF (N.GT.1) GO TO 10
      IFAX(1) = 0
      IF (N.LT.1) IFAX(1) = -99
      RETURN
   10 NN=N
      K=1
C     TEST FOR FACTORS OF 4
   20 IF (MOD(NN,4).NE.0) GO TO 30
      K=K+1
      IFAX(K)=4
      NN=NN/4
      IF (NN.EQ.1) GO TO 80
      GO TO 20
C     TEST FOR EXTRA FACTOR OF 2
   30 IF (MOD(NN,2).NE.0) GO TO 40
      K=K+1
      IFAX(K)=2
      NN=NN/2
      IF (NN.EQ.1) GO TO 80
C     TEST FOR FACTORS OF 3
   40 IF (MOD(NN,3).NE.0) GO TO 50
      K=K+1
      IFAX(K)=3
      NN=NN/3
      IF (NN.EQ.1) GO TO 80
      GO TO 40
C     NOW FIND REMAINING FACTORS
   50 L=5
      MAX = SQRT(FLOAT(NN))
      INC=2
C     INC ALTERNATELY TAKES ON VALUES 2 AND 4
   60 IF (MOD(NN,L).NE.0) GO TO 70
      K=K+1
      IFAX(K)=L
      NN=NN/L
      IF (NN.EQ.1) GO TO 80
      GO TO 60
   70 IF (L.GT.MAX) GO TO 75
      L=L+INC
      INC=6-INC
      GO TO 60
   75 K = K+1
      IFAX(K) = NN
   80 IFAX(1)=K-1
C     IFAX(1) NOW CONTAINS NUMBER OF FACTORS
      RETURN
      END
      SUBROUTINE CFTRIG(N,TRIGS)
C Modification 90/11/29 (BTD):
C      DIMENSION TRIGS(1)
      DIMENSION TRIGS(*)
      PI=2.0*ASIN(1.0)
      DEL=(PI+PI)/FLOAT(N)
      L=N+N
      DO 10 I=1,L,2
      ANGLE=0.5*FLOAT(I-1)*DEL
      TRIGS(I)=COS(ANGLE)
      TRIGS(I+1)=SIN(ANGLE)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE VPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA)
      DIMENSION A(N),B(N),C(N),D(N),TRIGS(N)
C
C     SUBROUTINE "VPASSM" - MULTIPLE VERSION OF "VPASSA"
C     PERFORMS ONE PASS THROUGH DATA
C     AS PART OF MULTIPLE COMPLEX (INVERSE) FFT ROUTINE
C     A IS FIRST REAL INPUT VECTOR
C     B IS FIRST IMAGINARY INPUT VECTOR
C     C IS FIRST REAL OUTPUT VECTOR
C     D IS FIRST IMAGINARY OUTPUT VECTOR
C     TRIGS IS PRECALCULATED TABLE OF SINES & COSINES
C     INC1 IS ADDRESSING INCREMENT FOR A AND B
C     INC2 IS ADDRESSING INCREMENT FOR C AND D
C     INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S
C     INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S
C     LOT IS THE NUMBER OF VECTORS
C     N IS LENGTH OF VECTORS
C     IFAC IS CURRENT FACTOR OF N
C     LA IS PRODUCT OF PREVIOUS FACTORS
C
      DATA SIN36/0.587785252292473/,COS36/0.809016994374947/,
     *     SIN72/0.951056516295154/,COS72/0.309016994374947/,
     *     SIN60/0.866025403784437/
C
c*** diagnostic
c      write(0,*)'in vpassm, n=',n
c***
      M=N/IFAC
      IINK=M*INC1
      JINK=LA*INC2
      JUMP=(IFAC-1)*JINK
      IBASE=0
      JBASE=0
      IGO=IFAC-1
      IF (IGO.GT.4) RETURN
      GO TO (10,50,90,130),IGO
C
C     CODING FOR FACTOR 2
C
   10 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      DO 20 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 15 IJK=1,LOT
      C(JA+J)=A(IA+I)+A(IB+I)
      D(JA+J)=B(IA+I)+B(IB+I)
      C(JB+J)=A(IA+I)-A(IB+I)
      D(JB+J)=B(IA+I)-B(IB+I)
      I=I+INC3
      J=J+INC4
   15 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   20 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 40 K=LA1,M,LA
      KB=K+K-2
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      DO 30 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 25 IJK=1,LOT
      C(JA+J)=A(IA+I)+A(IB+I)
      D(JA+J)=B(IA+I)+B(IB+I)
      C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))
      D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))
      I=I+INC3
      J=J+INC4
   25 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   30 CONTINUE
      JBASE=JBASE+JUMP
   40 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 3
C
   50 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      DO 60 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 55 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
      C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))
      C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))
      D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))
      D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))
      I=I+INC3
      J=J+INC4
   55 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   60 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 80 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      DO 70 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 65 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
      C(JB+J)=
     *    C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      D(JB+J)=
     *    S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))))
     *   +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
      C(JC+J)=
     *    C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      D(JC+J)=
     *    S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))))
     *   +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
      I=I+INC3
      J=J+INC4
   65 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
   70 CONTINUE
      JBASE=JBASE+JUMP
   80 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 4
C
   90 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      ID=IC+IINK
      JD=JC+JINK
      DO 100 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 95 IJK=1,LOT
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
      C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
      D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))
      C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))
      C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))
      D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))
      D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))
      I=I+INC3
      J=J+INC4
   95 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  100 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 120 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      KD=KC+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      C3=TRIGS(KD+1)
      S3=TRIGS(KD+2)
      DO 110 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 105 IJK=1,LOT
      C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
      D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
      C(JC+J)=
     *    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
     *   -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
      D(JC+J)=
     *    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
     *   +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
      C(JB+J)=
     *    C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))
     *   -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
      D(JB+J)=
     *    S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)))
     *   +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
      C(JD+J)=
     *    C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))
     *   -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
      D(JD+J)=
     *    S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)))
     *   +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
      I=I+INC3
      J=J+INC4
  105 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  110 CONTINUE
      JBASE=JBASE+JUMP
  120 CONTINUE
      RETURN
C
C     CODING FOR FACTOR 5
C
  130 IA=1
      JA=1
      IB=IA+IINK
      JB=JA+JINK
      IC=IB+IINK
      JC=JB+JINK
      ID=IC+IINK
      JD=JC+JINK
      IE=ID+IINK
      JE=JD+JINK
      DO 140 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 135 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
      C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *  -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
      C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *  +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
      D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *  +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
      D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *  -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
      C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *  -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
      C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *  +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
      D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *  +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
      D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *  -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
      I=I+INC3
      J=J+INC4
  135 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  140 CONTINUE
      IF (LA.EQ.M) RETURN
      LA1=LA+1
      JBASE=JBASE+JUMP
      DO 160 K=LA1,M,LA
      KB=K+K-2
      KC=KB+KB
      KD=KC+KB
      KE=KD+KB
      C1=TRIGS(KB+1)
      S1=TRIGS(KB+2)
      C2=TRIGS(KC+1)
      S2=TRIGS(KC+2)
      C3=TRIGS(KD+1)
      S3=TRIGS(KD+2)
      C4=TRIGS(KE+1)
      S4=TRIGS(KE+2)
      DO 150 L=1,LA
      I=IBASE
      J=JBASE
CDIR$ IVDEP
      DO 145 IJK=1,LOT
      C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
      D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
      C(JB+J)=
     *    C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      D(JB+J)=
     *    S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      C(JE+J)=
     *    C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      D(JE+J)=
     *    S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I)))
     *      +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))))
     *   +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I)))
     *      -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
      C(JC+J)=
     *    C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      D(JC+J)=
     *    S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      C(JD+J)=
     *    C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      D(JD+J)=
     *    S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I)))
     *      +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))))
     *   +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I)))
     *      -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
      I=I+INC3
      J=J+INC4
  145 CONTINUE
      IBASE=IBASE+INC1
      JBASE=JBASE+INC2
  150 CONTINUE
      JBASE=JBASE+JUMP
  160 CONTINUE
      RETURN
      END
      SUBROUTINE C3DFFT(CZ,I1,I2,I3,I4,I5,I6,I7)
C***********************************************************************
C    Note: This module should NOT be linked into DDSCAT when
C    running on a Convex with the vector library -- in that case
C    link with the -lveclib flag and use the native Convex C3DFFT
C    routine.
C    On other systems, however, you need to link in this routine
C    to avoid an error message from the compiler.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      COMPLEX CZ(1)
      INTEGER I1,I2,I3,I4,I5,I6,I7
      CALL ERRMSG('FATAL','C3DFFT',' CONVEX option but dummy routine!')
      RETURN
      END

      SUBROUTINE WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT,
     &                    CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET,
     &                    MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH,
     &                    NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV,
     &                    NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ,
     &                    MXCOMP,BETMID,BETMXD,THTMID,THTMXD,
     &                    PHIMID,PHIMXD,ICTHM,IPHIM,TOL,WAVEA,AEFFA,
     &                    THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD,
     &                    ITHETA,IBETA,IPHI,IORI,AEFF,WAVE,XX,AK1,BETAD,
     &                    THETAD,PHID,AKR,CXE01R,CXE02R,CXRFR,CXEPS,
     &                    QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB,
     &                    QTRQSC,SM,SMORI,PHIN,THETAN,TIMERS)
C***********************************************************************
C Dummy version of subroutine WRITENET
C (genuine subroutine WRITENET writes NetCDF portable binary file
C to file cnetfile 'dd.cdf')
C
C This dummy version is provided in order that use of DDSCAT is
C possible even if NetCDF software is not installed on system.
C
C History:
C 96.11.15 (PJF): First version of NetCDF  (Wolff's request)
C                 5a6 ---> 5a7 Converted to subroutine.
C 96.11.20 (BTD): dummy version created from real WRITENET
C end history
C Copyright (C) 1996 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C     .. Scalar Arguments ..
      REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD,
     &     THTMID,THTMXD,TOL,WAVE,XX
      INTEGER IBETA,ICTHM,IDNC,IORI,IORTH,IPHI,IPHIM,IRAD,ITHETA,IWAV,
     &        MXBETA,MXCOMP,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA,
     &        MXTHET,MXTIMERS,MXWAV,NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI,
     &        NRAD,NSCAT,NTHETA,NTIMERS,NWAV,NX,NY,NZ
      CHARACTER CALPHA*6,CMDFFT*6,CMDTRQ*6,CNETFLAG*6,CSHAPE*6,
     &          CNETFILE*14,CSTAMP*22,CDESCR*67
C     .. Array Arguments ..
      COMPLEX CXE01R(3),CXE02R(3),CXEPS(MXCOMP),CXRFR(MXCOMP)
      REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),G(2),PHI(MXPHI),
     &     PHIN(MXSCA),QABS(2),QAV(17),QBKSCA(2),QEXT(2),QPHA(2),
     &     QSCAT(2),QTRQAB(3,2),QTRQSC(3,2),SM(MXSCA,4,4),
     &     SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA),
     &     TIMERS(MXTIMERS),WAVEA(MXWAV)
      INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3)
C
C This subroutine should never be called, so if it is, bail out...
C
      CALL ERRMSG('FATAL','DUMMYWRITENET',
     &   ' NetCDF is not installed. Change to NOTCDF in ddscat.par!')
      RETURN
      END

