equations.f90

Go to the documentation of this file.

!********************************************************************************
!********************************************************************************
MODULE equations

  USE constitutive
  
  USE geometry, ONLY : pi
  USE parameters, ONLY : verbose_level
  USE parameters, ONLY : n_eqns , n_vars
  USE parameters, ONLY : n_cry , n_gas , n_mom

  USE parameters, ONLY : idx_p1 , idx_p2 , idx_u1 , idx_u2 , idx_t ,            &
       idx_xd_first , idx_xd_last , idx_alfa_first , idx_alfa_last ,            &
       idx_beta_first , idx_beta_last
  
  USE parameters, ONLY : idx_mix_mass_eqn , idx_vol1_eqn , idx_mix_mom_eqn ,    &
       idx_rel_vel_eqn , idx_mix_engy_eqn , idx_dis_gas_eqn_first ,             &
       idx_dis_gas_eqn_last , idx_ex_gas_eqn_first , idx_ex_gas_eqn_last ,      &
       idx_cry_eqn_first , idx_cry_eqn_last
  
  IMPLICIT none

  LOGICAL :: isothermal

  REAL*8 :: fixed_temp

  LOGICAL :: lateral_degassing_flag

  LOGICAL :: lateral_degassing


  REAL*8 :: alfa2_lat_thr

  LOGICAL :: ext_water

  LOGICAL :: inst_vaporization

  CHARACTER*30 :: aquifer_type

  REAL*8 :: total_water_influx

  REAL*8 :: min_z_influx

  REAL*8 :: delta_z_influx

  REAL*8 :: max_z_influx

  REAL*8 :: t_w

  COMPLEX*16 :: t_boiling

  REAL*8 :: lambda_w

  REAL*8 :: cv_w

  REAL*8 :: visc_w

CONTAINS


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

  SUBROUTINE phys_var_qp(r_qp,c_qp)

    USE complexify 
    IMPLICIT none

    REAL*8, INTENT(IN), OPTIONAL :: r_qp(n_eqns)
    COMPLEX*16, INTENT(IN), OPTIONAL :: c_qp(n_eqns)

    COMPLEX*16 :: qp(n_eqns)

    INTEGER :: i , j

    IF ( present(c_qp) ) THEN

       qp = c_qp

    ELSE

       qp = dcmplx( r_qp , 0.d0 )

    END IF

    p_1 = qp(idx_p1)
    p_2 = qp(idx_p2)
    u_1 = qp(idx_u1)
    u_2 = qp(idx_u2)
    t = qp(idx_t)

    ! dissolved and axsolved gas variables
    x_d_md(1:n_gas) = qp(idx_xd_first:idx_xd_last)

    alfa_g(1:n_gas) = qp(idx_alfa_first:idx_alfa_last)

    alfa_2 = sum( alfa_g(1:n_gas) )

    alfa_g_2(1:n_gas) = alfa_g(1:n_gas) / alfa_2

    alfa_1 = 1.d0 - alfa_2

    ! crystal variables
    IF ( n_mom .LE. 1 ) THEN
       
       beta(1:n_cry) = qp(idx_beta_first:idx_beta_last)

    ELSE

       DO i = 1,n_cry
          
          DO j = 0,n_mom-1
             
             mom_cry(i,j) = qp(idx_cry_eqn_first+n_mom*(i-1)+j) 

          END DO

          beta(i) = pi / 6.d0 * mom_cry(i,3)
          
       END DO

    END IF
       
    ! eval_densities requires: beta, x_d_md , p_1 , p_2 , T
    CALL eval_densities

    alfarho_2 = alfa_2 * rho_2

    x_c_1(1:n_cry) = rho_c(1:n_cry) * beta(1:n_cry) / rho_1
    x_md_1 = dcmplx(1.d0,0.d0) - sum(x_c_1)

    x_d_1(1:n_gas) = x_d_md(1:n_gas) * x_md_1
    x_m_1 = dcmplx(1.d0,0.d0) - sum(x_d_1) - sum(x_c_1)

    alfa_d_1(1:n_gas) = rho_1 * x_d_1(1:n_gas) / rho_d(1:n_gas) 
    alfa_m_1 = rho_1 * x_m_1 / rho_m

    rho_mix = alfa_1 * rho_1 + alfa_2 * rho_2

    x_1 = alfa_1 * rho_1 / rho_mix
    x_2 = alfa_2 * rho_2 / rho_mix

    x_g(1:n_gas) = alfa_g(1:n_gas) * rho_g(1:n_gas) / rho_mix

    x_g_2(1:n_gas) = x_g(1:n_gas) / x_2

    x_d(1:n_gas) = x_d_1(1:n_gas) * x_1
    x_m = x_m_1 * x_1
    x_c(1:n_cry) = x_c_1(1:n_cry) * x_1

    u_mix = x_1 * u_1 + x_2 * u_2

    rhob_m = rho_mix * x_m
    rhob_c(1:n_cry) = rho_mix * x_c(1:n_cry)

    cv_mix = x_m * cv_m + sum( x_c(1:n_cry) * cv_c(1:n_cry) )                   &
         + sum( x_d(1:n_gas) * cv_d(1:n_gas) )                                  &
         + sum( x_g(1:n_gas) * cv_g(1:n_gas) )

  END SUBROUTINE phys_var_qp

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

  SUBROUTINE eval_qp2( r_qp , qp2 )

    IMPLICIT none

    REAL*8, INTENT(IN) :: r_qp(n_vars)
    REAL*8, INTENT(OUT) :: qp2(1+n_cry+n_gas+n_gas+4)

    COMPLEX*16 :: qp(n_vars)

    INTEGER :: i

    DO i = 1,n_vars

       qp(i) = dcmplx( r_qp(i), 0.d0 )

    END DO

    CALL phys_var_qp( c_qp = qp )
    CALL eos

    CALL f_beta_eq

    CALL f_xdis_eq

    qp2(1) = REAL(rho_1)
    qp2(1+1:1+n_gas) = REAL(rho_g(1:n_gas))
    qp2(1+n_gas+1:1+n_gas+n_cry) = REAL(beta_eq(1:n_cry))
    qp2(1+n_gas+n_cry+1:1+n_gas+n_cry+n_gas) = REAL( x_d_md_eq(1:n_gas) )

    CALL mixture_viscosity

    qp2(1+n_gas+n_cry+n_gas+1) = REAL( visc_mix )

    CALL f_viscmelt

    qp2(1+n_gas+n_cry+n_gas+2) = REAL( visc_melt )

    CALL f_theta

    qp2(1+n_gas+n_cry+n_gas+3) = REAL( theta )

    CALL f_bubbles

    qp2(1+n_gas+n_cry+n_gas+4) = REAL( visc_rel_bubbles )

  END SUBROUTINE eval_qp2

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

  SUBROUTINE eval_fluxes_qp(c_qp,r_qp,c_flux,r_flux)

    USE complexify
    USE geometry, ONLY : radius

    IMPLICIT none

    COMPLEX*16, INTENT(IN), OPTIONAL :: c_qp(n_vars)
    COMPLEX*16, INTENT(OUT), OPTIONAL :: c_flux(n_eqns)
    REAL*8, INTENT(IN), OPTIONAL :: r_qp(n_vars)
    REAL*8, INTENT(OUT), OPTIONAL :: r_flux(n_eqns)

    COMPLEX*16 :: qp(n_eqns)
    COMPLEX*16 :: flux(n_eqns)

    INTEGER :: i , j


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

       qp = c_qp

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

       qp = dcmplx( r_qp , 0.d0 )

    ELSE

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

    END IF

    CALL phys_var_qp( c_qp = qp )
    CALL eos

    flux(1:n_eqns) = dcmplx(0.d0,0.d0)

    !---- Mixture Density -------------------------------------------------------
    flux(idx_mix_mass_eqn) = rho_mix * u_mix * radius**2

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'Mixture mass flux'
       WRITE(*,*) flux(idx_mix_mass_eqn)

    END IF
       
    !---- Volumetric Fraction First Phase ---------------------------------------
    IF ( p_relax_model .EQ. 'single' ) THEN

       flux(idx_vol1_eqn) = dcmplx( 0.d0 , 0.d0 ) 

    ELSE

       flux(idx_vol1_eqn) = rho_mix * u_mix * alfa_1 * radius**2

    END IF

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'Volume Fraction1 flux'
       WRITE(*,*) flux(idx_vol1_eqn)

    END IF
    
    !---- Mixture Momentum ------------------------------------------------------
    flux(idx_mix_mom_eqn) = ( alfa_1 * rho_1 * u_1**2 + alfa_2 * rho_2 * u_2**2 &
         + alfa_1 * p_1 + alfa_2 * p_2 ) * radius**2

    IF ( verbose_level .GE. 3 ) THEN
       
       WRITE(*,*) 'Mixture momentum flux'
       WRITE(*,*) flux(idx_mix_mom_eqn)
       
    END IF
           
    !---- Relative Velocity -----------------------------------------------------
    IF ( drag_funct_model .EQ. 'single_velocity' ) THEN

       flux(idx_rel_vel_eqn) = dcmplx( 0.d0 , 0.d0 ) 

    ELSE

       flux(idx_rel_vel_eqn) = ( 0.5d0 * ( u_1 * u_1 - u_2 * u_2 )              &
            + ( mu_1 - mu_2 ) ) * radius**2

    END IF

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'Relative velocity flux'
       WRITE(*,*) flux(idx_rel_vel_eqn)

    END IF
       
    !---- Total Energy ----------------------------------------------------------
    IF ( isothermal ) THEN

       flux(idx_mix_engy_eqn) = dcmplx( 0.d0 , 0.d0 )

    ELSE

       flux(idx_mix_engy_eqn) = ( alfa_1 * rho_1 * u_1 * ( e_1 + p_1/rho_1      &
            + 0.5d0 * u_1 * u_1 ) + alfa_2 * rho_2 * u_2 * ( e_2 + p_2/rho_2    &
            + 0.5d0 * u_2* u_2 ) - rho_mix * x_1 * x_2 * ( u_1 - u_2 )          &
            * ( s_1 - s_2 ) * t ) * radius**2
    END IF

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'Mixture energy flux'
       WRITE(*,*) flux(idx_mix_engy_eqn)

    END IF

    !----- Dissolved Gas Phases -------------------------------------------------
    flux(idx_dis_gas_eqn_first:idx_dis_gas_eqn_last) =                          &
         ( x_d_md(1:n_gas)*alfa_1 *                                             &
         ( rho_1 - sum( beta(1:n_cry) * rho_c(1:n_cry) ) ) * u_1 ) * radius**2

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'Dissolved gas mass fraction fluxes'
       WRITE(*,*) flux(idx_ex_gas_eqn_first:idx_ex_gas_eqn_last)

    END IF

    !---- Mass Fraction Exsolved Gas Phases -------------------------------------
    flux(idx_ex_gas_eqn_first:idx_ex_gas_eqn_last) = alfa_g(1:n_gas) *          &
         rho_g(1:n_gas) * u_2 * radius**2

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'Mass fraction exsolved gas'
       WRITE(*,*) flux(idx_ex_gas_eqn_first:idx_ex_gas_eqn_last)

    END IF

    !----- Crystal Phases -------------------------------------------------------

    IF ( n_mom .LE. 1 ) THEN
       
       flux(idx_cry_eqn_first:idx_cry_eqn_last) = alfa_1 * rho_c(1:n_cry) *     &
            beta(1:n_cry) * u_1 * radius**2

       IF ( verbose_level .GE. 3 ) THEN
       
          WRITE(*,*) 'Crystal volume fraction'
          WRITE(*,*) flux(idx_cry_eqn_first:idx_cry_eqn_last)
          
       END IF

    ELSE

       DO i=1,n_cry

          DO j=0,n_mom-1

             flux(idx_cry_eqn_first+n_mom*(i-1)+j) = alfa_1 * mom_cry(i,j) *    &
                  u_1 * radius**2

          END DO

       END DO
       
    END IF
       
    IF ( present(c_qp) .AND. present(c_flux) ) THEN

       c_flux = flux

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

       r_flux = REAL( flux )

    END IF


  END SUBROUTINE eval_fluxes_qp

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

  SUBROUTINE eval_nonhyperbolic_terms_qp( c_qp , c_nh_term_impl , r_qp ,        &
       r_nh_term_impl )

    USE complexify 
    USE geometry, ONLY : radius
    IMPLICIT NONE

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

    COMPLEX*16 :: qp(n_eqns)

    COMPLEX*16 :: nh_term_impl(n_eqns)

    COMPLEX*16 :: relaxation_term(n_eqns)
    COMPLEX*16 :: forces_term(n_eqns)
    COMPLEX*16 :: source_term(n_eqns)

    INTEGER :: i


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

       qp = c_qp

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

       DO i = 1,n_eqns

          qp(i) = dcmplx( r_qp(i) , 0.d0 )

       END DO

    ELSE

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

    END IF

    CALL phys_var_qp( c_qp = qp )
    CALL eos

    !--------------- Evaluate the relaxation terms -----------------------------
    CALL eval_relaxation_terms( relaxation_term )

    !--------------- Evaluate the forces terms ---------------------------------
    CALL eval_forces_terms( forces_term )

    !--------------- Evaluate the forces terms ---------------------------------
    CALL eval_source_terms( source_term )

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*) 'relaxation_term'
       WRITE(*,*) relaxation_term / ( radius**2 )
    
       WRITE(*,*) 'forces_term'
       WRITE(*,*) forces_term / ( radius **2 )
    
       WRITE(*,*) 'source_term'
       WRITE(*,*) source_term / ( radius **2 )
    
    END IF

    nh_term_impl = relaxation_term + forces_term + source_term

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

       c_nh_term_impl = nh_term_impl

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

       r_nh_term_impl = REAL( nh_term_impl )

    END IF

  END SUBROUTINE eval_nonhyperbolic_terms_qp

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

  SUBROUTINE eval_relaxation_terms( relaxation_term )

    USE complexify 
    USE geometry, ONLY : radius
    IMPLICIT none

    COMPLEX*16, INTENT(OUT) :: relaxation_term(n_eqns)

    COMPLEX*16 :: pressure_relaxation 
    COMPLEX*16 :: velocity_relaxation
    !COMPLEX*16 :: frag_relaxation

    INTEGER :: i , j
    INTEGER :: idx

    idx = 0

    relaxation_term(1:n_eqns) = dcmplx(0.d0,0.d0)

    ! relaxation term for the mixture density -----------------------------------
    relaxation_term(idx_mix_mass_eqn) = dcmplx(0.d0,0.d0)


    ! relaxation term for the volume fraction equation --------------------------
    CALL press_relax_term( pressure_relaxation )
    relaxation_term(idx_vol1_eqn) = pressure_relaxation * radius**2


    ! relaxation term for the mixture momentum ----------------------------------
    relaxation_term(idx_mix_mom_eqn) = dcmplx(0.d0,0.d0)


    ! relaxation term for relative velocity -------------------------------------
    CALL vel_relax_term( velocity_relaxation )
    relaxation_term(idx_rel_vel_eqn) = velocity_relaxation * radius**2

    ! relaxation term for the mixture energy ------------------------------------
    relaxation_term(idx_mix_engy_eqn) = dcmplx(0.d0,0.d0)

    ! relaxation term for dissolved gas -----------------------------------------
    CALL f_xdis_eq
    relaxation_term(idx_dis_gas_eqn_first:idx_dis_gas_eqn_last) =               &
         - ( x_d_md(1:n_gas) - x_d_md_eq(1:n_gas) ) * ( 1.d0 - alfa_2 )         &
         * ( rho_1 - sum( beta(1:n_cry) * rho_c(1:n_cry) ) ) / tau_d(1:n_gas)   &
         * radius ** 2
    
    DO i=1,n_gas

       IF ( REAL( x_d_md(i) ) .LT. REAL( x_d_md_eq(i) ) ) then

          relaxation_term(idx_dis_gas_eqn_first+i-1) = dcmplx(0.d0,0.d0)

       END IF

    END DO

    ! relaxation term for exsolved gas ------------------------------------------
    CALL f_xdis_eq

    relaxation_term(idx_ex_gas_eqn_first:idx_ex_gas_eqn_last) =                 &
         ( x_d_md(1:n_gas) - x_d_md_eq(1:n_gas) ) * ( 1.d0 - alfa_2)            &
         * ( rho_1 - sum( beta(1:n_cry) * rho_c(1:n_cry) ) ) / tau_d(1:n_gas)   &
         * radius **2
    
    DO i=1,n_gas

       IF ( REAL( x_d_md(i) ) .LT. REAL( x_d_md_eq(i) ) ) then

          relaxation_term(idx_ex_gas_eqn_first+i-1) = dcmplx(0.d0,0.d0)

       END IF

    END DO

    ! relaxation term for crystallization ---------------------------------------

    IF ( n_mom .LE. 1 ) THEN

       CALL f_beta_eq

       relaxation_term(idx_cry_eqn_first:idx_cry_eqn_last) =                    &
            - ( 1.d0 - alfa_2 ) * rho_c(1:n_cry)                                &
            * ( beta(1:n_cry) - beta_eq(1:n_cry) ) / tau_c(1:n_cry) * radius ** 2

    ELSE

       CALL f_growth_rate
       CALL f_nucleation_rate
       
       DO i = 1,n_cry

          DO j = 0,n_mom-1

             !relaxation_term(idx_cry_eqn_first+n_mom*(i-1)+j) =                 &
             !     function of growth_rate(i) and nucleation_rate(i)
                  
          END DO

       END DO

       
       
    END IF
       
  END SUBROUTINE eval_relaxation_terms

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

  SUBROUTINE eval_forces_terms( force_term )

    USE geometry, ONLY : radius
    USE complexify 
    IMPLICIT none

    COMPLEX*16, INTENT(OUT) :: force_term(n_eqns)

    COMPLEX*16 :: visc_force_1 , visc_force_2
    COMPLEX*16 :: visc_force_1_rel , visc_force_2_rel

    REAL*8 :: gas_wall_drag

    INTEGER :: i

    INTEGER :: idx

    idx = 0

    force_term(1:n_eqns) = dcmplx(0.d0,0.d0)

    ! Term for the mixture density ----------------------------------------------
    force_term(idx_mix_mass_eqn) = dcmplx(0.d0,0.d0)

    ! Term for the first phase volume fraction ----------------------------------
    force_term(idx_vol1_eqn) = dcmplx(0.d0,0.d0)

    ! Term for the mixture momentum equation ------------------------------------
    force_term(idx_mix_mom_eqn) = - rho_mix * grav * radius**2.d0

    CALL mixture_viscosity

    visc_force_1 = - 8.d0 * visc_mix * u_1

    ! Turbulent gas-wall friction (Degruyter et al. 2012)

    gas_wall_drag = 0.03d0

    visc_force_2 = - gas_wall_drag / 4.d0 * radius * rho_2 * cdabs( u_2 ) * u_2 

    visc_force_1 = visc_force_1 * ( 1.d0 - frag_eff )
    visc_force_2 = visc_force_2 * frag_eff

    !visc_force_1 =  visc_force_1 * t
    !visc_force_2 =  visc_force_2 * ( 1.0 - t )


    force_term(idx_mix_mom_eqn) = force_term(idx_mix_mom_eqn) + visc_force_2    &
         + visc_force_1 

    ! Term for the relative velocity equation -----------------------------------
    force_term(idx_rel_vel_eqn) = dcmplx(0.0,0.0)

    IF ( drag_funct_model .EQ. 'single_velocity' ) THEN

       force_term(idx_rel_vel_eqn) = dcmplx( 0.d0 , 0.d0 ) 

    ELSE

       visc_force_1_rel = visc_force_1 / ( ( 1.d0 - alfa_2 ) * rho_1 ) 

       visc_force_2_rel = - visc_force_2 / ( alfa_2 * rho_2 )

       force_term(idx_rel_vel_eqn) = force_term(idx_rel_vel_eqn)                &
            + visc_force_2_rel + visc_force_1_rel

    END IF

    ! Term for the mixture energy -----------------------------------------------
    IF ( isothermal ) THEN

       force_term(idx_mix_engy_eqn) = dcmplx(0.d0,0.d0)

    ELSE

       force_term(idx_mix_engy_eqn) = - rho_mix * u_mix * grav * radius**2

       force_term(idx_mix_engy_eqn) = force_term(idx_mix_engy_eqn)              &
            + 2.d0 * visc_force_2 * u_2 + 2.d0 * visc_force_1 * u_1

    END IF

    ! Terms for the exsolved gas phases -----------------------------------------
    DO i=idx_dis_gas_eqn_first,idx_dis_gas_eqn_last

       force_term(i) = dcmplx(0.d0,0.d0)

    END DO
    
    ! Terms for the dissolved gas phases ----------------------------------------
    DO i = idx_ex_gas_eqn_first,idx_ex_gas_eqn_last

       force_term(i) = dcmplx(0.d0,0.d0)

    END DO

    ! Terms for the crystal phases ----------------------------------------------
    DO i= idx_cry_eqn_first,idx_cry_eqn_first

       force_term(i) = dcmplx(0.d0,0.d0)

    END DO


  END SUBROUTINE eval_forces_terms


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

  SUBROUTINE eval_source_terms( source_term )

    USE geometry, ONLY : radius
    USE complexify 
    IMPLICIT none

    COMPLEX*16, INTENT(OUT) :: source_term(n_eqns)

    COMPLEX*16 :: q_lat

    COMPLEX*16 :: water_mass_flux

    COMPLEX*16 :: heat_flux

    INTEGER :: i

    INTEGER :: idx

    idx = 0

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

    IF ( lateral_degassing ) THEN
       
       CALL lithostatic_pressure
       
       IF ( p_2 .GE. p_lith ) THEN
          
          q_lat = ( rho_2 * alfa_2 * k_cr * ( p_2 - p_lith ) ) /                &
               ( visc_2 * radius )
          
       ELSE
          
          q_lat = dcmplx(0.d0,0.d0)
          
       END IF
       
    END IF


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

    source_term(1:n_eqns) = dcmplx(0.d0,0.d0)

    ! --- TOTAL MASS SOURCE TERM ------------------------------------------------
    IF ( lateral_degassing ) THEN

       source_term(idx_mix_mass_eqn) = - 2.d0 * q_lat * radius

    ELSE

       source_term(idx_mix_mass_eqn) = dcmplx(0.d0,0.d0)

    END IF

    IF ( ext_water ) THEN

       rho_w = 971.80  ! 80C. For shallow aquifers, the effect of P is expected
                       ! to be less important
       
       IF ( ( zeta_lith .GT. min_z_influx ) .AND.                               &
            ( zeta_lith .LT. min_z_influx + delta_z_influx ) ) THEN
          
          IF ( total_water_influx .GT. 0.d0 ) THEN
             
             water_mass_flux = dcmplx( total_water_influx/delta_z_influx , 0.d0)
             
          ELSE
             
             IF ( aquifer_type .EQ. 'unconfined' ) THEN
                
                CALL hydrostatic_pressure
  
                IF ( p_hydro .GE. REAL(p_1) ) then
                   
                   visc_w = 2.414d-5 * 10.d0 ** ( 247.8d0 / ( t_w - 140.d0 ) )
                   
                   water_mass_flux = ( 2.d0 * radius * 3.14d0) * rho_w * k_cr / &
                        visc_w * ( p_hydro - p_1 ) / radius
                   
                ELSE
                   
                   water_mass_flux = dcmplx( 0.d0 , 0.d0 )
                   
                END IF

             ELSEIF ( aquifer_type .EQ. 'confined') THEN

                CALL lithostatic_pressure
         
                IF ( p_lith .GE. REAL(p_1) ) then
                   
                   visc_w = 2.414d-5 * 10.d0 ** ( 247.8d0 / ( t_w - 140.d0 ) )
                   
                   water_mass_flux = ( 2.d0 * radius * 3.14d0) * rho_w * k_cr / &
                        visc_w * ( p_lith - p_1 ) / radius
                   
                ELSE
                   
                   water_mass_flux = dcmplx( 0.d0 , 0.d0 )
                   
                END IF
                
             END IF
             
          END IF
          
       ELSE
          
          water_mass_flux = dcmplx(0.d0,0.d0)
          
       END IF
       
       source_term(idx_mix_mass_eqn) = source_term(idx_mix_mass_eqn)            &
            + water_mass_flux
       
    END IF
    
    ! --- FIRST PHASE VOLUME FRACTION -------------------------------------------
    source_term(idx_vol1_eqn) = dcmplx(0.d0,0.d0)

    ! --- Mixture Momentum ------------------------------------------------------
    IF ( lateral_degassing ) THEN

       source_term(idx_mix_mom_eqn) = - 2.d0 * q_lat * radius * u_2

    ELSE

       source_term(idx_mix_mom_eqn) = dcmplx(0.d0,0.d0)

    END IF

    ! --- Relative Velocity -----------------------------------------------------
    source_term(idx_rel_vel_eqn) = dcmplx(0.d0,0.d0)

    ! --- MIXTURE ENERGY --------------------------------------------------------
    IF ( isothermal ) THEN

       source_term(idx_mix_engy_eqn) = dcmplx(0.d0,0.d0)

    ELSE

       IF ( lateral_degassing ) THEN

          source_term(idx_mix_engy_eqn) = - 2.d0 * q_lat * radius * alfa_2 *    &
               ( t * cv_2 + 0.5d0 * u_2*u_2 )

       ELSE

          source_term(idx_mix_engy_eqn) = dcmplx(0.d0,0.d0)

       END IF

       IF ( ext_water ) THEN

          IF ( ( zeta_lith .GT. min_z_influx ) .AND.                            &
               ( zeta_lith .LT. min_z_influx + delta_z_influx ) ) THEN

             IF ( inst_vaporization ) THEN

                heat_flux = - water_mass_flux * ( cv_d(1) * ( t_boiling-t_w )   &
                     + cv_2 * ( - t_boiling ) + lambda_w )

             ELSE

                heat_flux = water_mass_flux * cv_d(1) * t_w

             END IF

             source_term(idx_mix_engy_eqn) = source_term(idx_mix_engy_eqn)       &
                  + heat_flux

          END IF

       END IF

    END IF

    ! --- DISSOLVED GAS BULK DENSITY --------------------------------------------

    ! H2O source due to inlet
    IF ( ext_water .AND. .NOT.inst_vaporization ) THEN

       source_term(idx_dis_gas_eqn_first) = water_mass_flux

    ELSE

       source_term(idx_dis_gas_eqn_first) = dcmplx(0.d0,0.d0)

    END IF

    ! Other gas phases (from 2 to n_gas)
    DO i = idx_dis_gas_eqn_first+1,idx_dis_gas_eqn_last

       source_term(i) = dcmplx(0.d0,0.d0)

    END DO

    ! --- EXSOLVED GAS MASS FRACTIONS -------------------------------------------

    DO i = idx_ex_gas_eqn_first,idx_ex_gas_eqn_last

       IF ( lateral_degassing ) THEN

          source_term(i) = - 2.d0 * q_lat * radius

       ELSE

          source_term(i) = dcmplx(0.d0,0.d0)

       END IF

    END DO

    ! --- CRYSTALS BULK DENSITY -------------------------------------------------

    DO i = idx_cry_eqn_first,idx_cry_eqn_last

       source_term(i) = dcmplx(0.d0,0.d0)

    END DO

  END SUBROUTINE eval_source_terms


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

  SUBROUTINE eval_explicit_forces( expl_forces_term )

    USE geometry, ONLY : radius

    IMPLICIT NONE

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

    INTEGER :: i 

    INTEGER :: idx

    idx = 0

    expl_forces_term(1:n_eqns) = 0.d0


    ! --- TOTAL MASS TERM -------------------------------------------------------
    expl_forces_term(idx_mix_mass_eqn) = 0.d0

    ! --- FIRST PHASE VOLUME FRACTION -------------------------------------------
    expl_forces_term(idx_vol1_eqn) = 0.d0

    ! --- Mixture Momentum ------------------------------------------------------
    expl_forces_term(idx_mix_mom_eqn) = dreal( rho_mix * grav * radius ** 2 )

    ! --- Relative Velocity -----------------------------------------------------
    expl_forces_term(idx_rel_vel_eqn) = 0.d0

    ! --- MIXTURE ENERGY --------------------------------------------------------
    IF ( isothermal ) THEN

       expl_forces_term(idx_mix_engy_eqn) = 0.d0

    ELSE

       expl_forces_term(idx_mix_engy_eqn) = dreal( rho_mix * u_mix * grav *     &
            radius **2 )

    END IF

    ! --- DISSOLVED GAS BULK DENSITY --------------------------------------------
    DO i = idx_dis_gas_eqn_first,idx_dis_gas_eqn_last

       expl_forces_term(i) = 0.d0

    END DO

    ! --- EXSOLVED GAS MASS FRACTIONS -------------------------------------------
    DO i = idx_ex_gas_eqn_first,idx_ex_gas_eqn_last

       expl_forces_term(i) = 0.d0

    END DO

    ! --- CRYSTALS BULK DENSITY ------------------------------------------------
    DO i = idx_cry_eqn_first,idx_cry_eqn_last

       expl_forces_term(i) = 0.d0

    END DO

  END SUBROUTINE eval_explicit_forces


END MODULE equations