module mod_norms
    use mod_parameters
    use mod_variables
    use mod_write_f
    implicit none
    
contains
    

subroutine compute_dissipation
    implicit none
    integer :: i, k, j, m, jmi
    real(8) :: sjm2_0, sjm1_0, sjm0_0, sjp1_0, sjp2_0, dr, integral(1:numoflayers)
    real(8) :: sjm2_1, sjm1_1, sjm0_1, sjp1_1, sjp2_1



    layer(:)%dissipation = 0.0d0
    do i = 1, numoflayers
        layer(i)%spectrum = 0.0d0
    end do
    
    do j = jstart, jfinal
        do m = 0, min(j,2), 2
            jmi = jm2(j,m) - 1
            do i = 1, numoflayers
                integral(i)=0.0d0
                do k = 1, layer(i)%n - 1
                    dr   = layer(i)%r_i(k+1) - layer(i)%r_i(k)
                    sjm2_0 =  (layer(i)%r_i(k  )**2)*real(layer(i)%str(k  , 6*jmi+2)*dconjg(layer(i)%str(k  , 6*jmi+2))) / (2.0d0*layer(i)%eta_i(k))
                    sjm1_0 =  (layer(i)%r_i(k  )**2)*real(layer(i)%str(k  , 6*jmi+3)*dconjg(layer(i)%str(k  , 6*jmi+3))) / (2.0d0*layer(i)%eta_i(k))
                    sjm0_0 =  (layer(i)%r_i(k  )**2)*real(layer(i)%str(k  , 6*jmi+4)*dconjg(layer(i)%str(k  , 6*jmi+4))) / (2.0d0*layer(i)%eta_i(k))
                    sjp1_0 =  (layer(i)%r_i(k  )**2)*real(layer(i)%str(k  , 6*jmi+5)*dconjg(layer(i)%str(k  , 6*jmi+5))) / (2.0d0*layer(i)%eta_i(k))
                    sjp2_0 =  (layer(i)%r_i(k  )**2)*real(layer(i)%str(k  , 6*jmi+6)*dconjg(layer(i)%str(k  , 6*jmi+6))) / (2.0d0*layer(i)%eta_i(k))

                    sjm2_1 =  (layer(i)%r_i(k+1)**2)*real(layer(i)%str(k+1, 6*jmi+2)*dconjg(layer(i)%str(k+1, 6*jmi+2))) / (2.0d0*layer(i)%eta_i(k+1))
                    sjm1_1 =  (layer(i)%r_i(k+1)**2)*real(layer(i)%str(k+1, 6*jmi+3)*dconjg(layer(i)%str(k+1, 6*jmi+3))) / (2.0d0*layer(i)%eta_i(k+1))
                    sjm0_1 =  (layer(i)%r_i(k+1)**2)*real(layer(i)%str(k+1, 6*jmi+4)*dconjg(layer(i)%str(k+1, 6*jmi+4))) / (2.0d0*layer(i)%eta_i(k+1))
                    sjp1_1 =  (layer(i)%r_i(k+1)**2)*real(layer(i)%str(k+1, 6*jmi+5)*dconjg(layer(i)%str(k+1, 6*jmi+5))) / (2.0d0*layer(i)%eta_i(k+1))
                    sjp2_1 =  (layer(i)%r_i(k+1)**2)*real(layer(i)%str(k+1, 6*jmi+6)*dconjg(layer(i)%str(k+1, 6*jmi+6))) / (2.0d0*layer(i)%eta_i(k+1))
                    integral(i) = integral(i) + 0.5d0 * dr * (sjm2_0 + sjm1_0 + sjm0_0 + sjp1_0 + sjp2_0) &
                                              + 0.5d0 * dr * (sjm2_1 + sjm1_1 + sjm0_1 + sjp1_1 + sjp2_1) 
                end do
                integral(i)=integral(i)* dble(2-dirac(m,0))
                layer(i)%dissipation=layer(i)%dissipation       + integral(i) 
                layer(i)%spectrum(j)   = layer(i)%spectrum(j)   + integral(i)
            end do
        end do
    end do
end subroutine compute_dissipation


subroutine velocity_norm
    implicit none
    integer    :: il, ir, j, m, jm1, jm0, jp1, nn
    real(8)    :: integral(1:numoflayers)
    real(8)    :: vvjm1_0, vvj00_0, vvjp1_0
    real(8)    :: vvjm1_1, vvj00_1, vvjp1_1, dr


    do il = 1 , numoflayers
        layer(il)%spectrum_v = zero
    end do
    do j = jstart, jfinal
        do m = 0, min(2, j), 2
            jm1 = jml2(j, m, j-1)
            jm0 = jml2(j, m, j  )
            jp1 = jml2(j, m, j+1)
            do il = 1, numoflayers
                integral(il) = zero
                nn = layer(il)%n
                do ir = 1, nn
                    dr = layer(il)%r_l(ir+1) - layer(il)%r_l(ir)
                    vvjm1_0 = layer(il)%r_l(ir) * real(layer(il)%vel(ir, jm1)*dconjg(layer(il)%vel(ir, jm1)))
                    vvj00_0 = layer(il)%r_l(ir) * real(layer(il)%vel(ir, jm0)*dconjg(layer(il)%vel(ir, jm0)))
                    vvjp1_0 = layer(il)%r_l(ir) * real(layer(il)%vel(ir, jp1)*dconjg(layer(il)%vel(ir, jp1)))

                    vvjm1_1 = layer(il)%r_l(ir+1) * real(layer(il)%vel(ir+1, jm1)*dconjg(layer(il)%vel(ir+1, jm1)))
                    vvj00_1 = layer(il)%r_l(ir+1) * real(layer(il)%vel(ir+1, jm0)*dconjg(layer(il)%vel(ir+1, jm0)))
                    vvjp1_1 = layer(il)%r_l(ir+1) * real(layer(il)%vel(ir+1, jp1)*dconjg(layer(il)%vel(ir+1, jp1)))

                    integral(il) = integral(il) + dr * 0.5d0 * (vvjm1_0 + vvjm1_1 + vvj00_0 + vvj00_1 + vvjp1_0 + vvjp1_1)
                end do
                integral(il)=integral(il)* dble(2-dirac(m,0))
                layer(il)%spectrum_v(j) = layer(il)%spectrum_v(j) + integral(il)

            end do
        end do
    end do
end subroutine velocity_norm


subroutine check_spectrum(it2)
    implicit none
    integer :: i, change_status, it2
    real(8) :: trunc_3

    call write_spectrum(it2)
    change_status = 0
    old_jmax = jmax
    if(old_jmax == Max_jmax) goto 8
    do i = 1, numoflayers
        trunc_3 = log10(maxval(layer(i)%spectrum)) - 3.0d0
        if((jmax.ge.3)) then
            if(maxval(log10(layer(i)%spectrum(jmax-5:jmax))).gt.trunc_3) then
                change_status = 1
                exit
            end if
        end if
    end do
    if(change_status.eq.1)then
        jmax = jmax + 5
        jfinal = jmax
    end if
8 continue
end subroutine check_spectrum

subroutine check_stability(p)
    implicit none
    integer :: p, i
    real(8) :: delta_k22   (1:9)
    
    do i = 1, 10-1
        delta_k22(i) = abs(rk22_p(p-i+1) -rk22_p(p-i))
        delta_k22(i) = delta_k22(i)*100.0d0/rk22_p(p-i)
    end do
    if (maxval(delta_k22) <= 0.1d0) signal=0

end subroutine check_stability

subroutine check_signal
    implicit none
    
    if(signal==1)then
        open(400, file=trim(dir_output)//"/signal")
        read(400, *) signal
        close(400)
    end if
end subroutine check_signal

end module mod_norms