← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/code-saturne/oneiric

« back to all changes in this revision

Viewing changes to src/cfbl/cfbsc2.f90

  • Committer: Bazaar Package Importer
  • Author(s): Sylvestre Ledru
  • Date: 2009-11-02 23:21:17 UTC
  • Revision ID: james.westby@ubuntu.com-20091102232117-9brxj2l5e33ie45a
Tags: upstream-2.0.0.beta2
Import upstream version 2.0.0.beta2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/cfbl/cfbsc2.f90
 
1
!-------------------------------------------------------------------------------
 
2
 
 
3
!     This file is part of the Code_Saturne Kernel, element of the
 
4
!     Code_Saturne CFD tool.
 
5
 
 
6
!     Copyright (C) 1998-2008 EDF S.A., France
 
7
 
 
8
!     contact: saturne-support@edf.fr
 
9
 
 
10
!     The Code_Saturne Kernel is free software; you can redistribute it
 
11
!     and/or modify it under the terms of the GNU General Public License
 
12
!     as published by the Free Software Foundation; either version 2 of
 
13
!     the License, or (at your option) any later version.
 
14
 
 
15
!     The Code_Saturne Kernel is distributed in the hope that it will be
 
16
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
 
17
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
!     GNU General Public License for more details.
 
19
 
 
20
!     You should have received a copy of the GNU General Public License
 
21
!     along with the Code_Saturne Kernel; if not, write to the
 
22
!     Free Software Foundation, Inc.,
 
23
!     51 Franklin St, Fifth Floor,
 
24
!     Boston, MA  02110-1301  USA
 
25
 
 
26
!-------------------------------------------------------------------------------
 
27
 
 
28
                  subroutine cfbsc2                               &
 
29
!================
 
30
 
 
31
 ( idbia0 , idbra0 ,                                              &
 
32
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
 
33
   nnod   , lndfac , lndfbr , ncelbr ,                            &
 
34
   nvar   , nscal  , nphas  ,                                     &
 
35
   nideve , nrdeve , nituse , nrtuse ,                            &
 
36
   ivar   , iconvp , idiffp , nswrgp , imligp , ircflp ,          &
 
37
   ischcp , isstpp , inc    , imrgra , iccocg , iifbru ,          &
 
38
   ipp    , iwarnp ,                                              &
 
39
   blencp , epsrgp , climgp , extrap ,                            &
 
40
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , ifrusb ,          &
 
41
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
 
42
   idevel , ituser , ia     ,                                     &
 
43
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
44
   pvar   , coefap , coefbp , cofafp , cofbfp ,                   &
 
45
   flumas , flumab , viscf  , viscb  ,                            &
 
46
   smbrp  ,                                                       &
 
47
   dpdx   , dpdy   , dpdz   , dpdxa  , dpdya  , dpdza  ,          &
 
48
   rdevel , rtuser , ra     )
 
49
 
 
50
!===============================================================================
 
51
! FONCTION :
 
52
! ---------
 
53
 
 
54
! CALCUL DU BILAN EXPLICITE DE LA VARIABLE PVAR (VITESSE,SCALAIRES)
 
55
 
 
56
!                   -- .                  ----->        -->
 
57
! SMBRP = SMBRP -  \   m   PVAR +Visc   ( grad PVAR )  . n
 
58
!                  /__  ij    ij     ij             ij    ij
 
59
!                   j
 
60
 
 
61
! ATTENTION : SMBRP DEJA INITIALISE AVANT L'APPEL A BILSCA
 
62
!            IL CONTIENT LES TERMES SOURCES EXPLICITES, ETC....
 
63
 
 
64
! BLENCP = 1 : PAS D'UPWIND EN DEHORS DU TEST DE PENTE
 
65
! BLENCP = 0 : UPWIND
 
66
! ISCHCP = 1 : CENTRE
 
67
! ISCHCP = 0 : SECOND ORDER
 
68
 
 
69
!-------------------------------------------------------------------------------
 
70
! Arguments
 
71
!__________________.____._____.________________________________________________.
 
72
!    nom           !type!mode !                   role                         !
 
73
!__________________!____!_____!________________________________________________!
 
74
! idbia0           ! e  ! <-- ! numero de la 1ere case libre dans ia           !
 
75
! idbra0           ! e  ! <-- ! numero de la 1ere case libre dans ra           !
 
76
! ndim             ! e  ! <-- ! dimension de l'espace                          !
 
77
! ncelet           ! e  ! <-- ! nombre d'elements halo compris                 !
 
78
! ncel             ! e  ! <-- ! nombre d'elements actifs                       !
 
79
! nfac             ! e  ! <-- ! nombre de faces internes                       !
 
80
! nfabor           ! e  ! <-- ! nombre de faces de bord                        !
 
81
! nfml             ! e  ! <-- ! nombre de familles d entites                   !
 
82
! nprfml           ! e  ! <-- ! nombre de proprietese des familles             !
 
83
! nnod             ! e  ! <-- ! nombre de sommets                              !
 
84
! lndfac           ! e  ! <-- ! longueur du tableau nodfac (optionnel          !
 
85
! lndfbr           ! e  ! <-- ! longueur du tableau nodfbr (optionnel          !
 
86
! ncelbr           ! e  ! <-- ! nombre d'elements ayant au moins une           !
 
87
!                  !    !     ! face de bord                                   !
 
88
! nvar             ! e  ! <-- ! nombre total de variables                      !
 
89
! nscal            ! e  ! <-- ! nombre total de scalaires                      !
 
90
! nphas            ! e  ! <-- ! nombre de phases                               !
 
91
! nideve nrdeve    ! e  ! <-- ! longueur de idevel rdevel                      !
 
92
! nituse nrtuse    ! e  ! <-- ! longueur de ituser rtuser                      !
 
93
! ivar             ! e  ! <-- ! numero de la variable                          !
 
94
! iconvp           ! e  ! <-- ! indicateur = 1 convection, 0 sinon             !
 
95
! idiffp           ! e  ! <-- ! indicateur = 1 diffusion , 0 sinon             !
 
96
! nswrgp           ! e  ! <-- ! nombre de sweep pour reconstruction            !
 
97
!                  !    !     !             des gradients                      !
 
98
! imligp           ! e  ! <-- ! methode de limitation du gradient              !
 
99
!                  !    !     !  < 0 pas de limitation                         !
 
100
!                  !    !     !  = 0 a partir des gradients voisins            !
 
101
!                  !    !     !  = 1 a partir du gradient moyen                !
 
102
! ircflp           ! e  ! <-- ! indicateur = 1 rec flux, 0 sinon               !
 
103
! ischcp           ! e  ! <-- ! indicateur = 1 centre , 0 2nd order            !
 
104
! isstpp           ! e  ! <-- ! indicateur = 1 sans test de pente              !
 
105
!                  !    !     !            = 0 avec test de pente              !
 
106
! inc              ! e  ! <-- ! indicateur = 0 resol sur increment             !
 
107
!                  !    !     !              1 sinon                           !
 
108
! imrgra           ! e  ! <-- ! indicateur = 0 gradrc 97                       !
 
109
!                  ! e  ! <-- !            = 1 gradmc 99                       !
 
110
! iccocg           ! e  ! <-- ! indicateur = 1 pour recalcul de cocg           !
 
111
!                  !    !     !              0 sinon                           !
 
112
! iifbru           ! e  ! <-- ! pointeur flux de bord rusanov                  !
 
113
! ipp              ! e  ! <-- ! numero de variable pour post                   !
 
114
! iwarnp           ! e  ! <-- ! niveau d'impression                            !
 
115
! blencp           ! r  ! <-- ! 1 - proportion d'upwind                        !
 
116
! epsrgp           ! r  ! <-- ! precision relative pour la                     !
 
117
!                  !    !     !  reconstruction des gradients 97               !
 
118
! climgp           ! r  ! <-- ! coef gradient*distance/ecart                   !
 
119
! extrap           ! r  ! <-- ! coef extrap gradient                           !
 
120
! ifacel           ! te ! <-- ! elements voisins d'une face interne            !
 
121
! (2, nfac)        !    !     !                                                !
 
122
! ifabor           ! te ! <-- ! element  voisin  d'une face de bord            !
 
123
! (nfabor)         !    !     !                                                !
 
124
! ifmfbr           ! te ! <-- ! numero de famille d'une face de bord           !
 
125
! (nfabor)         !    !     !                                                !
 
126
! ifmcel           ! te ! <-- ! numero de famille d'une cellule                !
 
127
! (ncelet)         !    !     !                                                !
 
128
! iprfml           ! te ! <-- ! proprietes d'une famille                       !
 
129
! nfml  ,nprfml    !    !     !                                                !
 
130
! ifrusb(nfabor    ! te ! <-- ! indicateur flux de rusanov                     !
 
131
! ipnfac           ! te ! <-- ! position du premier noeud de chaque            !
 
132
!   (lndfac)       !    !     !  face interne dans nodfac (optionnel)          !
 
133
! nodfac           ! te ! <-- ! connectivite faces internes/noeuds             !
 
134
!   (nfac+1)       !    !     !  (optionnel)                                   !
 
135
! ipnfbr           ! te ! <-- ! position du premier noeud de chaque            !
 
136
!   (lndfbr)       !    !     !  face de bord dans nodfbr (optionnel)          !
 
137
! nodfbr           ! te ! <-- ! connectivite faces de bord/noeuds              !
 
138
!   (nfabor+1)     !    !     !  (optionnel)                                   !
 
139
! idevel(nideve    ! te ! <-- ! tab entier complementaire developemt           !
 
140
! ituser(nituse    ! te ! <-- ! tab entier complementaire utilisateur          !
 
141
! ia(*)            ! te ! --- ! macro tableau entier                           !
 
142
! xyzcen           ! tr ! <-- ! point associes aux volumes de control          !
 
143
! (ndim,ncelet     !    !     !                                                !
 
144
! surfac           ! tr ! <-- ! vecteur surface des faces internes             !
 
145
! (ndim,nfac)      !    !     !                                                !
 
146
! surfbo           ! tr ! <-- ! vecteur surface des faces de bord              !
 
147
! (ndim,nfabor)    !    !     !                                                !
 
148
! cdgfac           ! tr ! <-- ! centre de gravite des faces internes           !
 
149
! (ndim,nfac)      !    !     !                                                !
 
150
! cdgfbo           ! tr ! <-- ! centre de gravite des faces de bord            !
 
151
! (ndim,nfabor)    !    !     !                                                !
 
152
! xyznod           ! tr ! <-- ! coordonnes des noeuds (optionnel)              !
 
153
! (ndim,nnod)      !    !     !                                                !
 
154
! volume           ! tr ! <-- ! volume d'un des ncelet elements                !
 
155
! (ncelet          !    !     !                                                !
 
156
! pvar (ncelet     ! tr ! <-- ! variable resolue (instant precedent)           !
 
157
! coefap, b        ! tr ! <-- ! tableaux des cond lim pour p                   !
 
158
!   (nfabor)       !    !     !  sur la normale a la face de bord              !
 
159
! cofafp, b        ! tr ! <-- ! tableaux des cond lim pour le flux de          !
 
160
!   (nfabor)       !    !     !  diffusion de p                                !
 
161
! flumas(nfac)     ! tr ! <-- ! flux de masse aux faces internes               !
 
162
! flumab(nfabor    ! tr ! <-- ! flux de masse aux faces de bord                !
 
163
! viscf (nfac)     ! tr ! <-- ! visc*surface/dist aux faces internes           !
 
164
!                  !    !     !  pour second membre                            !
 
165
! viscb (nfabor    ! tr ! <-- ! visc*surface/dist aux faces de bord            !
 
166
!                  !    !     !  pour second membre                            !
 
167
! smbrp(ncelet     ! tr ! <-- ! bilan au second membre                         !
 
168
! dpdx,y,z         ! tr ! --- ! tableau de travail pour le grad de p           !
 
169
!    (ncelet)      !    !     !                                                !
 
170
! dpdxa,ya,za      ! tr ! --- ! tableau de travail pour le grad de p           !
 
171
!    (ncelet)      !    !     !  avec decentrement amont                       !
 
172
! rdevel(nrdeve    ! tr ! <-- ! tab reel complementaire developemt             !
 
173
! rtuser(nrtuse    ! tr ! <-- ! tab reel complementaire utilisateur            !
 
174
! ra(*)            ! tr ! --- ! macro tableau reel                             !
 
175
!__________________!____!_____!________________________________________________!
 
176
 
 
177
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
178
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
179
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
180
!            --- tableau de travail
 
181
!===============================================================================
 
182
 
 
183
implicit none
 
184
 
 
185
!===============================================================================
 
186
!     DONNEES EN COMMON
 
187
!===============================================================================
 
188
 
 
189
include "paramx.h"
 
190
include "pointe.h"
 
191
include "vector.h"
 
192
include "entsor.h"
 
193
include "period.h"
 
194
include "parall.h"
 
195
 
 
196
!===============================================================================
 
197
 
 
198
! Arguments
 
199
 
 
200
integer          idbia0 , idbra0
 
201
integer          ndim   , ncelet , ncel   , nfac   , nfabor
 
202
integer          nfml   , nprfml
 
203
integer          nnod   , lndfac , lndfbr , ncelbr
 
204
integer          nvar   , nscal  , nphas
 
205
integer          nideve , nrdeve , nituse , nrtuse
 
206
integer          ivar   , iconvp , idiffp , nswrgp , imligp
 
207
integer          ircflp , ischcp , isstpp
 
208
integer          inc    , imrgra , iccocg , iifbru
 
209
integer          iwarnp , ipp
 
210
double precision blencp , epsrgp , climgp, extrap
 
211
 
 
212
integer          ifacel(2,nfac) , ifabor(nfabor)
 
213
integer          ifmfbr(nfabor) , ifmcel(ncelet)
 
214
integer          iprfml(nfml,nprfml)
 
215
integer          ifrusb(nfabor)
 
216
integer          ipnfac(nfac+1), nodfac(lndfac)
 
217
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
 
218
integer          idevel(nideve), ituser(nituse)
 
219
integer          ia(*)
 
220
 
 
221
double precision xyzcen(ndim,ncelet)
 
222
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
 
223
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
 
224
double precision xyznod(ndim,nnod), volume(ncelet)
 
225
double precision pvar (ncelet), coefap(nfabor), coefbp(nfabor)
 
226
double precision                cofafp(nfabor), cofbfp(nfabor)
 
227
double precision flumas(nfac), flumab(nfabor)
 
228
double precision viscf (nfac), viscb (nfabor)
 
229
double precision smbrp(ncelet)
 
230
double precision dpdx (ncelet),dpdy (ncelet),dpdz (ncelet)
 
231
double precision dpdxa(ncelet),dpdya(ncelet),dpdza(ncelet)
 
232
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
 
233
 
 
234
! VARIABLES LOCALES
 
235
 
 
236
character*80     chaine
 
237
character*8      cnom
 
238
integer          idebia, idebra
 
239
integer          ifac,ii,jj,infac,iel,iupwin, iij, iii, iok
 
240
integer          itenso, idimte, iphydp
 
241
integer          iiu(nphsmx),iiv(nphsmx),iiw(nphsmx)
 
242
integer          iitytu(nphsmx)
 
243
integer          iir11(nphsmx),iir22(nphsmx),iir33(nphsmx)
 
244
integer          iir12(nphsmx),iir13(nphsmx),iir23(nphsmx)
 
245
double precision pfac,pfacd,pip,pjp,flui,fluj,flux
 
246
double precision difx,dify,difz,djfx,djfy,djfz,pif,pjf
 
247
double precision testi,testj,testij
 
248
double precision dpxf,dpyf,dpzf
 
249
double precision dcc, ddi, ddj, tesqck
 
250
double precision dijpfx, dijpfy, dijpfz
 
251
double precision diipfx, diipfy, diipfz
 
252
double precision djjpfx, djjpfy, djjpfz
 
253
double precision diipbx, diipby, diipbz
 
254
double precision pond, dist, surfan
 
255
double precision pfac1, pfac2, pfac3, unsvol
 
256
 
 
257
!===============================================================================
 
258
 
 
259
!===============================================================================
 
260
! 1.  INITIALISATION
 
261
!===============================================================================
 
262
 
 
263
idebia = idbia0
 
264
idebra = idbra0
 
265
 
 
266
chaine = nomvar(ipp)
 
267
cnom   = chaine(1:8)
 
268
 
 
269
if(iwarnp.ge.2) then
 
270
  if (ischcp.eq.1) then
 
271
    WRITE(NFECRA,1000)CNOM,'    CENTRE ',                         &
 
272
                                              (1.d0-blencp)*100.d0
 
273
  else
 
274
    WRITE(NFECRA,1000)CNOM,' 2ND ORDER ',                         &
 
275
                                              (1.d0-blencp)*100.d0
 
276
  endif
 
277
endif
 
278
 
 
279
iupwin = 0
 
280
if(blencp.eq.0.d0) iupwin = 1
 
281
 
 
282
!===============================================================================
 
283
! 2.  CALCUL DU BILAN AVEC TECHNIQUE DE RECONSTRUCTION
 
284
!===============================================================================
 
285
 
 
286
! ======================================================================
 
287
! ---> CALCUL DU GRADIENT DE P
 
288
! ======================================================================
 
289
!    DPDX sert a la fois pour la reconstruction des flux et pour le test
 
290
!    de pente. On doit donc le calculer :
 
291
!        - quand on a de la diffusion et qu'on reconstruit les flux
 
292
!        - quand on a de la convection SOLU
 
293
!        - quand on a de la convection, qu'on n'est pas en upwind pur
 
294
!          et qu'on reconstruit les flux
 
295
!        - quand on a de la convection, qu'on n'est pas en upwind pur
 
296
!          et qu'on n'a pas shunte le test de pente
 
297
 
 
298
if( (idiffp.ne.0 .and. ircflp.eq.1) .or.                          &
 
299
    (iconvp.ne.0 .and. iupwin.eq.0 .and.                          &
 
300
       (ischcp.eq.0 .or. ircflp.eq.1 .or. isstpp.eq.0)) ) then
 
301
 
 
302
  iphydp = 0
 
303
  call grdcel                                                     &
 
304
  !==========
 
305
 ( idebia , idebra ,                                              &
 
306
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
 
307
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
 
308
   nideve , nrdeve , nituse , nrtuse ,                            &
 
309
   ivar   , imrgra , inc    , iccocg , nswrgp , imligp ,  iphydp ,&
 
310
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
 
311
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
 
312
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
 
313
   idevel , ituser , ia     ,                                     &
 
314
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
315
   dpdxa  , dpdxa  , dpdxa  ,                                     &
 
316
   pvar   , coefap , coefbp ,                                     &
 
317
   dpdx   , dpdy   , dpdz   ,                                     &
 
318
!        ------   ------   ------
 
319
   dpdxa  , dpdya  , dpdza  ,                                     &
 
320
   rdevel , rtuser , ra     )
 
321
 
 
322
else
 
323
  do iel = 1, ncelet
 
324
    dpdx(iel) = 0.d0
 
325
    dpdy(iel) = 0.d0
 
326
    dpdz(iel) = 0.d0
 
327
  enddo
 
328
endif
 
329
 
 
330
 
 
331
! ======================================================================
 
332
! ---> CALCUL DU GRADIENT DECENTRE DPDXA, DPDYA, DPDZA POUR TST DE PENTE
 
333
! ======================================================================
 
334
 
 
335
do iel = 1, ncelet
 
336
  dpdxa(iel) = 0.d0
 
337
  dpdya(iel) = 0.d0
 
338
  dpdza(iel) = 0.d0
 
339
enddo
 
340
 
 
341
if( iconvp.gt.0.and.iupwin.eq.0.and.isstpp.eq.0 ) then
 
342
 
 
343
  if (ivecti.eq.1) then
 
344
 
 
345
!CDIR NODEP
 
346
    do ifac = 1, nfac
 
347
 
 
348
      ii = ifacel(1,ifac)
 
349
      jj = ifacel(2,ifac)
 
350
 
 
351
      difx = cdgfac(1,ifac) - xyzcen(1,ii)
 
352
      dify = cdgfac(2,ifac) - xyzcen(2,ii)
 
353
      difz = cdgfac(3,ifac) - xyzcen(3,ii)
 
354
      djfx = cdgfac(1,ifac) - xyzcen(1,jj)
 
355
      djfy = cdgfac(2,ifac) - xyzcen(2,jj)
 
356
      djfz = cdgfac(3,ifac) - xyzcen(3,jj)
 
357
 
 
358
      pif = pvar(ii) +difx*dpdx(ii)+dify*dpdy(ii)+difz*dpdz(ii)
 
359
      pjf = pvar(jj) +djfx*dpdx(jj)+djfy*dpdy(jj)+djfz*dpdz(jj)
 
360
 
 
361
      pfac = pjf
 
362
      if( flumas(ifac ).gt.0.d0 ) pfac = pif
 
363
 
 
364
      pfac1 = pfac*surfac(1,ifac )
 
365
      pfac2 = pfac*surfac(2,ifac )
 
366
      pfac3 = pfac*surfac(3,ifac )
 
367
 
 
368
      dpdxa(ii) = dpdxa(ii) +pfac1
 
369
      dpdya(ii) = dpdya(ii) +pfac2
 
370
      dpdza(ii) = dpdza(ii) +pfac3
 
371
 
 
372
      dpdxa(jj) = dpdxa(jj) -pfac1
 
373
      dpdya(jj) = dpdya(jj) -pfac2
 
374
      dpdza(jj) = dpdza(jj) -pfac3
 
375
 
 
376
    enddo
 
377
 
 
378
  else
 
379
 
 
380
! VECTORISATION NON FORCEE
 
381
    do ifac = 1, nfac
 
382
 
 
383
      ii = ifacel(1,ifac)
 
384
      jj = ifacel(2,ifac)
 
385
 
 
386
      difx = cdgfac(1,ifac) - xyzcen(1,ii)
 
387
      dify = cdgfac(2,ifac) - xyzcen(2,ii)
 
388
      difz = cdgfac(3,ifac) - xyzcen(3,ii)
 
389
      djfx = cdgfac(1,ifac) - xyzcen(1,jj)
 
390
      djfy = cdgfac(2,ifac) - xyzcen(2,jj)
 
391
      djfz = cdgfac(3,ifac) - xyzcen(3,jj)
 
392
 
 
393
      pif = pvar(ii) +difx*dpdx(ii)+dify*dpdy(ii)+difz*dpdz(ii)
 
394
      pjf = pvar(jj) +djfx*dpdx(jj)+djfy*dpdy(jj)+djfz*dpdz(jj)
 
395
 
 
396
      pfac = pjf
 
397
      if( flumas(ifac ).gt.0.d0 ) pfac = pif
 
398
 
 
399
      pfac1 = pfac*surfac(1,ifac )
 
400
      pfac2 = pfac*surfac(2,ifac )
 
401
      pfac3 = pfac*surfac(3,ifac )
 
402
 
 
403
      dpdxa(ii) = dpdxa(ii) +pfac1
 
404
      dpdya(ii) = dpdya(ii) +pfac2
 
405
      dpdza(ii) = dpdza(ii) +pfac3
 
406
 
 
407
      dpdxa(jj) = dpdxa(jj) -pfac1
 
408
      dpdya(jj) = dpdya(jj) -pfac2
 
409
      dpdza(jj) = dpdza(jj) -pfac3
 
410
 
 
411
    enddo
 
412
 
 
413
  endif
 
414
 
 
415
  if (ivectb.eq.1) then
 
416
 
 
417
!CDIR NODEP
 
418
    do ifac = 1, nfabor
 
419
      ii = ifabor(ifac )
 
420
      iii = idiipb-1+3*(ifac-1)
 
421
      diipbx = ra(iii+1)
 
422
      diipby = ra(iii+2)
 
423
      diipbz = ra(iii+3)
 
424
      pfac = inc*coefap(ifac )                                    &
 
425
            +coefbp(ifac )*(pvar(ii)+diipbx*dpdx(ii)              &
 
426
                    +diipby*dpdy(ii)+diipbz*dpdz(ii) )
 
427
      dpdxa(ii) = dpdxa(ii) +pfac*surfbo(1,ifac )
 
428
      dpdya(ii) = dpdya(ii) +pfac*surfbo(2,ifac )
 
429
      dpdza(ii) = dpdza(ii) +pfac*surfbo(3,ifac )
 
430
    enddo
 
431
 
 
432
  else
 
433
 
 
434
! VECTORISATION NON FORCEE
 
435
    do ifac =1,nfabor
 
436
      ii = ifabor(ifac )
 
437
      iii = idiipb-1+3*(ifac-1)
 
438
      diipbx = ra(iii+1)
 
439
      diipby = ra(iii+2)
 
440
      diipbz = ra(iii+3)
 
441
      pfac = inc*coefap(ifac )                                    &
 
442
            +coefbp(ifac )*(pvar(ii)+diipbx*dpdx(ii)              &
 
443
                    +diipby*dpdy(ii)+diipbz*dpdz(ii) )
 
444
      dpdxa(ii) = dpdxa(ii) +pfac*surfbo(1,ifac )
 
445
      dpdya(ii) = dpdya(ii) +pfac*surfbo(2,ifac )
 
446
      dpdza(ii) = dpdza(ii) +pfac*surfbo(3,ifac )
 
447
    enddo
 
448
 
 
449
  endif
 
450
 
 
451
  do iel = 1, ncel
 
452
    unsvol = 1.d0/volume(iel)
 
453
    dpdxa(iel) = dpdxa(iel)*unsvol
 
454
    dpdya(iel) = dpdya(iel)*unsvol
 
455
    dpdza(iel) = dpdza(iel)*unsvol
 
456
  enddo
 
457
 
 
458
!     TRAITEMENT DU PARALLELISME
 
459
 
 
460
  if(irangp.ge.0) then
 
461
    call parcom (dpdxa)
 
462
    !==========
 
463
    call parcom (dpdya)
 
464
    !==========
 
465
    call parcom (dpdza)
 
466
    !==========
 
467
  endif
 
468
 
 
469
! TRAITEMENT DE LA PERIODICITE
 
470
 
 
471
!  On echange pour la translation
 
472
!   pour la rotation, on prend le gradient simple (pas de temps precedent)
 
473
 
 
474
  if(iperio.eq.1) then
 
475
 
 
476
!        Pour les rotations
 
477
!          avec la vitesse et les tensions de Reynolds,
 
478
!          on utilise la valeur du gradient simple (PERING) a defaut de mieux.
 
479
!        Dans les autres cas, on echange DPDXA
 
480
 
 
481
!        On recupere d'abord certains COMMON necessaires a PERING
 
482
 
 
483
    call pergra                                                   &
 
484
    !==========
 
485
  ( nphsmx , nphas  ,                                             &
 
486
    iiu    , iiv    , iiw    ,                                    &
 
487
    iitytu ,                                                      &
 
488
    iir11  , iir22  , iir33  , iir12  , iir13  , iir23  )
 
489
 
 
490
    call pering                                                   &
 
491
    !==========
 
492
  ( nphas  , ivar   ,                                             &
 
493
    idimte , itenso , iperot , iguper , igrper ,                  &
 
494
    iiu    , iiv    , iiw    , iitytu ,                           &
 
495
    iir11  , iir22  , iir33  , iir12  , iir13  , iir23  ,         &
 
496
    dpdxa  , dpdya  , dpdza  ,                                    &
 
497
    ra(idudxy) , ra(idrdxy)  )
 
498
 
 
499
    call percom                                                   &
 
500
    !==========
 
501
  ( idimte , itenso ,                                             &
 
502
    dpdxa  , dpdxa  , dpdxa ,                                     &
 
503
    dpdya  , dpdya  , dpdya ,                                     &
 
504
    dpdza  , dpdza  , dpdza )
 
505
  endif
 
506
 
 
507
endif
 
508
 
 
509
 
 
510
! ======================================================================
 
511
! ---> ASSEMBLAGE A PARTIR DES FACETTES FLUIDES
 
512
! ======================================================================
 
513
 
 
514
infac = 0
 
515
 
 
516
if(ncelet.gt.ncel) then
 
517
  do iel = ncel+1, ncelet
 
518
    smbrp(iel) = 0.d0
 
519
  enddo
 
520
endif
 
521
 
 
522
 
 
523
!  --> FLUX UPWIND PUR
 
524
!  =====================
 
525
 
 
526
if(iupwin.eq.1) then
 
527
 
 
528
  if (ivecti.eq.1) then
 
529
 
 
530
!CDIR NODEP
 
531
    do ifac = 1, nfac
 
532
 
 
533
      ii = ifacel(1,ifac)
 
534
      jj = ifacel(2,ifac)
 
535
 
 
536
      iij = idijpf-1+3*(ifac-1)
 
537
      dijpfx = ra(iij+1)
 
538
      dijpfy = ra(iij+2)
 
539
      dijpfz = ra(iij+3)
 
540
 
 
541
      pond   = ra(ipond-1+ifac)
 
542
 
 
543
! ON RECALCULE A CE NIVEAU II' ET JJ'
 
544
 
 
545
      diipfx = cdgfac(1,ifac) - (xyzcen(1,ii)+                    &
 
546
               (1.d0-pond) * dijpfx)
 
547
      diipfy = cdgfac(2,ifac) - (xyzcen(2,ii)+                    &
 
548
               (1.d0-pond) * dijpfy)
 
549
      diipfz = cdgfac(3,ifac) - (xyzcen(3,ii)+                    &
 
550
               (1.d0-pond) * dijpfz)
 
551
      djjpfx = cdgfac(1,ifac) -  xyzcen(1,jj)+                    &
 
552
                   pond  * dijpfx
 
553
      djjpfy = cdgfac(2,ifac) -  xyzcen(2,jj)+                    &
 
554
                   pond  * dijpfy
 
555
      djjpfz = cdgfac(3,ifac) -  xyzcen(3,jj)+                    &
 
556
                   pond  * dijpfz
 
557
 
 
558
      dpxf = 0.5d0*(dpdx(ii) + dpdx(jj))
 
559
      dpyf = 0.5d0*(dpdy(ii) + dpdy(jj))
 
560
      dpzf = 0.5d0*(dpdz(ii) + dpdz(jj))
 
561
 
 
562
!     reconstruction uniquement si IRCFLP = 1
 
563
      pip = pvar(ii)                                              &
 
564
           + ircflp*(dpxf*diipfx+dpyf*diipfy+dpzf*diipfz)
 
565
      pjp = pvar(jj)                                              &
 
566
           + ircflp*(dpxf*djjpfx+dpyf*djjpfy+dpzf*djjpfz)
 
567
 
 
568
      flui = 0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
569
      fluj = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
570
 
 
571
      pif = pvar(ii)
 
572
      pjf = pvar(jj)
 
573
      infac = infac+1
 
574
 
 
575
      flux = iconvp*( flui*pif +fluj*pjf )                        &
 
576
           + idiffp*viscf(ifac)*( pip -pjp )
 
577
 
 
578
      smbrp(ii) = smbrp(ii) -flux
 
579
      smbrp(jj) = smbrp(jj) +flux
 
580
 
 
581
    enddo
 
582
 
 
583
  else
 
584
 
 
585
! VECTORISATION NON FORCEE
 
586
    do ifac = 1, nfac
 
587
 
 
588
      ii = ifacel(1,ifac)
 
589
      jj = ifacel(2,ifac)
 
590
 
 
591
      iij = idijpf-1+3*(ifac-1)
 
592
      dijpfx = ra(iij+1)
 
593
      dijpfy = ra(iij+2)
 
594
      dijpfz = ra(iij+3)
 
595
 
 
596
      pond   = ra(ipond-1+ifac)
 
597
 
 
598
! ON RECALCULE A CE NIVEAU II' ET JJ'
 
599
 
 
600
      diipfx = cdgfac(1,ifac) - (xyzcen(1,ii)+                    &
 
601
               (1.d0-pond) * dijpfx)
 
602
      diipfy = cdgfac(2,ifac) - (xyzcen(2,ii)+                    &
 
603
               (1.d0-pond) * dijpfy)
 
604
      diipfz = cdgfac(3,ifac) - (xyzcen(3,ii)+                    &
 
605
               (1.d0-pond) * dijpfz)
 
606
      djjpfx = cdgfac(1,ifac) -  xyzcen(1,jj)+                    &
 
607
                   pond  * dijpfx
 
608
      djjpfy = cdgfac(2,ifac) -  xyzcen(2,jj)+                    &
 
609
                   pond  * dijpfy
 
610
      djjpfz = cdgfac(3,ifac) -  xyzcen(3,jj)+                    &
 
611
                   pond  * dijpfz
 
612
 
 
613
      dpxf = 0.5d0*(dpdx(ii) + dpdx(jj))
 
614
      dpyf = 0.5d0*(dpdy(ii) + dpdy(jj))
 
615
      dpzf = 0.5d0*(dpdz(ii) + dpdz(jj))
 
616
 
 
617
      pip = pvar(ii)                                              &
 
618
           + ircflp*(dpxf*diipfx+dpyf*diipfy+dpzf*diipfz)
 
619
      pjp = pvar(jj)                                              &
 
620
           + ircflp*(dpxf*djjpfx+dpyf*djjpfy+dpzf*djjpfz)
 
621
 
 
622
      flui = 0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
623
      fluj = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
624
 
 
625
      pif = pvar(ii)
 
626
      pjf = pvar(jj)
 
627
      infac = infac+1
 
628
 
 
629
      flux = iconvp*( flui*pif +fluj*pjf )                        &
 
630
           + idiffp*viscf(ifac)*( pip -pjp )
 
631
 
 
632
      smbrp(ii) = smbrp(ii) -flux
 
633
      smbrp(jj) = smbrp(jj) +flux
 
634
 
 
635
    enddo
 
636
 
 
637
  endif
 
638
 
 
639
 
 
640
!  --> FLUX SANS TEST DE PENTE
 
641
!  ============================
 
642
 
 
643
elseif(isstpp.eq.1) then
 
644
 
 
645
  if (ivecti.eq.1) then
 
646
 
 
647
    iok = 0
 
648
!CDIR NODEP
 
649
    do ifac = 1, nfac
 
650
 
 
651
      ii = ifacel(1,ifac)
 
652
      jj = ifacel(2,ifac)
 
653
 
 
654
      iij = idijpf-1+3*(ifac-1)
 
655
 
 
656
      dijpfx = ra(iij+1)
 
657
      dijpfy = ra(iij+2)
 
658
      dijpfz = ra(iij+3)
 
659
 
 
660
      pond   = ra(ipond-1+ifac)
 
661
 
 
662
! ON RECALCULE A CE NIVEAU II' ET JJ'
 
663
 
 
664
 
 
665
      diipfx = cdgfac(1,ifac) - (xyzcen(1,ii)+                    &
 
666
               (1.d0-pond) * dijpfx)
 
667
      diipfy = cdgfac(2,ifac) - (xyzcen(2,ii)+                    &
 
668
               (1.d0-pond) * dijpfy)
 
669
      diipfz = cdgfac(3,ifac) - (xyzcen(3,ii)+                    &
 
670
               (1.d0-pond) * dijpfz)
 
671
      djjpfx = cdgfac(1,ifac) -  xyzcen(1,jj)+                    &
 
672
                   pond  * dijpfx
 
673
      djjpfy = cdgfac(2,ifac) -  xyzcen(2,jj)+                    &
 
674
                   pond  * dijpfy
 
675
      djjpfz = cdgfac(3,ifac) -  xyzcen(3,jj)+                    &
 
676
                   pond  * dijpfz
 
677
 
 
678
      dpxf = 0.5d0*(dpdx(ii) + dpdx(jj))
 
679
      dpyf = 0.5d0*(dpdy(ii) + dpdy(jj))
 
680
      dpzf = 0.5d0*(dpdz(ii) + dpdz(jj))
 
681
 
 
682
 
 
683
      pip = pvar(ii)                                              &
 
684
           + ircflp*(dpxf*diipfx+dpyf*diipfy+dpzf*diipfz)
 
685
      pjp = pvar(jj)                                              &
 
686
           + ircflp*(dpxf*djjpfx+dpyf*djjpfy+dpzf*djjpfz)
 
687
 
 
688
      flui = 0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
689
      fluj = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
690
 
 
691
 
 
692
!         CENTRE
 
693
!        --------
 
694
 
 
695
      if (ischcp.eq.1) then
 
696
 
 
697
        pif = pond*pip +(1.d0-pond)*pjp
 
698
        pjf = pif
 
699
 
 
700
 
 
701
!         SECOND ORDER
 
702
!        --------------
 
703
 
 
704
      elseif(ischcp.eq.0) then
 
705
 
 
706
        difx = cdgfac(1,ifac) - xyzcen(1,ii)
 
707
        dify = cdgfac(2,ifac) - xyzcen(2,ii)
 
708
        difz = cdgfac(3,ifac) - xyzcen(3,ii)
 
709
        djfx = cdgfac(1,ifac) - xyzcen(1,jj)
 
710
        djfy = cdgfac(2,ifac) - xyzcen(2,jj)
 
711
        djfz = cdgfac(3,ifac) - xyzcen(3,jj)
 
712
 
 
713
!     on laisse la reconstruction de PIF et PJF meme si IRCFLP=0
 
714
!     sinon cela revient a faire de l'upwind
 
715
        pif = pvar(ii)                                            &
 
716
             +difx*dpdx(ii)+dify*dpdy(ii)+difz*dpdz(ii)
 
717
        pjf = pvar(jj)                                            &
 
718
             +djfx*dpdx(jj)+djfy*dpdy(jj)+djfz*dpdz(jj)
 
719
 
 
720
      else
 
721
        write(nfecra,9000)ischcp
 
722
        iok = 1
 
723
      endif
 
724
 
 
725
 
 
726
!        BLENDING
 
727
!       ----------
 
728
 
 
729
      pif = blencp*pif+(1.d0-blencp)*pvar(ii)
 
730
      pjf = blencp*pjf+(1.d0-blencp)*pvar(jj)
 
731
 
 
732
 
 
733
!        FLUX
 
734
!       ------
 
735
 
 
736
      flux = iconvp*( flui*pif +fluj*pjf )                        &
 
737
           + idiffp*viscf(ifac)*( pip -pjp )
 
738
 
 
739
 
 
740
!        ASSEMBLAGE
 
741
!       ------------
 
742
 
 
743
      smbrp(ii) = smbrp(ii) -flux
 
744
      smbrp(jj) = smbrp(jj) +flux
 
745
 
 
746
    enddo
 
747
!        Le call csexit ne doit pas etre dans la boucle, sinon
 
748
!         le VPP la devectorise.
 
749
    if(iok.ne.0) then
 
750
      call csexit (1)
 
751
    endif
 
752
 
 
753
  else
 
754
 
 
755
    iok = 0
 
756
! VECTORISATION NON FORCEE
 
757
    do ifac = 1, nfac
 
758
 
 
759
      ii = ifacel(1,ifac)
 
760
      jj = ifacel(2,ifac)
 
761
 
 
762
      iij = idijpf-1+3*(ifac-1)
 
763
 
 
764
      dijpfx = ra(iij+1)
 
765
      dijpfy = ra(iij+2)
 
766
      dijpfz = ra(iij+3)
 
767
 
 
768
      pond   = ra(ipond-1+ifac)
 
769
 
 
770
! ON RECALCULE II' ET JJ'
 
771
 
 
772
      diipfx = cdgfac(1,ifac) - (xyzcen(1,ii)+                    &
 
773
               (1.d0-pond) * dijpfx)
 
774
      diipfy = cdgfac(2,ifac) - (xyzcen(2,ii)+                    &
 
775
               (1.d0-pond) * dijpfy)
 
776
      diipfz = cdgfac(3,ifac) - (xyzcen(3,ii)+                    &
 
777
               (1.d0-pond) * dijpfz)
 
778
      djjpfx = cdgfac(1,ifac) -  xyzcen(1,jj)+                    &
 
779
                   pond  * dijpfx
 
780
      djjpfy = cdgfac(2,ifac) -  xyzcen(2,jj)+                    &
 
781
                   pond  * dijpfy
 
782
      djjpfz = cdgfac(3,ifac) -  xyzcen(3,jj)+                    &
 
783
                   pond  * dijpfz
 
784
 
 
785
      dpxf = 0.5d0*(dpdx(ii) + dpdx(jj))
 
786
      dpyf = 0.5d0*(dpdy(ii) + dpdy(jj))
 
787
      dpzf = 0.5d0*(dpdz(ii) + dpdz(jj))
 
788
 
 
789
      pip = pvar(ii)                                              &
 
790
           + ircflp*(dpxf*diipfx+dpyf*diipfy+dpzf*diipfz)
 
791
      pjp = pvar(jj)                                              &
 
792
           + ircflp*(dpxf*djjpfx+dpyf*djjpfy+dpzf*djjpfz)
 
793
 
 
794
      flui = 0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
795
      fluj = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
796
 
 
797
 
 
798
!         CENTRE
 
799
!        --------
 
800
 
 
801
      if (ischcp.eq.1) then
 
802
 
 
803
        pif = pond*pip +(1.d0-pond)*pjp
 
804
        pjf = pif
 
805
 
 
806
 
 
807
!         SECOND ORDER
 
808
!        --------------
 
809
 
 
810
      elseif(ischcp.eq.0) then
 
811
 
 
812
        difx = cdgfac(1,ifac) - xyzcen(1,ii)
 
813
        dify = cdgfac(2,ifac) - xyzcen(2,ii)
 
814
        difz = cdgfac(3,ifac) - xyzcen(3,ii)
 
815
        djfx = cdgfac(1,ifac) - xyzcen(1,jj)
 
816
        djfy = cdgfac(2,ifac) - xyzcen(2,jj)
 
817
        djfz = cdgfac(3,ifac) - xyzcen(3,jj)
 
818
 
 
819
!     on laisse la reconstruction de PIF et PJF meme si IRCFLP=0
 
820
!     sinon cela revient a faire de l'upwind
 
821
        pif = pvar(ii)                                            &
 
822
             +difx*dpdx(ii)+dify*dpdy(ii)+difz*dpdz(ii)
 
823
        pjf = pvar(jj)                                            &
 
824
             +djfx*dpdx(jj)+djfy*dpdy(jj)+djfz*dpdz(jj)
 
825
 
 
826
      else
 
827
        write(nfecra,9000)ischcp
 
828
        iok = 1
 
829
      endif
 
830
 
 
831
 
 
832
!        BLENDING
 
833
!       ----------
 
834
 
 
835
      pif = blencp*pif+(1.d0-blencp)*pvar(ii)
 
836
      pjf = blencp*pjf+(1.d0-blencp)*pvar(jj)
 
837
 
 
838
 
 
839
!        FLUX
 
840
!       ------
 
841
 
 
842
      flux = iconvp*( flui*pif +fluj*pjf )                        &
 
843
           + idiffp*viscf(ifac)*( pip -pjp )
 
844
 
 
845
 
 
846
!        ASSEMBLAGE
 
847
!       ------------
 
848
 
 
849
      smbrp(ii) = smbrp(ii) -flux
 
850
      smbrp(jj) = smbrp(jj) +flux
 
851
 
 
852
    enddo
 
853
!        Le call csexit ne doit pas etre dans la boucle, sinon
 
854
!         le VPP la devectorise (pour la boucle non vectorisee forcee,
 
855
!         ce n'est pas grave, c'est juste pour recopier exactement
 
856
!         la precedente)
 
857
    if(iok.ne.0) then
 
858
      call csexit (1)
 
859
    endif
 
860
 
 
861
  endif
 
862
 
 
863
 
 
864
 
 
865
 
 
866
!  --> FLUX AVEC TEST DE PENTE (separe pour vectorisation)
 
867
!  =============================
 
868
 
 
869
else
 
870
 
 
871
  if (ivecti.eq.1) then
 
872
 
 
873
    iok = 0
 
874
!CDIR NODEP
 
875
    do ifac = 1, nfac
 
876
 
 
877
      ii = ifacel(1,ifac)
 
878
      jj = ifacel(2,ifac)
 
879
 
 
880
      iij = idijpf-1+3*(ifac-1)
 
881
 
 
882
      dijpfx = ra(iij+1)
 
883
      dijpfy = ra(iij+2)
 
884
      dijpfz = ra(iij+3)
 
885
 
 
886
      pond   = ra(ipond-1+ifac)
 
887
      dist   = ra(idist-1+ifac)
 
888
      surfan = ra(isrfan-1+ifac)
 
889
 
 
890
! ON RECALCULE II' ET JJ'
 
891
 
 
892
      diipfx = cdgfac(1,ifac) - (xyzcen(1,ii)+                    &
 
893
               (1.d0-pond) * dijpfx)
 
894
      diipfy = cdgfac(2,ifac) - (xyzcen(2,ii)+                    &
 
895
               (1.d0-pond) * dijpfy)
 
896
      diipfz = cdgfac(3,ifac) - (xyzcen(3,ii)+                    &
 
897
               (1.d0-pond) * dijpfz)
 
898
      djjpfx = cdgfac(1,ifac) -  xyzcen(1,jj)+                    &
 
899
                   pond  * dijpfx
 
900
      djjpfy = cdgfac(2,ifac) -  xyzcen(2,jj)+                    &
 
901
                   pond  * dijpfy
 
902
      djjpfz = cdgfac(3,ifac) -  xyzcen(3,jj)+                    &
 
903
                   pond  * dijpfz
 
904
 
 
905
      dpxf = 0.5d0*(dpdx(ii) + dpdx(jj))
 
906
      dpyf = 0.5d0*(dpdy(ii) + dpdy(jj))
 
907
      dpzf = 0.5d0*(dpdz(ii) + dpdz(jj))
 
908
 
 
909
      pip = pvar(ii)                                              &
 
910
           + ircflp*(dpxf*diipfx+dpyf*diipfy+dpzf*diipfz)
 
911
      pjp = pvar(jj)                                              &
 
912
           + ircflp*(dpxf*djjpfx+dpyf*djjpfy+dpzf*djjpfz)
 
913
 
 
914
      flui = 0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
915
      fluj = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
916
 
 
917
 
 
918
!         TEST DE PENTE
 
919
!        ---------------
 
920
 
 
921
      testi = dpdxa(ii)*surfac(1,ifac) +dpdya(ii)*surfac(2,ifac)  &
 
922
            + dpdza(ii)*surfac(3,ifac)
 
923
      testj = dpdxa(jj)*surfac(1,ifac) +dpdya(jj)*surfac(2,ifac)  &
 
924
            + dpdza(jj)*surfac(3,ifac)
 
925
      testij= dpdxa(ii)*dpdxa(jj)    +dpdya(ii)*dpdya(jj)         &
 
926
            + dpdza(ii)*dpdza(jj)
 
927
 
 
928
      if( flumas(ifac).gt.0.d0) then
 
929
        dcc = dpdx(ii)*surfac(1,ifac) +dpdy(ii)*surfac(2,ifac)    &
 
930
            + dpdz(ii)*surfac(3,ifac)
 
931
        ddi = testi
 
932
        ddj = ( pvar(jj)-pvar(ii) )/dist *surfan
 
933
      else
 
934
        dcc = dpdx(jj)*surfac(1,ifac) +dpdy(jj)*surfac(2,ifac)    &
 
935
            + dpdz(jj)*surfac(3,ifac)
 
936
        ddi = ( pvar(jj)-pvar(ii) )/dist *surfan
 
937
        ddj = testj
 
938
      endif
 
939
      tesqck = dcc**2 -(ddi-ddj)**2
 
940
 
 
941
 
 
942
!         UPWIND
 
943
!        --------
 
944
 
 
945
!MO          IF( (TESTI*TESTJ).LE.0.D0 .OR. TESTIJ.LE.0.D0 ) THEN
 
946
      if( tesqck.le.0.d0 .or. testij.le.0.d0 ) then
 
947
 
 
948
        pif = pvar(ii)
 
949
        pjf = pvar(jj)
 
950
        infac = infac+1
 
951
 
 
952
      else
 
953
 
 
954
 
 
955
!         CENTRE
 
956
!        --------
 
957
 
 
958
        if (ischcp.eq.1) then
 
959
 
 
960
          pif = pond*pip +(1.d0-pond)*pjp
 
961
          pjf = pif
 
962
 
 
963
 
 
964
!         SECOND ORDER
 
965
!        --------------
 
966
 
 
967
        elseif(ischcp.eq.0) then
 
968
 
 
969
          difx = cdgfac(1,ifac) - xyzcen(1,ii)
 
970
          dify = cdgfac(2,ifac) - xyzcen(2,ii)
 
971
          difz = cdgfac(3,ifac) - xyzcen(3,ii)
 
972
          djfx = cdgfac(1,ifac) - xyzcen(1,jj)
 
973
          djfy = cdgfac(2,ifac) - xyzcen(2,jj)
 
974
          djfz = cdgfac(3,ifac) - xyzcen(3,jj)
 
975
 
 
976
!     on laisse la reconstruction de PIF et PJF meme si IRCFLP=0
 
977
!     sinon cela revient a faire de l'upwind
 
978
          pif = pvar(ii)                                          &
 
979
              + difx*dpdx(ii)+dify*dpdy(ii)+difz*dpdz(ii)
 
980
          pjf = pvar(jj)                                          &
 
981
              + djfx*dpdx(jj)+djfy*dpdy(jj)+djfz*dpdz(jj)
 
982
 
 
983
        else
 
984
          write(nfecra,9000)ischcp
 
985
          iok = 1
 
986
        endif
 
987
 
 
988
      endif
 
989
 
 
990
 
 
991
!        BLENDING
 
992
!       ----------
 
993
 
 
994
      pif = blencp*pif+(1.d0-blencp)*pvar(ii)
 
995
      pjf = blencp*pjf+(1.d0-blencp)*pvar(jj)
 
996
 
 
997
 
 
998
!        FLUX
 
999
!       ------
 
1000
 
 
1001
      flux = iconvp*( flui*pif +fluj*pjf )                        &
 
1002
           + idiffp*viscf(ifac)*( pip -pjp )
 
1003
 
 
1004
 
 
1005
!        ASSEMBLAGE
 
1006
!       ------------
 
1007
 
 
1008
      smbrp(ii) = smbrp(ii) -flux
 
1009
      smbrp(jj) = smbrp(jj) +flux
 
1010
 
 
1011
    enddo
 
1012
!        Le call csexit ne doit pas etre dans la boucle, sinon
 
1013
!         le VPP la devectorise.
 
1014
    if(iok.ne.0) then
 
1015
      call csexit (1)
 
1016
    endif
 
1017
 
 
1018
  else
 
1019
 
 
1020
    iok = 0
 
1021
! VECTORISATION NON FORCEE
 
1022
    do ifac = 1, nfac
 
1023
 
 
1024
      ii = ifacel(1,ifac)
 
1025
      jj = ifacel(2,ifac)
 
1026
 
 
1027
      iij = idijpf-1+3*(ifac-1)
 
1028
 
 
1029
      dijpfx = ra(iij+1)
 
1030
      dijpfy = ra(iij+2)
 
1031
      dijpfz = ra(iij+3)
 
1032
 
 
1033
      pond   = ra(ipond-1+ifac)
 
1034
      dist   = ra(idist-1+ifac)
 
1035
      surfan = ra(isrfan-1+ifac)
 
1036
 
 
1037
! ON RECALCULE II' ET JJ'
 
1038
 
 
1039
      diipfx = cdgfac(1,ifac) - (xyzcen(1,ii)+                    &
 
1040
               (1.d0-pond) * dijpfx)
 
1041
      diipfy = cdgfac(2,ifac) - (xyzcen(2,ii)+                    &
 
1042
               (1.d0-pond) * dijpfy)
 
1043
      diipfz = cdgfac(3,ifac) - (xyzcen(3,ii)+                    &
 
1044
               (1.d0-pond) * dijpfz)
 
1045
      djjpfx = cdgfac(1,ifac) -  xyzcen(1,jj)+                    &
 
1046
                   pond  * dijpfx
 
1047
      djjpfy = cdgfac(2,ifac) -  xyzcen(2,jj)+                    &
 
1048
                   pond  * dijpfy
 
1049
      djjpfz = cdgfac(3,ifac) -  xyzcen(3,jj)+                    &
 
1050
                   pond  * dijpfz
 
1051
 
 
1052
      dpxf = 0.5d0*(dpdx(ii) + dpdx(jj))
 
1053
      dpyf = 0.5d0*(dpdy(ii) + dpdy(jj))
 
1054
      dpzf = 0.5d0*(dpdz(ii) + dpdz(jj))
 
1055
 
 
1056
      pip = pvar(ii)                                              &
 
1057
           + ircflp*(dpxf*diipfx+dpyf*diipfy+dpzf*diipfz)
 
1058
      pjp = pvar(jj)                                              &
 
1059
           + ircflp*(dpxf*djjpfx+dpyf*djjpfy+dpzf*djjpfz)
 
1060
 
 
1061
      flui = 0.5d0*( flumas(ifac) +abs(flumas(ifac)) )
 
1062
      fluj = 0.5d0*( flumas(ifac) -abs(flumas(ifac)) )
 
1063
 
 
1064
 
 
1065
!         TEST DE PENTE
 
1066
!        ---------------
 
1067
 
 
1068
      testi = dpdxa(ii)*surfac(1,ifac) +dpdya(ii)*surfac(2,ifac)  &
 
1069
            + dpdza(ii)*surfac(3,ifac)
 
1070
      testj = dpdxa(jj)*surfac(1,ifac) +dpdya(jj)*surfac(2,ifac)  &
 
1071
            + dpdza(jj)*surfac(3,ifac)
 
1072
      testij= dpdxa(ii)*dpdxa(jj)    +dpdya(ii)*dpdya(jj)         &
 
1073
            + dpdza(ii)*dpdza(jj)
 
1074
 
 
1075
      if( flumas(ifac).gt.0.d0) then
 
1076
        dcc = dpdx(ii)*surfac(1,ifac) +dpdy(ii)*surfac(2,ifac)    &
 
1077
            + dpdz(ii)*surfac(3,ifac)
 
1078
        ddi = testi
 
1079
        ddj = ( pvar(jj)-pvar(ii) )/dist *surfan
 
1080
      else
 
1081
        dcc = dpdx(jj)*surfac(1,ifac) +dpdy(jj)*surfac(2,ifac)    &
 
1082
            + dpdz(jj)*surfac(3,ifac)
 
1083
        ddi = ( pvar(jj)-pvar(ii) )/dist *surfan
 
1084
        ddj = testj
 
1085
      endif
 
1086
      tesqck = dcc**2 -(ddi-ddj)**2
 
1087
 
 
1088
 
 
1089
!         UPWIND
 
1090
!        --------
 
1091
 
 
1092
!MO          IF( (TESTI*TESTJ).LE.0.D0 .OR. TESTIJ.LE.0.D0 ) THEN
 
1093
      if( tesqck.le.0.d0 .or. testij.le.0.d0 ) then
 
1094
 
 
1095
        pif = pvar(ii)
 
1096
        pjf = pvar(jj)
 
1097
        infac = infac+1
 
1098
 
 
1099
      else
 
1100
 
 
1101
 
 
1102
!         CENTRE
 
1103
!        --------
 
1104
 
 
1105
        if (ischcp.eq.1) then
 
1106
 
 
1107
          pif = pond*pip +(1.d0-pond)*pjp
 
1108
          pjf = pif
 
1109
 
 
1110
 
 
1111
!         SECOND ORDER
 
1112
!        --------------
 
1113
 
 
1114
        elseif(ischcp.eq.0) then
 
1115
 
 
1116
          difx = cdgfac(1,ifac) - xyzcen(1,ii)
 
1117
          dify = cdgfac(2,ifac) - xyzcen(2,ii)
 
1118
          difz = cdgfac(3,ifac) - xyzcen(3,ii)
 
1119
          djfx = cdgfac(1,ifac) - xyzcen(1,jj)
 
1120
          djfy = cdgfac(2,ifac) - xyzcen(2,jj)
 
1121
          djfz = cdgfac(3,ifac) - xyzcen(3,jj)
 
1122
 
 
1123
!     on laisse la reconstruction de PIF et PJF meme si IRCFLP=0
 
1124
!     sinon cela revient a faire de l'upwind
 
1125
          pif = pvar(ii)                                          &
 
1126
              + difx*dpdx(ii)+dify*dpdy(ii)+difz*dpdz(ii)
 
1127
          pjf = pvar(jj)                                          &
 
1128
              + djfx*dpdx(jj)+djfy*dpdy(jj)+djfz*dpdz(jj)
 
1129
 
 
1130
        else
 
1131
          write(nfecra,9000)ischcp
 
1132
          iok = 1
 
1133
        endif
 
1134
 
 
1135
      endif
 
1136
 
 
1137
 
 
1138
!        BLENDING
 
1139
!       ----------
 
1140
 
 
1141
      pif = blencp*pif+(1.d0-blencp)*pvar(ii)
 
1142
      pjf = blencp*pjf+(1.d0-blencp)*pvar(jj)
 
1143
 
 
1144
 
 
1145
!        FLUX
 
1146
!       ------
 
1147
 
 
1148
      flux = iconvp*( flui*pif +fluj*pjf )                        &
 
1149
           + idiffp*viscf(ifac)*( pip -pjp )
 
1150
 
 
1151
 
 
1152
!        ASSEMBLAGE
 
1153
!       ------------
 
1154
 
 
1155
      smbrp(ii) = smbrp(ii) -flux
 
1156
      smbrp(jj) = smbrp(jj) +flux
 
1157
 
 
1158
    enddo
 
1159
!        Le call csexit ne doit pas etre dans la boucle, sinon
 
1160
!         le VPP la devectorise (pour la boucle non vectorisee forcee,
 
1161
!         ce n'est pas grave, c'est juste pour recopier exactement
 
1162
!         la precedente)
 
1163
    if(iok.ne.0) then
 
1164
      call csexit (1)
 
1165
    endif
 
1166
 
 
1167
  endif
 
1168
 
 
1169
endif
 
1170
 
 
1171
 
 
1172
 
 
1173
if(iwarnp.ge.2) then
 
1174
  write(nfecra,1100)cnom,infac,nfac
 
1175
endif
 
1176
 
 
1177
 
 
1178
! ======================================================================
 
1179
! ---> ASSEMBLAGE A PARTIR DES FACETTES DE BORD
 
1180
! ======================================================================
 
1181
 
 
1182
!     Lorsque IIFBRU.GT.0, on ne prend pas en compte le flux convectif
 
1183
!       sur les faces pour lesquelles on dispose par ailleurs d'un
 
1184
!       flux de C.L.
 
1185
 
 
1186
!     Lorsque IIFBRU.LE.0, on adopte le sch�ma standard (i.e. on prend
 
1187
!       en compte le flux convectitf).
 
1188
 
 
1189
!     Dans les deux cas, si on prend en compte le flux convectif, il
 
1190
!       pourrait etre utile de vraiment utiliser la condition imposee
 
1191
!       a la limite, i.e. de ne pas utiliser un sch�ma upwind (sinon,
 
1192
!       on ne verra pas certaines C.L.). Actuellement, avec les conditions
 
1193
!       aux limites propos�es, il n'y a que 3 cas pour lesquels on
 
1194
!       prend en compte le flux convectif ici : paroi, symetrie et
 
1195
!       sortie supersonique. Dans les deux premiers, le flux est nul.
 
1196
!       Pour le dernier, un traitement upwind correspond a utiliser
 
1197
!       effectivement la valeur de bord. On peut donc conserver
 
1198
!       le code tel quel.
 
1199
 
 
1200
!     On n'impose pas le flux convectif sur les faces pour lesquelles
 
1201
!       il sera impos� par les C.L.
 
1202
if(iifbru.gt.0) then
 
1203
 
 
1204
  if (ivectb.eq.1) then
 
1205
 
 
1206
!CDIR NODEP
 
1207
    do ifac = 1, nfabor
 
1208
 
 
1209
      ii = ifabor(ifac)
 
1210
 
 
1211
      iii = idiipb-1+3*(ifac-1)
 
1212
      diipbx = ra(iii+1)
 
1213
      diipby = ra(iii+2)
 
1214
      diipbz = ra(iii+3)
 
1215
 
 
1216
      flui = 0.5d0*( flumab(ifac) +abs(flumab(ifac)) )
 
1217
      fluj = 0.5d0*( flumab(ifac) -abs(flumab(ifac)) )
 
1218
 
 
1219
      pip = pvar(ii)                                              &
 
1220
       +ircflp*(dpdx(ii)*diipbx+dpdy(ii)*diipby+dpdz(ii)*diipbz)
 
1221
 
 
1222
      pfac  = inc*coefap(ifac) +coefbp(ifac)*pip
 
1223
      pfacd = inc*cofafp(ifac) +cofbfp(ifac)*pip
 
1224
 
 
1225
!            FLUX = ICONVP*( FLUI*PVAR(II) +FLUJ*PFAC )
 
1226
!     &           + IDIFFP*VISCB(IFAC)*( PIP -PFACD )
 
1227
      flux = iconvp*( flui*pvar(ii) +fluj*pfac )                  &
 
1228
                   *dble(1-ifrusb(ifac))                          &
 
1229
           + idiffp*viscb(ifac)*( pip -pfacd )
 
1230
      smbrp(ii) = smbrp(ii) -flux
 
1231
 
 
1232
    enddo
 
1233
 
 
1234
  else
 
1235
 
 
1236
    do ifac = 1, nfabor
 
1237
 
 
1238
      ii = ifabor(ifac)
 
1239
 
 
1240
      iii = idiipb-1+3*(ifac-1)
 
1241
      diipbx = ra(iii+1)
 
1242
      diipby = ra(iii+2)
 
1243
      diipbz = ra(iii+3)
 
1244
 
 
1245
      flui = 0.5d0*( flumab(ifac) +abs(flumab(ifac)) )
 
1246
      fluj = 0.5d0*( flumab(ifac) -abs(flumab(ifac)) )
 
1247
 
 
1248
      pip = pvar(ii)                                              &
 
1249
       +ircflp*(dpdx(ii)*diipbx+dpdy(ii)*diipby+dpdz(ii)*diipbz)
 
1250
 
 
1251
      pfac  = inc*coefap(ifac) +coefbp(ifac)*pip
 
1252
      pfacd = inc*cofafp(ifac) +cofbfp(ifac)*pip
 
1253
 
 
1254
!            FLUX = ICONVP*( FLUI*PVAR(II) +FLUJ*PFAC )
 
1255
!     &           + IDIFFP*VISCB(IFAC)*( PIP -PFACD )
 
1256
      flux = iconvp*( flui*pvar(ii) +fluj*pfac )                  &
 
1257
                   *dble(1-ifrusb(ifac))                          &
 
1258
           + idiffp*viscb(ifac)*( pip -pfacd )
 
1259
      smbrp(ii) = smbrp(ii) -flux
 
1260
 
 
1261
    enddo
 
1262
 
 
1263
  endif
 
1264
 
 
1265
!     On ne dispose pas de flux issu des C.L. : traitement std
 
1266
else
 
1267
 
 
1268
  if (ivectb.eq.1) then
 
1269
 
 
1270
!CDIR NODEP
 
1271
    do ifac = 1, nfabor
 
1272
 
 
1273
      ii = ifabor(ifac)
 
1274
 
 
1275
      iii = idiipb-1+3*(ifac-1)
 
1276
      diipbx = ra(iii+1)
 
1277
      diipby = ra(iii+2)
 
1278
      diipbz = ra(iii+3)
 
1279
 
 
1280
      flui = 0.5d0*( flumab(ifac) +abs(flumab(ifac)) )
 
1281
      fluj = 0.5d0*( flumab(ifac) -abs(flumab(ifac)) )
 
1282
 
 
1283
      pip = pvar(ii)                                              &
 
1284
       +ircflp*(dpdx(ii)*diipbx+dpdy(ii)*diipby+dpdz(ii)*diipbz)
 
1285
 
 
1286
      pfac  = inc*coefap(ifac) +coefbp(ifac)*pip
 
1287
      pfacd = inc*cofafp(ifac) +cofbfp(ifac)*pip
 
1288
 
 
1289
      flux = iconvp*( flui*pvar(ii) +fluj*pfac )                  &
 
1290
           + idiffp*viscb(ifac)*( pip -pfacd )
 
1291
      smbrp(ii) = smbrp(ii) -flux
 
1292
 
 
1293
    enddo
 
1294
 
 
1295
  else
 
1296
 
 
1297
    do ifac = 1, nfabor
 
1298
 
 
1299
      ii = ifabor(ifac)
 
1300
 
 
1301
      iii = idiipb-1+3*(ifac-1)
 
1302
      diipbx = ra(iii+1)
 
1303
      diipby = ra(iii+2)
 
1304
      diipbz = ra(iii+3)
 
1305
 
 
1306
      flui = 0.5d0*( flumab(ifac) +abs(flumab(ifac)) )
 
1307
      fluj = 0.5d0*( flumab(ifac) -abs(flumab(ifac)) )
 
1308
 
 
1309
      pip = pvar(ii)                                              &
 
1310
       +ircflp*(dpdx(ii)*diipbx+dpdy(ii)*diipby+dpdz(ii)*diipbz)
 
1311
 
 
1312
      pfac  = inc*coefap(ifac) +coefbp(ifac)*pip
 
1313
      pfacd = inc*cofafp(ifac) +cofbfp(ifac)*pip
 
1314
 
 
1315
      flux = iconvp*( flui*pvar(ii) +fluj*pfac )                  &
 
1316
           + idiffp*viscb(ifac)*( pip -pfacd )
 
1317
      smbrp(ii) = smbrp(ii) -flux
 
1318
 
 
1319
    enddo
 
1320
 
 
1321
  endif
 
1322
 
 
1323
endif
 
1324
 
 
1325
!--------
 
1326
! FORMATS
 
1327
!--------
 
1328
 
 
1329
 1000 format(1X,A8,' : CONVECTION EN ',A11,                             &
 
1330
                               ' BLENDING A ',F4.0,' % D''UPWIND')
 
1331
 1100 format(1X,A8,' : ',I10,' FACES UPWIND SUR ',                      &
 
1332
                               I10,' FACES INTERNES ')
 
1333
 9000 format(                                                           &
 
1334
'@                                                            ',/,&
 
1335
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
 
1336
'@                                                            ',/,&
 
1337
'@ @@ ATTENTION : ARRET DANS cfbsc2                           ',/,&
 
1338
'@    =========                                               ',/,&
 
1339
'@     APPEL DE cfbsc2 POUR ',A8 ,' AVEC ISCHCP = ',I10        ,/,&
 
1340
'@                                                            ',/,&
 
1341
'@  Le calcul ne peut pas etre execute.                       ',/,&
 
1342
'@                                                            ',/,&
 
1343
'@  Contacter l''assistance.                                  ',/,&
 
1344
'@                                                            ',/,&
 
1345
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
 
1346
'@                                                            ',/)
 
1347
 
 
1348
!----
 
1349
! FIN
 
1350
!----
 
1351
 
 
1352
return
 
1353
 
 
1354
end