tgmodel.c File Reference

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

Include dependency graph for tgmodel.c:

Go to the source code of this file.

Functions
void	trace_gas_model (int *jday, double *time, double *newminrl, double *ammonium, double nitrate[], int *texture, double *sand, double *silt, double *clay, double *afiel, double *bulkd, double *maxt, double *ppt, double *snow, double *avgwfps, double *stormf, double *basef, double frlechd[], double stream[], double *inorglch, double *critflow, double wfluxout[], double *newCO2, double *NOflux, double *Nn2oflux, double *Dn2oflux, double *Dn2flux, double *CH4, int *isdecid, int *isagri, double *grass_lai, double *tree_lai, double *NOabsorp_grass, double *NOabsorp_tree, double *nit_amt, double *nreduce, double dN2lyr[], double dN2Olyr[], SITEPAR_SPT sitepar, LAYERPAR_SPT layers, SOIL_SPT soil)

Function Documentation

void trace_gas_model	(	int *	jday,
		double *	time,
		double *	newminrl,
		double *	ammonium,
		double	nitrate[],
		int *	texture,
		double *	sand,
		double *	silt,
		double *	clay,
		double *	afiel,
		double *	bulkd,
		double *	maxt,
		double *	ppt,
		double *	snow,
		double *	avgwfps,
		double *	stormf,
		double *	basef,
		double	frlechd[],
		double	stream[],
		double *	inorglch,
		double *	critflow,
		double	wfluxout[],
		double *	newCO2,
		double *	NOflux,
		double *	Nn2oflux,
		double *	Dn2oflux,
		double *	Dn2flux,
		double *	CH4,
		int *	isdecid,
		int *	isagri,
		double *	grass_lai,
		double *	tree_lai,
		double *	NOabsorp_grass,
		double *	NOabsorp_tree,
		double *	nit_amt,
		double *	nreduce,
		double	dN2lyr[],
		double	dN2Olyr[],
		SITEPAR_SPT	sitepar,
		LAYERPAR_SPT	layers,
		SOIL_SPT	soil
	)

    {

      /* Local Variables */

      int    debug = 0;
      int    ilyr;
      int *soilclassification;
/*      double netmn_to_no3 = 0.20; */
      double netmn_to_no3 = 0.0;
      double turnovfrac = 0.02;
      double newNH4;
      double newNO3;
      double nh4_to_no3;
      double krainNO;
      double potential_NOflux;
      double dDO;
      double stdbulkd;
      double stdfieldc;
      double NO_N2O_ratio;
      double NH4_to_NO;
      double npool_sum;
      double canopy_reduction;
      double NOabsorp;
      double total_lai;

//      extern SITEPAR_SPT sitepar;

      *Nn2oflux = 0.0;
      *NOflux = 0.0;
      *Dn2oflux = 0.0;
      *Dn2flux = 0.0;

      /* Compute fraction of new mineralization that is converted to NH4 */
      /* and NO3 */

      if (debug) {
        printf("newminrl = %6.4lf\n", *newminrl);
      }

      if (*newminrl <= 0.0) {

        /* Immobilization */
        /* Distribute N loss proportionally between ammonium and nitrate   */
        /* layers.  There is no check that these N pools won't go negative */
        /* once immobilization is accounted for.  It is assumed that the   */
        /* immobilization calculation by the decomp model is moderated by  */
        /* the supply of minerl N.                                         */

        npool_sum = (*ammonium > 0.0) ? *ammonium : 0.0;
        for (ilyr=0; ilyr < MAXLYR; ilyr ++) {
          npool_sum += (nitrate[ilyr] > 0.0) ? nitrate[ilyr] : 0.0;
        }
        if (*ammonium > 0.0) {
          *ammonium += *newminrl * (*ammonium / npool_sum);
        }
        for (ilyr=0; ilyr < MAXLYR; ilyr ++) {
          if (nitrate[ilyr] > 0.0) {
            nitrate[ilyr] += *newminrl * (nitrate[ilyr] / npool_sum);
          }
        }
        newNH4 = 0.0;
        newNO3 = 0.0;
      } else {
        /* Mineralization */
        newNH4 = *newminrl * (1.0 - netmn_to_no3);
        newNO3 = *newminrl * netmn_to_no3;
      }

      if (debug) {
        printf("newNH4 = %6.4lf\n", newNH4);
        printf("newNO3 = %6.4lf\n", newNO3);
      }

      *ammonium += newNH4;

      /* Compute the amount of NH4 that is converted to NO3 due to */
      /* nitrification */

     nitrify(ammonium, &nh4_to_no3, maxt, nreduce,sitepar,layers,soil);
      *nit_amt = nh4_to_no3;

      if (debug) {
        printf("texture = %1d\n", *texture);
        printf("nh4_to_no3 = %6.4lf\n", nh4_to_no3);
        printf("maxt = %6.4lf\n", *maxt);
      }

      /* Compute fraction of new NO3 that is converted to N2O and NO */

      krainNO = nox_pulse(ppt, snow);

      getsoilprop(sand, silt, clay, &stdbulkd, &stdfieldc, texture,&soilclassification);

      /* Use standard field capacity and bulk density according */
      /* to the soil class in the texture triangle -mdh 10/26/99 */
/*      dDO = diffusiv(afiel(1), bulkd, *avgwfps) */
      /* No, change back to soils.in field capacity and bulk density. */
      /* -mdh 6/20/00 */
/*      dDO = diffusiv(&stdfieldc, &stdbulkd, avgwfps); */
      dDO = diffusiv(afiel, bulkd, avgwfps);

      newNO3 += nh4_to_no3;

      if (newNO3 > 1.0E-30) {
        *Nn2oflux = newNO3 * turnovfrac * sitepar->N2Oadjust;
        newNO3 -= *Nn2oflux; 

        /* Another update to NO flux calculation -mdh 10/26/99 */

/*        NO_N2O_ratio = 15.23 + (35.45*atan(0.676*PI*(10*dDO-1.86)))/PI; */
        NO_N2O_ratio = 8.0 + (18.0*atan(0.75*PI*(10*dDO-1.86)))/PI;
        /* If this is an agricultural system adjust the NO to N2O ratio */
        /* cak - 01/28/03 */
        if (*isagri) {
/*          NO_N2O_ratio *= 0.2; */
          NO_N2O_ratio *= 0.5;
        }
        potential_NOflux = NO_N2O_ratio * *Nn2oflux * krainNO;

        if (potential_NOflux <= newNO3) {
          *NOflux = potential_NOflux;
          newNO3 -= *NOflux;
        } else {
          /* take N out of ammonimum to get max NOflux possible */
          NH4_to_NO = min(*ammonium, (potential_NOflux-newNO3));
          *NOflux = newNO3 + NH4_to_NO;
          *ammonium -= NH4_to_NO;
          newNO3 = 0;
        }

        if (*NOflux < 1.0E-30) {
          *NOflux = 0.0;
        }

      } else {
        NO_N2O_ratio = 0.0;
      }

      /* Compute the N2O flux (Dn2oflux) and N2 flux (Dn2flux) due to */
      /* denitrification */

denitrify(newCO2, &newNO3, nitrate, wfluxout, critflow, frlechd,
                stream, basef, stormf, inorglch, Dn2oflux, Dn2flux,
                stdfieldc, stdbulkd, dN2lyr, dN2Olyr, jday,sitepar,layers);


      /* Now compute NOflux from denitrification (new calculation */
      /* -mdh 6/1/00 */
/*      potential_NOflux = NO_N2O_ratio * *Dn2oflux * krainNO; */
      /* For denitrification, krainNO is >= 1.0 -mdh 6/22/00 */

      potential_NOflux = NO_N2O_ratio * *Dn2oflux * min(1.0, krainNO);

      if (potential_NOflux <= *ammonium) {
        /* Take all N out of ammonimum pool */
        *NOflux += potential_NOflux;
        *ammonium -= potential_NOflux;
      } else {
        /* Take N out of available ammonium, then convert some Dn2oflux to */
        /* NOflux */
        *NOflux += *ammonium;
        potential_NOflux -= *ammonium;
        *ammonium = 0.0;
        if (potential_NOflux <= *Dn2oflux) {
          *NOflux += potential_NOflux;
          *Dn2oflux -= potential_NOflux;
        }
      }

      /* Compute the amount of the soil NO flux that is absorped by the canopy, cak - 09/23/03 */
      total_lai = *grass_lai + *tree_lai;
      if (total_lai > 0.0) {
        canopy_reduction = 0.0077 * pow(total_lai,2) + -0.13 * total_lai + 0.99;
        /* We need to retain the soil flux value */

        /* This previous code did not maintain and N balance -Melannie 9/26/2011 */
        /* NOsoil = *NOflux; */
        /* *NOflux *= canopy_reduction; */
        /* NOabsorp = NOsoil - *NOflux; */

        NOabsorp = *NOflux * (1 - canopy_reduction);
        if (NOabsorp > 0.0)
        {
            *NOabsorp_grass = NOabsorp * (*grass_lai / total_lai);
            *NOabsorp_tree = NOabsorp * (*tree_lai / total_lai);
            *NOflux -= NOabsorp;
        }

        /* NO absorped by canopy goes to crop storage and forest storage */
        /* Reset NOflux using the retained value */
        /* *NOflux = NOsoil; */
      }

      if (*NOflux < 1.0E-30) {
        *NOflux = 0.0;
      }
      if (*Nn2oflux < 1.0E-30) {
        *Nn2oflux = 0.0;
      }
      if (*Dn2oflux < 1.0E-30) {
        *Dn2oflux = 0.0;
      }
      if (*Dn2flux < 1.0E-30) {
        *Dn2flux = 0.0;
      }
      
      /* Calculate methane oxidation */
      methane_oxidation(CH4, isdecid, isagri,layers,soil);
  //    wrtsoiln(time, jday,ammonium,nitrate);
      return;
    }