constitutive.f90

Go to the documentation of this file.

!********************************************************************************
!********************************************************************************
MODULE 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
  
  IMPLICIT none

  !--------- Constants for the equations of state -------------------------------
  REAL*8 :: cv_2
  REAL*8 :: cv_m
  REAL*8, ALLOCATABLE :: cv_c(:)
  REAL*8, ALLOCATABLE :: cv_d(:)
  REAL*8, ALLOCATABLE :: cv_g(:)

  !~   REAL*8 :: C0_2        !< exsolved gas sound speed at atmospheric conditions
  REAL*8 :: c0_m
  REAL*8, ALLOCATABLE :: c0_c(:)
  REAL*8, ALLOCATABLE :: c0_d(:)
  REAL*8, ALLOCATABLE :: c0_g(:)

  REAL*8 :: gamma_2
  REAL*8 :: gamma_m
  REAL*8, ALLOCATABLE :: gamma_c(:)
  REAL*8, ALLOCATABLE :: gamma_d(:)
  REAL*8, ALLOCATABLE :: gamma_g(:)

  REAL*8 :: rho0_2
  REAL*8 :: rho0_m
  REAL*8, ALLOCATABLE :: rho0_c(:)
  REAL*8, ALLOCATABLE :: rho0_d(:)
  REAL*8, ALLOCATABLE :: rho0_g(:)

  REAL*8 :: t0_2
  REAL*8 :: t0_m
  REAL*8, ALLOCATABLE :: t0_c(:)
  REAL*8, ALLOCATABLE :: t0_d(:)
  REAL*8, ALLOCATABLE :: t0_g(:)

  !~   REAL*8 :: p0_2       !< exsolved gas reference pressure
  REAL*8 :: p0_m
  REAL*8, ALLOCATABLE :: p0_c(:)
  REAL*8, ALLOCATABLE :: p0_d(:)
  REAL*8, ALLOCATABLE :: p0_g(:)

  REAL*8 :: bar_e_2
  REAL*8 :: bar_e_m
  REAL*8, ALLOCATABLE :: bar_e_c(:)
  REAL*8, ALLOCATABLE :: bar_e_d(:)
  REAL*8, ALLOCATABLE :: bar_e_g(:)

  REAL*8 :: bar_p_m
  REAL*8, ALLOCATABLE :: bar_p_c(:)
  REAL*8, ALLOCATABLE :: bar_p_d(:)
  REAL*8, ALLOCATABLE :: bar_p_g(:)

  !~   REAL*8 :: s0_2      !< exsolved gas reference entropy
  REAL*8 :: s0_m
  REAL*8, ALLOCATABLE :: s0_c(:)
  REAL*8, ALLOCATABLE :: s0_d(:)
  REAL*8, ALLOCATABLE :: s0_g(:)

  !-----------------------------------------------------------------------------------------!
  REAL*8, ALLOCATABLE :: pc_g(:)
  REAL*8, ALLOCATABLE :: tc_g(:)
  REAL*8, ALLOCATABLE :: a_g(:)
  REAL*8, ALLOCATABLE :: b_g(:)

  CHARACTER*20 :: gas_law
  !----------------------------------------------------------------------------------------!


  COMPLEX*16 :: cv_1

  COMPLEX*16 :: e_mix

  COMPLEX*16 :: e_1
  COMPLEX*16 :: e_2
  COMPLEX*16 :: e_m
  COMPLEX*16, ALLOCATABLE :: e_c(:)
  COMPLEX*16, ALLOCATABLE :: e_d(:)
  COMPLEX*16, ALLOCATABLE :: e_g(:)

  COMPLEX*16 :: p_1
  COMPLEX*16 :: p_2

  COMPLEX*16 :: s_1
  COMPLEX*16 :: s_2
  COMPLEX*16 :: s_m
  COMPLEX*16, ALLOCATABLE :: s_c(:)
  COMPLEX*16, ALLOCATABLE :: s_d(:)
  COMPLEX*16, ALLOCATABLE :: s_g(:)

  COMPLEX*16 :: mu_1
  COMPLEX*16 :: mu_2
  COMPLEX*16 :: mu_m
  COMPLEX*16, ALLOCATABLE :: mu_c(:)
  COMPLEX*16, ALLOCATABLE :: mu_d(:)
  COMPLEX*16, ALLOCATABLE :: mu_g(:)

  COMPLEX*16 :: rho_1
  COMPLEX*16 :: rho_2
  COMPLEX*16 :: rho_m
  COMPLEX*16, ALLOCATABLE :: rho_c(:)
  COMPLEX*16, ALLOCATABLE :: rho_d(:)
  COMPLEX*16, ALLOCATABLE :: rho_g(:)
  COMPLEX*16 :: rho_md

  COMPLEX*16, ALLOCATABLE :: rhob_c(:)
  COMPLEX*16 :: rhob_m

  COMPLEX*16 :: alfa_1
  COMPLEX*16 :: alfa_2
  COMPLEX*16, ALLOCATABLE :: alfa_g(:)

  COMPLEX*16 :: alfarho_2

  COMPLEX*16 :: alfa_m_1
  COMPLEX*16, ALLOCATABLE :: alfa_d_1(:)
  COMPLEX*16, ALLOCATABLE :: alfa_g_2(:)


  COMPLEX*16, ALLOCATABLE :: beta(:)
  COMPLEX*16, ALLOCATABLE :: beta_eq(:)

  COMPLEX*16, ALLOCATABLE :: mom_cry(:,:)

  COMPLEX*16, ALLOCATABLE :: growth_rate(:)
  COMPLEX*16, ALLOCATABLE :: nucleation_rate(:)
  
  COMPLEX*16 :: u_1
  COMPLEX*16 :: u_2

  COMPLEX*16 :: x_1
  COMPLEX*16 :: x_2
  COMPLEX*16, ALLOCATABLE :: x_c(:)
  COMPLEX*16 :: x_m
  COMPLEX*16, ALLOCATABLE :: x_d(:)
  COMPlEX*16, ALLOCATABLE :: x_g(:)

  COMPLEX*16, ALLOCATABLE :: x_d_md(:)
  COMPLEX*16, ALLOCATABLE :: x_d_md_eq(:)


  COMPLEX*16, ALLOCATABLE :: x_c_1(:)
  COMPLEX*16 :: x_m_1
  COMPLEX*16 :: x_md_1
  COMPLEX*16, ALLOCATABLE :: x_d_1(:)
  COMPLEX*16, ALLOCATABLE :: x_g_2(:)

  COMPLEX*16 :: t
  COMPLEX*16 :: rho_mix
  COMPLEX*16 :: u_mix
  COMPLEX*16 :: cv_mix
  COMPLEX*16 :: visc_mix

  COMPLEX*16 :: s

  COMPLEX*16 :: c_1
  COMPLEX*16 :: c_2

  INTEGER :: n_drag_models

  CHARACTER (LEN=30), DIMENSION(20) :: available_drag_models
  
  CHARACTER*30 :: drag_funct_model

  COMPLEX*16 :: drag_funct

  REAL*8 :: drag_funct_coeff

  ! parameters of the Forchheimer model

  REAL*8 :: bubble_number_density
  REAL*8 :: log10_bubble_number_density

  REAL*8 :: tortuosity_factor

  REAL*8 :: throat_bubble_ratio

  ! friction coefficient
  REAL*8 :: friction_coefficient

  ! drag coefficient
  REAL*8 :: c_d 

  ! ash particl size
  REAL*8 :: r_a

  CHARACTER*20 :: p_relax_model

  COMPLEX*16 :: tau_p

  REAL*8 :: tau_p_coeff

  REAL*8, ALLOCATABLE :: tau_c(:)

  REAL*8, ALLOCATABLE :: beta0(:)

  REAL*8, ALLOCATABLE :: beta_max(:)

  CHARACTER*20 :: crystallization_model

  REAL*8 :: frag_eff

  INTEGER :: fragmentation_model

  REAL*8 :: frag_thr

  REAL*8 :: tau_frag_exp

  COMPLEX*16 :: tau_frag

  REAL*8 :: tau_frag_coeff

  REAL*8, ALLOCATABLE :: tau_d(:)

  CHARACTER*20 :: exsol_model

  REAL*8, ALLOCATABLE :: solub(:)

  REAL*8, ALLOCATABLE :: solub_exp(:)

  REAL*8, ALLOCATABLE :: x_ex_dis_in(:)

  REAL*8 :: grav

  REAL*8 :: k_cr

  REAL*8 :: rho_cr

  REAL*8 :: visc_2

  COMPLEX*16 :: visc_melt

  COMPLEX*16 :: visc_1

  COMPLEX*16 :: visc_rel_bubbles

  INTEGER :: n_theta_models 

  CHARACTER (LEN=30), DIMENSION(20) :: available_theta_models
  
  CHARACTER*30 :: theta_model

  COMPLEX*16 :: theta

  REAL*8 :: theta_fixed

  REAL*8 :: c1 , c2 , c3

  REAL*8 :: p_lith

  REAL*8 :: p_hydro

  REAL*8 :: zeta_lith

  INTEGER :: n_bubble_models

  CHARACTER (LEN=30), DIMENSION(20) :: available_bubble_models
  
  CHARACTER*20 :: bubbles_model

  LOGICAL :: explosive

  COMPLEX*16 :: permkc 
  COMPLEX*16 :: inv_permkc  

  REAL*8 :: alfa_switch
  REAL*8 :: a_2nd , b_2nd , perm0

  INTEGER :: n_visc_melt_models

  CHARACTER (LEN=30), DIMENSION(20) :: available_visc_melt_models
  
  CHARACTER*30 :: visc_melt_model

  
  REAL*8, DIMENSION(12) :: wt_init

  REAL*8 :: rho_w

  REAL*8 :: xa,xb,xc


CONTAINS

  SUBROUTINE initialize_models

    IMPLICIT NONE

    n_drag_models = 6
    
    available_drag_models(1) = 'constant'
    available_drag_models(2) = 'eval'
    available_drag_models(3) = 'darcy'
    available_drag_models(4) = 'forchheimer'
    available_drag_models(5) = 'drag'
    available_drag_models(6) = 'single_velocity'


    n_theta_models = 11

    available_theta_models(1) = 'Lejeune_and_Richet1995'
    available_theta_models(2) = 'Dingwell1993'
    available_theta_models(3) = 'Melnik_and_Sparks1999'
    available_theta_models(4) = 'Costa2005'
    available_theta_models(5) = 'Melnik_and_Sparks2005'
    available_theta_models(6) = 'Vona_et_al2011'
    available_theta_models(7) = 'Vona_et_al2011_mod'
    available_theta_models(8) = 'Vona_et_al2013_eq19'
    available_theta_models(9) = 'Vona_et_al2013_eq20'
    available_theta_models(10) = 'Vona_et_al2013_eq21'
    available_theta_models(11) = 'Fixed_value'

    n_bubble_models = 11

    available_bubble_models(1) = 'none'
    available_bubble_models(2) = 'Costa2007'
    available_bubble_models(3) = 'Einstein'
    available_bubble_models(4) = 'Quane-Russel'
    available_bubble_models(5) = 'Eilers'
    available_bubble_models(6) = 'Sibree'
    available_bubble_models(7) = 'Taylor'
    available_bubble_models(8) = 'Mackenzie'
    available_bubble_models(9) = 'DucampRaj'
    available_bubble_models(10) = 'BagdassarovDingwell'
    available_bubble_models(11) = 'Rahaman'


    n_visc_melt_models = 6

    available_visc_melt_models(1) = 'Hess_and_Dingwell1996'
    available_visc_melt_models(2) = 'Romano_et_al2003'
    available_visc_melt_models(3) = 'Giordano_et_al2008'
    available_visc_melt_models(4) = 'Giordano_et_al2009'
    available_visc_melt_models(5) = 'Di_Genova_et_al2013_eqn_3,5'
    available_visc_melt_models(6) = 'Di_Genova_et_al2013_eqn_4,5'
    
  END SUBROUTINE initialize_models
  
  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE allocate_phases_parameters

    ALLOCATE( cv_c(n_cry) )
    ALLOCATE( c0_c(n_cry) )
    ALLOCATE( gamma_c(n_cry) )
    ALLOCATE( rho0_c(n_cry) )
    ALLOCATE( t0_c(n_cry) )
    ALLOCATE( p0_c(n_cry) )
    ALLOCATE( bar_e_c(n_cry) )
    ALLOCATE( bar_p_c(n_cry) )
    ALLOCATE( s0_c(n_cry) )
    ALLOCATE( e_c(n_cry) )
    ALLOCATE( s_c(n_cry) )
    ALLOCATE( mu_c(n_cry) )
    ALLOCATE( rho_c(n_cry) )
    ALLOCATE( rhob_c(n_cry) )
    ALLOCATE( beta(n_cry) )
    ALLOCATE( beta_eq(n_cry) )
    ALLOCATE( x_c(n_cry) )
    ALLOCATE( x_c_1(n_cry) )
    ALLOCATE( tau_c(n_cry) )
    ALLOCATE( beta0(n_cry) )
    ALLOCATE( beta_max(n_cry) )

    IF ( n_mom .GT. 1 ) THEN

       ALLOCATE( mom_cry(1:n_cry,0:n_mom-1) )
       ALLOCATE( growth_rate(n_cry) )
       ALLOCATE( nucleation_rate(n_cry) )

    END IF

    ALLOCATE( cv_d(n_gas) )
    ALLOCATE( c0_d(n_gas) )
    ALLOCATE( gamma_d(n_gas) )
    ALLOCATE( rho0_d(n_gas) )
    ALLOCATE( t0_d(n_gas) )
    ALLOCATE( p0_d(n_gas) )
    ALLOCATE( bar_e_d(n_gas) )
    ALLOCATE( bar_p_d(n_gas) )
    ALLOCATE( s0_d(n_gas) )
    ALLOCATE( e_d(n_gas) )
    ALLOCATE( s_d(n_gas) )
    ALLOCATE( mu_d(n_gas) )
    ALLOCATE( rho_d(n_gas) )
    ALLOCATE( alfa_d_1(n_gas) )

    ALLOCATE( x_d(n_gas) )
    ALLOCATE( x_d_md(n_gas) )
    ALLOCATE( x_d_md_eq(n_gas) )
    ALLOCATE( x_d_1(n_gas) )
    ALLOCATE( tau_d(n_gas) )
    ALLOCATE( solub(n_gas) )
    ALLOCATE( solub_exp(n_gas) )

    ALLOCATE( x_ex_dis_in(n_gas) )

    ALLOCATE( cv_g(n_gas) )
    ALLOCATE( c0_g(n_gas) )
    ALLOCATE( gamma_g(n_gas) )
    ALLOCATE( rho0_g(n_gas) )
    ALLOCATE( t0_g(n_gas) )
    ALLOCATE( p0_g(n_gas) )
    ALLOCATE( bar_e_g(n_gas) )
    ALLOCATE( bar_p_g(n_gas) )
    ALLOCATE( s0_g(n_gas) )
    ALLOCATE( e_g(n_gas) )
    ALLOCATE( s_g(n_gas) )
    ALLOCATE( mu_g(n_gas) )
    ALLOCATE( rho_g(n_gas) )
    ALLOCATE( x_g(n_gas) )
    ALLOCATE( x_g_2(n_gas) )
    ALLOCATE( alfa_g(n_gas) )
    ALLOCATE( alfa_g_2(n_gas) )

    ALLOCATE( pc_g(n_gas) )
    ALLOCATE( tc_g(n_gas) )
    ALLOCATE( a_g(n_gas) )
    ALLOCATE( b_g(n_gas) )

  END SUBROUTINE allocate_phases_parameters


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

  SUBROUTINE eos

    USE complexify 
    IMPLICIT none

    e_c(1:n_cry) = dcmplx(bar_e_c(1:n_cry),0.d0) + dcmplx(cv_c(1:n_cry),0.d0)   &
         * t + dcmplx(bar_p_c(1:n_cry),0.d0) / rho_c(1:n_cry)

    e_m = dcmplx(bar_e_m,0.d0) + dcmplx(cv_m,0.d0) * t +                        &
         dcmplx(bar_p_m,0.d0) / rho_m

    e_d(1:n_gas) = dcmplx(bar_e_d(1:n_gas),0.d0) + dcmplx(cv_d(1:n_gas),0.d0)   &
         * t + dcmplx(bar_p_d(1:n_gas),0.d0) / rho_d(1:n_gas)

    e_1 = sum( x_c_1(1:n_cry) * e_c(1:n_cry) ) + x_m_1 * e_m +                  &
         sum( x_d_1(1:n_gas) * e_d(1:n_gas) )

    e_g(1:n_gas) = dcmplx(bar_e_g(1:n_gas),0.d0) + dcmplx(cv_g(1:n_gas),0.d0)   &
         * t - dcmplx(a_g(1:n_gas),0.0) * rho_g(1:n_gas)

    e_2 = sum( x_g_2(1:n_gas) * e_g(1:n_gas) )

    s_c(1:n_cry) = s0_c(1:n_cry) + cv_c(1:n_cry) * log( t / t0_c(1:n_cry) *     &
         ( rho0_c(1:n_cry) / rho_c(1:n_cry) )**( gamma_c(1:n_cry) - 1.d0 ) )

    s_m = s0_m + cv_m * log( t / t0_m * ( rho0_m / rho_m )**( gamma_m - 1.d0 ) )

    s_d(1:n_gas) = s0_d(1:n_gas) + cv_d(1:n_gas) * log( t / t0_d(1:n_gas) *     &
         ( rho0_d(1:n_gas) / rho_d(1:n_gas) )**( gamma_d(1:n_gas) - 1.d0 ) )

    s_1 = sum( x_c_1(1:n_cry) * s_c(1:n_cry) ) + x_m_1 * s_m +                  &
         sum( x_d_1(1:n_gas) * s_d(1:n_gas) )

    s_g(1:n_gas) = s0_g(1:n_gas) + cv_g(1:n_gas) * log( t / t0_g(1:n_gas) *     &
         ( ( rho0_g(1:n_gas) / rho_g(1:n_gas) ) * ( dcmplx(1.d0,0.d0) -         &
         b_g(1:n_gas) * rho_g(1:n_gas) ) ) ** ( gamma_g(1:n_gas) - 1.d0 ) )

    s_2 = sum( x_g_2(1:n_gas) * s_g(1:n_gas) )

    mu_m = e_m + p_1/rho_m - t * s_m

    mu_c(1:n_cry) = e_c(1:n_cry) + p_1/rho_c(1:n_cry) - t * s_c(1:n_cry)

    mu_d(1:n_gas) = e_d(1:n_gas) + p_1/rho_d(1:n_gas) - t * s_d(1:n_gas)

    mu_1 = sum( x_c_1(1:n_cry) * mu_c(1:n_cry) ) + x_m_1 * mu_m +               &
         sum( x_d_1(1:n_gas) * mu_d(1:n_gas) )

    mu_g(1:n_gas) = e_g(1:n_gas) + p_2/rho_g(1:n_gas) - t * s_g(1:n_gas)

    mu_2 = sum( x_g_2(1:n_gas) * mu_g(1:n_gas) )

    !mu_1 = e_1 + p_1/rho_1 - T * s_1
    !mu_2 = e_2 + p_2/rho_2 - T * s_2

    s = x_1 * s_1 + x_2 * s_2

  END SUBROUTINE eos

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

  SUBROUTINE sound_speeds(C_mix,mach) 

    USE complexify 
    IMPLICIT none

    REAL*8, INTENT(OUT) :: C_mix , mach

    COMPLEX*16 :: C2_c(1:n_cry) , C2_m , C2_d(1:n_gas) , C2_1 , C2_2
    COMPLEX*16 :: C2_g(1:n_gas)

    COMPLEX*16 :: K_1 , K_2 , K_mix
    COMPLEX*16 :: K_c(1:n_cry) , K_m , K_d(1:n_gas) , K_g(1:n_gas)

    !REAL*8 :: C_mix_Wallis

    c2_c(1:n_cry) = c0_c(1:n_cry)**2.d0 * ( rho_c(1:n_cry) / rho0_c(1:n_cry) )  &
         ** ( gamma_c(1:n_cry) - dcmplx(1.0,0.d0) ) * cdexp( ( s_c(1:n_cry) -   &
         s0_c(1:n_cry) ) / cv_c(1:n_cry) )

    c2_m = c0_m**2.d0 * ( rho_m / rho0_m ) **                                   &
         ( gamma_m - dcmplx(1.0,0.d0) ) * cdexp( ( s_m - s0_m ) / cv_m )

    c2_d(1:n_gas) = c0_d(1:n_gas)**2.d0 * ( rho_d(1:n_gas) / rho0_d(1:n_gas) )  &
         ** ( gamma_d(1:n_gas) - dcmplx(1.0,0.d0) ) * cdexp( ( s_d(1:n_gas) -   &
         s0_d(1:n_gas) ) / cv_d(1:n_gas) )

    k_c(1:n_cry) = 1.d0 / ( rho_c(1:n_cry) * c2_c(1:n_cry) )
    k_m = 1.d0 / ( rho_m * c2_m )
    k_d(1:n_gas) = 1.d0 / ( rho_d(1:n_gas) * c2_d(1:n_gas) )

    k_1 = alfa_m_1*k_m + sum( beta(1:n_cry)*k_c(1:n_cry) )                      &
         + sum( alfa_d_1(1:n_gas) * k_d(1:n_gas) )

    c2_1 = 1.d0 / ( k_1 * rho_1 )

    c_1 = cdsqrt( c2_1 )

    c2_g(1:n_gas) = cv_g(1:n_gas) * ( gamma_g(1:n_gas) - dcmplx(1.0,0.d0) )     &
         * gamma_g(1:n_gas) * t / ( dcmplx(1.0,0.d0) - b_g(1:n_gas) *           &
         rho_g(1:n_gas) )**2.0 - dcmplx(2.0,0.d0) * a_g(1:n_gas) * rho_g(1:n_gas)

    k_g(1:n_gas) = 1.d0 / ( rho_g(1:n_gas) * c2_g(1:n_gas) )

    k_2 = sum( alfa_g_2(1:n_gas) * k_g(1:n_gas) ) 

    c2_2 = 1.d0 / ( k_2 * rho_2 )

    c_2 = cdsqrt( c2_2 )

    k_mix = alfa_1*k_1 + alfa_2*k_2

    c_mix = dreal( 1.d0 / cdsqrt( k_mix * rho_mix ) )

    ! Sound speed of the mixture (Wallis, 1969)
    ! C_mix = REAL( C_2 / SQRT(x_2) * ( x_2 + (1.d0 - x_2 ) * rho_2 / rho_1 ) )

    mach = dreal( u_mix / c_mix )

    IF ( verbose_level .GE. 2 ) THEN

       WRITE(*,*) 'K_1,k_2',k_1,k_2
       WRITE(*,*) 'K_mix,rho_mix',k_mix,rho_mix
       WRITE(*,*) 'u_mix,C_mix',u_mix,c_mix

    END IF


  END SUBROUTINE sound_speeds

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

  SUBROUTINE eval_densities

    IMPLICIT none

    COMPLEX*16 :: a,b,c
    COMPLEX*16 :: y1,y2,y3

    COMPLEX*16 :: coeff1 , coeff2 , coeff3 

    INTEGER :: i

    rho_c(1:n_cry) = ( p_1 + bar_p_c(1:n_cry) ) / ( t * cv_c(1:n_cry) *         &
         ( gamma_c(1:n_cry) - dcmplx(1.d0,0.d0) ) )

    rho_m = ( p_1 + bar_p_m ) / ( t * cv_m * ( gamma_m - dcmplx(1.d0,0.d0) ) )
    rho_d(1:n_gas) = ( p_1 + bar_p_d(1:n_gas) ) / ( t * cv_d(1:n_gas) *         &
         ( gamma_d(1:n_gas) - dcmplx(1.d0,0.d0) ) )

    rho_md = dcmplx(1.d0,0.d0) / ( sum( x_d_md(1:n_gas) / rho_d(1:n_gas) )      &
         + ( dcmplx(1.d0,0.d0) - sum( x_d_md(1:n_gas) ) ) / rho_m )

    rho_1 = sum( beta(1:n_cry) * rho_c(1:n_cry) ) + ( dcmplx(1.d0,0.d0) -       &
         sum( beta(1:n_cry) ) ) * rho_md

    IF ( gas_law .EQ. 'IDEAL' )  THEN

       rho_g(1:n_gas) = ( p_2 ) / ( t * cv_g(1:n_gas) *                         &
            ( gamma_g(1:n_gas) - dcmplx(1.d0,0.d0) ) )

    ELSEIF (  gas_law .EQ. 'VDW' ) THEN

       DO i=1,n_gas

          a = a_g(i) * dcmplx(-1.d0,0.d0)
          b = cv_g(i) * ( gamma_g(i) - dcmplx(1.d0,0.d0) ) * t  + p_2 * b_g(i)
          c = - p_2

          coeff1 = a / ( a_g(i) * b_g(i) )
          coeff2 = b / ( a_g(i) * b_g(i) )
          coeff3 = c / ( a_g(i) * b_g(i) )

          CALL solve_cubic( coeff1 , coeff2 , coeff3 , y1 , y2 , y3 )

          rho_g(i) = y1

       END DO

    END IF

    rho_2 = sum( alfa_g_2(1:n_gas) * rho_g(1:n_gas) ) 

  END SUBROUTINE eval_densities

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

  SUBROUTINE solve_cubic(a,b,c,y1,y2,y3)

    USE complexify
    IMPLICIT NONE

    COMPLEX*16, INTENT(IN) :: a,b,c
    COMPLEX*16, INTENT(OUT) :: y1,y2,y3

    COMPLEX*16 :: Q,R,D_temp
    COMPLEX*16 :: phi , arg , u , v

    q = ( 3.d0*b - a**2) / dcmplx(9.d0,0.d0)
    r = ( 9.d0*a*b - 27.d0*c - 2.d0*a**3 ) / dcmplx(54.d0,0.d0)

    d_temp = q**3 + r**2

    IF ( d_temp .GE. 0 ) THEN

       arg = r + cdsqrt(d_temp)

       u = sign(1.d0,REAL(arg))*(abs(arg))**(1.D0/3.D0)

       arg = r - cdsqrt(d_temp)
       v = sign(1.d0,REAL(arg))*(abs(arg))**(1.D0/3.D0)

       y1 = u + v - a / dcmplx(3.d0,0.d0)

    ELSE

       phi = acos(r/cdsqrt(-(q**3)))

       y1 = 2.d0 * cdsqrt(-q) * cos(phi/ dcmplx(3.d0,0.d0) )                    &
            - a / dcmplx(3.d0,0.d0)

       y2 = 2.d0 * cdsqrt(-q) * cos((phi+2.d0*pi)/ dcmplx(3.d0,0.d0) )          &
            - a / dcmplx(3.d0,0.d0)

       y3 = 2.d0 * cdsqrt(-q) * cos((phi-2.d0*pi)/ dcmplx(3.d0,0.d0) )          &
            - a/ dcmplx(3.d0,0.d0)

    END IF

  END SUBROUTINE solve_cubic


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

  SUBROUTINE f_xdis_eq

    USE complexify 

    IMPLICIT NONE


    COMPLEX*16 :: aa , bb , cc , pp
    COMPLEX*16 :: P_bar, log10Ds, Ds, Dcl
    REAL*8 :: beta_S(4)
    REAL*8 :: melt_volume

    COMPLEX*16 :: melt_mass
    COMPLEX*16 :: co2_exs_mass,h2o_exs_mass,S_exs_mass,Cl_exs_mass
    COMPLEX*16 :: S_dis_mass,Cl_dis_mass
    COMPLEX*16 :: gas_mass, Ratio_exs_dis
    COMPLEX*16 :: S_mass,Cl_mass


    IF ( REAL(p_2) .LE. 0.D0 ) then

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

    ELSE

       SELECT CASE ( exsol_model )

       CASE DEFAULT


       CASE ( 'Henry' )

          ! Henry's law

          x_d_md_eq(1:n_gas) = solub(1:n_gas) * (alfa_g_2(1:n_gas) * p_2 )      &
               ** solub_exp(1:n_gas)

       CASE ( 'Zhang' )

          ! Zhang

          aa = 0.4874d0 - 608.0d0 / t + 489530.d0 / ( t * t )

          bb = -0.06062d0 + 135.6d0 / t - 69200.d0 / ( t * t )

          cc = 0.00253d0 - 4.154d0 / t + 1509.0d0 / ( t * t )

          pp = p_2 * 1.0d-6

          x_d_md_eq(1:n_gas) = 1.0d-2 * ( aa * cdsqrt(pp) + bb * pp + cc        &
               * pp ** 1.5d0 )

       CASE ( 'H20-CO2-S-CL' )

          ! Henry's law (s is global variable and is a constant for rhyolite)
          
          melt_volume=1.0 !%m3
          
          melt_mass=melt_volume*rho_md
          
          !water exsolved mass
          h2o_exs_mass=melt_mass*x_g(1)
          
          !co2 exsolved mass
          co2_exs_mass=melt_mass*x_g(2)
          
          !S exsolved mass
          s_exs_mass=melt_mass*x_g(3)
          
          !Cl exsolved mass
          cl_exs_mass=melt_mass*x_g(4)
          
          !Total volatile mass
          gas_mass=co2_exs_mass+h2o_exs_mass+s_exs_mass+cl_exs_mass 
          
          !Ratio_exs_dis is ratio of all gas mass to melt mass
          ratio_exs_dis=gas_mass/melt_mass 
          
          !Pressure in bar
          p_bar = p_2 / 1e5;
          
          !fitting coefficient
          beta_s = [4.8095e+02, -4.5850e+02, 3.0240e-03, -2.7519e-07] 
          
          !Solubility for S and Cl            
          log10ds = beta(1) * p_bar**(-0.0075) + beta(2) + beta(3) * p_bar      &
               + beta(4) * p_bar**(2.0)
          ds = 10**(log10ds) 
          dcl = dcmplx(2.d0, 0.d0)
          
          !Equilibrium dissolved water    
          x_d_md_eq(1) = solub(1) * (alfa_g_2(1) * p_2) ** solub_exp(1)
          !Equilibrium dissolved co2      
          x_d_md_eq(2) = solub(2) * (alfa_g_2(2) * p_2) ** solub_exp(2)

          !Total mass of S
          s_mass=melt_mass*x_ex_dis_in(3)
          
          !Dissolved mass of S 
          s_dis_mass=s_mass - (ds*ratio_exs_dis/(ds*ratio_exs_dis + 1.0))*s_mass
          
          !Equilibrium dissolved S
          x_d_md_eq(3)=s_dis_mass/melt_mass 
          
          !Total mass of S
          cl_mass=melt_mass*x_ex_dis_in(4)
          
          !Dissolved mass of Cl 
          cl_dis_mass=cl_mass- (dcl*ratio_exs_dis/(dcl* ratio_exs_dis + 1.0))   &
               *cl_mass
          
          !Equilibrium dissolved Cl
          x_d_md_eq(4)=cl_dis_mass/melt_mass 
          
       END SELECT

    END IF

  END SUBROUTINE f_xdis_eq

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

  SUBROUTINE f_beta_eq

    USE complexify 
    IMPLICIT NONE

    COMPLEX*16 :: x_d_md_tot ,x_d_md_wt_tot 
    COMPLEX*16 :: p_1_bar, T_celsius
    COMPLEX*16 :: crystal_mass_fraction(1:n_cry)

    INTEGER :: j

    p_1_bar = p_1 / 1.0e5
    t_celsius = t - 273.d0

    x_d_md_tot = sum( x_d_md(1:n_gas) )
    x_d_md_wt_tot = x_d_md_tot * 100.d0

    SELECT CASE ( crystallization_model )

    CASE DEFAULT

    CASE ( 'Vitturi2010' )

       DO j=1,n_cry
          
          !----------------------------------------------------------------------
          beta_eq(j)=beta0(j) + 0.55d0*( 0.58815d0*( p_1/1.d6 )**( -0.5226d0 ) )
          !----------------------------------------------------------------------
          
          beta_eq(j) = max( beta0(j) + 1d-15, beta_eq(j) )
          beta_eq(j) = min( beta_max(j) , beta_eq(j) )
          
       END DO

    END SELECT

  END SUBROUTINE f_beta_eq

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

  SUBROUTINE f_growth_rate

    USE complexify 
    IMPLICIT NONE

    integer :: i
    
    DO i = 1,n_cry

       growth_rate(i) = 0.d0
       
    END DO
    
  END SUBROUTINE f_growth_rate


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

  SUBROUTINE f_nucleation_rate

    USE complexify 
    IMPLICIT NONE

    integer :: i
    
    DO i = 1,n_cry

       nucleation_rate(i) = 0.d0
       
    END DO

    
  END SUBROUTINE f_nucleation_rate
    
  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE lithostatic_pressure

    USE geometry, ONLY : zn

    USE complexify 
    IMPLICIT NONE

    p_lith = ( zn - zeta_lith ) * rho_cr * grav

  END SUBROUTINE lithostatic_pressure

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

  SUBROUTINE hydrostatic_pressure

    USE geometry

    USE complexify 
    IMPLICIT NONE

    p_hydro = ( zn - zeta_lith ) * rho_w * grav

  END SUBROUTINE hydrostatic_pressure


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

  SUBROUTINE press_relax_term( pressure_relaxation )

    USE complexify
    IMPLICIT NONE

    COMPLEX*16, INTENT(OUT) :: pressure_relaxation

    CALL mixture_viscosity

    SELECT CASE (  p_relax_model )

    CASE DEFAULT

       tau_p = dcmplx( 1.d0 , 0.d0 )

    CASE ( 'eval' ) 

       IF ( ( explosive ) .AND. ( frag_eff .GT. 0.d0 ) ) THEN

          tau_p = ( visc_mix ** ( 1.d0 - frag_eff ) * visc_2 ** frag_eff )      &
               / ( alfa_1 * alfa_2 * rho_mix )

       ELSE

          tau_p = visc_mix / ( alfa_1 * alfa_2 * rho_mix )

       END IF


    CASE ( 'eval2' ) 


       tau_p = ( 1.d0 - frag_eff ) * visc_mix /                                 &
            ( (alfa_1 ** alfa_1) * rho_mix ) + 0.001 * frag_eff


    CASE ( 'constant' )

       tau_p = dcmplx( 1.d0 , 0.d0 )

    CASE ( 'single_pressure' )

       tau_p = dcmplx( 1.d0 , 0.d0 )

    END SELECT

    tau_p = tau_p_coeff * tau_p

    pressure_relaxation = - ( p_2 - p_1 ) / tau_p

  END SUBROUTINE press_relax_term

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

  SUBROUTINE vel_relax_term( velocity_relaxation )

    USE complexify
    IMPLICIT NONE

    COMPLEX*16, INTENT(OUT) :: velocity_relaxation

    COMPLEX*16 :: Rey 
    COMPLEX*16 :: diam
    COMPLEX*16 :: radius_bubble
    COMPLEX*16 :: effusive_drag , explosive_drag

    COMPLEX*16 :: k_1 , k_2

    COMPLEX*16 :: throat_radius
    
    REAL*8 :: frag_transition, t

    frag_transition = frag_thr + 0.05d0

    t = min(1.0d0,max(0.0d0, ( REAL(alfa_2) - frag_thr ) /                      &
         (frag_transition - frag_thr )))
         
    SELECT CASE ( drag_funct_model )

    CASE DEFAULT

       effusive_drag = dcmplx(1.d0,0.d0)

       explosive_drag = effusive_drag

    CASE ( 'eval' )

       ! permeability model

       CALL f_permkc

       diam = ( alfa_2 / ( 4.0 / 3.0 * pi *  bubble_number_density              &
            * ( alfa_1 ) ) ) ** ( 1.d0 / 3.d0 )

       rey = alfa_2 * rho_2 * diam * cdabs( u_2 - u_1 ) / visc_2

       effusive_drag = visc_2 * inv_permkc * ( 1.d0 + 1.d-2 * rey / alfa_1 )

       explosive_drag = effusive_drag

       IF ( verbose_level .GE. 3 ) THEN

          WRITE(*,*) 'Rey',rey
          WRITE(*,*) 'visc_2',visc_2
          WRITE(*,*) 'rho_2',rho_2
          WRITE(*,*) 'inv_permkc',inv_permkc

       END IF

    CASE ( 'Klug_and_Cashman' )

       ! permeability model

       CALL f_permkc

       diam = ( alfa_2 / ( 4.0 / 3.0 * pi *  bubble_number_density              &
            * ( alfa_1 ) ) ) ** ( 1.d0 / 3.d0 )

       rey = alfa_2 * rho_2 * diam * cdabs( u_2 - u_1 ) / visc_2

       effusive_drag = visc_2 * inv_permkc * ( 1.d0 + 1.d-2 * rey / alfa_1 )

       explosive_drag = 3.d0 * c_d / ( 8.d0 * r_a ) * rho_2 * cdabs( u_2 - u_1 ) 

    CASE ( 'darcy' )

       ! Darcy formulation from Eq. 16 Degruyter et al. 2012 

       radius_bubble = ( alfa_2 / ( 4.0 / 3.0 * pi *  bubble_number_density     &
            * ( alfa_1 ) ) ) ** ( 1.d0 / 3.d0 )

       throat_radius = radius_bubble * throat_bubble_ratio

       k_1 = 0.125d0 * throat_radius ** 2.d0 * alfa_2 ** tortuosity_factor

       effusive_drag = visc_2 / k_1

       explosive_drag = 3.d0 * c_d / ( 8.d0 * r_a ) * rho_2 * cdabs( u_2 - u_1 ) 

    CASE ( 'forchheimer' )

       ! Forchheimer (Eq. 16 Degruyter et al. 2012)

       IF ( REAL(alfa_2) .GT. 0.0001 ) then

          radius_bubble = ( alfa_2 / ( 4.0 / 3.0 * pi *  bubble_number_density  &
               * ( alfa_1 ) ) ) ** ( 1.d0 / 3.d0 )

          throat_radius = radius_bubble * throat_bubble_ratio

          k_1 = 0.125d0 * throat_radius ** 2.d0 * alfa_2 ** tortuosity_factor

          k_2 = throat_radius / friction_coefficient * alfa_2 ** ( ( 1.d0 +     &
               3.d0 * tortuosity_factor ) / 2.d0 )

          effusive_drag = visc_2 / k_1 + rho_2 / k_2 * cdabs( u_2 - u_1 )

          explosive_drag = 3.d0 * c_d / ( 8.d0 * r_a ) * rho_2                  &
               * cdabs( u_2 - u_1 ) 

       ELSE

          effusive_drag = dcmplx(1.0e20,0.0)

          explosive_drag = dcmplx(0.0,0.0)

       END IF


    CASE ( 'forchheimer_wt' )

       ! Forchheimer (Eq. 16 Degruyter et al. 2012) without transit. zone (phi_t)

       IF ( REAL(alfa_2) .GT. 0.0001 ) then

          radius_bubble = ( alfa_2 / ( 4.0 / 3.0 * pi *  bubble_number_density  &
               * ( alfa_1 ) ) ) ** ( 1.d0 / 3.d0 )

          throat_radius = radius_bubble * throat_bubble_ratio

          k_1 = 0.125d0 * throat_radius ** 2.d0 * alfa_2 ** tortuosity_factor

          k_2 = throat_radius / friction_coefficient * alfa_2 ** ( ( 1.d0 +     &
               3.d0 * tortuosity_factor ) / 2.d0 )

          effusive_drag = visc_2 / k_1 + rho_2 / k_2 * cdabs( u_2 - u_1 )

          explosive_drag = 3.d0 * c_d / ( 8.d0 * r_a ) * rho_2                  &
               * cdabs( u_2 - u_1 ) 

       ELSE

          effusive_drag = dcmplx(1.0e20,0.0)

          explosive_drag = dcmplx(0.0,0.0)

       END IF

    CASE ('drag')

       radius_bubble = ( alfa_2 / ( 4.0 / 3.0 * pi * ( alfa_1 ) ) )             &
            ** ( 1.d0 / 3.d0 )

       rey = 2.d0 * radius_bubble * cdabs( u_2 - u_1 ) / visc_1

       c_d = dreal( 24.0d0 / rey * ( 1.0d0 + 1.0d0 / 8.0d0 * rey**0.72) )

       effusive_drag = 3.d0 * c_d / ( 8.d0 * radius_bubble ) * rho_1            &
            * cdabs( u_2 - u_1 )

       explosive_drag = effusive_drag

    CASE ( 'constant' )

       effusive_drag = dcmplx(1.d0,0.d0)

       explosive_drag = effusive_drag

    CASE ( 'single_velocity' )

       effusive_drag = dcmplx(1.d0,0.d0)

       explosive_drag = effusive_drag

    END SELECT

    IF ( drag_funct_model .EQ. 'forchheimer' ) THEN

       drag_funct = effusive_drag ** ( 1.d0 - t ) *                             &
               ( explosive_drag + 1.d-10 ) ** t
            
    ELSE
    
       drag_funct = effusive_drag ** ( 1.d0 - frag_eff ) *                      &
          ( explosive_drag + 1.d-10 ) ** frag_eff

    END IF

    drag_funct = drag_funct_coeff * drag_funct

    velocity_relaxation = - drag_funct * ( u_1 - u_2 ) * ( rho_mix              &
         / ( rho_1 * rho_2 ) )

    IF ( drag_funct_model .EQ. 'eval' ) THEN

       velocity_relaxation = - drag_funct * ( u_1-u_2 ) / ( x_1 * x_2 * rho_mix )

    ELSE

       velocity_relaxation = - drag_funct * ( u_1 - u_2 ) * ( rho_mix           &
            / ( rho_1 * rho_2 ) )
       
    END IF

    IF ( verbose_level .GE. 3 ) THEN

       WRITE(*,*)  drag_funct , ( u_1 - u_2 ) , x_1 , x_2 , rho_mix 
       WRITE(*,*)  - drag_funct * ( u_1 - u_2 ) / ( x_1 * x_2 * rho_mix ) 
       WRITE(*,*) 'velocity_relaxation',velocity_relaxation

    END IF

  END SUBROUTINE vel_relax_term

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

  SUBROUTINE f_permkc

    IMPLICIT NONE

    REAL*8 :: permkc_alfa_switch
    REAL*8 :: d_permkc_alfa_switch
    REAL*8 :: log_base

    alfa_switch = 0.03

    log_base = 1.d0 / dlog10(dexp(1.d0))

    permkc_alfa_switch = (perm0*10.d0**(-10.2d0*(100.d0*alfa_switch)            &
         **(0.014d0/alfa_switch)) )

    d_permkc_alfa_switch =  log_base * permkc_alfa_switch *                     &
         dlog10(permkc_alfa_switch / perm0) * ( 0.014d0/alfa_switch**2) *       &
         ( 1 - dlog(100*alfa_switch) )

    a_2nd = ( permkc_alfa_switch - d_permkc_alfa_switch * alfa_switch)          &
         / ( alfa_switch**2 - 2.d0 * alfa_switch**2.d0 )

    b_2nd = d_permkc_alfa_switch - 2.d0 * a_2nd*alfa_switch

    IF ( REAL(alfa_2) .LT. alfa_switch ) then

       permkc = a_2nd * alfa_2*alfa_2 + b_2nd * alfa_2

    ELSE

       permkc = perm0*10.d0**(-10.2d0*(100.d0*alfa_2) ** ( 0.014d0 / alfa_2 ) )

    ENDIF

    inv_permkc = 1.d0 / permkc

  END SUBROUTINE f_permkc


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

  SUBROUTINE mixture_viscosity

    USE complexify 
    IMPLICIT NONE

    CALL f_viscliq

    CALL f_bubbles

    visc_mix = visc_1 * visc_rel_bubbles

  END SUBROUTINE mixture_viscosity

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

  SUBROUTINE f_viscliq

    USE complexify 
    IMPLICIT NONE

    ! Calculate viscosity of the melt (bubbles and crystal-free)
    CALL f_viscmelt

    ! Calculate relative visocisty due to crystals
    CALL f_theta

    ! Calculate viscosity of the melt+crystal
    visc_1 = visc_melt * theta

  END SUBROUTINE f_viscliq

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

  SUBROUTINE f_viscmelt

    USE complexify 
    IMPLICIT NONE

    COMPLEX*16 :: w
    COMPLEX*16 :: visc1exp , visc2exp , visc3exp 

    COMPLEX*16, DIMENSION(12) :: wt

    COMPLEX*16, DIMENSION(12) :: norm_wt, xmf
    REAL*8, DIMENSION(12) :: mw
    REAL*8 :: A
    COMPLEX*16 :: gfw , B , C
    REAL*8, DIMENSION(10) :: bb
    COMPLEX*16, DIMENSION(10) :: bcf
    COMPLEX*16, DIMENSION(7) :: cc , ccf
    COMPLEX*16 :: siti , tial , fmm , nak 
    COMPLEX*16 :: b1 , b2 , b3 , b4 , b5 , b6 , b7 , b11 , b12 , b13 , c1 , c2 ,&
         c3 , c4 , c5 , c6 , c11 

    COMPLEX*16 :: x_d_md_tot

    INTEGER :: i

    x_d_md_tot = sum( x_d_md(1:n_gas) )

    w = x_d_md_tot * 100.d0

    SELECT CASE ( visc_melt_model )

    CASE DEFAULT

    CASE ( 'Hess_and_Dingwell1996')

       !  melt viscosity of rhyolite reaches a maximum at 1.e12

       IF ( REAL(w) .LT. 0.03 ) then

          visc_melt = dcmplx( 1.0d12 , 0.d0 )

       ELSE

          visc1exp = (-3.545d0 + 0.833d0 * cdlog(w))

          visc2exp = (9601.0d0 - 2368.0d0 * cdlog(w))

          visc3exp = visc2exp / ( t - ( 195.70d0 + 32.250d0 * cdlog(w) ) )

          visc_melt = 10.0**(visc1exp+visc3exp)

       END IF

    CASE ( 'Romano_et_al2003')

       !  Campi-Flegrei - Romano et al. (2003)

       IF ( REAL(w) .LT. 0.03 ) then

          visc_melt = dcmplx( 1.0d12 , 0.d0 )

       ELSE

          visc1exp = ( -5.8996d0 -0.2857d0 * cdlog(w))

          visc2exp = ( 10775.0d0 - 394.83d0 * cdlog(w))

          visc3exp = visc2exp / ( t - ( 148.71d0 -21.65d0 * cdlog(w) ) )

          visc_melt = 10.0**(visc1exp+visc3exp)

       END IF

    CASE ( 'Giordano_et_al2008' )

       ! -------------------------- Giordano et al. (2008) ----------------------
       ! Call for dissolved water in te melt from the model and add it to the   
       ! composition of glass in wt

       DO i = 1,12

          wt(i) = dcmplx( wt_init(i) , 0.d0 )

       END DO


       wt(11) = w

       ! Create the normalized wt. % distribution

       norm_wt = 100.d0 * ( wt / sum(wt) )

       ! Change to molar fractions

       mw = [ 60.0843 , 79.8658 , 101.961276 , 71.8444 , 70.937449 , 40.3044 ,  &
            56.0774 , 61.97894 , 94.1960 , 141.9446 , 18.01528 , 18.9984 ]
       gfw = 100.d0 / sum( norm_wt / mw )
       xmf = ( norm_wt / mw ) * gfw


       ! Model coefficients

       bb  = [159.56, -173.34, 72.13, 75.69, -38.98, -84.08, 141.54, -2.43,     &
            -0.91, 17.62]
       cc  = [2.75, 15.72, 8.32, 10.2, -12.29, -99.54, 0.3]

       ! Load composition-based matrix for multiplication against 
       ! model-coefficients

       siti = xmf(1) + xmf(2)
       tial = xmf(2) + xmf(3)
       fmm  = xmf(4) + xmf(5) + xmf(6)
       nak  = xmf(8) + xmf(9)
       b1  = siti
       b2  = xmf(3)
       b3  = xmf(4) + xmf(5) + xmf(10)
       b4  = xmf(6)
       b5  = xmf(7)
       b6  = xmf(8) + xmf(11) + xmf(12)
       b7  = xmf(11) + xmf(12) + cdlog( 1.d0 + xmf(11) )
       b11 = siti * fmm
       b12 = ( siti + xmf(3) + xmf(10) ) * ( nak + xmf(11) )
       b13 = xmf(3) * nak

       c1 = xmf(1)
       c2 = tial
       c3 = fmm
       c4 = xmf(7)
       c5 = nak
       c6 = cdlog( 1.d0 + xmf(11) + xmf(12) )
       c11 = xmf(3) + fmm + xmf(7) - xmf(10)
       c11 = c11 * ( nak + xmf(11) + xmf(12) )

       bcf =   [b1, b2, b3, b4, b5, b6, b7, b11, b12, b13]
       ccf =   [c1, c2, c3, c4, c5, c6, c11]   

       ! Model main parameters

       a  = -4.55d0
       b  = sum( bb * bcf )
       c  = sum( cc * ccf )

       ! Calculates viscosity

       visc_melt = a + b / ( t - c )
       visc_melt = 10.d0 ** visc_melt

    CASE ('Di_Genova_et_al2013_eqn_3,5')

       ! Viscosity for peralkaline rhyolites from Pantelleria Island 
       ! Di Genova et al. 2013 (Eqs. 3,5)

       a = -4.55d0

       ! -------- Giordano et al. 2009 ------------------------------------------
       b1 = 4278.17d0
       b2 = 8.6d0
       b = b1 + b2 * w

       c1 = 513.d0
       c2 = -245.3d0
       c = c1 + c2 * cdlog( 1.d0 + w )

       visc_melt = a + b / ( t - c )
       visc_melt = 10.d0 ** visc_melt

    CASE ('Di_Genova_et_al2013_eqn_4,5')

       ! Viscosity for peralkaline rhyolites from Pantelleria Island 
       ! Di Genova et al. 2013 (Eqs. 4,5)

       a = -4.55d0

       ! -------- New parametrization -------------------------------------------
       b3 = 10528.64d0
       b4 = -4672.21d0
       b = b3 + b4 * cdlog( 1.d0 + w )

       c3 = 172.27d0
       c4 = 89.75d0
       c = c3 + c4 * cdlog( 1.d0 + w )


       visc_melt = a + b / ( t - c )
       visc_melt = 10.d0 ** visc_melt


    CASE ( 'Giordano_et_al2009' )

       a = -4.55d0

       b1 = 6101.0d0
       b2 = -63.66d0

       b = b1 + b2 * 3.5d0 * w
       !B = b1 + b2 * w


       c1 = 567.0d0
       c2 = -160.3d0
       c = c1 + c2 * cdlog( 1.d0 + 3.5d0 * w  )
       !C = c1 + c2 * CDLOG( 1.D0 + w  )

       visc_melt = a + b / ( t - c )
       visc_melt = 10.d0 ** visc_melt

    END SELECT


  END SUBROUTINE f_viscmelt

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

  SUBROUTINE f_theta

    USE complexify 
    IMPLICIT NONE

    COMPLEX*16 :: arg_erf , t , errorf, var_phi

    REAL*8 :: omega , betas , theta0

    REAL*8 :: p , a1 , a2 , a3 , a4 , a5
    REAL*8 :: phi_star, csi, delta, Einstein_coeff

    SELECT CASE (theta_model)

    CASE DEFAULT

    CASE ( 'Lejeune_and_Richet1995')

       !---------- Einstein-Roscoe ( with constants from Lejeune & Richet '95)
       theta = ( 1.0d0 - ( sum( beta(1:n_cry) ) / 0.7d0 ) ) ** ( -3.4d0 )

    CASE ('Dingwell1993')

       !---------- Dingwell & al '93
       theta = ( 1.0d0 + 0.75d0 * ( ( sum(beta(1:n_cry) ) / 0.84d0 ) / ( 1.d0 - &
            ( sum( beta(1:n_cry) ) / .84d0 ) ) ) ) **2.d0

    CASE ('Fixed_value')

       !-----------Constant theta
       theta = theta_fixed    ! Melnik & Sparks '99

    CASE ('Costa2005')

       !-----------Costa
       c1 = 0.9995d0
       c2 = 0.4d0
       c3 = 1.d0


       p = 0.3275911d0
       a1 = 0.254829592d0 
       a2 = -0.284496736d0 
       a3 = 1.421413741d0
       a4 = -1.453152027d0
       a5 = 1.061405429d0

       arg_erf = 0.5d0 * dsqrt(pi) * sum( beta(1:n_cry) ) * ( 1.d0 + c2 /       &
            ( 1.d0 - sum( beta(1:n_cry) ) ) ** c3 )

       t = 1.d0 / ( 1.d0 + p * arg_erf )

       errorf = 1.d0 - ( a1 * t + a2 * t**2 + a3 * t**3 + a4 * t**4 + a5 * t**5)&
            * cdexp( - arg_erf ** 2 )

       theta = ( 1.d0 - c1 * errorf ) ** ( - 2.5d0 / c1 )

    CASE ('Melnik_and_Sparks2005')

       !-----------Melnik & Sparks 2005 Eq.16
       c1 = 1.4d0
       c2 = 8.6d0
       c3 = 0.69d0

       theta = c1 * 10.d0 ** ( atan( c2 * ( sum( beta(1:n_cry) ) - c3 ) ) + pi ) 

    CASE ('Melnik_and_Sparks1999')

       omega = 20.6d0
       betas = 0.62d0
       theta0 = 3.5d0

       theta = theta0 * 10.d0 ** ( 0.5d0 * pi + atan( omega * ( sum(            &
            beta(1:n_cry) ) - betas )))

    CASE ('Vona_et_al2011')

       !-----------Costa et al. 2009, Vona et al. 2011

       p = 0.3275911d0
       a1 = 0.254829592d0 
       a2 = -0.284496736d0 
       a3 = 1.421413741d0
       a4 = -1.453152027d0
       a5 = 1.061405429d0

       ! This is the value reported in Vona et al 2011, but according to Fig 9
       ! of that article, a best fitting is obtained with 0.29

       phi_star = 0.274d0
       csi = 0.0327d0
       delta = 13.d0 - 0.84d0
       einstein_coeff = 2.8d0


       var_phi = sum( beta(1:n_cry) * rho_c(1:n_cry) / rho_1 ) / phi_star

       arg_erf = 0.5d0 * dsqrt(pi) / (1.0d0 - csi) * var_phi * ( 1.d0 +         &
            var_phi ** 0.84d0 ) 

       t = 1.d0 / ( 1.d0 + p * arg_erf )

       errorf = 1.d0 - ( a1 * t + a2 * t**2.0d0 + a3 * t**3.0d0                 &
            + a4 * t**4.0d0 + a5 * t**5.0d0) * cdexp( - arg_erf ** 2.0d0 )

       theta =  theta_fixed * (1.d0 + var_phi ** delta ) /                      &
            ( ( 1.d0 - (1.d0 - csi) * errorf ) ** (einstein_coeff * phi_star) ) 


    CASE ('Vona_et_al2011_mod')

       !-----------Costa et al. 2009, Vona et al. 2011

       p = 0.3275911d0
       a1 = 0.254829592d0 
       a2 = -0.284496736d0 
       a3 = 1.421413741d0
       a4 = -1.453152027d0
       a5 = 1.061405429d0

       ! Modification of Vona et al. 2011 in order to better reproduce
       ! the viscosity at Stromboli.

       phi_star = 0.39d0
       csi = 0.03d0
       delta = 2.0d0 - 0.84d0
       einstein_coeff = 2.8d0


       var_phi = sum( beta(1:n_cry) * rho_c(1:n_cry) / rho_1 ) / phi_star

       arg_erf = 0.5d0 * dsqrt(pi) / (1.0d0 - csi) * var_phi * ( 1.d0 +         &
            var_phi ** 0.84d0 ) 

       t = 1.d0 / ( 1.d0 + p * arg_erf )

       errorf = 1.d0 - ( a1 * t + a2 * t**2.0d0 + a3 * t**3.0d0 + a4 * t**4.0d0 &
            + a5 * t**5.0d0) * cdexp( - arg_erf ** 2.0d0 )

       theta =  theta_fixed * (1.d0 + var_phi ** delta )/( ( 1.d0 - (1.d0 - csi)&
            * errorf ) ** (einstein_coeff * phi_star) ) 
            
    CASE ('Vona_et_al2013_eq19')

       theta = ( 1.d0 - sum(beta(1:n_cry)) / (1.d0 - alfa_2 ) ) ** ( - 5.d0     &
            / 2.d0) * ( 1 - alfa_2 ) ** (- 1.d0)

    CASE ('Vona_et_al2013_eq20')

       theta = ( 1.d0 - sum(beta(1:n_cry)) - alfa_2 )  ** ( - ( 5.d0 *          &
            sum(beta(1:n_cry)) + 2.d0 * alfa_2 ) / ( 2.d0 * ( sum(beta(1:n_cry))&
            + alfa_2 ) ) )

    CASE ('Vona_et_al2013_eq21')

       theta = ( 1.d0 - alfa_2 / (1.d0 - sum(beta(1:n_cry)) ) ) ** ( -1.0 )     &
           * ( 1 - sum(beta(1:n_cry)) ) ** (- 5.d0/2.d0)

    END SELECT

  END SUBROUTINE f_theta

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

  SUBROUTINE f_bubbles

    USE complexify 
    USE geometry, ONLY : radius
    IMPLICIT NONE

    REAL*8 :: Ca, gamma_Ca

    IF(  (.NOT. (theta_model .EQ. 'Vona_et_al2013_eq19' ) ) .AND.               & 
         (.NOT. (theta_model .EQ. 'Vona_et_al2013_eq20' ) ) .AND.               & 
         (.NOT. (theta_model .EQ. 'Vona_et_al2013_eq21' ) ) ) THEN

       SELECT CASE ( bubbles_model )

       CASE DEFAULT

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0)

       CASE ( 'none' )

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0)

       CASE ( 'Costa2007' )

          gamma_ca = 0.25d0

          ca =  r_a * visc_melt * (8.00d0 * (u_mix * pi * radius**2)            &
               / (3.0d0 * pi * radius**3.0d0)) / gamma_ca

          visc_rel_bubbles = (1.0d0/(1.0d0 + 25.0d0*ca*ca))                     & 
               * ( 1.0d0 / (1.0d0-alfa_2)                                       &
               + 25.0d0*ca*ca*((1.0d0-alfa_2)**(5.0d0/3.0d0)))


       CASE ( 'Einstein' )

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) / ( 1.d0 - alfa_2 )

       CASE ( 'Quane-Russel' )

          ! For Campi-Flegrei we use 0.63 as reported in 2004 Report (Task2.2)

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * cdexp( ( - 0.63 * alfa_2 )   &
               / ( 1.d0 - alfa_2 ) )

       CASE ( 'Eilers' )

          ! Eq. (17) Mader et al. 2013

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * (1.0d0 + (1.25d0 * alfa_2)   &
               / (1.0d0 - 1.29 * alfa_2) ) ** 2.0d0

       CASE ( 'Sibree' )

          ! Eq. (18) Mader et al. 2013

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * ( 1.0d0 / ( 1.0d0            &
               - (1.2 * alfa_2)** 0.33333d0 ) )

       CASE ( 'Taylor' )

          ! Eq. (16) Mader et al. 2013

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * (1.0d0 + alfa_2)

       CASE ( 'Mackenzie' )

          ! Eq. (19) Mader et al. 2013

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * (1.0d0 - (5.d0 / 3.d0)       &
               * alfa_2)

       CASE ( 'DucampRaj' )

          ! Eq. (21) Mader et al. 2013, using b = 3

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * cdexp( -3.d0 * ( alfa_2      &
               / ( 1.d0 - alfa_2 ) ) )

       CASE ( 'BagdassarovDingwell' )

          ! Eq. (22) Mader et al. 2013

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * ( 1.d0 / (1.d0 + 22.4d0      &
               * alfa_2 ) )

       CASE ( 'Rahaman' )

          ! Eq. (20) Mader et al. 2013

          visc_rel_bubbles = dcmplx(1.0d0,0.0d0) * cdexp( - 11.2d0 * alfa_2 )

       END SELECT

    ELSE

       visc_rel_bubbles = dcmplx(1.0d0,0.0d0)

    END IF

  END SUBROUTINE f_bubbles

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

  SUBROUTINE f_alfa(xtot,xmax,r_beta,r_rho_md,r_rho_2,r_alfa_2)

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: xtot(n_gas) , xmax(n_gas) 
    REAL*8, INTENT(IN) :: r_beta(n_cry)
    REAL*8, INTENT(IN) :: r_rho_md
    REAL*8, INTENT(IN) :: r_rho_2 
    REAL*8, INTENT(OUT) :: r_alfa_2(n_gas)


    REAL*8 :: r_x_g(n_gas)

    ! Mass fraction of the exsolved gas with respect to the crystal-free magma
    r_x_g(1:n_gas) = xtot(1:n_gas) - xmax(1:n_gas)

    r_alfa_2 = ( r_x_g * ( 1.d0 - sum( r_beta(1:n_cry) ) ) * r_rho_md ) /       &
         ( ( 1.d0 - xtot) * r_rho_2 + r_x_g * ( 1.d0 - sum( r_beta(1:n_cry) ) ) &
         * r_rho_md )

  END SUBROUTINE f_alfa

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

  SUBROUTINE f_alfa3(r_p_2,xtot,r_beta,r_rho_md,r_rho_g,r_alfa_g)

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: xtot(n_gas) 
    REAL*8, INTENT(IN) :: r_beta(n_cry)
    REAL*8, INTENT(IN) :: r_p_2
    REAL*8, INTENT(IN) :: r_rho_md
    REAL*8, INTENT(IN) :: r_rho_g(n_gas) 
    REAL*8, INTENT(OUT) :: r_alfa_g(n_gas)


    !REAL*8 :: r_x_g(n_gas)    !< exsolved gas mass fraction
    REAL*8 :: r_x_d(n_gas)
    REAL*8 :: r_alfa_g_2_max(n_gas),best_alfa_g_2(n_gas)
    REAL*8 :: r_alfa_g_2(n_gas), r_alfa_g_2_new(n_gas)
    REAL*8 :: h,error,best_error
    REAL*8 :: A(n_gas,n_gas),num(n_gas), den(n_gas)
    REAL*8 :: r_rho_md_beta
    INTEGER :: iter,i,j
    INTEGER :: pivot(n_gas)
    INTEGER :: ok

    r_alfa_g = 1.0e-5;

    r_alfa_g_2_max = (xtot / solub)**(1.0 / solub_exp) / r_p_2;


    r_rho_md_beta = ( 1.d0 - sum( r_beta(1:n_cry) ) ) * r_rho_md

    h = 1e-5;

    iter = 1
    error = 1.0
    r_alfa_g_2 = [0.0 , 1.0]

    best_error = 1.0
    best_alfa_g_2 = r_alfa_g_2

    DO WHILE ( r_alfa_g_2(1) .LT. r_alfa_g_2_max(1) - h)

       r_alfa_g_2 = [r_alfa_g_2(1) + h, 1.0 - (r_alfa_g_2(1) + h)]

       DO i = 1,n_gas

          r_x_d(i) = dreal( solub(i) * (r_alfa_g_2(i) * p_2)**(solub_exp(i)) )

          DO j = 1,n_gas
             a(i,j) = - xtot(i) * r_rho_g(j) + r_rho_md_beta * &
                  ( xtot(i) - r_x_d(i) ) 
          END DO

          a(i,i) = a(i,i) + r_rho_g(i)

          r_alfa_g(i) = r_rho_md_beta * ( xtot(i) - r_x_d(i) )

       END DO

       call dgesv(n_gas, 1, a , n_gas, pivot, r_alfa_g , n_gas, ok)

       r_alfa_g_2_new = r_alfa_g / sum( r_alfa_g )

       error = maxval(abs(r_alfa_g_2_new - r_alfa_g_2))

       DO i = 1,n_gas
          IF (r_alfa_g(i) .LT. 0.0) THEN
             error = 1.0
          END IF
       END DO

       DO i = 1,n_gas
          IF (r_x_d(i) .GT. xtot(i)) THEN
             error = 1.0
          END IF
       END DO

       IF ( error .LT. best_error) THEN
          best_error = error
          best_alfa_g_2 = r_alfa_g_2
       END IF

    END DO

    r_alfa_g_2 = best_alfa_g_2

    DO i = 1,n_gas

       r_x_d(i) = dreal(solub(i) * (r_alfa_g_2(i) * p_2)**(solub_exp(i)) )

       DO j = 1,n_gas
          a(i,j) = - xtot(i) * r_rho_g(j) + r_rho_md_beta *                     &
               ( xtot(i) - r_x_d(i) ) 
       END DO

       a(i,i) = a(i,i) + r_rho_g(i)

       r_alfa_g(i) = r_rho_md_beta * ( xtot(i) - r_x_d(i) )

    END DO

    call dgesv(n_gas, 1, a , n_gas, pivot, r_alfa_g , n_gas, ok)

    num = r_alfa_g * r_rho_g + ( 1.0 - sum( r_alfa_g ))* r_rho_md_beta * r_x_d
    den = sum(r_alfa_g * r_rho_g) + ( 1.0 - sum( r_alfa_g )) * r_rho_md_beta
    error = maxval(abs(xtot - num/den));

    IF( error .GT. 1e-010 ) THEN
       WRITE(*,*) 'Initial exsolved volatile content error!'
       WRITE(*,*) 'error='
       WRITE(*,*) error
       stop
    END IF


  END SUBROUTINE f_alfa3

END MODULE constitutive