constitutive_2d.f90

Go to the documentation of this file.

!********************************************************************************
!********************************************************************************
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
  USE parameters_2d, ONLY : sed_vol_perc

  IMPLICIT none

  LOGICAL, ALLOCATABLE :: implicit_flag(:)

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

  COMPLEX*16 :: h
  COMPLEX*16 :: u
  COMPLEX*16 :: v
  COMPLEX*16 :: t
  COMPLEX*16 :: alphas
  COMPLEX*16 :: rho_m

  REAL*8 :: grav

  REAL*8 :: mu
  REAL*8 :: xi
  
  REAL*8 :: tau

  REAL*8 :: t_env

  REAL*8 :: rad_coeff

  REAL*8 :: frict_coeff

  REAL*8 :: rho

  REAL*8 :: t_ref

  REAL*8 :: nu_ref

  REAL*8 :: visc_par

  REAL*8 :: emme

  REAL*8 :: c_p

  REAL*8 :: atm_heat_transf_coeff

  REAL*8 :: exp_area_fract

  REAL*8, PARAMETER :: sbconst = 5.67d-8

  REAL*8 :: emissivity

  REAL*8 :: enne

  REAL*8 :: t_ground

  REAL*8 :: thermal_conductivity

  !--- Lahars rheology model parameters

  REAL*8 :: alpha2

  REAL*8 :: beta2

  REAL*8 :: alpha1

  REAL*8 :: beta1

  REAL*8 :: kappa

  REAL*8 :: n_td
 
  REAL*8 :: rho_w

  REAL*8 :: rho_s

  REAL*8 :: settling_vel

  REAL*8 :: erosion_coeff


CONTAINS

  !******************************************************************************
  !
  !******************************************************************************

  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

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE phys_var(r_qj,c_qj)
    
    USE complexify
    USE parameters_2d, ONLY : eps_sing
    IMPLICIT none
    
    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

    ELSE

       qj = dcmplx(r_qj)

    END IF

    h = qj(1)

    IF ( solid_transport_flag ) THEN

       IF ( dble( h ) .GT. 1.d-25 ) THEN
          
          alphas = qj(4) / h
          
          IF ( temperature_flag ) t = qj(5) / h
          
       ELSE
          
          alphas = dsqrt(2.d0) * h * qj(4) / cdsqrt( h**4 + eps_sing )
          
          IF ( temperature_flag ) THEN
             
             t =  dsqrt(2.d0) * h * qj(5) / cdsqrt( h**4 + eps_sing )
             
          END IF
          
       END IF

       rho_m = ( dcmplx(1.d0,0.d0) - alphas ) * rho_w + alphas * rho_s 

    ELSE

       IF ( temperature_flag ) THEN

          IF ( dble( h ) .GT. 1.d-25 ) THEN
          
             t = qj(4) / h
          
          ELSE
          
             t =  dsqrt(2.d0) * h * qj(4) / cdsqrt( h**4 + eps_sing )
             
          END IF
          
       END IF

       alphas = dcmplx(0.d0,0.d0) 
       rho_m = dcmplx(1.d0,0.d0) 

    END IF

    IF ( dble( h ) .GT. eps_sing ** 0.25d0 ) THEN

       u = qj(2) / ( rho_m * h )

       v = qj(3) / ( rho_m * h )

    ELSE

       u = dsqrt(2.d0) * h * ( qj(2) / rho_m ) / cdsqrt( h**4 + eps_sing )

       v = dsqrt(2.d0) * h * ( qj(3) / rho_m ) / cdsqrt( h**4 + eps_sing )

    END IF

    RETURN

  END SUBROUTINE phys_var

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_local_speeds_x(qj,vel_min,vel_max)
  
    IMPLICIT none
    
    REAL*8, INTENT(IN)  :: qj(n_vars)

    REAL*8, INTENT(OUT) :: vel_min(n_vars) , vel_max(n_vars)

    CALL phys_var(r_qj = qj)
    
    vel_min(1:n_eqns) = dble(u) - dsqrt( grav * dble(h) )
    vel_max(1:n_eqns) = dble(u) + dsqrt( grav * dble(h) )
        
  END SUBROUTINE eval_local_speeds_x

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_local_speeds_y(qj,vel_min,vel_max)
    
    IMPLICIT none
    
    REAL*8, INTENT(IN)  :: qj(n_vars)
    REAL*8, INTENT(OUT) :: vel_min(n_vars) , vel_max(n_vars)
    
    CALL phys_var(r_qj = qj)
    
    vel_min(1:n_eqns) = dble(v) - dsqrt( grav * dble(h) )
    vel_max(1:n_eqns) = dble(v) + dsqrt( grav * dble(h) )
       
  END SUBROUTINE eval_local_speeds_y

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE qc_to_qp(qc,B,qp)
    
    IMPLICIT none
    
    REAL*8, INTENT(IN) :: qc(n_vars)
    REAL*8, INTENT(IN) :: B
    REAL*8, INTENT(OUT) :: qp(n_vars)
    
    CALL phys_var(r_qj = qc)
    
    ! It is important to have as reconstructed variable at the interfaces
    ! the total height h+B, instead of h. Only in this way, when the topography
    ! is not flat and the free surface is horizontal (steady condition), the  
    ! latter is kept flat by the reconstruction and the solution is stable.
    qp(1) = dble(h+b)
    qp(2) = dble(u)
    qp(3) = dble(v)

    IF ( solid_transport_flag ) THEN
       
       qp(4) = dble(alphas)
       
       IF ( temperature_flag ) qp(5) = dble(t)
       
    ELSE
       
       IF ( temperature_flag ) qp(4) = dble(t)
       
    END IF
    
    
  END SUBROUTINE qc_to_qp

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE qp_to_qc(qp,B,qc)
    
    USE complexify 
    IMPLICIT none
    
    REAL*8, INTENT(IN) :: qp(n_vars)
    REAL*8, INTENT(IN) :: B
    REAL*8, INTENT(OUT) :: qc(n_vars)
    
    REAL*8 :: r_h
    REAL*8 :: r_u
    REAL*8 :: r_v
    REAL*8 :: r_T
    REAL*8 :: r_rho_m
    REAL*8 :: r_alphas
 
    r_h = qp(1) - b
    r_u  = qp(2)
    r_v  = qp(3)


    qc(1) = r_h

    IF ( solid_transport_flag ) THEN
       
       r_alphas = qp(4)
       qc(4) = r_h * r_alphas 
       
       IF ( temperature_flag ) THEN
          
          r_t  = qp(5)
          qc(5) = r_h * r_t 
       
       END IF
       
       ! TO BE REMOVED
       ! r_alphas = sed_vol_perc/100.D0

       r_rho_m = r_alphas * rho_s + ( 1.d0 - r_alphas ) * rho_w

    ELSE

       IF ( temperature_flag ) THEN
          
          r_t  = qp(4)
          qc(4) = r_h * r_t 
       
       END IF

       r_alphas = 0.d0
       r_rho_m = 1.d0

    END IF

    qc(2) = r_h * r_rho_m * r_u
    qc(3) = r_h * r_rho_m * r_v


  END SUBROUTINE qp_to_qc

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE qc_to_qc2(qc,B,qc2)
    
    IMPLICIT none
    
    REAL*8, INTENT(IN) :: qc(n_vars)
    REAL*8, INTENT(IN) :: B
    REAL*8, INTENT(OUT) :: qc2(n_vars)

    qc2(1) = qc(1) + b
    qc2(2:n_vars) = qc(2:n_vars)

  END SUBROUTINE qc_to_qc2

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE qc2_to_qc(qc2,B,qc)
    
    IMPLICIT none
    
    REAL*8, INTENT(IN) :: qc2(n_vars)
    REAL*8, INTENT(IN) :: B
    REAL*8, INTENT(OUT) :: qc(n_vars)

    REAL*8 :: qc_temp(n_vars)

    REAL*8 :: r_h
    REAL*8 :: r_u
    REAL*8 :: r_v
    REAL*8 :: r_alphas
    REAL*8 :: r_T
    REAL*8 :: r_rho_m

    qc_temp(1) = qc2(1) - b
    qc_temp(2:n_vars) = qc2(2:n_vars)

    ! Desingularization
    CALL phys_var(r_qj = qc_temp)
          
    r_h = dble(h)
    r_u = dble(u)
    r_v = dble(v)
    r_alphas = dble(alphas)

    ! TO BE REMOVED
    ! r_alphas = sed_vol_perc/100.D0

    IF ( solid_transport_flag ) THEN

       r_rho_m = dble(rho_m)

    ELSE

       r_rho_m = 1.d0

    END IF

    qc(1) = r_h
    qc(2) = r_h * r_rho_m * r_u
    qc(3) = r_h * r_rho_m * r_v
    
    IF ( solid_transport_flag ) THEN
       
       qc(4) = r_h * r_alphas
       
       IF ( temperature_flag ) THEN
          
          r_t = dble(t)
          qc(5) = r_h * r_t 
          
       END IF
       
    ELSE
       
       IF ( temperature_flag ) THEN
          
          r_t = dble(t)
          qc(4) = r_h * r_t 
          
       END IF
       
    END IF
    
  END SUBROUTINE qc2_to_qc


  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE eval_fluxes(c_qj,r_qj,c_flux,r_flux,dir)
    
    USE complexify
    IMPLICIT none

    COMPLEX*16, INTENT(IN), OPTIONAL :: c_qj(n_vars)
    COMPLEX*16, INTENT(OUT), OPTIONAL :: c_flux(n_eqns)
    REAL*8, INTENT(IN), OPTIONAL :: r_qj(n_vars)
    REAL*8, INTENT(OUT), OPTIONAL :: r_flux(n_eqns)
    INTEGER, INTENT(IN) :: dir

    COMPLEX*16 :: qj(n_vars)
    COMPLEX*16 :: flux(n_eqns)
    COMPLEX*16 :: h_temp , u_temp , v_temp
    COMPLEX*16 :: alphas_temp , rho_m_temp

    INTEGER :: i 

    IF ( present(c_qj) .AND. present(c_flux) ) THEN

       qj = c_qj

    ELSEIF ( present(r_qj) .AND. present(r_flux) ) THEN

       DO i = 1,n_vars

          qj(i) = dcmplx( r_qj(i) )

       END DO

    ELSE

       WRITE(*,*) 'Constitutive, eval_fluxes: problem with arguments'
       stop

    END IF

    h_temp = qj(1)

    pos_thick:IF ( dble(h_temp) .NE. 0.d0 ) THEN

       IF ( solid_transport_flag ) THEN
          
          alphas_temp = qj(4) / h_temp
          
          ! TO BE REMOVED
          ! alphas_temp = DCMPLX(sed_vol_perc/100.D0,0.D0)
          
          rho_m_temp = alphas_temp * rho_s + ( dcmplx(1.d0,0.d0)                 &
               - alphas_temp ) * rho_w
          
       ELSE
          
          alphas_temp = dcmplx( 0.d0 , 0.d0 )
          rho_m_temp = dcmplx(1.d0,0.d0) 
          
       END IF
       
       IF ( dir .EQ. 1 ) THEN

          ! Velocity in x-direction
          u_temp = qj(2) / ( h_temp * rho_m_temp )

          ! Volumetric flux in x-direction: h*u
          flux(1) = qj(2) / rho_m_temp

          ! x-momentum flux in x-direction + hydrostatic pressure term
          flux(2) = rho_m_temp * h_temp * u_temp**2 +                           &
               0.5d0 * rho_m_temp * grav * h_temp**2

          ! y-momentum flux in x-direction: rho*h*v*u
          flux(3) = u_temp * qj(3)

          IF ( solid_transport_flag ) THEN

             ! Volumetric flux of solid in x-direction: h * alphas * u
             flux(4) = u_temp * qj(4)

             ! Solid flux can't be larger than total flux
             IF ( flux(4) / flux(1) .GT. 1.d0 ) flux(4) = flux(1)
             
             ! Temperature flux in x-direction: U*T*h
             IF ( temperature_flag ) flux(5) = u_temp * qj(5)
             
          ELSE

             ! Temperature flux in x-direction: U*T*h
             IF ( temperature_flag ) flux(4) = u_temp * qj(4)
              
          END IF

       ELSEIF ( dir .EQ. 2 ) THEN

          ! flux G (derivated wrt y in the equations)

          v_temp = qj(3) / ( h_temp * rho_m_temp )
          
          flux(1) = qj(3) / rho_m_temp
          
          flux(2) = v_temp * qj(2)
          
          flux(3) = rho_m_temp * h_temp * v_temp**2 +                           &
               0.5d0 * rho_m_temp * grav * h_temp**2
          
          IF ( solid_transport_flag ) THEN

             flux(4) = v_temp * qj(4)

             ! Solid flux can't be larger than total flux
             IF ( flux(4) / flux(1) .GT. 1.d0 ) flux(4) = flux(1)

             ! Temperature flux in y-direction: V*T*h
             IF ( temperature_flag ) flux(5) = v_temp * qj(5)
             
          ELSE

             ! Temperature flux in y-direction: V*T*h
             IF ( temperature_flag ) flux(4) = v_temp * qj(4)

          END IF

       END IF

    ELSE
       
       flux(1:n_eqns) = 0.d0
       
    ENDIF pos_thick
 
    !WRITE(*,*) 'flux',DBLE(flux)
    !WRITE(*,*) 'qj',DBLE(qj)
    !WRITE(*,*) 'Bj',Bj
    !WRITE(*,*) 'rho_m_temp , h_temp', REAL(rho_m_temp) , REAL(h_temp)
    !READ(*,*)

    IF ( present(c_qj) .AND. present(c_flux) ) THEN

       c_flux = flux

    ELSEIF ( present(r_qj) .AND. present(r_flux) ) THEN

       r_flux = dble( flux )

    END IF

  END SUBROUTINE eval_fluxes

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_nonhyperbolic_terms( c_qj , c_nh_term_impl , r_qj ,           &
       r_nh_term_impl )

    USE complexify 

    IMPLICIT NONE

    COMPLEX*16, INTENT(IN), OPTIONAL :: c_qj(n_vars)

    COMPLEX*16, INTENT(OUT), OPTIONAL :: c_nh_term_impl(n_eqns)
    REAL*8, INTENT(IN), OPTIONAL :: r_qj(n_vars)
    REAL*8, INTENT(OUT), OPTIONAL :: r_nh_term_impl(n_eqns)

    COMPLEX*16 :: qj(n_vars)

    COMPLEX*16 :: nh_term(n_eqns)

    COMPLEX*16 :: relaxation_term(n_eqns)

    COMPLEX*16 :: forces_term(n_eqns)

    INTEGER :: i

    COMPLEX*16 :: mod_vel
    
    COMPLEX*16 :: gamma

    REAL*8 :: radiative_coeff

    COMPLEX*16 :: radiative_term

    REAL*8 :: convective_coeff

    COMPLEX*16 :: convective_term

    COMPLEX*16 :: conductive_coeff , conductive_term

    REAL*8 :: thermal_diffusivity

    REAL*8 :: h_threshold

    !--- Lahars rheology model variables
    
    COMPLEX*8 :: fluid_visc

    COMPLEX*8 :: s_f

    COMPLEX*8 :: s_v

    COMPLEX*8 :: s_td


    IF ( temperature_flag ) THEN

       IF ( ( thermal_conductivity .GT. 0.d0 ) .OR. ( emme .GT. 0.d0 ) ) THEN

          h_threshold = 1.d-10

       ELSE

          h_threshold = 0.d0

       END IF
       
    END IF
       

    IF ( present(c_qj) .AND. present(c_nh_term_impl) ) THEN

       qj = c_qj

    ELSEIF ( present(r_qj) .AND. present(r_nh_term_impl) ) THEN

       DO i = 1,n_vars

          qj(i) = dcmplx( r_qj(i) )

       END DO

    ELSE

       WRITE(*,*) 'Constitutive, eval_fluxes: problem with arguments'
       stop

    END IF

    ! initialize and evaluate the relaxation terms
    relaxation_term(1:n_eqns) = dcmplx(0.d0,0.d0) 

    ! initialize and evaluate the forces terms
    forces_term(1:n_eqns) = dcmplx(0.d0,0.d0)

    IF (rheology_flag) THEN

       CALL phys_var(c_qj = qj)
    
       mod_vel = cdsqrt( u**2 + v**2 )
       
       ! Voellmy Salm rheology
       IF ( rheology_model .EQ. 1 ) THEN
       
          IF ( dble(mod_vel) .NE. 0.d0 ) THEN 
          
             ! IMPORTANT: grav3_surv is always negative 
             forces_term(2) = forces_term(2) - rho_m * ( u / mod_vel ) *        &
                  ( grav / xi ) * mod_vel ** 2
          
             forces_term(3) = forces_term(3) - rho_m * ( v / mod_vel ) *        &
                  ( grav / xi ) * mod_vel ** 2
          
          ENDIF
        
       ! Plastic rheology
       ELSEIF ( rheology_model .EQ. 2 ) THEN
       
          IF ( dble(mod_vel) .NE. 0.d0 ) THEN 
          
             forces_term(2) = forces_term(2) - rho_m * tau * (u/mod_vel)
          
             forces_term(3) = forces_term(3) - rho_m * tau * (v/mod_vel)

          ENDIF

       ! Temperature dependent rheology
       ELSEIF ( rheology_model .EQ. 3 ) THEN

          IF ( dble(h) .GT. h_threshold ) THEN
    
             ! Equation 6 from Costa & Macedonio, 2005
             gamma = 3.d0 * nu_ref / h * cdexp( - visc_par * ( t - t_ref ) )

          ELSE

             ! Equation 6 from Costa & Macedonio, 2005
             gamma = 3.d0 * nu_ref / h_threshold * cdexp( - visc_par            &
                  * ( t - t_ref ) )
             
          END IF
          
          IF ( dble(mod_vel) .NE. 0.d0 ) THEN 
          
             ! Last R.H.S. term in equation 2 from Costa & Macedonio, 2005
             forces_term(2) = forces_term(2) - rho_m * gamma * u
          
             ! Last R.H.S. term in equation 3 from Costa & Macedonio, 2005
             forces_term(3) = forces_term(3) - rho_m * gamma * v

          ENDIF
          
       ! Lahars rheology (O'Brien 1993, FLO2D)
       ELSEIF ( rheology_model .EQ. 4 ) THEN

          h_threshold = 1.d-20

          ! THIS IS ONLY USED WHEN SEDIMENT PERCENTAGE IS GIVEN AS INPUT AND
          ! IT IS CONSTANT WITHIN THE SIMULATION
          ! TO BE REMOVED
          ! alphas = DCMPLX(sed_vol_perc/100.D0,0.D0)

          ! Fluid viscosity
          fluid_visc = alpha1 * cdexp( beta1 * alphas )

          IF ( h .GT. h_threshold ) THEN
             
             ! Viscous slope component
             s_v = kappa * fluid_visc * mod_vel / ( 8.d0 * h**2 )
             
             ! Turbulent dispersive component
             s_td = rho_m * n_td**2 * mod_vel**2 / ( h**(4.d0/3.d0) )
          
          ELSE
             
             ! Viscous slope component
             s_v = kappa * fluid_visc * mod_vel / ( 8.d0 * h_threshold**2 )
             
             ! Turbulent dispersive components
             s_td = rho_m * n_td**2 * (mod_vel**2) / ( h_threshold**(4.d0/3.d0) )
             
          END IF
          
          ! Total implicit friction slope
          s_f = s_v + s_td
         
          IF ( mod_vel .GT. 0.d0 ) THEN

             forces_term(2) = forces_term(2) - grav * h * ( u / mod_vel ) * s_f

             forces_term(3) = forces_term(3) - grav * h * ( v / mod_vel ) * s_f
  
          END IF

       ELSEIF ( rheology_model .EQ. 5 ) THEN

          tau = 1.d-3 / ( 1.d0 + 10.d0 * h ) * mod_vel
          
          IF ( dble(mod_vel) .NE. 0.d0 ) THEN
             
             forces_term(2) = forces_term(2) - rho_m * tau * ( u / mod_vel )
             forces_term(3) = forces_term(3) - rho_m * tau * ( v / mod_vel )

          END IF
          
       ENDIF
              
    ENDIF

    IF ( temperature_flag ) THEN

       CALL phys_var(c_qj = qj)

       IF ( dble(h) .GT. 0.d0 ) THEN

          ! Equation 8 from Costa & Macedonio, 2005
          radiative_coeff = emissivity * sbconst * exp_area_fract / ( rho * c_p )
       
       ELSE

          radiative_coeff = 0.d0

       END IF

       IF ( dble(t) .GT. t_env ) THEN

          ! First R.H.S. term in equation 4 from Costa & Macedonio, 2005
          radiative_term = - radiative_coeff * ( t**4 - t_env**4 )

       ELSE

          radiative_term = dcmplx(0.d0,0.d0)

       END IF

       IF ( dble(h) .GT. 0.d0 ) THEN
       
          ! Equation 9 from Costa & Macedonio, 2005
          convective_coeff = atm_heat_transf_coeff * exp_area_fract             &
               / ( rho * c_p )

       ELSE

          convective_coeff = 0.d0

       END IF

       IF ( dble(t) .GT. t_env ) THEN

          ! Second R.H.S. term in equation 4 from Costa & Macedonio, 2005
          convective_term = - convective_coeff * ( t - t_env )

       ELSE

          convective_term =  dcmplx(0.d0,0.d0)

       END IF

       IF ( dble(h) .GT. h_threshold ) THEN
    
          thermal_diffusivity = thermal_conductivity / ( rho * c_p ) 

          ! Equation 7 from Costa & Macedonio, 2005
          conductive_coeff = enne * thermal_diffusivity / h

       ELSE

          conductive_coeff =  dcmplx(0.d0,0.d0)
          conductive_coeff = enne * thermal_diffusivity / dcmplx(h_threshold,0.d0)

       END IF

       ! Third R.H.S. term in equation 4 from Costa & Macedonio, 2005
       IF ( dble(t) .GT. t_ground ) THEN

          conductive_term = - conductive_coeff * ( t - t_ground )

       ELSE

           conductive_term = dcmplx(0.d0,0.d0)

        END IF

        IF ( solid_transport_flag ) THEN

           relaxation_term(5) = radiative_term + convective_term + conductive_term

        ELSE

           relaxation_term(4) = radiative_term + convective_term + conductive_term

        END IF

    END IF

    nh_term = relaxation_term + forces_term

    IF ( present(c_qj) .AND. present(c_nh_term_impl) ) THEN

       c_nh_term_impl = nh_term

    ELSEIF ( present(r_qj) .AND. present(r_nh_term_impl) ) THEN

       r_nh_term_impl = dble( nh_term )

    END IF 

  END SUBROUTINE eval_nonhyperbolic_terms

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_nh_semi_impl_terms( grav3_surf , c_qj , c_nh_semi_impl_term , &
       r_qj , r_nh_semi_impl_term )

    USE complexify 


    IMPLICIT NONE

    REAL*8, INTENT(IN) :: grav3_surf

    COMPLEX*16, INTENT(IN), OPTIONAL :: c_qj(n_vars)
    COMPLEX*16, INTENT(OUT), OPTIONAL :: c_nh_semi_impl_term(n_eqns)
    REAL*8, INTENT(IN), OPTIONAL :: r_qj(n_vars)
    REAL*8, INTENT(OUT), OPTIONAL :: r_nh_semi_impl_term(n_eqns)

    COMPLEX*16 :: qj(n_vars)

    COMPLEX*16 :: forces_term(n_eqns)

    INTEGER :: i

    COMPLEX*16 :: mod_vel
    
    REAL*8 :: h_threshold

    !--- Lahars rheology model variables
    
    COMPLEX*8 :: tau_y

    COMPLEX*8 :: s_y


    IF ( present(c_qj) .AND. present(c_nh_semi_impl_term) ) THEN

       qj = c_qj

    ELSEIF ( present(r_qj) .AND. present(r_nh_semi_impl_term) ) THEN

       DO i = 1,n_vars

          qj(i) = dcmplx( r_qj(i) )

       END DO

    ELSE

       WRITE(*,*) 'Constitutive, eval_fluxes: problem with arguments'
       stop

    END IF

    ! initialize and evaluate the forces terms
    forces_term(1:n_eqns) = dcmplx(0.d0,0.d0)

    IF (rheology_flag) THEN

       CALL phys_var(c_qj = qj)
    
       mod_vel = cdsqrt( u**2 + v**2 )

       ! Voellmy Salm rheology
       IF ( rheology_model .EQ. 1 ) THEN

          IF ( mod_vel .GT. 0.d0 ) THEN

             forces_term(2) = forces_term(2) - rho_m * ( u / mod_vel ) *        &
                  mu * h * ( - grav * grav3_surf )
             
             forces_term(3) = forces_term(3) - rho_m * ( v / mod_vel ) *        &
                  mu * h * ( - grav * grav3_surf )
             
          END IF

          ! Plastic rheology
       ELSEIF ( rheology_model .EQ. 2 ) THEN
          

       ! Temperature dependent rheology
       ELSEIF ( rheology_model .EQ. 3 ) THEN

          
       ! Lahars rheology (O'Brien 1993, FLO2D)
       ELSEIF ( rheology_model .EQ. 4 ) THEN

          h_threshold = 1.d-20

          ! THIS IS ONLY USED WHEN SEDIMENT PERCENTAGE IS GIVEN AS INPUT AND
          ! IT IS CONSTANT WITHIN THE SIMULATION
          ! alphas = DCMPLX( sed_vol_perc*1.D-2 , 0.D0 )

          ! Yield strength
          tau_y = alpha2 * cdexp( beta2 * alphas )

          IF ( h .GT. h_threshold ) THEN
             
             ! Yield slope component
             s_y = tau_y / h
                       
          ELSE
             
             ! Yield slope component
             s_y = tau_y / h_threshold
                          
          END IF
  
          IF ( mod_vel .GT. 0.d0 ) THEN

             forces_term(2) = forces_term(2) - grav * h * ( u / mod_vel ) * s_y

             forces_term(3) = forces_term(3) - grav * h * ( v / mod_vel ) * s_y
  
          END IF

       ELSEIF ( rheology_model .EQ. 5 ) THEN

          
       ENDIF
              
    ENDIF

    IF ( temperature_flag ) THEN


    END IF

    
    IF ( present(c_qj) .AND. present(c_nh_semi_impl_term) ) THEN

       c_nh_semi_impl_term = forces_term

    ELSEIF ( present(r_qj) .AND. present(r_nh_semi_impl_term) ) THEN

       r_nh_semi_impl_term = dble( forces_term )

    END IF 

  END SUBROUTINE eval_nh_semi_impl_terms

  
  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_expl_terms( Bprimej_x , Bprimej_y , source_xy , qj ,          &
       expl_term )

    USE parameters_2d, ONLY : source_flag , vel_source , t_source
    
    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bprimej_x
    REAL*8, INTENT(IN) :: Bprimej_y
    REAL*8, INTENT(IN) :: source_xy
    
    REAL*8, INTENT(IN) :: qj(n_eqns)
    REAL*8, INTENT(OUT) :: expl_term(n_eqns)

    REAL*8 :: visc_heat_coeff

    expl_term(1:n_eqns) = 0.d0

    CALL phys_var(r_qj = qj)

    IF ( source_flag ) expl_term(1) = - source_xy * vel_source
    
    expl_term(2) = grav * dble(rho_m) * dble(h) * bprimej_x
   
    expl_term(3) = grav * dble(rho_m) * dble(h) * bprimej_y

    IF ( temperature_flag .AND. source_flag ) THEN

       IF ( solid_transport_flag ) THEN
    
          expl_term(5) = - source_xy * vel_source * t_source

       ELSE

          expl_term(4) = - source_xy * vel_source * t_source

       END IF

       IF ( rheology_model .EQ. 3 ) THEN
              
          IF ( dble(h) .GT. 0.d0 ) THEN
             
             ! Equation 10 from Costa & Macedonio, 2005
             visc_heat_coeff = emme * nu_ref / ( c_p * dble(h) ) 
             
          ELSE
             
             visc_heat_coeff = 0.d0
             
          END IF

          IF ( solid_transport_flag ) THEN
                          
             ! Viscous heating
             ! Last R.H.S. term in equation 4 from Costa & Macedonio, 2005
             expl_term(5) = expl_term(5) - visc_heat_coeff * ( dble(u)**2       &
                  + dble(v)**2 ) * dexp( - visc_par * ( dble(t) - t_ref ) ) 
             
          ELSE

             ! Viscous heating
             ! Last R.H.S. term in equation 4 from Costa & Macedonio, 2005
             expl_term(4) = expl_term(4) - visc_heat_coeff * ( dble(u)**2       &
                  + dble(v)**2 ) * dexp( - visc_par * ( dble(t) - t_ref ) ) 
             
          END IF

       END IF
          
    END IF
           
  END SUBROUTINE eval_expl_terms


  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_erosion_dep_term( qj , erosion_term , deposition_term )
    
    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: qj(n_eqns)
    
    REAL*8, INTENT(OUT) :: erosion_term
    REAL*8, INTENT(OUT) :: deposition_term
    
    REAL*8 :: r_alphas


    !REAL*8 :: settling_single_vel

    !REAL*8 :: RZ_exp

    !REAL*8 :: Rey

    !REAL*8 :: kin_visc

    !REAL*8 :: diam

    !REAL*8 :: R

    REAL*8 :: mod_vel

    !diam = 1.D-3

    CALL phys_var(r_qj = qj)

    deposition_term = 0.d0
    erosion_term = 0.d0
    
    r_alphas = dble(alphas)

    !Rey = DSQRT( grav * diam**3 ) / kin_visc

    !RZ_exp = ( 4.7D0 + 0.41D0 * Rey**0.75D0 ) / ( 1.D0 + 0.175D0 * Rey**0.75D0 )

    !settling_single_vel = kin_visc / diam * 

    !settling_vel = ( 1.D0 - r_alphas )**(2.7D0 - 0.15D0*Rz_exp )

    deposition_term = r_alphas * settling_vel

    mod_vel = dsqrt( dble(u)**2 + dble(v)**2 )
  
    IF ( dble(h) .GT. 1.d-2) THEN
    
       !erosion_term = erosion_coeff * mod_vel * DBLE(h) * ( rho_s - DBLE(rho_m) )  &
       !     / ( rho_s - rho_w )
       
       erosion_term = erosion_coeff * mod_vel * dble(h) * ( 1.d0 - r_alphas )

    ELSE

       erosion_term = 0.d0

    END IF
    
    !WRITE(*,*) 'r_alphas, settling_vel',r_alphas, settling_vel
    !WRITE(*,*) 'eval_erosion_dep_term',erosion_term,deposition_term
    !READ(*,*)

    RETURN
  
  END SUBROUTINE eval_erosion_dep_term


  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_topo_term( qj , deposition_avg_term , erosion_avg_term ,      &
       eqns_term, topo_term )
    
    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: qj(n_eqns)
    REAL*8, INTENT(IN) :: deposition_avg_term
    REAL*8, INTENT(IN) :: erosion_avg_term

    REAL*8, INTENT(OUT):: eqns_term(n_eqns)
    REAL*8, INTENT(OUT):: topo_term
    

    CALL phys_var(r_qj = qj)

    eqns_term(1:n_eqns) = 0.d0

    ! free surface (topography+flow) equation
    eqns_term(1) = erosion_avg_term - deposition_avg_term

    ! x-momenutm equation
    eqns_term(2) = - dble(u) * rho_s * deposition_avg_term

    ! y-momentum equation
    eqns_term(3) = - dble(v) * rho_s * deposition_avg_term

    ! solid phase thickness equation
    eqns_term(4) = erosion_avg_term - deposition_avg_term
  
    ! topography term
    topo_term = - erosion_avg_term + deposition_avg_term

  END SUBROUTINE eval_topo_term


END MODULE constitutive_2d