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

!                      Code_Saturne version 5.0.7
!                      --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2018 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_physical_properties.f90
!>
!> \brief Definition of physical variable laws.
!>
!> See \subpage physical_properties for examples.
!>

!===============================================================================
!> usphyv
!>
!> \brief Definition of physical variable laws.
!>
!> \section Warning
!>
!> It is \b forbidden to modify turbulent viscosity \c visct here
!> (a specific subroutine is dedicated to that: \ref usvist)
!>
!> - icp = 0 must <b> have been specified </b>
!>    in \ref usipsu if we wish to define a variable specific heat
!>    cpro_cp (otherwise: memory overwrite).
!>
!> - the kivisl field integer key (scalar_diffusivity_id)
!>    must <b> have been specified </b>
!>    in \ref usipsu if we wish to define a variable viscosity
!>    \c viscls.
!>
!>
!> \remarks
!>  - This routine is called at the beginning of each time step
!>    Thus, <b> AT THE FIRST TIME STEP </b> (non-restart case), the only
!>    values initialized before this call are those defined
!>      - in the GUI or  \ref usipsu (cs_user_parameters.f90)
!>             - density    (initialized at \c ro0)
!>             - viscosity  (initialized at \c viscl0)
!>      - in the GUI or \ref cs_user_initialization
!>             - calculation variables (initialized at 0 by defaut
!>             or to the value given in the GUI or in \ref cs_user_initialization)
!>
!>  - We may define here variation laws for cell properties, for:
!>     - density:                                    rom    kg/m3
!>     - density at boundary faces:                  romb   kg/m3)
!>     - molecular viscosity:                        cpro_viscl  kg/(m s)
!>     - specific heat:                              cpro_cp     J/(kg degrees)
!>     - diffusivities associated with scalars:      cpro_vscalt kg/(m s)
!>
!> \b Warning: if the scalar is the temperature, cpro_vscalt corresponds
!> to its conductivity (Lambda) in W/(m K)
!>
!>
!> The types of boundary faces at the previous time step are available
!>   (except at the first time step, where arrays \c itypfb and \c itrifb have
!>   not been initialized yet)
!>
!> It is recommended to keep only the minimum necessary in this file
!>   (i.e. remove all unused example code)
!>
!>
!> \section usphyv_cell_id Cells identification
!>
!> Cells may be identified using the \ref getcel subroutine.
!> The syntax of this subroutine is described in the
!> \ref cs_user_boundary_conditions subroutine,
!> but a more thorough description can be found in the user guide.
!
!-------------------------------------------------------------------------------

subroutine usphyv &
 ( nvar   , nscal  ,                                              &
   mbrom  ,                                                       &
   dt     )

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

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

use paramx
use pointe
use numvar
use optcal
use cstphy
use entsor
use parall
use period
use ppppar
use ppthch
use ppincl
use field
use mesh
use cs_cf_bindings

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

implicit none

! Arguments

integer          nvar   , nscal

integer          mbrom

double precision dt(ncelet)

! Local variables

!< [loc_var_dec]
integer          ivart, iel
integer          ith, iscal, ifcvsl
double precision varam, varbm, varcm, vardm
double precision varal, varbl, varcl, vardl
double precision varac, varbc
double precision xvart, mu_ref, T_ref, S

double precision, dimension(:), pointer :: cpro_viscl, cpro_viscv
double precision, dimension(:), pointer :: cpro_vtmpk, cpro_vscal
double precision, dimension(:), pointer :: cpro_cp, cpro_cv, mix_mol_mas
double precision, dimension(:), pointer :: cvar_scalt
!< [loc_var_dec]

ivivar=1

!===============================================================================
! 1. Mandatory initializations
!===============================================================================

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

! Warning: the examples provided below are physically meaningless.
! =======

! These examples must be adapted by the user

! It is adviced to discard all the examples that are not necessary, so
! as to minimize the risk of error.

! List of examples
! ================

! Ex. 1: molecular viscosity varying with temperature
! Ex. 2: molecular volumetric viscosity varying with temperature
! Ex. 3: isobaric specific heat varying with temperature
! Ex. 4: molecular thermal conductivity varying with temperature
! Ex. 5: molecular diffusivity of user-defined scalars varying with temperature

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


!===============================================================================
! Ex. 1: molecular viscosity varying with temperature
! =====
!    The values of the molecular viscosity are provided as a function of
!    the temperature. All variables are evaluated at the cell centers.
!===============================================================================

!!< [example_1]
!!     To refer to the user-defined scalar number 2 instead, for example, use
!!     ivart = isca(2)
!
!ivart = isca(itempk)
!call field_get_val_s(ivarfl(ivart), cvar_scalt)
!
!! --- Molecular dynamic viscosity 'cpro_viscl'
!!     (physical properties at the cell centers)
!
!call field_get_val_s(iviscl, cpro_viscl)
!
!! --- User-defined coefficients for the selected law.
!     The values hereafter are provided as a mere example. They
!!     are physically meaningless.
!
!varam = -3.4016d-9
!varbm =  6.2332d-7
!varcm = -4.5577d-5
!vardm =  1.6935d-3
!! --- Molecular dynamic viscosity mu at the cell centers, kg/(m s)
!!     In this example, mu is provided as a function of the temperature T:
!!       mu(T)              =    T  *( T  *( am  * T +  bm  )+ cm  )+ dm
!!     that is:
!!       cpro_viscl(iel) =   xvart*(xvart*(varam*xvart+varbm)+varcm)+vardm
!
!do iel = 1, ncel
!  xvart = cvar_scalt(iel)
!  cpro_viscl(iel) = xvart*(xvart*(varam*xvart+varbm)+varcm)+vardm
!enddo
!!< [example_1]


!< [SuthTest_1] 
!ivart = isca(itempk)
!call field_get_val_s(ivarfl(isca(itempk)), cvar_scalt)
call field_get_val_s_by_name('temperature',cvar_scalt)

call field_get_val_s(iviscl, cpro_viscl)

mu_ref=1.710E-05 !Moyenne grossière effectuee avec la compo gas entré
                 !http://www.springer.com/978-981-287-211-1 Crane Company (1988)
                 !Flow of fluids through valves, fittings, and pipe. Technical
                 !Paper No. 410 (TP 410). P A-5
T_ref=290.37
S=122.8

do iel = 1, ncel
  xvart = cvar_scalt(iel)
  cpro_viscl(iel) = mu_ref*((xvart/T_ref)**(3/2))*(T_ref+S)/(xvart+S)
enddo



!< [SuthTest_1]


return
end subroutine usphyv
                      








!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     nvar          total number of variables
!> \param[in]     nscal         total number of scalars
!> \param[in]     mbrom         indicator of filling of romb array
!> \param[in]     dt            time step (per cell)
!_______________________________________________________________________________

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

!> usvist
!> \brief Modify turbulent viscosity
!>
!> This subroutine is called at beginning of each time step
!> after the computation of the turbulent viscosity
!> (physical quantities have already been computed in \ref usphyv).
!>
!> Turbulent viscosity \f$ \mu_T \f$ (kg/(m s)) can be modified.
!>
!> A modification of the turbulent viscosity can lead to very
!> significant differences betwwen solutions and even give wrong
!> results.
!>
!> This subroutine is therefore reserved to expert users.
!
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! 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
!> \param[in]     ncesmp        number of cells with mass source term
!> \param[in]     icepdc        head loss cell numbering
!> \param[in]     icetsm        numbering of cells with mass source term
!> \param[in]     itypsm        kind of mass source for each variable
!>                               (cf. \ref cs_user_mass_source_terms)
!> \param[in]     dt            time step (per cell)
!> \param[in]     ckupdc        work array for head loss terms
!> \param[in]     smacel        values of variables related to mass source
!>                              term. If ivar=ipr, smacel=mass flux
!_______________________________________________________________________________

subroutine usvist &
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
   icepdc , icetsm , itypsm ,                                     &
   dt     ,                                                       &
   ckupdc , smacel )

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

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

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



use dimens, only: ndimfb
use pointe
use cstnum
use lagran
use ppppar
use ppthch
use ppincl
use turbomachinery


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

!implicit none

! Arguments

!integer ::  nvar   , nscal
!integer :: ncepdp , ncesmp

!double precision dt(ncelet)

! Local variables
double precision :: x1, y1, xmax

double precision, dimension(:), pointer :: crom
double precision, dimension(:), pointer :: viscl, visct
double precision, dimension(:), pointer :: cvar_k, cvar_ep, cvar_phi
!
!!===============================================================================
!
call field_get_val_s(ivisct, visct)
!
call field_get_val_s(ivarfl(ik), cvar_k)
call field_get_val_s(ivarfl(iep), cvar_ep)
!

call field_get_val_s(icrom, crom)


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

x1=-0.1
y1=0.005
xmax=-1.17

do iel = 1, ncel

!--------- Laminar casei turbulent visco
      xk = cvar_k(iel)
      xe = cvar_ep(iel)
      visct(iel) = crom(iel)*cmu*xk**2/xe
!---------

!---------
!  if (xyzcen(1,iel).gt.x1 .and. xyzcen(2,iel).gt.y1) then
  if (xyzcen(1,iel).lt.x1) then
    visct_save = visct(iel)
    visct(iel) = 99*visct_save/(xmax-x1)*(xyzcen(1,iel)-x1)+visct_save
  endif
enddo

!write(nfecra,*) 'n_iter: ', n_iter
!===============================================================================

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

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

return
end subroutine usvist