mixture.f90

Go to the documentation of this file.

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

MODULE mixture_module

  USE envi_module, ONLY: alpha, t_a, gi

  IMPLICIT NONE

  REAL*8 :: C_vmix

  REAL*8 :: epsilon

  REAL*8 :: Nmag

  REAL*8 :: N
  
  REAL*8 :: init_N

  REAL*8 :: lambda

  REAL*8 :: rhogas

  REAL*8 :: beta_old

  REAL*8 :: beta_new

  REAL*8 :: Ri
  
  REAL*8 :: flow_regime

  REAL*8 :: gas_constmix

  REAL*8 :: C_s = 1617.D0       

  REAL*8 :: beta

  REAL*8 :: rhosol_ave
  
  REAL*8 :: rwvapour
       
  REAL*8 :: cpwvapour
  
  REAL*8 :: pi
  
  REAL*8 :: initial_velocity
  
  REAL*8 :: initial_density
  
  REAL*8 :: Cv_magmix
  
  REAL*8 :: solidvolumefraction
  REAL*8 :: gasvolumefraction
  REAL*8 :: final_solid_mass_flux
  REAL*8 :: initial_SMF
  
  REAL*8 :: deltabeta
  
  REAL*8 :: deltabeta_g
  
  REAL*8 :: Pdyn
  
  REAL*8 :: entrainment_rate

  SAVE

CONTAINS

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

  SUBROUTINE compute_mixture
    !
    USE envi_module, ONLY: alpha ,p, c_vair, gas_constair
    USE current_module, ONLY: u, r, h, t, solid_mass_flux, t0, tvent_flag ,     &
         oned_model, initial_mf
    USE particles_module, ONLY: rhosol, diam, iclass, fracsolid , v_s, s, c_d , &
         sumsed ,solidmassflux_fract, acc_rate
    !      
    IMPLICIT NONE

    INTEGER :: i        
    
    REAL*8 :: var
    
    ! mass fraction of magmatic gas in the mixture (with air and solid)
    REAL*8 :: nmag_mix

    ! mass fraction pf magmatic gas in the gas (volcanic+air)
    REAL*8 :: nmag_gas
    
    pi = 4.d0 * atan(1.d0)
    
    init_n = n

    nmag_mix = nmag * ( 1.d0 - n ) / (1.d0 - nmag)

    lambda = (n - nmag) / ( 1.d0 - nmag )

    c_vmix = nmag_mix * cpwvapour + lambda * c_vair                   &
         + ( 1.d0 - n ) * c_s 
         
    cv_magmix = nmag * cpwvapour + ( 1.d0 - nmag ) * c_s

    nmag_gas = nmag_mix / n

    gas_constmix = rwvapour * nmag_gas + gas_constair * ( 1.d0 - nmag_gas ) 
    
  

    IF (tvent_flag) THEN

       ! Calculating the initial temperature using temperature of erupted material 
       ! if not explicitly defined in input file 
       ! Initial temperature in the flow based on the eruption temp (eq. 14 in 
       ! Bursik and Woods 96)    

       t = ( ( 1.d0 - lambda ) * cv_magmix * t0 + lambda * c_vair * t_a ) /        &
            ( ( 1.d0 - lambda )* cv_magmix + lambda * c_vair )
        
    ELSE   

       t = t0

    END IF

    rhogas = p / ( gas_constmix * t )
    rhosol_ave = 1/(sum(fracsolid/rhosol))
    beta = 1.d0 / ( n / rhogas + ( 1.d0 - n ) / rhosol_ave)
    
    initial_density = beta
  
    ! Flow velocity
    u = sqrt((beta - alpha) * gi * h / (ri * beta))
    initial_velocity = u

    ! Defining entrainment conditions          
     epsilon = 0.075d0 / dsqrt( 1.0d0 + 718.0d0 * ri**2.4d0 )
     
    IF ( oned_model ) THEN

       var = 1.d0

    ELSE

       var = r

    END IF

    solid_mass_flux = (beta * u * h * var) * ( 1.d0 - n ) 
    initial_smf = (beta * u * h * var) * ( 1.d0 - n ) 
    
    initial_mf =  (beta * u * h * var) * 2 * pi
 
    ! provides the mass_flux of a particle of a given size
    DO i=1,iclass

       solidmassflux_fract(i) = solid_mass_flux*fracsolid(i)  

       ! Settling velocity for the particle fraction

       v_s(i) = dsqrt( (rhosol(i) * gi * diam(i)) / (c_d(i) * beta) )  

       ! SEDIMENTATION for the particle fraction

       s(i) = 1.d0 / (h * u) * solidmassflux_fract(i) * v_s(i)

    ENDDO

    sumsed = sum(s)

    RETURN
  END SUBROUTINE compute_mixture

  !----------------------------------------------------------------------
END MODULE mixture_module
!----------------------------------------------------------------------