!===============================================================================
! 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     ,                                                       &
   ckupdc , smacel ,                                              &
   crvexp , crvimp )

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

!    User subroutine.

!    Additional right-hand side source terms for velocity components equation
!    (Navier-Stokes)

!
! Usage
! -----
!
! The additional source term is decomposed into an explicit part (crvexp) and
! an implicit part (crvimp) that must be provided here.
! The resulting equation solved by the code for a velocity is:
!
!  rho*volume*du/dt + .... = crvimp*u + crvexp
!
! Note that crvexp and crvimp are defined after the Finite Volume integration
! over the cells, so they include the "volume" term. More precisely:
!   - crvexp is expressed in kg.m/s2
!   - crvimp is expressed in kg/s
!
! The crvexp and crvimp arrays are already initialized to 0 before entering the
! the routine. It is not needed to do it in the routine (waste of CPU time).
!
! For stability reasons, Code_Saturne will not add -crvimp directly to the
! diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is
! treated implicitely only if it strengthens the diagonal of the matrix.
! However, when using the second-order in time scheme, this limitation cannot
! be done anymore and -crvimp is added directly. The user should therefore test
! the negativity of crvimp by himself.
!
! When using the second-order in time scheme, one should supply:
!   - crvexp at time n
!   - crvimp at time n+1/2
!
!
! The selection of cells where to apply the source terms is based on a getcel
! command. For more info on the syntax of the getcel command, refer to the
! user manual or to the comments on the similar command getfbr in the routine
! cs_user_boundary_conditions.

!-------------------------------------------------------------------------------
! Arguments
!__________________.____._____.________________________________________________.
! name             !type!mode ! role                                           !
!__________________!____!_____!________________________________________________!
! nvar             ! i  ! <-- ! total number of variables                      !
! nscal            ! i  ! <-- ! total number of scalars                        !
! ncepdp           ! i  ! <-- ! number of cells with head loss terms           !
! ncssmp           ! i  ! <-- ! number of cells with mass source terms         !
! ivar             ! i  ! <-- ! index number of the current variable           !
! icepdc(ncepdp)   ! ia ! <-- ! index number of cells with head loss terms     !
! icetsm(ncesmp)   ! ia ! <-- ! index number of cells with mass source terms   !
! itypsm           ! ia ! <-- ! type of mass source term for each variable     !
!  (ncesmp,nvar)   !    !     !  (see ustsma)                                  !
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
! ckupdc(ncepdp,6) ! ra ! <-- ! head loss coefficient                          !
! smacel           ! ra ! <-- ! value associated to each variable in the mass  !
!  (ncesmp,nvar)   !    !     !  source terms or mass rate (see ustsma)        !
! crvexp           ! ra ! --> ! explicit part of the source term               !
! crvimp           ! ra ! --> ! implicit part of the source term               !
!__________________!____!_____!________________________________________________!

!     Type: i (integer), r (real), s (string), a (array), l (logical),
!           and composite types (ex: ra real array)
!     mode: <-- input, --> output, <-> modifies data, --- work array
!===============================================================================

!===============================================================================
! 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(ncepdp,6), smacel(ncesmp,nvar)
double precision crvexp(3,ncelet), crvimp(3,3,ncelet)

! Local variables

integer          iel
integer          ifac , iflmu, nlfac
double precision surtot, mflow, mflowa, dp

double precision, dimension(:), pointer :: imasfl
integer, allocatable, dimension(:) :: lstfac
!===============================================================================

! 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

!calculate surface at mass flow at the inlet
allocate (lstfac(nfac))
call getfac('2',nlfac,lstfac)

do iel = 1, nlfac
  ifac = lstfac(iel)
  surtot = surtot + surfan(ifac)
  mflow = mflow + imasfl(ifac)*surfac(2,ifac)/surfan(ifac)
enddo

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

write(nfecra,*) 'surtot', surtot

if((ntcabs-ntpabs).eq.1) then
  mflow  = 0.99*ro0*uref*surtot
  mflowa = 1.01*ro0*uref*surtot
  dp = 0.d0
else
  open(file='MassFlow.dat',unit=impusr(1))
  read(impusr(1),*) mflowa,dp
  close(impusr(1))
endif

!calculate the source term needed to impose the target mass flow rate
dp = dp+0.95d0*(mflow/surtot-ro0*uref)/dtref
!dp =dp+(2.0d0*(mflow/surtot-ro0*uref)-(mflowa/surtot-ro0*uref))&
!               /(2.0d0*dtref)


write(nfecra,*) 'target:',ro0*uref*surtot,&
                         'actual:', mflow,&
                         'previous:', mflowa,&
                         'dp:', dp


!store the calculated values of dp and mflow for the next time step
if(irangp.le.0) then
  if((ntcabs-ntpabs).eq.1) then
     open(file="MassFlow.dat",unit=impusr(1),status='new')
  else
     open(file="MassFlow.dat",unit=impusr(1))
  endif
  write(impusr(1),*) mflow, dp
  close(impusr(1))
endif

!write mass flow rate and dp at every time step
open(file='mflow_source.dat',unit=impusr(3))
if((ntcabs-ntpabs).eq.1) then
  write(impusr(3),*) 'iter ','mflow ','dp'
endif
write(impusr(3),*) ntcabs-ntpabs,mflow,dp
if(ntcabs.eq.ntmabs) then
  close(impusr(3))
endif

!impose source term
do iel=1,ncel
   crvexp(2,iel) = -volume(iel)*dp
enddo

deallocate(lstfac)

return
end subroutine ustsnv


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


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

!    User subroutine.

!    Additional right-hand side source terms for scalar equations (user
!     scalars and specific physics scalars).

!
! Usage
! -----
! The routine is called for each scalar, user or specific physisc. It is
! therefore necessary to test the value of the scalar number iscal to separate
! the treatments of the different scalars (if (iscal.eq.p) then ....).
!
! The additional source term is decomposed into an explicit part (crvexp) and
! an implicit part (crvimp) that must be provided here.
! The resulting equation solved by the code for a scalar f is:
!
!  rho*volume*df/dt + .... = crvimp*f + crvexp
!
!
! Note that crvexp and crvimp are defined after the Finite Volume integration
! over the cells, so they include the "volume" term. More precisely:
!   - crvexp is expressed in kg.[scal]/s, where [scal] is the unit of the scalar
!   - crvimp is expressed in kg/s
!
!
! The crvexp and crvimp arrays are already initialized to 0 before entering the
! the routine. It is not needed to do it in the routine (waste of CPU time).
!
! For stability reasons, Code_Saturne will not add -crvimp directly to the
! diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is
! treated implicitely only if it strengthens the diagonal of the matrix.
! However, when using the second-order in time scheme, this limitation cannot
! be done anymore and -crvimp is added directly. The user should therefore test
! the negativity of crvimp by himself.
!
! When using the second-order in time scheme, one should supply:
!   - crvexp at time n
!   - crvimp at time n+1/2
!
!
! The selection of cells where to apply the source terms is based on a getcel
! command. For more info on the syntax of the getcel command, refer to the
! user manual or to the comments on the similar command getfbr in the routine
! cs_user_boundary_conditions.

! WARNING: If scalar is the temperature, the resulting equation
!          solved by the code is:
!
!  rho*Cp*volume*dT/dt + .... = crvimp*T + crvexp
!
!
! Note that crvexp and crvimp are defined after the Finite Volume integration
! over the cells, so they include the "volume" term. More precisely:
!   - crvexp is expressed in W
!   - crvimp is expressed in W/K
!

!
! STEEP SOURCE TERMS
!===================
! In case of a complex, non-linear source term, say F(f), for scalar f, the
! easiest method is to implement the source term explicitely.
!
!   df/dt = .... + F(f(n))
!   where f(n) is the value of f at time tn, the beginning of the time step.
!
! This yields :
!   crvexp = volume*F(f(n))
!   crvimp = 0
!
! However, if the source term is potentially steep, this fully explicit
! method will probably generate instabilities. It is therefore wiser to
! partially implicit the term by writing:
!
!   df/dt = .... + dF/df*f(n+1) - dF/df*f(n) + F(f(n))
!
! This yields:
!   crvexp = volume*( F(f(n)) - dF/df*f(n) )
!   crvimp = volume*dF/df

!-------------------------------------------------------------------------------
! 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
!>                                     (see 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 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(ncepdp,6), smacel(ncesmp,nvar)
double precision crvexp(ncelet), crvimp(ncelet)
double precision crvexp2(ncelet), crvimp2(ncelet)

! Local variables

integer          iel, nlfac
integer          ifac, iflmas

double precision ubulk, surf, aT, surfq, tbulk, forcing
double precision, dimension(:), pointer :: imasfl
double precision, dimension(:,:), pointer :: vel
double precision, dimension(:), pointer :: tempc 
!double precision, allocatable, dimension (:) :: masfl, masflb
integer, allocatable, dimension(:) :: lstfac
!double precision :: epaiss = 5.d0, larg = 1.5d0, Long = 40.d0, Haut = 20.d0
!double precision :: qpoint = 1.d0

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

if (iscal.ne.iscalt) return

ubulk = 0.d0
surf  = 0.d0
tbulk = 0.d0

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

call field_get_key_int(ivarfl(iv), kimasf, iflmas)
call field_get_val_s(iflmas, imasfl)
call field_get_val_s(ivarfl(isca(iscalt)), tempc)

!calculate ubulk and surf at the inlet 
allocate (lstfac(nfac))
call getfac('2',nlfac,lstfac)

do iel = 1, nlfac
  ifac = lstfac(iel)
  ubulk = ubulk + imasfl(ifac)*surfac(2,ifac)/surfan(ifac)
  surf  = surf + surfan(ifac)
  tbulk = tbulk + vel(2,ifacel(2,ifac))*tempc(ifacel(2,ifac))*surfan(ifac)
enddo

if(irangp.ge.0) call parsom(ubulk)
if(irangp.ge.0) call parsom(surf)
if(irangp.ge.0) call parsom(tbulk)

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

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

tbulk = tbulk/(surf*ubulk)

write(nfecra,*) 'ubulk ustssc = ', ubulk,&
                'tbulk ustssc = ', tbulk,&
                't0 = ', T0,&
                'tref-tbulk =', T0-tbulk

if((ntcabs-ntpabs).eq.1) then
  forcing = 0.d0
else
  open(file='Tbulk.dat',unit=impusr(2))
  read(impusr(2),*) forcing
  close(impusr(2))
endif

!calculate a forcing term fo ensure that tbulk is equal to tref (which is T0)
forcing = forcing + 0.5d0*(493.15d0-tbulk)

!calculate the surface area of the heated walls
call getfbr('4 or 5', nlfac, lstfac)

surfq=0.d0
do iel=1,nlfac
  ifac=lstfac(iel)
!  surfmod=sqrt(surfbo(1,ifac)**2+surfbo(2,ifac)**2+surfbo(3,ifac))
  surfq=surfq+surfbn(ifac)
enddo  

if(irangp.ge.0) call parsom(surfq)
write(nfecra,*) 'heated wall surface = ', surfq
! calculation of the streamwise temperature gradient (insert heat flux and streamwise domain length)
aT = surfq*65.8d3/(ubulk*surf*cp0*ro0*1.044d-3)
write(nfecra,*)'cp0',cp0
write(nfecra,*)'ro0',ro0


!store the calculated values of forcing 
if(irangp.le.0) then
  if((ntcabs-ntpabs).eq.1) then
     open(file="Tbulk.dat",unit=impusr(2),status='new')
  else
     open(file="Tbulk.dat",unit=impusr(2))
  endif
  write(impusr(2),*) forcing
  close(impusr(2))
endif

!write tbulk, forcing at every time step
open(file='temp_source.dat',unit=impusr(4))
if((ntcabs-ntpabs).eq.1) then
  write(impusr(4),*) 'iter ','tbulk ','forcing'
endif
write(impusr(4),*) ntcabs-ntpabs,tbulk,forcing
if(ntcabs.eq.ntmabs) then
  close(impusr(4))
endif

!calculate and impose the source term for the temperature field
do iel=1,ncel
  crvexp(iel) = volume(iel)*(1.d0+forcing)*aT*vel(2,iel)
  !write(nfecra,*) 'vel(2,iel)=',vel(2,iel)
enddo

deallocate(lstfac)
return
end subroutine ustssc