diffusiv.c File Reference

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

Include dependency graph for diffusiv.c:

Go to the source code of this file.

Macros
#define	PARTDENS   2.65 /* Particle Density (g/cm3) */
#define	min(a, b)   (((a) < (b)) ? (a) : (b))
#define	max(a, b)   (((a) > (b)) ? (a) : (b))

Functions
double	diffusiv (double *A, double *bulkden, double *wfps)

Macro Definition Documentation

#define max	(		a,
			b
	)		(((a) > (b)) ? (a) : (b))

#define min	(		a,
			b
	)		(((a) < (b)) ? (a) : (b))

#define PARTDENS   2.65 /* Particle Density (g/cm3) */

Function Documentation

double diffusiv	(	double *	A,
		double *	bulkden,
		double *	wfps
	)

    {
      double porosity;
      double theta_A, theta_P, theta_V, my_theta_V;
      double tp1, tp2, tp3, tp4, tp5, tp6, tp7, tp8;
      double vfrac;
      double s_wat, sw_p;
      double pfc;
      double dDO;
      int   debug;

      /* Initialization */
      debug = 0;

      porosity = 1.0 - (*bulkden)/PARTDENS;
      if ((porosity - (*A)) < 0.0) {
        printf("ERROR:  Mis-match between bulk density and field capacity\n");
        printf("value in soils.in file.  Ending simulation!\n");
        exit(1);
      }

      pfc = (*wfps)*100 /((*A)/porosity);
 
      if (pfc >= 100.0) {
        vfrac = ((*wfps)*porosity - *A) / (porosity - (*A));
      } else {
        vfrac = 0.0;
      }
  
      if (pfc < 100.0) {
        theta_V = (pfc/100.0)*(*A);
      } else {
        theta_V = (*A) + (min(vfrac,1.0) * (porosity - (*A)));
      }

      my_theta_V = *wfps * (1.0 - *bulkden/PARTDENS);

      theta_P = (theta_V < (*A)) ? 0.0 : theta_V - (*A);
      theta_A = (theta_V > (*A)) ? (*A) : theta_V;

      if (debug) {
        printf("In SUB, vfrac = %f\n", vfrac);
        printf("In SUB, porosity - A = %f\n", (porosity - (*A)));
        printf("In SUB, (*A) = %f\n", (*A));
        printf("In SUB, (*wfps) = %f\n", (*wfps));
        printf("In SUB, pfc = %f\n", pfc);
        printf("In SUB, porosity = %f\n", porosity);
        printf("In SUB, theta_V = %f\n", theta_V);
        printf("In SUB, my_theta_V = %f\n", my_theta_V);
        printf("In SUB, theta_P = %f\n", theta_P);
        printf("In SUB, theta_A = %f\n", theta_A);
      }

      s_wat = min(1.0, theta_V/(*A));
      sw_p = min(1.0, theta_P/(porosity-(*A)));

      tp1 = pow((1.0 - s_wat), 2.0);
      tp2 = ((*A) - theta_A) / ((*A) + (1.0 - porosity));
      tp3 = pow(tp2,(0.5*tp2 + 1.16));
      tp4 = (1.0 - pow(((porosity-(*A))),
                              (0.5*((porosity-(*A))) + 1.16)));
      tp5 = (porosity-(*A))-theta_P;
      if (tp5 > 0.0) {
        tp6 = pow(tp5, (0.5*tp5 + 1.16));
      } else {
        tp6 = 0.0;
      }
      tp7 = pow((1.0-sw_p), 2.0);
      tp8 = max(0.0, (tp1*tp3*tp4*(tp5-tp6)) /
                      (1.0E-6 + (tp1*tp3*tp4) + tp5 - tp6) * 1.0E7);

      if (debug) printf("From SUB: %f %f %f %f %f %f %f %f\n",
                        tp1, tp2, tp3, tp4, tp5, tp6, tp7, tp8);
      dDO = max(0.0, (tp8/1.0E7 + tp7*tp6));
      if (debug) {
            printf("In SUB, dDO = %f\n", dDO);
      }

/*      dDO /= 0.20;
      dDO = min(dDO, 1.0); */

      return(dDO);
    }