module coefficients
  implicit none
  !
  !  SUBROUTINES
  !  1) get_real_centroids_xy
  !  2) get_metrics
  !  3) get_radius
  !  4) get_density_at_nodes
  !  5) get_density_at_faces
  !  6) get_u_coefficients
  !  7) get_u_source
  !  8) get_v_coefficients
  !  9) get_v_source
  ! 10) get_T_coefficients_and_source
  ! 11) get_p_coefficients
  ! 12) get_p_source
  ! 13) get_velocities_at_faces
  ! 14) get_internal_simplec_coefficients
  ! 15) get_pressure_density_correction_with_pl
  ! 16) get_u_v_at_real_nodes_with_pl
  ! 17) get_velocities_at_internal_faces_with_pl
  !
  ! Last update: 03 Oct 2011
  !
contains

  subroutine get_real_centroids_xy( opt, nx, ny, x, y & ! Input
       ,                            xp, yp )            ! Output
    implicit none
    integer, intent(in) :: opt ! (1=simple mean, 2=weighted mean)
    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) :: x(nx*ny) ! Coord. x of the northest corner of volume P
    real(8), intent(in) :: y(nx*ny) ! Coord. y of the northest corner of volume P
    real(8), intent(out) :: xp(nx*ny) ! Coord. x of the centroid of volume P
    real(8), intent(out) :: yp(nx*ny) ! Coord. y of the centroid of volume P
    !
    integer :: i, j, np, nps, npw, npsw
    real(8) :: A1, A2

    ! Centroids of real volumes

    xp = 0.d0
    yp = 0.d0

    if ( opt == 1 ) then ! Simple mean
       ! Centroids are given by the mean of the corners
       do j = 2, ny-1
          do i = 2, nx-1

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

             xp(np) = ( x(np) + x(npw) + x(npsw) + x(nps) ) / 4.d0
             yp(np) = ( y(np) + y(npw) + y(npsw) + y(nps) ) / 4.d0

          end do
       end do
    else ! weighted mean
       ! The rectangle is divided in two triangles of area A1 and A2.
       ! The centroids of each triangle is calculated by the mean of the corners.
       ! The centroid of the rectangle is given by the weighted mean of the centroids of each triangle.
       do j = 2, ny-1
          do i = 2, nx-1

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

             A1 = 0.5d0 * ( x(npsw) * (y(nps)-y(np) ) + x(nps) * (y(np) -y(npsw)) + x(np)  * (y(npsw)-y(nps)))
             A2 = 0.5d0 * ( x(npsw) * (y(np) -y(npw)) + x(np)  * (y(npw)-y(npsw)) + x(npw) * (y(npsw)-y(np) ))

             xp(np) = ( (x(npsw)+x(nps)+x(np)) * A1 + (x(npsw)+x(np)+x(npw)) * A2 ) / ( 3.d0 * (A1+A2) )
             yp(np) = ( (y(npsw)+y(nps)+y(np)) * A1 + (y(npsw)+y(np)+y(npw)) * A2 ) / ( 3.d0 * (A1+A2) )
          end do
       end do
    end if
  end subroutine get_real_centroids_xy


  subroutine get_metrics(nx, ny, x, y, xp, yp            & ! Input
       ,          xe, ye, xen, yen, xk, yk, xke, yke, Jp & ! Output
       ,          Je, Jn, alphae, gamman, betae, betan )   ! Output
    implicit none
    integer, intent(in)  :: nx, ny  ! Number of volumes in csi and eta directions (real+fictitious)
    real(8), dimension (nx*ny), intent(in)  :: x,  y  ! Coord. of the northest corner of volume P
    real(8), dimension (nx*ny), intent(in)  :: xp, yp ! Coord. of the centroid of volume P
    real(8), dimension (nx*ny), intent(out) :: xe     ! x_eta at face east of volume P
    real(8), dimension (nx*ny), intent(out) :: ye     ! y_eta at face east of volume P
    real(8), dimension (nx*ny), intent(out) :: xen    ! x_eta at face north of volume P
    real(8), dimension (nx*ny), intent(out) :: yen    ! y_eta at face north of volume P
    real(8), dimension (nx*ny), intent(out) :: xk     ! x_csi at face north of volume P
    real(8), dimension (nx*ny), intent(out) :: yk     ! y_csi at face north of volume P
    real(8), dimension (nx*ny), intent(out) :: xke    ! x_csi at face east of volume P
    real(8), dimension (nx*ny), intent(out) :: yke    ! y_csi at face east of volume P
    real(8), dimension (nx*ny), intent(out) :: Jp     ! Jacobian at the center of volume P
    real(8), dimension (nx*ny), intent(out) :: Je     ! Jacobian at the center of east face of volume P
    real(8), dimension (nx*ny), intent(out) :: Jn     ! Jacobian at the center of north face of volume P
    real(8), dimension (nx*ny), intent(out) :: Alphae ! (metric) Alpha at the center of east face of volume P
    real(8), dimension (nx*ny), intent(out) :: Gamman ! (metric) Gamma at the center of north face of volume P
    real(8), dimension (nx*ny), intent(out) :: Betae  ! (metric) Beta  at the center of east face of volume P
    real(8), dimension (nx*ny), intent(out) :: Betan  ! (metric) Beta  at the center of north face of volume P
    !
    integer :: i, j, npsw, nps, npw, np, npe, npn
    real(8) :: fw, fe, fs, fn, xkp, xep, ykp, yep
    !
    xe  = 0.d0
    ye  = 0.d0
    xen = 0.d0
    yen = 0.d0
    xk  = 0.d0
    yk  = 0.d0
    xke = 0.d0
    yke = 0.d0
    Jp  = 0.d0
    Je  = 0.d0
    Jn  = 0.d0
    alphae = 0.d0
    gamman = 0.d0
    betae  = 0.d0
    betan  = 0.d0

    ! Derivatives relatively to eta at the center of the east face

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

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

          xe(np) = x(np) - x(nps)
          ye(np) = y(np) - y(nps)

          alphae(np) = xe(np)**2 + ye(np)**2

       end do
    end do

    ! Derivatives relatively to csi at the center of the north face

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

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

          xk(np) = x(np) - x(npw)
          yk(np) = y(np) - y(npw)

          gamman(np) = xk(np)**2 + yk(np)**2

       end do
    end do

    ! Derivatives relatively to csi at the center of the east face (only inner real faces)

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

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

          xke(np) = xp(npe) - xp(np)
          yke(np) = yp(npe) - yp(np)

       end do
    end do

    ! Derivatives relatively to csi at the center of the east face (only west boundary of the domain)

    i = 2
    do j = 2, ny-1

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

       fw = ( x(npw) + x(npsw) ) / 2.d0
       fe = ( x(np)  + x(nps)  ) / 2.d0

       xke(npw) = -3.d0 * fw + 4.d0 * xp(np) -fe

       fw = ( y(npw) + y(npsw) ) / 2.d0
       fe = ( y(np)  + y(nps)  ) / 2.d0

       yke(npw) = -3.d0 * fw + 4.d0 * yp(np) -fe

    end do

    ! Derivatives relatively to csi at the center of the east face (only east boundary of the domain)

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

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

       fw = ( x(npw) + x(npsw) ) / 2.d0
       fe = ( x(np)  + x(nps)  ) / 2.d0

       xke(np) = fw - 4.d0 * xp(np) + 3.d0 * fe

       fw = ( y(npw) + y(npsw) ) / 2.d0
       fe = ( y(np)  + y(nps)  ) / 2.d0

       yke(np) = fw - 4.d0 * yp(np) + 3.d0 * fe

    end do

    ! Derivatives relatively to eta at the center of the north face (only inner real faces)

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

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

          xen(np) = xp(npn) - xp(np)
          yen(np) = yp(npn) - yp(np)

       end do
    end do

    ! Derivatives relatively to eta at the center of the north face (only south boundary of the domain)

    j = 2

    do i = 2, nx-1

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

       fs = ( x(nps) + x(npsw) ) / 2.d0
       fn = ( x(np)  + x(npw)  ) / 2.d0

       xen(nps) = -3.d0 * fs + 4.d0 * xp(np) -fn

       fs = ( y(nps) + y(npsw) ) / 2.d0
       fn = ( y(np)  + y(npw)  ) / 2.d0

       yen(nps) = -3.d0 * fs + 4.d0 * yp(np) -fn

    end do

    ! Derivatives relatively to eta at the center of the north face (only north boundary of the domain)

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

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

       fs = ( x(nps) + x(npsw) ) / 2.d0
       fn = ( x(np)  + x(npw)  ) / 2.d0

       xen(np) = fs - 4.d0 * xp(np) + 3.d0 * fn

       fs = ( y(nps) + y(npsw) ) / 2.d0
       fn = ( y(np)  + y(npw)  ) / 2.d0

       yen(np) = fs - 4.d0 * yp(np) + 3.d0 * fn

    end do

    ! Beta and J at the center of the east face (all real faces)

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

          np   = nx * (j-1) + i

          betae(np) = xke(np) * xe(np) + yke(np) * ye(np)

          Je(np) = 1.d0 / ( xke(np) * ye(np) - xe(np) * yke(np) )

       end do
    end do

    ! Beta and J at center of the north face (all real faces)

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

          np   = nx * (j-1) + i

          betan(np) = xk(np) * xen(np) + yk(np) * yen(np)

          Jn(np) = 1.d0 / ( xk(np) * yen(np) - xen(np) * yk(np) )

       end do
    end do

    ! Jacobian J at the center of all real volumes

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

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

          fw = ( x(npw) + x(npsw) ) / 2.d0
          fe = ( x(np)  + x(nps)  ) / 2.d0
          fs = ( x(nps) + x(npsw) ) / 2.d0
          fn = ( x(np)  + x(npw)  ) / 2.d0

          xkp = fe - fw ! (x_csi)_P
          xep = fn - fs ! (x_eta)_P

          fw = ( y(npw) + y(npsw) ) / 2.d0
          fe = ( y(np)  + y(nps)  ) / 2.d0
          fs = ( y(nps) + y(npsw) ) / 2.d0
          fn = ( y(np)  + y(npw)  ) / 2.d0

          ykp = fe - fw ! (y_csi)_P
          yep = fn - fs ! (y_eta)_P

          Jp(np) = 1.d0 / ( xkp * yep - xep * ykp )

       end do
    end do

    Je = dabs(Je)
    Jn = dabs(Jn)
    Jp = dabs(Jp)

  end subroutine get_metrics

  subroutine get_radius(coord, nx, ny, y, yp & ! Input
       ,                radius, re, rn, rp   ) ! Output
    implicit none
    integer, intent(in)  :: coord   ! Coordinate system ( 1=cylindrical, 0 = cartesian)
    integer, intent(in)  :: nx, ny  ! Number of volumes in csi and eta directions (real+fictitious)
    real(8), dimension (nx*ny), intent(in)  :: y  ! Coord. y of the northeast corner of volume P
    real(8), dimension (nx*ny), intent(in)  :: yp ! Coord. y of the centroid of volume P
    real(8), dimension (nx*ny), intent(out) :: radius ! Radius of northest corner of volume P
    real(8), dimension (nx*ny), intent(out) :: re ! Radius of the center of east face of volume P
    real(8), dimension (nx*ny), intent(out) :: rn ! Radius of the center of north face of volume P
    real(8), dimension (nx*ny), intent(out) :: rp ! Radius of the center of volume P

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

    ! Cartesian coordinate system

    radius = 1.d0
    re = 1.d0
    rn = 1.d0
    rp = 1.d0

    ! Cylindrical coordinate system

    if ( coord == 1 ) then

       ! Radius of northeast corner

       radius = y

       ! Radius at the center of the east face

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

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

             re(np) = ( y(np) + y(nps) ) / 2.d0

          end do
       end do

       ! Radius at the center of the north face

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

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

             rn(np) = ( y(np) + y(npw) ) / 2.d0

          end do
       end do

       ! Radius of the center of the volume

       rp = yp

       ! Radius  of the center of the volume (south fictitious)

       j = 1
       do i = 2, nx-1

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

          !rp(np) = rn(np) - rp(npn) + rn(np)
          rp(np) = rp(npn) - rp(npn+nx) + rp(npn)

       end do

       ! Radius  of the center of the volume (north fictitious)

       j = ny
       do i = 2, nx-1

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

          !rp(np) = rn(nps) - rp(nps) + rn(nps)
          rp(np) = rp(nps) - rp(nps-nx) + rp(nps)

       end do

       ! Radius  of the center of the volume (west fictitious)

       i = 1
       do j = 2, ny-1

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

          !rp(np) = re(np) - rp(npe) + re(np)
          rp(np) = rp(npe) - rp(npe+1) + rp(npe)

       end do

       ! Radius  of the center of the volume (east fictitious)

       i = nx
       do j = 2, ny-1

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

          !rp(np) = re(npw) - rp(npw) + re(npw)
          rp(np) = rp(npw) - rp(npw-1) + rp(npw)

       end do

       ! Radius  of the center of the volume (corner fictitious: SW)

       i = 1 ; j = 1

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

       rp(np) = rp(npn) - rp(npn+nx) + rp(npn)

       ! Radius  of the center of the volume (corner fictitious: SE)

       i = nx; j = 1

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

       rp(np) = rp(npn) - rp(npn+nx) + rp(npn)

       ! Radius  of the center of the volume (corner fictitious: NW)

       i = 1; j = ny

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

       rp(np) = rp(nps) - rp(nps-nx) + rp(nps)

       ! Radius  of the center of the volume (corner fictitious: NE)

       i = nx; j = ny

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

       rp(np) = rp(nps) - rp(nps-nx) + rp(nps)

    end if

  end subroutine get_radius


  subroutine get_density_at_nodes( nx, ny, Rg, p, T, ro) ! ro 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) :: Rg ! Perfect gas constant
    real(8), dimension(nx*ny), intent(in)  :: p  ! Pressure at center of volumes
    real(8), dimension(nx*ny), intent(in)  :: T  ! Temperature at center of volumes
    real(8), dimension(nx*ny), intent(out) :: ro ! Specific mass (absolute density) at center of volumes

    ro = p / ( Rg * T )

  end subroutine get_density_at_nodes


  subroutine get_density_at_faces( nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron 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) :: beta ! Constant of the UDS/CDS mixing scheme
    real(8), dimension(nx*ny), intent(in)  :: ro  ! Specific mass (absolute density) at center of volumes
    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), dimension(nx*ny), intent(out) :: roe ! Absolute density at east face
    real(8), dimension(nx*ny), intent(out) :: ron ! Absolute density at north face

    integer :: i, j, np, npe, npn
    real(8) :: ae, an ! Coefficients of UDS scheme

    ! Density at east face

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

          ae = dsign(0.5d0,Uce(np))
          roe(np) = (0.5d0+ae) * ro(np) + (0.5d0-ae) * ro(npe) + beta * ae * ( ro(npe)-ro(np) )

       end do
    end do

    ! Density at north face

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

          an = dsign(0.5d0,Vcn(np))
          ron(np) = (0.5d0+an) * ro(np) + (0.5d0-an) * ro(npn) + beta * an * ( ro(npn)-ro(np) )

       end do
    end do

  end subroutine get_density_at_faces


  subroutine get_u_coefficients( nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
       ,                        ye, yk, alphae, betae, betan, gamman       &
       ,                        vle, vln, roe, ron, roa, Uce, Vcn, a ) ! a 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)
    integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
    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) :: Je     ! Jacobian at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: ye     ! y_eta at face east of volume P
    real(8), dimension (nx*ny), intent(in) :: yk     ! y_csi at face north of volume P
    real(8), dimension (nx*ny), intent(in) :: Alphae ! (metric) Alpha at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Betae  ! (metric) Beta  at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Betan  ! (metric) Beta  at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: Gamman ! (metric) Gamma at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
    real(8), dimension (nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
    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) :: 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), dimension (nx*ny,9), intent(out) :: a   ! Coefficients of the linear system

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe
    real(8) :: fmw, fme, fms, fmn, mpa
    real(8) :: ae, aw, an, as
    real(8) :: S1p, S2p

    ! Contribution to the coefficients due to the advection term
    if ( modvis == 0 ) then ! Euler
       do j = 2, ny-1
          do i = 2, nx-1

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

             fme = roe(np ) * re(np ) * Uce(np )
             fmw = roe(npw) * re(npw) * Uce(npw)
             fmn = ron(np ) * rn(np ) * Vcn(np )
             fms = ron(nps) * rn(nps) * Vcn(nps)

             mpa = roa(np) * rp(np) / Jp(np)

             as = dsign( 0.5d0, Vcn(nps) )
             an = dsign( 0.5d0, Vcn(np)  )
             aw = dsign( 0.5d0, Uce(npw) )
             ae = dsign( 0.5d0, Uce(np)  )

             a(np,1) = 0.d0 ! SW
             a(np,3) = 0.d0 ! SE
             a(np,7) = 0.d0 ! NW
             a(np,9) = 0.d0 ! NE

             a(np,2) =-fms * (0.5d0+as) ! S
             a(np,4) =-fmw * (0.5d0+aw) ! W
             a(np,6) = fme * (0.5d0-ae) ! E
             a(np,8) = fmn * (0.5d0-an) ! N

             a(np,5) = mpa / dt - ( a(np,8) +  a(np,2) + a(np,4) + a(np,6) )

          end do
       end do

    else ! Navier-Stokes

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

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

             fme = roe(np ) * re(np ) * Uce(np )
             fmw = roe(npw) * re(npw) * Uce(npw)
             fmn = ron(np ) * rn(np ) * Vcn(np )
             fms = ron(nps) * rn(nps) * Vcn(nps)

             mpa = roa(np) * rp(np) / Jp(np)

             as = dsign( 0.5d0, Vcn(nps) )
             an = dsign( 0.5d0, Vcn(np)  )
             aw = dsign( 0.5d0, Uce(npw) )
             ae = dsign( 0.5d0, Uce(np)  )

             ! Contribution to the coefficients due to the advection and diffusion term

             a(np,1) =  ( vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  +       vln(nps) * rn(nps) * Jn(nps) * betan(nps) ) / 4.d0


             a(np,2) =  - fms * (0.5d0+as) & ! S
                  -       vln(nps) * rn(nps) * Jn(nps) * gamman(nps) &
                  +     ( vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       vle(np ) * re(np ) * Je(np ) * betae(np ) ) / 4.d0


             a(np,3) = -( vle(np ) * re(np ) * Je(np ) * betae(np )  &
                  +       vln(nps) * rn(nps) * Jn(nps) * betan(nps) ) / 4.d0


             a(np,4) = -  fmw * (0.5d0+aw) & ! W
                  -       vle(npw) * re(npw) * Je(npw) * alphae(npw) &
                  +     ( vln(nps) * rn(nps) * Jn(nps) * betan(nps)  &
                  -       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,6) =    fme * (0.5d0-ae) & ! E
                  -       vle(np ) * re(np ) * Je(np ) * alphae(np ) &
                  +     (-vln(nps) * rn(nps) * Jn(nps) * betan(nps)  &
                  +       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,7) =  (-vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,8) =    fmn * (0.5d0-an) & ! N
                  -       vln(np ) * rn(np ) * Jn(np ) * gamman(np ) &
                  -     ( vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       vle(np ) * re(np ) * Je(np ) * betae(np ) ) / 4.d0


             a(np,9) =  ( vle(np ) * re(np ) * Je(np ) * betae(np )  &
                  +       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,5) = mpa / dt - ( a(np,8) +  a(np,2) + a(np,4) + a(np,6) )


             ! Contribution to the coefficients due to the source term

             S1p =  rp(np) * Je(np ) * vle(np ) * ye(np ) ** 2 / 3.d0 &
                  + rp(np) * Je(npw) * vle(npw) * ye(npw) ** 2 / 3.d0

             S2p =  rp(np) * Jn(np ) * vln(np ) * yk(np ) ** 2 / 3.d0 &
                  + rp(np) * Jn(nps) * vln(nps) * yk(nps) ** 2 / 3.d0

             a(np,5) = a(np,5) + S1p + S2p

          end do
       end do

    end if

  end subroutine get_u_coefficients


  subroutine get_u_source( nx, ny, modvis, beta, dt, rp, re, rn   &
       ,                  xe, ye, xk, yk, xke, yke, xen, yen     &
       ,                  Jp, Je, Jn, roe, ron, roa, p, vle, vln &
       ,                  Uce, Vcn, ua, u, v, cup, sup, b ) ! The last 3 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)
    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) :: 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) :: xe     ! x_eta at center of face east
    real(8), dimension (nx*ny), intent(in) :: ye     ! y_eta at center of face east
    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) :: xke    ! x_csi at center of face east
    real(8), dimension (nx*ny), intent(in) :: yke    ! y_csi at center of face east
    real(8), dimension (nx*ny), intent(in) :: xen    ! x_eta at center of face north
    real(8), dimension (nx*ny), intent(in) :: yen    ! y_eta at center of face north
    real(8), dimension (nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
    real(8), dimension (nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face 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) :: roa    ! Absolute density at a time step before at center of vol. P
    real(8), dimension (nx*ny), intent(in) :: p      ! Pressure at center o volume P
    real(8), dimension (nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
    real(8), dimension (nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
    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), dimension (nx*ny), intent(in) :: ua     ! Cartesian velocity of a time step before
    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) :: cup   ! Term of deferred correction for u
    real(8), dimension (nx*ny), intent(out) :: sup   ! Viscous term for u
    real(8), dimension (nx*ny), intent(out) :: b     ! Source vector of the linear system

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe, npsw, npse, npnw, npne
    real(8) :: fmw, fme, fmn, fms, mpa
    real(8) :: as, an, aw, ae
    real(8) :: S1, S2, S3, S4, S5, S6

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

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

          fme = roe(np ) * re(np ) * Uce(np )
          fmw = roe(npw) * re(npw) * Uce(npw)
          fmn = ron(np ) * rn(np ) * Vcn(np )
          fms = ron(nps) * rn(nps) * Vcn(nps)

          mpa = roa(np) * rp(np) / Jp(np)

          as = dsign( 0.5d0, Vcn(nps) )
          an = dsign( 0.5d0, Vcn(np)  )
          aw = dsign( 0.5d0, Uce(npw) )
          ae = dsign( 0.5d0, Uce(np)  )

          ! Contribution to b due to advection (UDS)

          b(np) = mpa * ua(np) / dt

          ! Contribution to b due to advection (Deferred correction)


          cup(np) = - beta * &
               (  fme * ae * ( u(npe) - u(np ) ) &
               -  fmw * aw * ( u(np ) - u(npw) ) &
               +  fmn * an * ( u(npn) - u(np ) ) &
               -  fms * as * ( u(np ) - u(nps) ) &
               )

          b(np) = b(np) + cup(np)


          ! Contribution to b due to pressure term

          b(np) = b(np) +  0.5d0 * rp(np) * ( &
               + yk(np ) * ( p(npn) + p(np) ) &
               - yk(nps) * ( p(nps) + p(np) ) &
               - ye(np ) * ( p(npe) + p(np) ) &
               + ye(npw) * ( p(npw) + p(np) ) &
               )

       end do
    end do

    ! Contribution to b due to viscous term

    if ( modvis == 1 ) then
       do j = 2, ny-1
          do i = 2, nx-1

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



             S1 = rp(np) / 3.d0 * (  &

                  + Je(np ) * vle(np ) * &

                  ( ye(np )**2 * u(npe) - yke(np ) * ye(np ) * ( u(npn) + u(npne) - u(nps) - u(npse) ) / 4.d0 ) &

                  - Je(npw) * vle(npw) * &

                  (-ye(npw)**2 * u(npw) - yke(npw) * ye(npw) * ( u(npn) + u(npnw) - u(nps) - u(npsw) ) / 4.d0 ) )



             S2 = rp(np) / 3.d0 * ( &

                  + Jn(np ) * vln(np ) * &

                  ( yk(np )**2 * u(npn) - yk(np ) * yen(np ) * ( u(npe) + u(npne) - u(npw) - u(npnw) ) / 4.d0 ) &

                  - Jn(nps) * vln(nps) * &

                  (-yk(nps)**2 * u(nps) - yk(nps) * yen(nps) * ( u(npe) + u(npse) - u(npw) - u(npsw) ) / 4.d0 ) )



             S3 =   Je(np ) * re(np ) * xe(np ) * &

                  ( yke(np ) * (v(npn)+v(npne)-v(nps)-v(npse)) / 4.d0 - ye(np ) * (v(npe)-v(np)) ) &

                  - Je(npw) * re(npw) * xe(npw) * &

                  ( yke(npw) * (v(npn)+v(npnw)-v(nps)-v(npsw)) / 4.d0 - ye(npw) * (v(np)-v(npw)) )



             S4 =   Jn(np ) * rn(np ) * vln(np ) * xk(np ) * &

                  (yen(np ) * (v(npne)+v(npe)-v(npnw)-v(npw)) / 4.d0 - yk(np ) * (v(npn)-v(np)) ) &

                  - Jn(nps) * rn(nps) * vln(nps) * xk(nps) * &

                  (yen(nps) * (v(npse)+v(npe)-v(npsw)-v(npw)) / 4.d0 - yk(nps) * (v(np)-v(nps)) )




             S5 = - 2.d0 * rp(np) / 3.d0 * Je(np ) * vle(np ) / re(np ) * ye(np ) * &

                  ( xke(np ) * ( rp(npn)*v(npn) + rp(npne)*v(npne) - rp(nps)*v(nps) - rp(npse)*v(npse) ) / 4.d0 &

                  - xe(np) * ( rp(npe)*v(npe) - rp(np)*v(np))) &

                  + 2.d0 * rp(np) / 3.d0 * Je(npw) * vle(npw) / re(npw) * ye(npw) * &

                  ( xke(npw) * ( rp(npn)*v(npn) + rp(npnw)*v(npnw) - rp(nps)*v(nps) - rp(npsw)*v(npsw) ) / 4.d0 &

                  - xe(npw) * ( rp(np)*v(np) - rp(npw)*v(npw)) )

             if ( j == 2 ) then

                S6 = - 2.d0 * rp(np) / 3.d0 * Jn(np ) * vln(np ) / rn(np ) * yk(np ) * &

                     ( xen(np) * ( rp(npne)*v(npne) + rp(npe)*v(npe) - rp(npnw)*v(npnw) - rp(npw)*v(npw) ) / 4.d0 &

                     - xk(np) * ( rp(npn)*v(npn) - rp(np)*v(np) ) )

             else


                S6 = - 2.d0 * rp(np) / 3.d0 * Jn(np ) * vln(np ) / rn(np ) * yk(np ) * &

                     ( xen(np) * ( rp(npne)*v(npne) + rp(npe)*v(npe) - rp(npnw)*v(npnw) - rp(npw)*v(npw) ) / 4.d0 &

                     - xk(np) * ( rp(npn)*v(npn) - rp(np)*v(np) ) ) &

                     + 2.d0 * rp(np) / 3.d0 * Jn(nps) * vln(nps) / rn(nps) * yk(nps) * &

                     ( xen(nps) * ( rp(npe)*v(npe) + rp(npse)*v(npse) - rp(npw)*v(npw) - rp(npsw)*v(npsw) ) / 4.d0 &

                     - xk(nps) * ( rp(np)*v(np) - rp(nps)*v(nps) ) )
             end if

             sup(np) = S1 + S2 + S3 + S4 + S5 + S6

             b(np) = b(np) + sup(np)

          end do
       end do

    end if
  end subroutine get_u_source


  subroutine get_v_coefficients( nx, ny, coord, modvis, dt, rp, re, rn, Jp, Je, Jn &
       ,                        xe, xk, alphae, betae, betan, gamman       &
       ,                        vle, vln, vlp, roe, ron, roa, Uce, Vcn, a ) ! Last 1 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)
    integer, intent(in) :: coord  ! Coordinate system ( 1 = cylindrical, 0 = cartesian)
    integer, intent(in) :: modvis ! modvis = 0 -> Euler;  modvis = 1 -> Navier-Stokes
    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) :: Je     ! Jacobian at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: xe     ! y_eta at face east of volume P
    real(8), dimension (nx*ny), intent(in) :: xk     ! y_csi at face north of volume P
    real(8), dimension (nx*ny), intent(in) :: Alphae ! (metric) Alpha at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Betae  ! (metric) Beta  at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Betan  ! (metric) Beta  at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: Gamman ! (metric) Gamma at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
    real(8), dimension (nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
    real(8), dimension (nx*ny), intent(in) :: vlp    ! Laminar viscosity at 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) :: 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), dimension (nx*ny,9), intent(out) :: a   ! Coefficients of the linear system

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe
    real(8) :: fmw, fme, fms, fmn, mpa
    real(8) :: ae, aw, an, as
    real(8) :: S1p, S2p, S7p

    ! Contribution to the coefficients due to the advection term
    if ( modvis == 0 ) then ! Euler
       do j = 2, ny-1
          do i = 2, nx-1

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

             fme = roe(np ) * re(np ) * Uce(np )
             fmw = roe(npw) * re(npw) * Uce(npw)
             fmn = ron(np ) * rn(np ) * Vcn(np )
             fms = ron(nps) * rn(nps) * Vcn(nps)

             mpa = roa(np) * rp(np) / Jp(np)

             as = dsign( 0.5d0, Vcn(nps) )
             an = dsign( 0.5d0, Vcn(np)  )
             aw = dsign( 0.5d0, Uce(npw) )
             ae = dsign( 0.5d0, Uce(np)  )

             a(np,1) = 0.d0 ! SW
             a(np,3) = 0.d0 ! SE
             a(np,7) = 0.d0 ! NW
             a(np,9) = 0.d0 ! NE

             a(np,2) =-fms * (0.5d0+as) ! S
             a(np,4) =-fmw * (0.5d0+aw) ! W
             a(np,6) = fme * (0.5d0-ae) ! E
             a(np,8) = fmn * (0.5d0-an) ! N

             a(np,5) = mpa / dt - ( a(np,8) +  a(np,2) + a(np,4) + a(np,6) )

          end do
       end do

    else ! Navier-Stokes

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

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

             fme = roe(np ) * re(np ) * Uce(np )
             fmw = roe(npw) * re(npw) * Uce(npw)
             fmn = ron(np ) * rn(np ) * Vcn(np )
             fms = ron(nps) * rn(nps) * Vcn(nps)

             mpa = roa(np) * rp(np) / Jp(np)

             as = dsign( 0.5d0, Vcn(nps) )
             an = dsign( 0.5d0, Vcn(np)  )
             aw = dsign( 0.5d0, Uce(npw) )
             ae = dsign( 0.5d0, Uce(np)  )

             ! Contribution to the coefficients due to the advection and diffusion term

             a(np,1) =  ( vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  +       vln(nps) * rn(nps) * Jn(nps) * betan(nps) ) / 4.d0


             a(np,2) =  - fms * (0.5d0+as) & ! S
                  -       vln(nps) * rn(nps) * Jn(nps) * gamman(nps) &
                  +     ( vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       vle(np ) * re(np ) * Je(np ) * betae(np ) ) / 4.d0


             a(np,3) = -( vle(np ) * re(np ) * Je(np ) * betae(np )  &
                  +       vln(nps) * rn(nps) * Jn(nps) * betan(nps) ) / 4.d0


             a(np,4) = -  fmw * (0.5d0+aw) & ! W
                  -       vle(npw) * re(npw) * Je(npw) * alphae(npw) &
                  +     ( vln(nps) * rn(nps) * Jn(nps) * betan(nps)  &
                  -       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,6) =    fme * (0.5d0-ae) & ! E
                  -       vle(np ) * re(np ) * Je(np ) * alphae(np ) &
                  +     (-vln(nps) * rn(nps) * Jn(nps) * betan(nps)  &
                  +       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,7) =  (-vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,8) =    fmn * (0.5d0-an) & ! N
                  -       vln(np ) * rn(np ) * Jn(np ) * gamman(np ) &
                  -     ( vle(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       vle(np ) * re(np ) * Je(np ) * betae(np ) ) / 4.d0


             a(np,9) =  ( vle(np ) * re(np ) * Je(np ) * betae(np )  &
                  +       vln(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,5) = mpa / dt - ( a(np,8) +  a(np,2) + a(np,4) + a(np,6) )


             ! Contribution to the coefficients due to the source term

             S1p =  vle(np ) * re(np ) * Je(np ) * xe(np )**2 / 3.d0 &
                  + vle(npw) * re(npw) * Je(npw) * xe(npw)**2 / 3.d0

             S2p =  vln(np ) * rn(np ) * Jn(np ) * xk(np )**2 / 3.d0 &
                  + vln(nps) * rn(nps) * Jn(nps) * xk(nps)**2 / 3.d0

             S7p = coord * vlp(np) * 4.d0 / ( 3.d0 * rp(np) * Jp(np) )

             a(np,5) = a(np,5) + S1p + S2p + S7p

          end do
       end do

    end if

  end subroutine get_v_coefficients


  subroutine get_v_source( nx, ny, coord, modvis, beta, dt, rp, re, rn   &
       ,                  xe, ye, xk, yk, xke, yke, xen, yen     &
       ,                  Jp, Je, Jn, roe, ron, roa, p, vle, vln &
       ,                  Uce, Vcn, va, u, v, cvp, svp, b ) ! Last 3 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)
    integer, intent(in) :: coord  ! Coordinate system ( 1 = cylindrical, 0 = cartesian)
    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) :: 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) :: xe     ! x_eta at center of face east
    real(8), dimension (nx*ny), intent(in) :: ye     ! y_eta at center of face east
    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) :: xke    ! x_csi at center of face east
    real(8), dimension (nx*ny), intent(in) :: yke    ! y_csi at center of face east
    real(8), dimension (nx*ny), intent(in) :: xen    ! x_eta at center of face north
    real(8), dimension (nx*ny), intent(in) :: yen    ! y_eta at center of face north
    real(8), dimension (nx*ny), intent(in) :: Jp     ! Jacobian at the center of volume P
    real(8), dimension (nx*ny), intent(in) :: Je     ! Jacobian at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face 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) :: roa    ! Absolute density at a time step before at center of vol. P
    real(8), dimension (nx*ny), intent(in) :: p      ! Pressure at center o volume P
    real(8), dimension (nx*ny), intent(in) :: vle    ! Laminar viscosity at center of face east
    real(8), dimension (nx*ny), intent(in) :: vln    ! Laminar viscosity at center of face north
    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), dimension (nx*ny), intent(in) :: va     ! Cartesian velocity of a time step before
    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) :: cvp   ! Term of deferred correction for v
    real(8), dimension (nx*ny), intent(out) :: svp   ! Viscous term for v
    real(8), dimension (nx*ny), intent(out) :: b     ! Source vector of the linear system

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe, npsw, npse, npnw, npne
    real(8) :: fmw, fme, fmn, fms, mpa
    real(8) :: as, an, aw, ae
    real(8) :: S1, S2, S3, S4, S5, S6, S8

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

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

          fme = roe(np ) * re(np ) * Uce(np )
          fmw = roe(npw) * re(npw) * Uce(npw)
          fmn = ron(np ) * rn(np ) * Vcn(np )
          fms = ron(nps) * rn(nps) * Vcn(nps)

          mpa = roa(np) * rp(np) / Jp(np)

          as = dsign( 0.5d0, Vcn(nps) )
          an = dsign( 0.5d0, Vcn(np)  )
          aw = dsign( 0.5d0, Uce(npw) )
          ae = dsign( 0.5d0, Uce(np)  )

          ! Contribution to b due to advection (UDS)

          b(np) = mpa * va(np) / dt

          ! Contribution to b due to advection (Deferred correction)


          cvp(np) = - beta * &
               (  fme * ae * ( v(npe) - v(np ) ) &
               -  fmw * aw * ( v(np ) - v(npw) ) &
               +  fmn * an * ( v(npn) - v(np ) ) &
               -  fms * as * ( v(np ) - v(nps) ) &
               )

          b(np) = b(np) + cvp(np)


          ! Contribution to b due to pressure term

          b(np) = b(np) +  0.5d0 * rp(np) * ( &
               + xe(np ) * ( p(np) + p(npe) ) &
               - xe(npw) * ( p(np) + p(npw) ) &
               - xk(np ) * ( p(np) + p(npn) ) &
               + xk(nps) * ( p(np) + p(nps) ) &
               )

       end do
    end do

    ! Contribution to b due to viscous term

    if ( modvis == 1 ) then

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

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


             S1 =   vle(np ) * re(np ) * Je(np ) / 3.d0 * ( xe(np )**2 * v(npe) &
                  - xke(np ) * xe(np ) * ( v(npn) + v(npne) - v(nps) - v(npse) ) / 4.d0 ) &
                  - vle(npw) * re(npw) * Je(npw) / 3.d0 * (-xe(npw)**2 * v(npw) &
                  - xke(npw) * xe(npw) * ( v(npn) + v(npnw) - v(nps) - v(npsw) ) / 4.d0 )

             S2 =   vln(np ) * rn(np ) * Jn(np ) / 3.d0 * ( xk(np )**2 * v(npn) &
                  - xk(np ) * xen(np ) * ( v(npe) + v(npne) - v(npnw) - v(npw) ) / 4.d0 ) &
                  - vln(nps) * rn(nps) * Jn(nps) / 3.d0 * (-xk(nps)**2 * v(nps) &
                  - xk(nps) * xen(nps) * ( v(npe) + v(npse) - v(npsw) - v(npw) ) / 4.d0 )

             S3 =   rp(np) * vle(np ) * Je(np ) * ye(np ) * ( &
                  + xke(np ) * ( u(npn) + u(npne) - u(nps) - u(npse) ) / 4.d0 &
                  - xe(np ) * ( u(npe) - u(np) ) ) &
                  - rp(np) * vle(npw) * Je(npw) * ye(npw) * ( &
                  + xke(npw) * ( u(npn) + u(npnw) - u(nps) - u(npsw) ) / 4.d0 &
                  - xe(npw) * ( u(np) - u(npw) ) )

             S4 =   rp(np) * vln(np ) * Jn(np ) * yk(np ) * ( &
                  + xen(np ) * ( u(npne) + u(npe) - u(npnw) - u(npw) ) / 4.d0 &
                  - xk(np ) * ( u(npn) - u(np) ) ) &
                  - rp(np) * vln(nps) * Jn(nps) * yk(nps) * ( &
                  + xen(nps) * ( u(npse) + u(npe) - u(npsw) - u(npw) ) / 4.d0 &
                  - xk(nps) * ( u(np) - u(nps) ) )

             S5 = - 2.d0 / 3.d0 * rp(np) * vle(np ) * Je(np ) * xe(np ) * ( &
                  + yke(np ) * ( u(npn) + u(npne) - u(nps) - u(npse) ) / 4.d0 &
                  - ye(np ) * ( u(npe) - u(np) ) ) &
                  + 2.d0 / 3.d0 * rp(np) * vle(npw) * Je(npw) * xe(npw) * ( &
                  + yke(npw) * ( u(npn) + u(npnw) - u(nps) - u(npsw) ) / 4.d0 &
                  - ye(npw) * ( u(np) - u(npw) ) )

             S6 = - 2.d0 / 3.d0 * rp(np) * vln(np ) * Jn(np ) * xk(np ) * ( &
                  + yen(np ) * ( u(npe) + u(npne) - u(npw) - u(npnw) ) / 4.d0 &
                  - yk(np ) * ( u(npn) - u(np) ) ) &
                  + 2.d0 / 3.d0 * rp(np) * vln(nps) * Jn(nps) * xk(nps) * ( &
                  + yen(nps) * ( u(npe) + u(npse) - u(npw) - u(npsw) ) / 4.d0 &
                  - yk(nps) * ( u(np) - u(nps) ) )

             S8 = 2.d0/3.d0 * coord * v(np) * ( &
                  - vln(np) * xk(np) + vln(nps) * xk(nps) &
                  + vle(np) * xe(np) - vle(npw) * xe(npw) &
                  )

             svp(np) = S1 + S2 + S3 + S4 + S5 + S6 + S8

             b(np) = b(np) + svp(np)

          end do
       end do

    end if

  end subroutine get_v_source

  subroutine get_T_coefficients_and_source( nx, ny, coord, modvis, beta, dt, rp, re, rn       &
       ,                  xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, Jn, roe, ron  &
       ,                  roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, a, b ) ! Last 2 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)
    integer, intent(in) :: coord  ! Coordinate system ( 1 = cylindrical, 0 = cartesian)
    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) :: 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) :: xe     ! x_eta at center of face east
    real(8), dimension (nx*ny), intent(in) :: ye     ! y_eta at center of face east
    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) :: Alphae ! (metric) Alpha at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Betae  ! (metric) Beta  at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Betan  ! (metric) Beta  at the center of north face of volume P
    real(8), dimension (nx*ny), intent(in) :: Gamman ! (metric) Gamma at 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) :: Je     ! Jacobian at the center of east face of volume P
    real(8), dimension (nx*ny), intent(in) :: Jn     ! Jacobian at the center of north face 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) :: roa    ! Absolute density at a time step before at center of vol. P
    real(8), dimension (nx*ny), intent(in) :: p      ! Pressure at center o volume P
    real(8), dimension (nx*ny), intent(in) :: pa     ! Pressure at center o volume P at a time step before
    real(8), dimension (nx*ny), intent(in) :: cp     ! Specific heat at constant pressure at center o volume P
    real(8), dimension (nx*ny), intent(in) :: vlp    ! Laminar viscosity at center of volume P
    real(8), dimension (nx*ny), intent(in) :: ke     ! Thermal conductivity at center of face east
    real(8), dimension (nx*ny), intent(in) :: kn     ! Thermal conductivity at center of face north
    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), 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(in) :: T      ! Temperature of the last iteraction
    real(8), dimension (nx*ny), intent(in) :: Ta     ! Temperature at the time step before

    real(8), dimension (nx*ny,9), intent(out) :: a   ! Coefficients of the linear system
    real(8), dimension (nx*ny),   intent(out) :: b   ! Source vector of the linear system

    ! Auxiliary variables
    integer :: i, j
    integer :: np, nps, npn, npw, npe
    real(8) :: fmw, fme, fms, fmn, mpa, pup, pvp
    real(8) :: as, an, aw, ae
    real(8) :: S1, S2, S3, S4, S5

    if ( modvis == 0 ) then ! Euler

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

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

             fme = roe(np ) * re(np ) * Uce(np )
             fmw = roe(npw) * re(npw) * Uce(npw)
             fmn = ron(np ) * rn(np ) * Vcn(np )
             fms = ron(nps) * rn(nps) * Vcn(nps)

             mpa = roa(np) * rp(np) / Jp(np)

             as = dsign( 0.5d0, Vcn(nps) )
             an = dsign( 0.5d0, Vcn(np)  )
             aw = dsign( 0.5d0, Uce(npw) )
             ae = dsign( 0.5d0, Uce(np)  )

             ! Contribution to the COEFFICIENTS due to advection (UDS)

             a(np,1) = 0.d0 ! SW
             a(np,3) = 0.d0 ! SE
             a(np,7) = 0.d0 ! NW
             a(np,9) = 0.d0 ! NE

             a(np,2) =-fms * (0.5d0+as) * cp(np) ! S
             a(np,4) =-fmw * (0.5d0+aw) * cp(np) ! W
             a(np,6) = fme * (0.5d0-ae) * cp(np) ! E
             a(np,8) = fmn * (0.5d0-an) * cp(np) ! N

             a(np,5) = mpa  * cp(np) / dt - ( a(np,8) +  a(np,2) + a(np,4) + a(np,6) )

             ! Contribution to the SOURCE due to advection (UDS)

             b(np) = mpa * Ta(np) * cp(np) / dt

             ! Contribution to the SOURCE due to advection (Deferred correction)

             b(np) = b(np) - beta * cp(np) *        &
                  (  fme * ae * ( T(npe) - T(np ) ) &
                  -  fmw * aw * ( T(np ) - T(npw) ) &
                  +  fmn * an * ( T(npn) - T(np ) ) &
                  -  fms * as * ( T(np ) - T(nps) ) &
                  )

             ! Contribution to the SOURCE due to pressure term

             pup  =   0.5d0 *  rp(np)  * ( &
                  + yk(np ) * ( p(npn) + p(np) ) &
                  - yk(nps) * ( p(nps) + p(np) ) &
                  - ye(np ) * ( p(npe) + p(np) ) &
                  + ye(npw) * ( p(npw) + p(np) ) &
                  )

             pvp  =   0.5d0 *   rp(np) * ( &
                  + xe(np ) * ( p(np) + p(npe) ) &
                  - xe(npw) * ( p(np) + p(npw) ) &
                  - xk(np ) * ( p(np) + p(npn) ) &
                  + xk(nps) * ( p(np) + p(nps) ) &
                  )

             b(np) = b(np) + rp(np) / Jp(np) * ( p(np) - pa(np) ) / dt &
                  - u(np) * pup - v(np) * pvp

          end do
       end do

    else ! Navier-Stokes

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

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

             fme = roe(np ) * re(np ) * Uce(np )
             fmw = roe(npw) * re(npw) * Uce(npw)
             fmn = ron(np ) * rn(np ) * Vcn(np )
             fms = ron(nps) * rn(nps) * Vcn(nps)

             mpa = roa(np) * rp(np) / Jp(np)

             as = dsign( 0.5d0, Vcn(nps) )
             an = dsign( 0.5d0, Vcn(np)  )
             aw = dsign( 0.5d0, Uce(npw) )
             ae = dsign( 0.5d0, Uce(np)  )

             ! Contribution to the COEFFICIENTS due to advection-diffusion (UDS-CDS)

             a(np,1) =  ( ke(npw) * re(npw) * Je(npw) * betae(npw)  &
                  +       kn(nps) * rn(nps) * Jn(nps) * betan(nps) ) / 4.d0


             a(np,2) =  - fms * (0.5d0+as) * cp(np) & ! S
                  -       kn(nps) * rn(nps) * Jn(nps) * gamman(nps) &
                  +     ( ke(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       ke(np ) * re(np ) * Je(np ) * betae(np ) ) / 4.d0


             a(np,3) = -( ke(np ) * re(np ) * Je(np ) * betae(np )  &
                  +       kn(nps) * rn(nps) * Jn(nps) * betan(nps) ) / 4.d0


             a(np,4) = -  fmw * (0.5d0+aw) * cp(np) & ! W
                  -       ke(npw) * re(npw) * Je(npw) * alphae(npw) &
                  +     ( kn(nps) * rn(nps) * Jn(nps) * betan(nps)  &
                  -       kn(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,6) =    fme * (0.5d0-ae) * cp(np) & ! E
                  -       ke(np ) * re(np ) * Je(np ) * alphae(np ) &
                  +     (-kn(nps) * rn(nps) * Jn(nps) * betan(nps)  &
                  +       kn(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,7) =  (-ke(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       kn(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,8) =    fmn * (0.5d0-an) * cp(np) & ! N
                  -       kn(np ) * rn(np ) * Jn(np ) * gamman(np ) &
                  -     ( ke(npw) * re(npw) * Je(npw) * betae(npw)  &
                  -       ke(np ) * re(np ) * Je(np ) * betae(np ) ) / 4.d0


             a(np,9) =  ( ke(np ) * re(np ) * Je(np ) * betae(np )  &
                  +       kn(np ) * rn(np ) * Jn(np ) * betan(np ) ) / 4.d0


             a(np,5) = mpa  * cp(np) / dt - ( a(np,8) +  a(np,2) + a(np,4) + a(np,6) )


             ! Contribution to the SOURCE due to advection (UDS)

             b(np) = mpa * Ta(np) * cp(np) / dt

             ! Contribution to the SOURCE due to advection (Deferred correction)

             b(np) = b(np) - beta * cp(np) *        &
                  (  fme * ae * ( T(npe) - T(np ) ) &
                  -  fmw * aw * ( T(np ) - T(npw) ) &
                  +  fmn * an * ( T(npn) - T(np ) ) &
                  -  fms * as * ( T(np ) - T(nps) ) &
                  )

             ! Contribution to the SOURCE due to pressure term

             pup  =   0.5d0 *  rp(np)  * ( &
                  + yk(np ) * ( p(npn) + p(np) ) &
                  - yk(nps) * ( p(nps) + p(np) ) &
                  - ye(np ) * ( p(npe) + p(np) ) &
                  + ye(npw) * ( p(npw) + p(np) ) &
                  )

             pvp  =   0.5d0 *   rp(np) * ( &
                  + xe(np ) * ( p(np) + p(npe) ) &
                  - xe(npw) * ( p(np) + p(npw) ) &
                  - xk(np ) * ( p(np) + p(npn) ) &
                  + xk(nps) * ( p(np) + p(nps) ) &
                  )

             b(np) = b(np) + rp(np) / Jp(np) * ( p(np) - pa(np) ) / dt &
                  - u(np) * pup - v(np) * pvp

             ! Contribution to the SOURCE due to the viscous term

             S1 = 0.5d0 * vlp(np) * Jp(np) * ( &
                  + ye(np) * ( u(np) + u(npe) ) - ye(npw) * ( u(np) + u(npw) ) &
                  - yk(np) * ( u(np) + u(npn) ) + yk(nps) * ( u(np) + u(nps) ) &
                  )**2

             S2 = 0.5d0 * vlp(np) * Jp(np) * ( &
                  + xk(np) * ( v(np) + v(npn) ) - xk(nps) * ( v(np) + v(nps) ) &
                  - xe(np) * ( v(np) + v(npe) ) + xe(npw) * ( v(np) + v(npw) ) &
                  )**2

             S3 = 2.d0 * coord * vlp(np) / Jp(np) * ( v(np) / rp(np) )**2

             S4 = 0.25d0 * vlp(np) * Jp(np) * ( &
                  + ye(np ) * ( v(np) + v(npe) ) - xe(np ) * ( u(np) + u(npe) ) &
                  - ye(npw) * ( v(np) + v(npw) ) + xe(npw) * ( u(np) + u(npw) ) &
                  + xk(np ) * ( u(np) + u(npn) ) - yk(np ) * ( v(np) + v(npn) ) &
                  - xk(nps) * ( u(np) + u(nps) ) + yk(nps) * ( v(np) + v(nps) ) &
                  )**2

             S5 = - 2.d0 / 3.d0 * vlp(np) * Jp(np) * (      &
                  + Uce(np) - Uce(npw) + Vcn(np) - Vcn(nps) &
                  + coord * v(np) / ( rp(np) * Jp(np) ) )**2

             b(np) = b(np) + rp(np) * ( S1 + S2 + S3 + S4 + S5 )

          end do
       end do

    end if

  end subroutine get_T_coefficients_and_source


    ! Calculates the coefficients of the linear system for pressure correction
    ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
    ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
    ! u* and v* are obtained.
    subroutine get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn &
        ,                       ro, g, de, dn, a) ! Output: last one
        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) :: Uce    ! Contravariant velocity U at east face
        real(8), dimension (nx*ny), intent(in) :: Vcn    ! Contravariant velocity V 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) :: g      ! ro/p = 1/(Rg T) for ideal gases
        real(8), dimension (nx*ny), intent(in) :: de     ! SIMPLEC coefficient for Uce
        real(8), dimension (nx*ny), intent(in) :: dn     ! SIMPLEC coefficient for Vcn
        real(8), dimension (nx*ny,5), intent(out) :: a   ! Coefficients of the linear system

        integer :: i, j, np, nps, npn, npw, npe
        real(8) :: as, an, aw, ae
        real(8) :: roe, row, ron, ros
        
        do j = 2, ny-1
            do i = 2, nx-1

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

                as = dsign( 0.5d0, Vcn(nps) )
                an = dsign( 0.5d0, Vcn(np)  )
                aw = dsign( 0.5d0, Uce(npw) )
                ae = dsign( 0.5d0, Uce(np)  )

                roe = (0.5d0 + ae) * ro(np)  + (0.5d0 - ae) * ro(npe) ! Density on east face
                row = (0.5d0 + aw) * ro(npw) + (0.5d0 - aw) * ro(np)  ! Density on west face
                ron = (0.5d0 + an) * ro(np)  + (0.5d0 - an) * ro(npn) ! Density on north face
                ros = (0.5d0 + as) * ro(nps) + (0.5d0 - as) * ro(np)  ! Density on south face

                ! South
                a(np,1) = - rn(nps) * ( (0.5d0 + as) * Vcn(nps) * g(nps) + ros * dn(nps) )

                ! West
                a(np,2) = - re(npw) * ( (0.5d0 + aw) * Uce(npw) * g(npw) + row * de(npw) )

                ! Center
                a(np,3) = g(np) * ( rp(np)/(Jp(np)*dt)  &
                    + (0.5d0 + ae) * re(np)  * Uce(np)  &
                    - (0.5d0 - aw) * re(npw) * Uce(npw) &
                    + (0.5d0 + an) * rn(np)  * Vcn(np)  &
                    - (0.5d0 - as) * rn(nps) * Vcn(nps) ) &
                    + roe * re(np)  * de(np)  &
                    + row * re(npw) * de(npw) &
                    + ron * rn(np)  * dn(np)  &
                    + ros * rn(nps) * dn(nps)

                ! East
                a(np,4) = re(np) * ( (0.5d0 - ae) * Uce(np) * g(npe) - roe * de(np) )

                ! North
                a(np,5) = rn(np) * ( (0.5d0 - an) * Vcn(np) * g(npn) - ron * dn(np) )

            end do
        end do

    end subroutine get_p_coefficients

    ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
    ! rom, Ucem and Vcnm are those of the previous iteraction
    subroutine get_p_source(nx, ny, beta, dt, rp, re, rn, Jp, rom, ro, roa, Ucem, Uce, Vcnm, Vcn, b) ! Output: last one
        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) :: beta   ! Constant of the UDS/CDS mixing scheme
        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) :: rom    ! Absolute density at center of vol. P (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: ro     ! Absolute density at center of vol. P (incorrect ro*)
        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) :: Ucem   ! Contravariant velocity U at east face (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: Uce    ! Contravariant velocity U at east face (incorrect Uce*)
        real(8), dimension (nx*ny), intent(in) :: Vcnm   ! Contravariant velocity V at north face (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: Vcn    ! Contravariant velocity V at north face (incorrect Vcn*)

        real(8), dimension (nx*ny), intent(out) :: b     ! Source vector of the linear system

        integer :: i, j, np, nps, npn, npw, npe
        real(8) :: as, an, aw, ae ! UPS coefficients
        real(8) :: roe, row, ron, ros ! Incorrect density (ro*) on the faces
        real(8) :: roem, rowm, ronm, rosm ! Density on the faces (previous iteraction)

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

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

                as = dsign( 0.5d0, Vcnm(nps) )
                an = dsign( 0.5d0, Vcnm(np)  )
                aw = dsign( 0.5d0, Ucem(npw) )
                ae = dsign( 0.5d0, Ucem(np)  )

                roem = (0.5d0+ae) * rom(np)  + (0.5d0-ae) * rom(npe) ! Density on east face of the previous iteraction
                rowm = (0.5d0+aw) * rom(npw) + (0.5d0-aw) * rom(np)  ! Density on west face of the previous iteraction
                ronm = (0.5d0+an) * rom(np)  + (0.5d0-an) * rom(npn) ! Density on north face of the previous iteraction
                rosm = (0.5d0+as) * rom(nps) + (0.5d0-as) * rom(np)  ! Density on south face of the previous iteraction

                roe  = (0.5d0+ae) * ro(np)   + (0.5d0-ae) * ro(npe)  ! Incorrect density (ro*) on east face
                row  = (0.5d0+aw) * ro(npw)  + (0.5d0-aw) * ro(np)   ! Incorrect density (ro*) on west face
                ron  = (0.5d0+an) * ro(np)   + (0.5d0-an) * ro(npn)  ! Incorrect density (ro*) on north face
                ros  = (0.5d0+as) * ro(nps)  + (0.5d0-as) * ro(np)   ! Incorrect density (ro*) on south face

                b(np) = -rp(np) * ( ro(np) - roa(np) ) / ( Jp(np) * dt ) &
                    + roem * re(np)  * ( Ucem(np ) - Uce(np)  ) &
                    - rowm * re(npw) * ( Ucem(npw) - Uce(npw) ) &
                    + ronm * rn(np)  * ( Vcnm(np)  - Vcn(np)  ) &
                    - rosm * rn(nps) * ( Vcnm(nps) - Vcn(nps) ) &
                    - roe  * re(np)  *   Ucem(np ) &
                    + row  * re(npw) *   Ucem(npw) &
                    - ron  * rn(np)  *   Vcnm(np)  &
                    + ros  * rn(nps) *   Vcnm(nps) &
                    - beta * (  &
                    + ae * ( rom(npe) - rom(np) ) * re(np)  * Ucem(np)  &
                    + aw * ( rom(npw) - rom(np) ) * re(npw) * Ucem(npw) &
                    + an * ( rom(npn) - rom(np) ) * rn(np)  * Vcnm(np)  &
                    + as * ( rom(nps) - rom(np) ) * rn(nps) * Vcnm(nps) )
            end do
        end do
    end subroutine get_p_source


   subroutine get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk     & ! Input
       ,                             xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
       ,                             au, av, roa, p, u, v, uea, vea, una, vna   & ! Input
       ,                             ue, ve, un, vn, Uce, Vcn)                    ! 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) :: xe   ! x_eta at center of face east
    real(8), dimension (nx*ny),   intent(in) :: ye   ! y_eta at center of face east
    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) :: xen  ! x_eta at center of face north
    real(8), dimension (nx*ny),   intent(in) :: yen  ! y_eta at center of face north
    real(8), dimension (nx*ny),   intent(in) :: xke  ! x_csi at center of face east
    real(8), dimension (nx*ny),   intent(in) :: yke  ! y_csi at center of face east
    real(8), dimension (nx*ny),   intent(in) :: Jp   ! Jacobian at the center of volume P
    real(8), dimension (nx*ny),   intent(in) :: cup  ! Term of deferred correction for u
    real(8), dimension (nx*ny),   intent(in) :: cvp  ! Term of deferred correction for v
    real(8), dimension (nx*ny),   intent(in) :: sup  ! Viscous term for u
    real(8), dimension (nx*ny),   intent(in) :: svp  ! Viscous term for v
    real(8), dimension (nx*ny,9), intent(in) :: au   ! Coefficients of the linear system for u
    real(8), dimension (nx*ny,9), intent(in) :: av   ! Coefficients of the linear system for v
    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) :: p    ! Pressure at center o volume P
    real(8), dimension (nx*ny),   intent(in) :: u    ! Cartesian velocity of the present iteraction
    real(8), dimension (nx*ny),   intent(in) :: v    ! Cartesian velocity of the present iteraction
    real(8), dimension (nx*ny),   intent(in) :: uea  ! Cartesian velocity u at center of east face in the last time step
    real(8), dimension (nx*ny),   intent(in) :: vea  ! Cartesian velocity v at center of east face in the last time step
    real(8), dimension (nx*ny),   intent(in) :: una  ! Cartesian velocity u at center of north face in the last time step
    real(8), dimension (nx*ny),   intent(in) :: vna  ! Cartesian velocity v at center of north face in the last time step
    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, npsw, npse, nps, npw, npe, npnw, npn, npne
    integer :: npee, npsee, npnee, npnn, npnnw, npnne
    real(8) :: mpa, mea, mna, aux
    real(8) :: sumup, sumun, sumvp, sumvn, sumue, sumve
    real(8) :: pue, pve, pun, pvn

    ! Calculation of Uce at the inner faces of the real domain

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

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

          npee  = npe  + 1
          npsee = npse + 1
          npnee = npne + 1

          mpa = roa(np ) * rp(np ) / Jp(np )
          mea = roa(npe) * rp(npe) / Jp(npe)

          sumup = au(np,1) * u(npsw) + au(np,2) * u(nps ) + au(np,3) * u(npse) + au(np,4) * u(npw ) &
               +  au(np,6) * u(npe ) + au(np,7) * u(npnw) + au(np,8) * u(npn ) + au(np,9) * u(npne)

          sumue = au(npe,1) * u(nps ) + au(npe,2) * u(npse) + au(npe,3) * u(npsee) + au(npe,4) * u(np   ) &
               +  au(npe,6) * u(npee) + au(npe,7) * u(npn ) + au(npe,8) * u(npne ) + au(npe,9) * u(npnee)

          sumvp = av(np,1) * v(npsw) + av(np,2) * v(nps ) + av(np,3) * v(npse) + av(np,4) * v(npw ) &
               +  av(np,6) * v(npe ) + av(np,7) * v(npnw) + av(np,8) * v(npn ) + av(np,9) * v(npne)

          sumve = av(npe,1) * v(nps ) + av(npe,2) * v(npse) + av(npe,3) * v(npsee) + av(npe,4) * v(np   ) &
               +  av(npe,6) * v(npee) + av(npe,7) * v(npn ) + av(npe,8) * v(npne ) + av(npe,9) * v(npnee)

          aux = ( p(npn) + p(npne) - p(nps) - p(npse) ) / 4.d0

          pue = re(np) * ( yke(np) * aux + ye(np) * ( p(np) - p(npe) ) )

          pve = re(np) * (-xke(np) * aux + xe(np) * ( p(npe) - p(np) ) )

          ue(np) = ( ( mpa + mea ) * uea(np) / dt + cup(np) + cup(npe) + sup(np) + sup(npe) &
               - sumup - sumue + 2.d0 * pue ) / ( au(np,5) + au(npe,5) )

          ve(np) = ( ( mpa + mea ) * vea(np) / dt + cvp(np) + cvp(npe) + svp(np) + svp(npe) &
               - sumvp - sumve + 2.d0 * pve ) / ( av(np,5) + av(npe,5) )

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

       end do
    end do

    ! Calculation of Vcn at the inner faces of the real domain

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

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

          npnn  = npn  + nx
          npnnw = npnn - 1
          npnne = npnn + 1

          mpa = roa(np ) * rp(np ) / Jp(np )
          mna = roa(npn) * rp(npn) / Jp(npn)

          sumup = au(np,1) * u(npsw) + au(np,2) * u(nps ) + au(np,3) * u(npse) + au(np,4) * u(npw ) &
               +  au(np,6) * u(npe ) + au(np,7) * u(npnw) + au(np,8) * u(npn ) + au(np,9) * u(npne)

          sumun = au(npn,1) * u(npw ) + au(npn,2) * u(np   ) + au(npn,3) * u(npe ) + au(npn,4) * u(npnw ) &
               +  au(npn,6) * u(npne) + au(npn,7) * u(npnnw) + au(npn,8) * u(npnn) + au(npn,9) * u(npnne)

          sumvp = av(np,1) * v(npsw) + av(np,2) * v(nps ) + av(np,3) * v(npse) + av(np,4) * v(npw ) &
               +  av(np,6) * v(npe ) + av(np,7) * v(npnw) + av(np,8) * v(npn ) + av(np,9) * v(npne)

          sumvn = av(npn,1) * v(npw ) + av(npn,2) * v(np   ) + av(npn,3) * v(npe ) + av(npn,4) * v(npnw ) &
               +  av(npn,6) * v(npne) + av(npn,7) * v(npnnw) + av(npn,8) * v(npnn) + av(npn,9) * v(npnne)

          aux = ( p(npne) + p(npe) - p(npnw) - p(npw) ) / 4.d0

          pun = rn(np) * ( yk(np) * ( p(npn) - p(np) ) - yen(np) * aux )

          pvn = rn(np) * ( xk(np) * ( p(np) - p(npn) ) + xen(np) * aux )

          un(np) = ( ( mpa + mna ) * una(np) / dt + cup(np) + cup(npn) + sup(np) + sup(npn) &
               - sumup - sumun + 2.d0 * pun ) / ( au(np,5) + au(npn,5) )

          vn(np) = ( ( mpa + mna ) * vna(np) / dt + cvp(np) + cvp(npn) + svp(np) + svp(npn) &
               - sumvp - sumvn + 2.d0 * pvn ) / ( av(np,5) + av(npn,5) )

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

       end do
    end do

  end subroutine get_velocities_at_faces

  subroutine get_internal_simplec_coefficients( nx, ny, re, rn, xe, ye, xk, yk, au, av & ! Input
       ,                                                    due, dve, dun, dvn, de, dn ) ! 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)  :: 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)  :: xe  ! x_eta at center of face east
    real(8), dimension (nx*ny),   intent(in)  :: ye  ! y_eta at center of face east
    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,9), intent(in)  :: au  ! Coefficients of the linear system for u
    real(8), dimension (nx*ny,9), intent(in)  :: av  ! Coefficients of the linear system for v
    real(8), dimension (nx*ny),   intent(out) :: due ! Simplec coef. for the cartesian velocity u (east face)
    real(8), dimension (nx*ny),   intent(out) :: dve ! Simplec coef. for the cartesian velocity v (east face)
    real(8), dimension (nx*ny),   intent(out) :: dun ! Simplec coef. for the cartesian velocity u (north face)
    real(8), dimension (nx*ny),   intent(out) :: dvn ! Simplec coef. for the cartesian velocity v (north face)
    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)

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

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

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

          due(np) = 2.d0 * re(np) * ye(np) / ( sum(au(np,:)) + sum(au(npe,:)) )
          dve(np) = 2.d0 * re(np) * xe(np) / ( sum(av(np,:)) + sum(av(npe,:)) )

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

       end do
    end do

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

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

          dun(np) = 2.d0 * rn(np) * yk(np) / ( sum(au(np,:)) + sum(au(npn,:)) )
          dvn(np) = 2.d0 * rn(np) * xk(np) / ( sum(av(np,:)) + sum(av(npn,:)) )

          dn(np) = xk(np) * dvn(np) + yk(np) * dun(np)

       end do
    end do

  end subroutine get_internal_simplec_coefficients

  subroutine get_pressure_density_correction_with_pl( nx, ny, pl, g & ! Input
       ,                                                      ro, p ) ! Input and 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)    :: pl ! Pressure correction
    real(8), dimension (nx*ny), intent(in)    :: g  ! ro/p = 1/(Rg T) for ideal gases
    real(8), dimension (nx*ny), intent(inout) :: p  ! p of present iteraction is input. Corrected p is output.
    real(8), dimension (nx*ny), intent(inout) :: ro ! ro of present iteraction is input. Corrected ro is output.

    p = p + pl

    ro = ro + pl * g

    !ro = p * g

  end subroutine get_pressure_density_correction_with_pl

  subroutine get_u_v_at_real_nodes_with_pl( nx, ny, xe, ye, xk & ! Input
       ,                                    yk, rp, pl, au, av & ! Input
       ,                                    u, v )               ! Input and 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)    :: xe  ! x_eta at center of face east
    real(8), dimension (nx*ny),   intent(in)    :: ye  ! y_eta at center of face east
    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)    :: rp  ! Radius of the center of volume P
    real(8), dimension (nx*ny),   intent(in)    :: pl  ! Pressure correction
    real(8), dimension (nx*ny,9), intent(in)    :: au  ! Coefficients of the linear system for u
    real(8), dimension (nx*ny,9), intent(in)    :: av  ! Coefficients of the linear system for v
    real(8), dimension (nx*ny),   intent(inout) :: u   ! u of present iteraction is input. u corrected is output
    real(8), dimension (nx*ny),   intent(inout) :: v   ! v of present iteraction is input. v corrected is output

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe
    real(8) :: plup, plvp

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

          plup  =   0.5d0 *  rp(np)  * ( &
               +  yk(np ) * ( pl(npn) + pl(np) ) &
               -  yk(nps) * ( pl(nps) + pl(np) ) &
               -  ye(np ) * ( pl(npe) + pl(np) ) &
               +  ye(npw) * ( pl(npw) + pl(np) ) &
               )

          plvp  =   0.5d0 *   rp(np) * ( &
               +  xe(np ) * ( pl(np) + pl(npe) ) &
               -  xe(npw) * ( pl(np) + pl(npw) ) &
               -  xk(np ) * ( pl(np) + pl(npn) ) &
               +  xk(nps) * ( pl(np) + pl(nps) ) &
               )

          u(np) = u(np) + plup / sum( au(np,:) )
          v(np) = v(np) + plvp / sum( av(np,:) )

       end do
    end do

  end subroutine get_u_v_at_real_nodes_with_pl



  subroutine get_velocities_at_internal_faces_with_pl( nx, ny, xe, ye, xk, yk   & ! Input
       ,                                               due, dve, dun, dvn, pl   & ! Input
       ,                                               ue, ve, un, vn, Uce, Vcn ) ! Input and 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)    :: xe  ! x_eta at center of face east
    real(8), dimension (nx*ny), intent(in)    :: ye  ! y_eta at center of face east
    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)    :: due ! Simplec coef. for the cartesian velocity u (east face)
    real(8), dimension (nx*ny), intent(in)    :: dve ! Simplec coef. for the cartesian velocity v (east face)
    real(8), dimension (nx*ny), intent(in)    :: dun ! Simplec coef. for the cartesian velocity u (north face)
    real(8), dimension (nx*ny), intent(in)    :: dvn ! Simplec coef. for the cartesian velocity v (north face)
    real(8), dimension (nx*ny), intent(in)    :: pl  ! Pressure correction
    real(8), dimension (nx*ny), intent(inout) :: ue  ! Cartesian velocity u at center of east face
    real(8), dimension (nx*ny), intent(inout) :: ve  ! Cartesian velocity v at center of east face
    real(8), dimension (nx*ny), intent(inout) :: un  ! Cartesian velocity u at center of north face
    real(8), dimension (nx*ny), intent(inout) :: vn  ! Cartesian velocity v at center of north face
    real(8), dimension (nx*ny), intent(inout) :: Uce ! Contravariant velocity U at east face
    real(8), dimension (nx*ny), intent(inout) :: Vcn ! Contravariant velocity V at north face

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

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

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

          ue(np) = ue(np) + due(np) * ( pl(np) - pl(npe) )

          ve(np) = ve(np) + dve(np) * ( pl(npe) - pl(np) )

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

       end do
    end do

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

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

          un(np) = un(np) + dun(np) * ( pl(npn) - pl(np) )

          vn(np) = vn(np) + dvn(np) * ( pl(np) - pl(npn) )

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

       end do
    end do

  end subroutine get_velocities_at_internal_faces_with_pl

end module coefficients
