comment:
This is an IPAT-S script. IPAT-S is an open-source scripting language
for developing sustainability scenarios. It can be run using the free
IPAT-S software, available from:

	http://ipat-s.kb-creative.net/

For more information, please contact the author of IPAT-S,
Eric Kemp-Benedict (eric@kb-creative.net)
:comment

comment:
This script implements land use change in a scenario by calculating
an area transition matrix (ATM). The scenario examines future changes
in US land area based on historical patterns and future demographics.
It uses the historical area transition matrix for 1982-1997 from the
US National Resource Inventory (1997 NRI, revised Dec. 2000,
http://www.nrcs.usda.gov/technical/NRI/1997/summary_report/table5.html).

The ATM is solved for in the scenario using a linear goal program (LGP).
The scenario ATM is set as close to an historically-based target ATM
as possible, while meeting the requirements for agricultural and
developed land.

Because a transition matrix approach is used, rather than net changes
in land area, it is possible to keep track of the loss of existing
forest as well as additions of new forest land separately. This is
useful because it can take several decades for new forest to achieve
the same biomass density of standing forest.
:comment

base year 2000
scenario years 2005 to 2015 by 5

# LUC is "Land Use Category"
ditto 'Cropland' 'Pasture & Range' 'Forest' 'Developed' 'Other':
dimension LUC ''
# For ATM, define "to" and "from" LUCs using same categories
dimension toLUC ''
dimension fromLUC ''

# Key scenario variables to be solved for
summvar Area{LUC} ATM{toLUC, fromLUC}

##
## Scenario assumptions
##
ratio AreaIntensity{LUC}
ratio driverGrowth{LUC}

# Three LUCs are set exogenously
AreaIntensity{LUC='Cropland'} = gr(<-0.7%>)
AreaIntensity{LUC='Pasture & Range'} = gr(<-0.3%>)
AreaIntensity{LUC='Developed'} = gr(<1.0%>)

# For most LUCs, there is no driver: just specify annual
#   % change in area (as change in AreaIntensity)
driverGrowth = 1
# Population drives developed land
# Growth rates from GEO-3 dataset (from 2000 revision of the UN mid-range
#   population projections)
driverGrowth{LUC='Developed'} = gr(<0.886%, 0.827%, 0.805%>)

##
## Base year data (from 1997 US National Resource Inventory, rev. Dec 2000)
##
# Data are for 1997, here assigned to 2000 base year
Area.0{LUC='Cropland'} = 376997.9
Area.0{LUC='Pasture & Range'} = 525969.1
Area.0{LUC='Forest'} = 406955.2
Area.0{LUC='Developed'} = 98251.7
Area.0{LUC='Other'} = 535955.8

##
## Target ATM values (transition matrix for 1982-1997 period)
##
summvar targetATM{toLUC, fromLUC} histATM{toLUC, fromLUC}
# First, load historical ATM, 1982-1997 as base
ditto fromLUC='Cropland':
histATM{'', toLUC='Pasture & Range'} = 22928.6
histATM{'', toLUC='Forest'} = 5606.5
histATM{'', toLUC='Developed'} = 7097.5
histATM{'', toLUC='Other'} = 35056.1
ditto fromLUC='Pasture & Range':
histATM{'', toLUC='Cropland'} = 22314.5
histATM{'', toLUC='Forest'} = 17113
histATM{'', toLUC='Developed'} = 7511.3
histATM{'', toLUC='Other'} = 9495.8
ditto fromLUC='Forest':
histATM{'', toLUC='Cropland'} = 2037.1
histATM{'', toLUC='Pasture & Range'} = 6267
histATM{'', toLUC='Developed'} = 10279.2
histATM{'', toLUC='Other'} = 4411.6
ditto fromLUC='Developed':
histATM{'', toLUC='Cropland'} = 196.7
histATM{'', toLUC='Pasture & Range'} = 189.4
histATM{'', toLUC='Forest'} = 227
histATM{'', toLUC='Other'} = 14
ditto fromLUC='Other':
histATM{'', toLUC='Cropland'} = 2184.3
histATM{'', toLUC='Pasture & Range'} = 4273.3
histATM{'', toLUC='Forest'} = 3665.4
histATM{'', toLUC='Developed'} = 745

# Next, perform simple extrapolation for baseline "targetATM"
:: histATM * (y - y.0)/(1997 - 1982) -> targetATM

##
## Scenario
##

# Areas that are set exogenously in the scenario:
summvar setArea{LUC} setFilter{LUC}
setFilter{LUC='Cropland'} = 1
setFilter{LUC='Pasture & Range'} = 1
setFilter{LUC='Developed'} = 1

# Set value for "setArea" in the base year
:: setFilter * Area -> setArea

# Scenario for growth in "set areas"
:: >>driverGrowth * AreaIntensity-> setArea

##
## Land use allocation LP
##
summvar addArea{LUC} takeArea{LUC}
var deltaPlus{toLUC, fromLUC} deltaMinus{toLUC, fromLUC}

# Set up "identity" arrays
var toIdent{LUC, toLUC} fromIdent{LUC, fromLUC}
ditto 'Cropland':
   toIdent{LUC='', toLUC=''}     = 1
   fromIdent{LUC='', fromLUC=''} = 1
ditto 'Pasture & Range':
   toIdent{LUC='', toLUC=''}     = 1
   fromIdent{LUC='', fromLUC=''} = 1
ditto 'Forest':
   toIdent{LUC='', toLUC=''}     = 1
   fromIdent{LUC='', fromLUC=''} = 1
ditto 'Developed':
   toIdent{LUC='', toLUC=''}     = 1
   fromIdent{LUC='', fromLUC=''} = 1
ditto 'Other':
   toIdent{LUC='', toLUC=''}     = 1
   fromIdent{LUC='', fromLUC=''} = 1

# Store the base year area data
summvar byArea{LUC}
:: byvalue(Area) -> byArea

##### The LP itself
LP: scen
solve for Area ATM addArea takeArea \
          deltaPlus deltaMinus
# LGP objective function
min deltaPlus + deltaMinus

# Sum across rows and columns in ATM to get additions
#   and removals of land by LUC
toIdent * ATM - takeArea  = 0{LUC}
fromIdent * ATM - addArea = 0{LUC}

# Remove and add land use areas to base year areas
# (Start with Area = byArea + addArea - takeArea, then rearrange)
Area - addArea + takeArea = byArea

# Make sure set areas are equal to their set values
setFilter * Area = setArea

# Ensure total area is preserved, as a sum over all LUCs
addArea - takeArea = 0

# Implement LGP: make as close to "targetATM" as possible, in % terms
ATM + targetATM * deltaPlus  >= targetATM
ATM - targetATM * deltaMinus <= targetATM

:LP


##
## Report results
##
report 0.001 * Area as "Area (1000 acres)"
ditto {LUC='Forest'}:
report 0.001 * addArea{LUC='Forest'} as "New Forest (1000 Acres)"
report 0.001 * (Area'' - takeArea'') as "Remaining Forest (1000 acres)"