module user
   !
   !  SUBROUTINES
   !
   !  1) set_wall_temperature
   !  2) set_cp_and_gamma
   !  3) set_laminar_viscosity_at_nodes
   !  4) set_thermal_conductivity_at_nodes
   !  5) get_laminar_viscosity_at_faces
   !  6) get_thermal_conductivity_at_faces
   !  7) set_bcu
   !  8) set_bcv
   !  9) set_bcT
   ! 10) set_bcp
   ! 11) get_uin_vin_pin_Tin_Mw
   ! 12) get_plin_and_p_fictitious
   ! 13) get_u_v_at_fictitious_nodes_with_pl
   ! 14) get_boundary_simplec_coefficients
   ! 15) get_Uce_Vcn_at_boundary_faces
   ! 16) get_Uce_Vcn_at_boundary_faces_with_pl
   ! 16) get_isentropic_mass_flow
   ! 17) get_mach_area
   ! 18) get_initial_guess
   ! 19) get_boundary_nodes - ALTERADA PARA GEOMETRIAS GENÉRICAS, MAIORES QUE 2 PARTES
   ! 20) verificacao - Verifica o comando e chama calculo de 'a','b' para linear e 'cx' e 'cy' para arco
   ! 21) linear - Calcula coef. 'a' e 'b' da reta
   ! 22) arco - Calcula o 'cx' e 'cy', pontos de centro do arco de arredondamento
   ! 23) pontos - A cada x, essa rotina calcula o respectivo y, dado o comando
   !
   ! Last update: 06 Ago 2012
   !

   use coefficients

   implicit none

contains

   subroutine set_wall_temperature( nx, Twall ) ! Last one is output
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      real(8), dimension(nx), intent(out) :: Twall ! Wall temperature

      Twall = 300.d0

   end subroutine set_wall_temperature


   subroutine set_cp_and_gamma( nx, ny, Rg, cp, gcp) ! Last two are output
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      real(8), intent(in) :: Rg  ! Perfect gas constant
      real(8), dimension(nx*ny), intent(out) :: cp  ! Specific heat at const pressure
      real(8), dimension(nx*ny), intent(out) :: gcp ! gcp = gamma = Cp/Cv at center of CV P

      cp = 1004.5d0

      gcp = cp / ( cp - Rg )

   end subroutine set_cp_and_gamma


   subroutine set_laminar_viscosity_at_nodes( nx, ny, vlp) ! Last one is output
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(out) :: vlp ! Laminar viscosity at center of volume P

      vlp = 1.d-5 ! (Units Pa.s)

   end subroutine set_laminar_viscosity_at_nodes


   subroutine set_thermal_conductivity_at_nodes( nx, ny, kp) ! Last one is output
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(out) :: kp ! Thermal conductivity at center of volume P

      kp = 1.0d-3 ! (Units W/m.K)

   end subroutine set_thermal_conductivity_at_nodes


   ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
   ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
   subroutine get_laminar_viscosity_at_faces( nx, ny, vlp, vle, vln) ! Last two are output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(in)  :: vlp ! Laminar viscosity at center of volume P
      real(8), dimension (nx*ny), intent(out) :: vle ! Laminar viscosity at center of face east
      real(8), dimension (nx*ny), intent(out) :: vln ! Laminar viscosity at center of face north

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

      do j = 2, ny-1
         do i = 1, nx-1
            np   = nx * (j-1) + i
            npe  = np + 1

            vle(np) = 2.d0 * vlp(np) * vlp(npe) / ( vlp(np) + vlp(npe) )

         end do
      end do

      do i = 2, nx-1
         do j = 1, ny-1
            np  = nx * (j-1) + i
            npn = np + nx

            vln(np) = 2.d0 * vlp(np) * vlp(npn) / ( vlp(np) + vlp(npn) )

         end do
      end do

   end subroutine get_laminar_viscosity_at_faces


   ! CAUTION: In order to apply correctly this function, thermal conductivity kp must be known
   ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
   subroutine get_thermal_conductivity_at_faces( nx, ny, kp, ke, kn) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(in)  :: kp ! Thermal conductivity at center of volume P
      real(8), dimension (nx*ny), intent(out) :: ke ! Thermal conductivity at center of face east
      real(8), dimension (nx*ny), intent(out) :: kn ! Thermal conductivity at center of face north

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

      do j = 2, ny-1
         do i = 1, nx-1
            np   = nx * (j-1) + i
            npe  = np + 1

            ke(np) = 2.d0 * kp(np) * kp(npe) / ( kp(np) + kp(npe) )

         end do
      end do

      do i = 2, nx-1
         do j = 1, ny-1
            np  = nx * (j-1) + i
            npn = np + nx

            kn(np) = 2.d0 * kp(np) * kp(npn) / ( kp(np) + kp(npn) )

         end do
      end do

   end subroutine get_thermal_conductivity_at_faces

   subroutine set_bcu( nx, ny, modvis, xk, yk, u, v, au, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), dimension (nx*ny),   intent(in)  :: xk ! x_csi at center of face north
      real(8), dimension (nx*ny),   intent(in)  :: yk ! y_csi at center of face north
      real(8), dimension (nx*ny),   intent(in)  :: u  ! Cartesian velocity of the last iteraction
      real(8), dimension (nx*ny),   intent(in)  :: v  ! Cartesian velocity of the last iteraction
      real(8), dimension (nx*ny,9), intent(out) :: au ! Coefficients of the linear system for u
      real(8), dimension (nx*ny),   intent(out) :: b  ! Source vector of the linear system
      ! Auxiliary variables

      integer :: i, j, np, nps

      real(8) :: aux1, aux2, aux3, lambda


      ! West boundary (corners are not included)

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         au(np,:) =  0.d0

         au(np,6) = -1.d0

         au(np,5) =  1.d0

         b(np) = 0.d0

      end do

      ! East boundary (corners are not included)

      i = nx

      do j = 2, ny-1

         np   = nx * (j-1) + i

         au(np,:) =  0.d0

         au(np,4) = -1.d0

         au(np,5) =  1.d0

         b(np) = 0.d0

      end do

      ! South boundary ( SW and SE corners are included)

      j = 1

      do i = 1, nx

         np   = nx * (j-1) + i

         au(np,:) =  0.d0

         au(np,8) = -1.d0

         au(np,5) =  1.d0

         b(np) = 0.d0

      end do

      ! North boundary ( NW and NE corners are included)

      j = ny

      if ( modvis == 1 ) then ! Navier Stokes

         do i = 1, nx

            np   = nx * (j-1) + i

            au(np,:) = 0.d0

            au(np,2) = 1.d0

            au(np,5) = 1.d0

            b(np) = 0.d0

         end do

      else ! Euler

         do i = 2, nx-1

            np   = nx * (j-1) + i

            nps  = np - nx

            au(np,:) = 0.d0

            au(np,2) = 1.d0

            au(np,5) = 1.d0

            aux1 = sign(1.d0, u(nps) * xk(nps) + v(nps) * yk(nps) )

            aux2 =  u(nps) ** 2 +  v(nps) ** 2

            aux3 = xk(nps) ** 2 + yk(nps) ** 2

            lambda = aux1 * sqrt( aux2 / aux3 )

            b(np) = 2.d0 * lambda * xk(nps)

         end do

         ! NW corner

         i = 1

         np   = nx * (j-1) + i

         au(np,:) =  0.d0

         au(np,5) =  1.d0

         b(np) = 0.d0

         ! NE corner

         i = nx

         np   = nx * (j-1) + i

         au(np,:) =  0.d0

         au(np,5) =  1.d0

         b(np) = 0.d0

      end if

   end subroutine set_bcu

   subroutine set_bcv( nx, ny, modvis, xk, yk, u, v, av, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), dimension (nx*ny),   intent(in)  :: xk ! x_csi at center of face north
      real(8), dimension (nx*ny),   intent(in)  :: yk ! y_csi at center of face north
      real(8), dimension (nx*ny),   intent(in)  :: u  ! Cartesian velocity of the last iteraction
      real(8), dimension (nx*ny),   intent(in)  :: v  ! Cartesian velocity of the last iteraction
      real(8), dimension (nx*ny,9), intent(out) :: av ! Coefficients of the linear system for v
      real(8), dimension (nx*ny),   intent(out) :: b  ! Source vector of the linear system

      ! Auxiliary variables

      integer :: i, j, np, nps

      real(8) :: aux1, aux2, aux3, lambda

      ! West boundary (corners are not included)

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         av(np,:) = 0.d0

         av(np,6) = 1.d0

         av(np,5) = 1.d0

         b(np) = 0.d0

      end do

      ! East boundary (corners are not included)

      i = nx

      do j = 2, ny-1

         np   = nx * (j-1) + i

         av(np,:) =  0.d0

         av(np,4) = -1.d0

         av(np,5) =  1.d0

         b(np) = 0.d0

      end do

      ! South boundary ( SW and SE corners are included)

      j = 1

      do i = 1, nx

         np   = nx * (j-1) + i

         av(np,:) = 0.d0

         av(np,8) = 1.d0

         av(np,5) = 1.d0

         b(np) = 0.d0

      end do


      ! North boundary ( NW and NE corners are included)

      j = ny

      if ( modvis == 1 ) then ! Navier-Stokes

         do i = 1, nx

            np   = nx * (j-1) + i

            av(np,:) = 0.d0

            av(np,2) = 1.d0

            av(np,5) = 1.d0

            b(np) = 0.d0

         end do

      else ! Euler

         do i = 2, nx-1

            np   = nx * (j-1) + i

            nps  = np - nx

            av(np,:) =  0.d0

            av(np,2) =  1.d0

            av(np,5) =  1.d0

            aux1 = sign(1.d0, u(nps) * xk(nps) + v(nps) * yk(nps) )

            aux2 =  u(nps) ** 2 +  v(nps) ** 2

            aux3 = xk(nps) ** 2 + yk(nps) ** 2

            lambda = aux1 * sqrt( aux2 / aux3 )

            b(np) = 2.d0 * lambda * yk(nps)

         end do

         ! NW corner

         i = 1

         np   = nx * (j-1) + i

         av(np,:) =  0.d0

         av(np,5) =  1.d0

         b(np) = 0.d0

         ! NE corner

         i = nx

         np   = nx * (j-1) + i

         av(np,:) =  0.d0

         av(np,5) =  1.d0

         b(np) = 0.d0

      end if

   end subroutine set_bcv

   ! CAUTION: Twall must be difined in the fictitious volumes too
   subroutine set_bcT( nx, ny, ccTw, Tin, Twall, T, at, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      integer, intent(in) :: ccTw   ! ccTw = 0 -> adiabatic;  ccTw = 1 -> prescribed temperature
      real(8), dimension(ny),      intent(in) ::  Tin   ! Temperature in the entrance
      real(8), dimension(nx),      intent(in) ::  Twall ! Wall temperature
      real(8), dimension(nx*ny),   intent(in) ::  T     ! Temperature of the last iteraction
      real(8), dimension(nx*ny,9), intent(out) :: aT    ! Coefficients of the linear system
      real(8), dimension(nx*ny),   intent(out) :: b     ! Source vector of the linear system

      ! Auxiliary variables
      integer :: i, j, np

      ! West boundary
      i = 1
      do j = 2, ny-1
         np   = nx * (j-1) + i

         at(np,:) = 0.d0
         at(np,6) = 1.d0
         at(np,5) = 1.d0

         b(np) = 2.d0 * Tin(j)
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np   = nx * (j-1) + i

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

         b(np) = 0.d0
      end do

      ! South boundary ( SW and SE corners are included)
      j = 1
      do i = 1, nx
         np = nx * (j-1) + i

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

         b(np) = 0.d0
      end do

      ! North boundary ( NW and NE corners are included)
      j = ny
      if ( ccTw == 0 ) then ! Adiabatic
         do i = 1, nx
            np = nx * (j-1) + i

            at(np,:) =  0.d0
            at(np,2) = -1.d0
            at(np,5) =  1.d0

            b(np) = 0.d0
         end do
      else ! Prescribed temperature
         do i = 1, nx
            np = nx * (j-1) + i

            at(np,:) = 0.d0
            at(np,2) = 1.d0
            at(np,5) = 1.d0

            b(np) = 2.d0 * Twall(i)
         end do
      end if

   end subroutine set_bcT

   subroutine set_bcp( nx, ny, ap, b) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension(nx*ny,5), intent(out) :: ap    ! Coefficients of the linear system
      real(8), dimension(nx*ny),   intent(out) :: b     ! Source vector of the linear system

      ! Auxliary variables
      integer :: i, j, np
      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx * (j-1) + i

         ap(np,:) = 0.d0
         ap(np,4) = 1.d0
         ap(np,3) = 1.d0

         b(np) = 0.d0
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np   = nx * (j-1) + i

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

         b(np) = 0.d0
      end do

      ! South boundary
      j = 1
      do i = 2, nx-1
         np  = nx * (j-1) + i

         ap(np,:) =  0.d0
         ap(np,5) = -1.d0
         ap(np,3) =  1.d0

         b(np) = 0.d0
      end do

      ! North boundary
      j = ny
      do i = 2, nx-1
         np  = nx * (j-1) + i

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

         b(np) = 0.d0
      end do
   end subroutine set_bcp

   subroutine get_uin_vin_pin_Tin_Mw( nx, ny, gamma, Rg, po, T0, u & ! Input
      ,                             uin, vin, pin, Tin, Mw)        ! Output
      implicit none
      integer, intent(in) :: nx    ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny    ! Number of volumes in eta direction (real+fictitious)
      real(8), intent(in) :: gamma ! gamma = Cpo / Cvo in the chamber
      real(8), intent(in) :: Rg    ! Perfect gas constant
      real(8), intent(in) :: po    ! Stagnation pressure in the chamber
      real(8), intent(in) :: T0    ! Stagnation temperature in the chamber
      real(8), dimension (nx*ny), intent(in)  :: u   ! Cartesian velocity of the last iteraction
      real(8), dimension (ny),    intent(out) :: uin ! Velocity u in the entrance
      real(8), dimension (ny),    intent(out) :: vin ! Velocity v in the entrance
      real(8), dimension (ny),    intent(out) :: pin ! Pressure in the entrance
      real(8), dimension (ny),    intent(out) :: Tin ! Temperature in the entrance
      real(8), dimension (ny),    intent(out) :: Mw  ! Mach number in the entrance

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

      i = 1
      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         uin(j) = ( u(np) + u(npe) ) / 2.d0

         vin(j) = 0.d0 ! From boundary conditions

         Tin(j) = T0 - (gamma-1.d0) / ( 2.d0 * gamma * Rg ) * ( uin(j)**2 + vin(j)**2 )

         Mw(j) = dsqrt( ( uin(j)**2 + vin(j)**2 ) / ( gamma * Rg * Tin(j) ) )

         pin(j) = po * (1.d0 + 0.5d0 * (gamma-1.d0) * Mw(j)**2 ) ** ( -gamma / (gamma-1.d0) )

      end do

   end subroutine get_uin_vin_pin_Tin_Mw

   subroutine get_plin_and_p_fictitious( nx, ny, pina, pin, plin, p) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx    ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny    ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension(ny),    intent(in)  :: pina ! Pressure in the entrance at a time step before
      real(8), dimension(ny),    intent(in)  :: pin  ! Pressure in the entrance
      real(8), dimension(ny),    intent(out) :: plin ! Pressure correction in the entrance
      real(8), dimension(nx*ny), intent(out) :: p    ! Pressure at center of volumes
      ! Auxiliary variables
      integer :: i, j, np, npe

      ! West boundary

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         p(np) = 2.d0 * pin(j) - p(npe)

         plin(j) = pin(j) - pina(j)

      end do
   end subroutine get_plin_and_p_fictitious


   subroutine get_u_v_at_fictitious_nodes_with_pl( nx, ny, modvis, xk, yk, u, v) ! InOutput: last two entries
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), dimension (nx*ny), intent(in)    :: xk ! x_csi at center of face north
      real(8), dimension (nx*ny), intent(in)    :: yk ! y_csi at center of face north
      real(8), dimension (nx*ny), intent(inout) :: u  ! Cartesian velocity
      real(8), dimension (nx*ny), intent(inout) :: v  ! Cartesian velocity

      ! Auxiliary variables

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

      real(8) :: aux1, aux2, aux3, lambda

      ! West boundary

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         u(np) = u(npe)

         v(np) = -v(npe)

      end do


      ! East boundary

      i = nx

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npw  = np - 1

         u(np) = u(npw)

         v(np) = v(npw)

      end do


      ! South boundary (Corners SW and SE are included)

      j = 1

      do i = 1, nx

         np   = nx * (j-1) + i

         npn  = np + nx

         u(np) =  u(npn)

         v(np) = -v(npn)

      end do


      ! North boundary (Corners NW and NE are included)

      j = ny

      if ( modvis == 1 ) then ! Viscous

         do i = 1, nx

            np   = nx * (j-1) + i

            nps  = np - nx

            u(np) = -u(nps)

            v(np) = -v(nps)

         end do

      else ! Inviscid

         do i = 2, nx-1

            np   = nx * (j-1) + i

            nps  = np - nx

            aux1 = sign(1.d0, u(nps) * xk(nps) + v(nps) * yk(nps) )

            aux2 =  u(nps) ** 2 +  v(nps) ** 2

            aux3 = xk(nps) ** 2 + yk(nps) ** 2

            lambda = aux1 * sqrt( aux2 / aux3 )


            u(np) = -u(nps) + lambda * xk(nps) + lambda * xk(nps)

            v(np) = -v(nps) + lambda * yk(nps) + lambda * yk(nps)

         end do

      end if
   end subroutine get_u_v_at_fictitious_nodes_with_pl

   subroutine get_boundary_simplec_coefficients( nx, ny, de, dn) ! InOutput: de, dn
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(inout) :: de  ! SIMPLEC coefficient for Uce
      real(8), dimension (nx*ny), intent(inout) :: dn  ! SIMPLEC coefficient for Vcn

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

      ! West boundary

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         de(np) = de(npe)

      end do

      ! East boundary

      i = nx-1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npw  = np - 1

         de(np) = de(npw)

      end do

      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         dn(np) = 0.d0

      end do

      ! North boundary

      j = ny-1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         dn(np) = 0.d0

      end do

   end subroutine get_boundary_simplec_coefficients

   subroutine get_velocities_at_boundary_faces( nx, ny, ye, u, v, ue, ve, un, vn, Uce, Vcn) ! Output: last six
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (nx*ny), intent(in)  :: ye  ! y_eta at face east of volume P
      real(8), dimension (nx*ny), intent(in)  :: u   ! Cartesian velocity of the last iteraction
      real(8), dimension (nx*ny), intent(in)  :: v   ! Cartesian velocity of the last iteraction
      real(8), dimension (nx*ny), intent(out) :: ue  ! Cartesian velocity u at center of east face
      real(8), dimension (nx*ny), intent(out) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension (nx*ny), intent(out) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension (nx*ny), intent(out) :: vn  ! Cartesian velocity v at center of north face
      real(8), dimension (nx*ny), intent(out) :: Uce ! Contravariant velocity U at east face
      real(8), dimension (nx*ny), intent(out) :: Vcn ! Contravariant velocity V at north face

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

      ! West boundary

      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) = ( v(np) + v(npe) ) / 2.d0

         Uce(np) = ue(np) * ye(np)

      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)

      end do

      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         npn  = np + nx

         un(np) = ( u(np) + u(npn) ) / 2.d0

         vn(np) = ( v(np) + v(npn) ) / 2.d0

         Vcn(np) = 0.d0

      end do

      ! North boundary

      j = ny-1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         npn  = np + nx

         un(np) = ( u(np) + u(npn) ) / 2.d0

         vn(np) = ( v(np) + v(npn) ) / 2.d0

         Vcn(np) = 0.d0

      end do

   end subroutine get_velocities_at_boundary_faces

   subroutine get_Uce_Vcn_at_boundary_faces_with_pl( nx, ny, pl, de, Uce) ! InOutput: Last entry
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(in)    :: pl  ! Pressure correction
      real(8), dimension(nx*ny), intent(in)    :: de  ! SIMPLEC coefficient for Uce
      real(8), dimension(nx*ny), intent(inout) :: Uce ! Contravariant velocity U at east face

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

      ! West boundary

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         Uce(np) = Uce(np) + de(np) * ( pl(np) - pl(npe) )

      end do

      ! West boundary

      i = nx-1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         Uce(np) = Uce(np) + de(np) * ( pl(np) - pl(npe) )

      end do

      ! North and South Vcn do not neet to be changed

   end subroutine get_Uce_Vcn_at_boundary_faces_with_pl


   subroutine get_isentropic_mass_flow( po, T0, gamma, Rg, Sg , fm1D)
      implicit none
      real(8), intent(in)  :: po    ! Stagnation pressure
      real(8), intent(in)  :: T0    ! Stagnation temperature
      real(8), intent(in)  :: gamma ! Specific heat ratio
      real(8), intent(in)  :: Rg    ! Perfect gas constant
      real(8), intent(in)  :: Sg    ! Throttle area
      real(8), intent(out) :: fm1D  ! Mass flow rate

      fm1D = po * Sg * sqrt( gamma/ (Rg*T0) * (2.d0/(gamma+1.d0)) ** ((gamma+1.d0)/(gamma-1.d0)) )

   end subroutine get_isentropic_mass_flow


   subroutine get_mach_area( kf, ar, gamma, M )
      implicit none
      integer, intent(in)  :: kf    ! Kind of flow: 0 = subsonic, 1 = supersonic
      real(8), intent(in)  :: ar    ! Area ratio = local area / throttle area
      real(8), intent(in)  :: gamma ! Specific heat ratio
      real(8), intent(out) :: M     ! Mach number at the local area

      ! Auxiliary variables

      integer :: it

      real(8) :: M0, M1, f, x

      f(x) = ( 2.d0 * ( 1.d0 + (gamma-1.d0) * x**2 / 2.d0 ) / (gamma+1.d0) ) &
         ** ((gamma+1.d0)/(gamma-1.d0)) - ( x * ar ) ** 2

      if ( abs(ar-1.d0) < 1.d-15 ) then ! Sonic flow in the throttle

         M = 1.d0

      else

         if ( kf == 0 ) then ! Subsonic
            M0 = 1.d-4
            M1 = 1.d0
         else                ! Supersonic
            M0 = 1.d0
            M1 = 100.d0
         end if

         do it = 1, 1000

            M = 0.5d0 * ( M0 + M1 )

            if ( f(M) * f(M0) > 0 ) then

               M0 = M

            else

               M1 = M

            end if

            if ( abs(M1-M0) < 1.d-15 ) exit

         end do

      end if

   end subroutine get_mach_area


   subroutine get_initial_guess( nx, ny, ig, modvis, beta, po, T0, gamma  & ! Input
      ,                        Rg, Sg, ye, xk, yk, radius, rn, x, xp    & ! Input
      ,                        M1D, p1D, T1D, u1D, p, T, u, v, ue       & ! Output
      ,                        un, ve, vn, Uce, Vcn, uin, vin, pin, Tin, Mw & ! Output
      ,                        fm1D, Fd1D, Fpv1D, de, dn, ro, roe, ron  & ! Output
      ,                        a, ap, b ) ! Output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      integer, intent(in) :: ig
      integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
      real(8), intent(in) :: beta   ! Constant of the UDS/CDS mixing scheme
      real(8), intent(in) :: po     ! Stagnation pressure in the chamber
      real(8), intent(in) :: T0     ! Stagnation temperature in the chamber
      real(8), intent(in) :: gamma  ! gamma = Cpo / Cvo in the chamber
      real(8), intent(in) :: Rg     ! Perfect gas constant
      real(8), intent(in) :: Sg     ! Throttle area
      real(8), dimension(nx*ny), intent(in) :: ye     ! y_eta at face east of volume P
      real(8), dimension(nx*ny), intent(in) :: xk     ! x_csi at center of face north
      real(8), dimension(nx*ny), intent(in) :: yk     ! y_csi at face north of volume P
      real(8), dimension(nx*ny), intent(in) :: radius ! Radius of northest corner of volume P
      real(8), dimension(nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension(nx*ny), intent(in) :: x      ! Coord. x of the northest corner of volume P
      real(8), dimension(nx*ny), intent(in) :: xp     ! Coord. x of the centroid of volume P

      real(8), intent(out) :: fm1D  !
      real(8), intent(out) :: Fd1D  !
      real(8), intent(out) :: Fpv1D !

      real(8), dimension(nx),    intent(out) :: M1D !  Mach number. Isentropic flow
      real(8), dimension(nx),    intent(out) :: p1D !  Pressure. Isentropic flow
      real(8), dimension(nx),    intent(out) :: T1D !  Temperature. Isentropic flow
      real(8), dimension(nx),    intent(out) :: u1D !  Velocity. Isentropic flow

      real(8), dimension(nx*ny), intent(out) :: p   ! Pressure at center o volume P
      real(8), dimension(nx*ny), intent(out) :: T   ! Temperature of the last iteraction
      real(8), dimension(nx*ny), intent(out) :: u   ! Cartesian velocity of the last iteraction
      real(8), dimension(nx*ny), intent(out) :: v   ! Cartesian velocity of the last iteraction

      real(8), dimension(nx*ny), intent(out) :: ue  ! Cartesian velocity u at center of east face
      real(8), dimension(nx*ny), intent(out) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension(nx*ny), intent(out) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension(nx*ny), intent(out) :: vn  ! Cartesian velocity v at center of north face

      real(8), dimension(nx*ny), intent(out) :: Uce ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(out) :: Vcn ! Contravariant velocity V at north face

      real(8), dimension (ny),   intent(out) :: uin ! Velocity u in the entrance
      real(8), dimension (ny),   intent(out) :: vin ! Velocity v in the entrance
      real(8), dimension (ny),   intent(out) :: pin ! Pressure in the entrance
      real(8), dimension (ny),   intent(out) :: Tin ! Temperature in the entrance
      real(8), dimension (ny),   intent(out) :: Mw  ! Mach number in the entrance

      real(8), dimension(nx*ny), intent(out) :: de  ! Simplec coef. for the contravariant velocity U (east face)
      real(8), dimension(nx*ny), intent(out) :: dn  ! Simplec coef. for the contravariant velocity V (north face)

      real(8), dimension(nx*ny), intent(out) :: ro  ! Specific mass (absolute density) at center of volumes
      real(8), dimension(nx*ny), intent(out) :: roe ! Absolute density at east face
      real(8), dimension(nx*ny), intent(out) :: ron ! Absolute density at north face

      real(8), dimension(nx*ny,9), intent(out) :: a  ! Coefficients of the linear system for u, v, T
      real(8), dimension(nx*ny,5), intent(out) :: ap ! Coefficients of the linear system
      real(8), dimension(nx*ny),   intent(out) :: b  ! Source vector of the linear system


      ! Auxiliary variables
      real(8), parameter :: pi = dacos(-1.d0)

      integer :: i, j, np, npe, npw, kf
      real(8) :: Si, Me1D, pe1D, Te1D, ue1D, aux

      call get_isentropic_mass_flow( po, T0, gamma, Rg, Sg , fm1D)


      ! Calculating isentropic solution

      do i = 2, nx-1

         j = ny-1

         np = nx * (j-1) + i

         Si = pi * rn(np) ** 2

         kf = 0

         if ( i > ig ) kf = 1

         call get_mach_area( kf, Si/Sg, gamma, M1D(i) )

         aux = 1.d0 + 0.5d0 * (gamma-1.d0) * M1D(i)**2

         p1D(i) = po * aux ** (-gamma/(gamma-1.d0))

         T1D(i) = T0 / aux

         u1D(i) = M1D(i) * sqrt( gamma * Rg * T1D(i) )

         do j = 1, ny

            np = nx * (j-1) + i

            p(np) = p1D(i)

            T(np) = T1D(i)

            u(np) = u1D(i)

            if ( j > 1 .and. j < ny-1 ) then

               un(np) = u1D(i)

               Vcn(np) = -un(np) * yk(np) ! vn = 0

            end if

         end do

      end do

      ! Calculates u and v at fictitious nodes using boundary conditions
      call get_u_v_at_fictitious_nodes_with_pl( nx, ny, modvis, xk, xk, u, v) ! InOutput: last two entries

      ! Calculates gas proprieties in the entrance
      call get_uin_vin_pin_Tin_Mw( nx, ny, gamma, Rg, po, T0, u & ! Input
      ,                             uin, vin, pin, Tin, Mw)  ! Output

      ! Calculates p and T in the fictitious nodes of the entrance according to boundary conditions
      i = 1
      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         p(np) = 2.d0 * pin(j) - p(npe)

         T(np) = 2.d0 * Tin(j) - T(npe)

      end do

      ! Calculates p and T in the fictitious nodes of the exit according to boundary conditions
      i = nx
      do j = 2, ny-1

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

         p(np) = p(npw)
         T(np) = T(npw)

      end do

      ! Initial guess of velocities in east face of each CV
      ! and analytical solution in the entrance and exit boundaries

      do i = 1, nx-1

         j = ny-1

         np = nx * (j-1) + i

         Si = pi * radius(np) ** 2

         kf = 0

         if ( i > ig ) kf = 1

         call get_mach_area( kf, Si/Sg, gamma, Me1D )

         aux = 1.d0 + 0.5d0 * (gamma-1.d0) * Me1D**2

         pe1D = po * aux ** (-gamma/(gamma-1.d0))

         Te1D = T0 / aux

         ue1D = Me1D * sqrt( gamma * Rg * Te1D )

         do j = 2, ny-1

            np = nx * (j-1) + i

            ue(np) = ue1D

            Uce(np) = ue1D * ye(np) ! v = 0

         end do

         ! Entrance boundary
         if ( i == 1 ) then
            u1D(i) = ue1D
            p1D(i) = pe1D
            T1D(i) = Te1D
            M1D(i) = Me1D
         end if

         ! Exit boundary
         if ( i == nx-1 ) then
            u1D(nx) = ue1D
            p1D(nx) = pe1D
            T1D(nx) = Te1D
            M1D(nx) = Me1D

            Fd1D = fm1D * ue1D
            Fpv1D = pe1D * Si
         end if

      end do

      ! Velocities at boundary faces
      call get_velocities_at_boundary_faces( nx, ny, ye, u, v, ue, ve, un, vn, Uce, Vcn) ! Output: last six

      ! Simplec coefficients at boundary faces
      call get_boundary_simplec_coefficients( nx, ny, de, dn) ! InOutput: de, dn

      ! Specific mass at nodes
      call get_density_at_nodes( nx, ny, Rg, p, T, ro) ! ro is output

      ! Specific mass at faces
      call get_density_at_faces( nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

      a  = 0.d0
      ap = 0.d0
      b  = 0.d0

      ap(:,3) = 1.d0

   end subroutine get_initial_guess

   subroutine get_boundary_nodes( nx, ny, ig, Sg, rcg) ! Output: last three entries ! Modificada P/ ALTERAR O GRID POR ARQUIVO DE SAIDA
      implicit none
      integer, intent(in)  :: nx
      integer, intent(in)  :: ny
      integer, intent(out) :: ig  ! i=ig when the throttle cross section is in the east face of the CV
      real(8), intent(out) :: Sg  ! area of the throttle
      real(8), intent(out) :: rcg ! Throttle curvature radius

      real(8), allocatable, dimension(:) ::  Lp, Dx ! Conjunto dos comprimentos de cada uma das partes da tubeira e distância de cada parte
      integer, allocatable, dimension(:) ::  Npvc ! Conjunto do número de volumes de controle que cada parte da tubeira ocupa

      integer :: NXdv ! Variavel de refino da malha, se 1=malha-base, 2 = 2*malha-base, etc.
      integer :: volumes ! Número de volumes de controle na direção axial
      integer :: partes ! Número de partes que a tubeira foi repartida

      real(8) :: L  ! comprimento total
      real(8) :: Ri ! raio na entrada da tubeira (m)
      real(8) :: rth
      real(8) :: a3,b3


      real(8) :: x(nx), y
      real(8),parameter :: pi = dacos(-1.d0)

      integer :: i, j, soma
      integer :: sistema, localg

      character*128 :: comando
      real(8) :: a, b, x1,y1,x2,y2
      real(8) :: r=0.d0


      ! ----- Leitura dos dados da tubeira pelo arquivo contorno.txt
      open(40,file='./mach2d-input/contorno.txt');
      ! Sistema em metros ou em polegadas
      read(40,*) sistema ; if (sistema == 1)  open(13,file='./mach2d-input/contorno_espelho_metrico.txt')
      if(sistema==1) write(13,*) sistema, " :: Sistema de Unidades (0 = metros) (1 = polegadas)"

      ! Raio de Entrada da Tubeira
      read(40,*) Ri   ; if(sistema==1) write(13,'(1PE15.5e3,A)') Ri*(25.4d0/1000.d0) , " :: Raio de entrada da tubeira"

      ! Raio da garganta
      read(40,*) rth     ; if(sistema==1) write(13,'(1PE15.5e3,A)') rth*(25.4d0/1000.d0) , " :: Raio da garganta"

      ! Raio de arredondamento da garganta
      read(40,*) rcg     ; if(sistema==1) write(13,'(1PE15.5e3,A)') rcg*(25.4d0/1000.d0) , " :: Raio de arredondamento da garganta"

      ! O número do volume de controle que está a leste da garganta
      read(40,*) localg  ; if(sistema==1) write(13,'(I5,A)') localg  ,&
         " :: Localização da garganta em relação ao volume de controle (ela está a leste dele)"

      ! partes que a tubeira foi repartida
      read(40,*) partes  ; if(sistema==1) write(13,'(I4,A)') partes  , " :: Número de partes em que a tubeira foi repartida"

      allocate(Lp(partes), Dx(partes), Npvc(partes))
      do i=1, 8
         read(40,*)
      enddo

      if(sistema ==1) then
         write(13,'(A)') " :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::"
         write(13,'(A,A)') " :: Põe-se o comprimento de cada parte da tubeira e a quantidade resp", &
            "ectiva de volumes de controle que ela ocupa"
         write(13,'(A)') " :: Após isto o contorno propriamente dito, comandos:"
         write(13,'(A)') " :: 1) Line : Interpolação Linear, Argumentos: Yfinal"
         write(13,'(A,A)') " :: 2) arcr : Interpolação Circular, Argumentos: Yfinal, raio de arredon", &
            "damento (+ concavidade para baixo) (- concavidade para cima)"
         write(13,'(A)') " :: 3) arcc+ : Interpolação Circular, Argumentos: Centro em x e em y , Concavidade para baixo"
         write(13,'(A)') " :: 4) arcc- : Interpolação Circular, Argumentos: Centro em x e em y , Concavidade para cima"
         write(13,'(A)') " :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::"
      endif
      ! Encontrar o comprimento total da tubeira e o numero de volumes
      L = 0.d0
      volumes = 0
      do i=1, partes
         read(40,*) Lp(i), Npvc(i) ; if(sistema==1) write(13,*) Lp(i)*25.4d0/1000d0, Npvc(i)
         L = L + Lp(i)
         volumes = volumes + Npvc(i)
      enddo

      NXdv = (nx-2) / volumes
      ig = NXdv*localg+1 ! Localizacao da garganta

      ! Encontrar o comprimento dos volumes de controle em cada parte
      do i=1, partes
         Dx(i) = Lp(i)/(NXdv*Npvc(i))
      enddo

      ! ---- Adequação do Sistema de Unidades (se estiver em polegadas, transforma em metros)
      if (sistema==1) then
         Ri=Ri*(25.4d0/1000d0)
         rth=rth*(25.4d0/1000d0)
         rcg=rcg*(25.4d0/1000d0)
         L = L*(25.4d0/1000d0)
         do i=1, partes
            Lp(i) = Lp(i)*(25.4d0/1000d0)
            Dx(i) = Dx(i)*(25.4d0/1000d0)
         enddo
      endif

      Sg = pi*rth**2  ! menor área da garganta
      open(10,file='./mach2d-input/gridboundary.dat')


      ! *** CONTORNO SUL ***
      write(10,*) "south boundary"

      I = 1
      x(i) = 0.d0
      y = 0.d0
      write(10,*) x(i), y
      ! Parte(1)
      DO I = 2, NXdv*Npvc(1)+1
         x(i) = x(i-1) + dx(1)
         write(10,*) x(i), y
      end do
      soma = Npvc(1)

      ! Parte(2) a (partes-1)
      do j=2, partes-1
         DO I = NXdv*soma+2, NXdv*(soma+Npvc(j))+1
            x(i) = x(i-1) + dx(j)
            write(10,*) x(i), y
         end do
         soma = soma + Npvc(j)
      end do

      ! Parte(partes)
      DO I = NXdv*soma+2, NXdv*(soma+Npvc(partes))
         x(i) = x(i-1) + dx(partes)
         write(10,*) x(i), y
      end do

      I = NX-1
      X(i) = L
      Y = 0.0d0
      write(10,*) x(i), y

      ! *** CONTORNO NORTE ***
      open(61,file='./mach2d-input/geometria_completa.txt'); write(61,*) "sistema = ", sistema
      write(10,*) "north boundary"

      x1 = 0
      x2 = Lp(1)
      y1 = Ri

      read(40,*) comando

      if(sistema==1) write(13,'(A)') comando
      if(comando == "line" .or. comando == "Line") then
         read(40,*) y2 ; if(sistema==1) write(13,*) y2*25.4d0/1000.d0
      endif
      if(comando == "arcr" .or. comando == "Arcr") then
         read(40,*) y2, r ; if(sistema==1) write(13,*) y2*25.4d0/1000.d0, r*25.4d0/1000.d0
      endif
      if(comando == "arcc-" .or. comando == "Arcc-".or.comando == "arcc+" .or. comando == "Arcc+") then
         read(40,*) a, b ; if(sistema==1) write(13,*) a*25.4d0/1000.d0, b*25.4d0/1000.d0
      endif

      call verificacao(sistema,comando,x1,y1,x2,y2,r,a,b)
      ! Parte(1)

      DO I = 1, NXdv*Npvc(1)
         call pontos(comando,x(i),a,b,r,y)
         write(10,*) x(i), y
      end do

      soma = Npvc(1)

      ! Parte(2) a (partes-1)
      do j=2, partes-1
         read(40,*) comando
         y1=y2
         x1 = x2
         x2 = x2+Lp(j)

         if(sistema==1) write(13,'(A)') comando
         if(comando == "line" .or. comando == "Line") then
            read(40,*) y2 ; if(sistema==1) write(13,*) y2*25.4d0/1000.d0
         endif
         if(comando == "arcr" .or. comando == "Arcr") then
            read(40,*) y2, r ; if(sistema==1) write(13,*) y2*25.4d0/1000.d0, r*25.4d0/1000.d0
         endif
         if(comando == "arcc-" .or. comando == "Arcc-".or.comando == "arcc+" .or. comando == "Arcc+") then
            read(40,*) a, b ; if(sistema==1) write(13,*) a*25.4d0/1000.d0, b*25.4d0/1000.d0
         endif

         call verificacao(sistema,comando,x1,y1,x2,y2,r,a,b)

         do i=NXdv*soma+1,NXdv*(soma+Npvc(j))
            call pontos(comando,x(i),a,b,r,y)
            write(10,*) x(i), y
         enddo
         soma = soma + Npvc(j)
      enddo

      ! Parte(partes)

      read(40,*) comando
      y1=y2
      x1 = x2

      x2 = x2 + Lp(partes)

      if(sistema==1) write(13,'(A)') comando
      if(comando == "line" .or. comando == "Line") then
         read(40,*) y2 ; if(sistema==1) write(13,*) y2*25.4d0/1000.d0
      endif
      if(comando == "arcr" .or. comando == "Arcr") then
         read(40,*) y2, r ; if(sistema==1) write(13,*) y2*25.4d0/1000.d0, r*25.4d0/1000.d0
      endif
      if(comando == "arcc-" .or. comando == "Arcc-".or.comando == "arcc+" .or. comando == "Arcc+") then
         read(40,*) a, b ; if(sistema==1) write(13,*) a*25.4d0/1000.d0, b*25.4d0/1000.d0
      endif

      call verificacao(sistema,comando,x1,y1,x2,y2,r,a,b)
      soma = soma + Npvc(partes)
      Do i=NXdv*(soma-Npvc(partes))+1,NXdv*soma+1
         call pontos(comando,x(i),a,b,r,y)
         write(10,*) x(i), y
      enddo

      I = NX
      X(i) = L
      Y = 0.0d0
      write(10,*) x(i), y

      close(40) ! Arquivo da Geometria
      close(10) ! Arquivo para entrada do contorno no mach2d
      close(61) ! Arquivo de saida da geometria para comparações
      if(sistema==1)  close(13) ! Arquivo espelho da Geometria em metros (se estiver em polegadas)

      deallocate(Lp, Dx, Npvc)
   end subroutine get_boundary_nodes

   !------------------------------------
   subroutine verificacao(sistema,comando,x1,y1,x2,y2,r,a,b)
      integer, intent(in) :: sistema
      character*128, intent(inout) :: comando
      real(8), intent(inout) :: x1,y1,x2,y2,r
      real(8), intent(inout)  :: a, b
      logical :: linha=.false.
      logical :: arcr=.false.
      logical :: arcc=.false.

      if(sistema == 1) then
         y2=y2*(25.4d0/1000d0)
         r = r*(25.4d0/1000d0)
         a = a*(25.4d0/1000d0)
         b = b*(25.4d0/1000d0)
      endif

      if (comando == "line" .or. comando == "Line") linha = .true.
      if (comando == "arcr" .or. comando == "Arcr") arcr = .true.
      if (comando == "arcc-" .or. comando == "Arcc-".or.comando == "arcc+" .or. comando == "Arcc+") arcc = .true.

      if (linha) call linear(x1,y1,x2,y2,a,b)

      if (arcr .and. (r**2-((dsqrt((x2-x1)**2 + (y2-y1)**2))/2)**2) >= 0) call arco(x1,y1,x2,y2,r,a,b)
      if (arcr .and. (r**2-((dsqrt((x2-x1)**2 + (y2-y1)**2))/2)**2) < 0) then
         write(*,*) "ERRO - ARCO DE ARREDONDAMENTO INVALIDO !!!! - PASSANDO PARA RETA ..."
         comando = "line"
         call linear(x1,y1,x2,y2,a,b)
      endif
      if(arcc) then
         r = dsqrt((x1-a)**2+(y1-b)**2)

         ! Altera a concavidade
         if(comando == "arcc-" .or. comando == "Arcc-") r=-r

         y2 = dsin(dacos((x2-a)/r))*r+b
      endif

      arcc=.false.
      arcr=.false.
      linha=.false.

      if(sistema == 1) then
         write(61,'(/A,A)') "comando = ", trim(adjustl(comando))
         write(61,*) "x1 = ",x1*1000d0/25.4d0
         write(61,*) "x2 = ",x2*1000d0/25.4d0
         write(61,*) "y1 = ",y1*1000d0/25.4d0
         write(61,*) "y2 = ", y2*1000d0/25.4d0
         if(trim(adjustl(comando)) == 'line' .or. trim(adjustl(comando)) == 'Line') then
            write(61,*) "a = ",a*1000d0/25.4d0
            write(61,*) "b = ",b*1000d0/25.4d0
         else
            write(61,*) "r = ",r*1000d0/25.4d0
            write(61,*) "arc(cx) = ",a*1000d0/25.4d0
            write(61,*) "arc(cy) = ",b*1000d0/25.4d0
         endif
      else
         write(61,'(/A,A)') "comando = ", trim(adjustl(comando))
         write(61,*) "x1 = ",x1
         write(61,*) "x2 = ",x2
         write(61,*) "y1 = ",y1
         write(61,*) "y2 = ", y2
         write(61,*) "r = ",r
         if(trim(adjustl(comando)) == 'line' .or. trim(adjustl(comando)) == 'Line') then
            write(61,*) "a = ",a
            write(61,*) "b = ",b
         else
            write(61,*) "arc(cx) = ",a
            write(61,*) "arc(cy) = ",b
         endif
      endif

   end subroutine verificacao
   !---------------------------------
   subroutine linear(x1,y1,x2,y2,a,b) ! Encontra os valores de 'a' e 'b' para uma reta
      real(8), intent(inout) :: x1,y1,x2,y2
      real(8), intent(out) :: a, b
      a = (y1-y2)/(x1-x2)
      b = y1-a*x1
   end subroutine linear
   !------------------------------------
   subroutine arco(x1,y1,x2,y2,r,cx,cy) ! Encontra o centro dos arcos de concordância(a = CX e b = CY)
      real(8), intent(inout) :: x1, y1, x2, y2, r
      real(8), intent(out) :: cx, cy
      real(8) q, x3, y3

      x3 = (x1+x2)/2.d0
      y3 = (y1+y2)/2.d0
      q = dsqrt((x2-x1)**2 + (y2-y1)**2)

      if(r > 0) then
         cx = x3 + dsqrt(r**2-(q/2)**2)*(y2-y1)/q
         cy = y3 - dsqrt(r**2-(q/2)**2)*(x2-x1)/q
      endif
      if(r < 0) then
         cx = x3 - dsqrt(r**2-(q/2)**2)*(y2-y1)/q
         cy = y3 + dsqrt(r**2-(q/2)**2)*(x2-x1)/q
      endif
   end subroutine arco
   !------------------------------------
   subroutine pontos(comando,x,a,b,r,y)
      character*128, intent(in)  :: comando
      real(8), intent(in)  :: x, a, b,r
      real(8), intent(out)  :: y
      real(8)   :: cx, cy
      logical :: linha=.false.
      logical :: arcr=.false.
      logical :: arcc=.false.

      if (comando == "line" .or. comando == "Line") linha = .true.
      if (comando == "arcr" .or. comando == "Arcr") arcr = .true.
      if (comando == "arcc-" .or. comando == "Arcc-".or.comando == "arcc+" .or. comando == "Arcc+") arcc = .true.

      if(linha)  y = a*x+b

      if (arcr.or.arcc) then
         cx = a
         cy = b
         y = dsin(dacos((x-cx)/r))*r+cy
         !y= sqrt(r**2-(x-cx)**2)+cy
      endif
      arcc=.false.
      arcr=.false.
      linha=.false.
   end subroutine pontos
   !------------------------------------


   subroutine get_mass_conservation_check( nx,  ny,  dt,  rp,  re,  rn &
      ,                                     Jp, roe, ron,  ro, roa, Uce &
      ,                                    Vcn, RES) ! Last one is output
      implicit none
      integer, intent(in) :: nx     ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     ! Number of volumes in eta direction (real+fictitious)
      real(8), intent(in) :: dt     ! Time step
      real(8), dimension (nx*ny), intent(in) :: rp     ! Radius of the center of volume P
      real(8), dimension (nx*ny), intent(in) :: re     ! Radius of the center of east face of volume P
      real(8), dimension (nx*ny), intent(in) :: rn     ! Radius of the center of north face of volume P
      real(8), dimension (nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
      real(8), dimension (nx*ny), intent(in) :: roe    ! Absolute density at east face
      real(8), dimension (nx*ny), intent(in) :: ron    ! Absolute density at north face
      real(8), dimension (nx*ny), intent(in) :: ro     ! Absolute density at center of vol. P
      real(8), dimension (nx*ny), intent(in) :: roa    ! Absolute density at a time step before at center of vol. P
      real(8), dimension (nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face
      real(8), dimension (nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face
      real(8),                    intent(out) :: RES    ! Norm L_inf of the residual of the mass conservation


      ! Auxiliary variables

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

      real(8) :: aux, DM, Ms, Mn, Mw, Me




      RES = 0.d0

      do i = 2, nx-1

         do j = 2, ny-1

            np  = nx * (j-1) + i

            nps  = np - nx

            npw  = np - 1


            DM = ( ro(np) - roa(np) ) * rp(np) / ( Jp(np) * dt )

            Me = roe(np ) * re(np ) * Uce(np )

            Mw = roe(npw) * re(npw) * Uce(npw)

            Mn = ron(np ) * rn(np ) * Vcn(np )

            Ms = ron(nps) * rn(nps) * Vcn(nps)



            aux = abs( DM + Me - Mw + Mn - Ms )

            if ( aux > RES ) RES = aux

         end do

      end do

   end subroutine get_mass_conservation_check




   subroutine get_p_and_ro_extrapolation_to_fictitious(nx, ny, pin, g, p, ro) ! InOutput: p, Output: ro
      implicit none
      integer, intent(in) :: nx    ! Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny    ! Number of volumes in eta direction (real+fictitious)
      real(8), dimension (ny),       intent(in) :: pin ! Pressure in the entrance
      real(8), dimension (nx*ny),    intent(in) :: g   ! ro/p = 1/(Rg T) for ideal gases
      real(8), dimension (nx*ny), intent(inout) :: p   ! Pressure at center o volume P
      real(8), dimension (nx*ny),   intent(out) :: ro  ! Specific mass (absolute density) at center of volumes


      ! Inner variables

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


      ! South boundary

      j = 1

      do i = 2, nx-1

         np   = nx * (j-1) + i

         npn  = np + nx

         p(np) = p(npn)

         ro(np) = g(np) * p(np)

      end do



      ! North boundary

      j = ny

      do i = 2, nx-1

         np   = nx * (j-1) + i

         nps  = np - nx

         p(np) = p(nps)

         ro(np) = g(np) * p(np)

      end do




      ! West boundary

      i = 1

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npe  = np + 1

         p(np) = pin(j) - p(npe) + pin(j)

         ro(np) = g(np) * p(np)

      end do



      ! East boundary

      i = nx

      do j = 2, ny-1

         np   = nx * (j-1) + i

         npw  = np - 1

         p(np) = p(npw)

         ro(np) = g(np) * p(np)

      end do

   end subroutine get_p_and_ro_extrapolation_to_fictitious



end module user
