!> \brief Contains user-DEPENDENT subroutines and functions
!! such as initial and boundary conditions as well as fluid thermophysical properties.
!<

module user

   implicit none
contains

   !> \brief Generates the north and south boundaries of the grid based on a power-law model.
   subroutine get_grid_boundary_power_law( nx, ny, akn, aks, la, lb, lr, rb, x, y) ! Last two are output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: akn !< Exponent of the power law for the north boundary
      real(8), intent(in) :: aks !< Exponent of the power law for the south boundary
      real(8), intent(in) :: la  !< Length of the elliptical x semi-axis (m)
      real(8), intent(in) :: lb  !< Length of the elliptical y semi-axis (m)
      real(8), intent(in) :: lr  !< length of the body (m)
      real(8), intent(in) :: rb  !< base radius/semi-height of the body (m)

      real(8), dimension(nx*ny), intent(out) :: x !< Coord. x at the northeast corner of the volume P (m)
      real(8), dimension(nx*ny), intent(out) :: y !< Coord. y at the northeast corner of the volume P (m)

      ! Constants

      real(8), parameter :: pi = acos(-1.d0)

      ! Auxiliary variables
      integer :: i, j, np

      ! North boundary

      j = ny-1

      do i = 1, nx-1

         np   = nx * (j-1) + i

         x(np) = la * ( dble(i-1) / dble(nx-2) ) ** akn - (la-lr)

         y(np) = lb * sqrt( 1.d0 - ( (x(np)-lr) / la )**2.d0 )

      end do

      ! South boundary

      j = 1

      do i = 1, nx-1

         np   = nx * (j-1) + i

         x(np) = lr * ( dble(i-1) / dble(nx-2) ) ** aks

         y(np) = f_body_shape( x(np), lr, rb )

      end do

   end subroutine get_grid_boundary_power_law

   !> \brief Defines the function g(x) of the body shape
   real(8) function f_body_shape ( x, lr, rb)
      implicit none
      real(8), intent(in) :: x   !< x coordinate (m)
      real(8), intent(in) :: lr  !< length of the body (m)
      real(8), intent(in) :: rb  !< base radius/semi-height of the body (m)

      f_body_shape = rb * ( x / lr ) ** 1.d0

   end function f_body_shape


   !> \brief Calculates cp at the center of each real volume
   subroutine set_cp( nx, ny, cp) ! Last is output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(out) :: cp  !< Specific heat at const pressure (J/(kg.K))

      cp = 1003.8d0 ! (J/kg.K) dry air, sea level ??, T=275K

   end subroutine set_cp

   !> \brief Calculates gamma = Cp / Cv at the center of each real volume based on Cp and Rg.
   subroutine set_gamma( nx, ny, Rg, cp, gcp) ! Last is output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: Rg  !< Perfect gas constant (J/(kg.K))
      real(8), dimension(nx*ny), intent(in) :: cp  !< Specific heat at const pressure (J/(kg.K))
      real(8), dimension(nx*ny), intent(out) :: gcp !< gcp = gamma = Cp/Cv at center of CV P (dimensionless)

      gcp = cp / ( cp - Rg ) ! dimensionless

   end subroutine set_gamma



   !> \brief Calculates the coefficients and source of the linear system of u for the fictitious volumes
   !! based on the boundary conditions.
   subroutine set_bcu( nx, ny, UF, xk, yk, alphae, betae, u, v, Uce, Vcw, a, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: UF  !< Free stream speed (m/s)
      real(8), dimension (nx*ny),   intent(in)  :: xk !< x_csi at the center of north face  of volume P (m)
      real(8), dimension (nx*ny),   intent(in)  :: yk !< y_csi at the center of north face of volume P (m)
      real(8), dimension(nx*ny),    intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),    intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension (nx*ny),   intent(in)  :: u  !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension (nx*ny),   intent(in)  :: v  !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension(nx*ny),    intent(in)  :: Uce !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),       intent(in)  :: Vcw !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension (nx*ny,9), intent(out) :: a   !< Coefficients of the linear system for u (kg/s)
      real(8), dimension (nx*ny),   intent(out) :: b   !< Source of the linear system for u (N)

      ! Auxiliary variables
      integer :: i, j, np, npw, npe, nps, npn, npse, npsw, npne, npnw

      integer :: npnnw, npssw, npnne, npsse, npss, npnn

      real(8) :: aux, u_eta

      ! West boundary

      i = 1

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1
      npnn = npn + nx
      npnne= npnn + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      a(np,6) = -1.d0

      u_eta = 0.5d0 * ( - u(npnne) + 4.d0 * u(npne) - 3.d0 * u(npe) )

      b(np) = - betae(np) / alphae(np) * u_eta

      do j = 3, ny-2

         np = nx * (j-1) + i

         nps  = np - nx

         npn  = np + nx

         npse = nps + 1

         npne = npn + 1

         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,6) = -1.d0

         u_eta = 0.5d0 * ( u(npne) - u(npse) )

         b(np) = - betae(np) / alphae(np) * u_eta

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1
      npss = nps - nx
      npsse= npss + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      a(np,6) = -1.d0

      u_eta = 0.5d0 * ( u(npsse) - 4.d0 * u(npse) + 3.d0 * u(npe) )

      b(np) = - betae(np) / alphae(np) * u_eta






      ! East boundary

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      a(np,:) =  0.d0
      a(np,4) = -1.d0
      a(np,5) =  1.d0

      u_eta = 0.5d0 * ( - u(npnnw) + 4.d0 * u(npnw) - 3.d0 * u(npw) )

      b(np) = - Vcw(j) / Uce(npw) * u_eta


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         a(np,:) =  0.d0
         a(np,4) = -1.d0
         a(np,5) =  1.d0

         u_eta = 0.5d0 * ( u(npnw) - u(npsw) )

         b(np) = - Vcw(j) / Uce(npw) * u_eta

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      a(np,:) =  0.d0
      a(np,4) = -1.d0
      a(np,5) =  1.d0

      u_eta = 0.5d0 * ( 3.d0 * u(npw) - 4.d0 * u(npsw) + u(npssw) )

      b(np) = - Vcw(j) / Uce(npw) * u_eta



      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i
         npn  = np + nx

         a(np,:) =  0.d0
         a(np,8) =  1.d0
         a(np,5) =  1.d0

         aux = sqrt( ( u(npn)**2 + v(npn)**2 ) / ( xk(np)**2 + yk(np)**2 ) ) &

            * sign( 1.d0, u(npn) * xk(np) + v(npn) * yk(np) )

         b(np) = 2.d0 * aux * xk(np)

      end do

      ! North boundary

      j = ny

      do i = 2, nx-1

         np = nx * (j-1) + i

         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,2) = 1.d0

         b(np) = 2.d0 * UF

      end do

      ! Corners (extrapolation)

      ! SW
      i = 1; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( u(npn) + u(npe) + u(npne) ) / 3.d0


      !SE
      i = nx; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( u(npn) + u(npw) + u(npnw) ) / 3.d0

      ! NW
      i = 1; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( u(nps) + u(npe) + u(npse) ) / 3.d0


      ! NE
      i = nx; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( u(nps) + u(npw) + u(npsw) ) / 3.d0

   end subroutine set_bcu

   !> Calculates the coefficients and source of the linear system of v for the fictitious volumes
   !! based on the boundary conditions.
   subroutine set_bcv( nx, ny, xk, yk, u, v, Uce, Vcw, a, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension (nx*ny),   intent(in)  :: xk  !< x_csi at the center of north face  of volume P (m)
      real(8), dimension (nx*ny),   intent(in)  :: yk  !< y_csi at the center of north face of volume P (m)
      real(8), dimension (nx*ny),   intent(in)  :: u   !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension (nx*ny),   intent(in)  :: v   !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension(nx*ny),    intent(in)  :: Uce !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),       intent(in)  :: Vcw !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension (nx*ny,9), intent(out) :: a   !< Coefficients of the linear system for v (kg/s)
      real(8), dimension (nx*ny),   intent(out) :: b   !< Source of the linear system for v (N)

      ! Auxiliary variables

      integer :: i, j, np, npe, npw, npn, nps, npse, npsw, npnw, npne

      integer :: npssw, npnnw

      real(8) :: aux, v_eta

      ! West boundary

      i = 1

      do j = 2, ny-1

         np = nx * (j-1) + i

         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,6) = 1.d0

         b(np) = 0.d0

      end do




      ! East boundary

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      a(np,:) =  0.d0
      a(np,4) = -1.d0
      a(np,5) =  1.d0

      v_eta = 0.5d0 * ( - v(npnnw) + 4.d0 * v(npnw) - 3.d0 * v(npw) )

      b(np) = - Vcw(j) / Uce(npw) * v_eta


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         a(np,:) =  0.d0
         a(np,4) = -1.d0
         a(np,5) =  1.d0

         v_eta = 0.5d0 * ( v(npnw) - v(npsw) )

         b(np) = - Vcw(j) / Uce(npw) * v_eta

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      a(np,:) =  0.d0
      a(np,4) = -1.d0
      a(np,5) =  1.d0

      v_eta = 0.5d0 * ( 3.d0 * v(npw) - 4.d0 * v(npsw) + v(npssw) )

      b(np) = - Vcw(j) / Uce(npw) * v_eta




      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i
         npn  = np + nx

         a(np,:) = 0.d0
         a(np,5) = 1.d0
         a(np,8) = 1.d0

         aux = sqrt( ( u(npn)**2 + v(npn)**2 ) / ( xk(np)**2 + yk(np)**2 ) ) &

            * sign( 1.d0, u(npn) * xk(np) + v(npn) * yk(np) )

         b(np) = 2.d0 * aux * yk(np)

      end do

      ! North boundary

      j = ny

      do i = 2, nx-1

         np = nx * (j-1) + i

         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,2) = 1.d0

         b(np) = 0.d0

      end do

      ! Corners (extrapolation)

      ! SW
      i = 1; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( v(npn) + v(npe) + v(npne) ) / 3.d0


      !SE
      i = nx; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( v(npn) + v(npw) + v(npnw) ) / 3.d0

      ! NW
      i = 1; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( v(nps) + v(npe) + v(npse) ) / 3.d0


      ! NE
      i = nx; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( v(nps) + v(npw) + v(npsw) ) / 3.d0

   end subroutine set_bcv

   !> Calculates the coefficients  and source of the linear system of T for the fictitious volumes
   !! based on the boundary conditions.
   subroutine set_bcT( nx, ny, TF, Uce, Vcw, T, alphae, betae, betan, gamman, a, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: TF     !< Free stream temperature (K)
      real(8), dimension(nx*ny),   intent(in)  :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),      intent(in)  :: Vcw    !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension(nx*ny),   intent(in)  :: T      !< Temperature at center of volume P (K)
      real(8), dimension(nx*ny),   intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betan  !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: gamman !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), dimension(nx*ny,9), intent(out) :: a  !< Coefficients of the linear system for T (J/(K.s))
      real(8), dimension(nx*ny),   intent(out) :: b  !< Source of the linear system for T (J/s)

      ! Auxiliary variables
      integer :: i, j, np, npe, npw, nps, npn, npse, npsw, npne, npnw

      integer :: npnn, npnne, npss, npsse, npnnw, npssw

      real(8) :: Te, Tk

      ! West boundary

      i = 1

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1
      npnn = npn + nx
      npnne= npnn + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      a(np,6) = -1.d0

      Te = 0.5d0 * ( - T(npnne) + 4.d0 * T(npne) - 3.d0 * T(npe) )

      b(np) = - betae(np) / alphae(np) * Te

      do j = 3, ny-2

         np = nx * (j-1) + i

         nps  = np - nx

         npn  = np + nx

         npse = nps + 1

         npne = npn + 1

         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,6) = -1.d0

         Te = 0.5d0 * ( T(npne) - T(npse) )

         b(np) = - betae(np) / alphae(np) * Te

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1
      npss = nps - nx
      npsse= npss + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      a(np,6) = -1.d0

      Te = 0.5d0 * ( T(npsse) - 4.d0 * T(npse) + 3.d0 * T(npe) )

      b(np) = - betae(np) / alphae(np) * Te




      ! East boundary

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      a(np,:) =  0.d0
      a(np,4) = -1.d0
      a(np,5) =  1.d0

      Te = 0.5d0 * ( - T(npnnw) + 4.d0 * T(npnw) - 3.d0 * T(npw) )

      b(np) = - Vcw(j) / Uce(npw) * Te


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         a(np,:) =  0.d0
         a(np,4) = -1.d0
         a(np,5) =  1.d0

         Te = 0.5d0 * ( T(npnw) - T(npsw) )

         b(np) = - Vcw(j) / Uce(npw) * Te

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      a(np,:) =  0.d0
      a(np,4) = -1.d0
      a(np,5) =  1.d0

      Te = 0.5d0 * ( 3.d0 * T(npw) - 4.d0 * T(npsw) + T(npssw) )

      b(np) = - Vcw(j) / Uce(npw) * Te

      ! South boundary

      j = 1

      do i = 2, nx-1

         np = nx * (j-1) + i

         npn  = np + nx

         npw  = np - 1

         npe  = np + 1

         npnw = npn - 1

         npne = npn + 1


         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,8) = -1.d0

         Tk = 0.25d0 * ( T(npe) + T(npne) - T(npw) - T(npnw) )

         b(np) = -betan(np) / gamman(np) * Tk

      end do

      ! North boundary

      j = ny

      do i = 2, nx-1

         np = nx * (j-1) + i

         a(np,:) = 0.d0

         a(np,5) = 1.d0

         a(np,2) = 1.d0

         b(np) = 2.d0 * TF

      end do

      ! Corners (extrapolation)

      ! SW
      i = 1; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( T(npn) + T(npe) + T(npne) ) / 3.d0


      !SE
      i = nx; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( T(npn) + T(npw) + T(npnw) ) / 3.d0

      ! NW
      i = 1; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( T(nps) + T(npe) + T(npse) ) / 3.d0


      ! NE
      i = nx; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1

      a(np,:) = 0.d0

      a(np,5) = 1.d0

      b(np) = ( T(nps) + T(npw) + T(npsw) ) / 3.d0

   end subroutine set_bcT


   !> Calculates the coefficients and source of the linear system of the pressure deviation
   !! for the fictitious volumes based on the boundary conditions
   subroutine set_bcp( nx, ny, alphae, betae, betan, gamman, Uce, Vcw, p, a, b) ! Output: last two
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny),   intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betan  !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: gamman !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),      intent(in)  :: Vcw    !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension(nx*ny),   intent(in)  :: p      !< Pressure at center o volume P (Pa)
      real(8), dimension(nx*ny,5), intent(out) :: a      !< Coefficients of the linear system for pl (m.s)
      real(8), dimension(nx*ny),   intent(out) :: b      !< Source of the linear system for pl (kg/s)

      ! Auxliary variables

      integer :: i, j, np, npe, npw, npn, nps, npnw, npne, npse, npsw

      integer :: npnee, npnww, npnne, npsse, npnnw, npssw

      real(8) :: pk, pe

      ! West boundary

      i = 1

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1
      npnne= npne + nx

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      a(np,4) = -1.d0

      pe = ( -p(npnne) + 4.d0 * p(npne) - 3.d0 * p(npe) ) / 2.d0

      b(np) = p(npe) - p(np) - betae(np) / alphae(np) * pe


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npe  = np + 1
         npse = nps + 1
         npne = npn + 1

         a(np,:) = 0.d0

         a(np,3) = 1.d0

         a(np,4) = -1.d0

         pe = ( p(npne) - p(npse) ) / 2.d0

         b(np) = p(npe) - p(np) - betae(np) / alphae(np) * pe

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1
      npsse= npse - nx

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      a(np,4) = -1.d0

      pe = ( 3.d0 * p(npe) - 4.d0 * p(npse) + p(npsse) ) / 2.d0

      b(np) = p(npe) - p(np) - betae(np) / alphae(np) * pe




      ! East boundary

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      a(np,:) =  0.d0
      a(np,2) = -1.d0
      a(np,3) =  1.d0

      pe = ( -p(npnnw) + 4.d0 * p(npnw) - 3.d0 * p(npw) ) / 2.d0

      b(np) =  p(npw) - p(np) - Vcw(j) / Uce(npw) * pe


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         a(np,:) =  0.d0
         a(np,2) = -1.d0
         a(np,3) =  1.d0

         pe = ( p(npnw) - p(npsw) ) / 2.d0

         b(np) =  p(npw) - p(np) - Vcw(j) / Uce(npw) * pe

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      a(np,:) =  0.d0
      a(np,2) = -1.d0
      a(np,3) =  1.d0

      pe = ( 3.d0 * p(npw) - 4.d0 * p(npsw) + p(npssw) ) / 2.d0

      b(np) =  p(npw) - p(np) - Vcw(j) / Uce(npw) * pe




      ! South boundary

      j = 1

      i = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npne = npn + 1
      npnee = npne + 1

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      a(np,5) = -1.d0

      pk = ( - p(npnee) + 4.d0 * p(npne) - 3.d0 * p(npn) ) / 2.d0

      b(np) = p(npn) - p(np) - betan(np) / gamman(np) * pk


      do i = 3, nx-2

         np = nx * (j-1) + i
         npn  = np + nx
         npnw = npn - 1
         npne = npn + 1

         a(np,:) = 0.d0

         a(np,3) = 1.d0

         a(np,5) = -1.d0

         pk = ( p(npne) - p(npnw) ) / 2.d0

         b(np) = p(npn) - p(np) - betan(np) / gamman(np) * pk

      end do

      i = nx-1

      np   = nx * (j-1) + i
      npn  = np + nx
      npnw = npn - 1
      npnww= npnw - 1

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      a(np,5) = -1.d0

      pk = ( 3.d0 * p(npn) - 4.d0 * p(npnw) + p(npnww) ) / 2.d0

      b(np) = p(npn) - p(np) - betan(np) / gamman(np) * pk

      ! North boundary

      j = ny

      do i = 2, nx-1

         np = nx * (j-1) + i

         a(np,:) = 0.d0

         a(np,3) = 1.d0

         a(np,1) = 1.d0

         b(np) = 0.d0

      end do


      ! Corners

      ! SW
      i = 1; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      b(np) = 0.d0


      !SE
      i = nx; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      b(np) = 0.d0


      ! NW
      i = 1; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      b(np) = 0.d0


      ! NE
      i = nx; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1

      a(np,:) = 0.d0

      a(np,3) = 1.d0

      b(np) = 0.d0

   end subroutine set_bcp

   !> \brief Calculates the contravariant velocity V on the west face
   !! of the fictitious volumes of the east boundary
   subroutine get_Vcw(nx, ny, xe, xk, ue, ve, Vcw)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(in)  :: xe  !< x_eta at the center of east face of volume P (m)
      real(8), dimension (nx*ny), intent(in)  :: xk  !< x_csi at the center of north face  of volume P (m)
      real(8), dimension (nx*ny), intent(in)  :: ue  !< Cartesian velocity u at center of east face of volume P (m/s)
      real(8), dimension (nx*ny), intent(in)  :: ve  !< Cartesian velocity v at center of east face of volume P (m/s)
      real(8), dimension(ny),     intent(out) :: Vcw !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)

      ! Inner variables

      integer :: i, j, np

      i = nx-1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         Vcw(j) = ve(np) * xk(np) - ue(np) * xe(np)

      end do

   end subroutine get_Vcw


   !> Extrapolates the corrected velocities of the real volumes to the fictitious volumes
   !! based on the boundary conditions.
   subroutine get_u_v_extrapolation_to_fictitious( nx, ny, UF, xk, yk, alphae, betae, Uce, Vcw, u, v) ! Last two are inout
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: UF  !< Free stream speed (m/s)
      real(8), dimension (nx*ny),   intent(in)  :: xk   !< x_csi at the center of north face  of volume P (m)
      real(8), dimension (nx*ny),   intent(in)  :: yk   !< y_csi at the center of north face of volume P (m)
      real(8), dimension(nx*ny),    intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),    intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)   :: Uce  !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),      intent(in)   :: Vcw  !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension (nx*ny), intent(inout) :: u    !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension (nx*ny), intent(inout) :: v    !< Cartesian velocity at the center of volume P (m/s)

      ! Auxiliary variables

      integer :: i, j, np, npe, npw, nps, npn

      integer :: npssw, npnnw, npsw, npnw, npsse, npnne, npss, npnn, npse, npne

      real(8) :: aux, u_eta, v_eta

      ! West boundary ( only for u )

      i = 1

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1
      npnn = npn + nx
      npnne= npnn + 1

      u_eta = 0.5d0 * ( - u(npnne) + 4.d0 * u(npne) - 3.d0 * u(npe) )

      u(np) = u(npe) - betae(np) / alphae(np) * u_eta

      do j = 3, ny-2

         np = nx * (j-1) + i

         npe  = np + 1

         nps  = np - nx

         npn  = np + nx

         npse = nps + 1

         npne = npn + 1

         u_eta = 0.5d0 * ( u(npne) - u(npse) )

         u(np) = u(npe) - betae(np) / alphae(np) * u_eta

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1
      npss = nps - nx
      npsse= npss + 1

      u_eta = 0.5d0 * ( u(npsse) - 4.d0 * u(npse) + 3.d0 * u(npe) )

      u(np) = u(npe) - betae(np) / alphae(np) * u_eta



      ! West boundary ( only for v )

      i = 1

      do j = 2, ny-1


         np   = nx * (j-1) + i
         npe  = np + 1

         v(np) = - v(npe)

      end do




      ! East boundary

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      u_eta = 0.5d0 * ( - u(npnnw) + 4.d0 * u(npnw) - 3.d0 * u(npw) )
      v_eta = 0.5d0 * ( - v(npnnw) + 4.d0 * v(npnw) - 3.d0 * v(npw) )

      u(np) = u(npw) - Vcw(j) / Uce(npw) * u_eta
      v(np) = v(npw) - Vcw(j) / Uce(npw) * v_eta


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         u_eta = 0.5d0 * ( u(npnw) - u(npsw) )
         v_eta = 0.5d0 * ( v(npnw) - v(npsw) )

         u(np) = u(npw) - Vcw(j) / Uce(npw) * u_eta
         v(np) = v(npw) - Vcw(j) / Uce(npw) * v_eta

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      u_eta = 0.5d0 * ( 3.d0 * u(npw) - 4.d0 * u(npsw) + u(npssw) )
      v_eta = 0.5d0 * ( 3.d0 * v(npw) - 4.d0 * v(npsw) + v(npssw) )

      u(np) = u(npw) - Vcw(j) / Uce(npw) * u_eta
      v(np) = v(npw) - Vcw(j) / Uce(npw) * v_eta


      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i
         npn  = np + nx

         aux = sqrt( ( u(npn)**2 + v(npn)**2 ) / ( xk(np)**2 + yk(np)**2 ) ) &

            * sign( 1.d0, u(npn) * xk(np) + v(npn) * yk(np) )

         u(np) = 2.d0 * aux * xk(np) - u(npn)

         v(np) = 2.d0 * aux * yk(np) - v(npn)

      end do

      ! North boundary (far field)

      j = ny

      do i = 2, nx-1

         np   = nx * (j-1) + i

         nps  = np - nx

         u(np) = 2.d0 * UF - u(nps)

         v(np) = - v(nps)

      end do

      call get_extrapolation_to_corners(nx, ny, u)

      call get_extrapolation_to_corners(nx, ny, v)


   end subroutine get_u_v_extrapolation_to_fictitious



   !> Extrapolates the temperature of the real volumes to the fictitious volumes
   !! based on the boundary conditions.
   subroutine get_T_extrapolation_to_fictitious( nx, ny, TF, alphae, betae &
         , betan, gamman, Uce, Vcw, T) ! Last is inout
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: TF  !< Free stream temperature (K)
      real(8), dimension(nx*ny),   intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betan  !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: gamman !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),      intent(in)  :: Vcw    !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension(nx*ny), intent(inout) :: T !< Temperature at center of volume P (K)

      ! Auxiliary variables

      integer :: i, j, np, npe, npw, nps, npn, npnw, npne, npse

      integer :: npnn, npnne, npss, npsse, npsw, npssw, npnnw

      real(8) :: Te, Tk

      ! West boundary (frontal symmetry line)

      i = 1


      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1
      npnn = npn + nx
      npnne= npnn + 1

      Te = 0.5d0 * ( - T(npnne) + 4.d0 * T(npne) - 3.d0 * T(npe) )

      T(np) = T(npe) - betae(np) / alphae(np) * Te



      do j = 3, ny-2

         np   = nx * (j-1) + i

         nps  = np - nx

         npn  = np + nx

         npe  = np + 1

         npse = nps + 1

         npne = npn + 1

         Te = 0.5d0 * ( T(npne) - T(npse) )

         T(np) = T(npe) - betae(np) / alphae(np) * Te

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1
      npss = nps - nx
      npsse= npss + 1

      Te = 0.5d0 * ( T(npsse) - 4.d0 * T(npse) + 3.d0 * T(npe) )

      T(np) = T(npe) - betae(np) / alphae(np) * Te



      ! East face

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      Te = 0.5d0 * ( - T(npnnw) + 4.d0 * T(npnw) - 3.d0 * T(npw) )

      T(np) = T(npw) - Vcw(j) / Uce(npw) * Te


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         Te = 0.5d0 * ( T(npnw) - T(npsw) )

         T(np) = T(npw) - Vcw(j) / Uce(npw) * Te

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npn  = np + nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      Te = 0.5d0 * ( 3.d0 * T(npw) - 4.d0 * T(npsw) + T(npssw) )

      T(np) = T(npw) - Vcw(j) / Uce(npw) * Te

      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         npn  = np + nx

         npw  = np - 1

         npe  = np + 1

         npnw = npn - 1

         npne = npn + 1


         Tk = 0.25d0 * ( T(npe) + T(npne) - T(npw) - T(npnw) )

         T(np) = T(npn) - betan(np) / gamman(np) * Tk

      end do

      ! North boundary (far field)

      j = ny

      do i = 2, nx-1

         np   = nx * (j-1) + i

         nps  = np - nx

         T(np) = 2.d0 * TF - T(nps)

      end do

      call get_extrapolation_to_corners(nx, ny, T)

   end subroutine get_T_extrapolation_to_fictitious



   !> Extrapolates the pressure of the real volumes to the fictitious volumes
   !! based on the boundary conditions.
   subroutine get_p_extrapolation_to_fictitious( nx, ny, PF, alphae, betae, betan, gamman, Uce, Vcw, p) ! Last is inout
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: PF  !< Free stream pressure (Pa)
      real(8), dimension(nx*ny),   intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: betan  !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: gamman !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), dimension(nx*ny),   intent(in)  :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(ny),      intent(in)  :: Vcw    !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)
      real(8), dimension(nx*ny), intent(inout) :: p !< Pressure at center of volume P (K)

      ! Auxiliary variables

      integer :: i, j, np, npe, npw, nps, npn

      integer :: npne, npnw, npse, npsw, npnww, npnee, npnne, npsse, npnnw, npssw

      real(8) :: pk, pe

      ! West boundary

      i = 1

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1
      npnne= npne + nx

      pe = ( -p(npnne) + 4.d0 * p(npne) - 3.d0 * p(npe) ) / 2.d0

      p(np) = p(npe) - betae(np) / alphae(np) * pe


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npe  = np + 1
         npse = nps + 1
         npne = npn + 1

         pe = ( p(npne) - p(npse) ) / 2.d0

         p(np) = p(npe) - betae(np) / alphae(np) * pe

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1
      npsse= npse - nx

      pe = ( 3.d0 * p(npe) - 4.d0 * p(npse) + p(npsse) ) / 2.d0

      p(np) = p(npe) - betae(np) / alphae(np) * pe




      ! East boundary

      i = nx

      j = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1
      npnnw= npnw + nx

      pe = ( -p(npnnw) + 4.d0 * p(npnw) - 3.d0 * p(npw) ) / 2.d0

      p(np) =  p(npw) - Vcw(j) / Uce(npw) * pe


      do j = 3, ny-2

         np   = nx * (j-1) + i
         nps  = np - nx
         npn  = np + nx
         npw  = np - 1
         npsw = nps - 1
         npnw = npn - 1

         pe = ( p(npnw) - p(npsw) ) / 2.d0

         p(np) =  p(npw) - Vcw(j) / Uce(npw) * pe

      end do

      j = ny-1

      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1
      npssw= npsw - nx

      pe = ( 3.d0 * p(npw) - 4.d0 * p(npsw) + p(npssw) ) / 2.d0

      p(np) =  p(npw) - Vcw(j) / Uce(npw) * pe





      ! South boundary

      j = 1

      i = 2

      np   = nx * (j-1) + i
      npn  = np + nx
      npne = npn + 1
      npnee = npne + 1

      pk = ( - p(npnee) + 4.d0 * p(npne) - 3.d0 * p(npn) ) / 2.d0

      p(np) = p(npn) - betan(np) / gamman(np) * pk


      do i = 3, nx-2

         np = nx * (j-1) + i
         npn  = np + nx
         npnw = npn - 1
         npne = npn + 1

         pk = ( p(npne) - p(npnw) ) / 2.d0

         p(np) = p(npn) - betan(np) / gamman(np) * pk

      end do

      i = nx-1

      np   = nx * (j-1) + i
      npn  = np + nx
      npnw = npn - 1
      npnww= npnw - 1

      pk = ( 3.d0 * p(npn) - 4.d0 * p(npnw) + p(npnww) ) / 2.d0

      p(np) = p(npn) - betan(np) / gamman(np) * pk





      ! North boundary (far field)

      j = ny

      do i = 2, nx-1

         np   = nx * (j-1) + i

         nps  = np - nx

         p(np) = 2.d0 * PF - p(nps)

      end do

      call get_extrapolation_to_corners(nx, ny, p)

   end subroutine get_p_extrapolation_to_fictitious


   !> \brief Extrapolation of a field to the fictitious volumes of the corners
   subroutine get_extrapolation_to_corners(nx, ny, f)
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(inout) :: f !< Field to be extrapolated to the corners

      ! Auxiliary variables

      integer :: i, j, np, npe, npw, npn, nps, npne, npnw, npse, npsw


      ! Corners (extrapolation)

      ! SW
      i = 1; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npe  = np + 1
      npne = npn + 1

      f(np) = ( f(npn) + f(npe) + f(npne) ) / 3.d0

      !SE
      i = nx; j = 1
      np   = nx * (j-1) + i
      npn  = np + nx
      npw  = np - 1
      npnw = npn - 1

      f(np) = ( f(npn) + f(npw) + f(npnw) ) / 3.d0

      ! NW
      i = 1; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npe  = np + 1
      npse = nps + 1

      f(np) = ( f(nps) + f(npe) + f(npse) ) / 3.d0


      ! NE
      i = nx; j = ny
      np   = nx * (j-1) + i
      nps  = np - nx
      npw  = np - 1
      npsw = nps - 1

      f(np) = ( f(nps) + f(npw) + f(npsw) ) / 3.d0

   end subroutine



   !> Calculates the SIMPLEC coefficients at the boundary faces.
   subroutine get_boundary_simplec_coefficients( nx, ny, xe, ye &
      ,        due, dve, dun, dvn, de, dn) ! InOutput: last six
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(in) :: xe     !< x_eta at the center of east face of volume P (m)
      real(8), dimension (nx*ny), intent(in) :: ye     !< y_eta at the center of east face of volume P (m)
      real(8), dimension (nx*ny), intent(inout) :: due !< SIMPLEC coefficients for ue (m2.s/kg)
      real(8), dimension (nx*ny), intent(inout) :: dve !< SIMPLEC coefficients for ve (m2.s/kg)
      real(8), dimension (nx*ny), intent(inout) :: dun !< SIMPLEC coefficients for un (m2.s/kg)
      real(8), dimension (nx*ny), intent(inout) :: dvn !< SIMPLEC coefficients for vn (m2.s/kg)
      real(8), dimension (nx*ny), intent(inout) :: de  !< SIMPLEC coefficients for Uce (m3.s/kg)
      real(8), dimension (nx*ny), intent(inout) :: dn  !< SIMPLEC coefficients for Vcn (m3.s/kg)

      ! Auxiliary variables
      integer :: i, j, np, npe, npw, npn

      ! West boundary

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         due(np) = due(npe)

         dve(np) = 0.d0

         de(np) = 0.d0

      end do

      ! East boundary

      i = nx-1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npw  = np - 1

         due(np) = due(npw)

         dve(np) = dve(npw)

         de(np) = ye(np) * due(np) + xe(np) * dve(np)

      end do

      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         npn  = np + nx

         dun(np) = dun(npn)

         dvn(np) = dvn(npn)

         dn(np) = 0.d0

      end do

      ! North boundary

      j = ny-1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         dun(np) = 0.d0

         dvn(np) = 0.d0

         dn(np) = 0.d0

      end do

   end subroutine get_boundary_simplec_coefficients


   !> Calculates the velocities at boundary faces based on the boundary conditions
   subroutine get_velocities_at_boundary_faces( nx, ny, UF, xe, ye, xk, yk, u, v, ue, ve, un, vn, Uce, Vcn) ! Last six are output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: UF  !< Free stream speed (m/s)
      real(8), dimension(nx*ny), intent(in)  :: xe     !< x_eta at the center of east face of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: ye     !< y_eta at the center of east face of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: xk     !< x_csi at the center of north face  of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: yk     !< y_csi at the center of north face of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: u      !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension(nx*ny), intent(in)  :: v      !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: ue     !< Cartesian velocity u at center of east face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: ve     !< Cartesian velocity v at center of east face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: un     !< Cartesian velocity u at center of north face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: vn     !< Cartesian velocity v at center of north face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(nx*ny), intent(out) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)

      ! Axuliary variables

      integer :: i, j, np, npe, npn

      real(8) :: aux

      ! West boundary ( frontal symmetry line )

      i = 1

      do j = 2, ny-1

         np  = nx * (j-1) + i

         npe = np + 1

         ue(np) = ( u(np) + u(npe) ) / 2.d0

         ve(np) = 0.d0

         Uce(np) = 0.d0 ! There is no mass flow through the symmetry line

      end do

      ! East boundary

      i = nx-1

      do j = 2, ny-1

         np  = nx * (j-1) + i

         npe = np + 1

         ue(np) = ( u(np) + u(npe) ) / 2.d0

         ve(np) = ( v(np) + v(npe) ) / 2.d0

         Uce(np) = ue(np) * ye(np) - ve(np) * xe(np)

      end do

      ! South boundary (body surface)

      j = 1

      do i = 2, nx-1

         np  = nx * (j-1) + i

         npn = np + nx

         aux = sqrt( ( u(npn)**2 + v(npn)**2 ) / ( xk(np)**2 + yk(np)**2 ) ) &

            * sign( 1.d0, u(npn) * xk(np) + v(npn) * yk(np) )

         un(np) = aux * xk(np)

         vn(np) = aux * yk(np)

         Vcn(np) = 0.d0 ! There is no mass flow through the body surface

      end do

      ! North boundary (far field surface)

      j = ny-1

      do i = 2, nx-1

         np = nx * (j-1) + i

         un(np) = UF

         vn(np) = 0.d0

         Vcn(np) = - un(np) * yk(np)

      end do

   end subroutine get_velocities_at_boundary_faces


   !> Prescribes an initial conditions for some variables.
   subroutine get_initial_conditions( nx, ny, PF, TF, Rg, GF, MF, UF &
      , xe, ye, xk, yk, alphae, betae, betan, gamman &
      , p, pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, Vcw ) ! UF and last 15 are output

      implicit none
      integer, intent(in)  :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in)  :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in)  :: PF  !< Free stream pressure (Pa)
      real(8), intent(in)  :: TF  !< Free stream temperature (K)
      real(8), intent(in)  :: Rg  !< Perfect gas constant (J/(kg.K)) (J/kg.K)
      real(8), intent(in)  :: GF  !< GF = gamma = Cp / Cv (for the free stream)  (dimensionless)
      real(8), intent(in)  :: MF  !< Free stream Mach number (dimensionless)
      real(8), intent(out) :: UF  !< Free stream speed (m/s)
      real(8), dimension(nx*ny), intent(in)  :: xe     !< x_eta at the center of east face of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: ye     !< y_eta at the center of east face of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: xk     !< x_csi at the center of north face  of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: yk     !< y_csi at the center of north face of volume P (m)
      real(8), dimension(nx*ny), intent(in)  :: alphae !< (metric tensor component) alpha at the center of east face of volume P (m2)
      real(8), dimension(nx*ny), intent(in)  :: betae  !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), dimension(nx*ny), intent(in)  :: betan  !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), dimension(nx*ny), intent(in)  :: gamman !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), dimension(nx*ny), intent(out) :: p      !< Pressure at center o volume P (Pa)
      real(8), dimension(nx*ny), intent(out) :: pl     !< Pressure deviation at center o volume P (Pa)
      real(8), dimension(nx*ny), intent(out) :: T      !< Temperature at center of volume P (K)
      real(8), dimension(nx*ny), intent(out) :: ro     !< Specific mass (absolute density) at center of volume P (kg/m3)
      real(8), dimension(nx*ny), intent(out) :: roe    !< Specific mass at east face of volume P (kg/m3)
      real(8), dimension(nx*ny), intent(out) :: ron    !< Specific mass at north face of volume P (kg/m3)
      real(8), dimension(nx*ny), intent(out) :: u      !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: v      !< Cartesian velocity at the center of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: ue     !< Cartesian velocity u at center of east face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: ve     !< Cartesian velocity v at center of east face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: un     !< Cartesian velocity u at center of north face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: vn     !< Cartesian velocity v at center of north face of volume P (m/s)
      real(8), dimension(nx*ny), intent(out) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
      real(8), dimension(nx*ny), intent(out) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)
      real(8), dimension(ny),    intent(out) :: Vcw    !< Contravariant velocity V on west face of fictitious volumes of east boundary (m2/s)

      ! Auxiliary variables

      integer :: i, j, np

      real(8) :: ROF


      ! Free-stream speed
      UF = MF * sqrt( GF * Rg * TF )

      ! Free stream density
      ROF = PF / ( Rg * TF )


      ! Values of the variables at the real nodes


      u = UF

      v = 0.d0

      T = TF

      p = PF

      ro = ROF

      pl = 0.d0


      ! Values of the variables at fictitious nodes

      ! Extrapolation of u and v according the boundary conditions

      Uce = UF

      Vcw = 0.d0

      call get_u_v_extrapolation_to_fictitious( nx, ny, UF, xk, yk, alphae, betae, Uce, Vcw, u, v) ! Last two are inout

      ! Values of the variables at faces

      roe = ROF

      ron = ROF

      un = UF

      ue = UF

      vn = 0.d0

      ve = 0.d0


      ! U and V on the interface between real volumes

      do i = 2, nx-2
         do j = 2, ny-1

            np = nx * (j-1) + i

            Uce(np) = ue(np) * ye(np) - ve(np) * xe(np)

         end do
      end do

      do i = 2, nx-1
         do j = 2, ny-2

            np = nx * (j-1) + i

            Vcn(np) = vn(np) * xk(np) - un(np) * yk(np)

         end do
      end do

      ! u, v, U and V  on the boundary real faces

      call get_velocities_at_boundary_faces( nx, ny, UF, xe, ye, xk, yk &

      , u, v, ue, ve, un, vn, Uce, Vcn) ! Last six are output


      ! Calculates the contravariant velocity V on the west face
      ! of the fictitious volumes of the east boundary
      call get_Vcw(nx, ny, xe, xk, ue, ve, Vcw)


      ! Values of the variables at fictitious nodes

      ! Extrapolation of T according the boundary conditions

      call get_T_extrapolation_to_fictitious(nx, ny, TF, alphae, betae &
         , betan, gamman, Uce, Vcw, T) ! Last is inout


      ! Extrapolation of p according the boundary conditions

      call get_p_extrapolation_to_fictitious(nx, ny, PF, alphae, betae, betan, gamman, Uce, Vcw, p) ! Last is inout


      ! Extrapolation of ro according to the state equation

      ro = p / ( Rg * T )


   end subroutine get_initial_conditions

end module user
