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

! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2019 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.

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

!> \file cfmspr.f90
!> \brief Update the convective mass flux before the velocity prediction step.
!> It is the first step of the compressible algorithm at each time iteration.
!>
!> This function solves the continuity equation in pressure formulation and then
!> updates the density and the mass flux.
!>
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     nvar          total number of variables
!> \param[in]     nscal         total number of scalars
!> \param[in]     iterns        Navier-Stokes iteration number
!> \param[in]     ncepdp        number of cells with head loss
!> \param[in]     ncesmp        number of cells with mass source term
!> \param[in]     icepdc        index of cells with head loss
!> \param[in]     icetsm        index of cells with mass source term
!> \param[in]     itypsm        type of mass source term for each variable
!>                               (see uttsma.f90)
!> \param[in]     dt            time step (per cell)
!> \param[in]     vela          velocity value at time step beginning
!> \param[in]     ckupdc        work array for the head loss
!> \param[in]     smacel        variable value associated to the mass source
!>                               term (for ivar=ipr, smacel is the mass flux
!>                               \f$ \Gamma^n \f$)
!_______________________________________________________________________________

subroutine cfmspr &
 ( nvar   , nscal  , iterns , ncepdp , ncesmp ,                   &
   icepdc , icetsm , itypsm ,                                     &
   dt     , vela   ,                                              &
   ckupdc , smacel )

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

use paramx
use pointe, only:rvoid1
use numvar
use entsor
use optcal
use cstphy
use cstnum
use parall
use period
use ppppar
use ppthch
use ppincl
use mesh
use pointe, only: itypfb
use field
use cs_c_bindings
use cs_cf_bindings

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

implicit none

! Arguments

integer          nvar   , nscal, iterns
integer          ncepdp , ncesmp

integer          icepdc(ncepdp)
integer          icetsm(ncesmp), itypsm(ncesmp,nvar)

double precision dt(ncelet)
double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
double precision vela  (3  ,ncelet)

! Local variables

character(len=80) :: chaine
integer          ifac  , iel
integer          init  , inc   , iccocg, ii, jj
integer          iflmas, iflmab
integer          iphydp, icvflb
integer          nswrgp, imligp, iwarnp
integer          istatp, iconvp, idiffp, ndircp
integer          nswrsp, ircflp, ischcp, isstpp, iescap
integer          ivoid(1)
double precision epsrgp, climgp, extrap, blencp, epsilp
double precision epsrsp
double precision sclnor, hint

integer          imucpp, idftnp, iswdyp
integer          imvis1, f_id0
integer          f_id

double precision thetv, relaxp
double precision normp

double precision rvoid(1)

double precision, allocatable, dimension(:) :: viscf, viscb
double precision, allocatable, dimension(:) :: smbrs, rovsdt
double precision, allocatable, dimension(:) :: w1
double precision, allocatable, dimension(:) :: w7, w8, w9
double precision, allocatable, dimension(:) :: w10
double precision, allocatable, dimension(:) :: wflmas, wflmab, ivolfl, bvolfl
double precision, allocatable, dimension(:) :: wbfa, wbfb
double precision, allocatable, dimension(:) :: dpvar

double precision, allocatable, dimension(:) :: c2

double precision, dimension(:), pointer :: coefaf_p, coefbf_p
double precision, dimension(:), pointer :: imasfl, bmasfl
double precision, dimension(:), pointer :: rhopre, crom, brom
double precision, dimension(:), pointer :: cvar_pr, cvara_pr, cpro_cp, cpro_cv
double precision, dimension(:,:), pointer :: coefau
double precision, dimension(:,:,:), pointer :: coefbu
double precision, dimension(:), pointer :: cpro_divq
double precision, dimension(:), pointer :: cvar_fracv, cvar_fracm, cvar_frace

type(var_cal_opt) :: vcopt_p

!===============================================================================
!===============================================================================
! 1. INITIALISATION
!===============================================================================

! Allocate temporary arrays for the mass resolution
allocate(viscf(nfac), viscb(nfabor))
allocate(smbrs(ncelet), rovsdt(ncelet))
allocate(wflmas(nfac), wflmab(nfabor))
allocate(ivolfl(nfac), bvolfl(nfabor))

allocate(wbfa(nfabor), wbfb(nfabor))

! Allocate work arrays
allocate(w1(ncelet))
allocate(w7(ncelet), w8(ncelet), w9(ncelet))
allocate(w10(ncelet))
allocate(dpvar(ncelet))

call field_get_coefa_v(ivarfl(iu), coefau)
call field_get_coefb_v(ivarfl(iu), coefbu)

! --- Number of computational variable and post for pressure

f_id0 = -1

! --- Mass flux associated to energy
call field_get_key_int(ivarfl(isca(ienerg)), kimasf, iflmas)
call field_get_key_int(ivarfl(isca(ienerg)), kbmasf, iflmab)
call field_get_val_s(iflmas, imasfl)
call field_get_val_s(iflmab, bmasfl)

call field_get_val_s(icrom, crom)
call field_get_val_s(ibrom, brom)
call field_get_val_prev_s(icrom, rhopre)

call field_get_val_s(ivarfl(ipr), cvar_pr)
call field_get_val_prev_s(ivarfl(ipr), cvara_pr)


call field_get_label(ivarfl(ipr), chaine)

call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt_p)

if (ippmod(icompf).gt.1) then
  call field_get_val_s(ivarfl(isca(ifracv)), cvar_fracv)
  call field_get_val_s(ivarfl(isca(ifracm)), cvar_fracm)
  call field_get_val_s(ivarfl(isca(ifrace)), cvar_frace)
else
  cvar_fracv => rvoid1
  cvar_fracm => rvoid1
  cvar_frace => rvoid1
endif

if(vcopt_p%iwarni.ge.1) then
  write(nfecra,1000) chaine(1:8)
endif

call field_get_coefaf_s(ivarfl(ipr), coefaf_p)
call field_get_coefbf_s(ivarfl(ipr), coefbf_p)

if (icp.ge.0) then
  call field_get_val_s(icp, cpro_cp)
else
  cpro_cp => rvoid1
endif

if (icv.ge.0) then
  call field_get_val_s(icv, cpro_cv)
else
  cpro_cv => rvoid1
endif

! Computation of the boundary coefficients for the pressure gradient
! recontruction in accordance with the diffusion boundary coefficients
! (coefaf_p, coefbf_p)

do ifac = 1, nfabor
  iel = ifabor(ifac)
  hint = dt(iel) / distb(ifac)

  call set_neumann_scalar                 &
     ( wbfa(ifac)    , coefaf_p(ifac),    &
       wbfb(ifac)    , rvoid(1)      ,    &
       coefaf_p(ifac), hint )
enddo

!===============================================================================
! 2. SOURCE TERMS
!===============================================================================

! --> Initialization

do iel = 1, ncel
  smbrs(iel) = 0.d0
enddo
do iel = 1, ncel
  rovsdt(iel) = 0.d0
enddo


!     MASS SOURCE TERM
!     ================

if (ncesmp.gt.0) then
  do ii = 1, ncesmp
    iel = icetsm(ii)
    smbrs(iel) = smbrs(iel) + smacel(ii,ipr)*cell_f_vol(iel)
  enddo
endif


!     UNSTEADY TERM
!     =============

! --- Calculation of the square of sound velocity c2.
!     Pressure is an unsteady variable in this algorithm
!     Varpos has been modified for that.

allocate(c2(ncelet))

call cs_cf_thermo_c_square(cpro_cp, cpro_cv, cvar_pr, crom, &
                           cvar_fracv, cvar_fracm, cvar_frace, c2, ncel)

do iel = 1, ncel
  rovsdt(iel) = rovsdt(iel) + vcopt_p%istat*(cell_f_vol(iel)/(dt(iel)*c2(iel)))
enddo

!===============================================================================
! 3. "MASS FLUX" AND FACE "VISCOSITY" CALCULATION
!===============================================================================

! computation of the "convective flux" for the density

! volumic flux (u + dt f)
call cfmsfp                                                                    &
( nvar   , nscal  , iterns , ncepdp , ncesmp ,                                 &
  icepdc , icetsm , itypsm ,                                                   &
  dt     , vela   ,                                                            &
  ckupdc , smacel ,                                                            &
  ivolfl , bvolfl )

! mass flux at internal faces (upwind scheme for the density)
! (negative because added to RHS)
do ifac = 1, nfac
  ii = ifacel(1,ifac)
  jj = ifacel(2,ifac)
  wflmas(ifac) = -0.5d0*                                                       &
                 ( crom(ii)*(ivolfl(ifac)+abs(ivolfl(ifac)))                   &
                 + crom(jj)*(ivolfl(ifac)-abs(ivolfl(ifac))))
enddo

! mass flux at boundary faces
! (negative because added to RHS)
do ifac = 1, nfabor
  iel = ifabor(ifac)
  wflmab(ifac) = -brom(ifac)*bvolfl(ifac)
enddo

init = 0
call divmas(init,wflmas,wflmab,smbrs)

call field_get_id_try("predicted_vel_divergence", f_id)
if (f_id.ge.0) then
  call field_get_val_s(f_id, cpro_divq)
  do iel = 1, ncel
    cpro_divq(iel) = smbrs(iel)
  enddo
endif

! (Delta t)_ij is calculated as the "viscosity" associated to the pressure
imvis1 = 1

call viscfa                                                                    &
( imvis1 ,                                                                     &
  dt     ,                                                                     &
  viscf  , viscb  )

!===============================================================================
! 4. SOLVING
!===============================================================================

istatp = vcopt_p%istat
iconvp = vcopt_p%iconv
idiffp = vcopt_p%idiff
ndircp = vcopt_p%ndircl
nswrsp = vcopt_p%nswrsm
nswrgp = vcopt_p%nswrgr
imligp = vcopt_p%imligr
ircflp = vcopt_p%ircflu
ischcp = vcopt_p%ischcv
isstpp = vcopt_p%isstpc
iescap = 0
imucpp = 0
idftnp = vcopt_p%idften
iswdyp = vcopt_p%iswdyn
iwarnp = vcopt_p%iwarni
blencp = vcopt_p%blencv
epsilp = vcopt_p%epsilo
epsrsp = vcopt_p%epsrsm
epsrgp = vcopt_p%epsrgr
climgp = vcopt_p%climgr
extrap = vcopt_p%extrag
relaxp = vcopt_p%relaxv
thetv  = vcopt_p%thetav
icvflb = 0
normp = -1.d0

call codits                                                                    &
( idtvar , init   , ivarfl(ipr)     , iconvp , idiffp , ndircp ,               &
  imrgra , nswrsp , nswrgp , imligp , ircflp ,                                 &
  ischcp , isstpp , iescap , imucpp , idftnp , iswdyp ,                        &
  iwarnp , normp  ,                                                            &
  blencp , epsilp , epsrsp , epsrgp , climgp , extrap ,                        &
  relaxp , thetv  ,                                                            &
  cvara_pr        , cvara_pr        ,                                          &
  wbfa   , wbfb   ,                                                            &
  coefaf_p        , coefbf_p        ,                                          &
  wflmas          , wflmab          ,                                          &
  viscf  , viscb  , viscf  , viscb  , rvoid  ,                                 &
  rvoid  , rvoid  ,                                                            &
  icvflb , ivoid  ,                                                            &
  rovsdt , smbrs  , cvar_pr         , dpvar  ,                                 &
  rvoid  , rvoid  )

!===============================================================================
! 5. PRINTINGS AND CLIPPINGS
!===============================================================================

! --- User intervention for a finer management of the bounds and possible
!       corrective treatement.

call cs_cf_check_pressure(cvar_pr, ncel)

! --- Explicit balance (see codits : the increment is withdrawn)

if (vcopt_p%iwarni.ge.2) then
  do iel = 1, ncel
    smbrs(iel) = smbrs(iel)                                                    &
                 - vcopt_p%istat*(cell_f_vol(iel)/dt(iel))                     &
                   *(cvar_pr(iel)-cvara_pr(iel))                               &
                   * max(0,min(vcopt_p%nswrsm-2,1))
  enddo
  sclnor = sqrt(cs_gdot(ncel,smbrs,smbrs))
  write(nfecra,1200) chaine(1:8) ,sclnor
endif

!===============================================================================
! 6. COMMUNICATION OF P
!===============================================================================

if (irangp.ge.0.or.iperio.eq.1) then
  call synsca(cvar_pr)
endif

!===============================================================================
! 7. ACOUSTIC MASS FLUX CALCULATION AT THE FACES
!===============================================================================

! mass flux = [dt (grad P).n] + [rho (u + dt f)]

! Computation of [dt (grad P).n] by itrmas
init   = 1
inc    = 1
iccocg = 1
iphydp = 0
! festx,y,z    = rvoid
! viscf, viscb = arithmetic mean at faces
! viscelx,y,z  = dt
! This flux is stored as the mass flux of the energy

call itrmas                                                                    &
( f_id0  , init   , inc    , imrgra , iccocg , nswrgp , imligp ,               &
  iphydp , 0      , iwarnp ,                                                   &
  epsrgp , climgp , extrap ,                                                   &
  rvoid  ,                                                                     &
  cvar_pr,                                                                     &
  wbfa   , wbfb   ,                                                            &
  coefaf_p        , coefbf_p        ,                                          &
  viscf  , viscb  ,                                                            &
  dt     ,                                                                     &
  imasfl, bmasfl)

! Incrementation of the flux with [rho (u + dt f)].n = wflmas
! (added with a negative sign since wflmas,wflmab was used above
! in the right hand side).
do ifac = 1, nfac
  imasfl(ifac) = imasfl(ifac) - wflmas(ifac)
enddo
do ifac = 1, nfabor
  bmasfl(ifac) = bmasfl(ifac) - wflmab(ifac)
enddo

!===============================================================================
! 8. UPDATING OF THE DENSITY
!===============================================================================

if (igrdpp.gt.0) then

  do iel = 1, ncel
    ! Backup of the current density values
    rhopre(iel) = crom(iel)
    ! Update of density values
    crom(iel) = crom(iel)+(cvar_pr(iel)-cvara_pr(iel))/c2(iel)
  enddo

!===============================================================================
! 9. DENSITY COMMUNICATION
!===============================================================================

  if (irangp.ge.0.or.iperio.eq.1) then
    call synsca(crom)
    call synsca(rhopre)
  endif

endif

deallocate(c2)
deallocate(viscf, viscb)
deallocate(smbrs, rovsdt)
deallocate(wflmas, wflmab)
deallocate(ivolfl, bvolfl)
deallocate(w1)
deallocate(w7, w8, w9)
deallocate(w10)
deallocate(dpvar)
deallocate(wbfa, wbfb)

!--------
! FORMATS
!--------

 1000 format(/,                                                                 &
'   ** RESOLUTION FOR THE VARIABLE ',A8                        ,/,              &
'      ---------------------------                            ',/)
 1200 format(1X,A8,' : EXPLICIT BALANCE = ',E14.5)

!----
! END
!----

return
end subroutine
