    !****************************************************************************
    !
    !  PROGRAM: Beams
    !
    !  PURPOSE:  Entry point for the console application.
    !
    !****************************************************************************
    program Beams
    implicit none
    ! Variables
    integer, parameter :: nt = 12
    integer, parameter :: nl = 18
    integer, parameter :: nk = 21
    integer, parameter :: mt = 6
    integer, parameter :: ml = 3
    integer, parameter :: mn = 2
    Double Precision A(nt,mt)
    Double Precision S(mt,mt)
    Double Precision SAT(mt,nt)
    Double Precision ASAT(nt,nt)
    Double Precision SR(mt,ml)
    Double Precision CORD(mt,ml)
    Integer JTYPE(mt,mt)
    Double Precision GK(nl,nk)
    Integer LSN(ml)
    Double Precision BETA(mt,mn)
    Double Precision GAMMAA(mt)
    integer MCON(mt,mn)
    Double Precision SMA(mt,mn)
    Double Precision SVT(mt)
    Double Precision AREA(mt)
    Double Precision G(mt)
    Double Precision E(mt)
    Double Precision X
    Double Precision Y
    Double Precision Z
    Double Precision LENA
    Integer JFN
    Integer JCT
    Integer NM
    Integer I
    Integer j
    integer flag
    integer ndf
    integer k
    integer ISW
    integer J1
    integer J2
    integer M1
    double precision AK1
    double precision VFACTOR
    double precision V4
    double precision V2
    double precision RL
    double precision GA
    double precision SG
    double precision CG
    double precision XL
    Integer L
    Integer LX
    integer N
    Integer LL
    Integer LJ
    Integer JB
    Integer JE
    Integer KB
    Integer JN2
    Integer JJJ
    Integer JF
    Integer KKK
    Integer KF
    Integer KP
    Integer JJ
    Integer KK
    Integer LSC
    Integer KE
    Integer LB
    Integer LE
    Integer JN1
    Integer LN
    Integer JN
    Integer KJ
    Integer IP1
    Double Precision TEMP
    Integer KL
    Integer ID
    Integer IC
    Integer J3
    Integer II
    Integer III
    Integer NLS
    Integer M
    Integer Ig
    Integer Jg
    write(*,140)nt
140 Format( //,'------------------------------------------------------------------------------------------------------',/&
        '            Solver for Beams',/,&
        '            Code from Harrison 1973 Samples',/,&
        '            Amendments J.M. Nichols 2015',/,&
        '            Problem Size :: ',i4,//,&
        '------------------------------------------------------------------------------------------------------',/)
    open(1, file='a.inp')
    open(2, file='a.rpt', status = 'UNKNOWN')
    flag = 1
10  READ(1,*, err = 1000, end = 1000)JFN
    write(*,15)JFN
15  Format( '    Frame Number for the analysis :: ',i4)
    IF(JFN .le. 0) then
        write(*,*)'Frame Number has to be positive'
20      STOP
    endif
    flag = 2
    READ(1,*, err = 1000, end = 1000)JCT,NM
    write(*,25)JCT,NM
25  Format(  '    Number of Junctions :: ',I4,/,&
        '    Number of Members   :: ',I4,//&
        '------------------------------------------------------------------------------------------------------',//,&
        '                  Frame Details and Design',//,&
        '   i         X          Y         Z          Displacement        Rotation ',/,&
        '               Coordinates                    Restraint          Restraint',//)
    flag = 3
    DO 40 I=1,JCT
        READ(1,*, err = 1000, end = 1000)(CORD(I,J),J=1,3),(JTYPE(I,J),J=1,6)
        write(*,35)i,(CORD(I,J),J=1,3),(JTYPE(I,J),J=1,6)
40  continue
35  Format( i4,2x,'|',3F10.3,'  |',3I6,'  |',3I6)
    flag = 4
    write(*,45)
45  Format( //,'    Zero is fixed point ',//&
        '------------------------------------------------------------------------------------------------------',//,&
        '                  Member Details ',//,&
        '   i    End   End         Moment of I           Area    Youngs    Shear           BETA          Gamma',/&
        '         1     2       Z       Y      X                 Modulus   Modulus        Values',//)
    DO 50 I=1,NM
        READ(1,*, err = 1000, end = 1000)(MCON(I,J),J=1,2),(SMA(I,J),J=1,2),SVT(I),AREA(I),E(I),G(I),(BETA(I,J),J=1,2),GAMMAA(I)
        Write(*,55)i,(MCON(I,J),J=1,2),(SMA(I,J),J=1,2),SVT(I),Area(I),E(I),G(I),(BETA(I,J),J=1,2),GAMMAA(I)
50  continue
55  Format(I4,2x,I4,2x,I4,' | ',F7.5,2x,F7.5,2x,F7.5,' |',1x, F6.3,2x,E9.2,1x,E9.2,1x,E9.2,1x,E9.2,1x, f5.1)
    DO 51 I=1,6
        DO 51 J=1,6
51  S(I,J)=0.
    NDF=0
    DO 52 J=1,JCT
        DO 52 K=1,6
52  NDF=NDF+JTYPE(J,K)
    DO 53 J=1,NDF
        DO 53 K=1,NDF
53  GK(J,K)=0.
    ISW=-1
54  DO 160 I=1,NM
        J1=ABS(MCON(I,1))
        J2=ABS(MCON(I,2))
        IF(J2-J1)58,59,59
58      M1=J1
        J1=J2
        J2=M1
59      X=CORD(J2,1)-CORD(J1,1)
        Y=CORD(J2,2)-CORD(J1,2)
        Z=CORD(J2,3)-CORD(J1,3)
        LENA=SQRT(X*X+Y*Y+Z*Z)
        S(1,1)=E(I)*AREA(I)/LENA
        AK1=4.0*E(I)*SMA(I,1)/LENA
        VFACTOR = BETA(I,1 )*AREA(I)*G(I)*LENA*LENA/(E(I)*SMA(I,1))
        V4=AK1*(VFACTOR+3.0)/(VFACTOR+12.0)
        V2 = AK1 *0.5* (VFACTOR - 6.0) /(VFACTOR +12.0)
        S(2,2)=V4
        S(2,3)=V2
        S(3,2) =V2
        S(3,3)=V4
        AK1 = 4.0*E(I)*SMA(I,2)/LENA
        VFACTOR=BETA(I,2)*AREA(I)*G(I)*LENA*LENA/(E(I)*SMA(I,2))
        V4 = AK1 * (VFACTOR + 3.0)/(VFACTOR +12.0)
        V2=AK1 *0.5*(VFACTOR - 6.0)/(VFACTOR +12.0)
        S(4,4)=V4
        S(4,5)=V2
        S(5,4) =V2
        S(5,5)=V4
        S(6,6)=G(I)*SVT(I)/LENA
        Call PrintS(S,mt,mt,I)
        DO 60 J= 1,12
            DO 60 K=1,6
60      SAT(K,J)=0.
        RL=1.0/LENA
        SAT(1,1)=1.
        SAT(2,2)=-RL
        SAT(3,2)=-RL
        SAT(4,3)=RL
        SAT(5,3)=RL
        SAT(6,4)=-1.
        SAT(4,5) =-1.
        SAT(2,6) =-1.
        SAT(1,7)=-1.
        SAT(2,8)=RL
        SAT(3,8)=RL
        SAT(4,9)=-RL
        SAT(5,9)=-RL
        SAT(6,10)=1.
        SAT(5,11)=-1.
        SAT(3,12)=-1.
        call PrintSAT(SAT,MT,NT,I)
        DO 70 J= 1,12
            DO 70 K=1,12
70      ASAT(J,K)=0.
        GA=GAMMAA(I)*3.14159265/180.0
        SG=SIN(GA)
        CG=COS(GA)
        IF(X+Y) 80,75,80
75      DO 76 J=1,4
            K=3*J-2
            ASAT(K+2,K)=-1.
            ASAT(K,K+1)=SG
            ASAT(K,K+2)=CG
            ASAT(K+1,K+1)=-CG
76      ASAT(K+1,K+2)=SG
        GO TO 86
80      XL=SQRT(X*X+Y*Y)
        Do 85 J=1,4
            K=3*J-2
            ASAT(K,K)=-X/LENA
            ASAT(K,K+1)=Y*CG/XL+X*Z*SG/(LENA*XL)
            ASAT(K,K+2)=-Y*SG/XL+X*Z*CG/(LENA*XL)
            ASAT(K+1,K)=-Y/LENA
            ASAT(K+1,K+1)=(-X*CG/XL)+((Y*Z*SG)/(LENA*XL))
            ASAT(K+1,K+2) = (X*SG/XL)+((Y*Z*CG)/(LENA*XL))
            ASAT(K+2,K)=-Z/LENA
            ASAT(K+2,K+1)=-XL*SG/LENA
85      ASAT(K+2,K+2)= -XL*CG/LENA
86      CONTINUE
        call PrintA(A,nt,mt,I)
        DO 90 J=1,12
            DO 90 K=1,6
                A(j,k) = 0.0
                DO 90 L=1,12
90      A(J,K) = A(J,K)+ASAT(J,L)*SAT(K,L)
        DO 100 J =1,6
            DO 100 K=1,12
                SAT(J,K)=0.
                DO 100 L=1,6
100     SAT(J,K) = SAT(J,K) + S(J,L)*A(K,L)
        IF(ISW)105,105,900
105     DO 110 J=1,12
            DO 110 K=1,12
                ASAT(J,K)=0.
                DO 110 L=1,6
110     ASAT(J,K) = ASAT(J,K) + A(J,L) * SAT(L,K)
        LSN(1)=J1
        LSN(2)=J2
        DO 114 N=1,2
            LL=0
            LJ=LSN(N)-1
            IF(LJ)114,114,112
112         DO 113 J=1,LJ
                DO 113 K=1,6
113         LL=LL+JTYPE(J,K)
114     SAT(N,N)=LL
        DO 160 N=1,4
            JB=6*N-12*((N+1)/2)+7
            JE=JB+5
            KB=6*((N+1)/3-N/4)+1
            KE=KB+5
            LB=JE/6
            LE=KE/6
            JN1=LSN(LB)
            JN2=LSN(LE)
            JJJ=0
            DO 160 J=JB,JE
                JF=J-JB+1
                JJJ=JJJ+JTYPE(JN1,JF)
                KKK=0
                DO 160 K=KB,KE
                    KF=K-KB+1
                    KKK=KKK+JTYPE(JN2,KF)
                    IF(JTYPE(JN1 ,JF)*JTYPE(JN2,KF))160,160,118
118                 JJ=SAT(LB,LB)+JJJ
                    KK=SAT(LE,LE)+KKK
                    GK(JJ,KK)=GK(JJ,KK)+ASAT(J,K)
160 CONTINUE
    
    CAll PrintGK(GK,nl,nk)
2011 format (12(1x,E12.5))
    write(*,*)'hereA'
    write(*,*)'Here ISW ',ISW
    IF(ISW)164,1040,1040
164 L=NDF
    LSC=L
    Flag = 5
460 READ(1,*, ERR = 1000, END = 1000)KK,LN
    LSC=LSC+1
    IF(KK)570,470,470
470 DO 480 I=1,L
480 GK(I,LSC)=0.0
    I = LSC-L
    LSN(i)=KK
    DO 520 I=1,LN
        FLAG = 6
        READ(1,*, ERR = 1000, END = 1000)JN,(SR(J,1),J=1,6)
        LL=0
        LJ=JN-1
        IF(LJ) 493,493,492
492     DO 500 J=1,LJ
            DO 500 K=1,6
500     LL=LL+JTYPE(J,K)
493     DO 520 K=1,6
            IF(JTYPE(JN,K))520,520,510
510         LL=LL+1
            GK(LL,LSC)=SR(K,1)
520 CONTINUE
    GO TO 460
570 KJ=LSC-1
    NLS=KJ-L
    DO 680 I=1,L
        IP1=I+1
        TEMP=ABS(GK(I,I))
        K=I
        DO 590 J=I,L
            IF(ABS(GK(J,I))-TEMP)590,590,585
585         K=J
            TEMP=ABS(GK(J,I))
590     CONTINUE
        IF(K-1)600,620,600
600     DO 610 J=I,KJ
            TEMP=GK(I,J)
            GK(I,J)=GK(K,J)
610     GK(K,J)=TEMP
620     IF(GK(I,I))665,1060,665
665     TEMP=1.0/GK(I,I)
        DO 650 J=I,KJ
650     GK(I,J)=GK(I,J)*TEMP
        DO 680 J=1,L
            IF(I-J)660,680,660
660         TEMP=GK(J,I)
            DO 670 K=IP1,KJ
670         GK(J,K)=GK(J,K)-TEMP*GK(I,K)
680 CONTINUE
    DO 820 I=1,NLS
        KL=L+I
        write(*,700)LSN(I)
700     FORMAT(//'-------------------------------------------------------------------------------------------------------',//&
            20x,36H Deformations caused by Load Set No.,I3//)
        LL=0
        write(*,*)  ' NODE X-MOVEMENT Y-MOVEMENT Z-MOVEMENT X-ROTAT10N Y-ROTATION Z-ROTATION'
        DO 820 J=1,JCT
            DO 740 K=1,6
                IF(JTYPE(J,K))730,730,720
720             LL=LL+1
                SR(K,1)=GK(LL,KL)
                GO TO 740
730             SR(K,1)=-0.
740         CONTINUE
            Write(*,810)J,(SR(K,1),K= 1,6)
810         FORMAT(/I3,1X,3F11.6,1X,3F11.6)
820 CONTINUE
    ISW=1
    GO TO 54
900 ID=0
    M=0
    L=NDF
903 JJ=J1-1
    IF(JJ)912,912,905
905 DO 910 J=1,JJ
        DO 910 K=1,6
910 ID = ID + JTYPE(J,K)
912 DO 965 IC=1,6
        IF(JTYPE(J1,IC))915,915,920
915     DO 916 N=1,NLS
916     ASAT(M + IC,N) =0.0
        GO TO 965
920     DO 925 N=1,NLS
925     ASAT(M+IC,N)=GK(ID+1,L+N)
        ID=ID+1
965 CONTINUE
    IF(M)980,980,990
980 M=6
    J3=J1
    J1=J2
    ID=0
    GO TO 903
990 DO 1100 II=1,6
        DO 1100 J=1,NLS
            SR(II,J)=0.
            DO 1100 K=1,12
1100 SR(II,J)=SR(II,J)+SAT(II,K)*ASAT(K,J)
    Write(*,1010)I
1010 FORMAT(//11H MEMBER NO.,I3/)
    Write(*,*)' L/SET AX.TENSION MAB-Z-AX1S MBA-Z-AXIS MAB-Y-AXIS MBA-Y-AXIS   TORQUE'
    DO 1020 J=1,NLS
1020 PRINT 1030,LSN(J),(SR(II,J),II=1,6)
1030 FORMAT(/I4,2X,5F11.4,' | ',F10.4)
    N=4
    J=12
    write(*,*)'here'
    GO TO 160
1040 PRINT 1050,JFN
1050 FORMAT(//33H ANALYSIS COMPLETED FOR FRAME NO.,I3/)
    CAll PrintGK(GK,nl,nk)
2010 format (12(1x,E12.5))
    GO TO 10
1060 PRINT 1070,JFN
1070 FORMAT(/48H ZERO DIVISION IN EQUATION SOLUTION IN FRAME No.,I3)
    Go to 10
    stop
1000 call error(flag)
    end program Beams

    subroutine error(i)
    implicit none
    integer i
    write(*,10)i
10  format(    'Error in read ',i4)
    return
    end

    subroutine PrintS(S,n,m,NM)

    implicit none
    integer m
    integer n
    integer i,j,NM
    double precision S(n,m)
    !write(*,1000)NM
    write(2,1000)NM
    do 10 i = 1, n
        !   write(*,100)(S(i,j),j=1,6)
        write(2,100)(S(i,j),j=1,6)
10  end do
100 Format(6F15.3)
1000 Format(//'-----------------------------------------------------------------------------------',//,&
        '  Print S Matrix ',/,'  Member :: ',I4,/)

    return
    end subroutine PrintS

    subroutine PrintSAT(S,n,m,NM)
    implicit none
    integer m
    integer n
    integer i,j,NM
    double precision S(n,m)
    !write(*,1000)NM
    write(2,1000)NM
    do 10 i = 1, n
        !write(*,100)(S(i,j),j=1,12)
        write(2,100)(S(i,j),j=1,12)
10  end do
100 Format(12F10.3)
1000 Format(//'-----------------------------------------------------------------------------------',//,&
        '  Print SAT  Member :: ',I4,/)
    return
    end subroutine PrintSAT
    
    !-----------------------------------------------------------------------------------------------------------------------
    !
    !   Print A Matrix
    !
    !----------------------------------------------------------------------------------------------------------------------

    subroutine PrintA(S,n,m,NM)
    implicit none
    integer m
    integer n
    integer i,j,NM
    double precision S(n,m)
    ! write(*,1000)NM
    write(2,1000)NM
    do 10 i = 1, n
        !write(*,100)(S(i,j),j=1,6)
        write(2,100)(S(i,j),j=1,6)
10  end do
100 Format(6F12.3)
    1000 Format(//'-----------------------------------------------------------------------------------',//,&
    '  A Matrix',/,&    
    '  Print A  Member :: ',I4,/)
    return
    end subroutine PrintA
    
    
    !-----------------------------------------------------------------------------------------------------------------------
    !
    !   Print A Matrix
    !
    !----------------------------------------------------------------------------------------------------------------------

    
    Subroutine PrintGK(GK,m,n)
    
    implicit none
    double precision GK(m,n)
    integer m,n,ig,jg
    
    write(2,100)
    100 Format(//,'--------------------------------------------------------------------------------------------------------------------',//&
                 '  Stiffenss Matrix K',//)
    
    do 2001 ig =1,12
        write(2,2011)(GK(ig,jg),jg=1,12)
2001 end do
     
2011 format (12(1x,F11.0))
     
     return
     end subroutine