-
Notifications
You must be signed in to change notification settings - Fork 0
/
PROFV.f90
169 lines (141 loc) · 3.98 KB
/
PROFV.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
!****************************************************************************
!
! Subroutine program: PROFV
! Date: 2024.09.30
! Detail: Velocity solver by momentum equation (U,V).
!
!****************************************************************************
Subroutine PROFV
Use NUM_VAR
Implicit None
Integer I,J,K,L,M,N,K2,LIM,IL,IR
Real(Kind=8) AAAA(KBM,KBM),AA(3*KBM-2),BBBB(KBM),BB(KBM)
Real(Kind=8) GRAD2X(IJM),GRAD2Y(IJM),BU,BV,UD
LIM = 1
!============================================================================
! Coefficient matrix of the equations
!============================================================================
Do K = 1,KBM
Do K2 = 1,KBM
AAAA(K,K2) = 0
Enddo
Enddo
Do K = 1,3*KBM-2
AA(K) = 0
Enddo
Do K = 1,KBM
BB(K) = 0
BBBB(K) = 0
Enddo
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(I,J,K,AAAA,BBBB,AA,BB,BU,BV,IL,IR,UD)
!$OMP DO
Do I = 1,IJM
If (CCM(I) == 1) Then
UAVE(I) = 0
VAVE(I) = 0
Do K = 1,KBM
UAVE(I) = UAVE(I) + U(I,K) * DZ(K)
VAVE(I) = VAVE(I) + V(I,K) * DZ(K)
Enddo
Endif
Enddo
!$OMP END DO
!$OMP DO
Do I = 1,IJE
If (CFM(I) == 1) Then
IL = INDEX_EDGE(I,1,1)
IR = INDEX_EDGE(I,1,2)
UD = 0.5 * (UAVE(IL) + UAVE(IR)) * (CXY(IR,1) - CXY(IL,1)) + &
0.5 * (VAVE(IL) + VAVE(IR)) * (CXY(IR,2) - CXY(IL,2))
If (UD < 0.0) Then
INDEX_EDGE(I,1,1) = IR
INDEX_EDGE(I,1,2) = IL
Endif
Endif
Enddo
!$OMP END DO
!$OMP DO
Do I = 1,IJM
GRAD2X(I) = 0
GRAD2Y(I) = 0
If (CCM(I) == 1) Then
Do J = 1,CELL_POLYGEN(I)
If (CFM(CELL_SIDE(I,J,1)) == 1) Then
GRAD2X(I) = GRAD2X(I) + WIX(I,J) * (ELF(CELL_SIDE(I,J,2)) - ELF(I))
GRAD2Y(I) = GRAD2Y(I) + WIY(I,J) * (ELF(CELL_SIDE(I,J,2)) - ELF(I))
Endif
Enddo
Endif
Enddo
!$OMP END DO
!$OMP DO
Do I = 1,IJM
IF (CCM(I) == 1) Then
!----------------------------------------------------------------------------
Do K = 2,KBM
AAAA(K,K-1) = -1 * DTI * ((KM(I,K-1) + KM(I,K)) / 2 + UMOL) / &
DC(I)**2 / DZZ(K-1)
AAAA(K-1,K) = AAAA(K,K-1)
Enddo
Do K = 2,KBM-1
AAAA(K,K) = DZ(K) - AAAA(K,K-1) - AAAA(K,K+1)
Enddo
AAAA(1,1) = DZ(1) - AAAA(1,2)
AAAA(KBM,KBM) = DZ(KBM) - AAAA(KBM,KBM-1) + DTI * CBC(I) * &
Sqrt(U(I,KBM)**2 + V(I,KBM)**2) / DC(I)**2
Do K = 1,KBM
BBBB(K) = GRAV * DC(I) * DTI * THETA * DZ(K)
Enddo
BU = GRAD2X(I)
BV = GRAD2Y(I)
!============================================================================
! Velocity solver based on forward elimination and back substitution
!============================================================================
AA(1) = AAAA(1,1)
AA(2) = AAAA(1,2)
Do K = 2,KBM-1
AA(2*(K-1)+K-1) = AAAA(K,K-1)
AA(2*(K-1)+K) = AAAA(K,K)
AA(2*(K-1)+K+1) = AAAA(K,K+1)
Enddo
AA(3*KBM-3) = AAAA(KBM,KBM-1)
AA(3*KBM-2) = AAAA(KBM,KBM)
Do K = 1,KBM
BB(K) = -1 * BBBB(K) * BU
Enddo
Call MAT_TRD(AA,KBM,3*KBM-2,BB,LIM)
Do K = 1,KBM
U(I,K) = BB(K)
Enddo
AA(1) = AAAA(1,1)
AA(2) = AAAA(1,2)
Do K = 2,KBM-1
AA(2*(K-1)+K-1) = AAAA(K,K-1)
AA(2*(K-1)+K) = AAAA(K,K)
AA(2*(K-1)+K+1) = AAAA(K,K+1)
Enddo
AA(3*KBM-3) = AAAA(KBM,KBM-1)
AA(3*KBM-2) = AAAA(KBM,KBM)
Do K = 1,KBM
BB(K) = -1 * BBBB(K) * BV
Enddo
Call MAT_TRD(AA,KBM,3*KBM-2,BB,LIM)
Do K = 1,KBM
V(I,K) = BB(K)
Enddo
!----------------------------------------------------------------------------
Endif
Enddo
!$OMP END DO
!$OMP DO
Do I = 1,IJM
If (CCM(I) == 1) Then
Do K = 1,KBM
U(I,K) = U(I,K) + USTAR(I,K)
V(I,K) = V(I,K) + VSTAR(I,K)
Enddo
Endif
Enddo
!$OMP END DO
!$OMP END PARALLEL
End Subroutine PROFV