!===============================================================================
! User source terms definition.
!
! 1) Momentum equation (segegrated solver)
! 2) Momentum equation (coupled solver)
! 3) Species transport
! 4) 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-2012 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.

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

subroutine ustsnv &
!================

 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
   ivar   ,                                                       &
   icepdc , icetsm , itypsm ,                                     &
   dt     , rtpa   , propce , propfa , propfb ,                   &
   ckupdc , smacel ,                                              &
   crvexp , crvimp )

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

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

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

implicit none

! Arguments

integer          nvar   , nscal
integer          ncepdp , ncesmp
integer          ivar

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

double precision dt(ncelet), rtpa(ncelet,*)
double precision propce(ncelet,*)
double precision propfa(nfac,*), propfb(nfabor,*)
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
double precision crvexp(3,ncelet), crvimp(3,3,ncelet)

! Local variables

character*80     chaine
integer          iel, ipcrom, ipp
double precision ckp, qdm, retau, mu, delta, dpdz, rho

integer, allocatable, dimension(:) :: lstelt

!===============================================================================
! 1. Initialization
!===============================================================================

allocate(lstelt(ncel))

ipp    = ipprtp(ivar)

if (iwarni(ivar).ge.1) then
  chaine = nomvar(ipp)
  write(nfecra,1000) chaine(1:8)
endif

retau = 180.d0
delta = 1.d0

do iel = 1, ncel
  rho   = propce(iel,ipproc(irom)) ! expected value is 1.17862d0
  mu    = propce(iel,ipproc(iviscl)) ! expected value is 1.83337d-5
  dpdz  = ( 1/(rho*delta) )*( (retau*mu)/(delta) )**2
  crvexp(1,iel) =   volume(iel)*dpdz
enddo

!--------
! Formats
!--------

 1000 format(' User source termes for variable ',A8,/)

!----
! End
!----

deallocate(lstelt)

return
end subroutine ustsnv








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

subroutine ustssc &
!================

 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
   iscal  ,                                                       &
   icepdc , icetsm , itypsm ,                                     &
   dt     , rtpa   , rtp    , propce , propfa , propfb ,          &
   ckupdc , smacel ,                                              &
   crvexp , crvimp )

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

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

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

implicit none

! Arguments
integer          nvar   , nscal
integer          ncepdp , ncesmp
integer          iscal

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

double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
double precision propce(ncelet,*)
double precision propfa(nfac,*), propfb(nfabor,*)
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
double precision crvexp(ncelet), crvimp(ncelet)

! Local variables

character*80     chaine
integer          ivar, iiscvr, ipcrom, iel
integer          ilelt, nlelt
integer, allocatable, dimension(:) :: lstelt


! Local variables

integer          icel
double precision xx, masfl, perim, htf, dtdx, rho, u, source_sum

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

if (iscal.ne.iscalt) return
if (ntcabs.le.1) return
if (.false.) return

ivar = isca(iscal)

source_sum = 0.d0
masfl = 0
do icel = 1, ncel
  xx = xyzcen(1,icel)
  if (abs(xx).lt.0.1) then
    masfl = masfl + propfa( icel,ipprof( ifluma(iu) ) )
 endif
enddo

if (irangp.ge.0) then
  call parsom(masfl)
endif

! write(unit=nfecra, fmt=*) "masfl: ",masfl
! write(unit=nfecra, fmt=*) "cp: ",cp0

perim = 2.d0*(3.2d0) ! normally 2*(6.4+3.2)
htf   = -0.2d0
dtdx = (perim*htf)/(masfl*cp0) 

do iel = 1, ncel
   rho   = propce(iel,ipproc(irom))
!    rho   = 1
   u   = propce(iel,iu)
   crvexp(iel) = -volume(iel)*rho*u*cp0*dtdx
   source_sum = source_sum + crvexp(iel)
enddo

if (irangp.ge.0) then
  call parsom(source_sum)
endif

open(impusr(2),file='source.dat',status='old',position='append',action='write')

write(impusr(2),99) ttcabs, source_sum
close(impusr(2))

99    format(f14.1, g17.3)


return
end subroutine ustssc