initialize_atmos.c

超强的大尺度水文模拟工具

    }
  }

  /**************************************************
    calculate the subdaily and daily vapor pressure 
    and vapor pressure deficit
  **************************************************/

  if(!param_set.TYPE[VP].SUPPLIED) {
    for (rec = 0; rec < global_param.nrecs; rec++) {
      atmos[rec].vp[NR] = vp[rec/stepspday];
      atmos[rec].vpd[NR] = svp(atmos[rec].air_temp[NR]) - atmos[rec].vp[NR];

      if(atmos[rec].vpd[NR]<0) {
	atmos[rec].vpd[NR]=0;
	atmos[rec].vp[NR]=svp(atmos[rec].air_temp[NR]);
      }
      
      for (i = 0; i < NF; i++) {
	atmos[rec].vp[i]  = atmos[rec].vp[NR];
	atmos[rec].vpd[i]  = (svp(atmos[rec].air_temp[i]) - atmos[rec].vp[i]);

	if(atmos[rec].vpd[i]<0) {
	  atmos[rec].vpd[i]=0;
	  atmos[rec].vp[i]=svp(atmos[rec].air_temp[i]);
	}
      
      }
    }
  }
  else {
    if(param_set.FORCE_DT[param_set.TYPE[VP].SUPPLIED-1] == 24) {
      /* daily vp provided */
      rec = 0;
      for (day = 0; day < Ndays; day++) {
	for (i = 0; i < stepspday; i++) {
	  sum = 0;
	  for (j = 0; j < NF; j++) {
	    atmos[rec].vp[j] = forcing_data[VP][day];
	    atmos[rec].vpd[j] = (svp(atmos[rec].air_temp[j]) 
				 - atmos[rec].vp[j]);
	    sum += atmos[rec].vp[j];
	  }
	  if(NF > 1) {
	    atmos[rec].vp[NR] = sum / (float)NF;
	    atmos[rec].vpd[NR] = (svp(atmos[rec].air_temp[NR]) 
				  - atmos[rec].vp[NR]);
	    rec++;
	  }
	}
      }
    }
    else {
      /* sub-daily vp provided */
      idx = 0;
      for(rec = 0; rec < global_param.nrecs; rec++) {
	sum = 0;
	for(i = 0; i < NF; i++) {
	  atmos[rec].vp[i] = forcing_data[VP][idx];


	  atmos[rec].vp[i] = ((float)(int)(atmos[rec].vp[i]*1000 + 0.5)/1000);


	  atmos[rec].vpd[i] = (svp(atmos[rec].air_temp[i]) 
			       - atmos[rec].vp[i]);
	  sum += atmos[rec].vp[i];
	  idx++;
	}
	if(NF > 1) {
	  atmos[rec].vp[NR] = sum / (float)NF;
	  atmos[rec].vpd[NR] = (svp(atmos[rec].air_temp[NR]) 
				    - atmos[rec].vp[NR]);
	}
      }
    }
  }

  /****************************************************************************
    calculate the daily and sub-daily longwave.  There is a separate case for
    the full energy and the water balance modes.  For water balance mode we 
    need to calculate the net longwave for the daily timestep and the incoming
    longwave for the SNOW_STEPs, for the full energy balance mode we always
    want the incoming longwave. 
  ****************************************************************************/

  if ( !param_set.TYPE[LONGWAVE].SUPPLIED ) {
    /** Incoming longwave radiation not supplied **/
    for (rec = 0; rec < global_param.nrecs; rec++) {
      if( NF > 1 ) {
	for (i = 0; i < NF; i++) {
	  calc_longwave(&(atmos[rec].longwave[i]), tskc[rec/stepspday],
			atmos[rec].air_temp[i], atmos[rec].vp[i]);
	}
	calc_netlongwave(&(atmos[rec].longwave[NR]), tskc[rec/stepspday],
			 atmos[rec].air_temp[NR], atmos[rec].vp[NR]);
      }
      else {
	calc_longwave(&(atmos[rec].longwave[NR]), tskc[rec/stepspday],
		      atmos[rec].air_temp[NR], atmos[rec].vp[NR]);
      }
    }
  }
  else if(param_set.FORCE_DT[param_set.TYPE[LONGWAVE].SUPPLIED-1] == 24) {
    /* daily incoming longwave radiation provided */
    rec = 0;
    for (day = 0; day < Ndays; day++) {
      for (i = 0; i < stepspday; i++) {
	sum = 0;
	for (j = 0; j < NF; j++) {
	  atmos[rec].longwave[j] = forcing_data[LONGWAVE][day];
	  sum += atmos[rec].longwave[j];
	}
	if(NF>1) atmos[rec].longwave[NR] = sum / (float)NF - STEFAN_B 
		   * atmos[rec].air_temp[NR] * atmos[rec].air_temp[NR] 
		   * atmos[rec].air_temp[NR] * atmos[rec].air_temp[NR];
	rec++;
      }
    }
  }
  else {
    /* sub-daily incoming longwave radiation provided */
    idx = 0;
    for(rec = 0; rec < global_param.nrecs; rec++) {
      sum = 0;
      for(i = 0; i < NF; i++) {
	atmos[rec].longwave[i] = forcing_data[LONGWAVE][idx];
	sum += atmos[rec].longwave[i];
	idx++;
      }
      if(NF>1) atmos[rec].longwave[NR] = sum / (float)NF - STEFAN_B 
		 * atmos[rec].air_temp[NR] * atmos[rec].air_temp[NR] 
		 * atmos[rec].air_temp[NR] * atmos[rec].air_temp[NR];
    }
  }

  /********************
    set the windspeed 
  ********************/

  if (!param_set.TYPE[WIND].SUPPLIED) {
    /* no wind data provided, use default constant */
    for (rec = 0; rec < global_param.nrecs; rec++) {
      for (i = 0; i < NF; i++) {
	atmos[rec].wind[i] = 1.5;
      }
      atmos[rec].wind[NR] = 1.5;	
    }
  }
  else {
    if(param_set.FORCE_DT[param_set.TYPE[WIND].SUPPLIED-1] == 24) {
      /* daily wind provided */
      rec = 0;
      for (day = 0; day < Ndays; day++) {
	for (i = 0; i < stepspday; i++) {
	  sum = 0;
	  for (j = 0; j < NF; j++) {
	    if(forcing_data[WIND][day] < options.MIN_WIND_SPEED)
	      atmos[rec].wind[j] = options.MIN_WIND_SPEED;
	    else 
	      atmos[rec].wind[j] = forcing_data[WIND][day];
	    sum += atmos[rec].wind[j];
	  }
	  if(NF>1) atmos[rec].wind[NR] = sum / (float)NF;
	  if(global_param.dt == 24) {
	    if(forcing_data[WIND][day] < options.MIN_WIND_SPEED)
	      atmos[rec].wind[j] = options.MIN_WIND_SPEED;
	    else 
	      atmos[rec].wind[NR] = forcing_data[WIND][day];
	  }
	  rec++;
	}
      }
    }
    else {
      /* sub-daily wind speed provided */
      idx = 0;
      for(rec = 0; rec < global_param.nrecs; rec++) {
	sum = 0;
	for(i = 0; i < NF; i++) {
	  if(forcing_data[WIND][idx] <options.MIN_WIND_SPEED)
	    atmos[rec].wind[i] = options.MIN_WIND_SPEED;
	  else
	    atmos[rec].wind[i] = forcing_data[WIND][idx];
	  sum += atmos[rec].wind[i];
	  idx++;
	}
	if(NF>1) atmos[rec].wind[NR] = sum / (float)NF;
      }
    }
  }

  /*************************************************
    Store atmospheric density if provided (kg/m^3)
  *************************************************/

  if(param_set.TYPE[DENSITY].SUPPLIED) {
    if(param_set.FORCE_DT[param_set.TYPE[DENSITY].SUPPLIED-1] == 24) {
      /* daily density provided */
      rec = 0;
      for (day = 0; day < Ndays; day++) {
	for (i = 0; i < stepspday; i++) {
	  sum = 0;
	  for (j = 0; j < NF; j++) {
	    atmos[rec].density[j] = forcing_data[DENSITY][day];
	    sum += atmos[rec].density[j];
	  }
	  if(NF>1) atmos[rec].density[NR] = sum / (float)NF;
	  rec++;
	}
      }
    }
    else {
      /* sub-daily density provided */
      idx = 0;
      for(rec = 0; rec < global_param.nrecs; rec++) {
	sum = 0;
	for(i = 0; i < NF; i++) {
	  atmos[rec].density[i] = forcing_data[DENSITY][idx];
	  sum += atmos[rec].density[i];
	  idx++;
	}
	if(NF>1) atmos[rec].density[NR] = sum / (float)NF;
      }
    }
  }

  /**************************************
    Estimate Atmospheric Pressure (kPa) 
  **************************************/

  if(!param_set.TYPE[PRESSURE].SUPPLIED) {
    if(!param_set.TYPE[DENSITY].SUPPLIED) {
      /* set pressure to constant value */
      for (rec = 0; rec < global_param.nrecs; rec++) {
	atmos[rec].pressure[NR] = 95.5;
	for (i = 0; i < NF; i++) {
	  atmos[rec].pressure[i] = atmos[rec].pressure[NR];
	}
      }
    }
    else {
      /* use observed densities to estimate pressure */
      for (rec = 0; rec < global_param.nrecs; rec++) {
	atmos[rec].pressure[NR] = (275.0 + atmos[rec].air_temp[NR])
	  *atmos[rec].density[NR]/3.486;
	for (i = 0; i < NF; i++) {
	  atmos[rec].pressure[i] = (275.0 + atmos[rec].air_temp[i])
	    *atmos[rec].density[i]/3.486;
	}
      }
    }
  }
  else {
    /* observed atmospheric pressure supplied */
    if(param_set.FORCE_DT[param_set.TYPE[PRESSURE].SUPPLIED-1] == 24) {
      /* daily pressure provided */
      rec = 0;
      for (day = 0; day < Ndays; day++) {
	for (i = 0; i < stepspday; i++) {
	  sum = 0;
	  for (j = 0; j < NF; j++) {
	    atmos[rec].pressure[j] = forcing_data[PRESSURE][day];
	    sum += atmos[rec].pressure[j];
	  }
	  if(NF>1) atmos[rec].pressure[NR] = sum / (float)NF;
	  rec++;
	}
      }
    }
    else {
      /* sub-daily pressure provided */
      idx = 0;
      for(rec = 0; rec < global_param.nrecs; rec++) {
	sum = 0;
	for(i = 0; i < NF; i++) {
	  atmos[rec].pressure[i] = forcing_data[PRESSURE][idx];
	  sum += atmos[rec].pressure[i];
	  idx++;
	}
	if(NF>1) atmos[rec].pressure[NR] = sum / (float)NF;
      }
    }
  }

  /********************************************************
    Estimate Atmospheric Density if not provided (kg/m^3)
  ********************************************************/

  if(!param_set.TYPE[DENSITY].SUPPLIED) {
    for (rec = 0; rec < global_param.nrecs; rec++) {
      atmos[rec].density[NR] = 3.486*atmos[rec].pressure[NR]/
	(275.0 + atmos[rec].air_temp[NR]);
      for (i = 0; i < NF; i++) {
	atmos[rec].density[i] = 3.486*atmos[rec].pressure[i]/
	  (275.0 + atmos[rec].air_temp[i]);
      }
    }
  }

  /****************************************************
    Determine if Snow will Fall During Each Time Step
  ****************************************************/

#if !OUTPUT_FORCE
  min_Tfactor = Tfactor[0];
  for (band = 1; band < options.SNOW_BAND; band++) {
    if (Tfactor[band] < min_Tfactor)
      min_Tfactor = Tfactor[band];
  }
  for (rec = 0; rec < global_param.nrecs; rec++) {
    atmos[rec].snowflag[NR] = FALSE;
    for (i = 0; i < NF; i++) {
      if ((atmos[rec].air_temp[i] + min_Tfactor) < global_param.MAX_SNOW_TEMP
	  &&  atmos[rec].prec[i] > 0) {
	atmos[rec].snowflag[i] = TRUE;
	atmos[rec].snowflag[NR] = TRUE;
      }
      else
	atmos[rec].snowflag[i] = FALSE;
    }
  }
#endif
 
  // Free temporary parameters
  free(hourlyrad);
  free(prec);
  free(tair);
  free(tmax);
  free(tmaxhour);
  free(tmin);
  free(tminhour);
  free(tskc);
  free(vp);

  for(i=0;i<N_FORCING_TYPES;i++) 
    if (param_set.TYPE[i].SUPPLIED) 
      free((char *)forcing_data[i]);
  free((char *)forcing_data);

  // If COMPUTE_TREELINE is TRUE and the treeline computation hasn't
  // specifically been turned off for this cell (by supplying avgJulyAirTemp
  // and setting it to -999), calculate which snowbands are above the
  // treeline, based on average July air temperature.
  if (options.COMPUTE_TREELINE) {
    if ( !(options.JULY_TAVG_SUPPLIED && avgJulyAirTemp == -999) ) {
      if ( options.SNOW_BAND ) {
        compute_treeline( atmos, dmy, avgJulyAirTemp, Tfactor, AboveTreeLine );
      }
    }
  }

#if OUTPUT_FORCE
  // If OUTPUT_FORCE is set to TRUE in user_def.h then the full
  // forcing data array is dumped into a new set of files.
  write_forcing_file(atmos, global_param.nrecs, outfiles);
#endif /* OUTPUT_FORCE */

}