!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	Timeshenko分層彈塑性梁單元材料非線性計算						  !
!	 材料采用彈塑性，混凝土材料同樣采用改種材料						  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



$DEBUG
!	NPOIN(結點總數),NELEM(單元總數),NBOUN(邊界點總數),NLAYR(分層數),NLOAD(),
!				NPROP(完全確定材料需要的參數),NNODE(每個單元的結點數),IINCS(),
!     IITER,KRESL,NCHEK,TOLER,NALGO(用以識別解法的指示變量),NSVAB,NDOFN(每個結點的自由度數),
!     NINCS(施加荷載增量的數目),NEVAB,NITER,NOUTP,FACTO,PVALU

!     IEVAB(結點自由度變量)


!	PROPS(5,25)(材料數組),COORD(26)(結點坐標),LNODS(25,2)(單元拓撲信息),IFPRE(52),
!     FIXED(52),TLOAD(25,4),RLOAD(25,4),ELOAD(25,4),
!     MATNO(25)(每個單元對應的材料號),STRES(25,2),PLAST(250),XDISP(52),
!     TDISP(26,2),TREAC(26,2),ASTIF(52,52),ASLOD(52),
!     REACT(52),FRESV(1352),PEFIX(52),ESTIF(4,4),
!     STRSL(250,2)
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	

	MODULE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::NPOIN,NELEM,NBOUN,NLAYR,NLOAD,NPROP,NNODE,IINCS,
     $			IITER,KRESL,NCHEK,NALGO,NSVAB,NDOFN,
     $			NINCS,NEVAB,NITER,NOUTP,NMATS
      DOUBLEPRECISION::FACTO,PVALU,TOLER


	DOUBLEPRECISION,ALLOCATABLE::PROPS(:,:),COORD(:),TLOAD(:,:) 
	DOUBLEPRECISION,ALLOCATABLE::RLOAD(:,:),ELOAD(:,:)
	DOUBLEPRECISION,ALLOCATABLE::STRES(:,:),PLAST(:),XDISP(:)
	DOUBLEPRECISION,ALLOCATABLE::TDISP(:,:),TREAC(:,:)
	DOUBLEPRECISION,ALLOCATABLE::ASTIF(:,:),ASLOD(:),REACT(:)
	DOUBLEPRECISION,ALLOCATABLE::FRESV(:),PEFIX(:)
	DOUBLEPRECISION,ALLOCATABLE::STRSL(:,:)
	DOUBLEPRECISION::ESTIF(4,4)
	INTEGER,ALLOCATABLE::LNODS(:,:),IFPRE(:),FIXED(:),MATNO(:)

	END MODULE CONFIG_ARRAY

!*============================================================*!　
!*												 　			*!
!*============================================================*! 	


	PROGRAM BEAM_FENGCHENG
	USE CONFIG_ARRAY
	IMPLICIT NONE
	CHARACTER PN*40,FN_1*12,FN_5*12,FN_6*12,FN_7*12
	WRITE(*,'(/A)')' 鋼筋混凝土平面Timeshenko分層梁單元'
      write(*,'(A)')'	         承載能力計算程序'
	WRITE(*,'(/A)') ' 輸入計算問題名(PN):'
	READ(*,'(A)')PN

	CALL FNAME(PN,'.DAT',FN_5)
	WRITE(*,'(/2A)') '  輸入數據文件名為 ',FN_5
	CALL FNAME(PN,'.OUT',FN_6)
	WRITE(*,'(/2A)') '  輸入數據文件名為 ',FN_6
	CALL FNAME(PN,'.PFX',FN_7)
	WRITE(*,'(/2A)') '  荷載撓度關系輸出數據文件名為 ',FN_7
 	OPEN(1,STATUS='NEW',FORM='UNFORMATTED')
	OPEN(5,FILE=FN_5)
	OPEN(6,FILE=FN_6)
	OPEN(7,FILE=FN_7)
	CALL DATAIN
	CALL INITAL
	DO IINCS=1,NINCS
		CALL INCLOD
		DO IITER=1,NITER
			WRITE(*,*)'荷載步:' ,IINCS,'	迭代步:', IITER
			CALL NONAL
			IF (KRESL.EQ.1) CALL STIFBL
			CALL ASSEMB
			IF(KRESL.EQ.1)CALL GREDUC
			IF(KRESL.EQ.2) CALL RESOLV
			CALL BAKSUB
			CALL REFORBL
			CALL CONUND
			IF (NCHEK.EQ.0) GOTO 20
			IF(IITER.EQ.1.AND.NOUTP.EQ.1)CALL RESULT
			IF(NOUTP.EQ.2)CALL RESULT
		ENDDO
		WRITE(6,900)
900		FORMAT(1HO,5X,'SOLUTION NOT COVERGED')
		STOP
20		CALL RESULT
	ENDDO

	CLOSE(5)
	CLOSE(6)
	CLOSE(7)
	DEALLOCATE(PROPS,COORD,TLOAD )
	DEALLOCATE(RLOAD,ELOAD)
	DEALLOCATE(STRES,PLAST,XDISP)
	DEALLOCATE(TDISP,TREAC)
	DEALLOCATE(ASTIF,ASLOD,REACT)
	DEALLOCATE(FRESV,PEFIX)
	DEALLOCATE(STRSL)
	DEALLOCATE(LNODS,IFPRE,FIXED,MATNO)

	STOP
	END	PROGRAM BEAM_FENGCHENG

!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE DATAIN
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::IMATS,JMATS,IPROP,IELEM,JELEM,INODE
	INTEGER::IPOIN,JPOIN,ISVAB,IBOUN,NODFX,IDOFN,NPOSN,IEVAB
	INTEGER::ICODE(2)
	DOUBLEPRECISION::VALUE(2) 
!	CHARACTER TITLE*140 
!
!	READ(5,*) TITLE
!	WRITE(6,*) TITLE
	READ(5,*) NPOIN,NELEM,NBOUN,NMATS,NPROP,NNODE,NINCS,
     $NALGO,NDOFN,NLAYR
	WRITE(6,905)NPOIN,NELEM,NBOUN,NMATS,NPROP,NNODE,
     $NINCS,NALGO,NDOFN,NLAYR
905	FORMAT(//1X,'NPOIN=',I5,3X,'NELEM=',I5,3X,'NBOUN=',I5,3X,
     $		    'NMATS=',I5//1X,'NPROP=',I5,3X,'NNODE=',I5,3X,
     $			'NINCS=',I5,3X,'NALGO=',I5//1X,'NDOFN=',I5,3X,
     $			'NLAYR=',I5)
	NEVAB=NDOFN*NNODE
	NSVAB=NDOFN*NPOIN

	ALLOCATE (PROPS(NMATS,NPROP))
	ALLOCATE(LNODS(NELEM,NNODE),MATNO(NELEM))
	ALLOCATE(COORD(NPOIN))
	ALLOCATE(IFPRE(NSVAB),PEFIX(NSVAB),FIXED(NSVAB))
	ALLOCATE(RLOAD(NELEM,NEVAB))
	ALLOCATE(PLAST(NLAYR*NELEM),STRSL(NLAYR*NELEM,2))
	ALLOCATE(STRES(NELEM,NDOFN))
	ALLOCATE(ELOAD(NELEM,NEVAB),TLOAD(NELEM,NEVAB))
	ALLOCATE(TDISP(NPOIN,NDOFN),TREAC(NPOIN,NDOFN))
	ALLOCATE(ASLOD(NSVAB),ASTIF(NSVAB,NSVAB))
	ALLOCATE(FRESV((NSVAB-1)*NSVAB/2))
	ALLOCATE(REACT(NSVAB),XDISP(NSVAB)) 



	WRITE(6,910)
910	FORMAT(5X,'MATERIAL PROPERTIES')
 
	DO IMATS=1,NMATS
		READ(5,*) JMATS,(PROPS(JMATS,IPROP),IPROP=1,NPROP)
		WRITE(6,915)JMATS,(PROPS(JMATS,IPROP),IPROP=1,NPROP)
	ENDDO
915	FORMAT(I10,4F15.5/10X,4F15.5/10X,4F15.5/10X,4F15.5/10X,4F15.5)
	WRITE(6,920)
920	FORMAT(3X,'EL NODFS MAT.')

	DO IELEM=1,NELEM
	 READ(5,*) JELEM,(LNODS(JELEM,INODE),INODE=1,NNODE),MATNO(JELEM)
	 WRITE(6,925)JELEM,(LNODS(JELEM,INODE),INODE=1,NNODE),MATNO(JELEM)
	ENDDO
925	FORMAT(4I5)
	WRITE(6,930)
930	FORMAT(5X,'NODE',5X,'COORD.')

	DO IPOIN=1,NPOIN
		READ(5,*)JPOIN,COORD(JPOIN)
		WRITE(6,935)JPOIN,COORD(JPOIN)
	ENDDO
935	FORMAT(I10,F15.5)
	
	DO ISVAB=1,NSVAB
		IFPRE(ISVAB)=0.
		PEFIX(ISVAB)=0.0
	ENDDO

	IF(NDOFN.EQ.1)WRITE(6,940)
940	FORMAT(1X,'RES.NODE',2X,'CODE',3X,'PRES.VALUES')
	IF(NDOFN.EQ.2)WRITE(6,945)
945	FORMAT(1X,'RES.NODE',2X,'CODE',3X,'PRES.VALUES',2X,
     $		'CODE',3X,'PRES.VALUES')

	DO IBOUN=1,NBOUN
		READ(5,*) NODFX,(ICODE(IDOFN),VALUE(IDOFN),IDOFN=1,NDOFN)
		WRITE(6,950)NODFX,(ICODE(IDOFN),VALUE(IDOFN),IDOFN=1,NDOFN)
950	FORMAT(I10,2(I5,F15.5))
		NPOSN=(NODFX-1)*NDOFN
		DO IDOFN=1,NDOFN
			NPOSN=NPOSN+1
			IFPRE(NPOSN)=ICODE(IDOFN)
 			PEFIX(NPOSN)=VALUE(IDOFN)
		ENDDO
	ENDDO
	WRITE(6,955)
955	FORMAT(3X,'ELEMENT',10X,'NODAL LOADS')
	DO IELEM=1,NELEM
		DO IEVAB=1,NEVAB
			RLOAD(IELEM,IEVAB)=0.0
		ENDDO
	ENDDO
	
	DO IELEM=1,NELEM
 		READ(5,*)JELEM,(RLOAD(JELEM,IEVAB),IEVAB=1,NEVAB)
	ENDDO

	DO IELEM=1,NELEM
		WRITE(6,960)IELEM,(RLOAD(IELEM,IEVAB),IEVAB=1,NEVAB)
	ENDDO
960	FORMAT(I10,5F15.5)

	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE INITAL
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::IXiaBiao,IELEM,IDOFN,IEVAB,IPOIN
	DO IXiaBiao=1,NLAYR*NELEM
		PLAST(IXiaBiao)=0.0
		STRSL(IXiaBiao,1)=0.0
		STRSL(IXiaBiao,2)=0.0
	ENDDO
	DO IELEM=1,NELEM
!		PLAST(IELEM)=0.0
		DO IDOFN=1,NDOFN
			STRES(IELEM,IDOFN)=0.0
		ENDDO
		DO IEVAB=1,NEVAB
			ELOAD(IELEM,IEVAB)=0.0
			TLOAD(IELEM,IEVAB)=0.0
		ENDDO
	ENDDO
	DO IPOIN=1,NPOIN
		DO IDOFN=1,NDOFN
			TDISP(IPOIN,IDOFN)=0.0
			TREAC(IPOIN,IDOFN)=0.0
		ENDDO
	ENDDO
	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE NONAL
	USE CONFIG_ARRAY
	IMPLICIT NONE

	INTEGER::ISVAB
	KRESL=2
	IF(NALGO.EQ.1)KRESL=1
	IF(NALGO.EQ.2)KRESL=1
	IF(NALGO.EQ.3.AND.IINCS.EQ.1.AND.IITER.EQ.1)KRESL=1
	IF(NALGO.EQ.4.AND.IITER.EQ.1)KRESL=1
	IF(NALGO.EQ.5.AND.IITER.EQ.2)KRESL=1
	IF(IITER.EQ.1.OR.NALGO.EQ.1)GOTO 20
	DO ISVAB=1,NSVAB
		FIXED(ISVAB)=0.0
	ENDDO
	RETURN
20	DO ISVAB=1,NSVAB
		FIXED(ISVAB)=PEFIX(ISVAB)*FACTO
	ENDDO
	RETURN
	END	

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     荷載增量
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE INCLOD
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::IELEM,IEVAB
	READ(5,*)NITER,NOUTP,FACTO,TOLER
	WRITE(6,905)IINCS,NITER,NOUTP,FACTO,TOLER
905	FORMAT(5X,'IINCS=',I5,3X,'NITER= ',I5,3X,' NOUTP=',I5,	
     $	3X,'FACTO=',E14.6,3X,'TOLER=',E14.6)
	DO IELEM=1,NELEM
		DO IEVAB=1,NEVAB
			ELOAD(IELEM,IEVAB)=ELOAD(IELEM,IEVAB)+
     $					RLOAD(IELEM,IEVAB)*FACTO
			TLOAD(IELEM,IEVAB)=TLOAD(IELEM,IEVAB)+
     $					RLOAD(IELEM,IEVAB)*FACTO
		ENDDO
	ENDDO
	RETURN
	END

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     剛度矩陣的組集
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE ASSEMB 
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::ISVAB,JSVAB,IELEM,INODE,NODEI,NODEJ,NROWS,NROWE
	INTEGER::JNODE,NCOLS,NCOLE,I,J,IDOFN,JDOFN
 	REWIND 1
	DO ISVAB=1,NSVAB
		ASLOD(ISVAB)=0.0
	ENDDO

	IF(KRESL.EQ.2) GOTO 30
	DO ISVAB=1,NSVAB
		DO JSVAB=1,NSVAB
			ASTIF(ISVAB,JSVAB)=0.0
		ENDDO
	ENDDO
30	CONTINUE
C
C *** ASSEMBLE THE ELEMNT LAODS
C
	DO IELEM=1,NELEM
		DO I=1,4
			DO J=1,4
				ESTIF(I,J)=0.0	
			ENDDO
		ENDDO

			READ(1)ESTIF

 		DO INODE=1,NNODE
 			NODEI=LNODS(IELEM,INODE)
 			DO IDOFN=1,NDOFN
 				NROWS=(NODEI-1)*NDOFN+IDOFN
 				NROWE=(INODE-1)*NDOFN+IDOFN
 				ASLOD(NROWS)=ASLOD(NROWS)+ELOAD(IELEM,NROWE)
C
C *** ASSEMBLE THE ELEMNT STIFFNESS MATRICES 
C
 				IF(KRESL.EQ.2)GOTO 40
 				DO JNODE=1,NNODE
 					NODEJ=LNODS(IELEM,JNODE)
 					DO  JDOFN=1,NDOFN
 						NCOLS=(NODEJ-1)*NDOFN+JDOFN
 						NCOLE=(JNODE-1)*NDOFN+JDOFN
 						ASTIF(NROWS,NCOLS)=ASTIF(NROWS,NCOLS)
     $					+ESTIF(NROWE,NCOLE)
					ENDDO
				ENDDO
			ENDDO
40		ENDDO
	ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!調試語句
!	OPEN(2,FILE='zg.dat')
!	DO I=1,NROWS
!			WRITE(2,111)(ASTIF(I,J),J=1,NCOLS)
!	ENDDO
!111	FORMAT(5X,22F25.5)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!調試語句
	RETURN
	END
!*============================================================*!　
!*												 　			*!
!*============================================================*! 	
	SUBROUTINE STIFBL
	USE CONFIG_ARRAY
	IMPLICIT NONE
	INTEGER::IELEM,LPROP,NODE1,NODE2,ISTIF,JSTIF
	DOUBLEPRECISION::EIVAL,SVALU,VALU1,VALU2,VALU3,VALU4,HARDS
	DOUBLEPRECISION::ELENG
 	REWIND 1
	DO IELEM=1,NELEM
		LPROP=MATNO(IELEM)
		CALL LAYER(IELEM,EIVAL,SVALU)
	WRITE(*,*)EIVAL,SAVLU
		HARDS=PROPS(LPROP,4)
		NODE1=LNODS(IELEM,1)
		NODE2=LNODS(IELEM,2)
		ELENG=ABS(COORD(NODE2)-COORD(NODE1))
!		IF(PLAST(IELEM).NE.0) EIVAL=EIVAL*(1.0-EIVAL/(EIVAL+HARDS))
		VALU1=0.5*SVALU
		VALU2=SVALU/ELENG
		VALU3=EIVAL/ELENG
		VALU4=0.25*SVALU*ELENG
		ESTIF(1,1)=VALU2
		ESTIF(1,2)=VALU1
		ESTIF(1,3)=-VALU2
		ESTIF(1,4)=VALU1
		ESTIF(2,2)=VALU3+VALU4