s6rcoor3.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    s6rcoor3      ../engine/source/elements/thickshell/solide6c/s6rcoor3.F
!||--- called by ------------------------------------------------------
!||    s6cforc3      ../engine/source/elements/thickshell/solide6c/s6cforc3.F
!||--- calls      -----------------------------------------------------
!||    s6cortho3     ../engine/source/elements/thickshell/solide6c/s6cortho3.F
!||    s6proj3       ../engine/source/elements/thickshell/solide6c/s6proj3.f
!||    vrrota3       ../engine/source/elements/thickshell/solide6c/vrrota3.F
!||--- uses       -----------------------------------------------------
!||    element_mod   ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE s6rcoor3(
     .   X  ,IXS , V  , W  ,GAMA0,GAMA,
     .   X1 , X2 , X3 , X4 , X5 ,  X6   ,
     .   Y1 , Y2 , Y3 , Y4 , Y5 ,  Y6   ,
     .   Z1 , Z2 , Z3 , Z4 , Z5 ,  Z6   ,
     .   VX1, VX2, VX3, VX4, VX5, VX6   , 
     .   VY1, VY2, VY3, VY4, VY5, VY6   ,
     .   VZ1, VZ2, VZ3, VZ4, VZ5, VZ6   ,
     .   VD2,VIS,OFFG,OFF,SAV,RHO,RHOO  ,
     .   R11, R12, R13, R21, R22, R23, R31, R32, R33,
     .   NC1,NC2,NC3,NC4,NC5,NC6,NGL,MXT,NGEO,
     .   IOUTPRT, VGAX, VGAY, VGAZ, VGA2,DI,
     .   NEL, XGAX, XGAY, XGAZ, XGXA2, XGYA2, XGZA2,
     .   XGXYA, XGYZA, XGZXA, IPARG,GAMA_R)
      use element_mod , only : nixs
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "vect01_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NEL
C     REAL
      my_real
     .   X(3,*),V(3,*),W(3,*), VIS(*),
     .   X1(*), X2(*), X3(*), X4(*), X5(*), X6(*), 
     .   Y1(*), Y2(*), Y3(*), Y4(*), Y5(*), Y6(*), 
     .   Z1(*), Z2(*), Z3(*), Z4(*), Z5(*), Z6(*), 
     .  VX1(*), VX2(*), VX3(*), VX4(*), VX5(*), VX6(*), 
     .  VY1(*), VY2(*), VY3(*), VY4(*), VY5(*), VY6(*), 
     .  VZ1(*), VZ2(*), VZ3(*), VZ4(*), VZ5(*), VZ6(*), 
     .  VD2(*), OFFG(*), OFF(*),  RHO(*), RHOO(*),
     .   R11(MVSIZ),R12(MVSIZ),R13(MVSIZ),
     .   R21(MVSIZ),R22(MVSIZ),R23(MVSIZ),
     .   r31(mvsiz),r32(mvsiz),r33(mvsiz),
     .   rx(mvsiz) , ry(mvsiz) , rz(mvsiz) ,
     .   sx(mvsiz) , sy(mvsiz) , sz(mvsiz) ,
     .   tx(mvsiz) , ty(mvsiz) , tz(mvsiz) ,
     .   gama0(nel,6),gama(mvsiz,6),
     .   vgax(*), vgay(*), vgaz(*), vga2(*),di(mvsiz,6),
     .   xgax(*), xgay(*), xgaz(*),
     .   xgxa2(mvsiz),xgya2(mvsiz),xgza2(mvsiz),
     .   xgxya(mvsiz),xgyza(mvsiz),xgzxa(mvsiz),gama_r(nel,6)
      double precision
     .  sav(nel,15)
      INTEGER NC1(*), NC2(*), NC3(*), NC4(*),
     .        NC5(*), NC6(*),  MXT(*), NGL(*),NGEO(*)
      INTEGER IXS(NIXS,*),IOUTPRT,IPARG(*)
CMasParINCLUDE 'scoor3.intmap.inc'
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
C     REAL
      my_real
     .   XL,YL,ZL
      my_real
     .   xx,yy,zz,xy,xz,yz,rtr(6),abc,xxyz2,zzxy2,yyxz2,deta   
      my_real
     .   off_l
C-----------------------------------------------
      off_l  = zero
C
      DO i=lft,llt
      vis(i)=zero
      ngeo(i)=ixs(10,i)
      ngl(i)=ixs(11,i)
      mxt(i)=ixs(1,i)
      nc1(i)=ixs(2,i)
      nc2(i)=ixs(3,i)
      nc3(i)=ixs(4,i)
      nc4(i)=ixs(6,i)
      nc5(i)=ixs(7,i)
      nc6(i)=ixs(8,i)
      rhoo(i)=rho(i)
      ENDDO
C----------------------------
C     NODAL COORDINATES     |
C----------------------------
      DO i=lft,llt
        x1(i)=x(1,nc1(i))
        y1(i)=x(2,nc1(i))
        z1(i)=x(3,nc1(i))
        x2(i)=x(1,nc2(i))
        y2(i)=x(2,nc2(i))
        z2(i)=x(3,nc2(i))
        x3(i)=x(1,nc3(i))
        y3(i)=x(2,nc3(i))
        z3(i)=x(3,nc3(i))
        x4(i)=x(1,nc4(i))
        y4(i)=x(2,nc4(i))
        z4(i)=x(3,nc4(i))
        x5(i)=x(1,nc5(i))
        y5(i)=x(2,nc5(i))
        z5(i)=x(3,nc5(i))
        x6(i)=x(1,nc6(i))
        y6(i)=x(2,nc6(i))
        z6(i)=x(3,nc6(i))
        off(i) = min(one,abs(offg(i)))
        off_l  = min(off_l,offg(i))
      ENDDO
C-----------
      DO i=lft,llt
       vx1(i)=v(1,nc1(i))
       vy1(i)=v(2,nc1(i))
       vz1(i)=v(3,nc1(i))
       vx2(i)=v(1,nc2(i))
       vy2(i)=v(2,nc2(i))
       vz2(i)=v(3,nc2(i))
       vx3(i)=v(1,nc3(i))
       vy3(i)=v(2,nc3(i))
       vz3(i)=v(3,nc3(i))
       vx4(i)=v(1,nc4(i))
       vy4(i)=v(2,nc4(i))
       vz4(i)=v(3,nc4(i))
       vx5(i)=v(1,nc5(i))
       vy5(i)=v(2,nc5(i))
       vz5(i)=v(3,nc5(i))
       vx6(i)=v(1,nc6(i))
       vy6(i)=v(2,nc6(i))
       vz6(i)=v(3,nc6(i))
      ENDDO
      IF(off_l<zero)THEN
        DO i=lft,llt
          IF(offg(i)<zero)THEN
            vx1(i)=zero
            vy1(i)=zero
            vz1(i)=zero
            vx2(i)=zero
            vy2(i)=zero
            vz2(i)=zero
            vx3(i)=zero
            vy3(i)=zero
            vz3(i)=zero
            vx4(i)=zero
            vy4(i)=zero
            vz4(i)=zero
            vx5(i)=zero
            vy5(i)=zero
            vz5(i)=zero
            vx6(i)=zero
            vy6(i)=zero
            vz6(i)=zero
          ENDIF
        ENDDO
      ENDIF
C-----------
C     prepares outputs by part.
C-----------
      IF(ioutprt/=0)THEN
       DO i=lft,llt
        vgax(i)=vx1(i)+vx2(i)+vx3(i)+vx4(i)+vx5(i)+vx6(i)
        vgay(i)=vy1(i)+vy2(i)+vy3(i)+vy4(i)+vy5(i)+vy6(i)
        vgaz(i)=vz1(i)+vz2(i)+vz3(i)+vz4(i)+vz5(i)+vz6(i)
        vga2(i)=vx1(i)*vx1(i)+vx2(i)*vx2(i)+vx3(i)*vx3(i)+vx4(i)*vx4(i)
     1         +vx5(i)*vx5(i)+vx6(i)*vx6(i)
     2         +vy1(i)*vy1(i)+vy2(i)*vy2(i)+vy3(i)*vy3(i)+vy4(i)*vy4(i)
     3         +vy5(i)*vy5(i)+vy6(i)*vy6(i)
     4         +vz1(i)*vz1(i)+vz2(i)*vz2(i)+vz3(i)*vz3(i)+vz4(i)*vz4(i)
     5         +vz5(i)*vz5(i)+vz6(i)*vz6(i)
       ENDDO
       IF(iparg(80)==1) THEN
         DO i=lft,llt
          xgax(i)=x1(i)+x2(i)+x3(i)+x4(i)+x5(i)+x6(i)
          xgay(i)=y1(i)+y2(i)+y3(i)+y4(i)+y5(i)+y6(i)
          xgaz(i)=z1(i)+z2(i)+z3(i)+z4(i)+z5(i)+z6(i)
          xgxa2(i)=x1(i)**2+x2(i)**2+x3(i)**2+x4(i)**2
     1            +x5(i)**2+x6(i)**2
          xgya2(i)=y1(i)**2+y2(i)**2+y3(i)**2+y4(i)**2
     1            +y5(i)**2+y6(i)**2
          xgza2(i)=z1(i)**2+z2(i)**2+z3(i)**2+z4(i)**2
     1            +z5(i)**2+z6(i)**2
          xgxya(i)=x1(i)*y1(i)+x2(i)*y2(i)+x3(i)*y3(i)+x4(i)*y4(i)
     1            +x5(i)*y5(i)+x6(i)*y6(i)
          xgyza(i)=y1(i)*z1(i)+y2(i)*z2(i)+y3(i)*z3(i)+y4(i)*z4(i)
     1            +y5(i)*z5(i)+y6(i)*z6(i)
          xgzxa(i)=z1(i)*x1(i)+z2(i)*x2(i)+z3(i)*x3(i)+z4(i)*x4(i)
     1            +z5(i)*x5(i)+z6(i)*x6(i)
         ENDDO
       ENDIF
      ENDIF
C  -----------
      DO i=lft,llt
        xl=one_over_6*(x1(i)+x2(i)+x3(i)+x4(i)+x5(i)+x6(i))
        yl=one_over_6*(y1(i)+y2(i)+y3(i)+y4(i)+y5(i)+y6(i))
        zl=one_over_6*(z1(i)+z2(i)+z3(i)+z4(i)+z5(i)+z6(i))
        x1(i)=x1(i)-xl
        y1(i)=y1(i)-yl
        z1(i)=z1(i)-zl
        x2(i)=x2(i)-xl
        y2(i)=y2(i)-yl
        z2(i)=z2(i)-zl
        x3(i)=x3(i)-xl
        y3(i)=y3(i)-yl
        z3(i)=z3(i)-zl
        x4(i)=x4(i)-xl
        y4(i)=y4(i)-yl
        z4(i)=z4(i)-zl
        x5(i)=x5(i)-xl
        y5(i)=y5(i)-yl
        z5(i)=z5(i)-zl
        x6(i)=x6(i)-xl
        y6(i)=y6(i)-yl
        z6(i)=z6(i)-zl
      ENDDO
C-----------
C     convected frame 
C-----------
       CALL s6cortho3(
     1   x1,      x2,      x3,      x4,
     2   x5,      x6,      y1,      y2,
     3   y3,      y4,      y5,      y6,
     4   z1,      z2,      z3,      z4,
     5   z5,      z6,      r11,     r12,
     6   r13,     r21,     r22,     r23,
     7   r31,     r32,     r33,     rx,
     8   ry,      rz,      sx,      sy,
     9   sz,      tx,      ty,      tz,
     a   nel)
C------stocker [R] in %GAMA_R--------------     
       gama_r(lft:llt,1) = r11(lft:llt)                             
       gama_r(lft:llt,2) = r21(lft:llt)                             
       gama_r(lft:llt,3) = r31(lft:llt)                             
       gama_r(lft:llt,4) = r12(lft:llt)                             
       gama_r(lft:llt,5) = r22(lft:llt)                             
       gama_r(lft:llt,6) = r32(lft:llt)                             
C-----------
C     Change to the convected frame.
C-----------
C   X=RX' <=> X'=t(R)X chgt de base.
      IF(ismstr<=4.AND.jlag>0) THEN
       DO i=lft,llt
       IF(offg(i)>one)THEN
        x1(i)=sav(i,1)
        y1(i)=sav(i,2)
        z1(i)=sav(i,3)
        x2(i)=sav(i,4)
        y2(i)=sav(i,5)
        z2(i)=sav(i,6)
        x3(i)=sav(i,7)
        y3(i)=sav(i,8)
        z3(i)=sav(i,9)
        x4(i)=sav(i,10)
        y4(i)=sav(i,11)
        z4(i)=sav(i,12)
        x5(i)=sav(i,13)
        y5(i)=sav(i,14)
        z5(i)=sav(i,15)
        x6(i)=zero
        y6(i)=zero
        z6(i)=zero
        off(i) = offg(i) -one
        xl=one_over_6*(x1(i)+x2(i)+x3(i)+x4(i)+x5(i)+x6(i))
        yl=one_over_6*(y1(i)+y2(i)+y3(i)+y4(i)+y5(i)+y6(i))
        zl=one_over_6*(z1(i)+z2(i)+z3(i)+z4(i)+z5(i)+z6(i))
        x1(i)=x1(i)-xl
        y1(i)=y1(i)-yl
        z1(i)=z1(i)-zl
        x2(i)=x2(i)-xl
        y2(i)=y2(i)-yl
        z2(i)=z2(i)-zl
        x3(i)=x3(i)-xl
        y3(i)=y3(i)-yl
        z3(i)=z3(i)-zl
        x4(i)=x4(i)-xl
        y4(i)=y4(i)-yl
        z4(i)=z4(i)-zl
        x5(i)=x5(i)-xl
        y5(i)=y5(i)-yl
        z5(i)=z5(i)-zl
        x6(i)=x6(i)-xl
        y6(i)=y6(i)-yl
        z6(i)=z6(i)-zl
       ELSE
        xl=r11(i)*x1(i)+r21(i)*y1(i)+r31(i)*z1(i)
        yl=r12(i)*x1(i)+r22(i)*y1(i)+r32(i)*z1(i)
        zl=r13(i)*x1(i)+r23(i)*y1(i)+r33(i)*z1(i)
        x1(i)=xl
        y1(i)=yl
        z1(i)=zl
        xl=r11(i)*x2(i)+r21(i)*y2(i)+r31(i)*z2(i)
        yl=r12(i)*x2(i)+r22(i)*y2(i)+r32(i)*z2(i)
        zl=r13(i)*x2(i)+r23(i)*y2(i)+r33(i)*z2(i)
        x2(i)=xl
        y2(i)=yl
        z2(i)=zl
        xl=r11(i)*x3(i)+r21(i)*y3(i)+r31(i)*z3(i)
        yl=r12(i)*x3(i)+r22(i)*y3(i)+r32(i)*z3(i)
        zl=r13(i)*x3(i)+r23(i)*y3(i)+r33(i)*z3(i)
        x3(i)=xl
        y3(i)=yl
        z3(i)=zl
        xl=r11(i)*x4(i)+r21(i)*y4(i)+r31(i)*z4(i)
        yl=r12(i)*x4(i)+r22(i)*y4(i)+r32(i)*z4(i)
C          ZL=R13(I)*X4(I)+R23(I)*Y4(I)+R33(I)*Z4(I)
        x4(i)=xl
        y4(i)=yl
        z4(i)=-z1(i)
        xl=r11(i)*x5(i)+r21(i)*y5(i)+r31(i)*z5(i)
        yl=r12(i)*x5(i)+r22(i)*y5(i)+r32(i)*z5(i)
C          ZL=R13(I)*X5(I)+R23(I)*Y5(I)+R33(I)*Z5(I)
        x5(i)=xl
        y5(i)=yl
        z5(i)=-z2(i)
        xl=r11(i)*x6(i)+r21(i)*y6(i)+r31(i)*z6(i)
        yl=r12(i)*x6(i)+r22(i)*y6(i)+r32(i)*z6(i)
C          ZL=R13(I)*X6(I)+R23(I)*Y6(I)+R33(I)*Z6(I)
        x6(i)=xl
        y6(i)=yl
        z6(i)=-z3(i)
C       sav in local frame.
C       SAV(I,1)=X1(I)-X6(I)
C       SAV(I,2)=Y1(I)-Y6(I)
C       SAV(I,3)=Z1(I)-Z6(I)
C       SAV(I,4)=X2(I)-X6(I)
C       SAV(I,5)=Y2(I)-Y6(I)
C       SAV(I,6)=Z2(I)-Z6(I)
C       SAV(I,7)=X3(I)-X6(I)
C       SAV(I,8)=Y3(I)-Y6(I)
C       SAV(I,9)=Z3(I)-Z6(I)
C       SAV(I,10)=X4(I)-X6(I)
C       SAV(I,11)=Y4(I)-Y6(I)
C       SAV(I,12)=Z4(I)-Z6(I)
C       SAV(I,13)=X5(I)-X6(I)
C       SAV(I,14)=Y5(I)-Y6(I)
C       SAV(I,15)=Z5(I)-Z6(I)
        off(i) = offg(i)
       ENDIF
       ENDDO
C
      ELSE
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     x1,      y1,      z1,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     x2,      y2,      z2,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     x3,      y3,      z3,
     4   nel)
       DO i=lft,llt
        xl=r11(i)*x4(i)+r21(i)*y4(i)+r31(i)*z4(i)
        yl=r12(i)*x4(i)+r22(i)*y4(i)+r32(i)*z4(i)
C          ZL=R13(I)*X4(I)+R23(I)*Y4(I)+R33(I)*Z4(I)
        x4(i)=xl
        y4(i)=yl
        z4(i)=-z1(i)
        xl=r11(i)*x5(i)+r21(i)*y5(i)+r31(i)*z5(i)
        yl=r12(i)*x5(i)+r22(i)*y5(i)+r32(i)*z5(i)
C          ZL=R13(I)*X5(I)+R23(I)*Y5(I)+R33(I)*Z5(I)
        x5(i)=xl
        y5(i)=yl
        z5(i)=-z2(i)
        xl=r11(i)*x6(i)+r21(i)*y6(i)+r31(i)*z6(i)
        yl=r12(i)*x6(i)+r22(i)*y6(i)+r32(i)*z6(i)
C          ZL=R13(I)*X6(I)+R23(I)*Y6(I)+R33(I)*Z6(I)
        x6(i)=xl
        y6(i)=yl
        z6(i)=-z3(i)
        off(i) = min(one,offg(i))
       ENDDO
C
       ENDIF
C-----------
C     transfer of velocities to the convected (or orthotropic) frame.
C-----------
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx1,     vy1,     vz1,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx2,     vy2,     vz2,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx3,     vy3,     vz3,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx4,     vy4,     vz4,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx5,     vy5,     vz5,
     4   nel)
      CALL vrrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx6,     vy6,     vz6,
     4   nel)
C-----------
      DO i=lft,llt
        vd2(i) = zero
      ENDDO
C-----PROJECTION----
      DO i=lft,llt
        xx = x1(i)*x1(i)+x2(i)*x2(i)+x3(i)*x3(i)
     1      +x4(i)*x4(i)+x5(i)*x5(i)+x6(i)*x6(i)
        yy = y1(i)*y1(i)+y2(i)*y2(i)+y3(i)*y3(i)
     1      +y4(i)*y4(i)+y5(i)*y5(i)+y6(i)*y6(i)
        xy = x1(i)*y1(i)+x2(i)*y2(i)+x3(i)*y3(i)
     1      +x4(i)*y4(i)+x5(i)*y5(i)+x6(i)*y6(i)
        xz = x1(i)*z1(i)+x2(i)*z2(i)+x3(i)*z3(i)
     1      +x4(i)*z4(i)+x5(i)*z5(i)+x6(i)*z6(i)
        yz = y1(i)*z1(i)+y2(i)*z2(i)+y3(i)*z3(i)
     1      +y4(i)*z4(i)+y5(i)*z5(i)+y6(i)*z6(i)
        zz = z1(i)*z1(i)+z2(i)*z2(i)+z3(i)*z3(i)
     1      +z4(i)*z4(i)+z5(i)*z5(i)+z6(i)*z6(i)
        rtr(1)= yy+zz
        rtr(2)= xx+zz
        rtr(3)= xx+yy
        rtr(4)= -xy
        rtr(5)= -xz
        rtr(6)= -yz
C       
        abc = rtr(1)*rtr(2)*rtr(3)
        xxyz2 = rtr(1)*rtr(6)*rtr(6)
        yyxz2 = rtr(2)*rtr(5)*rtr(5)
        zzxy2 = rtr(3)*rtr(4)*rtr(4)
        deta = abc + two*rtr(4)*rtr(5)*rtr(6)-xxyz2-yyxz2-zzxy2
        IF (deta<em20) THEN
          deta=one
        ELSE
         deta=one/deta
        ENDIF
        di(i,1) = (abc-xxyz2)*deta/rtr(1)
        di(i,2) = (abc-yyxz2)*deta/rtr(2)
        di(i,3) = (abc-zzxy2)*deta/rtr(3)
        di(i,4) = (rtr(5)*rtr(6)-rtr(4)*rtr(3))*deta
        di(i,5) = (rtr(6)*rtr(4)-rtr(5)*rtr(2))*deta
        di(i,6) = (rtr(4)*rtr(5)-rtr(6)*rtr(1))*deta
      ENDDO
      CALL s6proj3(
     1   x1,      x2,      x3,      x4,
     2   x5,      x6,      y1,      y2,
     3   y3,      y4,      y5,      y6,
     4   z1,      z2,      z3,      z4,
     5   z5,      z6,      vx1,     vx2,
     6   vx3,     vx4,     vx5,     vx6,
     7   vy1,     vy2,     vy3,     vy4,
     8   vy5,     vy6,     vz1,     vz2,
     9   vz3,     vz4,     vz5,     vz6,
     a   di,      nel)
      RETURN
      END