error775

I have a program with error775.please help me to run it.what should I do?

PROGRAM p72
!------------------------------------------------------------------------- ! Program 7.2 Plane or axisymmetric analysis of steady seepage using ! 4-node rectangular quadrilaterals. Mesh numbered ! in x(r)- or y(z)- direction. !-------------------------------------------------------------------------

IMPLICIT NONE INTEGER,PARAMETERiwp=SELECTED_REAL_KIND(15) INTEGERfixed_freedoms,i,iel,k,loaded_nodes,nci,ndim=2,nels,neq,nip=4, & nod=4,nn,np_types,nxe,nye,determinant REAL(iwp)det,one=1.0_iwp,penalty=1.0e20_iwp,zero=0.0_iwp
CHARACTER(LEN=15)dir,element='quadrilateral',type_2d !-----------------------dynamic arrays------------------------------------ INTEGER,ALLOCATABLEetype(:),g_num(:,:),kdiag(:),node(:),num(:) REAL(iwp),ALLOCATABLEcoord(:,:),der(:,:),deriv(:,:),disps(:),fun(:), & gc(:),g_coord(:,:),jac(:,:),kay(:,:),kc(:,:),kv(:),kvh(:),loads(:), & points(:,:),prop(:,:),value(:),weights(:),x_coords(:),y_coords(:) !-----------------------input and initialisation-------------------------- OPEN(10,FILE='fe95.dat') OPEN(11,FILE='fe95.res') READ(10,)type_2d,dir,nxe,nye,np_types CALL mesh_size(element,nod,nels,nn,nxe,nye) neq=nn ALLOCATE(points(nip,ndim),g_coord(ndim,nn),coord(nod,ndim), & jac(ndim,ndim),weights(nip),der(ndim,nod),deriv(ndim,nod), & kc(nod,nod),num(nod),g_num(nod,nels),kay(ndim,ndim),etype(nels), & x_coords(nxe+1),y_coords(nye+1),prop(ndim,np_types),gc(ndim),fun(nod), & kdiag(neq),loads(0:neq),disps(0:neq)) READ(10,)prop etype=1 IF(np_types>1)READ(10,)etype READ(10,)x_coords,y_coords
!-----------------------loop the elements to find global arrays sizes----- kdiag=0 elements_1: DO iel=1,nels CALL geom_rect(element,iel,x_coords,y_coords,coord,num,dir) g_num(:,iel)=num g_coord(:,num)=TRANSPOSE(coord) CALL fkdiag(kdiag,num) END DO elements_1 CALL mesh(g_coord,g_num,12) DO i=2,neq kdiag(i)=kdiag(i)+kdiag(i-1) END DO ALLOCATE(kv(kdiag(neq)),kvh(kdiag(neq))) WRITE(11,'(2(A,I5))') & 'There are',neq,' equations and the skyline storage is',kdiag(neq) CALL sample(element,points,weights) !-----------------------global conductivity matrix assembly--------------- kv=zero gc=one elements_2: DO iel=1,nels kay=zero DO i=1,ndim kay(i,i)=prop(i,etype(iel)) END DO num=g_num(:,iel) coord=TRANSPOSE(g_coord(:,num)) kc=zero gauss_pts_1: DO i=1,nip CALL shape_der(der,points,i) CALL shape_fun(fun,points,i) jac=MATMUL(der,coord) det=determinant(jac) CALL invert(jac) deriv=MATMUL(jac,der) IF(type_2d=='axisymmetric')gc=MATMUL(fun,coord) kc=kc+MATMUL(MATMUL(TRANSPOSE(deriv),kay),deriv)detweights(i)gc(1) END DO gauss_pts_1 CALL fsparv(kv,kc,num,kdiag) END DO elements_2 kvh=kv !-----------------------specify boundary values--------------------------- loads=zero READ(10,)loaded_nodes IF(loaded_nodes/=0)READ(10,)(k,loads(k),i=1,loaded_nodes) READ(10,)fixed_freedoms IF(fixed_freedoms/=0)THEN ALLOCATE(node(fixed_freedoms),value(fixed_freedoms)) READ(10,*)(node(i),value(i),i=1,fixed_freedoms) kv(kdiag(node))=kv(kdiag(node))+penalty loads(node)=kv(kdiag(node))value END IF !-----------------------equation solution--------------------------------- CALL sparin(kv,kdiag) CALL spabac(kv,loads,kdiag)
!-----------------------retrieve nodal net flow rates--------------------- CALL linmul_sky(kvh,loads,disps,kdiag) WRITE(11,'(/A)')' Node Total Head Flow rate' DO k=1,nn WRITE(11,'(I5,2E12.4)')k,loads(k),disps(k) END DO disps(0)=zero WRITE(11,'(/A)')' Inflow Outflow' WRITE(11,'(5X,2E12.4)') & SUM(disps,MASK=disps>zero),SUM(disps,MASK=disps<zero) READ(10,)nci IF(nod==4)CALL contour(loads,g_coord,g_num,nci,13) STOP

CONTAINS SUBROUTINE mesh_size(element,nod,nels,nn,nxe,nye,nze) ! ! This subroutine returns the number of elements (nels) and the number ! of nodes (nn) in a 2-d geometry-created mesh. ! IMPLICIT NONE CHARACTER(LEN=15),INTENT(IN)::element INTEGER,INTENT(IN)nod,nxe,nye INTEGER,INTENT(IN),OPTIONALnze INTEGER,INTENT(OUT)::nels,nn IF(element=='triangle')THEN nels=nxenye2 IF(nod3)nn=(nxe+1)*(nye+1) IF(nod6)nn=(2nxe+1)(2nye+1) IF(nod==10)nn=(3nxe+1)(3nye+1) IF(nod15)nn=(4nxe+1)(4*nye+1) ELSE IF(element'quadrilateral')THEN nels=nxenye IF(nod==4)nn=(nxe+1)(nye+1) IF(nod8)nn=(2nxe+1)(nye+1)+(nxe+1)*nye IF(nod9)nn=(2nxe+1)(2nye+1) ELSE IF(element=='hexahedron')THEN nels=nxenyenze IF(nod==8)nn=(nxe+1)(nye+1)(nze+1) IF(nod==14)nn=4nxenyenze+2*(nxenye+nyenze+nzenxe)+nxe+nye+nze+1 IF(nod==20)nn=((2nxe+1)(nze+1)+(nxe+1)nze)(nye+1)+ & (nxe+1)(nze+1)*nye END IF RETURN END SUBROUTINE mesh_size END

SUBROUTINE geom_rect(element,iel,x_coords,y_coords,coord,num,dir) ! ! This subroutine forms the coordinates and connectivity for a ! rectangular mesh of rt. angled triangular elements (3, 6, 10 or 15-node) ! or quadrilateral elements (4, 8 or 9-node) counting in the ! x- or y-dir. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN)::x_coords(:),y_coords(:) REAL(iwp),INTENT(OUT)::coord(:,:) CHARACTER(LEN=15),INTENT(IN)::element CHARACTER(LEN=1),INTENT(IN)::dir INTEGER,INTENT(IN)::iel INTEGER,INTENT(OUT)num(:) INTEGERip,iq,jel,fac1,nod,nxe,nye REAL(iwp)::pt5=0.5_iwp,two=2.0_iwp,d3=3.0_iwp nxe=UBOUND(x_coords,1)-1 nod=UBOUND(num,1) IF(element=='triangle')THEN nye=(UBOUND(y_coords,1)-1)2 IF(dir=='x'.OR.dir=='r')THEN jel=2nxe*((iel-1)/(2nxe)) ip=(iel-jel+1)/2 iq=2((iel-1)/(2nxe)+1)-1+((iel/2)2)/iel ELSE
jel=(iel-1)/nye ip=jel+1 iq=iel-nyejel END IF SELECT CASE(nod) CASE(3) IF(MOD(iq,2)/=0)THEN IF(dir=='x'.OR.dir=='r')THEN num(1)=(nxe+1)(iq-1)/2+ip num(2)=num(1)+1
num(3)=(nxe+1)(iq+1)/2+ip ELSE num(1)=(ip-1)(nye+2)/2+(iq+1)/2 num(2)=num(1)+(nye+2)/2 num(3)=num(1)+1 END IF ! coord(1,1)=x_coords(ip) coord(1,2)=y_coords((iq+1)/2) coord(2,1)=x_coords(ip+1)
coord(2,2)=y_coords((iq+1)/2) coord(3,1)=x_coords(ip)
coord(3,2)=y_coords((iq+3)/2) ELSE IF(dir=='x'.OR.dir=='r')THEN num(1)=(nxe+1)iq/2+ip+1
num(2)=num(1)-1
num(3)=(nxe+1)(iq-2)/2+ip+1 ELSE num(1)=ip*(nye+2)/2+(iq+2)/2 num(2)=(ip-1)(nye+2)/2+(iq+1)/2+1 num(3)=num(1)-1 END IF ! coord(1,1)=x_coords(ip+1) coord(1,2)=y_coords((iq+2)/2) coord(2,1)=x_coords(ip)
coord(2,2)=y_coords((iq+2)/2) coord(3,1)=x_coords(ip+1) coord(3,2)=y_coords(iq/2) END IF CASE(6) IF(MOD(iq,2)/=0)THEN IF(dir=='x'.OR.dir=='r')THEN num(1)=(iq-1)(2nxe+1)+2ip-1 num(2)=num(1)+1 num(3)=num(1)+2 num(4)=(iq-1)(2nxe+1)+2nxe+2ip+1 num(5)=(iq+1)(2nxe+1)+2ip-1 num(6)=num(4)-1 ELSE num(1)=2(nye+1)(ip-1)+iq num(2)=2(nye+1)(ip-1)+nye+1+iq num(3)=2(nye+1)ip+iq num(4)=num(2)+1 num(5)=num(1)+2 num(6)=num(1)+1 END IF ! coord(1,1)=x_coords(ip) coord(1,2)=y_coords((iq+1)/2) coord(3,1)=x_coords(ip+1)
coord(3,2)=y_coords((iq+1)/2) coord(5,1)=x_coords(ip)
coord(5,2)=y_coords((iq+3)/2) ELSE IF(dir=='x'.OR.dir=='r')THEN num(1)=iq(2nxe+1)+2ip+1 num(2)=num(1)-1 num(3)=num(1)-2 num(4)=(iq-2)(2nxe+1)+2nxe+2ip+1 num(5)=(iq-2)(2nxe+1)+2ip+1 num(6)=num(4)+1 ELSE num(1)=2(nye+1)ip+iq+1 num(2)=2(nye+1)(ip-1)+nye+iq+2 num(3)=2(nye+1)(ip-1)+iq+1 num(4)=num(2)-1 num(5)=num(1)-2 num(6)=num(1)-1 END IF ! coord(1,1)=x_coords(ip+1) coord(1,2)=y_coords((iq+2)/2) coord(3,1)=x_coords(ip)
coord(3,2)=y_coords((iq+2)/2) coord(5,1)=x_coords(ip+1) coord(5,2)=y_coords(iq/2) END IF coord(2,:)=pt5(coord(1,:)+coord(3,:)) coord(4,:)=pt5*(coord(3,:)+coord(5,:)) coord(6,:)=pt5*(coord(5,:)+coord(1,:)) CASE(10) IF(MOD(iq,2)/=0)THEN IF(dir=='x'.OR.dir=='r')THEN num(1)=(iq-1)/2*(3nxe+1)3+3ip-2 num(2)=num(1)+1 num(3)=num(1)+2 num(4)=num(1)+3 num(5)=(iq-1)/2(3nxe+1)3+3nxe+1+3ip num(6)=(iq-1)/2*(3nxe+1)3+6nxe+2+3ip-1 num(7)=(iq-1)/2*(3nxe+1)3+9nxe+3+3ip-2 num(8)=num(6)-1 num(9)=num(5)-2 num(10)=num(9)+1 ELSE num(1)=(9*(nye-2)/2+12)(ip-1)+3(iq-1)/2+1 num(2)=(9*(nye-2)/2+12)(ip-1)+3(nye-2)/2+4+3*(iq-1)/2+1 num(3)=(9*(nye-2)/2+12)(ip-1)+3(nye-2)+8+3*(iq-1)/2+1 num(4)=(9*(nye-2)/2+12)(ip-1)+9(nye-2)/2+12+3*(iq-1)/2+1 num(5)=num(3)+1 num(6)=num(2)+2 num(7)=num(1)+3 num(8)=num(1)+2 num(9)=num(1)+1 num(10)=num(2)+1 END IF ! coord(1,1)=x_coords(ip) coord(2,1)=x_coords(ip)+(x_coords(ip+1)-x_coords(ip))/d3 coord(3,1)=x_coords(ip)+two*(x_coords(ip+1)-x_coords(ip))/d3 coord(4,1)=x_coords(ip+1) coord(4,2)=y_coords((iq+1)/2) coord(5,2)=y_coords((iq+1)/2)+ & (y_coords((iq+3)/2)-y_coords((iq+1)/2))/d3 coord(6,2)=y_coords((iq+1)/2)+ & two*(y_coords((iq+3)/2)-y_coords((iq+1)/2))/d3 coord(7,2)=y_coords((iq+3)/2) ELSE IF(dir=='x'.OR.dir=='r')THEN num(1)=(iq-2)/2*(3nxe+1)3+9nxe+3+3ip+1 num(2)=num(1)-1 num(3)=num(1)-2 num(4)=num(1)-3 num(5)=(iq-2)/2*(3nxe+1)3+6nxe+2+3ip-1 num(6)=(iq-2)/2*(3nxe+1)3+3nxe+1+3ip num(7)=(iq-2)/2*(3nxe+1)3+3ip+1 num(8)=num(6)+1 num(9)=num(5)+2 num(10)=num(9)-1 ELSE num(1)=(9(nye-2)/2+12)(ip-1)+9(nye-2)/2+12+3iq/2+1 num(2)=(9(nye-2)/2+12)(ip-1)+3(nye-2)+8+3iq/2+1 num(3)=(9(nye-2)/2+12)(ip-1)+3(nye-2)/2+4+3iq/2+1 num(4)=(9(nye-2)/2+12)(ip-1)+3iq/2+1 num(5)=num(3)-1 num(6)=num(2)-2 num(7)=num(1)-3 num(8)=num(1)-2 num(9)=num(1)-1 num(10)=num(2)-1 END IF ! coord(1,1)=x_coords(ip+1) coord(2,1)=x_coords(ip+1)-(x_coords(ip+1)-x_coords(ip))/d3 coord(3,1)=x_coords(ip+1)-two*(x_coords(ip+1)-x_coords(ip))/d3 coord(4,1)=x_coords(ip) coord(4,2)=y_coords((iq+2)/2) coord(5,2)=y_coords((iq+2)/2)-(y_coords((iq+2)/2)-y_coords(iq/2))/d3 coord(6,2)=y_coords((iq+2)/2)- & two*(y_coords((iq+2)/2)-y_coords(iq/2))/d3 coord(7,2) =y_coords(iq/2) END IF coord(5,1)=coord(3,1) coord(6,1)=coord(2,1) coord(7,1)=coord(1,1) coord(8,1)=coord(1,1) coord(9,1)=coord(1,1) coord(10,1)=coord(2,1) coord(1,2)=coord(4,2) coord(2,2)=coord(4,2) coord(3,2)=coord(4,2) coord(8,2)=coord(6,2) coord(9,2)=coord(5,2) coord(10,2)=coord(5,2) CASE(15) IF(MOD(iq,2)/=0)THEN IF(dir=='x'.OR.dir=='r')THEN fac1=4*(4nxe+1)(iq-1)/2 num(1)=fac1+4ip-3 num(2)=num(1)+1 num(3)=num(1)+2 num(4)=num(1)+3 num(5)=num(1)+4 num(6)=fac1+ 4nxe+1+4ip num(7)=fac1+ 8nxe+1+4ip num(8)=fac1+12nxe+1+4ip num(9)=fac1+16nxe+1+4ip num(10)=num(8)-1 num(11)=num(7)-2 num(12)=num(6)-3 num(13)=num(12)+1 num(14)=num(12)+2 num(15)=num(11)+1 ELSE fac1=4(2nye+1)(ip-1)+2iq-1 num(1)=fac1 num(2)=fac1+2nye+1 num(3)=fac1+4nye+2 num(4)=fac1+6nye+3 num(5)=fac1+8nye+4 num(6)=fac1+6nye+4 num(7)=fac1+4nye+4 num(8)=fac1+2nye+4 num(9)=fac1+4 num(10)=fac1+3 num(11)=fac1+2 num(12)=fac1+1 num(13)=fac1+2nye+2 num(14)=fac1+4nye+3 num(15)=fac1+2nye+3
END IF ! coord(1,1)=x_coords(ip) coord(1,2)=y_coords((iq+1)/2) coord(5,1)=x_coords(ip+1)
coord(5,2)=y_coords((iq+1)/2) coord(9,1)=x_coords(ip)
coord(9,2)=y_coords((iq+3)/2) ELSE IF(dir=='x'.OR.dir=='r')THEN fac1=4(4nxe+1)(iq-2)/2 num(1)=fac1+16nxe+5+4ip num(2)=num(1)-1 num(3)=num(1)-2 num(4)=num(1)-3 num(5)=num(1)-4 num(6)=fac1+12nxe+1+4ip num(7)=fac1+8nxe+1+4ip num(8)=fac1+4nxe+1+4ip num(9)=fac1+4ip+1 num(10)=num(8)+1 num(11)=num(7)+2 num(12)=num(6)+3 num(13)=num(12)-1 num(14)=num(12)-2 num(15)=num(11)-1 ELSE fac1=4(2nye+1)(ip-1)+2iq+8nye+5 num(1)=fac1 num(2)=fac1-2nye-1 num(3)=fac1-4nye-2 num(4)=fac1-6nye-3 num(5)=fac1-8nye-4 num(6)=fac1-6nye-4
num(7)=fac1-4nye-4 num(8)=fac1-2nye-4 num(9)=fac1-4 num(10)=fac1-3 num(11)=fac1-2 num(12)=fac1-1 num(13)=fac1-2nye-2
num(14)=fac1-4nye-3 num(15)=fac1-2nye-3 END IF ! coord(1,1)=x_coords(ip+1) coord(1,2)=y_coords((iq+2)/2) coord(5,1)=x_coords(ip)
coord(5,2)=y_coords((iq+2)/2) coord(9,1)=x_coords(ip+1) coord(9,2)=y_coords(iq/2) END IF coord(3,:)=pt5*(coord(1,:)+coord(5,:)) coord(7,:)=pt5*(coord(5,:)+coord(9,:)) coord(11,:)=pt5*(coord(9,:)+coord(1,:)) coord(2,:)=pt5*(coord(1,:)+coord(3,:)) coord(4,:)=pt5*(coord(3,:)+coord(5,:)) coord(6,:)=pt5*(coord(5,:)+coord(7,:)) coord(8,:)=pt5*(coord(7,:)+coord(9,:)) coord(10,:)=pt5*(coord(9,:)+coord(11,:)) coord(12,:)=pt5*(coord(11,:)+coord(1,:)) coord(15,:)=pt5*(coord(7,:)+coord(11,:)) coord(14,:)=pt5*(coord(3,:)+coord(7,:)) coord(13,:)=pt5*(coord(2,:)+coord(15,:)) CASE DEFAULT WRITE(11,'(a)')'Wrong number of nodes for triangular element' STOP END SELECT ELSE nye=UBOUND(y_coords,1)-1 IF(dir=='x'.OR.dir=='r')THEN iq=(iel-1)/nxe+1 ip=iel-(iq-1)nxe ELSE ip=(iel-1)/nye+1 iq=iel-(ip-1)nye END IF SELECT CASE(nod) CASE(4) IF(dir=='x'.OR.dir=='r')THEN num(1)=iq(nxe+1)+ip num(2)=(iq-1)(nxe+1)+ip num(3)=num(2)+1 num(4)=num(1)+1 ELSE num(1)=(ip-1)(nye+1)+iq+1 num(2)=num(1)-1 num(3)=ip(nye+1)+iq num(4)=num(3)+1 END IF ! coord(1:2,1)=x_coords(ip) coord(3:4,1)=x_coords(ip+1) coord(1,2)=y_coords(iq+1) coord(2:3,2)=y_coords(iq) coord(4,2)=coord(1,2) CASE(8) IF(dir=='x'.OR.dir=='r')THEN num(1)=iq*(3nxe+2)+2ip-1 num(2)=iq*(3nxe+2)+ip-nxe-1 num(3)=(iq-1)(3nxe+2)+2ip-1 num(4)=num(3)+1 num(5)=num(4)+1 num(6)=num(2)+1 num(7)=num(1)+2 num(8)=num(1)+1 ELSE num(1)=(ip-1)(3nye+2)+2iq+1 num(2)=num(1)-1 num(3)=num(1)-2 num(4)=(ip-1)(3nye+2)+2nye+iq+1 num(5)=ip*(3nye+2)+2iq-1 num(6)=num(5)+1 num(7)=num(5)+2 num(8)=num(4)+1 END IF ! coord(1:3,1)=x_coords(ip) coord(5:7,1)=x_coords(ip+1) coord(4,1)=pt5*(coord(3,1)+coord(5,1)) coord(8,1)=pt5*(coord(7,1)+coord(1,1)) coord(1,2)=y_coords(iq+1) coord(7:8,2)=y_coords(iq+1) coord(3:5,2)=y_coords(iq) coord(2,2)=pt5*(coord(1,2)+coord(3,2)) coord(6,2)=pt5*(coord(5,2)+coord(7,2)) CASE(9) IF(dir=='x'.OR.dir=='r')THEN num(1)=iq*(4nxe+2)+2ip-1 num(2)=iq*(4nxe+2)+2ip-nxe-4 num(3)= (iq-1)(4nxe+2)+2ip-1 num(4)=num(3)+1 num(5)=num(4)+1 num(6)=num(2)+2 num(7)=num(1)+2 num(8)=num(1)+1 num(9)=num(2)+1 ELSE num(1)=(ip-1)2(2nye+1)+2iq+1 num(2)=num(1)-1 num(3)=num(1)-2 num(4)=(ip-1)2(2nye+1)+2nye+2iq num(5)=ip2(2nye+1)+2iq-1 num(6)=num(5)+1 num(7)=num(5)+2 num(8)=num(4)+2 num(9)=num(4)+1 END IF ! coord(1:3,1)=x_coords(ip) coord(5:7,1)=x_coords(ip+1) coord(4,1)=pt5*(coord(3,1)+coord(5,1)) coord(8,1)=pt5*(coord(7,1)+coord(1,1)) coord(1,2)=y_coords(iq+1) coord(7:8,2)=y_coords(iq+1) coord(3:5,2)=y_coords(iq) coord(2,2)=pt5*(coord(1,2)+coord(3,2)) coord(6,2)=pt5*(coord(5,2)+coord(7,2)) coord(9,:)=pt5*(coord(4,:)+coord(8,:)) CASE DEFAULT WRITE(11,'(a)')'Wrong number of nodes for quadrilateral element' STOP END SELECT END IF RETURN END SUBROUTINE geom_rect

SUBROUTINE fkdiag(kdiag,g) ! ! This subroutine computes the skyline profile. ! IMPLICIT NONE INTEGER,INTENT(IN)::g(:) INTEGER,INTENT(OUT)kdiag(:) INTEGERidof,i,iwp1,j,im,k idof=SIZE(g) DO i=1,idof iwp1=1 IF(g(i)/=0)THEN DO j=1,idof IF(g(j)/=0)THEN im=g(i)-g(j)+1 IF(im>iwp1)iwp1=im END IF END DO k=g(i) IF(iwp1>kdiag(k))kdiag(k)=iwp1 END IF END DO RETURN END SUBROUTINE fkdiag

SUBROUTINE mesh(g_coord,g_num,ips) ! ! This subroutine produces a PostScript output file '.msh' displaying ! the undeformed finite element mesh. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN)::g_coord(:,:) INTEGER,INTENT(IN)::g_num(:,:),ips REAL(iwp)xmin,xmax,ymin,ymax,width,height,scale=72,sxy,xo,yo,x,y, & pt5=0.5_iwp,opt5=1.5_iwp,fpt5=5.5_iwp,d8=8.0_iwp,ept5=8.5_iwp, & d11=11.0_iwp INTEGERi,ii,j,jj,nn,nod,nel OPEN(ips,FILE='fe95.msh') ! ! compute size of mesh ! nn=UBOUND(g_coord,2) xmin=g_coord(1,1) xmax=g_coord(1,1) ymin=g_coord(2,1) ymax=g_coord(2,1) DO i=2,nn IF(g_coord(1,i)<xmin)xmin=g_coord(1,i)
IF(g_coord(1,i)>xmax)xmax=g_coord(1,i)
IF(g_coord(2,i)<ymin)ymin=g_coord(2,i)
IF(g_coord(2,i)>ymax)ymax=g_coord(2,i)
END DO width =xmax-xmin height=ymax-ymin ! ! allow 1.5' margin minimum on each side of figure ! IF(height.GE.d11/ept5width)THEN ! ! height governs the scale ! sxy=scaled8/height xo=scalept5*(ept5-d8width/height) yo=scaleopt5 ELSE ! ! width governs the scale ! sxy=scalefpt5/width xo=scaleopt5 yo=scalept5(d11-fpt5height/width) END IF ! ! start PostScript output ! WRITE(ips,'(a)')'%!PS-Adobe-1.0' WRITE(ips,'(a)')'%%DocumentFonts: none' WRITE(ips,'(a)')'%%Pages: 1' WRITE(ips,'(a)')'%%EndComments' WRITE(ips,'(a)')'/m def' WRITE(ips,'(a)')'/l def' WRITE(ips,'(a)')'/s def' WRITE(ips,'(a)')'/c def' WRITE(ips,'(a)')'%%EndProlog' WRITE(ips,'(a)')'%%Page: 0 1' WRITE(ips,'(a)')'gsave' WRITE(ips,'(2f9.2,a)') xo, yo, ' translate' WRITE(ips,'(f9.2,a)') 0.5, ' setlinewidth' ! ! draw the mesh ! nod=UBOUND(g_num,1) nel=UBOUND(g_num,2) IF(nod9)nod=8 IF(nod10)nod=9 IF(nod15)nod=12 DO i=1,nel ii=g_num(1,i) IF(ii0)CYCLE x=sxy(g_coord(1,ii)-xmin) y=sxy*(g_coord(2,ii)-ymin) WRITE(ips,'(2f9.2,a)')x,y,' m' DO j=2,nod jj=g_num(j,i) x=sxy*(g_coord(1,jj)-xmin) y=sxy*(g_coord(2,jj)-ymin) WRITE(ips,'(2f9.2,a)') x, y,' l' END DO WRITE(ips,'(a)')'c s' END DO ! ! close output file ! WRITE(ips,'(a)')'grestore' WRITE(ips,'(a)')'showpage' CLOSE(ips) ! RETURN END SUBROUTINE mesh

SUBROUTINE sample(element,s,wt) ! ! This subroutine returns the local coordinates and weighting coefficients ! of the integrating points. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(OUT)s(:,:) REAL(iwp),INTENT(OUT),OPTIONALwt(:) CHARACTER(),INTENT(IN)element INTEGERnip REAL(iwp)::root3,r15,w(3),v(9),b,c root3=1.0_iwp/SQRT(3.0_iwp) r15=0.2_iwpSQRT(15.0_iwp) nip=UBOUND(s,1) w=(/5.0_iwp/9.0_iwp,8.0_iwp/9.0_iwp,5.0_iwp/9.0_iwp/) v=(/5.0_iwp/9.0_iwpw,8.0_iwp/9.0_iwpw,5.0_iwp/9.0_iwpw/) SELECT CASE(element) CASE('line') SELECT CASE(nip) CASE(1) s(1,1)=0.0_iwp wt(1) =2.0_iwp CASE(2) s(1,1)=-0.577350269189626_iwp s(2,1)= 0.577350269189626_iwp wt(1) = 1.000000000000000_iwp wt(2) = 1.000000000000000_iwp CASE(3) s(1,1)=-0.774596669241484_iwp s(2,1)= 0.000000000000000_iwp s(3,1)= 0.774596669241484_iwp wt(1) = 0.555555555555556_iwp wt(2) = 0.888888888888889_iwp wt(3) = 0.555555555555556_iwp CASE(4) s(1,1)=-0.861136311594053_iwp s(2,1)=-0.339981043584856_iwp s(3,1)= 0.339981043584856_iwp s(4,1)= 0.861136311594053_iwp wt(1) = 0.347854845137454_iwp wt(2) = 0.652145154862546_iwp wt(3) = 0.652145154862546_iwp wt(4) = 0.347854845137454_iwp CASE(5) s(1,1)=-0.906179845938664_iwp s(2,1)=-0.538469310105683_iwp s(3,1)= 0.000000000000000_iwp s(4,1)= 0.538469310105683_iwp s(5,1)= 0.906179845938664_iwp wt(1) = 0.236926885056189_iwp wt(2) = 0.478628670499366_iwp wt(3) = 0.568888888888889_iwp wt(4) = 0.478628670499366_iwp wt(5) = 0.236926885056189_iwp CASE(6) s(1,1)=-0.932469514203152_iwp s(2,1)=-0.661209386466265_iwp s(3,1)=-0.238619186083197_iwp s(4,1)= 0.238619186083197_iwp s(5,1)= 0.661209386466265_iwp s(6,1)= 0.932469514203152_iwp wt(1) = 0.171324492379170_iwp wt(2) = 0.360761573048139_iwp wt(3) = 0.467913934572691_iwp wt(4) = 0.467913934572691_iwp wt(5) = 0.360761573048139_iwp wt(6) = 0.171324492379170_iwp CASE DEFAULT WRITE(,)'Wrong number of integrating points for a line' END SELECT CASE('triangle') SELECT CASE(nip) CASE(1) s(1,1)= 0.333333333333333_iwp s(1,2)= 0.333333333333333_iwp wt(1) = 0.500000000000000_iwp CASE(3) s(1,1)= 0.500000000000000_iwp s(1,2)= 0.500000000000000_iwp s(2,1)= 0.500000000000000_iwp s(2,2)= 0.000000000000000_iwp s(3,1)= 0.000000000000000_iwp s(3,2)= 0.500000000000000_iwp wt(1:3)=0.333333333333333_iwp wt=0.5_iwpwt CASE(4) s(1,1)= 0.6_iwp s(1,2)= 0.2_iwp s(2,1)= 0.2_iwp s(2,2)= 0.6_iwp s(3,1)= 0.2_iwp s(3,2)= 0.2_iwp s(4,1)= 0.333333333333333_iwp s(4,2)= 0.333333333333333_iwp wt(1:3)= -0.520833333333333_iwp wt(4)=-0.5625_iwp wt=0.5_iwpwt CASE(6) s(1,1)= 0.816847572980459_iwp s(1,2)= 0.091576213509771_iwp s(2,1)= 0.091576213509771_iwp s(2,2)= 0.816847572980459_iwp s(3,1)= 0.091576213509771_iwp s(3,2)= 0.091576213509771_iwp s(4,1)= 0.108103018168070_iwp s(4,2)= 0.445948490915965_iwp s(5,1)= 0.445948490915965_iwp s(5,2)= 0.108103018168070_iwp s(6,1)= 0.445948490915965_iwp s(6,2)= 0.445948490915965_iwp wt(1:3)=0.109951743655322_iwp wt(4:6)=0.223381589678011_iwp wt=0.5_iwpwt CASE(7) s(1,1)= 0.333333333333333_iwp s(1,2)= 0.333333333333333_iwp s(2,1)= 0.797426985353087_iwp s(2,2)= 0.101286507323456_iwp s(3,1)= 0.101286507323456_iwp s(3,2)= 0.797426985353087_iwp s(4,1)= 0.101286507323456_iwp s(4,2)= 0.101286507323456_iwp s(5,1)= 0.470142064105115_iwp s(5,2)= 0.059715871789770_iwp s(6,1)= 0.059715871789770_iwp s(6,2)= 0.470142064105115_iwp s(7,1)= 0.470142064105115_iwp s(7,2)= 0.470142064105115_iwp wt(1) = 0.225000000000000_iwp wt(2:4)=0.125939180544827_iwp wt(5:7)=0.132394152788506_iwp wt=0.5_iwpwt CASE(12) s(1,1)= 0.873821971016996_iwp s(1,2)= 0.063089014491502_iwp s(2,1)= 0.063089014491502_iwp s(2,2)= 0.873821971016996_iwp s(3,1)= 0.063089014491502_iwp s(3,2)= 0.063089014491502_iwp s(4,1)= 0.501426509658179_iwp s(4,2)= 0.249286745170910_iwp s(5,1)= 0.249286745170910_iwp s(5,2)= 0.501426509658179_iwp s(6,1)= 0.249286745170910_iwp s(6,2)= 0.249286745170910_iwp s(7,1) =0.053145049844817_iwp s(7,2) =0.310352451033784_iwp s(8,1) =0.310352451033784_iwp s(8,2) =0.053145049844817_iwp s(9,1) =0.053145049844817_iwp s(9,2) =0.636502499121398_iwp s(10,1)=0.310352451033784_iwp s(10,2)=0.636502499121398_iwp s(11,1)=0.636502499121398_iwp s(11,2)=0.053145049844817_iwp s(12,1)=0.636502499121398_iwp s(12,2)=0.310352451033784_iwp wt(1:3)=0.050844906370207_iwp wt(4:6)=0.116786275726379_iwp wt(7:12)=0.082851075618374_iwp wt=0.5_iwpwt CASE(16) s(1,1)=0.333333333333333_iwp s(1,2)=0.333333333333333_iwp s(2,1)=0.658861384496478_iwp s(2,2)=0.170569307751761_iwp s(3,1)=0.170569307751761_iwp s(3,2)=0.658861384496478_iwp s(4,1)=0.170569307751761_iwp s(4,2)=0.170569307751761_iwp s(5,1)=0.898905543365938_iwp s(5,2)=0.050547228317031_iwp s(6,1)=0.050547228317031_iwp s(6,2)=0.898905543365938_iwp s(7,1)=0.050547228317031_iwp s(7,2)=0.050547228317031_iwp s(8,1)=0.081414823414554_iwp s(8,2)=0.459292588292723_iwp s(9,1)=0.459292588292723_iwp s(9,2)=0.081414823414554_iwp s(10,1)=0.459292588292723_iwp s(10,2)=0.459292588292723_iwp s(11,1)=0.008394777409958_iwp s(11,2)=0.263112829634638_iwp s(12,1)=0.008394777409958_iwp s(12,2)=0.728492392955404_iwp s(13,1)=0.263112829634638_iwp s(13,2)=0.008394777409958_iwp s(14,1)=0.263112829634638_iwp s(14,2)=0.728492392955404_iwp s(15,1)=0.728492392955404_iwp s(15,2)=0.008394777409958_iwp s(16,1)=0.728492392955404_iwp s(16,2)=0.263112829634638_iwp wt(1)=0.144315607677787_iwp wt(2:4)=0.103217370534718_iwp wt(5:7)=0.032458497623198_iwp wt(8:10)=0.095091634267284_iwp wt(11:16)=0.027230314174435_iwp wt=0.5_iwpwt CASE DEFAULT WRITE(,)'wrong number of integrating points for a triangle' END SELECT CASE('quadrilateral') SELECT CASE(nip) CASE(1) s(1,1)=0.0_iwp s(1,2)=0.0_iwp wt(1)=4.0_iwp CASE(4) s(1,1)=-root3 s(1,2)= root3 s(2,1)= root3 s(2,2)= root3 s(3,1)=-root3 s(3,2)=-root3 s(4,1)= root3 s(4,2)=-root3 wt=1.0_iwp CASE(9) s(1:7:3,1)=-r15 s(2:8:3,1)=0.0_iwp s(3:9:3,1)=r15 s(1:3,2) =r15 s(4:6,2) =0.0_iwp s(7:9,2) =-r15 wt= v CASE(16) s(1:13:4,1)=-0.861136311594053_iwp s(2:14:4,1)=-0.339981043584856_iwp s(3:15:4,1)= 0.339981043584856_iwp s(4:16:4,1)= 0.861136311594053_iwp s(1:4,2) = 0.861136311594053_iwp s(5:8,2) = 0.339981043584856_iwp s(9:12,2) =-0.339981043584856_iwp s(13:16,2) =-0.861136311594053_iwp wt(1) = 0.121002993285602_iwp wt(4) = wt(1) wt(13) = wt(1) wt(16) = wt(1) wt(2) = 0.226851851851852_iwp wt(3) = wt(2) wt(5) = wt(2) wt(8) = wt(2) wt(9) = wt(2) wt(12) = wt(2) wt(14) = wt(2) wt(15) = wt(2) wt(6) = 0.425293303010694_iwp wt(7) = wt(6) wt(10) = wt(6) wt(11) = wt(6) CASE(25) s(1:21:5,1)= 0.906179845938664_iwp s(2:22:5,1)= 0.538469310105683_iwp s(3:23:5,1)= 0.0_iwp s(4:24:5,1)=-0.538469310105683_iwp s(5:25:5,1)=-0.906179845938664_iwp s( 1: 5,2) = 0.906179845938664_iwp s( 6:10,2) = 0.538469310105683_iwp s(11:15,2) = 0.0_iwp s(16:20,2) =-0.538469310105683_iwp s(21:25,2) =-0.906179845938664_iwp wt(1) =0.056134348862429_iwp wt(2) =0.113400000000000_iwp wt(3) =0.134785072387521_iwp wt(4) =0.113400000000000_iwp wt(5) =0.056134348862429_iwp wt(6) =0.113400000000000_iwp wt(7) =0.229085404223991_iwp wt(8) =0.272286532550750_iwp wt(9) =0.229085404223991_iwp wt(10)=0.113400000000000_iwp wt(11)=0.134785072387521_iwp wt(12)=0.272286532550750_iwp wt(13)=0.323634567901235_iwp wt(14)=0.272286532550750_iwp wt(15)=0.134785072387521_iwp wt(16)=0.113400000000000_iwp wt(17)=0.229085404223991_iwp wt(18)=0.272286532550750_iwp wt(19)=0.229085404223991_iwp wt(20)=0.113400000000000_iwp wt(21)=0.056134348862429_iwp wt(22)=0.113400000000000_iwp wt(23)=0.134785072387521_iwp wt(24)=0.113400000000000_iwp wt(25)=0.056134348862429_iwp CASE DEFAULT WRITE(,)'wrong number of integrating points for a quadrilateral' END SELECT CASE('tetrahedron') ! for tetrahedra weights multiplied by 1/6 SELECT CASE(nip) CASE(1) s(1,1)=0.25_iwp s(1,2)=0.25_iwp s(1,3)=0.25_iwp wt(1)=1.0_iwp/6.0_iwp CASE(4) s(1,1)=0.58541020_iwp s(1,2)=0.13819660_iwp s(1,3)=s(1,2) s(2,2)=s(1,1) s(2,3)=s(1,2) s(2,1)=s(1,2) s(3,3)=s(1,1) s(3,1)=s(1,2) s(3,2)=s(1,2) s(4,1)=s(1,2) s(4,2)=s(1,2) s(4,3)=s(1,2) wt(1:4)=0.25_iwp/6.0_iwp CASE(5) s(1,1)=0.25_iwp s(1,2)=0.25_iwp s(1,3)=0.25_iwp s(2,1)=0.5_iwp s(2,2)=1.0_iwp/6.0_iwp s(2,3)=s(2,2) s(3,2)=0.5_iwp s(3,3)=1.0_iwp/6.0_iwp s(3,1)=s(3,3) s(4,3)=0.5_iwp s(4,1)=1.0_iwp/6.0_iwp s(4,2)=s(4,1) s(5,1)=1.0_iwp/6.0_iwp s(5,2)=s(5,1) s(5,3)=s(5,1) wt(1)=-0.8_iwp wt(2)=9.0_iwp/20.0_iwp wt(3:5)=wt(2) wt=wt/6.0_iwp CASE DEFAULT WRITE(,)'wrong number of integrating points for a tetrahedron' END SELECT CASE('hexahedron') SELECT CASE(nip) CASE(1) s(1,1:3)=0.0_iwp wt(1)=8.0_iwp CASE(8) s(1,1)= root3 s(1,2)= root3 s(1,3)= root3 s(2,1)= root3 s(2,2)= root3 s(2,3)=-root3 s(3,1)= root3 s(3,2)=-root3 s(3,3)= root3 s(4,1)= root3 s(4,2)=-root3 s(4,3)=-root3 s(5,1)=-root3 s(5,2)= root3 s(5,3)= root3 s(6,1)=-root3 s(6,2)=-root3 s(6,3)= root3 s(7,1)=-root3 s(7,2)= root3 s(7,3)=-root3 s(8,1)=-root3 s(8,2)=-root3 s(8,3)=-root3 wt=1.0_iwp CASE(14) b=0.795822426_iwp c=0.758786911_iwp wt(1:6)=0.886426593_iwp wt(7:14)=0.335180055_iwp s(1,1)=-b s(2,1)=b s(3,2)=-b s(4,2)=b s(5,3)=-b s(6,3)=b s(7:,:)=c s(7,1)=-c s(7,2)=-c s(7,3)=-c s(8,2)=-c s(8,3)=-c s(9,1)=-c s(9,3)=-c s(10,3)=-c s(11,1)=-c s(11,2)=-c s(12,2)=-c s(13,1)=-c CASE(15) b=1.0_iwp c =0.674199862_iwp wt(1) =1.564444444_iwp wt(2:7) =0.355555556_iwp wt(8:15)=0.537777778_iwp s(2,1)=-b s(3,1)=b s(4,2)=-b s(5,2)=b s(6,3)=-b s(7,3)=b s(8:,:)=c s(8,1)=-c s(8,2)=-c s(8,3)=-c s(9,2)=-c s(9,3)=-c s(10,1)=-c s(10,3)=-c s(11,3)=-c s(12,1)=-c s(12,2)=-c s(13,2)=-c s(14,1)=-c CASE(27) wt=(/5.0_iwp/9.0_iwpv,8.0_iwp/9.0_iwpv,5.0_iwp/9.0_iwpv/) s(1:7:3,1)=-r15 s(2:8:3,1)=0.0_iwp s(3:9:3,1)=r15 s(1:3,3)=r15 s(4:6,3)=0.0_iwp s(7:9,3)=-r15 s(1:9,2)=-r15 s(10:16:3,1)=-r15 s(11:17:3,1)=0.0_iwp s(12:18:3,1)=r15 s(10:12,3)=r15 s(13:15,3)=0.0_iwp s(16:18,3)=-r15 s(10:18,2)=0.0_iwp s(19:25:3,1)=-r15 s(20:26:3,1)=0.0_iwp s(21:27:3,1)=r15 s(19:21,3)= r15 s(22:24,3)=0.0_iwp s(25:27,3)=-r15 s(19:27,2)= r15 CASE DEFAULT WRITE(,)'wrong number of integrating points for a hexahedron' END SELECT CASE DEFAULT WRITE(,)'not a valid element type' END SELECT RETURN END SUBROUTINE sample

SUBROUTINE shape_der(der,points,i) ! ! This subroutine produces derivatives of shape functions withe respect ! to local coordinates. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) INTEGER,INTENT(IN)::i REAL(iwp),INTENT(IN)::points(:,:) REAL(iwp),INTENT(OUT)::der(:,:) REAL(iwp)eta,xi,zeta,xi0,eta0,zeta0,etam,etap,xim,xip,c1,c2,c3 REAL(iwp)t1,t2,t3,t4,t5,t6,t7,t8,t9,x2p1,x2m1,e2p1,e2m1,zetam,zetap REAL,PARAMETERzero=0.0_iwp,pt125=0.125_iwp,pt25=0.25_iwp,pt5=0.5_iwp, & pt75=0.75_iwp,one=1.0_iwp,two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp,& d6=6.0_iwp,d8=8.0_iwp,d9=9.0_iwp,d10=10.0_iwp,d11=11.0_iwp, & d12=12.0_iwp,d16=16.0_iwp,d18=18.0_iwp,d27=27.0_iwp,d32=32.0_iwp, & d36=36.0_iwp,d54=54.0_iwp,d64=64.0_iwp,d128=128.0_iwp INTEGERxii(20),etai(20),zetai(20),l,ndim,nod ndim=UBOUND(der,1) nod= UBOUND(der,2) SELECT CASE(ndim) CASE(1) ! one dimensional elements xi=points(i,1) SELECT CASE(nod) CASE(2) der(1,1)=-pt5 der(1,2)= pt5 CASE(3) t1=-one-xi t2=-xi
t3=one-xi der(1,1)=-(t3+t2)/two
der(1,2)=(t3+t1)
der(1,3)=-(t2+t1)/two
CASE(4) t1=-one-xi t2=-one/d3-xi t3=one/d3-xi t4=one-xi der(1,1)=-(t3t4+t2t4+t2t3)d9/d16
der(1,2)=(t3t4+t1t4+t1t3)d27/d16 der(1,3)=-(t2t4+t1t4+t1t2)d27/d16 der(1,4)=(t2t3+t1t3+t1t2)d9/d16
CASE(5) t1=-one-xi t2=-pt5-xi t3=-xi t4=pt5-xi t5=one-xi der(1,1)=-(t3t4t5+t2t4t5+t2t3t5+t2t3t4)two/d3
der(1,2)=(t3t4t5+t1t4t5+t1t3t5+t1t3t4)d8/d3 der(1,3)=-(t2t4t5+t1t4t5+t1t2t5+t1t2t4)d4 der(1,4)=(t2t3t5+t1t3t5+t1t2t5+t1t2t3)d8/d3 der(1,5)=-(t2t3t4+t1t3t4+t1t2t4+t1t2t3)two/d3 CASE DEFAULT WRITE(,)'wrong number of nodes in shape_der'
END SELECT CASE(2) ! two dimensional elements xi=points(i,1) eta=points(i,2) c1=xi c2=eta c3=one-c1-c2 etam=pt25(one-eta) etap=pt25*(one+eta) xim= pt25*(one-xi) xip= pt25*(one+xi) x2p1=twoxi+one x2m1=twoxi-one e2p1=twoeta+one e2m1=twoeta-one SELECT CASE(nod) CASE(3) der(1,1)=one der(1,3)=zero der(1,2)=-one der(2,1)=zero der(2,3)=one der(2,2)=-one CASE(6) der(1,1)=d4c1-one der(1,6)=d4c2 der(1,5)=zero
der(1,4)=-d4c2 der(1,3)=-(d4c3-one) der(1,2)=d4*(c3-c1) der(2,1)=zero der(2,6)=d4c1 der(2,5)=d4c2-one der(2,4)=d4*(c3-c2) der(2,3)=-(d4c3-one)
der(2,2)=-d4c1 CASE(10)
der(1,1)=(d27c1**2-d18c1+two)/two der(1,9)=(d9*(d6c1-one)c2)/two der(1,8)=(d9(d3c2-one)c2)/two der(1,7)=zero der(1,6)=-(d9(d3c2-one)c2)/two der(1,5)= (d9(d6c1+d6c2-d5)c2)/two der(1,4)=-(d27c1**2+d54c1c2-d36c1+d27c2**2-d36c2+d11)/two der(1,3)= (d9*(d9c1**2+d12c1c2-d10c1+d3c2**2-d5c2+two))/two der(1,2)=-(d9*(d9c1**2+d6c1c2-d8c1-c2+one))/two der(1,10)=-d27*(((c2-one)+c1)+c1)c2 der(2,1)=zero der(2,9)= (d9(d3c1-one)c1)/two der(2,8)= (d9(d6c2-one)c1)/two der(2,7)=(d27c22-d18c2+two)/two der(2,6)=-(d9((c1+c2-one)(d6c2-one)+(d3c2-one)c2))/two der(2,5)= (d9(d3c12+d12c1c2-d5c1+d9c22-d10c2+two))/two der(2,4)=-(d27c12+d54c1c2-d36c1+d27c2**2-d36c2+d11)/two der(2,3)= (d9(d6c1+d6c2-d5)c1)/two der(2,2)=-(d9(d3c1-one)c1)/two der(2,10)=-d27(((c2-one)+c1)+c2)c1 CASE(15)
t1=c1-pt25
t2=c1-pt5 t3=c1-pt75
t4=c2-pt25 t5=c2-pt5
t6=c2-pt75 t7=c3-pt25
t8=c3-pt5 t9=c3-pt75 der(1,1)=d32/d3(t2t3*(t1+c1)+c1t1(t3+t2)) der(1,12)=d128/d3c2(t2*(t1+c1)+c1t1) der(1,11)=d64c2t4(t1+c1) der(1,10)=d128/d3c2t4t5
der(1,9)=zero der(1,8)=-d128/d3c2t4t5 der(1,7)=-d64c2t4*(t7+c3) der(1,6)=-d128/d3c2(t8*(t7+c3)+c3t7) der(1,5)=-d32/d3(t8t9(t7+c3)+c3t7(t8+t9)) der(1,4)=d128/d3*(c3t7t8-c1*(t8*(t7+c3)+c3t7)) der(1,3)=d64(c3t7(t1+c1)-c1t1(t7+c3)) der(1,2)=d128/d3*(c3*(t2*(t1+c1)+c1t1)-c1t1t2) der(1,13)=d128c2*(c3*(t1+c1)-c1t1) der(1,15)=d128c2t4(c3-c1) der(1,14)=d128c2(c3t7-c1(t7+c3)) der(2,1)=zero der(2,12)=d128/d3c1t1t2 der(2,11)=d64c1t1(t4+c2) der(2,10)=d128/d3c1(t5*(t4+c2)+c2t4) der(2,9)=d32/d3(t5t6(t4+c2)+c2t4(t6+t5)) der(2,8)=d128/d3*((c3*(t5*(t4+c2)+c2t4))-c2t4t5) der(2,7)=d64(c3t7(t4+c2)-c2t4(t7+c3)) der(2,6)=d128/d3*(c3t7t8-c2*(t8*(t7+c3)+c3t7)) der(2,5)=-d32/d3(t8t9(t7+c3)+c3t7(t8+t9)) der(2,4)=-d128/d3c1(t8*(t7+c3)+c3t7) der(2,3)=-d64c1t1(t7+c3)
der(2,2)=-d128/d3c1t1t2 der(2,13)=d128c1t1(c3-c2) der(2,15)=d128c1(c3*(t4+c2)-c2t4) der(2,14)=d128c1*(c3t7-c2(c3+t7))
CASE (4)
der(1,1)=-etam der(1,2)=-etap der(1,3)=etap der(1,4)=etam der(2,1)=-xim der(2,2)=xim der(2,3)=xip der(2,4)=-xip CASE(8) der(1,1)=etam*(twoxi+eta) der(1,2)=-d8etametap der(1,3)=etap(twoxi-eta) der(1,4)=-d4etapxi der(1,5)=etap(twoxi+eta) der(1,6)=d8etapetam der(1,7)=etam(twoxi-eta) der(1,8)=-d4etamxi der(2,1)=xim(xi+twoeta) der(2,2)=-d4ximeta der(2,3)=xim(twoeta-xi) der(2,4)=d8ximxip der(2,5)=xip(xi+twoeta) der(2,6)=-d4xipeta der(2,7)=xip(twoeta-xi) der(2,8)=-d8ximxip
CASE(9) etam=eta-one etap=eta+one xim=xi-one xip=xi+one der(1,1)=pt25x2m1etaetam
der(1,2)=-pt5x2m1etapetam der(1,3)=pt25x2m1etaetap
der(1,4)=-xietaetap der(1,5)=pt25x2p1etaetap
der(1,6)=-pt5x2p1etapetam der(1,7)=pt25x2p1etaetam
der(1,8)=-xietaetam der(1,9)=twoxietapetam
der(2,1)=pt25xixime2m1 der(2,2)=-xiximeta
der(2,3)=pt25xixime2p1 der(2,4)=-pt5xipxime2p1
der(2,5)=pt25xixipe2p1 der(2,6)=-xixipeta
der(2,7)=pt25xixipe2m1 der(2,8)=-pt5xipxime2m1
der(2,9)=twoxipximeta CASE DEFAULT WRITE(,)'wrong number of nodes in shape_der'
END SELECT CASE(3) ! d3 dimensional elements xi=points(i,1) eta=points(i,2) zeta=points(i,3) etam=one-eta xim=one-xi zetam=one-zeta etap=eta+one xip=xi+one zetap=zeta+one SELECT CASE(nod) CASE(4) der(1:3,1:4)=zero der(1,1)=one der(2,2)=one
der(3,3)=one der(1,4)=-one der(2,4)=-one der(3,4)=-one
CASE(8) der(1,1)=-pt125etamzetam
der(1,2)=-pt125etamzetap der(1,3)= pt125etamzetap
der(1,4)= pt125etamzetam der(1,5)=-pt125etapzetam
der(1,6)=-pt125etapzetap der(1,7)= pt125etapzetap
der(1,8)= pt125etapzetam der(2,1)=-pt125ximzetam
der(2,2)=-pt125ximzetap der(2,3)=-pt125xipzetap
der(2,4)=-pt125xipzetam der(2,5)= pt125ximzetam
der(2,6)= pt125ximzetap der(2,7)= pt125xipzetap
der(2,8)= pt125xipzetam der(3,1)=-pt125ximetam
der(3,2)= pt125ximetam der(3,3)= pt125xipetam
der(3,4)=-pt125xipetam der(3,5)=-pt125ximetap
der(3,6)= pt125ximetap der(3,7)= pt125xipetap
der(3,8)=-pt125xipetap
CASE(14) ! type 6 element der(1,1)= (twoxieta+twoxizeta+d4xi+etazeta+eta+zeta) & (eta-one)(zeta-one)/d8 der(1,2)=-(twoxieta-twoxizeta+d4xi-etazeta+eta-zeta) & (eta-one)(zeta+one)/d8 der(1,3)=-(twoxieta-twoxizeta+d4xi+etazeta-eta+zeta) & (eta-one)(zeta+one)/d8 der(1,4)= (twoxieta+twoxizeta+d4xi-etazeta-eta-zeta) & (eta-one)(zeta-one)/d8 der(1,5)= -(eta-one)(zeta+one)(zeta-one)xi der(1,6)=-(eta+one)(eta-one)(zeta+one)(zeta-one)/two der(1,7)= (eta+one)(eta-one)(zeta+one)xi der(1,8)= (eta+one)(eta-one)(zeta+one)(zeta-one)/two der(1,9)= -(eta+one)(eta-one)(zeta-one)xi
der(1,10)= (twoxieta-twoxizeta-d4xi+etazeta+eta-zeta)* & (eta+one)(zeta-one)/d8 der(1,11)=-(twoxieta+twoxizeta-d4xi-etazeta+eta+zeta) & (eta+one)(zeta+one)/d8 der(1,12)=-(twoxieta+twoxizeta-d4xi+etazeta-eta-zeta) & (eta+one)(zeta+one)/d8 der(1,13)= (twoxieta-twoxizeta-d4xi-etazeta-eta+zeta) & (eta+one)(zeta-one)/d8 der(1,14)= (eta+one)(zeta+one)(zeta-one)xi der(2,1)= (twoxieta+xizeta+xi+twoetazeta+d4eta+zeta)* & (xi-one)(zeta-one)/d8
der(2,2)=-(twoxieta-xizeta+xi-twoetazeta+d4eta-zeta) & (xi-one)(zeta+one)/d8 der(2,3)=-(twoxieta-xizeta+xi+twoetazeta-d4eta+zeta) & (xi+one)(zeta+one)/d8 der(2,4)= (twoxieta+xizeta+xi-twoetazeta-d4eta-zeta) & (xi+one)(zeta-one)/d8 der(2,5)=-(xi+one)(xi-one)(zeta+one)(zeta-one)/two der(2,6)= -(xi-one)(zeta+one)(zeta-one)eta der(2,7)= (xi+one)(xi-one)(zeta+one)eta der(2,8)= (xi+one)(zeta+one)(zeta-one)eta der(2,9)= -(xi+one)(xi-one)(zeta-one)eta der(2,10)= (twoxieta-xizeta-xi+twoetazeta+d4eta-zeta)* & (xi-one)(zeta-one)/d8 der(2,11)=-(twoxieta+xizeta-xi-twoetazeta+d4eta+zeta) & (xi-one)(zeta+one)/d8 der(2,12)=-(twoxieta+xizeta-xi+twoetazeta-d4eta-zeta) & (xi+one)(zeta+one)/d8 der(2,13)= (twoxieta-xizeta-xi-twoetazeta-d4eta+zeta) & (xi+one)(zeta-one)/d8 der(2,14)= (xi+one)(xi-one)(zeta+one)(zeta-one)/two der(3,1)= (xieta+twoxizeta+xi+twoetazeta+eta+d4zeta)* & (xi-one)(eta-one)/d8 der(3,2)=-(xieta-twoxizeta+xi-twoetazeta+eta-d4zeta) & (xi-one)(eta-one)/d8 der(3,3)=-(xieta-twoxizeta+xi+twoetazeta-eta+d4zeta) & (xi+one)(eta-one)/d8 der(3,4)= (xieta+twoxizeta+xi-twoetazeta-eta-d4zeta) & (xi+one)(eta-one)/d8 der(3,5)= -(xi+one)(xi-one)(eta-one)zeta
der(3,6)= -(xi-one)(eta+one)(eta-one)zeta
der(3,7)= (xi+one)(xi-one)(eta+one)(eta-one)/two der(3,8)= (xi+one)(eta+one)(eta-one)zeta der(3,9)=-(xi+one)(xi-one)(eta+one)(eta-one)/two der(3,10)= (xieta-twoxizeta-xi+twoetazeta+eta-d4zeta)* & (xi-one)(eta+one)/d8 der(3,11)=-(xieta+twoxizeta-xi-twoetazeta+eta+d4zeta) & (xi-one)(eta+one)/d8 der(3,12)=-(xieta+twoxizeta-xi+twoetazeta-eta-d4zeta) & (xi+one)(eta+one)/d8 der(3,13)= (xieta-twoxizeta-xi-twoetazeta-eta+d4zeta) & (xi+one)(eta+one)/d8 der(3,14)= (xi+one)(xi-one)(eta+one)zeta CASE(20) xii=(/-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0/) etai=(/-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1/) zetai=(/-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1/) DO l=1,20 xi0=xixii(l) eta0=etaetai(l) zeta0=zetazetai(l) IF(l4.OR.l8.OR.l16.OR.l20)THEN der(1,l)=-pt5xi*(one+eta0)(one+zeta0) der(2,l)=pt25etai(l)(one-xixi)(one+zeta0) der(3,l)=pt25zetai(l)(one-xixi)(one+eta0) ELSE IF(l>=9.AND.l<=12)THEN der(1,l)=pt25xii(l)(one-etaeta)(one+zeta0) der(2,l)=-pt5eta*(one+xi0)(one+zeta0) der(3,l)=pt25zetai(l)(one+xi0)(one-etaeta) ELSE IF(l2.OR.l6.OR.l14.OR.l18) THEN der(1,l)=pt25xii(l)(one+eta0)(one-zetazeta) der(2,l)=pt25etai(l)(one+xi0)(one-zetazeta) der(3,l)=-pt5zeta*(one+xi0)(one+eta0) ELSE der(1,l)=pt125xii(l)(one+eta0)(one+zeta0)* & (twoxi0+eta0+zeta0-one) der(2,l)=pt125etai(l)(one+xi0)(one+zeta0)* & (xi0+twoeta0+zeta0-one) der(3,l)=pt125zetai(l)(one+xi0)(one+eta0)* & (xi0+eta0+twozeta0-one) END IF END DO CASE DEFAULT WRITE(,)'wrong number of nodes in shape_der'
END SELECT CASE DEFAULT WRITE(,*)'wrong number of dimensions in shape_der' END SELECT RETURN END SUBROUTINE shape_der

SUBROUTINE shape_fun(fun,points,i) ! ! This subroutine computes the values of the shape functions. ! to local coordinates ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) INTEGER,INTENT(in)::i REAL(iwp),INTENT(IN)::points(:,:) REAL(iwp),INTENT(OUT)::fun(:) REAL(iwp)::eta,xi,etam,etap,xim,xip,zetam,zetap,c1,c2,c3
REAL(iwp)t1,t2,t3,t4,t5,t6,t7,t8,t9 REAL(iwp)zeta,xi0,eta0,zeta0 INTEGERxii(20),etai(20),zetai(20),l,ndim,nod REAL,PARAMETERpt125=0.125_iwp,pt25=0.25_iwp,pt5=0.5_iwp,pt75=0.75_iwp, & one=1.0_iwp,two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d8=8.0_iwp,d9=9.0_iwp, & d16=16.0_iwp,d27=27.0_iwp,d32=32.0_iwp,d64=64.0_iwp,d128=128.0_iwp ndim=UBOUND(points,2) nod=UBOUND(fun,1)
SELECT CASE(ndim) CASE(1) ! one dimensional case xi=points(i,1) SELECT CASE(nod) CASE(2) t1=-one-xi t2= one-xi fun(1)=t2/two fun(2)=-t1/two CASE(3) t1=-one-xi t2=-xi t3=one-xi fun(1)=t2t3/two fun(2)=-t1t3 fun(3)=t1t2/two CASE(4) t1=-one-xi t2=-one/d3-xi t3=one/d3-xi t4=one-xi fun(1)=t2t3t4d9/d16
fun(2)=-t1t3t4d27/d16 fun(3)=t1t2t4d27/d16 fun(4)=-t1t2t3d9/d16 CASE(5) t1=-one -xi t2=-pt5-xi t3=-xi t4=pt5-xi t5=one-xi fun(1)=t2t3t4t5two/d3 fun(2)=-t1t3t4t5d8/d3 fun(3)=t1t2t4t5d4 fun(4)=-t1t2t3t5d8/d3 fun(5)=t1t2t3t4two/d3 CASE DEFAULT WRITE(,)'wrong number of nodes in shape_fun' END SELECT CASE(2) ! two dimensional case c1=points(i,1) c2=points(i,2) c3=one-c1-c2 xi=points(i,1) eta=points(i,2) etam=pt25(one-eta) etap=pt25*(one+eta) xim=pt25*(one-xi) xip=pt25*(one+xi) SELECT CASE(nod) CASE(3) fun = (/c1,c3,c2/)
CASE(6) fun(1)=(twoc1-one)c1 fun(2)=d4c3c1 fun(3)=(twoc3-one)c3 fun(4)=d4c2c3
fun(5)=(twoc2-one)c2 fun(6)=d4c1c2 CASE(10) fun(1)= ((d3c1-one)(d3c1-two)c1)/two fun(2)= -(d9(d3c1-one)(c1+c2-one)c1)/two fun(3)= (d9(d3c1+d3c2-two)(c1+c2-one)c1)/two fun(4)=-((d3c1+d3c2-one)(d3c1+d3c2-two)(c1+c2-one))/two
fun(5)= (d9(d3c1+d3c2-two)(c1+c2-one)c2)/two fun(6)= -(d9(c1+c2-one)(d3c2-one)c2)/two fun(7)= ((d3c2-one)(d3c2-two)c2)/two fun(8)= (d9(d3c2-one)c1c2)/two fun(9)= (d9*(d3c1-one)c1c2)/two fun(10)=-d27((c2-one)+c1)c1c2 CASE(15) t1=c1-pt25
t2=c1-pt5 t3=c1-pt75
t4=c2-pt25 t5=c2-pt5
t6=c2-pt75 t7=c3-pt25
t8=c3-pt5 t9=c3-pt75 fun(1)=d32/d3c1t1t2t3
fun(2)=d128/d3c3c1t1t2 fun(3)=d64c3c1t1t7
fun(4)=d128/d3c3c1t7t8 fun(5)=d32/d3c3t7t8t9
fun(6)=d128/d3c2c3t7t8 fun(7)=d64c2c3t4t7
fun(8)=d128/d3c2c3t4t5 fun(9)=d32/d3c2t4t5t6
fun(10)=d128/d3c1c2t4t5 fun(11)=d64c1c2t1t4
fun(12)=d128/d3c1c2t1t2 fun(13)=d128c1c2t1c3
fun(15)=d128c1c2c3t4 fun(14)=d128c1c2c3t7
CASE(4) fun=(/d4ximetam,d4ximetap,d4xipetap,d4xipetam/) CASE(8) fun=(/d4etamxim*(-xi-eta-one),d32etamximetap, & d4etapxim(-xi+eta-one),d32ximxipetap, & d4etapxip(xi+eta-one), d32etapxipetam, & d4xipetam(xi-eta-one), d32ximxipetam/) CASE(9) etam=eta-one etap=eta+one xim=xi-one xip=xi+one fun=(/pt25xiximetaetam,-pt5xiximetapetam, & pt25xiximetaetap,-pt5xipximetaetap, & pt25xixipetaetap,-pt5xixipetapetam, & pt25xixipetaetam,-pt5xipximetaetam, & xipximetapetam/) CASE DEFAULT WRITE(,)'wrong number of nodes in shape_fun' END SELECT CASE(3) ! d3 dimensional case xi=points(i,1) eta=points(i,2) zeta=points(i,3) etam=one-eta xim=one-xi
zetam=one-zeta etap=eta+one xip=xi+one
zetap=zeta+one SELECT CASE(nod) CASE(4) fun(1)=xi
fun(2)=eta fun(3)=zeta fun(4)=one-fun(1)-fun(2)-fun(3) CASE(8) fun=(/pt125ximetamzetam,pt125ximetamzetap, & pt125xipetamzetap,pt125xipetamzetam, & pt125ximetapzetam,pt125ximetapzetap, & pt125xipetapzetap,pt125xipetapzetam/) CASE(14) ! type 6 element fun(1) = (xieta+xizeta+twoxi+etazeta+twoeta+twozeta+two)* & (xi-one)(eta-one)(zeta-one)/d8 fun(2) =-(xieta-xizeta+twoxi-etazeta+twoeta-twozeta+two)* & (xi-one)(eta-one)(zeta+one)/d8 fun(3) =-(xieta-xizeta+twoxi+etazeta-twoeta+twozeta-two)* & (xi+one)(eta-one)(zeta+one)/d8 fun(4) = (xieta+xizeta+twoxi-etazeta-twoeta-twozeta-two)* & (xi+one)(eta-one)(zeta-one)/d8 fun(5) =-(xi+one)(xi-one)(eta-one)(zeta+one)(zeta-one)/two fun(6) =-(xi-one)(eta+one)(eta-one)(zeta+one)(zeta-one)/two fun(7) = (xi+one)(xi-one)(eta+one)(eta-one)(zeta+one)/two fun(8) = (xi+one)(eta+one)(eta-one)(zeta+one)(zeta-one)/two fun(9) =-(xi+one)(xi-one)(eta+one)(eta-one)(zeta-one)/two fun(10)= (xieta-xizeta-twoxi+etazeta+twoeta-twozeta-two)* & (xi-one)(eta+one)(zeta-one)/d8 fun(11)=-(xieta+xizeta-twoxi-etazeta+twoeta+twozeta-two)* & (xi-one)(eta+one)(zeta+one)/d8 fun(12)=-(xieta+xizeta-twoxi+etazeta-twoeta-twozeta+two)* & (xi+one)(eta+one)(zeta+one)/d8 fun(13)= (xieta-xizeta-twoxi-etazeta-twoeta+twozeta+two)* & (xi+one)(eta+one)(zeta-one)/d8 fun(14)= (xi+one)(xi-one)(eta+one)(zeta+one)(zeta-one)/two CASE(20) xii=(/-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0/) etai=(/-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1/) zetai=(/-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1/) DO l=1,20 xi0=xixii(l) eta0=etaetai(l) zeta0=zetazetai(l) IF(l4.OR.l8.OR.l16.OR.l20)THEN fun(l)=pt25(one-xixi)(one+eta0)(one+zeta0) ELSE IF(l>=9.AND.l<=12)THEN fun(l)=pt25(one+xi0)(one-etaeta)(one+zeta0) ELSE IF(l2.OR.l6.OR.l14.OR.l18)THEN fun(l)=pt25(one+xi0)(one+eta0)(one-zetazeta) ELSE fun(l)=pt125(one+xi0)(one+eta0)(one+zeta0)(xi0+eta0+zeta0-2) END IF END DO CASE DEFAULT WRITE(,)'wrong number of nodes in shape_fun' END SELECT CASE DEFAULT WRITE(,*)'wrong number of dimensions in shape_fun' END SELECT RETURN END SUBROUTINE shape_fun

SUBROUTINE invert(matrix) ! ! This subroutine inverts a small square matrix onto itself. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN OUT)::matrix(:,:) REAL(iwp)det,j11,j12,j13,j21,j22,j23,j31,j32,j33,con INTEGERndim,i,k ndim=UBOUND(matrix,1) IF(ndim2)THEN det=matrix(1,1)*matrix(2,2)-matrix(1,2)*matrix(2,1) j11=matrix(1,1) matrix(1,1)=matrix(2,2) matrix(2,2)=j11 matrix(1,2)=-matrix(1,2) matrix(2,1)=-matrix(2,1) matrix=matrix/det ELSE IF(ndim3)THEN det=matrix(1,1)*(matrix(2,2)*matrix(3,3)-matrix(3,2)matrix(2,3)) det=det-matrix(1,2)(matrix(2,1)*matrix(3,3)-matrix(3,1)matrix(2,3)) det=det+matrix(1,3)(matrix(2,1)*matrix(3,2)-matrix(3,1)*matrix(2,2)) j11=matrix(2,2)*matrix(3,3)-matrix(3,2)*matrix(2,3) j21=-matrix(2,1)*matrix(3,3)+matrix(3,1)*matrix(2,3) j31=matrix(2,1)*matrix(3,2)-matrix(3,1)*matrix(2,2) j12=-matrix(1,2)*matrix(3,3)+matrix(3,2)*matrix(1,3) j22=matrix(1,1)*matrix(3,3)-matrix(3,1)*matrix(1,3) j32=-matrix(1,1)*matrix(3,2)+matrix(3,1)*matrix(1,2) j13=matrix(1,2)*matrix(2,3)-matrix(2,2)*matrix(1,3) j23=-matrix(1,1)*matrix(2,3)+matrix(2,1)*matrix(1,3) j33=matrix(1,1)*matrix(2,2)-matrix(2,1)*matrix(1,2) matrix(1,1)=j11 matrix(1,2)=j12 matrix(1,3)=j13 matrix(2,1)=j21 matrix(2,2)=j22 matrix(2,3)=j23 matrix(3,1)=j31 matrix(3,2)=j32 matrix(3,3)=j33 matrix=matrix/det ELSE DO k=1,ndim con=matrix(k,k) matrix(k,k)=1.0_iwp matrix(k,:)=matrix(k,:)/con DO i=1,ndim IF(i/=k)THEN con=matrix(i,k) matrix(i,k)=0.0_iwp matrix(i,:)=matrix(i,:)-matrix(k,:)*con END IF END DO END DO END IF RETURN END SUBROUTINE invert

SUBROUTINE fsparv(kv,km,g,kdiag) ! ! This subroutine assembles element matrices into a symmetric skyline ! global matrix. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) INTEGER,INTENT(IN)::g(:),kdiag(:) REAL(iwp),INTENT(IN)::km(:,:) REAL(iwp),INTENT(OUT)kv(:) INTEGERi,idof,k,j,iw,ival idof=UBOUND(g,1) DO i=1,idof k=g(i) IF(k/=0)THEN DO j=1,idof IF(g(j)/=0)THEN iw=k-g(j) IF(iw>=0)THEN ival=kdiag(k)-iw kv(ival)=kv(ival)+km(i,j) END IF END IF END DO END IF END DO RETURN END SUBROUTINE fsparv

SUBROUTINE sparin(kv,kdiag) ! ! This subroutine performs Cholesky factorisation on a symmetric ! skyline global matrix. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN OUT)::kv(:) INTEGER,INTENT(IN)kdiag(:) INTEGERn,i,ki,l,kj,j,ll,m,k REAL(iwp)::x n=UBOUND(kdiag,1)
kv(1)=SQRT(kv(1)) DO i=2,n ki=kdiag(i)-i l=kdiag(i-1)-ki+1 DO j=l,i x=kv(ki+j) kj=kdiag(j)-j IF(j/=1)THEN ll=kdiag(j-1)-kj+1 ll=max(l,ll) IF(ll/=j)THEN m=j-1 DO k=ll,m x=x-kv(ki+k)*kv(kj+k) END DO END IF END IF kv(ki+j)=x/kv(kj+j) END DO kv(ki+i)=SQRT(x) END DO RETURN END SUBROUTINE sparin

SUBROUTINE spabac(kv,loads,kdiag) ! ! This subroutine performs Cholesky forward and back-substitution ! on a symmetric skyline global matrix. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN)::kv(:) REAL(iwp),INTENT(IN OUT)::loads(0:) INTEGER,INTENT(IN)kdiag(:) INTEGERn,i,ki,l,m,j,it,k REAL(iwp)::x n=UBOUND(kdiag,1) loads(1)=loads(1)/kv(1) DO i=2,n ki=kdiag(i)-i l=kdiag(i-1)-ki+1 x=loads(i) IF(l/=i)THEN m=i-1 DO j=l,m x=x-kv(ki+j)loads(j) END DO END IF loads(i)=x/kv(ki+i) END DO DO it=2,n i=n+2-it ki=kdiag(i)-i x=loads(i)/kv(ki+i) loads(i)=x l=kdiag(i-1)-ki+1 IF(l/=i)THEN m=i-1 DO k=l,m loads(k)=loads(k)-xkv(ki+k) END DO END IF END DO loads(1)=loads(1)/kv(1) RETURN END SUBROUTINE spabac

SUBROUTINE linmul_sky(kv,disps,loads,kdiag) ! ! This subroutine forms the product of symmetric matrix stored as ! a skyline and a vector. ! IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN)::kv(:),disps(0:) REAL(iwp),INTENT(OUT)::loads(0:) INTEGER,INTENT(IN)kdiag(:) INTEGERn,i,j,low,lup,k REAL(iwp)::x,zero=0.0_iwp n=UBOUND(disps,1) DO i=1,n x=zero lup=kdiag(i) IF(i1)low=lup IF(i/=1)low=kdiag(i-1)+1 DO j=low,lup x=x+kv(j)*disps(i+j-lup) END DO loads(i)=x IF(i1)CYCLE
lup=lup-1 DO j=low,lup k=i+j-lup-1 loads(k)=loads(k)+kv(j)*disps(i)
END DO END DO RETURN END SUBROUTINE linmul_sky

SUBROUTINE contour(loads,g_coord,g_num,ned,ips) ! ! This subroutine produces a PostScript output file '.con' displaying ! a contour map of loads over the finite element mesh. ! IMPLICIT NONE INTEGER,PARAMETERiwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN) loads(0:),g_coord(:,:) INTEGER,INTENT(IN)g_num(:,:),ned,ips REAL(iwp)xmin,xmax,ymin,ymax,width,height,scale=72,sxy,xo,yo REAL(iwp)pmin,pmax,ratio,x12,y12,x23,y23,x34,y34,x41,y41,x1,y1 REAL(iwp)pt5=0.5_iwp,elvn=11.0_iwp,ept5=8.5_iwp,eight=8.0_iwp, & onept5=1.5_iwp,fpt5=5.5_iwp LOGICALs12,s23,s34,s41,draw INTEGER,ALLOCATABLEcorner(:,:) REAL(iwp),ALLOCATABLEcont(:) INTEGERi,j,k,l,nn,nels,nci,nod,ns,i1,i2,j1,j2 ! OPEN(ips,FILE='fe95.con') ! ! compute size of mesh ! nn=UBOUND(g_coord,2) nod= UBOUND(g_num,1) nels=UBOUND(g_num,2) xmin=g_coord(1,1) xmax=g_coord(1,1) ymin=g_coord(2,1) ymax=g_coord(2,1) DO i=2,nn IF(g_coord(1,i)<xmin)xmin=g_coord(1,i)
IF(g_coord(1,i)>xmax)xmax=g_coord(1,i)
IF(g_coord(2,i)<ymin)ymin=g_coord(2,i)
IF(g_coord(2,i)>ymax)ymax=g_coord(2,i)
END DO width =xmax-xmin height=ymax-ymin ! ! allow 1.5' margin minimum on each side of figure ! ! portrait mode ! IF(height.GE.elvn/ept5width)THEN ! ! height governs the scale ! sxy=scaleeight/height xo=scalept5*(ept5-eightwidth/height) yo=scaleonept5 ELSE ! ! width governs the scale ! sxy=scalefpt5/width xo=scaleonept5 yo=scalept5(elvn-fpt5height/width) END IF ! ! find range of potentials and contour values ! nci=ned+1 pmin=MINVAL(loads(1:)) pmax=MAXVAL(loads(1:)) ALLOCATE(cont(nci)) cont(1)=pmin DO i=2,nci cont(i)=cont(i-1)+(pmax-pmin)/ned END DO ! ! start PostScript output ! WRITE(ips,'(a)')'%!PS-Adobe-1.0' WRITE(ips,'(a)')'%%DocumentFonts: none' WRITE(ips,'(a)')'%%Pages: 1' WRITE(ips,'(a)')'%%EndComments' WRITE(ips,'(a)')'/m def' WRITE(ips,'(a)')'/l def' WRITE(ips,'(a)')'/s def' WRITE(ips,'(a)')'%%EndProlog' WRITE(ips,'(a)')'%%Page: 0 1' WRITE(ips,'(a)')'gsave' ! WRITE(ips,'(2f9.2,a)') xo, yo, ' translate' WRITE(ips,'(f9.2,a)') 1.0, ' setlinewidth' ! ! draw the mesh outline ! IF(nod3.OR.nod6.OR.nod10.OR.nod15)ns=3 IF(nod4.OR.nod8.OR.nod9)ns=4 ALLOCATE(corner(ns,2)) IF(nod 3)corner=RESHAPE((/1,2,3,2,3,1/),(/3,2/)) IF(nod== 6)corner=RESHAPE((/1,3,5,3,5,1/),(/3,2/)) IF(nod10)corner=RESHAPE((/1,4,7,4,7,1/),(/3,2/)) IF(nod15)corner=RESHAPE((/1,5,9,5,9,1/),(/3,2/)) IF(nod== 4)corner=RESHAPE((/1,2,3,4,2,3,4,1/),(/4,2/)) IF(nod== 8)corner=RESHAPE((/1,3,5,7,3,5,7,1/),(/4,2/)) IF(nod== 9)corner=RESHAPE((/1,3,5,7,3,5,7,1/),(/4,2/)) DO i=1,nels DO j=1,ns draw=.TRUE. i1=g_num(corner(j,1),i) i2=g_num(corner(j,2),i) DO k=1,nels DO l=1,ns j1=g_num(corner(l,1),k) j2=g_num(corner(l,2),k) IF((i1j2).AND.(i2j1))THEN draw=.FALSE. EXIT END IF END DO IF(.NOT.draw)EXIT END DO IF(draw)THEN x1=sxy(g_coord(1,i1)-xmin) y1=sxy*(g_coord(2,i1)-ymin) WRITE(ips,'(2f9.2,a)') x1, y1,' m' x1=sxy*(g_coord(1,i2)-xmin) y1=sxy*(g_coord(2,i2)-ymin) WRITE(ips,'(2f9.2,a)') x1, y1,' l' WRITE(ips,'(a)')' s' END IF END DO END DO ! ! check intersection of contours with each element ! WRITE(ips,'(f9.2,a)') 0.5, ' setlinewidth' DO i=2,nci-1 DO j=1,nels s12=.FALSE. s23=.FALSE. s34=.FALSE. s41=.FALSE. IF((loads(g_num(1,j))<=cont(i).AND.loads(g_num(2,j))>cont(i)).OR. & (loads(g_num(2,j))<=cont(i).AND.loads(g_num(1,j))>cont(i))) & s12=.TRUE. IF((loads(g_num(2,j))<=cont(i).AND.loads(g_num(3,j))>cont(i)).OR. & (loads(g_num(3,j))<=cont(i).AND.loads(g_num(2,j))>cont(i))) & s23=.TRUE. IF((loads(g_num(3,j))<=cont(i).AND.loads(g_num(4,j))>cont(i)).OR. & (loads(g_num(4,j))<=cont(i).AND.loads(g_num(3,j))>cont(i))) & s34=.TRUE. IF((loads(g_num(4,j))<=cont(i).AND.loads(g_num(1,j))>cont(i)).OR. & (loads(g_num(1,j))<=cont(i).AND.loads(g_num(4,j))>cont(i))) & s41=.TRUE. IF(s12)THEN ratio=(cont(i)-loads(g_num(1,j)))/ & (loads(g_num(2,j))-loads(g_num(1,j))) x12=sxy*(g_coord(1,g_num(1,j))+ & ratio*(g_coord(1,g_num(2,j))-g_coord(1,g_num(1,j)))-xmin) y12=sxy*(g_coord(2,g_num(1,j))+ & ratio*(g_coord(2,g_num(2,j))-g_coord(2,g_num(1,j)))-ymin) END IF IF(s23)THEN ratio=(cont(i)-loads(g_num(2,j)))/ & (loads(g_num(3,j))-loads(g_num(2,j))) x23=sxy*(g_coord(1,g_num(2,j))+ & ratio*(g_coord(1,g_num(3,j))-g_coord(1,g_num(2,j)))-xmin) y23=sxy*(g_coord(2,g_num(2,j))+ & ratio*(g_coord(2,g_num(3,j))-g_coord(2,g_num(2,j)))-ymin) END IF IF(s34)THEN ratio=(cont(i)-loads(g_num(3,j)))/ & (loads(g_num(4,j))-loads(g_num(3,j))) x34=sxy*(g_coord(1,g_num(3,j))+ & ratio*(g_coord(1,g_num(4,j))-g_coord(1,g_num(3,j)))-xmin) y34=sxy*(g_coord(2,g_num(3,j))+ & ratio*(g_coord(2,g_num(4,j))-g_coord(2,g_num(3,j)))-ymin) END IF IF(s41)THEN ratio=(cont(i)-loads(g_num(4,j)))/ & (loads(g_num(1,j))-loads(g_num(4,j))) x41=sxy*(g_coord(1,g_num(4,j))+ & ratio*(g_coord(1,g_num(1,j))-g_coord(1,g_num(4,j)))-xmin) y41=sxy*(g_coord(2,g_num(4,j))+ & ratio*(g_coord(2,g_num(1,j))-g_coord(2,g_num(4,j)))-ymin) END IF ! ! draw contours ! IF(s12)THEN WRITE(ips,'(2f9.2,a)')x12,y12,' m' IF(s23)WRITE(ips,'(2f9.2,a)')x23,y23,' l s' IF(s34)WRITE(ips,'(2f9.2,a)')x34,y34,' l s' IF(s41)WRITE(ips,'(2f9.2,a)')x41,y41,' l s' END IF IF(s23)THEN WRITE(ips,'(2f9.2,a)')x23,y23,' m' IF(s34)WRITE(ips,'(2f9.2,a)')x34,y34,' l s' IF(s41)WRITE(ips,'(2f9.2,a)')x41,y41,' l s' END IF IF(s34.AND.s41)THEN WRITE(ips,'(2f9.2,a)')x34,y34,' m' WRITE(ips,'(2f9.2,a)')x41,y41,' l s' END IF END DO END DO ! close output file WRITE(ips,'(a)')'grestore' WRITE(ips,'(a)')'showpage' CLOSE(ips) ! RETURN END SUBROUTINE contour END PROGRAM p72

Can you give us a clue: on which line do you get the error?

Perhaps you could include the program as code and not as normal test. This makes reading it easier. 1.) The begin of the code block starts with [code] 2.) The code is ended by [/code]

Are you developing this code, as the error 775 is probably due to using an END rather than an END DO statement.

Your use of ALLOCATE is beyond my experience, although it may be valid. I would not have used the ALLOCATE statements like: .... x_coords(nxe+1),y_coords(nye+1) ... or ALLOCATE(kv(kdiag(neq)),kvh(kdiag(neq)))

and preferred each array in a seperate allocate statement

Your use of arrays appears luxurious by my old practices, but can work for a smallish problem.

!-----------------------dynamic arrays------------------------------------ 
INTEGER,ALLOCATABLE::  &
  etype(:),     &  ! nels
  g_num(:,:),   &  ! nod,nels
  kdiag(:),     &  ! neq
  node(:),      &  ! fixed_freedoms
  num(:)           ! nod

REAL(iwp),ALLOCATABLE::  &
  coord(:,:),   &  ! nod,ndim
  der(:,:),     &  ! ndim,nod
  deriv(:,:),   &  ! ndim,nod
  disps(:),     &  ! 0:neq
  fun(:),       &  ! nod
  gc(:),        &  ! ndim
  g_coord(:,:), &  ! ndim,nn
  jac(:,:),     &  ! ndim,ndim
  kay(:,:),     &  ! ndim,ndim
  kc(:,:),      &  ! nod,nod
  kv(:),        &  ! kdiag(neq)   ?? I'd use a variable nstif = kdiag(neq)
  kvh(:),       &  ! kdiag(neq)   ??
  loads(:),     &  ! 0:neq
  points(:,:),  &  ! nip,ndim
  prop(:,:),    &  ! ndim,np_types
  value(:),     &  ! fixed_freedoms
  weights(:),   &  ! nip
  x_coords(:)   &  ! nxe+1        ?? I'd use a variable nxe1 = nxe+1
  y_coords(:)      ! nye+1        ??

I'd also import this into excel and calculate the byte size of each array to check where the memory might disappear.

John