constitutive.f90

Go to the documentation of this file.

!********************************************************************************
!********************************************************************************
MODULE constitutive

  USE parameters, ONLY : n_eqns , n_vars

  IMPLICIT none

  LOGICAL, ALLOCATABLE :: implicit_flag(:)

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

  !--------- Constants for the equations of state -----------------------------


  COMPLEX*16 :: h
  COMPLEX*16 :: u


  REAL*8 :: grav

CONTAINS

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

  SUBROUTINE init_problem_param

    USE parameters, ONLY : n_nh

    ALLOCATE( implicit_flag(n_eqns) )

    implicit_flag(1:n_eqns) = .false.
    implicit_flag(2) = .true.

    n_nh = count( implicit_flag )

  END SUBROUTINE init_problem_param

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

  SUBROUTINE phys_var(Bj,r_qj,c_qj)
    
    USE complexify
    USE parameters, 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 )

    
    ! Correction for small values of h (Eq. 2.17 K&P 2007)
    IF ( REAL( h ) .GT. eps_sing ** 0.25d0 ) then

       u = qj(2) / h 

    ELSE

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

    END IF

  END SUBROUTINE phys_var

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

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

    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
  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_local_speeds2(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) = REAL(u_temp) - dsqrt( grav * REAL(h_temp) ) 
    vel_max(1:n_eqns) = REAL(u_temp) + dsqrt( grav * REAL(h_temp) )


  END SUBROUTINE eval_local_speeds2



  !******************************************************************************
  !
  !******************************************************************************
  
  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)

    
  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      !> batimetry + height 
    REAL*8 :: r_u       !> velocity
    
    r_hb = qp(1)
    r_u = qp(2)
    
    qc(1) = r_hb
    qc(2) = ( r_hb - b ) * r_u

  END SUBROUTINE qp_to_qc

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE qc_to_qp2(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)
    qp(2) = REAL(u)

    
  END SUBROUTINE qc_to_qp2


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

  SUBROUTINE qp2_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      !> batimetry + height 
    REAL*8 :: r_u       !> velocity
    
    r_h = qp(1)
    r_u = qp(2)
    
    qc(1) = r_h + b
    qc(2) = r_h * r_u

  END SUBROUTINE qp2_to_qc

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE eval_fluxes(Bj,c_qj,r_qj,c_flux,r_flux)
    
    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)
    
    COMPLEX*16 :: qj(n_vars)
    COMPLEX*16 :: flux(n_eqns)

    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
   
    flux(1) = qj(2)
    
    ! check if h is small
    IF((qj(1)-bj).GT.10.d0**(-8))THEN

      flux(2) = qj(2)**2.d0 / ( qj(1) - bj ) + 0.5d0 * grav * ( qj(1) - bj ) ** 2.d0 

    ELSE

      flux(2) = 0.d0

    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 , c_qj , c_nh_term_impl ,   &
       r_qj , r_nh_term_impl )

    USE complexify 
    IMPLICIT NONE

    REAL*8, INTENT(IN) :: bj
    REAL*8, INTENT(IN) :: bprimej

    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 :: heat , drag
    COMPLEX*16 :: vel_term

    COMPLEX*16 :: forces_term(n_eqns)

    INTEGER :: i
    
    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) 

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

    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_explicit_forces( Bj , Bprimej , qj , expl_forces_term )

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: bj
    REAL*8, INTENT(IN) :: bprimej

    REAL*8, INTENT(IN) :: qj(n_eqns)
    REAL*8, INTENT(OUT) :: expl_forces_term(n_eqns)

    expl_forces_term(1:n_eqns) = 0.d0

    CALL phys_var(bj,r_qj = qj)

    expl_forces_term(2) = grav * h * bprimej

    ! friction term
    expl_forces_term(2) = expl_forces_term(2) + ( 0.001 / (1.d0+10.d0*h) ) * u

  END SUBROUTINE eval_explicit_forces

END MODULE constitutive