!     Last change:  DV    7 Apr 2005   11:54 am
PROGRAM test_8
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp)::coord(8,2),km(16,16),ym,pr,two=2.0_iwp
 INTEGER::i,j,nlen
 CHARACTER(LEN=15)::argv
 CALL getname(argv,nlen)
 OPEN(10,FILE=argv(1:nlen)//'.dat')
 OPEN(11,FILE=argv(1:nlen)//'.res')
 WRITE(11,'(a/)')"      PLANE STRAIN 8-NODE QUADRILATERAL"
 WRITE(11,'(a)')"      Nodal Coordinates"
 DO i=1,7,2
   READ(10,*)coord(i,:)
 END DO
 DO j=1,2
   coord(2,j)=(coord(1,j)+coord(3,j))/two
   coord(4,j)=(coord(3,j)+coord(5,j))/two
   coord(6,j)=(coord(5,j)+coord(7,j))/two
   coord(8,j)=(coord(7,j)+coord(1,j))/two
 END DO
 DO i=1,8
   WRITE(11,'(2F12.4)')coord(i,:)
 END DO
 WRITE(11,'(/a)')"      Young's Modulus and Poisson's Ratio"
 READ(10,*)ym,pr
 WRITE(11,'(/2F12.4)')ym,pr
 CALL stiff8(coord,ym,pr,km)
 WRITE(11,'(/a)')"      Stiffness Matrix"
 DO i=1,16
 DO j=1,i
  WRITE(11,'(2I3,F18.4)')i,j,km(i,j)
 END DO
 END DO
STOP
CONTAINS
!
SUBROUTINE stiff8(COORD,E,V,KM)
!
! This subroutine generates the "semi-analytical" stiffness matrix
! for an eight node quadrilateral in plane strain based on four
! Gauss-Legendre integrating points (nip=4).
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::COORD(8,2),E,V
 REAL(iwp),INTENT(OUT)::KM(16,16)
 REAL(iwp)::E1,E2,G,A2,F1,F2,S1
 REAL(iwp)::X51ORG,X57ORG,X73ORG,Y51ORG,Y57ORG,Y73ORG
 REAL(iwp)::X51,X73,X57,Y51,Y73,Y57,X1,X3,X5,X7,Y1,Y3,Y5,Y7
 REAL(iwp)::zero=0.0_iwp,one=1.0_iwp,two=2.0_iwp
 INTEGER::I,J
 X1=COORD(1,1)
 X3=COORD(3,1)
 X5=COORD(5,1)
 X7=COORD(7,1)
 Y1=COORD(1,2)
 Y3=COORD(3,2)
 Y5=COORD(5,2)
 Y7=COORD(7,2)
 E1=E*(one-V)/(one+V)/(one-two*V)
 E2=V*E1/(one-V)
 G=E/two/(one+V)
 X73=X7-X3
 X51=X5-X1
 X57=X5-X7
 Y73=Y7-Y3
 Y51=Y5-Y1
 Y57=Y5-Y7
 X73ORG=X73
 X51ORG=X51
 X57ORG=X57
 Y73ORG=Y73
 Y51ORG=Y51
 Y57ORG=Y57
 A2=Y51*X73-Y73*X51
 F1=two*Y73*X57-two*X73*Y57+A2
 F2=two*Y51*X57-two*X51*Y57+A2
 S1=E1*Y73**2+G*X73**2
 CALL FAM1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 S1=E1*Y73**2+G*X73**2
 CALL FAM9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 CALL FAM16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 CALL ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 DO I=1,15
   DO J=I+1,16
     KM(I,J)=KM(J,I)
   END DO
 END DO
RETURN
END SUBROUTINE stiff8

SUBROUTINE ORG(COORD,A2,X51,X57,X73,Y51,Y57,Y73,F1,F2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::COORD(:,:),A2
 REAL(iwp),INTENT(OUT)::X51,X57,X73,Y51,Y57,Y73,F1,F2
 REAL(iwp)::X1,X3,X5,X7,Y1,Y3,Y5,Y7,two=2.0_iwp
 X1=COORD(1,1)
 X3=COORD(3,1)
 X5=COORD(5,1)
 X7=COORD(7,1)
 Y1=COORD(1,2)
 Y3=COORD(3,2)
 Y5=COORD(5,2)
 Y7=COORD(7,2)
 X73=X7-X3
 X51=X5-X1
 X57=X5-X7
 Y73=Y7-Y3
 Y51=Y5-Y1
 Y57=Y5-Y7
 F1=two*Y73*X57-two*X73*Y57+A2
 F2=two*Y51*X57-two*X51*Y57+A2
RETURN
END SUBROUTINE ORG

SUBROUTINE FAM1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(1,1)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(5,5)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(9,9)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(13,13)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(14,14)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(10,10)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(6,6)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2 + G*X73**2
 KM(2,2)=DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM1

FUNCTION DAD1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD1)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,DD1
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d36=36.0_iwp
 DIAG=(d3*(two*A2-F1)*S1)/(two*(d3*A2**2-F1**2))
 A=two*X51+two*X57+X73
 B=two*Y51+two*Y57+Y73
 OPPO=(d3*A2*S1+A*(A*A2-two*X73*F2)*G+                                    &
       B*(B*A2-two*Y73*F2)*E1)/(d36*(d3*A2**2-F2**2))
 DD1=DIAG+OPPO
RETURN
END FUNCTION DAD1

SUBROUTINE FAM2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(2,1)=DAD2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(6,5)=DAD2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,9)=DAD2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,13)=DAD2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 RETURN
END SUBROUTINE FAM2

FUNCTION DAD2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD2)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,C,D,DD2
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d18=18.0_iwp
 DIAG=-d3*(E2+G)*(two*A2-F1)*X73*Y73/(two*(d3*A2**2-F1**2))
 C=(two*Y51-Y57+Y73)*X73-F1+d3*X57*Y73
 D=(two*Y51+two*Y57+d3*Y73)*X57+(two*Y51-Y57+two*Y73)*X73+                &
    two*(Y51+Y57)*X51-F1
 OPPO=-(E2+G)*(A2*D-C*F2)/(d18*(d3*A2**2-F2**2))
 DD2=DIAG+OPPO
RETURN
END FUNCTION DAD2

SUBROUTINE FAM3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(3,1)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(7,5)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,9)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,13)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,12)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,8)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(6,4)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,2)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(15,1)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(13,11)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(9,7)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(5,3)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(8,6)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(12,10)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(16,14)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(4,2)=DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM3

FUNCTION DAD3(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD3)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN):: X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD3
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp
 REAL(iwp)::d6=6.0_iwp,d7=7.0_iwp,d9=9.0_iwp,d10=10.0_iwp,d12=12.0_iwp
 A=(d12*A2-(d7*F1))*X51-((d9*A2-(d5*F1))*X57)+X73*F1
 B=(d12*A2-(d7*F1))*Y51-((d9*A2-(d5*F1))*Y57)+Y73*F1
 DIAG=-(A*G*X73+B*Y73*E1)/(d3*(d3*A2**2-F1**2))
 C=X51**2*A2+X51**2*F2+d4*X51*X73*A2-(d3*X51*X73*F2)+d5*X51*X57*          &
   A2-(X51*X57*F2)+d7*X73**2*A2-(d4*X73**2*F2)+d10*X73*X57*A2-(           &
   d6*X73*X57*F2)+d4*X57**2*A2-(two*X57**2*F2)
 D=Y51**2*A2+Y51**2*F2+d4*Y51*Y73*A2-(d3*Y51*Y73*F2)+d5*Y51*Y57*          &
   A2-(Y51*Y57*F2)+d7*Y73**2*A2-(d4*Y73**2*F2)+d10*Y73*Y57*A2-(           &
   d6*Y73*Y57*F2)+d4*Y57**2*A2-(two*Y57**2*F2)
 OPPO=-(C*G+D*E1)/(d9*(d3*A2**2-F2**2))
 DD3=DIAG+OPPO
RETURN
END FUNCTION DAD3

SUBROUTINE FAM4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(4,1)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(8,5)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(12,9)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,13)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,11)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,7)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(6,3)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,2)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(16,1)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(13,12)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(9,8)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(5,4)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(7,6)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(11,10)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(15,14)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(3,2)=DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM4

FUNCTION DAD4(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD4)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN):: X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,TOT,DE,DD4
 REAL(iwp)::zero=0.0_iwp,one=1.0_iwp,two=2.0_iwp,d3=3.0_iwp
 REAL(iwp)::d4=4.0_iwp,d8=8.0_iwp,d12=12.0_iwp,d24=24.0_iwp,sq3
 REAL(iwp)::d5=5.0_iwp
 sq3=SQRT(d3)
 TOT=zero
 DE=(sq3*A2-F1)/(d8*sq3)
 A=-((sq3-one)*Y73)/d8
 B=((sq3-one)*X73)/d8
 C=-(d3*sq3*Y51-(two*sq3*Y57)-(sq3*Y73)-(d5*Y51)+d4*Y57-Y73)/(d12*sq3)
 D=(d3*sq3*X51-(two*sq3*X57)-(sq3*X73)-(d5*X51)+d4*X57-X73)/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DE=(sq3*A2+F1)/(d8*sq3)
 A=((sq3+one)*Y73)/d8
 B=-((sq3+one)*X73)/d8
 C=-(d3*sq3*Y51-(two*sq3*Y57)-(sq3*Y73)+d5*Y51-                           &
      (d4*Y57)+Y73)/(d12*sq3)
 D=(d3*sq3*X51-(two*sq3*X57)-(sq3*X73)+d5*X51-(d4*X57)+X73)/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DE=(sq3*A2+F2)/(d8*sq3)
 A=-(sq3*Y73+two*Y51+two*Y57+Y73)/(d24*sq3)
 B=(sq3*X73+two*X51+two*X57+X73)/(d24*sq3)
 C=-(sq3*Y51-(two*sq3*Y57)-(d3*sq3*Y73)-Y51-(d4*Y57)-                     &
    (d5*Y73))/(d12*sq3)
 D=(sq3*X51-(two*sq3*X57)-(d3*sq3*X73)-X51-(d4*X57)-                      &
   (d5*X73))/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DE=(sq3*A2-F2)/(d8*sq3)
 A=-(sq3*Y73-(two*Y51)-(two*Y57)-Y73)/(d24*sq3)
 B=(sq3*X73-(two*X51)-(two*X57)-X73)/(d24*sq3)
 C=-(sq3*Y51-(two*sq3*Y57)-(d3*sq3*Y73)+Y51+d4*Y57+d5*Y73)/(d12*sq3)
 D=(sq3*X51-(two*sq3*X57)-(d3*sq3*X73)+X51+d4*X57+d5*X73)/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DD4=TOT
RETURN
END FUNCTION DAD4

SUBROUTINE FAM5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(5,1)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(9,5)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(13,9)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(13,1)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,10)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,6)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(6,2)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,2)=DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM5

FUNCTION DAD5(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD5)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN):: X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD5
 REAL(iwp):: two=2.0_iwp,d3=3.0_iwp,d6=6.0_iwp
 A=(X73-X57)*(d3*A2-F1)+(d3*A2-(two*F1))*X51
 B=(Y73-Y57)*(d3*A2-F1)+(d3*A2-(two*F1))*Y51
 DIAG=(A*G*X73+B*Y73*E1)/(d6*(d3*A2**2-F1**2))
 C=(X51+X57)*(d3*A2-F2)+(d3*A2-(two*F2))*X73
 D=(Y51+Y57)*(d3*A2-F2)+(d3*A2-(two*F2))*Y73
 OPPO=(C*G*X51+D*Y51*E1)/(d6*(d3*A2**2-F2**2))
 DD5=DIAG+OPPO
RETURN
END FUNCTION DAD5

SUBROUTINE FAM6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(6,1)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,5)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,9)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(13,2)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(14,1)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(5,2)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(9,6)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(13,10)=DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM6

FUNCTION DAD6(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD6)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,DD6
 REAL(iwp)::one=1.0_iwp,two=2.0_iwp,d3=3.0_iwp,d24=24.0_iwp,sq3
 sq3=SQRT(d3)
 DIAG=(-(((sq3*X51+X51-two*X57+two*X73)*E2*Y73+(sq3*Y51+Y51-two*Y57+      &
      two*Y73)*G*X73)*(sq3*(A2)-F1)*(sq3+one)-((sq3*X51-X51+two*X57-      &
      two*X73)*E2*Y73+(sq3*Y51-Y51+two*Y57-two*Y73)*G*X73)*(sq3*(A2)+F1)* &
      (sq3-one)))/(d24*(d3*A2**2-F1**2))
 OPPO=(-(((sq3*X73+two*X51+two*X57+X73)*G*Y51+(sq3*Y73+two*Y51+two*       &
      Y57+Y73)*E2*X51)*(sq3*A2-F2)*(sq3+one)-((sq3*X73-two*X51-two*X57-   &
      X73)*G*Y51+(sq3*Y73-two*Y51-two*Y57-Y73)*E2*X51)*(sq3*(A2)+F2)*     &
      (sq3-one)))/(d24*(d3*A2**2-F2**2))
 DD6=DIAG+OPPO
RETURN
END FUNCTION DAD6

SUBROUTINE FAM7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(7,1)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,5)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,9)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(13,3)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,8)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,4)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,6)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(12,2)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(11,1)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(13,7)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(9,3)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(15,5)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(12,6)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(16,10)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(14,4)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(8,2)=DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM7

FUNCTION DAD7(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD7)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD7
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp
 REAL(iwp)::d6=6.0_iwp,d7=7.0_iwp,d9=9.0_iwp,d10=10.0_iwp
 A=(d3*A2+F1)*X57-((d3*A2-(two*F1))*X73)+X51*F1
 B=(d3*A2+F1)*Y57-((d3*A2-(two*F1))*Y73)+Y51*F1
 DIAG=(A*G*X73+B*Y73*E1)/(d3*(d3*A2**2-F1**2))
 C=d5*X51**2*A2-(X51**2*F2)-(X51*X73*A2)+X51*X57*A2+X51*X57*F2-(          &
   d7*X73**2*A2)+d4*X73**2*F2-(d10*X73*X57*A2)+d6*X73*X57*F2-(d4*X57**2*  &
   A2)+two*X57**2*F2
 D=d5*Y51**2*A2-(Y51**2*F2)-(Y51*Y73*A2)+Y51*Y57*A2+Y51*Y57*F2-(          &
   d7*Y73**2*A2)+d4*Y73**2*F2-(d10*Y73*Y57*A2)+d6*Y73*Y57*F2-(d4*Y57**2*  &
   A2)+two*Y57**2*F2
 OPPO=-(C*G+D*E1)/(d9*(d3*A2**2-F2**2))
 DD7=DIAG+OPPO
RETURN
END FUNCTION DAD7

SUBROUTINE FAM8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(8,1)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(12,5)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,9)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(13,4)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,7)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(10,3)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,6)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,2)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(12,1)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(13,8)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(9,4)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(16,5)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(11,6)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(15,10)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(14,3)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(7,2)=DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM8

FUNCTION DAD8(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD8)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN):: X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,TOT,DE,DD8
 REAL(iwp)::zero=0.0_iwp,one=1.0_iwp,two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp
 REAL(iwp)::d5=5.0_iwp,d8=8.0_iwp,d12=12.0_iwp,d24=24.0_iwp,sq3
 sq3=SQRT(d3)
 TOT=zero
 DE=(sq3*A2-F1)/(d8*sq3)
 A=-((sq3-one)*Y73)/d8
 B=((sq3-one)*X73)/d8
 C=-(sq3*Y51+two*sq3*Y57+sq3*Y73+Y51+d4*Y57-Y73)/(d12*sq3)
 D=(sq3*X51+two*sq3*X57+sq3*X73+X51+d4*X57-X73)/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DE=(sq3*A2+F1)/(d8*sq3)
 A=((sq3+one)*Y73)/d8
 B=-((sq3+one)*X73)/d8
 C=-(sq3*Y51+two*sq3*Y57+sq3*Y73-Y51-(d4*Y57)+Y73)/(d12*sq3)
 D=(sq3*X51+two*sq3*X57+sq3*X73-X51-(d4*X57)+X73)/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DE=(sq3*A2+F2)/(d8*sq3)
 A=-(sq3*Y73+two*Y51+two*Y57+Y73)/(d24*sq3)
 B=(sq3*X73+two*X51+two*X57+X73)/(d24*sq3)
 C=(sq3*Y51-(two*sq3*Y57)-(d3*sq3*Y73)+d5*Y51-(d4*Y57)-(d5*Y73))/(d12*sq3)
 D=-(sq3*X51-(two*sq3*X57)-(d3*sq3*X73)+d5*X51-(d4*X57)-(d5*X73))/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DE=(sq3*A2-F2)/(d8*sq3)
 A=-(sq3*Y73-(two*Y51)-(two*Y57)-Y73)/(d24*sq3)
 B=(sq3*X73-(two*X51)-(two*X57)-X73)/(d24*sq3)
 C=(sq3*Y51-(two*sq3*Y57)-(d3*sq3*Y73)-(d5*Y51)+d4*Y57+d5*Y73)/(d12*sq3)
 D=-(sq3*X51-(two*sq3*X57)-(d3*sq3*X73)-(d5*X51)+d4*X57+d5*X73)/(d12*sq3)
 TOT=TOT+(A*D*E2+B*C*G)/DE
 DD8=TOT
RETURN
END FUNCTION DAD8

SUBROUTINE FAM9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 REAL(iwp)::S1
 KM(9,1)=DAD9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 S1=E1*Y73**2+G*X73**2
 KM(13,5)=DAD9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 S1=E1*Y73**2+G*X73**2
 KM(14,6)=DAD9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 S1=E1*Y73**2+G*X73**2
 KM(10,2)=DAD9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM9

FUNCTION DAD9(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD9)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN):: X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,C,D,DD9
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d18=18.0_iwp
 DIAG=(d3*A2*S1)/(two*(d3*A2**2-F1**2))
 C=((A2-F2)*X73+two*X57*A2)*X51-((two*A2-F2)*X73**2)-(two*(A2-F2)*        &
   X73*X57)+d4*X51**2*A2-(two*X57**2*A2)
 D=((A2-F2)*Y73+two*Y57*A2)*Y51-((two*A2-F2)*Y73**2)-(two*(A2-F2)*        &
   Y73*Y57)+d4*Y51**2*A2-(two*Y57**2*A2)
 OPPO=(C*G+D*E1)/(D18*(d3*A2**2-F2**2))
 DD9=DIAG+OPPO
RETURN
END FUNCTION DAD9

SUBROUTINE FAM10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(10,1)=DAD10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(14,5)=DAD10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(13,6)=DAD10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(9,2)=DAD10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM10

FUNCTION DAD10(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD10)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,C,DD10
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d8=8.0_iwp,d36=36.0_iwp
 DIAG=-(d3*(G+E2)*Y73*X73*A2)/(two*(d3*A2**2-F1**2))
 C=two*((two*A2-F2)*X73-((A2-F2)*X51))*Y73+d4*((A2-F2)*X73-(X51*A2)+      &
   X57*A2)*Y57+(F1+F2)*A2-(d8*Y51*X51*A2)-A2**2-(F1*F2)
 OPPO=((G+E2)*C)/(d36*(d3*A2**2-F2**2))
 DD10=DIAG+OPPO
RETURN
END FUNCTION DAD10

SUBROUTINE FAM11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(3,3)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(7,7)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,11)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,15)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(12,12)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(8,8)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(4,4)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,16)=DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM11

FUNCTION DAD11(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD11)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD11
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d7=7.0_iwp,d9=9.0_iwp
 REAL(IWP)::d11=11.0_iwp,d15=15.0_iwp,d19=19.0_iwp,d26=26.0_iwp
 A=two*((d19*A2-(d11*F1))*X57+(two*A2-F1)*X73)*X51-((d26*A2-(d15*F1))*    &
   X51**2)-(two*(d7*A2-(d4*F1))*X57**2)-((two*A2+F1)*X73**2)-(two*(A2-F1)*&
   X73*X57)
 B=-(two*((d19*A2-(d11*F1))*Y57+(two*A2-F1)*Y73)*Y51-((d26*A2-(d15*F1))*  &
     Y51**2)-(two*(d7*A2-(d4*F1))*Y57**2)-((two*A2+F1)*Y73**2)-(two*(A2-  &
     F1)*Y73*Y57))
 DIAG=-(two*(A*G-(B*E1)))/(d9*(d3*A2**2-F1**2))
 C=two*((two*A2-F2)*X73+(A2-F2)*X57)*X51-((d26*A2-(d15*F2))*X73**2)-(two* &
   (d19*A2-(d11*F2))*X73*X57)-(two*(d7*A2-(d4*F2))*X57**2)-((two*A2+F2)*  &
   X51**2)
 D=-(two*((two*A2-F2)*Y73+(A2-F2)*Y57)*Y51-((d26*A2-(d15*F2))*Y73**2)-    &
    (two*(d19*A2-(d11*F2))*Y73*Y57)-(two*(d7*A2-(d4*F2))*Y57**2)-((two*A2+&
    F2)*Y51**2))
 OPPO=-(two*(C*G-(D*E1)))/(d9*(d3*A2**2-F2**2))
 DD11=DIAG+OPPO
RETURN
END FUNCTION DAD11

SUBROUTINE FAM12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(4,3)=DAD12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(8,7)=DAD12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(12,11)=DAD12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,15)=DAD12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM12

FUNCTION DAD12(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD12)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,C,DD12
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d7=7.0_iwp,d9=9.0_iwp
 REAL(IWP)::d11=11.0_iwp,d15=15.0_iwp,d19=19.0_iwp,d26=26.0_iwp
 REAL(iwp)::d14=14.0_iwp,d22=22.0_iwp
 A=d4*((d19*A2-(d11*F1))*X51-((d7*A2-(d4*F1))*X57)-((A2-F1)*X73))*Y57-    &
   (two*((two*A2+F1)*X73-(two*(two*A2-F1)*X51))*Y73)-(two*(d26*A2-        &
   (d15*F1))*Y51*X51)+(d7*F1+d19*F2)*A2-(d14*A2**2)+F1**2-(d11*F1*F2)
 DIAG=((G+E2)*A)/(d9*(d3*A2**2-F1**2))
 C=two*((d26*A2-(d15*F2))*X73-(two*(two*A2-F2)*X51))*Y73+d4*((d19*A2-     &
   (d11*F2))*X73+(d7*A2-(d4*F2))*X57-((A2-F2)*X51))*Y57+two*(two*A2+F2)*  &
   Y51*X51+(d19*F1+d11*F2)*A2-(d22*A2**2)-(d11*F1*F2)+F2**2
 OPPO=-((G+E2)*C)/(d9*(d3*A2**2-F2**2))
 DD12=DIAG+OPPO
RETURN
END FUNCTION DAD12

 SUBROUTINE FAM13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(7,3)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,7)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,11)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,3)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(16,4)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(8,4)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(12,8)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(16,12)=DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM13

FUNCTION DAD13(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD13)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD13
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d7=7.0_iwp,d9=9.0_iwp
 REAL(iwp)::d5=5.0_iwp,d8=8.0_iwp,d10=10.0_iwp
 REAL(IWP)::d11=11.0_iwp,d15=15.0_iwp,d19=19.0_iwp,d26=26.0_iwp
 REAL(iwp)::d14=14.0_iwp,d22=22.0_iwp,d38=38.0_iwp
 A=(d5*X73-(two*X57))*(A2-F1)*X51+(two*A2-F1)*X51**2-(X73**2*A2)-         &
   (d10*X73*X57*A2)+two*X57**2*A2
 B=(d5*Y73-(two*Y57))*(A2-F1)*Y51+(two*A2-F1)*Y51**2-(Y73**2*A2)-         &
   (d10*Y73*Y57*A2)+two*Y57**2*A2
 DIAG=(two*(A*G+B*E1))/(d9*(d3*A2**2-F1**2))
 C=X51**2*A2+two*X51**2*F2+d19*X51*X73*A2-(d11*X51*X73*F2)+d14*           &
   X51*X57*A2-(d8*X51*X57*F2)-(d26*X73**2*A2)+d15*X73**2*F2-(d38*X73*     &
   X57*A2)+d22*X73*X57*F2-(d14*X57**2*A2)+d8*X57**2*F2
 D=Y51**2*A2+two*Y51**2*F2+d19*Y51*Y73*A2-(d11*Y51*Y73*F2)+d14*           &
   Y51*Y57*A2-(d8*Y51*Y57*F2)-(d26*Y73**2*A2)+d15*Y73**2*F2-(d38*Y73*     &
   Y57*A2)+d22*Y73*Y57*F2-(d14*Y57**2*A2)+d8*Y57**2*F2
 OPPO=(two*(C*G+D*E1))/(d9*(d3*A2**2-F2**2))
 DD13 = DIAG + OPPO
RETURN
END FUNCTION DAD13

SUBROUTINE FAM14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(8,3)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(12,7)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,11)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,4)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,8)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(7,4)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,3)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,12)=DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM14

FUNCTION DAD14(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD14)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD14
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp,d7=7.0_iwp
 REAL(iwp)::d9=9.0_iwp,d11=11.0_iwp,d15=15.0_iwp,d19=19.0_iwp
 REAL(iwp)::d22=22.0_iwp,d26=26.0_iwp,d38=38.0_iwp,d40=40.0_iwp
 REAL(iwp)::d6=6.0_iwp
 A=two*(d5*(A2-F1)*X51-(X73*A2))*Y73-(d4*((A2-F1)*X51+d5*X73*A2-(X57*A2))*&
   Y57)+two*(two*A2-F1)*Y51*X51-((d11*F1+F2)*A2)+two*A2**2+d3*F1**2+F1*F2
 B=-(d6*(d3*A2**2-F1**2))
 DIAG=-((G+E2)*A+B*E2)/(d9*(d3*A2**2-F1**2))
 C=d4*((X51-X57)*(d7*A2-(d4*F2))-((d19*A2-(d11*F2))*X73))*Y57-            &
   (two*((d26*A2-(d15*F2))*X73-((d19*A2-(d11*F2))*X51))*Y73)+two*(A2+two* &
   F2)*Y51*X51-((d19*F1+d11*F2)*A2)+D40*A2**2+d11*F1*F2-(d7*F2**2)
 D=d6*(d3*A2**2-F2**2)
 OPPO=-((G+E2)*C+D*E2)/(d9*(d3*A2**2-F2**2))
 DD14=DIAG+OPPO
RETURN
END FUNCTION DAD14

SUBROUTINE FAM15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(11,3)=DAD15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,7)=DAD15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS3(X51,X57,X73,Y51,Y57,Y73)
 KM(16,8)=DAD15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.TRUE.)
 KM(12,4)=DAD15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM15

FUNCTION DAD15(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD15)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,B,C,D,DD15
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp,d7=7.0_iwp
 REAL(iwp)::d9=9.0_iwp,d11=11.0_iwp,d15=15.0_iwp,d19=19.0_iwp
 REAL(iwp)::d22=22.0_iwp,d26=26.0_iwp,d38=38.0_iwp,d40=40.0_iwp
 A=two*((d5*A2-(two*F1))*X73+(A2-F1)*X57)*X51+(d4*A2+d3*F1)*X73**2-       &
   ((two*A2-F1)*X51**2)-(two*(A2-(d3*F1))*X73*X57)-(two*X57**2*A2)
 B=-(two*((d5*A2-(two*F1))*Y73+(A2-F1)*Y57)*Y51+(d4*A2+d3*F1)*Y73**2-     &
   ((two*A2-F1)*Y51**2)-(two*(A2-(d3*F1))*Y73*Y57)-(two*Y57**2*A2))
 DIAG=(two*(A*G-(B*E1)))/(d9*(d3*A2**2-F1**2))
 C=two*((d5*A2-(two*F2))*X73+(A2-(d3*F2))*X57)*X51+(d4*A2+d3*F2)*X51**2-  &
   ((two*A2-F2)*X73**2)-(two*(A2-F2)*X73*X57)-(two*X57**2*A2)
 D=two*((d5*A2-(two*F2))*Y73+(A2-(d3*F2))*Y57)*Y51+(d4*A2+d3*F2)*Y51**2-  &
   ((two*A2-F2)*Y73**2)-(two*(A2-F2)*Y73*Y57)-(two*Y57**2*A2)
 OPPO=(two*(C*G+D*E1))/(d9*(d3*A2**2-F2**2))
 DD15=DIAG+OPPO
RETURN
END FUNCTION DAD15

SUBROUTINE FAM16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G,KM)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp),INTENT(OUT)::KM(16,16)
 KM(12,3)=DAD16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(16,7)=DAD16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(11,4)=DAD16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
 CALL TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,.FALSE.)
 KM(15,8)=DAD16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)
RETURN
END SUBROUTINE FAM16

FUNCTION DAD16(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G)RESULT(DD16)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN):: X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,S1,E1,E2,G
 REAL(iwp)::DIAG,OPPO,A,C,DD16
 REAL(iwp)::two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp,d7=7.0_iwp
 REAL(iwp)::d9=9.0_iwp,d10=10.0_iwp,d15=15.0_iwp,d19=19.0_iwp
 REAL(iwp)::d22=22.0_iwp,d26=26.0_iwp,d38=38.0_iwp,d40=40.0_iwp
 A=d4*((A2-F1)*X51-((A2-(d3*F1))*X73)-(X57*A2))*Y57+two*(d10*X51          &
       *A2-(d4*X51*F1)+d4*X73*A2+d3*X73*F1)*Y73-(two*(two*A2-F1)*Y51*     &
        X51)-((d7*F1-F2)*A2)+d10*A2**2+d3*F1**2-(F1*F2)
 DIAG=-((G+E2)*A)/(d9*(d3*A2**2-F1**2))
 C=d4*((A2-F2)*X73-((A2-(d3*F2))*X51)+X57*A2)*Y57-(two*(d10*X51*A2-       &
  (d4*X51*F2)-(two*X73*A2)+X73*F2)*Y73)-(two*(d4*A2+d3*F2)*Y51*X51)+      &
  (F1+F2)*A2-(d10*A2**2)-(F1*F2)+d3*F2**2
 OPPO=((G+E2)*C)/(d9*(d3*A2**2-F2**2))
 DD16=DIAG+OPPO
RETURN
END FUNCTION DAD16

SUBROUTINE TRANS1(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,SWITCH)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 LOGICAL,INTENT(IN)::SWITCH
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2
 REAL(iwp),INTENT(OUT)::F1,F2
 REAL(iwp)::X51OLD,X57OLD,X73OLD,Y51OLD,Y57OLD,Y73OLD,two=2.0_iwp
 X51OLD=X51
 X57OLD=X57
 X73OLD=X73
 Y51OLD=Y51
 Y57OLD=Y57
 Y73OLD=Y73
 X51=X73OLD
 X57=X51OLD-X57OLD
 X73=-X51OLD
 Y51=Y73OLD
 Y57=Y51OLD-Y57OLD
 Y73=-Y51OLD
 IF(SWITCH)THEN
   F1=two*X73*Y57-two*Y73*X57+A2
   F2=two*X51*Y57-two*Y51*X57+A2
 ELSE
   F1=two*Y73*X57-two*X73*Y57+A2
   F2=two*Y51*X57-two*X51*Y57+A2
 ENDIF
RETURN
END SUBROUTINE TRANS1

SUBROUTINE TRANS2(X51,X57,X73,Y51,Y57,Y73,A2,F1,F2,SWITCH)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 LOGICAL,INTENT(IN)::SWITCH
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73,A2
 REAL(iwp),INTENT(OUT)::F1,F2
 REAL::X51OLD,X57OLD,X73OLD,Y51OLD,Y57OLD,Y73OLD,two=2.0_iwp
 X51OLD=X51
 X57OLD=X57
 X73OLD=X73
 Y51OLD=Y51
 Y57OLD=Y57
 Y73OLD=Y73
 X51=Y51OLD
 X57=Y57OLD+Y73OLD
 X73=-Y73OLD
 Y51=X51OLD
 Y57=X57OLD+X73OLD
 Y73=-X73OLD
 IF(SWITCH)THEN
   F1=two*X73*Y57-two*Y73*X57+A2
   F2=two*X51*Y57-two*Y51*X57+A2
 ELSE
   F1=two*Y73*X57-two*X73*Y57+A2
   F2=two*Y51*X57-two*X51*Y57+A2
 ENDIF
RETURN
END SUBROUTINE TRANS2

SUBROUTINE TRANS3(X51,X57,X73,Y51,Y57,Y73)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::X51,X57,X73,Y51,Y57,Y73
 REAL(iwp)::X51OLD,X57OLD,X73OLD,Y51OLD,Y57OLD,Y73OLD
 X51OLD=X51
 X57OLD=X57
 X73OLD=X73
 Y51OLD=Y51
 Y57OLD=Y57
 Y73OLD=Y73
 X51=Y51OLD
 X57=Y57OLD
 X73=Y73OLD
 Y51=X51OLD
 Y57=X57OLD
 Y73=X73OLD
RETURN
END SUBROUTINE TRANS3

END PROGRAM test_8


