constitutive_2d.f90

!> @author
!> Mattia de' Michieli Vitturi
!> \date 15/08/2011
!>
!> Modification :: add friction weakening law
!> @Author T. Frasson & A. Lucas
!> \data 15/07/2020


!********************************************************************************
!> \brief Constitutive equations
!********************************************************************************
MODULE constitutive_2d

  USE parameters_2d, ONLY : n_eqns , n_vars
  USE parameters_2d, ONLY : rheology_flag , rheology_model
  USE parameters_2d, ONLY : temperature_flag
  USE parameters_2d, ONLY : solid_transport_flag

  IMPLICIT none

  LOGICAL, ALLOCATABLE :: implicit_flag(:)

  CHARACTER(LEN=20) :: phase1_name
  CHARACTER(LEN=20) :: phase2_name


  COMPLEX*16 :: h      !< height [m]
  COMPLEX*16 :: u      !< velocity (x direction)
  COMPLEX*16 :: v      !< velocity (y direction)
  COMPLEX*16 :: T      !< temperature
  COMPLEX*16 :: xs     !< sediment concentration

  !> gravitational acceleration
  REAL*8 :: grav

  !> drag coefficients (Voellmy-Salm model)
  REAL*8 :: mu
  REAL*8 :: xi

  !> drag coefficients (plastic model)
  REAL*8 :: tau

  !> evironment temperature [K]
  REAL*8 :: T_env

  !> radiative coefficient
  REAL*8 :: rad_coeff

  !> friction coefficient
  REAL*8 :: frict_coeff

  !> fluid density [kg/m3]
  REAL*8 :: rho

  !> reference temperature [K]
  REAL*8 :: T_ref

  !> reference kinematic viscosity [m2/s]
  REAL*8 :: nu_ref

  !> viscosity parameter [K-1] (b in Table 1 Costa & Macedonio, 2005)
  REAL*8 :: visc_par

  !> velocity boundary layer fraction of total thickness
  REAL*8 :: emme

  !> specific heat [J kg-1 K-1]
  REAL*8 :: c_p

  !> atmospheric heat trasnfer coefficient [W m-2 K-1] (lambda in C&M, 2005)
  REAL*8 :: atm_heat_transf_coeff

  !> fractional area of the exposed inner core (f in C&M, 2005)
  REAL*8 :: exp_area_fract

  !> Stephan-Boltzmann constant [W m-2 K-4]
  REAL*8, PARAMETER :: SBconst = 5.67D-8

  !> emissivity (eps in Costa & Macedonio, 2005)
  REAL*8 :: emissivity

  !> thermal boundary layer fraction of total thickness
  REAL*8 :: enne

  !> temperature of lava-ground interface
  REAL*8 :: T_ground

  !> thermal conductivity [W m-1 K-1] (k in Costa & Macedonio, 2005)
  REAL*8 :: thermal_conductivity

  !--- Lahars rheology model parameters

  !> 1st parameter for yield strenght empirical relationship (O'Brian et al, 1993)
  REAL*8 :: alpha2

  !> 2nd parameter for yield strenght empirical relationship (O'Brian et al, 1993)
  REAL*8 :: beta2

  !> 1st parameter for fluid viscosity empirical relationship (O'Brian et al, 1993)
  REAL*8 :: alpha1

  !> 2nd parameter for fluid viscosity empirical relationship (O'Brian et al, 1993)
  REAL*8 :: beta1

  !> Empirical resistance parameter
  REAL*8 :: Kappa

  !> Mannings roughness coefficient ( units: T L^(-1/3) )
  REAL*8 :: n_td

  !> Specific weight of water
  REAL*8 :: gamma_w

  !> Specific weight of sediments
  REAL*8 :: gamma_s

  !> parametres Pouliquen
  REAL*8 :: beta
  REAL*8 :: theta1
  REAL*8 :: theta2
  REAL*8 :: hstop
  REAL*8 :: L
  REAL*8 :: d

  !> parametres Bingham
  REAL*8 :: eta
  REAL*8 :: thet
  REAL*8 :: tau_bing

  !> parametres weakening

  REAL*8 :: mu_0
  REAL*8 :: mu_w
  REAL*8 :: U_0



CONTAINS

  !******************************************************************************
  !> \brief Initialization of relaxation flags
  !
  !> This subroutine set the number and the flags of the non-hyperbolic
  !> terms.
  !> \date 07/09/2012
  !******************************************************************************

  SUBROUTINE init_problem_param

    USE parameters_2d, ONLY : n_nh

    ALLOCATE( implicit_flag(n_eqns) )

    implicit_flag(1:n_eqns) = .FALSE.
    implicit_flag(2) = .TRUE.
    implicit_flag(3) = .TRUE.

    IF ( solid_transport_flag ) THEN

       IF ( temperature_flag ) implicit_flag(5) = .TRUE.

    ELSE

       IF ( temperature_flag ) implicit_flag(4) = .TRUE.

    END IF

    n_nh = COUNT( implicit_flag )

  END SUBROUTINE init_problem_param

  !******************************************************************************
  !> \brief Physical variables
  !
  !> This subroutine evaluates from the conservative local variables qj
  !> the local physical variables  (\f$h+B, u, v, xs , T \f$).
  !> \param[in]    r_qj     real conservative variables
  !> \param[in]    c_qj     complex conservative variables
  !> @author
  !> Mattia de' Michieli Vitturi
  !> \date 15/08/2011
  !******************************************************************************

  SUBROUTINE phys_var(Bj,r_qj,c_qj)

    USE COMPLEXIFY
    USE parameters_2d, ONLY : eps_sing
    IMPLICIT none

    REAL*8, INTENT(IN) :: Bj
    REAL*8, INTENT(IN), OPTIONAL :: r_qj(n_vars)
    COMPLEX*16, INTENT(IN), OPTIONAL :: c_qj(n_vars)

    COMPLEX*16 :: qj(n_vars)

    IF ( present(c_qj) ) THEN

       qj = c_qj