!===============================================================================
! User source terms definition.
!
! 1) Momentum equation (coupled solver)
! 2) Species transport
! 3) Turbulence (k-epsilon, k-omega, Rij-epsilon, v2-f, Spalart-Allmaras)
!===============================================================================

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

!VERS

! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2017 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_source_terms.f90
!>
!> \brief Additional right-hand side source terms
!>
!> See \subpage cs_user_source_terms and \subpage cs_user_source_terms-scalar_in_a_channel
!> for examples.
!>
!> \brief Additional right-hand side source terms for velocity components equation
!> (Navier-Stokes)
!
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     nvar          total number of variables
!> \param[in]     nscal         total number of scalars
!> \param[in]     ncepdp        number of cells with head loss terms
!> \param[in]     ncesmp        number of cells with mass source terms
!> \param[in]     ivar          index number of the current variable
!> \param[in]     icepdc        index number of cells with head loss terms
!> \param[in]     icetsm        index number of cells with mass source terms
!> \param[in]     itypsm        type of mass source term for each variable
!>                               (see \ref cs_user_mass_source_terms)
!> \param[in]     dt            time step (per cell)
!> \param[in]     ckupdc        head loss coefficient
!> \param[in]     smacel        value associated to each variable in the mass
!>                               source terms or mass rate (see
!>                               \ref cs_user_mass_source_terms)
!> \param[out]    crvexp        explicit part of the source term
!> \param[out]    crvimp        implicit part of the source term
!_______________________________________________________________________________

subroutine ustsnv &
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
   ivar   ,                                                       &
   icepdc , icetsm , itypsm ,                                     &
   dt     ,                                                       &
   ckupdc , smacel ,                                              &
   crvexp , crvimp )

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

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

use paramx
use numvar
use entsor
use optcal
use cstphy
use parall
use period
use mesh
use field
use setup

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

implicit none

! Arguments

integer          nvar   , nscal
integer          ncepdp , ncesmp
integer          ivar

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

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

! Local variables

integer          iel

integer          ifac , iflmu
double precision err, surtot, mflow, ratio, sur
double precision bbuf(2)
double precision, save :: mflowa = 0.d0, dp = 0.d0

double precision, dimension(:), pointer :: imasfl

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

! Map field arrays

call field_get_key_int(ivarfl(ivar), kimasf, iflmu)
call field_get_val_s(iflmu, imasfl)

surtot = 0.d0
mflow  = 0.d0

do ifac = 1, nfac
  if(abs(cdgfac(1,ifac)-0.0d0).lt.0.0001d0) then
    surtot = surtot + surfan(ifac)
    mflow = mflow + imasfl(ifac)*surfac(1,ifac)/surfan(ifac)
  endif
enddo

if(irangp.ge.0)then
  call parsom(surtot)
  call parsom(mflow)
endif

write(nfecra,*) 'surtot', surtot

if (ntcabs.eq.1) then
  mflow  = 0.99*uref*surtot
  mflowa = 1.01*uref*surtot
else if (ntcabs.eq.ntpabs+1) then
  if(irangp.le.0) then
    open(file="MassFlow.dat",unit=impusr(1))
    read(impusr(1),*) mflowa, dp
    close(impusr(1))
    bbuf(1) = mflowa
    bbuf(2) = dp
  endif
  if(irangp.ge.0) then
    call parbcr(0, 2, bbuf)
    mflowa = bbuf(1)
    dp = bbuf(2)
  endif
endif

dp = dp+0.5d0*(2.d0*(mflow/surtot-1.d0)-(mflowa/surtot-1.d0))&
        /dtref

if(ntcabs.eq.ntmabs) then
  if(irangp.le.0) then
    open(file="MassFlow.dat",unit=impusr(1))
    write(impusr(1),*) mflow, dp
    close(impusr(1))
  endif
endif

voltot = 0.0d0
do iel=1,ncel
  voltot = voltot + volume(iel)
enddo
if (irangp.ge.0) then
  call parsom(voltot)
endif

write(nfecra,*) 'target:',1.d0*surtot,&
                         'actual:', mflow,&
                         'old:', mflowa

ratio = 0.8d0
err = (ratio*(1.d0*surtot-mflow) &
               + (1.d0-ratio)*(1.d0*surtot-mflowa) )/voltot

write(nfecra,*) 'erreur:' ,err

do iel=1,ncel
   !crvimp(1,1,iel) = err*volume(iel)
   crvexp(1,iel) = -volume(iel)*dp
enddo

mflowa = mflow

return
end subroutine ustsnv


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

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

!>    User subroutine.
!> \brief    Additional right-hand side source terms for scalar equations (user
!>     scalars and specific physics scalars).
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role
!______________________________________________________________________________!
!> \param[in]     nvar          total number of variables
!> \param[in]     nscal         total number of scalars
!> \param[in]     ncepdp        number of cells with head loss terms
!> \param[in]     ncesmp        number of cells with mass source terms
!> \param[in]     iscal         index number of the current scalar
!> \param[in]     icepdc        index number of cells with head loss terms
!> \param[in]     icetsm        index number of cells with mass source terms
!> \param[in]     itypsm        type of mass source term for each variable
!> \param[in]                    (see \ref cs_user_mass_source_terms)
!> \param[in]     dt            time step (per cell)
!> \param[in]     ckupdc        head loss coefficient
!> \param[in]     smacel        value associated to each variable in the mass
!> \param[in]                    source terms or mass rate (see \ref
!>                               cs_user_mass_source_terms)
!> \param[out]    crvexp        explicit part of the source term
!> \param[out]    crvimp        implicit part of the source term
!______________________________________________________________________________!


subroutine ustssc &
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
   iscal  ,                                                       &
   icepdc , icetsm , itypsm ,                                     &
   dt     ,                                                       &
   ckupdc , smacel ,                                              &
   crvexp , crvimp )

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

use paramx
use numvar
use entsor
use optcal
use cstphy
use parall
use period
use mesh
use field
use setup

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

implicit none

! Arguments

integer          nvar   , nscal
integer          ncepdp , ncesmp
integer          iscal

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

double precision dt(ncelet)
double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
double precision crvexp(ncelet), crvimp(ncelet)

! Local variables

integer          iel
integer          ifac, iflmas

double precision ubulk, surf, xx
double precision, dimension(:), pointer :: imasfl
double precision, dimension(:,:), pointer :: vel

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

if (iscal.ne.iscalt) return

ubulk = 0.d0
surf  = 0.d0

! Map field arrays
call field_get_val_v(ivarfl(iu), vel)

call field_get_key_int(ivarfl(iu), kimasf, iflmas)
call field_get_val_s(iflmas, imasfl)

do ifac = 1, nfac
  xx = cdgfac(1,ifac)
  if(abs(xx).lt.1.d-4) then
    ubulk = ubulk + imasfl(ifac)*surfac(1,ifac)/surfan(ifac)
    surf  = surf + surfac(1,ifac)
  endif
enddo
if(irangp.ge.0) call parsom(ubulk)
if(irangp.ge.0) call parsom(surf)

write(nfecra,*)'surtot ustssc:',surf

ubulk = abs(ubulk)/abs(surf)
ubulk = max(ubulk,1d-12)

write(nfecra,*) 'ustssc ubulk = ', ubulk
do iel = 1, ncel
  crvimp(iel) = 0.d0
  crvexp(iel) = -2.d0*volume(iel)*(1.d0-0.05d0/0.9d0)*qwall &
                  *vel(1,iel)/(0.9d0*ubulk)
enddo

return
end subroutine ustssc

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