C                  INTERNATIONAL AVS CENTER
C        (This disclaimer must remain at the top of all files)
C
C WARRANTY DISCLAIMER
C
C This module and the files associated with it are distributed free of charge.
C It is placed in the public domain and permission is granted for anyone to use,
C duplicate, modify, and redistribute it unless otherwise noted.  Some modules
C may be copyrighted.  You agree to abide by the conditions also included in
C the AVS Licensing Agreement, version 1.0, located in the main module
C directory located at the International AVS Center ftp site and to include
C the AVS Licensing Agreement when you distribute any files downloaded from
C that site.
C
C The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module provide absolutely
C NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
C the quality and performance of this software is with the user.  IN NO EVENT
C WILL The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module BE LIABLE TO
C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.
C
C This AVS module and associated files are public domain software unless
C otherwise noted.  Permission is hereby granted to do whatever you like with
C it, subject to the conditions that may exist in copyrighted materials. Should
C you wish to make a contribution toward the improvement, modification, or
C general performance of this module, please send us your comments:  why you
C liked or disliked it, how you use it, and most important, how it helps your
C work. We will receive your comments at avs@ncsc.org.
C
C Please send AVS module bug reports to avs@ncsc.org.
C
C  NOTE:  THIS MODULE AND SOURCE CODE IS FOR USE 
C  WITH THE AVS SOFTWARE ENVIRONMENT ONLY 
      SUBROUTINE VERTEX_NORMAL( NX, NY, P, I, J, FNX, FNY, FNZ )
*
      INTEGER  NX, NY, I, J
      REAL     P(NX,NY,3), FNX, FNY, FNZ
*
      REAL     AX, AY, AZ, BX, BY, BZ, LEN
*
*--   branch on grid region
*
*--   allocate x basis vector
      IF ( (I.GT.1) .AND. (I.LT.NX) ) THEN
*
*--     central area of the grid
        AX = P(I-1,J,1) - P(I+1,J,1)
        AY = P(I-1,J,2) - P(I+1,J,2)
        AZ = P(I-1,J,3) - P(I+1,J,3)
*
      ELSE IF (I.EQ.1) THEN
*
*--     left border
        AX = P(I,J,1) - P(I+1,J,1)
        AY = P(I,J,2) - P(I+1,J,2)
        AZ = P(I,J,3) - P(I+1,J,3)
*
      ELSE IF (I.EQ.NX) THEN
*
*--     right border
        AX = P(I-1,J,1) - P(I,J,1)
        AY = P(I-1,J,2) - P(I,J,2)
        AZ = P(I-1,J,3) - P(I,J,3)
*
      END IF
*
*--   allocate y basis vector
      IF ( (J.GT.1) .AND. (J.LT.NY) ) THEN
*
*--     central area of the grid
        BX = P(I,J-1,1) - P(I,J+1,1)
        BY = P(I,J-1,2) - P(I,J+1,2)
        BZ = P(I,J-1,3) - P(I,J+1,3)
*
      ELSE IF (J.EQ.1) THEN
*
*--     lower border
        BX = P(I,J,1) - P(I,J+1,1)
        BY = P(I,J,2) - P(I,J+1,2)
        BZ = P(I,J,3) - P(I,J+1,3)
*
      ELSE IF (J.EQ.NY) THEN
*
*--     upper border
        BX = P(I,J-1,1) - P(I,J,1)
        BY = P(I,J-1,2) - P(I,J,2)
        BZ = P(I,J-1,3) - P(I,J,3)
*
      END IF 
*
*--   vector product of basis vectors
      FNX = AY*BZ - AZ*BY
      FNY = AZ*BX - AX*BZ
      FNZ = AX*BY - AY*BX
*
*--   normalize
      LEN = SQRT( FNX*FNX + FNY*FNY + FNZ*FNZ )
      IF (LEN.GT.1E-9) THEN
        FNX = FNX / LEN
        FNY = FNY / LEN
        FNZ = FNZ / LEN
      ELSE
        PRINT *, 'ERROR : vertex normal has zero length'
        FNX = 0.0
        FNY = 0.0
        FNZ = 1.0
      END IF
*
      RETURN
      END