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

!                      Code_Saturne version 5.0.3
!                      --------------------------
! 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_parameters.f90
!>
!> \brief User subroutines for input of calculation parameters (Fortran modules).
!>        These subroutines are called in all cases.
!>
!>  See \subpage f_parameters for examples.
!>
!>   If the Code_Saturne GUI is used, this file is not required (but may be
!>   used to override parameters entered through the GUI, and to set
!>   parameters not accessible through the GUI).
!>
!>   Several routines are present in the file, each destined to defined
!>   specific parameters.
!>
!>   To modify the default value of parameters which do not appear in the
!>   examples provided, code should be placed as follows:
!>   - usipsu   for numerical and physical options
!>   - usipes   for input-output related options
!>
!>   As a convention, "specific physics" defers to the following modules only:
!>   pulverized coal, gas combustion, electric arcs.
!>
!>   In addition, specific routines are provided for the definition of some
!>   "specific physics" options.
!>   These routines are described at the end of this file and will be activated
!>   when the corresponding option is selected in the usppmo routine.
!-------------------------------------------------------------------------------

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

!> \brief User subroutine for selection of specific physics module

!> Define the use of a specific physics amongst the following:
!>   - combustion with gas / coal / heavy fuel oil
!>   - compressible flows
!>   - electric arcs
!>   - atmospheric modelling
!>   - radiative transfer
!>   - cooling towers modelling
!>
!>    Only one specific physics module can be activated at once.

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     ixmlpu        indicates if the XML file from the GUI is used
!>                              (1 : yes, 0 : no)
!______________________________________________________________________________!

subroutine usppmo &
 ( ixmlpu )

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

use paramx
use entsor
use cstphy
use ppppar
use ppthch
use ppincl
use ppcpfu
use coincl
use radiat
use cs_c_bindings

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

implicit none

! Arguments

integer ixmlpu

!< [usppmo]

!===============================================================================
! 1.  Choice for a specific physics
!===============================================================================

! --- cod3p: Diffusion flame with complete fast chemistry (3 points)
! ==========

!        if = -1   module not activated
!        if =  0   adiabatic model
!        if =  1   extended model with enthalpy source term

if (ixmlpu.eq.0) then

  ippmod(icod3p) = -1

endif

! --- coebu: Eddy-Break Up pre-mixed flame
! ==========

!        if = -1   module not activated
!        if =  0   reference Spalding model
!                   (adiabatic, homogeneous mixture fraction)
!        if =  1   extended model with enthalpy source term
!                   (homogeneous mixture fraction : perfect premix)
!        if =  2   extended model with mixture fraction transport
!                   (adiabatic, no variance of mixture fraction)
!        if =  3   extended model with enthalpy and mixture fraction transport
!                   (dilution, thermal losses, etc.)

if (ixmlpu.eq.0) then

  ippmod(icoebu) = -1

endif

! --- colwc: Libby-Williams pre-mixed flame
! ==========

!        if = -1   module not activated
!        if =  0   reference two-peak model with adiabatic condition
!        if =  1   extended two-peak model with enthapy source terms
!        if =  2   extended three-peak model, adiabatic
!        if =  3   extended three-peak model with enthalpy source terms
!        if =  4   extended four-peak model, adiabatic
!        if =  5   extended four-peak model with enthalpy source terms

if (ixmlpu.eq.0) then

  ippmod(icolwc) = -1

endif


! --- Soot model
! =================

!        if = -1   module not activated
!        if =  0   constant fraction of fuel Xsoot
!        if =  1   2 equations model of Moss et al.
if (.false.) then
  isoot = 0
  
  xsoot  = 0.1d0 ! ( if isoot = 0 )
  rosoot = 2000.d0 ! kg/m3
end if
! --- cfuel: Heavy fuel oil combustion
! ==========

!        Progressive evaporation (temperature gap)
!        Char residue
!        Sulphur tracking

!        if = -1   module not activated
!        if = 0    module activated

if (ixmlpu.eq.0) then

  ippmod(icfuel) = -1

endif

! --- coal :
! ==========
!
!     Pulverized coal combustion
!        Description of granulometry
!        Assumption of diffusion flame around particles
!         (extension of 3-point fast chemistry "D3P")
!        Between a mixture of gaseous fuels (volatiles matters, CO from char
!                                            oxydation)
!            and a mixture of oxidisers (air and water vapor)
!        Enthalpy for both mix and solid phase are solved
!
!        if = -1   module not activated
!        if = 0    module activated
!        if = 1    with drying

if (ixmlpu.eq.0) then

  ippmod(iccoal) = -1

endif

! Activate the drift: 0 (no activation),
!                     1 (transported particle velocity)
!                     2 (limit drop particle velocity)

if (.false.) i_comb_drift = 1


! --- cpl3c: Pulverized coal with Lagrangian reciprocal approach
! ==========

!        Not recently tested... at least outdated, may be obsolete

!        if = -1   module not activated
!        if = 0    module activated
!        if = 1    with drying (NOT functional)

if (ixmlpu.eq.0) then

  ippmod(icpl3c) = -1

endif

! --- compf: Compressible flows
! ==========

!        if = -1   module not activated
!        if = 0    module activated

if (ixmlpu.eq.0) then

  ippmod(icompf) = -1

endif

! --- eos: equation of state for compressible flows
! ========

!        if = 1    ideal gas with constant gamma
!        if = 2    stiffened gas
!        if = 3    ideal gas mix

if (ixmlpu.eq.0.and.ippmod(icompf).ge.0) then

  ieos = 1

endif

! --- eljou: Joule effect
! ==========

!        if = -1   module not activated
!        if = 1    Potentiel reel
!        if = 2    Potentiel complexe
!        if = 3    Potentiel reel     + CDL Transfo
!        if = 4    Potentiel complexe + CDL Transfo

if (ixmlpu.eq.0) then

  ippmod(ieljou) = -1

endif

! --- elarc: Electric arcs
! ==========

!        if = -1   module not activated
!        if = 1    electric potential
!        if = 2    electric potential and vector potential (hence 3D modelling)

if (ixmlpu.eq.0) then

  ippmod(ielarc) = 2

endif

! --- atmos: Atmospheric flows
! ==========

!        if = -1   module not activated
!        if = 0    standard modelling
!        if = 1    dry atmosphere
!        if = 2    humid atmosphere (experimental)

if (ixmlpu.eq.0) then

  ippmod(iatmos) = -1

endif

! --- aeros: Cooling towers
! ==========

!        if = -1   module not activated
!        if = 0    no model (NOT functional)
!        if = 1    Poppe's model
!        if = 2    Merkel's model

if (ixmlpu.eq.0) then

  ippmod(iaeros) = -1

endif

! --- igmix: Gas mixtures modelling
! ==========
!        if =-1 module not activated
!        if = 0  Air/Helium   gas mixtures
!        if = 1  Air/Hydrogen gas mixtures
!        if = 2  Air/Steam    gas mixtures
!        if = 3  Air/Helium/Steam gas mixtures
!        if = 4  Air/Hydrogen/Steam gas mixtures


if (.false.) ippmod(igmix) = 0


! Radiative transfer module (iirayo)
!--------------------------
!        if = 0: not activated (Default)
!        if = 1: DOM
!        if = 2: approximation P1 method

if (.false.) iirayo = 1

! --- richards model
! ==========

!        if = -1   module not activated
!        if =  1   module activated

if (.false.) ippmod(idarcy) = -1

!===============================================================================
! 2.  Specific options related to herebefore modules
!===============================================================================

! These options are defined here at the moment, this might change in the future

! --- Enthalpy-Temperature conversion law (for gas combustion modelling)

!       if = 0   user-specified
!       if = 1   tabulated by JANAF (default)

if (ixmlpu.eq.0) then

  indjon = 1

endif

! --- Kinetic model for CO <=> CO2

!         Compatible with coal and heavy fuel oil combustion

!         if = 0  unused (maximal conversion in turbulent model)
!         if = 1  transport of CO2 mass fraction
!         if = 2  transport of CO mass fraction

if (ixmlpu.eq.0) then

  ieqco2 = 0

endif

!===============================================================================
! 2.  Data file related to modules above
!===============================================================================

if (ixmlpu.eq.0) then

  ! Combustion

  if (     ippmod(icod3p).ge.0                                          &
      .or. ippmod(icoebu).ge.0 .or. ippmod(icolwc).ge.0) then

    if (indjon.eq.1) then
      ficfpp = 'dp_C3P'
    else
      ficfpp = 'dp_C3PSJ'
    endif

  endif

  ! Fuel combustion

  if (ippmod(icfuel).ge.0) then
    ficfpp = 'dp_FUE'
  endif

  ! Atmospheric flows

  if (ippmod(iatmos).ge.0) then
    ficmet = 'meteo'
  endif

 if ( ippmod(igmix).ge.0 ) then
   ! Specific condensation modelling

   ! wall condensation
   !      if = -1 module not activated
   !      if =  0 condensation source terms activated
   icondb = -1

   ! internal condensation
   !      if = -1 module not activated
   !      if =  0 condensation source terms with metal
   !                               structures activate
   icondv = -1
 endif

endif

!< [usppmo]

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

return
end subroutine usppmo


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

!> \brief User subroutine for input of model selection parameters.

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]      ixmlpu       indicates if the XML file from the GUI is used
!>                              used (1: yes, 0: no
!> \param[in, out] iturb        turbulence model
!> \param[in, out] itherm       thermal model
!> \param[in, out] iale         ale module
!> \param[in, out] ivofmt       vof method
!> \param[in, out] icavit       cavitation model
!______________________________________________________________________________!

subroutine usipph &
 ( ixmlpu, iturb , itherm, iale , ivofmt, icavit )

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

use entsor, only: nfecra ! No other module should appear here
use optcal, only: irijco ! No other module should appear here

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

implicit none

! Arguments

integer ixmlpu
integer iturb, itherm, iale, ivofmt, icavit

! Local variables

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

!>    In this subroutine, only the parameters which already appear may
!>    be set, to the exclusion of any other.
!>
!>    If we are not using the Code_Saturne GUI:
!>    All the parameters which appear in this subroutine must be set.
!>
!>    If we are using the Code_Saturne GUI:
!>    parameters protected by a test of the form:
!>
!>      if (ixmlpu.eq.0) then
!>         ...
!>      endif
!>
!>    should already have been defined using the GUI, so only
!>    experts should consider removing the test and adapting them here.

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

!< [usipph]

! --- Turbulence
!       0: Laminar
!      10: Mixing length
!      20: k-epsilon
!      21: k-epsilon (linear production)
!      30: Rij-epsilon, (standard LRR)
!      31: Rij-epsilon (SSG)
!      32: Rij-epsilon (EBRSM)
!      40: LES (Smagorinsky)
!      41: LES (Dynamic)
!      42: LES (WALE)
!      50: v2f (phi-model)
!      51: v2f (BL-v2/k)
!      60: k-omega SST
!      70: Spalart Allmaras
!  For 10, contact the development team before use

if (ixmlpu.eq.0) then

  iturb = 0 !WE DON'T USE TURBULENCE FOR ELECTRIC ARCS SIMULATION

endif


! Coupled solver for Rij components (when iturb=30, 31 or 32)

if (.false.) irijco = 1

! --- Thermal model
!      0: none
!      1: temperature
!      2: enthalpy
!      3: total energy (only for compressible module)
!
!  For temperature, the temperature scale may be set later using itpscl
!  (1 for Kelvin, 2 for Celsius).
!
!  Warning: When using specific physics, this value is
!           set automatically by the physics model.

if (.false.) then
  itherm = 2 !ENTHALPY MUST BE THERMAL SCALAR FOR ELECTRIC ARCS SIMULATION 
endif

! --- Cavitation module
!    - -1: module not activated
!    -  0: no vaporization/condensation model
!    -  1: Merkle's model
!
!  Specific cavitation module input parameters should be set usipsu
!  (see example in cs_user_parameters-cavitation.f90)
!


if (.false.) icavit = -1


! --- Activation of ALE (Arbitrary Lagrangian Eulerian) method


if (.false.) iale = 1

!< [usipph]

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


return
end subroutine usipph


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

!> \brief User subroutine for the input of additional user parameters.
!
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     nmodpp         number of active specific physics models
!______________________________________________________________________________!

subroutine usipsu &
 ( nmodpp )

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

use paramx
use cstnum
use dimens
use numvar
use optcal
use cstphy
use entsor
use parall
use period
use ihmpre
use albase
use ppppar
use ppthch
use ppincl
use coincl
use cpincl
use field
use cavitation
use post
use rotation
use cs_c_bindings
!use cs_var_cal_opt_t 

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

implicit none

! Arguments

integer nmodpp

! Local variables

logical       inoprv
integer       ii, jj, ivar, kscmin, kscmax, keydri, kbfid, kccmin, kccmax
integer       klimiter
integer       f_id, idim1, itycat, ityloc, iscdri, iscal, ifcvsl, b_f_id
double precision minent
type(var_cal_opt) :: vcopt

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

!>  This subroutine allows setting parameters
!>  which do not already appear in the other subroutines of this file.
!>
!>  It is possible to add or remove parameters.
!>  The number of physical properties and variables is known here.

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

!< [usipsu]

! Calculation options (optcal)
! ============================

! In case of restart, read auxiliary restart file ileaux (= 1) or not (0).

! By default, this file is read, but it may be useful to deactivate
! its use when restarting after a preprocessing stage possibly leading
! to a different number of faces (such as simply joining meshes on
! a different architecture or optimization level or with different options).

! Writing of auxiliary restart files may also be deactivated using: iecaux = 0

if (.false.) ileaux = 0


! --- Time stepping  (0 : uniform and constant
!                     1 : variable in time, uniform in space
!                     2 : variable in time and space
!                    -1 : steady algorithm)

if (.false.) idtvar = 0


! --- Duration
!       ntmabs = absolute number of the last time step required
!         if we have already run 10 time steps and want to
!         run 10 more, ntmabs must be set to 10 + 10 = 20

if (.false.) ntmabs = 10


! --- Reference time step
!     The example given below is probably not adapted to your case.


if (.false.) dtref  = 1d-7

! --- Maximum time step: dtmax
!     Set a value base on characteristic values of your case.
!      otherwise, the code will use a multiple of dtref by default.
!     Example with
!        Ld: "dynamic" length (for example, the domain length)
!        Ud: characteristic flow velocity
!        Lt: thermal length (for example, the domain height gravity-wise)
!        Delta_rho/rho: relative density difference
!        g: gravity acceleration

!     dtmax = min(Ld/Ud, sqrt(Lt/(g.Delta_rho/rho)))

! --- Handling of hydrostatic pressure
!     iphydr = 0 : ignore hydrostatic pressure (by default)
!              1 : with hydrotatic pressure computation to handle the balance
!                  between the pressure gradient and source terms (gravity and
!                  head losses)
!              2 : with hydrostatic pressure computation to handle the imbalance
!                  between the pressure gradient and gravity source term

if (.false.) iphydr = 1


! --- Algorithm to take into account the density variation in time
!
!     idilat = 0 : Boussinesq algorithm with constant density (not available)
!              1 : dilatable steady algorithm (default)
!              2 : dilatable unsteady algorithm
!              3 : low-Mach algorithm

if (.false.) idilat = 1

! --- Algorithm to take into account the thermodynamical pressure variation in time
!     (not used by default except if idilat = 3)

!     by default:
!     ----------
!      - the thermodynamic pressure (pther) is initialized with p0 = p_atmos
!      - the maximum thermodynamic pressure (pthermax) is initialized with -1
!        (no maximum by default, this term is used to model a venting effect when
!         a positive value is given by the user)
!      - a global leak can be set through a leakage surface sleak with a head
!      loss kleak of 2.9 (Idelcick)
if (.false.) then
  ipthrm = 0

  pthermax= -1.d0

  sleak = 0.d0
  kleak = 2.9d0
end if

! --- Temperature or enthalpy

!   When used without specific physics, if we have chosen to solve in temperature
!     (that is if itherm = 1), the fluid temperature is considered to be in
!     degrees Kelvin by default (be careful for boundary conditions an expression
!     of physical properties depending on temperature)t.

!     If we wish for the fluid solver to work with a temperature in degrees Celsius,
!     we must set itpscl = 2.

!     This is recommended for Syrthes Coupling, but not recommended for the
!     radiative model, as it is a source of user errors in this case:
!     Indeed, the boundary conditions for the fluid temperature will then be
!     in degrees Celsius, while the boundary conditions for radiation in
!     cs_user_radiative_transfer_bcs must still be in Kelvin.

if (nmodpp.eq.0) then
  itpscl = 2
endif


!   If a USER scalar behaves like a temperature (relative to Cp):
!     we set iscacp(isca) = 1.
!
!   Otherwise, we do not modify iscacp(isca)

if (.false.) then
  if (nscaus.gt.0) then
    do ii = 1, nscaus
      iscacp(isca(ii)) = 1
    enddo
  endif
end if

! --- Solver taking a scalar porosity into account:
!       0 No porosity taken into account (Standard)
!       1 Porosity taken into account
!

if (.false.)iporos = 1


! --- Calculation (restart) with frozen velocity field (1 yes, 0 no)


if (.false.) iccvfg = 1


! Physical constants (cstphy)
! ===========================

! --- gravity (g in m/s2, with the sign in the calculation coordinate axes).

gx = 0.d0
gy = 0.d0
gz = 9.81d0


! --- rotation of the reference frame (omega in rad/s)

!       If the rotation is not nul, then
!          icorio = 0: rotation is taken into account by rotating the mesh
!                      (simulation in the absolute frame)
!                 = 1: rotation is taken into account by Coriolis source terms
!                      (simulation in the relative frame)

if (.false.) then
  icorio = 0

  call rotation_define(0.d0, 0.d0, 0.d0,  &    ! rotation vector
                       0.d0, 0.d0, 0.d0)       ! invariant point
end if


! --- Reference fluid properties

!       ro0        : density in kg/m3
!       viscl0     : dynamic viscosity in kg/(m s)
!       cp0        : specific heat in J/(Kelvin kg)
!       t0         : reference temperature in Kelvin
!       p0         : total reference pressure in Pascal
!                    the calculation is based on a
!                    reduced pressure P*=Ptot-ro0*g.(x-xref)
!                    (except in compressible case)
!       xyzp0(3)   : coordinates of the reference point for
!                    the total pressure (where it is equal to p0)

!     In general, it is not necessary to furnish a reference point xyz0.
!       If there are outlets, the code will take the center of the
!       reference outlet face.
!       On the other hand, if we plan to explicitly fix Dirichlet conditions
!       for pressure, it is better to indicate to which reference the
!       values relate (for a better resolution of reduced pressure).


!     Other properties are given by default in all cases.

!     Nonetheless, we may note that:

!       In the standard case (no gas combustion, coal, electric arcs,
!                             compressibility):
!       ---------------------
!         ro0, viscl0 and cp0
!             are useful and represent either the fluid properties if they
!             are constant, either simple mean values for the initialization
!             if properties are variable and defined in usphyv.
!         t0  is not useful
!         p0  is useful but is not used in an equation of state. p0
!             is a reference value for the incompressible solver
!             which will serve to set the (possible) domain outlet pressure.
!             We may also take it as 0 or as a physical value in Pascals.

!       With the electric module:
!       ------------------------
!         ro0, viscl0 and cp0
!             are useful but simply represent mean initial values;
!             the density, molecular dynamic viscosity, and specific
!             heat are necessarily defined as fields (whether they are
!             physically variable or not): see cs_user_physical_properties
!             for the Joule effect
!             module and the electric arcs dp_ELE data file.
!         t0  is useful an must be in Kelvin (> 0) but represents a simple
!             initialization value.
!         p0  is useful bu is not used in the equation of state. p0
!             is a reference value for the incompressible solver which
!             will be used to calibrate the (possible) outlet pressure
!             of the domain. We may take it as zero or as a physical
!             value in Pascals.

!       With gas combustion:
!       --------------------
!         ro0 is not useful (it is automatically recalculated by the
!             law of ideal gases from t0 and p0).
!         viscl0 is indispensable: it is the molecular dynamic viscosity,
!             assumed constant for the fluid.
!         cp0 is indispensable: it is the heat capacity, assumed constant,
!             (modelization of source terms involving a local Nusselt in
!             the Lagrangian module, reference value allowing the
!             calculation of a radiative
!             (temperature, exchange coefficient) couple).
!         t0  is indispensible and must be in Kelvin (> 0).
!         p0  is indispensable and must be in Pascal (> 0).

!       With pulverized coal:
!       ---------------------
!         ro0 is not useful (it is automatically recalculated by the
!             law of ideal gases from t0 and p0).
!         viscl0 is indispensable: it is the molecular dynamic viscosity,
!             assumed constant for the fluid (its effect is expected to
!             be small compared to turbulent effects).
!         cp0 is indispensable: it is the heat capacity, assumed constant,
!             (modelization of source terms involving a local Nusselt in
!             the coal or Lagrangian module, reference value allowing the
!             calculation of a radiative
!             (temperature, exchange coefficient) couple).
!         t0  is indispensable and must be in Kelvin (> 0).
!         p0  is indispensable and must be in Pascal (> 0).

!       With compressibility:
!       ---------------------
!         ro0 is not useful, stricto sensu; nonetheless, as experience
!             shows that users often use this variable, it is required
!             to assign to it a strictly positive value (for example,
!             an initial value).
!         viscl0 is useful and represents the molecular dynamic viscosity,
!             when it is constant, or a value which will be used during
!             initializations (or in inlet turbulence conditions,
!             depending on the user choice.
!         cp0 is indispensable: it is the heat capacity, assumed constant
!             in the thermodynamics available by default
!         t0  is indispensable and must be in Kelvin (> 0).
!         p0  is indispensable and must be in Pascal (> 0).
!             With the thermodynamic law available by default,
!             t0 and p0 are used for the initialization of the density.
!         xyzp0 is not useful because the pressure variable directly
!             represents the total pressure.

ro0    = 1.24d0
viscl0 = 2.16d-5
cp0    = 751.0d0
t0 = 20.D0 + 273.15d0
p0 = 1.01325d5


! We only specify XYZ0 if we explicitely fix Dirichlet conditions
! for the pressure.

xyzp0(1) = 0.d0   ! I had it in Maher's case. May be useless.
xyzp0(2) = 0.d0   ! We only specify XYZ0 if we explicitely fix Dirichlet conditions for the pressure.
xyzp0(3) = 0.11d-2

! --- Minimum and maximum admissible values for each USER scalar:

!      Results are clipped at the end of each time step.

!      If min > max, we do not clip.

!      For a scalar jj representing the variance of another, we may
!        abstain from defining these values
!        (a default clipping is set in place).
!        This is the purpose of the test on iscavr(jj) in the example below.

!      For non-user scalars relative to specific physics (coal, combustion,
!        electric arcs: see usppmo) implicitly defined according to the
!        model, the information is automatically set elsewhere: we
!        do not set min or max values here.

call field_get_key_id("min_scalar_clipping", kscmin)
call field_get_key_id("max_scalar_clipping", kscmax)

! Thermal scalar:
if (iscalt.gt.0) then
  ! We define the min and max bounds of the thermal scalar. For electric arc thermal scalar is enthalpy.
  minent=1.4D4
  write(*,'(A,F10.4,A,E10.2)')'Enthalpy clipping: Min=',minent,' Max=',grand
  call field_set_key_double(ivarfl(isca(iscalt)), kscmin, -grand)
  call field_set_key_double(ivarfl(isca(iscalt)), kscmax, +grand)
endif
write(*,'(A,I4)')'Type of gradient reconstruction imrgra=',imrgra
! Loop on user scalars:
if (.false.) then
  do jj = 1, nscaus
    ! For scalars which are not variances
    if (iscavr(jj).le.0) then
      ! We define the min and max bounds
      call field_set_key_double(ivarfl(isca(jj)), kscmin, -grand)
      call field_set_key_double(ivarfl(isca(jj)), kscmax, +grand)
    endif
  enddo
  
  ! --- Convective scheme for user (and non-user) scalars
  
  ! ischcv is the type of convective scheme:
  !   - 0: second order linear upwind
  !   - 1: centered
  !   - 2: pure upwind gradient in SOLU
  
  ! isstpc is the slope test, Min/Max limiter or Roe and Sweby limiters
  !   - 0: swich on the slope test
  !   - 1: swich off the slope test (default)
  !   - 2: continuous limiter ensuring boundedness (beta limiter)
  !   - 3: NVD/TVD Scheme
  !        Then "limiter_choice" keyword must be set:
  !        * 0: Gamma
  !        * 1: SMART
  !        * 2: CUBISTA
  !        * 3: SUPERBEE
  !        * 4: MUSCL
  !        * 5: MINMOD
  !        * 6: CLAM
  !        * 7: STOIC
  !        * 8: OSHER
  !        * 9: WASEB
  !        * --- VOF scheme ---
  !        * 10: M-HRIC
  !        * 11: M-CICSAM
  
  ! Get the Key for the Sup and Inf for the convective scheme
  call field_get_key_id("min_scalar", kccmin)
  call field_get_key_id("max_scalar", kccmax)
  
  ! Thermal model:
  if (iscalt.gt.0) then
    ivar = isca(iscalt)
  
    call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
    vcopt%ischcv = 0
    vcopt%isstpc = 3
    call field_set_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
  
    ! Get the Key for the limiter choice of the studied scalar
    call field_get_key_id("limiter_choice", klimiter)
    call field_set_key_int(ivarfl(isca(iscalt)), klimiter, 3)! Set SUPERBEE
  
    ! Set the Value for the Sup and Inf of the studied scalar
    call field_set_key_double(ivarfl(ivar), kccmin, 0.d0)
    call field_set_key_double(ivarfl(ivar), kccmax, 1.d0)
  
  endif
  
  ! We loop on user scalars:
  do jj = 1, nscaus
    ivar = isca(jj)
  
    call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
    vcopt%ischcv = 0
    vcopt%isstpc = 2
    call field_set_key_struct_var_cal_opt(ivarfl(ivar), vcopt)
  
    ! Set the Value for the Sup and Inf of the studied scalar
    call field_set_key_double(ivarfl(ivar), kccmin, 0.d0)
    call field_set_key_double(ivarfl(ivar), kccmax, 1.d0)
  enddo
  
  
  ! --- Variable diffusivity field id (ifcvsl>=0) or constant
  !     diffusivity (ifcvsl=-1) for the thermal scalar and USER scalars.
  
  !     With ifcvsl = 0, the field will be added automatically, and later calls to
  !       field_get_key_int(ivarfl(isca(iscal)), kivisl, ifcvsl)
  !       will return its id.
  !     With ifcvsl > 0, the id of an existing, predifined field is given. This
  !       may allow sharing a diffusivity between multiple scalars.
  
  !     For user scalars iscal which represent the variance of another user
  !       scalar, the diffusivity of the variance of a scalar is assumed to
  !       have the same behavior as the diffusivity of this scalar,
  !       so values set here will be ignored.
  
  !     For non-user scalars relative to specific physics (coal, combustion,
  !       electric arcs: see usppmo) implicitly defined in the model,
  !       the diffusivity should not be modified here.
  
  !     Caution:   complete usphyv with the law defining the diffusivity
  !     ========   if and only if ifcvsl = 0 has been set here.
  
  
  ! For thermal scalar
  if (ippmod(icompf).ge.0) then
    ifcvsl = -1
    call field_set_key_int(ivarfl(isca(itempk)), kivisl, ifcvsl)
  else if (iscalt.gt.0) then
    ifcvsl = -1
    call field_set_key_int(ivarfl(isca(iscalt)), kivisl, ifcvsl)
  endif
  
  do iscal = 1, nscaus
    if (iscavr(iscal).le.0) then
      ifcvsl = -1
      call field_set_key_int(ivarfl(isca(iscal)), kivisl, ifcvsl)
    endif
  enddo
  
  
  ! --- Variable density field id (ifcvsl>=0) or bulk
  !     density (ifcvsl=-1) for USER scalars.
  
  !     With ifcvsl = 0, the field will be added automatically, and later calls to
  !       field_get_key_int(ivarfl(isca(iscal)), kromsl, ifcvsl)
  !       will return its id.
  !     With ifcvsl > 0, the id of an existing, predifined field is given. This
  !       may allow sharing a density between multiple scalars.
  
  !     For user scalars iscal which represent the variance of another user
  !       scalar, the density of the variance of a scalar is assumed to
  !       have the same behavior as the density of this scalar,
  !       so values set here will be ignored.
  
  !     Caution:   complete usphyv with the law defining the density
  !     ========   if and only if ifcvsl = 0 has been set here.
  
  do iscal = 1, nscaus
    if (iscavr(iscal).le.0) then
      ifcvsl = -1
      call field_set_key_int(ivarfl(isca(iscal)), kromsl, ifcvsl)
    endif
  enddo
  
  
  ! --- Turbulent flux model u'T' for the scalar T
  !     Algebraic Model
  !      0  SGDH
  !      10 GGDH
  !      11 EB-GGDH (Elliptic Blending)
  !      20 AFM
  !      21 EB-AFM (Elliptic Blending)
  !     Model with transport equations
  !      30 DFM
  !      31 EB-DFM (Elliptic Blending)
  
  ! GGDH for thermal scalar:
  if (iscalt.gt.0) iturt(iscalt) = 10
  
  ! GGDH for all the scalars:
  do jj = 1, nscaus
    iturt(jj) = 10
  enddo
  
  
  ! --- Reference velocity for turbulence initialization (m2/s)
  !       (useful only with turbulence)
  
  uref = 1.d0
  
  
  ! --- Reference length scale in meters for initialization
  !       of epsilon (and specific clipping of turbulence, but
  !       this is not the default option)
  !       Assign a value of the order of the largest dimension of the
  !       physical domain in which the flow may develop.
  !       If a negative value is set here, or no value set and the GUI not
  !       used, the cubic root of the domain will be used.
  !       (useful only for turbulence).
  
  almax = 0.5
  
  ! --- Scalar with a drift (key work "drift_scalar_model">0) or without drift
  !       ((key work "drift_scalar_model"=0, default option) for each USER scalar.
  !       - to specify that a scalar have a drift and need the drift computation:
  !       iscdri = ibset(iscdri, DRIFT_SCALAR_ADD_DRIFT_FLUX)
  !
  ! --- Then, for each scalar with a drift, add a flag to specify if
  !     specific terms have to be taken into account:
  !       - thermophoresis terms:
  !       iscdri = ibset(iscdri, DRIFT_SCALAR_THERMOPHORESIS)
  !       - turbophoresis terms:
  !       iscdri = ibset(iscdri, DRIFT_SCALAR_TURBOPHORESIS)
  !       - centrifugal force terms:
  !       iscdri = ibset(iscdri, DRIFT_SCALAR_CENTRIFUGALFORCE)
  
  
  ! Key id for drift scalar
  call field_get_key_id("drift_scalar_model", keydri)
  
  if (nscaus.ge.1) then
  
    iscdri = 1
    iscdri = ibset(iscdri, DRIFT_SCALAR_ADD_DRIFT_FLUX)
  
    if (.false.) then
      iscdri = ibset(iscdri, DRIFT_SCALAR_THERMOPHORESIS)
    endif
  
    if (.false.) then
      iscdri = ibset(iscdri, DRIFT_SCALAR_TURBOPHORESIS)
    endif
  
    if (.false.) then
      iscdri = ibset(iscdri, DRIFT_SCALAR_CENTRIFUGALFORCE)
    endif
  
    iscal = 1
    f_id = ivarfl(isca(iscal))
  
    ! Set the key word "drift_scalar_model" into the field structure
    call field_set_key_int(f_id, keydri, iscdri)
  
  endif
end if

! Postprocessing-related fields
! =============================

! Example: enforce existence of 'yplus', 'tplus' and 'tstar' fields, so that
!          yplus may be saved, or a local Nusselt number may be computed using
!          the post_boundary_nusselt subroutine.
!          When postprocessing of these quantities is activated, those fields
!          are present, but if we need to compute them in the
!          cs_user_extra_operations user subroutine without postprocessing them,
!          forcing the definition of these fields to save the values computed
!          for the boundary layer is necessary.
!Definition of the field with cell numbers !Rodion!

!field_set_key_int(f_id, keyvis, POST_ON_LOCATION + POST_BOUNDARY_NR + POST_MONITOR)
!choose which post-processing you need and put it as sum in the argument:
!parameter (POST_ON_LOCATION=1)
!parameter (POST_BOUNDARY_NR=2)
!parameter (POST_MONITOR=4)

itycat = FIELD_INTENSIVE + FIELD_PROPERTY  !mask of field property and category values
!FIELD_INTENSIVE=1 !An intensive property is a bulk property, meaning that it is a physical property of a system that does not depend on the system size or the amount of material in the system.
!FIELD_EXTENSIVE=2 !By contrast, an extensive property is additive for subsystems.
!FIELD_PROPERTY=16 
ityloc = 1 ! cells
idim1 = 1 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('CellNumber', f_id)
if (f_id.lt.0) then
  call field_create('CellNumber', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1) !CellNumber post treatment
  call field_set_key_int(f_id, keylog, 0)
endif

! Boundary condition for vector potential CALCULATED
ityloc = 3 ! boundary faces
idim1 = 3 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('PotVecBC', f_id)
if (f_id.lt.0) then
  call field_create('PotVecBC', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1) !PotVecBC post treatment
  call field_set_key_int(f_id, keylog, 1) !and in the log
endif
call field_get_id_try('PotVecBC', f_id)
write(*,'(A,I3)')'Field for potential vector boundary condition PotVecBC, Field Id ',f_id

! Boundary condition for vector potential Relaxed
ityloc = 3 ! boundary faces
idim1 = 3 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('PVBCRel', f_id)
if (f_id.lt.0) then
  call field_create('PVBCRel', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1) !PVBCRel post treatment
  call field_set_key_int(f_id, keylog, 1) !and in the log
endif
call field_get_id_try('PVBCRel', f_id)
write(*,'(A,I3)')'Potential vector boundary condition with relax PVBCRel, Field Id ',f_id

! Domain type for interface heat source calculation
ityloc = 1 ! cells
idim1 = 1 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('domain_type', f_id)
if (f_id.lt.0) then
  call field_create('domain_type', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1) ! domain_type post treatment
  call field_set_key_int(f_id, keylog, 0)
endif

ityloc = 1 ! cells
idim1 = 1 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('inb_from_fl', f_id)
if (f_id.lt.0) then
  call field_create('inb_from_fl', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1)  ! inb_from_fl post treatment
  call field_set_key_int(f_id, keylog, 0)
endif

ityloc = 1 ! cells
idim1 = 1 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('intrf_face', f_id)
if (f_id.lt.0) then
  call field_create('intrf_face', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1)  ! intrf_face post treatment
  call field_set_key_int(f_id, keylog, 0)
endif

ityloc = 1 ! cells
idim1 = 1 ! dimension
inoprv = .false. ! no previous time step values needed
call field_get_id_try('inb_heatflux', f_id)
if (f_id.lt.0) then
  call field_create('inb_heatflux', itycat, ityloc, idim1, inoprv, f_id)
  call field_set_key_int(f_id, keyvis, 1)  ! inb_heatflux post treatment
  call field_set_key_int(f_id, keylog, 0)
endif

!THIS IS FOR DEBUGGING
if (.false.) then
  ityloc = 1 ! boundary faces
  idim1 = 1 ! dimension
  inoprv = .false. ! no previous time step values needed
  call field_get_id_try('TempField1', f_id)
  if (f_id.lt.0) then
    call field_create('TempField1', itycat, ityloc, idim1, inoprv, f_id)
    call field_set_key_int(f_id, keyvis, 1) !TempField1 post treatment
    call field_set_key_int(f_id, keylog, 1) !and in the log
  endif
  call field_get_id_try('TempField1', f_id)
  
  ityloc = 1 ! boundary faces
  idim1 = 1 ! dimension
  inoprv = .false. ! no previous time step values needed
  call field_get_id_try('TempField2', f_id)
  if (f_id.lt.0) then
    call field_create('TempField2', itycat, ityloc, idim1, inoprv, f_id)
    call field_set_key_int(f_id, keyvis, 1) !TempField2 post treatment
    call field_set_key_int(f_id, keylog, 1) !and in the log
  endif
  call field_get_id_try('TempField2', f_id)
  
  ityloc = 1 ! boundary faces
  idim1 = 1 ! dimension
  inoprv = .false. ! no previous time step values needed
  call field_get_id_try('TempField3', f_id)
  if (f_id.lt.0) then
    call field_create('TempField3', itycat, ityloc, idim1, inoprv, f_id)
    call field_set_key_int(f_id, keyvis, 1) !TempField3 post treatment
    call field_set_key_int(f_id, keylog, 1) !and in the log
  endif
  call field_get_id_try('TempField3', f_id)
end if

if (.false.) then
  itycat = FIELD_INTENSIVE + FIELD_PROPERTY
  ityloc = 3 ! boundary faces
  idim1 = 1 ! dimension
  inoprv = .false. ! no previous time step values needed
  
  call field_get_id_try('yplus', f_id)
  if (f_id.lt.0) then
    call field_create('yplus', itycat, ityloc, idim1, inoprv, f_id)
    ! yplus postprocessed and in the log
    call field_set_key_int(f_id, keyvis, POST_ON_LOCATION + POST_MONITOR)
    call field_set_key_int(f_id, keylog, 1)
  endif
  
  call field_get_id_try('tplus', f_id)
  if (f_id.lt.0) then
    call field_create('tplus', itycat, ityloc, idim1, inoprv, f_id)
    ! tplus postreated and in the log
    call field_set_key_int(f_id, keyvis, POST_ON_LOCATION + POST_MONITOR)
    call field_set_key_int(f_id, keylog, 1)
  endif
  
  call field_get_id_try('tstar', f_id)
  if (f_id.lt.0) then
    call field_create('tstar', itycat, ityloc, idim1, inoprv, f_id)
    ! tstar postreated and in the log
    call field_set_key_int(f_id, keyvis, POST_ON_LOCATION + POST_MONITOR)
    call field_set_key_int(f_id, keylog, 1)
  endif
end if

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

return
end subroutine usipsu