-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathDEMSETTING.F
198 lines (137 loc) · 5.71 KB
/
DEMSETTING.F
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
c-----------------------------------------------------------------------------------c
c SUBROUTINE PROGRAM #
C DESINGED BY ChenJun #
C SHANGHAI JIAO TONG UNIVERSITY #
C SHANGHAI, CHINA #
c-----------------------------------------------------------------------------------#
c 子程序功能：
c 1------DEM变量初始化
c 2------颗粒初始信息读入（从文件）
c 3------DEM计算设置读入（从文件）
c-----------------------------------------------------------------------------------#
Subroutine DEMSETTING
Include './Include/OCERM_INF'
CALL DEMSETREADING
IF(I_TIMESAVE .EQ. 1)THEN
Open(1121,file='./Results/DEM-CPU-TIME-TOTAL.DAT')
ENDIF
c Write (1121,'(9X8A6)') ,' NNDEM',' TOTAL','IPPKPP','IPJPKP',
c & ' FCOUP',' BCDEM',' SOFT',' PMOVE','ACHIVE'
IF(I_TIMESAVE .EQ. 1)THEN
Open(1131,file='./Results/DEM-CPU-TIME-SELECT.DAT')
ENDIF
c Write (1131,'(9X8A6)') ,' NNDEM',' TOTAL','IPPKPP','IPJPKP',
c & ' FCOUP',' BCDEM',' SOFT',' PMOVE','ACHIVE'
!PRINT*,'INPUT NREAD:'
!READ(*,*)NREAD
I_READ_DEMSET = 0
I_READ_PARTICLE_SOFT = 0
NNDEM = 0
TTTDEM = 0.0
N_TEMP1 = 0
DO I = 1, N_PSLEV
IF(I .EQ. 1)THEN
NBX(I) = NBX1
NBY(I) = NBY1
NBZ(I) = NBZ1
DELBM(I) = DELBM1
IBCPMAX(I) = IBCPMAX1
C NBL1 = NBX1 * NBY1 * NBZ1
NBCP_INDEX(I) = 0
ELSEIF(I .EQ. 2)THEN
NBX(I) = NBX2
NBY(I) = NBY2
NBZ(I) = NBZ2
DELBM(I) = DELBM2
IBCPMAX(I) = IBCPMAX2
C NBL2 = NBX2 * NBY2 * NBZ2
NBCP_INDEX(I) = NBL1
ELSEIF(I .EQ. 3)THEN
NBX(I) = NBX3
NBY(I) = NBY3
NBZ(I) = NBZ3
DELBM(I) = DELBM3
IBCPMAX(I) = IBCPMAX3
C NBL3 = NBX3 * NBY3 * NBZ3
NBCP_INDEX(I) = NBL1 + NBL2
ELSEIF(I .EQ. 4)THEN
NBX(I) = NBX4
NBY(I) = NBY4
NBZ(I) = NBZ4
DELBM(I) = DELBM4
IBCPMAX(I) = IBCPMAX4
C NBL4 = NBX4 * NBY4 * NBZ4
NBCP_INDEX(I) = NBL1 + NBL2 + NBL3
ELSEIF(I .EQ. 5)THEN
NBX(I) = NBX5
NBY(I) = NBY5
NBZ(I) = NBZ5
DELBM(I) = DELBM5
IBCPMAX(I) = IBCPMAX5
C NBL5 = NBX5 * NBY5 * NBZ5
NBCP_INDEX(I) = NBL1 + NBL2 + NBL3 + NBL4
ELSE
PRINT*,'ERROR OF N_PSLEV IN OCERM_INF'
PAUSE
ENDIF
DO L = 1, NBZ(I)
DO K = 1, NBY(I)
DO J = 1, NBX(I)
N_TEMP = NBCP_INDEX(I) + (L-1)*NBX(I)*NBY(I) +
& NBX(I)*(K-1) + J !根据层级和所在层I,J,K确定网格全局编号
IBCP_INDEX(N_TEMP) = N_TEMP1 !确定当前网格颗粒存储起始位置
N_TEMP1 = N_TEMP1 + IBCPMAX(I) !当前网格颗粒存储终止位置，作为下一网格颗粒存储起始位置
ENDDO
ENDDO
ENDDO
ENDDO
!NBL = NBL1 + NBL2 + NBL3 + NBL4 + NBL5
CALL PARTICLE_INFO
C--- 床底单元边方向向量
DO I = 1, IJE
N1 = IEND_EDGE(I,1)
N2 = IEND_EDGE(I,2)
VECTOR_BE(I,1) = PXY(N2,1) - PXY(N1,1)
VECTOR_BE(I,2) = PXY(N2,2) - PXY(N1,2)
VECTOR_BE(I,3) = - HP(N2) + HP(N1)
VECTOR_BE(I,4) = VECTOR_BE(I,1)**2 + VECTOR_BE(I,2)**2 +
& VECTOR_BE(I,3)**2
ENDDO
C DO I = 1, IJM
C DO J = 1, CELL_POLYGEN(I)
C N1 = IEND_EDGE(I,1)
C N2 = IEND_EDGE(I,2)
C VECTOR_BE(I,1) = PXY(N2,1) - PXY(N1,1)
C VECTOR_BE(I,2) = PXY(N2,2) - PXY(N1,2)
C VECTOR_BE(I,3) = - HP(N2) + HP(N1)
C VECTOR_BE(I,4) = VECTOR_BE(I,1)**2 + VECTOR_BE(I,2)**2 +
C & VECTOR_BE(I,3)**2
C ENDDO
C--- 床底非与周围网格不共面的单元搜索
DO I = 1, IJM
NSYM_BC(I) = 0
IF(CELL_POLYGEN(I) .EQ. 3)THEN
NLS = 3
ELSEIF(CELL_POLYGEN(I) .EQ. 4)THEN
NLS = 2
ELSE
PRINT*,'ERROR OF CELL_POLYGEN(I) IN IPPKPPINITIAL.F'
PAUSE
ENDIF
DO J = 2, INL(I,NLS)
K = INE(I,J)
TEMP1 = SQRT(DEMCOSB(I,1)**2+DEMCOSB(I,2)**2+DEMCOSB(I,3)**2)
TEMP2 = SQRT(DEMCOSB(K,1)**2+DEMCOSB(K,2)**2+DEMCOSB(K,3)**2)
TEMP3 = DEMCOSB(I,1)*DEMCOSB(K,1)+DEMCOSB(I,2)*DEMCOSB(K,2)+
& DEMCOSB(I,3)*DEMCOSB(K,3)
TEMP4 = TEMP3 / (TEMP1 * TEMP2)
IF(ABS(ABS(TEMP4) - 1.0) .GT. 1.E-10)THEN
NSYM_BC(I) = 1
ENDIF
ENDDO
ENDDO
!DO I = 1, IJM
! IF(NSYM_BC(I) == 1)PRINT*,I
!ENDDO
!PRINT*,CXY(1198,1),CXY(1198,2)
End