module mod_solver
    use mod_variables
    use mod_rhs
    use mod_matrix
    implicit none
    
contains

subroutine solver(t2)
    implicit none
    integer  :: j, m
    real(8)  :: t2

    do j = jstart, jfinal
        call create_matrix(j)
        do m = 0, min(2,j), 2
            call generate_rhs(t2, j, m)
            call banbks(S_U, smtot, sph_lo_diag, sph_hi_diag, smtot, sph_band, S_L, sph_lo_diag, S_I, RHS%spheriodal)
            call banbks(T_U, tmtot, tor_lo_diag, tor_hi_diag, tmtot, tor_band, T_L, tor_lo_diag, T_I, RHS%torodial  )
            call separate_variables(j,m)
        end do
    end do
end subroutine solver

subroutine separate_variables(j,m)
    implicit none
    integer :: il, ir, jm1, jp1, j00, jms, j, m
    integer :: nn, smm, tmm
    
    jm1 = jml2(j, m, j-1)
    j00 = jml2(j, m, j  )
    jp1 = jml2(j, m, j+1)
    jms =  jm2(j, m)


! !==========VELOCITY==========
    smm = 0; tmm = 0
    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 1,nn
            layer(il)%vel(ir, jm1) = RHS%spheriodal(smm+7*(ir-1)+1)
            layer(il)%vel(ir, jp1) = RHS%spheriodal(smm+7*(ir-1)+2)
            layer(il)%vel(ir, j00) = RHS%torodial  (tmm+3*(ir-1)+1)
        end do
        smm = sum(layer(1:il)%sm); tmm = sum(layer(1:il)%tm)
        layer(il)%vel(nn+1, jm1) = RHS%spheriodal(smm-1)
        layer(il)%vel(nn+1, jp1) = RHS%spheriodal(smm  )
        layer(il)%vel(nn+1, j00) = RHS%torodial  (tmm)
    end do
! !==========DENSITY==========
    smm = 0
    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 1, nn
            layer(il)%rhoe(ir, jms) = RHS%spheriodal(smm + 7*(ir-1)+3)
        end do
        smm = sum(layer(1:il)%sm)
    end do
! !==========STRESS=========
    smm = 0; tmm = 0
    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 1, nn
            layer(il)%str(ir, 6*(jms-1)+1) = RHS%spheriodal(smm + 7*(ir-1)+4)
            layer(il)%str(ir, 6*(jms-1)+2) = RHS%spheriodal(smm + 7*(ir-1)+5)
            layer(il)%str(ir, 6*(jms-1)+4) = RHS%spheriodal(smm + 7*(ir-1)+6)
            layer(il)%str(ir, 6*(jms-1)+6) = RHS%spheriodal(smm + 7*(ir-1)+7)
            

            layer(il)%str(ir, 6*(jms-1)+3) = RHS%torodial  (tmm + 3*(ir-1)+2)
            layer(il)%str(ir, 6*(jms-1)+5) = RHS%torodial  (tmm + 3*(ir-1)+3)
        end do
        smm = sum(layer(1:il)%sm); tmm = sum(layer(1:il)%tm)
    end do
end subroutine separate_variables
    
end module mod_solver