cbadef.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/.
!||====================================================================
!||    cbadef     ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
      SUBROUTINE cbadef(JFT,JLT,NG,VCORE,AREA,CDET,VQN,VQ,VJFI,
     1                  VXYZ,RXYZ,VDEF,VNRM,VASTN,
     2                  HX,HY,VETA,VKSI,BM,BMF,BF,BC,TC,NPLAT,IPLAT,
     3                  ISROT,BRZ )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C        CALCULATION OF GENERALIZED DEFORMATIONS AT GAUSS POINTS AND [B]
C        ENTREES :  NEL,NG,VCORE,VQN,VXYZ,RXYZ,VNRM,VASM,VASTN
C        SORTIES : VDEF,BM (MEMBRANE),BMF(COUPLAGE MEM-FLEXION),BF (FLEXION),BC(C.T)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
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   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER NPG,NG,JFT,JLT,NPLAT,IPLAT(*),ISROT
      PARAMETER (NPG = 4)
      my_real 
     .   rxyz(mvsiz,8),vcore(mvsiz,12),vxyz(mvsiz,12),
     .   vqn(mvsiz,9,4),vksi(4,4),veta(4,4),
     .   bm(mvsiz,36),bmf(mvsiz,36),bf(mvsiz,24),bc(mvsiz,40),hx(mvsiz,4),hy(mvsiz,4)
      my_real 
     .   vnrm(mvsiz,12),vastn(mvsiz,16),vjfi(mvsiz,3,2,4),
     .   vq(mvsiz,3,3,4),vdef(mvsiz,8),dt1,area(*)
      my_real 
     .   cdet(*),tc(mvsiz,2,2),brz(mvsiz,4,4)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP, IUN
      MY_REAL 
     .   VPG(2,NPG),PG1,PG,THK,DETJ,DET,
     .   tfn(3,2),bcx,bcy,bxy(3),byx(3),
     .   v1(2),v2(2),rv1,rv2,
     .   c1,c2,vt1,vt2,bc1,bc2,vb1,
     .   vjf1(2,3),vjf(3,3),tbi(2,2),tbc(2,2),v11(4)
      my_real 
     .   a_1,c11,c12,c21,c22,cc,beta1,ksi1,ksiy1,beta2,ksi2,ksiy2
       parameter(pg=.577350269189626)
       parameter(pg1=-.577350269189626)
C--------------------------
C     INITIALISATION
C--------------------------
      DATA iun/1/
      DATA vpg/pg1,pg1,pg,pg1,pg,pg,pg1,pg/
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include "vectorize.inc"
       DO i=jft,nplat 
        ep=iplat(i)
C-------VCORE(12):BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
C-------V13(I)->VXYZ(j,1,I),V24(I)->VXYZ(j,2,I),VHI(I)->VXYZ(j,3,I)------------
        bm(ep,1)=vcore(ep,1)+hx(ep,ng)*vcore(ep,5)
        bm(ep,2)=vcore(ep,2)+hx(ep,ng)*vcore(ep,6)
        bm(ep,3)=hx(ep,ng)*fourth
        bm(ep,4)=-bm(ep,3)
        bm(ep,5)=vcore(ep,3)+hy(ep,ng)*vcore(ep,5)
        bm(ep,6)=vcore(ep,4)+hy(ep,ng)*vcore(ep,6)
        bm(ep,7)=hy(ep,ng)*fourth
        bm(ep,8)=-bm(ep,7)
C
        vdef(ep,1)=bm(ep,1)*vxyz(ep,1)+bm(ep,2)*vxyz(ep,4)
     1          +bm(ep,3)*vxyz(ep,7)
        vdef(ep,2)=bm(ep,5)*vxyz(ep,2)+bm(ep,6)*vxyz(ep,5)
     1          +bm(ep,7)*vxyz(ep,8)
C--------FLEXION-----
        vdef(ep,6)=bm(ep,1)*rxyz(ep,2)+bm(ep,2)*rxyz(ep,4)
     1           +bm(ep,3)*rxyz(ep,6)
        vdef(ep,7)=-(bm(ep,5)*rxyz(ep,1)+bm(ep,6)*rxyz(ep,3)
     1           +bm(ep,7)*rxyz(ep,5))
        vdef(ep,8)=-(bm(ep,1)*rxyz(ep,1)+bm(ep,2)*rxyz(ep,3)
     1            +bm(ep,3)*rxyz(ep,5))
     2            +bm(ep,5)*rxyz(ep,2)+bm(ep,6)*rxyz(ep,4)
     3            +bm(ep,7)*rxyz(ep,6)
C--------C.T-----
C-------VCORE(12):BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
       a_1 = 0.25/max(cdet(ep),em20)
       c11=(vcore(ep,10)+vcore(ep,12)*vpg(1,ng))*a_1
       c12=-(vcore(ep,8)+vcore(ep,12)*vpg(2,ng))*a_1
       c21=-(vcore(ep,9)+vcore(ep,11)*vpg(1,ng))*a_1
       c22=( vcore(ep,7)+vcore(ep,11)*vpg(2,ng))*a_1
       beta1=vcore(ep,12)+vcore(ep,8)*vpg(2,ng)
        ksi1=vcore(ep,12)+vcore(ep,10)*vpg(1,ng)
       beta2=vcore(ep,11)+vcore(ep,7)*vpg(2,ng)
        ksi2=vcore(ep,11)+vcore(ep,9)*vpg(1,ng)
C-------NODE--- 1
       bc(ep,1)=-c11-c12
       bc(ep,2)=-c21-c22
       bc(ep,3)= beta1*c11+ksi1*c12
       bc(ep,4)= beta1*c21+ksi1*c22
       bc(ep,5)= -beta2*c11-ksi2*c12
       bc(ep,6)= -beta2*c21-ksi2*c22
C-------NODE--- 2
       bc(ep,7)= c11-c12
       bc(ep,8)= c21-c22
       bc(ep,9)= beta1*c11-ksi1*c12
       bc(ep,10)=beta1*c21-ksi1*c22
       bc(ep,11)=-beta2*c11+ksi2*c12
       bc(ep,12)=-beta2*c21+ksi2*c22
C
       beta1=vcore(ep,8)+vcore(ep,12)*vpg(2,ng)
        ksi1=vcore(ep,10)+vcore(ep,12)*vpg(1,ng)
       beta2=vcore(ep,7)+vcore(ep,11)*vpg(2,ng)
        ksi2=vcore(ep,9)+vcore(ep,11)*vpg(1,ng)
C-----PART SYM---------
C-------NODE--- 1
       bc(ep,13)=c11*vpg(2,ng)+c12*vpg(1,ng)
       bc(ep,14)=c21*vpg(2,ng)+c22*vpg(1,ng)
       bc(ep,15)=-beta1*c11-ksi1*c12
       bc(ep,16)=-beta1*c21-ksi1*c22
       bc(ep,17)=c11*beta2+c12*ksi2
       bc(ep,18)=c21*beta2+c22*ksi2
C-------NODE--- 2
       bc(ep,19)=-bc(ep,13)
       bc(ep,20)=-bc(ep,14)
       bc(ep,21)=bc(ep,15)
       bc(ep,22)=bc(ep,16)
       bc(ep,23)=bc(ep,17)
       bc(ep,24)=bc(ep,18)
C
        vdef(ep,4)=bc(ep,1)*vxyz(ep,3)+bc(ep,7)*vxyz(ep,6)
     1            +bc(ep,13)*vxyz(ep,9)+bc(ep,3)*rxyz(ep,1)
     2            +bc(ep,9)*rxyz(ep,3)+bc(ep,15)*rxyz(ep,7)
     3            +bc(ep,5)*rxyz(ep,2)+bc(ep,11)*rxyz(ep,4)
     4            +bc(ep,17)*rxyz(ep,8)
        vdef(ep,5)=bc(ep,2)*vxyz(ep,3)+bc(ep,8)*vxyz(ep,6)
     1            +bc(ep,14)*vxyz(ep,9)+bc(ep,4)*rxyz(ep,1)
     2            +bc(ep,10)*rxyz(ep,3)+bc(ep,16)*rxyz(ep,7)
     3            +bc(ep,6)*rxyz(ep,2)+bc(ep,12)*rxyz(ep,4)
     4            +bc(ep,18)*rxyz(ep,8)
       ENDDO
       IF (isrot>0) THEN
#include "vectorize.inc"
        DO i=jft,nplat
         ep=iplat(i)
C----------no more constant shear-----------
         vdef(ep,3)=bm(ep,1)*vxyz(ep,2)+bm(ep,2)*vxyz(ep,5)
     1             +bm(ep,3)*vxyz(ep,8)
     1             +bm(ep,5)*vxyz(ep,1)+bm(ep,6)*vxyz(ep,4)
     1             +bm(ep,7)*vxyz(ep,7)
        END DO
       END IF !(ISROT>0) THEN
C------------ELEMENT GAUCH-------------------
#include "vectorize.inc"
       DO 150 i=nplat+1,jlt 
        ep=iplat(i)
C---------------------------------------------------
C  CALCULATION OF [FN]
C---------------------------------------------------
          tfn(1,1)=vksi(1,ng)*vqn(ep,7,1)+vksi(2,ng)*vqn(ep,7,2)
     1            +vksi(3,ng)*vqn(ep,7,3)+vksi(4,ng)*vqn(ep,7,4)
          tfn(2,1)=vksi(1,ng)*vqn(ep,8,1)+vksi(2,ng)*vqn(ep,8,2)
     1            +vksi(3,ng)*vqn(ep,8,3)+vksi(4,ng)*vqn(ep,8,4)
          tfn(3,1)=vksi(1,ng)*vqn(ep,9,1)+vksi(2,ng)*vqn(ep,9,2)
     1            +vksi(3,ng)*vqn(ep,9,3)+vksi(4,ng)*vqn(ep,9,4)
          tfn(1,2)=veta(1,ng)*vqn(ep,7,1)+veta(2,ng)*vqn(ep,7,2)
     1            +veta(3,ng)*vqn(ep,7,3)+veta(4,ng)*vqn(ep,7,4)
          tfn(2,2)=veta(1,ng)*vqn(ep,8,1)+veta(2,ng)*vqn(ep,8,2)
     1            +veta(3,ng)*vqn(ep,8,3)+veta(4,ng)*vqn(ep,8,4)
          tfn(3,2)=veta(1,ng)*vqn(ep,9,1)+veta(2,ng)*vqn(ep,9,2)
     1            +veta(3,ng)*vqn(ep,9,3)+veta(4,ng)*vqn(ep,9,4)
C--------------------------------------------------------
C  CALCUL [BN]=[F]^-1 [FN]
C--------------------------------------------------------
          tbi(2,2)=vjfi(ep,1,1,ng)*tfn(1,1)+vjfi(ep,2,1,ng)*tfn(2,1)
     1           + vjfi(ep,3,1,ng)*tfn(3,1)
          tbi(2,1)=vjfi(ep,1,2,ng)*tfn(1,1)+vjfi(ep,2,2,ng)*tfn(2,1)
     1           + vjfi(ep,3,2,ng)*tfn(3,1)
          tbi(1,2)=vjfi(ep,1,1,ng)*tfn(1,2)+vjfi(ep,2,1,ng)*tfn(2,2)
     1           + vjfi(ep,3,1,ng)*tfn(3,2)
          tbi(1,1)=vjfi(ep,1,2,ng)*tfn(1,2)+vjfi(ep,2,2,ng)*tfn(2,2)
     1           + vjfi(ep,3,2,ng)*tfn(3,2)
C
        thk =-(tbi(1,1)+tbi(2,2))
C       THK(2)=TBI(1,1)*TBI(2,2)-TBI(1,2)*TBI(2,1)
        tbi(1,2)=-tbi(1,2)
        tbi(2,1)=-tbi(2,1)
C---------------------------
C  CALCUL DE TC=VJFI*VQ
C---------------------------
          tc(ep,1,1)=vjfi(ep,1,1,ng)*vq(ep,1,1,ng)+vjfi(ep,2,1,ng)
     +             *vq(ep,2,1,ng)+ vjfi(ep,3,1,ng)*vq(ep,3,1,ng)
          tc(ep,2,1)=vjfi(ep,1,2,ng)*vq(ep,1,1,ng)+vjfi(ep,2,2,ng)
     +             *vq(ep,2,1,ng)+ vjfi(ep,3,2,ng)*vq(ep,3,1,ng)
          tc(ep,1,2)=vjfi(ep,1,1,ng)*vq(ep,1,2,ng)+vjfi(ep,2,1,ng)
     +             *vq(ep,2,2,ng)+ vjfi(ep,3,1,ng)*vq(ep,3,2,ng)
          tc(ep,2,2)=vjfi(ep,1,2,ng)*vq(ep,1,2,ng)+vjfi(ep,2,2,ng)
     +             *vq(ep,2,2,ng)+ vjfi(ep,3,2,ng)*vq(ep,3,2,ng)
C------------------------------
C  CALCUL DE TBC=TBI*TC
C------------------------------
          tbc(1,1)=tbi(1,1)*tc(ep,1,1)+tbi(1,2)*tc(ep,2,1)
          tbc(2,1)=tbi(2,1)*tc(ep,1,1)+tbi(2,2)*tc(ep,2,1)
          tbc(1,2)=tbi(1,1)*tc(ep,1,2)+tbi(1,2)*tc(ep,2,2)
          tbc(2,2)=tbi(2,1)*tc(ep,1,2)+tbi(2,2)*tc(ep,2,2)
C--------------------------
C     BUILD [B] MATRIX AND DEFORMATION
C--------------------------
C--------------------------
C     DEFORMATION IN-PLANE
C--------------------------
C      [BM](3,3*NPG),[BMF](3,3*NPG),
C      [BF](3,2*NPG),VDEF(8), 
C--------------------------
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C---------- LOOP OVER NODES---
C--------J=1         II =0 KK = 0
         vt1=vq(ep,1,1,ng)*vxyz(ep,1)+vq(ep,2,1,ng)*vxyz(ep,2)
     1        +vq(ep,3,1,ng)*vxyz(ep,3)
         vt2=vq(ep,1,2,ng)*vxyz(ep,1)+vq(ep,2,2,ng)*vxyz(ep,2)
     1        +vq(ep,3,2,ng)*vxyz(ep,3)
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(1,ng)*tc(ep,1,1)+veta(1,ng)*tc(ep,2,1)
         c2=vksi(1,ng)*tc(ep,1,2)+veta(1,ng)*tc(ep,2,2)
C--------TERMES ASSOCIEES AUX U,V,W DE [B1]=<T>BCI-------
         bc1=vksi(1,ng)*tbc(1,1)+veta(1,ng)*tbc(2,1)
         bc2=vksi(1,ng)*tbc(1,2)+veta(1,ng)*tbc(2,2)
C          I = 1; JJ = 3*KK =0
         IF (isrot>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          brz(i,1,1)=-bxy(1)+byx(1)
          brz(i,2,1)=-bxy(2)+byx(2)
          brz(i,3,1)=-bxy(3)+byx(3)
          bm(ep,3)=bxy(1)+byx(1)
          bm(ep,6)=bxy(2)+byx(2)
          bm(ep,9)=bxy(3)+byx(3)
          vdef(ep,3)= c1*vt2+c2*vt1
         ELSE
          bm(ep,9)=zero
         END IF
C---------CONTRIBUTION DU <T1> :EPS-X
          bm(ep,1)=vq(ep,1,1,ng)*c1
          bmf(ep,1)=thk*bm(ep,1)+vq(ep,1,1,ng)*bc1
C---------CONTRIBUTION DU <T2> :EPS-Y
          bm(ep,2)=vq(ep,1,2,ng)*c2
          bmf(ep,2)=thk*bm(ep,2)+vq(ep,1,2,ng)*bc2
C---------CONTRIBUTION DU <T> :EPS-XY
          bmf(ep,3)=thk*bm(ep,3)+vq(ep,1,1,ng)*bc2+vq(ep,1,2,ng)*bc1
C          I = 2; JJ = 3*(KK+1)=3
C---------:EPS-X :EPS-Y :EPS-XY
          bm(ep,4)=vq(ep,2,1,ng)*c1
          bmf(ep,4)=thk*bm(ep,4)+vq(ep,2,1,ng)*bc1
          bm(ep,5)=vq(ep,2,2,ng)*c2
          bmf(ep,5)=thk*bm(ep,5)+vq(ep,2,2,ng)*bc2
          bmf(ep,6)=thk*bm(ep,6)+vq(ep,2,1,ng)*bc2+vq(ep,2,2,ng)*bc1
C          I = 3; JJ = 3*(KK+2)=6
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,7)=vq(ep,3,1,ng)*c1
          bmf(ep,7)=thk*bm(ep,7)+vq(ep,3,1,ng)*bc1
          bm(ep,8)=vq(ep,3,2,ng)*c2
          bmf(ep,8)=thk*bm(ep,8)+vq(ep,3,2,ng)*bc2
          bmf(ep,9)=thk*bm(ep,9)+vq(ep,3,1,ng)*bc2+vq(ep,3,2,ng)*bc1
C
C---------DEF MEMBRANE------
          vdef(ep,1)= c1*vt1
          vdef(ep,2)= c2*vt2
C
C---------FLEXION ASSOCIEE AUX U,V,W ([BMF]=2H[B0]+[B01])------
C----------[B01]---------
          vdef(ep,6)= bc1*vt1
          vdef(ep,7)= bc2*vt2
          vdef(ep,8)= bc1*vt2+bc2*vt1
C---------- TERMES ASSOCIEES AU BETA ----------
C---------PRODUIT SCALAIRE <T(PG)>*<-T^2 T^1>----------
         v1(2)=vq(ep,1,1,ng)*vqn(ep,1,1)+vq(ep,2,1,ng)*vqn(ep,2,1)
     1           +vq(ep,3,1,ng)*vqn(ep,3,1)
         v1(1)=-(vq(ep,1,1,ng)*vqn(ep,4,1)+vq(ep,2,1,ng)*vqn(ep,5,1)
     1           +vq(ep,3,1,ng)*vqn(ep,6,1))
         v2(2)=vq(ep,1,2,ng)*vqn(ep,1,1)+vq(ep,2,2,ng)*vqn(ep,2,1)
     1           +vq(ep,3,2,ng)*vqn(ep,3,1)
         v2(1)=-(vq(ep,1,2,ng)*vqn(ep,4,1)+vq(ep,2,2,ng)*vqn(ep,5,1)
     1           +vq(ep,3,2,ng)*vqn(ep,6,1))
          rv1=v1(1)*rxyz(ep,1)+v1(2)*rxyz(ep,2)
          rv2=v2(1)*rxyz(ep,1)+v2(2)*rxyz(ep,2)
C----------TERMES ASSOCIEES AUX BETA DE [B01]--
C           JJ = 3*II =0
            bf(ep,1)=v1(1)*c1
            bf(ep,2)=v2(1)*c2
            bf(ep,3)=v1(1)*c2+v2(1)*c1
            bf(ep,4)=v1(2)*c1
            bf(ep,5)=v2(2)*c2
            bf(ep,6)=v1(2)*c2+v2(2)*c1
C---------DEF DE FLEXION ASSOCIEE AUX THETA1,THETA2------
          vdef(ep,6)= vdef(ep,6)+c1*rv1
          vdef(ep,7)= vdef(ep,7)+c2*rv2
          vdef(ep,8)= vdef(ep,8)+c1*rv2+c2*rv1
C
C--------J=2---II=(J-1)*2 =2  KK = 3*(J-1) =3
         vt1=vq(ep,1,1,ng)*vxyz(ep,4)+vq(ep,2,1,ng)*vxyz(ep,5)
     1        +vq(ep,3,1,ng)*vxyz(ep,6)
         vt2=vq(ep,1,2,ng)*vxyz(ep,4)+vq(ep,2,2,ng)*vxyz(ep,5)
     1        +vq(ep,3,2,ng)*vxyz(ep,6)
C--------- [B0], [B01]---------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(2,ng)*tc(ep,1,1)+veta(2,ng)*tc(ep,2,1)
         c2=vksi(2,ng)*tc(ep,1,2)+veta(2,ng)*tc(ep,2,2)
C--------TERMES ASSOCIEES AUX U,V,W DE [B1]=<T>BCI-------
         bc1=vksi(2,ng)*tbc(1,1)+veta(2,ng)*tbc(2,1)
         bc2=vksi(2,ng)*tbc(1,2)+veta(2,ng)*tbc(2,2)
C
C         I = 1; JJ = 3*KK =9
         IF (isrot>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          brz(i,1,2)=-bxy(1)+byx(1)
          brz(i,2,2)=-bxy(2)+byx(2)
          brz(i,3,2)=-bxy(3)+byx(3)
          bm(ep,12)=bxy(1)+byx(1)
          bm(ep,15)=bxy(2)+byx(2)
          bm(ep,18)=bxy(3)+byx(3)
          vdef(ep,3)= vdef(ep,3)+ c1*vt2+c2*vt1
         ELSE
          bm(ep,18)=zero
         END IF 
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,10)=vq(ep,1,1,ng)*c1
          bmf(ep,10)=thk*bm(ep,10)+vq(ep,1,1,ng)*bc1
          bm(ep,11)=vq(ep,1,2,ng)*c2
          bmf(ep,11)=thk*bm(ep,11)+vq(ep,1,2,ng)*bc2
          bmf(ep,12)=thk*bm(ep,12)+vq(ep,1,1,ng)*bc2+vq(ep,1,2,ng)*bc1
C
C         I = 2; JJ = 3*(KK+1)=12
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,13)=vq(ep,2,1,ng)*c1
          bmf(ep,13)=thk*bm(ep,13)+vq(ep,2,1,ng)*bc1
          bm(ep,14)=vq(ep,2,2,ng)*c2
          bmf(ep,14)=thk*bm(ep,14)+vq(ep,2,2,ng)*bc2
          bmf(ep,15)=thk*bm(ep,15)+vq(ep,2,1,ng)*bc2+vq(ep,2,2,ng)*bc1
C
C         I = 3; JJ = 3*(KK+2)=15
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,16)=vq(ep,3,1,ng)*c1
          bmf(ep,16)=thk*bm(ep,16)+vq(ep,3,1,ng)*bc1
          bm(ep,17)=vq(ep,3,2,ng)*c2
          bmf(ep,17)=thk*bm(ep,17)+vq(ep,3,2,ng)*bc2
          bmf(ep,18)=thk*bm(ep,18)+vq(ep,3,1,ng)*bc2+vq(ep,3,2,ng)*bc1
C
C---------DEF MEMBRANE------
          vdef(ep,1)= vdef(ep,1)+ c1*vt1
          vdef(ep,2)= vdef(ep,2)+ c2*vt2
C---------DEF DE FLEXION ASSOCIEE AUX U,V,W ([BMF]=2H[B0]+[B01])------
C----------[B01]---------
          vdef(ep,6)= vdef(ep,6)+ bc1*vt1
          vdef(ep,7)= vdef(ep,7)+ bc2*vt2
          vdef(ep,8)= vdef(ep,8)+ bc1*vt2+bc2*vt1
C---------- TERMES ASSOCIEES AU BETA ----------
C---------PRODUIT SCALAIRE <T(PG)>*<-T^2 T^1>----------
         v1(2)=vq(ep,1,1,ng)*vqn(ep,1,2)+vq(ep,2,1,ng)*vqn(ep,2,2)
     1           +vq(ep,3,1,ng)*vqn(ep,3,2)
         v1(1)=-(vq(ep,1,1,ng)*vqn(ep,4,2)+vq(ep,2,1,ng)*vqn(ep,5,2)
     1           +vq(ep,3,1,ng)*vqn(ep,6,2))
         v2(2)=vq(ep,1,2,ng)*vqn(ep,1,2)+vq(ep,2,2,ng)*vqn(ep,2,2)
     1           +vq(ep,3,2,ng)*vqn(ep,3,2)
         v2(1)=-(vq(ep,1,2,ng)*vqn(ep,4,2)+vq(ep,2,2,ng)*vqn(ep,5,2)
     1           +vq(ep,3,2,ng)*vqn(ep,6,2))
          rv1=v1(1)*rxyz(ep,2+1)+v1(2)*rxyz(ep,2+2)
          rv2=v2(1)*rxyz(ep,2+1)+v2(2)*rxyz(ep,2+2)
C----------TERMES ASSOCIEES AUX BETA DE [B01]--
C           JJ = 3*II=6
            bf(ep,7)=v1(1)*c1
            bf(ep,8)=v2(1)*c2
            bf(ep,9)=v1(1)*c2+v2(1)*c1
            bf(ep,10)=v1(2)*c1
            bf(ep,11)=v2(2)*c2
            bf(ep,12)=v1(2)*c2+v2(2)*c1
C---------DEF DE FLEXION ASSOCIEE AUX THETA1,THETA2------
          vdef(ep,6)= vdef(ep,6)+c1*rv1
          vdef(ep,7)= vdef(ep,7)+c2*rv2
          vdef(ep,8)= vdef(ep,8)+c1*rv2+c2*rv1
C
C--------J=3---II=(J-1)*2 =4 KK = 3*(J-1)=6
         vt1=vq(ep,1,1,ng)*vxyz(ep,7)+vq(ep,2,1,ng)*vxyz(ep,8)
     1        +vq(ep,3,1,ng)*vxyz(ep,9)
         vt2=vq(ep,1,2,ng)*vxyz(ep,7)+vq(ep,2,2,ng)*vxyz(ep,8)
     1        +vq(ep,3,2,ng)*vxyz(ep,9)
C--------- [B0], [B01]---------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(3,ng)*tc(ep,1,1)+veta(3,ng)*tc(ep,2,1)
         c2=vksi(3,ng)*tc(ep,1,2)+veta(3,ng)*tc(ep,2,2)
C--------TERMES ASSOCIEES AUX U,V,W DE [B1]=<T>BCI-------
         bc1=vksi(3,ng)*tbc(1,1)+veta(3,ng)*tbc(2,1)
         bc2=vksi(3,ng)*tbc(1,2)+veta(3,ng)*tbc(2,2)
C         I = 1;JJ = 3*KK=18
         IF (isrot>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          brz(i,1,3)=-bxy(1)+byx(1)
          brz(i,2,3)=-bxy(2)+byx(2)
          brz(i,3,3)=-bxy(3)+byx(3)
          bm(ep,21)=bxy(1)+byx(1)
          bm(ep,24)=bxy(2)+byx(2)
          bm(ep,27)=bxy(3)+byx(3)
          vdef(ep,3)= vdef(ep,3)+ c1*vt2+c2*vt1
         ELSE
          bm(ep,27)=zero
         END IF 
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,19)=vq(ep,1,1,ng)*c1
          bmf(ep,19)=thk*bm(ep,19)+vq(ep,1,1,ng)*bc1
          bm(ep,20)=vq(ep,1,2,ng)*c2
          bmf(ep,20)=thk*bm(ep,20)+vq(ep,1,2,ng)*bc2
          bmf(ep,21)=thk*bm(ep,21)+vq(ep,1,1,ng)*bc2+vq(ep,1,2,ng)*bc1
C
C         I = 2; JJ = 3*(KK+1)=21
          bm(ep,22)=vq(ep,2,1,ng)*c1
          bmf(ep,22)=thk*bm(ep,22)+vq(ep,2,1,ng)*bc1
          bm(ep,23)=vq(ep,2,2,ng)*c2
          bmf(ep,23)=thk*bm(ep,23)+vq(ep,2,2,ng)*bc2
          bmf(ep,24)=thk*bm(ep,24)+vq(ep,2,1,ng)*bc2+vq(ep,2,2,ng)*bc1
C
C         I = 3; JJ = 3*(KK+2)=24
          bm(ep,25)=vq(ep,3,1,ng)*c1
          bmf(ep,25)=thk*bm(ep,25)+vq(ep,3,1,ng)*bc1
          bm(ep,26)=vq(ep,3,2,ng)*c2
          bmf(ep,26)=thk*bm(ep,26)+vq(ep,3,2,ng)*bc2
          bmf(ep,27)=thk*bm(ep,27)+vq(ep,3,1,ng)*bc2+vq(ep,3,2,ng)*bc1
C
C---------DEF MEMBRANE------
          vdef(ep,1)= vdef(ep,1)+ c1*vt1
          vdef(ep,2)= vdef(ep,2)+ c2*vt2
C---------DEF DE FLEXION ASSOCIEE AUX U,V,W ([BMF]=2H[B0]+[B01])------
C----------[B01]---------
          vdef(ep,6)= vdef(ep,6)+ bc1*vt1
          vdef(ep,7)= vdef(ep,7)+ bc2*vt2
          vdef(ep,8)= vdef(ep,8)+ bc1*vt2+bc2*vt1
         v1(2)=vq(ep,1,1,ng)*vqn(ep,1,3)+vq(ep,2,1,ng)*vqn(ep,2,3)
     1           +vq(ep,3,1,ng)*vqn(ep,3,3)
         v1(1)=-(vq(ep,1,1,ng)*vqn(ep,4,3)+vq(ep,2,1,ng)*vqn(ep,5,3)
     1           +vq(ep,3,1,ng)*vqn(ep,6,3))
         v2(2)=vq(ep,1,2,ng)*vqn(ep,1,3)+vq(ep,2,2,ng)*vqn(ep,2,3)
     1           +vq(ep,3,2,ng)*vqn(ep,3,3)
         v2(1)=-(vq(ep,1,2,ng)*vqn(ep,4,3)+vq(ep,2,2,ng)*vqn(ep,5,3)
     1           +vq(ep,3,2,ng)*vqn(ep,6,3))
          rv1=v1(1)*rxyz(ep,4+1)+v1(2)*rxyz(ep,4+2)
          rv2=v2(1)*rxyz(ep,4+1)+v2(2)*rxyz(ep,4+2)
            bf(ep,13)=v1(1)*c1
            bf(ep,14)=v2(1)*c2
            bf(ep,15)=v1(1)*c2+v2(1)*c1
            bf(ep,16)=v1(2)*c1
            bf(ep,17)=v2(2)*c2
            bf(ep,18)=v1(2)*c2+v2(2)*c1
          vdef(ep,6)= vdef(ep,6)+c1*rv1
          vdef(ep,7)= vdef(ep,7)+c2*rv2
          vdef(ep,8)= vdef(ep,8)+c1*rv2+c2*rv1
C--------J=4---II=(J-1)*2  =6 KK = 3*(J-1)=9
         vt1=vq(ep,1,1,ng)*vxyz(ep,10)+vq(ep,2,1,ng)*vxyz(ep,11)
     1        +vq(ep,3,1,ng)*vxyz(ep,12)
         vt2=vq(ep,1,2,ng)*vxyz(ep,10)+vq(ep,2,2,ng)*vxyz(ep,11)
     1        +vq(ep,3,2,ng)*vxyz(ep,12)
C--------- [B0], [B01]---------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(4,ng)*tc(ep,1,1)+veta(4,ng)*tc(ep,2,1)
         c2=vksi(4,ng)*tc(ep,1,2)+veta(4,ng)*tc(ep,2,2)
C--------TERMES ASSOCIEES AUX U,V,W DE [B1]=<T>BCI-------
         bc1=vksi(4,ng)*tbc(1,1)+veta(4,ng)*tbc(2,1)
         bc2=vksi(4,ng)*tbc(1,2)+veta(4,ng)*tbc(2,2)
C         I = 1; JJ = 3*KK =27
         IF (isrot>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          brz(i,1,4)=-bxy(1)+byx(1)
          brz(i,2,4)=-bxy(2)+byx(2)
          brz(i,3,4)=-bxy(3)+byx(3)
          bm(ep,30)=bxy(1)+byx(1)
          bm(ep,33)=bxy(2)+byx(2)
          bm(ep,36)=bxy(3)+byx(3)
          vdef(ep,3)= vdef(ep,3)+ c1*vt2+c2*vt1
         ELSE
          bm(ep,36)=zero
         END IF 
          bm(ep,28)=vq(ep,1,1,ng)*c1
          bmf(ep,28)=thk*bm(ep,28)+vq(ep,1,1,ng)*bc1
          bm(ep,29)=vq(ep,1,2,ng)*c2
          bmf(ep,29)=thk*bm(ep,29)+vq(ep,1,2,ng)*bc2
C---------CONTRIBUTION DU <T> :EPS-XY
          bmf(ep,30)=thk*bm(ep,30)+vq(ep,1,1,ng)*bc2+vq(ep,1,2,ng)*bc1
C
C         I = 2; JJ = 3*(KK+1)=30
          bm(ep,31)=vq(ep,2,1,ng)*c1
          bmf(ep,31)=thk*bm(ep,31)+vq(ep,2,1,ng)*bc1
          bm(ep,32)=vq(ep,2,2,ng)*c2
          bmf(ep,32)=thk*bm(ep,32)+vq(ep,2,2,ng)*bc2
          bmf(ep,33)=thk*bm(ep,33)+vq(ep,2,1,ng)*bc2+vq(ep,2,2,ng)*bc1
C
C         I = 3; JJ = 3*(KK+2)=33
          bm(ep,34)=vq(ep,3,1,ng)*c1
          bmf(ep,34)=thk*bm(ep,34)+vq(ep,3,1,ng)*bc1
          bm(ep,35)=vq(ep,3,2,ng)*c2
          bmf(ep,35)=thk*bm(ep,35)+vq(ep,3,2,ng)*bc2
          bmf(ep,36)=thk*bm(ep,36)+vq(ep,3,1,ng)*bc2+vq(ep,3,2,ng)*bc1
C
C---------DEF MEMBRANE------
          vdef(ep,1)= vdef(ep,1)+ c1*vt1
          vdef(ep,2)= vdef(ep,2)+ c2*vt2
C---------DEF DE FLEXION ASSOCIEE AUX U,V,W ([BMF]=2H[B0]+[B01])------
C----------ONLY THE [B01] PART---------
          vdef(ep,6)= vdef(ep,6)+ bc1*vt1
          vdef(ep,7)= vdef(ep,7)+ bc2*vt2
          vdef(ep,8)= vdef(ep,8)+ bc1*vt2+bc2*vt1
C---------- TERMES ASSOCIEES AU BETA ----------
C---------PRODUIT SCALAIRE <T(PG)>*<-T^2 T^1>----------
         v1(2)=vq(ep,1,1,ng)*vqn(ep,1,4)+vq(ep,2,1,ng)*vqn(ep,2,4)
     1           +vq(ep,3,1,ng)*vqn(ep,3,4)
         v1(1)=-(vq(ep,1,1,ng)*vqn(ep,4,4)+vq(ep,2,1,ng)*vqn(ep,5,4)
     1           +vq(ep,3,1,ng)*vqn(ep,6,4))
         v2(2)=vq(ep,1,2,ng)*vqn(ep,1,4)+vq(ep,2,2,ng)*vqn(ep,2,4)
     1           +vq(ep,3,2,ng)*vqn(ep,3,4)
         v2(1)=-(vq(ep,1,2,ng)*vqn(ep,4,4)+vq(ep,2,2,ng)*vqn(ep,5,4)
     1           +vq(ep,3,2,ng)*vqn(ep,6,4))
          rv1=v1(1)*rxyz(ep,6+1)+v1(2)*rxyz(ep,6+2)
          rv2=v2(1)*rxyz(ep,6+1)+v2(2)*rxyz(ep,6+2)
C           JJ = 3*II=18
            bf(ep,19)=v1(1)*c1
            bf(ep,20)=v2(1)*c2
            bf(ep,21)=v1(1)*c2+v2(1)*c1
            bf(ep,22)=v1(2)*c1
            bf(ep,23)=v2(2)*c2
            bf(ep,24)=v1(2)*c2+v2(2)*c1
          vdef(ep,6)= vdef(ep,6)+thk*vdef(ep,1)+c1*rv1
          vdef(ep,7)= vdef(ep,7)+thk*vdef(ep,2)+c2*rv2
          vdef(ep,8)= vdef(ep,8)+thk*vdef(ep,3)+c1*rv2+c2*rv1
C
C--------------------------
C     DEFORMATION OUT-PLANE (C.T)
C--------------------------
C      [BC](2,5*NPG),VDEF(4,5) 
C--------------------------
C
        v11(1)=vksi(2,ng)
        v11(2)=vksi(3,ng)
        v11(3)=veta(4,ng)
        v11(4)=veta(3,ng)
C--- J=1---------------------
C-----------POINTEURS POUR BC(II,J1~J2)
C         J1=(KA1(1,J)-1)*5=0
C         J2=(KA1(2,J)-1)*5=5
C         II=KIND(J)       =1
C         JJ = 3*(J-1)    =0
C--------TERMES <NAI>-------
          c1=v11(1)*vnrm(ep,1)
          bc(ep,1)=-c1
          bc(ep,11)=c1
C
          c1=v11(1)*vnrm(ep,2)
          bc(ep,3)=-c1
          bc(ep,13)=c1
C
          c1=v11(1)*vnrm(ep,3)
          bc(ep,5)=-c1
          bc(ep,15)=c1
C--------TERMES <ASI>-------
          bc(ep,7)=v11(1)*vastn(ep,1)
          bc(ep,9)=v11(1)*vastn(ep,2)
          bc(ep,17)=v11(1)*vastn(ep,3)
          bc(ep,19)=v11(1)*vastn(ep,4)
C--- J=2---------------------
C-----------POINTEURS POUR BC(II,J1~J2)
C         J1=(KA1(1,J)-1)*5 =15
C         J2=(KA1(2,J)-1)*5 =10
C         II=KIND(J)        =1
C         JJ = 3*(J-1)     =3
C--------TERMES <NAI>-------
          c1=v11(2)*vnrm(ep,4)
          bc(ep,31)=-c1
          bc(ep,21)=c1
C
          c1=v11(2)*vnrm(ep,5)
          bc(ep,33)=-c1
          bc(ep,23)=c1
C
          c1=v11(2)*vnrm(ep,6)
          bc(ep,35)=-c1
          bc(ep,25)=c1
C--------TERMES <ASI>-------
          bc(ep,37)=v11(2)*vastn(ep,5)
          bc(ep,39)=v11(2)*vastn(ep,6)
          bc(ep,27)=v11(2)*vastn(ep,7)
          bc(ep,29)=v11(2)*vastn(ep,8)
C--- J=3---------------------
C-----------POINTEURS POUR BC(II,J1~J2)
C         J1=(KA1(1,J)-1)*5 =0
C         J2=(KA1(2,J)-1)*5 =15
C         II=KIND(J)        =2
C         JJ = 3*(J-1)     =6
C--------TERMES <NAI>-------
          c1=v11(3)*vnrm(ep,7)
          bc(ep,2)=-c1
          bc(ep,32)=c1
C
          c1=v11(3)*vnrm(ep,8)
          bc(ep,4)=-c1
          bc(ep,34)=c1
C
          c1=v11(3)*vnrm(ep,9)
          bc(ep,6)=-c1
          bc(ep,36)=c1
C--------TERMES <ASI>-------
          bc(ep,8)=v11(3)*vastn(ep,9)
          bc(ep,10)=v11(3)*vastn(ep,10)
          bc(ep,38)=v11(3)*vastn(ep,11)
          bc(ep,40)=v11(3)*vastn(ep,12)
C--- J=4---------------------
C-----------POINTEURS POUR BC(II,J1~J2)
C         J1=(KA1(1,J)-1)*5 =5
C         J2=(KA1(2,J)-1)*5 =10
C         II=KIND(J)        =2
C         JJ = 3*(J-1)     =9
C--------TERMES <NAI>-------
          c1=v11(4)*vnrm(ep,10)
          bc(ep,12)=-c1
          bc(ep,22)=c1
C
          c1=v11(4)*vnrm(ep,11)
          bc(ep,14)=-c1
          bc(ep,24)=c1
C
          c1=v11(4)*vnrm(ep,12)
          bc(ep,16)=-c1
          bc(ep,26)=c1
C--------TERMES <ASI>-------
          bc(ep,18)=v11(4)*vastn(ep,13)
          bc(ep,20)=v11(4)*vastn(ep,14)
          bc(ep,28)=v11(4)*vastn(ep,15)
          bc(ep,30)=v11(4)*vastn(ep,16)
C
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C----------J=1--II=(J-1)*3 =0 JJ=(J-1)*2 =0 KK=10*(J-1)=0
         bcx=
     1            bc(ep,1)*vxyz(ep,1)+bc(ep,3)*vxyz(ep,2)
     2           +bc(ep,5)*vxyz(ep,3)+bc(ep,7)*rxyz(ep,1)
     3           +bc(ep,9)*rxyz(ep,2)       
         bcy=
     1            bc(ep,2)*vxyz(ep,1)+bc(ep,4)*vxyz(ep,2)
     2           +bc(ep,6)*vxyz(ep,3)+bc(ep,8)*rxyz(ep,1)
     3           +bc(ep,10)*rxyz(ep,2)       
C----------J=2 II=(J-1)*3 =3 JJ=(J-1)*2 =2 KK=10*(J-1) =10
         bcx=bcx
     1           +bc(ep,11)*vxyz(ep,4)+bc(ep,13)*vxyz(ep,5)
     2           +bc(ep,15)*vxyz(ep,6)+bc(ep,17)*rxyz(ep,3)
     3           +bc(ep,19)*rxyz(ep,4)       
         bcy=bcy
     1           +bc(ep,12)*vxyz(ep,4)+bc(ep,14)*vxyz(ep,5)
     2           +bc(ep,16)*vxyz(ep,6)+bc(ep,18)*rxyz(ep,3)
     3           +bc(ep,20)*rxyz(ep,4)       
C----------J=3--  II=(J-1)*3=6 JJ=(J-1)*2=4  KK=10*(J-1)=20
         bcx=bcx
     1           +bc(ep,21)*vxyz(ep,7)+bc(ep,23)*vxyz(ep,8)
     2           +bc(ep,25)*vxyz(ep,9)+bc(ep,27)*rxyz(ep,5)
     3           +bc(ep,29)*rxyz(ep,6)       
         bcy=bcy
     1           +bc(ep,22)*vxyz(ep,7)+bc(ep,24)*vxyz(ep,8)
     2           +bc(ep,26)*vxyz(ep,9)+bc(ep,28)*rxyz(ep,5)
     3           +bc(ep,30)*rxyz(ep,6)       
C----------J=4-- II=(J-1)*3 =9 JJ=(J-1)*2 =6 KK=10*(J-1)=30
         bcx=bcx
     1           +bc(ep,31)*vxyz(ep,10)+bc(ep,33)*vxyz(ep,11)
     2           +bc(ep,35)*vxyz(ep,12)+bc(ep,37)*rxyz(ep,7)
     3           +bc(ep,39)*rxyz(ep,8)       
         bcy=bcy
     1           +bc(ep,32)*vxyz(ep,10)+bc(ep,34)*vxyz(ep,11)
     2           +bc(ep,36)*vxyz(ep,12)+bc(ep,38)*rxyz(ep,7)
     3           +bc(ep,40)*rxyz(ep,8)
         vdef(ep,4)=tc(ep,1,1)*bcx+tc(ep,2,1)*bcy
         vdef(ep,5)=tc(ep,1,2)*bcx+tc(ep,2,2)*bcy
C      
 150   CONTINUE
       RETURN
       END
!||====================================================================
!||    cbadefsh   ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
        SUBROUTINE cbadefsh(JFT  ,JLT  ,X13  ,X24  ,Y13  ,Y24  ,
     1                                 BM, VDEF, VXYZ , NPLAT, IPLAT)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C        CALCUL 'membrane shear traitement' use only PART CONSTANTE
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   D U M M Y   A R G U M E N T S
C-----------------------------------------------
       INTEGER JFT ,JLT,NPLAT,IPLAT(*)
       my_real 
     .     X13(*)   ,X24(*)    ,Y13(*)     ,Y24(*)    ,
     .     BM(MVSIZ,36),VDEF(MVSIZ,8) ,VXYZ(MVSIZ,3,4)  
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
        INTEGER I, EP
C
       DO EP=jft,nplat 
        i=iplat(ep)
          vdef(i,3)=y24(i)*vxyz(i,2,1)-y13(i)*vxyz(i,2,2)
     .             -x24(i)*vxyz(i,1,1)+x13(i)*vxyz(i,1,2)
       ENDDO
       DO ep=nplat+1,jlt 
        i=iplat(ep)
          bm(i,3)=-x24(i)
          bm(i,6)=y24(i)
          bm(i,12)=x13(i)
          bm(i,15)=-y13(i)
          bm(i,21)=-bm(i,3)
          bm(i,24)=-bm(i,6)
          bm(i,30)=-bm(i,12)
          bm(i,33)=-bm(i,15)
          vdef(i,3)=y24(i)*(vxyz(i,2,1)-vxyz(i,2,3))+
     1              y13(i)*(-vxyz(i,2,2)+vxyz(i,2,4))
     2             -x24(i)*(vxyz(i,1,1)-vxyz(i,1,3))
     3             +x13(i)*(vxyz(i,1,2)-vxyz(i,1,4))
       ENDDO
       RETURN
       END
!||====================================================================
!||    cbaderirz    ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3     ../engine/source/elements/shell/coqueba/cbaforc3.F
!||    cbake3       ../engine/source/elements/shell/coqueba/cbake3.F
!||--- calls      -----------------------------------------------------
!||    cbaderirz0   ../engine/source/elements/shell/coqueba/cbadef.F
!||====================================================================
        SUBROUTINE cbaderirz(JFT  ,JLT  ,AREA ,X13 ,X24   ,
     2                       Y13  ,Y24  ,BM0RZ,BMKRZ,BMERZ,
     3                       VCORE,NPLAT,IPLAT,ISMSTR)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT,NPLAT,IPLAT(*),ISMSTR
      MY_REAL 
     .   VCORE(MVSIZ,3,4),X13(*),X24(*),Y13(*),Y24(*),
     .   BM0RZ(MVSIZ,4,4),BMKRZ(MVSIZ,4,4),BMERZ(MVSIZ,4,4),
     .   AREA(*)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP
      my_real 
     .   MX13(MVSIZ),MX23(MVSIZ),MX34(MVSIZ),MY13(MVSIZ),
     .   MY23(MVSIZ),MY34(MVSIZ),A1(MVSIZ),
     .   X13L(MVSIZ),X24L(MVSIZ),Y13L(MVSIZ),Y24L(MVSIZ)
C-----------------------------------------------
#include "vectorize.inc"
C-------VCORE(12):BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
      DO i=jft,nplat
        ep =iplat(i)
        mx23(i) = vcore(ep,1,3)
        my23(i) = vcore(ep,2,3)
        mx34(i) = vcore(ep,3,3)
        my34(i) = vcore(ep,1,4)
        mx13(i) = vcore(ep,2,4)
        my13(i) = vcore(ep,3,4)
        x13l(i) = x13(ep)
        x24l(i) = x24(ep)
        y13l(i) = y13(ep)
        y24l(i) = y24(ep)
        a1(i) = area(ep)
      END DO !I=JFT,NPLAT
C       
      DO i=nplat+1,jlt
       ep =iplat(i)
       mx13(i) = (vcore(ep,1,1)+vcore(ep,1,3))*half
       my13(i) = (vcore(ep,2,1)+vcore(ep,2,3))*half
       mx23(i) = (vcore(ep,1,2)+vcore(ep,1,3))*half
       my23(i) = (vcore(ep,2,2)+vcore(ep,2,3))*half
       mx34(i) = (vcore(ep,1,3)+vcore(ep,1,4))*half
       my34(i) = (vcore(ep,2,3)+vcore(ep,2,4))*half
       a1(i) = area(ep)
       x13l(i) = x13(ep)
       x24l(i) = x24(ep)
       y13l(i) = y13(ep)
       y24l(i) = y24(ep)
      END DO !I=NPLAT+1,JLT
C      
        CALL cbaderirz0(jft ,jlt  ,a1  ,x13l  ,x24l   ,
     2                 y13l  ,y24l  ,mx13 ,mx23 ,mx34  ,
     3                 my13 ,my23 ,my34 ,bm0rz,bmkrz ,
     4                 bmerz,ismstr)
C
      RETURN
      END
!||====================================================================
!||    cbadefrz   ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
        SUBROUTINE cbadefrz(JFT ,JLT  ,AREA  ,RLZ   ,VDEF ,
     1                     VXYZ ,BM0RZ,BMKRZ,BMERZ  ,VRLZ ,
     2                     BMRZ ,BRZ  ,BM   ,NPLAT  ,IPLAT,
     3                     NG   )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT,NPLAT,IPLAT(*),NG
      my_real 
     .   RLZ(MVSIZ,4),VXYZ(MVSIZ,3,4),
     .   BM0RZ(MVSIZ,4,4),BMKRZ(MVSIZ,4,4),BMERZ(MVSIZ,4,4),
     .   BMRZ(MVSIZ,3,4),BM(MVSIZ,36)
      my_real 
     .   AREA(*),VDEF(MVSIZ,8),VRLZ(*),BRZ(MVSIZ,4,4)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I,J,EP,NPG
      PARAMETER (NPG = 4)
      my_real 
     .   aa4,bxv2,byv1,erz,vpg(2,npg),pg,pg1
       parameter(pg=.577350269189626)
       parameter(pg1=-.577350269189626)
C--------------------------
C     INITIALISATION
C--------------------------
      DATA vpg/pg1,pg1,pg,pg1,pg,pg,pg1,pg/
C--------pay attention QEPH: Bi*A--- RLZ: VRZ/A-- QBAT:Bi--RLZ:VRZ
      DO i=jft,jlt
       DO j=1,4
        bmrz(i,1,j) = bm0rz(i,1,j)+bmkrz(i,1,j)*vpg(1,ng)+
     .                bmerz(i,1,j)*vpg(2,ng)
        bmrz(i,2,j) = bm0rz(i,2,j)+bmkrz(i,2,j)*vpg(1,ng)+
     .                bmerz(i,2,j)*vpg(2,ng)
        bmrz(i,3,j) = bm0rz(i,3,j)+bmkrz(i,3,j)*vpg(1,ng)+
     .                bmerz(i,3,j)*vpg(2,ng)
C
        brz(i,4,j) = bm0rz(i,4,j)+bmkrz(i,4,j)*vpg(1,ng)+
     .               bmerz(i,4,j)*vpg(2,ng)
       ENDDO
      ENDDO
C      
#include "vectorize.inc"
      DO i=jft,jlt
C---------------
C  MEMBRANE CONSTANT
C---------------
        ep=iplat(i)
       vdef(ep,1)=vdef(ep,1)+bmrz(i,1,1)*rlz(i,1)+bmrz(i,1,2)*rlz(i,2)
     1                      +bmrz(i,1,3)*rlz(i,3)+bmrz(i,1,4)*rlz(i,4)
       vdef(ep,2)=vdef(ep,2)+bmrz(i,2,1)*rlz(i,1)+bmrz(i,2,2)*rlz(i,2)
     1                      +bmrz(i,2,3)*rlz(i,3)+bmrz(i,2,4)*rlz(i,4)
       vdef(ep,3)=vdef(ep,3)+bmrz(i,3,1)*rlz(i,1)+bmrz(i,3,2)*rlz(i,2)
     1                      +bmrz(i,3,3)*rlz(i,3)+bmrz(i,3,4)*rlz(i,4)
      ENDDO
#include "vectorize.inc"
      DO i=jft,nplat
        ep=iplat(i)
        bxv2= bm(ep,1)*vxyz(ep,2,1)+bm(ep,2)*vxyz(ep,2,2)
     1       +bm(ep,3)*vxyz(ep,2,3)
        byv1=bm(ep,5)*vxyz(ep,1,1)+bm(ep,6)*vxyz(ep,1,2)
     1       +bm(ep,7)*vxyz(ep,1,3)
        vrlz(i)=(bxv2-byv1+
     1           brz(i,4,1)*rlz(i,1)+brz(i,4,2)*rlz(i,2)
     1          +brz(i,4,3)*rlz(i,3)+brz(i,4,4)*rlz(i,4))*half
      ENDDO
C      
#include "vectorize.inc"
      DO i=nplat+1,jlt 
       ep=iplat(i)
       erz=zero
       DO j=1,4
        erz=erz+brz(i,1,j)*vxyz(ep,1,j)+
     1          brz(i,2,j)*vxyz(ep,2,j)+brz(i,3,j)*vxyz(ep,3,j)+
     1          brz(i,4,j)*rlz(i,j)
       ENDDO
        vrlz(i)=erz*half
      ENDDO
C
      RETURN
      END
!||====================================================================
!||    cbaderirz0   ../engine/source/elements/shell/coqueba/cbadef.f
!||--- called by ------------------------------------------------------
!||    cbaderirz    ../engine/source/elements/shell/coqueba/cbadef.f
!||--- calls      -----------------------------------------------------
!||    cbpatch      ../engine/source/elements/shell/coqueba/cbadef.F
!||====================================================================
        SUBROUTINE cbaderirz0(JFT ,JLT  ,AREA ,X13 ,X24   ,
     2                       Y13  ,Y24  ,MX13 ,MX23,MX34  ,
     3                       MY13 ,MY23 ,MY34 ,BM0RZ,BMKRZ,
     4                       BMERZ,ISMSTR )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT,ISMSTR
      my_real
     .   x13(*),x24(*),y13(*),y24(*),
     .   mx13(*),mx23(*),mx34(*),my13(*),my23(*),my34(*)
      my_real
     .   area(*),
     .   bm0rz(mvsiz,4,4),bmkrz(mvsiz,4,4),bmerz(mvsiz,4,4)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I,J
      my_real
     .   aa4,bxv2,byv1,a05,nxy,nyx,a05k(4),a05e(4),
     .   xij(mvsiz,4),yij(mvsiz,4),xil(mvsiz,4),yil(mvsiz,4)
C--------(14,23,32,41)/4A----------X13:X13/A-----------------------------
      DO i=jft,jlt
       a05=half/area(i)
       xij(i,1) = a05*(mx13(i)-mx34(i))
       xij(i,2) = half*x13(i)+a05*(mx23(i)-mx13(i))
       xij(i,3) = -xij(i,2)
       xij(i,4) = -xij(i,1)
       yij(i,1) = a05*(my13(i)-my34(i))
       yij(i,2) = half*y13(i)+a05*(my23(i)-my13(i))
       yij(i,3) = -yij(i,2)
       yij(i,4) = -yij(i,1)
C--------(12,21,34,43)/4---------------------------------------
       xil(i,1) = a05*(mx13(i)-mx23(i))
       xil(i,2) = -xil(i,1)
       xil(i,3) = a05*(mx13(i)-mx34(i))-half*x13(i)
       xil(i,4) = -xil(i,3)
       yil(i,1) = a05*(my13(i)-my23(i))
       yil(i,2) = -yil(i,1)
       yil(i,3) = a05*(my13(i)-my34(i))-half*y13(i)
       yil(i,4) = -yil(i,3)
      ENDDO
C
      DO i=jft,jlt
C--------NxI,x ------
C-------- MY34YIJ KSI - MY23YIL ETA
       bm0rz(i,1,1) = (-my34(i)*yij(i,1)+my23(i)*yil(i,1))
       bm0rz(i,1,2) = ( my34(i)*yij(i,2)+my23(i)*yil(i,2))
       bm0rz(i,1,3) = ( my34(i)*yij(i,3)-my23(i)*yil(i,3))
       bm0rz(i,1,4) = (-my34(i)*yij(i,4)-my23(i)*yil(i,4))
C-------- MY13YIJ KSI - 2MY34YIL
       bmkrz(i,1,1) = (-my13(i)*yij(i,1)-two*my34(i)*yil(i,1))
       bmkrz(i,1,2) = ( my13(i)*yij(i,2)-two*my34(i)*yil(i,2))
       bmkrz(i,1,3) = ( my13(i)*yij(i,3)-two*my34(i)*yil(i,3))
       bmkrz(i,1,4) = (-my13(i)*yij(i,4)-two*my34(i)*yil(i,4))
C-------- -MY13YIL ETA + 2MY23YIJ
       bmerz(i,1,1) = ( my13(i)*yil(i,1)+two*my23(i)*yij(i,1))
       bmerz(i,1,2) = ( my13(i)*yil(i,2)+two*my23(i)*yij(i,2))
       bmerz(i,1,3) = (-my13(i)*yil(i,3)+two*my23(i)*yij(i,3))
       bmerz(i,1,4) = (-my13(i)*yil(i,4)+two*my23(i)*yij(i,4))
C -------NyI,y -------------
C-------- MX34XIJ KSI - MX23XIL ETA
       bm0rz(i,2,1) = (-mx34(i)*xij(i,1)+mx23(i)*xil(i,1))
       bm0rz(i,2,2) = ( mx34(i)*xij(i,2)+mx23(i)*xil(i,2))
       bm0rz(i,2,3) = ( mx34(i)*xij(i,3)-mx23(i)*xil(i,3))
       bm0rz(i,2,4) = (-mx34(i)*xij(i,4)-mx23(i)*xil(i,4))
C-------- MX13XIJ KSI - 2MX34XIL
       bmkrz(i,2,1) = (-mx13(i)*xij(i,1)-two*mx34(i)*xil(i,1))
       bmkrz(i,2,2) = ( mx13(i)*xij(i,2)-two*mx34(i)*xil(i,2))
       bmkrz(i,2,3) = ( mx13(i)*xij(i,3)-two*mx34(i)*xil(i,3))
       bmkrz(i,2,4) = (-mx13(i)*xij(i,4)-two*mx34(i)*xil(i,4))
C-------- -MX13XIL ETA + 2MX23XIJ
       bmerz(i,2,1) = ( mx13(i)*xil(i,1)+two*mx23(i)*xij(i,1))
       bmerz(i,2,2) = ( mx13(i)*xil(i,2)+two*mx23(i)*xij(i,2))
       bmerz(i,2,3) = (-mx13(i)*xil(i,3)+two*mx23(i)*xij(i,3))
       bmerz(i,2,4) = (-mx13(i)*xil(i,4)+two*mx23(i)*xij(i,4))
C--------NxI,y
C-------- -MX34YIJ KSI + MX23YIL ETA
       bm0rz(i,3,1) = ( mx34(i)*yij(i,1)-mx23(i)*yil(i,1))
       bm0rz(i,3,2) = (-mx34(i)*yij(i,2)-mx23(i)*yil(i,2))
       bm0rz(i,3,3) = (-mx34(i)*yij(i,3)+mx23(i)*yil(i,3))
       bm0rz(i,3,4) = ( mx34(i)*yij(i,4)+mx23(i)*yil(i,4))
C-------- -MX13YIJ KSI + 2MX34YIL
       bmkrz(i,3,1) = ( mx13(i)*yij(i,1)+two*mx34(i)*yil(i,1))
       bmkrz(i,3,2) = (-mx13(i)*yij(i,2)+two*mx34(i)*yil(i,2))
       bmkrz(i,3,3) = (-mx13(i)*yij(i,3)+two*mx34(i)*yil(i,3))
       bmkrz(i,3,4) = ( mx13(i)*yij(i,4)+two*mx34(i)*yil(i,4))
C-------- MX13YIL ETA - 2MX23YIJ
       bmerz(i,3,1) = (-mx13(i)*yil(i,1)-two*mx23(i)*yij(i,1))
       bmerz(i,3,2) = (-mx13(i)*yil(i,2)-two*mx23(i)*yij(i,2))
       bmerz(i,3,3) = ( mx13(i)*yil(i,3)-two*mx23(i)*yij(i,3))
       bmerz(i,3,4) = ( mx13(i)*yil(i,4)-two*mx23(i)*yij(i,4))
C--------NyI,x
C-------- -MY34XIJ KSI + MY23XIJ ETA
       bm0rz(i,4,1) = ( my34(i)*xij(i,1)-my23(i)*xil(i,1))
       bm0rz(i,4,2) = (-my34(i)*xij(i,2)-my23(i)*xil(i,2))
       bm0rz(i,4,3) = (-my34(i)*xij(i,3)+my23(i)*xil(i,3))
       bm0rz(i,4,4) = ( my34(i)*xij(i,4)+my23(i)*xil(i,4))
C-------- -MY13XIJ KSI + 2MY34XIL
       bmkrz(i,4,1) = ( my13(i)*xij(i,1)+two*my34(i)*xil(i,1))
       bmkrz(i,4,2) = (-my13(i)*xij(i,2)+two*my34(i)*xil(i,2))
       bmkrz(i,4,3) = (-my13(i)*xij(i,3)+two*my34(i)*xil(i,3))
       bmkrz(i,4,4) = ( my13(i)*xij(i,4)+two*my34(i)*xil(i,4))
C-------- MY13XIL ETA - 2MY23XIJ
       bmerz(i,4,1) = (-my13(i)*xil(i,1)-two*my23(i)*xij(i,1))
       bmerz(i,4,2) = (-my13(i)*xil(i,2)-two*my23(i)*xij(i,2))
       bmerz(i,4,3) = ( my13(i)*xil(i,3)-two*my23(i)*xij(i,3))
       bmerz(i,4,4) = ( my13(i)*xil(i,4)-two*my23(i)*xij(i,4))
      ENDDO
C--------correction to pass the patch test
      CALL cbpatch(jft ,jlt  ,bm0rz,bmkrz,bmerz  ,area ,
     +             mx13 ,mx23,mx34  ,my13 ,my23 ,my34)
C
        a05= half
       IF (ismstr==10) a05=zero
      DO i=jft,jlt
C--------BM0RZ(I,3,1)->(NxI,y+NyI,x) -;BM0RZ(I,4,1)->(-NxI,y+NyI,x -2NI)-----
       a05k(1)= -a05
       a05k(2)= a05
       a05k(3)= a05
       a05k(4)= -a05
       a05e(1)= -a05
       a05e(2)= -a05
       a05e(3)= a05
       a05e(4)= a05
      DO j=1,4
       nxy=bm0rz(i,3,j)
       nyx=bm0rz(i,4,j)
       bm0rz(i,3,j)=nxy+nyx
       bm0rz(i,4,j)=-nxy+nyx-a05
       nxy=bmkrz(i,3,j)
       nyx=bmkrz(i,4,j)
       bmkrz(i,3,j)=nxy+nyx
       bmkrz(i,4,j)=-nxy+nyx-a05k(j)
       nxy=bmerz(i,3,j)
       nyx=bmerz(i,4,j)
       bmerz(i,3,j)=nxy+nyx
       bmerz(i,4,j)=-nxy+nyx-a05e(j)
      ENDDO
      ENDDO
C
      RETURN
      END
!||====================================================================
!||    cbpatch      ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaderirz0   ../engine/source/elements/shell/coqueba/cbadef.F
!||    czderirz     ../engine/source/elements/shell/coquez/czdef.F
!||    czderirzt    ../engine/source/elements/shell/coquez/czdef.F
!||====================================================================
        SUBROUTINE cbpatch(JFT  ,JLT  ,BM0RZ,BMKRZ,BMERZ  ,AREA ,
     3                     MX13 ,MX23,MX34  ,MY13 ,MY23 ,MY34)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT
      my_real
     .   BM0RZ(MVSIZ,4,4),BMKRZ(MVSIZ,4,4),BMERZ(MVSIZ,4,4),
     .   MX13(*),MX23(*),MX34(*),MY13(*),MY23(*),MY34(*),AREA(*)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, J, NPG, NG, IAC
      PARAMETER (NPG = 4)
      my_real
     .   bm0rzm(mvsiz,4,4),j1,j2,j0,jac(mvsiz,4),det,jk_a,je_a
      my_real
     .   vpg(2,npg),pg,pg1
       parameter(pg=.577350269189626)
       parameter(pg1=-.577350269189626)
C--------------------------
C     INITIALISATION
C--------------------------
      DATA vpg/pg1,pg1,pg,pg1,pg,pg,pg1,pg/
C-------------Correction of BM0RZ to pass Patch test- IAC: analytical compute
      iac = 1
            IF (iac>0) THEN
      DO i=jft,jlt
             jk_a=four_over_3*(mx23(i)*my13(i)-mx13(i)*my23(i))/area(i)
             je_a=four_over_3*(mx13(i)*my34(i)-mx34(i)*my13(i))/area(i)
C--------NxI,x ------
       bm0rz(i,1,1) = -bmkrz(i,1,1)*jk_a-bmerz(i,1,1)*je_a
       bm0rz(i,1,2) = -bmkrz(i,1,2)*jk_a-bmerz(i,1,2)*je_a
       bm0rz(i,1,3) = -bmkrz(i,1,3)*jk_a-bmerz(i,1,3)*je_a
       bm0rz(i,1,4) = -bmkrz(i,1,4)*jk_a-bmerz(i,1,4)*je_a
C -------NyI,y -------------
       bm0rz(i,2,1) = -bmkrz(i,2,1)*jk_a-bmerz(i,2,1)*je_a
       bm0rz(i,2,2) = -bmkrz(i,2,2)*jk_a-bmerz(i,2,2)*je_a
       bm0rz(i,2,3) = -bmkrz(i,2,3)*jk_a-bmerz(i,2,3)*je_a
       bm0rz(i,2,4) = -bmkrz(i,2,4)*jk_a-bmerz(i,2,4)*je_a
C--------NxI,y ------
       bm0rz(i,3,1) = -bmkrz(i,3,1)*jk_a-bmerz(i,3,1)*je_a
       bm0rz(i,3,2) = -bmkrz(i,3,2)*jk_a-bmerz(i,3,2)*je_a
       bm0rz(i,3,3) = -bmkrz(i,3,3)*jk_a-bmerz(i,3,3)*je_a
       bm0rz(i,3,4) = -bmkrz(i,3,4)*jk_a-bmerz(i,3,4)*je_a
C -------NyI,y -------------
       bm0rz(i,4,1) = -bmkrz(i,4,1)*jk_a-bmerz(i,4,1)*je_a
       bm0rz(i,4,2) = -bmkrz(i,4,2)*jk_a-bmerz(i,4,2)*je_a
       bm0rz(i,4,3) = -bmkrz(i,4,3)*jk_a-bmerz(i,4,3)*je_a
       bm0rz(i,4,4) = -bmkrz(i,4,4)*jk_a-bmerz(i,4,4)*je_a
      ENDDO
      RETURN
            END IF !(ICN>0) THEN
C    
      DO i=jft,jlt
        j1=(mx23(i)*my13(i)-mx13(i)*my23(i))*pg
        j2=-(mx13(i)*my34(i)-mx34(i)*my13(i))*pg
        j0=area(i)*fourth
        jac(i,1)=j0+j2-j1
        jac(i,2)=j0+j2+j1
        jac(i,3)=j0-j2+j1
        jac(i,4)=j0-j2-j1
      ENDDO
C
      DO i=jft,jlt
       DO j=1,4
        bm0rzm(i,1,j)=zero
        bm0rzm(i,2,j)=zero
        bm0rzm(i,3,j)=zero
        bm0rzm(i,4,j)=zero
       ENDDO
      ENDDO
C
      DO ng= 1,4
      DO i=jft,jlt
             det=jac(i,ng)/area(i)
       DO j=1,4
        bm0rzm(i,1,j) = bm0rzm(i,1,j)+ (bmkrz(i,1,j)*vpg(1,ng)+
     .                  bmerz(i,1,j)*vpg(2,ng))*det
        bm0rzm(i,2,j) = bm0rzm(i,2,j)+(bmkrz(i,2,j)*vpg(1,ng)+
     .                  bmerz(i,2,j)*vpg(2,ng))*det
        bm0rzm(i,3,j) = bm0rzm(i,3,j)+(bmkrz(i,3,j)*vpg(1,ng)+
     .                  bmerz(i,3,j)*vpg(2,ng))*det
        bm0rzm(i,4,j) = bm0rzm(i,4,j)+(bmkrz(i,4,j)*vpg(1,ng)+
     .                  bmerz(i,4,j)*vpg(2,ng))*det
C
       ENDDO
      ENDDO
      ENDDO
C    
      DO i=jft,jlt
       DO j=1,4
        bm0rz(i,1,j) = -bm0rzm(i,1,j)
        bm0rz(i,2,j) = -bm0rzm(i,2,j)
        bm0rz(i,3,j) = -bm0rzm(i,3,j)
        bm0rz(i,4,j) = -bm0rzm(i,4,j)
C
       ENDDO
      ENDDO
C
      RETURN
      END
!||====================================================================
!||    cbadef1    ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
      SUBROUTINE cbadef1(JFT,JLT,NG,VCORE,VXYZ,VDEF,
     1                  HX,HY,BM,NPLAT,IPLAT,ISROT)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C        CALCULATION OF GENERALIZED DEFORMATIONS AT GAUSS POINTS AND [B] (MEMBRANE ONLY)
C        ENTREES :  NEL,NG,VCORE,VXYZ
C        SORTIES : VDEF,BM (MEMBRANE)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
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   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER NG,JFT,JLT,NPLAT,IPLAT(*),ISROT
      MY_REAL 
     .   VCORE(MVSIZ,12),VXYZ(MVSIZ,12),BM(MVSIZ,36),HX(MVSIZ,4),HY(MVSIZ,4)
      MY_REAL 
     .   VDEF(MVSIZ,8)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include "vectorize.inc"
       DO I=jft,nplat 
        ep=iplat(i)
C-------VCORE(12):BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
C-------V13(I)->VXYZ(j,1,I),V24(I)->VXYZ(j,2,I),VHI(I)->VXYZ(j,3,I)------------
        bm(ep,1)=vcore(ep,1)+hx(ep,ng)*vcore(ep,5)
        bm(ep,2)=vcore(ep,2)+hx(ep,ng)*vcore(ep,6)
        bm(ep,3)=hx(ep,ng)*fourth
        bm(ep,4)=-bm(ep,3)
        bm(ep,5)=vcore(ep,3)+hy(ep,ng)*vcore(ep,5)
        bm(ep,6)=vcore(ep,4)+hy(ep,ng)*vcore(ep,6)
        bm(ep,7)=hy(ep,ng)*fourth
        bm(ep,8)=-bm(ep,7)
C
        vdef(ep,1)=bm(ep,1)*vxyz(ep,1)+bm(ep,2)*vxyz(ep,4)
     1          +bm(ep,3)*vxyz(ep,7)
        vdef(ep,2)=bm(ep,5)*vxyz(ep,2)+bm(ep,6)*vxyz(ep,5)
     1          +bm(ep,7)*vxyz(ep,8)
        vdef(ep,4)=zero
        vdef(ep,5)=zero
        vdef(ep,6)=zero
        vdef(ep,7)=zero
        vdef(ep,8)=zero
       ENDDO
       IF (isrot>0) THEN
#include "vectorize.inc"
        DO i=jft,nplat
         ep=iplat(i)
C----------no more constant shear-----------
         vdef(ep,3)=bm(ep,1)*vxyz(ep,2)+bm(ep,2)*vxyz(ep,5)
     1             +bm(ep,3)*vxyz(ep,8)
     1             +bm(ep,5)*vxyz(ep,1)+bm(ep,6)*vxyz(ep,4)
     1             +bm(ep,7)*vxyz(ep,7)
        END DO
       END IF !(ISROT>0) THEN
C       
       RETURN
       END
!||====================================================================
!||    cbadeft    ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
      SUBROUTINE cbadeft(JFT,JLT,VXYZ,RLZ,VDEF,
     2                   BM,NPLAT,IPLAT,ISROT,BMRZ )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
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   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT,NPLAT,IPLAT(*),ISROT
      my_real 
     .   rlz(mvsiz,4),vxyz(mvsiz,3,4),
     .   bm(mvsiz,9,4),vdef(mvsiz,8)
C     .   BM(9,4,*),BMF(9,4,*),BF(6,4,*),VDEF(8,*)
      my_real 
     .   bmrz(mvsiz,4,4)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C---VDEF: 1:8 exx eyy,exy,eyx,byx,bxx,byy,bxy
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP
C
      MY_REAL 
     .   a_1,c11,c12,c21,c22,cdemi
C--------------------------
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include "vectorize.inc"
       DO i=jft,nplat 
        ep=iplat(i)
C-------BM(8):BXI(2),HX(1),not used,BYI(2),HY(1),not used
C-------V13(I)->VXYZ(j,1,I),V24(I)->VXYZ(j,2,I),VHI(I)->VXYZ(j,3,I)------------
C
        vdef(ep,1)=bm(ep,1,1)*vxyz(ep,1,1)+bm(ep,2,1)*vxyz(ep,1,2)
     1            +bm(ep,3,1)*vxyz(ep,1,3)
        vdef(ep,2)=bm(ep,5,1)*vxyz(ep,2,1)+bm(ep,6,1)*vxyz(ep,2,2)
     1            +bm(ep,7,1)*vxyz(ep,2,3)
C--------FLEXION-----
C        VDEF(EP,6)=BM(EP,1,1)*RXYZ(EP,2,1)+BM(EP,2,1)*RXYZ(EP,2,2)
C     1            +BM(EP,3,1)*RXYZ(EP,2,3)
C        VDEF(EP,7)=-(BM(EP,5,1)*RXYZ(EP,1,1)+BM(EP,6,1)*RXYZ(EP,1,2)
C     1            +BM(EP,7,1)*RXYZ(EP,1,3))
C        VDEF(EP,5)=-(BM(EP,1,1)*RXYZ(EP,1,1)+BM(EP,2,1)*RXYZ(EP,1,2)
C     1            +BM(EP,3,1)*RXYZ(EP,1,3))
C        VDEF(EP,8)=BM(EP,5,1)*RXYZ(EP,2,1)+BM(EP,6,1)*RXYZ(EP,2,2)
C     3            +BM(EP,7,1)*RXYZ(EP,2,3)
       ENDDO
c       constant shear is removed 
#include "vectorize.inc"
        DO i=jft,nplat
         ep=iplat(i)
C---------eyx-----------
         vdef(ep,4)=bm(ep,1,1)*vxyz(ep,2,1)+bm(ep,2,1)*vxyz(ep,2,2)
     1             +bm(ep,3,1)*vxyz(ep,2,3)
C-------- exy---     
         vdef(ep,3)=bm(ep,5,1)*vxyz(ep,1,1)+bm(ep,6,1)*vxyz(ep,1,2)
     1             +bm(ep,7,1)*vxyz(ep,1,3)
        END DO
c       
C------------ELEMENT GAUCH-------------------
#include "vectorize.inc"
       DO ep=nplat+1,jlt 
        i=iplat(ep)
C          I = 1; JJ = 3*KK =0 terms VXYZ(3,j,I) shared
          cdemi = half*(bm(i,9,1)*vxyz(i,3,1)+bm(i,9,2)*vxyz(i,3,2)+
     +                   bm(i,9,3)*vxyz(i,3,3)+bm(i,9,4)*vxyz(i,3,4))
          vdef(i,4)=bm(i,6,1)*vxyz(i,2,1)+bm(i,6,2)*vxyz(i,2,2)+
     1              bm(i,6,3)*vxyz(i,2,3)+bm(i,6,4)*vxyz(i,2,4)+cdemi
          vdef(i,3)=bm(i,3,1)*vxyz(i,1,1)+bm(i,3,2)*vxyz(i,1,2)+
     1              bm(i,3,3)*vxyz(i,1,3)+bm(i,3,4)*vxyz(i,1,4)+cdemi
C---------DEF MEMBRANE------
          vdef(i,1)= 
     +bm(i,1,1)*vxyz(i,1,1)+bm(i,4,1)*vxyz(i,2,1)+bm(i,7,1)*vxyz(i,3,1)+
     +bm(i,1,2)*vxyz(i,1,2)+bm(i,4,2)*vxyz(i,2,2)+bm(i,7,2)*vxyz(i,3,2)+
     +bm(i,1,3)*vxyz(i,1,3)+bm(i,4,3)*vxyz(i,2,3)+bm(i,7,3)*vxyz(i,3,3)+
     +bm(i,1,4)*vxyz(i,1,4)+bm(i,4,4)*vxyz(i,2,4)+bm(i,7,4)*vxyz(i,3,4)
          vdef(i,2)= 
     +bm(i,2,1)*vxyz(i,1,1)+bm(i,5,1)*vxyz(i,2,1)+bm(i,8,1)*vxyz(i,3,1)+
     +bm(i,2,2)*vxyz(i,1,2)+bm(i,5,2)*vxyz(i,2,2)+bm(i,8,2)*vxyz(i,3,2)+
     +bm(i,2,3)*vxyz(i,1,3)+bm(i,5,3)*vxyz(i,2,3)+bm(i,8,3)*vxyz(i,3,3)+
     +bm(i,2,4)*vxyz(i,1,4)+bm(i,5,4)*vxyz(i,2,4)+bm(i,8,4)*vxyz(i,3,4)
C
C---------FLEXION ASSOCIEE AUX U,V,W ([BMF]=2H[B0]+[B01])------
C----------[B01]---------
C          VDEF(I,6)= BMF(1,1,I)*VXYZ(1,1,I)+
C     +                BMF(4,1,I)*VXYZ(2,1,I)+BMF(7,1,I)*VXYZ(3,1,I)+
C     +                BMF(1,2,I)*VXYZ(1,2,I)+
C     +                BMF(4,2,I)*VXYZ(2,2,I)+BMF(7,2,I)*VXYZ(3,2,I)+
C     +                BMF(1,3,I)*VXYZ(1,3,I)+
C     +                BMF(4,3,I)*VXYZ(2,3,I)+BMF(7,3,I)*VXYZ(3,3,I)+
C     +                BMF(1,4,I)*VXYZ(1,4,I)+
C     +                BMF(4,4,I)*VXYZ(2,4,I)+BMF(7,4,I)*VXYZ(3,4,I)
C          VDEF(I,7)= BMF(2,1,I)*VXYZ(1,1,I)+
C     +                BMF(5,1,I)*VXYZ(2,1,I)+BMF(8,1,I)*VXYZ(3,1,I)+
C     +                BMF(2,2,I)*VXYZ(1,2,I)+
C     +                BMF(5,2,I)*VXYZ(2,2,I)+BMF(8,2,I)*VXYZ(3,2,I)+
C     +                BMF(2,3,I)*VXYZ(1,3,I)+
C     +                BMF(5,3,I)*VXYZ(2,3,I)+BMF(8,3,I)*VXYZ(3,3,I)+
C     +                BMF(2,4,I)*VXYZ(1,4,I)+
C     +                BMF(5,4,I)*VXYZ(2,4,I)+BMF(8,4,I)*VXYZ(3,4,I)
C          VXYZ(3,j,I) shared
C          CDEMI = HALF*(BMF(9,1,I)*VXYZ(3,1,I)+BMF(9,2,I)*VXYZ(3,2,I)+
C     +                   BMF(9,3,I)*VXYZ(3,3,I)+BMF(9,4,I)*VXYZ(3,4,I))
C          VDEF(I,8)= BMF(3,1,I)*VXYZ(1,1,I)+BMF(3,2,I)*VXYZ(1,2,I)+
C     +               BMF(3,3,I)*VXYZ(1,3,I)+BMF(3,4,I)*VXYZ(1,4,I)+CDEMI
C          VDEF(I,5)= BMF(6,1,I)*VXYZ(2,1,I)+BMF(6,2,I)*VXYZ(2,2,I)+
C     +               BMF(6,3,I)*VXYZ(2,3,I)+BMF(6,4,I)*VXYZ(2,4,I)+CDEMI
C---------- TERMES ASSOCIEES AU BETA ----------
C----------TERMES ASSOCIEES AUX BETA DE [B01]--
C           JJ = 3*II =0
C---------DEF DE FLEXION ASSOCIEE AUX THETA1,THETA2------
C          VDEF(I,6)= VDEF(I,6)+
C     +           BF(1,1,I)*RXYZ(1,1,I)+BF(4,1,I)*RXYZ(2,1,I)+
C     +           BF(1,2,I)*RXYZ(1,2,I)+BF(4,2,I)*RXYZ(2,2,I)+
C     +           BF(1,3,I)*RXYZ(1,3,I)+BF(4,3,I)*RXYZ(2,3,I)+
C     +           BF(1,4,I)*RXYZ(1,4,I)+BF(4,4,I)*RXYZ(2,4,I)
C          VDEF(I,7)= VDEF(I,7)+
C     +           BF(2,1,I)*RXYZ(1,1,I)+BF(5,1,I)*RXYZ(2,1,I)+
C     +           BF(2,2,I)*RXYZ(1,2,I)+BF(5,2,I)*RXYZ(2,2,I)+
C     +           BF(2,3,I)*RXYZ(1,3,I)+BF(5,3,I)*RXYZ(2,3,I)+
C     +           BF(2,4,I)*RXYZ(1,4,I)+BF(5,4,I)*RXYZ(2,4,I)
C          VDEF(I,8)= VDEF(I,5+
C     +           BF(6,1,I)*RXYZ(2,1,I)+BF(6,2,I)*RXYZ(2,2,I)+
C     +           BF(6,3,I)*RXYZ(2,3,I)+BF(6,4,I)*RXYZ(2,4,I)
C          VDEF(I,5)= VDEF(I,8+
C     +           BF(3,1,I)*RXYZ(1,1,I)+BF(3,2,I)*RXYZ(1,2,I)+
C     +           BF(3,3,I)*RXYZ(1,3,I)+BF(3,4,I)*RXYZ(1,4,I)
C      
      END DO
      IF (isrot>0) THEN
#include "vectorize.inc"
       DO i=jft,jlt
C---------------
C  MEMBRANE CONSTANT
C---------------
        ep=iplat(i)
       vdef(ep,1)=vdef(ep,1)+bmrz(i,1,1)*rlz(ep,1)+bmrz(i,1,2)*
     1      rlz(ep,2)+bmrz(i,1,3)*rlz(ep,3)+bmrz(i,1,4)*rlz(ep,4)
       vdef(ep,2)=vdef(ep,2)+bmrz(i,2,1)*rlz(ep,1)+bmrz(i,2,2)*
     1      rlz(ep,2)+bmrz(i,2,3)*rlz(ep,3)+bmrz(i,2,4)*rlz(ep,4)
       vdef(ep,3)=vdef(ep,3)+ 
     1                bmrz(i,3,1)*rlz(ep,1)+bmrz(i,3,2)*rlz(ep,2)
     2               +bmrz(i,3,3)*rlz(ep,3)+bmrz(i,3,4)*rlz(ep,4)
       vdef(ep,4)=vdef(ep,4)+ 
     1                bmrz(i,4,1)*rlz(ep,1)+bmrz(i,4,2)*rlz(ep,2)
     2               +bmrz(i,4,3)*rlz(ep,3)+bmrz(i,4,4)*rlz(ep,4)
       ENDDO
      END IF !(ISROT>0) THEN
C      
       RETURN
       END
!||====================================================================
!||    cbadeft1   ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
      SUBROUTINE cbadeft1(JFT,JLT,NG,VCORE,VXYZ,VDEF,
     1                    HX,HY,BM,NPLAT,IPLAT,ISROT,
     2                    BMRZ,RXYZ,WXY )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
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   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER NG,JFT,JLT,NPLAT,IPLAT(*),ISROT
      MY_REAL 
     .   vcore(mvsiz,12),vxyz(mvsiz,3,4),bm(mvsiz,9,4),hx(mvsiz,4),hy(mvsiz,4)
      my_real 
     .   vdef(mvsiz,8),rxyz(mvsiz,4),bmrz(mvsiz,4,4),wxy(*)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C---VDEF: 1:8 exx eyy,exy,eyx,byx,bxx,byy,bxy
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP
C
      MY_REAL 
     .   A_1,C11,C12,C21,C22,CDEMI
C--------------------------
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
       DO EP=jft,jlt 
C--------FLEXION-----
        vdef(ep,6)=zero
        vdef(ep,7)=zero
        vdef(ep,8)=zero
        vdef(ep,5)=zero
        wxy(ep)=zero
       ENDDO
#include "vectorize.inc"
       DO i=jft,nplat 
        ep=iplat(i)
C-------BM(8):BXI(2),HX(1),not used,BYI(2),HY(1),not used
        bm(ep,1,1)=vcore(ep,1)+hx(ep,ng)*vcore(ep,5)
        bm(ep,2,1)=vcore(ep,2)+hx(ep,ng)*vcore(ep,6)
        bm(ep,3,1)=hx(ep,ng)*fourth
        bm(ep,4,1)=-bm(ep,3,1)
        bm(ep,5,1)=vcore(ep,3)+hy(ep,ng)*vcore(ep,5)
        bm(ep,6,1)=vcore(ep,4)+hy(ep,ng)*vcore(ep,6)
        bm(ep,7,1)=hy(ep,ng)*fourth
        bm(ep,8,1)=-bm(ep,7,1)
C-------V13(I)->VXYZ(j,1,I),V24(I)->VXYZ(j,2,I),VHI(I)->VXYZ(j,3,I)------------
C
        vdef(ep,1)=bm(ep,1,1)*vxyz(ep,1,1)+bm(ep,2,1)*vxyz(ep,1,2)
     1            +bm(ep,3,1)*vxyz(ep,1,3)
        vdef(ep,2)=bm(ep,5,1)*vxyz(ep,2,1)+bm(ep,6,1)*vxyz(ep,2,2)
     1            +bm(ep,7,1)*vxyz(ep,2,3)
C---------eyx-----------
         vdef(ep,4)=bm(ep,1,1)*vxyz(ep,2,1)+bm(ep,2,1)*vxyz(ep,2,2)
     1             +bm(ep,3,1)*vxyz(ep,2,3)
C-------- exy---     
         vdef(ep,3)=bm(ep,5,1)*vxyz(ep,1,1)+bm(ep,6,1)*vxyz(ep,1,2)
     1             +bm(ep,7,1)*vxyz(ep,1,3)
       ENDDO
      IF (isrot>0) THEN
#include "vectorize.inc"
       DO i=jft,nplat
C---------------
C  MEMBRANE CONSTANT
C---------------
       ep=iplat(i)
       vdef(ep,1)=vdef(ep,1)+bmrz(i,1,1)*rxyz(ep,1)+bmrz(i,1,2)*
     1         rxyz(ep,2)+bmrz(i,1,3)*rxyz(ep,3)+bmrz(i,1,4)*rxyz(ep,4)
       vdef(ep,2)=vdef(ep,2)+bmrz(i,2,1)*rxyz(ep,1)+bmrz(i,2,2)*
     1         rxyz(ep,2)+bmrz(i,2,3)*rxyz(ep,3)+bmrz(i,2,4)*rxyz(ep,4)
       vdef(ep,3)=vdef(ep,3)+ 
     1                bmrz(i,3,1)*rxyz(ep,1)+bmrz(i,3,2)*rxyz(ep,2)
     2               +bmrz(i,3,3)*rxyz(ep,3)+bmrz(i,3,4)*rxyz(ep,4)
       vdef(ep,4)=vdef(ep,4)+ 
     1                bmrz(i,4,1)*rxyz(ep,1)+bmrz(i,4,2)*rxyz(ep,2)
     2               +bmrz(i,4,3)*rxyz(ep,3)+bmrz(i,4,4)*rxyz(ep,4)
       ENDDO
      END IF !(isrot>0) THEN
C      
       RETURN
       END
!||====================================================================
!||    cbaderit1   ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3    ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
      SUBROUTINE cbaderit1(JFT,JLT,NG,VCORE,VQ,VJFI,
     2                  HX,HY,VETA,VKSI,BM,NPLAT,IPLAT,
     3                  ISROT)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C        CALCULATION OF GENERALIZED DEFORMATIONS AT GAUSS POINTS AND [B]
C        ENTREES :  NEL,NG,VCORE,VQN,VXYZ,RXYZ,VNRM,VASM,VASTN
C        SORTIES : VDEF,BM (MEMBRANE),BMF(COUPLAGE MEM-FLEXION),BF (FLEXION),BC(C.T)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C-----------------------------------------------
C   I M P L I C I T   T Y P E S
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER NPG,NG,JFT,JLT,NPLAT,IPLAT(*),ISROT
      PARAMETER (NPG = 4)
      my_real 
     .   vcore(mvsiz,12),vksi(4,4),veta(4,4),
     .   bm(mvsiz,36),hx(mvsiz,4),hy(mvsiz,4)
      my_real 
     .   vjfi(mvsiz,3,2,4),vq(mvsiz,3,3,4)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP, IUN
      MY_REAL 
     .   VPG(2,NPG),PG1,PG,THK,DETJ,DET,
     .   TFN(3,2),BCX,BCY,BXY(3),BYX(3),
     .   V1(2),V2(2),RV1,RV2,
     .   C1,C2,VT1,VT2,BC1,BC2,VB1,TC(MVSIZ,2,2),
     .   VJF1(2,3),VJF(3,3),TBI(2,2),TBC(2,2),V11(4)
      MY_REAL 
     .   A_1,C11,C12,C21,C22,CC,BETA1,KSI1,KSIY1,BETA2,KSI2,KSIY2
       parameter(pg=.577350269189626)
       parameter(pg1=-.577350269189626)
C--------------------------
C     INITIALISATION
C--------------------------
      DATA iun/1/
      DATA vpg/pg1,pg1,pg,pg1,pg,pg,pg1,pg/
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include "vectorize.inc"
       DO i=jft,nplat 
        ep=iplat(i)
C-------VCORE(12):BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
        bm(ep,1)=vcore(ep,1)+hx(ep,ng)*vcore(ep,5)
        bm(ep,2)=vcore(ep,2)+hx(ep,ng)*vcore(ep,6)
        bm(ep,3)=hx(ep,ng)*fourth
        bm(ep,4)=-bm(ep,3)
        bm(ep,5)=vcore(ep,3)+hy(ep,ng)*vcore(ep,5)
        bm(ep,6)=vcore(ep,4)+hy(ep,ng)*vcore(ep,6)
        bm(ep,7)=hy(ep,ng)*fourth
        bm(ep,8)=-bm(ep,7)
       ENDDO
C------------ELEMENT GAUCH-------------------
#include "vectorize.inc"
       DO 150 i=nplat+1,jlt 
        ep=iplat(i)
C---------------------------
C  CALCUL DE TC=VJFI*VQ
C---------------------------
          tc(ep,1,1)=vjfi(ep,1,1,ng)*vq(ep,1,1,ng)+vjfi(ep,2,1,ng)
     +             *vq(ep,2,1,ng)+ vjfi(ep,3,1,ng)*vq(ep,3,1,ng)
          tc(ep,2,1)=vjfi(ep,1,2,ng)*vq(ep,1,1,ng)+vjfi(ep,2,2,ng)
     +             *vq(ep,2,1,ng)+ vjfi(ep,3,2,ng)*vq(ep,3,1,ng)
          tc(ep,1,2)=vjfi(ep,1,1,ng)*vq(ep,1,2,ng)+vjfi(ep,2,1,ng)
     +             *vq(ep,2,2,ng)+ vjfi(ep,3,1,ng)*vq(ep,3,2,ng)
          tc(ep,2,2)=vjfi(ep,1,2,ng)*vq(ep,1,2,ng)+vjfi(ep,2,2,ng)
     +             *vq(ep,2,2,ng)+ vjfi(ep,3,2,ng)*vq(ep,3,2,ng)
C--------------------------
C     BUILD [B] MATRIX AND DEFORMATION
C--------------------------
C--------------------------
C     DEFORMATION IN-PLANE
C--------------------------
C      [BM](3,3*NPG),[BMF](3,3*NPG),
C--------------------------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(1,ng)*tc(ep,1,1)+veta(1,ng)*tc(ep,2,1)
         c2=vksi(1,ng)*tc(ep,1,2)+veta(1,ng)*tc(ep,2,2)
C          I = 1; JJ = 3*KK =0
c         IF (ISROT>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          bm(ep,3)=bxy(1)+byx(1)
          bm(ep,6)=bxy(2)+byx(2)
          bm(ep,9)=bxy(3)+byx(3)
C---------CONTRIBUTION DU <T1> :EPS-X
          bm(ep,1)=vq(ep,1,1,ng)*c1
C---------CONTRIBUTION DU <T2> :EPS-Y
          bm(ep,2)=vq(ep,1,2,ng)*c2
C---------:EPS-X :EPS-Y :EPS-XY
          bm(ep,4)=vq(ep,2,1,ng)*c1
          bm(ep,5)=vq(ep,2,2,ng)*c2
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,7)=vq(ep,3,1,ng)*c1
          bm(ep,8)=vq(ep,3,2,ng)*c2
C
C--------J=2---II=(J-1)*2 =2  KK = 3*(J-1) =3
C--------- [B0], [B01]---------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(2,ng)*tc(ep,1,1)+veta(2,ng)*tc(ep,2,1)
         c2=vksi(2,ng)*tc(ep,1,2)+veta(2,ng)*tc(ep,2,2)
C         I = 1; JJ = 3*KK =9
c         IF (ISROT>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          bm(ep,12)=bxy(1)+byx(1)
          bm(ep,15)=bxy(2)+byx(2)
          bm(ep,18)=bxy(3)+byx(3)
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,10)=vq(ep,1,1,ng)*c1
          bm(ep,11)=vq(ep,1,2,ng)*c2
C         I = 2; JJ = 3*(KK+1)=12
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,13)=vq(ep,2,1,ng)*c1
          bm(ep,14)=vq(ep,2,2,ng)*c2
C         I = 3; JJ = 3*(KK+2)=15
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,16)=vq(ep,3,1,ng)*c1
          bm(ep,17)=vq(ep,3,2,ng)*c2
C---------DEF DE FLEXION ASSOCIEE AUX THETA1,THETA2------
C--------J=3---II=(J-1)*2 =4 KK = 3*(J-1)=6
C--------- [B0], [B01]---------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(3,ng)*tc(ep,1,1)+veta(3,ng)*tc(ep,2,1)
         c2=vksi(3,ng)*tc(ep,1,2)+veta(3,ng)*tc(ep,2,2)
c         IF (ISROT>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          bm(ep,21)=bxy(1)+byx(1)
          bm(ep,24)=bxy(2)+byx(2)
          bm(ep,27)=bxy(3)+byx(3)
C---------:EPS-X:EPS-Y:EPS-XY
          bm(ep,19)=vq(ep,1,1,ng)*c1
          bm(ep,20)=vq(ep,1,2,ng)*c2
C         I = 2; JJ = 3*(KK+1)=21
          bm(ep,22)=vq(ep,2,1,ng)*c1
          bm(ep,23)=vq(ep,2,2,ng)*c2
C         I = 3; JJ = 3*(KK+2)=24
          bm(ep,25)=vq(ep,3,1,ng)*c1
          bm(ep,26)=vq(ep,3,2,ng)*c2
C---------- TERMES ASSOCIEES AU BETA ----------
C--------TERMES DE [B0]=<T>CI-------
         c1=vksi(4,ng)*tc(ep,1,1)+veta(4,ng)*tc(ep,2,1)
         c2=vksi(4,ng)*tc(ep,1,2)+veta(4,ng)*tc(ep,2,2)
c         IF (ISROT>0) THEN
          bxy(1)=vq(ep,1,1,ng)*c2
          bxy(2)=vq(ep,2,1,ng)*c2
          bxy(3)=vq(ep,3,1,ng)*c2
          byx(1)=vq(ep,1,2,ng)*c1
          byx(2)=vq(ep,2,2,ng)*c1
          byx(3)=vq(ep,3,2,ng)*c1
          bm(ep,30)=bxy(1)+byx(1)
          bm(ep,33)=bxy(2)+byx(2)
          bm(ep,36)=bxy(3)+byx(3)
          bm(ep,28)=vq(ep,1,1,ng)*c1
          bm(ep,29)=vq(ep,1,2,ng)*c2
C         I = 2; JJ = 3*(KK+1)=30
          bm(ep,31)=vq(ep,2,1,ng)*c1
          bm(ep,32)=vq(ep,2,2,ng)*c2
C         I = 3; JJ = 3*(KK+2)=33
          bm(ep,34)=vq(ep,3,1,ng)*c1
          bm(ep,35)=vq(ep,3,2,ng)*c2
C
 150   CONTINUE
       RETURN
       END
!||====================================================================
!||    cbaderirzt   ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3     ../engine/source/elements/shell/coqueba/cbaforc3.F
!||====================================================================
      SUBROUTINE cbaderirzt(JFT,JLT,NG,BM0RZ,BMKRZ,BMERZ,BMRZ)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C        CALCULATION OF GENERALIZED DEFORMATIONS AT GAUSS POINTS AND [B]
C        ENTREES :  JFT,JLT,NG,BM0RZ,BMKRZ,BMERZ
C        SORTIES : BMRZ
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C-----------------------------------------------
C   I M P L I C I T   T Y P E S
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER NG,JFT,JLT
      MY_REAL 
     .   BM0RZ(MVSIZ,4,4),BMKRZ(MVSIZ,4,4),BMERZ(MVSIZ,4,4),
     .   BMRZ(MVSIZ,4,4)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER NPG,I,J
      PARAMETER (NPG = 4)
      my_real 
     .   vpg(2,npg),pg1,pg,nxyp,nyxn
       parameter(pg=.577350269189626)
       parameter(pg1=-.577350269189626)
C--------------------------
C     INITIALISATION
C--------------------------
      DATA vpg/pg1,pg1,pg,pg1,pg,pg,pg1,pg/
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C---+BMRZ(3,J,I): NXY ,BMRZ(4,J,I):NYX
      DO i=jft,jlt
       DO j=1,4
        bmrz(i,1,j) = bm0rz(i,1,j)+bmkrz(i,1,j)*vpg(1,ng)+
     .                bmerz(i,1,j)*vpg(2,ng)
        bmrz(i,2,j) = bm0rz(i,2,j)+bmkrz(i,2,j)*vpg(1,ng)+
     .                bmerz(i,2,j)*vpg(2,ng)
        nxyp        = bm0rz(i,3,j)+bmkrz(i,3,j)*vpg(1,ng)+
     .                bmerz(i,3,j)*vpg(2,ng)
        nyxn        = bm0rz(i,4,j)+bmkrz(i,4,j)*vpg(1,ng)+
     .                bmerz(i,4,j)*vpg(2,ng)
        bmrz(i,4,j) = half*(nxyp + nyxn)
        bmrz(i,3,j) = half*(nxyp - nyxn)
       ENDDO
      ENDDO
       RETURN
       END
!||====================================================================
!||    cbadeftw   ../engine/source/elements/shell/coqueba/cbadef.F
!||--- called by ------------------------------------------------------
!||    cbaforc3   ../engine/source/elements/shell/coqueba/cbaforc3.f
!||====================================================================
      SUBROUTINE cbadeftw(JFT,JLT,VXYZ,RXYZ,
     2                    BM,BMF,BF,NPLAT,IPLAT,
     3                    WXY )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
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   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER JFT,JLT,NPLAT,IPLAT(*)
      MY_REAL 
     .   RXYZ(MVSIZ,2,4),VXYZ(MVSIZ,3,4),
     .   BM(MVSIZ,9,4),BMF(MVSIZ,9,4),BF(MVSIZ,6,4),WXY(*)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
C---VDEF: 1:8 exx eyy,exy,eyx,byx,bxx,byy,bxy
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I, EP
C
      MY_REAL 
     .   A_1,C11,C12,C21,C22,CDEMI,KXY,KYX
C--------------------------
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#include "vectorize.inc"
       DO I=jft,nplat 
        ep=iplat(i)
        kyx =-(bm(ep,1,1)*rxyz(ep,1,1)+bm(ep,2,1)*rxyz(ep,1,2)
     1            +bm(ep,3,1)*rxyz(ep,1,3))
        kxy =bm(ep,5,1)*rxyz(ep,2,1)+bm(ep,6,1)*rxyz(ep,2,2)
     3            +bm(ep,7,1)*rxyz(ep,2,3)
        wxy(ep) = kxy - kyx
       ENDDO
C------------ELEMENT GAUCH-------------------
#include "vectorize.inc"
      DO i=nplat+1,jlt 
        ep=iplat(i)
C          
        cdemi =half*(bmf(ep,9,1)*vxyz(ep,3,1)+bmf(i,9,2)*vxyz(ep,3,2)+
     +                 bmf(ep,9,3)*vxyz(ep,3,3)+bmf(i,9,4)*vxyz(ep,3,4))
          kxy = bmf(ep,3,1)*vxyz(ep,1,1)+bmf(ep,3,2)*vxyz(ep,1,2)+
     +          bmf(ep,3,3)*vxyz(ep,1,3)+bmf(ep,3,4)*vxyz(ep,1,4)+cdemi
          kyx = bmf(ep,6,1)*vxyz(ep,2,1)+bmf(ep,6,2)*vxyz(ep,2,2)+
     +          bmf(ep,6,3)*vxyz(ep,2,3)+bmf(ep,6,4)*vxyz(ep,2,4)+cdemi
          kxy = kxy +
     +          bf(ep,6,1)*rxyz(ep,2,1)+bf(ep,6,2)*rxyz(ep,2,2)+
     +          bf(ep,6,3)*rxyz(ep,2,3)+bf(ep,6,4)*rxyz(ep,2,4)
          kyx = kyx +
     +          bf(ep,3,1)*rxyz(ep,1,1)+bf(ep,3,2)*rxyz(ep,1,2)+
     +          bf(ep,3,3)*rxyz(ep,1,3)+bf(ep,3,4)*rxyz(ep,1,4)
C      
        wxy(ep) = kxy - kyx
      END DO
C      
       RETURN
       END