!> \brief Contains subroutines for data post-processing.

module postp
    use newton_interpolation

    character (len = *), parameter, private :: text_editor     = "emacs -nw " !< Text editor name
    character (len = *), parameter, private :: graph_viewer    = "evince "    !< Graph viewer name
    character (len = *), parameter, private :: graph_generator = "gnuplot "   !< Plotter name

contains

    !> \brief Generates the post processed data
    subroutine post_processing(nx, ny, it, sem_a, sem_g, w_g, w_cam, lid &
            ,  sim_id, tcpu, RAM, res, lr, rb, x, y, xe, ye, xk, yk, Jp, Rg, cp &
            ,  gcp, Uce, Vcn, de, dn, pl, bp, xp, yp, u, v, p, T, ro, Cdfi)

        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)
        integer, intent(in) :: it   !< Iteractions counter for time evolution cycle
        integer, intent(in) :: sem_a !< 1 = do not open result files, 0 = open
        integer, intent(in) :: sem_g !< 0 = visualize the plot, 1 = do not visualize
        integer, intent(in) :: w_g   !< Frequency of writing data for graphics
        integer, intent(in) :: w_cam !< 1 = write the fields, 0 = do not
        integer, intent(in) :: lid   !< lid = listing file id
        character (len = *), intent(in) :: sim_id     !< Simulation identification
        real(8), intent(in) :: tcpu  !< Total cpu time (s)
        real(8), intent(in) :: RAM    !< RAM memory (MB)
        real(8), intent(in) :: res    !< Sum of all residuals of the linear systems for u, v, T and p.
        real(8), intent(in) :: lr    !< length of the body (m)
        real(8), intent(in) :: rb    !< base radius/semi-height of the body (m)
        real(8), dimension (nx*ny), intent(in) :: x   !< Coord. x at the northeast corner of the volume P (m)
        real(8), dimension (nx*ny), intent(in) :: y   !< Coord. y at the northeast corner of the 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), intent(in) :: Rg     !< Perfect gas constant (J/(kg.K))
        real(8), dimension (nx*ny), intent(in) :: cp    !< Specific heat at const pressure (J/(kg.K))
        real(8), dimension (nx*ny), intent(in) :: gcp   !< gcp = gamma = Cp/Cv at center of CV P (dimensionless)
        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) :: 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(in) :: bp    !< Source of the linear system for pl (kg/s)
        real(8), dimension (nx*ny), intent(inout) :: xp !< Coord. x of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(inout) :: yp !< Coord. y of the centroid of volume P (m)
        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)
        real(8), dimension (nx*ny), intent(inout) :: p  !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(inout) :: T  !< Temperature at center of volume P (K)
        real(8), dimension (nx*ny), intent(inout) :: ro !< Specific mass (absolute density) at center of volume P (kg/m3)
        real(8), intent(in) :: Cdfi !< Pressure foredrag coefficient (dimensionless)

        ! Auxiliary variables
        real(8), dimension (nx*ny) :: M    ! Mach number at the center of the volumes

        ! Writes grid parameters
        call write_grid_parameters( nx, ny, lid, x, y) ! No output

        ! Writes body parameters
        call write_body_parameters( lid, lr, rb) ! Output: last entry

        ! Writes grid boundary nodes
        call write_grid_boundary( nx, ny, w_g, sim_id, x, y)

        ! Plots grid boundaries
        call plotter04( sem_g, sim_id)

        ! Writes grid vertices
        call write_grid(nx, ny, w_g, sim_id, x, y)

        ! Plots the grid
        call plotter05( sem_g, sim_id)


        ! Calculates the boundary values and store them in the fictitious
        call post_processing_boundaries( nx, ny, x, y, xp, yp, u, v, p, T ) ! InOutput: last six entries

        call write_fields_on_boundaries( nx, ny, sim_id, xp, yp, u, v, p, T ) ! Output: none

        ! Generates the plots of u, v, T and p over the boundaries.
        call plotter07( sim_id)


        ! Calculates the Mach number field
        M = dsqrt ( (u**2 + v**2) / (gcp * Rg * T) )

        ! Calculates the specific mass field
        ro = p / ( Rg * T )

        ! Write fields

        !call write_main_fields(nx, ny, sim_id, xp, yp, p, ro, T, u, v, M)

        !call plotter06( sem_g, sim_id)

        if ( w_cam == 1 ) then

            write(lid,09)
            09     format(//,10x,'*** FIELDS ***', //, &
                10x,'The fictitious nodes have the boundary values.', //, &
                10x,'Horizon (from left to right): boundary west to east.', //, &
                10x,'Vertical (from top to bottom): boundary north to south.' )

            write(lid,10)
            10     format(//,30x,'M: Mach number (dimensionless)')
            call write_field(nx, ny, lid, xp, yp, M)

            write(lid,11)
            11     format(//,30x,'u: nodal cartesian velocity u (m/s)')
            call write_field(nx, ny, lid, xp, yp, u)

            write(lid,12)
            12     format(//,30x,'v: nodal cartesian velocity v (m/s)')
            call write_field(nx, ny, lid, xp, yp, v)

            write(lid,13)
            13     format(//,30x,'p: pressure (Pa)')
            call write_field(nx, ny, lid, xp, yp, p)

            write(lid,14)
            14     format(//,30x,'T: Temperature (K)')
            call write_field(nx, ny, lid, xp, yp, T)

            write(lid,15)
            15     format(//,30x,'ro: nodal specific mass (kg/m3)')
            call write_field(nx, ny, lid, xp, yp, ro)

            write(lid,55)
            55     format(//,30x,'cp: specific heat at constant pressure (J/kg.K)')
            call write_field(nx, ny, lid, xp, yp, cp)

            write(lid,58)
            58     format(//,30x,'gcp: specific heat ratio (dimensionless)')
            call write_field(nx, ny, lid, xp, yp, gcp)

            write(lid,16)
            16     format(//,30x,'Uce: contravariant velocity U at east face (m2/s)')
            call write_field(nx, ny, lid, xp, yp, Uce)

            write(lid,17)
            17     format(//,30x,'Vcn: contravariant velocity V at north face (m2/s)')
            call write_field(nx, ny, lid, xp, yp, Vcn)

            write(lid,18)
            18     format(//,30x,'pl: pressure variation (Pa)')
            call write_field(nx, ny, lid, xp, yp, pl)

            write(lid,19)
            19     format(//,30x,'bp: source of pl (kg/s)')
            call write_field(nx, ny, lid, xp, yp, bp)

            write(lid,21)
            21     format(//,30x,'X: x coordinate at northeast corner (m)')
            call write_field(nx, ny, lid, xp, yp, x)

            write(lid,22)
            22     format(//,30x,'Y: y coordinate at northeast corner (m)')
            call write_field(nx, ny, lid, xp, yp, y)

            write(lid,23)
            23     format(//,30x,'Jp: Jacobian at the center of volume (1/m2)')
            call write_field(nx, ny, lid, xp, yp, Jp)

            write(lid,24)
            24     format(//,30x,'xe: metric x eta at east face (m)')
            call write_field(nx, ny, lid, xp, yp, xe)

            write(lid,25)
            25     format(//,30x,'ye: metric y eta at east face (m)')
            call write_field(nx, ny, lid, xp, yp, ye)

            write(lid,26)
            26     format(//,30x,'xk: metric x ksi at north face (m)')
            call write_field(nx, ny, lid, xp, yp, xk)

            write(lid,27)
            27     format(//,30x,'yk: metric y ksi at north face (m)')
            call write_field(nx, ny, lid, xp, yp, yk)

            write(lid,28)
            28     format(//,30x,'de: coefficient d of SIMPLEC at east face (m3.s/kg)')
            call write_field(nx, ny, lid, xp, yp, de)

            write(lid,29)
            29     format(//,30x,'dn: coefficient d of SIMPLEC at north face (m3.s/kg)')
            call write_field(nx, ny, lid, xp, yp, dn)

        end if

        ! Print main results
        call write_main_results( nx, ny, lid, it, lr, tcpu, RAM, res, Cdfi, xp, u, v, p, T, ro)


        ! Plots the residual as a function of the iteractions
        call plotter02( sem_g, sim_id)

        ! show the listing file of iteractions
        if ( sem_a == 0 ) call system(text_editor // './mach2d_output/mach2d_' // trim(adjustl(sim_id)) // '_residual.dat')

        ! close the main listing file
        close (lid)

        ! show main listing file
        if ( sem_a == 0 ) call system(text_editor // './mach2d_output/mach2d_' // trim(adjustl(sim_id)) // '.txt')

    end subroutine post_processing


    !> \brief Plots the residual as a function of the iteractions.
    subroutine plotter02( sem_g, sim_id)
        implicit none
        integer, intent(in) :: sem_g !< 0 = visualize the plot, 1 = do not visualize
        character (len = *), intent(in) :: sim_id    !< Simulation identification

        ! Auxiliary variables
        character (len=250) :: str1, str2, str3

        !
        ! Generating gnuplot file
        !

        str1 = "./gnu/mach2d_SIM_ID_residual.gnu"

        str2 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) //"_residual.gnu"

        str3 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) //"_residual.eps"

        ! Replaces all SIM_ID string in the original file by sim_id
        call file_replace( str1, str2, "SIM_ID", trim(adjustl(sim_id)) )


        call system( graph_generator // trim(adjustl(str2)) )

        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) )

    end subroutine plotter02

    !> \brief Generates the boundary grid plot.
    subroutine plotter04( sem_g, sim_id)
        implicit none
        integer, intent(in) :: sem_g !< 0 = visualize the plot, 1 = do not visualize
        character (len = *), intent(in) :: sim_id    !< Simulation identification

        ! Auxiliary variables
        character (len=250) :: str1, str2, str3

        !
        ! Generating gnuplot file
        !

        ! Template file
        str1 = "./gnu/mach2d_SIM_ID_boundary.gnu"

        ! Case file
        str2 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // '_boundary.gnu'

        ! Output file
        str3 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // '_boundary.eps'

        ! Replaces all SIM_ID string in the original file by sim_id
        call file_replace( str1, str2, "SIM_ID", trim(adjustl(sim_id)) )

        call system( graph_generator // trim(adjustl(str2)) )

        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) )

    end subroutine plotter04

    !> \brief Generates the grid plot.
    subroutine plotter05( sem_g, sim_id)
        implicit none
        integer, intent(in) :: sem_g !< 0 = visualize the plot, 1 = do not visualize
        character (len = *), intent(in) :: sim_id    !< Simulation identification

        ! Auxiliary variables
        character (len=250) :: str1, str2, str3

        !
        ! Generating gnuplot file
        !

        ! Template file
        str1 = "./gnu/mach2d_SIM_ID_grid.gnu"

        ! Case file
        str2 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // '_grid.gnu'

        ! Output file
        str3 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // '_grid.eps'

        ! Replaces all SIM_ID string in the original file by sim_id
        call file_replace( str1, str2, "SIM_ID", trim(adjustl(sim_id)) )

        call system( graph_generator // trim(adjustl(str2)) )

        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) )

    end subroutine plotter05


    !> \brief Generates the plot of the following fields: p, ro, T, u, v, M.
    subroutine plotter06( sem_g, sim_id)
        implicit none
        integer, intent(in) :: sem_g !< 0 = visualize the plot, 1 = do not visualize
        character (len = *), intent(in) :: sim_id    !< Simulation identification

        ! Auxiliary variables
        character (len=250) :: str1, str2, str3

        !
        ! Generating gnuplot file
        !

        ! Template file
        str1 = "./gnu/mach2d_SIM_ID_main_fields.gnu"

        ! Case file
        str2 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) // '_main_fields.gnu'

        ! Replaces all SIM_ID string in the original file by sim_id
        call file_replace( str1, str2, "SIM_ID", trim(adjustl(sim_id)) )

        ! Auxiliary string
        str3 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id))

        call system( graph_generator // trim(adjustl(str2)) )

        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) // "_p_field.eps" )
        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) // "_ro_field.eps" )
        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) // "_T_field.eps" )
        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) // "_u_field.eps" )
        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) // "_v_field.eps" )
        if ( sem_g == 0 ) call system( graph_viewer // trim(adjustl(str3)) // "_M_field.eps" )

    end subroutine plotter06

    !> \brief Generates the plots of u, v, T and p over the boundaries.
    subroutine plotter07( sim_id)
        implicit none
        character (len = *), intent(in) :: sim_id    !< Simulation identification

        ! Auxiliary variables
        character (len=250) :: str1, str2

        !
        ! Generating gnuplot file
        !

        str1 = "./gnu/mach2d_SIM_ID_main_fields_boundaries.gnu"

        str2 = "./mach2d_output/mach2d_" // trim(adjustl(sim_id)) //"_main_fields_boundaries.gnu"

        ! Replaces all SIM_ID string in the original file by sim_id
        call file_replace( str1, str2, "SIM_ID", trim(adjustl(sim_id)) )

        call system( graph_generator // trim(adjustl(str2)) )

    end subroutine plotter07

    !> \brief Prints body parameters.
    subroutine write_body_parameters( unit, lr, rb) ! Output: last entry
        implicit none
        integer, intent(in) :: unit !< Unit to which the results will be printed
        real(8), intent(in) :: lr   !< length of the body (m)
        real(8), intent(in) :: rb   !< base radius/semi-height of the body (m)

        write(unit,*)
        write(unit,*)
        write(unit,*) '*** Body parameters ***'
        write(unit,*)
        write(unit,"(ES23.16,A)") lr, " = lr: Body length (m)"
        write(unit,"(ES23.16,A)") rb, " = rb: Body base radius (m)"
        write(unit,"(ES23.16,A)") lr/(2*rb), " = lr / (2 rb): Body length to base diameter ratio (dimensionless)"


    end subroutine write_body_parameters

    !> \brief Print grid parameters: min and max of dx and dy.
    subroutine write_grid_parameters( nx, ny, unit, x, y)
        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)
        integer, intent(in) :: unit !< Unit to which the results will be printed
        real(8), dimension (nx*ny), intent(in) :: x   !< Coord. x at the northeast corner of the volume P (m)
        real(8), dimension (nx*ny), intent(in) :: y   !< Coord. y at the northeast corner of the volume P (m)

        ! Auxiliary variables
        integer :: i, j, np
        real(8) :: dy, dx, dxmax, dxmin, dymax, dymin

        dxmin = maxval(abs(x))
        dxmax = minval(abs(x))
        dymin = maxval(abs(y))
        dymax = minval(abs(y))

        do j = 1, ny-1

            do i = 2, nx-1

                np = (j-1)*nx + i

                dx = dabs ( x(np) - x(np-1) )

                if ( dx < dxmin ) dxmin = dx

                if ( dx > dxmax ) dxmax = dx

            end do

        end do

        do i = 1, nx-1

            do j = 2, ny-1

                np = (j-1)*nx + i

                dy = dabs ( y(np) - y(np-nx) )

                if ( dy < dymin ) dymin = dy

                if ( dy > dymax ) dymax = dy

            end do

        end do

        write(unit,1) dxmin, dxmax, dxmin/dxmax, dymin, dymax, dymin/dymax

        1   format(//,'*** Grid parameters ***', //, &
            1pe22.9,' = dxmin: min(size of volume controls) in the x direction (m)',/, &
            1pe22.9,' = dxmax: max(size of volume controls) in the x direction (m)',/, &
            1pe22.9,' = dxmin / dxmax',/, &
            1pe22.9,' = dymin: min(size of volume controls) in the y direction (m)',/, &
            1pe22.9,' = dymax: max(size of volume controls) in the y direction (m)',/, &
            1pe22.9,' = dymin / dymax')

    end subroutine write_grid_parameters

    !> \brief Prints the grid boundaries.
    subroutine write_grid_boundary( nx, ny, w_g, sim_id, x, y)
        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)
        integer, intent(in) :: w_g   !< Frequency of writing data for graphics
        character (len = *), intent(in) :: sim_id     !< Simulation identification
        real(8), dimension (nx*ny), intent(in) :: x   !< Coord. x at the northeast corner of the volume P (m)
        real(8), dimension (nx*ny), intent(in) :: y   !< Coord. y at the northeast corner of the volume P (m)

        ! Auxiliary variables
        integer :: i, j, np

        ! writes south boundary

        open(7,file="./mach2d_output/mach2d_" // trim(adjustl(sim_id))//'_boundary_south.dat')

        write(7,1)

        j = 1
        do i = 1, nx-1
            np = (j-1)*nx + i
            if ( i==1 .or. i==nx-1 .or. mod(i,w_g)==0 ) &
                write(7,2) i, j, np, x(np), y(np)
        end do

        close(7)

        ! writes north boundary

        open(7,file="./mach2d_output/mach2d_" // trim(adjustl(sim_id))//'_boundary_north.dat')

        write(7,1)

        j = ny-1
        do i = 1, nx-1
            np = (j-1)*nx + i
            if ( i==1 .or. i==nx-1 .or. mod(i,w_g)==0 ) &
                write(7,2) i, j, np, x(np), y(np)
        end do

        close(7)

        ! writes west boundary

        open(7,file="./mach2d_output/mach2d_" // trim(adjustl(sim_id))//'_boundary_west.dat')

        write(7,1)

        i = 1
        do j = 1, ny-1
            np = (j-1)*nx + i
            if ( j==1 .or. j==ny-1 .or. mod(j,w_g)==0 ) &
                write(7,2) i, j, np, x(np), y(np)
        end do

        close(7)

        ! writes east boundary

        open(7,file="./mach2d_output/mach2d_" // trim(adjustl(sim_id))//'_boundary_east.dat')

        write(7,1)

        i = nx-1
        do j = 1, ny-1
            np = (j-1)*nx + i
            if ( j==1 .or. j==ny-1 .or. mod(j,w_g)==0 ) &
                write(7,2) i, j, np, x(np), y(np)
        end do

        close(7)

        1   format('#', t5,'i', t10,'j', t19,'np', t24,'x (m)', t48,'y (m)')

        2   format(2(i5),i10,2(1pe24.15))

    end subroutine write_grid_boundary

    !> \brief Prints the grid.
    subroutine write_grid(nx, ny, w_g, sim_id, x, y)
        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)
        integer, intent(in) :: w_g   !< Frequency of writing data for graphics
        character (len = *), intent(in) :: sim_id     !< Simulation identification
        real(8), dimension (nx*ny), intent(in) :: x   !< Coord. x at the northeast corner of the volume P (m)
        real(8), dimension (nx*ny), intent(in) :: y   !< Coord. y at the northeast corner of the volume P (m)

        ! Auxiliary variables
        integer :: i, j, np


        open(7,file= "./mach2d_output/mach2d_" // trim(adjustl(sim_id))//'_grid.dat')

        ! write the lines between west and east boundaries

        do j = 1, ny-1

            write(7,*)

            do i = 1, nx-1
                np = (j-1)*nx + i
                if ( i==1 .or. i==nx-1 .or. mod(i,w_g)==0 ) &
                    write(7,2) x(np), y(np)
            end do

        end do

        ! write the lines between south and north boundaries

        do i = 1, nx-1

            write(7,*)

            do j = 1, ny-1
                np = (j-1)*nx + i
                if ( j==1 .or. j==ny-1 .or. mod(j,w_g)==0 ) &
                    write(7,2) x(np), y(np)
            end do

        end do

        close(7)

        2   format(2(1pe24.15))

    end subroutine write_grid

    !> \brief Prints the main results.
    subroutine write_main_results( nx, ny, unit, it, lr, tcpu, RAM, res, Cdfi, xp, u, v, p, T, ro)
        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)
        integer, intent(in) :: unit !< Unit to which the results will be printed
        integer, intent(in) :: it   !< Iteractions counter for time evolution cycle
        real(8), intent(in) :: lr    !< length of the body (m)
        real(8), intent(in) :: tcpu  !< Total cpu time (s)
        real(8), intent(in) :: RAM    !< RAM memory (MB)
        real(8), intent(in) :: res    !< Sum of all residuals of the linear systems for u, v, T and p.
        real(8), intent(in) :: Cdfi   !< Pressure foredrag coefficient (dimensionless)
        real(8), dimension (nx*ny), intent(in) :: xp !< Coord. x of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: u  !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny), intent(in) :: v  !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny), intent(in) :: 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(in) :: ro !< Specific mass (absolute density) at center of volume P (kg/m3)

        ! Parameters

        integer, parameter :: ns = 20 ! Number points along the body surface where fields will be printed

        ! Auxiliary variables

        integer :: i
        integer :: resu, resv, resp, resT, resro

        real(8) :: dx
        real(8) :: xc(ns)
        real(8) :: uc, vc, pc, Tc, roc


        write(unit,*)
        write(unit,*)
        write(unit,*) "*** Main results ***"
        write(unit,*)
        write(unit,"(I23,A)")   nx-2,   " = nx-2:  Number of real volumes in the ksi direction"
        write(unit,"(I23,A)")   ny-2,   " = ny-2:  Number of real volumes in the eta direction"
        write(unit,"(I23,A)")     it,   " =   it:  Number of iteractions took by the time cycle"
        write(unit,"(ES23.16,A)") tcpu, " = tcpu:  Total cpu time (s)"
        write(unit,"(ES23.16,A)")  RAM, " =  RAM:  Memory (MB)"
        write(unit,"(ES23.16,A)")  res, " =  res:  Sum of all residuals of the linear systems for u, v, T and p."
        write(unit,*)
        write(unit,*) "*** Global variables ***"
        write(unit,*)
        write(unit,"(ES23.16,A)") Cdfi, " = Cdfi:  Pressure foredrag coefficient (dimensionless)"

        write(unit,*)
        write(unit,*) "*** Local variables over the body surface (quadratic interpolation) ***"
        write(unit,*)

        dx = lr / dble(ns)

        do i = 1, ns

            xc(i) = dx * ( dble(i) - 0.5d0 )

        end do

        write(unit,"(6(1X,A21))") "x", "u", "v", "p", "T", "ro"

        do i = 1, ns

            call get_ni_f0_vec_d2( nx, xp(1:nx), u(1:nx),  xc(i), uc,  resu)
            call get_ni_f0_vec_d2( nx, xp(1:nx), v(1:nx),  xc(i), vc,  resv)
            call get_ni_f0_vec_d2( nx, xp(1:nx), p(1:nx),  xc(i), pc,  resp)
            call get_ni_f0_vec_d2( nx, xp(1:nx), T(1:nx),  xc(i), Tc,  resT)
            call get_ni_f0_vec_d2( nx, xp(1:nx), ro(1:nx), xc(i), roc, resro)

            write(unit,"(6(1X,ES21.14))") xc(i), uc, vc, pc, Tc, roc

        end do

        write(unit,*)
        write(unit,*) "*** Local variables over the body surface (linear interpolation) ***"
        write(unit,*)

        write(unit,"(6(1X,A21))") "x", "u", "v", "p", "T", "ro"

        do i = 1, ns

            call get_ni_f0_vec_d1( nx, xp(1:nx), u(1:nx),  xc(i), uc,  resu)
            call get_ni_f0_vec_d1( nx, xp(1:nx), v(1:nx),  xc(i), vc,  resv)
            call get_ni_f0_vec_d1( nx, xp(1:nx), p(1:nx),  xc(i), pc,  resp)
            call get_ni_f0_vec_d1( nx, xp(1:nx), T(1:nx),  xc(i), Tc,  resT)
            call get_ni_f0_vec_d1( nx, xp(1:nx), ro(1:nx), xc(i), roc, resro)

            write(unit,"(6(1X,ES21.14))") xc(i), uc, vc, pc, Tc, roc

        end do

    end subroutine write_main_results

    !> \brief Prints some field.
    subroutine write_field(nx, ny, unit, xp, yp, b)
        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)
        integer, intent(in) :: unit !< Unit in which the field will be written
        real(8), dimension (nx*ny), intent(in) :: xp  !< Coord. x of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: yp  !< Coord. y of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: b   !< Field to be printed

        ! Auxiliary variables
        integer :: i, j, k, nt, maxi, i1, i2, np
        integer :: nloc(1)

        maxi = 10

        nt = nx / maxi

        do k = 1, nt

            i1 = (k-1) * maxi + 1

            i2 = k * maxi

            write(unit,"(A4,3X,10(I14))") "i = ", ( i, i = i1, i2 )

            do j = ny, 1, -1
                write(unit,"(A4,I3,10(ES14.6))") "j = ", j, (b(np), np = nx * (j-1) + i1, nx * (j-1) + i2)
            end do

            write(unit,*)

        end do

        i1 = nt * maxi + 1

        i2 = nx

        if ( i1 <= i2 ) then

            write(unit,"(A4,3X,10(I14))") "i = ", ( i, i = i1, i2 )

            do j = ny, 1, -1
                write(unit,"(A4,I3,10(ES14.6))") "j = ", j, (b(np), np = nx * (j-1) + i1, nx * (j-1) + i2)
            end do

        end if

        nloc = minloc(b)

        write(unit,*)
        write(unit,*) "Minimum of the field"
        write(unit,*)
        write(unit,"((ES14.6,A9,5X),2(ES14.6,A5,5X))") b(nloc(1)), " = Minval", xp(nloc(1)), " = xp", yp(nloc(1)), " = yp"


        nloc = maxloc(b)

        write(unit,*)
        write(unit,*) "Maximum of the field"
        write(unit,*)

        write(unit,"((ES14.6,A9,5X),2(ES14.6,A5,5X))") b(nloc(1)), " = Maxval", xp(nloc(1)), " = xp", yp(nloc(1)), " = yp"

        write(unit,*)
    end subroutine write_field


    !> \brief Prints the main fields.
    subroutine write_main_fields(nx, ny, sim_id, xp, yp, p, ro, T, u, v, M)
        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)
        character (len = 100), intent(in) :: sim_id !< Simulation identification
        real(8), dimension (nx*ny), intent(in) :: xp !< Coord. x of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: yp !< Coord. y of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: p  !< Pressure at center o volume P (Pa)
        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) :: T  !< Temperature at center of volume P (K)
        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) :: M  !< Mach number field

        ! Auxiliary variables

        integer :: i, j, np
        character (len = 100) :: str1

        str1 = "./mach2d_output/mach2d_"      &
            // trim(adjustl(sim_id)) // "_" &
            // "main_fields.dat"

        open(10, file = str1)

        write(10,"(8(X,A11))") "# x", "y", "p", "ro", "T", "u", "v", "M"

        do j = 1, ny
            do i = 1, nx
                np = (j-1)*nx + i
                write(10,"(8(X,ES11.4))") xp(np), yp(np), p(np), ro(np), T(np), u(np), v(np), M(np)
            end do
            write(10,*)
        end do

        close(10)

    end subroutine write_main_fields


    !> \brief Prints the fields u, v, p and T on the boundaries.
    subroutine write_fields_on_boundaries( nx, ny, sim_id, xp, yp, u, v, p, T ) ! Output: none
        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)
        character (len = 100), intent(in) :: sim_id !< Simulation identification
        real(8), dimension (nx*ny), intent(in) :: xp  !< Coord. x of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: yp  !< Coord. y of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: u   !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny), intent(in) :: v   !< Cartesian velocity at the center of volume P (m/s)
        real(8), dimension (nx*ny), intent(in) :: p   !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(in) :: T   !< Temperature at center of volume P (K)

        ! Auxiliary variables
        integer :: i, j, np
        character (len = 100) :: str1

        ! west boundary (frontal symmetry line)

        str1 = "./mach2d_output/mach2d_"      &
            // trim(adjustl(sim_id)) // "_" &
            // "main_fields_west_boundary.dat"

        open(10,file = str1)

        write(10,"(8(X,A23))") "# x", "y", "u", "v", "p", "T"

        i = 1
        do j = 1, ny

            np = (j-1)*nx + i

            write(10,"(6(X,ES23.16))") xp(np), yp(np), u(np), v(np), p(np), T(np)

        end do

        close(10)

        ! east boundary (exit)

        str1 = "./mach2d_output/mach2d_"      &
            // trim(adjustl(sim_id)) // "_" &
            // "main_fields_east_boundary.dat"

        open(10,file = str1)

        write(10,"(8(X,A23))") "# x", "y", "u", "v", "p", "T"

        i = nx

        do j = 1, ny

            np = (j-1)*nx + i

            write(10,"(6(X,ES23.16))") xp(np), yp(np), u(np), v(np), p(np), T(np)

        end do

        close(10)


        ! lower boundary (body wall)

        str1 = "./mach2d_output/mach2d_"      &
            // trim(adjustl(sim_id)) // "_" &
            // "main_fields_south_boundary.dat"

        open(10,file = str1)

        write(10,"(8(X,A23))") "# x", "y", "u", "v", "p", "T"

        j = 1

        do i = 1, nx

            np   = (j-1)*nx + i

            write(10,"(6(X,ES23.16))") xp(np), yp(np), u(np), v(np), p(np), T(np)

        end do

        close(10)

        ! upper boundary (far field)

        str1 = "./mach2d_output/mach2d_"      &
            // trim(adjustl(sim_id)) // "_" &
            // "main_fields_north_boundary.dat"

        open(10,file = str1)

        write(10,"(8(X,A23))") "# x", "y", "u", "v", "p", "T"

        j = ny

        do i = 1, nx

            np  = (j-1)*nx + i

            write(10,"(6(X,ES23.16))") xp(np), yp(np), u(np), v(np), p(np), T(np)

        end do

        close(10)

    end subroutine


    !> \brief Calculates the fields u, v, p and T over the boundaries and store them in the fictitious volumes.
    subroutine post_processing_boundaries( nx, ny, x, y, xp, yp, u, v, p, T ) ! InOutput: last six entries
        implicit none
        integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny), intent(in)  :: x     !< Coord. x at the northeast corner of the volume P (m)
        real(8), dimension (nx*ny), intent(in)  :: y     !< Coord. y at the northeast corner of the volume P (m)
        real(8), dimension (nx*ny), intent(inout) :: xp  !< Coord. x of the centroid of volume P (m)
        real(8), dimension (nx*ny), intent(inout) :: yp  !< Coord. y of the centroid of volume P (m)
        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)
        real(8), dimension (nx*ny), intent(inout) :: p   !< Pressure at center o volume P (Pa)
        real(8), dimension (nx*ny), intent(inout) :: T   !< Temperature at center of volume P (K)

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

        ! west boundary (frontal symmetry line)
        i = 1
        do j = 2, ny-1
            np = (j-1)*nx + i
            nps = np - nx
            npe  = np + 1
            u(np) = ( u(np) + u(npe) ) / 2.d0
            v(np) = ( v(np) + v(npe) ) / 2.d0
            p(np) = ( p(np) + p(npe) ) / 2.d0
            T(np) = ( T(np) + T(npe) ) / 2.d0
            xp(np) = ( x(nps) + x(np) ) / 2.d0
            yp(np) = ( y(nps) + y(np) ) / 2.d0
        end do

        ! east boundary (exit)
        i = nx
        do j = 2, ny-1
            np = (j-1)*nx + i
            npw  = np - 1
            u(np) = ( u(npw) + u(np) ) / 2
            v(np) = ( v(npw) + v(np) ) / 2
            p(np) = ( p(npw) + p(np) ) / 2
            T(np) = ( T(npw) + T(np) ) / 2
            xp(np) = ( x(npw-nx) + x(npw) ) / 2
            yp(np) = ( y(npw-nx) + y(npw) ) / 2
        end do

        ! lower boundary (body wall)
        j = 1
        do i = 2, nx-1
            np   = (j-1)*nx + i
            npw  = np - 1
            npn  = np + nx
            u(np) = ( u(npn) + u(np) ) / 2.d0
            v(np) = ( v(npn) + v(np) ) / 2.d0
            p(np) = ( p(npn) + p(np) ) / 2.d0
            T(np) = ( T(npn) + T(np) ) / 2.d0
            xp(np) = ( x(npw) + x(np) ) / 2.d0
            yp(np) = ( y(npw) + y(np) ) / 2.d0
        end do

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

        u(np) = ( u(npe) + u(npn) + u(npne) ) / 3.d0
        v(np) = ( v(npe) + v(npn) + v(npne) ) / 3.d0
        p(np) = ( p(npe) + p(npn) + p(npne) ) / 3.d0
        T(np) = ( T(npe) + T(npn) + T(npne) ) / 3.d0
        xp(np) = x(np)
        yp(np) = y(np)

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

        u(np) = ( u(npw) + u(npn) + u(npnw) ) / 3.d0
        v(np) = ( v(npw) + v(npn) + v(npnw) ) / 3.d0
        p(np) = ( p(npw) + p(npn) + p(npnw) ) / 3.d0
        T(np) = ( T(npw) + T(npn) + T(npnw) ) / 3.d0
        xp(np) = x(npw)
        yp(np) = y(npw)

        ! upper boundary (far field)
        j = ny
        do i = 2, nx-1
            np  = (j-1)*nx + i
            nps = np - nx
            u(np) = (u(np) + u(nps)) / 2
            v(np) = (v(np) + v(nps)) / 2
            p(np) = (p(np) + p(nps)) / 2
            T(np) = (T(np) + T(nps)) / 2
            xp(np) = ( x(nps) + x(nps-1) ) / 2
            yp(np) = ( y(nps) + y(nps-1) ) / 2
        end do

        ! NW corner
        j = ny
        i = 1
        np   = (j-1)*nx + i
        nps  = np - nx
        npe  = np + 1
        npse = nps + 1
        u(np) = ( u(npe) + u(nps) + u(npse) ) / 3.d0
        v(np) = ( v(npe) + v(nps) + v(npse) ) / 3.d0
        p(np) = ( p(npe) + p(nps) + p(npse) ) / 3.d0
        T(np) = ( T(npe) + T(nps) + T(npse) ) / 3.d0
        xp(np) = x(nps)
        yp(np) = y(nps)


        ! NE corner
        j = ny
        i = nx
        np   = (j-1)*nx + i
        nps  = np - nx
        npw  = np - 1
        npsw = nps - 1
        u(np) = ( u(npw) + u(nps) + u(npsw) ) / 3.d0
        v(np) = ( v(npw) + v(nps) + v(npsw) ) / 3.d0
        p(np) = ( p(npw) + p(nps) + p(npsw) ) / 3.d0
        T(np) = ( T(npw) + T(nps) + T(npsw) ) / 3.d0
        xp(np) = x(npsw)
        yp(np) = y(npsw)


    end subroutine post_processing_boundaries

    !> \brief Calculates the pressure foredrag coefficient.
    subroutine get_cdfi(nx, ny, coord, Rg, PF, TF, UF, rb, yk, rn, p, Cdfi)
        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)
        integer, intent(in) :: coord !< Kind of coord. system ( 1=cylindrical, 0 = cartesian)
        real(8), intent(in) :: Rg  !< Perfect gas constant (J/(kg.K))
        real(8), intent(in) :: PF  !< Free stream pressure (Pa)
        real(8), intent(in) :: TF  !< Free stream temperature (K)
        real(8), intent(in) :: UF  !< Free stream speed (m/s)
        real(8), intent(in) :: rb  !< base radius/semi-height of the body (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) :: rn !< Radius at the center of north face of volume P (m)
        real(8), dimension (nx*ny), intent(in) :: p  !< Pressure at center o volume P (Pa)
        real(8), intent(out) :: Cdfi !< Pressure foredrag coefficient (dimensionless)

        ! Auxiliary variables

        integer :: i, j, np, npn

        real(8) :: ROF, QF, pn


        ! Far field density
        ROF = PF / ( Rg * TF )

        ! Far field dynamic pressure
        QF = ROF * UF ** 2 / 2.d0

        Cdfi = 0.d0

        j = 1

        do i = 2, nx-1

            np   = nx * (j-1) + i

            npn  = np + nx

            pn = 0.5d0 * ( p(np) + p(npn) )

            Cdfi = Cdfi + ( pn - PF ) * rn(np) * yk(np)

        end do

        Cdfi = Cdfi * 2.d0 ** coord / ( rb ** (coord+1) * QF )

    end subroutine get_cdfi



    !> \brief Finds all the occurences of str1 in file1, replaces them by str2 and save the results in file2
    subroutine file_replace( file1, file2, str1, str2)
        implicit none
        character (len=*), intent(in) :: file1 !< Input file
        character (len=*), intent(in) :: file2 !< Output file
        character (len=*), intent(in) :: str1  !< Find
        character (len=*), intent(in) :: str2  !< Replace

        ! Auxiliary variables
        integer :: io
        character (len=1000) :: str

        if ( file1 /= file2 ) then

            open(10, file = file1)
            open(11, file = file2)

            do

                read(10,"(A)",iostat=io) str

                if ( io /= 0 ) exit

                call replace(str,str1,str2)

                write(11,*) trim(adjustl(str))

            end do

            close(10)
            close(11)

        else

            open(10, file = file1)
            open(11, file = trim(adjustl(file2))//".tmp")

            do

                read(10,"(A)",iostat=io) str

                if ( io /= 0 ) exit

                call replace(str,str1,str2)

                write(11,*) trim(adjustl(str))

            end do

            close(10)
            close(11)

            call rename( trim(adjustl(file2))//".tmp", file1 )

        end if

    end subroutine file_replace

    !> \brief Finds all the occurrences of str2 in str1 and replaces them by str3
    subroutine replace(str1, str2, str3)
        implicit none
        character (len=*), intent(inout) :: str1 !< Original string
        character (len=*), intent(in)    :: str2 !< Find
        character (len=*), intent(in)    :: str3 !< Replace

        ! Auxiliary variables

        integer :: idx !< First position where str2 was found
        integer :: n1  !< Length of str1
        integer :: n2  !< Length of str2

        n1 = len(str1)

        n2 = len(str2)

        do

            idx = index(str1,str2)

            if ( idx /= 0 ) then

                str1 = str1(1:idx-1) // str3 // str1(idx+n2:n1)

            else

                exit

            end if

        end do

    end subroutine replace


end module postp

