SHAPE Subroutine

                   S H A P E   P R O G R A M

    F O R  A R B I T R A R Y  V O R O N O I  P O L Y H E D R A


                                                       N.Stefanou

 In order to improve the efficency of the code and to
 check more, this program was changed in summer 1998 by N.Stefanou
 ATTENTION: BUG is removed! 
            In the old version IPAN=-1 is wrong and has to be
            set to IPAN=0.

 THIS PROGRAM CALCULATES THE ANGULAR MOMENTUM COMPONENTS OF THE 
 SHAPEFUNCTION FOR AN ARBITRARY VORONOI POLYHEDRON.
 A REAL SPHERICAL HARMONIC BASIS IS USED FOR THE DECOMPOSITION.
 ON INPUT WE GIVE :

    LMAX          :  MAXIMUM ANGULAR MOMENTUM

    DLT           :  DEFINES THE STEP FOR GAUSS-LEGENDRE CALC.

    NFACE         :  NUMBER OF FACES OF THE POLYHEDRON
    KEYPAN        :  KEY TO DEFINE  THE  RADIAL  MESH.  IF KEYPAN=
                     THE DEFAULT  RADIAL  DIVISION  OF PANNELS GIVE
                     IN DATA STATEMENT IS USED.OTHERWISE THE  NUMBE
                     OF  RADIAL  MESH  POINTS PER PANNEL  (NM(IPAN)
                     IS READ IN INPUT
                  ** IN THIS VERSION THE MESH IS DETERMINED
                     BY SUBROUTINE MESH0.
    Z(I)          :  COEFFICIENTS OF THE EQUATION OF A FACE
                     Z(1)*X + Z(2)*Y + Z(3)*Z  =  1
    NVERT         :  NUMBER OF VERTICES OF A FACE
    V(I,IVERT)    :  COORDINATES OF THE VERTICES OF A FACE
    NEWSCH(IFACE) :  INTEGER   PARAMETER TO CALCULATE   (=1)
                     THE CONTRIBUTION OF THE CORRESPONDING
                     PYRAMID TO THE SHAPEFUNCTIONS  .  IF
                     NEWSCH.NE.1 THE CONTRIBUTION IS TAKEN
                     EQUAL TO THAT OF THE PREVIOUS PYRAMID


 IN ORDER TO SAVE MEMORY WE STORE IN LOCAL TEMPORARY FILES IN  UNIT
 30+1 , 30+2 , ... , 30+NFACE THE TRANSFORMATION MATRICES ASSOCIATE
 WITH THE ROTATION OF EACH PYRAMID. THE TEMPORARY DIRECT ACCESS FIL
 IN UNIT 10 CONTAINS THE CALCULATED COMPONENTS OF THE SHAPEFUNCTIO

              ...........I N P U T  C A R D...(Bcc/fcc)

                                   IF not (Bcc/fcc) change main prg

63   32   30    7   21   15   15   15   15   15 \     
15   17   15   15   23   15   15    0    0    0  \ This is the
 0   17   15   15   23   15   15    0    0    0  / radial mesh info
 0    0    0    0    0    0    0    0    0    0 /
              .........................................


 THE DEFINITION OF REAL SPHERICAL HARMONICS IS NOT THE STANDARD  ON
 REFERED IN THE PAPER:
 N.STEFANOU,H.AKAI AND R.ZELLER,COMPUTER PHYS.COMMUN. 60 (1990) 231
 IF YOU WANT TO HAVE ANGULAR MOMENTUM COMPONENTS IN THE STANDARD BA
 SIS CHANGE THE FOLLOWING STATEMENTS IN THE ROUTINES :
 IN CCOEF      ISI=1                        ---->    ISI=1-2*MOD(M,
 IN DREAL      IF(MOD(M+MP),2).EQ.0) D=-D   ---->    DELETE THE LIN

 .. PARAMETER STATEMENTS ..