!-------------------------------------------------------------------------------

!                      Code_Saturne version 3.0.1
!                      --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2013 EDF S.A.
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 2 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.

!-------------------------------------------------------------------------------

!===============================================================================
! Purpose:
! -------

!> \file cs_user_postprocess_var.f90
!> \brief Output additional variables on a postprocessing mesh.
!>
!> Several "automatic" postprocessing meshes may be defined:
!> - The volume mesh (ipart=-1) if 'ichrvl' = 1
!> - The boundary mesh (ipart=-2) if 'ichrbo' = 1
!> - SYRTHES coupling surface (ipart < -2) if 'ichrsy' = 1
!> - Cooling tower exchange zone meshes (ipart < -2) if 'ichrze' = 1
!>
!> Additional meshes (cells or faces) may also be defined through the GUI or
!> using the cs_user_postprocess_meshes() function from the
!> cs_user_postprocess.c file.
!>
!> This subroutine is called once for each post-processing mesh
!> (with a different value of 'ipart') for each time step at which output
!> on this mesh is active.
!
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     ipart         number of the post-processing mesh (< 0 or > 0)
!> \param[in]     nvar          total number of variables
!> \param[in]     nscal         total number of scalars
!> \param[in]     nvlsta        number of Lagrangian statistical variables
!> \param[in]     ncelps        number of cells in post-processing mesh
!> \param[in]     nfacps        number of interior faces in post-process. mesh
!> \param[in]     nfbrps        number of boundary faces in post-process. mesh
!> \param[in]     itypps        global presence flag (0 or 1) for cells (1),
!>                              interior faces (2), or boundary faces (3) in
!>                              post-processing mesh
!> \param[in]     lstcel        list of cells in post-processing mesh
!> \param[in]     lstfac        list of interior faces in post-processing mesh
!> \param[in]     lstfbr        list of boundary faces in post-processing mesh
!> \param[in]     dt            time step (per cell)
!> \param[in]     rtp, rtpa     calculated variables at cell centers
!>                               (at current and previous time steps)
!> \param[in]     propce        physical properties at cell centers
!> \param[in]     propfa        physical properties at interior face centers
!> \param[in]     propfb        physical properties at boundary face centers
!> \param[in]     statis        statistic values (Lagrangian)
!_______________________________________________________________________________

subroutine usvpst &
 ( ipart  ,                                                       &
   nvar   , nscal  , nvlsta ,                                     &
   ncelps , nfacps , nfbrps ,                                     &
   itypps ,                                                       &
   lstcel , lstfac , lstfbr ,                                     &
   dt     , rtpa   , rtp    , propce , propfa , propfb ,          &
   statis )

!===============================================================================

!===============================================================================

!===============================================================================
! Module files
!===============================================================================

use paramx
use cstnum
use pointe
use entsor
use optcal
use numvar
use parall
use period
use mesh
use field
use post
use user_module
!===============================================================================

implicit none

! Arguments

integer          ipart
integer          nvar,   nscal , nvlsta
integer          ncelps, nfacps, nfbrps

integer          itypps(3)
integer          lstcel(ncelps), lstfac(nfacps), lstfbr(nfbrps)

double precision dt(ncelet), rtpa(ncelet,*), rtp(ncelet,*)
double precision propce(ncelet,*)
double precision propfa(nfac,*), propfb(nfabor,*)
double precision statis(ncelet,nvlsta)

! Local variables

character*32     namevr

integer          ntindp
integer          iel, ifac, iloc, ivar
integer          idimt, ii , jj
logical          ientla, ivarpr
integer          imom1, imom2, ipcmo1, ipcmo2, idtcm
double precision pnd
double precision rvoid(1)

double precision, dimension(:), allocatable :: scel, sfac, sfbr
double precision, dimension(:,:), allocatable :: vcel, vfac, vfbr
double precision, dimension(:), pointer :: coefap, coefbp

integer          intpst
data             intpst /0/
save             intpst

!===============================================================================


!===============================================================================
! Increment call counter once per time step (possibly used in some tests)
!===============================================================================

if (ipart .eq. -1) then
  intpst = intpst + 1
endif

!===============================================================================
! 1. Handle variables to output
!    MUST BE FILLED IN by the user at indicated places
!===============================================================================

! The ipart argument matches a post-processing maehs id (using the EnSight
! vocabulary; the MED and CGNS equivalents are "mesh" and "base" respectively).
! The user will have defined post-processing meshes using the GUI or the
! cs_user_postprocess_meshes() function from the cs_user_postprocess.c
! file.

! This subroutine is called once for each post-processing mesh
! (with a different value of 'ipart') for each time step at which output
! on this mesh is active. For each mesh and for all variables we wish to
! post-process here, we must define certain parameters and pass them to
! the 'post_write_var' subroutine, which is in charge of the actual output.
! These parameters are:

! namevr <-- variable name
! idimt  <-- variable dimension
!            (1: scalar, 3: vector, 6: symmetric tensor, 9: tensor)
! ientla <-- when idimt >1, this flag specifies if the array containing the
!            variable values is interlaced when ientla = .true.
!            (x1, y1, z1, x2, y2, z2, x3, y3, z3...), or non-interlaced when
!            ientla = .false. (x1, x2, x3,...,y1, y2, y3,...,z1, z2, z3,...).
! ivarpr <-- specifies if the array containing the variable is defined on
!            the "parent" mesh or locally.
!            Even if the 'ipart' post-processing mesh contains all the
!            elements of its parent mesh, their numbering may be different,
!            especially when different element types are present.
!            A local array passed as an argument to 'post_write_var' is built
!            relative to the numbering of the 'ipart' post-processing mesh.
!            To post-process a variable contained for example in the 'user'
!            array, it should first be re-ordered, as shown here:
!              do iloc = 1, ncelps
!                iel = lstcel(iloc)
!                scel(iloc) = user(iel)
!              enddo
!            An alternative option is provided, to avoid unnecessary copies:
!            an array defined on the parent mesh, such our 'user' example,
!            may be passed directly to 'post_write_var', specifying that values
!            are defined on the parent mesh instead of the post-processing mesh,
!            by setting the 'ivarpr' argument of 'post_write_var' to .true..

! Note: be cautious with variable name lengths.

! We allow up to 32 characters here, but names may be truncted depending on the
! output format.

! The name length is not limited internally, so in case of 2 variables whoses
! names differ only after the truncation character, the corresponding names will
! both appear in the ".case" file; simply renaming one of the field descriptors
! in this text file will correct the output.

! Whitespace at the beginning or the end of a line is truncated automatically.
! Depending on the format used, prohibited characters (under EnSight, characters
! (  ) ] [ + - @           ! # * ^ $ / as well as white spaces and tabulations
! are automatically replaced by the _ character.

! Examples:

!   For post-processing mesh 2, we output the velocity, pressure, and prescribed
!   temperature at boundary faces (as well as 0 on possible interior faces)

!   For post-processing mesh 1, we output all the variables usually
!   post-processed, using a more compact coding.

!   Examples given here correspond to the meshes defined in
!   cs_user_postprocess.c

!===============================================================================
! Examples of volume variables on the main volume mesh (ipart = -1)
!===============================================================================

if (ipart .eq. -1) then

  ! Output of k=1/2(R11+R22+R33) for the Rij-epsilon model
  ! ------------------------------------------------------

  if (itytur .eq. 3) then

    allocate(scel(ncelps))

    do iloc = 1, ncelps
      iel = lstcel(iloc)
      scel(iloc) = 0.5d0*(  rtp(iel,ir11)  &
                          + rtp(iel,ir22)  &
                          + rtp(iel,ir33))
    enddo

    idimt = 1        ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
    ientla = .true.  ! dimension 1 here, so no effect
    ivarpr = .false. ! defined on the work array, not on the parent

    ! Output values; as we have no face values, we can pass a
    ! trivial array rvoid for those.
    call post_write_var(ipart, 'Turb energy', idimt, ientla, ivarpr,  &
                        ntcabs, ttcabs, scel, rvoid, rvoid)

    deallocate(scel)

  endif


  ! Output of a combination of moments
  ! ----------------------------------

  ! We assume in this example that we have 2 temporal means (moments):
  !   <u>  for imom=1
  !   <uu> for imom=2
  ! We seek to plot <u'u'>=<uu>-<U>**2

!===============================================================================
! Examples of volume variables on the boundary mesh (ipart = -2)
!===============================================================================

else if (ipart .eq. -2) then

  ! Output of the density at the boundary
  ! -------------------------------------

  idimt = 1        ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
  ientla = .true.  ! dimension 1 here, so no effect
  ivarpr = .true.  ! we use the propfb array defined on the parent mesh

  ! Output values; as we have no cell or interior face values, we can pass a
  ! trivial array for those.
  !call post_write_var(ipart, 'Density at boundary', idimt, ientla, ivarpr,    &
  !                    ntcabs, ttcabs, rvoid, rvoid, propfb(1,ipprob(irom)))
  call post_write_var(ipart, 'Wall_Shear_Stress', idimt, ientla, ivarpr,    &
                      ntcabs, ttcabs, rvoid, rvoid, usertau)
  call post_write_var(ipart, 'Friction_Velcoity', idimt, ientla, ivarpr,    &
                      ntcabs, ttcabs, rvoid, rvoid, userutau)


!===============================================================================
! Examples of volume variables on user meshes 1 or 2
!===============================================================================

else if (ipart.eq.1 .or. ipart.eq.2) then

  ! Output of the velocity
  ! ----------------------

  ! Compute variable values on interior faces.
  ! In this example, we use a simple linear interpolation.
  ! For parallel calculations, if neighbors are used, they must be synchronized
  ! first. This also applies for periodicity.

  if (irangp.ge.0.or.iperio.eq.1) then
    call synvec(rtp(1,iu), rtp(1,iv), rtp(1,iw))
    !==========
  endif

  allocate(vfac(3,nfacps), vfbr(3,nfbrps))

  do iloc = 1, nfacps

    ifac = lstfac(iloc)
    ii = ifacel(1, ifac)
    jj = ifacel(2, ifac)
    pnd = pond(ifac)

    vfac(1,iloc) = pnd  * rtp(ii,iu) + (1.d0 - pnd) * rtp(jj,iu)
    vfac(2,iloc) = pnd  * rtp(ii,iv) + (1.d0 - pnd) * rtp(jj,iv)
    vfac(3,iloc) = pnd  * rtp(ii,iw) + (1.d0 - pnd) * rtp(jj,iw)

  enddo

  ! Compute variable values on boundary faces.
  ! In this example, we use a simple copy of the adjacent cell value.

  do iloc = 1, nfbrps

    ifac = lstfbr(iloc)
    ii = ifabor(ifac)

    vfbr(1,iloc) = rtp(ii, iu)
    vfbr(2,iloc) = rtp(ii, iv)
    vfbr(3,iloc) = rtp(ii, iw)

  enddo

  idimt = 3        ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
  ientla = .true.  ! interleaved
  ivarpr = .false. ! defined on the work array, not on the parent

  ! Output values; as we have no cell values, we can pass a
  ! trivial array for those.
  call post_write_var(ipart, 'Interpolated velocity', idimt, ientla, ivarpr,  &
                      ntcabs, ttcabs, rvoid, vfac, vfbr)

  deallocate(vfac, vfbr)

  ! Output of the pressure
  ! ----------------------

  ! Variable number
  ivar = ipr

  ! Compute variable values on interior faces.
  ! In this example, we use a simple linear interpolation.
  ! For parallel calculations, if neighbors are used, they must be synchronized
  ! first. This also applies for periodicity.

  if (irangp.ge.0.or.iperio.eq.1) then
    call synsca(rtp(1,ivar))
    !==========
  endif

  allocate(sfac(nfacps), sfbr(nfbrps))

  do iloc = 1, nfacps

    ifac = lstfac(iloc)
    ii = ifacel(1, ifac)
    jj = ifacel(2, ifac)
    pnd = pond(ifac)

    sfac(iloc) =           pnd  * rtp(ii, ivar)  &
                 + (1.d0 - pnd) * rtp(jj, ivar)

  enddo

  ! Compute variable values on boundary faces.
  ! In this example, we use a simple copy of the adjacent cell value.

  do iloc = 1, nfbrps

    ifac = lstfbr(iloc)
    ii = ifabor(ifac)

    sfbr(iloc) = rtp(ii, ivar)

  enddo

  idimt = 1        ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
  ientla = .true.  ! dimension 1 here, so no effect
  ivarpr = .false. ! defined on the work array, not on the parent

  ! Output values; as we have no cell values, we can pass a
  ! trivial array for those.
  call post_write_var(ipart, 'Interpolated pressure', idimt, ientla, ivarpr,  &
                      ntcabs, ttcabs, rvoid, sfac, sfbr)

  deallocate(sfac, sfbr)

  ! The examples below illustrate how to output a same variable in different
  ! ways (interlaced or not, using an indirection or not).


  ! Output of the centers of gravity, interlaced
  ! --------------------------------

  if (intpst.eq.1) then

    allocate(vfac(3,nfacps), vfbr(3,nfbrps))

    do iloc = 1, nfacps

      ifac = lstfac(iloc)

      vfac(1,iloc) = cdgfac(1, ifac)
      vfac(2,iloc) = cdgfac(2, ifac)
      vfac(3,iloc) = cdgfac(3, ifac)

    enddo

    ! Compute variable values on boundary faces

    do iloc = 1, nfbrps

      ifac = lstfbr(iloc)

      vfbr(1, iloc) = cdgfbo(1, ifac)
      vfbr(2, iloc) = cdgfbo(2, ifac)
      vfbr(3, iloc) = cdgfbo(3, ifac)

    enddo

    ! We assign a negative time step and output this variable once only
    ! to avoid duplicating it at each output time (assuming a fixed mesh).
    ntindp = -1

    idimt = 3        ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
    ientla = .true.  ! interleaved
    ivarpr = .false. ! defined on the work array, not on the parent

    ! Output values; as we have no cell values, we can pass a
    ! trivial array for those.
    call post_write_var(ipart, 'face cog (interlaced)', idimt,               &
                        ientla, ivarpr,                                      &
                        ntindp, ttcabs, rvoid, vfac, vfbr)

    deallocate(vfac, vfbr)

  endif

  ! Output of the centers of gravity, non-interlaced, time independent
  ! --------------------------------

  if (intpst.eq.1) then

    allocate(vfac(nfacps, 3), vfbr(nfbrps, 3))

    do iloc = 1, nfacps

      ifac = lstfac(iloc)

      vfac(iloc,1) = cdgfac(1, ifac)
      vfac(iloc,2) = cdgfac(2, ifac)
      vfac(iloc,3) = cdgfac(3, ifac)

    enddo

    ! Compute variable values on boundary faces

    do iloc = 1, nfbrps

      ifac = lstfbr(iloc)

      vfbr(iloc,1) = cdgfbo(1, ifac)
      vfbr(iloc,2) = cdgfbo(2, ifac)
      vfbr(iloc,3) = cdgfbo(3, ifac)

    enddo

    ! We assign a negative time step and output this variable once only
    ! to avoid duplicating it at each output time (assuming a fixed mesh).
    ntindp = -1

    idimt = 3         ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
    ientla = .false.  ! not interleaved
    ivarpr = .false.  ! defined on the work array, not on the parent

    ! Output values; as we have no cell values, we can pass a
    ! trivial array for those.
    call post_write_var(ipart, 'face cog (non interlaced)', idimt,           &
                        ientla, ivarpr,                                      &
                        ntindp, ttcabs, rvoid, vfac, vfbr)

    deallocate(vfac, vfbr)

  endif

  ! Output of the centers of gravity, with indirection (parent-based)
  ! --------------------------------

  if (intpst.eq.1) then

    ! We assign a negative time step and output this variable once only
    ! to avoid duplicating it at each output time (assuming a fixed mesh).
    ntindp = -1

    idimt = 3        ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
    ientla = .true.  ! interleaved
    ivarpr = .true.  ! defined on the parent

    ! Output values; as we have no cell values, we can pass a
    ! trivial array for those.
    call post_write_var(ipart, 'face cog (parent)', idimt, ientla, ivarpr,   &
                        ntindp, ttcabs, rvoid, cdgfac, cdgfbo)

  endif

!===============================================================================
! Examples of volume variables on user meshes 3 or 4
!===============================================================================

else if (ipart.ge.3 .and. ipart.le.4) then

  ! Output of the velocity
  ! ----------------------

  ! Compute variable values on interior faces.
  ! In this example, we use a simple linear interpolation.
  ! For parallel calculations, if neighbors are used, they must be synchronized
  ! first. This also applies for periodicity.

  if (irangp.ge.0.or.iperio.eq.1) then
    call synvec(rtp(1,iu), rtp(1,iv), rtp(1,iw))
    !==========
  endif

  allocate(vfac(3,nfacps), vfbr(3,nfbrps))

  do iloc = 1, nfacps

    ifac = lstfac(iloc)
    ii = ifacel(1, ifac)
    jj = ifacel(2, ifac)
    pnd = pond(ifac)

    vfac(1,iloc) =            pnd  * rtp(ii, iu)   &
                    + (1.d0 - pnd) * rtp(jj, iu)
    vfac(2,iloc) =            pnd  * rtp(ii, iv)   &
                    + (1.d0 - pnd) * rtp(jj, iv)
    vfac(3,iloc) =            pnd  * rtp(ii, iw)   &
                    + (1.d0 - pnd) * rtp(jj, iw)

  enddo

  ! Compute variable values on boundary faces.
  ! In this example, we use a simple copy of the adjacent cell value.

  do iloc = 1, nfbrps

    ifac = lstfbr(iloc)
    ii = ifabor(ifac)

    vfbr(1,iloc) = rtp(ii, iu)
    vfbr(2,iloc) = rtp(ii, iv)
    vfbr(3,iloc) = rtp(ii, iw)

  enddo

  idimt = 3         ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
  ientla = .true.   ! interleaved
  ivarpr = .false.  ! defined on the work array

  ! Output values; as we have no cell values, we can pass a
  ! trivial array for those.
  call post_write_var(ipart, 'Velocity', idimt, ientla, ivarpr,              &
                      ntcabs, ttcabs, rvoid, vfac, vfbr)

  deallocate(vfac, vfbr)

  ! Output of the pressure
  ! ----------------------

  ! Variable number
  ivar = ipr

  ! Compute variable values on interior faces.
  ! In this example, we use a simple linear interpolation.
  ! For parallel calculations, if neighbors are used, they must be synchronized
  ! first. This also applies for periodicity.

  if (irangp.ge.0.or.iperio.eq.1) then
    call synsca(rtp(1,ivar))
    !==========
  endif

  allocate(sfac(nfacps), sfbr(nfbrps))

  do iloc = 1, nfacps

    ifac = lstfac(iloc)
    ii = ifacel(1, ifac)
    jj = ifacel(2, ifac)
    pnd = pond(ifac)

    sfac(iloc)  =           pnd  * rtp(ii, ivar)   &
                  + (1.d0 - pnd) * rtp(jj, ivar)

  enddo

  ! Compute variable values on boundary faces.
  ! In this example, we use a simple copy of the adjacent cell value.

  do iloc = 1, nfbrps

    ifac = lstfbr(iloc)
    ii = ifabor(ifac)

    sfbr(iloc) = rtp(ii, ivar)

  enddo

  idimt = 1         ! 1: scalar, 3: vector, 6/9: symm/non-symm tensor
  ientla = .true.   ! interleaved
  ivarpr = .false.  ! defined on the work array

  ! Output values; as we have no cell values, we can pass a
  ! trivial array for those.
  call post_write_var(ipart, 'Pressure', idimt, ientla, ivarpr,              &
                      ntcabs, ttcabs, rvoid, sfac, sfbr)

  deallocate(sfac, sfbr)

endif ! end of test on post-processing mesh number

return

end subroutine usvpst