Upload src

Large_canyon/cs_user_physical_properties.f90

+!-------------------------------------------------------------------------------
+
+!VERS
+
+! This file is part of Code_Saturne, a general-purpose CFD tool.
+!
+! Copyright (C) 1998-2019 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.
+
+!-------------------------------------------------------------------------------
+
+!===============================================================================
+!> \file cs_user_physical_properties.f90
+!>
+!> \brief Definition of physical variable laws.
+!>
+!> See \subpage physical_properties for examples.
+!>
+
+!===============================================================================
+!> \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 (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.
+!
+!-------------------------------------------------------------------------------
+
+!-------------------------------------------------------------------------------
+! 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)
+!_______________________________________________________________________________
+
+subroutine usphyv &
+ ( nvar   , nscal  ,                                              &
+   mbrom  ,                                                       &
+   dt     )
+
+!===============================================================================
+
+!===============================================================================
+! Module files
+!===============================================================================
+
+use atincl
+use paramx
+use pointe
+use numvar
+use optcal
+use cstphy
+use cstnum
+use entsor
+use parall
+use period
+use ppppar
+use ppthch
+use ppincl
+use field
+use mesh
+use lagran
+use cs_c_bindings
+
+!===============================================================================
+
+implicit none
+
+! Arguments
+
+integer          nvar   , nscal
+
+integer          mbrom
+
+double precision dt(ncelet)
+
+! Local variables
+
+integer          ivart, iel
+
+double precision xvart, rhum, rscp, pp, zent
+double precision lrhum, lrscp
+double precision qsl, deltaq
+double precision qliq, qwt, tliq, dum
+
+double precision, dimension(:), pointer :: brom, crom
+double precision, dimension(:), pointer :: cvar_vart, cvar_totwt
+double precision, dimension(:), pointer :: cpro_tempc, cpro_liqwt
+double precision, dimension(:), pointer :: cpro_beta
+
+logical activate
+
+!===============================================================================
+! 0. INITIALISATIONS A CONSERVER
+!===============================================================================
+
+activate = .false.
+
+! Initialize variables to avoid compiler warnings
+
+ivart = -1
+
+if (idilat.eq.0) then
+  call field_get_val_s_by_name("thermal_expansion", cpro_beta)
+endif
+
+! This routine computes the density and the thermodynamic temperature.
+! The computations require the pressure profile which is here taken from
+! the meteo file. If no meteo file is used, the user should
+! give the laws for RHO and T in cs_user_physical_properties.f90
+
+if (imeteo.ne.0) return
+
+!===============================================================================
+
+!   Positions des variables, coefficients
+!   -------------------------------------
+
+! --- Numero de variable thermique
+!       (et de ses conditions limites)
+!       (Pour utiliser le scalaire utilisateur 2 a la place, ecrire
+!          IVART = ISCA(2)
+
+if (iscalt.gt.0) then
+  ivart = isca(iscalt)
+else
+  write(nfecra,9010) iscalt
+  call csexit (1)
+endif
+
+! --- Masse volumique
+
+call field_get_val_s(icrom, crom)
+call field_get_val_s(ibrom, brom)
+
+call field_get_val_s(itempc,cpro_tempc)
+call field_get_val_s(ivarfl(ivart), cvar_vart)
+
+! From potential temperature, compute:
+! - Temperature in Celsius
+! - Density
+! ----------------------
+
+! Computes the perfect gaz constants according to the physics
+
+rhum = rair
+rscp = rair/cp0
+
+lrhum = rair
+lrscp = rair/cp0
+
+do iel = 1, ncel
+
+  xvart = cvar_vart(iel) !  The thermal scalar is potential temperature
+
+  zent = xyzcen(3,iel)
+
+  pp = p0
+
+  ! Temperature in Celsius in cell centers:
+  ! ---------------------------------------
+  ! law: T = theta * (p/ps) ** (Rair/Cp0)
+
+  cpro_tempc(iel) = xvart*(pp/ps)**lrscp
+  cpro_tempc(iel) = cpro_tempc(iel) - tkelvi
+
+  !   Density in cell centers:
+  !   ------------------------
+  !   law:    RHO       =   P / ( Rair * T(K) )
+
+  if (idilat.eq.0) then
+    crom(iel) = ro0
+    ! "delta rho = - beta rho0 delta theta" gives
+    ! "beta = 1 / theta"
+    cpro_beta(iel) = 1.d0 / xvart
+  else
+    crom(iel) = pp/(lrhum*xvart)*(ps/pp)**lrscp
+  endif
+
+enddo
+
+!===============================================================================
+! FORMATS
+!----
+
+9010 format(                                                           &
+'@                                                            ',/,&
+'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
+'@                                                            ',/,&
+'@ @@ ATTENTION : ARRET LORS DU CALCUL DES GRANDEURS PHYSIQUES',/,&
+'@    =========                                               ',/,&
+'@    APPEL A csexit DANS LE SOUS PROGRAMME atphyv            ',/,&
+'@                                                            ',/,&
+'@    La variable dont dependent les proprietes physiques ne  ',/,&
+'@      semble pas etre une variable de calcul.               ',/,&
+'@    En effet, on cherche a utiliser la temperature alors que',/,&
+'@      ISCALT = ',I10                                         ,/,&
+'@    Le calcul ne sera pas execute.                          ',/,&
+'@                                                            ',/,&
+'@    Verifier le codage de cs_user_physical_properties,',      /,&
+'@     (et le test lors de la definition de IVART).'           ,/,&
+'@    Verifier la definition des variables de calcul dans     ',/,&
+'@      usipsu. Si un scalaire doit jouer le role de la       ',/,&
+'@      temperature, verifier que ISCALT a ete renseigne.     ',/,&
+'@                                                            ',/,&
+'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
+'@                                                            ',/)
+
+!----
+! FIN
+!----
+
+return
+
+end subroutine usphyv