/* s_c_frac.aml
/* Script for producing three VIC layers of Sand, and three layers
/* of clay fractions, corresponding to the layers defined for Ksat.
/* The original grid of soil fracitons was obtained from Penn State, 
/* and is based on processed STATSGO soils data. 
/*
/* this is run with &run s_c_frac from the arc prompt.
/* The required input file that must be in the active workspace are:
/*   mo_soilfract1
/* mo_soilfract1 contains columns MUID, Rockdepm, Sand_L1, Silt_L1,
/* Clay_L1, Sand_l2,..., Clay_L11, Sand_L11

/* First, remove old version of the output files

&if [exists sand_frac_l1 -grid] &then
   kill sand_frac_l1 all
&if [exists sand_frac_l2 -grid] &then
   kill sand_frac_l2 all
&if [exists sand_frac_l3 -grid] &then
   kill sand_frac_l3 all
&if [exists clay_frac_l1 -grid] &then
   kill clay_frac_l1 all
&if [exists clay_frac_l2 -grid] &then
   kill clay_frac_l2 all
&if [exists clay_frac_l3 -grid] &then
   kill clay_frac_l3 all

grid

/* The values in the original database are in percent, so divide by 100
/* Again, since layers 10 and 11 have no data for the Missouri Basin,
/* they are excluded from these calculations

s_pct_l1 = mean ( MO_SOILFRACT1.SAND_L1, MO_SOILFRACT1.SAND_L2 )
s_pct_l2 = mean ( MO_SOILFRACT1.SAND_L3, MO_SOILFRACT1.SAND_L4, ~
  MO_SOILFRACT1.SAND_L5, MO_SOILFRACT1.SAND_L6, MO_SOILFRACT1.SAND_L7)
s_pct_l3 = mean ( MO_SOILFRACT1.SAND_L8, MO_SOILFRACT1.SAND_L9 )

c_pct_l1 = mean ( MO_SOILFRACT1.CLAY_L1, MO_SOILFRACT1.CLAY_L2 )
c_pct_l2 = mean ( MO_SOILFRACT1.CLAY_L3, MO_SOILFRACT1.CLAY_L4, ~
  MO_SOILFRACT1.CLAY_L5, MO_SOILFRACT1.CLAY_L6, MO_SOILFRACT1.CLAY_L7)
c_pct_l3 = mean ( MO_SOILFRACT1.CLAY_L8, MO_SOILFRACT1.CLAY_L9 )

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

sand_pct_l1 = aggregate (s_pct_l1, 15, mean, expand, data)
sand_pct_l2 = aggregate (s_pct_l2, 15, mean, expand, data)
sand_pct_l3 = aggregate (s_pct_l3, 15, mean, expand, data)

kill s_pct_l1 all
kill s_pct_l2 all
kill s_pct_l3 all

clay_pct_l1 = aggregate (c_pct_l1, 15, mean, expand, data)
clay_pct_l2 = aggregate (c_pct_l2, 15, mean, expand, data)
clay_pct_l3 = aggregate (c_pct_l3, 15, mean, expand, data)

kill c_pct_l1 all
kill c_pct_l2 all
kill c_pct_l3 all

/* change percent to fraction

sand_frac_l1 = (sand_pct_l1 * 0.01)
kill sand_pct_l1 all

/* sand_frac_l1 is the VIC sand fraction grid for layer 1.

/* display result for debugging
display 9999
mape sand_frac_l1
grids sand_frac_l1

/* now continue with sand layers 2 and 3, and clay layers 1, 2, 3

sand_frac_l2 = (sand_pct_l2 * 0.01)
sand_frac_l3 = (sand_pct_l3 * 0.01)
clay_frac_l1 = (clay_pct_l1 * 0.01)
clay_frac_l2 = (clay_pct_l2 * 0.01)
clay_frac_l3 = (clay_pct_l3 * 0.01)

kill sand_pct_l2 all
kill sand_pct_l3 all
kill clay_pct_l1 all
kill clay_pct_l2 all
kill clay_pct_l3 all