module mod_variables
    use mod_parameters
    implicit none
    integer :: numoflayers
    integer :: layer_ocean

    !!Switches
    integer                :: SWITCH_self_grav
    integer                :: SWITCH_initial_run
    integer                :: SWITCH_inertia
    integer                :: SWITCH_coriolis
    integer                :: SWITCH_compressible
    integer                :: SWITCH_velocity
    integer                :: SWITCH_mmtides
    character(1)           :: SWITCH_solid_layer   !! Elastic: e; Maxwell: m; Andrade: a
    character(1)           :: SWITCH_ocean_status  !! Viscous: v; Hydrostatic: h

    !! Spherical harmonic degrees
    integer                :: jstart, jfinal, jmax, old_jmax
    integer                :: total_scalar_harm
    integer                :: total_vector_harm

    !! Time stepping
    integer                :: num_time_steps
    integer                :: pc
    integer                :: signal
    real(8)                :: time_step
    real(8)                :: t
    real(8)                :: min_dr
    real(8)                :: weight_sg

    !! The parameters to be read from the file
    real(8)                :: ecc
    real(8)                :: ang_vel
    real(8)                :: rotation_period
    real(8)                :: final_period
    real(8)                :: moi
    real(8)                :: Ocean_viscosity

    !! Radial profiles from files
    real(8), allocatable  :: temp_radius(:), temp_density(:), temp_shear_m(:), temp_bulk_m(:), temp_visc(:)
    integer               :: N_initial

    
    !!Layers
    real(8), allocatable  :: rad(:)
    type :: spherical_layers
        integer :: n, sm, tm
        !!radial profiles
        real(8), allocatable    ::  r_i     (:),     r_l(:), &
                                    rho_i   (:),   rho_l(:), &
                                    rho_r_i (:), rho_r_l(:), &
                                    g_i     (:),     g_l(:), &
                                    mass_i  (:),  mass_l(:), &
                                    mu_i    (:)            , &
                                    eta_i   (:)            , &
                                    K_i     (:)
        !! unknown variables
        complex(8), allocatable :: rhoe(:,:), str(:,:), vel(:,:), pvel(:,:)
        !! Biproduct of the unknown variablse
        complex(8), allocatable :: dis(:,:), ur_i(:,:), pot(:,:), grad_pot(:,:), disl(:,:)
        complex(8), allocatable ::  fc(:,:),  fc0(:,:), temp_array(:,:), fc2(:,:), fc20(:,:), fc1(:,:)
        complex(8), allocatable :: maxwell_term(:,:), andrade_term(:,:), str_a(:,:,:)
        complex(8), allocatable :: lap_density(:,:), rhoe_l(:,:)
        complex(8), allocatable :: divergence_u20(:)
        real(8), allocatable    :: andrade_CC(:), spectrum(:), temp_array2(:), spectrum_v(:)
        real(8)                 :: dissipation
        complex(8), allocatable :: div_tensor(:,:), div_v(:,:), grad_v(:,:)
    end type

    type(spherical_layers), allocatable :: layer(:)


    !! Declaration of the matrix variables
    integer                 :: smtot, tmtot
    real(8)   , allocatable :: S_U(:,:), T_U(:,:), S_L(:,:), T_L(:,:)
    integer   , allocatable :: S_I(:), T_I(:)
    type :: rhs_array
        complex(8), allocatable :: spheriodal(:)
        complex(8), allocatable :: torodial(:)
    end type

    type(rhs_array) :: RHS

    !! input section
    character(200) :: file_radial_profile
    character(200) :: dir_output
    character(200) :: body_name
    !! Abdrade parameters
    real(8) :: alpha, zeta

    !!Love numbers
    real(8) :: k20, real_k22, imag_k22
    real(8) :: h20, real_h22, imag_h22
    real(8) :: l20, real_l22, imag_l22
    real(8), allocatable   :: kh(:), khp(:,:), kh_per(:,:)
    real(8), allocatable   :: gp20(:), rgp22(:), igp22(:)
    real(8), allocatable   :: ur20(:), rur22(:), iur22(:)
    real(8), allocatable   :: ut20(:), rut22(:), iut22(:)
    real(8)                :: k20_p(1:Max_num_period), rk22_p(1:Max_num_period), ik22_p(1:Max_num_period)
    real(8)                :: h20_p(1:Max_num_period), rh22_p(1:Max_num_period), ih22_p(1:Max_num_period)
    real(8)                :: l20_p(1:Max_num_period), rl22_p(1:Max_num_period), il22_p(1:Max_num_period)
    real(8)                :: dissip_p(1:Max_num_period)
    real(8)                :: avrg_dissip
    real(8)                :: gp20_0, gp20_1, gp20_2, rgp22_0, rgp22_1, rgp22_2, igp22_0, igp22_1, igp22_2
    real(8)                :: peak_time_1, peak_time_2, pt_d_1, pt_d_2
    real(8)                :: dissipation_0, dissipation_1, dissipation_2

    !!For opening and writing files
    integer :: iu(1:20)
    integer :: status
    integer :: WRITE_density, WRITE_velocity, WRITE_stress, WRITE_displacement
    integer :: checkpoint, start_output, write_output, start_simulation, start_cp

    !! Harmonics to grid
    type(spherical_layers), allocatable :: layer_grid(:)
    real(8), allocatable :: theta(:), phi(:)
    real(8), allocatable :: rho_phi(:,:) !! density for given phi
    complex(8), allocatable :: velocity_ocean_m2(:,:)

contains
pure integer function jml2(j,m,l)
    implicit none
    integer, intent(in) :: j, m, l

    if(j.eq.0) then
        jml2 = 1
        return
    end if

    if(j.eq.1) then
        jml2 = l + 2
        return
    end if

    if(m.eq.0) then
        jml2 = (2*j - 2)*3 + l - j
    else
        jml2 = (2*j - 1)*3 + l - j
    end if
end function jml2

pure integer function jm2(j,m)
    implicit none
    integer, intent(in) :: j, m

    if(j.eq.0) then
        jm2 = 1
        return
    end if
    if(j.eq.1) then
        jm2 = 2
        return
    end if
    if(m.eq.0) then
        jm2 = 2*j-1 
    else
        jm2 = 2*j
    end if    
end function jm2

pure integer function dirac(a, b)
    implicit none
    integer, intent(in) :: a, b
    if (a.eq.b) then
        dirac = 1
    else
        dirac = 0
    end if
end function 
end module mod_variables