C     --------------------------------------------------------
C     THIS PROGRAM PROVIDES THE ICE IH STRUCTURE WITH ZERO
C     DIPOLE MOMENT AND VERY SAMLL MULTIPOLE MOMENT
C     --------------------------------------------------------

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 576)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD)
      COMMON /DIRECTION/ MTYPE(192)
      DIMENSION VECTOR(2,3,2)

      SQ2 = SQRT(2.)
      SQ3 = SQRT(3.)
      SQ32 = SQ3/2.

      NM = 192
      NP = NM*3      
      AA = 5.
      R  = SQ3/SQ2/2.*AA
      XLN = 4*AA
      YLN = 3*SQ3*AA
      ZLN = 4*SQ2/SQ3*AA
      BOND = 1.
      
      VECTOR(1,1,1) = 0.
      VECTOR(1,1,2) = 2*SQ2/3.
          
      VECTOR(1,2,1) = SQ2/SQ3
      VECTOR(1,2,2) = -SQ2/3.
      
      VECTOR(1,3,1) = -SQ2/SQ3
      VECTOR(1,3,2) = -SQ2/3.
      
      VECTOR(2,1,1) = -SQ2/SQ3
      VECTOR(2,1,2) = SQ2/3.
      
      VECTOR(2,2,1) = SQ2/SQ3
      VECTOR(2,2,2) = SQ2/3.
      
      VECTOR(2,3,1) = 0.
      VECTOR(2,3,2) = -2*SQ2/3.

      BASEX = 0.
      BASEY = 0.
      BASEZ = ZLN
      IP    = 1
      IM    = 1

      DO 100 I=1,8

      IH = MOD(I-1,4)+1
      IF (IH.EQ.1.OR.IH.EQ.4) THEN
      IQ = 1
      ELSE 
      IQ = 2
      ENDIF
      IF (IH.EQ.1.OR.IH.EQ.3) THEN
      PH = 1.
      ELSE 
      PH = -1.
      ENDIF

      DO 200 J=1,6
      DO 200 K=1,4

      IF (IH.EQ.1.OR.IH.EQ.4) THEN
      X(IP) = BASEX+(K-1)*AA+(1-MOD(J,2))*0.5*AA
      ELSE
      X(IP) = BASEX+(K-1)*AA-(1-MOD(J,2))*0.5*AA
      ENDIF

      Y(IP) = BASEY+(J-1)*SQ32*AA
      Z(IP) = BASEZ
      
      MT = MTYPE(IM)
   
      IF (MT.LT.0) THEN
      X(IP+1) = X(IP) 
      Y(IP+1) = Y(IP)
      Z(IP+1) = Z(IP)+PH*BOND
      MT = -MT
      X(IP+2) = X(IP)+VECTOR(IQ,MT,1)*BOND
      Y(IP+2) = Y(IP)+VECTOR(IQ,MT,2)*BOND
      Z(IP+2) = Z(IP)-PH/3.*BOND
     
      ELSE
      MT1 = MOD(MT,3)+1
      X(IP+1) = X(IP)+VECTOR(IQ,MT1,1)*BOND 
      Y(IP+1) = Y(IP)+VECTOR(IQ,MT1,2)*BOND
      Z(IP+1) = Z(IP)-PH/3.*BOND
      MT2 = MOD(MT+1,3)+1
      X(IP+2) = X(IP)+VECTOR(IQ,MT2,1)*BOND
      Y(IP+2) = Y(IP)+VECTOR(IQ,MT2,2)*BOND
      Z(IP+2) = Z(IP)-PH/3.*BOND
      ENDIF

      IP = IP+3
      IM = IM+1
200   CONTINUE

      IF (IH.EQ.1) THEN
      BASEX = BASEX+SQ2/SQ3*R
      BASEY = BASEY-SQ2/3.*R
      BASEZ = BASEZ-1./3.*R
      
      ELSEIF (IH.EQ.3) THEN
      BASEX = BASEX-SQ2/SQ3*R
      BASEY = BASEY+SQ2/3.*R
      BASEZ = BASEZ-1./3.*R

      ELSE
      BASEZ = BASEZ-R
      ENDIF

100   CONTINUE
           
      DO 400 I=1,NP
      IF (X(I).GT.XLN) X(I) = X(I)-XLN     
      IF (X(I).LT.0.)  X(I) = X(I)+XLN     
      IF (Y(I).GT.YLN) Y(I) = Y(I)-YLN     
      IF (Y(I).LT.0.)  Y(I) = Y(I)+YLN     
      IF (Z(I).GT.ZLN) Z(I) = Z(I)-ZLN     
      IF (Z(I).LT.0.)  Z(I) = Z(I)+ZLN     
400   CONTINUE

      DO 500 I=1,NP-2,3
      PRINT *,'NM=',(I-1)/3+1

      DO 600 J=1,NP
      RX = X(I)-X(J)
      RY = Y(I)-Y(J)
      RZ = Z(I)-Z(J)
      IF (RX.GT.XLN/2.) RX=RX-XLN
      IF (RX.LT.-XLN/2.) RX=RX+XLN
      IF (RY.GT.YLN/2.) RY=RY-YLN
      IF (RY.LT.-YLN/2.) RY=RY+YLN
      IF (RZ.GT.ZLN/2.) RZ=RZ-ZLN
      IF (RZ.LT.-ZLN/2.) RZ=RZ+ZLN
      R1 = SQRT(RX*RX+RY*RY+RZ*RZ)
      IF (R1.LT.3.2.AND.R1.NE.0) print *,R1
600   CONTINUE 
      PRINT *,' '
      PRINT *,' '
500   CONTINUE

      STOP
      END
C     ---------------------------------------------------------------

C     ---------------------------------------------------------------
      BLOCK DATA
      COMMON /DIRECTION/ MTYPE(192)

      DATA MTYPE/ 3,-1,1,-2,-1,3,-3,2,-1,1,-2,-2,1,3,2,-3,2,-3,
     1            -3,3,-2,1,-1,2,
     2            1,2,-2,-2,3,3,-3,2,-1,-2,1,-2,-1,1,-3,2,2,-1,
     3            3,-3,3,-3,1,-1,
     4            -1,-3,3,1,-2,-3,2,-1,1,3,-2,1,3,-2,1,-1,-2,2,
     5            -1,2,-3,3,-3,2,
     6            -1,1,-3,3,-2,2,2,1,3,-1,-2,3,3,-2,-1,-3,1,-1,
     7            -3,-3,2,-2,2,1,
     8            3,-3,2,-2,2,-3,-3,-3,-1,1,1,-2,-2,3,2,1,-1,1,
     9            3,3,-2,2,-1,-1,
     A            -3,-1,3,-2,-1,2,-1,2,-1,-3,1,1,-2,1,-3,2,2,-2,
     B            -2,1,3,-3,3,3,
     C            2,2,-1,2,1,-2,2,-2,3,3,-3,-3,3,-1,3,-3,-1,1,1,
     D            -1,-3,1,-2,-2,
     E            -3,2,-3,2,3,-1,1,-2,2,-1,3,1,-1,-3,-2,2,-3,1,3,
     F            -1,1,-2,3,-2/
      END
C     ----------------------------------------------------------------