1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Module Library: //
9 // sim_landscape_evolution //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // salem.cpp //
14 // //
15 // Michael Bock, Olaf Conrad //
16 // (C) 2017 //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49
50
51 ///////////////////////////////////////////////////////////
52 // //
53 // //
54 // //
55 ///////////////////////////////////////////////////////////
56
57 //---------------------------------------------------------
58 #include "salem.h"
59
60
61 ///////////////////////////////////////////////////////////
62 // //
63 // Climate //
64 // //
65 ///////////////////////////////////////////////////////////
66
67 //---------------------------------------------------------
CSaLEM_Climate(void)68 CSaLEM_Climate::CSaLEM_Climate(void)
69 {
70 m_pTrend = NULL;
71 m_pAnnual = NULL;
72 }
73
74 //---------------------------------------------------------
Destroy(void)75 void CSaLEM_Climate::Destroy(void)
76 {
77 m_pTrend = NULL;
78 m_pAnnual = NULL;
79 }
80
81 //---------------------------------------------------------
Add_Parameters(CSG_Parameters & Parameters,const CSG_String & Parent)82 bool CSaLEM_Climate::Add_Parameters(CSG_Parameters &Parameters, const CSG_String &Parent)
83 {
84 if( !Parent.is_Empty() )
85 {
86 Parameters.Add_Node(Parent, Parent, _TL("Climate"), _TL(""));
87 }
88
89 //-----------------------------------------------------
90 Parameters.Add_Table(Parent,
91 "TREND" , _TL("Long-term Temperature Signal"),
92 _TL("Long-term temperature signal, used as adjustment for annual scenarios, i.e. to let their original values appear cooler or warmer."),
93 PARAMETER_INPUT
94 );
95
96 Parameters.Add_Table_Field("TREND",
97 "TREND_YEAR" , _TL("Year"),
98 _TL("Time expected as 1000 years before present (ka BP).")
99 );
100
101 Parameters.Add_Table_Field("TREND",
102 "TREND_T" , _TL("Temperature"),
103 _TL("Temperature expected as degree Celsius.")
104 );
105
106 Parameters.Add_Double("TREND",
107 "TREND_T_OFFSET", _TL("Temperature Offset"),
108 _TL("Temperature offset (degree Celsius)."),
109 31.0
110 );
111
112 //-----------------------------------------------------
113 Parameters.Add_Table(Parent,
114 "ANNUAL" , _TL("Annual Climate"),
115 _TL(""),
116 PARAMETER_INPUT
117 );
118
119 Parameters.Add_Table_Field("ANNUAL",
120 "ANNUAL_T" , _TL("Mean Temperature"),
121 _TL("")
122 );
123
124 Parameters.Add_Table_Field("ANNUAL",
125 "ANNUAL_TMIN" , _TL("Minimum Temperature"),
126 _TL("")
127 );
128
129 Parameters.Add_Table_Field("ANNUAL",
130 "ANNUAL_TMAX" , _TL("Maximum Temperature"),
131 _TL("")
132 );
133
134 Parameters.Add_Table_Field("ANNUAL",
135 "ANNUAL_P" , _TL("Precipitation"),
136 _TL("")
137 );
138
139 Parameters.Add_Choice("ANNUAL",
140 "ANNUAL_T_UNIT" , _TL("Temperature Unit"),
141 _TL(""),
142 CSG_String::Format("%s|%s|",
143 _TL("Celsius"),
144 _TL("Kelvin")
145 ), 0
146 );
147
148 Parameters.Add_Double(Parent,
149 "T_LAPSE" , _TL("Temperature Lapse Rate"),
150 _TL("Temperature lapse rate as degree Celsius per 100 meter."),
151 0.6, 0.0, true
152 );
153
154 Parameters.Add_Bool("T_LAPSE",
155 "T_LAPSE_CELL" , _TL("Temperature Height Correction"),
156 _TL("Cellwise temperature correction applying specified temperature lapse rate to surface elevation."),
157 true
158 );
159
160 //-----------------------------------------------------
161 return( true );
162 }
163
164 //---------------------------------------------------------
Set_Parameters(CSG_Parameters & Parameters)165 bool CSaLEM_Climate::Set_Parameters(CSG_Parameters &Parameters)
166 {
167 m_pTrend = Parameters("TREND" )->asTable();
168 m_pAnnual = Parameters("ANNUAL")->asTable();
169
170 if( !m_pTrend || !m_pTrend ->is_Valid() || m_pTrend ->Get_Count() < 1
171 || !m_pAnnual || !m_pAnnual->is_Valid() || m_pAnnual->Get_Count() < 12 )
172 {
173 return( false );
174 }
175
176 m_fTrend_Year = Parameters("TREND_YEAR" )->asInt();
177 m_fTrend_T = Parameters("TREND_T" )->asInt();
178
179 m_fAnnual_T = Parameters("ANNUAL_T" )->asInt();
180 m_fAnnual_Tmin = Parameters("ANNUAL_TMIN" )->asInt();
181 m_fAnnual_Tmax = Parameters("ANNUAL_TMAX" )->asInt();
182 m_fAnnual_P = Parameters("ANNUAL_P" )->asInt();
183
184 m_TLapse = Parameters("T_LAPSE" )->asDouble() / 100.0;
185 m_TLapse_bCell = Parameters("T_LAPSE_CELL" )->asBool();
186
187 m_T_Offset = Parameters("TREND_T_OFFSET")->asDouble() - (Parameters("ANNUAL_T_UNIT")->asInt() == 1 ? 273.15 : 0.0);
188
189 m_pTrend->Set_Index(m_fTrend_Year, TABLE_INDEX_Ascending);
190
191 m_Scenario = 0;
192
193 return( true );
194 }
195
196 //---------------------------------------------------------
Set_Year(int Year_BP)197 bool CSaLEM_Climate::Set_Year(int Year_BP)
198 {
199 if( !m_pTrend || !m_pTrend->is_Valid() || m_pTrend->Get_Count() <= 0 )
200 {
201 return( false );
202 }
203
204 //-----------------------------------------------------
205 int i = 0;
206
207 double kYear_BP = -0.001 * Year_BP;
208
209 while( i < m_pTrend->Get_Count() && m_pTrend->Get_Record_byIndex(i)->asDouble(m_fTrend_Year) <= kYear_BP ) // kYears BP (ascending)
210 {
211 i++;
212 }
213
214 //-----------------------------------------------------
215 if( i >= m_pTrend->Get_Count() )
216 {
217 m_TTrend = m_pTrend->Get_Record_byIndex(i - 1)->asDouble(m_fTrend_T);
218 }
219 else if( i <= 0 )
220 {
221 m_TTrend = m_pTrend->Get_Record_byIndex( 0)->asDouble(m_fTrend_T);
222 }
223 else
224 {
225 CSG_Table_Record *pGE = m_pTrend->Get_Record_byIndex(i);
226 CSG_Table_Record *pLT = m_pTrend->Get_Record_byIndex(i - 1);
227
228 double dY = pGE->asDouble(m_fTrend_Year) - pLT->asDouble(m_fTrend_Year);
229 double dT = pGE->asDouble(m_fTrend_T ) - pLT->asDouble(m_fTrend_T );
230
231 m_TTrend = pLT->asDouble(m_fTrend_T) + (kYear_BP - pLT->asDouble(m_fTrend_Year)) * dT / dY;
232 }
233
234 m_TTrend += m_T_Offset;
235
236 return( true );
237 }
238
239 //---------------------------------------------------------
Set_Month(int Month)240 bool CSaLEM_Climate::Set_Month(int Month)
241 {
242 if( Month == 0 )
243 {
244 m_Scenario = (m_Scenario + 1) % (m_pAnnual->Get_Count() / 12); // number of scenario years
245 }
246
247 CSG_Table_Record *pMonth = m_pAnnual->Get_Record(12 * m_Scenario + Month % 12);
248
249 m_T = pMonth->asDouble(m_fAnnual_T ) + m_TTrend;
250 m_Tmin = pMonth->asDouble(m_fAnnual_Tmin) + m_TTrend;
251 m_Tmax = pMonth->asDouble(m_fAnnual_Tmax) + m_TTrend;
252 m_P = pMonth->asDouble(m_fAnnual_P );
253
254 return( true );
255 }
256
257
258 ///////////////////////////////////////////////////////////
259 // //
260 // Bedrock //
261 // //
262 ///////////////////////////////////////////////////////////
263
264 //---------------------------------------------------------
CSaLEM_Bedrock(void)265 CSaLEM_Bedrock::CSaLEM_Bedrock(void)
266 {
267 m_pRocks = NULL;
268
269 m_Weathering[0] = NULL;
270 m_Weathering[1] = NULL;
271 }
272
273 //---------------------------------------------------------
Destroy(void)274 void CSaLEM_Bedrock::Destroy(void)
275 {
276 m_pRocks = NULL;
277
278 if( m_Weathering[0] ) { delete[](m_Weathering[0]); m_Weathering[0] = NULL; }
279 if( m_Weathering[1] ) { delete[](m_Weathering[1]); m_Weathering[1] = NULL; }
280 }
281
282 //---------------------------------------------------------
Add_Parameters(CSG_Parameters & Parameters,const CSG_String & Parent)283 bool CSaLEM_Bedrock::Add_Parameters(CSG_Parameters &Parameters, const CSG_String &Parent)
284 {
285 if( !Parent.is_Empty() )
286 {
287 Parameters.Add_Node(Parent, Parent, _TL("Bedrock"), _TL(""));
288 }
289
290 Parameters.Add_Grid_List(Parent, "ROCK_LAYERS", _TL("Lithology"), _TL(""), PARAMETER_INPUT_OPTIONAL);
291
292 CSG_Table t;
293
294 t.Add_Field(_TL("Frost" ), SG_DATATYPE_String);
295 t.Add_Field(_TL("Chemical"), SG_DATATYPE_String);
296
297 #define ADD_FORMULAS(a, b) { CSG_Table_Record *p = t.Add_Record(); p->Set_Value(0, a); p->Set_Value(1, b); }
298
299 ADD_FORMULAS("0.0250 * (0.00175 * R + T - Tmax) / (-Tamp * cos(S))", "0.0006 * (exp(-4*R) - exp(-6*R)) * (P/0.0002)");
300 ADD_FORMULAS("0.0250 * (0.00175 * R + T - Tmax) / (-Tamp * cos(S))", "0.0006 * (exp(-4*R) - exp(-6*R)) * (P/0.0002)");
301 ADD_FORMULAS("0.0125 * (0.03750 * R + T - Tmax) / (-Tamp * cos(S))", "0.0050 * (exp(-4*R) - exp(-6*R)) * (P/0.0002)");
302 ADD_FORMULAS("0.0075 * (0.07500 * R + T - Tmax) / (-Tamp * cos(S))", "0.0075 * (exp(-4*R) - exp(-6*R)) * (P/0.0002)");
303 ADD_FORMULAS("0.0250 * (0.03500 * R + T - Tmax) / (-Tamp * cos(S))", "0.0050 * (exp(-4*R) - exp(-6*R)) * (P/0.0002)");
304 ADD_FORMULAS("0.0200 * (0.08750 * R + T - Tmax) / (-Tamp * cos(S))", "0.0060 * (exp(-4*R) - exp(-6*R)) * (P/0.0002)");
305
306 Parameters.Add_FixedTable(Parent,
307 "WEATHERING", _TL("Weathering Formulas"),
308 _TL("Variables that can be used in formulas are 'T', 'Tmin', 'Tmax', 'Tamp' (mean, minimum, maximum, range of temperature), 'P' (precipitation), 'R' (regolith thickness), 'S' (slope gradient)."),
309 &t
310 );
311
312 Parameters.Add_Node(Parent,
313 "DEFAULTS" , _TL("Defaults"),
314 _TL("")
315 );
316
317 Parameters.Add_String("DEFAULTS",
318 "FROST" , _TL("Frost Weathering"),
319 _TL(""),
320 "0.0250 * (0.00175 * R + T - Tmax) / (-Tamp * cos(S))"
321 );
322
323 Parameters.Add_String("DEFAULTS",
324 "CHEMICAL" , _TL("Chemical Weathering"),
325 _TL(""),
326 "0.0002 * exp(-5.0 * R)"
327 );
328
329 return( true );
330 }
331
332 //---------------------------------------------------------
Set_Parameters(CSG_Parameters & Parameters)333 bool CSaLEM_Bedrock::Set_Parameters(CSG_Parameters &Parameters)
334 {
335 Destroy();
336
337 m_pRocks = Parameters("ROCK_LAYERS")->asGridList();
338
339 //-----------------------------------------------------
340 int n = m_pRocks->Get_Grid_Count();
341
342 m_Weathering[0] = new CSG_Formula[1 + n];
343 m_Weathering[1] = new CSG_Formula[1 + n];
344
345 m_Weathering[0][n].Set_Formula(Get_Weathering_Formula(Parameters("FROST" )->asString(), "0")); // default
346 m_Weathering[1][n].Set_Formula(Get_Weathering_Formula(Parameters("CHEMICAL")->asString(), "0")); // default
347
348 CSG_Table &Formulas = *Parameters("WEATHERING")->asTable();
349
350 for(int i=0; i<n; i++)
351 {
352 if( i < Formulas.Get_Count() )
353 {
354 m_Weathering[0][i].Set_Formula(Get_Weathering_Formula(Formulas[i].asString(0), m_Weathering[0][n].Get_Formula()));
355 m_Weathering[1][i].Set_Formula(Get_Weathering_Formula(Formulas[i].asString(1), m_Weathering[1][n].Get_Formula()));
356 }
357 else
358 {
359 m_Weathering[0][i].Set_Formula(m_Weathering[0][n].Get_Formula());
360 m_Weathering[1][i].Set_Formula(m_Weathering[1][n].Get_Formula());
361 }
362 }
363
364 return( true );
365 }
366
367
368 ///////////////////////////////////////////////////////////
369 // //
370 ///////////////////////////////////////////////////////////
371
372 //---------------------------------------------------------
Get_Bedrock_Index(int x,int y,double z) const373 int CSaLEM_Bedrock::Get_Bedrock_Index(int x, int y, double z) const
374 {
375 int iMax = -1;
376 double zMax;
377
378 for(int i=0; i<m_pRocks->Get_Grid_Count(); i++)
379 {
380 if( !m_pRocks->Get_Grid(i)->is_NoData(x, y) )
381 {
382 double iz = m_pRocks->Get_Grid(i)->asDouble(x, y);
383
384 if( z <= iz && (iMax < 0 || iz < zMax) )
385 {
386 iMax = i;
387 zMax = iz;
388 }
389 }
390 }
391
392 return( iMax < 0 ? m_pRocks->Get_Grid_Count() : iMax );
393 }
394
395 //---------------------------------------------------------
Get_Bedrock_Name(int x,int y,double z) const396 CSG_String CSaLEM_Bedrock::Get_Bedrock_Name(int x, int y, double z) const
397 {
398 int i = Get_Bedrock_Index(x, y, z);
399
400 if( i < m_pRocks->Get_Grid_Count() )
401 {
402 return( m_pRocks->Get_Grid(i)->Get_Name() );
403 }
404
405 return( _TL("unknown") );
406 }
407
408
409 ///////////////////////////////////////////////////////////
410 // //
411 ///////////////////////////////////////////////////////////
412
413 //---------------------------------------------------------
414 enum
415 {
416 WV_Tmin , // Temperature minimum
417 WV_Tmax , // Temperature maximum
418 WV_Tamp , // Temperature amplitude
419 WV_T , // Temperature mean
420 WV_P , // Precipitation
421 WV_S , // Slope
422 WV_R , // Regolith cover thickness
423 WV_Count
424 };
425
426 //---------------------------------------------------------
Get_Weathering_Formula(const CSG_String & Formula,const CSG_String & Default)427 CSG_String CSaLEM_Bedrock::Get_Weathering_Formula(const CSG_String &Formula, const CSG_String &Default)
428 {
429 const char Var_Code[] = "abcdefghijklmnopqrstuvwxyz";
430 // 01234567890134567890123456
431 // V 012345678901234567
432
433 const char Var_Key[WV_Count][8] =
434 {
435 "Tmin", // Temperature minimum
436 "Tmax", // Temperature maximum
437 "Tamp", // Temperature amplitude
438 "T" , // Temperature mean
439 "P" , // Precipitation mean
440 "S" , // Slope
441 "R" // Regolith cover thickness
442 };
443
444 //-----------------------------------------------------
445 CSG_String s = Formula; s.Trim(true); s.Trim(false);
446
447 for(int i=0; i<WV_Count; i++)
448 {
449 s.Replace(Var_Key[i], Var_Code[i]);
450 }
451
452 CSG_Formula f;
453
454 if( !f.Set_Formula(s) )
455 {
456 f.Get_Error(s);
457
458 return( Default );
459 }
460
461 return( s );
462 }
463
464
465 ///////////////////////////////////////////////////////////
466 // //
467 ///////////////////////////////////////////////////////////
468
469 //---------------------------------------------------------
Set_Weathering(double dTime,CSaLEM_Climate & Climate,CSG_Grid & Surface,CSG_Grid & Slope,CSG_Grid & Regolith)470 bool CSaLEM_Bedrock::Set_Weathering(double dTime, CSaLEM_Climate &Climate, CSG_Grid &Surface, CSG_Grid &Slope, CSG_Grid &Regolith)
471 {
472 //-----------------------------------------------------
473 dTime *= 0.01 / 12.;
474
475 for(int Month=0; Month<12; Month++)
476 {
477 Climate.Set_Month(Month);
478
479 #pragma omp parallel for
480 for(int y=0; y<Surface.Get_NY(); y++)
481 {
482 CSG_Vector Values(WV_Count);
483
484 Values[WV_T ] = Climate.Get_T ();
485 Values[WV_Tmin] = Climate.Get_Tmin();
486 Values[WV_Tmax] = Climate.Get_Tmax();
487 Values[WV_Tamp] = Climate.Get_Tamp();
488 Values[WV_P ] = Climate.Get_P ();
489
490 for(int x=0; x<Surface.Get_NX(); x++)
491 {
492 Values[WV_S] = Slope .asDouble(x, y);
493 Values[WV_R] = Regolith.asDouble(x, y);
494
495 if( Climate.Get_TLapse_Cellwise() )
496 {
497 double dT = Climate.Get_TLapse() * Surface.asDouble(x, y);
498
499 Values[WV_T ] = Climate.Get_T () - dT;
500 Values[WV_Tmin] = Climate.Get_Tmin() - dT;
501 Values[WV_Tmax] = Climate.Get_Tmax() - dT;
502 }
503
504 if( Values[WV_Tmax] > 0.0 ) // no significant weathering under permanent frost !
505 {
506 int i = Get_Bedrock_Index(x, y, Surface.asDouble(x, y) - Regolith.asDouble(x, y));
507
508 double dR = 0.0;
509
510 if( Values[WV_Tmin] < 0.0 ) { dR += m_Weathering[0][i].Get_Value(Values); } // frost change
511 if( Values[WV_T ] > 0.0 ) { dR += m_Weathering[1][i].Get_Value(Values); } // chemical weathering
512
513 if( dR > 0.0 )
514 {
515 Regolith.Add_Value(x, y, dTime * dR);
516 }
517 }
518 }
519 }
520 }
521
522 return( true );
523 }
524
525
526 ///////////////////////////////////////////////////////////
527 // //
528 // Tracers //
529 // //
530 ///////////////////////////////////////////////////////////
531
532 //---------------------------------------------------------
533 enum
534 {
535 TRACER_TID = 0,
536 TRACER_ROCKTYPE,
537 TRACER_ORIGIN_X,
538 TRACER_ORIGIN_Y,
539 TRACER_HEIGHT,
540 TRACER_DEPTH,
541 TRACER_DISTANCE,
542 TRACER_T_WEATHERED,
543 TRACER_T_MOV_FIRST,
544 TRACER_T_MOV_LAST
545 };
546
547 //---------------------------------------------------------
CSaLEM_Tracers(void)548 CSaLEM_Tracers::CSaLEM_Tracers(void)
549 {
550 m_Candidates.Create(SHAPE_TYPE_Point, _TL("Tracers"), NULL, SG_VERTEX_TYPE_XYZ);
551
552 m_Candidates.Add_Field("TID" , SG_DATATYPE_Int );
553 m_Candidates.Add_Field("ROCKTYPE" , SG_DATATYPE_String);
554 m_Candidates.Add_Field("ORIGIN_X" , SG_DATATYPE_Double);
555 m_Candidates.Add_Field("ORIGIN_Y" , SG_DATATYPE_Double);
556 m_Candidates.Add_Field("HEIGHT" , SG_DATATYPE_Double);
557 m_Candidates.Add_Field("DEPTH" , SG_DATATYPE_Double);
558 m_Candidates.Add_Field("DISTANCE" , SG_DATATYPE_Double);
559 m_Candidates.Add_Field("T_WEATHERED", SG_DATATYPE_Int );
560 m_Candidates.Add_Field("T_MOV_FIRST", SG_DATATYPE_Int );
561 m_Candidates.Add_Field("T_MOV_LAST" , SG_DATATYPE_Int );
562 }
563
564 //---------------------------------------------------------
Destroy(void)565 void CSaLEM_Tracers::Destroy(void)
566 {
567 m_Candidates.Del_Records();
568
569 m_Bedrock.Destroy();
570
571 if( m_Trim == 1 && m_pPoints )
572 {
573 for(int i=m_Trim_Points.Get_Count()-1; i>=0; i--) // remove those tracers and paths that have left the scene
574 {
575 m_pPoints->Add_Shape(m_Trim_Points.Get_Shape(i)); m_Trim_Points.Del_Shape(i);
576
577 if( m_pLines )
578 {
579 m_pLines->Add_Shape(m_Trim_Lines.Get_Shape(i)); m_Trim_Lines.Del_Shape(i);
580 }
581 }
582
583 m_Trim_Points.Destroy();
584 m_Trim_Lines .Destroy();
585 }
586
587 m_pPoints = NULL;
588 m_pLines = NULL;
589 }
590
591 //---------------------------------------------------------
Add_Parameters(CSG_Parameters & Parameters,const CSG_String & Parent)592 bool CSaLEM_Tracers::Add_Parameters(CSG_Parameters &Parameters, const CSG_String &Parent)
593 {
594 if( !Parent.is_Empty() )
595 {
596 Parameters.Add_Node(Parent, Parent, _TL("Tracers"), _TL(""));
597 }
598
599 Parameters.Add_Shapes(Parent,
600 "POINTS" , _TL("Tracer"),
601 _TL(""),
602 PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Point
603 );
604
605 Parameters.Add_Shapes(Parent,
606 "LINES" , _TL("Tracer Paths"),
607 _TL(""),
608 PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Line
609 );
610
611 Parameters.Add_Choice("",
612 "TRIM" , _TL("Trim"),
613 _TL("Remove tracers and associated paths that moved out of the investigated area (performance gain)."),
614 CSG_String::Format("%s|%s|%s|",
615 _TL("no"),
616 _TL("temporarily"),
617 _TL("permanently")
618 ), 1
619 );
620
621 Parameters.Add_Double("",
622 "DIR_RAND" , _TL("Randomize Direction"),
623 _TL("Uncertainty in routing direction expressed as standard deviation (degree)."),
624 0.0, 0.0, true
625 );
626
627 Parameters.Add_Int("",
628 "H_DENSITY" , _TL("Horizontal Density"),
629 _TL("Horizontal tracer density in cells."),
630 1, 1, true
631 );
632
633 Parameters.Add_Bool("H_DENSITY",
634 "H_RANDOM" , _TL("Randomize"),
635 _TL(""),
636 true
637 );
638
639 Parameters.Add_Double("",
640 "V_DENSITY" , _TL("Vertical Density"),
641 _TL("Vertical tracer density in meter."), 0.5, 0.01, true
642 );
643
644 Parameters.Add_Bool("V_DENSITY",
645 "V_RANDOM" , _TL("Randomize"),
646 _TL(""),
647 true
648 );
649
650 return( true );
651 }
652
653 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)654 int CSaLEM_Tracers::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
655 {
656 if( pParameters->Cmp_Identifier("TRACERS") )
657 {
658 if( pParameter->Cmp_Identifier("POINTS") )
659 {
660 pParameters->Set_Enabled("LINES" , pParameter->asDataObject() != NULL);
661 pParameters->Set_Enabled("TRIM" , pParameter->asDataObject() != NULL);
662 pParameters->Set_Enabled("DIR_RAND" , pParameter->asDataObject() != NULL);
663 pParameters->Set_Enabled("H_DENSITY", pParameter->asDataObject() != NULL);
664 pParameters->Set_Enabled("V_DENSITY", pParameter->asDataObject() != NULL);
665 }
666 }
667
668 return( 0 );
669 }
670
671 //---------------------------------------------------------
Set_Parameters(CSG_Parameters & Parameters,CSG_Grid * pSurface)672 bool CSaLEM_Tracers::Set_Parameters(CSG_Parameters &Parameters, CSG_Grid *pSurface)
673 {
674 Destroy();
675
676 if( !pSurface || !pSurface->is_Valid() )
677 {
678 return( false );
679 }
680
681 if( (m_pPoints = Parameters("POINTS")->asShapes()) != NULL )
682 {
683 m_pPoints->Create(SHAPE_TYPE_Point, _TL("Tracers"), &m_Candidates, SG_VERTEX_TYPE_XYZ);
684
685 if( (m_pLines = Parameters("LINES")->asShapes()) != NULL )
686 {
687 m_pLines->Create(SHAPE_TYPE_Line, _TL("Tracer Paths"), NULL, SG_VERTEX_TYPE_XYZ);
688 m_pLines->Add_Field("TID" , SG_DATATYPE_Int);
689 m_pLines->Add_Field("ROCKTYPE", SG_DATATYPE_String);
690 }
691
692 m_Trim = Parameters("TRIM" )->asInt ();
693 m_dDirection = Parameters("DIR_RAND" )->asDouble() * M_DEG_TO_RAD;
694 m_hDensity = Parameters("H_DENSITY")->asInt ();
695 m_hRandom = Parameters("H_RANDOM" )->asBool ();
696 m_vDensity = Parameters("V_DENSITY")->asDouble();
697 m_vRandom = Parameters("V_RANDOM" )->asBool ();
698
699 m_Bedrock.Create(*pSurface);
700
701 if( m_Trim == 1 )
702 {
703 m_Trim_Points.Create(*m_pPoints);
704
705 if( m_pLines )
706 {
707 m_Trim_Lines.Create(*m_pLines);
708 }
709 }
710 }
711
712 return( true );
713 }
714
715 //---------------------------------------------------------
Add_Tracers(double Time,CSG_Grid & Surface,CSG_Grid & Regolith,const CSaLEM_Bedrock & Bedrock)716 int CSaLEM_Tracers::Add_Tracers(double Time, CSG_Grid &Surface, CSG_Grid &Regolith, const CSaLEM_Bedrock &Bedrock)
717 {
718 int nAdded = 0;
719
720 if( m_pPoints )
721 {
722 for(int y=m_hDensity/2; y<Surface.Get_NY(); y+=m_hDensity)
723 {
724 if( y == 0 || y >= Surface.Get_NY() - 1 ) continue;
725
726 for(int x=m_hDensity/2; x<Surface.Get_NX(); x+=m_hDensity)
727 {
728 if( x == 0 || x >= Surface.Get_NX() - 1 || Surface.is_NoData(x, y) ) continue;
729
730 double z, zBedrock = Surface.asDouble(x, y) - Regolith.asDouble(x, y);
731
732 while( zBedrock < (z = m_Bedrock.asDouble(x, y)) )
733 {
734 m_Bedrock.Set_Value(x, y, z - m_vDensity);
735
736 z -= m_vDensity * (m_vRandom ? CSG_Random::Get_Uniform(0, 1) : 0.5);
737
738 CSG_Shape *pTracer = m_Candidates.Add_Shape();
739
740 TSG_Point p = Surface.Get_System().Get_Grid_to_World(x, y);
741
742 if( m_hRandom )
743 {
744 p.x += CSG_Random::Get_Uniform(-0.5, 0.5) * Surface.Get_Cellsize();
745 p.y += CSG_Random::Get_Uniform(-0.5, 0.5) * Surface.Get_Cellsize();
746 }
747
748 pTracer->Set_Point(p, 0); pTracer->Set_Z(z, 0);
749 pTracer->Set_Value(TRACER_ROCKTYPE , Bedrock.Get_Bedrock_Name(x, y, z));
750 pTracer->Set_Value(TRACER_ORIGIN_X , p.x);
751 pTracer->Set_Value(TRACER_ORIGIN_Y , p.y);
752 pTracer->Set_Value(TRACER_HEIGHT , z);
753 pTracer->Set_Value(TRACER_DEPTH , Surface.asDouble(x, y) - z);
754 pTracer->Set_Value(TRACER_DISTANCE , 0.0);
755 pTracer->Set_Value(TRACER_T_WEATHERED, Time);
756
757 nAdded++;
758 }
759 }
760 }
761 }
762
763 return( nAdded );
764 }
765
766 //---------------------------------------------------------
Do_Move(double Depth,double dActive)767 inline bool CSaLEM_Tracers::Do_Move(double Depth, double dActive)
768 {
769 if( dActive < 0.0 || Depth > dActive ) return( false );
770
771 if( Depth <= 0.0 ) return( true );
772
773 return( CSG_Random::Get_Uniform(0, 1) < (Depth / dActive) );
774 }
775
776 //---------------------------------------------------------
Set_Tracers(double Time,double k,CSG_Grid & Surface,CSG_Grid Gradient[3],CSG_Grid & dHin,CSG_Grid & dHout)777 bool CSaLEM_Tracers::Set_Tracers(double Time, double k, CSG_Grid &Surface, CSG_Grid Gradient[3], CSG_Grid &dHin, CSG_Grid &dHout)
778 {
779 if( !m_pPoints )
780 {
781 return( false );
782 }
783
784 CSG_Grid_System System(Surface.Get_System());
785
786 int i;
787
788 //-----------------------------------------------------
789 for(i=m_Candidates.Get_Count()-1; i>=0; i--) // candidates' first move ?
790 {
791 CSG_Shape *pTracer = m_Candidates.Get_Shape(i);
792
793 TSG_Point p = pTracer->Get_Point(0);
794
795 double Height;
796
797 if( Surface.Get_Value(p, Height) )
798 {
799 double Depth = Height - pTracer->asDouble(TRACER_HEIGHT);
800
801 if( Depth <= dHout.Get_Value(p) ) // we take dHout as rough estimation of potentially active layer thickness
802 {
803 pTracer = m_pPoints->Add_Shape(pTracer);
804
805 pTracer->Set_Value(TRACER_TID, m_pPoints->Get_Count());
806 pTracer->Set_Value(TRACER_T_MOV_FIRST, Time);
807
808 m_Candidates.Del_Shape(i);
809
810 if( m_pLines )
811 {
812 CSG_Shape *pLine = m_pLines->Add_Shape();
813
814 pLine->Add_Point(pTracer->Get_Point(0));
815 pLine->Set_Z(pTracer->Get_Z(0), pLine->Get_Point_Count() - 1);
816 pLine->Set_Value(0, pTracer->asInt (TRACER_TID ));
817 pLine->Set_Value(1, pTracer->asString(TRACER_ROCKTYPE));
818 }
819 }
820 }
821 }
822
823 //-----------------------------------------------------
824 #pragma omp parallel for private(i)
825 for(i=0; i<m_pPoints->Get_Count(); i++)
826 {
827 CSG_Shape *pTracer = m_pPoints->Get_Shape(i);
828
829 TSG_Point p = pTracer->Get_Point(0);
830
831 double Height;
832
833 if( Surface.Get_Value(p, Height) )
834 {
835 double Depth = Height - pTracer->asDouble(TRACER_HEIGHT);
836
837 if( !Do_Move(Depth, dHout.Get_Value(p)) ) // tracer will not move because it's deeper than the active layer is thick
838 {
839 pTracer->Set_Value(TRACER_DEPTH, Depth);
840 }
841 else
842 {
843 //-----------------------------------------
844 { // aspect driven routing
845 double s, sin_a, cos_a;
846
847 if( Gradient[0].Get_Value(p, s) && Gradient[1].Get_Value(p, sin_a) && Gradient[2].Get_Value(p, cos_a) )
848 {
849 double tanS = tan(s);
850
851 if( m_dDirection > 0.0 ) // directional uncertainty in dependence of slope gradient
852 {
853 double a = CSG_Random::Get_Gaussian(atan2(sin_a, cos_a), exp(-tanS) * m_dDirection);
854
855 sin_a = sin(a);
856 cos_a = cos(a);
857 }
858
859 double d = tanS * k * System.Get_Cellsize();
860
861 p.x += sin_a * d;
862 p.y += cos_a * d;
863 }
864 }
865
866 //-----------------------------------------
867 if( p.x != pTracer->Get_Point(0).x || p.y != pTracer->Get_Point(0).y ) // tracer has moved
868 {
869 pTracer->Add_Value(TRACER_DISTANCE, SG_Get_Distance(p, pTracer->Get_Point(0)));
870 pTracer->Set_Value(TRACER_T_MOV_LAST, Time);
871 pTracer->Set_Point(p, 0);
872
873 if( Surface.Get_Value(p, Height) ) // returns false once the tracer was moved out of the investigated area
874 {
875 double in = dHin .Get_Value(p);
876 double out = dHout.Get_Value(p);
877
878 Depth = in * pow(CSG_Random::Get_Uniform(0, 1), 3); // bury tracer using some randomness
879 Height += (in - out);
880
881 pTracer->Set_Z ( Height - Depth, 0);
882 pTracer->Set_Value(TRACER_HEIGHT, Height - Depth);
883 pTracer->Set_Value(TRACER_DEPTH , Depth);
884 }
885
886 if( m_pLines )
887 {
888 CSG_Shape *pLine = m_pLines->Get_Shape(i);
889 pLine->Add_Point(pTracer->Get_Point(0));
890 pLine->Set_Z(Height, pLine->Get_Point_Count() - 1);
891 }
892 }
893 }
894 }
895 }
896
897 //-----------------------------------------------------
898 if( m_Trim )
899 {
900 for(i=m_pPoints->Get_Count()-1; i>=0; i--) // remove those tracers and paths that have left the scene
901 {
902 if( !Surface.is_InGrid_byPos(m_pPoints->Get_Shape(i)->Get_Point(0)) )
903 {
904 if( m_Trim == 1 )
905 {
906 m_Trim_Points.Add_Shape(m_pPoints->Get_Shape(i));
907 }
908
909 m_pPoints->Del_Shape(i);
910
911 if( m_pLines )
912 {
913 if( m_Trim == 1 )
914 {
915 m_Trim_Lines.Add_Shape(m_pLines->Get_Shape(i));
916 }
917
918 m_pLines->Del_Shape(i);
919 }
920 }
921 }
922 }
923
924 return( true );
925 }
926
927
928 ///////////////////////////////////////////////////////////
929 // //
930 // //
931 // //
932 ///////////////////////////////////////////////////////////
933
934 //---------------------------------------------------------
CSaLEM(void)935 CSaLEM::CSaLEM(void)
936 {
937 //-----------------------------------------------------
938 Set_Name (_TL("SaLEM"));
939
940 Set_Author ("M.Bock, O.Conrad (c) 2017");
941
942 Set_Description (_TW(
943 "This is the implementation of a Soil and Landscape Evolution Model (SaLEM) "
944 "for the spatiotemporal investigation of soil parent material evolution following "
945 "a lithologically differentiated approach. The model needs a digital elevation "
946 "model and (paleo-)climatic data for the simulation of weathering, erosion and "
947 "transport processes. Weathering is controlled by user defined functions in "
948 "dependence of climate conditions, local slope, regolith cover and outcropping "
949 "bedrock lithology. Lithology can be supplied as a set of grids, of which each "
950 "grid represents the top elevation of the underlying bedrock type. "
951 ));
952
953 Add_Reference(
954 "Bock, M., Conrad, O., Guenther, A., Gehrt, E., Baritz, R., and Boehner, J.", "2018",
955 "SaLEM (v1.0) - the Soil and Landscape Evolution Model (SaLEM) for simulation of regolith depth in periglacial environments",
956 "Geosci. Model Dev., 11, 1641-1652.",
957 SG_T("https://doi.org/10.5194/gmd-11-1641-2018")
958 );
959
960 Add_Reference(
961 "Alley, R.", "2000",
962 "The Younger Dryas cold interval as viewed from central Greenland",
963 "Quaternary Science Reviews 19: 213-226."
964 );
965
966 Add_Reference(
967 "Tucker, G.E. & Slingerland, R.", "1994",
968 "Erosional dynamics, flexural isostasy, and long-lived escarpments: A numerical modeling study",
969 "Journal of Geophysical Research 99: 12229-12243."
970 );
971
972 //-----------------------------------------------------
973 Parameters.Add_Grid("",
974 "SURFACE_T0" , _TL("Initial Surface Elevation"),
975 _TL(""),
976 PARAMETER_INPUT
977 );
978
979 Parameters.Add_Grid_or_Const("",
980 "REGOLITH_T0" , _TL("Initial Regolith Thickness"),
981 _TL("Initial sediment cover (m)"),
982 0.0, 0.0, true
983 );
984
985 Parameters.Add_Grid_or_Const("",
986 "ALLOCHTHONE" , _TL("Allochthone Input"),
987 _TL("Allochthone input in meter per year (e.g. of aeolian origin, 'Loess')."),
988 0.0, 0.0, true
989 );
990
991 Parameters.Add_Grid("",
992 "SURFACE" , _TL("Surface Elevation"),
993 _TL(""),
994 PARAMETER_OUTPUT
995 );
996
997 Parameters.Add_Grid("",
998 "REGOLITH" , _TL("Regolith Thickness"),
999 _TL(""),
1000 PARAMETER_OUTPUT
1001 );
1002
1003 Parameters.Add_Grid("",
1004 "DIFFERENCE" , _TL("Difference"),
1005 _TL(""),
1006 PARAMETER_OUTPUT
1007 );
1008
1009 //-----------------------------------------------------
1010 m_Climate.Add_Parameters(*Parameters.Add_Parameters("", "CLIMATE", _TL("Climate"), _TL(""))->asParameters(), "");
1011 m_Bedrock.Add_Parameters(*Parameters.Add_Parameters("", "BEDROCK", _TL("Bedrock"), _TL(""))->asParameters(), "");
1012 m_Tracers.Add_Parameters(*Parameters.Add_Parameters("", "TRACERS", _TL("Tracers"), _TL(""))->asParameters(), "");
1013
1014 //-----------------------------------------------------
1015 Parameters.Add_Node("", "TIME", _TL("Simulation Time"), _TL(""));
1016
1017 Parameters.Add_Int("TIME", "TIME_START", _TL("Start [Years BP]"), _TL(""), 50000);
1018 Parameters.Add_Int("TIME", "TIME_STOP" , _TL("Stop [Years BP]" ), _TL(""), 10000);
1019 Parameters.Add_Int("TIME", "TIME_STEP" , _TL("Step [Years]" ), _TL(""), 100, 1, true);
1020
1021 //-----------------------------------------------------
1022 Parameters.Add_Node("", "DIFFUSIVE", _TL("Diffusive Hillslope Processes"), _TL(""));
1023
1024 Parameters.Add_Double("DIFFUSIVE",
1025 "DIFFUSIVE_KD", _TL("Diffusivity Coefficient Kd [m^2/a]"),
1026 _TL(""),
1027 0.01, 0.0, true
1028 );
1029
1030 Parameters.Add_Choice("DIFFUSIVE",
1031 "DIFFUSIVE_NEIGHBOURS", _TL("Neighbourhood"),
1032 _TL(""),
1033 CSG_String::Format("%s|%s|",
1034 _TL("Neumann"),
1035 _TL("Moore")
1036 ), 1
1037 );
1038
1039 //-----------------------------------------------------
1040 if( has_GUI() )
1041 {
1042 Parameters.Add_Int("", "UPDATE", _TL("Display Update Frequency"),
1043 _TL("The frequency (simulation time steps) with which associated maps will be updated. Set to zero to switch off."),
1044 1, 0, true
1045 );
1046
1047 Parameters.Add_Bool ("UPDATE" , "UPDATE_VEC", _TL("Vectors" ), _TL(""), true);
1048
1049 Parameters.Add_Bool ("UPDATE" , "UPDATE_ADJ", _TL("Adjust Regolith"), _TL(""), true);
1050 Parameters.Add_Double("UPDATE_ADJ", "UPDATE_MIN", _TL("Minimum" ), _TL(""), 0.0, 0.0, true);
1051 Parameters.Add_Double("UPDATE_ADJ", "UPDATE_MAX", _TL("Maximum" ), _TL(""), 2.0, 0.0, true);
1052 }
1053 }
1054
1055
1056 ///////////////////////////////////////////////////////////
1057 // //
1058 ///////////////////////////////////////////////////////////
1059
1060 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)1061 int CSaLEM::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
1062 {
1063 m_Tracers.On_Parameters_Enable(pParameters, pParameter);
1064
1065 if( pParameter->Cmp_Identifier("UPDATE") )
1066 {
1067 pParameters->Set_Enabled("UPDATE_ADJ", pParameter->asInt() > 0);
1068 pParameters->Set_Enabled("UPDATE_VEC", pParameter->asInt() > 0);
1069 }
1070
1071 if( pParameter->Cmp_Identifier("UPDATE_ADJ") )
1072 {
1073 pParameters->Set_Enabled("UPDATE_MIN", pParameter->asBool());
1074 pParameters->Set_Enabled("UPDATE_MAX", pParameter->asBool());
1075 }
1076
1077 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
1078 }
1079
1080
1081 ///////////////////////////////////////////////////////////
1082 // //
1083 ///////////////////////////////////////////////////////////
1084
1085 //---------------------------------------------------------
On_Execute(void)1086 bool CSaLEM::On_Execute(void)
1087 {
1088 if( !Initialize() )
1089 {
1090 Error_Set(_TL("initialization failed"));
1091
1092 Finalize();
1093
1094 return( false );
1095 }
1096
1097 //-----------------------------------------------------
1098 int time_start = -Parameters("TIME_START")->asInt();
1099 int time_stop = -Parameters("TIME_STOP" )->asInt();
1100
1101 m_dTime = Parameters("TIME_STEP")->asInt();
1102
1103 //-----------------------------------------------------
1104 int Update = Parameters("UPDATE") ? Parameters("UPDATE")->asInt() * m_dTime : 0;
1105
1106 if( Update > 0 && Parameters("UPDATE_ADJ")->asBool() )
1107 {
1108 DataObject_Update(m_pRegolith,
1109 Parameters("UPDATE_MIN")->asDouble(),
1110 Parameters("UPDATE_MAX")->asDouble(), true
1111 );
1112 }
1113
1114 //-----------------------------------------------------
1115 for(m_Time=time_start; m_Time<=time_stop && Set_Progress(m_Time-time_start, time_stop-time_start); m_Time+=m_dTime)
1116 {
1117 Process_Set_Text("%s: %d", _TL("Years BP"), -m_Time);
1118
1119 SG_UI_Progress_Lock(true);
1120
1121 //-------------------------------------------------
1122 Set_Gradient ();
1123
1124 Set_Allochthone();
1125
1126 Set_Weathering ();
1127
1128 Set_Diffusive (); // diffusive regolith creep
1129
1130 //-------------------------------------------------
1131 if( Update > 0 && ((m_Time - time_start) % Update) == 0 )
1132 {
1133 if( Parameters("UPDATE_ADJ")->asBool() )
1134 {
1135 DataObject_Update(m_pRegolith,
1136 Parameters("UPDATE_MIN")->asDouble(),
1137 Parameters("UPDATE_MAX")->asDouble()
1138 );
1139 }
1140 else
1141 {
1142 DataObject_Update(m_pRegolith);
1143 }
1144
1145 if( Parameters("UPDATE_VEC")->asBool() )
1146 {
1147 DataObject_Update(m_Tracers.Get_Tracers());
1148 DataObject_Update(m_Tracers.Get_Paths ());
1149 }
1150 }
1151
1152 SG_UI_Progress_Lock(false);
1153 }
1154
1155 //-----------------------------------------------------
1156 Finalize();
1157
1158 return( true );
1159 }
1160
1161
1162 ///////////////////////////////////////////////////////////
1163 // //
1164 ///////////////////////////////////////////////////////////
1165
1166 //---------------------------------------------------------
Initialize(void)1167 bool CSaLEM::Initialize(void)
1168 {
1169 m_pSurface = Parameters("SURFACE" )->asGrid();
1170 m_pRegolith = Parameters("REGOLITH")->asGrid();
1171
1172 //-----------------------------------------------------
1173 m_pSurface->Assign(Parameters("SURFACE_T0")->asGrid());
1174
1175 if( Parameters("REGOLITH_T0")->asGrid() == NULL )
1176 {
1177 m_pRegolith->Assign(Parameters("REGOLITH_T0")->asDouble());
1178 }
1179 else if( Parameters("REGOLITH_T0")->asGrid() != m_pRegolith )
1180 {
1181 m_pRegolith->Assign(Parameters("REGOLITH_T0")->asGrid());
1182 }
1183
1184 //-----------------------------------------------------
1185 if( !m_Climate.Set_Parameters(*Parameters("CLIMATE")->asParameters())
1186 || !m_Bedrock.Set_Parameters(*Parameters("BEDROCK")->asParameters())
1187 || !m_Tracers.Set_Parameters(*Parameters("TRACERS")->asParameters(), m_pSurface) )
1188 {
1189 return( false );
1190 }
1191
1192 return( true );
1193 }
1194
1195 //---------------------------------------------------------
Finalize(void)1196 bool CSaLEM::Finalize(void)
1197 {
1198 m_Climate.Destroy();
1199 m_Bedrock.Destroy();
1200 m_Tracers.Destroy();
1201
1202 m_Gradient[0].Destroy();
1203 m_Gradient[1].Destroy();
1204 m_Gradient[2].Destroy();
1205
1206 //-----------------------------------------------------
1207 if( Parameters("DIFFERENCE")->asGrid() )
1208 {
1209 CSG_Grid *pDifference = Parameters("DIFFERENCE")->asGrid();
1210 CSG_Grid *pSurface_T0 = Parameters("SURFACE_T0")->asGrid();
1211
1212 #pragma omp parallel for
1213 for(sLong i=0; i<Get_NCells(); i++)
1214 {
1215 if( m_pSurface->is_NoData(i) )
1216 {
1217 pDifference->Set_NoData(i);
1218 }
1219 else
1220 {
1221 pDifference->Set_Value(i, m_pSurface->asDouble(i) - pSurface_T0->asDouble(i));
1222 }
1223 }
1224 }
1225
1226 //-----------------------------------------------------
1227 return( true );
1228 }
1229
1230
1231 ///////////////////////////////////////////////////////////
1232 // //
1233 ///////////////////////////////////////////////////////////
1234
1235 //---------------------------------------------------------
Set_Gradient(void)1236 bool CSaLEM::Set_Gradient(void)
1237 {
1238 if( !Get_System().is_Equal(m_Gradient[0].Get_System()) )
1239 {
1240 m_Gradient[0].Create(Get_System()); // slope gradient [radians]
1241 m_Gradient[1].Create(Get_System()); // sine of aspect
1242 m_Gradient[2].Create(Get_System()); // cosine of aspect
1243 }
1244
1245 #pragma omp parallel for
1246 for(int y=0; y<m_pSurface->Get_NY(); y++)
1247 {
1248 for(int x=0; x<m_pSurface->Get_NX(); x++)
1249 {
1250 double s, a;
1251
1252 bool bOkay = m_pSurface->Get_Gradient(x, y, s, a);
1253
1254 //---------------------------------------------
1255 if( bOkay ) // the following lines make sure that aspect always points downslope (if possible)
1256 {
1257 int i = (int)(0.5 + a * 8. / M_PI_360);
1258
1259 int ix = Get_xTo(i, x);
1260 int iy = Get_yTo(i, y);
1261
1262 if( m_pSurface->is_InGrid(ix, iy) && m_pSurface->asDouble(x, y) < m_pSurface->asDouble(ix, iy) )
1263 {
1264 i = m_pSurface->Get_Gradient_NeighborDir(x, y);
1265
1266 if( (bOkay = i >= 0) == true )
1267 {
1268 a = i * M_PI_360 / 8.;
1269 s = atan((m_pSurface->asDouble(x, y) - m_pSurface->asDouble(ix, iy)) / Get_Length(i));
1270 }
1271 }
1272 }
1273
1274 //---------------------------------------------
1275 if( bOkay )
1276 {
1277 m_Gradient[0].Set_Value(x, y, s );
1278 m_Gradient[1].Set_Value(x, y, sin(a));
1279 m_Gradient[2].Set_Value(x, y, cos(a));
1280 }
1281 else
1282 {
1283 m_Gradient[0].Set_NoData(x, y);
1284 m_Gradient[1].Set_NoData(x, y);
1285 m_Gradient[2].Set_NoData(x, y);
1286 }
1287 }
1288 }
1289
1290 return( true );
1291 }
1292
1293
1294 ///////////////////////////////////////////////////////////
1295 // //
1296 ///////////////////////////////////////////////////////////
1297
1298 //---------------------------------------------------------
Set_Allochthone(void)1299 bool CSaLEM::Set_Allochthone(void)
1300 {
1301 CSG_Grid *pAllochthone = Parameters("ALLOCHTHONE")->asGrid ();
1302 double Allochthone = Parameters("ALLOCHTHONE")->asDouble();
1303
1304 if( pAllochthone || Allochthone > 0.0 )
1305 {
1306 #pragma omp parallel for
1307 for(int y=0; y<Get_NY(); y++)
1308 {
1309 for(int x=0; x<Get_NX(); x++)
1310 {
1311 double d = pAllochthone && !pAllochthone->is_NoData(x, y) ? pAllochthone->asDouble(x, y) : Allochthone;
1312
1313 if( d > 0.0 )
1314 {
1315 d *= m_dTime;
1316
1317 m_pRegolith->Add_Value(x, y, d);
1318 m_pSurface ->Add_Value(x, y, d);
1319 }
1320 }
1321 }
1322 }
1323
1324 return( true );
1325 }
1326
1327
1328 ///////////////////////////////////////////////////////////
1329 // //
1330 ///////////////////////////////////////////////////////////
1331
1332 //---------------------------------------------------------
Set_Weathering(void)1333 bool CSaLEM::Set_Weathering(void)
1334 {
1335 m_Climate.Set_Year(m_Time); // Years BP
1336
1337 m_Bedrock.Set_Weathering(m_dTime, m_Climate, *m_pSurface, m_Gradient[0], *m_pRegolith);
1338
1339 m_Tracers.Add_Tracers(m_Time, *m_pSurface, *m_pRegolith, m_Bedrock);
1340
1341 return( true );
1342 }
1343
1344
1345 ///////////////////////////////////////////////////////////
1346 // //
1347 ///////////////////////////////////////////////////////////
1348
1349 //---------------------------------------------------------
1350 // Reimplements GOLEM::DiffuseExplIrreg
1351 //---------------------------------------------------------
Set_Diffusive(void)1352 bool CSaLEM::Set_Diffusive(void)
1353 {
1354 double k = Parameters("DIFFUSIVE_KD")->asDouble() * m_dTime / Get_Cellarea(); // Diffusivity coefficient Kd [m^2/a]
1355 int y, iStep = Parameters("DIFFUSIVE_NEIGHBOURS")->asInt() == 1 ? 1 : 2;
1356
1357 CSG_Grid dHin (Get_System()); // dHin.Assign(0.0);
1358 CSG_Grid dHout(Get_System());
1359
1360 //-----------------------------------------------------
1361 for(y=0; y<Get_NY(); y++)
1362 {
1363 for(int x=0; x<Get_NX(); x++)
1364 {
1365 if( !m_pSurface->is_NoData(x, y) )
1366 {
1367 int i;
1368
1369 double qs[8], dHout_Sum = 0.0, z = m_pSurface->asDouble(x, y);
1370
1371 for(i=0; i<8; i+=iStep)
1372 {
1373 int ix = Get_xTo(i, x);
1374 int iy = Get_yTo(i, y);
1375
1376 qs[i] = 0.0;
1377
1378 if( m_pSurface->is_InGrid(ix, iy) )
1379 {
1380 double dz = z - m_pSurface->asDouble(ix, iy);
1381
1382 if( dz > 0.0 )
1383 {
1384 dHout_Sum += (qs[i] = dz * k / Get_UnitLength(i));
1385 }
1386 }
1387 else if( m_pSurface->is_InGrid(ix = Get_xFrom(i, x), iy = Get_yFrom(i, y)) )
1388 {
1389 double dz = m_pSurface->asDouble(ix, iy) - z;
1390
1391 if( dz > 0.0 )
1392 {
1393 dHout_Sum += dz * k / Get_UnitLength(i);
1394 }
1395 }
1396 }
1397
1398 double scale = dHout_Sum > m_pRegolith->asDouble(x, y) ? m_pRegolith->asDouble(x, y) / dHout_Sum : 1.0;
1399
1400 dHout.Set_Value(x, y, scale * dHout_Sum);
1401
1402 for(i=0; i<8; i+=iStep)
1403 {
1404 if( qs[i] > 0.0 )
1405 {
1406 dHin.Add_Value(Get_xTo(i, x), Get_yTo(i, y), scale * qs[i]);
1407 }
1408 }
1409 }
1410 }
1411 }
1412
1413 //-----------------------------------------------------
1414 m_Tracers.Set_Tracers(m_Time, Parameters("DIFFUSIVE_KD")->asDouble() * m_dTime, *m_pSurface, m_Gradient, dHin, dHout);
1415
1416 //-----------------------------------------------------
1417 #pragma omp parallel for private(y)
1418 for(y=0; y<Get_NY(); y++)
1419 {
1420 for(int x=0; x<Get_NX(); x++)
1421 {
1422 if( !m_pSurface->is_NoData(x, y) )
1423 {
1424 double R = m_pRegolith->asDouble(x, y);
1425 double dH = dHin.asDouble(x, y) - dHout.asDouble(x, y);
1426
1427 if( R + dH < 0.0 )
1428 {
1429 dH = -R;
1430 }
1431
1432 m_pSurface ->Add_Value(x, y, dH);
1433 m_pRegolith->Add_Value(x, y, dH);
1434 }
1435 }
1436 }
1437
1438 //-----------------------------------------------------
1439 return( true );
1440 }
1441
1442
1443 ///////////////////////////////////////////////////////////
1444 // //
1445 // Relevant Original GOLEM Implementations //
1446 // //
1447 ///////////////////////////////////////////////////////////
1448
1449 //---------------------------------------------------------
1450 // int Dt, // Time step size (years)
1451 // mw, // Depth at which weathering rate decays to 1/e
1452 // kd, // Hillslope transport (diffusivity) coefficient
1453
1454 //---------------------------------------------------------
1455 // double elev [XGRID + 2][YGRID + 2]; // cell elevations
1456 // double chansed [XGRID + 2][YGRID + 2]; // Sediment cover thickness
1457 // double rock [XGRID + 2][YGRID + 2]; // Thickness of uppermost rock layer
1458 // int rockType[XGRID + 2][YGRID + 2]; // Lithology of uppermost rock layer
1459 // int area [XGRID + 2][YGRID + 2]; // Drainage area in cells
1460
1461 //---------------------------------------------------------
1462 // Lithology structure: Version 5.10 adds fields, Kmn and tc. Kmn is the mean,
1463 // unmodified K value---that is, unmodified by changes in precip frequency.
1464 // If precip frequency/duration varies, the "effective" K value will also
1465 // vary. Thus, Kmn stores the mean, and K the value at any particular
1466 // time, as computed by the CalcPrecip function. If precip doesn't vary
1467 // in time, K and Kmn will be identical. tc is a bedrock erosion threshold,
1468 // the units of which depend on exponents mb and nb. It represents that
1469 // value of Q^mb S^nb at which erosion just begins to be initiated.
1470 // struct Lithology // Lithology record
1471 // {
1472 // double Kmn; // Erodability by flowing water (mean)
1473 // double K; // Erodability by flowing water (time-varying)
1474 // double T; // Initial thickness
1475 // double W; // Bedrock weathering rate
1476 // double Scr; // Failure slope
1477 // double tc; // Erosion threshold
1478 // }*layer;
1479
1480 //---------------------------------------------------------
1481 // This function calculates the addition to soil thickness
1482 // during a time step, using an analytical solution to the
1483 // equation dC/dt = W * exp( -C ). With large soil thick-
1484 // nesses, numerical errors in computing e raised to a
1485 // large number can actually result in decreases in the
1486 // computed soil thickness (and increases in the rock
1487 // layer thickness). To avoid this error, the
1488 // computation is only performed for soils less than 100m
1489 // thick.
1490 // Modified for version 5.10 to allow weathering only
1491 // where drainage area, in cells, is less than a specified
1492 // minimum threshold.
1493 //---------------------------------------------------------
1494 /*/Weather()
1495 {
1496 int i, j;
1497 double oldsed, kwdivmw;
1498
1499 for (i = 1; i <= XGRID; i++)
1500 for (j = 1; j <= YGRID; j++)
1501 {
1502 if (area[i][j] < permChannelA && chansed[i][j] <= 100.0)
1503 {
1504 kwdivmw = layer[rockType[i][j]].W / mw;
1505 oldsed = chansed[i][j];
1506 chansed[i][j] = mw * log(exp(chansed[i][j] / mw) + kwdivmw);
1507 rock[i][j] -= (chansed[i][j] - oldsed);
1508 if (rock[i][j] < 0.0)
1509 {
1510 rockType[i][j]++;
1511 rock[i][j] += layer[rockType[i][j]].T;
1512 }
1513 }
1514 }
1515 }/*/
1516
1517 //---------------------------------------------------------
1518 /*/void DiffuseExpl()
1519 {
1520 double dhin[XGRID + 2][YGRID + 2], dhout[XGRIDP1][YGRIDP1], qsn, qss, qse, qsw, dhtot, scale, k;
1521 int i, j, jn, js;
1522
1523 k = kd * Dt / cellarea;
1524
1525 for (i = 1; i <= XGRID; i++)
1526 for (j = 1; j <= YGRID; j++)
1527 dhin[i][j] = 0.0;
1528 for (i = 1; i <= XGRID; i++)
1529 for (j = 1; j <= YGRID; j++)
1530 {
1531 scale = 1.0;
1532 jn = (j == YGRID) ? 1 : j + 1;
1533 js = (j == 1) ? YGRID : j - 1;
1534 qsn = (elev[i][j] > elev[i][jn]) ? k * (elev[i][j] - elev[i][jn]) : 0.0;
1535 qss = (elev[i][j] > elev[i][js]) ? k * (elev[i][j] - elev[i][js]) : 0.0;
1536 qse = (elev[i][j] > elev[i + 1][j]) ? k * (elev[i][j] - elev[i + 1][j]) : 0.0;
1537 qsw = (elev[i][j] > elev[i - 1][j]) ? k * (elev[i][j] - elev[i - 1][j]) : 0.0;
1538 if ((dhtot = qsn + qss + qse + qsw) > chansed[i][j])
1539 scale = chansed[i][j] / dhtot;
1540 dhout[i][j] = scale * dhtot;
1541 dhin[i][jn] += scale * qsn;
1542 dhin[i][js] += scale * qss;
1543 dhin[i + 1][j] += scale * qse;
1544 dhin[i - 1][j] += scale * qsw;
1545
1546 }
1547
1548 // Adjust elevations, sediment reservoirs and denudation rates
1549 for (j = 1; j <= YGRID; j++)
1550 {
1551 for (i = 1; i <= XGRID; i++)
1552 {
1553 chansed[i][j] += dhin[i][j] - dhout[i][j];
1554 elev[i][j] += dhin[i][j] - dhout[i][j];
1555 }
1556
1557 // Accumulate influx at boundaries in sed_yield vbl
1558 sed_yield += dhin[0][j] + dhin[YGRIDP1][j];
1559 }
1560 }*/
1561
1562
1563 ///////////////////////////////////////////////////////////
1564 // //
1565 // //
1566 // //
1567 ///////////////////////////////////////////////////////////
1568
1569 //---------------------------------------------------------
1570