-
Notifications
You must be signed in to change notification settings - Fork 0
/
momx2.src
115 lines (115 loc) · 3.61 KB
/
momx2.src
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
#include "zeus2d.def"
c=======================================================================
c///////////////////////// SUBROUTINE MOMX2 \\\\\\\\\\\\\\\\\\\\\\\\\\
c
subroutine momx2(mflx2,s1,s2)
c
c PURPOSE: Transports the 1- and 2-components of vector variables
c in the 2-direction. Currently transported are:
c 1- and 2-components of momentum density
c The consistent transport algorithm is used, including the effects of
c grid compression. The transported fluxes are thus given by the mass
c fluxes times the time centered area of the control volume faces times
c the interpolated velocities.
c Interpolations are performed in X2INTFC and X2INTZC.
c
c INPUT ARGUMENTS:
c mflx2 = mass flux in 2-direction
c s1 = momentum density in 1-direction
c s2 = momentum density in 2-direction
c
c OUTPUT ARGUMENTS:
c s1 = "half"-updated momentum density in 1-direction
c s2 = "half"-updated momentum density in 2-direction
c
c EXTERNALS: X2INTFC
c X2INTZC
c
c LOCALS:
c-----------------------------------------------------------------------
implicit NONE
#include "param.h"
#include "root.h"
#include "grid.h"
#include "field.h"
#include "scratch.h"
REAL mflx2(in,jn),s1(in,jn),s2(in,jn)
c
integer i,j
REAL dflx2
REAL atwid(in),vel2(jn),tv1(jn),v1twid2(jn),s1flx2(jn)
. ,tv2(jn),v2twid2(jn),s2flx2(jn)
equivalence (atwid,wi0),(vel2,wj1),(v1twid2,wj2)
. ,(tv1,s1flx2,wj3),(v2twid2,wj4),(tv2,s2flx2,wj5)
c
external x2intfc,x2intzc
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////////
c=======================================================================
c Check for 1-D problem in the 1-direction
c
if (nx2z .le. 1) return
c================= TRANSPORT 1-COMPONENT IN 2-DIRECTION ==============
c Compute time-centered are factors
c
do 10 i=is+1,ie
atwid(i) = g31a(i)*dx1b(i)/dvl1b(i)
10 continue
c
c Interpolate v1 to interfaces in the 2-direction. vel2 is
c the relative fluid velocity at interpolation points.
c
do 100 i=is+1,ie
do 20 j=ji(i),jop1(i)
vel2(j) = 0.5*((v2(i-1,j) - vg2(j)) + (v2(i,j) - vg2(j)))
20 continue
do 30 j=jim2(i),jop2(i)
tv1(j) = v1(i,j)
30 continue
call x2intzc(tv1 ,vel2,wj0,i,g2a,iords1 ,0,v1twid2 )
c
do 40 j=ji(i),jop1(i)
dflx2 = 0.5*(mflx2(i-1,j) + mflx2(i,j))*atwid(i)*g32ah(j)
s1flx2(j) = dflx2*v1twid2(j)
40 continue
c
c Perform advection using fluxes. Note timestep dt is hidden in the
c fluxes.
c
do 50 j=ji(i),jo(i)
s1(i,j) = (s1(i,j)*dvl2a(j)-(s1flx2(j+1)-s1flx2(j)))/dvl2an(j)
50 continue
100 continue
c
c=============== TRANSPORT 2-COMPONENTS IN 2-DIRECTION ===============
c Compute time centered area factor
c
do 110 i=is,ie
atwid(i) = g31b(i)*dx1a(i)/dvl1a(i)
110 continue
c
c Interpolate v2 to zone centers in the 2-direction.
c
do 200 i=is,ie
do 120 j=ji(i),jo(i)
vel2(j) = 0.5*(v2(i,j) - vg2(j) + v2(i,j+1) - vg2(j+1))
120 continue
do 130 j=jim1(i),jop2(i)
tv2(j) = v2(i,j)
130 continue
call x2intfc(tv2 ,vel2,i,g2b,iords2 ,istps2 ,v2twid2 )
c
do 140 j=ji(i),jo(i)
dflx2 = 0.5*(mflx2(i,j) + mflx2(i,j+1))*atwid(i)*g32bh(j)
s2flx2(j) = dflx2*v2twid2(j)*g2b(i)
140 continue
c
c Perform advection using fluxes. Note timestep dt is hidden in the
c fluxes.
c
do 150 j=jip1(i),jo(i)
s2(i,j) = (s2(i,j)*dvl2b(j)-(s2flx2(j)-s2flx2(j-1)))/dvl2bn(j)
150 continue
200 continue
c
return
end