/* poros.aml
/* Script for producing three VIC layers of porosity values.
/* The original grid of soil porosity was obtained from Penn State, 
/* and is based on processed STATSGO soils data with values for 11 layers. 
/* The Penn State/STATSGO data is in fraction, which is the same as needed
/* by the VIC input
/*
/* this is run with &run poros from the arc prompt.
/* The required input file that must be in the active workspace are:
/*     mo_poros
/* mo_poros contains columns for L1_POROS, L2_POROS,...,L11_POROS with the
/* data.

/* First, remove old version of the output files

&if [exists poros_l1 -grid] &then
   kill poros_l1 all
&if [exists poros_l2 -grid] &then
   kill poros_l2 all
&if [exists poros_l3 -grid] &then
   kill poros_l3 all

grid

/* Again, since layers 10 and 11 have no data for the Missouri Basin,
/* they are excluded from these calculations

poros_l1_pix = mean ( MO_POROS.L1_POROS, MO_POROS.L2_POROS )
poros_l2_pix = mean ( MO_POROS.L3_POROS, MO_POROS.L4_POROS, ~
  MO_POROS.L5_POROS, MO_POROS.L6_POROS, MO_POROS.L7_POROS)
poros_l3_pix = mean ( MO_POROS.L8_POROS, MO_POROS.L9_POROS )

/* aggregate from 30 arc sec. to cell size of 1/8 degree

poros_l1 = aggregate (poros_l1_pix, 15, mean, expand, data)
poros_l2 = aggregate (poros_l2_pix, 15, mean, expand, data)
poros_l3 = aggregate (poros_l3_pix, 15, mean, expand, data)

kill poros_l1_pix all
kill poros_l2_pix all
kill poros_l3_pix all

/* poros_l1 is the VIC porosity grid for layer 1.

/* display result for debugging
display 9999
mape poros_l1