module data
  !
  !  SUBROUTINES
  !
  !  1) get_parameters
  !  2) write_parameters
  !  3) allocate_variables
  !  4) initialize_variables
  !  5) date_time
  !
  ! Last update: 03 Oct 2011
  !
  implicit none
  !
  ! Global variables
  !
  integer :: nx     ! Number of real + ficititious volumes in the csi direction
  integer :: ny     ! Number of real + ficititious volumes in the eta direction
  integer :: nxy    ! Total number of real + fictitious volumes
  integer :: kg     ! Kind of grid
  integer :: coord  ! Coordinate system ( 1=cylindrical, 0 = cartesian)
  integer :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
  integer :: ccTw   ! ccTw = 0 -> adiabatic;  ccTw = 1 -> prescribed temperature
  integer :: ig     ! Number of volume (x direction) of the throttle
  integer :: it1    ! dt = dt1  for it  <= it1
  integer :: it2    ! dt = dt2  for it2 <= it
  integer :: itb1   ! beta = beta1 for   it <= itb1
  integer :: itb2   ! beta = beta2 for itb2 <= it
  integer :: ito    ! First iteraction (ito = 1 at first iteraction and other value for reload)
  integer :: itmax  ! Maximum number of iteractions of the time evolution
  integer :: imax   ! Maximum number of iteractions for the pressure correction
  integer :: nitm_u ! Maximum number of iteractions for the TDMA method for u, v and T
  integer :: nitm_p ! Maximum number of iteractions for the TDMA method for p
  integer :: reload ! Upload backup data and continue computation (0 = no, 1 = yes)
  integer :: wbkp   ! Frequency of saving backup data
  integer :: wlf    ! Frequency of printing in the listing file
  integer :: sem_a  ! 1 = do not open result files, 0 = open
  integer :: sem_g  ! 0 = visualize the plot, 1 = do not visualize
  integer :: w_g    ! Frequency of writing data for graphics
  integer :: w_cam  ! 1 = write the fields, 0 = do not

  integer, parameter :: lid = 100 ! Listing file id
  integer, parameter :: rid = 101 ! Residual file id

  real(8) :: tcpu   ! Total CPU time
  real(8) :: tcpuo  ! Accumulated CPU time before interruption
  real(8) :: tcpu1  ! Starts CPU time
  real(8) :: tcpu2  ! Finishes CPU time

  real(8) :: beta   ! Constant of the UDS/CDS mixing scheme
  real(8) :: beta1  ! Initial constant of the UDS/CDS mixing scheme
  real(8) :: beta2  ! Final constant of the UDS/CDS mixing scheme
  real(8) :: dt     ! Time step
  real(8) :: dt1    ! Initial time step
  real(8) :: dt2    ! Final time step
  real(8) :: a1     ! Initial step for the geometric progression grid

  real(8) :: gamma  ! gamma = Cpo / Cvo in the chamber
  real(8) :: Rg     ! Perfect gas constant
  real(8) :: po     ! Stagnation pressure in the chamber
  real(8) :: T0     ! Stagnation temperature in the chamber
  real(8) :: pr     ! atmospheric pressure at the sea level (Pa)
  real(8) :: go     ! gravitational acceleration at the sea level
  real(8) :: Sg     ! Throttle area
  real(8) :: rcg    ! Throttle radius of curvature

  real(8) :: fmi    ! Mass flow rate at entrance
  real(8) :: fme    ! Mass flow rate at exit
  real(8) :: Fd     ! Thrust

  real(8) :: fm1D   ! Mass flow rate (Isentropic 1D flow)
  real(8) :: Fd1D   ! Dynamic thrust (Isentropic 1D flow)
  real(8) :: Fpv1D  ! Pressure thrust into vacuum (Isentropic 1D flow)

  real(8) :: norm   ! Norm of the linear systems
  real(8) :: norm1  ! Norm of the linear systems of the first iteraction

  character (len = 100) :: sim_id ! Simulation identification
  character (len = 100) :: input_file_parameters ! Input parameters data file
  character (len=20) :: date ! System date
  character (len=20) :: time ! System time

  real(8), allocatable, dimension(:) :: x ! X coordinate of the corner ne of each volume (real+ficititious)
  real(8), allocatable, dimension(:) :: y ! Y coordinate of the corner ne of each volume (real+ficititious)
  real(8), allocatable, dimension(:) :: xp ! Centroid X coordinate of each volume (real+ficititious)
  real(8), allocatable, dimension(:) :: yp ! Centroid Y coordinate of each volume (real+ficititious)

  real(8), allocatable, dimension(:) :: xe     ! x_eta at face east of volume P
  real(8), allocatable, dimension(:) :: ye     ! y_eta at face east of volume P
  real(8), allocatable, dimension(:) :: xen    ! x_eta at face north of volume P
  real(8), allocatable, dimension(:) :: yen    ! y_eta at face north of volume P
  real(8), allocatable, dimension(:) :: xk     ! x_csi at face north of volume P
  real(8), allocatable, dimension(:) :: yk     ! y_csi at face north of volume P
  real(8), allocatable, dimension(:) :: xke    ! x_csi at face east of volume P
  real(8), allocatable, dimension(:) :: yke    ! y_csi at face east of volume P
  real(8), allocatable, dimension(:) :: Jp     ! Jacobian at the center of volume P
  real(8), allocatable, dimension(:) :: Je     ! Jacobian at the center of east face of volume P
  real(8), allocatable, dimension(:) :: Jn     ! Jacobian at the center of north face of volume P
  real(8), allocatable, dimension(:) :: Alphae ! (metric) Alpha at the center of east face of volume P
  real(8), allocatable, dimension(:) :: Gamman ! (metric) Gamma at the center of north face of volume P
  real(8), allocatable, dimension(:) :: Betae  ! (metric) Beta  at the center of east face of volume P
  real(8), allocatable, dimension(:) :: Betan  ! (metric) Beta  at the center of north face of volume P

  real(8), allocatable, dimension(:) :: radius ! Radius of northest corner of volume P
  real(8), allocatable, dimension(:) :: re     ! Radius of the center of east face of volume P
  real(8), allocatable, dimension(:) :: rn     ! Radius of the center of north face of volume P
  real(8), allocatable, dimension(:) :: rp     ! Radius of the center of volume P

  real(8), allocatable, dimension(:) :: p      ! Pressure at center of volumes
  real(8), allocatable, dimension(:) :: pa     ! Pressure at center of volumes at a time step before
  real(8), allocatable, dimension(:) :: T      ! Temperature at center of volumes
  real(8), allocatable, dimension(:) :: ro     ! Specific mass (absolute density) at center of volumes
  real(8), allocatable, dimension(:) :: Ta     ! Temperature at the time step before
  real(8), allocatable, dimension(:) :: roa    ! Absolute density at a time step before at center of vol. P
  real(8), allocatable, dimension(:) :: roe    ! Absolute density at east face
  real(8), allocatable, dimension(:) :: ron    ! Absolute density at north face

  real(8), allocatable, dimension(:) :: u      ! Cartesian velocity u
  real(8), allocatable, dimension(:) :: v      ! Cartesian velocity v
  real(8), allocatable, dimension(:) :: ua     ! Cartesian velocity of a time step before
  real(8), allocatable, dimension(:) :: va     ! Cartesian velocity of a time step before
  real(8), allocatable, dimension(:) :: ue     ! Cartesian velocity u at center of east face
  real(8), allocatable, dimension(:) :: ve     ! Cartesian velocity v at center of east face
  real(8), allocatable, dimension(:) :: un     ! Cartesian velocity u at center of north face
  real(8), allocatable, dimension(:) :: vn     ! Cartesian velocity v at center of north face
  real(8), allocatable, dimension(:) :: uea    ! Cartesian velocity u at center of east face in the last time step
  real(8), allocatable, dimension(:) :: vea    ! Cartesian velocity v at center of east face in the last time step
  real(8), allocatable, dimension(:) :: una    ! Cartesian velocity u at center of north face in the last time step
  real(8), allocatable, dimension(:) :: vna    ! Cartesian velocity v at center of north face in the last time step

  real(8), allocatable, dimension(:) :: Uce    ! Contravariant velocity U at east face
  real(8), allocatable, dimension(:) :: Vcn    ! Contravariant velocity V at north face

  real(8), allocatable, dimension(:) :: Ucem   ! Contravariant velocity U at east face (previous iteraction)
  real(8), allocatable, dimension(:) :: Vcnm   ! Contravariant velocity V at north face (previous iteraction)


  real(8), allocatable, dimension(:) :: vle    ! Laminar viscosity at center of face east
  real(8), allocatable, dimension(:) :: vln    ! Laminar viscosity at center of face north
  real(8), allocatable, dimension(:) :: vlp    ! Laminar viscosity at center of volume P
  real(8), allocatable, dimension(:) :: cp     ! Specific heat at constant pressure at center o volume P
  real(8), allocatable, dimension(:) :: gcp    ! gcp = gamma = Cp/Cv at center of CV P
  real(8), allocatable, dimension(:) :: kp     ! Thermal conductivity at center of volume P
  real(8), allocatable, dimension(:) :: ke     ! Thermal conductivity at center of face east
  real(8), allocatable, dimension(:) :: kn     ! Thermal conductivity at center of face north

  real(8), allocatable, dimension(:) :: cup    ! Term of deferred correction for u
  real(8), allocatable, dimension(:) :: cvp    ! Term of deferred correction for v
  real(8), allocatable, dimension(:) :: sup    ! Viscous term for u
  real(8), allocatable, dimension(:) :: svp    ! Viscous term for v

  real(8), allocatable, dimension(:) :: pl     ! Pressure correction
  real(8), allocatable, dimension(:) :: due    ! SIMPLEC coefficient for ue
  real(8), allocatable, dimension(:) :: dve    ! SIMPLEC coefficient for ve
  real(8), allocatable, dimension(:) :: dun    ! SIMPLEC coefficient for un
  real(8), allocatable, dimension(:) :: dvn    ! SIMPLEC coefficient for vn
  real(8), allocatable, dimension(:) :: de     ! SIMPLEC coefficient for Uce
  real(8), allocatable, dimension(:) :: dn     ! SIMPLEC coefficient for Vcn
  real(8), allocatable, dimension(:) :: bu     ! Source of the linear system for u
  real(8), allocatable, dimension(:) :: b      ! Source of the linear system for v, T
  real(8), allocatable, dimension(:) :: bp     ! Source of the linear system for pl

  real(8), allocatable, dimension(:,:) :: au   ! Coefficients of the linear system for u
  real(8), allocatable, dimension(:,:) :: av   ! Coefficients of the linear system for v
  real(8), allocatable, dimension(:,:) :: at   ! Coefficients of the linear system for T
  real(8), allocatable, dimension(:,:) :: ap   ! Coefficients of the linear system for pl

  real(8), allocatable, dimension(:) :: Twall  ! Wall temperature
  real(8), allocatable, dimension(:) :: Tin    ! Temperature in the entrance
  real(8), allocatable, dimension(:) :: plin   ! Pressure correction in the entrance
  real(8), allocatable, dimension(:) :: uin    ! Cartesian velocity in the entrance
  real(8), allocatable, dimension(:) :: vin    ! Cartesian velocity in the entrance
  real(8), allocatable, dimension(:) :: Mw     ! Mach number in the entrance
  real(8), allocatable, dimension(:) :: pin    ! Pressure in the entrance
  real(8), allocatable, dimension(:) :: pina   ! Pressure in the entrance at a time step before

  real(8), allocatable, dimension(:) :: M1D    ! Mach number distribution. Isentropic flow 1D.
  real(8), allocatable, dimension(:) :: p1D    ! Pressure distribution. Isentropic flow 1D.
  real(8), allocatable, dimension(:) :: T1D    ! Temperature distribution. Isentropic flow 1D.
  real(8), allocatable, dimension(:) :: u1D    ! Velocity distribution. Isentropic flow 1D.

  ! Variables used in the eddy viscosity calculation
  integer :: modtur ! turbulence model option (0=laminar, 1=Baldwin-Lomax)
  real(8), allocatable, dimension(:) :: vtp    ! Eddy viscosity at center of volume P
  real(8), allocatable, dimension(:) :: rona   ! Ro of the previous iteraction level
  real(8), allocatable, dimension(:) :: ro_av  ! Average rho of two iteraction levels
  real(8), allocatable, dimension(:) :: ron_av ! Average rho_n of two iteraction levels
  real(8), allocatable, dimension(:) :: u_av   ! Average u of two iteraction levels
  real(8), allocatable, dimension(:) :: v_av   ! Average v of two iteraction levels

contains

  subroutine get_parameters
    implicit none
    !
    call date_time(date, time)
    !
    open(10,file="./mach2d-input/mach2d-input.dat")
    ! Reading file name from which parameters will be read
    read(10,*) input_file_parameters
    close(10)

    open(10, file = "./mach2d-input/" // input_file_parameters)
    read(10,*) sim_id ! Simulation identification  (up to 100 characters)
    read(10,*) nx     ! Number of real volumes in the csi direction
    read(10,*) ny     ! Number of real volumes in the eta direction
    read(10,*) kg     ! Kind of grid (1 = uniform, 2 = Geometric Progression, 3=power law)
    read(10,*) a1     ! Initial step for the geometric progression grid
    read(10,*) coord  ! Coordinate system ( 1=cylindrical, else cartesian)
    read(10,*) Rg     ! Perfect gas constant
    read(10,*) gamma  ! gamma = Cpo / Cvo in the chamber (Specific heat ratio)
    read(10,*) po     ! Stagnation pressure in the chamber
    read(10,*) T0     ! Stagnation temperature in the chamber
    read(10,*) pr     ! atmospheric pressure at the sea level (Pa)
    read(10,*) go     ! gravitational acceleration at the sea level
    read(10,*) beta1  ! Initial constant of the UDS/CDS mixing scheme
    read(10,*) beta2  ! Final constant of the UDS/CDS mixing scheme
    read(10,*) itb1   ! beta = beta1 for   it <= itb1
    read(10,*) itb2   ! beta = beta2 for itb2 <= it
    read(10,*) modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
    read(10,*) modtur ! turbulence model option (0=laminar, 1=Baldwin-Lomax)
    read(10,*) ccTw   ! ccTw = 0 -> adiabatic;  ccTw = 1 -> prescribed temperature
    read(10,*) it1    ! dt = dt1  for it  <= it1
    read(10,*) it2    ! dt = dt2  for it2 <= it
    read(10,*) dt1    ! Initial time step
    read(10,*) dt2    ! Final time step
    read(10,*) itmax  ! Maximum number of iteractions of the time evolution
    read(10,*) imax   ! Maximum number of iteractions for the pressure correction
    read(10,*) nitm_u ! Maximum number of iteractions for the TDMA method for u, v and T
    read(10,*) nitm_p ! Maximum number of iteractions for the TDMA method for p
    read(10,*) reload ! Upload backup data and continue computation (0 = no, 1 = yes)
    read(10,*) wbkp   ! Frequency of saving backup data
    read(10,*) wlf    ! Frequency of printing in the listing file
    read(10,*) sem_a  ! 1 = do not open result files, 0 = open
    read(10,*) sem_g  ! 0 = visualize the plot, 1 = do not visualize
    read(10,*) w_g    ! Frequency of writing data for graphics
    read(10,*) w_cam  ! 1 = write the fields, 0 = do not
    close(10)
    !
    nx = nx + 2 ! Real + Fictitious
    ny = ny + 2 ! Real + Fictitious
    nxy = nx * ny

    beta = beta1

  end subroutine get_parameters


  subroutine write_parameters(fid)
    implicit none
    integer, intent(in) :: fid

    write(fid,*)
    write(fid,*) "Date: ", date
    write(fid,*) "Time: ", time
    write(fid,*)
    write(fid,*) "          PARAMETERS         "
    write(fid,*)
    write(fid,"(A14,' ....: ',A)") trim(adjustl(sim_id)) , " Simulation identification  (up to 100 characters)"
    write(fid,"(I14,' ....: ',A)")    nx     , " nx     - Number of real+ficititious volumes in the csi direction"
    write(fid,"(I14,' ....: ',A)")    ny     , " ny     - Number of real+ficititious volumes in the eta direction"
    write(fid,"(I14,' ....: ',A)")    kg     , " kg     - Kind of grid (1 = uniform, 2 = Geometric Progression, 3=power law)"
    write(fid,"(ES14.7,' ....: ',A)") a1     , " a1     - Initial step for the geometric progression grid (m)"
    write(fid,"(I14,' ....: ',A)")    coord  , " coord  - Coordinate system ( 1=cylindrical, else cartesian)"
    write(fid,"(ES14.7,' ....: ',A)") Rg     , " Rg     - Perfect gas constant (J/(kg.K))"
    write(fid,"(ES14.7,' ....: ',A)") gamma  , " gamma  -  Cpo / Cvo in the chamber (Specific heat ratio)"
    write(fid,"(ES14.7,' ....: ',A)") po     , " po     - Stagnation pressure in the chamber (Pa)"
    write(fid,"(ES14.7,' ....: ',A)") T0     , " T0     - Stagnation temperature in the chamber (K)"
    write(fid,"(ES14.7,' ....: ',A)") pr     , " pr     - atmospheric pressure at the sea level (Pa)"
    write(fid,"(ES14.7,' ....: ',A)") go     , " go     - gravitational acceleration at the sea level (m/s2)"
    write(fid,"(ES14.7,' ....: ',A)") beta   , " beta   - Constant of the UDS/CDS mixing scheme"
    write(fid,"(ES14.7,' ....: ',A)") beta1  , " beta1  - Initial constant of the UDS/CDS mixing scheme"
    write(fid,"(ES14.7,' ....: ',A)") beta2  , " beta2  - Final constant of the UDS/CDS mixing scheme"
    write(fid,"(I14,' ....: ',A)")    itb1   , " itb1   - beta = beta1 for   it <= itb1"
    write(fid,"(I14,' ....: ',A)")    itb2   , " itb2   - beta = beta2 for itb2 <= it"
    write(fid,"(I14,' ....: ',A)")    modvis , " modvis - Viscosity model ( 0 = Euler,  1 = Navier-Stokes)"
    write(fid,"(I14,' ....: ',A)")    modtur , " modtur - turbulence model option (0=laminar, 1=Baldwin-Lomax)"
    write(fid,"(I14,' ....: ',A)")    ccTw   , " ccTw   - Temp. bc (0 = adiabatic,  1 = prescribed temperature)"
    write(fid,"(I14,' ....: ',A)")    it1    , " dt = dt1  for it  <= it1"
    write(fid,"(I14,' ....: ',A)")    it2    , " dt = dt2  for it2 <= it"
    write(fid,"(ES14.7,' ....: ',A)") dt1    , " Initial time step (s)"
    write(fid,"(ES14.7,' ....: ',A)") dt2    , " Final time step (s)"
    write(fid,"(I14,' ....: ',A)")    itmax  , " itmax  - Maximum number of iteractions of the time evolution"
    write(fid,"(I14,' ....: ',A)")    imax   , " imax   - Maximum number of iteractions for the pressure correction"
    write(fid,"(I14,' ....: ',A)")    nitm_u , " nitm_u - Maximum number of iteractions for the TDMA method for u, v and T"
    write(fid,"(I14,' ....: ',A)")    nitm_p , " nitm_p - Maximum number of iteractions for the TDMA method for p"
    write(fid,"(I14,' ....: ',A)")    reload , " reload - Upload backup data and continue computation (0 = no, 1 = yes)"
    write(fid,"(I14,' ....: ',A)")    wbkp   , " wbkp   - Frequency of saving backup data"
    write(fid,"(I14,' ....: ',A)")    wlf    , " wlf    - Frequency of printing in the listing file"
    write(fid,"(I14,' ....: ',A)")    sem_a  , " sem_a  - 1 = do not open result files, 0 = open"
    write(fid,"(I14,' ....: ',A)")    sem_g  , " sem_g  - 0 = visualize the plot, 1 = do not visualize"
    write(fid,"(I14,' ....: ',A)")    w_g    , " w_g    - Frequency of writing data for graphics"
    write(fid,"(I14,' ....: ',A)")    w_cam  , " w_cam  - 1 = write the fields, 0 = do not"
    write(fid,*)

  end subroutine write_parameters


  subroutine allocate_variables

    implicit none

    allocate( Tin(ny), pin(ny), pina(ny), plin(ny), uin(ny), vin(ny), Mw(ny) )

    allocate( Twall(nx), M1D(nx), p1D(nx), T1D(nx), u1D(nx) )

    allocate( x(nxy),       y(nxy),       xp(nxy),      yp(nxy),        xe(nxy),        ye(nxy)     &
         ,    xen(nxy),     yen(nxy),     xk(nxy),      yk(nxy),        xke(nxy),       yke(nxy)    &
         ,    Jp(nxy),      Je (nxy),     Jn(nxy),      alphae(nxy),    gamman(nxy),    betae(nxy)  &
         ,    betan(nxy),   radius(nxy),  re(nxy),      rn(nxy),        rp(nxy),        p(nxy)      &
         ,    T(nxy),       ro(nxy),      roe(nxy),     ron(nxy),       Uce(nxy),       Vcn(nxy)    &
         ,    roa(nxy),     vle(nxy),     vln(nxy),     u(nxy),         v(nxy),         ua(nxy)     &
         ,    va(nxy),      cup(nxy),     sup(nxy),     vlp(nxy),       cvp(nxy),       svp(nxy)    &
         ,    cp(nxy),      ke(nxy),      kn(nxy),      Ta(nxy),        pl(nxy),        b(nxy)      &
         ,    bp(nxy),      due(nxy),     dve(nxy),     dun(nxy),       dvn(nxy),       de(nxy)     &
         ,    dn(nxy),      ue(nxy),      ve(nxy),      un(nxy),        vn(nxy),        uea(nxy)    &
         ,    vea(nxy),     una(nxy),     vna(nxy),     gcp(nxy),       kp(nxy),        pa(nxy)     &
         ,    bu(nxy),      Ucem(nxy),    Vcnm(nxy),    vtp(nxy),       ro_av(nxy),     ron_av(nxy) &
         ,    u_av(nxy),    v_av(nxy),    rona(nxy) &
         )

    allocate( au(nxy,9), av(nxy,9), at(nxy,9), ap(nxy,9) )

  end subroutine allocate_variables


  subroutine initialize_variables
    implicit none

    ! Auxiliary variable
    logical :: lexist

    call get_parameters

    ! Openning listing file

    inquire(file = "./mach2d-output/mach2d-" // trim(adjustl(sim_id)) // ".lst", exist = lexist)

    if ( lexist .and. reload == 1 ) then
        open(lid, file = "./mach2d-output/mach2d-" // trim(adjustl(sim_id)) // ".lst", POSITION='APPEND')

        write(lid,*)
        write(lid,*) " ==================================================================================== "
        write(lid,*)
        write(lid,*) "           RELOADED. CONTINUING THE CALCULATIONS..."
        write(lid,*)

    else
        open(lid, file = "./mach2d-output/mach2d-" // trim(adjustl(sim_id)) // ".lst")
    end if

    write(lid,*)
    write(lid,*) "LISTING FILE OF MACH2D"
    write(lid,*)

    call write_parameters(lid)

    call allocate_variables

    ! Initializing variables

    ito = 1

    tcpu = 0.d0

    tcpuo = 0.d0

    x = 0.d0
    y = 0.d0
    xp = 0.d0
    yp = 0.d0

    xe  = 0.d0
    ye  = 0.d0
    xen = 0.d0
    yen = 0.d0
    xk  = 0.d0
    yk  = 0.d0
    xke = 0.d0
    yke = 0.d0
    Jp  = 0.d0
    Je  = 0.d0
    Jn  = 0.d0
    alphae = 0.d0
    gamman = 0.d0
    betae  = 0.d0
    betan  = 0.d0

    radius = 0.d0
    re = 0.d0
    rn = 0.d0
    rp = 0.d0

    p   = 0.d0
    T   = 0.d0
    ro  = 0.d0
    roa = 0.d0
    roe = 0.d0
    ron = 0.d0

    Uce = 0.d0
    Vcn = 0.d0

    cp  = 0.d0
    gcp = 0.d0
    kp  = 0.d0
    ke  = 0.d0
    kn  = 0.d0
    vle = 0.d0
    vln = 0.d0
    vlp = 0.d0

    u   = 0.d0
    v   = 0.d0
    ua  = 0.d0
    va  = 0.d0

    cup = 0.d0
    cvp = 0.d0
    sup = 0.d0
    svp = 0.d0

    pl  = 0.d0
    due = 0.d0
    dve = 0.d0
    dun = 0.d0
    dvn = 0.d0
    de  = 0.d0
    dn  = 0.d0
    ue  = 0.d0
    ve  = 0.d0
    un  = 0.d0
    vn  = 0.d0
    uea  = 0.d0
    vea  = 0.d0
    una  = 0.d0
    vna  = 0.d0

    au  = 0.d0
    av  = 0.d0
    at  = 0.d0
    ap  = 0.d0
    bp  = 0.d0
    bu  = 0.d0
    b   = 0.d0

    Twall = 0.d0
    Tin   = 0.d0
    pin   = 0.d0
    pina  = 0.d0
    plin  = 0.d0
    uin   = 0.d0
    vin   = 0.d0
    Mw    = 0.d0


    u1D = 0.d0
    M1D = 0.d0
    p1D = 0.d0
    T1D = 0.d0

    ! Variables used in the eddy viscosity calculation
    vtp    = 0.d0
    ro_av  = 0.d0
    ron_av = 0.d0
    u_av   = 0.d0
    v_av   = 0.d0
    rona   = 0.d0

  end subroutine initialize_variables


  subroutine date_time(date,time)
    implicit none
    character (len=20), intent(out) :: date
    character (len=20), intent(out) :: time
    ! Auxiliary variabless
    integer(4)   :: var(8) ! date and time
    character*20 :: vardate,vartime,varzone ! date and time
    character*2  :: aux1,aux2
    character*4  :: aux3
    character*50 :: aux

    call date_and_time(vardate,vartime,varzone,var)

    write(aux,*) var(3)
    aux1 = trim(adjustl(aux))
    write(aux,*) var(2)
    aux2 = trim(adjustl(aux))
    write(aux,*) var(1)
    aux3 = trim(adjustl(aux))
    date = '('//trim(aux1)//'/'//trim(aux2)//'/'//aux3//')'

    write(aux,*) var(5)
    aux1 = trim(adjustl(aux))
    write(aux,*) var(6)
    aux2 = trim(adjustl(aux))
    write(aux,*) var(7)
    aux3 = trim(adjustl(aux))
    time = trim(aux1)//':'//trim(aux2)//':'//aux3

  end subroutine date_time

end module data
