! 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; Derivata dp
!       SUB-AB: A(ng.ngl); b(ngl)
!		ELIM (...): eliminare Gauss
! 
! w=dU3

subroutine NewOp(t, PrintIt, kod_Elim)
use work

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

real(8) u3eq(ngl) ! pt. mediere

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_h3, gama_h2, delta_h, &
				teta_beta_h3, teta1_gama_h2, teta2_delta_h
logical PrintIt


! Calculeaza "U~":
uEnd =  u0 +u10*h +u20*h*h/2. +u30*beta_h3 ! u~
u1End= u10 +u20*h +u30*gama_h2  ! u'~
u2End= u20 +u30*delta_h  ! u"~
u3End= u30
       
! 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 + teta_beta_h3 *w
X1 = u1End + teta1_gama_h2*w
X2 = u2End + teta2_delta_h*w         

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

kod =0; kodPrint =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  + teta_beta_h3  *w
u1End = u1End  + teta1_gama_h2 *w
u2End = u2End  + teta2_delta_h *w
u3End = u3End  + w


IF(kodm /=0) THEN
! Mediere:
! M*U3 = dP - jG*U1 - jG1*U2 - jF*U1
	call FGP_J(t, uEnd, u1End)

DO l =1,NGL
b1 =0d0
do j =1,ngl
b1 = b1 + jG(l,j)*U1End(j) + jG1(l,j)*U2End(j) + jF(l,j)*U1End(j)
enddo

b(l) =dP(l) -b1
ENDDO

kod =0; kodPrint =0;
call gen_M()
call ELIM (ngl, M, b, prag, kod, kodPrint)	!
u3eq =b

u3End =niu*u3End +(1-niu)*u3eq

ENDIF

500 RETURN
END

