subroutine CoeffExtSplit(dtq,omgxq,omgyq,omgzq,A,u12,u13,u23,l21,l31,l32,Vr,v1,v2,v3)
!----------------------------------------------------------------------!
!  input: dtq, omgxq, omgyq, omgzq                                     !
!                                                                      !
!----------------------------------------------------------------------!

	implicit none
        real(16), intent(in) :: omgxq,omgyq,omgzq
        real(16) :: dtq
        real(8), intent(out) :: l21,l31,l32,u12,u13,u23, v1, v2,v3
        real(8), intent(out) :: A(3,3),Vr(3,3)
        real(16) :: lq21,lq31,lq32,uq12,uq13,uq23, vq1, vq2,vq3
        real(16) :: Aq(3,3),Vrq(3,3),Theta1_xi(3,3), Theta2_xi_mi1(3,3), Theta2_eta(3,3),W1(3,3),W(3,3),I3(3,3)
        real(16) :: m,c,s,H11,H12,H13,H21,H22,H23,H31,H32,H33 
     
        m = qsqrt(omgxq**2.q0+omgyq**2.q0+omgzq**2.q0)
        c = (1.q0-cos(m*dtq))/(m*dtq)**2.q0
        s = sin(m*dtq)/(m*dtq)
 
        vq1 = 1.q0/2.q0*(omgyq**2.q0+omgzq**2.q0)*c/(cos(m*dtq)+omgxq**2.q0*c*dtq**2.q0)
        vq2 = 1.q0/2.q0*(omgxq**2.q0-omgzq**2.q0)*c/(cos(m*dtq)+omgzq**2.q0*c*dtq**2.q0)
        vq3 = -1.q0/2.q0*(omgxq**2.q0+omgyq**2.q0)*c
        lq21 = -(omgzq*s+omgxq*omgyq*c*dtq)/(cos(m*dtq)+omgxq**2.q0*c*dtq**2.q0)
        lq31 = (omgyq*s+omgxq*omgzq*c*dtq)/(cos(m*dtq)+omgzq**2.q0*c*dtq**2.q0)
        lq32 = (-omgxq*s+omgyq*omgzq*c*dtq)/(cos(m*dtq)+omgzq**2.q0*c*dtq**2.q0)
        uq12 = (omgzq*s-omgxq*omgyq*c*dtq)*(cos(m*dtq)+omgxq**2.q0*c*dtq**2.q0)/(cos(m*dtq)+omgzq**2.q0*c*dtq**2.q0)
        uq13 = -(omgyq*s+omgxq*omgzq*c*dtq)*(cos(m*dtq)+omgzq**2.q0*c*dtq**2.q0)
        uq23 = (omgxq*s+omgyq*omgzq*c*dtq)*(cos(m*dtq)+omgzq**2.q0*c*dtq**2.q0)/(cos(m*dtq)+omgxq**2.q0*c*dtq**2.q0)

                   I3 = reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
            Theta1_xi = reshape((/1.q0, 0.q0, 0.q0, dtq*uq12, 1.q0, 0.q0, dtq**2.q0*uq12*uq23+dtq*uq13, dtq*uq23, 1.q0 /),(/3,3/))       
        Theta2_xi_mi1 = reshape((/1.q0, -dtq*uq12, -dtq*uq13, 0.q0,1.q0,-dtq*uq23, 0.q0, 0.q0,1.q0 /),(/3,3/))
           Theta2_eta = reshape((/ 1.q0, dtq*lq21, dtq*lq31, 0.q0, 1.q0, dtq*lq32, 0.q0, 0.q0,1.q0/),(/3,3/))
            
        H11 = cos(m*dtq)+omgxq**2.q0*(1.q0-cos(m*dtq))/m**2.q0
        H12 = omgzq*sin(m*dtq)/m+omgxq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H13 = -omgyq*sin(m*dtq)/m+ omgxq*omgzq*(1.q0-cos(m*dtq))/m**2.q0
        H21 = -omgzq*sin(m*dtq)/m+ omgxq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H22 = cos(m*dtq)+omgyq**2.q0*(1.q0-cos(m*dtq))/m**2.q0
        H23 = omgxq*sin(m*dtq)/m+omgzq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H31 = omgyq*sin(m*dtq)/m+ omgxq*omgzq*(1.q0-cos(m*dtq))/m**2.q0
        H32 = -omgxq*sin(m*dtq)/m+omgzq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H33 = cos(m*dtq)+omgzq**2.q0*(1.q0-cos(m*dtq))/m**2.q0
          W = reshape((/H11, H21,H31,H12,H22,H32,H13,H23,H33/),(/3,3/))
        
        dtq = -dtq
        H11 = cos(m*dtq)+omgxq**2.q0*(1.q0-cos(m*dtq))/m**2.q0
        H12 = omgzq*sin(m*dtq)/m+omgxq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H13 = -omgyq*sin(m*dtq)/m+ omgxq*omgzq*(1.q0-cos(m*dtq))/m**2.q0
        H21 = -omgzq*sin(m*dtq)/m+ omgxq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H22 = cos(m*dtq)+omgyq**2.q0*(1.q0-cos(m*dtq))/m**2.q0
        H23 = omgxq*sin(m*dtq)/m+omgzq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H31 = omgyq*sin(m*dtq)/m+ omgxq*omgzq*(1.q0-cos(m*dtq))/m**2.q0
        H32 =-omgxq*sin(m*dtq)/m+omgzq*omgyq*(1.q0-cos(m*dtq))/m**2.q0
        H33 = cos(m*dtq)+omgzq**2.q0*(1.q0-cos(m*dtq))/m**2.q0
         W1 = reshape((/H11, H21,H31,H12,H22,H32,H13,H23,H33/),(/3,3/))

        dtq = -dtq
        Vrq = 1.q0/2.q0/dtq**2.q0*(matmul(matmul(W1,Theta1_xi),Theta2_eta)-I3)
         Aq = 1.q0/2.q0*matmul(matmul(Transpose(Theta2_xi_mi1),W),Transpose(Theta2_eta))
        l21 = real(lq21, 8)
        l31 = real(lq31, 8)
        l32 = real(lq32, 8)
        u12 = real(uq12, 8)
        u13 = real(uq13, 8)
        u23 = real(uq23, 8)
        u12 = real(uq12, 8)
         v1 = real(vq1, 8)
         v2 = real(vq2, 8)
         v3 = real(vq3, 8) 
         Vr = real(Vrq,8)
          A = real(Aq,8)

        end
