!> \brief Contains user-INDEPENDENT subroutines related to the calculation
!! of coefficients and source of the linear systems.
!<

module coefficients
    implicit none


contains

    !> \brief Calculates the specific mass at nodes using p, T and the state equation.
    subroutine get_density_at_nodes( nx, ny, Rg, p, T, ro) ! ro is output
        implicit none
        integer, intent(in) :: nx !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: Rg !< Perfect gas constant (J/(kg.K))
        real(8), dimension(nx*ny), intent(in)  :: p  !< Pressure at center o volume P (Pa)
        real(8), dimension(nx*ny), intent(in)  :: T  !< Temperature at center of volume P (K)
        real(8), dimension(nx*ny), intent(out) :: ro !< Specific mass (absolute density) at center of volume P (kg/m3)

        ro = p / ( Rg * T )

    end subroutine get_density_at_nodes


    !> \brief Interpolates the specific mass at center of volume faces using the specific mass at nodes
    !! based on the scheme UDS with deferred correction to CDS.
    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 the csi direction (real+fictitious)
        integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: beta !< UDS/CDS mixing constant (0=UDS, 1=CDS)
        real(8), dimension(nx*ny), intent(in)  :: ro  !< Specific mass (absolute density) at center of volume P (kg/m3)
        real(8), dimension(nx*ny), intent(in)  :: Uce !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension(nx*ny), intent(in)  :: Vcn !< Contravariant velocity V at north face of volume P (m2/s)
        real(8), dimension(nx*ny), intent(out) :: roe !< Specific mass at east face of volume P (kg/m3)
        real(8), dimension(nx*ny), intent(out) :: ron !< Specific mass at north face of volume P (kg/m3)

        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

    !> \brief Calculates the coefficients of the linear system for u for the real volumes.
    subroutine get_u_coefficients( nx, ny, dt, rp, re, rn, Jp &
            ,                        roe, ron, roa, Uce, Vcn, a ) ! Last one is output
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: roa    !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)

        real(8), dimension (nx*ny,9), intent(out) :: a   !< Coefficients of the linear system for u (kg/s)

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

        ! Contribution to the coefficients due to the advection term
        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

    end subroutine get_u_coefficients

    !> \brief Calculates the source of the linear system for u for the real volumes.
    subroutine get_u_source( nx, ny, beta, dt, rp, re, rn, ye, yk &
            ,             Jp, roe, ron, roa, p, Uce, Vcn, ua, u, cup, b ) ! The last two are output
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: beta   !< UDS/CDS mixing constant (0=UDS, 1=CDS)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: ye     !< y_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: yk     !< y_csi at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: roa    !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: p      !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: ua     !< u of the previous time step (m/s)
        real(8), dimension (nx*ny), intent(in) :: u      !< Cartesian velocity at the center of volume P (m/s)

        real(8), dimension (nx*ny), intent(out) :: cup   !< Term of deferred correction for u (N)
        real(8), dimension (nx*ny), intent(out) :: b     !< Source of the linear system for u (N)

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

        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

    end subroutine get_u_source

    !>  \brief Calculates the coefficients of the linear system for v  for the real volumes.
    subroutine get_v_coefficients( nx, ny, dt, rp, re, rn, Jp &
            ,                        roe, ron, roa, Uce, Vcn, a ) ! Last 1 is output
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: roa    !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)

        real(8), dimension (nx*ny,9), intent(out) :: a   !< Coefficients of the linear system for v (kg/s)

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

        ! Contribution to the coefficients due to the advection term
        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

    end subroutine get_v_coefficients

    !> \brief Calculates the source of the linear system for v for the real volumes.
    subroutine get_v_source( nx, ny, beta, dt, rp, re, rn, xe, xk &
            ,          Jp, roe, ron, roa, p, Uce, Vcn, va, v, cvp, b ) ! Last 3 are output
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: beta   !< UDS/CDS mixing constant (0=UDS, 1=CDS)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: xe     !< x_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: xk     !< x_csi at the center of north face  of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: roa    !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: p      !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: va     !< v of the previous time step (m/s)
        real(8), dimension (nx*ny), intent(in) :: v      !< Cartesian velocity at the center of volume P (m/s)

        real(8), dimension (nx*ny), intent(out) :: cvp   !< Term of deferred correction for v (N)
        real(8), dimension (nx*ny), intent(out) :: b     !< Source of the linear system for v (N)

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

        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

    end subroutine get_v_source


    !> \brief Calculates the coefficients and source of the linear system for T for the real volumes
    subroutine get_T_coefficients_and_source( nx, ny, beta, dt, rp, re, rn  &
            ,                                      xe, ye, xk, yk, Jp, roe, ron  &
            ,                      roa, p, pa, cp, Uce, Vcn, u, v, T, Ta, a, b ) ! Last 2 are output
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: beta   !< UDS/CDS mixing constant (0=UDS, 1=CDS)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: xe     !< x_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: ye     !< y_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: xk     !< x_csi at the center of north face  of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: yk     !< y_csi at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: roa    !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: p      !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(in) :: pa     !< Pressure of the previous time step (Pa)
        real(8), dimension (nx*ny), intent(in) :: cp     !< Specific heat at const pressure (J/(kg.K))
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(in) :: u      !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny), intent(in) :: v      !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny), intent(in) :: T      !< Temperature at center of volume P (K)
        real(8), dimension (nx*ny), intent(in) :: Ta     !< Temperature of the previous time step (K)

        real(8), dimension (nx*ny,9), intent(out) :: a   !< Coefficients of the linear system for T (J/(K.s))
        real(8), dimension (nx*ny),   intent(out) :: b   !< Source of the linear system for T (J/s)

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

        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

    end subroutine get_T_coefficients_and_source


    !> \brief 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 &
            ,                       roe, ron, g, de, dn, a) ! Output: last one
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: g      !< ro/p = 1/(Rg T) for ideal gases (kg/J)
        real(8), dimension (nx*ny), intent(in) :: de     !< SIMPLEC coefficients for Uce (m3.s/kg)
        real(8), dimension (nx*ny), intent(in) :: dn     !< SIMPLEC coefficients for Vcn (m3.s/kg)

        real(8), dimension (nx*ny,5), intent(out) :: a   !< Coefficients of the linear system for pl (m.s)

        integer :: i, j, np, nps, npn, npw, npe
        real(8) :: as, an, aw, ae
        real(8) :: roem, rowm, ronm, rosm

        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)  )

                roem = roe(np)   ! Density on east face
                rowm = roe(npw)  ! Density on west face
                ronm = ron(np)   ! Density on north face
                rosm = ron(nps)  ! Density on south face

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

                ! West
                a(np,2) = - re(npw) * ( (0.5d0 + aw) * Uce(npw) * g(npw) + rowm * 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) ) &
                    + roem * re(np)  * de(np)  &
                    + rowm * re(npw) * de(npw) &
                    + ronm * rn(np)  * dn(np)  &
                    + rosm * rn(nps) * dn(nps)

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

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

            end do
        end do

    end subroutine get_p_coefficients

    !> \brief Calculates the source of the linear system of the pressure deviation.
    !! 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, dt, rp, re, rn, Jp, roe, ron, roem, ronm, ro, roa, Ucem, Uce, Vcnm, Vcn, b) ! Output: last one
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny), intent(in) :: rp     !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: re     !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: rn     !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: Jp     !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny), intent(in) :: roe    !< Specific mass at east face of volume P (kg/m3) (incorrect roe*)
        real(8), dimension (nx*ny), intent(in) :: ron    !< Specific mass at north face of volume P (kg/m3) (incorrect ron*)
        real(8), dimension (nx*ny), intent(in) :: roem   !< Specific mass at east face of volume P (kg/m3) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: ronm   !< Specific mass at north face of volume P (kg/m3) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: ro     !< Specific mass (absolute density) at center of volume P (kg/m3) (incorrect ro*)
        real(8), dimension (nx*ny), intent(in) :: roa    !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny), intent(in) :: Ucem   !< Contravariant velocity U at east face of volume P (m2/s) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: Uce    !< Contravariant velocity U at east face of volume P (m2/s) (incorrect Uce*)
        real(8), dimension (nx*ny), intent(in) :: Vcnm   !< Contravariant velocity V at north face of volume P (m2/s) (previous iteraction)
        real(8), dimension (nx*ny), intent(in) :: Vcn    !< Contravariant velocity V at north face of volume P (m2/s) (incorrect Vcn*)

        real(8), dimension (nx*ny), intent(out) :: b     !< Source of the linear system for pl (kg/s)

        integer :: i, j, np, nps, npn, npw, npe
        real(8) :: as, an, aw, ae ! UDS coefficients

        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)  )

                b(np) = -rp(np) * ( ro(np) - roa(np) ) / ( Jp(np) * dt ) &
                    - roem(np)  * re(np)  * Uce(np)  &
                    + roem(npw) * re(npw) * Uce(npw) &
                    - ronm(np)  * rn(np)  * Vcn(np)  &
                    + ronm(nps) * rn(nps) * Vcn(nps) &
                    + ( roem(np)  - roe(np)  ) * re(np)  * Ucem(np ) &
                    - ( roem(npw) - roe(npw) ) * re(npw) * Ucem(npw) &
                    + ( ronm(np)  - ron(np)  ) * rn(np)  * Vcnm(np)  &
                    - ( ronm(nps) - ron(nps) ) * rn(nps) * Vcnm(nps)

            end do
        end do
    end subroutine get_p_source


    !> \brief Calculates the velocities at the interfaces between real volumes.
    !! The interpolation scheme applied here is the one appropriate for the SIMPLEC method.
    subroutine get_velocities_at_internal_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk     & ! Input
        ,                             xen, yen, xke, yke, Jp, cup, cvp           & ! 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 the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), intent(in) :: dt     !< Time step (s)
        real(8), dimension (nx*ny),   intent(in) :: rp   !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: re   !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: rn   !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: xe   !< x_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: ye   !< y_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: xk   !< x_csi at the center of north face  of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: yk   !< y_csi at the center of north face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: xen  !< x_eta at the center of north face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: yen  !< y_eta at the center of north face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: xke  !< x_csi at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: yke  !< y_csi at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in) :: Jp   !< Jacobian at the center of volume P (1/m2)
        real(8), dimension (nx*ny),   intent(in) :: cup  !< Term of deferred correction for u (N)
        real(8), dimension (nx*ny),   intent(in) :: cvp  !< Term of deferred correction for v (N)
        real(8), dimension (nx*ny,9), intent(in) :: au   !< Coefficients of the linear system for u (kg/s)
        real(8), dimension (nx*ny,9), intent(in) :: av   !< Coefficients of the linear system for v (kg/s)
        real(8), dimension (nx*ny),   intent(in) :: roa  !< Specific mass of the previous time step (kg/m3)
        real(8), dimension (nx*ny),   intent(in) :: p    !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny),   intent(in) :: u    !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny),   intent(in) :: v    !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny),   intent(in) :: uea  !< ue of the previous time step (m/s)
        real(8), dimension (nx*ny),   intent(in) :: vea  !< ve of the previous time step (m/s)
        real(8), dimension (nx*ny),   intent(in) :: una  !< un of the previous time step (m/s)
        real(8), dimension (nx*ny),   intent(in) :: vna  !< vn of the previous time step (m/s)
        real(8), dimension (nx*ny),   intent(out) :: ue  !< Cartesian velocity u at center of east face of volume P (m/s)
        real(8), dimension (nx*ny),   intent(out) :: ve  !< Cartesian velocity v at center of east face of volume P (m/s)
        real(8), dimension (nx*ny),   intent(out) :: un  !< Cartesian velocity u at center of north face of volume P (m/s)
        real(8), dimension (nx*ny),   intent(out) :: vn  !< Cartesian velocity v at center of north face of volume P (m/s)
        real(8), dimension (nx*ny),   intent(out) :: Uce !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny),   intent(out) :: Vcn !< Contravariant velocity V at north face of volume P (m2/s)

        ! 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) &
                    - sumup - sumue + 2.d0 * pue ) / ( au(np,5) + au(npe,5) )

                ve(np) = ( ( mpa + mea ) * vea(np) / dt + cvp(np) + cvp(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) &
                    - sumup - sumun + 2.d0 * pun ) / ( au(np,5) + au(npn,5) )

                vn(np) = ( ( mpa + mna ) * vna(np) / dt + cvp(np) + cvp(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_internal_faces

    !> \brief Calculates SIMPLEC coefficients at the interface between real volumes.
    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 the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny),   intent(in)  :: re  !< Radius at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)  :: rn  !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)  :: xe  !< x_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)  :: ye  !< y_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)  :: xk  !< x_csi at the center of north face  of volume P (m)
        real(8), dimension (nx*ny),   intent(in)  :: yk  !< y_csi at the center of north face of volume P (m)
        real(8), dimension (nx*ny,9), intent(in)  :: au  !< Coefficients of the linear system for u (kg/s)
        real(8), dimension (nx*ny,9), intent(in)  :: av  !< Coefficients of the linear system for v (kg/s)
        real(8), dimension (nx*ny),   intent(out) :: due !< SIMPLEC coefficients for ue (m2.s/kg)
        real(8), dimension (nx*ny),   intent(out) :: dve !< SIMPLEC coefficients for ve (m2.s/kg)
        real(8), dimension (nx*ny),   intent(out) :: dun !< SIMPLEC coefficients for un (m2.s/kg)
        real(8), dimension (nx*ny),   intent(out) :: dvn !< SIMPLEC coefficients for vn (m2.s/kg)
        real(8), dimension (nx*ny),   intent(out) :: de  !< SIMPLEC coefficients for Uce (m3.s/kg)
        real(8), dimension (nx*ny),   intent(out) :: dn  !< SIMPLEC coefficients for Vcn (m3.s/kg)

        ! 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


    !> \brief Corrects pressure and specific mass at all nodes (real+fictitious) using the
    !! pressure deviation pl.
    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 the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny), intent(in)    :: pl !< Pressure deviation at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(in)    :: g  !< ro/p = 1/(Rg T) for ideal gases (kg/J)
        real(8), dimension (nx*ny), intent(inout) :: p  !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(inout) :: ro !< Specific mass (absolute density) at center of volume P (kg/m3)

        p = p + pl

        ro = ro + pl * g

        !ro = p * g

    end subroutine get_pressure_density_correction_with_pl


    !> \brief Corrects u and v at all real nodes with the pressure deviation pl.
    !! Input (u*,v*) ->  Output (u,v).
    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 the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny),   intent(in)    :: xe  !< x_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)    :: ye  !< y_eta at the center of east face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)    :: xk  !< x_csi at the center of north face  of volume P (m)
        real(8), dimension (nx*ny),   intent(in)    :: yk  !< y_csi at the center of north face of volume P (m)
        real(8), dimension (nx*ny),   intent(in)    :: rp  !< Radius at the center of volume P (m)
        real(8), dimension (nx*ny),   intent(in)    :: pl  !< Pressure deviation at center o volume P (Pa)
        real(8), dimension (nx*ny,9), intent(in)    :: au  !< Coefficients of the linear system for u (kg/s)
        real(8), dimension (nx*ny,9), intent(in)    :: av  !< Coefficients of the linear system for v (kg/s)
        real(8), dimension (nx*ny),   intent(inout) :: u   !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny),   intent(inout) :: v   !< Cartesian velocity at the center of volume P (m/s)

        ! Auxiliary variables
        integer :: i, j, np, 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



    !> \brief Corrects the velocities ue, ve, un, vn, Uce, Vcn at all faces using the pressure
    !! deviation pl. Input (ue*, ve*, un*, vn*, Uce*, Vcn*) -> Output (ue, ve, un, vn, Uce, Vcn)
    subroutine get_velocity_correction_at_faces_with_pl( nx, ny     & ! Input
        ,                            due, dve, dun, dvn, de, dn, pl  & ! Input
        ,                            ue, ve, un, vn, Uce, Vcn )        ! InOutput
        implicit none
        integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny), intent(in)    :: due !< SIMPLEC coefficients for ue (m2.s/kg)
        real(8), dimension (nx*ny), intent(in)    :: dve !< SIMPLEC coefficients for ve (m2.s/kg)
        real(8), dimension (nx*ny), intent(in)    :: dun !< SIMPLEC coefficients for un (m2.s/kg)
        real(8), dimension (nx*ny), intent(in)    :: dvn !< SIMPLEC coefficients for vn (m2.s/kg)
        real(8), dimension (nx*ny), intent(in)    :: de  !< SIMPLEC coefficients for Uce (m3.s/kg)
        real(8), dimension (nx*ny), intent(in)    :: dn  !< SIMPLEC coefficients for Vcn (m3.s/kg)
        real(8), dimension (nx*ny), intent(in)    :: pl  !< Pressure deviation at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(inout) :: ue  !< Cartesian velocity u at center of east face of volume P (m/s)
        real(8), dimension (nx*ny), intent(inout) :: ve  !< Cartesian velocity v at center of east face of volume P (m/s)
        real(8), dimension (nx*ny), intent(inout) :: un  !< Cartesian velocity u at center of north face of volume P (m/s)
        real(8), dimension (nx*ny), intent(inout) :: vn  !< Cartesian velocity v at center of north face of volume P (m/s)
        real(8), dimension (nx*ny), intent(inout) :: Uce !< Contravariant velocity U at east face of volume P (m2/s)
        real(8), dimension (nx*ny), intent(inout) :: Vcn !< Contravariant velocity V at north face of volume P (m2/s)

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

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

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

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

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

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

            end do
        end do

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

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

                Vcn(np) = Vcn(np) + dn(np) * ( pl(np) - pl(npn) )

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

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

            end do
        end do

    end subroutine get_velocity_correction_at_faces_with_pl

end module coefficients
