forked from csiro-coasts/EMS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
trvdiffsettl.c
451 lines (400 loc) · 13.6 KB
/
trvdiffsettl.c
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
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
/*
*
* ENVIRONMENTAL MODELLING SUITE (EMS)
*
* File: model/lib/sediment/trvdiffsettl.c
*
* Description:
* Calculate vertical diffusion and settling of tracer
*
* Copyright:
* Copyright (c) 2018. Commonwealth Scientific and Industrial
* Research Organisation (CSIRO). ABN 41 687 119 230. All rights
* reserved. See the license file for disclaimer and full
* use/redistribution conditions.
*
* $Id: trvdiffsettl.c 5848 2018-06-29 05:01:15Z riz008 $
*
*/
#ifdef __cplusplus
#include "stdafx.h"
extern "C" {
#endif
#include "sediments.h"
static int implicit_vadv_vdiff(double dt, int nz, double *dz, double *dzold,
double *por, double *porold, double *var,
double *dvar, double *Kz, double *w,
double botinflux, double topinflux,
double botoutfluxcoef,
double topoutfluxcoef, int kb, int kt,
double *Splus, double *Sminus, double *dzface);
void vdiff_sedim_wc(sediment_t *sediment, sed_column_t *sm,
sed_tracer_t *tracer,
sed_tracer_t *tracercar, double botinflux,
double topinflux, double botoutfluxcoef,
double topoutfluxcoef, double *dvar)
{
sed_params_t *param = sediment->msparam;
long k, kt, kb;
double top, bot;
double *Kzij; /* diffusion coefficients for 1 column */
double *w;
double *var = sm->tr_wc[tracer->n];
double s = 1.; /* Kz scale factor */
double *dummy;
int ncar = tracercar->n;
/* Allocate temporary arrays */
Kzij = d_alloc_1d(param->nz + 1);
w = d_alloc_1d(param->nz + 1);
dummy = d_alloc_1d(param->nz + 1);
kb = sm->botk_wc;
kt = sm->topk_wc;
top = sm->topz_wc;
bot = sm->botz_wc;
/* nothing to do if dry */
/* if( top <= bot ) */
/* continue; */
if ( tracer->diffuse == 0 )
s=0;
/* Get diffusion coefficients for this column */
Kzij[kt + 1] = 0;
Kzij[kt] = s * sm->Kz_wc[kt];
for (k = kt - 1; k > kb; k--)
Kzij[k] = s * sm->Kz_wc[k];
Kzij[kb] = 0.;
/* sediment velocity */
for (k = kb; k <= kt + 1; k++)
w[k] = -sm->gridvel_wc[k] + sm->svel_wc[ncar][k];
for (k = kb; k <= kt; k++)
dummy[k] = 1.;
/* MH 07/2012: included instability handling */
if (implicit_vadv_vdiff(param->dt, param->nz, sm->dz_wc, sm->dzold_wc,
dummy, dummy, var, dvar, Kzij, w,
botinflux, topinflux, botoutfluxcoef, topoutfluxcoef,
kb, kt, NULL, NULL,sm->dzface_wc))
i_set_error(sediment->hmodel, sm->col_number, LFATAL, "trvdiffsettl:vdiff_sedim_wc: Error encountered in implicit_vadv_vdiff().\n");
/* END MH */
/* free temporary arrays */
d_free_1d(Kzij);
d_free_1d(w);
d_free_1d(dummy);
}
/******************************************************************************/
void vdiff_dissolved_wc(sediment_t *sediment, sed_column_t *sm, sed_tracer_t *tracer,
double botinflux, double topinflux,
double botoutfluxcoef, double topoutfluxcoef,
double *dvar)
{
sed_params_t *param = sediment->msparam;
long k, kt, kb;
double top, bot;
double *Kzij; /* diffusion coefficients for 1 column */
double *w;
double *var = sm->tr_wc[tracer->n];
double s = 1.; /* Kz scale factor */
double *dummy;
/* Allocate temporary arrays */
Kzij = d_alloc_1d(param->nz + 1);
w = d_alloc_1d(param->nz + 1);
dummy = d_alloc_1d(param->nz + 1);
kb = sm->botk_wc;
kt = sm->topk_wc;
top = sm->topz_wc;
bot = sm->botz_wc;
/* nothing to do if dry */
/* if( top <= bot ) */
/* continue; */
if ( tracer->diffuse == 0 )
s=0;
/* Get diffusion coefficients for this column */
Kzij[kt + 1] = 0;
Kzij[kt] = s * sm->Kz_wc[kt];
for (k = kt - 1; k > kb; k--)
Kzij[k] = s * sm->Kz_wc[k];
Kzij[kb] = 0.;
for (k = kb; k <= kt + 1; k++)
w[k] = sm->watvel_wc[k];
for (k = kb; k <= kt; k++)
dummy[k] = 1.;
/* MH 07/2012: included instability handling */
if (implicit_vadv_vdiff(param->dt, param->nz, sm->dz_wc, sm->dzold_wc,
sm->por_wc, sm->porold_wc, var, dvar, Kzij, w,
botinflux, topinflux, botoutfluxcoef, topoutfluxcoef,
kb, kt, NULL, NULL,sm->dzface_wc))
i_set_error(sediment->hmodel, sm->col_number, LFATAL, "trvdiffsettl:vdiff_dissolved_wc: Error encountered in implicit_vadv_vdiff().\n");
/* END MH */
/* free temporary arrays */
d_free_1d(Kzij);
d_free_1d(w);
d_free_1d(dummy);
}
/******************************************************************************/
void vdiff_sedim_sed(sediment_t *sediment, sed_column_t *sm, sed_tracer_t *tracer,
double botinflux, double topinflux,
double botoutfluxcoef, double topoutfluxcoef,
double *dvar)
{
sed_params_t *param = sediment->msparam;
long k, kt, kb;
double top, bot;
double *Kzij; /* diffusion coefficients for 1 column */
double *var;
double *w;
double s = 1.; /* Kz scale factor */
int n = tracer->n;
double *dummy;
/* Allocate temporary arrays */
Kzij = d_alloc_1d(param->sednz + 1);
w = d_alloc_1d(param->sednz + 1);
var = d_alloc_1d(param->sednz + 1);
dummy = d_alloc_1d(param->sednz + 1);
kb = sm->botk_sed;
kt = sm->topk_sed;
top = sm->topz_sed;
bot = sm->botz_sed;
/* nothing to do if dry */
/* if( top <= bot ) */
/* continue; */
/* Get diffusion coefficients and vertical velocities for this column */
Kzij[kt + 1] = 0;
Kzij[kt] = s * sm->partic_kz[kt];
for (k = kt - 1; k > kb; k--)
Kzij[k] = s * sm->partic_kz[k];
Kzij[kb] = 0.;
/* relative velocity of bottom sediment */
w[kt + 1] = 0.;
for (k = kb; k <= kt; k++)
w[k] = -sm->gridvel_sed[k] + sm->svel_consolid[k];
w[kb] = 0.;
for (k = kb; k <= kt; k++)
var[k] = sm->tr_sed[n][k];
for (k = kb; k <= kt; k++)
dummy[k] = 1.;
/* MH 07/2012: included instability handling */
if (implicit_vadv_vdiff(param->dt, param->sednz, sm->dz_sed, sm->dzold_sed,
dummy, dummy, var, dvar, Kzij, w,
botinflux, topinflux, botoutfluxcoef, topoutfluxcoef,
kb, kt, NULL, NULL,sm->dzface_sed))
i_set_error(sediment->hmodel, sm->col_number, LFATAL, "trvdiffsettl:vdiff_sedim_sed: Error encountered in implicit_vadv_vdiff().\n");
/* END MH */
/* free temporary arrays */
d_free_1d(Kzij);
d_free_1d(w);
d_free_1d(var);
d_free_1d(dummy);
}
/******************************************************************************/
void vdiff_dissolved_sed(sediment_t *sediment, sed_column_t *sm, sed_tracer_t *tracer,
double botinflux, double topinflux,
double botoutfluxcoef, double topoutfluxcoef,
double *dvar)
{
sed_params_t *param = sediment->msparam;
long k, kt, kb;
double top, bot;
double *Kzij; /* diffusion coefficients for 1 column */
double *var;
double *w;
double s = 1.; /* Kz scale factor */
int n = tracer->n;
/* Allocate temporary arrays */
Kzij = d_alloc_1d(param->sednz + 1);
w = d_alloc_1d(param->sednz + 1);
var = d_alloc_1d(param->sednz + 1);
kb = sm->botk_sed;
kt = sm->topk_sed;
top = sm->topz_sed;
bot = sm->botz_sed;
/* nothing to do if dry */
/* if( top <= bot ) */
/* continue; */
/* Get diffusion coefficients and vertical velocities for this column */
Kzij[kt + 1] = 0;
Kzij[kt] = s * sm->dissol_kz[kt];
for (k = kt - 1; k > kb; k--)
Kzij[k] = s * sm->dissol_kz[k];
Kzij[kb] = 0.;
for (k = kb; k <= kt; k++)
var[k] = sm->tr_sed[n][k];
/* relative velocity of bed water */
for (k = kb; k <= kt + 1; k++)
w[k] = sm->watvel_sed[k];
/* MH 07/2012: included instability handling */
if (implicit_vadv_vdiff(param->dt, param->sednz, sm->dz_sed, sm->dzold_sed,
sm->por_sed, sm->porold_sed, var, dvar, Kzij, w,
botinflux, topinflux, botoutfluxcoef, topoutfluxcoef,
kb, kt, NULL, NULL,sm->dzface_sed))
i_set_error(sediment->hmodel, sm->col_number, LFATAL, "trvdiffsettl:vdiff_dissolved_sed: Error encountered in implicit_vadv_vdiff().\n");
/* END MH */
/* free temporary arrays */
d_free_1d(Kzij);
d_free_1d(w);
d_free_1d(var);
}
/*********************************************************************
File: FIvdiffsettl.c
Created: Wed Jul 14 18:03:11 EST 1993
Author:
Purpose: A general routine to implement fully
implicit vertical mixing for a single
column (i,j) of any of
the model 3d variables. The calling
routine should itself deal with trivial
cases where there is only 1 wet layer.
Arguments: var - pointer to the 1d array of input values
dvar - output 1d array of changes (ie new value
is var+dvar)
Kz - 1d array of mixing coefficients at
layer interfaces
fb - flux at bottom (+ve up)
ft - flux at water surface (+ve up)
kb - bottom k index
kt - top k index
Splus - Positive part of source term (see below)
Sminus - Negative part of source term (see below)
Returns: void
Revisions:
11/01/1999 SJW
Temporary arrays are now locally allocated.
Added parameters and code to allow sources/sinks.
Source linearisation and positivity treatment
according to Patankar 1980. Source is make up of
two parts:
S = Splus - Sminus
where Splus and Sminus are both positive.
07/2012 MH
Made the function type int, returning 1 if it fails.
$Id: trvdiffsettl.c 5848 2018-06-29 05:01:15Z riz008 $
*********************************************************************/
static int implicit_vadv_vdiff(double dt, int nz, double *dz, double *dzold, double *por, double *porold, double *var, /* Array of variable values */
double *dvar, /* Array of changes in values */
double *Kz, /* Diffusivity values */
double *w,
double botinflux, /* Flux out of bottom */
double topinflux, /* Flux out of top */
double botoutfluxcoef, /* Flux out of bottom */
double topoutfluxcoef, /* Flux out of top */
int kb,
int kt, /* Vertical indices for water bottom and top */
double *Splus, /* Positive part of source term */
double *Sminus, /* Negative part of source term */
double *dzface /*UR added, calculate the dface ones, not for every tracer */
)
{
int k = 0;
double dzdt, dzdtold;
double div;
#if HAVE_ALLOCA
double *Cm1 = (double *)alloca((nz + 1) * sizeof(double));
double *C = (double *)alloca((nz + 1) * sizeof(double));
double *Cp1 = (double *)alloca((nz + 1) * sizeof(double));
double *rhs = (double *)alloca((nz + 1) * sizeof(double));
double *sol = (double *)alloca((nz + 1) * sizeof(double));
double *ud = (double *)alloca((nz + 1) * sizeof(double));
/* double *dzface = (double *)alloca((nz + 1) * sizeof(double)); */
#else
double *Cm1 = d_alloc_1d(nz + 1);
double *C = d_alloc_1d(nz + 1);
double *Cp1 = d_alloc_1d(nz + 1);
double *rhs = d_alloc_1d(nz + 1);
double *sol = d_alloc_1d(nz + 1);
double *ud = d_alloc_1d(nz + 1);
/* double *dzface = d_alloc_1d(nz + 1); */
#endif
/* Single layer - the calling routine should deal with this */
if (kt <= kb) {
/* MH 07/2012: included instability handling */
sedtag(LWARN,"sed:trvdiffsettl:implicit_vadv_vdiff"," less than 2 layers\n");
return(1);
/*exit(1);*/
/* END MH */
}
/* Calculate dz values at interfaces */
/* for (k = kb + 1; k <= kt; k++)
dzface[k] = (dz[k - 1] + dz[k]) / 2.;
*/
/* Set up tri-diagonal set of equations */
/* Bottom layer */
dzdt = por[kb] * dz[kb] / dt;
dzdtold = porold[kb] * dzold[kb] / dt;
Cm1[kb] = 0.0;
Cp1[kb] = -Kz[kb + 1] / dzface[kb + 1];
C[kb] = dzdt - Cp1[kb];
rhs[kb] = dzdtold * var[kb] + botinflux;
if (w[kb + 1] < 0.)
Cp1[kb] += w[kb + 1];
else
C[kb] += w[kb + 1];
C[kb] -= botoutfluxcoef;
/* Mid-water layers */
for (k = kb + 1; k < kt; k++) {
dzdt = por[k] * dz[k] / dt;
dzdtold = porold[k] * dzold[k] / dt;
Cm1[k] = -Kz[k] / dzface[k];
Cp1[k] = -Kz[k + 1] / dzface[k + 1];
C[k] = dzdt - Cm1[k] - Cp1[k];
rhs[k] = dzdtold * var[k];
if (w[k + 1] < 0.)
Cp1[k] += w[k + 1];
else
C[k] += w[k + 1];
if (w[k] < 0.)
C[k] -= w[k];
else
Cm1[k] -= w[k];
}
/* Surface layer */
dzdt = por[kt] * dz[kt] / dt;
dzdtold = porold[kt] * dzold[kt] / dt;
Cm1[kt] = -Kz[kt] / dzface[kt];
Cp1[kt] = 0.0;
C[kt] = dzdt - Cm1[kt];
rhs[kt] = dzdtold * var[kt] - topinflux;
if (w[kt] < 0.)
C[kt] -= w[kt];
else
Cm1[kt] -= w[kt];
C[kt] += topoutfluxcoef;
/* Add positive part of source terms, if specified */
if (Splus)
for (k = kb; k <= kt; k++)
rhs[k] += dz[k] * Splus[k];
/* Add negative part of source terms, if specified */
if (Sminus)
for (k = kb; k <= kt; k++)
C[k] += dz[k] * Sminus[k] / var[k];
/* Solve tridiagonal system */
div = C[kb];
sol[kb] = rhs[kb] / div;
for (k = kb + 1; k <= kt; k++) {
ud[k] = Cp1[k - 1] / div;
div = C[k] - Cm1[k] * ud[k];
if (div == 0.0) {
/* MH 07/2012: included instability handling */
sedtag(LWARN,"sed:trvdiffsettl:implicit_vadv_vdiff","(solvetri): zero divisor\n");
return(1);
/*exit(1);*/
/* END MH */
}
sol[k] = (rhs[k] - Cm1[k] * sol[k - 1]) / div;
}
dvar[kt] = sol[kt] - var[kt];
for (k = kt - 1; k >= kb; k--) {
sol[k] -= ud[k + 1] * sol[k + 1];
dvar[k] = sol[k] - var[k];
}
#if !HAVE_ALLOCA
/* Free temporary storage */
d_free_1d(Cm1);
d_free_1d(C);
d_free_1d(Cp1);
d_free_1d(rhs);
d_free_1d(sol);
d_free_1d(ud);
/*2007jul d_free_1d(dzface); */
#endif
return(0);
}
#ifdef __cplusplus
}
#endif