-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathnewx2.src
More file actions
124 lines (124 loc) · 4.13 KB
/
Copy pathnewx2.src
File metadata and controls
124 lines (124 loc) · 4.13 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include "zeus2d.def"
c=======================================================================
c////////////////////////// SUBROUTINE NEWX2 \\\\\\\\\\\\\\\\\\\\\\\\\
c
subroutine newx2
c
c PURPOSE: Computes "new" x2 grid variables (grid variables at advanced
c time) to be used in TRANSPRT. Grid values are calculated for
c j=js-2 to je+2, except for dvl2a (j=js,je+2) and dvl2b (j=js+1,je+2).
c Note similarity of the expressions used to those in the grid
c generator routine GRIDI.
c
c EXTERNALS: [none]
c
c LOCALS:
c-----------------------------------------------------------------------
implicit NONE
#include "param.h"
#include "grid.h"
#include "root.h"
integer j
REAL qa,qb,qc,qd
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\//////////////////////////////////
c=======================================================================
x2an(js-2) = x1an(js-2) + vg2(js-2)*dt
do 10 j=js-1,je+2
x2an(j ) = x2a (j) + vg2(j)*dt
dx2an(j-1) = x2an(j) - x2an(j-1)
10 continue
dx2an(je+2) = (dx2an(je+1)/dx2an(je))*dx2an(je+1)
c
dx2bn(js-2) = dx2an(js-2)
x2bn(js-2) = x2an(js-1) - 0.5*dx2an(js-2)
do 20 j=js-1,je+2
x2bn(j) = x2an(j) + 0.5*dx2an(j)
dx2bn(j) = x2bn(j) - x2bn(j-1)
20 continue
c
do 30 j=js-2,je+2
#ifdef XY
g32a h (j) = 1.0
g32b h (j) = 1.0
g4 a h (j) = 0.5*(x2a(j) + x2an(j))
g4 b h (j) = 0.5*(x2b(j) + x2bn(j))
g32a n (j) = 1.0
g32b n (j) = 1.0
g4 a n (j) = x2an(j)
g4 b n (j) = x2bn(j)
#endif
#ifdef RZ
g32a h (j) = 0.5*(x2a(j) + x2an(j))
g32b h (j) = 0.5*(x2b(j) + x2bn(j))
g4 a h (j) = 0.125*(x2a(j) + x2an(j))**2
g4 b h (j) = 0.125*(x2b(j) + x2bn(j))**2
g32a n (j) = x2an(j)
g32b n (j) = x2bn(j)
g4 a n (j) = 0.5*x2an(j)**2
g4 b n (j) = 0.5*x2bn(j)**2
#endif
#ifdef RT
g32a h (j) = sin(0.5*(x2a(j) + x2an(j)))
g32b h (j) = sin(0.5*(x2b(j) + x2bn(j)))
g4 a h (j) = -cos(0.5*(x2a(j) + x2an(j)))
g4 b h (j) = -cos(0.5*(x2b(j) + x2bn(j)))
g32a n (j) = sin(x2a n (j))
g32b n (j) = sin(x2b n (j))
g4 a n (j) = -cos(x2a n (j))
g4 b n (j) = -cos(x2b n (j))
#endif
30 continue
c
c New volume factors
c
do 40 j=js-2,je+1
dvl2ah(j ) = g4ah(j+1) - g4ah(j)
dvl2an(j ) = g4an(j+1) - g4an(j)
40 continue
do 50 j=js-2,je+1
dvl2bh(j+1) = g4bh(j+1) - g4bh(j)
dvl2bn(j+1) = g4bn(j+1) - g4bn(j)
50 continue
c
c New PPA coefficients
c
do 60 j=js,je+1
qa = dx2b(j+1) + dx2b(j)
qb = 1.0/(dx2b(j-1) + qa)
ppafc2(1,j) = qb/(dx2b(j+1) + dx2b(j))
ppafc2(2,j) = qb/(dx2b(j-1) + dx2b(j))
ppafc2(3,j) = (2.0*dx2b(j-1)+dx2b(j))*dx2b(j)*ppafc2(1,j)
ppafc2(4,j) = (2.0*dx2b(j+1)+dx2b(j))*dx2b(j)*ppafc2(2,j)
60 continue
do 70 j=js,je
qa = dx2b(j+1) + dx2b(j)
qc = dx2b(j+2) + dx2b(j-1) + qa
ppafc2(5,j) = dx2b(j )/qa
ppafc2(6,j) = dx2b(j+1)/qa
ppafc2(8,j) = dx2b(j )*(dx2b(j )+dx2b(j-1))
& /((2.0*dx2b(j)+dx2b(j+1))*qc)
ppafc2(7,j) = dx2b(j+1)*(dx2b(j+1)+dx2b(j+2))
& /((2.0*dx2b(j+1)+dx2b(j))*qc)
ppafc2(7,j) = 2.0/qa*(dx2b(j+1)*ppafc2(8,j)-dx2b(j)*ppafc2(7,j))
70 continue
do 80 j=js-1,je+1
qa = dx2a(j-1) + dx2a(j)
qb = 1.0/(dx2a(j+1) + qa)
ppazc2(1,j) = qb/(dx2a(j+1) + dx2a(j))
ppazc2(2,j) = qb/(dx2a(j-1) + dx2a(j))
ppazc2(3,j) = (2.0*dx2a(j-1)+dx2a(j))*dx2a(j)*ppazc2(1,j)
ppazc2(4,j) = (2.0*dx2a(j+1)+dx2a(j))*dx2a(j)*ppazc2(2,j)
80 continue
do 90 j=js,je+1
qa = dx2a(j-1) + dx2a(j)
qc = dx2a(j+1) + dx2a(j-2) + qa
ppazc2(7,j) = dx2a(j-1)*(dx2a(j-1)+dx2a(j-2))
& /((2.0*dx2a(j-1)+dx2a(j ))*qc)
ppazc2(8,j) = dx2a(j )*(dx2a(j )+dx2a(j+1))
& /((2.0*dx2a(j )+dx2a(j-1))*qc)
qd = 2.0/qa*(dx2a(j)*ppazc2(7,j)-dx2a(j-1)*ppazc2(8,j))
ppazc2(5,j) = dx2a(j-1)/qa + qd
ppazc2(6,j) = dx2a(j )/qa - qd
90 continue
return
end