-
Notifications
You must be signed in to change notification settings - Fork 0
/
advect.f90
200 lines (162 loc) · 6.73 KB
/
advect.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
subroutine advect(var,uvarx)
! ---------------------------------------------
USE header
! convecn(m,n,dtime) means use level m and write s,T to level n
! computes Cx,Cy,Cz at the cell centers.
! uflx,vflx,wflx,sflx,Tflx are the fluxes defined at
! cell faces (using QUICK).
!
INTEGER :: i,i2,j,k,n,m,step
REAL(kind=rc_kind), PARAMETER :: Eighth=0.125d0,Half=0.5
REAL(kind=rc_kind), PARAMETER :: courinv=2.d0
REAL(kind=rc_kind) :: absvel,left,right,ctr,fc,upos,uneg
REAL(kind=rc_kind) :: dtimel
REAL(kind=rc_kind), dimension( 0:NI+1,0:NJ+1, 0:NK+1) :: var ! var=(s,T,u,v,w,Tr(it,:,:,:))
REAL(kind=rc_kind), dimension( 0:NI+1,0:NJ+1, 0:NK+1) :: uvarx ! uvarx is the source term, divergence of the advective fluxes
REAL(kind=rc_kind), dimension( 0:NI ,0:NJ , 0:NK ) :: varflx ! varflx is a temporary variable to compute fluxes at the faces
REAL(kind=rc_kind) :: dvarx(0:NI+1),dvary(0:NJ),dvarz(0:NK+1) ! dvarX is d(var)/dX at the faces (linear difference).
!------------------------------------------
! computation of the source term in x.
do k=1,NK
do j=1,NJ
dvarx(0)=var(1,j,k)-var(NI,j,k)
do i=1,NI-1
dvarx(i)= var(i+1,j,k) - var(i,j,k)
end do
dvarx(NI)=dvarx(0)
dvarx(NI+1)=dvarx(1)
var(NI+1,j,k)= var(1,j,k)
var(0,j,k)= var(NI,j,k)
do i=1,NI
if (uf(i,j,k).gt.0.d0) then
left= Eighth*(dvarx(i) -dvarx(i-1))
ctr= Half* (var(i+1,j,k) +var(i,j,k))
fc= ctr -left
call ultim(var(i-1,j,k),var(i,j,k),var(i+1,j,k),dvarx(i),dvarx(i-1),courinv,fc)
varflx(i,j,k)= uf(i,j,k)*fc
else
if(i==NI) then
i2=0;
else
i2=i;
endif
right= Eighth*(dvarx(i2+1) -dvarx(i2))
ctr= Half* (var(i2+1,j,k) +var(i,j,k))
fc= ctr - right
call ultim(var(i2+2,j,k),var(i2+1,j,k),var(i,j,k),dvarx(i),dvarx(i2+1),courinv,fc)
varflx(i,j,k)= uf(i,j,k)*fc
end if
end do
varflx(0,j,k)= varflx(NI,j,k)
end do
end do
do i=1,NI
uvarx(i,1:NJ,1:NK)=varflx(i,1:NJ,1:NK) -varflx(i-1,1:NJ,1:NK)
end do
! computation of the source term in x done.
!------------------------------------------
!------------------------------------------
! computation of the source term in y.
do k=1,NK
do i=1,NI
dvary(0)= 0.d0;
dvary(NJ)= 0.d0;
do j=1,NJ-1
dvary(j)= var(i,j+1,k) - var(i,j,k)
end do
var(i,0,k)= var(i,1,k)
var(i,NJ+1,k)= var(i,NJ,k)
do j=1,NJ-1
if (vf(i,j,k).ge.0.d0) then
left= Eighth*(dvary(j) -dvary(j-1))
ctr= Half* (var(i,j+1,k) +var(i,j,k))
fc= ctr -left
call ultim(var(i,j-1,k),var(i,j,k),var(i,j+1,k),dvary(j),dvary(j-1),courinv,fc)
varflx(i,j,k)= vf(i,j,k)*fc
else
right= Eighth*(dvary(j+1) -dvary(j))
ctr= Half* (var(i,j+1,k) +var(i,j,k))
fc= ctr - right
call ultim(var(i,j+2,k),var(i,j+1,k),var(i,j,k),dvary(j),dvary(j+1),courinv,fc)
varflx(i,j,k)= vf(i,j,k)*fc
end if
end do
varflx(i,0,k)=0.d0;
varflx(i,NJ,k)=0.d0;
end do
end do
do j=1,NJ
uvarx(1:NI,j,1:NK)= (varflx(1:NI,j,1:NK) -varflx(1:NI,j-1,1:NK) ) + uvarx(1:NI,j,1:NK)
end do
! computation of the source term in y done.
!------------------------------------------
!------------------------------------------
! computation of the source term in z.
do j=1,NJ
do i=1,NI
var(i,j,0)= var(i,j,1)
dvarz(0)= 0.d0
do k=1,NK-1
dvarz(k)= var(i,j,k+1) - var(i,j,k)
end do
dvarz(NK:NK+1)= 0.d0 ! Top boundary - Zero gradient
var(i,j,NK+1)= var(i,j,NK)
do k=1,NK
if (wf(i,j,k).ge.0.d0) then
left= Eighth*(dvarz(k) -dvarz(k-1))
ctr= Half* (var(i,j,k+1) +var(i,j,k))
fc= ctr -left
call ultim(var(i,j,k-1),var(i,j,k),var(i,j,k+1),dvarz(k),dvarz(k-1),courinv,fc)
varflx(i,j,k)= wf(i,j,k)*fc
else
if (k.eq.NK) then
varflx(i,j,k)= wf(i,j,k)*Half*(var(i,j,k+1) +var(i,j,k))
else
right= Eighth*(dvarz(k+1) -dvarz(k))
ctr= Half* (var(i,j,k+1) +var(i,j,k))
fc= ctr - right
call ultim(var(i,j,k+2),var(i,j,k+1),var(i,j,k),dvarz(k),dvarz(k+1),courinv,fc)
varflx(i,j,k)= wf(i,j,k)*fc
end if
end if
end do
varflx(i,j,0)= 0.d0 ! Solid Boundary
end do
end do
do k=1,NK
uvarx(1:NI,1:NJ,k)= ( varflx(1:NI,1:NJ,k) -varflx(1:NI,1:NJ,k-1) ) + uvarx(1:NI,1:NJ,k)
enddo
! computation of the source term in z done.
!------------------------------------------
return
END
!----------------------------------------------------------------------
! SUBROUTINE ultim
! Implementation of the ULTIMATE limiter based on upstrm,dnstrm,ctr
! points and d(phi)/dx at the up- and dnstrm locations
! Returns the "ulitimate limited" face value of the variable
! that should be multiplied by uf,vf or wf to give the flux.
SUBROUTINE ultim(up,ct,dn,gdn,gup,courinv,fc)
USE header, only : rc_kind
REAL(kind=rc_kind) :: up,ct,dn,gup,gdn,courinv,fc
REAL(kind=rc_kind) :: ref,del,adel,acurv
del= dn -up
adel= abs(del)
acurv= abs(gdn -gup)
if (acurv.ge.adel) then
fc= ct
return
else
ref= up + (ct -up)*courinv
if (del.gt.0) then
fc= max(ct,fc)
ref= min(ref,dn)
fc= min(fc,ref)
else
fc= min(ct,fc)
ref= max(ref,dn)
fc= max(ref,fc)
endif
end if
return
END