lag_direct.F File Reference

#include "implicit_f.inc"
#include "lagmult.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	lag_direct (iadll, lll, jll, xll, ltsm, v, vr, a, ar, ms, in, nc_ini, ncl)

Function/Subroutine Documentation

lag_direct()

subroutine lag_direct	(	integer, dimension(*)	iadll,
		integer, dimension(*)	lll,
		integer, dimension(*)	jll,
			xll,
			ltsm,
			v,
			vr,
			a,
			ar,
			ms,
			in,
		integer	nc_ini,
		integer	ncl )

Definition at line 34 of file lag_direct.F.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "lagmult.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NC_INI, NCL, IADLL(*), LLL(*), JLL(*)
C     REAL
      my_real
     .   ltsm(6,*),xll(*),ms(*),in(*),v(3,*),vr(3,*),a(3,*),ar(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,K,L,IC,JC,IK
      my_real
     .        hloc(mxdlen,mxdlen),rloc(mxdlen),s,hij
C======================================================================|
      IF (ncl>mxdlen) THEN
        CALL ancmsg(msgid=111,anmode=aninfo,
     .            i1=ncl)
        CALL arret(2)
      ENDIF
C---  Local H matrix
      DO k=1,ncl
        rloc(k) = zero
        ic = nc_ini + k
        DO ik=iadll(ic),iadll(ic+1)-1
          i = lll(ik)
          j = jll(ik)
          IF (j>3) THEN
            ltsm(j,i) = xll(ik)/in(i)
          ELSE
            ltsm(j,i) = xll(ik)/ms(i)
          ENDIF
        ENDDO
        DO l = 1,k
          jc = nc_ini + l
          hij = zero
          DO ik=iadll(jc),iadll(jc+1)-1
            hij = hij + xll(ik)*ltsm(jll(ik),lll(ik))
          ENDDO
          hloc(l,k)  = hij
        ENDDO
        DO ik=iadll(ic),iadll(ic+1)-1
          ltsm(jll(ik),lll(ik)) = zero
        ENDDO
      ENDDO
      DO k = 2,ncl
        DO l = 1,k
          hloc(k,l) = hloc(l,k)
        ENDDO
      ENDDO
C---  second membre
      DO k = 1,ncl
        ic = nc_ini + k
        DO ik=iadll(ic),iadll(ic+1)-1
          i = lll(ik)
          j = jll(ik)
          IF (j>3) THEN
            j = j-3
            rloc(k) = rloc(k) + xll(ik)*(vr(j,i)/dt12+ar(j,i))
          ELSE
            rloc(k) = rloc(k) + xll(ik)*(v(j,i)/dt12+a(j,i))
          ENDIF
        ENDDO
      ENDDO
C---  Factorise H (Full Cholesky)
      DO j=1,ncl
              IF (hloc(j,j)<=zero) THEN
          CALL ancmsg(msgid=112,anmode=aninfo,
     .            i1=j)
          CALL arret(2)
              ENDIF
        hloc(j,j) = sqrt(hloc(j,j))
        DO k=1,j-1
          DO i=j+1,ncl
                  hloc(i,j) = hloc(i,j) - hloc(i,k)*hloc(j,k)
          ENDDO
        ENDDO
        DO i=j+1,ncl
          hloc(i,j) = hloc(i,j)/hloc(j,j)
          hloc(i,i) = hloc(i,i) - hloc(i,j)*hloc(i,j)
        ENDDO
      ENDDO
C---  back subst Ly = r,
      DO i=1,ncl
        s = rloc(i)
        DO j=1,i-1
          s = s - hloc(i,j)*rloc(j)
        ENDDO
              rloc(i) = s / hloc(i,i)
      ENDDO
C---  back subst Lz = y
      DO i=ncl,1,-1
        s = rloc(i)
        DO j=i+1,ncl
          s = s - hloc(j,i)*rloc(j)
        ENDDO
              rloc(i) = s / hloc(i,i)
      ENDDO
C---  update accelerations
      DO k=1,ncl
        ic = nc_ini + k
        DO ik=iadll(ic),iadll(ic+1)-1
          i = lll(ik)
          j = jll(ik)
          IF(j>3) THEN
            j = j-3
            ar(j,i) = ar(j,i) - xll(ik)*rloc(k)/in(i)
          ELSE
            a(j,i)  = a(j,i)  - xll(ik)*rloc(k)/ms(i)
          ENDIF
        ENDDO
      ENDDO
C---
      RETURN