inpout.f90

Go to the documentation of this file.

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

MODULE inpout

  ! -- Variables for the namelist TRANSIENT_PARAMETERS
  USE parameters, ONLY : verbose_level
  USE parameters, ONLY : n_cry , n_gas , n_eqns , n_vars , n_mom

  ! -- Variables for the namelist NEWRUN_PARAMETERS
  USE geometry, ONLY : z0 , zn , radius_fixed, radius_min, radius_max,          &
                       radius_z, radius_z_sig, radius_model , comp_cells ,      &
                       update_radius

  ! -- Variables for the namelist STEADY_BOUNDARY_CONDITIONS
  USE init, ONLY : t_in , p1_in , delta_p_in , p_out , u1_in
  USE constitutive, ONLY : x_ex_dis_in
  USE parameters, ONLY : shooting
  USE parameters, ONLY : eps_conv

  ! -- Variables for the namelist EXSOLVED_GAS_PARAMETERS
  USE constitutive, ONLY : rho0_g , cv_g , gamma_g , t0_g ,                     &
       bar_e_g , visc_2 , perm0, pc_g, &
       tc_g, a_g, b_g, s0_g, gas_law

  USE equations, ONLY : lateral_degassing_flag , alfa2_lat_thr

  ! -- Variables for the namelist DISSOLVED_GAS_PARAMETERS
  USE constitutive, ONLY : rho0_d , c0_d , cv_d , gamma_d , p0_d , t0_d ,       &
       bar_e_d , bar_p_d , s0_d , exsol_model , solub , solub_exp

  ! -- Variables for the namelist CRYSTALS_PARAMETERS
  USE constitutive, ONLY : rho0_c , c0_c , cv_c , gamma_c , p0_c , t0_c ,       &
       bar_e_c , bar_p_c , s0_c , beta0 , beta_max, crystallization_model

  ! -- Variables for the namelist MELT_PARAMETERS
  USE constitutive, ONLY : rho0_m , c0_m , cv_m , gamma_m , p0_m , t0_m ,       &
       bar_e_m , bar_p_m , s0_m

  ! -- Variables for the namelist VISCOSITY_PARAMETERS
  USE constitutive, ONLY : visc_melt_model , bubbles_model , theta_model ,      &
       theta_fixed

  ! -- Variables for the namelist TEMPERATURE_PARAMETERS
  USE equations, ONLY : isothermal , fixed_temp

  ! -- Variables for the namelist FRAGMENTATION_PARAMETERS
  USE constitutive, ONLY : explosive , fragmentation_model , frag_thr

  ! -- Variables for the namelist EXTERNAL_WATER_PARAMETERS
  USE equations, ONLY : ext_water , total_water_influx , min_z_influx ,      &
       delta_z_influx , t_w , inst_vaporization , aquifer_type

  ! -- Variables for the namelist SOURCE_PARAMETERS
  USE constitutive, ONLY : grav

  ! -- Variables for the namelist COUNTRY_ROCK_PARAMETERS
  USE constitutive, ONLY : rho_cr, k_cr

  ! -- Variables for the namelist RELAXATION_PARAMETERS
  USE constitutive, ONLY : drag_funct_model , drag_funct_coeff , p_relax_model ,&
       tau_p_coeff , tau_c , tau_d

  ! -- Variables for the namelist FORCHHEIMER_PARAMETERS
  USE constitutive, ONLY : bubble_number_density , tortuosity_factor ,          &
       throat_bubble_ratio , friction_coefficient , c_d , r_a, xa, xb, xc,       &
       log10_bubble_number_density 

  !  USE parameters, ONLY : atmospheric_pressure, chocked_flow

  ! -- Variables for the card meltcomposition
  USE constitutive, ONLY : wt_init

  IMPLICIT NONE

  CHARACTER(LEN=50) :: run_name
  CHARACTER(LEN=50) :: bak_name
  CHARACTER(LEN=50) :: input_file
  CHARACTER(LEN=50) :: output_q_file
  CHARACTER(LEN=50) :: output_p_file
  CHARACTER(LEN=50) :: steady_p_file
  CHARACTER(LEN=50) :: steady_q_file
  CHARACTER(LEN=50) :: exit_file

  INTEGER, PARAMETER :: input_unit = 7
  INTEGER, PARAMETER :: backup_unit = 8
  INTEGER, PARAMETER :: output_q_unit = 10
  INTEGER, PARAMETER :: output_p_unit = 11
  INTEGER, PARAMETER :: steady_p_output_unit = 13
  INTEGER, PARAMETER :: steady_q_output_unit = 14
  INTEGER, PARAMETER :: dakota_unit = 15
  INTEGER, PARAMETER :: dakota_unit2 = 17
  INTEGER, PARAMETER :: exit_unit = 16


  LOGICAL :: close_units

  REAL*8 :: log10_k_cr

  ! -- Variables for the namelist RELAXATION_PARAMETERS
  REAL*8 :: log10_drag_funct_coeff , log10_tau_p_coeff

  REAL*8, ALLOCATABLE :: log10_tau_c(:)
  REAL*8, ALLOCATABLE :: log10_tau_d(:)

  namelist / run_parameters / run_name , verbose_level

  namelist / geometry_parameters / z0 , zn , radius_model , radius_fixed ,      &
       radius_min, radius_max, radius_z, radius_z_sig,  comp_cells

  namelist / phases_parameters / n_gas , n_cry , n_mom

  namelist / steady_boundary_conditions / t_in , p1_in , delta_p_in ,           &
       x_ex_dis_in , p_out , u1_in , eps_conv, shooting

  namelist / exsolved_gas_parameters / gas_law, pc_g , tc_g , cv_g , gamma_g ,  &
       rho0_g , t0_g , bar_e_g , s0_g, visc_2 , lateral_degassing_flag ,        &
       alfa2_lat_thr , perm0
  namelist / dissolved_gas_parameters / rho0_d , c0_d , cv_d , gamma_d , p0_d , &
       t0_d , bar_e_d , bar_p_d , s0_d , exsol_model , solub , solub_exp

  namelist / crystals_parameters / rho0_c , c0_c , cv_c , gamma_c , p0_c ,      &
       t0_c , bar_e_c , bar_p_c , s0_c , beta0 , beta_max, crystallization_model

  namelist / melt_parameters / rho0_m , c0_m , cv_m , gamma_m , p0_m , t0_m ,   &
       bar_e_m , bar_p_m , s0_m

  namelist / viscosity_parameters / visc_melt_model , bubbles_model ,           &
       theta_model , theta_fixed

  namelist / temperature_parameters / isothermal , fixed_temp

  namelist / fragmentation_parameters / explosive , fragmentation_model ,       &
       frag_thr

  namelist / external_water_parameters / ext_water , total_water_influx ,       &
       min_z_influx , delta_z_influx , t_w , inst_vaporization , aquifer_type

  namelist / source_parameters /  grav

  namelist / country_rock_parameters / rho_cr , log10_k_cr
 
  namelist / relaxation_parameters / drag_funct_model , log10_drag_funct_coeff ,&
       p_relax_model , log10_tau_p_coeff, log10_tau_c , log10_tau_d

  namelist / bubbles_parameters / bubble_number_density

  namelist / forchheimer_parameters / log10_bubble_number_density ,             &
       tortuosity_factor , throat_bubble_ratio , friction_coefficient , c_d ,   &
       r_a

  namelist / permeability_parameters / xa, xb, xc


CONTAINS

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

  SUBROUTINE init_param

    IMPLICIT none

    LOGICAL :: lexist

    INTEGER :: n_gas_init

    INTEGER :: n_cry_init

    ! -- Variables for the namelist DISSOLVED_GAS_PARAMETERS
    REAL*8 :: rho0_d_init, C0_d_init , cv_d_init , gamma_d_init , p0_d_init , &
         T0_d_init , bar_e_d_init , bar_p_d_init , s0_d_init , solub_init ,   &
         solub_exp_init

    ! -- Variables for the namelist CRYSTALS_PARAMETERS
    REAL*8 :: rho0_c_init , C0_c_init , cv_c_init , gamma_c_init , p0_c_init ,&
         T0_c_init , bar_e_c_init , bar_p_c_init , s0_c_init , beta0_init ,   &
         beta_max_init

    REAL*8 :: log10_tau_c_init , log10_tau_d_init

    REAL*8 :: x_ex_dis_in_init, xa, xb, xc

    REAL*8 :: Pc_g_init , Tc_g_init , cv_g_init , gamma_g_init , rho0_g_init ,&
         T0_g_init , bar_e_g_init, s0_g_init

    CHARACTER*20 :: gas_law_init 

    namelist / phases_parameters_init / n_gas_init , n_cry_init

    namelist / steady_boundary_conditions_init / t_in , p1_in , delta_p_in ,  &
         x_ex_dis_in_init , p_out , u1_in , eps_conv, shooting

    namelist / exsolved_gas_parameters_init / gas_law_init, pc_g_init ,       &
         tc_g_init , cv_g_init , gamma_g_init , rho0_g_init , t0_g_init ,     &
         bar_e_g_init , s0_g_init, visc_2 , lateral_degassing_flag ,          &
         alfa2_lat_thr , perm0

    namelist / dissolved_gas_parameters_init / rho0_d_init ,                  &
         c0_d_init , cv_d_init , gamma_d_init , p0_d_init , t0_d_init ,       &
         bar_e_d_init , bar_p_d_init , s0_d_init , exsol_model ,              &
         solub_init , solub_exp_init

    namelist / crystals_parameters_init / rho0_c_init ,                       &
         c0_c_init , cv_c_init , gamma_c_init , p0_c_init , t0_c_init ,       &
         bar_e_c_init , bar_p_c_init , s0_c_init , beta0_init , beta_max_init


    namelist / relaxation_parameters_init / drag_funct_model ,                &
         log10_drag_funct_coeff , p_relax_model , log10_tau_p_coeff ,         &
         log10_tau_c_init , log10_tau_d_init


    !-- Inizialization of the Variables for the namelist RUN_PARAMETERS
    run_name = 'test'
    verbose_level = 4 

    !-- Inizialization of the Variables for the namelist geometry_parameters
    z0 = 0.d0
    zn = 5200.d0
    radius_fixed = 0.d0
    radius_min = 0.d0
    radius_max = 0.d0
    radius_z = 0.d0
    radius_z_sig = 0.d0
    radius_model = 'fixed'
    comp_cells = 500

    !-- Inizialization of the Variables for the namelist 
    !-- steady_boundary_condition_parameters
    t_in = 1323.d0
    p1_in = 135000000.d0
    delta_p_in = 0.d0
    x_ex_dis_in_init = 5.0d-2
    p_out = 101300
    u1_in = 0.d0
    shooting = .true.
    eps_conv = 1.d-5

    ! exsolved gas parameters
    gas_law_init = 'VDW'
    pc_g_init = 22064000.d0
    tc_g_init = 647.d0
    cv_g_init = 1571.d0
    gamma_g_init = 1.324d0
    rho0_g_init = 0.52924425d0
    t0_g_init = 373.0
    s0_g_init = 2373.0
    visc_2 = 1.5d-2
    lateral_degassing_flag = .false.
    alfa2_lat_thr = 1.1d0
    perm0 = 5.0d-3

    n_gas_init = 1
    n_cry_init = 1

    n_mom = 1
    
    n_vars = 5 + 2 * n_gas_init + n_cry_init * n_mom
    n_eqns = n_vars

    ! dissolved gas parameters
    rho0_d_init = 1000.d0
    c0_d_init = 2000.d0
    cv_d_init = 1200.d0
    gamma_d_init = 2.3d0
    p0_d_init = 1.0d8
    s0_d_init = 0.d0
    exsol_model = 'Henry'
    solub_init = 4.11d-06 
    solub_exp_init = 0.5d0

    ! crystals parameters
    rho0_c_init = 2600.d0
    c0_c_init = 2000.d0
    cv_c_init = 1200.d0
    gamma_c_init = 2.3d0
    p0_c_init = 1.0d8
    s0_c_init = 0.d0
    beta0_init = 0.0d0
    beta_max_init = 0.6d0
    crystallization_model = 'Vitturi2010'

    ! MELT_PARAMETERS
    rho0_m = 2300.0000000000000  
    c0_m = 2000.0000000000000    
    cv_m = 1200.0000000000000    
    gamma_m = 2.3    
    p0_m = 100000000.00000000  
    s0_m = 0.0000000000000000  

    !-- Inizialization of the Variables for the namelist viscosity_parameters
    visc_melt_model = 'Hess_and_Dingwell1996'
    bubbles_model = 'none'
    theta_model = 'Fixed_value'
    theta_fixed = 50.d0

    !-- Inizialization of the Variables for the namelist temperature_parameters
    isothermal = .false.
    fixed_temp = 1323.d0

    !-- Inizialization of the Variables for the namelist 
    !-- fragmentation_parameters
    explosive = .true.
    fragmentation_model = 1
    frag_thr = 0.60d+0
    !tau_frag_coeff = 1.D0
    !tau_frag_exp = 5.0D0

    !-- Inizialization of the Variables for the namelist 
    !-- external_water_parameters 
    ext_water = .false. 
    total_water_influx = 0.d0
    min_z_influx = 0.d0
    delta_z_influx = 0.d0
    t_w = 0.d0
    inst_vaporization = .false.
    aquifer_type = "unconfined"
    
    !-- Inizialization of the Variables for the namelist source_parameters
    grav = 9.81d0


    !-- Inizialization of the Variables for the namelist country_rock_parameters
    rho_cr = 2600
    log10_k_cr = -12


    !-- Initialization of the variables for the namelist relaxation_parameters
    drag_funct_model = 'forchheimer'
    log10_drag_funct_coeff = 1.0d0
    p_relax_model = 'constant'
    log10_tau_p_coeff = -8.d0
    log10_tau_d_init = -4.0d0
    log10_tau_c_init = -4.0d0

    drag_funct_coeff = 10.d0 ** log10_drag_funct_coeff
    tau_p_coeff = 10.d0 ** log10_tau_p_coeff
    tau_d(1:n_gas) = 10.d0 ** log10_tau_d
    tau_c(1:n_cry) = 10.d0 ** log10_tau_c(1:n_cry)

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

    log10_bubble_number_density = 15.0d0   
    bubble_number_density = 10.0d0 ** log10_bubble_number_density
    tortuosity_factor = 3.5 
    throat_bubble_ratio = 0.1 
    friction_coefficient = 10.d0

    c_d = 0.8d0
    r_a = 1.d-3

    xa = 1.0
    xb = 1.0
    xc = 1.0

    input_file = 'conduit_solver.inp'

    INQUIRE (file=input_file,exist=lexist)

    IF ( .NOT. lexist ) THEN

       OPEN(input_unit,file=input_file,status='NEW')

       WRITE(input_unit, run_parameters )

       WRITE(input_unit, geometry_parameters )

       WRITE(input_unit, phases_parameters_init )

       WRITE(input_unit, steady_boundary_conditions_init )

       WRITE(input_unit, exsolved_gas_parameters_init )

       WRITE(input_unit, dissolved_gas_parameters_init )

       WRITE(input_unit, crystals_parameters_init )

       WRITE(input_unit, melt_parameters )

       WRITE(input_unit, viscosity_parameters )

       WRITE(input_unit, temperature_parameters )

       WRITE(input_unit, fragmentation_parameters )

       WRITE(input_unit, external_water_parameters )

       WRITE(input_unit, source_parameters )

       IF ( lateral_degassing_flag ) THEN

          WRITE(input_unit, country_rock_parameters )

       END IF

       WRITE(input_unit, relaxation_parameters_init )

       IF ( drag_funct_model .EQ. 'eval' ) THEN

          WRITE(input_unit, bubbles_parameters )

       ELSEIF ( ( drag_funct_model .EQ. 'darcy' ) .OR.                          &
            drag_funct_model .EQ. 'forchheimer' .OR.                            &
            drag_funct_model .EQ. 'forchheimer_wt') THEN

          WRITE(input_unit, forchheimer_parameters )

       END IF

       CLOSE(input_unit)

       WRITE(*,*) 'Input file not found'
       WRITE(*,*) 'A new one with default values has been created'
       stop

    ELSE


    END IF

  END SUBROUTINE init_param

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

  SUBROUTINE read_param

    ! external subroutines
    USE constitutive, ONLY : allocate_phases_parameters

    ! external variables
    USE constitutive, ONLY : n_drag_models , available_drag_models
    USE constitutive, ONLY : n_theta_models , available_theta_models
    USE constitutive, ONLY : n_bubble_models, available_bubble_models
    USE constitutive, ONLY : n_visc_melt_models , available_visc_melt_models
    USE constitutive, ONLY : t0_c , bar_p_c !, bar_e_c
    USE constitutive, ONLY : t0_m , bar_p_m !, bar_e_m

    USE init, ONLY : beta_in , xd_md_in
    
    IMPLICIT none

    INTEGER :: i
    LOGICAL :: tend1
    CHARACTER(LEN=80) :: card

    LOGICAL :: check_model
    
    INTEGER :: ios

    OPEN(input_unit,file=input_file,status='old')


    ! ------- READ run_parameters NAMELIST --------------------------------------

    READ(input_unit, run_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist RUN_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF
    
    READ(input_unit,geometry_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist GEOMETRY_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF
    
    IF ( radius_model == 'fixed' ) THEN

       IF ( radius_min .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_min not used'
       IF ( radius_max .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_max not used'
       IF ( radius_z .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_z not used'
       IF ( radius_z_sig .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_z_sig not used'

    ELSEIF ( radius_model == 'linear' ) THEN

       IF ( radius_fixed .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_fixed not used'
       IF ( radius_z .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_z not used'
       IF ( radius_z_sig .NE. 0.d0 ) WRITE(*,*) 'WARNING: radius_z_sig not used'

    END IF
              

    ! ------- READ phases_parameters NAMELIST -----------------------------------
    READ(input_unit, phases_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist PHASES_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    IF ( n_mom .LE. 1 ) THEN

       WRITE(*,*) 'Solving for crystal volume fraction only'

    ELSE

       WRITE(*,*) 'Solving for ',n_mom,' moments for each crystal phase'

    END IF
    
    n_vars = 5 + 2 * n_gas + n_cry * n_mom
    n_eqns = n_vars

    ALLOCATE( beta_in(n_cry) )
    ALLOCATE( xd_md_in(n_gas) )

    CALL allocate_phases_parameters

    READ(input_unit,steady_boundary_conditions , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist STEADY_BOUNDARY_CONDITIONS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    IF ( u1_in .GT. 0.d0 ) THEN

       IF ( shooting .EQV. .true.) THEN
          WRITE(*,*) 'Shooting technique to search for inlet velocity'
       ELSE
          WRITE(*,*) 'Single run: no shooting'
       END IF

    ELSE

       shooting = .true.
       WRITE(*,*) 'Shooting technique to search for inlet velocity'

    END IF


    ! ------- READ exsolved_gas_parameters NAMELIST ----------------------------
    READ(input_unit, exsolved_gas_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist EXSOLVED_GAS_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    IF ( gas_law .EQ. 'IDEAL' )  THEN

       a_g(1:n_gas) = 0.d0
       b_g(1:n_gas) = 0.d0

    ELSEIF (  gas_law .EQ. 'VDW' ) THEN

       a_g(1:n_gas) = 27.0 / 64.0 * ( cv_g(1:n_gas)*(gamma_g(1:n_gas) - 1)      &
            * tc_g(1:n_gas))**2.0 / pc_g(1:n_gas)

       b_g(1:n_gas) = 1.0 / 8.0 * ( cv_g(1:n_gas)*(gamma_g(1:n_gas) - 1)        &
            * tc_g(1:n_gas)) / pc_g(1:n_gas)

    ELSE 

       WRITE(*,*) ''
       WRITE(*,*) 'Wrong gas law chosen.'
       WRITE(*,*) 'Please choose between:'
       WRITE(*,*) ''
       WRITE(*,*) 'IDEAL'
       WRITE(*,*) 'VDW'
       WRITE(*,*) ''
       CALL abort

    END IF

    ! ------- READ dissolved_gas_parameters NAMELIST ---------------------------
    READ(input_unit, dissolved_gas_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist DISSOLVED_GAS_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    t0_d(1:n_gas) = c0_d(1:n_gas) **2.d0 / ( cv_d(1:n_gas) * gamma_d(1:n_gas)   &
         * ( gamma_d(1:n_gas) - 1.d0 ) )

    bar_p_d(1:n_gas) = ( rho0_d(1:n_gas) * c0_d(1:n_gas)**2.d0 -                &
         gamma_d(1:n_gas) * p0_d(1:n_gas) ) / gamma_d(1:n_gas)


    ! ------- READ crystals_parameters NAMELIST ---------------------------------
    READ(input_unit, crystals_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist CRYSTALS_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    t0_c(1:n_cry) = c0_c(1:n_cry) **2.d0 / ( cv_c(1:n_cry) * gamma_c(1:n_cry)   &
         * ( gamma_c(1:n_cry) - 1.d0 ) )

    bar_p_c(1:n_cry) = ( rho0_c(1:n_cry) * c0_c(1:n_cry) ** 2.d0 -              &
         gamma_c(1:n_cry) * p0_c(1:n_cry) ) / gamma_c(1:n_cry)


    IF (.NOT. (crystallization_model .EQ. 'Vitturi2010' ) ) THEN

       WRITE(*,*) ''
       WRITE(*,*) 'Wrong crystallization model chosen.'
       WRITE(*,*) 'Please choose between:'
       WRITE(*,*) ''
       WRITE(*,*) 'Vitturi2010'
       WRITE(*,*) ''

       CALL abort

    END IF

    IF ( (crystallization_model .EQ. 'Vitturi2010' ) .AND.                      & 
         (.NOT.(n_cry .EQ. 1 ))  ) THEN

       WRITE(*,*) ''
       WRITE(*,*) 'Wrong number of crystal components inserted'
       WRITE(*,*) 'Using Vitturi2010, n_cry has to be 1'
       WRITE(*,*) ''

       CALL abort

    END IF

    ! ------- READ melt_parameters NAMELIST -------------------------------------
    READ(input_unit, melt_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist MELT_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    t0_m = c0_m **2.d0 / ( cv_m * gamma_m * ( gamma_m - 1.d0 ) )
    bar_p_m = ( rho0_m * c0_m ** 2.d0 - gamma_m * p0_m ) / gamma_m

    ! ------- READ viscosity_parameters NAMELIST --------------------------------
    READ(input_unit, viscosity_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist VISCOSITY_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    check_model = .false.

    DO i=1,n_visc_melt_models

       IF ( trim(visc_melt_model) .EQ. trim(available_visc_melt_models(i)) ) THEN
          
          check_model = .true.

       END IF

    END DO

    IF ( .NOT.check_model ) THEN

       WRITE(*,*) 'Wrong drag_funct_model chosen.'
       WRITE(*,*) 'Please choose between:'
       
       DO i=1,n_visc_melt_models
          
          WRITE(*,*) available_visc_melt_models(i)

       END DO

       stop

    END IF


    IF ( visc_melt_model .EQ. 'Giordano_et_al2008' ) THEN
       
       tend1 = .false.
       
       WRITE(*,*) 'search melt composition'
       
       melt_comp_search: DO
          
          READ(input_unit,*, end = 300 ) card
          
          IF( trim(card) == 'MELTCOMPOSITION' ) THEN
             
             EXIT melt_comp_search

          END IF
          
       END DO melt_comp_search
       
       READ(input_unit,*) 
       READ(input_unit,*) wt_init(1:12)
       
       IF ( verbose_level .GE. 1 ) WRITE(*,*) wt_init(1:12)

       GOTO 310
300    tend1 = .true.
310    CONTINUE

       rewind(input_unit)

    END IF

    check_model = .false.

    DO i=1,n_theta_models

       IF ( trim(theta_model) .EQ. trim(available_theta_models(i)) ) THEN
          
          check_model = .true.

       END IF

    END DO

    IF ( .NOT.check_model ) THEN

       WRITE(*,*) 'Wrong theta_model chosen.'
       WRITE(*,*) 'Please choose between:'
       
       DO i=1,n_theta_models
          
          WRITE(*,*) available_theta_models(i)

       END DO

       stop

    END IF

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

        WRITE(*,*) 'WARNING: bubbles_model not used.' 
        WRITE(*,*) 'Effect of bubbles is included by using Vona 2013 equations'  

    END IF

    check_model = .false.

    DO i=1,n_bubble_models

       IF ( trim(bubbles_model) .EQ. trim(available_bubble_models(i)) ) THEN
          
          check_model = .true.

       END IF

    END DO

    IF ( .NOT.check_model ) THEN

       WRITE(*,*) 'Wrong bubbles_model chosen.'
       WRITE(*,*) 'Please choose between:'
       
       DO i=1,n_bubble_models
          
          WRITE(*,*) available_bubble_models(i)

       END DO

       stop

    END IF


    ! ------- READ temperature_parameters NAMELIST ------------------------------
    READ(input_unit, temperature_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist TEMPERATURE_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    IF ( isothermal ) THEN

       t_in = fixed_temp

    END IF

    ! ------- READ fragmentation_parameters NAMELIST ----------------------------
    READ(input_unit, fragmentation_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist FRAGMENTATION_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    IF (fragmentation_model .NE. 1 ) THEN

       WRITE(*,*) ''
       WRITE(*,*) 'Wrong fragmentation model chosen.'
       WRITE(*,*) 'Please choose between:'
       WRITE(*,*) ''
       WRITE(*,*) '1 (Volume_Fraction)'
       WRITE(*,*) ''

       CALL abort

    END IF

    ! ------- READ external_water_parameters NAMELIST ---------------------------
    READ(input_unit, external_water_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist EXTERNAL_WATER_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    IF (ext_water) THEN

             IF ( (.NOT. (aquifer_type .EQ. 'confined' ) ) .AND.                & 
            (.NOT. (aquifer_type .EQ. 'unconfined' ) ) ) THEN

          WRITE(*,*) ''
          WRITE(*,*) 'Wrong aquifer type chosen.'
          WRITE(*,*) 'Please choose between:'
          WRITE(*,*) ''
          WRITE(*,*) 'confined'
          WRITE(*,*) 'unconfined'
          WRITE(*,*) ''

          CALL abort

       END IF
    
    END IF

    ! ------- READ source_parameters NAMELIST -----------------------------------
    READ(input_unit, source_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist SOURCE_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    ! ------- READ country_rock_parameters NAMELIST -----------------------------
    IF ( lateral_degassing_flag ) THEN
       
       READ(input_unit, country_rock_parameters , iostat = ios )

       IF ( ios .NE. 0 ) THEN
          
          WRITE(*,*) 'IOSTAT=',ios
          WRITE(*,*) 'ERROR: problem with namelist COUNTRY_ROCK_PARAMETERS'
          WRITE(*,*) 'Please check the input file'
          stop
          
       ELSE
          
          rewind(input_unit)
          
       END IF
              
    END IF

    k_cr = 10 ** log10_k_cr

    ! ------- READ relaxation_parameters NAMELIST -------------------------------
    ALLOCATE( log10_tau_d(n_gas) )
    ALLOCATE( log10_tau_c(n_cry) )

    READ(input_unit, relaxation_parameters , iostat = ios )

    IF ( ios .NE. 0 ) THEN
       
       WRITE(*,*) 'IOSTAT=',ios
       WRITE(*,*) 'ERROR: problem with namelist RELAXATION_PARAMETERS'
       WRITE(*,*) 'Please check the input file'
       stop
       
    ELSE
       
       rewind(input_unit)
       
    END IF

    
    drag_funct_coeff = 10.d0 ** log10_drag_funct_coeff
    tau_p_coeff = 10.d0 ** log10_tau_p_coeff
    tau_d(1:n_gas) = 10.d0 ** log10_tau_d(1:n_gas)
    tau_c(1:n_cry) = 10.d0 ** log10_tau_c(1:n_cry)


    check_model = .false.

    DO i=1,n_drag_models

       IF ( trim(drag_funct_model) .EQ. trim(available_drag_models(i)) ) THEN
          
          check_model = .true.

       END IF

    END DO

    IF ( .NOT.check_model ) THEN

       WRITE(*,*) 'Wrong drag_funct_model chosen.'
       WRITE(*,*) 'Please choose between:'
       
       DO i=1,n_drag_models
          
          WRITE(*,*) available_drag_models(i)

       END DO

       stop

    END IF
       
    IF ( ( drag_funct_model .EQ. 'eval' ) .OR.                                  &
         ( drag_funct_model .EQ. 'Klug_and_Cashman' ) .OR.                      &
         ( drag_funct_model .EQ. 'drag' ) ) THEN
       
       READ(input_unit, bubbles_parameters , iostat = ios )

       IF ( ios .NE. 0 ) THEN
          
          WRITE(*,*) 'IOSTAT=',ios
          WRITE(*,*) 'ERROR: problem with namelist BUBBLES_PARAMETERS'
          WRITE(*,*) 'Please check the input file'
          stop
          
       ELSE
          
          rewind(input_unit)
          
       END IF
              
    ELSEIF ( ( drag_funct_model .EQ. 'darcy' ) .OR.                             & 
         ( drag_funct_model .EQ. 'forchheimer') .OR.                            & 
         ( drag_funct_model .EQ. 'forchheimer_wt')) THEN
       
       READ(input_unit, forchheimer_parameters , iostat = ios )

       IF ( ios .NE. 0 ) THEN
          
          WRITE(*,*) 'IOSTAT=',ios
          WRITE(*,*) 'ERROR: problem with namelist FORCHHEIMER_PARAMETERS'
          WRITE(*,*) 'Please check the input file'
          stop
          
       ELSE
       
          rewind(input_unit)
          
       END IF
              
       bubble_number_density = 10.0d0 ** log10_bubble_number_density
       
    ELSEIF ( ( drag_funct_model .EQ. 'forchheimer_mod' ) .OR.                   &
         ( drag_funct_model .EQ. 'forchheimer_mod2' ) .OR.                      &
         ( drag_funct_model .EQ. 'forchheimer_mod3' ) ) THEN

       READ(input_unit, permeability_parameters , iostat = ios )

       IF ( ios .NE. 0 ) THEN
          
          WRITE(*,*) 'IOSTAT=',ios
          WRITE(*,*) 'ERROR: problem with namelist PERMEABILITY_PARAMETERS'
          WRITE(*,*) 'Please check the input file'
          stop
          
       ELSE
          
          rewind(input_unit)
          
       END IF       
       
    END IF

    CLOSE( input_unit )

    bak_name = trim(run_name)//'.bak'

    OPEN(backup_unit,file=bak_name,status='unknown')

    WRITE(backup_unit, run_parameters )

    WRITE(backup_unit, geometry_parameters )

    WRITE(backup_unit, phases_parameters )

    WRITE(backup_unit, steady_boundary_conditions )

    WRITE(backup_unit, exsolved_gas_parameters )

    WRITE(backup_unit, dissolved_gas_parameters )

    WRITE(backup_unit, crystals_parameters )

    WRITE(backup_unit, melt_parameters )

    WRITE(backup_unit, viscosity_parameters )

    WRITE(backup_unit, temperature_parameters )

    WRITE(backup_unit, fragmentation_parameters )

    WRITE(backup_unit, external_water_parameters )

    WRITE(backup_unit, source_parameters )

    IF ( lateral_degassing_flag ) THEN
       
       WRITE(backup_unit, country_rock_parameters )
       
    END IF

    WRITE(backup_unit, relaxation_parameters )


    IF ( ( drag_funct_model .EQ. 'eval' ) .OR.                                  &
         ( drag_funct_model .EQ. 'Klug_and_Cashman' ) .OR.                      &
         ( drag_funct_model .EQ. 'drag' ) ) THEN
       
       WRITE(backup_unit, bubbles_parameters )
              
    ELSEIF ( ( drag_funct_model .EQ. 'darcy' ) .OR.                             & 
         ( drag_funct_model .EQ. 'forchheimer') .OR.                            & 
         ( drag_funct_model .EQ. 'forchheimer_wt')) THEN
       
       WRITE(backup_unit, forchheimer_parameters )
       
    ELSEIF ( ( drag_funct_model .EQ. 'forchheimer_mod' ) .OR.                   &
         ( drag_funct_model .EQ. 'forchheimer_mod2' ) .OR.                      &
         ( drag_funct_model .EQ. 'forchheimer_mod3' ) ) THEN

       WRITE(backup_unit, permeability_parameters  )
       
    END IF

    
    IF ( visc_melt_model .EQ. 'Giordano_et_al2008' ) THEN
       
       WRITE(backup_unit,*) '''MELTCOMPOSITION'''
       WRITE(backup_unit,*) 'SiO2  TiO2 Al2O3  FeO   MnO   MgO   CaO   Na2O  K2O   P2O5  H2O   F2O-1'     
       WRITE(backup_unit,107) wt_init(1:12)
       
107    FORMAT(12(f5.2,1x))
       
    END IF
    

    CLOSE(backup_unit)

  END SUBROUTINE read_param

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

  SUBROUTINE output_steady(zeta,qp,radius)

    ! External variables
    USE parameters, ONLY : n_vars , n_cry , n_gas

    USE parameters, ONLY : idx_p1 , idx_p2 , idx_u1 , idx_u2 , idx_t ,          &
         idx_alfa_first , idx_alfa_last , idx_beta_first , idx_beta_last

    USE geometry, ONLY : zeta_exit

    ! External subroutine
    USE constitutive, ONLY : sound_speeds
    USE constitutive, ONLY : frag_eff
    USE equations, ONLY : eval_qp2

    IMPLICIT NONE

    REAL*8, INTENT(IN) :: zeta
    REAL*8 :: qp(n_vars)
    REAL*8, INTENT(IN) :: radius

    REAL*8 :: qp2(1+n_cry+n_gas+n_gas+4)
    REAL*8 :: mass_flow_rate
    REAL*8 :: pi

    REAL*8 :: C_mix, mach, mix_velocity

    REAL*8 :: rho1, rho2(n_gas), rho_mix

    REAL*8 :: gasTotVolFrac

    INTEGER :: i,j

    CALL eval_qp2( qp, qp2 )

    IF ( zeta .EQ. z0 ) THEN

       steady_p_file = trim(run_name)//'_p.std'

       OPEN(steady_p_output_unit,file=steady_p_file,status='UNKNOWN')

       OPEN(steady_q_output_unit,file='test.out',status='UNKNOWN')
       !WRITE(steady_p_output_unit,*) radius_fixed
       !WRITE(steady_p_output_unit,*) radius_min
       !WRITE(steady_p_output_unit,*) radius_max
       !WRITE(steady_p_output_unit,*) radius_z
       !WRITE(steady_p_output_unit,*) radius_z_sig
       !WRITE(steady_p_output_unit,*) radius_model
       !WRITE(steady_p_output_unit,*) comp_cells

    END IF



    DO i = 1,n_vars

       IF ( abs(qp(i)) .LT. 1d-20 ) THEN
          qp(i) = 0.d0
       END IF

    END DO

    DO i = 1,1+n_cry+n_gas+n_gas+4

       IF ( abs(qp2(i)) .LT. 1d-20 ) THEN
          qp2(i) = 0.d0
       END IF

    END DO

    WRITE(steady_p_output_unit, *) (zeta)
    WRITE(steady_p_output_unit, 1006) (qp(i), i=1,n_vars)
    WRITE(steady_p_output_unit, 1006) (qp2(i), i=1,1+n_cry+n_gas+n_gas+4)
    WRITE(steady_p_output_unit, *) radius

1006 FORMAT(4e20.12) 

    WRITE(steady_q_output_unit,1007,advance="no") zeta

    DO i = 1,n_vars

       WRITE(steady_q_output_unit,1007,advance="no") qp(i)

    END DO

    DO i = 1,1+n_cry+n_gas+n_gas+4

       WRITE(steady_q_output_unit,1007,advance="no") qp2(i)

    END DO
       
    WRITE(steady_q_output_unit,1007) radius

1007 FORMAT(1x,e15.8)

    IF ( zeta .EQ. zeta_exit ) THEN

       rho1 = qp2(1)

       gastotvolfrac = 0.0

       DO j=1,n_gas

          rho2(j) = qp2(1+j)
          gastotvolfrac = gastotvolfrac + qp(j)

       END DO

       rho_mix = (1.0d0 - gastotvolfrac) * rho1

       DO j=1,n_gas

          rho_mix = rho_mix + qp(j) * rho2(j)

       END DO

       mix_velocity = (1.0d0 - gastotvolfrac) * rho1 / rho_mix * qp(idx_u1)

       DO j=1,n_gas

          mix_velocity = mix_velocity + qp(j) * rho2(j) / rho_mix * qp(idx_u2)

       END DO

       CLOSE(steady_p_output_unit)
       CLOSE(steady_q_output_unit)

       OPEN( dakota_unit,file='MAMMA.out',status='UNKNOWN')

       pi = 4.d0 * atan(1.d0)

       WRITE(dakota_unit,*) 'Total Gas volume fraction', gastotvolfrac
       WRITE(dakota_unit,*) 'Pressure 1',qp(idx_p1)
       WRITE(dakota_unit,*) 'Pressure 2',qp(idx_p2)
       WRITE(dakota_unit,*) 'Liquid/particles velocity',qp(idx_u1)
       WRITE(dakota_unit,*) 'Gas velocity',qp(idx_u2)
       WRITE(dakota_unit,*) 'Exit Temperature',qp(idx_t)
       WRITE(dakota_unit,*) 'Crystals volume fraction',sum(qp(idx_beta_first:   &
            idx_beta_last))

       CALL sound_speeds(c_mix,mach) 
       CALL update_radius(zeta)

       WRITE(dakota_unit,*) 'Mach number',mach
       WRITE(dakota_unit,*) 'Exit mixture velocity', mix_velocity
       WRITE(dakota_unit,*) 'Mass flow rate', pi * radius * radius &
            * mix_velocity * rho_mix
       WRITE(dakota_unit,*) 'Fragmentation', frag_eff

       CLOSE( dakota_unit )

       OPEN( dakota_unit2,file='dakota.out',status='UNKNOWN')

       pi = 4.d0 * atan(1.d0)

       WRITE(dakota_unit2,*) gastotvolfrac
       WRITE(dakota_unit2,*) qp(idx_p1)
       WRITE(dakota_unit2,*) qp(idx_p2)
       WRITE(dakota_unit2,*) qp(idx_u1)
       WRITE(dakota_unit2,*) qp(idx_u2)
       WRITE(dakota_unit2,*) qp(idx_t)
       WRITE(dakota_unit2,*) mach
       WRITE(dakota_unit2,*) mix_velocity
       WRITE(dakota_unit2,*) pi * radius * radius * mix_velocity * rho_mix
       WRITE(dakota_unit2,*) frag_eff

       CLOSE( dakota_unit2 )

       exit_file = trim(run_name)//'_exit.std'

       OPEN( exit_unit,file=exit_file,status='UNKNOWN')

       WRITE(exit_unit,*) 'Mass_flow_rate',mass_flow_rate
       WRITE(exit_unit,*) 'Gas volume fraction',1.d0-sum(qp(idx_alfa_first:     &
            idx_alfa_last))
       WRITE(exit_unit,*) 'Pressure 1',qp(idx_p1)
       WRITE(exit_unit,*) 'Pressure 2',qp(idx_p2)
       WRITE(exit_unit,*) 'Crystals volume fraction',sum(qp(idx_beta_first:     &
            idx_beta_last))
       WRITE(exit_unit,*) 'Liquid/particles velocity',qp(idx_u1)
       WRITE(exit_unit,*) 'Gas velocity',qp(idx_u2)
       WRITE(exit_unit,*) 'Mach number',mach

       CLOSE( exit_unit )

    END IF

  END SUBROUTINE output_steady

  !------------------------------------------------------------------------------
  !
  !------------------------------------------------------------------------------

  CHARACTER*4 FUNCTION lettera(k)
    IMPLICIT NONE
    CHARACTER ones,tens,hund,thou
    !
    INTEGER :: k
    !
    INTEGER :: iten, ione, ihund, ithou
    !
    ithou=int(k/1000)
    ihund=int((k-(ithou*1000))/100)
    iten=int((k-(ithou*1000)-(ihund*100))/10)
    ione=k-ithou*1000-ihund*100-iten*10
    ones=char(ione+48)
    tens=char(iten+48)
    hund=char(ihund+48)
    thou=char(ithou+48)
    lettera=thou//hund//tens//ones
    !
    RETURN
  END FUNCTION lettera

END MODULE inpout