! Author: jnshuti
! Created: 2025-10-03 14:09:49
! Last modified: 2025-10-03 14:10:10 jnshuti
! model for L-matrix of NO3 radical 

module diab_mod
  use accuracy_constants, only: dp, idp
  use dim_parameter, only: ndiab, nstat, ntot,qn,pst
  implicit none
  private
  public :: diab 
  
  contains 
    subroutine diab(lx,ly,lz,n,x1,x2,p)
        implicit none 
        real(dp), intent(out),dimension(ndiab,ndiab):: lx,ly,lz 
        real(dp), intent(in), dimension(qn):: x1,x2 
        real(dp), intent(in),dimension(:):: p 
        integer(idp),intent(in):: n

        call Lx_diab(lx,x1,x2,p)
        call Ly_diab(ly,x1,x2,p)
        call Lz_diab(lz,x1,x2,p)

    end subroutine diab 


    subroutine Lx_diab(E,q,t,p)
        implicit none 
        real(dp),dimension(ndiab,ndiab), intent(out):: E 
        real(dp),dimension(:),intent(in):: q,t
        real(dp),dimension(:),intent(in):: p 
        real(dp):: xs,ys,xb,yb,a,b 
        real(dp):: v3_vec(8), v2(6)
        integer(idp):: i,j,id



        ! check the dimension of the matrix 
        if (size(E,1) .ne. ndiab) then 
            write(*,*) " Error in Lx_diab: wrong dimension of L matrix ", size(E,1)
            stop
        endif

        ! rewrite the coordinate array q into symmetry adapted coordinates
        call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
        v2(1)=xs**2-ys**2
        v2(2)=xb**2-yb**2
        v2(3)=xs*xb-ys*yb
        v2(4)=2*xs*ys 
        v2(5)=2*xb*yb 
        v2(6)=xs*yb+xb*ys

        e = 0.0_dp 
        id = 1
        e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 param
        id=id+1   ! 2
        e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb  ! 2 p
        e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
        id =id+1 ! 3
        e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb   ! 2 p
        e(5,5)=e(5,5)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb

        
        id=id+1 ! 4
        ! order 2
        e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 5 p
                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
        id =id+1 ! 5
        e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
        e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
        id =id+1  ! 6
        e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
        e(5,5)=e(5,5)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
        ! W and Z term of E1
        ! order 0
        id=id+1  ! 7
        e(2,2)=e(2,2)+p(pst(1,id))
        e(3,3)=e(3,3)-p(pst(1,id))
        !e(2,3)=e(2,3)

        ! order 1
        id=id+1  ! 8 ! 2 param
        e(2,2)=e(2,2)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb
        e(3,3)=e(3,3)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
        e(2,3)=e(2,3)- p(pst(1,id))*ys -p(pst(1,id)+1)*yb 
        ! order 2
        id=id+1 ! 9  ! 3p 
        do i=1,3
        e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i)
        e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i)
        e(2,3)=e(2,3)+ p(pst(1,id)+(i-1))*v2(i+3)
        enddo 
        ! order 3
        

        ! try the testing of higher order terms 
        !e(2,3)=e(2,3)- p(pst(1,id))*ys*ss +p(pst(1,id)+1)*ss*2*xs*ys 
        



        ! W and Z for E2
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        id=id+1  ! 10
        e(4,4)=e(4,4)+p(pst(1,id))
        e(5,5)=e(5,5)-p(pst(1,id))
        e(4,5)=e(4,5)

        ! order 1
        id=id+1  !112 param 15
        e(4,4)=e(4,4)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb  
        e(5,5)=e(5,5)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
        e(4,5)=e(4,5)- p(pst(1,id))*ys-p(pst(1,id)+1)*yb 
        ! order 2
        id=id+1 ! 12 ! 3p 
        do i=1,3
        e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i)
        e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i)
        e(4,5)=e(4,5)+ p(pst(1,id)+(i-1))*v2(i+3)
        enddo 

        ! make the dipole E = b* E

        e = b * e 

        ! E1 X E2 
        ! WW and ZZ 
        

        id =id+1 ! 13
        e(2,4)=e(2,4)+p(pst(1,id))*b
        e(3,5)=e(3,5)-p(pst(1,id))*b

        ! ORDER 1
        id=id+1 ! 14  ! 6 parama
        e(2,4)=e(2,4)+b*((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*xb)
        e(3,5)=e(3,5)+b*((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
        e(2,5)=e(2,5)+b*((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
        e(3,4)=e(3,4)+b*((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
        ! order 2
        id=id+1   ! 15

        do i=1,3 !  param
        e(2,4)=e(2,4)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
        e(3,5)=e(3,5)+b*(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
        e(2,5)=e(2,5)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
        e(3,4)=e(3,4)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i+3)
        enddo 
        
        ! pseudo A2 & E1
        ! ##################################################
        !###################################################
        ! order 0
        id=id+1 ! 1 param ! 16

        e(1,3)=e(1,3)+b*(p(pst(1,id)))
        
        ! order 1 
        id = id +1 ! 17 
        e(1,2)=e(1,2)-b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
        e(1,3)=e(1,3)+b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
        
        ! order 2 
        id=id+1 ! 18

        e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
        e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) &
                        +p(pst(1,id)+2)*(xs*xb-ys*yb))

        ! COUPLING OF A2 WITH E2

        !##########################################################################################################
                    
        ! order 0
        id =id+1 !19
        e(1,5)=e(1,5)+p(pst(1,id))
        
        ! order 1 
        id = id +1 ! 20 
        e(1,4)=e(1,4)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
        e(1,5)=e(1,5)+(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
        
        ! order 2 
        id=id+1 ! 21

        e(1,4)=e(1,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)
        e(1,5)=e(1,5)+p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) &
                        +p(pst(1,id)+2)*(xs*xb-ys*yb)

        e(1,4:5) = b* e(1,4:5)
        

        call copy_2_lower_triangle(e)
    end subroutine Lx_diab

        ! Ly matrix
    subroutine Ly_diab(e,q,t,p)
        implicit none 
        real(dp),dimension(ndiab,ndiab), intent(out):: e
        real(dp),dimension(:),intent(in):: q,t
        real(dp),dimension(:),intent(in):: p 
        real(dp):: xs,ys,xb,yb,a,b 
        real(dp):: v2(6)
        integer(idp):: i,j,id

        ! check the dimension of the matrix 
        if (size(e,1) .ne. ndiab) then 
            write(*,*) " Error in Ly_diab: wrong dimension of L matrix ", size(e,1)
            stop
        endif

        ! rewrite the coordinate array q into symmetry adapted coordinates
        call rewrite_coord(q,a,xs,ys,xb,yb,b,1) 
        v2(1)=xs**2-ys**2
        v2(2)=xb**2-yb**2
        v2(3)=xs*xb-ys*yb
        v2(4)=2*xs*ys 
        v2(5)=2*xb*yb 
        v2(6)=xs*yb+xb*ys

        e = 0.0_dp 
        
        ! V-term 

        id=1  ! 1
        ! order 1
        e(1,1)=e(1,1)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
        id=id+1  ! 2
        e(2,2)=e(2,2)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
        e(3,3)=e(3,3)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
        id =id+1  ! 3
        e(4,4)=e(4,4)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
        e(5,5)=e(5,5)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb

        id=id+1  ! 4b*(
        e(1,1)=e(1,1)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
        id =id+1 ! 5
        e(2,2)=e(2,2)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
        e(3,3)=e(3,3)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
        id=id+1 ! 6
        e(4,4)=e(4,4)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
        e(5,5)=e(5,5)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
         
        
        !  W and Z of E1
        ! order 0
        id=id+1 ! 7
        e(2,3)=e(2,3)+p(pst(1,id))
        ! order 1
        id=id+1 ! 8
        e(2,2)=e(2,2)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb 
        e(3,3)=e(3,3)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb 
        e(2,3)=e(2,3)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
        ! order 2
        id=id+1 ! 9
        do i=1,3
        e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i+3)
        e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i+3)
        e(2,3)=e(2,3)-p(pst(1,id)+(i-1))*v2(i) 
        enddo 
        

        !! W and Z of E2
        !!!!!!!!!!!!!!!!!!!!!!!!!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!

        ! order 0
        id=id+1 ! 10
        e(4,5)=e(4,5)+p(pst(1,id))
        ! order 1
        id=id+1 ! 11
        e(4,4)=e(4,4)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb 
        e(5,5)=e(5,5)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb 
        e(4,5)=e(4,5)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
        ! order 2
        id=id+1 ! 12
        do i=1,3
        e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i+3)
        e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i+3)
        e(4,5)=e(4,5)-p(pst(1,id)+(i-1))*v2(i) 
        enddo 

        ! PSEUDO JAHN-TELLER E1 AND E2

        e = b* e 

        !ORDER 0 
        id=id+1 ! 13

        e(2,5)=e(2,5)+p(pst(1,id))
        e(3,4)=e(3,4)+p(pst(1,id))
        ! order 1
        
        id=id+1 ! 14
        e(2,4)=e(2,4)+((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
        e(3,5)=e(3,5)+((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*yb)
        e(2,5)=e(2,5)+((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
        e(3,4)=e(3,4)+((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(+p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*xb)

        ! order 2
        id=id+1 ! 15 

        e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
        e(3,5)=e(3,5)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
        e(2,5)=e(2,5)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i)
        e(3,4)=e(3,4)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)

        ! no order 3
        !!!!!!!!!!!!!!!!

        ! the coupling A2 & E1
        ! #####################
        ! order 0

        id=id+1 ! 16 
        e(1,2)=e(1,2)+(p(pst(1,id)))
        ! order 1

        id=id+1 ! 17 
        e(1,2)=e(1,2)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
        e(1,3)=e(1,3)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)

        ! order 2
        id=id+1 !18 
        e(1,2)=e(1,2)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                        +p(pst(1,id)+2)*(xs*xb-ys*yb))
        e(1,3)=e(1,3)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ &
                        +p(pst(1,id)+2)*(xs*yb+xb*ys))

        ! COUPLING OF A2 WITH E2
        !#######################################################################################
        !###############################################################################
        ! order 0 

        id = id+1 !19
        e(1,4)=e(1,4)+p(pst(1,id))
        ! order 1

        id=id+1 ! 20
        e(1,4)=e(1,4)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
        e(1,5)=e(1,5)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)

        ! order 2
        id=id+1 ! 21 
        e(1,4)=e(1,4)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                        +p(pst(1,id)+2)*(xs*xb-ys*yb))
        e(1,5)=e(1,5)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ &
                        p(pst(1,id)+2)*(xs*yb+xb*ys))

        !write(*,*)'idy=',id
        e(1:4,5) = b * e(1:4,5)


        call copy_2_lower_triangle(e)

    end subroutine Ly_diab
         ! Lz matrix
    subroutine Lz_diab(e,q,t,p)
        implicit none 
        real(dp),dimension(ndiab,ndiab), intent(out):: e 
        real(dp),dimension(:),intent(in):: q,t
        real(dp),dimension(:),intent(in):: p 
        real(dp):: xs,ys,xb,yb,a,b 
        real(dp):: v2(6)
        integer(idp):: i,j,id 

        ! check the dimension of the matrix 
        if (size(e,1) .ne. ndiab) then 
            write(*,*) " Error in Lz_diab: wrong dimension of e matrix ", size(e,1)
            stop
        endif
        call rewrite_coord(q,a,xs,ys,xb,yb,b,1)


        e = 0.0_dp
        ! id for lz 
        id =  22 ! has to be 
        ! the diagonal terms 

        ! the v-term is 0th order and 3rd order.
        ! There is no zeroth order for diagonal

        ! w and z of E''
        ! order 1 
        id = id  
        e(2,2) = e(2,2) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb  
        e(3,3) = e(3,3) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
        e(2,3) = e(2,3) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb

        ! order 2
        id = id +1 
        do i =1,3 
            e(2,2) = e(2,2) + p(pst(1,id)+(i-1))*v2(i+3) 
            e(3,3) = e(3,3) - p(pst(1,id)+(i-1))*v2(i+3)
            e(2,3) = e(2,3) + p(pst(1,id)+(i-1))*v2(i)
        enddo 

        ! W and Z of E'
        ! order 1 
        
        id = id +1 
        e(4,4) = e(4,4) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb  
        e(5,5) = e(5,5) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
        e(4,5) = e(4,5) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb

        ! order 2
        id = id +1 
         do i =1,3 
            e(4,4) = e(4,4) + p(pst(1,id)+(i-1))*v2(i+3) 
            e(5,5) = e(5,5) - p(pst(1,id)+(i-1))*v2(i+3)
            e(4,5) = e(4,5) + p(pst(1,id)+(i-1))*v2(i)
        enddo 
        ! the coupling 
        ! Pseudo of E' and E'' 
        ! it must have odd power of b

        id = id +1 
        ! order 0 
        e(2,4) = e(2,4) 
        e(3,5) = e(3,5) 
        e(2,5) = e(2,5) + b*(p(pst(1,id)))
        e(3,4) = e(3,4) - b*(p(pst(1,id)))

        ! order 1 
        id = id +1 
        e(2,4) = e(2,4) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) 
        e(3,5) = e(3,5) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
        e(2,5) = e(2,5) - b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
        e(3,4) = e(3,4) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
        ! order 2 
        id = id +1 
        do i=1,3 
            e(2,4) = e(2,4) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
            e(3,5) = e(3,5) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
            e(2,5) = e(2,5) + b*(p(pst(1,id)+(i-1)))*v2(i)
            e(3,4) = e(3,4) - b*(p(pst(1,id)+(i-1)))*v2(i)
        enddo 
        ! no third order 

        ! the coupling between A2' and E''
        ! order 1 
        id = id +1 
        e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)*xb))
        e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)*yb))
        
        id = id +1 
        ! order 2 
        do i=1,3
            e(1,2) = e(1,2) + b*(p(pst(1,id)+(i-1)))*v2(i)
            e(1,3) = e(1,3) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
        enddo 

        ! the coupling of A2' and E'

        ! order 1 
        id = id +1 
        e(1,2) = e(1,2) + (p(pst(1,id))*xs + p(pst(1,id)*xb))
        e(1,3) = e(1,3) - (p(pst(1,id))*ys + p(pst(1,id)*yb))
        
        id = id +1 
        ! order 2 
        do i=1,3
            e(1,2) = e(1,2) + (p(pst(1,id)+(i-1)))*v2(i)
            e(1,3) = e(1,3) + (p(pst(1,id)+(i-1)))*v2(i+3)
        enddo 
        call copy_2_lower_triangle(e)

 

    end subroutine Lz_diab

    subroutine rewrite_coord(q,a,xs,ys,xb,yb,b,start)
        implicit none 
        real(dp),dimension(:),intent(in):: q
        real(dp),intent(out):: xs,ys,xb,yb,a,b
        integer(idp),intent(in):: start
        integer(idp):: i,j

        a= q(start)
        xs = q(start+1)
        ys = q(start+2)
        xb = q(start+3)
        yb = q(start+4)
        b = q(start+5)
    end subroutine rewrite_coord

    subroutine copy_2_lower_triangle(mat)
        real(dp), intent(inout) :: mat(:, :)
        integer :: m, n
    !     write lower triangle of matrix symmetrical
        do n=1,size(mat,1)
            do m=n,size(mat,1)
                mat(m,n)=mat(n,m)
            enddo 
        enddo
   end subroutine copy_2_lower_triangle


end module diab_mod