

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!===(!!!) Cos(Pi*x), Sin(Pi*x)   '../f.txt'  .      !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	
	program START

	use AVDef    
	use DFLIB
	
	integer, parameter:: dil=150,dl=64,del=31 

      integer, parameter::n=(dil)*4,k=n/2,nl=n/2,np=n/2
	  	common Ph0,Ph1,Ph2,Ph3,Ph4,Pha0T,oleole
	Common PhN1,PhN2,PhN3,PhN1T
	real z(N*(N-K+1)+(N-K+1)*(N-K+2)/2+1)
     real r1(N*(N-K+1)),z1(N*(N-K+1)),y(n+1),a1(N*(N+1)),oleole(N*(N+1))
     real om((N-K+1)*(N-K+2)/2)
     real rew(dil+4,31),ww(dl,dil+4),ww1(dil+4,dl)
	 real Aleft(K*(n+1)),Aright(K*(N+1)),AMat1(n/2,n+1),Amat2(n/2,n+1)
	 dimension itp(n)
	integer num1
	real xx(-5:dil+4)
	real ess(dil),x,ds,esd,ess1(n/4),ess2(n/4)
	real Mm(del,dl)
      real sina(dl)
		real suk,vnu,tempv
		character key
	real Pha0T(n/4,n/4)
	real Ph0(n/4,n/4),Ph1(n/4,n/4)
	real Ph2(n/4,n/4),Ph3(n/4,n/4),Ph4(n/4,n/4)
	real PhN1(n/4,n/4),PhN1T(n/4,n/4)
	real PhN2(n/4,n/4),PhN3(n/4,n/4),PhNN(n/4,n/4)
	open (6,file='dior.txt',status='unknown')
	open (7,file='dior1.txt',status='unknown')
		open (8,file='dior2.txt',status='unknown')

!****************************** ITP *************************************	
	

 itp=(/(i,i=1,n)/)	 !   
!======        ========
 do i=n/3 + 1,n/2
	itp(i) = itp(i) + n/6;
 enddo

 do i=n/2 + 1,2*n/3
	itp(i) = itp(i) - n/6;
 enddo
!======================================================================

!======    =======================================
!do i=n/6 + 1,n/3
!  itp(i) = itp(i) + 4*n/6;
!enddo
 
!do i=5*n/6 + 1,n
!  itp(i) = itp(i) - 4*n/6;
!enddo
!======================================================================

!************************************************************************



!===================================================
AMat1=0.; AMat2=0.; oleole=0.;

Aleft=0.
Aright=0.
!--     ------------------------------
open(102,file='../BqLeft.txt',status='unknown')
      do i=1,n/2
	  do j=1,n 
		read(102,*)AMat1(i,j)
		
	  enddo
   enddo
close(102)

open(102,file='../Bq.txt',status='unknown')
      do i=1,n/2
	  do j=1,n 
		read(102,*)AMat2(i,j)
		
	  enddo
   enddo
close(102)



!--      . Left ~ z = 0,  Right ~ z = c.
open(103,file='../f0.txt',status='unknown')
	  do i=1,n/2
		read(103,*)AMat1(i,n + 1)
	  enddo
close(103)

open(104,file='../fc.txt',status='unknown')
	  do i=1,n/2
		read(104,*)AMat2(i,n + 1)
	  enddo
close(104)

do j=1,n+1 
do i=1,n/2
Aleft((j-1)*N/2+i)=AMat1(i,j)
Aright((j-1)*N/2+i)=AMat2(i,j)
enddo
enddo
!===================================================

!--   --------------------------------------
	open(105,file='../A.txt',status='unknown')
      do j=1,n 
	  do i=1,n
		read(105,*)oleole((j-1)*N+i)
	  enddo
	  enddo
	close(105)

!--     -------------------------
	open(106,file='../f.txt',status='unknown')
	  do i=1,n
		read(106,*)oleole(N*N+i)
	  enddo
	close(106)
!---------------------------------------------------------
!=========================================================
!---------------------------
open (666,file='Aleft.txt',status='unknown')
write (666,1001)Aleft
close(666)

open (667,file='Aright.txt',status='unknown')
write (667,1001)Aright
close(667)

1001  FORMAT(1(E12.5,2X))
!---------------------------

PhN1=0.;PhN2=0.;PhN3=0.



      N1=(N-K+1)*N
      M1=(N-K+1)*(N-K+2)/2
      NZ=N1+M1+1
      NY=N+1
      NN=N*NY
      CALL DO(N,K,X,H,N1,M1,NZ,NY,NN,NL,NP, &
     Z,R1,Z1,Y,A1,OM,Aleft,Aright,ITP)
  
  close(6)

  open(6,file='dior.txt',status='unknown')
	
    

print *, "Complete"
key = GETCHARQQ()
 
END


	  SUBROUTINE FORWAR(N,M,X,H,IP1,IP2,IP3,NA,A)
      real(8) A(NA),h,x
	  IP1=1; IP2=0;   IP3=1
      
	  H=0.01

	  IF(ABS(0.-X).LE.H/2.) IP2=1
	  IF(ABS(0.1-X).LE.H/2.) IP2=1
	  IF(ABS(0.2-X).LE.H/2.) IP2=1
	  IF(ABS(0.3-X).LE.H/2.) IP2=1
	  IF(ABS(0.4-X).LE.H/2.) IP2=1
	  IF(ABS(0.5-X).LE.H/2.) IP2=1
	  IF(ABS(0.6-X).LE.H/2.) IP2=1
	  IF(ABS(0.7-X).LE.H/2.) IP2=1
	  IF(ABS(0.8-X).LE.H/2.) IP2=1
	  IF(ABS(0.9-X).LE.H/2.) IP2=1
	  IF(ABS(1.-X).LE.H/2.) IP2=1

      IF(ABS(1.-X).LE.H/2.) IP3=0
      RETURN
      END

      SUBROUTINE RESULT(N,NI,A)
      real(8) A(N)
		!write (6,1)a(2:n)
		do i=2,n
		   write (6,1)a(i)
		1  FORMAT(600(E12.5,1X))
		enddo

      RETURN
      END
 

     
    SUBROUTINE MATRIX(NN,X,NA,A)
	
	integer, parameter:: dil=150,n=(dil)*4
	common Ph0,Ph1,Ph2,Ph3,Ph4,Pha0T,oleole
	Common PhN1,PhN2,PhN3,PhN1T
      real A(NA),oleole(N*(N+1))
	real h,x,vn_st,vn_ts 
	real tmp4(n/4,n/4)
	real qdz(n/4),qqdz(n/4),qf(n/4)
	real dd(4,4,n/4,n/4),dd1((n/4)*4,(n/4)*4)
	real Pha0T(n/4,n/4),Phb0T(n/4,n/4)
	real Ph0(n/4,n/4),Ph1(n/4,n/4)
	real Ph2(n/4,n/4),Ph3(n/4,n/4),Ph4(n/4,n/4)
	real PhN1(n/4,n/4),PhN1T(n/4,n/4)
	real PhN2(n/4,n/4),PhN3(n/4,n/4)
	real Test(n/4 * n/4)
	real E1, E2, G12, niu1, niu2, D11, D12, D22, D66 , hhh
	real aaa, bbb, qqq


		a=oleole

!		open(106,file='../f.txt',status='unknown')
!		  do i=1,n
!			read(106,*)a(N*N+i)
!		  enddo
!		close(106)

		!open (668,file='A.txt',status='unknown')
		!write (668,1002)A
		!close(668)

1002  FORMAT(1(E12.5,2X))
 111  FORMAT(1(E12.5,2X)) 
      RETURN
	 
		contains
			subroutine multmv(nu,au,b,c)
			real(8) au(nu,nu),c(nu),b(nu)
			integer i,j,k
	
				do i=1,nu
					c(i)=0.
					do ku=1,nu
						c(i)=c(i)+au(i,ku)*b(ku)
					enddo
				end do
			end subroutine multmv

     END subroutine


