/*************************************************************************
  Program:     soilresp.c
  --------
  Description:
  ----------- 
	Output soil respiration.
***************************************************************************
  CCRS (EMS/Applications Division)
  Written by:   J. Liu      
  Last update:	May 1998

*****************************************************************************/	
#include "beps.h" 
 
void soilresp(b,g,x,z) 
    double b[], g[], x[], z[];
{  
    /* output */
    double soilresp;	/* soil respiration in umol/m2/s */

    double thita;		/* average water content for root zone */
    double temp_soil;	/* soil temperture at the depth of 5 cm */

    double a1;			/* one-half field capacity */
    double a2;			/* one-half saturation */
    double a3;			/* optiaml respiration rate at 10 C */
    double a4=2.0;    	/* temperature sensitivity parameter */

    int short lc_p;		/* landcover */

    a1=0.5*b[35];
    a2=0.5*b[36];
    thita=x[2]/z[22];
    temp_soil=z[7];

    lc_p=(int)z[23];
    switch(lc_p)
    {
	case 7: case 8: case 9: case 10:  /*mixe*/ 
		a3=4.1;
		break;
	case 6: /*deci*/ 
		a3=3.5;
		break;
	case 1: case 2: case 3: case 4: case 5:   /*tran coni*/ 
		a3=4.1;
		break;
	case 11: case 12:	/*tund*/ 
		a3=1.7;
		break;
	case 23: case 24: case 25: case 26: case 27: case 28:	/*agri*/ 
		a3=2.5;
		break;
	case 16:	/*gras*/ 
		a3=1.7;
	default:
		a3=1.7;
    }

    /* Ref: Bonan 1995, Journal of Geophysical Research 100:2817-2831 */
    a1=thita/(a1+thita);
    a2=a2/(a2+thita);
    a4=1000/95*86.5*exp(308.56*(1/(1.4872*10.2+46.02)-1/(temp_soil+46.02)))/31536000000;
    a3=a3/4.1;
    //a4=pow(a4, (temp_soil-10)/10);  // J. Liu's
	// Model comparison: a4=pow(a4, (temp_soil-20));
    
	/* for hourly mode, mg C/m2/s) */
    //soilresp=a1*a2*a3*a4;  // J.Liu's
    soilresp=a3*a4;
//	g[42]=RPERIOD*soilresp*12/1000000000;   //J.Liu's 
//	g[42]=soilresp*1000; 
	g[42]=soilresp*RPERIOD*24; 

/*	printf("res_soil in mol/s: %10.8f",soilresp*1000);
    printf("res_soil: %10.8f",g[42]*1000); */

    return;
} 
