-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsediment_source.f90
275 lines (219 loc) · 9.67 KB
/
sediment_source.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
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
! ----------------------------------------------------------------
! file: sediment_source.f90
! ----------------------------------------------------------------
! ----------------------------------------------------------------
! Battelle Memorial Institute
! Pacific Northwest Laboratory
! ----------------------------------------------------------------
! ----------------------------------------------------------------
! Created August 23, 2000 by William A. Perkins
! Last Change: 2014-06-09 15:49:14 d3g096
! ----------------------------------------------------------------
! ----------------------------------------------------------------
! MODULE sediment_source
! ----------------------------------------------------------------
MODULE sediment_source
USE utility
USE config
USE block_module
USE scalars, ONLY: max_species
IMPLICIT NONE
CHARACTER (LEN=80), PRIVATE, SAVE :: rcsid = "$Id$"
TYPE sediment_source_block_rec
! these are rates, units: (mass)/ft^2/s
TYPE (block_var), POINTER :: bv_deposition
DOUBLE PRECISION, POINTER :: deposition(:,:)
TYPE (block_var), POINTER :: bv_erosion
DOUBLE PRECISION, POINTER :: erosion(:,:)
END TYPE sediment_source_block_rec
TYPE sediment_source_rec
INTEGER :: ifract ! sediment fraction index
DOUBLE PRECISION :: erode, eshear
DOUBLE PRECISION :: setvel, dshear
DOUBLE PRECISION :: pdens, d50
TYPE (sediment_source_block_rec), POINTER :: block(:)
END TYPE sediment_source_rec
! total number of sediment fractions
! simulated (may be needed by bed more
! than here)
INTEGER, PUBLIC, SAVE :: sediment_fractions = 0
! map sediment fraction index to
! scalar index
INTEGER, PUBLIC, ALLOCATABLE :: sediment_scalar_index(:)
CONTAINS
! ----------------------------------------------------------------
! SUBROUTINE sediment_source_initialize
! ----------------------------------------------------------------
SUBROUTINE sediment_source_initialize()
IMPLICIT NONE
INTEGER :: iblk
ALLOCATE(sediment_scalar_index(sediment_fractions))
END SUBROUTINE sediment_source_initialize
! ----------------------------------------------------------------
! TYPE (SEDIMENT_SOURCE_REC) FUNCTION sediment_parse_options
! ----------------------------------------------------------------
TYPE (SEDIMENT_SOURCE_REC) FUNCTION sediment_parse_options(options)
IMPLICIT NONE
POINTER sediment_parse_options
CHARACTER (LEN=*) :: options(:)
CHARACTER (LEN=1024) :: msg
INTEGER :: nopt
INTEGER :: i, iblk
i = 1
nopt = UBOUND(options, 1)
ALLOCATE(sediment_parse_options)
sediment_parse_options%erode = 0.0
sediment_parse_options%eshear = -1.0
sediment_parse_options%setvel = 0.0
sediment_parse_options%dshear = 0.0
sediment_parse_options%pdens = 0.0
sediment_parse_options%d50 = 0.0
! allocate storage for erosion and
! deposition rates
ALLOCATE(sediment_parse_options%block(max_blocks))
DO iblk = 1, max_blocks
sediment_parse_options%block(iblk)%bv_deposition => &
&block_var_allocate("deposition", block(iblk)%varbase, .TRUE.)
sediment_parse_options%block(iblk)%deposition => &
&sediment_parse_options%block(iblk)%bv_deposition%current
sediment_parse_options%block(iblk)%bv_erosion => &
&block_var_allocate("erosion", block(iblk)%varbase, .TRUE.)
sediment_parse_options%block(iblk)%erosion => &
&sediment_parse_options%block(iblk)%bv_erosion%current
sediment_parse_options%block(iblk)%deposition = 0
CALL block_var_put(sediment_parse_options%block(iblk)%bv_deposition)
sediment_parse_options%block(iblk)%erosion = 0
CALL block_var_put(sediment_parse_options%block(iblk)%bv_erosion)
END DO
DO WHILE ((LEN_TRIM(options(i)) .GT. 0) .AND. (i .LE. nopt))
SELECT CASE (options(i))
! median partical diameter, feet
CASE ('D50')
IF ((i + 1 .GT. nopt) .OR. (LEN_TRIM(options(i+1)) .LE. 0)) THEN
WRITE(msg, 100) 'D50'
CALL error_message(msg, fatal=.TRUE.)
END IF
READ(options(i+1), *) sediment_parse_options%d50
i = i + 1
! erodibility coefficient, (mass)/ft^2/s
CASE ('ERODIBILTY')
IF ((i + 1 .GT. nopt) .OR. (LEN_TRIM(options(i+1)) .LE. 0)) THEN
WRITE(msg, 100) 'ERODIBILTY'
CALL error_message(msg, fatal=.TRUE.)
END IF
READ(options(i+1), *) sediment_parse_options%erode
i = i + 1
! critical shear stress for erosion, lbf/ft^2
CASE ('ESHEAR')
IF ((i + 1 .GT. nopt) .OR. (LEN_TRIM(options(i+1)) .LE. 0)) THEN
WRITE(msg, 100) 'ESHEAR'
CALL error_message(msg, fatal=.TRUE.)
END IF
READ(options(i+1), *) sediment_parse_options%eshear
i = i + 1
! critical shear stress for deposition, lbf/ft^2
CASE ('DSHEAR')
IF ((i + 1 .GT. nopt) .OR. (LEN_TRIM(options(i+1)) .LE. 0)) THEN
WRITE(msg, 100) 'DSHEAR'
CALL error_message(msg, fatal=.TRUE.)
END IF
READ(options(i+1), *) sediment_parse_options%dshear
i = i + 1
! settling velocity for deposition, ft/s
CASE ('SETVEL')
IF ((i + 1 .GT. nopt) .OR. (LEN_TRIM(options(i+1)) .LE. 0)) THEN
WRITE(msg, 100) 'SETVEL'
CALL error_message(msg, fatal=.TRUE.)
END IF
READ(options(i+1), *) sediment_parse_options%setvel
i = i + 1
! solids density, mass/ft^3
CASE ('DENSITY')
IF ((i + 1 .GT. nopt) .OR. (LEN_TRIM(options(i+1)) .LE. 0)) THEN
WRITE(msg, 100) 'SETVEL'
CALL error_message(msg, fatal=.TRUE.)
END IF
READ(options(i+1), *) sediment_parse_options%pdens
i = i + 1
CASE DEFAULT
WRITE(msg, *) 'GEN scalar option "', &
&TRIM(options(i)), '" not understood and ignored'
CALL error_message(msg, fatal=.TRUE.)
END SELECT
i = i + 1
END DO
! in order for a sediment species to
! be meaningful, all of the parameters
! should be non zero
IF (sediment_parse_options%setvel .LE. 0.0) THEN
WRITE (msg, *) 'Invalid or unspecified value for SETVEL'
CALL error_message(msg, fatal=.TRUE.)
END IF
IF (sediment_parse_options%eshear .LE. 0.0) THEN
WRITE (msg, *) 'Invalid or unspecified value for ESHEAR'
CALL error_message(msg, fatal=.TRUE.)
END IF
IF (sediment_parse_options%dshear .LE. 0.0) THEN
WRITE (msg, *) 'Invalid or unspecified value for DSHEAR'
CALL error_message(msg, fatal=.TRUE.)
END IF
IF (sediment_parse_options%erode .LT. 0.0) THEN
WRITE (msg, *) 'Invalid or unspecified value for ERODIBILTY'
CALL error_message(msg, fatal=.TRUE.)
END IF
IF (sediment_parse_options%pdens .LE. 0.0) THEN
WRITE (msg, *) 'Invalid or unspecified value for DENSITY'
CALL error_message(msg, fatal=.TRUE.)
END IF
IF (sediment_parse_options%d50 .LE. 0.0) THEN
WRITE (msg,*) 'Bad D50 value for SED species: ', sediment_parse_options%d50
CALL error_message(msg, fatal=.TRUE.)
END IF
100 FORMAT('additional argument missing for ', A10, ' keyword')
END FUNCTION sediment_parse_options
! ----------------------------------------------------------------
! DOUBLE PRECISION FUNCTION sediment_erosion
! ----------------------------------------------------------------
DOUBLE PRECISION FUNCTION sediment_erosion(rec, iblk, i, j)
IMPLICIT NONE
INCLUDE 'bed_functions.inc'
TYPE (sediment_source_rec) :: rec
INTEGER, INTENT(IN) :: iblk, i, j
DOUBLE PRECISION :: shear
sediment_erosion = 0.0
shear = block(iblk)%shear(i,j)
IF (shear > rec%eshear) THEN
sediment_erosion = &
&MIN((shear/rec%eshear - 1.0)*rec%erode,&
&bed_max_erosion(rec%ifract, iblk, i, j, delta_t))
END IF
END FUNCTION sediment_erosion
! ----------------------------------------------------------------
! DOUBLE PRECISION FUNCTION sediment_deposition
! ----------------------------------------------------------------
DOUBLE PRECISION FUNCTION sediment_deposition(rec, iblk, i, j, sconc)
IMPLICIT NONE
TYPE (sediment_source_rec) :: rec
INTEGER, INTENT(IN) :: iblk, i, j
DOUBLE PRECISION, INTENT(IN) :: sconc
DOUBLE PRECISION :: shear
sediment_deposition = 0.0
shear = block(iblk)%shear(i,j)
IF (shear < rec%dshear) THEN
sediment_deposition = (1.0 - shear/rec%dshear)*rec%setvel*sconc
sediment_deposition = MAX(sediment_deposition, 0.0d0)
END IF
END FUNCTION sediment_deposition
! ----------------------------------------------------------------
! DOUBLE PRECISION FUNCTION sediment_source_term
! ----------------------------------------------------------------
DOUBLE PRECISION FUNCTION sediment_source_term(rec, iblk, i, j, sconc)
IMPLICIT NONE
TYPE(sediment_source_rec), INTENT(IN) :: rec
INTEGER, INTENT(IN) :: iblk, i, j
DOUBLE PRECISION :: sconc
sediment_source_term = &
&rec%block(iblk)%erosion(i, j) - &
&rec%block(iblk)%deposition(i, j)
END FUNCTION sediment_source_term
END MODULE sediment_source