inpout.f90

Go to the documentation of this file.

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

MODULE inpout

  ! -- Variables for the namelist RUN_PARAMETERS
  USE parameters, ONLY : solver_scheme , max_dt , t_start , t_end ,              &
       dt_output , cfl, limiter , theta , reconstr_coeff , reconstr_variables ,  &
       interfaces_relaxation , n_rk , batimetry_function_flag , riemann_flag

  USE solver, ONLY : verbose_level

  ! -- Variables for the namelist NEWRUN_PARAMETERS
  USE geometry, ONLY : x0 , xn , comp_cells
  USE geometry, ONLY : batimetry_profile , n_batimetry_profile
  USE init, ONLY : riemann_interface

  ! -- Variables for the namelist LEFT_STATE
  USE init, ONLY : hb_l , u_l

  ! -- Variables for the namelist RIGHT_STATE
  USE init, ONLY : hb_r , u_r

  ! -- Variables for the namelists LEFT/RIGHT_BOUNDARY_CONDITIONS
  USE parameters, ONLY : bc

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

  IMPLICIT NONE

  CHARACTER(LEN=40) :: run_name
  CHARACTER(LEN=40) :: bak_name
  CHARACTER(LEN=40) :: input_file
  CHARACTER(LEN=40) :: output_file
  CHARACTER(LEN=40) :: restart_file
  CHARACTER(LEN=40) :: batimetry_file
  CHARACTER(LEN=40) :: dakota_file

  INTEGER, PARAMETER :: input_unit = 7
  INTEGER, PARAMETER :: backup_unit = 8
  INTEGER, PARAMETER :: output_unit = 9
  INTEGER, PARAMETER :: restart_unit = 10
  INTEGER, PARAMETER :: dakota_unit = 11


  INTEGER :: output_idx 

  LOGICAL :: restart

  ! -- Variables for the namelists LEFT_BOUNDARY_CONDITIONS
  TYPE(bc) :: hB_bcL , u_bcL

  ! -- Variables for the namelists RIGHT_BOUNDARY_CONDITIONS
  TYPE(bc) :: hB_bcR , u_bcR


  namelist / run_parameters / run_name , restart , batimetry_function_flag ,    &
       riemann_flag , max_dt , t_start , t_end , dt_output , solver_scheme ,    &
       cfl , limiter , theta , reconstr_coeff , reconstr_variables , n_rk ,     &
       verbose_level

  namelist / restart_parameters / restart_file

  namelist / newrun_parameters / x0 , xn , comp_cells , riemann_interface

  namelist / left_state / hb_l , u_l

  namelist / right_state / hb_r , u_r

  namelist / left_boundary_conditions / hb_bcl , u_bcl

  namelist / right_boundary_conditions / hb_bcr , u_bcr

  namelist / source_parameters / grav

CONTAINS

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

  SUBROUTINE init_param

    USE parameters , ONLY : n_vars

    IMPLICIT none

    LOGICAL :: lexist

    INTEGER :: i

    !-- Inizialization of the Variables for the namelist RUN_PARAMETERS
    run_name = 'default'
    restart = .false.
    batimetry_function_flag=.false.
    riemann_flag=.true.
    max_dt = 1.d-3
    t_start = 0.0
    t_end = 5.0d-4
    dt_output = 1.d-4
    solver_scheme = 'LxF'
    n_rk = 2
    verbose_level = 0


    cfl = 0.5
    limiter(1:n_vars) = 0
    theta=1.0
    reconstr_coeff = 1.0
    reconstr_variables = 0
    !-- Inizialization of the Variables for the namelist restart parameters
    restart_file = ''

    !-- Inizialization of the Variables for the namelist newrun_parameters
    x0 = 0.d0
    xn = 1.d0
    comp_cells = 500
    riemann_interface = 0.5d0

    !-- Inizialization of the Variables for the namelist left_state
    hb_l = 1.d0
    u_l = 0.d0

    !-- Inizialization of the Variables for the namelist right_state
    hb_r = 0.5d0
    u_r = 0.d0

    !-- Inizialization of the Variables for the namelist left boundary conditions

    hb_bcl%flag = 1 
    hb_bcl%value = 0.d0 

    u_bcl%flag = 1 
    u_bcl%value = 0.d0 

    !-- Inizialization of the Variables for the namelist right boundary conditions

    hb_bcr%flag = 1 
    hb_bcr%value = 0.d0 

    u_bcr%flag = 1 
    u_bcr%value = 0.d0 

    !-- Inizialization of the Variables for the namelist source_parameters
    grav = -9.81d0

    input_file = 'shallow_water.inp'

    INQUIRE (file=input_file,exist=lexist)

    IF (lexist .EQV. .false.) THEN

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

       WRITE(input_unit, run_parameters )
       WRITE(input_unit, restart_parameters )
       WRITE(input_unit, newrun_parameters )
       WRITE(input_unit, left_state )
       WRITE(input_unit, right_state )
       WRITE(input_unit, left_boundary_conditions )
       WRITE(input_unit, right_boundary_conditions )
       WRITE(input_unit, source_parameters )

       n_batimetry_profile = 2

       ALLOCATE( batimetry_profile(2,n_batimetry_profile) )

       batimetry_profile(1,1) = 0.d0
       batimetry_profile(2,1) = 1.d0

       batimetry_profile(2,2) = 0.d0
       batimetry_profile(2,2) = 10.d0


       WRITE(input_unit,*) '''BATIMETRY_PROFILE'''
       WRITE(input_unit,*) n_batimetry_profile

       DO i = 1, n_batimetry_profile

          WRITE(input_unit,108) batimetry_profile(1:2,i)

108       FORMAT(2(1x,e14.7))


       END DO

       CLOSE(input_unit)

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

    ELSE

    END IF

    ! output file index
    output_idx = 0

  END SUBROUTINE init_param

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

  SUBROUTINE read_param

    USE parameters, ONLY : bcl , bcr

    IMPLICIT none

    REAL*8 :: max_cfl

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

    INTEGER :: i

    CHARACTER(LEN=4) :: idx_string

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


    ! ------- READ run_parameters NAMELIST -----------------------------------
    READ(input_unit, run_parameters )

    IF ( ( solver_scheme .NE. 'LxF' ) .AND. ( solver_scheme .NE. 'KT' ) .AND. &
         ( solver_scheme .NE. 'GFORCE' ) ) THEN

       WRITE(*,*) 'WARNING: no correct solver scheme selected',solver_scheme
       WRITE(*,*) 'Chose between: LxF, GFORCE or KT'
       stop

    END IF

    IF  ( ( solver_scheme.EQ.'LxF' ) .OR. ( solver_scheme.EQ.'GFORCE' ) ) THEN 

       max_cfl = 1.0

    ELSEIF ( solver_scheme .EQ. 'KT' ) THEN

       max_cfl = 0.5

    END IF


    IF ( ( cfl .GT. max_cfl ) .OR. ( cfl .LT. 0.d0 ) ) THEN

       WRITE(*,*) 'WARNING: wrong value of cfl ',cfl
       WRITE(*,*) 'Choose a value between 0.0 and ',max_cfl
       READ(*,*)

    END IF


    WRITE(*,*) 'Limiters',limiter

    IF ( ( maxval(limiter) .GT. 3 ) .OR. ( minval(limiter) .LT. 0 ) ) THEN

       WRITE(*,*) 'WARNING: wrong limiter ',limiter
       WRITE(*,*) 'Choose among: none, minmod,superbee,van_leer'
       stop         

    END IF

    IF ( ( reconstr_coeff .GT. 1.0d0 ) .OR. ( reconstr_coeff .LT. 0.d0 ) ) THEN

       WRITE(*,*) 'WARNING: wrong value of reconstr_coeff ',reconstr_coeff
       WRITE(*,*) 'Change the value between 0.0 and 1.0 in the input file'
       READ(*,*)

    END IF

    IF ( ( reconstr_variables .GT. 2 ) .OR. ( reconstr_variables .LT. 0 ) ) THEN

       WRITE(*,*)'WARNING: wrong value of reconstr_variables ',reconstr_variables
       WRITE(*,*)'Choose the value among: 0, 1 and 2'
       READ(*,*)

    END IF


    IF ( restart ) THEN

       READ(input_unit,restart_parameters)

    ELSE

       READ(input_unit,newrun_parameters)

       IF ( riemann_flag ) THEN
       
          READ(input_unit,left_state)
          READ(input_unit,right_state)

       END IF

    END IF

    READ(input_unit,left_boundary_conditions)

    bcl(1) = hb_bcl 
    bcl(2) = u_bcl 

    READ(input_unit,right_boundary_conditions)

    bcr(1) = hb_bcr 
    bcr(2) = u_bcr 

    READ(input_unit, source_parameters )

    tend1 = .false.

    WRITE(*,*) 'search batimetry_profile'

    batimetry_profile_search: DO

       READ(input_unit,*, END = 200 ) card

       IF( trim(card) == 'BATIMETRY_PROFILE' ) THEN

          EXIT batimetry_profile_search

       END IF

    END DO batimetry_profile_search

    READ(input_unit,*) n_batimetry_profile

    IF ( verbose_level .GE. 1 ) WRITE(*,*) 'n_batimetry_profile' ,            &
         n_batimetry_profile

    ALLOCATE( batimetry_profile(2,n_batimetry_profile) )

    DO i = 1, n_batimetry_profile

       READ(input_unit,*) batimetry_profile(1:2,i)

       IF ( verbose_level .GE. 1 ) WRITE(*,*) i,batimetry_profile(1:2,i)

    END DO

    goto 210
200 tend1 = .true.
210 CONTINUE


    CLOSE( input_unit )

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

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

    WRITE(backup_unit, run_parameters )

    IF ( restart ) THEN

       WRITE(backup_unit,restart_parameters)

    ELSE

       WRITE(backup_unit,newrun_parameters)

       IF ( riemann_flag ) THEN

          WRITE(backup_unit,left_state)
          WRITE(backup_unit,right_state)

       END IF

    END IF

    WRITE(backup_unit,left_boundary_conditions)
    WRITE(backup_unit,right_boundary_conditions)

    WRITE(backup_unit, source_parameters )

    WRITE(backup_unit,*) '''BATIMETRY_PROFILE'''
    WRITE(backup_unit,*) n_batimetry_profile

    DO i = 1, n_batimetry_profile

       WRITE(backup_unit,107) batimetry_profile(1:2,i)

107    FORMAT(2(1x,e14.7))


    END DO


    CLOSE(backup_unit)


  END SUBROUTINE read_param


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

  SUBROUTINE read_solution

    USE geometry, ONLY : init_grid
    USE solver, ONLY : allocate_solver_variables

    USE geometry , ONLY : comp_cells , x0 , dx
    USE parameters, ONLY : n_vars
    USE parameters, ONLY : t_start
    USE solver, ONLY : q

    IMPLICIT none

    INTEGER :: j
    INTEGER :: i

    CHARACTER(LEN=30) :: string

    OPEN(restart_unit,file=restart_file,status='old')

    READ(restart_unit,1001) x0,dx,comp_cells,t_start
1001 FORMAT(e18.8,'    x0', /, e18.8,'    dx', /, i18,'    cells', /, e18.8,'    t', /)


    CALL init_grid

    CALL allocate_solver_variables


    DO i = 1,n_vars

       DO j = 1,comp_cells

          ! Exponents with more than 2 digits cause problems reading
          ! into matlab... reset tiny values to zero:
          IF ( dabs(q(i,j)) .LT. 1d-99) q(i,j) = 0.d0

       END DO

       READ(restart_unit,1003) (q(i,j), j=1,n_vars)
1003   FORMAT(4e20.12)

    END DO

    j = scan(restart_file, '.' , .true. )

    string = trim(restart_file(j+2:j+5))

    READ( string,* ) output_idx

  END SUBROUTINE read_solution

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

  SUBROUTINE output_solution(t)

    ! external procedures
    USE constitutive, ONLY : qc_to_qp

    ! external variables
    USE constitutive, ONLY : h , u
    USE geometry , ONLY : comp_cells , x0 , dx , x_comp , b_cent
    USE parameters, ONLY : n_vars
    USE parameters, ONLY : t_output , dt_output
    USE solver, ONLY : q

    IMPLICIT none

    REAL*8, INTENT(IN) :: t

    CHARACTER(LEN=4) :: idx_string

    REAL*8 :: qp(n_vars)

    INTEGER :: j
    INTEGER :: i

    output_idx = output_idx + 1

    idx_string = lettera(output_idx-1)

    output_file = trim(run_name)//'.q'//idx_string

    WRITE(*,*) 'WRITING ',output_file

    OPEN(output_unit,file=output_file,status='unknown',form='formatted')

    WRITE(output_unit,1002) x0,dx,comp_cells,t

    DO j = 1,comp_cells

       DO i = 1,n_vars

          ! Exponents with more than 2 digits cause problems reading
          ! into matlab... reset tiny values to zero:
          IF ( dabs(q(i,j)) .LT. 1d-99) q(i,j) = 0.d0

       END DO

       WRITE(output_unit,*) (q(i,j), i=1,n_vars)

    END DO

    WRITE(output_unit,*) ' '
    WRITE(output_unit,*) ' '

    CLOSE(output_unit)

    output_file = trim(run_name)//'.p'//idx_string

    WRITE(*,*) 'WRITING ',output_file

    OPEN(output_unit,file=output_file,status='unknown',form='formatted')

    WRITE(output_unit,1002) x0,dx,comp_cells,t

    DO j = 1,comp_cells

       DO i = 1,n_vars

          ! Exponents with more than 2 digits cause problems reading
          ! into matlab... reset tiny values to zero:
          IF ( dabs(q(i,j)) .LT. 1d-99) q(i,j) = 0.d0

       END DO

       CALL qc_to_qp(q(:,j),b_cent(j),qp(:))

       IF ( REAL(h) .LT. 1d-99) h = 0.d0
       IF ( abs(REAL(u)) .LT. 1d-99) u = 0.d0


       WRITE(output_unit,1005) x_comp(j), REAL(h) , REAL(u) , b_cent(j) ,           &
            REAL(h) + b_cent(j)

    END DO

    WRITE(output_unit,*) ' '
    WRITE(output_unit,*) ' '

    CLOSE(output_unit)

1002 FORMAT(e18.8,'    x0', /, e18.8,'    dx', /, i18,'    cells', /, e18.8,'    t', /)


1005 FORMAT(5e20.12)


    t_output = t + dt_output

  END SUBROUTINE output_solution


  SUBROUTINE close_units

    IMPLICIT NONE

  END SUBROUTINE close_units

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

  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

  SUBROUTINE output_dakota

    IMPLICIT NONE

    dakota_file = 'shallow_water.out'

    OPEN(dakota_unit,file=dakota_file,status='unknown',form='formatted')

    CLOSE(dakota_unit)

  END SUBROUTINE output_dakota


END MODULE inpout