* Water inflows are based on streamflows for long run period of record * Management institution is current (2001) law of the river * Written by Jim F. Booker, University of Colorado (1998-2001) * * Compact equations added on 10/16/00. * New ag demand calcs and figures from 11/14/00 added. * Non-linearities removed: Vu,Bosque,Evap 12/01/00. * Modified full RGR project release of 790,000 going into a drought; implementation 12/00. * Conejos flows increased to include Los Pinos and San Antonio 10/22/00. * Conveyance function from Lobatos to Embudo added 10/22/00; Future * non-constant additions will require changing CONVEY equation as linear * term is keyed to RG source flows. Embudo to Otowi also added as constant. * New interpretation of MRGCD consumption - 10/23/00. * * Basic hydrology: * MRG Water Assembly Figures, Oct. 1999 version * supplemented by SSPA MRG Water Supply Study (June, 2000), plus * Lower Rio Grande general Hydro Relationships, Aug. 1997 (case 96-888) * ********************************************************************** * Built for Regional Research Project * Institutional Adjustments for Coping with Severe and Sustained * Drought in the Rio Grande Basin; * Sponsored by US Geological Survey, * Water Resources Research Institutes, * and Agricultural Experiment Stations * of Colorado, New Mexico, and Texas ********************************************************************** *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * INITIAL CONTROL * OPTION RESLIM=100; $EOLCOM > *$ONSYMXREF include reference maps for now $OFFUPPER OPTION NLP=MINOS OPTION LIMROW=10,LIMCOL=10; *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * SET DEFINITIONS * SETS **************************************************************** i Flows--location of all flow nodes on Rio Grande CO to TX **************************************************************** / * river(i) * / Lobatos Embudo Chamita Otowi Cocgage Acacia Marcial EBgage ElPaso Quitman * / * inflow(i) * / RG Conejos Chama Jemez Puerco Salado * / * div(i) * / SLV MRG ABQ SOC EBID MX EP EPAG * / * return(i) * / SLV-r MRG-r ABQ-r SOC-r EBID-r MX-r EP-r EPAG-r * / * bosque(i) * / MRG-b SOC-b * / * convey(i) * / Lobatos-c Embudo-c Acacia-c Marcial-c EBgage-c * / * gw(i) * / SLV-g MRG-g ABQ-g SOC-g EBID-g MX-g EP-g EPAG-g * / * evap(i) * / ChamR-e CochR-e EBR-e * / * rel(i) * / ChamR-r CochR-r EBR-r / ; SETS * All subsets of SET i * Different nodes have different things happening at them * examples include places on the river, inflow locations, * diversion locations, coveyance fn locations, groundwater interaction spots,etc * [insert cells B6 to B### from worksheet ii below, without leading "*"s.] * [remove "+" from line 100 or so] river(i) / Lobatos Embudo Chamita Otowi Cocgage Acacia Marcial EBgage ElPaso Quitman / inflow(i) / RG Conejos Chama Jemez Puerco Salado / div(i) / SLV MRG ABQ SOC EBID MX EP EPAG / return(i) / SLV-r MRG-r ABQ-r SOC-r EBID-r MX-r EP-r EPAG-r / bosque(i) / MRG-b SOC-b / convey(i) / Lobatos-c Embudo-c Acacia-c Marcial-c EBgage-c / gw(i) / SLV-g MRG-g ABQ-g SOC-g EBID-g MX-g EP-g EPAG-g / evap(i) / ChamR-e CochR-e EBR-e / rel(i) / ChamR-r CochR-r EBR-r / ; PARAMETER startyear Defines order of starting year * for the run, given the structure of the input data file. ; startyear=2; SETS ***************************************************************** t time ***************************************************************** / 1*44/ >static set of time periods - can be expanded to many tsub(t) >dynamic subset of t ***************************************************************** year Calendar year ***************************************************************** /1941*1996/ >based on inflows data file ***************************************************************** u Stocks--location of all nodes on Rio Grande CO to TX ***************************************************************** / ChamR,CochR,EBR >reservoir stock nodes / ***************************************************************** policy range of institutions to be used * Used to capture results for different policies in a single run. Currently inactive. /law law of the river *bank1 administrative water bank with price and quantity limits bank2 administrative water bank with quantity limits only instream minimum instream flow requirements for minnow protection market market allocations within compact limits / type(policy) >dynamic subset of institutions. Used to choose type for a given run. ***************************************************************** * Institutional characteristics. Used primarily for reporting summary results. ag(div) /SLV,MRG,SOC,EBID,EPAG/ muni(div) /ABQ,EP/ ; ALIAS (t,tp); ALIAS (i,ip); > lets river flow node table be row or column ALIAS (river,riverp); *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * * READS MATRIX DEFINING MAPPING SET FOR BASIN FLOW OBJECTS * *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * Defines network locations and relationships TABLE ii(i,ip) *defined 6/13/00 Lobatos Embudo Chamita Otowi Cocgage Acacia Marcial EBgage ElPaso Quitman *river(i) River flows at gages */ Lobatos 1 Embudo 1 Chamita 1 Otowi 1 Cocgage 1 Acacia 1 Marcial 1 EBgage 1 ElPaso 1 Quitman */ *inflow(i) Native inflows to gages */ RG 1 Conejos 1 Chama 1 Jemez 1 Puerco 1 Salado 1 */ *div(i) Surface diversions and-or consumptive uses and-or groundwater basins */ SLV -1 MRG -1 ABQ -1 SOC -1 EBID -1 MX -1 EP -1 EPAG -1 */ *return(i) Direct return flows from diversions and-or use */ SLV-r 1 MRG-r 1 ABQ-r 1 SOC-r 1 EBID-r 1 MX-r 1 EP-r 1 EPAG-r 1 */ *bosque(i) Depletions by bosque vegetation */ MRG-b -1 SOC-b -1 */ *convey(i) Conveyance function additions aka fudge factors for unmeasured gains and losses */ Lobatos-c 1 Embudo-c 1 Acacia-c 1 Marcial-c 1 EBgage-c 1 */ *gw(i) Groundwater return flow from net seepage */ SLV-g 1 MRG-g 1 ABQ-g 1 SOC-g 1 EBID-g 1 MX-g 1 EP-g 1 EPAG-g 1 */ *evap(i) Reservoir evaporation */ ChamR-e CochR-e EBR-e */ *rel(i) Net release from reservoir storage */ ChamR-r 1 CochR-r 1 EBR-r 1 */ + SLV MRG ABQ SOC EBID MX EP EPAG * surface return flows come from ... SLV-r 1 MRG-r 1 ABQ-r 1 SOC-r 1 EBID-r 1 MX-r 1 EP-r 1 EPAG-r 1 * groundwater to river flows come from ... SLV-g 1 MRG-g 1 ABQ-g 1 SOC-g 1 EBID-g 1 MX-g 1 EP-g 1 EPAG-g 1 ; display i; option decimals=0;display ii; parameter iu(i,u) /ChamR-r.ChamR 1 ChamR-e.ChamR 1 CochR-r.CochR 1 CochR-e.CochR 1 EBR-r.EBR 1 EBR-e.EBR 1 /; display iu; *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * * DATA * *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * Constructed 12/01/00 faw * Inflows from river gages, long run averages for period of record TABLE inflows(year,*) RG Conejos Chama Jemez Puerco Salado Salado_only total Mogote Los_Pinos San_Antonio * Year 1941 659800 345760 439000 32238 1574919 239400 87500 18860 1942 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1943 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1944 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1945 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1946 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1947 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1948 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1949 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1950 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1951 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1952 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1953 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1954 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1955 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1956 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1957 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1958 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1959 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1960 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1961 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1962 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1963 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1964 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1965 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1966 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1967 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1968 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1969 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1970 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1971 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1972 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1973 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1974 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1975 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1976 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1977 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1978 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1979 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1980 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1981 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1982 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1983 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1984 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1985 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 1986 659800 345760 439000 45170 32238 40515 12436 239400 87500 18860 1987 659800 345760 439000 45170 32238 40515 12436 239400 87500 18860 1988 659800 345760 439000 45170 32238 40515 12436 239400 87500 18860 1989 659800 345760 439000 45170 32238 40515 12436 239400 87500 18860 1990 659800 345760 439000 45170 239400 87500 18860 1991 659800 345760 439000 45170 239400 87500 18860 1992 659800 345760 439000 45170 239400 87500 18860 1993 659800 345760 439000 45170 239400 87500 18860 1994 659800 345760 439000 45170 239400 87500 18860 1995 659800 345760 439000 45170 239400 87500 18860 1996 659800 345760 439000 45170 239400 87500 18860 ; * Aves 659800 345760 439000 45170 32238 40515 12436 1574919 239400 87500 18860 PARAMETER * This section converts USGS gaged flows (in af) by year into kaf flows by period t. * Note that Chama flows are those at Chamita, and thus already have SJC water included. source(inflow,t) Inflows SJC(t) San Juan - Chama water ; * Puts correct structure of table of inflows for the set t of years which will be * modeled in the given run. * All flows converted to 1,000 af (kaf) here. source(inflow,t)= SUM(year$ (ORD(t) EQ (ORD(year)-startyear+1)), inflows(year,inflow)/1000 ); source('Chama',t)= source('Chama',t)+46; * Adds evaporation loss to Chamita flows, because * model later takes out evaporation. * Choice of 46 kaf is calculated evap at average 1973-99 aggregate level of 548 kaf. * The following accounts for historical SJ-C flows embedded in Chamita flows (recently), * while adjusting for future conditions which include the presence of SJ-C flows. SJC(t)=0 + 55 $ (ORD(t) GE startyear+30); * For years 1972 on * For future conditions: * Take out embedded SJC flows from Chamita gaged flows. source('Chama',t)=source('Chama',t)-SJC(t); * Define future SJC flows. SJC(t)=96.2; * SJC(t)=75.855; * Add back future SJC flows to flows at Chamita. source('Chama',t)=source('Chama',t)+SJC(t); display source, SJC; *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * historical credit/debit, 1942-85 * Schedule of NM credits/debits in deliveries to Texas * from Papadopulos report, 2000 TABLE credits(year,*) nmcredit cocredit cobalance * year 1941 110 1942 -50 > spill year 1943 -60 1944 -80 1945 -10 1946 40 1947 -70 1948 -110 1949 10 1950 20 1951 -70 -24 1952 -120 -152 1953 -20 -18 1954 -20 -60 1955 20 -55 1956 -50 -63 1957 50 -156 1958 10 14 1959 -30 -21 1960 50 -59 1961 50 -54 1962 50 -86 1963 -10 -22 1964 -70 -77 1965 -30 -129 -940 1966 20 22.2 -917.8 1967 40 22.2 -895.6 1968 80 22.2 -873.3 1969 120 22.2 -851.1 1970 30 22.2 -828.9 1971 40 22.2 -806.7 1972 150 22.2 -784.4 1973 -80 22.2 -762.2 1974 50 22.2 -740.0 1975 60 22.2 -717.8 1976 -30 22.2 -695.6 1977 -20 22.2 -673.3 1978 -60 22.2 -651.1 1979 -100 22.2 -628.9 1980 -20 22.2 -606.7 1981 -50 22.2 -584.4 1982 30 22.2 -562.2 1983 50 28 -540 1984 130 109 -512 1985 472 95 > spill year 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 ; *- The following section of historical debits and credits is for calibration work only; *- the past record of credits and debits is not taken to preduct future conditions. PARAMETERS dnmpar(t) NM historical debits - positive number only cnmpar(t) NM historical credits - postive number only ccopar(t) CO historical credit-debits ; cnmpar(t)=SUM(year$ (ORD(t) EQ (ORD(year)-startyear+1)), credits(year,'nmcredit') ); dnmpar(t) $ (cnmpar(t) LE 0) = -cnmpar(t); cnmpar(t) $ (cnmpar(t) LE 0) = 0; ccopar(t)= SUM(year$ (ORD(t) EQ (ORD(year)-startyear+1)), credits(year,'cocredit') ); display dnmpar,cnmpar,ccopar; * ---------------------------------------------------------------------- * Elephant Butte end of year storage, in a-f. * Used for model calibration * ---------------------------------------------------------------------- TABLE ebstorage(year,*) storage * Year 1941 1937700 1942 1780500 1943 1250600 1944 1290600 1945 1130100 1946 584300 1947 435500 1948 492000 1949 621000 1950 335400 1951 36500 1952 376900 1953 110600 1954 97600 1955 155000 1956 32900 1957 776100 1958 988800 1959 586400 1960 439600 1961 360900 1962 390300 1963 112000 1964 87300 1965 517200 1966 344000 1967 267100 1968 333600 1969 528200 1970 324500 1971 177000 1972 301600 1973 794200 1974 402500 1975 617200 1976 315800 1977 181400 1978 182600 1979 943500 1980 1222000 1981 725300 1982 916700 1983 1322700 1984 1610500 1985 1955700 1986 2090100 1987 2043400 1988 2034700 1989 1675300 1990 1369600 1991 1713900 1992 1855100 1993 1907500 1994 2038700 1995 2040200 1996 1718100 ; PARAMETERS ebhistoric(t) Elephant Butte historical storage volume ebmodeled(t) Elephant Butte modeled storage - for calibration + error reporting delta(t) ; ebmodeled(t)=0; delta(t)=0; ebhistoric(t)=SUM(year$ (ORD(t) EQ (ORD(year)-startyear+1)),ebstorage(year,'storage')/1000 ); display startyear,ebhistoric; * ---------------------------------------------------------------------- PARAMETERS * The following section calculates lag coefficients used to calculate * limits to pumping, bosque depletions, and groundwater return flows, * based on net seepage from points of use which add to river surface flows. * ---------------------------------------------------------------------- * Pumping limits and bosque depletion limits, based on river flows. * Calculation of each uses identical linear lag calculations. * Provides upper bound after several low flow years. n(i) sumtime(div) sumgamma(i) pbase(div) flow independent pump limit pmax(div) full pumping limit pfull(div) flow needed to keep aquifer replenished over time bbase(bosque) flow independent bosque depletion bmax(bosque) maximum bosque depletion bfull(bosque) flow needed to keep bosque use replenished over time gwbase(gw) groundwater base flows prop(gw) proportion of gw flow ending up in surface water r0(gw) integrated unit groundwater withdrawal return flow over time alpha(gw) x0(i) divisor for depletion limits gamma(i,t) lagged impact coefficients ; TABLE pumppar(div,*) Pumping limits * base - base level of pumping always allowed (from deep, non-tributary sources). * max - maximum total pumping due to either physical or institutional factors. * xmax - average river flow needed to sustain maximimum pumping rate. * lagtime - linear lag time over which to calculate this average river flow. * The purpose here is to determine the degree to which pumping would be scaled back under * sustained low flow conditions. The parameter gamma is used, in conjunction with * modeled river flows, to determine the maximum level of pumping in any given time * period. Coefficients gamma(t) are used with modeled river flows to determine this * pumping limit (see equation Plimit1). base max xhalf lagtime SLV 0 200 100 3 ABQ 110 230 100 10 EBID 50 152 1 7 EP 100 150 1 10 ; n(div)=pumppar(div,'lagtime'); * maximum length of lagged impacts n(div)=MIN(n(div),CARD(t)); * limit lag to model time steps * Calculate gammas which define the lag structure for the flow dependent pumping limit. gamma(div,t)$(ORD(t) LE n(div))=1/n(div); * use simple average over n periods x0(div)=pumppar(div,'xhalf')/log(2); x0(div)$ (pumppar(div,'xhalf') EQ 0) = 1 ; * prevents division by zero for non-pumping uses pbase(div)=pumppar(div,'base'); * deep aquifer pumping unaffected by flows. pmax(div)=pumppar(div,'max'); * maximum possible pumping. ****************************************** TABLE bosqpar(bosque,*) Bosque depletion parameters * base -Base level of bosque depletion, even with 'zero' flow. * max - Maximimum level of depletion * xhalf - Average river flow at which depletion/limit is halved * lagtime - linear lag time over which to calculate average river flow. base max xhalf lagtime MRG-b 105 180 1000 3 SOC-b 40 75 900 1 ; * Calculate gammas which define the lag structure for the flow dependent pumping limit. n(bosque)=bosqpar(bosque,'lagtime'); n(bosque)=MIN(n(bosque),CARD(t)); * limit lag to model time steps gamma(bosque,t)$(ORD(t) LE n(bosque))=1/n(bosque); * use simple average over n periods x0(bosque)=bosqpar(bosque,'xhalf')/log(2); bbase(bosque)=bosqpar(bosque,'base'); * base bosque useage not affected by flows. bmax(bosque)=bosqpar(bosque,'max'); * maximum possible pumping. ***************************************************************************** * --------------------------------------------------------------------------- TABLE gwpar(gw,*) Groundwater return flow * base - base level of groundwater return flow * prop - total proportion of net seepage which ultimately is return flow * (Includes discharge to drainage ditches) * lagtime - linear lag time over which seepage returns to river. base prop lagtime SLV-g 0 0 3 MRG-g 60 1 3 ABQ-g 0 .60 38 SOC-g 0 1 3 EBID-g 0 .975 4 EP-g 0 0 15 ; * Calculate gammas which define the lag structure return flows from net seepage. n(gw)=gwpar(gw,'lagtime'); n(gw)=MIN(n(gw),CARD(t)); >limit lag to models time steps option decimals=2; display n; gamma(gw,'1')=n(gw); LOOP(gw, LOOP(t $ (ORD(t) LE n(gw)), >start by counting down from N gamma(gw,t+1)=gamma(gw,t)-1 ); ); sumgamma(gw)=SUM(t,gamma(gw,t)); >find summ gamma(gw,t)$sumgamma(gw)=gamma(gw,t)/sumgamma(gw); >normalize gamma(gw,t)$ (n(gw) EQ 0)=0; >setss non-specified basins to zero ******** These parameters will define actual groundwater return flow levels gwbase(gw)=gwpar(gw,'base'); > deep aquifer return flows unaffected by seepage. prop(gw)=gwpar(gw,'prop'); > maximum possible pumping. gamma(gw,t)=prop(gw)*gamma(gw,t); >setss gammas to reflect proportion returned to surface water. ***********CHECKS option decimals=3; display gamma; * By definition, gammas should sum proportion returned to surface water. * This is just a check to make sure they do. sumgamma(i)=SUM(t,gamma(i,t)); display sumgamma; * ---------------------------------------------------------------------------- * PARAMETERS * * The following section defines parameters for reservoir evaporation, and conveyance functions. * * * Evaporation is given by EVAP = b1*X + b2*X*X. * Surface area could be calculated using pan evaporation data; * methodology shaky for Chama because of hybrid nature of * reservoir. * b0(i,t) constant term b1(i,t) linear term b2(i,t) quadratic term b3(i,t) additional coefficient - use as needed evapest(t) EB evaporation estimate ; * Reservoir evaporation TABLE evappar(evap,*) b0 b1 b2 pan * (kaf) (unitless) (/kaf) (inches) ChamR-e 0.0995 -2.676E-05 56.0 CochR-e 0.1365 -9.496E-05 92.98 EBR-e 0.2378 -3.562E-05 111.17 ; b0(evap,t)=evappar(evap,'b0'); b1(evap,t)=evappar(evap,'b1'); b2(evap,t)=evappar(evap,'b2'); * Adjustments Elephant butte - seems much too high. Max is over 300 kaf/yr. Reduce by 2/3: b1('EBR-e',t)=0.775*b1('EBR-e',t); b2('EBR-e',t)=0.775*b2('EBR-e',t); * Net conveyance function parameters: current implementation. * form is Y(river+1) = b0 + b1*X('RG') where Y(river) is the added downstream flow, * net of all otherwise measured or accounted for additions and depletions. * these adjust for otherwise unexplained gains or losses * These can also be used for model calibration. TABLE conveypar(convey,*) b0 b1 Lobatos-c 95.994 0.2698 * b1 for Lobatos is multiplied by DelNorte flows. Embudo-c 35.12 Marcial-c EBgage-c -89.97 ; b0(convey,t)= conveypar(convey,'b0'); >modified in calibration methods forr ebgage-c. b1(convey,t)= conveypar(convey,'b1'); *----------------------------------------------------- PARAMETER * return flows and seepage * Implemented for each user. returnprop(i,t) seepprop(i,t) rate(i) consumptive use fraction ; TABLE divchar(div,*) retp seepp SLV 0 0.1 >arbitrary MRG 0.5570 0.1899 SOC 0.5570 0.1899 EBID 0.2047 0.3386 EPAG 0.2047 0.3386 MX 0 0 ABQ 0.4118 0.0588 EP 0.4118 0.0588 ; returnprop(return,t)=SUM(div,ii(return,div)*divchar(div,'retp') ); seepprop(div,t)=divchar(div,'seepp'); rate(div)=1-SUM(return,ii(return,div)*returnprop(return,'1') )-seepprop(div,'1'); display returnprop,seepprop,rate; * Reservoir data PARAMETER zmax(u) maximum storage zterminal(u) final storage - terminal condition beta0 constant - governs reduction from 790 release rule beta1 EBR coeff beta2 Salado coeff ; zmax('ChamR')=1000; zmax('CochR')=200; zmax('EBR')=1956; zterminal(u)=0.5*zmax(u); beta0 = -118.00; beta1 = 0.14; beta2 = -1.55; PARAMETERS * initial conditions z0(u); >initial levels at each stock node z0('ChamR')=548; z0('CochR')=53; *z0('EBR')=1938; z0('EBR')=SUM(year$ (ORD(year) EQ startyear-1),ebstorage(year,'storage')/1000); display z0; *INSTITUTIONAL AND ECONOMIC PARAMETERS PARAMETERS lawdiv(div,t) used to capture diversion levels under law of the river prior(div,t) priorities for diversions abqsjc(t) estimate of effective SJC diversion right. div0(div,t) maximum diversion per acre or household divmax(div,t) maximum total diversion use0(div,t) max consumptive use per acre or household usemax(div,t) max total consumptive use kacres(div,t) calculated number of acres in thousands khhlds(div,t) number of households in thousands khhldqvlt(muni,t) equivalent number of hhlds for all urban uses ; * The following are institutional limits. prior('SLV',t)=1; prior('MRG',t)=100; prior('SOC',t)=100; prior('ABQ',t)=100; prior('EBID',t)=2; prior('EPAG',t)=2; prior('EP',t)=2; prior('MX',t)=1.5; abqsjc(t)=0; >initialize TABLE divlimits(div,*) Limits to diversions * Consumptive use max Surface diversion max usemax divmax SLV 9999 9999 > SLV are arbitrary - controlled by Compact MRG 130 9999 *ABQ 14 0 >historic; grows with linear trend ABQ 90 SOC 65 9999 EBID 279 495 > use derived by running model with pumping capped at 115. MX 160 9999 EP 72.3 63 EPAG 148 324; usemax(div,t)=divlimits(div,'usemax'); divmax(div,t)=divlimits(div,'divmax'); lawdiv(div,t)=usemax(div,t); * Municipal historic diversions * usemax('ABQ',t)=rate('ABQ')*( divlimits('ABQ','usemax')+((133-14)/35)*(ORD(t)-1) ); * Municipal future diversions account for population growth in Albuquerque and El Paso usemax('ABQ',t)=rate('ABQ')*((divlimits('ABQ','usemax')/rate('ABQ'))+2.40*(ORD(t)-1) ); usemax('EP',t) =rate('EP') *((divlimits('EP','usemax') /rate('EP')) +3.76*(ORD(t)-1) ); * divmax('ABQ',t) $ (ORD(t) GE 10) = 97; divmax('EP',t)= divlimits('EP','divmax')+3.76*(ORD(t)-1); divmax('EPag',t)= divlimits('EPag','divmax')-3.76*(ORD(t)-1) *divmax(muni,t)=divmax(muni,'1')+ * (usemax(muni,t)-usemax(muni,'1'))/rate(muni); display usemax,divmax; TABLE betas(ag,*) Per acre benefit function parameters for ag diversions * Max diversion is that used in LP modeling - corresponds to that delivered at farm, not * diversion from river. * Water quantity here is in units of deliveries, not consumptive use * Constant linear squared linear max LP div kacres * tot wat tot wat gw only beta0 beta1 beta2 beta3 div0 kac * ($/ac) ($/ac-af($/ac-af($/ac-af2 (af/ac) (kac) SLV 0 265.50 -28.94 -11.01 2.92 600.000 >CSU AG ECON Mark Sperov thesis generalized w regs MRG -30.993 67.00 -5.92 -999 4.00 48.781 >Table 3.7, Papadolous report SOC -30.993 67.00 -5.92 -999 4.00 15.130 >with above, in ratio 50:160 from usemax EBID 137.54 94.20 -2.50 -76.42 6.00 78.899 >USBR Sum Stats 1991 + TAMU models John Ellis EPAG 0 193.00 -21.50 -999 4.00 49.637 >same ; * Calculate coefficients for estimating total value function - takes the per acre * ag benefit functions and translates to benefit over all acres and water usage. * New method added 11/14/00. div0(ag,t)=betas(ag,'div0'); kacres(ag,t)=betas(ag,'kac'); use0(ag,t)=usemax(ag,t)/kacres(ag,t); display use0; *Coefficients b give total benefit in $1,000 as a function of total consumptive use * at diversion div, and the groundwater withdrawal at div. * Note: marginals are now in $1,000/kaf, i.e. in $/af for convenience. b0(ag,t)=kacres(ag,t)*betas(ag,'beta0'); b1(ag,t)=(div0(ag,t)/use0(ag,t))*betas(ag,'beta1'); b2(ag,t)=((div0(ag,t)/use0(ag,t))**2)*betas(ag,'beta2')/kacres(ag,t); b3(ag,t)=(div0(ag,t)/use0(ag,t))*betas(ag,'beta3'); display div0,use0,kacres,rate; *-------------------------------------- TABLE mibetas(muni,*) Per household benefits for M&I diversions * Data based on Frank Ward interpretation of published M&I work by Tom McGuckin. * Constant linear squared linear in gw diversion houselholds beta0 beta1 beta2 beta3 div0 khhlds * ($/hhld) ($/af) ($ hhld/af^2)($/af) (af/hhld) (khhld) ABQ 0 10843 -9627 -10 0.54135 107.000 EP 0 9507 -9392 -100 0.48979 120.553 ; * water quantity is in units of deliveries, not consumptive use * Calculate coefficients for estimating total value function - takes the per household * muni benefit functions and translates to benefit over all households. div0(muni,t)=mibetas(muni,'div0'); use0(muni,t)=rate(muni)*div0(muni,t); * Alternative household numbers khhldqvlt(muni,t)= usemax(muni,t)/use0(muni,t); khhlds(muni,t) = mibetas(muni,'khhlds'); * Growth in number of households to be added here. * Coefficients b give total benefit in $1,000 as a function of total consumptive use * at diversion div, and the groundwater withdrawal at div. * Note: marginals are now in $1,000/kaf, i.e. in $/af for convenience. b0(muni,t)=mibetas(muni,'beta0')*khhldqvlt(muni,t); b1(muni,t)=mibetas(muni,'beta1')/rate(muni); b2(muni,t)=mibetas(muni,'beta2')/(khhldqvlt(muni,t)*rate(muni)**2); b3(muni,t)=mibetas(muni,'beta3')/rate(muni); display kacres,khhlds,khhldqvlt; option decimals=3; display rate,betas,mibetas, b0,b1,b2,b3; option decimals=0; *------------------------ * Recreation benefits * * Data from Richard Cole and Frank Wards published research based on Riofish, 1980-1990 * Chama composite calculation by J. Booker * * Form is Benefit = a + bZ + cZ**2 PARAMETERS a(u) Constant coefficient b(u) linear coefficient c(u) quadratic coefficient ; TABLE recpars(u,*) a b c ChamR 4246.98 7.359 -0.002092785 CochR 256.62 4.104 -0.002875613 EBR 379.82 2.210 -0.000503852 ; a(u)=recpars(u,'a'); b(u)=recpars(u,'b'); c(u)=recpars(u,'c'); *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * * VARIABLE DEFINITIONS * *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . FREE VARIABLES X(i,t) water flows - inflows and outflows SEEP(div,t) gross seepage from diversions NS(div,t) net seepage from diversions and pumping DCO(t) Colorado and New Mexico debits - or credits - from under or over delivery *DNM(t) ADCO(t) Accrued debit - credit ADNM(t) TRANSFER transfer of diversions - compared to law of the river DSUB total diversions over a subset of t USUB total use over a subset of t Vi(div,t) total economic benefit function for consumptive use Vu(u,t) total economic benefit function for reservoir recreation ViSUB total of all consumptive use benefits over a subset of t VuSUB total of all recreation benefits over subset of t VSUB total of all benefits over a subset of t BANK2SUB ; POSITIVE VARIABLES Z(u,t) water stocks - reservoirs USE(div,t) consumptive use by diversions PUMP(div,t) pumping at diversions DNM(t) NM compact debits CNM(t) NM compact credits NMEXCESS(t) NM compact deliveries in excess of compact and credits SPILL(t) EBID unable to use all available surface water EBSHORT(t) Project release shortfall EBSPILL(t) EB spill BSHORT(t) bosque violation - essential system flexibility BEXCESS(t) bosque excess depletions *DEFICIT(river,t) violation of minimum streamflow constraint *TRANSVIOL(div,t) excess transfers - penalty set to 1000 dollars per af ; *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * * EQUATION DEFINITIONS and STATEMENTS * *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . EQUATIONS *=================================================================== *... Hydrologic balance: flows balance(river,t) mass balance equation for gaged flows river0(river,t) non-negative river flows Sources(inflow,t) set source nodes *=================================================================== *... Hydrologic balance: reservoirs reservoirs(u,t) stock mass balance accounting evapcalc(evap,t) reservoir evaporation calculation *net release equations to be added here *=================================================================== *... Hydrologic flow adjustments gwflow(gw,t) find impacts of net seepage to aquifers from diversions minus pumping bosqueuse(bosque,t) bosque depletions - fixed conveyance(convey,t) calculate net conveyance function *=================================================================== *... Disposition of offstream diversions between uses and various losses returns(return,t) relate return flows to diversion levels seepage(div,t) find gross seepage netseep(div,t) net seepage - difference between pumping and gross seepage cnsmption(div,t) consumptive use by each diversion *=================================================================== *... Set maximums and mins for off stream uses divmin(div,t) set minimum level on diversions - usually zero divlimit(div,t) set maximum diversions uselimit(div,t) set maximum consumptive use pumplimit1(div,t) maximum pumping based on flow - aquifer conditions pumplimit2(div,t) minimum pumping level for El Paso *=================================================================== *... Reservoir limits resmax(u,t) reservoir capacity resterm(u,t) reservoir terminal condition EBflow(t) Annual flow from project storage *=================================================================== *... Compact constraints COxNM(t) CO delivery requirement at Lobatos NMxNMTX(t) NM delivery to project water COcredit(t) annual credit limit COdebit(t) annual debit limit COaccrued1(t) accrued debit-credit definition COaccrued2(t) accrued debit limit NMcredit(t) annual credit limit NMdebit(t) annual debit limit NMaccrued1(t) accrued debit-credit definition NMaccrued2(t) accrued debit limit *=================================================================== *... More institutional stuff NMdistrib(t) distribution between MRG and SOC ABQdistrib(t) distribution of between Abq and MRG TXdistrib(t) distribution between EP and EPag ABQdiv(t) diversion limit for ABQ - uses SJC and net returns NMxTX(t) distribution of project water between TX and EBID MXmin(t) delivery to Mexico *TrnsfrDEF(div,t) Quantity of diversion transfered from law of the river base *Trnsfrlimp(div,t) Limit on size of transfers - currently set to 20% *Trnsfrlimn(div,t) this takes care of -ve transfers *=================================================================== * Endogenous economic variables Benefiti(div,t) total benefit function for consumptive uses Benefitu(u,t) total benefit function for reservoir recreation *=================================================================== *... Objective function definitions OBJdiv total diversions over subset tsub OBJuse total use over a subset of time t OBJiben total benefit in consumptive uses only OBJben total benefit over a subset of time t OBJbank2 total benefit with penalty for excess water transfer ; *=================================================================== * MODEL EQUATIONS *=================================================================== *... Hydrologic balance: flows Balance(river,t)$tsub(t).. X(river,t)=E=SUM(ip, X(ip,t)*ii(ip,river)); *- Mass balance of all inflows and outflows at each model node [designated by the set river]. *- Possible flows present at a model node include inflows, diversions, surface return flows, *- groundwater return flows and losses, bosque (riparian vegetation) depletions, *- reservoir evaporation losses, changes to reservoir storage levels net of evaporation, *- and other uncategorized conveyance gains or losses. River0(river,t)$tsub(t).. X(river,t)=G=0; *- Non-negativity constraint on river flows. Sources(inflow,t)$tsub(t).. X(inflow,t)=E=source(inflow,t); * Definition of inflows. *=================================================================== *... Hydrologic balance: stocks reservoirs(u,t)$tsub(t).. Z(u,t)=E=z0(u)$(ORD(t) EQ 1) +Z(u,t-1)-SUM(i,iu(i,u)*X(i,t)); *- Intertemporal mass balance of reservoir stocks, for each reservoir, updates *- reservoir contents by adjusting for reservoir evaporation losses (in the *- previous period), and gains in storage from retained river flows, net of evaporation. *evapcalc(evap,t)$tsub(t).. X(evap,t)=E=b0(evap,t)+b1(evap,t)*SUM(u,iu(evap,u)*Z(u,t-1)) * +b2(evap,t)*SUM(u,iu(evap,u)*Z(u,t-1)*Z(u,t-1)); evapcalc(evap,t)$tsub(t).. X(evap,t)=E= 46 $ (ORD(evap) EQ 1) + 7 $ (ORD(evap) EQ 2) + evapest(t) $ (ORD(evap) EQ 3); *- Reservoir evaporation is calculated as a quadratic function of reservoir volume in *- previous time periods. The function is represents the relationship between volume and *- surface area, and assumes a constant evaporation rate per unit area. *=================================================================== *... Hydrologic flow adjustments *... Groundwater return flows from net seepage Gwflow(gw,t)$tsub(t).. X(gw,t)=E=gwbase(gw)+ SUM(tp $ (ORD(tp) LE ORD(t)), SUM(div,ii(gw,div)*gamma(gw,tp)*NS(div,t-(ORD(tp)-1)) ) ) +SUM(tp $ (ORD(tp) GT ORD(t)), SUM(div,ii(gw,div)*gamma(gw,tp)*NS(div,'1') ) ); *- Net gains or losses from groundwater return flows are a function of the lagged seepage *- from, or depletion to shallow tributary aquifers. Net seepage, the difference between *- percolation associated with water use, and pumping *- depletions in the same aquifer, is used together with the lag structure to calculate the *- net effect on river flows in any given time period. * *- The lag is a simple linearly declining function of net seepage. The lag time may vary *- from just the current year (no lag) to the full number of model time steps (years). *- The proportion of net seepage impacting river flows over the full lag ranges from *- zero to one. For lags longer than the number of time steps to the first modeled period (e.g. *- a five year lag in model year 3), the net seepage in period one is used as a proxy *- for the missing periods. bosqueuse(bosque,t)$tsub(t).. X(bosque,t)=E=bbase(bosque)+ bmax(bosque)*(1-EXP( -(SUM(tp $ (ORD(tp) LE ORD(t)), SUM(river,-ii(bosque,river)*gamma(bosque,tp)*X(river,t-(ORD(tp)-1)) ) ) +SUM(tp $ (ORD(tp) GT ORD(t)), SUM(river,-ii(bosque,river)*gamma(bosque,tp)*X(river,'1') ) ) )/x0(bosque) ) ) +(-BSHORT(t)+BEXCESS(t)) $ (ORD(bosque) EQ 1); *- Bosque use, or riparian depletions, are represented as an increasing function of lagged river flows. *- The purpose of the the functional form is to capture bosque use of shallow, river flow dependent, *- groundwater. *- Specifically, X=base + bmax(1-exp(-alpha*sum(lagged flows)), where base defines a river independent *- constant bosque depletion, and the second term utilizes a simple average of recent (e.g. three years) *- river flows. The maximum bosque use in any one year is then base+bmax. *[Simple linear formulation is below] *bosqueuse(bosque,t)$tsub(t).. X(bosque,t)=E=bbase(bosque)+bmax(bosque)/2 * +(-BSHORT(t)+BEXCESS(t)) $ (ORD(bosque) EQ 1); *... Net Conveyance Function. Defined explicitly for now for RG. conveyance(convey,t)$tsub(t).. X(convey,t) =E= b0(convey,t)+ b1(convey,t)*source('RG',t) ; *- A conveyance function is defined to account for otherwise unmeasured *- gains and depletions. It is a linear function of river flow. *=================================================================== *... Division of sources between use and various losses cnsmption(div,t)$tsub(t).. USE(div,t)=E=X(div,t)+PUMP(div,t)-SEEP(div,t) -SUM(return,ii(return,div)*X(return,t)); *- Consumptive use is defined as the difference between surface diversions plus pumped *- groundwater, and surface return flows plus deep percolation. The consumptive *- use defines the quantity of flows which are lost (through evapotranspiration, or simple *- evaporation) to any future use by the system. Returns(return,t)$tsub(t).. X(return,t)=E=SUM(div,returnprop(return,t)* (X(div,t)+PUMP(div,t))*ii(return,div) ); *- Surface return flows result from both surface diversions and groundwater pumping. Seepage(div,t)$tsub(t).. SEEP(div,t)=E=seepprop(div,t)*(X(div,t)+PUMP(div,t)); *... Net seepage calculation netseep(div,t)$tsub(t).. NS(div,t)=E=SEEP(div,t)-PUMP(div,t); *- Net seepage is the difference between percolation to shallow aquifers, and *- pumping from these aquifers. *=========================================== *... Set maximums and mins for off stream uses divmin(div,t)$tsub(t).. X(div,t) =G= 0; divlimit(div,t)$tsub(t).. X(div,t) =L= divmax(div,t); uselimit(div,t)$tsub(t).. USE(div,t) =L= usemax(div,t); *- Limits are defined from the total consumptive use, and surface diversions for each user. pumplimit1(div,t)$tsub(t).. PUMP(div,t)=L= pbase(div)+(pmax(div)-pbase(div))*(1-EXP( -(SUM(tp $ (ORD(tp) LE ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*X(river,t-(ORD(tp)-1)) ) ) +SUM(tp $ (ORD(tp) GT ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*X(river,'1') ) ) )/x0(div) ) ) ; *- Groundwater pumping limits reflect both available infrastructure, and the short-term *- possiblity of substantial drawdown or depletion of shallow aquifers during drought. *- The latter effect is captured through a pumping limit which is a decreasing function of *- lagged river flows. *- The purpose of the the functional form is to capture decreasing ability to pump from *- shallow, river flow dependent, groundwater. *- Specifically, Plimit=pbase + (pmax-pbase)*(1-exp(-alpha*sum(lagged flows)), where pbase *- defines the pumping capacity from flow independent, very large deep groundwawater reserves, *- and the second term defines the pumping capacity from the shallow aquifer. *- Lagged flows utilize a simple average of recent (e.g. three years) *- river flows. pumplimit2(div,t)$tsub(t).. PUMP('EP',t)=G=51.127+22.136; >existing level of pumping; *- El Paso is required to maintain an base level of groundwater pumping no lower *- than the absolute level of 1999 pumping. Increasing future water demands are satisfied *- largely from increased utilization of surface water. *=================================================================== *... Reservoir operating rules resmax(u,t)$tsub(t).. Z(u,t)=L=zmax(u); resterm(u,t)$(ORD(u) EQ 2).. Z(u,t)=G=zterminal(u); >Keep Cochiti above minimum *resterm(u,t)$(ORD(t) EQ CARD(t)).. * Z(u,t)$(ORD(t) EQ CARD(t))=G=zterminal(u); EBflow(t)$tsub(t).. X('EBgage',t)=E= 790 + beta0 + beta1 * Z('EBR',t) + beta2 * X('Salado',t)- EBSHORT(t) + EBSPILL(t); *- Penalty functions, EBSHORT and EBSPILL are used for departures of river flows *- below Elephant Butte Reservoir from a level of 790 thousand a-f per year. * This says that EBgage flow = 790 - 118 +0.14 EB storage -1.55 Salado flow * Reduction in flows below 790 based on reservoir level & Salado flows * This means that districts reduce Rio Grande project releases below * full releases as EB volume declines in drought periods, as a * measure for carryover storage * *=================================================================== *... Compact constraints. Modified 10/16/00 COxNM(t)$tsub(t).. X('Lobatos',t)=E=-10+0.03292*X('RG',t)+0.1149*X('Conejos',t) +0.00040325*X('RG',t)*X('RG',t) +0.00083979*X('Conejos',t)*X('Conejos',t)-DCO(t); *- Colorado's obligation to New Mexico under the Rio Grande Compact, as laid out in the *- Compact delivery schedules, is captured by quadratic functions defining the obligation *- given the Rio Grande *- and Conejos supply indices, respectively. Departures from the schedule result in debits *- or credits, DCO, charged to Colorado. For this report, Colorado debits and credits are assumed *- to be zero in all years. NMxNMTX(t)$tsub(t).. X('EBgage',t)-(X('EBR-r',t)+X('EBR-e',t))=E=0.56025*(X('Otowi',t)-SJC(t)) +0.00010839*(X('Otowi',t)-SJC(t))*(X('Otowi',t)-SJC(t)) -DNM(t)+CNM(t)+NMEXCESS(t); *- New Mexico's obligation to water users below Elephant Butte Reservoir is captured by a *- quadratic function defining required flows given the Otowi supply index. Departures from *- the schedule results in debits or credits, DNM and CNM, respectively, charged to New Mexico. *- Any flows in excess of those which accrue as credits under the Compact are captured by NMEXCESS. COcredit(t)$tsub(t).. -DCO(t)=L=0; COdebit(t)$tsub(t).. DCO(t)=E=0; * Colorado credit/debit set to zero for this study COaccrued1(t)$tsub(t).. ADCO(t)=E=0.85*ADCO(t-1)+DCO(t); COaccrued2(t)$tsub(t).. ADCO(t)=L=99999; *NMcredit(t)$tsub(t).. CNM(t)=E=cnmpar(t); *NMdebit(t)$tsub(t).. DNM(t)=E=dnmpar(t); * Above is for calibration work NMcredit(t)$tsub(t).. CNM(t)=L=999; NMdebit(t)$tsub(t).. DNM(t)=L=999; * NM credit/debit limits are altered below using *.lo and *.up bounds. NMaccrued1(t)$tsub(t).. ADNM(t)=E=0.85*ADNM(t-1)+DNM(t)-CNM(t); *- Accrued debits and credits are subject to system losses, including evaporation which would have *- occurred had the debit or credit not been incurred. We do not attempt to calculate *- such losses precisely, but estimate them at 15% annually. NMaccrued2(t)$tsub(t).. ADNM(t)=L=200; *- New Mexico, under the Compact, may accrue debits up to a total of 200 thousand a-f. *=================================================================== *... More institutional stuff NMdistrib(t)$tsub(t).. X('MRG',t)=E=(usemax('MRG',t)/usemax('SOC',t))*X('SOC',t); *- Use within MRGCD is assumed to be reduced proportionally when necessary to meet Compact *- obligations. While this neglects the reality of senior Native and acequia rights, it *- captures the reality that the dominant uses (by quantity) within the irrigation district *- are likely to be treated simlarly in times of water shortage. ABQdistrib(t)$tsub(t).. X('Abq',t)=L=abqsjc(t) +( (97-abqsjc(t))*(X('MRG',t)*rate('MRG')/usemax('MRG',t)) ) $ (ORD(t) GE 10); *- Upper bound on ABQ use. Requires that it not exceed the proportion of full supply which *- MRG enjoys in any given year. TXdistrib(t)$tsub(t).. X('EPag',t)=E=(divmax('EPag',t)/divmax('EP',t))*X('EP',t); *- Surface water used by El Paso is provided out of the allocation of the EPWID #1. Water users *- within the district are subject to the same allocation, and hence El Paso municipal use *- of surface water is cut back proportionally to remaining agricultural uses in times of less *- than full allocations. ABQdiv(t)$(tsub(t) AND ORD(t) GE 10).. X('Abq',t)=G=48.2+X('Abq-g',t)+X('Abq-r',t); *- Albuquerque surface diversions of its San Juan - Chama rights are limited in the model to *- those diversions leading to a net river depletion by Albuquerque equal to these rights. *- Return flows accruing to the river increase diversion rights, while the estimated depletions *- to the river resulting from pumping reduce the diversion right. * NMxTX(t)$tsub(t).. X('EBID',t)=E=0.567742*(1.3377994*X('EBgage',t)-X('MX',t)-89.97)-SPILL(t); * Alternative - may avoid underuse problems by TX. NMxTX(t)$tsub(t).. X('EBID',t)=E=(0.567742/0.432258)* (X('EPag',t)+X('EP',t)); *- Approximates the DII delivery rule by requiring proportional diversions from EBID and state of Texas *- water users. Estimated system losses of 89.97 thousand a-f annually from the DII rule *- are incorporated as a conveeyance loss below Elephant Butte Reservoir. MXmin(t)$tsub(t).. X('MX',t)=E=60; *- For model simplicity, and because of the potential issues raised with any future *- delivery reductions to Mexico, a constant 60 thousand acre-feet annual delivery is assumed. *TrnsfrDEF(div,t)$tsub(t).. TRANSFER(div,t)=E=lawdiv(div,t)-X(div,t); *Trnsfrlimp(div,t)$tsub(t).. TRANSFER(div,t)=L=0.20*lawdiv(div,t)+TRANSVIOL(div,t); *Trnsfrlimn(div,t)$tsub(t).. TRANSFER(div,t)=G=-0.20*lawdiv(div,t)-TRANSVIOL(div,t); *=================================================================== *... Total benefits Benefiti(div,t)$tsub(t).. Vi(div,t)=E=b0(div,t)+b1(div,t)*USE(div,t) +b2(div,t)*USE(div,t)*USE(div,t) +b3(div,t)*PUMP(div,t)*rate(div); *- Total benefits of water use are quadratic function of total consumptive use, minus *- the net added cost per unit of consumptive use derived from pumped groundwater *- rather than surface water. NOT ENDOGENOUS IN 9.30. Benefitu(u,t)$tsub(t).. Vu(u,t)=E=a(u)+b(u)*Z(u,t)+c(u)*Z(u,t)*Z(u,t); *- Recreation benefits are estimated as a quadratic function of reservoir volume. *- Implicit in the function is the dependence of benefits on reservoir surface area. *- NOT ENDOGENOUS *=================================================================== *... Objective function definitions OBJdiv.. DSUB=E= SUM(div,SUM(tsub,prior(div,tsub)*X(div,tsub)) ) +1.5*SUM(tsub,X('Quitman',tsub)) +1*SUM(div,SUM(tsub,PUMP(div,tsub)) ) +0.1*SUM(tsub,Z('EBR',tsub)) -5*SUM(tsub,DNM(tsub)+CNM(tsub) ) ; *- The objective function is constructed as a weighted sum of surface diversions, *- roundwater pumping, system outflow, and reservoir storage. Use of New Mexico debits *- and credits reduces the objective function value, but by less than the diversions *- which are enabled or limited, respectively. OBJuse.. USUB=E=SUM(div,SUM(tsub,prior(div,tsub)*USE(div,tsub)) ) +0.2*SUM(tsub,Z('EBR',tsub)) -0.1*SUM(tsub,DNM(tsub)+CNM(tsub) ) ; OBJiben.. ViSUB=E= SUM(div,SUM(tsub,Vi(div,tsub))) -1E3*SUM(tsub,DNM(tsub)+CNM(tsub) ) * -5E4*SUM(tsub,CNM(tsub) ) -5E3*SUM(bosque,SUM(tsub,BSHORT(tsub)+BEXCESS(tsub) )) -5E4*SUM(tsub,NMEXCESS(tsub) ) * -1E4*SUM(tsub,EBSHORT(tsub)+EBSPILL(tsub)) ; *- An objective function maximizing benefits includes benefits from all offstream consumptive *- uses, and from reservoir recreation. The objective also addresses the use *- of credits and debits, and flow levels below Elephant Butte Reservoir. This *- is accomplished through use of an artificial cost to debits and credits, and to deviations *- of flow levels below Elephant Butte from a standard level of 790 thousand a-f. OBJben.. VSUB=E= ViSUB + SUM(u,SUM(tsub,Vu(u,tsub))) ; OBJbank2.. BANK2SUB=E=SUM(div,SUM(tsub,Vi(div,tsub))) + SUM(u,SUM(tsub,Vu(u,tsub))); * -1E6*SUM(div,SUM(tsub,TRANSVIOL(div,tsub))) * -1E7*SUM(river,SUM(tsub,DEFICIT(river,tsub))); *================================================================== *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * MODELS *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MODEL USE1 / balance,river0,Sources, reservoirs,evapcalc, gwflow,bosqueuse,conveyance, cnsmption,returns,seepage,netseep, divmin,divlimit,uselimit, pumplimit1,pumplimit2, resmax,EBflow, * resterm, COxNM,NMxNMTX, COcredit,COdebit,COaccrued1,COaccrued2, NMcredit,NMdebit,NMaccrued1,NMaccrued2, NMdistrib,TXdistrib,ABQdistrib,ABQdiv,NMxTX,MXmin, * TrnsfrDEF,Trnsfrlimp,Trnsfrlimn, Benefiti,Benefitu, OBJdiv,OBJuse,OBJiben *,OBJben *,OBJbank2 /; MODEL BENI0/ balance,river0,Sources, reservoirs,evapcalc, gwflow,bosqueuse,conveyance, cnsmption,returns,seepage,netseep, divmin,divlimit,uselimit, pumplimit1,pumplimit2, resmax,EBflow, * resterm, COxNM,NMxNMTX, COcredit,COdebit,COaccrued1,COaccrued2, NMcredit,NMdebit,NMaccrued1,NMaccrued2, NMdistrib,TXdistrib,ABQdiv,NMxTX,MXmin, * TrnsfrDEF,Trnsfrlimp,Trnsfrlimn, Benefiti,Benefitu, OBJdiv,OBJuse,OBJiben,OBJben *,OBJbank2 /; ******************************** * Model development fixed values X.up(rel,t)=999; >Reservoirs release up to this amt in ea period EBSHORT.up(t)=700; EBSPILL.up(t)=0; EBSPILL.up(t)$ (ORD(t) EQ 1 OR ORD(t) EQ CARD(t))=899; X.up('EBgage',t)=790 + 600 $ (ORD(t) EQ 1); Z.fx('chamR',t)=z0('chamR'); Z.fx('cochR',t)=z0('cochR'); Z.lo('EBR',t)=33; NMEXCESS.up(t)=0; CNM.up(t)=300; DNM.up(t)=150; *- Compact credit and debit limits for New Mexico are 200 and 150 thousand a-f annually, *- respectively. These limits are enforced in all years (excepting year 32, during the *- recovery from drought.) USE.fx('Abq',t)=usemax('Abq',t); *- Albuquerque utilizes groundwater pumping, together with its use of San Juan - Chama water, *- and native water to maintain full water use in all years. X.up('Abq',t)=97; *USE.lo('EP',t)=80 $ (ORD(t) GE 10); *PUMP.fx('EP',t)=51.127+22.136; X.lo('EPag',t)=15; *X.up('EPag',t)=(usemax('EPag',t)/rate('EPag'))-3.76*(ORD(t)-1)+1; PUMP.lo('EBID',t)=115-10; SPILL.fx(t)=0; BSHORT.up(t)=0; BEXCESS.fx(t)=0; *DEFICIT.fx(river,t)=0; *TRANSVIOL.fx(div,t)=0; *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * SOLUTIONS *. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . * Define output parameters to capture each solution. PARAMETER flowout(i,t,policy) storout(u,t,policy) useout(div,t,policy) benout(div,t,policy) debitout(*,t,policy) tranout(div,t,policy) ordyear ordt flag(t) EBgagfix; * Purpose of the following is to solve the model sequentially, year by year. Note that the * model can also be solved dynamically - but this implies perfect future knowledge! * Starting values: USE.l(div,t)=usemax(div,t); PUMP.l(div,t)=pmax(div); X.l(div,t)=usemax(div,t)/rate(div)-pmax(div)/2; X.l('EBgage',t)=790; X.l('EBID',t)=divmax('EBID',t); BSHORT.l(t)=0;BEXCESS.l(t)=0; DCO.l(t)=0;DNM.l(t)=0; CNM.l(t)=0; evapest('1')=255; tsub(t)=NO; >initialize *tsub(t)=YES; DISPLAY TSUB; LOOP(year $ (ORD(year) GE startyear AND (ORD(year)-startyear+1) LE CARD(t)), tsub(t)$ (ORD(t) LE (ORD(year)-startyear+1)) = YES ; ordyear=ORD(year) DISPLAY TSUB,ordyear; * DISPLAY source; * Show solution in .lst file only for first and last period. IF (ordyear EQ startyear, OPTION SOLPRINT=ON; ELSEIF ((ORD(year)-startyear+1) EQ CARD(t)), OPTION SOLPRINT=ON; ELSE OPTION SOLPRINT=OFF; ); * Set ABQ diversion of SJC water before solving model. *abqsjc(t)$ (ORD(t) GE 10)=48.2/(1-returnprop('Abq-g',t))-X.l('Abq-g',t-1); abqsjc(t)$ (ORD(t) GE 10)=48.2+X.l('Abq-g',t-1)+X.l('Abq-r',t-1); display abqsjc; * Total benefits solution, including recreation values, no instream constraint, Compact compliant. SOLVE USE1 USING NLP MAXIMIZING DSUB; * SOLVE BENI0 USING NLP MAXIMIZING ViSUB; ********************************************************************* * Fix all variables before moving on to next time step. X.fx(i,tsub)=X.l(i,tsub);SEEP.fx(div,tsub)=SEEP.l(div,tsub);NS.fx(div,tsub)=NS.l(div,tsub); DCO.fx(tsub)=DCO.l(tsub);DNM.fx(tsub)=DNM.l(tsub); CNM.fx(tsub)=CNM.l(tsub); ADCO.fx(tsub)=ADCO.l(tsub);ADNM.fx(tsub)=ADNM.l(tsub); BSHORT.fx(tsub)=BSHORT.l(tsub); BEXCESS.fx(tsub)=BEXCESS.l(tsub); EBSHORT.fx(tsub)=EBSHORT.l(tsub); EBSPILL.fx(tsub)=EBSPILL.l(tsub); *NMEXCESS.fx(tsub)=NMEXCESS.l(tsub); *Vi.fx(div,tsub)=Vi.l(div,tsub);Vu.fx(u,tsub)=Vu.l(u,tsub); Z.fx(u,tsub)=Z.l(u,tsub);USE.fx(div,tsub)=USE.l(div,tsub);PUMP.fx(div,tsub)=PUMP.l(div,tsub); *DEFICIT.fx(river,tsub)=DEFICIT.l(river,tsub); *TRANSVIOL.fx(div,tsub)=TRANSVIOL.l(div,tsub); * Prepare for next solve evapest(t+1)= b0('EBR-e',t+1)+b1('EBR-e',t+1)*Z.l('EBR',t) +b2('EBR-e',t+1)*Z.l('EBR',t)*Z.l('EBR',t); ); > END OF TIME STEP LOOP. * ---------------- DISPLAYS FOLLOW ----------------------------------------------- * *Define attributes of uses: SETS state /CO,NM1,NM2,TX/ sector /ag,muni,rec,env/; * Attributes TABLE attloc(i,state) * State locations * CO NM1 NM2 TX *river(i) River flows at gages */ Lobatos 1 Embudo 1 Chamita 1 Otowi 1 Marcial 1 EBgage 1 Quitman 1 */ *inflow(i)Inflows to gages */ RG 1 Conejos 1 Chama 1 Jemez 1 Puerco 1 Salado 1 */ *div(i) Surface diversions */ SLV 1 Abq 1 MRG 1 SOC 1 EBID 1 MX EP 1 EPAG 1 */ *return(i)Direct return flows from diversions */ SLV-r 1 Abq-r 1 MRG-r 1 SOC-r 1 EBID-r 1 EP-r 1 EPAG-r 1 */ *bosque(i)Depletions by bosque vegetation */ MRG-b 1 SOC-b 1 */ *convey(i)Conveyance function additions aka fudge factors for unmeasured gains and losses */ Lobatos-c 1 Embudo-c 1 Acacia-c 1 Marcial-c 1 EBgage-c 1 */ *gw(i) Groundwater return flow from net seepage */ SLV-g 1 Abq-g 1 MRG-g 1 SOC-g 1 EBID-g 1 EP-g 1 EPAG-g 1 */ ; TABLE attsector(i,sector) Sector definitions * div ag muni rec env SLV 1 Abq 1 MRG 1 SOC 1 EBID 1 MX EP 1 EPAG 1 ; * * FLOW MASS BALANCE AT EACH RIVER NODE - use as diagnostic ************************** parameter priver1(t,river,*) FLOW MASS BALANCE CHECK; priver1(t,river,'river')=x.l(river,t); priver1(t,river,'inflow')=sum(inflow,ii(inflow,river)*x.l(inflow,t)); priver1(t,river,'convey')=sum(convey,ii(convey,river)*x.l(convey,t)); priver1(t,river,'bq')=sum(bosque,ii(bosque,river)*x.l(bosque,t)); priver1(t,river,'rel')=sum(rel,ii(rel,river)*x.l(rel,t))+ sum(evap,ii(evap,river)*x.l(evap,t)); priver1(t,river,'div-ret')=sum(div,ii(div,river)*x.l(div,t))+ sum(return,ii(return,river)*x.l(return,t)); priver1(t,river,'gw')=sum(gw,ii(gw,river)*x.l(gw,t)); priver1(t,river,'CHK+1')=1-priver1(t,river,'river')+ SUM(riverp,ii(riverp,river)*priver1(t,riverp,'river'))+ priver1(t,river,'inflow')+ priver1(t,river,'convey')+ priver1(t,river,'bq')+ priver1(t,river,'rel')+ priver1(t,river,'div-ret')+ priver1(t,river,'gw'); OPTION priver1:0:1:1; display priver1; *************************** parameter puse(t,div,*) MASS BALANCE AT USE NODES & CHECK; puse(t,div,'div')=x.l(div,t); puse(t,div,'pump')=pump.l(div,t); puse(t,div,'return')=sum(return,ii(return,div)*x.l(return,t)); puse(t,div,'seep')=seep.l(div,t); puse(t,div,'use')=USE.l(div,t); puse(t,div,'CHK-USE')=puse(t,div,'div')+ puse(t,div,'pump')- puse(t,div,'return')- puse(t,div,'seep'); OPTION puse:0:1:1; display puse; ************************** parameter pres(t,u,*) RESERVOIR MASS BALANCE; pres(t,u,'vol-ini')=z0(u)$(ORD(t) EQ 1)+z.l(u,t-1); pres(t,u,'vol-fin')=z.l(u,t); pres(t,u,'rel')=sum(rel,iu(rel,u)*x.l(rel,t)); pres(t,u,'evap')=sum(evap,iu(evap,u)*x.l(evap,t)); pres(t,u,'rec-ben')=a(u)+b(u)*Z.l(u,t)+c(u)*Z.l(u,t)*Z.l(u,t); OPTION pres:0:1:1; display pres; ************************** parameter pecon(t,div,*) ECONOMIC BENEFIT REPORTING V1(div,t) calculation of total benefit; V1(div,t) = b0(div,t) + b1(div,t) * USE.l(div,t) + b2(div,t) * USE.l(div,t)**2 + b3(div,t) * PUMP.l(div,t); pecon(t,div,'usemax')=usemax(div,t); pecon(t,div,'use')=USE.l(div,t); pecon(t,div,'pumpuse')=PUMP.l(div,t)*rate(div); pecon(t,div,'totben-mil')= v1(div,t)/1000; pecon(t,div,'avgben')=v1(div,t)/usemax(div,t); pecon(t,div,'marg-surf')=b1(div,t)+2*b2(div,t)*USE.l(div,t); pecon(t,div,'marg-grnd')=b1(div,t)+2*b2(div,t)*USE.l(div,t)+b3(div,t); OPTION pecon:0:1:1; display pecon; ************************** parameter annual1(t,*) ANNUAL SUMMARY 1 ; annual1(t,'Inflow')=SUM(inflow,source(inflow,t)); annual1(t,'EB-stor')=Z.l('EBR',t); annual1(t,'Marcial')=X.l('Marcial',t); annual1(t,'EBgage')=X.l('EBgage',t); annual1(t,'Use')=SUM(div,USE.l(div,t)); annual1(t,'MRGuse')=USE.l('MRG',t)+USE.l('SOC',t)+X.l('Abq',t)*rate('Abq'); annual1(t,'bq')=sum(bosque,x.l(bosque,t)); annual1(t,'Pumping')=SUM(div,PUMP.l(div,t)); OPTION annual1:0:1:1; display annual1; ************************** parameter annual2(t,*) ANNUAL SUMMARY 2: Index Flows ; annual2(t,'Conejs')=X.l('Conejos',t); annual2(t,'RioG')=X.l('RG',t); annual2(t,'Lab-sched')=-10+0.03292*X.l('RG',t)+0.1149*X.l('Conejos',t) +0.00040325*X.l('RG',t)*X.l('RG',t) +0.00083979*X.l('Conejos',t)*X.l('Conejos',t); annual2(t,'Lab-Actual')=X.l('Lobatos',t); annual2(t,'Otwi-ind')=X.l('Otowi',t)-SJC(t); annual2(t,'EB-effsup')=X.l('EBgage',t)-(X.l('EBR-r',t)+X.l('EBR-e',t)); annual2(t,'EB-sched')=0.56025*(X.l('Otowi',t)-SJC(t)) +0.00010839*(X.l('Otowi',t)-SJC(t))*(X.l('Otowi',t)-SJC(t)); OPTION annual2:0:1:1; display annual2; ************************** parameter annual3(t,*) ANNUAL SUMMARY 3: Sources users - NM ; annual3(t,'Inflow')=SUM(inflow,source(inflow,t)); annual3(t,'EBgage')=X.l('EBgage',t); annual3(t,'EB-stor')=Z.l('EBR',t); annual3(t,'ABQuse')=USE.l('ABQ',t); annual3(t,'ABQdiv')=X.l('ABQ',t); annual3(t,'ABQpump')=PUMP.l('ABQ',t); annual3(t,'MRGuse')=USE.l('MRG',t); annual3(t,'MRG-bq')=X.l('MRG-b',t); OPTION annual3:0:1:1; display annual3; ************************** parameter annual4(t,*) ANNUAL SUMMARY 4: Sources users - TX ; annual4(t,'Inflow')=SUM(inflow,source(inflow,t)); annual4(t,'EB-stor')=Z.l('EBR',t); annual4(t,'EBgage')=X.l('EBgage',t); annual4(t,'EBIDuse')=USE.l('EBID',t); annual4(t,'EBIDpump')=PUMP.l('EBID',t); annual4(t,'EPAGuse')=USE.l('EPAG',t); annual4(t,'EPdiv')=X.l('EP',t); annual4(t,'EPpump')=PUMP.l('EP',t); OPTION annual4:0:1:1; display annual4; ************************** parameter annual5(t,*) ANNUAL SUMMARY 5: Calibration ; annual5(t,'year')=1941+startyear+ORD(t)-2; annual5(t,'Otowi')=X.l('Otowi',t); annual5(t,'EBgage')=X.l('EBgage',t); annual5(t,'Bexcess')=-BSHORT.l(t)+BEXCESS.l(t); annual5(t,'NMexcess')=NMEXCESS.l(t); annual5(t,'EB-hist')=ebhistoric(t); annual5(t,'EBmodeled')=ebmodeled(t); annual5(t,'Quitman')=X.l('Quitman',t); OPTION annual5:0:1:1; display annual5; ******************* * MORE DISPLAYS * begin out_tabl.gms * ************************** parameter avg_water(*,div) Annual Average Water Supply By User avg_swater(*,div) Annual Average surface water by user avg_gwater(*,div) Annual Average ground water by user avg_vol(*, u) Annual Average Reservoir Volume ; avg_water('1942-85',div)=SUM(t,X.l(div,t)+PUMP.l(div,t))/44; avg_water('prop1',div)=SUM(t, X.l(div,t)+PUMP.l(div,t))/SUM(t,usemax(div,t)/rate(div)); avg_water('1950-59',div)=SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18),X.l(div,t)+PUMP.l(div,t))/10; avg_water('prop2',div)=SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18), X.l(div,t)+PUMP.l(div,t))/SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18),usemax(div,t)/rate(div)); avg_water('1960-69',div)=SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28),X.l(div,t)+PUMP.l(div,t))/10; avg_water('prop3',div)=SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28), X.l(div,t)+PUMP.l(div,t))/SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28),usemax(div,t)/rate(div)); avg_water('1970-79',div)=SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38),X.l(div,t)+PUMP.l(div,t))/10; avg_water('prop4',div)=SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38), X.l(div,t)+PUMP.l(div,t))/SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38),usemax(div,t)/rate(div)); avg_water('prop1','SLV')=9999; avg_water('prop2','SLV')=9999; avg_water('prop3','SLV')=9999; avg_water('prop4','SLV')=9999; OPTION avg_water:2:1:1; display avg_water; avg_swater('1942-85',div)=SUM(t,X.l(div,t) )/44; avg_swater('1950-59',div)=SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18),X.l(div,t) )/10; avg_swater('1960-69',div)=SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28),X.l(div,t) )/10; avg_swater('1970-79',div)=SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38),X.l(div,t) )/10; OPTION avg_swater:2:1:1; display avg_swater; avg_gwater('1942-85',div)=SUM(t, PUMP.l(div,t))/44; avg_gwater('1950-59',div)=SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18), PUMP.l(div,t))/10; avg_gwater('1960-69',div)=SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28), PUMP.l(div,t))/10; avg_gwater('1970-79',div)=SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38), PUMP.l(div,t))/10; OPTION avg_gwater:2:1:1; display avg_gwater; avg_vol('1942-85', u) = SUM(t, z.l(u,t)) / 44; avg_vol('1950-59', u) = SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18), z.l(u,t)) / 10; avg_vol('1960-69', u) = SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28), z.l(u,t)) / 10; avg_vol('1970-79', u) = SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38), z.l(u,t)) / 10; OPTION avg_vol:2:1:1; display avg_vol; SET div_econ(div) /SLV, MRG, ABQ, SOC, EBID, EP, EPAG/ ; parameter prop_ben(*,div_econ) Prop (0-1) econ div benefits of sw + gw use compared to MAX sw capac (all yrs) abs_ben(*,div_econ) Ave Absolute Econ Div Benefits from sw + gw use abs_vben(*) Ave Absolute Economic Volume Benefits from reservoirs only vimax(div_econ,t) Absolute econ div benefits of sw + gw use if MAX sw capacity used (1 yr) c_fact(div_econ) Correction Factor economic benefits /SLV 2.08 MRG 1 ABQ 1 SOC 1 EBID 0.75 EP 1 EPAG 1 / ; vimax(div_econ,t)= b0(div_econ,t) + b1(div_econ,t) * usemax(div_econ,t) +b2(div_econ,t) * usemax(div_econ,t) * usemax(div_econ,t); * leaves out pumping costs - hard to specify minimum pumping level generally prop_ben('1942-85',div_econ)=SUM(t, Vi.l(div_econ,t))/SUM(t,vimax(div_econ,t)); prop_ben('1950-59',div_econ)=SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18), Vi.l(div_econ,t))/SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18),vimax(div_econ,t)); prop_ben('1960-69',div_econ)=SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28), Vi.l(div_econ,t))/SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28),vimax(div_econ,t)); prop_ben('1970-79',div_econ)=SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38), Vi.l(div_econ,t))/SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38),vimax(div_econ,t)); OPTION prop_ben:3:1:1; display prop_ben; abs_ben('1942-85',div_econ)= c_fact(div_econ) * SUM(t, Vi.l(div_econ,t)) /44; abs_ben('1950-59',div_econ)= c_fact(div_econ) * SUM(t $(ORD(t) GE 9 AND ORD(t) LE 18), Vi.l(div_econ,t)) /10; abs_ben('1960-69',div_econ)= c_fact(div_econ) * SUM(t $(ORD(t) GE 19 AND ORD(t) LE 28), Vi.l(div_econ,t)) /10; abs_ben('1970-79',div_econ)= c_fact(div_econ) * SUM(t $(ORD(t) GE 29 AND ORD(t) LE 38), Vi.l(div_econ,t)) /10; OPTION abs_ben:0:1:1; display abs_ben; abs_vben('1942-85')=SUM((t,u), Vu.l(u,t)) /44; abs_vben('1950-59')=SUM((t,u) $(ORD(t) GE 9 AND ORD(t) LE 18), Vu.l(u,t)) /10; abs_vben('1960-69')=SUM((t,u) $(ORD(t) GE 19 AND ORD(t) LE 28), Vu.l(u,t)) /10; abs_vben('1970-79')=SUM((t,u) $(ORD(t) GE 29 AND ORD(t) LE 38), Vu.l(u,t)) /10; OPTION abs_vben:0; display abs_vben; * end of out_tabl.gms; parameters ravg(div,t) average river flows given lag n curt(div,t) curtailment ratio rtest(river,t) river flows; ravg(div,t)= SUM(tp $ (ORD(tp) LE ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*x.l(river,t-(ORD(tp)-1)) ) ) +SUM(tp $ (ORD(tp) GT ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*x.l(river,'1') ) ); curt(div,t)= (1-EXP( -(SUM(tp $ (ORD(tp) LE ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*x.l(river,t-(ORD(tp)-1)) ) ) +SUM(tp $ (ORD(tp) GT ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*x.l(river,'1') ) ) )/x0(div) ) ); rtest(river,t)=X.l(river,t); ************** parameter pump1test(div); LOOP(t, pump1test(div)= pbase(div)+(pmax(div)-pbase(div))*(1-EXP( -(SUM(tp $ (ORD(tp) LE ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*X.l(river,t-(ORD(tp)-1)) ) ) +SUM(tp $ (ORD(tp) GT ORD(t)), SUM(river,-ii(div,river)*gamma(div,tp)*X.l(river,'1') ) ) )/x0(div) ) ) ; ); display pump1test; display evapest; * THE END