)=SS(I,J)+DD(I,K)*BB(K,J)Р DO 20 I=1,6Р DO 20 J=1,6Р AKE(I,J)=0.0Р DO 20 K=1,3Р 20 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTР RETURN Р?ENDРcРc to form banded global stiffness matrix?Р?SUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE,Р & PR,PT,AK)Р DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)Р N_DOF=2*N_NODEР DO 40 I=1,N_DOFР DO 40 J=1,N_BANDР 40 AK(I,J)=0Р DO 50 IE=1,N_ELEР CALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)Р DO 50 I=1,3Р DO 50 II=1,2Р IH=2*(I-1)+IIР IDH=2*(IJK_ELE(IE,I)-1)+IIР DO 50 J=1,3Р DO 50 JJ=1,2Р IL=2*(J-1)+JJР IZL=2*(IJK_ELE(IE,J)-1)+JJР IDL=IZL-IDH+1Р IF(IDL.LE.0) THEN Р ELSEР AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL) Р END IFР 50?CONTINUEР?RETURNР?ENDРcРc to calculate the area of element Р SUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)