File: compute_gaseous_chemistry.f90

package info (click to toggle)
code-saturne 4.3.3%2Brepack-1
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 77,992 kB
  • sloc: ansic: 281,257; f90: 122,305; python: 56,490; makefile: 3,915; xml: 3,285; cpp: 3,183; sh: 1,139; lex: 176; yacc: 101; sed: 16
file content (203 lines) | stat: -rw-r--r-- 5,791 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
!-------------------------------------------------------------------------------

! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2016 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 compute_gaseous_chemistry.f90
!> \brief Calls the rosenbrock resolution for atmospheric chemistry
!
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     dt            time step (per cell)
!_______________________________________________________________________________

subroutine compute_gaseous_chemistry ( dt )

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

use paramx
use numvar
use entsor
use optcal
use cstphy
use cstnum
use parall
use period
use pointe
use ppppar
use ppthch
use ppincl
use mesh
use field
use dimens
use atchem
use siream

implicit none

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

! Arguments

double precision dt(ncelet)

! Local Variables

integer iel,ii
double precision rom

!  Variables used for Rosenbrock resolution
double precision  dlconc(nespg)
double precision  source(nespg)
double precision  dchema(nespg)
double precision  conv_factor(nespg) ! conversion factors for reaction rates
double precision rk(nrg)
double precision dtc
integer ncycle
double precision dtrest

double precision, dimension(:), pointer :: crom
type(pmapper_double_r1), dimension(:), allocatable :: cvar_espg, cvara_espg

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

allocate(cvar_espg(nespg), cvara_espg(nespg))

call field_get_val_s(icrom, crom)

! Arrays of pointers containing the fields values for each species
! (loop on cells outside loop on species)
do ii = 1, nespg
  call field_get_val_s(ivarfl(isca(ii)), cvar_espg(ii)%p)
  call field_get_val_prev_s(ivarfl(isca(ii)), cvara_espg(ii)%p)
enddo

do iel = 1, ncel

  ! time step
  dtc = dt(iel)

  ! density
  rom = crom(iel)

  ! Filling working array rk
  do ii = 1, nrg
    rk(ii) = reacnum((ii-1)*ncel+iel)
  enddo

  do ii = 1, nespg
    conv_factor(chempoint(ii)) = rom*navo*(1.0d-12)/dmmk(ii)
    source(ii) = 0.0d0
  enddo

  if ((isepchemistry.eq.1).or.(ntcabs.eq.1)) then
    ! -----------------------------
    ! -- splitted Rosenbrock solver
    ! -----------------------------

    ! Filling working array dlconc with values at current time step
    do ii = 1, nespg
      dlconc(chempoint(ii)) = cvar_espg(ii)%p(iel)
    enddo

  else
    ! -----------------------------
    ! -- semi-coupled Rosenbrock solver
    ! -----------------------------

    ! Filling working array dlconc with values at previous time step
    do ii = 1, nespg
      dlconc(chempoint(ii)) = cvara_espg(ii)%p(iel)
    enddo

    ! Computation of C(Xn)
    if (ichemistry.eq.1) then
      call fexchem_1 (nespg,nrg,dlconc,rk,source,conv_factor,dchema)
    else if (ichemistry.eq.2) then
      call fexchem_2 (nespg,nrg,dlconc,rk,source,conv_factor,dchema)
    else if (ichemistry.eq.3) then
      if (iaerosol.eq.1) then
        call fexchem_siream (nespg,nrg,dlconc,rk,source,conv_factor,dchema)
      else
        call fexchem_3 (nespg,nrg,dlconc,rk,source,conv_factor,dchema)
      endif
    else if (ichemistry.eq.4) then
      call fexchem (nespg,nrg,dlconc,rk,source,conv_factor,dchema)
    endif

    ! Explicit contribution from dynamics as a source term:
    ! (X*-Xn)/dt(dynamics) - C(Xn). See usatch.f90
    ! The first nespg user scalars are supposed to be chemical species
    do ii = 1, nespg
      source(chempoint(ii)) = (cvar_espg(ii)%p(iel)-cvara_espg(ii)%p(iel))/dtc  &
                            - dchema(chempoint(ii))
    enddo

  endif ! End test isepchemistry

  ! Rosenbrock resoluion

  ! The maximum time step used for chemistry resolution is dtchemmax
  if (dtc.le.dtchemmax) then
    call roschem (dlconc,source,source,conv_factor,dtc,rk,rk)
  else
    ncycle = int(dtc/dtchemmax)
    dtrest = mod(dtc,dtchemmax)
    do ii = 1, ncycle
      call roschem (dlconc,source,source,conv_factor,dtchemmax,rk,rk)
    enddo
    call roschem (dlconc,source,source,conv_factor,dtrest,rk,rk)
  endif

  ! Update of values at current time step
  do ii = 1, nespg
    cvar_espg(ii)%p(iel) = dlconc(chempoint(ii))
  enddo

enddo

deallocate(cvar_espg, cvara_espg)

! TODO: clipping or not ?
! Clipping
do ii = 1, nespg
  call clpsca (ii)
enddo

!--------
! FORMATS
!--------
!----
! END
!----

return
end subroutine compute_gaseous_chemistry