Upload src

Large_canyon/cs_user_modules.f90

+!-------------------------------------------------------------------------------
+
+!                      Code_Saturne version 6.1-alpha
+!                      --------------------------
+! 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.
+
+!-------------------------------------------------------------------------------
+
+!===============================================================================
+! Purpose:
+! -------
+
+!> \file cs_user_modules.f90
+!>
+!> \brief User-defined module: it allows to create any user array.
+!>
+!> See \subpage cs_user_modules for examples.
+!>
+!> This file is compiled before all other user Fortran files.
+!> To ensure this, it must not be renamed.
+!>
+!> The user may define an arbitrary number of modules here, even though
+!> only one is defined in the example.
+!
+!> \cond DOXYGEN_SHOULD_SKIP_THIS
+
+!-------------------------------------------------------------------------------
+
+module user_module
+
+   use parall, only : irangp, nrangp, parbcr, parsom, parmax, parmin
+   use cstnum, only : pi
+
+   implicit none
+
+   ! User variables
+   double precision, parameter :: longueur_rue = 1139.54d0
+
+   ! variables pour arbres
+   character(len=500) :: zone_arbre
+   double precision volume_arbre
+   double precision surface_arbre
+   double precision surface_projetee_arbre
+   double precision ratio_volume_ab
+   double precision nombre_arbre, tree_configuration
+   double precision tree_height, tree_radius, tree_position_y
+   double precision alpha
+   double precision, parameter :: microg = 1.d-3 ! from µg to mg
+   double precision, parameter :: lz = 1d0
+
+   contains
+
+   subroutine calcul_volume_zone(zone, volume_zone)
+
+      use mesh ! cell_f_vol
+      ! getcel ?
+
+      implicit none
+
+      ! Arguments
+      character(len=*), intent(in) :: zone
+      double precision, intent(out) :: volume_zone
+      ! Local variables
+      integer :: nlelt, ilelt
+      integer, dimension(ncel) :: lstelt
+
+      volume_zone = 0.d0
+      call getcel(zone, nlelt, lstelt)
+      do ilelt = 1, nlelt
+         volume_zone = volume_zone + cell_f_vol(lstelt(ilelt))
+      enddo
+      if (nrangp.gt.1) then
+         call parsom(volume_zone)
+      endif
+
+      return
+   end subroutine calcul_volume_zone
+
+   subroutine calcul_surface_projetee(zone, surface_projetee)
+
+      use mesh ! cell_f_vol
+      
+      ! getcel ?
+
+      implicit none
+
+      ! Arguments
+      character(len=*), intent(in) :: zone
+      double precision, intent(out) :: surface_projetee
+      ! Local variables
+      integer :: nlelt, ilelt
+      integer, dimension(ncel) :: lstelt
+      double precision :: xmin, xmax, dx
+
+      xmin = 999999999.
+      xmax = - xmin
+      call getcel(zone, nlelt, lstelt)
+      do ilelt = 1, nlelt
+         ! Use coordinates of center of the cell
+         xmin = min(xmin, xyzcen(1,lstelt(ilelt)))
+         xmax = max(xmax, xyzcen(1,lstelt(ilelt)))
+      enddo
+      if (nrangp.gt.1) then
+         call parmax(xmax)
+         call parmin(xmin)
+      endif
+
+      ! Estimate the coordinate of the boundary of the cell
+      ! This is valid only for a 2D mesh with constant cell size
+      dx = sqrt(cell_f_vol(1)/lz)
+      xmin = xmin - dx/2.
+      xmax = xmax + dx/2.
+      surface_projetee = (xmax - xmin) * lz
+
+      return
+   end subroutine calcul_surface_projetee
+
+   !
+   ! Initializer la zone des arbres
+   subroutine initialize_arbre()
+
+      use cs_c_bindings ! notebook_parameter_value_by_name
+      use entsor
+      implicit none
+
+      ! Local variables
+      character(len=300) :: zone_1_arbre
+      character(len=3) :: char_radius
+      character(len=4) :: char_height
+      character(len=4) :: char_longueur_rue
+      character(len=4) :: char_position_y
+      double precision lai
+
+      ! Get tree height
+      tree_height = notebook_parameter_value_by_name("hauteur_arbres")
+      if (tree_height .ge. 10.0) then
+         write(char_height,fmt='(F4.1)') tree_height
+      else
+         write(char_height,fmt='(F3.1)') tree_height
+      endif
+
+      ! Get tree position on Y axis
+      tree_position_y = notebook_parameter_value_by_name("position_y_arbres")
+      if (tree_position_y .ge. 10.0) then
+         write(char_position_y,fmt='(F4.1)') tree_position_y
+      else
+         write(char_position_y,fmt='(F3.1)') tree_position_y
+      endif
+
+      ! Get tree radius
+      tree_radius = notebook_parameter_value_by_name("rayon_arbres")
+      write(char_radius,fmt='(F3.1)') tree_radius
+
+      ! Get longueur rue
+      write(char_longueur_rue,fmt='(F4.1)') lz
+
+      ! Define zone arbre
+      tree_configuration = notebook_parameter_value_by_name("configuration_arbres")
+      if (tree_configuration == 1) then
+         if (nombre_arbre == 2) then
+             zone_1_arbre = "sphere[34.5, "//char_position_y//", "//char_height//", " &
+             //char_radius//"]"
+             zone_arbre = trim(zone_1_arbre) &
+             //" or sphere[48.0, "//char_position_y//", "//char_height//", " &
+             //char_radius//"]"
+         elseif (nombre_arbre == 1) then
+             zone_arbre = "sphere[41.25, "//char_position_y//", "//char_height//", " &
+             //char_radius//"]"
+             zone_1_arbre = trim(zone_arbre)
+         endif
+      	
+      elseif (tree_configuration == 2) then
+         if (nombre_arbre == 2) then
+             zone_1_arbre = "cylinder[34.5, 0.0, "//char_height &
+             //", 34.5, "//char_longueur_rue//", "//char_height//", "//char_radius//"]"
+             zone_arbre = trim(zone_1_arbre) &
+             //" or " &
+             //"cylinder[48.0, 0.0, "//char_height &
+             //", 48.0, "//char_longueur_rue//", "//char_height//", "//char_radius//"]"
+         elseif (nombre_arbre == 1) then
+             zone_arbre = "cylinder[41.25, 0.0, "//char_height &
+             //", 41.25, "//char_longueur_rue//", "//char_height//", "//char_radius//"]"
+             zone_1_arbre = trim(zone_arbre)
+         endif
+      endif
+
+      call calcul_volume_zone(trim(zone_arbre), volume_arbre)
+      call calcul_surface_projetee(trim(zone_1_arbre), surface_projetee_arbre)
+      surface_projetee_arbre = surface_projetee_arbre * nombre_arbre
+      surface_arbre = volume_arbre / lz
+
+      ! Compute alpha
+      lai = notebook_parameter_value_by_name("LAI") ! no unit
+      alpha = lai * surface_projetee_arbre / volume_arbre ! LAI * surface projetee / volume
+
+      ! Some output in the listing
+      write(nfecra,*) "**********************"
+      write(nfecra,*) "*** Zone arbre     ***"
+      write(nfecra,*) "**********************"
+      write(nfecra,*) "LAI : ", lai
+      write(nfecra,*) "Zone : ", trim(zone_arbre)
+      write(nfecra,*) "Zone_1_arbre : ", trim(zone_1_arbre)
+      write(nfecra,*) "Volume : ", volume_arbre
+      write(nfecra,*) "Surface projetée : ", surface_projetee_arbre
+      write(nfecra,*) "alpha : ", alpha
+      write(nfecra,*) "**********************"
+      write(nfecra,*) "*** Fin zone arbre ***"
+      write(nfecra,*) "**********************"
+      write(nfecra,*) ""
+
+      return
+   end subroutine initialize_arbre
+
+   !
+   ! Fonction pour initialiser les termes sources des voitures et des arbres
+   !
+   subroutine initialize_chem_source_term(meca_voiture, meca_arbre, emission_voiture, emission_arbre, depot_arbre)
+
+      use cs_c_bindings ! notebook_parameter_value_by_name
+
+      implicit none
+
+      logical, intent(in) :: meca_voiture, meca_arbre
+      logical, intent(in) :: emission_voiture, emission_arbre, depot_arbre
+
+      nombre_arbre = notebook_parameter_value_by_name("nombre_arbre")
+
+      if (meca_arbre .or. emission_arbre .or. depot_arbre) then
+         call initialize_arbre()
+      endif
+
+      return
+   end subroutine initialize_chem_source_term
+
+end module user_module
+
+!-------------------------------------------------------------------------------
+
+!> (DOXYGEN_SHOULD_SKIP_THIS) \endcond
+