module mod_time_integration
    use mod_parameters
    use mod_variables
    use mod_allocation
    use mod_displacement
    use mod_matrix
    use mod_solver
    use mod_potential
    use mod_norms
    use mod_coriolis
    use mod_write_f
    use mod_math
    implicit none
    
contains
    
subroutine time_integration
    implicit none
    integer :: it, ip, iter, fiter
    real(8) :: delta_k20_1, delta_k20_0

    t = time_step; fiter = 1
    it = 1; ip = 1; pc = 1
    call open_files_time_series
    print*, '#Time loop has started...'
    do
        call solver(t)
        if(SWITCH_ocean_status == 'v') then
            if((SWITCH_coriolis==1).and.(it>num_time_steps*10+5)) then
                call assign_coriolis(it)
            end if
        end if
        call displacement
        call assign_pvel_boundary
        if(pc>0) then
            call induced_potential
            if((SWITCH_self_grav==1)) call self_gravity
            weight_sg = 1.0d0 !min(1.0d0, 0.1d0 + (it-6*num_time_steps)*0.01d0)
            fiter=1
        end if
    
        if(mod(it, checkpoint)==0) call write_checkpoint(t, it)
        
        call love_numbers
        call compute_dissipation
        call velocity_norm
        call peaks(t, it)
        call assign_love_numbers(ip)

        call write_time_series(t)
        if(mod(it, num_time_steps) == 0) then
            call check_signal
            if(pc>10) call check_stability(pc-1)
            call love_numbers_per(pc)
            ip = 0; pc = pc + 1;
        end if
        if(mod(it, 5)==0)then
            if(status == 1) call check_spectrum(it)
        end if
        if(pc >= start_output) then
            if(mod(it, write_output)==0) then
                if(WRITE_density  == 1) call write_density_f(t, it)
                if(WRITE_velocity == 1) call write_velocity_f(t, it)
                if(WRITE_stress   == 1) call write_stress_f(t, it)
            end if
        end if
        if((pc==nint(final_period)+1).or.(signal == 0)) exit
        if(it > 250) call determine_dt
        t = t + time_step; it = it + 1; ip = ip + 1
        if((SWITCH_coriolis == 1).and.(it > num_time_steps*0)) then
            if(status == 0) call convert2jmax
        end if
    end do
    call close_files
end subroutine time_integration

subroutine time_integration_rerun
    implicit none
    integer :: it, ip, pc_start

    ! t = t + time_step;
    it = start_simulation+1; ip = 1; 
    t = dble(it) * time_step
    pc = nint(t/rotation_period)+1

    pc_start = pc
    !!for the rerun only:
    call radial_displacement
    call induced_potential !! should be followed by the computation of the Coriolis term from vel and pvel
    call assign_coriolis_rerun
    call open_files_time_series_rerun
    print*, '#Time loop has started...'
    do
        call solver(t)
        if((SWITCH_coriolis==1)) call assign_coriolis(it)
        call displacement
        
        call induced_potential
        call love_numbers
        call compute_dissipation
        call assign_love_numbers(ip)
        call peaks(t, it)
        call write_time_series_rerun(t)
        if(mod(it, num_time_steps) == 0) then
            call check_signal
            if(pc>10) call check_stability(pc-1)
            if(pc >= pc_start + 1) call love_numbers_per(pc)
            ip = 0; pc = pc + 1;
        end if
        if(mod(it, 200)==0)then
            if(status == 1) call check_spectrum(it)
        end if
        if(pc >= start_output) then
            if(mod(it, write_output)==0) then
                if(WRITE_density  == 1) call write_density_f(t, it)
                if(WRITE_velocity == 1) call write_velocity_f(t, it)
                if(WRITE_stress   == 1) call write_stress_f(t, it)
            end if
        end if
        if(mod(it, checkpoint)==0) call write_checkpoint(t, it)
        if((pc==nint(final_period)+1).or.(signal == 0)) exit
        call assign_pvel
        t = t + time_step; it = it + 1; ip = ip + 1
        if((SWITCH_coriolis == 1).and.(it > num_time_steps*10)) then
            if(status == 0) call convert2jmax
        end if
    end do
    call close_files
    
end subroutine time_integration_rerun


subroutine assign_pvel
    implicit none
    integer :: il
    do il = 1, numoflayers
        layer(il)%pvel(:,:) = layer(il)%vel(:,:)
    end do
end subroutine assign_pvel


subroutine assign_pvel_boundary
    implicit none
    integer :: il, nn

    do il = 1, numoflayers
        nn = layer(il)%n
        layer(il)%pvel(nn  , :) = layer(il)%vel(nn  , :)
        layer(il)%pvel(nn+1, :) = layer(il)%vel(nn+1, :)
    end do
end subroutine assign_pvel_boundary

subroutine assign_love_numbers(ipp)
    implicit none
    integer :: ipp

    gp20 (ipp) = k20
    rgp22(ipp) = real_k22
    igp22(ipp) = imag_k22

    khp(:,ipp)  = kh(:)
    

    ur20 (ipp) = h20
    rur22(ipp) = real_h22
    iur22(ipp) = imag_h22

    ut20 (ipp) = l20
    rut22(ipp) = real_l22
    iut22(ipp) = imag_l22
end subroutine assign_love_numbers


subroutine love_numbers_per(pp)
    implicit none
    integer :: pp, nn, ii, j
    real(8) :: forcing_potential
    
    nn = layer(numoflayers)%n
    forcing_potential = ang_vel * ang_vel * ecc * dsqrt(18.0d0 * pi/10.0d0) * rad(numoflayers)**2
    k20_p(pp) = maxval(gp20)/forcing_potential
    h20_p(pp) = maxval(ur20)/(forcing_potential/layer(numoflayers)%g_i(nn))
    l20_p(pp) = maxval(ut20)/(forcing_potential*dsqrt(30.0d0/5.0d0)/layer(numoflayers)%g_i(nn))

    do j = 1, 3
        kh_per(3*(j-1)+1, pp) = maxval(khp(3*(j-1)+1,:))/forcing_potential !!4,0
    end do
    
    forcing_potential = ang_vel * ang_vel * ecc * dsqrt(27.0d0 * pi/10.0d0) * rad(numoflayers)**2
    rk22_p(pp) = maxval(rgp22)/forcing_potential
    rh22_p(pp) = maxval(rur22)/(forcing_potential/layer(numoflayers)%g_i(nn))
    rl22_p(pp) = maxval(rut22)/(forcing_potential*dsqrt(30.0d0/5.0d0)/layer(numoflayers)%g_i(nn))

    do j = 1, 3
        kh_per(3*(j-1)+2, pp) = maxval(khp(3*(j-1)+2,:))/forcing_potential !!4,2, real
    end do

    forcing_potential = ang_vel * ang_vel * ecc * dsqrt(48.0d0 * pi/10.0d0) * rad(numoflayers)**2
    ik22_p(pp) = maxval(igp22)/forcing_potential
    ih22_p(pp) = maxval(iur22)/(forcing_potential/layer(numoflayers)%g_i(nn))
    il22_p(pp) = maxval(iut22)/(forcing_potential*dsqrt(30.0d0/5.0d0)/layer(numoflayers)%g_i(nn))
    do j = 1,3
        kh_per(3*(j-1)+3, pp) = maxval(khp(3*(j-1)+3,:))/forcing_potential !!4,2, imag
    end do

    dissip_p(pp) = avrg_dissip / dble(num_time_steps)
    ! if(0 == mod(pp, 3)) then
        open(iu(5), file=trim(dir_output)//"/love_numbers_period.dat", status="old", position="append", action="write")
        do ii = pp, pp
            write(iu(5), '(8(F20.5), 3(F20.5), 9(F20.5))') dble(ii), k20_p(ii), rk22_p(ii), ik22_p(ii), h20_p(ii), rh22_p(ii), ih22_p(ii), dissip_p(ii)*W_to_GW, l20_p(ii), rl22_p(ii), il22_p(ii), &
            kh_per(1, ii), kh_per(2,ii), kh_per(3,ii), kh_per(3+1, ii), kh_per(3+2,ii), kh_per(3+3,ii), kh_per(6+1, ii), kh_per(6+2,ii), kh_per(6+3,ii)
        end do
        close(iu(5))
    ! end if
    ! if ((signal==0).and.(0 /= mod(pp,3))) then
    !     open(iu(5), file=trim(dir_output)//"/love_numbers_period.dat", status="old", position="append", action="write")
    !     do ii = pp - mod(pp,3)+1, pp
    !         write(iu(5), '(8(F20.10), 3(F20.10), 9(F20.10))') dble(ii), k20_p(ii), rk22_p(ii), ik22_p(ii), h20_p(ii), rh22_p(ii), ih22_p(ii), dissip_p(ii)*W_to_GW, l20_p(ii), rl22_p(ii), il22_p(ii), &
    !         kh_per(1, ii), kh_per(2,ii), kh_per(3,ii), kh_per(3+1, ii), kh_per(3+2,ii), kh_per(3+3,ii), kh_per(6+1, ii), kh_per(6+2,ii), kh_per(6+3,ii)
    !     end do
    !     close(iu(5))
    ! end if
    avrg_dissip = 0.5d0 * layer(layer_ocean)%dissipation
    
end subroutine love_numbers_per

subroutine assign_coriolis(it2)
    implicit none
    integer :: lo, it2

    lo = layer_ocean !! layer ocean
    if(SWITCH_coriolis == 1)then
        call coriolis
        if(it2.ge.2) then
            layer(lo)%fc = dzero
            layer(lo)%fc  = 1.5d0 * layer(lo)%fc1 - 0.5d0 * layer(lo)%fc0
        end if
        layer(lo)%fc0 = layer(lo)%fc1
    end if
    
end subroutine assign_coriolis

subroutine assign_coriolis_rerun
    implicit none
    integer :: lo
    lo = layer_ocean !! layer ocean
    call coriolis
    call coriolis_pvel
    layer(lo)%fc = dzero
    layer(lo)%fc  = 1.5d0 * layer(lo)%fc1 - 0.5d0 * layer(lo)%fc0
    layer(lo)%fc0 = layer(lo)%fc1
end subroutine assign_coriolis_rerun
end module mod_time_integration