      PROGRAM CALLTARGET
      IMPLICIT NONE
C
C Program CALLTARGET is used to generate target arrays using target
C generation routines employed by DDSCAT.
C
C B.T.Draine, Princeton Univ. Obs., 91.01.01
C History:
C 91.05.24 (BTD): Modified to allow entry target names in lower case
C 91.07.16 (BTD): Replaced DO...ENDDO with DO #...# CONTINUE
C 92.12.16 (BTD): Added IDVSHP to argument list of TARGET
C 93.01.07 (BTD): Added options TWOELL and TWOAEL
C 94.05.15 (BTD): Added option CONELL
C 95.12.11 (BTD): Added option BLOCKS
C 96.01.04 (BTD): Added option DW1996
C 96.01.25 (BTD): Added option TWOSPH
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ here and in TARGET
C                 Added calls to SIZER
C 96.01.29 (BTD): Support selection of principal axis calculation
C                 in TAR2SP via SHPAR(6)
C 97.04.25 (BTD): Added option BLK_AN
C 97.04.30 (BTD): Removed code offering choice on target axis for
C                 HEXGON option (target axis now hardwired within
C                 TARHEX)
C 97.12.26 (BTD): Introduced lattice spacing vector DX(1-3).
C                 Added code to input DX interactively.
C                 Added DX to argument list in each of 15 separate
C                 calls to TARGET.
C 02.02.12 (BTD): Added options NSPHER and PRISM3
C end history
C
C Copyright (C) 1993,1994,1995,1996,1997,2002
C                B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C Adjustable parameters:
C    MXNAT=max.number of dipoles in target
      INTEGER MXNAT
      PARAMETER(MXNAT=200000)
C Derived parameters:
      INTEGER MXN3
      PARAMETER(MXN3=3*MXNAT)
C Variables:
      CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6
      INTEGER IDVSHP,IDVOUT,IOSHP,J,JJ,NAT,NAT3
      INTEGER LXYZ(3)
      INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3)
      REAL ERR,F,PI,RGYR2,RXY,RXZ,RYX,RYZ,RZX,RZY,
     &     X,XCM,Y,YCM,Z,ZCM
      REAL A1(3),A2(3),DX(3),SHPAR(6)
      DOUBLE PRECISION DXL,DYL,DZL,X2,Y2,Z2
C
C***********************************************************************
      DATA IDVOUT/0/,IOSHP/15/
      PI=4.*ATAN(1.)
C Define variables CFLSHP and IDVSHP to please flint:
      CFLSHP='none'
      IDVSHP=999
 1000 WRITE(IDVOUT,*)'What shape? (Enter choice in quotes) Choices:'
      WRITE(IDVOUT,*)'  BLK_AN = construct from anisotropic cubic ',
     &               'blocks, input from "blocks.par"'
      WRITE(IDVOUT,*)'  BLOCKS = construct from cubic blocks, input ',
     &               'from file "blocks.par"'
      WRITE(IDVOUT,*)'  CONELL = two concentric ellipsoids'
      WRITE(IDVOUT,*)'  CYLNDR = cylinder'
      WRITE(IDVOUT,*)'  DW1996 = 13 cube target used by Draine & ',
     &               'Weingartner (1996)'
      WRITE(IDVOUT,*)'  ELLIPS = ellipsoid'
      WRITE(IDVOUT,*)'  HEXGON = hexagonal prism'
      WRITE(IDVOUT,*)'  RCTNGL = rectangular solid'
      WRITE(IDVOUT,*)'  SPHERE = sphere'
      WRITE(IDVOUT,*)'  TETRAH = regular tetrahedron'
      WRITE(IDVOUT,*)'  THRAEL = three collinear touching anisotropic ',
     &               'ellipsoids'
      WRITE(IDVOUT,*)'  THRELL = three collinear touching ellipsoids'
      WRITE(IDVOUT,*)'  TWOAEL = two touching anisotropic ellipsoids'
      WRITE(IDVOUT,*)'  TWOELL = two touching isotropic ellipsoids'
      WRITE(IDVOUT,*)'  TWOSPH = two touching spheroids (with angle ',
     &               'phi between symm. axes'
      WRITE(IDVOUT,*)'  NSPHER = multisphere (union of spheres)'
      WRITE(IDVOUT,*)'  PRISM3 = triangular prism'
      READ(*,*)CSHAPE
      WRITE(IDVOUT,*)'Specify lattice spacing, e.g.,'
      WRITE(IDVOUT,*)'  1 1 1   (for uniform spacing)'
      WRITE(IDVOUT,*)' .5 1 1   (for double resolution in x direction)'
      READ(*,*)DX(1),DX(2),DX(3)
C
C Normalize DX so that DX(1)*DX(2)*DX(3)=1.
C
      X=(DX(1)*DX(2)*DX(3))**(1./3.)
      DX(1)=DX(1)/X
      DX(2)=DX(2)/X
      DX(3)=DX(3)/X
      WRITE(IDVOUT,*)'Renormalized lattice spacings:'
      WRITE(IDVOUT,*)' dx/d =',DX(1)
      WRITE(IDVOUT,*)' dy/d =',DX(2)
      WRITE(IDVOUT,*)' dz/d =',DX(3)
      IF(CSHAPE.EQ.'TETRAH'.OR.CSHAPE.EQ.'tetrah')THEN
         CSHAPE='TETRAH'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter tetrahedron side length/d ',
     &                     ' (0 to stop)'
            READ(*,*)SHPAR(1)
            SHPAR(2)=SHPAR(1)
            SHPAR(3)=SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(*,6000)SHPAR(1),NAT
         ENDDO
         STOP
      ELSEIF(CSHAPE.EQ.'SPHERE'.OR.CSHAPE.EQ.'sphere')THEN
         CSHAPE='SPHERE'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter sphere diameter/d (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)stop
            SHPAR(2)=SHPAR(1)
            SHPAR(3)=SHPAR(1)
            CSHAPE='ELLIPS'
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            RGYR2=0.
            XCM=0.
            YCM=0.
            ZCM=0.
            DO JJ=1,NAT
               X=IXYZ(JJ,1)
               Y=IXYZ(JJ,2)
               Z=IXYZ(JJ,3)
               XCM=XCM+X
               YCM=YCM+Y
               ZCM=ZCM+Z
               RGYR2=RGYR2+X*X+Y*Y+Z*Z+0.25
            ENDDO
            XCM=XCM/REAL(NAT)
            YCM=YCM/REAL(NAT)
            ZCM=ZCM/REAL(NAT)
            RGYR2=RGYR2/REAL(NAT)-XCM**2-YCM**2-ZCM**2
            F=(5./3.)*(4.*PI/REAL(3.*NAT))**(2./3.)*RGYR2
            WRITE(*,6100)SHPAR(1),SHPAR(2),SHPAR(3),NAT,F
         ENDDO
         STOP
      ELSEIF(CSHAPE.EQ.'ELLIPS'.OR.CSHAPE.EQ.'ellips')THEN
         CSHAPE='ELLIPS'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter ellipsoid diameters Dx/d,Dy/d,Dz/d',
     &                     ' (0,0,0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0.)stop
            CSHAPE='ELLIPS'
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            XCM=0.
            YCM=0.
            ZCM=0.
            DO JJ=1,NAT
               X=IXYZ(JJ,1)
               Y=IXYZ(JJ,2)
               Z=IXYZ(JJ,3)
               XCM=XCM+X
               YCM=YCM+Y
               ZCM=ZCM+Z
            ENDDO
            XCM=XCM/REAL(NAT)
            YCM=YCM/REAL(NAT)
            ZCM=ZCM/REAL(NAT)
            X2=0.
            Y2=0.
            Z2=0.
            DO JJ=1,NAT
               DXL=IXYZ(JJ,1)-XCM
               DYL=IXYZ(JJ,2)-YCM
               DZL=IXYZ(JJ,3)-ZCM
               X2=X2+DXL*DXL
               Y2=Y2+DYL*DYL
               Z2=Z2+DZL*DZL
            ENDDO
            X2=X2/NAT+1./12.
            Y2=Y2/NAT+1./12.
            Z2=Z2/NAT+1./12.
            ERR=20.*PI*SQRT(5.*X2*Y2*Z2)/(3.*NAT)-1.
            RXY=SQRT(X2/Y2)
            RYZ=SQRT(Y2/Z2)
            RZX=SQRT(Z2/X2)
            RYX=1./RXY
            RZY=1./RYZ
            RXZ=1./RZX
            X2=2.*SQRT(5.*X2)
            Y2=2.*SQRT(5.*Y2)
            Z2=2.*SQRT(5.*Z2)
            WRITE(*,6200)SHPAR(1),SHPAR(2),SHPAR(3),NAT,X2,Y2,Z2,ERR,
     &                   RYX,RZX,RYZ,RXY,RXZ,RZY
         ENDDO
         STOP
      ELSEIF(CSHAPE.EQ.'CYLNDR'.OR.CSHAPE.EQ.'cylndr')THEN
         CSHAPE='CYLNDR'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter cylinder length/d (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            WRITE(IDVOUT,*)'Enter cylinder diameter/d'
            READ(*,*)SHPAR(2)
            WRITE(IDVOUT,*)'Enter 1,2,3 for axis along x,y,z'
            READ(*,*)SHPAR(3)
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
         STOP
      ELSEIF(CSHAPE.EQ.'RCTNGL'.OR.CSHAPE.EQ.'rctngl')THEN
         CSHAPE='RCTNGL'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter x-length/d (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            WRITE(IDVOUT,*)'Enter y-length/d'
            READ(*,*)SHPAR(2)
            WRITE(IDVOUT,*)'Enter z-length/d'
            READ(*,*)SHPAR(3)
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
      ELSEIF(CSHAPE.EQ.'HEXGON'.OR.CSHAPE.EQ.'hexgon')THEN
         CSHAPE='HEXGON'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter length (along symmetry axis)',
     &                     ' (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            WRITE(IDVOUT,*)'Enter (hexagon diameter=2*side width)/d'
            READ(*,*)SHPAR(2)
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
      ELSEIF(CSHAPE.EQ.'TWOELL'.OR.CSHAPE.EQ.'twoell')THEN
         CSHAPE='TWOELL'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions/d (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
      ELSEIF(CSHAPE.EQ.'TWOAEL'.OR.CSHAPE.EQ.'twoael')THEN
         CSHAPE='TWOAEL'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions/d (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
      ELSEIF(CSHAPE.EQ.'THRELL'.OR.CSHAPE.EQ.'threll')THEN
         CSHAPE='THRELL'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions/d (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
      ELSEIF(CSHAPE.EQ.'THRAEL'.OR.CSHAPE.EQ.'thrael')THEN
         CSHAPE='THRAEL'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions/d (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
            CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            WRITE(IDVOUT,*)'NAT=',NAT
         ENDDO
      ELSEIF(CSHAPE.EQ.'CONELL'.OR.CSHAPE.EQ.'conell')THEN
         CSHAPE='CONELL'
         DO J=1,100
            WRITE(IDVOUT,*)'Enter length of outer ellipsoid in x,y,z',
     &                     'directions/d  (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0)STOP
            WRITE(IDVOUT,*)'Enter length of inner ellipsoid in x,y,z',
     &                     'directions/d'
            READ(*,*)SHPAR(4),SHPAR(5),SHPAR(6)
            IF(SHPAR(4).LE.SHPAR(1).AND.
     &         SHPAR(5).LE.SHPAR(2).AND.
     &         SHPAR(6).LE.SHPAR(3))THEN
               CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &              MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
               CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
            ELSE
               WRITE(IDVOUT,*)'Input error: inner ellipsoid is not ',
     &                        'contained within outer ellipsoid'
            ENDIF
         ENDDO
      ELSEIF(CSHAPE.EQ.'BLOCKS'.OR.CSHAPE.EQ.'blocks')THEN
         CSHAPE='BLOCKS'
         CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &               MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
         CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
      ELSEIF(CSHAPE.EQ.'BLK_AN'.OR.CSHAPE.EQ.'blk_an')THEN
         CSHAPE='BLK_AN'
         CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &               MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
         CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
      ELSEIF(CSHAPE.EQ.'DW1996'.OR.CSHAPE.EQ.'dw1996')THEN
         CSHAPE='DW1996'
         WRITE(0,*)'Enter block size (lattice units)'
         READ(*,*)SHPAR(1)
         CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &               MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
         CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
      ELSEIF(CSHAPE.EQ.'TWOSPH'.OR.CSHAPE.EQ.'twosph')THEN
         CSHAPE='TWOSPH'
         WRITE(0,*)'Enter length a_1/d and diam. b_1/d of 1st spheroid'
         READ(*,*)SHPAR(1),SHPAR(2)
         WRITE(0,*)'Enter length a_2/d and diam. b_2/d of 2nd spheroid'
         READ(*,*)SHPAR(3),SHPAR(4)
         WRITE(0,*)'Enter angle phi (deg) between axes a_1 and a_2'
         READ(*,*)SHPAR(5)
         WRITE(0,*)'Enter PRINAX=0 for a_1=(1,0,0),a_2=(0,1,0) in TF'
         WRITE(0,*)'            =1 for a_1,a_2=principal axes'
         READ(*,*)SHPAR(6)
         CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &               MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
         CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
       ELSEIF(CSHAPE.EQ.'NSPHER'.OR.CSHAPE.EQ.'nspher')THEN
         CSHAPE='NSPHER'
         WRITE(0,*)'Enter DIAMX = target extent in x-direction'
         READ(*,*)SHPAR(1)
         WRITE(0,*)'Enter PRINAX=0 for a_1=(1,0,0),a_2=(0,1,0) in TF'
         WRITE(0,*)'            =1 for a_1,a_2 = principal axes'
         READ(*,*)SHPAR(2)
         WRITE(0,*)'Enter name of sphere parameter file (e.g., ',
     &             'tarnsp.par [in quotes])'
         READ(*,*)CFLSHP
         CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &               MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
         CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
      ELSEIF(CSHAPE.EQ.'PRISM3'.OR.CSHAPE.EQ.'prism3')THEN
         CSHAPE='PRISM3'
         WRITE(0,*)'Enter a/d = first triangle side/d'
         READ(*,*)SHPAR(1)
         WRITE(0,*)'Enter b/a , c/a, L/a'
         READ(*,*)SHPAR(2),SHPAR(3),SHPAR(4)
         CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,
     &               MXNAT,SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT,MXN3,NAT3)
         CALL SIZER(MXNAT,NAT,IXYZ,LXYZ)
      ELSE

         WRITE(IDVOUT,*)'Do not know how to make shape = ',CSHAPE
         GOTO 1000
      ENDIF
 6000 FORMAT(' AX=',F7.4,' NAT=',I7)
 6100 FORMAT(' AX=',F7.4,' AY=',F7.4,' AZ=',F7.4,' NAT=',I7,' F=',F10.6)
 6200 FORMAT('    XV =',F7.4,'    YV =',F7.4,'    ZV =',F7.4,
     &       ' NAT=',I8,/,
     &       ' 2a_eff=',F7.4,' 2b_eff=',F7.4,' 2c_eff=',F7.4,
     &       ' err=',F8.7,/,
     &       '   y:x =',F7.4,'   z:x =',F7.4,'   y:z =',F7.4,/,
     &       '   x:y =',F7.4,'   x:z =',F7.4,'   z:y =',F7.4)
      END


      SUBROUTINE CXFFTW(CX,MX,MY,MZ,ISIGN)
      IMPLICIT NONE
C
C Arguments:
C
      INTEGER ISIGN,MX,MY,MZ
      COMPLEX CX(MX,MY,MZ)
C***********************************************************************
C Given:
C     CX = 3d array (complex)
C
C Returns:
C     CX = 3d complex fourier transform of input array CX
C
C Purpose:
C Interface to the FFTW (fastest FFT code in the West ?) 
C Required routines:
C FFTW package compiled in SINGLE PRECISION. 
C FFTW is available from http://www.fftw.org
C It was written in C by Matteo Frigo and Steven G. Johnson of MIT.
C To link you would typically do the following
C f77 TSTFFT.o -L/usr/local/lib -lsfftw (thus the library is called libsfftw.a)
C 
C
C Advantages. (a) The code works for arbitrary N and is comparable in speed 
C                 at least on an SGI to the GPFA code of Temperton 
C                 (more tests needed)
C             (b) The code has several parallel processing hooks including MPI
C Disadvantages
C             Written in C (for those of us who are FORTRAN programers)
C             Not always faster than the GPFA of Temperton. Sometimes slower. 
C             
C History:
C 00.06.21 (PJF): Original version
C 00.06.22 (BTD): minor cosmetic changes
C end history
C Copyright (C) 2000 
C               B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C
C WARNING: NOT TESTED IN DDSCAT YET
C NOT tested on vector machines like CRAY C90
C***********************************************************************
C
C Local variables:
C
      INTEGER MXOLD,MYOLD,MZOLD
      DATA MXOLD,MYOLD,MZOLD/0,0,0/
      SAVE MXOLD,MYOLD,MZOLD

C include file (this comes with FFTW)

      INTEGER FFTW_FORWARD,FFTW_BACKWARD
      PARAMETER (FFTW_FORWARD=-1,FFTW_BACKWARD=1)

      INTEGER FFTW_REAL_TO_COMPLEX,FFTW_COMPLEX_TO_REAL
      PARAMETER (FFTW_REAL_TO_COMPLEX=-1,FFTW_COMPLEX_TO_REAL=1)

      INTEGER FFTW_ESTIMATE,FFTW_MEASURE
      PARAMETER (FFTW_ESTIMATE=0,FFTW_MEASURE=1)

      INTEGER FFTW_OUT_OF_PLACE,FFTW_IN_PLACE,FFTW_USE_WISDOM
      PARAMETER (FFTW_OUT_OF_PLACE=0)
      PARAMETER (FFTW_IN_PLACE=8,FFTW_USE_WISDOM=16)

      INTEGER FFTW_THREADSAFE
      PARAMETER (FFTW_THREADSAFE=128)

C     Constants for the MPI wrappers (not used here yet
      INTEGER FFTW_TRANSPOSED_ORDER, FFTW_NORMAL_ORDER
      INTEGER FFTW_SCRAMBLED_INPUT, FFTW_SCRAMBLED_OUTPUT
      PARAMETER(FFTW_TRANSPOSED_ORDER=1, FFTW_NORMAL_ORDER=0)
      PARAMETER(FFTW_SCRAMBLED_INPUT=8192)
      PARAMETER(FFTW_SCRAMBLED_OUTPUT=16384)
C define pointers (whatever it is it is needed by FFTW)
      INTEGER P_FWD, P_RVS
      INTEGER N_DIM(3)
      SAVE P_FWD, P_RVS
C
C once per each mx, my, mz call
C
      IF(MX.NE.MXOLD .OR. MY.NE.MYOLD .OR. MZ.NE.MZOLD) THEN
         MXOLD=MX
         MYOLD=MY
         MZOLD=MZ
         N_DIM(1)=MX
         N_DIM(2)=MY
         N_DIM(3)=MZ
C destroy may not work for the first call, may neeed fix 
C (there is nothing to destroy)
         CALL FFTWND_F77_DESTROY_PLAN(P_FWD)
         CALL FFTWND_F77_DESTROY_PLAN(P_RVS)
         CALL FFTWND_F77_CREATE_PLAN(P_FWD,3,N_DIM,FFTW_FORWARD, 
     $                               FFTW_MEASURE+FFTW_IN_PLACE)
         CALL FFTWND_F77_CREATE_PLAN(P_RVS,3,N_DIM,FFTW_BACKWARD, 
     $                               FFTW_MEASURE+FFTW_IN_PLACE)
      ENDIF

C it seems that ISIGN=1 of Temperton is equivalent to FFTW_BACKWARD in FFTW
C I did not test normalization for ISIGN=-1 (can be different from DDSCAT convention) PJF
      IF(ISIGN.EQ.1) THEN
         CALL FFTWND_F77_ONE(P_RVS,CX,0)
      ELSE
         CALL FFTWND_F77_ONE(P_FWD,CX,0)
      ENDIF
C
      RETURN
      END



















      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C***********************************************************************
C This version of timeit uses the system call "etime" available under
C SunOS and some other bsd-like Unix systems (e.g., Convex unix).
C
C Copyright (C) 1993,1994,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C 94.06.20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing in CMSGTM='noprint'
C Arguments:
      CHARACTER CMSGTM*(*)
      REAL DTIME
C
C External system calls:
      REAL ETIME
      EXTERNAL ETIME
C
C Local variables:
      REAL OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C
C External subroutines:
      EXTERNAL WRIMSG
C
C Local variables to be saved:
      SAVE CSTA,OLDTIM
C
C Data statements:
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         DTIME=ETIME(TARRAY)-OLDTIM
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            IF(DTIME.GT.1.E4)THEN
               WRITE(CMSGNM,FMT='(F9.0, A)')DTIME,' CPU time (in secs)'
            ELSE
               WRITE(CMSGNM,FMT='(F9.3, A)')DTIME,' CPU time (in secs)'
            ENDIF
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END

      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_convex
C
C This version of timeit uses the system call "etime" available under
C Convex unix and other bsd-like Unix systems.
C
C 94.06.20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
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:
      CHARACTER CMSGTM* (*)
C Local variables:
      REAL DTIME,OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C External Subroutines:
      EXTERNAL WRIMSG
C External Functions:
      REAL ETIME
      EXTERNAL ETIME
C***********************************************************************
      DATA CSTA/'ONE'/
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         DTIME=ETIME(TARRAY)-OLDTIM
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_cray
C
C CRAY version. Call to second is executed.
C 94.06.20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C     .. Scalar Arguments ..
      CHARACTER CMSGTM* (*)
C     ..
C     .. Local Scalars ..
      REAL DTIME, OLDTIM
      CHARACTER CSTA*3, CMSGNM*70
C
C     ..
C
      save csta, time1
C     .. Data statements ..
      DATA CSTA/'ONE'/
C     ..
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         call second(time1)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         CALL SECOND(TIME2)
         DTIME=TIME2-TIME1
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C timeit_hp
C
C This version uses the HP-AUX system calls SYSCONF and TIMES
C 94.06.20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
C
C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C External system calls:
$ALIAS sysconf = 'sysconf' (%val)
      INTEGER SYSCONF,TIMES
      EXTERNAL SYSCONF,TIMES
C
C Arguments:
      CHARACTER CMSGTM*(*)
C
C Local variables:
      INTEGER*2 NAME
      INTEGER CLK_TCK,I
      INTEGER TMS(4)
      REAL DTIME, OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
C
C External subroutines:
      EXTERNAL WRIMSG

C Data statements:
      DATA CSTA/'ONE'/
C     ..
      IF(CSTA.EQ.'ONE')THEN
         NAME=2
         CLK_TCK=SYSCONF(NAME)
         OLDTIM=REAL(TMS(1)+TMS(2))/REAL(CLK_TCK)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         I=TIMES(TMS)
         DTIME=REAL(TMS(1)+TMS(2))/REAL(CLK_TCK)-OLDTIM
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_ibm6000
C 94/06/20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
C
C This version of timeit uses the system call "mclock" available under
C RISC6000/320 IBM AIX3.1
C
C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      CHARACTER CMSGTM* (*)
C
C Local variables:
      REAL DTIME, OLDTIM
      CHARACTER CSTA*3, CMSGNM*70
      REAL TARRAY(2)
C
C External Subroutines ..
      EXTERNAL WRIMSG
C
C External Functions ..
      REAL ETIME
      EXTERNAL ETIME
C
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         I1=MCLOCK()
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         I2=MCLOCK()
         IALL=I2-I1
         DTIME=FLOAT(IALL)/100.
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C***********************************************************************
C This version of timeit uses the system call "etime" available under
C SunOS and some other bsd-like Unix systems (e.g., DEC OSF).
C
C Copyright (C) 1993,1994,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C 94.06.20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing in CMSGTM='noprint'
C Arguments:
      CHARACTER CMSGTM*(*)
C
C External system calls:
      REAL ETIME
      EXTERNAL ETIME
C
C Local variables:
      REAL DTIME,OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C
C External subroutines:
      EXTERNAL WRIMSG
C
C Local variables to be saved:
      SAVE CSTA,OLDTIM
C
C Data statements:
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         DTIME=ETIME(TARRAY)-OLDTIM
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            IF(DTIME.GT.1.E4)THEN
               WRITE(CMSGNM,FMT='(F9.0, A)')DTIME,' CPU time (in secs)'
            ELSE
               WRITE(CMSGNM,FMT='(F9.3, A)')DTIME,' CPU time (in secs)'
            ENDIF
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END

      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C***********************************************************************
C This version of timeit uses the SGI IRIX system call "etime"
C
C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C 94.06.20 (PJF) add dtime to the formal parameters
C 95.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
C Arguments:
      CHARACTER CMSGTM*(*)
C
C External system calls:
      REAL ETIME
      EXTERNAL ETIME
C
C Local variables:
      REAL DTIME,OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C
C External subroutines:
      EXTERNAL WRIMSG
C
C Local variables to be saved:
      SAVE CSTA,OLDTIM
C
C Data statements:
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         DTIME=ETIME(TARRAY)-OLDTIM
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            IF(DTIME.GT.1.E4)THEN
               WRITE(CMSGNM,FMT='(F9.0, A)')DTIME,' CPU time (in secs)'
            ELSE
               WRITE(CMSGNM,FMT='(F9.3, A)')DTIME,' CPU time (in secs)'
            ENDIF
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END


      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_titan
C
C This version of TIMEIT is for Titan computers.
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 94/06/20 (PJF) add dtime to the formal parameters
C***********************************************************************
C     .. Scalar Arguments ..
      CHARACTER CMSGTM* (*)
C     ..
C     .. Local Scalars ..
      REAL DTIME, OLDTIM, TIME
      INTEGER IAD
      CHARACTER CSTA*3, CMSGNM*70
C     ..
C     .. External Subroutines ..
      EXTERNAL WRIMSG
C     ..
C     .. External Functions ..
      REAL CPUTIM
      EXTERNAL CPUTIM
C     ..
C     .. Data statements ..
      DATA CSTA/'ONE'/
C     .. 
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=CPUTIM(0.)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         TIME=CPUTIM(0.)
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
            WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_vms
C 94.06.20 (PJF) add dtime to the formal parameters
C 96.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
C This version of TIMEIT is for VMS systems.
C
C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
c***********************************************************************
C     .. Scalar Arguments ..
      CHARACTER CMSGTM* (*)
C     ..
C     .. Local Scalars ..
      REAL DTIME, OLDTIM, TIME
      INTEGER IAD
      CHARACTER CSTA*3, CMSGNM*70
C     ..
C     .. External Subroutines ..
      EXTERNAL LIB$INIT_TIMER, LIB$SHOW_TIMER, WRIMSG
C     ..
C     .. Data statements ..
      DATA CSTA/'ONE'/
C     ..
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         IAD=0
         CALL LIB$INIT_TIMER(IAD)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            CALL LIB$SHOW_TIMER(IAD,2)
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      PROGRAM TSTFFT
      IMPLICIT NONE
C Parameters:
      INTEGER MX,MY,MZ,NF235MX
        PARAMETER (MX=200, MY=200, MZ=200, NF235MX=45)
C Local variables:
      INTEGER IN,N,NX,NY,NZ
      INTEGER NF235(NF235MX)
      COMPLEX A(MX,MY,MZ),C(MX,MY,MZ),CW(MX,MY,MZ)
C***********************************************************************
C Purpose:
C Tests different 3D FFT implementations (CPU time).
C Required routines:
C cfft99f.o  cxfft3.o  cxfft3n.o cxfftw.o  
C errmsg.o  fourx.o 
C gpfa.o  timeit_sgi.o  tst.o  wrimsg.o gen.o
C
C History:
C 94.06.20 (PJF): Original version
C 94.12.18 (BTD): Minor changes in output formatting.
C                 Explicit declaration of all variables.
C 95.07.12 (BTD): Minor changes in formatting.
C 96.11.14 (PJF): getset not needed
C 00.07.05 (PJF): modified to compare Brenner, GPFA, and FFTW
C 03.07.15 (BTD): modified to also do non-prime-factorable cases
C Copyright (C) 1994,1995,2000,2003
C               B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C
C***********************************************************************
C Elements of NF235 are numbers.le.4096  which can be factored by 2,3,5:
       DATA NF235/2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,
     &    45,48,50,54,60,64,72,75,80,81,90,96,100,108,120,125,128,135,
     &    144,150,160,162,180,192,200/
       WRITE(0,*)'Calculation time for different 3-D FFT methods'
       WRITE(0,*)'=========================',
     &           '=========================='
       WRITE(0,*)' '
       open(unit=20,file='err.dat',status='unknown')
       open(unit=21,file='res.dat',status='unknown')
       WRITE(21,6021)
 6021 FORMAT(' CPU time (sec) for different 3-D FFT methods',/,
     &   ' Machine=2.0 GHz Xeon, 256kB L2 cache, pgf77 compiler ',/,
     &   ' parameter LVR = 64 ',/,
     &   ' LVR="length of vector registers" in gpfa2f,gpfa3f,gpfa5f',/,
     &   '============================================================',
     &   /,
     &   '                  Brenner    Temperton ',
     &   'Frigo & Johnson',/,
     &   '   NX  NY  NZ     (FOURX)     (GPFA)   ',
     &   '    (FFTW)')
      rewind(20)
C
C IN=33 corresponds to 100x100x100 FFT. 
C On 2 GHz Xeon, this requires approx 1.6 sec for Brenner, 
C                                     0.23 sec for FFTW, 
C                                     0.18 sec for GPFA
C IN=40 corresponds to 150x150x150 FFT.
C On 2 GHz Xeon, this requires approx 6.9 sec for Brenner,
C                                     0.87 sec for FFTW,
C                                     0.68 sec for GPFA 
      DO IN=18,45
         IF(NF235(IN-1).LT.NF235(IN)-1)THEN
            NX=NINT(.5*(NF235(IN-1)+NF235(IN)))
            NY=NX
            NZ=NX
            CALL CASE(A,C,CW,NX,NY,NZ,NF235,NF235MX)
         ENDIF
         Nx=NF235(IN)
         Ny=NX
         Nz=NX
         CALL CASE(A,C,CW,NX,NY,NZ,NF235,NF235MX)
      ENDDO
      STOP
      END

      SUBROUTINE CASE(A,C,CW,NX,NY,NZ,NF235,NF235MX)
      IMPLICIT NONE

C Arguments:

      INTEGER NX,NY,NZ,NF235MX
      INTEGER NF235(NF235MX)
      COMPLEX A(NX,NY,NZ),C(NX,NY,NZ),CW(NX,NY,NZ)

C Local variables:

      INTEGER I,J,JX,JY,JZ,K,L,NRPT
      INTEGER NN(3)
      REAL T0,T1,T2,T3
      COMPLEX C1, C2, C3, C4, C5, C6
C***********************************************************************
C Purpose: interface to 3D FFT codes (for use in FFT tests)
C History:
C 94.06.20 (PJF): Original version
C 94.12.19 (BTD): Explicit declaration of all variables.
C 03.07.15 (BTD): Modified to all use of NX,NY,NZ not factorizable
C                 as 2^a3^b5^c
C 03.07.15 (BTD): Modified to repeat for more accurate timings
C Copyright (C) 1994,2003, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
***********************************************************************C
C--------- BRENNR
         NN(1)=NX
         NN(2)=NY
         NN(3)=NZ

         CALL GEN(A,NX,NY,NZ)

         NRPT=1+INT(1E5/(NX*NY*NZ))

         CALL TIMEIT('BREN',T0)
         DO L=1,NRPT
            DO K=1,NZ
               DO J=1,NY
                  DO I=1,NX
                     C(I,J,K)=A(I,J,K)
                  ENDDO
               ENDDO
            ENDDO
            CALL FOURX(C,NN,3,1)
         ENDDO
         CALL TIMEIT('BREN',T1)
         T1=T1/REAL(NRPT)

         C1=C(1,1,1)
         C2=C(2,1,1)

C--------- GPFA
C initialize

c check whether NX,NY,NZ are of form 2^i3^j5^k
         JX=0
         JY=0
         JZ=0
         DO J=1,NF235MX
            IF(NX.LE.NF235(J).AND.JX.EQ.0)JX=NF235(J)
            IF(NY.LE.NF235(J).AND.JY.EQ.0)JY=NF235(J)
            IF(NZ.LE.NF235(J).AND.JZ.EQ.0)JZ=NF235(J)
         ENDDO

         CALL CXFFT3n(C,JX,JY,JZ,1)
         CALL GEN(A,NX,NY,NZ)

         NRPT=1+INT(1E6/(JX*JY*JZ))

         CALL TIMEIT('GPFA',T0)
         DO L=1,NRPT
            DO K=1,JZ
               DO J=1,JY
                  DO I=1,JX
                     C(I,J,K)=A(I,J,K)
                  ENDDO
               ENDDO
            ENDDO
            CALL CXFFT3n(C,JX,JY,JZ,1)
         ENDDO
         CALL TIMEIT('GPFA',T2)
         T2=T2/REAL(NRPT)
         C5=C(1,1,1)
         C6=C(2,1,1)

C--------- FFTW
C initialize

         CALL CXFFTW(C,NX,NY,NZ,1)
         CALL GEN(A,NX,NY,NZ)

         NRPT=1+INT(1E6/(NX*NY*NZ))

         CALL TIMEIT('FFTW',T0)
         DO L=1,NRPT
            DO K=1,NZ
               DO J=1,NY
                  DO I=1,NX
                     C(I,J,K)=A(I,J,K)
                  ENDDO
               ENDDO
            ENDDO
            CALL CXFFTW(C,NX,NY,NZ,1)
         ENDDO
         CALL TIMEIT('FFTW',T3)
         T3=T3/REAL(NRPT)
         C3=C(1,1,1)
         C4=C(2,1,1)

C--------- Print
         WRITE(21, '(1X,3I4,3F12.3)') NX,NY,NZ,T1,T2,T3 
         WRITE(20,100) NX, C1,C2
         WRITE(20,100) NX, C3,C4
         WRITE(20,100) NX, C5,C6
 100     FORMAT(1X,I5,4G15.4)
         RETURN
         END

      subroutine gen(c,nx,ny,nz)
C Arguments:
      INTEGER NX,NY,NZ
      complex c(nx,ny,nz)
C Local variables:
      INTEGER I1,I2,I3
      REAL F
C Purpose:
C generates a 3D matrix for use in 3D FFT CPU timmings
      do 103 i1=1, nx
         do 102 i2=1, ny
            do 101 i3=1, nz
c              f=i1*i2*i3
               f=sin(float(i1))+sin(float(i2))+sin(float(i3))
               c(i1,i2,i3)=cmplx(f*1.,f*0.5)
 101        CONTINUE
 102     CONTINUE
 103  continue
      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,DX,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 Subroutine  "writenet" writes NetCDF portable binary  file
C to file cnetfile ('dd.cdf')
C
C To use netCDF option 
c (a) set "include" statement in this routine to proper location of 
c     the file "netcdf.inc"  
c (b) set "ddscat.par"  variable "cnetflag" to 'ALLCDF' or 'ORICDF', 
c (c) modify Makefile to link with NetCDF library, by appropriately
c     defining the strings LIBNETCDF and LINKNETCDF
c     (see comments in the Makefile)
C
C Original version created by P.J.Flatau, University of California, SD
C History:
C 96.11.15 (PJF): First version of NetCDF  (Wolff's request)
C                 5a6 ---> 5a7 Converted to subroutine.
C 96.11.21 (PJF): Add IF(CNEFLAG.EQ.'ALLCDF')THEN and ENDIF
C 96.11.21 (BTD): Removed second occurrence of MXSCA in argument list
C                 here and call from WRITESCA
C 96.11.21 (BTD): added call to NCCLOS if called with CNETFLAG='NCCLOS'
C                 in order to get call to NCCLOS out of DDSCAT.f
C 98.01.11 (BTD): added DX(3) to argument list
C                 added variable IDDX
C 03.10.30 (BTD): corrected typo (missing & in col 6)
C end history
C Copyright (C) 1996,1998 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C IMPORTANT warning:
C NetCDF "include" statement. "include" is NOT standard FORTRAN 77.
C For location of "netcdf.inc" ask your system administrator.
C Flatau's system:
C      include '/usr/local/netcdf-2.4.2/include/netcdf.inc'
C Draine's system:
      include '/u/flatau/netcdf-2.4.3/include/netcdf.inc'
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),DX(3),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     .. Local Scalars ..
      INTEGER ID17,ID2,ID3,IDA1,IDA2,IDAEFF,IDAEFFA,IDAK1,IDBETA,
     &        IDBETMID,IDBETMXD,IDC22,IDC6,IDC67,IDCALPHA,IDCDESCR,
     &        IDCMDFFT,IDCMDTRQ,IDCSHAPE,IDCSTAMP,IDCXRFR,IDDX,IDG,
     &        IDICOMP,IDICTHM,IDIORTH,IDIPHIM,IDIX0,IDIY0,IDIZ0,IDNAT,
     &        IDNAT0,IDNAT3,IDNBETA,IDNCOMP,IDNORI,IDNPHI,IDNRAD,
     &        IDNSCAT,IDNTHETA,IDNTIMERS,IDNWAV,IDNX,IDNY,IDNZ,IDONE,
     &        IDPHI,IDPHIMID,IDPHIMXD,IDPHIN,IDQABS,IDQAV,IDQBKSCA,
     &        IDQEXT,IDQPHA,IDQSCAT,IDQTRQAB,IDQTRQSC,IDSM11,IDSM12,
     &        IDSM13,IDSM14,IDSM21,IDSM22,IDSM23,IDSM24,IDSM31,IDSM32,
     &        IDSM33,IDSM34,IDSM41,IDSM42,IDSM43,IDSM44,IDSMORI11,
     &        IDSMORI12,IDSMORI13,IDSMORI14,IDSMORI21,IDSMORI22,
     &        IDSMORI23,IDSMORI24,IDSMORI31,IDSMORI32,IDSMORI33,
     &        IDSMORI34,IDSMORI41,IDSMORI42,IDSMORI43,IDSMORI44,
     &        I,IDTHETA,IDTHETAN,DTHTMID,IDTHTMXD,IDTOL,IDWAVE,IDWAVEA,
     &        IDXX,IENTRY,IERR
C (defined in netcdf.inc)      &        NCCHAR,NCCLOB,NCFLOAT,NCLONG,NCSHORT
C     .. Local Arrays ..
      INTEGER ICOUNT(7),IDDIM(7),ISTART(7)
C     .. External Functions ..
C (defined in netcdf.inc)       INTEGER NCCRE,NCDDEF,NCVDEF
C (defined in netcdf.inc)      EXTERNAL NCCRE,NCDDEF,NCVDEF
C     .. External Subroutines ..
      EXTERNAL NCAPTC,NCENDF,NCVPT,NCVPT1,NCVPTC,WRIMSG
C     .. Save statement ..
      SAVE IENTRY,IDQEXT,IDQABS,IDQSCAT,IDG,IDQBKSCA,IDQPHA,IDQAV,
     &     IDQTRQAB,IDQTRQSC,IDSM11,IDSM12,IDSM13,IDSM14,IDSM21,IDSM22,
     &     IDSM23,IDSM24,IDSM31,IDSM32,IDSM33,IDSM34,IDSM41,IDSM42,
     &     IDSM43,IDSM44,IDSMORI11,IDSMORI12,IDSMORI13,IDSMORI14,
     &     IDSMORI21,IDSMORI22,IDSMORI23,IDSMORI24,IDSMORI31,IDSMORI32,
     &     IDSMORI33,IDSMORI34,IDSMORI41,IDSMORI42,IDSMORI43,IDSMORI44,
     &     IDWAVE,IDAEFF,IDXX,IDAK1,IDCXRFR
C     .. Data statements ..

      DATA IENTRY/0/

C 1) IENTRY=0 defines first entry to this routine. All NetCDF
C    dimensions will be  defined in this step. This part of the
C     code is never executed again.NetCDF "id" variables are saved.
C 2) iori=0. Routine is called for iwav, irad, itheta, ibeta, iphi
C    These 5 variables can be split to two groups:
C    ("iwav", "irad") and ("itheta", "iphi", "ibeta").
C    The second group defines orientations of target and is also
C    defined by variable "iori". If "iori=0"  the code
C    finished all specified ntheta, nbeta, nphi orientations
C    and it is time to write orientational averages.
C 3) If iori .ne. 0 we cycle over orientations and iwav, irad.
C    Many-dimensional arrays are stored on disk.

      IF(CNETFLAG.EQ.'NCCLOS')THEN
         CALL NCCLOS(IDNC,IERR)
         RETURN
      ENDIF
      IF (IENTRY.EQ.0) THEN
          CALL WRIMSG('WRITENET',' First entry to NetCDF write')
          IENTRY = 1
C this dd.cdf should be "cnetfile"
          IDNC = NCCRE('dd.cdf',NCCLOB,IERR)
C NetCDF dimensions
C character variables
          IDC6 = NCDDEF(IDNC,'c6',6,IERR)
          IDC67 = NCDDEF(IDNC,'c67',67,IERR)
          IDC22 = NCDDEF(IDNC,'c22',22,IERR)
          IDCALPHA = NCVDEF(IDNC,'calpha',NCCHAR,1,IDC6,IERR)
          IDCDESCR = NCVDEF(IDNC,'cdescr',NCCHAR,1,IDC67,IERR)
          IDCMDFFT = NCVDEF(IDNC,'cmdfft',NCCHAR,1,IDC6,IERR)
          IDCMDTRQ = NCVDEF(IDNC,'cmdtrq',NCCHAR,1,IDC6,IERR)
          IDCSHAPE = NCVDEF(IDNC,'cshape',NCCHAR,1,IDC6,IERR)
          IDCSTAMP = NCVDEF(IDNC,'cstamp',NCCHAR,1,IDC22,IERR)
C single variables
          IDONE = NCDDEF(IDNC,'one',1,IERR)
          IDIORTH = NCVDEF(IDNC,'iorth',NCLONG,1,IDONE,IERR)
          IDBETMID = NCVDEF(IDNC,'betmid',NCFLOAT,1,IDONE,IERR)
          IDBETMXD = NCVDEF(IDNC,'betmxd',NCFLOAT,1,IDONE,IERR)
          IDTHTMID = NCVDEF(IDNC,'thtmid',NCFLOAT,1,IDONE,IERR)
          IDTHTMXD = NCVDEF(IDNC,'thtmxd',NCFLOAT,1,IDONE,IERR)
          IDPHIMID = NCVDEF(IDNC,'phimid',NCFLOAT,1,IDONE,IERR)
          IDPHIMXD = NCVDEF(IDNC,'phimxd',NCFLOAT,1,IDONE,IERR)
          IDICTHM = NCVDEF(IDNC,'icthm',NCLONG,1,IDONE,IERR)
          IDIPHIM = NCVDEF(IDNC,'iphim',NCLONG,1,IDONE,IERR)
          IDTOL = NCVDEF(IDNC,'tol',NCFLOAT,1,IDONE,IERR)
C dimensions arrays
          IDNX = NCDDEF(IDNC,'nx',NX,IERR)
          IDNY = NCDDEF(IDNC,'ny',NY,IERR)
          IDNZ = NCDDEF(IDNC,'nz',NZ,IERR)
          IDNAT0 = NCDDEF(IDNC,'nat0',NAT0,IERR)
          IDNAT = NCDDEF(IDNC,'nat',NAT,IERR)
          IDNAT3 = NCDDEF(IDNC,'nat3',NAT3,IERR)
          IDNSCAT = NCDDEF(IDNC,'nscat',NSCAT,IERR)
          IDNCOMP = NCDDEF(IDNC,'ncomp',NCOMP,IERR)
          IDNORI = NCDDEF(IDNC,'nori',NORI,IERR)
          IDNWAV = NCDDEF(IDNC,'nwav',NWAV,IERR)
          IDNRAD = NCDDEF(IDNC,'nrad',NRAD,IERR)
          IDNTHETA = NCDDEF(IDNC,'ntheta',NTHETA,IERR)
          IDNPHI = NCDDEF(IDNC,'nphi',NPHI,IERR)
          IDNBETA = NCDDEF(IDNC,'nbeta',NBETA,IERR)
          IDNTIMERS = NCDDEF(IDNC,'ntimers',NTIMERS,IERR)
          ID2 = NCDDEF(IDNC,'two',2,IERR)
          ID3 = NCDDEF(IDNC,'three',3,IERR)
          ID17 = NCDDEF(IDNC,'seventeen',17,IERR)
C arrays
          IDWAVEA = NCVDEF(IDNC,'wavea',NCFLOAT,1,IDNWAV,IERR)
          CALL NCAPTC(IDNC,IDWAVEA,'long_name',NCCHAR,16,
     &                'wavelength array',IERR)
          IDAEFFA = NCVDEF(IDNC,'aeffa',NCFLOAT,1,IDNRAD,IERR)
          CALL NCAPTC(IDNC,IDAEFFA,'long_name',NCCHAR,22,
     &                'effective radius array',IERR)
          IDTHETA = NCVDEF(IDNC,'theta',NCFLOAT,1,IDNTHETA,IERR)
          CALL NCAPTC(IDNC,IDTHETA,'long_name',NCCHAR,38,
     &                'theta (target orientation angle array)',IERR)
          IDBETA = NCVDEF(IDNC,'beta',NCFLOAT,1,IDNBETA,IERR)
          CALL NCAPTC(IDNC,IDBETA,'long_name',NCCHAR,37,
     &                'beta (target orientation angle array)',IERR)
          IDPHI = NCVDEF(IDNC,'phi',NCFLOAT,1,IDNPHI,IERR)
          CALL NCAPTC(IDNC,IDPHI,'long_name',NCCHAR,36,
     &                'phi (target orientation angle array)',IERR)

          IDA1 = NCVDEF(IDNC,'a1',NCFLOAT,1,ID3,IERR)
          CALL NCAPTC(IDNC,IDA1,'long_name',NCCHAR,24,
     &                'target definition vector',IERR)
          IDA2 = NCVDEF(IDNC,'a2',NCFLOAT,1,ID3,IERR)
          CALL NCAPTC(IDNC,IDA2,'long_name',NCCHAR,24,
     &                'target definition vector',IERR)

          IDDX = NCVDEF(IDNC,'dx',NCFLOAT,1,ID3,IERR)
          CALL NCAPTC(IDNC,IDDX,'long_name',NCCHAR,22,
     &                'lattice spacing vector',IERR)

          IDICOMP = NCVDEF(IDNC,'icomp',NCSHORT,1,IDNAT3,IERR)
          CALL NCAPTC(IDNC,IDICOMP,'long_name',NCCHAR,39,
     &                'dipole composition definition (0-empty)',IERR)
C let as split ixyz0
          IDIX0 = NCVDEF(IDNC,'ix0',NCSHORT,1,IDNAT,IERR)
          IDIY0 = NCVDEF(IDNC,'iy0',NCSHORT,1,IDNAT,IERR)
          IDIZ0 = NCVDEF(IDNC,'iz0',NCSHORT,1,IDNAT,IERR)
C declare "running" variables dimensions
cc         WRITE(IOBIN) IWAV, IRAD, IORI, (timers(i), i=1,ntimers)
cc         WRITE(iobin) AEFF,WAVE,XX, AK1, BETAD,THETAD,PHID,
cc     $         AKR, CXE01R, CXE02R,
cc     $         (CXRFR(i), i=1, ncomp), (cxeps(i),i=1, ncomp)
cc         WRITE(iobin) QEXT, QABS, QSCAT, G, QBKSCA, QPHA, QAV,
cc     $               QTRQAB, QTRQSC
c          idone=ncddef(idnc,'one',1,ierr)
c          idiorth=ncvdef(idnc,'iorth',nclong,1, idone,ierr)
C
          IDDIM(1) = IDNWAV
          IDWAVE = NCVDEF(IDNC,'wave',NCFLOAT,1,IDDIM,IERR)
          IDDIM(1) = IDNRAD
          IDAEFF = NCVDEF(IDNC,'aeff',NCFLOAT,1,IDDIM,IERR)
          IDDIM(1) = IDNWAV
          IDDIM(2) = IDNRAD
          IDXX = NCVDEF(IDNC,'xx',NCFLOAT,2,IDDIM,IERR)
          IDAK1 = NCVDEF(IDNC,'ak1',NCFLOAT,2,IDDIM,IERR)
C cxrfr is complex but NetCDF doesn't support complex:
          IDDIM(1) = IDNWAV
          IDDIM(2) = ID2
          IDDIM(3) = IDNCOMP
          IDCXRFR = NCVDEF(IDNC,'cxrfr',NCFLOAT,3,IDDIM,IERR)
          IDDIM(1) = IDNWAV
          IDDIM(2) = IDNRAD
          IDDIM(3) = IDNSCAT
          IDSMORI11 = NCVDEF(IDNC,'smori11',NCFLOAT,3,IDDIM,IERR)
          IDSMORI12 = NCVDEF(IDNC,'smori12',NCFLOAT,3,IDDIM,IERR)
          IDSMORI13 = NCVDEF(IDNC,'smori13',NCFLOAT,3,IDDIM,IERR)
          IDSMORI14 = NCVDEF(IDNC,'smori14',NCFLOAT,3,IDDIM,IERR)
          IDSMORI21 = NCVDEF(IDNC,'smori21',NCFLOAT,3,IDDIM,IERR)
          IDSMORI22 = NCVDEF(IDNC,'smori22',NCFLOAT,3,IDDIM,IERR)
          IDSMORI23 = NCVDEF(IDNC,'smori23',NCFLOAT,3,IDDIM,IERR)
          IDSMORI24 = NCVDEF(IDNC,'smori24',NCFLOAT,3,IDDIM,IERR)
          IDSMORI31 = NCVDEF(IDNC,'smori31',NCFLOAT,3,IDDIM,IERR)
          IDSMORI32 = NCVDEF(IDNC,'smori32',NCFLOAT,3,IDDIM,IERR)
          IDSMORI33 = NCVDEF(IDNC,'smori33',NCFLOAT,3,IDDIM,IERR)
          IDSMORI34 = NCVDEF(IDNC,'smori34',NCFLOAT,3,IDDIM,IERR)
          IDSMORI41 = NCVDEF(IDNC,'smori41',NCFLOAT,3,IDDIM,IERR)
          IDSMORI42 = NCVDEF(IDNC,'smori42',NCFLOAT,3,IDDIM,IERR)
          IDSMORI43 = NCVDEF(IDNC,'smori43',NCFLOAT,3,IDDIM,IERR)
          IDSMORI44 = NCVDEF(IDNC,'smori44',NCFLOAT,3,IDDIM,IERR)
C real 6-7 dimensional cases
          IDDIM(1) = IDNWAV
          IDDIM(2) = IDNRAD
          IDDIM(3) = IDNTHETA
          IDDIM(4) = IDNBETA
          IDDIM(5) = IDNPHI
          IDDIM(6) = ID2
          IDQEXT = NCVDEF(IDNC,'qext',NCFLOAT,6,IDDIM,IERR)
          IDQABS = NCVDEF(IDNC,'qabs',NCFLOAT,6,IDDIM,IERR)
          IDQSCAT = NCVDEF(IDNC,'qscat',NCFLOAT,6,IDDIM,IERR)
          IDG = NCVDEF(IDNC,'g',NCFLOAT,6,IDDIM,IERR)
          IDQBKSCA = NCVDEF(IDNC,'qbksca',NCFLOAT,6,IDDIM,IERR)
          IDQPHA = NCVDEF(IDNC,'qpha',NCFLOAT,6,IDDIM,IERR)
          IDDIM(6) = ID17
          IDQAV = NCVDEF(IDNC,'qav',NCFLOAT,6,IDDIM,IERR)
          IDDIM(6) = ID3
          IDDIM(7) = ID2
          IDQTRQAB = NCVDEF(IDNC,'qtrqab',NCFLOAT,7,IDDIM,IERR)
          IDQTRQSC = NCVDEF(IDNC,'qtrqsc',NCFLOAT,7,IDDIM,IERR)
          IDDIM(6) = IDNSCAT
          IF(CNETFLAG.EQ.'ALLCDF')THEN
             IDSM11 = NCVDEF(IDNC,'sm11',NCFLOAT,6,IDDIM,IERR)
             IDSM12 = NCVDEF(IDNC,'sm12',NCFLOAT,6,IDDIM,IERR)
             IDSM13 = NCVDEF(IDNC,'sm13',NCFLOAT,6,IDDIM,IERR)
             IDSM14 = NCVDEF(IDNC,'sm14',NCFLOAT,6,IDDIM,IERR)
             IDSM21 = NCVDEF(IDNC,'sm21',NCFLOAT,6,IDDIM,IERR)
             IDSM22 = NCVDEF(IDNC,'sm22',NCFLOAT,6,IDDIM,IERR)
             IDSM23 = NCVDEF(IDNC,'sm23',NCFLOAT,6,IDDIM,IERR)
             IDSM24 = NCVDEF(IDNC,'sm24',NCFLOAT,6,IDDIM,IERR)
             IDSM31 = NCVDEF(IDNC,'sm31',NCFLOAT,6,IDDIM,IERR)
             IDSM32 = NCVDEF(IDNC,'sm32',NCFLOAT,6,IDDIM,IERR)
             IDSM33 = NCVDEF(IDNC,'sm33',NCFLOAT,6,IDDIM,IERR)
             IDSM34 = NCVDEF(IDNC,'sm34',NCFLOAT,6,IDDIM,IERR)
             IDSM41 = NCVDEF(IDNC,'sm41',NCFLOAT,6,IDDIM,IERR)
             IDSM42 = NCVDEF(IDNC,'sm42',NCFLOAT,6,IDDIM,IERR)
             IDSM43 = NCVDEF(IDNC,'sm43',NCFLOAT,6,IDDIM,IERR)
             IDSM44 = NCVDEF(IDNC,'sm44',NCFLOAT,6,IDDIM,IERR)
          ENDIF
          IDPHIN = NCVDEF(IDNC,'phin',NCFLOAT,1,IDNSCAT,IERR)
          CALL NCAPTC(IDNC,IDPHIN,'long_name',NCCHAR,16,
     &                'scattering angle',IERR)
          IDTHETAN = NCVDEF(IDNC,'thetan',NCFLOAT,1,IDNSCAT,IERR)
          CALL NCAPTC(IDNC,IDTHETAN,'long_name',NCCHAR,16,
     &                'scattering angle',IERR)
          CALL NCENDF(IDNC,IERR)
C write characters
          CALL NCVPTC(IDNC,IDCALPHA,1,6,CALPHA,6,IERR)
          CALL NCVPTC(IDNC,IDCDESCR,1,67,CDESCR,67,IERR)
          CALL NCVPTC(IDNC,IDCMDFFT,1,6,CMDFFT,6,IERR)
          CALL NCVPTC(IDNC,IDCMDTRQ,1,6,CMDTRQ,6,IERR)
          CALL NCVPTC(IDNC,IDCSHAPE,1,6,CSHAPE,6,IERR)
          CALL NCVPTC(IDNC,IDCSTAMP,1,22,CSTAMP,22,IERR)
C write single variables:
          CALL NCVPT1(IDNC,IDIORTH,1,IORTH,IERR)
          CALL NCVPT1(IDNC,IDBETMID,1,BETMID,IERR)
          CALL NCVPT1(IDNC,IDBETMXD,1,BETMXD,IERR)
          CALL NCVPT1(IDNC,IDTHTMID,1,THTMID,IERR)
          CALL NCVPT1(IDNC,IDTHTMXD,1,THTMXD,IERR)
          CALL NCVPT1(IDNC,IDPHIMID,1,PHIMID,IERR)
          CALL NCVPT1(IDNC,IDPHIMXD,1,PHIMXD,IERR)
          CALL NCVPT1(IDNC,IDICTHM,1,ICTHM,IERR)
          CALL NCVPT1(IDNC,IDIPHIM,1,IPHIM,IERR)
          CALL NCVPT1(IDNC,IDTOL,1,TOL,IERR)
C write arrays
          CALL NCVPT(IDNC,IDWAVEA,1,NWAV,WAVEA,IERR)
          CALL NCVPT(IDNC,IDAEFFA,1,NRAD,AEFFA,IERR)
          CALL NCVPT(IDNC,IDTHETA,1,NTHETA,THETA,IERR)
          CALL NCVPT(IDNC,IDBETA,1,NBETA,BETA,IERR)
          CALL NCVPT(IDNC,IDPHI,1,NPHI,PHI,IERR)
          CALL NCVPT(IDNC,IDA1,1,3,A1,IERR)
          CALL NCVPT(IDNC,IDA2,1,3,A2,IERR)
          CALL NCVPT(IDNC,IDDX,1,3,DX,IERR)
          CALL NCVPT(IDNC,IDICOMP,1,NAT3,ICOMP,IERR)
C print, idix0
          CALL NCVPT(IDNC,IDIX0,1,NAT,IXYZ0(1,1),IERR)
          CALL NCVPT(IDNC,IDIY0,1,NAT,IXYZ0(1,2),IERR)
          CALL NCVPT(IDNC,IDIZ0,1,NAT,IXYZ0(1,3),IERR)
          CALL NCVPT(IDNC,IDPHIN,1,NSCAT,PHIN,IERR)
          CALL NCVPT(IDNC,IDTHETAN,1,NSCAT,THETAN,IERR)
C endif      "if(ientry .eq. 0) then"
      END IF


      IF (IORI.NE.0 .AND. CNETFLAG.EQ.'ALLCDF') THEN
C NetCDF option --- running variables
C first 5 dimensions are always the same:
C define "lower left"  corner of the (many dimensional) box

C some stuff will be written over and over
          ISTART(1) = IWAV
          ICOUNT(1) = 1
          CALL NCVPT(IDNC,IDWAVE,ISTART,ICOUNT,WAVE,IERR)
          ISTART(1) = IRAD
          ICOUNT(1) = 1
          CALL NCVPT(IDNC,IDAEFF,ISTART,ICOUNT,AEFF,IERR)
          ISTART(1) = IWAV
          ISTART(2) = IRAD
          ICOUNT(1) = 1
          ICOUNT(2) = 1
          CALL NCVPT(IDNC,IDXX,ISTART,ICOUNT,XX,IERR)
          CALL NCVPT(IDNC,IDAK1,ISTART,ICOUNT,AK1,IERR)
          ISTART(1) = IWAV
          ISTART(2) = 1
          ISTART(3) = 1
          ICOUNT(1) = 1
          ICOUNT(2) = 2
          ICOUNT(3) = NCOMP
          CALL NCVPT(IDNC,IDCXRFR,ISTART,ICOUNT,CXRFR,IERR)
C
          ISTART(1) = IWAV
          ISTART(2) = IRAD
          ISTART(3) = ITHETA
          ISTART(4) = IBETA
          ISTART(5) = IPHI
          ISTART(6) = 1
          ISTART(7) = 1
C define how many elements to write
          ICOUNT(1) = 1
          ICOUNT(2) = 1
          ICOUNT(3) = 1
          ICOUNT(4) = 1
          ICOUNT(5) = 1
          ICOUNT(6) = 2
          CALL NCVPT(IDNC,IDQEXT,ISTART,ICOUNT,QEXT,IERR)
          CALL NCVPT(IDNC,IDQABS,ISTART,ICOUNT,QABS,IERR)
          CALL NCVPT(IDNC,IDQSCAT,ISTART,ICOUNT,QSCAT,IERR)
          CALL NCVPT(IDNC,IDG,ISTART,ICOUNT,G,IERR)
          CALL NCVPT(IDNC,IDQBKSCA,ISTART,ICOUNT,QBKSCA,IERR)
          CALL NCVPT(IDNC,IDQPHA,ISTART,ICOUNT,QPHA,IERR)
          ICOUNT(6) = 17
          CALL NCVPT(IDNC,IDQAV,ISTART,ICOUNT,QAV,IERR)
          ICOUNT(6) = 3
          ICOUNT(7) = 2
          CALL NCVPT(IDNC,IDQTRQAB,ISTART,ICOUNT,QTRQAB,IERR)
          CALL NCVPT(IDNC,IDQTRQSC,ISTART,ICOUNT,QTRQSC,IERR)
C Muller matrix (IDL handles only 7 dimensions) - so 16 ncvpt's
          ICOUNT(6) = NSCAT
          CALL NCVPT(IDNC,IDSM11,ISTART,ICOUNT,SM(1,1,1),IERR)
          CALL NCVPT(IDNC,IDSM12,ISTART,ICOUNT,SM(1,1,2),IERR)
          CALL NCVPT(IDNC,IDSM13,ISTART,ICOUNT,SM(1,1,3),IERR)
          CALL NCVPT(IDNC,IDSM14,ISTART,ICOUNT,SM(1,1,4),IERR)
          CALL NCVPT(IDNC,IDSM21,ISTART,ICOUNT,SM(1,2,1),IERR)
          CALL NCVPT(IDNC,IDSM22,ISTART,ICOUNT,SM(1,2,2),IERR)
          CALL NCVPT(IDNC,IDSM23,ISTART,ICOUNT,SM(1,2,3),IERR)
          CALL NCVPT(IDNC,IDSM24,ISTART,ICOUNT,SM(1,2,4),IERR)
          CALL NCVPT(IDNC,IDSM31,ISTART,ICOUNT,SM(1,3,1),IERR)
          CALL NCVPT(IDNC,IDSM32,ISTART,ICOUNT,SM(1,3,2),IERR)
          CALL NCVPT(IDNC,IDSM33,ISTART,ICOUNT,SM(1,3,3),IERR)
          CALL NCVPT(IDNC,IDSM34,ISTART,ICOUNT,SM(1,3,4),IERR)
          CALL NCVPT(IDNC,IDSM41,ISTART,ICOUNT,SM(1,4,1),IERR)
          CALL NCVPT(IDNC,IDSM42,ISTART,ICOUNT,SM(1,4,2),IERR)
          CALL NCVPT(IDNC,IDSM43,ISTART,ICOUNT,SM(1,4,3),IERR)
          CALL NCVPT(IDNC,IDSM44,ISTART,ICOUNT,SM(1,4,4),IERR)
C endif   "if(iori. ne. 0  .and. cnetflag.eq.'ALLCDF') then"
      END IF

C Here only for iori=0. Write smori to netCDF:
      IF (IORI.EQ.0) THEN
          ISTART(1) = IWAV
          ISTART(2) = IRAD
          ISTART(3) = 1
          ICOUNT(1) = 1
          ICOUNT(2) = 1
          ICOUNT(3) = NSCAT
          CALL NCVPT(IDNC,IDSMORI11,ISTART,ICOUNT,SMORI(1,1,1),IERR)
          CALL NCVPT(IDNC,IDSMORI12,ISTART,ICOUNT,SMORI(1,1,2),IERR)
          CALL NCVPT(IDNC,IDSMORI13,ISTART,ICOUNT,SMORI(1,1,3),IERR)
          CALL NCVPT(IDNC,IDSMORI14,ISTART,ICOUNT,SMORI(1,1,4),IERR)
          CALL NCVPT(IDNC,IDSMORI21,ISTART,ICOUNT,SMORI(1,2,1),IERR)
          CALL NCVPT(IDNC,IDSMORI22,ISTART,ICOUNT,SMORI(1,2,2),IERR)
          CALL NCVPT(IDNC,IDSMORI23,ISTART,ICOUNT,SMORI(1,2,3),IERR)
          CALL NCVPT(IDNC,IDSMORI24,ISTART,ICOUNT,SMORI(1,2,4),IERR)
          CALL NCVPT(IDNC,IDSMORI31,ISTART,ICOUNT,SMORI(1,3,1),IERR)
          CALL NCVPT(IDNC,IDSMORI32,ISTART,ICOUNT,SMORI(1,3,2),IERR)
          CALL NCVPT(IDNC,IDSMORI33,ISTART,ICOUNT,SMORI(1,3,3),IERR)
          CALL NCVPT(IDNC,IDSMORI34,ISTART,ICOUNT,SMORI(1,3,4),IERR)
          CALL NCVPT(IDNC,IDSMORI41,ISTART,ICOUNT,SMORI(1,4,1),IERR)
          CALL NCVPT(IDNC,IDSMORI42,ISTART,ICOUNT,SMORI(1,4,2),IERR)
          CALL NCVPT(IDNC,IDSMORI43,ISTART,ICOUNT,SMORI(1,4,3),IERR)
          CALL NCVPT(IDNC,IDSMORI44,ISTART,ICOUNT,SMORI(1,4,4),IERR)
c endif    "if(iori.eq.0) then
      END IF

      RETURN

      END

