← Back to branch summary

~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to users/cfbl/uscfth.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-24 00:00:08 UTC
  • mfrom: (6.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20111124000008-2vo99e38267942q5
Tags: 2.1.0-3
Install a missing file

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/cfbl/uscfth.f90
2
2
 
3
3
!VERS
4
4
 
5
 
 
6
 
!     This file is part of the Code_Saturne Kernel, element of the
7
 
!     Code_Saturne CFD tool.
8
 
 
9
 
!     Copyright (C) 1998-2009 EDF S.A., France
10
 
 
11
 
!     contact: saturne-support@edf.fr
12
 
 
13
 
!     The Code_Saturne Kernel is free software; you can redistribute it
14
 
!     and/or modify it under the terms of the GNU General Public License
15
 
!     as published by the Free Software Foundation; either version 2 of
16
 
!     the License, or (at your option) any later version.
17
 
 
18
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
19
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
20
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21
 
!     GNU General Public License for more details.
22
 
 
23
 
!     You should have received a copy of the GNU General Public License
24
 
!     along with the Code_Saturne Kernel; if not, write to the
25
 
!     Free Software Foundation, Inc.,
26
 
!     51 Franklin St, Fifth Floor,
27
 
!     Boston, MA  02110-1301  USA
 
5
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
6
!
 
7
! Copyright (C) 1998-2011 EDF S.A.
 
8
!
 
9
! This program is free software; you can redistribute it and/or modify it under
 
10
! the terms of the GNU General Public License as published by the Free Software
 
11
! Foundation; either version 2 of the License, or (at your option) any later
 
12
! version.
 
13
!
 
14
! This program is distributed in the hope that it will be useful, but WITHOUT
 
15
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 
16
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 
17
! details.
 
18
!
 
19
! You should have received a copy of the GNU General Public License along with
 
20
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
 
21
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
28
22
 
29
23
!-------------------------------------------------------------------------------
30
24
 
31
25
subroutine uscfth &
32
26
!================
33
27
 
34
 
 ( idbia0 , idbra0 ,                                              &
35
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
36
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
37
 
   nvar   , nscal  , nphas  ,                                     &
38
 
   iccfth , imodif , iphas  ,                                     &
39
 
   nideve , nrdeve , nituse , nrtuse ,                            &
40
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
41
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
42
 
   idevel , ituser , ia     ,                                     &
43
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
28
 ( nvar   , nscal  ,                                              &
 
29
   iccfth , imodif ,                                              &
44
30
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
45
31
   coefa  , coefb  ,                                              &
46
 
   sorti1 , sorti2 , gamagr , xmasmr ,                            &
47
 
   rdevel , rtuser , ra     )
 
32
   sorti1 , sorti2 , gamagr , xmasmr )
48
33
 
49
34
!===============================================================================
50
35
! Purpose:
201
186
!__________________.____._____.________________________________________________.
202
187
!    nom           !type!mode !                   role                         !
203
188
!__________________!____!_____!________________________________________________!
204
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
205
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
206
 
! ndim             ! i  ! <-- ! spatial dimension                              !
207
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
208
 
! ncel             ! i  ! <-- ! number of cells                                !
209
 
! nfac             ! i  ! <-- ! number of interior faces                       !
210
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
211
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
212
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
213
 
! nnod             ! i  ! <-- ! number of vertices                             !
214
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
215
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
216
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
217
189
! nvar             ! i  ! <-- ! total number of variables                      !
218
190
! nscal            ! i  ! <-- ! total number of scalars                        !
219
 
! nphas            ! i  ! <-- ! number of phases                               !
220
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
221
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
222
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
223
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
224
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
225
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
226
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
227
 
!  (nfml, nprfml)  !    !     !                                                !
228
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
229
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
230
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
231
 
! nodfac(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
232
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary developpement !
233
 
! ituser(nituse    ! ia ! <-> ! user-reserved integer work array               !
234
 
! ia(*)            ! ia ! --- ! main integer work array                        !
235
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
236
 
!  (ndim, ncelet)  !    !     !                                                !
237
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
238
 
!  (ndim, nfac)    !    !     !                                                !
239
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
240
 
!  (ndim, nfavor)  !    !     !                                                !
241
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
242
 
!  (ndim, nfac)    !    !     !                                                !
243
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
244
 
!  (ndim, nfabor)  !    !     !                                                !
245
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
246
 
!  (ndim, nnod)    !    !     !                                                !
247
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
248
191
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
249
192
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
250
193
!  (ncelet, *)     !    !     !  (at current and preceding time steps)         !
259
202
!                  !    !     !   (first value only used for perfect gas)      !
260
203
! xmasmr(*)        ! ra ! --> ! molar mass of the components of the gas        !
261
204
!                  !    !     !   (unused if iccfth.lt.0)                      !
262
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary developpement    !
263
 
! rtuser(nituse    ! ra ! <-> ! user-reserved real work array                  !
264
 
! ra(*)            ! ra ! --- ! main real work array                           !
265
205
!__________________!____!_____!________________________________________________!
266
206
 
267
207
!     Type: i (integer), r (real), s (string), a (array), l (logical),
269
209
!     mode: <-- input, --> output, <-> modifies data, --- work array
270
210
!===============================================================================
271
211
 
 
212
!===============================================================================
 
213
! Module files
 
214
!===============================================================================
 
215
 
 
216
use paramx
 
217
use numvar
 
218
use optcal
 
219
use cstphy
 
220
use cstnum
 
221
use parall
 
222
use pointe
 
223
use entsor
 
224
use ppppar
 
225
use ppthch
 
226
use ppincl
 
227
use mesh
 
228
 
 
229
!===============================================================================
 
230
 
272
231
implicit none
273
232
 
274
 
!===============================================================================
275
 
! Common blocks
276
 
!===============================================================================
277
 
 
278
 
include "paramx.h"
279
 
include "numvar.h"
280
 
include "optcal.h"
281
 
include "cstphy.h"
282
 
include "cstnum.h"
283
 
include "pointe.h"
284
 
include "entsor.h"
285
 
include "parall.h"
286
 
include "ppppar.h"
287
 
include "ppthch.h"
288
 
include "ppincl.h"
289
 
 
290
 
!===============================================================================
291
 
 
292
233
! Arguments
293
234
 
294
 
integer          idbia0 , idbra0
295
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
296
 
integer          nfml   , nprfml
297
 
integer          nnod   , lndfac , lndfbr , ncelbr
298
 
integer          nvar   , nscal  , nphas
299
 
integer          iccfth   , imodif , iphas
300
 
integer          nideve , nrdeve , nituse , nrtuse
301
 
 
302
 
integer          ifacel(2,nfac) , ifabor(nfabor)
303
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
304
 
integer          iprfml(nfml,nprfml)
305
 
integer          ipnfac(nfac+1), nodfac(lndfac)
306
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
307
 
integer          idevel(nideve), ituser(nituse), ia(*)
308
 
 
309
 
double precision xyzcen(ndim,ncelet)
310
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
311
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
312
 
double precision xyznod(ndim,nnod), volume(ncelet)
 
235
integer          nvar   , nscal
 
236
integer          iccfth   , imodif
 
237
 
313
238
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
314
239
double precision propce(ncelet,*),propfa(nfac,*),propfb(nfabor,*)
315
240
double precision coefa(nfabor,*), coefb(nfabor,*)
316
241
 
317
242
double precision sorti1(*), sorti2(*), gamagr(*), xmasmr(*)
318
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
319
243
 
320
244
! Local variables
321
245
 
322
 
integer          idebia, idebra
323
 
integer          iiph   , ifac0
 
246
integer          ifac0
324
247
integer          ierr
325
248
integer          iel    , ifac   , ivar
326
 
integer          ipriph , irhiph , itkiph , ieniph
327
 
integer          iuiph  , iviph  , iwiph
 
249
integer          irh    , itk    , ien
328
250
integer          iclp   , iclr   , iclt   , icle
329
251
integer          iclu   , iclv   , iclw
330
252
integer          iutile
357
279
!    No user input required.
358
280
!===============================================================================
359
281
 
360
 
! Memory pointers
361
 
idebia = idbia0
362
 
idebra = idbra0
363
 
 
364
282
! Error indicator (stop if non zero)
365
283
ierr   = 0
366
284
 
367
285
! Rank of the variables in their associated arrays
368
286
if(iccfth.ge.0.or.iccfth.le.-2) then
369
 
  ipriph = ipr(iphas)
370
 
  irhiph = isca(irho  (iphas))
371
 
  itkiph = isca(itempk(iphas))
372
 
  ieniph = isca(ienerg(iphas))
373
 
  iuiph = iu(iphas)
374
 
  iviph = iv(iphas)
375
 
  iwiph = iw(iphas)
376
 
  iclp = iclrtp(ipriph,icoef)
377
 
  iclr = iclrtp(irhiph,icoef)
378
 
  iclt = iclrtp(itkiph,icoef)
379
 
  icle = iclrtp(ieniph,icoef)
380
 
  iclu = iclrtp(iuiph,icoef)
381
 
  iclv = iclrtp(iviph,icoef)
382
 
  iclw = iclrtp(iwiph,icoef)
 
287
  irh = isca(irho  )
 
288
  itk = isca(itempk)
 
289
  ien = isca(ienerg)
 
290
  iclp = iclrtp(ipr,icoef)
 
291
  iclr = iclrtp(irh,icoef)
 
292
  iclt = iclrtp(itk,icoef)
 
293
  icle = iclrtp(ien,icoef)
 
294
  iclu = iclrtp(iu,icoef)
 
295
  iclv = iclrtp(iv,icoef)
 
296
  iclw = iclrtp(iw,icoef)
383
297
endif
384
298
 
385
299
! For calculation of values at the cell centers,
397
311
! --> ieos = 2: Perfect gas with variable Gamma
398
312
! --> ieos = 3: Van Der Waals
399
313
 
400
 
do iiph = 1, nphas
401
 
  ieos(iiph) = 1
402
 
enddo
 
314
ieos = 1
403
315
 
404
316
 
405
317
! Warning: once the thermodynamic law has been chosen,
410
322
! 2. Perfect gas
411
323
!===============================================================================
412
324
 
413
 
if(ieos(iphas).eq.1) then
 
325
if(ieos.eq.1) then
414
326
 
415
327
!===============================================================================
416
328
! 2.1. Parameters to be completed by the user
421
333
! --- Molar mass of the gas (kg/mol)
422
334
 
423
335
  if(iccfth.ge.0) then
424
 
 
425
 
    if(iphas.eq.1) then
426
 
 
427
 
      xmasml = 28.8d-3
428
 
 
429
 
    endif
430
 
 
 
336
    xmasml = 28.8d-3
431
337
  endif
432
338
 
433
339
!===============================================================================
445
351
    !   constant is not saved. A ''save'' instruction and a test would
446
352
    !   be sufficient to avoid computing gamagp at each call if necessary.
447
353
 
448
 
    gamagp = 1.d0 + rr/(xmasml*cp0(iphas)-rr)
 
354
    gamagp = 1.d0 + rr/(xmasml*cp0-rr)
449
355
 
450
356
    if(gamagp.lt.1.d0) then
451
357
      write(nfecra,1010) gamagp
469
375
    !   the user subroutine ''usini1''. The value for the isochoric
470
376
    !   specific heat Cv0 is calculated in a subsequent section (from Cp0)
471
377
 
472
 
    icp(iphas) = 0
473
 
    icv(iphas) = 0
 
378
    icp = 0
 
379
    icv = 0
474
380
 
475
381
 
476
382
! --- Default initializations (before uscfxi)
480
386
 
481
387
  elseif(iccfth.eq.0) then
482
388
 
483
 
    cv0(iphas) = cp0(iphas) - rr/xmasml
 
389
    cv0 = cp0 - rr/xmasml
484
390
 
485
391
    if ( isuite .eq. 0 ) then
486
392
      do iel = 1, ncel
487
 
        rtp(iel,irhiph) = p0(iphas)*xmasml/(rr*t0(iphas))
488
 
        rtp(iel,ieniph) = cv0(iphas)*t0(iphas)
 
393
        rtp(iel,irh) = p0*xmasml/(rr*t0)
 
394
        rtp(iel,ien) = cv0*t0
489
395
      enddo
490
396
    endif
491
397
 
502
408
 
503
409
    ierr = 0
504
410
    do iel = 1, ncel
505
 
      if(rtp(iel,irhiph).le.0.d0) then
506
 
        rtp(iel,irhiph) = epzero
 
411
      if(rtp(iel,irh).le.0.d0) then
 
412
        rtp(iel,irh) = epzero
507
413
        ierr = ierr + 1
508
414
      endif
509
415
    enddo
526
432
 
527
433
    ierr = 0
528
434
    do iel = 1, ncel
529
 
      enint = rtp(iel,ieniph)                                     &
530
 
               - 0.5d0*( rtp(iel,iuiph)**2                        &
531
 
                       + rtp(iel,iviph)**2                        &
532
 
                       + rtp(iel,iwiph)**2 )
 
435
      enint = rtp(iel,ien)                                     &
 
436
               - 0.5d0*( rtp(iel,iu)**2                        &
 
437
                       + rtp(iel,iv)**2                        &
 
438
                       + rtp(iel,iw)**2 )
533
439
      if(enint.le.0.d0) then
534
 
        rtp(iel,ieniph) = epzero                                  &
535
 
               + 0.5d0*( rtp(iel,iuiph)**2                        &
536
 
                       + rtp(iel,iviph)**2                        &
537
 
                       + rtp(iel,iwiph)**2 )
 
440
        rtp(iel,ien) = epzero                                  &
 
441
               + 0.5d0*( rtp(iel,iu)**2                        &
 
442
                       + rtp(iel,iv)**2                        &
 
443
                       + rtp(iel,iw)**2 )
538
444
        ierr = ierr + 1
539
445
      endif
540
446
    enddo
554
460
    ! Verification of the values of the density
555
461
    ierr = 0
556
462
    do iel = 1, ncel
557
 
      if(rtp(iel,irhiph).le.0.d0) then
558
 
        write(nfecra,3010)rtp(iel,irhiph),iel
 
463
      if(rtp(iel,irh).le.0.d0) then
 
464
        write(nfecra,3010)rtp(iel,irh),iel
559
465
      endif
560
466
    enddo
561
467
    ! Stop if a negative value is detected (since the density has been
567
473
 
568
474
    do iel = 1, ncel
569
475
      ! Temperature
570
 
      sorti1(iel) = xmasml*rtp(iel,ipriph)/(rr*rtp(iel,irhiph))
 
476
      sorti1(iel) = xmasml*rtp(iel,ipr)/(rr*rtp(iel,irh))
571
477
      ! Total energy
572
 
      sorti2(iel) = cv0(iphas)*sorti1(iel)                        &
573
 
           + 0.5d0*( rtp(iel,iuiph)**2 + rtp(iel,iviph)**2        &
574
 
                                       + rtp(iel,iwiph)**2 )
 
478
      sorti2(iel) = cv0*sorti1(iel)                        &
 
479
           + 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 )
575
480
    enddo
576
481
 
577
482
    ! Transfer to the array rtp
578
483
    if(imodif.gt.0) then
579
484
      do iel = 1, ncel
580
 
        rtp(iel,itkiph) = sorti1(iel)
581
 
        rtp(iel,ieniph) = sorti2(iel)
 
485
        rtp(iel,itk) = sorti1(iel)
 
486
        rtp(iel,ien) = sorti2(iel)
582
487
      enddo
583
488
    endif
584
489
 
590
495
    ! Verification of the values of the temperature
591
496
    ierr = 0
592
497
    do iel = 1, ncel
593
 
      if(rtp(iel,itkiph).le.0.d0) then
594
 
        write(nfecra,2010)rtp(iel,itkiph),iel
 
498
      if(rtp(iel,itk).le.0.d0) then
 
499
        write(nfecra,2010)rtp(iel,itk),iel
595
500
      endif
596
501
    enddo
597
502
    ! Stop if a negative value is detected (since the temperature has been
603
508
 
604
509
    do iel = 1, ncel
605
510
      ! Density
606
 
      sorti1(iel) = xmasml*rtp(iel,ipriph)/(rr*rtp(iel,itkiph))
 
511
      sorti1(iel) = xmasml*rtp(iel,ipr)/(rr*rtp(iel,itk))
607
512
      ! Total energy
608
 
      sorti2(iel) = cv0(iphas)*rtp(iel,itkiph)                    &
609
 
           + 0.5d0*( rtp(iel,iuiph)**2 + rtp(iel,iviph)**2        &
610
 
                                       + rtp(iel,iwiph)**2 )
 
513
      sorti2(iel) = cv0*rtp(iel,itk)                    &
 
514
           + 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 )
611
515
    enddo
612
516
 
613
517
    ! Transfer to the array rtp
614
518
    if(imodif.gt.0) then
615
519
      do iel = 1, ncel
616
 
        rtp(iel,irhiph) = sorti1(iel)
617
 
        rtp(iel,ieniph) = sorti2(iel)
 
520
        rtp(iel,irh) = sorti1(iel)
 
521
        rtp(iel,ien) = sorti2(iel)
618
522
      enddo
619
523
    endif
620
524
 
626
530
    do iel = 1, ncel
627
531
      ! Internal energy (to avoid the need to divide by the temperature
628
532
      ! to compute density)
629
 
      enint = rtp(iel,ieniph)                                     &
630
 
               - 0.5d0*( rtp(iel,iuiph)**2                        &
631
 
                       + rtp(iel,iviph)**2                        &
632
 
                       + rtp(iel,iwiph)**2 )
 
533
      enint = rtp(iel,ien)                                     &
 
534
               - 0.5d0*( rtp(iel,iu)**2                        &
 
535
                       + rtp(iel,iv)**2                        &
 
536
                       + rtp(iel,iw)**2 )
633
537
      ! Density
634
 
      sorti1(iel) = rtp(iel,ipriph) / ( (gamagp-1.d0) * enint )
 
538
      sorti1(iel) = rtp(iel,ipr) / ( (gamagp-1.d0) * enint )
635
539
      ! Temperature
636
540
      sorti2(iel) = xmasml * (gamagp-1.d0) * enint / rr
637
541
    enddo
639
543
    ! Transfer to the array rtp
640
544
    if(imodif.gt.0) then
641
545
      do iel = 1, ncel
642
 
        rtp(iel,irhiph) = sorti1(iel)
643
 
        rtp(iel,itkiph) = sorti2(iel)
 
546
        rtp(iel,irh) = sorti1(iel)
 
547
        rtp(iel,itk) = sorti2(iel)
644
548
      enddo
645
549
    endif
646
550
 
651
555
 
652
556
    do iel = 1, ncel
653
557
      ! Pressure
654
 
      sorti1(iel) = rtp(iel,irhiph)*rtp(iel,itkiph)*rr/xmasml
 
558
      sorti1(iel) = rtp(iel,irh)*rtp(iel,itk)*rr/xmasml
655
559
      ! Total energy
656
 
      sorti2(iel) = cv0(iphas)*rtp(iel,itkiph)                    &
657
 
           + 0.5d0*( rtp(iel,iuiph)**2 + rtp(iel,iviph)**2        &
658
 
                                       + rtp(iel,iwiph)**2 )
 
560
      sorti2(iel) = cv0*rtp(iel,itk)                    &
 
561
           + 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 )
659
562
    enddo
660
563
 
661
564
    ! Transfer to the array rtp
662
565
    if(imodif.gt.0) then
663
566
      do iel = 1, ncel
664
 
        rtp(iel,ipriph) = sorti1(iel)
665
 
        rtp(iel,ieniph) = sorti2(iel)
 
567
        rtp(iel,ipr) = sorti1(iel)
 
568
        rtp(iel,ien) = sorti2(iel)
666
569
      enddo
667
570
    endif
668
571
 
674
577
    do iel = 1, ncel
675
578
      ! Internal energy (to avoid the need to divide by the temperature
676
579
      ! to compute density)
677
 
      enint = rtp(iel,ieniph)                                     &
678
 
               - 0.5d0*( rtp(iel,iuiph)**2                        &
679
 
                       + rtp(iel,iviph)**2                        &
680
 
                       + rtp(iel,iwiph)**2 )
 
580
      enint = rtp(iel,ien)                                     &
 
581
               - 0.5d0*( rtp(iel,iu)**2                        &
 
582
                       + rtp(iel,iv)**2                        &
 
583
                       + rtp(iel,iw)**2 )
681
584
      ! Pressure
682
 
      sorti1(iel) = (gamagp-1.d0) * rtp(iel,irhiph) * enint
 
585
      sorti1(iel) = (gamagp-1.d0) * rtp(iel,irh) * enint
683
586
      ! Temperature
684
587
      sorti2(iel) = xmasml * (gamagp-1.d0) * enint / rr
685
588
    enddo
687
590
    ! Transfer to the array rtp
688
591
    if(imodif.gt.0) then
689
592
      do iel = 1, ncel
690
 
        rtp(iel,ipriph) = sorti1(iel)
691
 
        rtp(iel,itkiph) = sorti2(iel)
 
593
        rtp(iel,ipr) = sorti1(iel)
 
594
        rtp(iel,itk) = sorti2(iel)
692
595
      enddo
693
596
    endif
694
597
 
707
610
    if(iutile.eq.1) then
708
611
      ierr = 0
709
612
      do iel = 1, ncel
710
 
        if(rtp(iel,irhiph).le.0.d0) then
711
 
          write(nfecra,4010)rtp(iel,irhiph),iel
 
613
        if(rtp(iel,irh).le.0.d0) then
 
614
          write(nfecra,4010)rtp(iel,irh),iel
712
615
        endif
713
616
      enddo
714
617
      if(ierr.eq.1) then
717
620
    endif
718
621
 
719
622
    do iel = 1, ncel
720
 
      sorti1(iel) = gamagp * rtp(iel,ipriph) / rtp(iel,irhiph)
 
623
      sorti1(iel) = gamagp * rtp(iel,ipr) / rtp(iel,irh)
721
624
    enddo
722
625
 
723
626
 
734
637
    if(iutile.eq.1) then
735
638
      ierr = 0
736
639
      do iel = 1, ncel
737
 
        if(rtp(iel,irhiph).lt.0.d0) then
738
 
          write(nfecra,4020)rtp(iel,irhiph),iel
 
640
        if(rtp(iel,irh).lt.0.d0) then
 
641
          write(nfecra,4020)rtp(iel,irh),iel
739
642
        endif
740
643
      enddo
741
644
      if(ierr.eq.1) then
744
647
    endif
745
648
 
746
649
    do iel = 1, ncel
747
 
      sorti1(iel) = rtp(iel,irhiph)**gamagp
 
650
      sorti1(iel) = rtp(iel,irh)**gamagp
748
651
    enddo
749
652
 
750
653
 
765
668
    !     density is <= 0, the calculation will simply fail)
766
669
    ierr = 0
767
670
    do iel = 1, ncel
768
 
      if(rtp(iel,irhiph).le.0.d0) then
769
 
        write(nfecra,4030)rtp(iel,irhiph),iel
 
671
      if(rtp(iel,irh).le.0.d0) then
 
672
        write(nfecra,4030)rtp(iel,irh),iel
770
673
      endif
771
674
    enddo
772
675
    if(ierr.eq.1) then
774
677
    endif
775
678
 
776
679
    do iel = 1, ncel
777
 
      sorti1(iel) = rtp(iel,ipriph) / (rtp(iel,irhiph)**gamagp)
 
680
      sorti1(iel) = rtp(iel,ipr) / (rtp(iel,irh)**gamagp)
778
681
    enddo
779
682
 
780
683
 
807
710
    ! Calculation of the Mach number at the boundary face, using the
808
711
    !   cell center velocity projected on the vector normal to the boundary
809
712
    xmach =                                                       &
810
 
         ( rtp(iel,iuiph)*surfbo(1,ifac)                          &
811
 
         + rtp(iel,iviph)*surfbo(2,ifac)                          &
812
 
         + rtp(iel,iwiph)*surfbo(3,ifac) ) / ra(isrfbn+ifac-1)    &
813
 
         / sqrt( gamagp * rtp(iel,ipriph) / rtp(iel,irhiph) )
 
713
         ( rtp(iel,iu)*surfbo(1,ifac)                          &
 
714
         + rtp(iel,iv)*surfbo(2,ifac)                          &
 
715
         + rtp(iel,iw)*surfbo(3,ifac) ) / surfbn(ifac)         &
 
716
         / sqrt( gamagp * rtp(iel,ipr) / rtp(iel,irh) )
814
717
 
815
718
    ! Pressure
816
719
 
886
789
    ! Calculation of the Mach number at the boundary face, using the
887
790
    !   cell center velocity projected on the vector normal to the boundary
888
791
    xmachi =                                                      &
889
 
         ( rtp(iel,iuiph)*surfbo(1,ifac)                          &
890
 
         + rtp(iel,iviph)*surfbo(2,ifac)                          &
891
 
         + rtp(iel,iwiph)*surfbo(3,ifac) ) / ra(isrfbn+ifac-1)    &
892
 
         / sqrt( gamagp * rtp(iel,ipriph) / rtp(iel,irhiph) )
 
792
         ( rtp(iel,iu)*surfbo(1,ifac)                          &
 
793
         + rtp(iel,iv)*surfbo(2,ifac)                          &
 
794
         + rtp(iel,iw)*surfbo(3,ifac) ) / surfbn(ifac)         &
 
795
         / sqrt( gamagp * rtp(iel,ipr) / rtp(iel,irh) )
893
796
    xmache =                                                      &
894
797
         ( coefa(ifac,iclu)*surfbo(1,ifac)                        &
895
798
         + coefa(ifac,iclv)*surfbo(2,ifac)                        &
896
 
         + coefa(ifac,iclw)*surfbo(3,ifac) ) /ra(isrfbn+ifac-1)   &
897
 
         / sqrt( gamagp * rtp(iel,ipriph) / rtp(iel,irhiph) )
 
799
         + coefa(ifac,iclw)*surfbo(3,ifac) ) /surfbn(ifac)        &
 
800
         / sqrt( gamagp * rtp(iel,ipr) / rtp(iel,irh) )
898
801
    dxmach = xmachi - xmache
899
802
 
900
803
    ! Pressure: rarefaction wave (Rusanov)
901
804
    if(dxmach.le.0.d0) then
902
805
 
903
806
      if(dxmach.gt.2.d0/(1.d0-gamagp)) then
904
 
        coefa(ifac,iclp) = rtp(iel,ipriph)*                       &
 
807
        coefa(ifac,iclp) = rtp(iel,ipr)*                       &
905
808
             ( (1.d0 + (gamagp-1.d0)*0.50d0*dxmach)               &
906
809
               ** (2.d0*gamagp/(gamagp-1.d0))    )
907
810
      elseif(dxmach.le.2.d0/(1.d0-gamagp) ) then
910
813
 
911
814
      ! Pressure: shock (Rusanov)
912
815
    else
913
 
      coefa(ifac,iclp) = rtp(iel,ipriph)*                         &
 
816
      coefa(ifac,iclp) = rtp(iel,ipr)*                         &
914
817
           (  1.d0 + gamagp*dxmach                                &
915
818
           *( (gamagp+1.d0)*0.25d0*dxmach                         &
916
819
           + sqrt(1.d0 + (gamagp+1.d0)**2/16.d0*dxmach**2) )  )
917
820
    endif
918
821
 
919
822
    ! This choice overrides the previous Rusanov choice
920
 
    coefa(ifac,iclp) = rtp(iel,ipriph)
 
823
    coefa(ifac,iclp) = rtp(iel,ipr)
921
824
 
922
825
    ! Total energy
923
826
    coefa(ifac,icle) =                                            &
960
863
    iel  = ifabor(ifac)
961
864
 
962
865
    ! Rarefaction case
963
 
    if(coefa(ifac,iclp).le.rtp(iel,ipriph)) then
 
866
    if(coefa(ifac,iclp).le.rtp(iel,ipr)) then
964
867
 
965
868
      ! Density
966
 
      coefa(ifac,iclr) = rtp(iel,irhiph)                          &
967
 
           * (coefa(ifac,iclp)/rtp(iel,ipriph))**(1.d0/gamagp)
 
869
      coefa(ifac,iclr) = rtp(iel,irh)                          &
 
870
           * (coefa(ifac,iclp)/rtp(iel,ipr))**(1.d0/gamagp)
968
871
 
969
872
      ! Velocity
970
 
      coefa(ifac,iclu) = rtp(iel,iuiph)                           &
971
 
           + 2.d0/(gamagp-1.d0)                                   &
972
 
           * sqrt(gamagp*rtp(iel,ipriph)/rtp(iel,irhiph))         &
973
 
           * (1.d0-(coefa(ifac,iclp)/rtp(iel,ipriph)              &
974
 
                        )**((gamagp-1.d0)/(2.d0*gamagp)))         &
975
 
           * surfbo(1,ifac)/ra(isrfbn+ifac-1)
976
 
 
977
 
      coefa(ifac,iclv) = rtp(iel,iviph)                           &
978
 
           + 2.d0/(gamagp-1.d0)                                   &
979
 
           * sqrt( gamagp*rtp(iel,ipriph)/rtp(iel,irhiph))        &
980
 
           * (1.d0-(coefa(ifac,iclp)/rtp(iel,ipriph)              &
981
 
                        )**((gamagp-1.d0)/(2.d0*gamagp)))         &
982
 
           * surfbo(2,ifac)/ra(isrfbn+ifac-1)
983
 
 
984
 
      coefa(ifac,iclw) = rtp(iel,iwiph)                           &
985
 
           + 2.d0/(gamagp-1.d0)                                   &
986
 
           * sqrt( gamagp*rtp(iel,ipriph)/rtp(iel,irhiph))        &
987
 
           * (1.d0-(coefa(ifac,iclp)/rtp(iel,ipriph)              &
 
873
      coefa(ifac,iclu) = rtp(iel,iu)                           &
 
874
           + 2.d0/(gamagp-1.d0)                                   &
 
875
           * sqrt(gamagp*rtp(iel,ipr)/rtp(iel,irh))         &
 
876
           * (1.d0-(coefa(ifac,iclp)/rtp(iel,ipr)              &
 
877
                        )**((gamagp-1.d0)/(2.d0*gamagp)))         &
 
878
           * surfbo(1,ifac)/surfbn(ifac)
 
879
 
 
880
      coefa(ifac,iclv) = rtp(iel,iv)                           &
 
881
           + 2.d0/(gamagp-1.d0)                                   &
 
882
           * sqrt( gamagp*rtp(iel,ipr)/rtp(iel,irh))        &
 
883
           * (1.d0-(coefa(ifac,iclp)/rtp(iel,ipr)              &
 
884
                        )**((gamagp-1.d0)/(2.d0*gamagp)))         &
 
885
           * surfbo(2,ifac)/surfbn(ifac)
 
886
 
 
887
      coefa(ifac,iclw) = rtp(iel,iw)                           &
 
888
           + 2.d0/(gamagp-1.d0)                                   &
 
889
           * sqrt( gamagp*rtp(iel,ipr)/rtp(iel,irh))        &
 
890
           * (1.d0-(coefa(ifac,iclp)/rtp(iel,ipr)              &
988
891
                        )**((gamagp-1.d0)/(2.d0/gamagp)))         &
989
 
           * surfbo(3,ifac)/ra(isrfbn+ifac-1)
 
892
           * surfbo(3,ifac)/surfbn(ifac)
990
893
 
991
894
      ! Total energy
992
895
      coefa(ifac,icle) =                                          &
998
901
    else
999
902
 
1000
903
      ! Density
1001
 
      coefa(ifac,iclr) = rtp(iel,irhiph)                          &
 
904
      coefa(ifac,iclr) = rtp(iel,irh)                          &
1002
905
           * ( (gamagp+1.d0)*coefa(ifac,iclp)                     &
1003
 
             + (gamagp-1.d0)*rtp(iel,ipriph) )                    &
 
906
             + (gamagp-1.d0)*rtp(iel,ipr) )                    &
1004
907
           / ( (gamagp-1.d0)*coefa(ifac,iclp)                     &
1005
 
             + (gamagp+1.d0)*rtp(iel,ipriph) )
 
908
             + (gamagp+1.d0)*rtp(iel,ipr) )
1006
909
 
1007
910
      ! Velocity
1008
 
      coefa(ifac,iclu) = rtp(iel,iuiph)                           &
1009
 
           - (coefa(ifac,iclp)-rtp(iel,ipriph))                   &
1010
 
           * sqrt(2.d0/                                           &
1011
 
                  (rtp(iel,irhiph)                                &
1012
 
                   *((gamagp+1.d0)*coefa(ifac,iclp)               &
1013
 
                    +(gamagp-1.d0)*rtp(iel,ipriph) )))            &
1014
 
           * surfbo(1,ifac)/ra(isrfbn+ifac-1)
1015
 
 
1016
 
      coefa(ifac,iclv) = rtp(iel,iviph)                           &
1017
 
           - (coefa(ifac,iclp)-rtp(iel,ipriph))                   &
1018
 
           * sqrt(2.d0/                                           &
1019
 
                  (rtp(iel,irhiph)                                &
1020
 
                   *((gamagp+1.d0)*coefa(ifac,iclp)               &
1021
 
                    +(gamagp-1.d0)*rtp(iel,ipriph) )))            &
1022
 
           * surfbo(2,ifac)/ra(isrfbn+ifac-1)
1023
 
 
1024
 
      coefa(ifac,iclw) = rtp(iel,iwiph)                           &
1025
 
           - (coefa(ifac,iclp)-rtp(iel,ipriph))                   &
1026
 
           * sqrt(2.d0/                                           &
1027
 
                  (rtp(iel,irhiph)                                &
1028
 
                   *((gamagp+1.d0)*coefa(ifac,iclp)               &
1029
 
                    +(gamagp-1.d0)*rtp(iel,ipriph) )))            &
1030
 
           * surfbo(3,ifac)/ra(isrfbn+ifac-1)
 
911
      coefa(ifac,iclu) = rtp(iel,iu)                           &
 
912
           - (coefa(ifac,iclp)-rtp(iel,ipr))                   &
 
913
           * sqrt(2.d0/                                           &
 
914
                  (rtp(iel,irh)                                &
 
915
                   *((gamagp+1.d0)*coefa(ifac,iclp)               &
 
916
                    +(gamagp-1.d0)*rtp(iel,ipr) )))            &
 
917
           * surfbo(1,ifac)/surfbn(ifac)
 
918
 
 
919
      coefa(ifac,iclv) = rtp(iel,iv)                           &
 
920
           - (coefa(ifac,iclp)-rtp(iel,ipr))                   &
 
921
           * sqrt(2.d0/                                           &
 
922
                  (rtp(iel,irh)                                &
 
923
                   *((gamagp+1.d0)*coefa(ifac,iclp)               &
 
924
                    +(gamagp-1.d0)*rtp(iel,ipr) )))            &
 
925
           * surfbo(2,ifac)/surfbn(ifac)
 
926
 
 
927
      coefa(ifac,iclw) = rtp(iel,iw)                           &
 
928
           - (coefa(ifac,iclp)-rtp(iel,ipr))                   &
 
929
           * sqrt(2.d0/                                           &
 
930
                  (rtp(iel,irh)                                &
 
931
                   *((gamagp+1.d0)*coefa(ifac,iclp)               &
 
932
                    +(gamagp-1.d0)*rtp(iel,ipr) )))            &
 
933
           * surfbo(3,ifac)/surfbn(ifac)
1031
934
 
1032
935
      ! Total energy
1033
936
      coefa(ifac,icle) =                                          &
1054
957
 
1055
958
    ! Energie totale
1056
959
    coefa(ifac,icle) =                                            &
1057
 
         cv0(iphas)*coefa(ifac,iclt)                              &
 
960
         cv0*coefa(ifac,iclt)                              &
1058
961
         + 0.5d0*( coefa(ifac,iclu)**2                            &
1059
962
                 + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2 )
1060
963
 
1072
975
 
1073
976
    ! Total energy
1074
977
    coefa(ifac,icle) =                                            &
1075
 
         cv0(iphas)*coefa(ifac,iclt)                              &
 
978
         cv0*coefa(ifac,iclt)                              &
1076
979
         + 0.5d0*( coefa(ifac,iclu)**2                            &
1077
980
                 + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2 )
1078
981
 
1108
1011
                                       *coefa(ifac,iclt)
1109
1012
 
1110
1013
    ! Total energy
1111
 
    coefa(ifac,icle) = cv0(iphas) * coefa(ifac,iclt)              &
 
1014
    coefa(ifac,icle) = cv0 * coefa(ifac,iclt)              &
1112
1015
         + 0.5d0*( coefa(ifac,iclu)**2                            &
1113
1016
                 + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2 )
1114
1017
 
1143
1046
 
1144
1047
! This section requires further checking and testing
1145
1048
 
1146
 
elseif(ieos(iphas).eq.2) then
 
1049
elseif(ieos.eq.2) then
1147
1050
 
1148
1051
!===============================================================================
1149
1052
 
1164
1067
! --- Ex. 1: Perfect gas containing 3 components
1165
1068
!     Molar mass, gamma
1166
1069
 
1167
 
    ! Phase
1168
 
  if(iphas.eq.1) then
1169
 
 
1170
1070
    ! Molar mass of the components (kg/mol)
1171
1071
    cstgr(1)  = 18.d-3
1172
1072
    cstgr(2)  = 32.d-3
1183
1083
 
1184
1084
      ! Calculation of the equivalent gamma of the mixture at cell centers
1185
1085
      do iel = 1, ncel
1186
 
        gamagr(iel) = propce(iel,ipproc(icp(iphas)))              &
1187
 
           / ( propce(iel,ipproc(icp(iphas))) - rr/xmasmr(iel) )
 
1086
        gamagr(iel) = propce(iel,ipproc(icp))              &
 
1087
           / ( propce(iel,ipproc(icp)) - rr/xmasmr(iel) )
1188
1088
      enddo
1189
1089
 
1190
1090
    endif
1191
1091
 
1192
1092
  endif
1193
1093
 
1194
 
  endif
1195
 
 
1196
1094
!-------------------------------------------------------------------------------
1197
1095
 
1198
1096
! End of the examples
1219
1117
 
1220
1118
  if(iccfth.eq.-1) then
1221
1119
 
1222
 
    icp(iphas) = 1
1223
 
    cp0(iphas) = epzero
1224
 
    icv(iphas) = 1
1225
 
    cv0(iphas) = epzero
 
1120
    icp = 1
 
1121
    cp0 = epzero
 
1122
    icv = 1
 
1123
    cv0 = epzero
1226
1124
 
1227
1125
 
1228
1126
! Default initializations
1230
1128
  elseif(iccfth.eq.0) then
1231
1129
 
1232
1130
    do iel = 1, ncel
1233
 
      propce(iel,ipproc(icp(iphas))) = cp0(iphas)
1234
 
      propce(iel,ipproc(icv(iphas))) =                            &
1235
 
           cp0(iphas) - rr/xmasmr(iel)
1236
 
      rtp(iel,irhiph) = p0(iphas)*xmasmr(iel)/rr/t0(iphas)
1237
 
      rtp(iel,ieniph) =                                           &
1238
 
           propce(iel,ipproc(icv(iphas)))*t0(iphas)
 
1131
      propce(iel,ipproc(icp)) = cp0
 
1132
      propce(iel,ipproc(icv)) =                            &
 
1133
           cp0 - rr/xmasmr(iel)
 
1134
      rtp(iel,irh) = p0*xmasmr(iel)/rr/t0
 
1135
      rtp(iel,ien) = propce(iel,ipproc(icv))*t0
1239
1136
    enddo
1240
1137
 
1241
1138
 
1247
1144
 
1248
1145
      ! Temperature
1249
1146
      sorti1(iel) =                                               &
1250
 
           xmasmr(iel)/rr*rtp(iel,ipriph)/rtp(iel,irhiph)
 
1147
           xmasmr(iel)/rr*rtp(iel,ipr)/rtp(iel,irh)
1251
1148
 
1252
1149
      ! Total energy
1253
 
      sorti2(iel) = propce(iel,ipproc(icv(iphas)))*sorti1(iel)    &
1254
 
    + 0.5d0*( rtp(iel,iuiph)**2                                   &
1255
 
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 )
 
1150
      sorti2(iel) = propce(iel,ipproc(icv))*sorti1(iel)    &
 
1151
    + 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 )
1256
1152
 
1257
1153
    enddo
1258
1154
 
1259
1155
    ! Transfer to the array rtp
1260
1156
    if(imodif.gt.0) then
1261
1157
      do iel = 1, ncel
1262
 
        rtp(iel,itkiph) = sorti1(iel)
1263
 
        rtp(iel,ieniph) = sorti2(iel)
 
1158
        rtp(iel,itk) = sorti1(iel)
 
1159
        rtp(iel,ien) = sorti2(iel)
1264
1160
      enddo
1265
1161
    endif
1266
1162
 
1273
1169
 
1274
1170
      ! Density
1275
1171
      sorti1(iel) =                                               &
1276
 
           xmasmr(iel)/rr*rtp(iel,ipriph)/rtp(iel,itkiph)
 
1172
           xmasmr(iel)/rr*rtp(iel,ipr)/rtp(iel,itk)
1277
1173
 
1278
1174
      ! Total energy
1279
1175
      sorti2(iel) =                                               &
1280
 
           propce(iel,ipproc(icv(iphas)))*rtp(iel,itkiph)         &
1281
 
    + 0.5d0*( rtp(iel,iuiph)**2                                   &
1282
 
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 )
 
1176
           propce(iel,ipproc(icv))*rtp(iel,itk)         &
 
1177
    + 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 )
1283
1178
 
1284
1179
    enddo
1285
1180
 
1286
1181
    ! Transfer to the array rtp
1287
1182
    if(imodif.gt.0) then
1288
1183
      do iel = 1, ncel
1289
 
        rtp(iel,irhiph) = sorti1(iel)
1290
 
        rtp(iel,ieniph) = sorti2(iel)
 
1184
        rtp(iel,irh) = sorti1(iel)
 
1185
        rtp(iel,ien) = sorti2(iel)
1291
1186
      enddo
1292
1187
    endif
1293
1188
 
1300
1195
 
1301
1196
      ! Density
1302
1197
      sorti1(iel) =                                               &
1303
 
           rtp(iel,ipriph)/(gamagr(iel)-1.d0)/( rtp(iel,ieniph)   &
1304
 
  - 0.5d0*( rtp(iel,iuiph)**2                                     &
1305
 
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 ) )
 
1198
           rtp(iel,ipr)/(gamagr(iel)-1.d0)/( rtp(iel,ien)   &
 
1199
  - 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 ) )
1306
1200
 
1307
1201
      ! Temperature
1308
 
      sorti2(iel) = xmasmr(iel)/rr*rtp(iel,ipriph)/sorti1(iel)
 
1202
      sorti2(iel) = xmasmr(iel)/rr*rtp(iel,ipr)/sorti1(iel)
1309
1203
 
1310
1204
    enddo
1311
1205
 
1312
1206
    ! Transfer to the array rtp
1313
1207
    if(imodif.gt.0) then
1314
1208
      do iel = 1, ncel
1315
 
        rtp(iel,irhiph) = sorti1(iel)
1316
 
        rtp(iel,itkiph) = sorti2(iel)
 
1209
        rtp(iel,irh) = sorti1(iel)
 
1210
        rtp(iel,itk) = sorti2(iel)
1317
1211
      enddo
1318
1212
    endif
1319
1213
 
1326
1220
 
1327
1221
      ! Pressure
1328
1222
      sorti1(iel) =                                               &
1329
 
           rtp(iel,irhiph)*rr/xmasmr(iel)*rtp(iel,itkiph)
 
1223
           rtp(iel,irh)*rr/xmasmr(iel)*rtp(iel,itk)
1330
1224
 
1331
1225
      ! Total energy
1332
1226
      sorti2(iel) =                                               &
1333
 
           propce(iel,ipproc(icv(iphas)))*rtp(iel,itkiph)         &
1334
 
    + 0.5d0*( rtp(iel,iuiph)**2                                   &
1335
 
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 )
 
1227
           propce(iel,ipproc(icv))*rtp(iel,itk)         &
 
1228
    + 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 )
1336
1229
 
1337
1230
    enddo
1338
1231
 
1339
1232
    ! Transfer to the array rtp
1340
1233
    if(imodif.gt.0) then
1341
1234
      do iel = 1, ncel
1342
 
        rtp(iel,ipriph) = sorti1(iel)
1343
 
        rtp(iel,ieniph) = sorti2(iel)
 
1235
        rtp(iel,ipr) = sorti1(iel)
 
1236
        rtp(iel,ien) = sorti2(iel)
1344
1237
      enddo
1345
1238
    endif
1346
1239
 
1353
1246
 
1354
1247
      ! Pressure
1355
1248
      sorti1(iel) =                                               &
1356
 
           (gamagr(iel)-1.d0)*rtp(iel,irhiph)*( rtp(iel,ieniph)   &
1357
 
  - 0.5d0*( rtp(iel,iuiph)**2                                     &
1358
 
           + rtp(iel,iviph)**2 + rtp(iel,iwiph)**2 ) )
 
1249
           (gamagr(iel)-1.d0)*rtp(iel,irh)*( rtp(iel,ien)   &
 
1250
  - 0.5d0*( rtp(iel,iu)**2 + rtp(iel,iv)**2 + rtp(iel,iw)**2 ) )
1359
1251
 
1360
1252
      ! Temperature
1361
 
      sorti2(iel) = xmasmr(iel)/rr*sorti1(iel)/rtp(iel,irhiph)
 
1253
      sorti2(iel) = xmasmr(iel)/rr*sorti1(iel)/rtp(iel,irh)
1362
1254
 
1363
1255
    enddo
1364
1256
 
1365
1257
    ! Transfer to the array rtp
1366
1258
    if(imodif.gt.0) then
1367
1259
      do iel = 1, ncel
1368
 
        rtp(iel,ipriph) = sorti1(iel)
1369
 
        rtp(iel,itkiph) = sorti2(iel)
 
1260
        rtp(iel,ipr) = sorti1(iel)
 
1261
        rtp(iel,itk) = sorti2(iel)
1370
1262
      enddo
1371
1263
    endif
1372
1264
 
1379
1271
    do iel = 1, ncel
1380
1272
 
1381
1273
      ! Verification of the positivity of the pressure
1382
 
      if(rtp(iel,ipriph).lt.0.d0) then
1383
 
        write(nfecra,1110) iel , rtp(iel,ipriph)
 
1274
      if(rtp(iel,ipr).lt.0.d0) then
 
1275
        write(nfecra,1110) iel , rtp(iel,ipr)
1384
1276
        ierr = 1
1385
1277
 
1386
1278
      ! Verification of the positivity of the density
1387
 
      elseif(rtp(iel,irhiph).le.0.d0) then
1388
 
        write(nfecra,1120) iel , rtp(iel,irhiph)
 
1279
      elseif(rtp(iel,irh).le.0.d0) then
 
1280
        write(nfecra,1120) iel , rtp(iel,irh)
1389
1281
        ierr = 1
1390
1282
 
1391
1283
      else
1392
1284
 
1393
1285
        ! Computation
1394
1286
        sorti1(iel) =                                             &
1395
 
             gamagr(iel) * rtp(iel,ipriph) / rtp(iel,irhiph)
 
1287
             gamagr(iel) * rtp(iel,ipr) / rtp(iel,irh)
1396
1288
 
1397
1289
      endif
1398
1290
 
1410
1302
    do iel = 1, ncel
1411
1303
 
1412
1304
      ! Verification of the positivity of the density
1413
 
      if(rtp(iel,irhiph).lt.0.d0) then
1414
 
        write(nfecra,1220) iel , rtp(iel,irhiph)
 
1305
      if(rtp(iel,irh).lt.0.d0) then
 
1306
        write(nfecra,1220) iel , rtp(iel,irh)
1415
1307
        ierr = 1
1416
1308
 
1417
1309
      else
1418
1310
 
1419
1311
        ! Computation
1420
 
        sorti1(iel) = rtp(iel,irhiph)**gamagr(iel)
 
1312
        sorti1(iel) = rtp(iel,irh)**gamagr(iel)
1421
1313
 
1422
1314
      endif
1423
1315
 
1433
1325
 
1434
1326
    do iel = 1, ncel
1435
1327
 
1436
 
      sorti1(iel) = propce(iel,ipproc(icp(iphas)))-rr/xmasmr(iel)
 
1328
      sorti1(iel) = propce(iel,ipproc(icp))-rr/xmasmr(iel)
1437
1329
 
1438
1330
    enddo
1439
1331
 
1450
1342
    do iel = 1, ncel
1451
1343
 
1452
1344
      ! Verification of the positivity of the pressure
1453
 
      if(rtp(iel,ipriph).lt.0.d0) then
1454
 
        write(nfecra,1310) iel , rtp(iel,ipriph)
 
1345
      if(rtp(iel,ipr).lt.0.d0) then
 
1346
        write(nfecra,1310) iel , rtp(iel,ipr)
1455
1347
        ierr = 1
1456
1348
 
1457
1349
      ! Verification of the positivity of the density
1458
 
      elseif(rtp(iel,irhiph).le.0.d0) then
1459
 
        write(nfecra,1320) iel , rtp(iel,irhiph)
 
1350
      elseif(rtp(iel,irh).le.0.d0) then
 
1351
        write(nfecra,1320) iel , rtp(iel,irh)
1460
1352
        ierr = 1
1461
1353
 
1462
1354
      else
1463
1355
 
1464
1356
        ! Computation
1465
1357
        sorti1(iel) =                                             &
1466
 
             rtp(iel,ipriph) / (rtp(iel,irhiph)**gamagr(iel))
 
1358
             rtp(iel,ipr) / (rtp(iel,irh)**gamagr(iel))
1467
1359
 
1468
1360
      endif
1469
1361
 
1504
1396
 
1505
1397
    ! Calculation of the Mach number at the boundary face, using the
1506
1398
    !   cell center velocity projected on the vector normal to the boundary
1507
 
    xmach = ( rtp(iel,iuiph)*surfbo(1,ifac)                       &
1508
 
           + rtp(iel,iviph)*surfbo(2,ifac)                        &
1509
 
           + rtp(iel,iwiph)*surfbo(3,ifac) ) / ra(isrfbn+ifac-1)  &
1510
 
         / sqrt( gamagr(iel)*rtp(iel,ipriph)/rtp(iel,irhiph) )
 
1399
    xmach = ( rtp(iel,iu)*surfbo(1,ifac)                       &
 
1400
           + rtp(iel,iv)*surfbo(2,ifac)                        &
 
1401
           + rtp(iel,iw)*surfbo(3,ifac) ) / surfbn(ifac)       &
 
1402
         / sqrt( gamagr(iel)*rtp(iel,ipr)/rtp(iel,irh) )
1511
1403
 
1512
1404
    coefa(ifac,iclp) = 0.d0
1513
1405
 
1529
1421
            *( (gamagr(iel)+1.d0)/4.d0*xmach                      &
1530
1422
           + sqrt(1.d0 + (gamagr(iel)+1.d0)**2/16.d0*xmach**2) )
1531
1423
      coefb(ifac,iclt) = coefb(ifac,iclp)/(1.d0-coefb(ifac,iclp)) &
1532
 
          / rtp(iel,ipriph) * ( rtp(iel,irhiph)                   &
1533
 
              * (rtp(iel,iuiph)**2                                &
1534
 
                +rtp(iel,iviph)**2+rtp(iel,iwiph)**2)             &
1535
 
              + rtp(iel,ipriph) *(1.d0-coefb(ifac,iclp)) )
 
1424
          / rtp(iel,ipr) * ( rtp(iel,irh)                         &
 
1425
              * (rtp(iel,iu)**2+rtp(iel,iv)**2+rtp(iel,iw)**2)    &
 
1426
              + rtp(iel,ipr) *(1.d0-coefb(ifac,iclp)) )
1536
1427
    endif
1537
1428
 
1538
1429
    ! Total energy: 'internal energy - Cv T'
1552
1443
 
1553
1444
    ! Calculation of the Mach number at the boundary face, using the
1554
1445
    !   cell center velocity projected on the vector normal to the boundary
1555
 
    xmachi = ( rtp(iel,iuiph)*surfbo(1,ifac)                      &
1556
 
         + rtp(iel,iviph)*surfbo(2,ifac)                          &
1557
 
         + rtp(iel,iwiph)*surfbo(3,ifac) )/ra(isrfbn+ifac-1)      &
1558
 
         / sqrt(gamagr(iel)*rtp(iel,ipriph)/rtp(iel,irhiph))
 
1446
    xmachi = ( rtp(iel,iu)*surfbo(1,ifac)                      &
 
1447
         + rtp(iel,iv)*surfbo(2,ifac)                          &
 
1448
         + rtp(iel,iw)*surfbo(3,ifac) )/surfbn(ifac)           &
 
1449
         / sqrt(gamagr(iel)*rtp(iel,ipr)/rtp(iel,irh))
1559
1450
    xmache = ( coefa(ifac,iclu)*surfbo(1,ifac)                    &
1560
1451
         + coefa(ifac,iclv)*surfbo(2,ifac)                        &
1561
 
         + coefa(ifac,iclw)*surfbo(3,ifac) )/ra(isrfbn+ifac-1)    &
1562
 
         / sqrt(gamagr(iel)*rtp(iel,ipriph)/rtp(iel,irhiph))
 
1452
         + coefa(ifac,iclw)*surfbo(3,ifac) )/surfbn(ifac)         &
 
1453
         / sqrt(gamagr(iel)*rtp(iel,ipr)/rtp(iel,irh))
1563
1454
    dxmach = xmachi - xmache
1564
1455
 
1565
1456
    ! Pressure: rarefaction wave
1566
1457
    if(dxmach.le.0.d0) then
1567
1458
 
1568
1459
      if(dxmach.gt.2.d0/(1.d0-gamagr(iel))) then
1569
 
        coefa(ifac,iclp) = rtp(iel,ipriph)*                       &
 
1460
        coefa(ifac,iclp) = rtp(iel,ipr)*                       &
1570
1461
             ( (1.d0 + (gamagr(iel)-1.d0)*0.50d0*dxmach)          &
1571
1462
               ** (2.d0*gamagr(iel)/(gamagr(iel)-1.d0))  )
1572
1463
      elseif(dxmach.le.2.d0/(1.d0-gamagr(iel)) ) then
1575
1466
 
1576
1467
    ! Pressure: shock
1577
1468
    else
1578
 
      coefa(ifac,iclp) = rtp(iel,ipriph)*                         &
 
1469
      coefa(ifac,iclp) = rtp(iel,ipr)*                         &
1579
1470
           (  1.d0 + gamagr(iel)*dxmach                           &
1580
1471
           *( (gamagr(iel)+1.d0)*0.25d0*dxmach                    &
1581
1472
           + sqrt(1.d0 + (gamagr(iel)+1.d0)**2/16.d0              &
1583
1474
    endif
1584
1475
 
1585
1476
    ! This choice overrides the previous Rusanov choice
1586
 
    coefa(ifac,iclp) = rtp(iel,ipriph)
 
1477
    coefa(ifac,iclp) = rtp(iel,ipr)
1587
1478
 
1588
1479
    ! Total energy
1589
1480
    coefa(ifac,icle) =                                            &
1600
1491
 
1601
1492
    ! Calculation of the Mach number at the boundary face, using the
1602
1493
    !   cell center velocity projected on the vector normal to the boundary
1603
 
    xmach = ( rtp(iel,iuiph)*surfbo(1,ifac)                       &
1604
 
           + rtp(iel,iviph)*surfbo(2,ifac)                        &
1605
 
           + rtp(iel,iwiph)*surfbo(3,ifac) ) / ra(isrfbn+ifac-1)  &
1606
 
         / sqrt(gamagr(iel)*rtp(iel,ipriph)/rtp(iel,irhiph))
 
1494
    xmach = ( rtp(iel,iu)*surfbo(1,ifac)                       &
 
1495
           + rtp(iel,iv)*surfbo(2,ifac)                        &
 
1496
           + rtp(iel,iw)*surfbo(3,ifac) ) / surfbn(ifac)       &
 
1497
         / sqrt(gamagr(iel)*rtp(iel,ipr)/rtp(iel,irh))
1607
1498
 
1608
1499
    ! Supersonic outlet: Dirichlet for all variables
1609
1500
    if(xmach.ge.1.d0) then
1613
1504
 
1614
1505
      ! Entropy
1615
1506
      coefa(ifac,iclt) =                                          &
1616
 
           rtp(iel,ipriph)/rtp(iel,irhiph)**gamagr(iel)
 
1507
           rtp(iel,ipr)/rtp(iel,irh)**gamagr(iel)
1617
1508
 
1618
1509
    ! Subsonic outlet
1619
1510
    elseif(xmach.lt.1.d0 .and. xmach.ge.0.d0) then
1620
1511
 
1621
1512
      ! Rarefaction:
1622
 
      if(coefa(ifac,iclp).le.rtp(iel,ipriph)) then
 
1513
      if(coefa(ifac,iclp).le.rtp(iel,ipr)) then
1623
1514
 
1624
1515
        ! Density
1625
 
        coefa(ifac,iclr) = rtp(iel,irhiph)                        &
1626
 
             * (coefa(ifac,iclp)/rtp(iel,ipriph))                 &
 
1516
        coefa(ifac,iclr) = rtp(iel,irh)                        &
 
1517
             * (coefa(ifac,iclp)/rtp(iel,ipr))                 &
1627
1518
                **(1.d0/gamagr(iel))
1628
1519
 
1629
1520
        ! Velocity
1630
 
        coefa(ifac,iclu) = rtp(iel,iuiph)                         &
1631
 
             + 2.d0/(gamagr(iel)-1.d0)                            &
1632
 
 * sqrt( gamagr(iel) * rtp(iel,ipriph) / rtp(iel,irhiph) )        &
1633
 
             * ( 1.d0                                             &
1634
 
 - (coefa(ifac,iclp)/rtp(iel,ipriph))                             &
1635
 
               **((gamagr(iel)-1.d0)/2.d0/gamagr(iel)) )          &
1636
 
 * surfbo(1,ifac) / ra(isrfbn+ifac-1)
1637
 
 
1638
 
        coefa(ifac,iclv) = rtp(iel,iviph)                         &
1639
 
             + 2.d0/(gamagr(iel)-1.d0)                            &
1640
 
 * sqrt( gamagr(iel) * rtp(iel,ipriph) / rtp(iel,irhiph) )        &
1641
 
             * ( 1.d0                                             &
1642
 
 - (coefa(ifac,iclp)/rtp(iel,ipriph))                             &
1643
 
               **((gamagr(iel)-1.d0)/2.d0/gamagr(iel)) )          &
1644
 
 * surfbo(2,ifac) / ra(isrfbn+ifac-1)
1645
 
 
1646
 
        coefa(ifac,iclw) = rtp(iel,iwiph)                         &
1647
 
             + 2.d0/(gamagr(iel)-1.d0)                            &
1648
 
 * sqrt( gamagr(iel) * rtp(iel,ipriph) / rtp(iel,irhiph) )        &
1649
 
             * ( 1.d0                                             &
1650
 
 - (coefa(ifac,iclp)/rtp(iel,ipriph))                             &
1651
 
               **((gamagr(iel)-1.d0)/2.d0/gamagr(iel)) )          &
1652
 
 * surfbo(3,ifac) / ra(isrfbn+ifac-1)
 
1521
        coefa(ifac,iclu) = rtp(iel,iu)                         &
 
1522
             + 2.d0/(gamagr(iel)-1.d0)                            &
 
1523
 * sqrt( gamagr(iel) * rtp(iel,ipr) / rtp(iel,irh) )        &
 
1524
             * ( 1.d0                                             &
 
1525
 - (coefa(ifac,iclp)/rtp(iel,ipr))                             &
 
1526
               **((gamagr(iel)-1.d0)/2.d0/gamagr(iel)) )          &
 
1527
 * surfbo(1,ifac) / surfbn(ifac)
 
1528
 
 
1529
        coefa(ifac,iclv) = rtp(iel,iv)                         &
 
1530
             + 2.d0/(gamagr(iel)-1.d0)                            &
 
1531
 * sqrt( gamagr(iel) * rtp(iel,ipr) / rtp(iel,irh) )        &
 
1532
             * ( 1.d0                                             &
 
1533
 - (coefa(ifac,iclp)/rtp(iel,ipr))                             &
 
1534
               **((gamagr(iel)-1.d0)/2.d0/gamagr(iel)) )          &
 
1535
 * surfbo(2,ifac) / surfbn(ifac)
 
1536
 
 
1537
        coefa(ifac,iclw) = rtp(iel,iw)                         &
 
1538
             + 2.d0/(gamagr(iel)-1.d0)                            &
 
1539
 * sqrt( gamagr(iel) * rtp(iel,ipr) / rtp(iel,irh) )        &
 
1540
             * ( 1.d0                                             &
 
1541
 - (coefa(ifac,iclp)/rtp(iel,ipr))                             &
 
1542
               **((gamagr(iel)-1.d0)/2.d0/gamagr(iel)) )          &
 
1543
 * surfbo(3,ifac) / surfbn(ifac)
1653
1544
 
1654
1545
        ! Total energy
1655
1546
        coefa(ifac,icle) = coefa(ifac,iclp)                       &
1666
1557
      else
1667
1558
 
1668
1559
        ! Density
1669
 
        coefa(ifac,iclr) = rtp(iel,irhiph)                        &
 
1560
        coefa(ifac,iclr) = rtp(iel,irh)                        &
1670
1561
 * ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
1671
 
   + (gamagr(iel)-1.d0)*rtp(iel,ipriph) )                         &
 
1562
   + (gamagr(iel)-1.d0)*rtp(iel,ipr) )                         &
1672
1563
 / ( (gamagr(iel)-1.d0)*coefa(ifac,iclp)                          &
1673
 
   + (gamagr(iel)+1.d0)*rtp(iel,ipriph) )
 
1564
   + (gamagr(iel)+1.d0)*rtp(iel,ipr) )
1674
1565
 
1675
1566
        ! Velocity
1676
 
        coefa(ifac,iclu) = rtp(iel,iuiph)                         &
1677
 
 - (coefa(ifac,iclp)-rtp(iel,ipriph))*sqrt(2.d0/rtp(iel,irhiph)   &
1678
 
 / ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
1679
 
   + (gamagr(iel)-1.d0)*rtp(iel,ipriph) ))                        &
1680
 
 * surfbo(1,ifac) / ra(isrfbn+ifac-1)
1681
 
 
1682
 
        coefa(ifac,iclv) = rtp(iel,iviph)                         &
1683
 
 - (coefa(ifac,iclp)-rtp(iel,ipriph))*sqrt(2.d0/rtp(iel,irhiph)   &
1684
 
 / ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
1685
 
   + (gamagr(iel)-1.d0)*rtp(iel,ipriph) ))                        &
1686
 
 * surfbo(2,ifac) / ra(isrfbn+ifac-1)
1687
 
 
1688
 
        coefa(ifac,iclw) = rtp(iel,iwiph)                         &
1689
 
 - (coefa(ifac,iclp)-rtp(iel,ipriph))*sqrt(2.d0/rtp(iel,irhiph)   &
1690
 
 / ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
1691
 
   + (gamagr(iel)-1.d0)*rtp(iel,ipriph) ))                        &
1692
 
 * surfbo(3,ifac) / ra(isrfbn+ifac-1)
 
1567
        coefa(ifac,iclu) = rtp(iel,iu)                         &
 
1568
 - (coefa(ifac,iclp)-rtp(iel,ipr))*sqrt(2.d0/rtp(iel,irh)   &
 
1569
 / ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
 
1570
   + (gamagr(iel)-1.d0)*rtp(iel,ipr) ))                        &
 
1571
 * surfbo(1,ifac) / surfbn(ifac)
 
1572
 
 
1573
        coefa(ifac,iclv) = rtp(iel,iv)                         &
 
1574
 - (coefa(ifac,iclp)-rtp(iel,ipr))*sqrt(2.d0/rtp(iel,irh)   &
 
1575
 / ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
 
1576
   + (gamagr(iel)-1.d0)*rtp(iel,ipr) ))                        &
 
1577
 * surfbo(2,ifac) / surfbn(ifac)
 
1578
 
 
1579
        coefa(ifac,iclw) = rtp(iel,iw)                         &
 
1580
 - (coefa(ifac,iclp)-rtp(iel,ipr))*sqrt(2.d0/rtp(iel,irh)   &
 
1581
 / ( (gamagr(iel)+1.d0)*coefa(ifac,iclp)                          &
 
1582
   + (gamagr(iel)-1.d0)*rtp(iel,ipr) ))                        &
 
1583
 * surfbo(3,ifac) / surfbn(ifac)
1693
1584
 
1694
1585
        ! Total energy
1695
1586
        coefa(ifac,icle) = coefa(ifac,iclp)                       &
1723
1614
                                        /coefa(ifac,iclr)
1724
1615
 
1725
1616
    ! Total energy
1726
 
    coefa(ifac,icle) = propce(iel,ipproc(icv(iphas)))             &
 
1617
    coefa(ifac,icle) = propce(iel,ipproc(icv))             &
1727
1618
               * coefa(ifac,iclt) + 0.5d0*( coefa(ifac,iclu)**2   &
1728
1619
                    + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2)
1729
1620
 
1740
1631
                                       /coefa(ifac,iclt)
1741
1632
 
1742
1633
    ! Total energy
1743
 
    coefa(ifac,icle) = propce(iel,ipproc(icv(iphas)))             &
 
1634
    coefa(ifac,icle) = propce(iel,ipproc(icv))             &
1744
1635
               * coefa(ifac,iclt) + 0.5d0*( coefa(ifac,iclu)**2   &
1745
1636
                    + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2)
1746
1637
 
1774
1665
                                       *coefa(ifac,iclt)
1775
1666
 
1776
1667
    ! Total energy
1777
 
    coefa(ifac,icle) = propce(iel,ipproc(icv(iphas)))             &
 
1668
    coefa(ifac,icle) = propce(iel,ipproc(icv))             &
1778
1669
               * coefa(ifac,iclt) + 0.5d0*( coefa(ifac,iclu)**2   &
1779
1670
                    + coefa(ifac,iclv)**2 + coefa(ifac,iclw)**2)
1780
1671
 
1952
1843
'@     by setting a minimum value for the density variable in',/, &
1953
1844
'@     the GUI or in the user subroutine ''usini1'' (set the ',/, &
1954
1845
'@     scamin value associated to the variable ',/,               &
1955
 
'@     isca(irho(iphas)).',/,                                     &
 
1846
'@     isca(irho).',/,                                     &
1956
1847
'@',/,                                                            &
1957
1848
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
1958
1849
'@',/)
1959
 
 8100 format (                                                          &
 
1850
 8100 format (                                                    &
1960
1851
'@',/,                                                            &
1961
1852
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
1962
1853
'@',/,                                                            &
1978
1869
! gamma variable option will have been fixed
1979
1870
 
1980
1871
 
1981
 
 1110 format(                                                           &
 
1872
 1110 format(                                                     &
1982
1873
'@',/,                                                            &
1983
1874
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
1984
1875
'@',/,                                                            &
1994
1885
'@',/,                                                            &
1995
1886
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
1996
1887
'@',/)
1997
 
 1120 format(                                                           &
 
1888
 1120 format(                                                     &
1998
1889
'@',/,                                                            &
1999
1890
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2000
1891
'@',/,                                                            &
2010
1901
'@',/,                                                            &
2011
1902
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2012
1903
'@',/)
2013
 
 1220 format(                                                           &
 
1904
 1220 format(                                                     &
2014
1905
'@',/,                                                            &
2015
1906
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2016
1907
'@',/,                                                            &
2026
1917
'@',/,                                                            &
2027
1918
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2028
1919
'@',/)
2029
 
 1310 format(                                                           &
 
1920
 1310 format(                                                     &
2030
1921
'@',/,                                                            &
2031
1922
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2032
1923
'@',/,                                                            &
2042
1933
'@',/,                                                            &
2043
1934
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2044
1935
'@',/)
2045
 
 1320 format(                                                           &
 
1936
 1320 format(                                                     &
2046
1937
'@',/,                                                            &
2047
1938
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
2048
1939
'@',/,                                                            &