constitutive_2d.f90

Go to the documentation of this file.

!********************************************************************************
!********************************************************************************
MODULE constitutive_2d

  USE variables, ONLY : wp
  USE parameters_2d, ONLY : tolh
  USE parameters_2d, ONLY : n_eqns , n_vars , c_d

  IMPLICIT none

  LOGICAL, ALLOCATABLE :: implicit_flag(:)

  REAL(wp) :: u_atm_umbl , v_atm_umbl
  REAL(wp) :: n
  REAL(wp) :: grav , drho_dz , rho_nbl

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.

    n_nh = count( implicit_flag )

    RETURN
    
  END SUBROUTINE init_problem_param

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

  SUBROUTINE r_phys_var(r_qj,r_h,r_u,r_v)

    USE parameters_2d, ONLY : eps_sing
    IMPLICIT none

    REAL(wp), INTENT(IN) :: r_qj(n_vars)
    REAL(wp), INTENT(OUT) :: r_h
    REAL(wp), INTENT(OUT) :: r_u
    REAL(wp), INTENT(OUT) :: r_v

    r_h = r_qj(1)

    ! velocity components
    IF ( r_qj(1) .GT. eps_sing ) THEN

       r_u = r_qj(2) / r_qj(1)
       r_v = r_qj(3) / r_qj(1)

    ELSE

       r_u = sqrt(2.0_wp) * r_qj(1) * r_qj(2) / sqrt( r_qj(1)**4 + eps_sing**4 )
       r_v = sqrt(2.0_wp) * r_qj(1) * r_qj(3) / sqrt( r_qj(1)**4 + eps_sing**4 )

    END IF

    RETURN

  END SUBROUTINE r_phys_var

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

  SUBROUTINE c_phys_var(c_qj,h,u,v)

    USE complexify
    USE parameters_2d, ONLY : eps_sing
    IMPLICIT none

    COMPLEX(wp), INTENT(IN) :: c_qj(n_vars)
    COMPLEX(wp), INTENT(OUT) :: h
    COMPLEX(wp), INTENT(OUT) :: u
    COMPLEX(wp), INTENT(OUT) :: v

    COMPLEX(wp) :: inv_cqj1

    IF ( REAL(c_qj(1)) .GT. eps_sing ) then

       inv_cqj1 = 1.0_wp / c_qj(1)

    ELSE

       inv_cqj1 = cmplx(0.0_wp,0.0_wp,wp)

    END IF

    h = c_qj(1)

    ! velocity components
    IF ( REAL( c_qj(1) ) .GT. eps_sing ) then

       u = c_qj(2) * inv_cqj1
       v = c_qj(3) * inv_cqj1

    ELSE

       u = sqrt(2.0_wp) * c_qj(1) * c_qj(2) / sqrt( c_qj(1)**4 + eps_sing**4 )
       v = sqrt(2.0_wp) * c_qj(1) * c_qj(3) / sqrt( c_qj(1)**4 + eps_sing**4 )

    END IF

    RETURN

  END SUBROUTINE c_phys_var

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

  SUBROUTINE mixt_var(qpj)

    IMPLICIT none

    REAL(wp), INTENT(IN) :: qpj(n_vars+2)
    REAL(wp) :: r_u
    REAL(wp) :: r_v
    REAL(wp) :: r_h

    r_h = qpj(1)

    IF ( qpj(1) .LE. 0.0_wp ) THEN

       r_u = 0.0_wp
       r_v = 0.0_wp

       RETURN

    END IF

    r_u = qpj(n_vars+1)
    r_v = qpj(n_vars+2)

    RETURN

  END SUBROUTINE mixt_var

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

  SUBROUTINE qc_to_qp(qc,qp)

    IMPLICIT none

    REAL(wp), INTENT(IN) :: qc(n_vars)
    REAL(wp), INTENT(OUT) :: qp(n_vars+2)

    REAL(wp) :: r_h
    REAL(wp) :: r_u
    REAL(wp) :: r_v

    CALL r_phys_var(qc,r_h,r_u,r_v)
    
    qp(1) = r_h
    qp(2) = r_h*r_u
    qp(3) = r_h*r_v
    
    qp(n_vars+1) = r_u
    qp(n_vars+2) = r_v
    
    RETURN

  END SUBROUTINE qc_to_qp

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

  SUBROUTINE qp_to_qc(qp,qc)

    IMPLICIT none

    REAL(wp), INTENT(IN) :: qp(n_vars+2)
    REAL(wp), INTENT(OUT) :: qc(n_vars)

    REAL(wp) :: r_u
    REAL(wp) :: r_v
    REAL(wp) :: r_h
    REAL(wp) :: r_hu
    REAL(wp) :: r_hv

    r_h = qp(1)

    IF ( r_h .GT. 0.0_wp ) THEN

       r_hu = qp(2)
       r_hv = qp(3)

       r_u = qp(n_vars+1)
       r_v = qp(n_vars+2)

    ELSE

       r_hu = 0.0_wp
       r_hv = 0.0_wp

       r_u = 0.0_wp
       r_v = 0.0_wp

       qc(1:n_vars) = 0.0_wp
       RETURN

    END IF

    qc(1) = r_h 
    qc(2) = r_hu
    qc(3) = r_hv

    RETURN

  END SUBROUTINE qp_to_qc

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

  SUBROUTINE qp_to_qp2(qpj,qp2j)

    IMPLICIT none

    REAL(wp), INTENT(IN) :: qpj(n_vars)
    REAL(wp), INTENT(OUT) :: qp2j(3)

    qp2j(1) = qpj(1)

    IF ( qpj(1) .LE. 0.0_wp ) THEN

       qp2j(2) = 0.0_wp
       qp2j(3) = 0.0_wp

    ELSE

       qp2j(2) = qpj(2)/qpj(1)
       qp2j(3) = qpj(3)/qpj(1)

    END IF

    RETURN

  END SUBROUTINE qp_to_qp2

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

  SUBROUTINE eval_local_speeds_x(qpj,vel_min,vel_max)

    IMPLICIT none

    REAL(wp), INTENT(IN) :: qpj(n_vars+2)

    REAL(wp), INTENT(OUT) :: vel_min(n_vars) , vel_max(n_vars)

    REAL(wp) :: r_h
    REAL(wp) :: r_u

    r_h = qpj(1)
    r_u = qpj(n_vars+1)

    vel_min(1:n_eqns) = r_u - 0.5_wp * n * r_h
    vel_max(1:n_eqns) = r_u + 0.5_wp * n * r_h 

    RETURN

  END SUBROUTINE eval_local_speeds_x

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

  SUBROUTINE eval_local_speeds_y(qpj,vel_min,vel_max)

    IMPLICIT none

    REAL(wp), INTENT(IN)  :: qpj(n_vars+2)
    REAL(wp), INTENT(OUT) :: vel_min(n_vars) , vel_max(n_vars)

    REAL(wp) :: r_h
    REAL(wp) :: r_v

    r_h = qpj(1)
    r_v = qpj(n_vars+2)

    vel_min(1:n_eqns) = r_v - 0.5_wp * n * r_h
    vel_max(1:n_eqns) = r_v + 0.5_wp * n * r_h

    RETURN

  END SUBROUTINE eval_local_speeds_y

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

  SUBROUTINE eval_fluxes(qcj,qpj,dir,flux)

    IMPLICIT none

    REAL(wp), INTENT(IN) :: qcj(n_vars)
    REAL(wp), INTENT(IN) :: qpj(n_vars+2)
    INTEGER, INTENT(IN) :: dir

    REAL(wp), INTENT(OUT) :: flux(n_eqns)

    REAL(wp) :: r_h
    REAL(wp) :: r_u
    REAL(wp) :: r_v

    pos_thick:IF ( qcj(1) .GT. 0.0_wp ) THEN

       r_h = qpj(1)
       r_u = qpj(n_vars+1)
       r_v = qpj(n_vars+2)

       IF ( dir .EQ. 1 ) THEN

          ! Volume flux in x-direction: u * h
          flux(1) = r_u * qcj(1)

          ! x-momentum flux in x-direction + hydrostatic pressure term
          flux(2) = r_u * qcj(2) + n**2 * r_h**3 / 12.0_wp

          ! y-momentum flux in x-direction: u * (h * v)
          flux(3) = r_u * qcj(3)

       ELSEIF ( dir .EQ. 2 ) THEN

          ! Volume flux in y-direction: v*h
          flux(1) = r_v * qcj(1)

          ! y-momentum flux in y-direction: v * (h * v)
          flux(2) = r_v * qcj(2)

          ! y-momentum flux in y-direction + hydrostatic pressure term
          flux(3) = r_v * qcj(3) + n**2 * r_h**3 / 12.0_wp

       END IF

    ELSE

       flux(1:n_eqns) = 0.0_wp

    ENDIF pos_thick

    RETURN

  END SUBROUTINE eval_fluxes

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

  SUBROUTINE eval_nonhyperbolic_terms( cell_fract_jk , dx_rel_jk , dy_rel_jk ,  &
       c_qj , c_nh_term_impl , r_qj , r_nh_term_impl )

    USE parameters_2d, ONLY: r_source , u_source, v_source, vol_flux_source

    USE complexify 

    IMPLICIT NONE

    COMPLEX(wp), INTENT(IN), OPTIONAL :: c_qj(n_vars)
    COMPLEX(wp), INTENT(OUT), OPTIONAL :: c_nh_term_impl(n_eqns)
    REAL(wp), INTENT(IN), OPTIONAL :: r_qj(n_vars)
    REAL(wp), INTENT(OUT), OPTIONAL :: r_nh_term_impl(n_eqns)

    REAL(wp), INTENT(IN) :: cell_fract_jk
    REAL(wp), INTENT(IN) :: dx_rel_jk
    REAL(wp), INTENT(IN) :: dy_rel_jk

    COMPLEX(wp) :: h
    COMPLEX(wp) :: u
    COMPLEX(wp) :: v
 
    COMPLEX(wp) :: qj(n_vars)
    COMPLEX(wp) :: forces_term(n_eqns)

    COMPLEX(wp) :: mod_vel , delta_u , delta_v

    REAL(wp) :: pi_g

    REAL(wp) :: h_dot

    INTEGER :: i

    pi_g = 4.0_wp * atan(1.0_wp)

    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) = cmplx( r_qj(i),0.0_wp,wp )

       END DO

    ELSE

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

    END IF

    ! initialize and evaluate the forces terms
    forces_term(1:n_eqns) = cmplx(0.0_wp,0.0_wp,wp)

    CALL c_phys_var(qj,h,u,v)

    delta_u = u - u_atm_umbl
    delta_v = v - v_atm_umbl

    mod_vel = sqrt( delta_u**2 + delta_v**2 )

    forces_term(2) = forces_term(2) - c_d * delta_u * mod_vel 
    
    forces_term(3) = forces_term(3) - c_d * delta_v * mod_vel 

    h_dot = cell_fract_jk * vol_flux_source / (pi_g * r_source**2)
    
    forces_term(1) = h_dot
    
    forces_term(2) = forces_term(2) + ( u_source + dx_rel_jk * h_dot ) * h_dot
    
    forces_term(3) = forces_term(3) + ( v_source + dy_rel_jk * h_dot ) * h_dot
    
    IF ( present(c_qj) .AND. present(c_nh_term_impl) ) THEN

       c_nh_term_impl = forces_term

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

       r_nh_term_impl = REAL( forces_term )

    END IF

    RETURN

  END SUBROUTINE eval_nonhyperbolic_terms

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

  SUBROUTINE eval_source_bdry( time, vect_x , vect_y , source_bdry )

    USE parameters_2d, ONLY : h_source , vel_source , time_param

    IMPLICIT NONE

    REAL(wp), INTENT(IN) :: time
    REAL(wp), INTENT(IN) :: vect_x
    REAL(wp), INTENT(IN) :: vect_y
    REAL(wp), INTENT(OUT) :: source_bdry(n_vars)

    REAL(wp) :: t_rem
    REAL(wp) :: t_coeff
    REAL(wp) :: pi_g

    IF ( time .GE. time_param(4) ) THEN

       ! The exponents of t_coeff are such that Ri does not depend on t_coeff
       source_bdry(1) = 0.0_wp
       source_bdry(2) = 0.0_wp
       source_bdry(3) = 0.0_wp
       source_bdry(n_vars+1) = 0.0_wp
       source_bdry(n_vars+2) = 0.0_wp

       RETURN

    END IF

    t_rem = mod( time , time_param(1) )

    pi_g = 4.0_wp * atan(1.0_wp) 

    t_coeff = 0.0_wp

    IF ( time_param(3) .EQ. 0.0_wp ) THEN

       IF ( t_rem .LE. time_param(2) ) t_coeff = 1.0_wp

    ELSE

       IF ( t_rem .LT. time_param(3) ) THEN

          t_coeff = 0.5_wp * ( 1.0_wp - cos( pi_g * t_rem / time_param(3) ) )

       ELSEIF ( t_rem .LE. ( time_param(2) - time_param(3) ) ) THEN

          t_coeff = 1.0_wp

       ELSEIF ( t_rem .LE. time_param(2) ) THEN

          t_coeff = 0.5_wp * ( 1.0_wp + cos( pi_g * ( ( t_rem - time_param(2) ) &
               / time_param(3) + 1.0_wp ) ) )

       END IF

    END IF

    ! The exponents of t_coeff are such that Ri does not depend on t_coeff
    source_bdry(1) = t_coeff * h_source
    source_bdry(2) = t_coeff**1.5_wp * h_source * vel_source * vect_x
    source_bdry(3) = t_coeff**1.5_wp * h_source * vel_source * vect_y
    source_bdry(n_vars+1) = t_coeff**0.5_wp * vel_source * vect_x
    source_bdry(n_vars+2) = t_coeff**0.5_wp * vel_source * vect_y 

    RETURN

  END SUBROUTINE eval_source_bdry

END MODULE constitutive_2d