      SUBROUTINE L2DFCV (ELK,ELM,ELF,GEOM,NODEL,KSW,RO,CP,A,TELEM,DT)
      DOUBLE PRECISION ELK(*),ELF(*),ELM(*),GEOM(*)
      INTEGER NODEL, KSW
      DOUBLE PRECISION RO,CP,TELEM(*),DT
      DOUBLE PRECISION ZERO,PII2
      PARAMETER (ZERO = 0.D0, PII2 = 2.D0*3.141592653589793D0)
      INTEGER DIM, DOFNOD, MNODEL, MNQP, MDOFEL
      PARAMETER (DIM = 2, DOFNOD = 1)
      PARAMETER (MNODEL = 16, MNQP = 9, MDOFEL = MNODEL*DOFNOD)
      DOUBLE PRECISION BTAB(MDOFEL*MDOFEL), B(DIM*MDOFEL)
      DOUBLE PRECISION A(DIM*DIM), BT(MDOFEL*DIM), AB(DIM*MDOFEL) 
      DOUBLE PRECISION JAC(DIM,DIM), JACIN(DIM,DIM)
      DOUBLE PRECISION SFUN(MNODEL), LDER(DIM*MNODEL)
      DOUBLE PRECISION SFUN2(MNODEL*MNODEL),XK(DIM)
      DOUBLE PRECISION WG(MNQP), XI(MNQP), ETA(MNQP)
      DOUBLE PRECISION DETK, AMUL, FXK
      DOUBLE PRECISION VELO(DIM),SFUGR(MDOFEL*MDOFEL),
     +                 VGRAD(MDOFEL)
      INTEGER DOFEL, NQP, I, K

      DOFEL = DOFNOD*NODEL
      CALL GAUSQ4(WG,XI,ETA,NQP)
      CALL VFILL (DOFEL*DOFEL,ELK,ZERO)
      CALL VFILL (DOFEL*DOFEL,ELM,ZERO)
      CALL VFILL (DOFEL,ELF,ZERO)

      DO 100 K=1,NQP
         CALL SQUAD8(SFUN,LDER,XI(K),ETA(K))
         CALL MXM (SFUN,1,GEOM,NODEL,XK,DIM)
         CALL MXM (LDER,DIM,GEOM,NODEL,JAC,DIM)
         CALL INV2X2 (JAC,JACIN,DETK)
         CALL MXM (JACIN,DIM,LDER,DIM,B,NODEL)

         CALL MATRAN (B,BT,DIM,DOFEL)
         CALL MXM (A,DIM,B,DIM,AB,DOFEL)
         CALL MXM (BT,DOFEL,AB,DIM,BTAB,DOFEL)
         CALL MXM (SFUN,DOFEL,SFUN,1,SFUN2,DOFEL)
         IF (KSW.EQ.1) THEN
            AMUL = WG(K)*DETK
         ELSE
            AMUL = WG(K)*PII2*XK(1)*DETK
         ENDIF


         CALL SAXPY (DOFEL*DOFEL,AMUL*DT,BTAB,1,ELK,1)


         CALL SAXPY (DOFEL*DOFEL,AMUL*RO*CP,SFUN2,1,ELM,1)



         FXK=0.0D0
         DO 500 I=1,NODEL
            FXK=FXK+TELEM(I)*SFUN(I)
 500     CONTINUE

         IF (DABS(FXK).GT.ZERO)
     &      CALL SAXPY (DOFEL,RO*CP*AMUL*FXK,SFUN,1,ELF,1)


 100  CONTINUE

      RETURN
      END
