      SUBROUTINE L2NAVI (ELK,ELF,GEOM,VISCOS,DENSIT,DT,ITVELO,EPS,NODEL,
     + 		ITTEMP,GRAVX,GRAVY,HEXPCO,REFTEM)
      DOUBLE PRECISION ELK(*),ITTEMP(*),ITVELO(*),ELF(*),GEOM(*)
      DOUBLE PRECISION VISCOS,DENSIT,DT,EPS,GRAVX,GRAVY,HEXPCO,REFTEM
      INTEGER NODEL
      DOUBLE PRECISION ZERO,PII2
      INTEGER DOFNOD, MNODEL, DOFEL, MDOFEL, NS, DIM, MQP
      PARAMETER ( DOFNOD=2, MDOFEL=16, MNODEL=8, MQP=4 )
      PARAMETER (NS=4,DIM=2,ZERO=0.D0,PII2=2.D0*3.141592653589793D0)
      DOUBLE PRECISION BTB(MDOFEL*MDOFEL), B(NS*MDOFEL)
      DOUBLE PRECISION GDER(DIM,MNODEL), SFUN(MNODEL)
      DOUBLE PRECISION FUN(DIM,MDOFEL),BT(MDOFEL*NS)
      DOUBLE PRECISION BFOR(DIM), JAC(DIM,DIM), JACIN(DIM,DIM)
      DOUBLE PRECISION BFUN(MDOFEL)
      DOUBLE PRECISION LDER(DIM,MNODEL),XK(DIM),SFUN2(MDOFEL*MDOFEL),
     +     SFUMU(2*MDOFEL),SFUMUT(MDOFEL*2),UZZERO(2*MDOFEL),
     +     SFUNR(MNODEL)
      DOUBLE PRECISION WG(MQP), XI(MQP), ETA(MQP)
      DOUBLE PRECISION DIVDIV(MDOFEL*MDOFEL), DIV(MDOFEL)

      DOUBLE PRECISION VGRAD(MNODEL), MVGRAD(DIM,MDOFEL)
      DOUBLE PRECISION VELO(DIM), SFVGR(MDOFEL*MDOFEL)
      DOUBLE PRECISION DETK, AMUL,ELM(MDOFEL,MDOFEL)
      DOUBLE PRECISION PDER(MDOFEL),PDERIN(MDOFEL),PRES(MDOFEL*MDOFEL),
     +     PFUN(MDOFEL),PFUNI(MDOFEL),MEASEL,g(8)
      INTEGER I,J,K,IND,NQP
      DOFEL = NODEL*DOFNOD
      CALL VFILL (DOFEL*DOFEL,ELK,ZERO)
      CALL VFILL (DOFEL*DOFEL,ELM,ZERO)
      CALL VFILL (DOFEL,ELF,ZERO)

      CALL GAUSQ4 (WG,XI,ETA,NQP)

      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,GDER,NODEL)
         CALL B2STOK (B,GDER,SFUMU,SFUN,SFUN,NODEL)
         CALL MATRAN (B,BT,NS,DOFEL)
         CALL MATRAN (SFUMU,SFUMUT,DIM,DOFEL)
         CALL MXM (BT,DOFEL,B,NS,BTB,DOFEL)
         CALL MXM (SFUMUT,DOFEL,SFUMU,DIM,SFUN2,DOFEL)

         AMUL = WG(K)*DETK*XK(1)


         DO 10 I=1,DOFEL*DOFEL
            ELK(I) = ELK(I) + DT*VISCOS*AMUL*BTB(I)
 10      CONTINUE                 


         CALL SAXPY (DOFEL*DOFEL,AMUL*DENSIT,SFUN2,1,ELK,1)
         CALL SAXPY (DOFEL*DOFEL,AMUL*DENSIT,SFUN2,1,ELM,1)


                  
         BFOR(1) = ZERO
         BFOR(2) = ZERO-DENSIT*GRAVY

         DO 400 I=1,NODEL
            BFOR(1) = BFOR(1)+DENSIT*GRAVX*
     +           (HEXPCO*(ITTEMP(I)-REFTEM))*SFUN(I)
            BFOR(2) = BFOR(2)+DENSIT*GRAVY*
     +           (HEXPCO*(ITTEMP(I)-REFTEM))*SFUN(I)
 400     CONTINUE

         CALL VFILL (DOFEL*DIM,FUN,ZERO)
         DO 80 I=1,NODEL
            FUN(1,I) = SFUN(I)
            FUN(2,I+NODEL) = SFUN(I)
 80      CONTINUE     

         IF (BFOR(1).NE.0.D0.OR.BFOR(2).NE.0.D0) THEN
            CALL MXM (FUN,DOFEL,BFOR,DIM,BFUN,1)
            CALL SAXPY (DOFEL,DT*AMUL,BFUN,1,ELF,1)
         END IF


         AMUL=VISCOS*DT*WG(K)*DETK/XK(1)
         CALL VFILL(2*MDOFEL,UZZERO,ZERO)

         CALL B2STOK (B,GDER,SFUMU,SFUN,UZZERO,NODEL)
         CALL MATRAN (SFUMU,SFUMUT,DIM,DOFEL)
         CALL MXM (SFUMUT,DOFEL,SFUMU,DIM,SFUN2,DOFEL)
         CALL SAXPY (DOFEL*DOFEL,AMUL,SFUN2,1,ELK,1)

 100  CONTINUE


      DO 111 I=1,DOFEL
         DO 112 J=1,DOFEL
            ELF(I)=ELF(I)+ELM(I,J)*ITVELO(J) 
 112     CONTINUE
 111  CONTINUE                 


      
      CALL GAUSQ4 (WG,XI,ETA,NQP)

      DO 200 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,DIV,NODEL)

         CALL SAXPY (NODEL,1.D0/XK(1),SFUN,1,DIV(1),2)

         AMUL = DT*WG(1)*DETK*XK(1)/EPS

         CALL MXM (DIV,DOFEL,DIV,1,DIVDIV,DOFEL)
         CALL SAXPY (DOFEL*DOFEL,AMUL,DIVDIV,1,ELK,1)
 200  CONTINUE

      RETURN
      END
