! Subrutina de integrare a unui sistem de ec. neliniare cu NEWOP
! Ordinul sistemului = ngl
! MU'' +G(U,U') +F(U) =P(t)
! Cu sau fara mediere pt. U3
! Subrutinele pentru metoda Newton:
!       FGP_J: f, g, p; Jacobienii: jf; jg, jg1;
!       SUB-AB: A(ng.ngl); b(ngl)
!		ELIM (...): eliminare Gauss
! 
! w=dU2

subroutine Newmark(t, PrintIt, kod_Elim)
use work

implicit real(8) (a-h,o-z)


real(8) X(ngl), X1(ngl), X2(ngl)
real(8) norma
real(8) w(ngl), delta_w(ngl)
real(8) a(ngl,ngl), b(ngl)

common /coef/	beta_h2, gama_h

logical PrintIt


! Calculeaza "U_bar":
uEnd =  u0 +u10*h +u20*h*h/2. ! u_bar
u1End= u10 +u20*h  ! u'_bar
u2End= u20
       
! Rezolvare cu met. NEWTON: w = dU3
! J(W)*delta_w=-F(W)
! w=w+delta_w

w =0
nit=0

90 nit=nit+1
! Argumentele pt. F, JF, JG: X=U~+W*..., X1=U'~+W*..., X2=U''~+W*...;
X  =  uEnd + beta_h2 *w
X1 = u1End + gama_h*w
X2 = u2End + w         

	call FGP_J(t, X, X1)
	call sub_ab(X, X1, X2, a, b)	! stabileste a =JF; b =Termenii liberi.

kod =0
	call ELIM (ngl, A, b, prag, kod_Elim, kodPrint)	!

if(kod_Elim <0) RETURN

	delta_w = b(:)
	w = w + delta_w

if (norma(delta_w) <=eps .or. nit >=lnit) then
goto 80
else
goto 90
endif

80 continue
	if(PrintIt) then
	write(9,'(1pE12.4, i3, 1pE13.4)') t, nit, norma(delta_w)
	endif

uEnd  =  uEnd  + beta_h2*w
u1End = u1End  + gama_h*w
u2End = u2End  + w


500 RETURN
END

