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

  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 :: xs

  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 :: gamma_w

  REAL*8 :: gamma_s

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(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

    ELSE

       qj = dcmplx(r_qj)

    END IF

    h = qj(1) - dcmplx( bj , 0.d0 )

    IF ( REAL( h ) .GT. eps_sing ** 0.25D0 ) then

       u = qj(2) / h 

       v = qj(3) / h

    ELSE

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

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

    END IF

    IF ( solid_transport_flag ) THEN

       IF ( REAL( h ) .GT. 1.D-25 ) then
          
          xs = qj(4) / h
          
          IF ( temperature_flag ) t = qj(5) / h
          
       ELSE
          
          xs = 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

    ELSE

       IF ( temperature_flag ) THEN

          IF ( REAL( 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
     
    END IF
    
  END SUBROUTINE phys_var

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

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

    CALL phys_var(bj,r_qj = qj)
    
    vel_min(1:n_eqns) = REAL(u) - DSQRT( grav * REAL(h) )
    vel_max(1:n_eqns) = REAL(u) + DSQRT( grav * REAL(h) )
        
  END SUBROUTINE eval_local_speeds_x

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

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

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

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

    REAL*8 :: h_temp , u_temp

    h_temp = qj(1) - bj

    IF ( h_temp .NE. 0.d0 ) THEN

       u_temp = qj(2) / h_temp

    ELSE

       u_temp = 0.d0

    END IF

    vel_min(1:n_eqns) = u_temp - dsqrt( grav * h_temp )
    vel_max(1:n_eqns) = u_temp + dsqrt( grav * h_temp )

  END SUBROUTINE eval_local_speeds2_x

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

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

    REAL*8 :: h_temp , v_temp

    h_temp = qj(1) - bj

    IF ( h_temp .NE. 0.d0 ) THEN

       v_temp = qj(3) / h_temp

    ELSE

       v_temp = 0.d0

    END IF

    vel_min(1:n_eqns) = v_temp - dsqrt( grav * h_temp )
    vel_max(1:n_eqns) = v_temp + dsqrt( grav * h_temp )
    
  END SUBROUTINE eval_local_speeds2_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(b,r_qj = qc)
        
    qp(1) = REAL(h+b)
    qp(2) = REAL(u)
    qp(3) = REAL(v)

    IF ( solid_transport_flag ) THEN
       
       qp(4) = REAL(xs)
       
       IF ( temperature_flag ) qp(5) = REAL(t)
       
    ELSE
       
       IF ( temperature_flag ) qp(4) = REAL(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_hB
    REAL*8 :: r_u
    REAL*8 :: r_v
    REAL*8 :: r_xs
    REAL*8 :: r_T

    r_hb = qp(1)
    r_u  = qp(2)
    r_v  = qp(3)

    qc(1) = r_hb
    qc(2) = ( r_hb - b ) * r_u
    qc(3) = ( r_hb - b ) * r_v

    IF ( solid_transport_flag ) THEN
       
       r_xs = qp(4)
       qc(4) = ( r_hb - b ) * r_xs 
       
       IF ( temperature_flag ) THEN
          
          r_t  = qp(5)
          qc(5) = ( r_hb - b ) * r_t 
       
       END IF

    ELSE

       IF ( temperature_flag ) THEN
          
          r_t  = qp(4)
          qc(4) = ( r_hb - b ) * r_t 
       
       END IF

    END IF

  END SUBROUTINE qp_to_qc

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

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

    REAL*8 :: r_hB
    REAL*8 :: r_u
    REAL*8 :: r_v
    REAL*8 :: r_xs
    REAL*8 :: r_T

    ! Desingularization
    CALL phys_var(b,r_qj = qrec)
          
    r_hb = REAL(h) + B
    r_u = REAL(u)
    r_v = REAL(v)

    qc(1) = r_hb
    qc(2) = REAL(h) * r_u
    qc(3) = REAL(h) * r_v
    
    IF ( solid_transport_flag ) THEN
       
       r_xs = REAL(xs)
       qc(4) = REAL(h) * r_xs
       
       IF ( temperature_flag ) THEN
          
          r_t = REAL(t)
          qc(5) = REAL(h) * r_T 
          
       END IF
       
    ELSE
       
       IF ( temperature_flag ) THEN
          
          r_t = REAL(t)
          qc(4) = REAL(h) * r_T 
          
       END IF
       
    END IF
    
  END SUBROUTINE qrec_to_qc


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

    REAL*8, INTENT(IN) :: Bj
    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

    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


    IF ( dir .EQ. 1 ) THEN

    ! flux F (derivated wrt x in the equations)

       flux(1) = qj(2)

       h_temp = qj(1) - bj

       IF ( REAL(h_temp) .NE. 0.D0 ) then

          u_temp = qj(2) / h_temp
          
          flux(2) = h_temp * u_temp**2 + 0.5d0 * grav * h_temp**2  

          flux(3) = u_temp * qj(3)

          IF ( solid_transport_flag ) THEN

             flux(4) = u_temp * qj(4)

             ! 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

       ELSE

          flux(2:n_eqns) = 0.d0

       ENDIF

    ELSEIF ( dir .EQ. 2 ) THEN

       ! flux G (derivated wrt y in the equations)

       flux(1) = qj(3)

       h_temp = qj(1) - bj

       IF(REAL(h_temp).NE.0.d0)then

          v_temp = qj(3) / h_temp

          flux(2) = v_temp * qj(2)
          
          flux(3) = h_temp * v_temp**2 + 0.5d0 * grav * h_temp**2

          IF ( solid_transport_flag ) THEN

             flux(4) = v_temp * qj(4)

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

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

          END IF

       ELSE

          flux(2:n_eqns) = 0.d0

       ENDIF

    ELSE

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

       stop

    ENDIF

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

       c_flux = flux

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

       r_flux = REAL( flux )

    END IF

  END SUBROUTINE eval_fluxes

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

  SUBROUTINE eval_nonhyperbolic_terms( Bj , Bprimej_x , Bprimej_y , grav3_surf ,&
       c_qj , c_nh_term_impl , r_qj , r_nh_term_impl )

    USE complexify 
    USE parameters_2d, ONLY : sed_vol_perc

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bj
    REAL*8, INTENT(IN) :: Bprimej_x, Bprimej_y
    REAL*8, INTENT(IN) :: grav3_surf

    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 :: tau_y

    COMPLEX*8 :: fluid_visc

    COMPLEX*8 :: sed_vol_fract_cmplx

    COMPLEX*8 :: gamma_m

    COMPLEX*8 :: s_f

    COMPLEX*8 :: s_y

    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(bj,c_qj = qj)
    
       mod_vel = cdsqrt( u**2 + v**2 )
       
       ! Voellmy Salm rheology
       IF ( rheology_model .EQ. 1 ) THEN
       
          IF ( REAL(mod_vel) .NE. 0.D0 ) then 
          
             ! IMPORTANT: grav3_surv is always negative 
             forces_term(2) = forces_term(2) -  ( u / mod_vel ) *               &
                  ( grav / xi ) * mod_vel ** 2
          
             forces_term(3) = forces_term(3) -  ( v / mod_vel ) *               &
                  ( grav / xi ) * mod_vel ** 2
          
          ENDIF
        
       ! Plastic rheology
       ELSEIF ( rheology_model .EQ. 2 ) THEN
       
          IF ( REAL(mod_vel) .NE. 0.D0 ) then 
          
             forces_term(2) = forces_term(2) - tau * (u/mod_vel)
          
             forces_term(3) = forces_term(3) - tau * (v/mod_vel)

          ENDIF

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

          IF ( REAL(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 ( REAL(mod_vel) .NE. 0.D0 ) then 
          
             ! Last R.H.S. term in equation 2 from Costa & Macedonio, 2005
             forces_term(2) = forces_term(2) - gamma * u
          
             ! Last R.H.S. term in equation 3 from Costa & Macedonio, 2005
             forces_term(3) = forces_term(3) - gamma * v

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

          h_threshold = 1.d-20

          sed_vol_fract_cmplx = dcmplx(sed_vol_perc/100.d0,0.d0)

          ! Convert from mass fraction to volume fraction
          ! sed_vol_fract_cmplx = xs * gamma_w / ( xs * gamma_w +               &
          !                 (  DCMPLX(1.D0,0.D0) - xs ) * gamma_s )

          !IF ( xs .NE. 0.D0 ) THEN
             
             !WRITE(*,*) 'xs',xs
             !WRITE(*,*) 'sed_vol_fract',DBLE(sed_vol_fract_cmplx)
             !READ(*,*) 
             
          !END IF


          ! Mixture density
          gamma_m = ( dcmplx(1.d0,0.d0) - sed_vol_fract_cmplx ) * gamma_w       &
               + sed_vol_fract_cmplx * gamma_s 

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

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


          IF ( h .GT. h_threshold ) THEN
             
             ! Yield slope component
             s_y = tau_y / ( gamma_m * h )
             
             ! Viscous slope component
             s_v = kappa * fluid_visc * mod_vel / ( 8.d0 * gamma_m * h**2 )
             
             
             ! Turbulent dispersive component
             s_td = n_td**2 * mod_vel**2 / ( h**(4.d0/3.d0) )
          
             ! WRITE(*,*) 's_terms',REAL(s_y),REAL(s_v),REAL(s_td)
             ! WRITE(*,*) ' u', REAL(u)

          ELSE
             
             ! Yield slope component
             s_y = tau_y / ( gamma_m * h_threshold )
             
             ! Viscous slope component
             s_v = kappa * fluid_visc * mod_vel / ( 8.d0 * gamma_m *            &
                  h_threshold**2 )
             
             ! Turbulent dispersive components
             s_td = 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


          !WRITE(*,*) 's_terms',DBLE(s_y), &
          !     DBLE(Kappa * fluid_visc / ( 8.D0 * gamma_m * h**2 )), &
          !     DBLE(n_td**2 / ( h**(4.D0/3.D0) ))

          !WRITE(*,*) 'grav*h',grav*h

          !WRITE(*,*) 'eval_nh',DBLE(u),DBLE(forces_term(2))

       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) - tau * ( u / mod_vel )
             forces_term(3) = forces_term(3) - tau * ( v / mod_vel )

          END IF
          
       ENDIF
              
    ENDIF

    IF ( temperature_flag ) THEN

       CALL phys_var(bj,c_qj = qj)

       IF ( REAL(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 ( REAL(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 ( REAL(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 ( REAL(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 ( REAL(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 ( REAL(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 = REAL( nh_term )

    END IF 

  END SUBROUTINE eval_nonhyperbolic_terms

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

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

    USE complexify 
    USE parameters_2d, ONLY : sed_vol_perc

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: Bj
    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 :: sed_vol_fract_cmplx

    COMPLEX*8 :: gamma_m

    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(bj,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) -  ( u / mod_vel ) *               &
                  mu * h * ( - grav * grav3_surf )
             
             forces_term(3) = forces_term(3) -  ( 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

          sed_vol_fract_cmplx = dcmplx( sed_vol_perc*1.d-2 , 0.d0 )

          ! Convert from mass fraction to volume fraction
          ! sed_vol_fract_cmplx = xs * gamma_w / ( xs * gamma_w +               &
          !                 (  DCMPLX(1.D0,0.D0) - xs ) * gamma_s )

          !IF ( xs .NE. 0.D0 ) THEN
             
             !WRITE(*,*) 'xs',xs
             !WRITE(*,*) 'sed_vol_fract',DBLE(sed_vol_fract_cmplx)
             !READ(*,*) 
             
          !END IF


          ! Mixture density
          gamma_m = ( dcmplx(1.d0,0.d0) - sed_vol_fract_cmplx ) * gamma_w       &
               + sed_vol_fract_cmplx * gamma_s 

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

          IF ( h .GT. h_threshold ) THEN
             
             ! Yield slope component
             s_y = tau_y / ( gamma_m * h )
                       
          ELSE
             
             ! Yield slope component
             s_y = tau_y / ( gamma_m * 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( Bj, 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) :: Bj
    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(bj,r_qj = qj)

    IF ( source_flag ) expl_term(1) = - source_xy * vel_source
    
    expl_term(2) = grav * REAL(h) * Bprimej_x
   
    expl_term(3) = grav * REAL(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 ( REAL(h) .GT. 0.D0 ) then
             
             ! Equation 10 from Costa & Macedonio, 2005
             visc_heat_coeff = emme * nu_ref / ( c_p * REAL(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 * ( REAL(u)**2          &
                  + REAL(v)**2 ) * dexp( - visc_par * ( REAL(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 * ( REAL(u)**2          &
                  + REAL(v)**2 ) * dexp( - visc_par * ( REAL(T) - t_ref ) ) 
             
          END IF

       END IF
          
    END IF
           
  END SUBROUTINE eval_expl_terms

END MODULE constitutive_2d