calc_root_fraction.c

植被分类和植被覆盖度计算程序

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

static char vcid[] = "$Id: calc_root_fraction.c,v 4.1 2000/05/16 21:07:16 vicadmin Exp $";

void calc_root_fractions(veg_con_struct  *veg_con,
			 soil_con_struct  *soil_con)
/**********************************************************************
  calc_root_fraction.c    Keith Cherkauer      September 24, 1998

  This routine computes the fraction of roots in each soil layer based
  on the root zone distribution defined in the vegetation parameter
  file.  Roots are assumed to be linearly distributed within each
  root zone.

**********************************************************************/
{
  extern option_struct options;

  char   ErrStr[MAXSTRING];
  int    Nveg;
  int    veg;
  int    layer;
  int    zone;
  int    i;
  float  sum_depth;
  float  sum_fract;
  float  dum;
  double Zstep;
  double Zsum;
  double Lstep;
  double Lsum;
  double Zmin_fract;
  double Zmin_depth;
  double Zmax;

  Nveg      = veg_con[0].vegetat_type_num;

  for(veg=0;veg<Nveg;veg++) {
    sum_depth  = 0;
    sum_fract  = 0;
    layer      = 0;
    Lstep      = soil_con->depth[layer];
    Lsum       = Lstep;
    Zsum       = 0;
    zone       = 0;
    
    while(zone<options.ROOT_ZONES) {
      Zstep = (double)veg_con[veg].zone_depth[zone];
      if((Zsum + Zstep) <= Lsum && Zsum >= Lsum - Lstep) {
	/** CASE 1: Root Zone Completely in Soil Layer **/
	sum_fract += veg_con[veg].zone_fract[zone];
      }
      else {
	/** CASE 2: Root Zone Partially in Soil Layer **/
	if(Zsum < Lsum - Lstep) {
	  /** Root zone starts in previous soil layer **/
	  Zmin_depth = Lsum - Lstep;
	  Zmin_fract = linear_interp(Zmin_depth,Zsum,Zsum+Zstep,0,
				     veg_con[veg].zone_fract[zone]);
	}
	else {
	  /** Root zone starts in current soil layer **/
	  Zmin_depth = Zsum;
	  Zmin_fract = 0.;
	}
	if(Zsum + Zstep <= Lsum) {
	  /** Root zone ends in current layer **/
	  Zmax = Zsum + Zstep;
	}
	else {
	  /** Root zone extends beyond bottom of current layer **/
	  Zmax = Lsum;
	}
	sum_fract += linear_interp(Zmax,Zsum,Zsum+Zstep,0,
				   veg_con[veg].zone_fract[zone]) - Zmin_fract;
      }

      /** Update Root Zone and Soil Layer **/
      if(Zsum + Zstep < Lsum) {
	Zsum += Zstep;
	zone ++;
      }
      else if(Zsum + Zstep == Lsum) {
	Zsum += Zstep;
	zone ++;
	if(layer<options.Nlayer) {
	  veg_con[veg].root[layer] = sum_fract;
	  sum_fract = 0.;
	}
	layer++;
	if(layer<options.Nlayer) {
	  Lstep  = soil_con->depth[layer];
	  Lsum  += Lstep;
	}
	else if(layer==options.Nlayer) {
	  Lstep  = Zsum + Zstep - Lsum;
	  if(zone<options.ROOT_ZONES-1) {
	    for(i=zone+1;i<options.ROOT_ZONES;i++) {
	      Lstep += veg_con[veg].zone_depth[i];
	    }
	  }
	  Lsum  += Lstep;
	}
      }
      else if(Zsum + Zstep > Lsum) {
	if(layer<options.Nlayer) {
	  veg_con[veg].root[layer] = sum_fract;
	  sum_fract = 0.;
	}
	layer++;
	if(layer<options.Nlayer) {
	  Lstep  = soil_con->depth[layer];
	  Lsum  += Lstep;
	}
	else if(layer==options.Nlayer) {
	  Lstep  = Zsum + Zstep - Lsum;
	  if(zone<options.ROOT_ZONES-1) {
	    for(i=zone+1;i<options.ROOT_ZONES;i++) {
	      Lstep += veg_con[veg].zone_depth[i];
	    }
	  }
	  Lsum  += Lstep;
	}
      }
	
    }

    if(sum_fract > 0 && layer >= options.Nlayer) {
      veg_con[veg].root[options.Nlayer-1] += sum_fract;
    }
    else if(sum_fract > 0) {
      veg_con[veg].root[layer] += sum_fract;
    }

    dum=0.;
    for (layer=0;layer<options.Nlayer;layer++) {
      if(veg_con[veg].root[layer] < 1.e-4) veg_con[veg].root[layer] = 0.;
      dum += veg_con[veg].root[layer];
    }
    if(dum == 0.0){
      sprintf(ErrStr,"Root fractions sum equals zero: %f , Vege Class: %d\n",
	      dum, veg_con[veg].veg_class);
      nrerror(ErrStr);
    }
    for (layer=0;layer<options.Nlayer;layer++) {
      veg_con[veg].root[layer] /= dum;
    }

  }

}