|
| 1 | +#include "depeffects.h" |
| 2 | + |
| 3 | +#include <math.h> |
| 4 | + |
| 5 | +#include "common/context.h" |
| 6 | + |
| 7 | +#include "events.h" |
| 8 | +#include "state.h" |
| 9 | + |
| 10 | +double getClippedWaterFrac(double water, double whc) { |
| 11 | + return fmin(fmax(water / whc, 0.0), 1.0); |
| 12 | +} |
| 13 | + |
| 14 | +double calcAnaerobicIndex(double water, double whc) { |
| 15 | + double f_whc = getClippedWaterFrac(water, whc); |
| 16 | + double f_a = params.fAnoxia; |
| 17 | + |
| 18 | + // Anaerobic index (oxygen limitation proxy) |
| 19 | + return fmin(fmax((f_whc - f_a) / (1 - f_a), 0), 1); |
| 20 | +} |
| 21 | + |
| 22 | +double calcRespMoistEffect(double water, double whc) { |
| 23 | + double moistEffect; |
| 24 | + |
| 25 | + if ((!ctx.waterHResp) || (climate->tsoil < 0)) { |
| 26 | + // if not waterHResp, soil moisture does not affect heterotrophic |
| 27 | + // respiration; also, ignore moisture effects in frozen soils |
| 28 | + // :: from [2], snowpack addition |
| 29 | + moistEffect = 1.0; |
| 30 | + } else { |
| 31 | + double f_whc = getClippedWaterFrac(water, whc); |
| 32 | + if (!ctx.anaerobic) { |
| 33 | + // :: from [1], first part of eq (A20), with added exponent |
| 34 | + // Original formulation from [1], based on PnET is: |
| 35 | + // moistEffect = water / whc |
| 36 | + // which matches here with params.soilRespMoistEffect=1 (the default |
| 37 | + // value) |
| 38 | + // |
| 39 | + // [TAG:UNKNOWN_PROVENANCE] soilRespMoistEffect |
| 40 | + // Note: older versions of sipnet note this as "using PnET formulation", |
| 41 | + // but we have been unable to verify that the exponent comes from PnET |
| 42 | + moistEffect = pow(f_whc, params.soilRespMoistEffect); |
| 43 | + } else { |
| 44 | + // Unimodal moisture response: suppressed under dry conditions, maximal |
| 45 | + // at intermediate moisture, reduced under saturated/anoxic conditions |
| 46 | + |
| 47 | + // Aerobic water availability (dry limitation) |
| 48 | + double D_aer = fmin(fmax(f_whc / params.fAnoxia, 0), 1); |
| 49 | + // Anaerobic index (oxygen limitation proxy) |
| 50 | + double A = calcAnaerobicIndex(water, whc); |
| 51 | + // Uni-modal moisture response |
| 52 | + // Note that this will not go above 1, since: |
| 53 | + // - when f_whc <= f_a, A = 0, and moistEffect = D_aer |
| 54 | + // - when f_whc > f_a, D_aer caps at 1, and this becomes |
| 55 | + // moistEffect = 1 - (1 - adr)A, and adr <= 1 |
| 56 | + moistEffect = (1 - A) * D_aer + params.anaerobicDecompRate * A; |
| 57 | + } |
| 58 | + } |
| 59 | + |
| 60 | + return moistEffect; |
| 61 | +} |
| 62 | + |
| 63 | +double calcMethaneMoistEffect(double water, double whc) { |
| 64 | + // Anaerobic index (oxygen limitation proxy) |
| 65 | + double A = calcAnaerobicIndex(water, whc); |
| 66 | + |
| 67 | + return pow(A, params.anaerobicTransExp); |
| 68 | +} |
| 69 | + |
| 70 | +double calcTempEffect(double tsoil) { |
| 71 | + // :: from [1], D_temp calc as part of eq (A20) |
| 72 | + return pow(params.soilRespQ10, tsoil / 10); |
| 73 | +} |
| 74 | + |
| 75 | +double calcTillageEffect(void) { return 1 + eventTrackers.d_till_mod; } |
| 76 | + |
| 77 | +double calcCNEffect(double kCN, double poolC, double poolN) { |
| 78 | + if (!ctx.nitrogenCycle) { |
| 79 | + return 1.0; |
| 80 | + } |
| 81 | + |
| 82 | + // CN ratio dependency term, using k_CN term that indicates CN value |
| 83 | + // at which the dependency is 1/2 |
| 84 | + double cn = calcRatio(poolC, poolN); |
| 85 | + return kCN / (kCN + cn); |
| 86 | +} |
| 87 | + |
| 88 | +double calcVolatilizationMoistEffect(double water, double whc) { |
| 89 | + // Anaerobic index (oxygen limitation proxy) |
| 90 | + double A = calcAnaerobicIndex(water, whc); |
| 91 | + |
| 92 | + // 0.05 represents the baseline aerobic volatilization |
| 93 | + // The factor of 3.8 makes the max = 1, as with other dependency functions |
| 94 | + return 0.05 + 3.8 * A * (1 - A); |
| 95 | +} |
0 commit comments