module mod_preprocess
    use mod_parameters
    use mod_variables
    use mod_allocation
    use mod_read_f
    use mod_write_f
    implicit none
    
contains


subroutine prep
    implicit none
    
    call initial_read
    call number_of_layers
    call alloc_rads
    call assign_layers
    call find_ocean_layer
    call complete__
    call SWITCH_imp
    call mid_layers
    call radial_grads
    call hydrostatic_core
    call compute_mass
    call compute_grav_acc
    call compute_MOI
    call make_dir
    call write_summary
    call write_radial_profiles
    call alloc_NM
    call alloc_RHS
    call alloc_vars
    call write_initial_signal
    call read_initial_cond
    call determine_mindr

end subroutine prep

subroutine initial_read
    implicit none
    
    call read_input
    call read_radial_profile
    if(rotation_period<100) rotation_period = rotation_period * day2sec
    if(temp_radius(1)>1d5) call reverse_temps
    if(log10(maxval(temp_shear_m(N_initial-5:N_initial))) < 6) temp_shear_m = temp_shear_m * GPa_to_Pa
    if(log10(maxval(temp_bulk_m (N_initial-5:N_initial))) < 6) temp_bulk_m  = temp_bulk_m  * GPa_to_Pa
end subroutine initial_read

subroutine number_of_layers
    implicit none
    integer :: il, ir, lnn(1:150), i
    real(8) :: rr(1:150)
    logical :: cond1
    
    il=0; ir=1;
    do i = 1, N_initial-1
        cond1 = abs(temp_radius(i+1) - temp_radius(i)) < 1d-8
        if(cond1) then
            il = il + 1
            rr(il) = temp_radius(i)
            lnn(il) = ir
            ir = 0
        end if
        ir = ir + 1
    end do
    numoflayers = il + 1;            !if(numoflayers/=5) stop "Number of layers do not correspond to Ganymede"
    lnn(il+1)   = ir
    rr(il+1)    = temp_radius(N_initial)

    allocate(layer(numoflayers), rad(1:numoflayers))
    do il = 1, numoflayers
        rad(il)     = rr(il)
        layer(il)%n = lnn(il)
    end do

end subroutine number_of_layers


subroutine assign_layers
    implicit none
    integer :: i, il, ir, nn

    i = 1
    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 1, nn
            layer(il)%r_i(ir)   = temp_radius  (i)
            layer(il)%rho_i(ir) = temp_density (i)
            layer(il)%K_i(ir)   = temp_bulk_m  (i)
            layer(il)%mu_i(ir)  = temp_shear_m (i)
            layer(il)%eta_i(ir) = temp_visc    (i)
            i = i + 1
        end do
    end do
end subroutine assign_layers


subroutine complete__
    implicit none
    call complete_time
    call complete_nr
end subroutine complete__

subroutine complete_time
    implicit none
    
    time_step = rotation_period / dble(num_time_steps)
    ang_vel   = 2.0d0 * Pi / rotation_period
end subroutine complete_time

subroutine complete_nr
    implicit none
    integer :: il

    do il = 1, numoflayers
        layer(il)%sm = layer(il)%n * 7 + 2
        layer(il)%tm = layer(il)%n * 3 + 1
    end do
    smtot = sum(layer(1:numoflayers)%sm)
    tmtot = sum(layer(1:numoflayers)%tm)

    total_scalar_harm = jm2 (Max_jmax,2)
    total_vector_harm = jml2(Max_jmax,2,Max_jmax+1)
    
    jstart = 2; jfinal = 2; status = 0
end subroutine complete_nr


subroutine reverse_temps
    implicit none
    call reverse_array(temp_radius  (:))
    call reverse_array(temp_density (:))
    call reverse_array(temp_bulk_m  (:))
    call reverse_array(temp_shear_m (:))
    call reverse_array(temp_visc    (:))
end subroutine reverse_temps


subroutine mid_layers
    implicit none
    integer :: il

    do il = 1, numoflayers
        call interface2layer(layer(il)%r_i  , layer(il)%r_l  )
        call interface2layer(layer(il)%rho_i, layer(il)%rho_l)
    end do
    layer(1)%r_l(1) = 0.0d0
end subroutine mid_layers

subroutine radial_grads
    implicit none
    integer :: il, ir, nn
    real(8) :: dr, drho

    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 1, nn
            dr   = layer(il)%r_l  (ir+1) - layer(il)%r_l  (ir)
            drho = layer(il)%rho_l(ir+1) - layer(il)%rho_l(ir)
            layer(il)%rho_r_i(ir) = drho / dr
        end do
    end do

    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 2, nn
            dr   = layer(il)%r_i  (ir  ) - layer(il)%r_i  (ir-1)
            drho = layer(il)%rho_i(ir  ) - layer(il)%rho_i(ir-1)
            layer(il)%rho_r_i(ir) = drho / dr
        end do
    end do
end subroutine radial_grads

subroutine SWITCH_imp
    implicit none
    integer :: il


    !! The rocky parts are always incompressible
    layer(1)%K_i = 1d30
    layer(2)%K_i = 1d30

    if(SWITCH_ocean_status == 'h') then
        layer(layer_ocean)%mu_i  = 1d2
        layer(layer_ocean)%eta_i = 1d30
    else
        layer(layer_ocean)%mu_i  = 1d30
        layer(layer_ocean)%eta_i = Ocean_viscosity
    end if

    if(SWITCH_solid_layer == 'e') then
        do il = 2, numoflayers
            if(il/=layer_ocean) layer(il)%eta_i = 1d30
        end do
    else
        do il = 2, numoflayers
            if(il/=layer_ocean) layer(il)%eta_i = 10d0**layer(2)%eta_i
        end do
    end if
end subroutine SWITCH_imp


subroutine find_ocean_layer
    implicit none
    integer :: il
    
    do il = 2, numoflayers
        if(((minval(layer(il)%eta_i)) < 8.0d0).and.((maxval(layer(il)%mu_i)) < 10.0d0)) then
            layer_ocean = il
        end if
    end do
end subroutine find_ocean_layer

subroutine hydrostatic_core
    implicit none
    
    layer(1)%mu_i  = 1d2
    layer(1)%eta_i = 1d30
end subroutine hydrostatic_core

subroutine compute_mass
    implicit none
    integer :: il, ir, nn
    real(8) :: f, vol

    f   = 4.0d0 * Pi / 3.0d0
    vol = f * layer(1)%r_i(1)
    layer(1)%mass_i(1) = layer(1)%rho_i(1)*vol
    do il = 1, numoflayers
        nn = layer(il)%n
        if(il >= 2) layer(il)%mass_i(1) = layer(il-1)%mass_i(layer(il-1)%n)
        do ir = 2, layer(il)%n
            vol = f * (layer(il)%r_i(ir)**3 - layer(il)%r_i(ir-1)**3)
            layer(il)%mass_i(ir) = layer(il)%mass_i(ir-1) + vol *layer(il)%rho_i(ir)
        end do
        call interface2layer(layer(il)%mass_i, layer(il)%mass_l)
    end do
    layer(1)%mass_l(1) = 0.0d0
end subroutine compute_mass

subroutine compute_grav_acc
    implicit none
    integer :: il,ir, sp

    do il = 1, numoflayers
        do ir = 1, layer(il)%n
            layer(il)%g_i(ir) = layer(il)%mass_i(ir) * Grav_Const / layer(il)%r_i(ir)**2
        end do
    end do
    do il = 1, numoflayers
        if(il==1)sp=2
        if(il/=1)sp=1
        do ir = sp, layer(il)%n+1
            layer(il)%g_l(ir) = layer(il)%mass_l(ir) * Grav_Const / layer(il)%r_l(ir)**2
        end do
    end do
    layer(1)%g_l(1) = 0.0d0
end subroutine compute_grav_acc

subroutine compute_MOI
    implicit none
    integer :: i, ir
    real(8) :: fac, dr
    
    fac=8.0d0*pi/(3.0d0*(rad(numoflayers)**2)*layer(numoflayers)%mass_i(layer(numoflayers)%n))

    moi = 0.0d0
    do i = 1, numoflayers
        do ir = 1, layer(i)%n-1
            dr = layer(i)%r_i(ir+1) - layer(i)%r_i(ir)
            moi = moi + (layer(i)%rho_i(ir+1)*layer(i)%r_i(ir+1)**4 + layer(i)%rho_i(ir)*layer(i)%r_i(ir)**4)*0.5d0*dr
        end do
    end do
    moi = moi*fac
end subroutine compute_MOI

subroutine determine_mindr
    implicit none
    integer :: ir, il, nn

    min_dr = 1d20
    do il = 1, numoflayers
        nn = layer(il)%n
        do ir = 1, nn
            if ( abs(layer(il)%r_l(ir+1) - layer(il)%r_l(ir)) < min_dr) then
                min_dr = abs(layer(il)%r_l(ir+1) - layer(il)%r_l(ir))
            end if
        end do
    end do
end subroutine determine_mindr

end module mod_preprocess