!> \brief Contains global parameters and variables and subroutines related to them.

module data

   implicit none
   !
   ! GRID PARAMETERS AND VARIABLES
   !
   integer :: nx    !< Number of volumes in the csi direction (real+fictitious)
   integer :: ny    !< Number of volumes in the eta direction (real+fictitious)
   integer :: nxy   !< nxy = nx * ny
   integer :: kg    !< Kind of grid (1=uniform, 2=geometric progression, 3=power law)
   integer :: kcm   !< Kind of centroid mean (1=simple mean, 2=weighted mean)
   integer :: coord !< Kind of coord. system ( 1=cylindrical, 0 = cartesian)
   real(8) :: la    !< Length of the elliptical x semi-axis (m)
   real(8) :: lb    !< Length of the elliptical y semi-axis (m)
   real(8) :: a1    !< Width of the volume closer to the wall (m)
   real(8) :: akn   !< Exponent of the power law for the north boundary
   real(8) :: aks   !< Exponent of the power law for the south boundary


   real(8), allocatable, dimension(:) :: x   !< Coord. x at the northeast corner of the volume P (m)
   real(8), allocatable, dimension(:) :: y   !< Coord. y at the northeast corner of the volume P (m)
   real(8), allocatable, dimension(:) :: xp  !< Coord. x of the centroid of volume P (m)
   real(8), allocatable, dimension(:) :: yp  !< Coord. y of the centroid of volume P (m)
   real(8), allocatable, dimension(:) :: xe  !< x_eta at the center of east face of volume P (m)
   real(8), allocatable, dimension(:) :: ye  !< y_eta at the center of east face of volume P (m)
   real(8), allocatable, dimension(:) :: xen !< x_eta at the center of north face of volume P (m)
   real(8), allocatable, dimension(:) :: yen !< y_eta at the center of north face of volume P (m)
   real(8), allocatable, dimension(:) :: xk  !< x_csi at the center of north face  of volume P (m)
   real(8), allocatable, dimension(:) :: yk  !< y_csi at the center of north face of volume P (m)
   real(8), allocatable, dimension(:) :: xke !< x_csi at the center of east face of volume P (m)
   real(8), allocatable, dimension(:) :: yke !< y_csi at the center of east face of volume P (m)
   real(8), allocatable, dimension(:) :: Jp  !< Jacobian at the center of volume P (1/m2)
   real(8), allocatable, dimension(:) :: Je  !< Jacobian at the center of east face of volume P (1/m2)
   real(8), allocatable, dimension(:) :: Jn  !< Jacobian at the center of north face of volume P (1/m2)
   real(8), allocatable, dimension(:) :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
   real(8), allocatable, dimension(:) :: gamman !< (metric tensor component) gamma at the center of north face of volume P (m2)
   real(8), allocatable, dimension(:) :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
   real(8), allocatable, dimension(:) :: betan  !< (metric tensor component) beta  at the center of north face of volume P (m2)
   real(8), allocatable, dimension(:) :: re  !< Radius at the center of east face of volume P (m)
   real(8), allocatable, dimension(:) :: rn  !< Radius at the center of north face of volume P (m)
   real(8), allocatable, dimension(:) :: rp  !< Radius at the center of volume P (m)

   !
   ! NUMERIC PARAMETERS AND VARIABLES
   !
   integer :: it     !< Iteractions counter for time evolution cycle
   integer :: iti    !< Iteractions counter for iterative correction cycle
   integer :: it1    !< Number of iteractions up to which dt = dt1
   integer :: it2    !< Number of iteractions from which dt = dt2
   integer :: itb1   !< Number of iteractions up to which beta = beta1
   integer :: itb2   !< Number of iteractions from which beta = beta2
   integer :: itm    !< Iteractions counter for the mass cycle
   integer :: imax   !< Maximum number of iteractions for mass cycle
   integer :: itmax  !< Maximum number of iteractions for time evolution cycle
   integer :: itimax !< Maximum number of iteractions for the correction cycle
   integer :: nitm_u !< Maximum number of iteractions for solving the linear systems for u, v and T
   integer :: nitm_p !< Maximum number of iteractions for solving the linear system for p
   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 :: it_stop !< Number of iteractions before break time evolution cycle
   integer :: nit_res !< Number of iteractions to calculate the mean of the residuals
   integer :: res_check !< Checks if the it_stop was defined. 0 = false, 1 = true.
   real(8) :: tol_u  !< Tolerance in the MSI for solving the linear systems for u, v and T
   real(8) :: tol_p  !< Tolerance in the MSI for solving the linear system for p
   real(8) :: beta   !< UDS/CDS mixing constant (0=UDS, 1=CDS)
   real(8) :: beta1  !< Initial beta (UDS/CDS mixing constant (0=UDS, 1=CDS))
   real(8) :: beta2  !< Final beta (UDS/CDS mixing constant (0=UDS, 1=CDS))
   real(8) :: dt     !< Time step (s)
   real(8) :: dt1    !< Initial time step (s)
   real(8) :: dt2    !< Final time step (s)
   real(8) :: res_u  !< Relative residual of the linear system for u
   real(8) :: res_v  !< Relative residual of the linear system for v
   real(8) :: res_T  !< Relative residual of the linear system for T
   real(8) :: res_p  !< Residual of the linear system for p
   real(8) :: res    !< Sum of all residuals of the linear systems for u, v, T and p.
   real(8) :: res_mean !< Residual mean
   real(8) :: tol_res !< Tolerance for the sum of residuals
   real(8) :: norm   !< Norm of the residuals of all linear systems
   real(8) :: tcpu1  !< First time measurement (s)
   real(8) :: tcpu2  !< Second time measurement (s)
   real(8) :: tcpu   !< Total cpu time (s)
   real(8) :: RAM    !< RAM memory (MB)


   real(8), allocatable, dimension(:) :: bu  !< Source of the linear system for u (N)
   real(8), allocatable, dimension(:) :: bv  !< Source of the linear system for v (N)
   real(8), allocatable, dimension(:) :: bt  !< Source of the linear system for T (J/s)
   real(8), allocatable, dimension(:) :: bp  !< Source of the linear system for pl (kg/s)
   real(8), allocatable, dimension(:) :: pl  !< Pressure deviation at center o volume P (Pa)
   real(8), allocatable, dimension(:) :: due !< SIMPLEC coefficients for ue (m2.s/kg)
   real(8), allocatable, dimension(:) :: dve !< SIMPLEC coefficients for ve (m2.s/kg)
   real(8), allocatable, dimension(:) :: dun !< SIMPLEC coefficients for un (m2.s/kg)
   real(8), allocatable, dimension(:) :: dvn !< SIMPLEC coefficients for vn (m2.s/kg)
   real(8), allocatable, dimension(:) :: de  !< SIMPLEC coefficients for Uce (m3.s/kg)
   real(8), allocatable, dimension(:) :: dn  !< SIMPLEC coefficients for Vcn (m3.s/kg)

   real(8), allocatable, dimension(:,:) :: au !< Coefficients of the linear system for u (kg/s)
   real(8), allocatable, dimension(:,:) :: av !< Coefficients of the linear system for v (kg/s)
   real(8), allocatable, dimension(:,:) :: at !< Coefficients of the linear system for T (J/(K.s))
   real(8), allocatable, dimension(:,:) :: ap !< Coefficients of the linear system for pl (m.s)
   real(8), allocatable, dimension(:,:) :: dl9 !< Lower matrix of the MSI method for 9 diagonals
   real(8), allocatable, dimension(:,:) :: du9 !< Upper matrix of the MSI method for 9 diagonals
   real(8), allocatable, dimension(:,:) :: dl5 !< Lower matrix of the MSI method for 5 diagonals
   real(8), allocatable, dimension(:,:) :: du5 !< Upper matrix of the MSI method for 5 diagonals


   !
   ! GEOMETRIC PARAMETERS
   !
   real(8) :: lr  !< length of the body (m)
   real(8) :: rb  !< base radius/semi-height of the body (m)

   !
   ! FILE ID NUMBERS
   !
   integer, parameter :: tfid = 10  !< Temporary file id
   integer, parameter :: lid  = 100 !< Listing file id
   integer, parameter :: rid  = 101 !< Residual file id


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

   !
   ! GAS PROPERTIES AND VARIABLES
   !

   real(8) :: GF !< GF = gamma = Cp / Cv (for the free stream)  (dimensionless)
   real(8) :: Rg !< Perfect gas constant (J/(kg.K))
   real(8) :: PF !< Free stream pressure (Pa)
   real(8) :: TF !< Free stream temperature (K)
   real(8) :: MF !< Free stream Mach number (dimensionless)
   real(8) :: UF !< Free stream speed (m/s)

   real(8) :: Cdfi !< Pressure foredrag coefficient (dimensionless)

   real(8), allocatable, dimension(:) :: p    !< Pressure at center o volume P (Pa)
   real(8), allocatable, dimension(:) :: pa   !< Pressure of the previous time step (Pa)
   real(8), allocatable, dimension(:) :: T    !< Temperature at center of volume P (K)
   real(8), allocatable, dimension(:) :: Ta   !< Temperature of the previous time step (K)
   real(8), allocatable, dimension(:) :: ro   !< Specific mass (absolute density) at center of volume P (kg/m3)
   real(8), allocatable, dimension(:) :: roa  !< Specific mass of the previous time step (kg/m3)
   real(8), allocatable, dimension(:) :: roe  !< Specific mass at east face of volume P (kg/m3)
   real(8), allocatable, dimension(:) :: ron  !< Specific mass at north face of volume P (kg/m3)
   real(8), allocatable, dimension(:) :: g    !< ro / p (kg/J)
   real(8), allocatable, dimension(:) :: u    !< Cartesian velocity at the center of volume P (m/s)
   real(8), allocatable, dimension(:) :: ua   !< u of the previous time step (m/s)
   real(8), allocatable, dimension(:) :: v    !< Cartesian velocity at the center of volume P (m/s)
   real(8), allocatable, dimension(:) :: va   !< v of the previous time step (m/s)
   real(8), allocatable, dimension(:) :: ue   !< Cartesian velocity u at center of east face of volume P (m/s)
   real(8), allocatable, dimension(:) :: un   !< Cartesian velocity u at center of north face of volume P (m/s)
   real(8), allocatable, dimension(:) :: ve   !< Cartesian velocity v at center of east face of volume P (m/s)
   real(8), allocatable, dimension(:) :: vn   !< Cartesian velocity v at center of north face of volume P (m/s)
   real(8), allocatable, dimension(:) :: Uce  !< Contravariant velocity U at east face of volume P (m2/s)
   real(8), allocatable, dimension(:) :: Vcn  !< Contravariant velocity V at north face of volume P (m2/s)
   real(8), allocatable, dimension(:) :: cp   !< Specific heat at const pressure (J/(kg.K))
   real(8), allocatable, dimension(:) :: gcp  !< gcp = gamma = Cp/Cv at center of CV P (dimensionless)
   real(8), allocatable, dimension(:) :: cup  !< Term of deferred correction for u (N)
   real(8), allocatable, dimension(:) :: cvp  !< Term of deferred correction for v (N)
   real(8), allocatable, dimension(:) :: uea  !< ue of the previous time step (m/s)
   real(8), allocatable, dimension(:) :: vea  !< ve of the previous time step (m/s)
   real(8), allocatable, dimension(:) :: una  !< un of the previous time step (m/s)
   real(8), allocatable, dimension(:) :: vna  !< vn of the previous time step (m/s)
   real(8), allocatable, dimension(:) :: Ucea !< Uce of the previous time step (m2/s)
   real(8), allocatable, dimension(:) :: Vcna !< Vcn of the previous time step (m2/s)
   real(8), allocatable, dimension(:) :: Vcw  !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)


contains

   !> \brief Reads the parameters from an input data file, which name is in ./mach2d_input/mach2d_input.txt.
   !! Input and output variables are not explicited.
   !!
   subroutine get_parameters
      implicit none
      !
      call date_time(date, time)
      !
      open(10,file="./mach2d_input/mach2d_input.txt")
      ! 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,*) la     ! length of the elliptical x semi-axis (m)
      read(10,*) lb     ! length of the elliptical y semi-axis (m)
      read(10,*) lr     ! length of the body (m)
      read(10,*) rb     ! base radius/semi-height of the body (m)
      read(10,*) kg     ! Kind of grid (1=uniform, 2=geometric progression, 3=power law)
      read(10,*) kcm    ! Kind of centroid mean (1=simple mean, 2=weighted mean)
      read(10,*) coord  ! Kind of coord. system ( 1=cylindrical, 0 = cartesian)
      read(10,*) a1     ! width of the volume closer to the wall (m)
      read(10,*) akn    ! Exponent of the power law for the north boundary
      read(10,*) aks    ! Exponent of the power law for the south boundary
      read(10,*) itmax  ! Maximum number of iterations for time cycle
      read(10,*) itimax ! Maximum number of iteractions for the correction cycle
      read(10,*) imax   ! Maximum number of iteractions for mass cycle
      read(10,*) it1    ! Number of iteractions up to which dt = dt1
      read(10,*) it2    ! Number of iteractions from which dt = dt2
      read(10,*) dt1    ! Initial time step (s)
      read(10,*) dt2    ! Final time step (s)
      read(10,*) nitm_u ! Maximum number of iteractions for solving the linear systems for u, v and T
      read(10,*) nitm_p ! Maximum number of iteractions for solving the linear system for p
      read(10,*) nit_res! Number of iteractions to calculate the mean of the residuals
      read(10,*) tol_u  ! Tolerance in the MSI for solving the linear systems for u, v and T
      read(10,*) tol_p  ! Tolerance in the MSI for solving the linear system for p
      read(10,*) tol_res! Tolerance for the sum of residuals
      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
      read(10,*) itb1   ! Number of iteractions up to which beta = beta1
      read(10,*) itb2   ! Number of iteractions from which beta = beta2
      read(10,*) beta1  ! Initial beta (UDS/CDS mixing constant (0=UDS, 1=CDS))
      read(10,*) beta2  ! Final beta (UDS/CDS mixing constant (0=UDS, 1=CDS))
      read(10,*) GF     ! GF = gamma = Cp / Cv (for the free stream)
      read(10,*) Rg     ! Perfect gas constant (J/(kg.K))
      read(10,*) PF     ! Free stream pressure (Pa)
      read(10,*) TF     ! Free stream temperature (K)
      read(10,*) MF     ! Free stream Mach number

      close(10)
      !
      nx = nx + 2 ! Real + Fictitious
      ny = ny + 2 ! Real + Fictitious

      nxy = nx * ny

   end subroutine get_parameters


   !> \brief Writes the parameters, read from the input data file, to the file with id number fid.
   !! Input and output variables are not explicited.
   subroutine write_parameters(fid)
      implicit none
      integer, intent(in) :: fid !< File id number.

      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,"(ES14.7,' ....: ',A)") la     , " la     - length of the elliptical x semi-axis (m)"
      write(fid,"(ES14.7,' ....: ',A)") lb     , " lb     - length of the elliptical y semi-axis (m)"
      write(fid,"(ES14.7,' ....: ',A)") lr     , " lr     - length of the body (m)"
      write(fid,"(ES14.7,' ....: ',A)") rb     , " rb     - base radius/semi-height of the body (m)"
      write(fid,"(I14,' ....: ',A)")    kg     , " kg     - Kind of grid (1=uniform, 2=geometric progression, 3=power law)"
      write(fid,"(I14,' ....: ',A)")    kcm    , " kcm    - Kind of centroid mean (1=simple mean, 2=weighted mean)"
      write(fid,"(I14,' ....: ',A)")    coord  , " coord  - Kind of coord. system ( 1=cylindrical, 0 = cartesian)"
      write(fid,"(ES14.7,' ....: ',A)") a1     , " a1     - width of the volume closer to the wall (m)"
      write(fid,"(ES14.7,' ....: ',A)") akn    , " akn    - Exponent of the power law for the north boundary"
      write(fid,"(ES14.7,' ....: ',A)") aks    , " aks    - Exponent of the power law for the south boundary"
      write(fid,"(I14,' ....: ',A)")    itmax  , " itmax  - Maximum number of iterations for time cycle"
      write(fid,"(I14,' ....: ',A)")    itimax , " itimax - Maximum number of iterations for the correction cycle"
      write(fid,"(I14,' ....: ',A)")    imax   , " imax   - Maximum number of iteractions for mass cycle"
      write(fid,"(I14,' ....: ',A)")    it1    , " it1    - Number of iteractions up to which dt = dt1"
      write(fid,"(I14,' ....: ',A)")    it2    , " it2    - Number of iteractions from which dt = dt2"
      write(fid,"(ES14.7,' ....: ',A)") dt1    , " dt1    - Initial time step (s)"
      write(fid,"(ES14.7,' ....: ',A)") dt2    , " dt2    - Final time step (s)"
      write(fid,"(I14,' ....: ',A)")    nitm_u , " nitm_u - Max. # of iteractions for solving the linear systems for u, v and T"
      write(fid,"(I14,' ....: ',A)")    nitm_p , " nitm_p - Max. # of iteractions for solving the linear system for p"
      write(fid,"(I14,' ....: ',A)")    nit_res, " nit_res- Number of iteractions to calculate the mean of the residuals"
      write(fid,"(ES14.7,' ....: ',A)") tol_u  , " tol_u  - Tolerance in the MSI for solving the linear systems for u, v and T"
      write(fid,"(ES14.7,' ....: ',A)") tol_p  , " tol_p  - Tolerance in the MSI for solving the linear system for p"
      write(fid,"(ES14.7,' ....: ',A)") tol_res, " tol_res- Tolerance for the sum of residuals"
      write(fid,"(I14,' ....: ',A)")    wlf    , " Frequency of printing in the listing file"
      write(fid,"(I14,' ....: ',A)")    sem_a  , " 1 = do not open result files, 0 = open"
      write(fid,"(I14,' ....: ',A)")    sem_g  , " 0 = visualize the plot, 1 = do not visualize"
      write(fid,"(I14,' ....: ',A)")    w_g    , " Frequency of writing data for graphics"
      write(fid,"(I14,' ....: ',A)")    w_cam  , " 1 = write the fields, 0 = do not"
      write(fid,"(I14,' ....: ',A)")    itb1   , " itb1   - Number of iteractions up to which beta = beta1"
      write(fid,"(I14,' ....: ',A)")    itb2   , " itb2   - Number of iteractions from which beta = beta2"
      write(fid,"(ES14.7,' ....: ',A)") beta1  , " beta1  - Initial beta (UDS/CDS mixing constant (0=UDS, 1=CDS))"
      write(fid,"(ES14.7,' ....: ',A)") beta2  , " beta2  - Final beta (UDS/CDS mixing constant (0=UDS, 1=CDS))"
      write(fid,"(ES14.7,' ....: ',A)") GF     , " GF     - GF = gamma = Cp / Cv (for the free stream)"
      write(fid,"(ES14.7,' ....: ',A)") Rg     , " Rg     - Perfect gas constant (J/(kg.K))"
      write(fid,"(ES14.7,' ....: ',A)") PF     , " PF     - Free stream pressure (Pa)"
      write(fid,"(ES14.7,' ....: ',A)") TF     , " TF     - Free stream temperature (K)"
      write(fid,"(ES14.7,' ....: ',A)") MF     , " MF     - Free stream Mach number"
      write(fid,*)

   end subroutine write_parameters


   !> \brief Allocates the vectors.
   !! Input and output variables are not explicited.
   subroutine allocate_variables


      allocate( Vcw(ny) )

      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)               &
         ,      re(nxy),       rn(nxy),        rp(nxy),         p(nxy),          T(nxy)  &
         ,      ro(nxy),       roe(nxy),       ron(nxy),        u(nxy),          v(nxy)   &
         ,      ue(nxy),       un(nxy),        ve(nxy),         vn(nxy),         Uce(nxy) &
         ,      Vcn(nxy),      cp(nxy),        gcp(nxy),        ua(nxy),         va(nxy)  &
         ,      cup(nxy),      cvp(nxy),       g(nxy)                                     &
         ,      bu(nxy),       bv(nxy),        roa(nxy),        due(nxy),        dve(nxy) &
         ,      dun(nxy),      dvn(nxy),       de(nxy),         dn(nxy),         uea(nxy) &
         ,      vea(nxy),      una(nxy),       vna(nxy),        Ucea(nxy),       Vcna(nxy) &
         ,      bp(nxy),       pl(nxy),        Ta(nxy),         bt(nxy),         pa(nxy)   )


      allocate( au(nxy,9), av(nxy,9), at(nxy,9), ap(nxy,5), dl9(nxy,5), du9(nxy,4), dl5(nxy,4), du5(nxy,3) )


   end subroutine allocate_variables

   !> \brief Reads the parameters from an input file, allocates the vectors and initialise them.
   !! Input and output variables are not explicited.
   !!
   subroutine initialize_variables
      implicit none

      ! Reading data from input data file

      call get_parameters

      ! Openning listing file

      open(lid, file = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // ".txt")

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

      call write_parameters(lid)

      call allocate_variables

      ! Initializing variables

      it_stop = 0
      res_check = 0
      res_mean = 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

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

      p   = 0.d0
      T   = 0.d0
      ro  = 0.d0
      roe = 0.d0
      ron = 0.d0
      g   = 0.d0
      u   = 0.d0
      v   = 0.d0
      ue  = 0.d0
      un  = 0.d0
      ve  = 0.d0
      vn  = 0.d0
      Uce = 0.d0
      Vcn = 0.d0
      Vcw = 0.d0

      cp  = 0.d0
      gcp = 0.d0

      ua  = 0.d0
      va  = 0.d0
      cup = 0.d0
      cvp = 0.d0
      bu  = 0.d0
      bv  = 0.d0
      au  = 0.d0
      av  = 0.d0
      roa = 0.d0

      due = 0.d0
      dve = 0.d0
      dun = 0.d0
      dvn = 0.d0
      de  = 0.d0
      dn  = 0.d0

      uea = 0.d0
      vea = 0.d0
      una = 0.d0
      vna = 0.d0

      Ucea = 0.d0
      Vcna = 0.d0

      ap = 0.d0
      bp = 0.d0
      pl = 0.d0

      pa = 0.d0
      Ta = 0.d0
      at = 0.d0
      bt = 0.d0

      dl9 = 0.d0
      du9 = 0.d0
      dl5 = 0.d0
      du5 = 0.d0

   end subroutine initialize_variables

   !> \brief Gets the system date and time.
   subroutine date_time(date,time)
      implicit none
      character (len=20), intent(out) :: date !< system date
      character (len=20), intent(out) :: time !< system 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

   !> \brief Gets the current process memory in Mb
   real(8) function get_RAM(sim_id)
      implicit none
      character (len = *), intent(in) :: sim_id !< Simulation identification

      ! Inner variables

      character (len=20)  :: PID
      character (len=100) :: fout
      real(8) :: RAM


      write(PID,*) getPID()

      fout = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // "_memory.dat"

      call system("ps -o drs -p " // trim(adjustl(PID)) // " | tail -n 1 > " // trim(adjustl(fout)) )

      open(10,file=fout)

      read(10,*) RAM

      close(10)

      get_RAM = RAM / 1024.d0

   end function

end module data
