methane.c File Reference

#include <R.h>
#include <math.h>
#include <Rmath.h>
#include <Rinternals.h>
#include "soilwater.h"
#include "n2o_model.h"

Include dependency graph for methane.c:

Go to the source code of this file.

Macros
#define	CH4DEPTH   15.0

Functions
void	methane_oxidation (double *CH4, int *isdecid, int *isagri, LAYERPAR_SPT layers, SOIL_SPT soil)

Macro Definition Documentation

#define CH4DEPTH   15.0

Function Documentation

void methane_oxidation	(	double *	CH4,
		int *	isdecid,
		int *	isagri,
		LAYERPAR_SPT	layers,
		SOIL_SPT	soil
	)

    { 
      int    ilyr;
      double  bulkdensity;
      double  fieldcapacity;
      double soiltemp;
      double soilwater;
      double wfps;
      double CH4max;
      double  Dopt;
      double Wmin;
      double Wmax;
      double  Wopt;
      double agri_adjust;
      double temp_adjust;
      double watr_adjust;
      double wfps_adjust;
      double  percentlayer;
      double  temp;

/*
      extern LAYERPAR_SPT layers;
      extern SOIL_SPT soil;
*/
      /* Compute a weighted average for soil temperature, field capacity, */
      /* bulk density, water filled pore space, and volumetric soil water */
      /* content in top 15 cm of soil profile */
      ilyr = 0;
      bulkdensity = 0.0;
      fieldcapacity = 0.0;
      soiltemp = 0.0;
      soilwater = 0.0;
      wfps = 0.0;
      for (ilyr = 0; ilyr < layers->numlyrs; ilyr++) {
        if (layers->dpthmn[ilyr] < CH4DEPTH) {
          if (layers->dpthmx[ilyr] <= CH4DEPTH) {
            bulkdensity += layers->bulkd[ilyr] * layers->width[ilyr] / 
                           CH4DEPTH;
            fieldcapacity += layers->fieldc[ilyr] * layers->width[ilyr] /
                             CH4DEPTH;
            soiltemp += soil->soiltavg[ilyr] * layers->width[ilyr] / CH4DEPTH;
            soilwater += layers->wfps[ilyr] * layers->width[ilyr] / CH4DEPTH;
            wfps += layers->wfps[ilyr] * layers->width[ilyr] / CH4DEPTH;
          } else if ((layers->dpthmx[ilyr] - layers->dpthmn[ilyr]) > 0.0) {
            percentlayer = (CH4DEPTH - layers->dpthmn[ilyr]) / 
                           (layers->dpthmx[ilyr] - layers->dpthmn[ilyr]);
            bulkdensity += layers->bulkd[ilyr] * layers->width[ilyr] /
                           CH4DEPTH * percentlayer;
            fieldcapacity += layers->fieldc[ilyr] * layers->width[ilyr] /
                             CH4DEPTH * percentlayer;
            soiltemp += soil->soiltavg[ilyr] * layers->width[ilyr] /
                        CH4DEPTH * percentlayer;
            soilwater += layers->wfps[ilyr] * layers->width[ilyr] / CH4DEPTH *
                         percentlayer;
            wfps += layers->wfps[ilyr] * layers->width[ilyr] / CH4DEPTH *
                    percentlayer;
          }
        }
      }
      /* Convert from water filled pore space to volumetric water */
      soilwater = soilwater * (1.0 - (bulkdensity / PARTDENS));
      soilwater *= 100.0;

      /* CH4 oxidation for a deciduous system */
      if (*isdecid) {
        CH4max = 40.0 - 18.3 * bulkdensity;
        temp_adjust = 0.0209 * soiltemp + 0.845;
        /* Use bounded value for wfps_adjust if wfps falls below a critical */
        /* value, cak - 11/12/02 */
        if (wfps <= 0.05) {
          wfps_adjust = 0.1;
        } else {
          wfps_adjust = pow((10.0 * wfps - 0.5) / (1.84 - 0.5), 0.13);
          wfps_adjust *= pow((10.0 * wfps - 55) / (1.84 - 55),
                             (0.13 * (55 - 1.84)) / (1.84 - 0.5));
          wfps_adjust = max(0.1, wfps_adjust);
        }
        *CH4 = CH4max * wfps_adjust * temp_adjust;

      } else {
        /* CH4 oxidation for a grassland/coniferous/tropical system */
        Wmin = 3.0 * fieldcapacity - 0.28;
        Wopt = 6.3 * fieldcapacity - 0.58;
        Wmax = 10.6 * fieldcapacity + 1.9;
        temp = Wopt * 0.1 / (1.0 - (bulkdensity / PARTDENS));
        Dopt = diffusiv(&fieldcapacity, &bulkdensity, &temp);
        CH4max = 53.8 * Dopt + 0.58;
        if ((0.1*soilwater < Wmin) || 0.1*soilwater > Wmax) {
          watr_adjust = 0.1;
        } else {
          watr_adjust = pow(((0.1 * soilwater - Wmin) / (Wopt - Wmin)), 0.4) *
                        pow(((0.1 * soilwater - Wmax) / (Wopt - Wmax)), 
                             ((0.4 * (Wmax - Wopt)) / (Wopt - Wmin))); 
          watr_adjust = max(0.1, watr_adjust);
        }
        if (*isagri) {
          if (Dopt < 0.15) {
            agri_adjust = 0.9;
          } else if (Dopt > 0.28) {
            agri_adjust = 0.28;
          } else {
            agri_adjust = -4.6 * Dopt + 1.6;
          }
        } else {
          agri_adjust = 1.0;
        }
        temp_adjust = (soiltemp * max(0.11, Dopt) * 0.095) + 0.9;
        *CH4 = CH4max * watr_adjust * temp_adjust * agri_adjust;
      }

      return;
    }