indus.gms : Indus Agricultural Model

Description

This is an example model of the indus basin model family. These models
have been used by the Government of Pakistan to analyze investment
design and water policy issues.


Reference

  • Duloy, J H, and O'Mara, G T, Issues of Efficiency and Interdependence in Water Resource Investments: Lessons from the Indus Basin of Pakistan. Tech. rep., The World Bank, 1984.

Large Model of Type : LP


Category : GAMS Model library


Main file : indus.gms

$title Indus Agricultural Model (INDUS,SEQ=90)

$onText
This is an example model of the indus basin model family. These models
have been used by the Government of Pakistan to analyze investment
design and water policy issues.


Duloy, J H, and O'Mara, G T, Issues of Efficiency and Interdependence
in Water Resource Investments: Lessons from the Indus Basin of Pakistan.
Tech. rep., The World Bank, 1984.

Keywords: linear programming, investment planning, water policy
$offText

$sTitle Set Definition
Set
   c       'crop types'
           / basmati  'basmati rice'
             cotton
             berseem  'fodder crop'
             gram
             irri     'irri rice'
             sorghum
             maiz
             mustard
             sc-mill  'sugar cane for mill'
             sc-gur   'sugar cane processed on the farm'
             wheat                                       /
   cf(c)   'fodder crops' / berseem, sorghum /
   cfr(c)  'crops with yields given by optimal fertilizer response'
           / basmati, cotton, irri, sc-mill, sc-gur, wheat /
   cw(c)   / wheat  , cotton /
   cri(c)  'rice crops' / basmati, irri /
   cc(c)   'non-consumption crops'
           / gram, irri, maiz, mustard, sc-gur, wheat /
   g       'polygons'           / poly-17+19, poly-18 /
   l       'types of livestock'
           / bullock 'bullock pair', buffalo, cattle 'cattle cow' /
   q       'livestock commodities'
           / buff-milk 'buffaloes milk', catl-milk 'cattle cows milk', meat /
   sea     'seasons'                           / kharif, rabi         /
   ps      'price schedules'                   / 76-77,  80-81, 85-86 /
   psr(ps) 'price schedules used for this run' / 76-77                /
   i       'canals'                            / panjnad, abbasia     /
   sg(g)   'polygon with saline groundwater'   / poly-18              /
   fsg(g)  'polygon with fresh groundwater'    / poly-17+19           /

* cropping technology can be constructed from three codes i.e technology code, sequence codes and
* water stress level. these codes are given in set t, s and w.
   t       'technology' / bullock , semi-mech /
   s       'sequence'
           / standard  'standard sequence', el-plant  'early planting'
             la-plant  'late planting'    , qk-harv   'quick harvest'  /
   w       'water stress level'
           / standard  'no stress'   , heavy    'heavy stress'
             light     'light stress', january  'stress in january' /
   m       'months'
           / jan, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec /
   y       'year' / 1960*1975 /
   cnf(c)  'non fodder crops';

cnf(c)  = yes;
cnf(cf) =  no;

Parameter fc 'maunds to pounds conversion factor';
fc = 82.286;

$sTitle Crop Data
Table lando(c,t,m) 'land requirements for crops other than wheat and cotton by month'
                                         jan   feb  mar  apr  may  jun   jul   aug   sep   oct   nov   dec
  (basmati,irri).(bullock,semi-mech)                               1     1     1     1     1     1
   berseem.bullock                       1     1    1    1    1                       .5   1     1     1
   berseem.semi-mech                     1     1    1    1    1                       .25  1     1     1
   gram.bullock                          1     1    1                                 .5   1     1     1
   gram.semi-mech                        1     1    1                                 .25  1     1     1
   sorghum.(bullock,semi-mech)                       .5   .5   .5  1     1     1      .5    .5
   maiz.bullock                                                     .5   1     1     1     1     1      .5
   maiz.semi-mech                                                        1     1     1     1     1      .5
   mustard.(bullock,semi-mech)           1     1    1                                      1     1     1
  (sc-mill,sc-gur).(bullock,semi-mech)   1     1    1    1    1    1     1     1     1     1     1     1  ;

Table landcw(cw,t,s,w,m) 'land requirements for cotton and wheat by month'
                                               jan   feb  mar   apr   may   jun   jul   aug   sep   oct   nov   dec
   cotton.(bullock,semi-mech).standard.
                             (standard)                               1     1     1     1     1     1     1      .5
   cotton.bullock.el-plant.standard                              .5   1     1     1     1     1     1     1      .5
   cotton.semi-mech.la-plant.standard                                       1     1     1     1     1     1      .5
   wheat.bullock.standard.(standard,light,
                           heavy,january)      1     1    1     1                                    .5   1     1
   wheat.bullock.la-plant.(standard,light,
                           heavy,january)      1     1    1     1                                         1     1
   wheat.bullock.qk-harv.(standard,light,
                          heavy,january)       1     1    1      .5                                  .5   1     1
   wheat.semi-mech.standard.(standard,light,
                             heavy,january)    1     1    1     1                                         1     1
   wheat.semi-mech.la-plant.(standard,
                             light,heavy)      1     1    1     1                                          .5   1
   wheat.semi-mech.la-plant.january            1     1    1     1                                               1
   wheat.semi-mech.qk-harv.(standard,light,
                            heavy,january)     1     1    1      .5                                       1     1  ;

Table bso(c,t,m) 'bullock req. other than cotton and wheat in bullock pair hrs per month'
                        jan   feb   mar   apr   may   jun   jul   aug   sep   oct   nov   dec
   basmati.bullock                                   22.0  17.2                     2.0  15.6
   basmati.semi-mech                                                                     13.6
   berseem.bullock      2.0   4.0   4.0   2.0   2.0                    10.1  13.01  7.6   2.0
   gram.bullock                     7.0                                10.8   5.6
   irri.bullock                                      18.9  19.2                     1.5  18.4
   irri.semi-mech                                                                        16.4
   sorghum.bullock                 16.0   1.0   1.0  16.0   1.5  15.0   1.0   0.5
   maiz.bullock                                      10.8   4.5  14.2                     5.0
   maiz.semi-mech                                                 3.0
   mustard.bullock      1.0   1.0   1.0                                      10.8  10.1   1.0
   sc-mill.bullock     16.85 15.1  15.0  12.0   8.0   1.75                         10.0  12.5
   sc-mill.semi-mech          2.5               1.75  1.75
   sc-gur.bullock      17.85 11.6   8.2         1.75  1.75                         11.5  13.5
   sc-gur.semi-mech    11.0   6.5   2.5         1.75  1.75                         11.5  11.0;

Table bscw(cw,t,s,w,m) 'bullock requirements for cotton and wheat in bullock pair hours per month'
                                              jan   feb  mar    apr   may   jun   jul   aug   sep   oct   nov   dec
   cotton.bullock.standard.standard                                 17.09  15.2               1.0   1.0   1.0   1.0
   cotton.bullock.el-plant.standard                            16.0  4.0   12.3               1.0   1.0   1.0   1.0
   wheat.bullock.standard.standard                              8.4  8.4                           18.6  20.6
   wheat.bullock.standard.light                                 7.1  7.1                           18.6  20.6
   wheat.bullock.standard.(heavy,january)                       5.5  5.5                           18.6  20.6
   wheat.bullock.la-plant.standard                              7.4  7.4                                 39.2
   wheat.bullock.la-plant.light                                 6.3  6.3                                 39.2
   wheat.bullock.la-plant.(heavy,january)                       4.8  4.8                                 39.2
   wheat.bullock.qk-harv.standard                              16.8                                18.6  20.6
   wheat.bullock.qk-harv.light                                 14.3                                18.6  20.6
   wheat.bullock.qk-harv.(heavy,january)                       10.9                                18.6  20.6
   wheat.semi-mech.standard.standard                            5.9  5.9
   wheat.semi-mech.standard.light                               5.0  5.0
   wheat.semi-mech.standard.(heavy,january)                     3.8  3.8
   wheat.semi-mech.la-plant.standard                            5.2  5.2
   wheat.semi-mech.la-plant.light                               4.4  4.4
   wheat.semi-mech.la-plant.(heavy,january)                     3.4  3.4
   wheat.semi-mech.qk-harv.standard                            11.8
   wheat.semi-mech.qk-harv.light                               10.
   wheat.semi-mech.qk-harv.(heavy,january)                      7.7                                                ;

Table lbco(c,t,m) 'labor requirements for crops other than wheat and cotton in man hours per month'
                         jan   feb  mar   apr   may   jun   jul   aug   sep   oct   nov   dec
   basmati.bullock                                   29.1  88.8  65.9   5.6   5.6  47.9  17.6
   basmati.semi-mech                                  8.4  71.6  65.9   5.6   5.6  47.4  16.1
   berseem.bullock      32.3  41.5 41.8  29.4  18.2                    10.1  15.5  23.5  29.1
   berseem.semi-mech    31.1  39.8 40.3  28.7  17.6                     7.5  13.2  21.6  28.4
   gram.bullock                    21.7                                10.8   8.9   0.7   2.5
   gram.semi-mech                  20.0                                 1.6   4.4   0.7   2.5
   irri.bullock                                      22.9 122.3  35.9   8.4   6.4  43.9  20.4
   irri.semi-mech                                     5.3 105.5  29.9   8.4   6.4  43.4  17.9
   sorghum.bullock                 18.5   2.0   3.0  18.7   4.0  18.0   3.0   1.0
   sorghum.semi-mech                6.0   1.3   4.5   6.5   4.2   5.5   2.0   1.5
   maiz.bullock                                      10.8   4.5  16.7  44.2   2.5  26.8  27.0
   maiz.semi-mech                                           2.4   8.4  44.2   2.5  26.0  26.8
   mustard.bullock      21.7  20.5  6.0                                      10.8  13.6  13.9
   mustard.semi-mech    21.7  20.2  6.0                                       1.6   5.1  13.9
   sc-mill.bullock      90.0  85.0 95.0  72.1  30.0   5.05  3.0   3.0   3.0   1.5  85.0  95.0
   sc-mill.semi-mech    86.5  81.5 90.0  64.0  34.5   5.05  3.0   3.0   3.0   1.5  80.5  90.5
   sc-gur.bullock      159.   80.4 54.4  22.5   4.75  5.05  3.0   3.0   3.0   1.5 150.3 148.5
   sc-gur.semi-mech    151.   80.6 48.9  22.5   4.75  5.05  3.0   3.0   3.0   1.5 148.3 142.3;

Table lbcw(cw,t,s,w,m) 'labor requirements for wheat and cotton in man hours per month'
                                              jan   feb   mar   apr   may   jun   jul   aug   sep   oct   nov   dec
   cotton.bullock.standard.standard                                  35.6  18.4   2.5   7.6  13.2  41.3  57.4  22.1
   cotton.bullock.el-plant.standard                            26.8   6.7  20.5   2.5   7.6  13.2  41.3  57.4  22.1
   cotton.semi-mech.standard.standard                                 5.0   9.4   2.5   7.6  12.7  40.8  56.9  21.6
   cotton.semi-mech.la-plant.standard                                      14.4   2.5   7.6  12.7  40.8  56.9  21.6
   wheat.bullock.standard.standard            4.3   3.9   3.9  64.4  23.8                          18.6  23.9   3.9
   wheat.bullock.standard.light               4.3   3.9   3.9  54.7  20.2                          18.6  23.9   3.9
   wheat.bullock.standard.(heavy,january)     4.3   3.9   3.9  41.9  15.5                          18.6  23.9   3.9
   wheat.bullock.la-plant.standard            4.3   3.9   3.9  51.4  18.9                                45.2   3.9
   wheat.bullock.la-plant.light               4.3   3.9   3.9  43.4  16.1                                45.2   3.9
   wheat.bullock.la-plant.(heavy,january)     4.3   3.9   3.9  33.2  12.3                                45.2   3.9
   wheat.bullock.qk-harv.standard             4.3   3.9   3.9  88.1                                18.6  23.9   3.9
   wheat.bullock.qk-harv.light                4.3   3.9   3.9  74.9                                18.6  23.9   3.9
   wheat.bullock.qk-harv.(heavy,january)      4.3   3.9   3.9  57.3                                18.6  23.9   3.9
   wheat.semi-mech.standard.standard          4.3   3.9   3.9  60.1  19.5                                15.6   3.9
   wheat.semi-mech.standard.light             4.3   3.9   3.9  51.1  16.6                                15.6   3.9
   wheat.semi-mech.standard.(heavy,january)   4.3   3.9   3.9  39.1  12.7                                15.6   3.9
   wheat.semi-mech.la-plant.standard          4.3   3.9   3.9  52.8  17.2                                 8.0  11.5
   wheat.semi-mech.la-plant.light             4.3   3.9   3.9  44.9  14.6                                 8.0  11.5
   wheat.semi-mech.la-plant.(heavy,january)   4.3   3.9   3.9  34.3  11.2                                 8.0  11.5
   wheat.semi-mech.qk-harv.standard           4.3   3.9   3.9  79.6                                      15.6   3.9
   wheat.semi-mech.qk-harv.light              4.3   3.9   3.9  67.7                                      15.6   3.9
   wheat.semi-mech.qk-harv.(heavy,january)    4.3   3.9   3.9  44.0                                      15.6   3.9;

Table tr(c,t,s,w,m) 'tractor requirements in tractor hours per acre per month'
                                              jan   feb   mar   apr   may   jun   jul   aug   sep   oct   nov   dec
   basmati.semi-mech.standard.standard                                      2.7   1.7                     1.5   0.5
   cotton.semi-mech.standard.standard                                 2.5   2.2               0.5   0.5   0.5   0.5
   cotton.semi-mech.la-plant.standard                                 4.7               0.5   0.5   0.5   0.5
   berseem.semi-mech.standard.standard        1.5   1.5   1.5   1.5   1.5                     1.3   2.0   2.3   1.5
   gram.semi-mech.standard.standard                       3.7                                 1.3   0.8
   irri.semi-mech.standard.standard                                         2.6   2.2                     1.0   0.5
   sorghum.semi-mech.standard.standard                    2.5   0.5   1.8   1.0   2.2   1.5   1.0   0.5
   maiz.semi-mech.standard.standard                                               1.8   1.6                     2.5
   mustard.semi-mech.standard.standard        1.0   1.0   1.0                                       1.3   1.3   1.0
   sc-mill.semi-mech.standard.standard        6.4   6.6   5.2   4.5   4.0                                 5.8   5.8
   sc-gur.semi-mech.standard.standard         2.4   0.6   0.7                                                   0.3
   wheat.semi-mech.standard.standard                            0.7   0.8                                 7.6
   wheat.semi-mech.standard.light                               0.6   0.7                                 7.6
   wheat.semi-mech.standard.(heavy,january)                     0.5   0.5                                 7.6
   wheat.semi-mech.la-plant.standard                            0.6   0.7                                 4.0   3.6
   wheat.semi-mech.la-plant.light                               0.5   0.6                                 4.0   3.6
   wheat.semi-mech.la-plant.(heavy,january)                     0.4   0.4                                 4.0   3.6
   wheat.semi-mech.qk-harv.standard                             1.5                                       7.6
   wheat.semi-mech.qk-harv.light                                1.3                                       7.6
   wheat.semi-mech.qk-harv.(heavy,january)                      1.0                                       7.6      ;

Table wcro(c,t,m) 'water requirements in acre feet per acre for crops othe than wheat and cotton'
                                           jan   feb  mar   apr   may   jun   jul   aug   sep   oct   nov   dec
  (basmati,irri).(bullock,semi-mech)                                   .374  .842  .935  .935  .468
   berseem.(bullock,semi-mech)            .2    0.3  .45   .55   .25                     .20   .25   .15   .15
   gram.(bullock,semi-mech)               .1     .05                                     .25   .15   .10   .10
   sorghum.(bullock,semi-mech)                       .10   .15   .25   .30   .40   .30   .20
   maiz.bullock                                                        .10   .10   .30   .30   .40   .20
   maiz.semi-mech                                                            .15   .30   .30   .40   .25
   mustard.(bullock,semi-mech)            .2     .15                                           .30   .20   .20
  (sc-mill,sc-gur).(bullock,semi-mech)    .094   .14 .234  .374  .374  .327  .374  .468  .468  .281  .374  .187;

Table wcrcw(cw,t,s,w,m) 'water requirements in acre feet per acre for wheat and cotton'
                                                     jan   feb   mar   apr   may   jun   jul   aug   sep   oct   nov   dec
   cotton.(bullock,semi-mech).standard.(standard)                           .187  .281  .327  .374  .468  .327  .094
   cotton.bullock.el-plant.standard                                   .187  .094  .374  .327  .374  .468  .327  .094
   cotton.semi-mech.la-plant.standard                                       .374  .327  .374  .468  .327  .094
   wheat.bullock.(standard,qk-harv).standard        .234  .374  .468                                      .187  .187  .187
   wheat.bullock.(standard,qk-harv).light           .234  .374                                            .187  .187  .187
   wheat.bullock.(standard,qk-harv).heavy           .281                                                  .187  .187  .187
   wheat.bullock.(standard,qk-harv).january               .281                                            .187  .187  .187
   wheat.bullock.la-plant.standard                  .234  .374  .468                                            .281  .187
   wheat.bullock.la-plant.light                     .234  .374                                                  .281  .187
   wheat.bullock.la-plant.heavy                     .281                                                        .281  .187
   wheat.bullock.la-plant.january                         .281                                                  .281  .187
   wheat.semi-mech.(standard,qk-harv).standard      .234  .374  .468                                            .281  .187
   wheat.semi-mech.(standard,qk-harv).light         .234  .374                                                  .281  .187
   wheat.semi-mech.(standard,qk-harv).heavy         .281                                                        .281  .187
   wheat.semi-mech.(standard,qk-harv).january             .281                                                  .281  .187
   wheat.semi-mech.la-plant.standard                .234  .374  .468                                            .234  .234
   wheat.semi-mech.la-plant.light                   .234  .374                                                  .234  .234
   wheat.semi-mech.la-plant.heavy                   .281                                                        .234  .234
   wheat.semi-mech.la-plant.january                       .281                                                  .234  .234;

Table crio(c,t,s,w,*) 'crop inputs and outputs'
* yield and straw-yld is in maunds/acre
* cash-input  working capital required net of fertilizer.
* nit-input and pho-input  are nutrients required in kg per acre.
* adjustments for fertilizer response made below
                                                                   yield   straw-yld   cash-input   nit-input   pho-input
   basmati.(bullock,semi-mech).standard.standard                    9.00       21.00        28.83       18.80        6.90
   cotton.(bullock,semi-mech).standard.standard                    11.00                    35.34       17.40        5.20
   cotton.bullock.el-plant.standard                                11.00                    35.34       17.40        5.20
   cotton.semi-mech.la-plant.standard                              11.00                    35.34       17.40        5.20
   berseem.(bullock,semi-mech).standard.standard                  471.00                    23.20       11.50        2.07
   gram.(bullock,semi-mech).standard.standard                       5.50        8.30        25.36        1.60
   irri.bullock.standard.standard                                  16.70       27.00        28.70       25.40        7.90
   sorghum.(bullock,semi-mech).standard.standard                  192.10                   153.08        8.72         .95
   maiz.(bullock,semi-mech).standard.standard                      10.00       30.00        35.05       21.70        5.40
   mustard.(bullock,semi-mech).standard.standard                    6.00        3.00        16.80
   sc-mill.(bullock,semi-mech).standard.standard                  380.00       50.50       270.93       24.80        7.20
   sc-gur.(bullock,semi-mech).standard.standard                    30.40       50.50       270.93       24.80        7.20
   wheat.(bullock,semi-mech).(standard,qk-harv).standard           20.20       27.30        54.25       19.20        7.40
   wheat.(bullock,semi-mech).(standard,qk-harv).light              17.20       23.20        54.31       16.30        6.30
   wheat.(bullock,semi-mech).(standard,qk-harv).(heavy,january)    13.10       17.70        54.16       12.50        4.80
   wheat.(bullock,semi-mech).la-plant.standard                     17.80       24.00        54.18       16.90        6.50
   wheat.(bullock,semi-mech).la-plant.light                        15.10       20.40        54.08       14.40        5.50
   wheat.(bullock,semi-mech).la-plant.(heavy,january)              11.50       15.60        54.18       11.00        4.20;

Table weed(c,t,s,w,sea) 'weed yield from crops maunds per acre'
                                                                  kharif   rabi
  (basmati,irri,cotton).(bullock,semi-mech).standard.standard       60.0
   cotton.bullock.el-plant.standard                                 60.0
   cotton.semi-mech.la-plant.standard                               60.0
  (sc-mill,sc-gur).(bullock,semi-mech).standard.standard            60.0   30.0
   wheat.(bullock,semi-mech).(standard,qk-harv).standard                   30.0
   wheat.(bullock,semi-mech).(standard,qk-harv).light                      25.5
   wheat.(bullock,semi-mech).(standard,qk-harv).(heavy,january)            19.5
   wheat.(bullock,semi-mech).la-plant.standard                             26.4
   wheat.(bullock,semi-mech).la-plant.light                                22.4
   wheat.(bullock,semi-mech).la-plant.(heavy,january)                      17.1;

Table sconv(*,sea,c) 'tdn and dp conversion factor from crop straw'
                      basmati  cotton  berseem  gram   irri   sorghum   maiz  mustard  sc-mill sc-gur  wheat
   tdn.kharif                                     .6             .14             .5                     .5
   tdn.rabi              .5               .14          .5                 .5              .17    .17
   dp.kharif                                      .1             .017            .005                   .005
   dp.rabi               .005             .025         .005               .1              .005   .005       ;

Parameter
   land(c,t,s,w,m)   'land requirements by crop technology and month in acres'
   bpr(c,t,s,w,m)    'bullock power requirements by crop technology and month in bullock pair hours'
   labor(c,t,s,w,m)  'labor requirements by crop technology and month in man-hours'
   watreq(c,t,s,w,m) 'water requiremetns by crop technology and monthin acre feet';

land(c,t,"standard","standard",m)   = lando(c,t,m);
bpr(c,t,"standard","standard",m)    = bso(c,t,m);
labor(c,t,"standard","standard",m)  = lbco(c,t,m);
watreq(c,t,"standard","standard",m) = wcro(c,t,m);

land(cw,t,s,w,m)   = landcw(cw,t,s,w,m);
bpr(cw,t,s,w,m)    = bscw(cw,t,s,w,m);
labor(cw,t,s,w,m)  = lbcw(cw,t,s,w,m);
watreq(cw,t,s,w,m) = wcrcw(cw,t,s,w,m);

Table dev60(c,y) 'revenue deviation (1960 rupees)'
                       1960     1961     1962     1963      1964     1965     1966      1967
  (basmati,irri)      23.8     20.84   -22.76   -12.2     -15.6    -34.17   -13.56    -11.42
   cotton                              -15.97   -17.11     57.01    -6.33    -1.24    -41.37
  (berseem,sorghum)   22.38   -12.03   -28.12   -34.9     -21.6     11.7     17.49     17.99
   gram               11.03     2.04     2.25     8.35    -13.42    10.86   -19.79     -8.49
   maiz               25.19   -17.17   -19.34    24.41    -16.36    -3.21     1.21    -16.17
   mustard            -9.82   -13.58    -5.53     1.51    -12.77    16.01    23.22     16.92
  (sc-mill,sc-gur)                      52.24    18.16   -104.2    -11.1      4.4    -195.3
   wheat              12.87     7.58   -11.77    -9.44       .57    -1.58   -33.53     28.4

   +                   1968     1969     1970     1971      1972     1973     1974      1975
  (basmati,irri)      95.66   -16.06    -5.81    -6.94      2.74     1.83    -8.03      1.69
   cotton            -63.62   -42.76   -28.96    36.79     30.73    40.42    78.63    -99.33
  (berseem,sorghum)   16.45    40.79    10.16    -2.43     20.17    -1.45   -23.06    -33.62
   gram              -15.5      8.02    17.6    -14.77     -2.88     9.5    -17.55     22.73
   maiz               26.15    18.63     -.38   -28.43     12.44    -6.79   -51.35     51.18
   mustard           -21.64    14.64    34.12    22.87    -37.17   -17.17   -25.34     13.72
  (sc-mill,sc-gur)   -25.52   151.48   184.69    -9.53   -137.6    189.2    -21.49   -132.4
   wheat               5.3     -1.74    18.49   -17.75      14.1    -9.48   -20.31     18.3 ;

Scalar kl 'normalizing contant for risk';
kl = sqrt(3.1445926*card(y)/(card(y)-1)/2);

display kl;

Parameter r(g) 'risk aversion coefficient';
r(g) = 5.5;

$sTitle Price Data, Fertilizer Response and Livestock Data
Table pri(*,ps,c) 'crop price yield change working capital increase'
* crop prices are rs/lb,  working capital is in rupees
                            basmati   cotton   berseem   gram   irri   sorghum   maiz  mustard  sc-mill sc-gur  wheat
   price.76-77                1.086    1.544             .450   .485             .450     .720     .070   .650    .5
   price.(80-81,85-86)        1.5      1.86              .563   .760             .563     .9       .1    1.0      .63
   yldchng.(80-81,85-86)               1.5
   wcapchng.(80-81,85-86)     10.34  210.35       13.6  12.   10.34       13.6  5.4      8.0     19.93  19.93    6.25;

Set p1 / nitrogen, phosphorus, prdef, protein, drinvt, twinvt, trinvt, twopc, tropc, dropc /;

Table pri1(p1,ps) 'fertilizer tubewell tractor and protein prices'
* fertilizer prices are rupees per kg
* prdef price index deflator
                     76-77         80-81         85-86
   nitrogen         2.96          4.12          6.47
   phosphorus       1.80          2.86          4.26
   prdef            3.1285        4.0671        4.0671
   protein          2.25         3.0             na
   drinvt       88904        124465        124465
   twinvt        5300          7420          7420
   trinvt       15890         22246         22246
   twopc           75           154           216
   tropc           20            33.75         41
   dropc           15.872        32.59         32.59  ;

Table liveprc(q,ps) 'prices of livestock commodities'
                     76-77   80-81   85-86
   buff-milk         1.250   1.562   1.562
   catl-milk         0.950   1.188   1.188
   meat              2.5     3.125   3.125;

Table wageps(m,ps)  'wage rates rs per man hour'
                                      76-77   80-81   85-86
   (jan,feb,mar,jun,jul,aug,sep,dec)   0.75   0.975   0.975
   (apr,may,oct,nov)                   1.5    1.95    1.95 ;

* price setup for this run
Parameter
   pbc(c)      'buying price, crop comodities rupees per maund'
   psc(c)      'selling price for crop comodities rupees per maund'
   pbq(q)      'buying price, livestock comodities'
   psq(q)      'livestock commodity selling prices'
*              (rupees per litre for buff-milk and catl-milk, rupees per pound meat)
   dev(c,y)    'revenue deviation'
   wcapchng(c) 'capital change by crop'
   yldchng(c)  'yield change by crop'
   misc(p1)    'miscellenious prices'
   wage(m)     'wage rates rs per man hour by month'
   pp          'purchase price of prorein concentrate re per lbs';

Scalar
   cdrwell     'annualized cost of public drainage well rupees per well'
   trcap       'tractor capacity in tractor hours per month'         /  250   /
   twcap       'private tubewell capacity acre feet per month'       /   44.6 /
   drcap       'drainage well capacity in acre feet per year'        / 2138.4 /;

loop(psr,
   psc(c)      = pri("price",psr,c);
   pp          = pri1("protein",psr);
   psq(q)      = liveprc(q,psr);
   yldchng(c)  = pri("yldchng",psr,c);
   wcapchng(c) = pri("wcapchng",psr,c);
   misc(p1)    = pri1(p1,psr);
   dev(c,y)    = dev60(c,y)*pri1("prdef",psr);
   wage(m)     = wageps(m,psr);
);

pbc(c)  = psc(c)*1.5;
pbq(q)  = psq(q)*2.0;
cdrwell = drcap*misc("dropc") + misc("drinvt");

* fertilizer response calculations

* underlying equation for fertilizer response calculations is  y = y0*(f/f0)**gamma  where f0 anf y0
* are the fertilizer application and yield at the base price f and y are fertilizer  application for a
* given price and corresponding yield.

Parameter gammafrt(c) 'fertilizer response elasticities'
                      / basmati  .153, cotton  .0795, irri   .246
                        sc-mill  .179, sc-gur  .179 , wheat  .101 /;

Set
   fresp         / with, without /
   sfresp(fresp) / with          /;

* to incorporate fertilizer response in the model change member of set sfresp to "with" and to delete
* the fertilizer response change sfrest to "without" .

Parameter
   nit(c,t,s,w)            'nitrogen input kg per acre'
   yc(c,t,s,w)             'crop yield maunds per acre'
   straw(c,t,s,w)          'straw yield maunds per acre'
   frtfac(cfr,t,s,w,fresp) 'fertilizer response factor'
   rep3(t,s,w,c,*)         'report on yield fertilizer response factor working capital'
   cwcaptl(c,t,s,w)        'total working capital by scenario crop and tecnology in rs per acre'
   tech(c,t,s,w)           'land use indicator';

frtfac(cfr,t,s,w,"with")$(crio(cfr,t,s,w,"nit-input") <> 0) =
   min(3, (gammafrt(cfr)*psc(cfr)*crio(cfr,t,s,w,"yield")*fc /misc("nitrogen")/crio(cfr,t,s,w,"nit-input"))**(1/(1-gammafrt(cfr))));

frtfac(cfr,t,s,w,"without") = 1.0;
nit(c,t,s,w)   = crio(c,t,s,w,"nit-input");
yc(c,t,s,w)    = crio(c,t,s,w,"yield");
straw(c,t,s,w) = crio(c,t,s,w,"straw-yld");

loop(sfresp,
   nit(cfr,t,s,w)$(frtfac(cfr,t,s,w,sfresp)   <> 0) =  crio(cfr,t,s,w,"nit-input")*frtfac(cfr,t,s,w,sfresp);
   yc(cfr,t,s,w)$ (frtfac(cfr,t,s,w,sfresp)   <> 0) = (crio(cfr,t,s,w,"yield")+yldchng(cfr))*frtfac(cfr,t,s,w,sfresp)**gammafrt(cfr);
   straw(cfr,t,s,w)$(frtfac(cfr,t,s,w,sfresp) <> 0) =  crio(cfr,t,s,w,"straw-yld")*frtfac(cfr,t,s,w,sfresp)**gammafrt(cfr);
   rep3(t,s,w,cfr,"frtfac") = frtfac(cfr,t,s,w,sfresp);
);
cwcaptl(c,t,s,w)$(crio(c,t,s,w,"cash-input") <> 0) =   crio(c,t,s,w,"cash-input") + misc("nitrogen")
                                                     * nit(c,t,s,w) + misc("phosphorus")
                                                     * crio(c,t,s,w,"pho-input") + wcapchng(c);

rep3(t,s,w,c,"yield")      = yc(c,t,s,w);
rep3(t,s,w,c,"nitrogen")   = nit(c,t,s,w);
rep3(t,s,w,c,"phosphorus") = crio(c,t,s,w,"pho-input");
rep3(t,s,w,c,"straw")      = straw(c,t,s,w);
rep3(t,s,w,c,"cwcaptl")    = cwcaptl(c,t,s,w);

* nutrients availablity for livestock from fodder crops, crop straw and weeds
Parameter
   tdy(c,t,s,w,sea) 'crop tdn yield'
   dpy(c,t,s,w,sea) 'crop dp yield'
   wtd(c,t,s,w,sea) 'tdn yield from weeds'
   wdp(c,t,s,w,sea) 'dp yield from weeds'
   gfd(g,sea)       'green fodder available from grazing by season'
   gdp(g,sea)       'protein available from grazing by season';

tdy(cnf,t,s,w,sea) = sconv("tdn",sea,cnf)*straw(cnf,t,s,w)*fc;
dpy(cnf,t,s,w,sea) = sconv("dp",sea,cnf)*straw(cnf,t,s,w)*fc;
tdy(cf,t,s,w,sea)  = sconv("tdn",sea,cf)*yc(cf,t,s,w)*fc;
dpy(cf,t,s,w,sea)  = sconv("dp",sea,cf)*yc(cf,t,s,w)*fc;
wtd(c,t,s,w,sea)   = sconv("tdn","rabi","berseem")*weed(c,t,s,w,sea)*fc;
wdp(c,t,s,w,sea)   = sconv("dp","kharif","sorghum")*weed(c,t,s,w,sea)*fc;
gfd(g,sea) = 0.0;
gdp(g,sea) = 0.0;
tech(c,t,s,w) = sum(m, land(c,t,s,w,m));

Table livio(l,*) 'livestock comodity yields and output'
* milk  production is in liters/year per animal.
* meat  production is in lbs/year per animal, in case of bullock per bullock pair.
* tn    seasonal tdn requirements (lbs per season)
* pr    seasonal dp requirements  (lbs per season)
             buff-milk  catl-milk   meat     tn    pr
   buffalo         750              35.2   2300   210
   cattle                     550   26.0   1500   135
   bullock                          55.4   2800   256;

Scalar repco  'reproductive coefficient' / 1.25 /;

Parameter
   lbq(l,m)   'livestock labor requirement         (man hours per month)'
   bp(m)      'bullock pair draft power capacity  (pair-hours per month)'
   lwcaptl(l) 'livestock miscellaneous cash requirement         (rupees)'
   yq(l,q)    'livestock comodity yields'
   gr         'required proportion of green fodder in total fodder';

lbq("bullock",m) = 30.1;
lbq("buffalo",m) = 33.6;
lbq("cattle",m)  = 25.1;

bp(m)      = 96;
bp("may")  = 77;
bp("jun")  = 77;
lwcaptl(l) = 0.0;
yq(l,q)    = livio(l,q);
gr         = 0.3 ;

* report on crop input output coefficients
Set r1 / land, labor, bullock, tractor, watreq /;

Parameter
   rep1(t,s,w,c,r1,*) 'report on crop input output'
   rep2(t,s,w,c,*)    'report on crop weed fodder yield';

rep1(t,s,w,c,"land",m)    = land(c,t,s,w,m);
rep1(t,s,w,c,"bullock",m) = bpr(c,t,s,w,m);
rep1(t,s,w,c,"tractor",m) = tr(c,t,s,w,m);
rep1(t,s,w,c,"labor",m)   = labor(c,t,s,w,m);
rep1(t,s,w,c,"watreq",m)  = watreq(c,t,s,w,m);
rep2(t,s,w,c,sea)         = weed(c,t,s,w,sea);

rep1(t,s,w,c,"land","total")    = sum(m, land(c,t,s,w,m));
rep1(t,s,w,c,"bullock","total") = sum(m, bpr(c,t,s,w,m));
rep1(t,s,w,c,"tractor","total") = sum(m, tr(c,t,s,w,m));
rep1(t,s,w,c,"labor","total")   = sum(m, labor(c,t,s,w,m));
rep1(t,s,w,c,"watreq","total")  = sum(m, watreq(c,t,s,w,m));
rep2(t,s,w,c,"total")           = sum(sea, weed(c,t,s,w,sea));

display rep1, rep2, rep3;

$sTitle Population, Labor and Consumption Parameters
* population and labor availability
Set
   labs 'labor sex'       / men, women /
   int  'labor intensity' / fulltime, parttime, occasional /
   hht  'house hold type' / fh 'farm household', lh 'livestock labor household' /
   dis  'districts'       / rahim-kh /;

Parameter fpop(dis) 'farm population by districts';
fpop("rahim-kh") = 846897;

* source: agricultural census of population 1972.
Scalar lstd 'standard labor limit (hours per month)' / 200 /;

Parameter
   plab(hht)        'proportional labor participation by household type' / fh 1.0, lh .49 /
   flab(g,hht)      'available farm labor by household type thousand man-hrs per month'
   flabps(ps,g,hht) 'available farm labor by hh type thousand man-hrs per month by scenario'
   popadj(ps)       'farm population increase'
   pop(g)           'farm population in thousands'
   popps(g,ps)      'farm population in thousands by scenarios';

* popadj is estimated using increase in total farm population reported in world bank report no. 2018-pak
* over level of 1972 agricultural census .

popadj("76-77") = 1.1048;
popadj("80-81") = 1.1965;
popadj("85-86") = na;

Table alab(dis,hht,labs,int) 'available household labor by district'
                              fulltime   parttime   occasional
   rahim-kh.fh.men              217336      11446         8202
   rahim-kh.fh.women             65647      84106        10039
   rahim-kh.lh.men               20528       7207          782
   rahim-kh.lh.women              5030      11299          783;

Table labint(labs,int) 'mapping from labor class to standard intensity'
           fulltime     parttime   occasional
   men           .9           .5          .1
   women         .6           .3          .06;

Table pmap(dis,g) 'mapping from distrits to polygon'
                   poly-17+19     poly-18
   rahim-kh                .7          .3;

flabps(ps,g,hht) = sum((labs, int,dis),   alab(dis,hht,labs,int)*labint(labs,int)
                                        * pmap(dis,g)*plab(hht)*lstd*popadj(ps))/1000;

popps(g,ps) = sum(dis, fpop(dis)*pmap(dis,g)*popadj(ps))/1000;

loop(psr,
   pop(g)      = popps(g,psr);
   flab(g,hht) =  flabps(psr,g,hht);
);
display pop, flab;

* consumption parameters
Table conspc(cc,*) 'crop commodity consumption parameters'
* gamc  minimum per capita commodity consumption crops
* eec   expenditure elasticity crop
* bwc   budget weight crop
               gamc    eec     bwc
   gram       12.60    .5    .0246
   irri       62.46    .3    .0264
   maiz       20.00    .3    .0097
   mustard     1.73    .55   .0301
   sc-gur     19.50   1.65   .0270
   wheat     270.57    .4    .1995;

Table conspl(q,*) 'livestock products consumption parameters'
* gamq  minimum per capita commodity consumption livestock
* eeq   expenditure elasticity livestock
* bwq   budget weight livestock
                   gamq   eeq     bwq
   buff-milk      16.83   2.1   .0761
   catl-milk      10.76   2.1   .0507
   meat            5.92   1.1   .0174;

Scalar
   fs  'share of food in minimum per capita budget'     / .601/
   mu  'marginal propensity to consume from income'     / .931/
   sub 'minimum per capita consumption budget in rupees';

Parameter
   mpcc(cc)   'marginal propensity to consume out of income, crop'
   mpcq(q)    'marginal propensity to consume out of income, livestock'
   alphc(g,c) 'autonomous consumption of crop comodity in thousand lbs'
   alphq(g,q) 'autonomous consumption of livestock comodities in thousand lbs or liters'
   betac(g,c) 'induced consumption, crop'
   betaq(g,q) 'induced consumption, livestock';

mpcc(cc) =  conspc(cc,"eec")*conspc(cc,"bwc");
mpcq(q)  =  conspl(q,"eeq")*conspl(q,"bwq");
sub      = (sum(cc, psc(cc)*conspc(cc,"gamc")) + sum(q, psq(q)*conspl(q,"gamq")))/fs;

betac(g,cc) =  mu*mpcc(cc)/psc(cc);
betaq(g,q ) =  mu*mpcq(q)/psq(q);
alphc(g,cc) = (conspc(cc,"gamc") - betac(g,cc)*sub)*pop(g);
alphq(g,q ) = (conspl(q,"gamq")  - betaq(g,q )*sub)*pop(g);

display sub, mpcc, mpcq, betac, betaq,alphc, alphq;

$sTitle Water Data
Table wal(*,i,g,m) 'water diversions in maf'
                                  jan       feb       mar       apr       may       jun
   post.panjnad.poly-17+19   .1187200  .1352400  .1411200  .2305800  .3306800  .4071200
   post.panjnad.poly-18      .0508800  .0579600  .0604800  .0988200  .1417200  .1744800
   post.abbasia.poly-17+19   .0145800  .0235800  .0288000  .0399600  .0489600  .0547200
   expj.panjnad.poly-17+19   .0567875  .0652750  .0708750  .1423625  .2977625  .3733625
   expj.panjnad.poly-18      .0597998  .1041390  .1342239  .1075394  .1276125  .2249383
   expj.abbasia.poly-17+19   .0110250  .0140625  .0131625  .0239625  .0412875  .0516375

   +                              jul       aug       sep       oct       nov       dec
   post.panjnad.poly-17+19   .3844400  .3991400  .3924200  .3869600  .1877400  .1173200
   post.panjnad.poly-18      .1647600  .1710600  .1681800  .1658400  .0804600  .0502800
   post.abbasia.poly-17+19   .0477000  .0489600  .0477000  .0478800  .0343800  .0268200
   expj.panjnad.poly-17+19   .4321625  .4032875  .3858040  .2609250  .1193500  .0601125
   expj.panjnad.poly-18      .2229836  .2677125  .2851960  .1765128  .1056247  .0616534
   expj.abbasia.poly-17+19   .0564750  .0522000  .0486000  .0326250  .0219375  .0129375;

Parameter wdiv(i,g,m) 'water diversions used for this run';
wdiv(i,g,m) = wal("expj",i,g,m);

Table height(g,*) 'land surface and water table height from sea level in feet'
                   hs      hwr
   poly-17+19     290   280.58
   poly-18        285   279   ;

Table wpara1(*,*,m) 'evaporation rain and river seepage'
* pevap    unit pan evaporation             (in acre feet per acre)
* rain     mean monthly rainfall on total land area   (in af/month)
* rivseep  seepage to aquifer from river reaches               (af)
                            jan      feb      mar      apr      may      jun       jul       aug     sep    oct     nov     dec
   pevap.poly-17+19        .207     .317     .544     .780     .977    1.111     1.026      .913    .799   .556    .336    .439
   pevap.poly-18           .200     .314     .529     .771     .968    1.112     1.026      .913    .800   .543    .329    .457
   rain.poly-17+19       21.51    24.76    29.32    18.24    13.69    18.24    171.3     143.3    14.99    .978  10.43   11.73
   rain.poly-18           9.96    11.31    13.12     8.15     6.34     8.15     75.12     63.80    6.79    .453   4.98    5.43
   rivseep.(poly-17+19)   1000    1000     1000     1000                                  2600     4000    5000   4100    2000 ;

Table areac(g,*) 'area characteristics'
* twg     installed tubewell capacity polygon     (acre feet)
* trg1    installed tractor capacity polygon  (tractor hours)
* at      total area polygon                          (acres)
* alg     irrigated land in polgon                    (acres)
* sra     area within transportation radius of sugar mills
* scoeff  storage coefficient
                   twg     trg1        at        alg      sra   scoeff
   poly-17+19   227000   274000   1173000   862665.2   511000      .25
   poly-18               117000    543000   369713.6   392000      .25;

* following calculations are to compute water loss coefficints, we, ws and wl
Set
   u  'water transit mode from barrage to root zone'
      / lc 'link canal'   , mb 'main and branch canals', dm 'distributories and minors'
        wc 'water courses', fd 'fields'                , pg 'private tubewell'          /
   ud 'components of canal system'
      / lc 'link canal'  ,  mc 'main canal'
        bc 'branch canal',  di 'distributory', mi 'minors'   /
   dsp(g,i) 'downstream polygon canal map' / poly-18.panjnad /
   usp(g)   'upstream polygon'             / poly-17+19      /
   n        / 1*6 /;

Table wpara(g,*) 'canal data'
* wamb     wetted area main and branch canals units acres
* wadm     wetted area distributory and minor canals units acres
* factor   adjustment factor distributory and minor canals
* factor1  adjustment factor main and branch canals
* mb       main and branch canal capacity        (kaf/m for all months)
* dm       capacity of distributories and minors (kaf/m for all months)
                  wamb   wadm   factor   factor1    mb    dm
   poly-17+19   2624.0   2464   3.4366    2.4335   800   542
   poly-18        46.0    452   3.9044    9.0671   201   201;

Table wpara3(*,*,m) 'parameters for water course and  field efficiencies'
* hc     height of canal surface above sea level (in feet)
* cfw    monthly adjustment factor for operational losses
* gamma  seasonal evaporation loss factor
* zeta   seasonal seepage loss factor
                       jan     feb     mar      apr      may     jun    jul      aug     sep     oct      nov      dec
   zeta.wtrcrs      -.075   -.075   -.075   -.075    -.0025    .070   .070     .07    -.0025  -.075   -.075    -.075
   zeta.fld         -.048   -.048   -.048   -.048    -.0063    .035   .035     .035   -.0063  -.048   -.048    -.048
   gamma.wtrcrs      .025    .025    .025    .025     .0525    .08    .08      .08     .0525   .025    .025     .025
   gamma.fld         .0475   .0475   .0475   .0475    .0563    .065   .065     .065    .0563   .0475   .0475    .0475
   hc.poly-17+19     290.1   290.2   290.2    290.6    291.5   292.0  292.3    292.2   292.0   291.5    290.6    290.1
   hc.poly-18        285.1   285.1   285.2    285.6    286.4   286.9  287.2    287.0   286.9   286.4    285.5    285.1
   cfw.poly-17+19        4       4       4        4        5       6      6        6       5       4        4        4
   cfw.poly-18           4       4       4        4        5       6      6        6       5       4        4        4;

Table lw(*,g,ud) 'length(in k ft) and width of canals (in ft)'
                       lc    mc    bc     di     mi
   length.poly-17+19        300   775   2135   1450
   length.poly-18                  50    300    300
   width.poly-17+19         150    80     25      7
   width.poly-18                   45     25      7;

Parameter
   pemax(g)   'maximum of pan evaporation over the year in feet'
   alpha      'scale factor watercourse and field evaporation'
   rpe(g,m)   'ratio pan evap. to maximum pan evaporation in the year'
   wa(g)      'wetted area along canal units acres'
   phi1mb(g)  'seepage proportionality constant for main and branch canals'
   phi1dm(g)  'seepage proportionality canstant for distributories and minors'
   we(u,g,m)  'evaporation loss coefficient in transit mode'
   ws(u,g,m)  'water seepage coefficient within polygon'
   dephc(g,m) 'depth from water level in the canals to watertable'
   wl(u,g,m)  'water loss coefficient within polygon'
   wsu        'seepage coefficient upstream polygon'
   weu        'evaporation coefficient upstream polygon'
   wlu        'water loss coefficient upstream polygon'
   beta       'loss level constant'
   per        'permeability coefficient'
   nrc        'normalizing constant to convert cusecs to af per month';

alpha = 0.10;
beta  = 0.25;
per   = 2*10**(-5);
nrc   = 59.504;

* evaporation loss coefficients
pemax(g)     = smax(m,wpara1("pevap",g,m));
rpe(g,m)     = wpara1("pevap",g,m)/pemax(g);
wa(g)        = wpara(g,"wamb") + wpara(g,"wadm");
we("mb",g,m) = wpara(g,"wamb")*wpara1("pevap",g,m)*(wpara(g,"factor1")*10**(-6));
we("dm",g,m) = wpara(g,"wadm")*wpara1("pevap",g,m)*(wpara(g,"factor")*10**(-6));
we("wc",g,m) = alpha * rpe(g,m) + wpara3("gamma","wtrcrs",m);
we("fd",g,m) = alpha * rpe(g,m) + wpara3("gamma","fld",m);
we("pg",g,m) = 0.5*we("wc",g,m);

* seepage coefficients
Table funcction (*,n) 'function of canal width'
            1     2     3     4      5      6
   input    0    10    30   100    300    500
   output   0   .36   .67   .90   1.30   1.90;

Parameter
   step 'grid point for function'
   ftn  'function value'
   ftna 'interpolated function output'
   ftni 'function input'
   ftno 'function output';

ftni(n) = funcction("input",n);
ftno(n) = funcction("output",n);
step(g,ud,n) = (ftni(n) <= lw("width",g,ud)) and (ftni(n+1) > lw("width",g,ud));
ftna(g,ud,n) = (ftno(n) + (lw("width",g,ud) -ftni(n))*(ftno(n+1)-ftno(n))/(ftni(n+1)-ftni(n)))$step(g,ud,n);
ftn(g,ud)    = sum(n, ftna(g,ud,n));
dephc(g,m)   = wpara3("hc",g,m)-height(g,"hwr");
phi1mb(g)    = per*nrc*(lw("length",g,"mc")*ftn(g,"mc") + lw("length",g,"bc")*ftn(g,"bc"));
phi1dm(g)    = per*nrc*(lw("length",g,"di")*ftn(g,"di") + lw("length",g,"mi")*ftn(g,"mi"));
ws("mb",g,m) = wpara3("cfw",g,m)*phi1mb(g)*dephc(g,m)/wpara(g,"mb");
ws("dm",g,m) = wpara3("cfw",g,m)*phi1dm(g)*dephc(g,m)/wpara(g,"dm");
ws("wc",g,m) = beta - alpha*rpe(g,m) + wpara3("zeta","wtrcrs",m);
ws("pg",g,m) = 0.5*ws("wc",g,m);
ws("fd",g,m) = beta - alpha*rpe(g,m) + wpara3("zeta","fld",m);
wl(u,g,m)    = ws(u,g,m) + we(u,g,m);

wsu("mb",i,g,m)$dsp(g,i) = sum(usp, .824*ws("mb",usp,m));
weu("mb",i,g,m)$dsp(g,i) = sum(usp, .824*we("mb",usp,m));
wlu("mb",i,g,m)          = wsu("mb",i,g,m) + weu("mb",i,g,m);

* delivery efficiencies and water availability at the root zone
Parameter
   wcdeleff(g,m)  'water course delivery eff. from water course head to root zone'
   twdeleff(g,m)  'delivery eff. from private tubwell to root zone'
   cnldeleff(g,m) 'delivery eff. from canal head to root zone'
   wr(g,m)        'canal water netted to the root zone in kaf';

wcdeleff(g,m)  = (1-wl("wc",g,m))*(1-wl("fd",g,m));
twdeleff(g,m)  = (1-wl("pg",g,m))*(1-wl("fd",g,m));
cnldeleff(g,m) = (1.0-wl("mb",g,m))*(1-wl("dm",g,m))*wcdeleff(g,m);
wr(g,m)        = sum(i,wdiv(i,g,m)*(1-wlu("mb",i,g,m)))*cnldeleff(g,m)*1000;

* sub-irrigation and crop water requirements net of sub-irr. and rain fall
Set v / 1*14 /;

Table fg(*,v) 'fireman gardner curve'
               1   2   3     4     5     6     7      8       9    10    11    12   13    14
   depth       0   1   2   2.5   2.7   3.0   3.4      4   5.167    6     8     10   20   100
   pro-pan     1   1   1     1    .8    .6    .4   .317    .2    .133  .064   .02    0     0;

Parameter
   fgi(v)          'input point for fg function'
   fgo(v)          'output point for fg function'
   dep2(g)         'mean depth to watertable from canal wetted area feet'
   dep(g)          'depth to watertable feet'
   fgc1(g,v)       'initial fg function value canal fg evap'
   fgc(g)          'interpolated fg function value fg evap'
   fga1(g,v)       'initial fg function value aquifer fg evap'
   fga(g)          'interpolated fg function value aquifer fg evap'
   wn(g,c,t,s,w,m) 'crop water req net of rainfall and subirrigation'
   efr(g,m)        'effective rainfall in feet'
   subirr(g,m)     'water supplied by capillary action from the aquifer'
   drc             'proportion of rainfall which is run-off'
   the1            'proportion of nonconsumptive evapotranspiration'
   test            'interpolation tests'
   test2           'interpolation tests';

drc  = .15;
the1 = 0.6;

fgi(v)     = fg("depth",v);
fgo(v)     = fg("pro-pan",v);
dep(g)     = height(g,"hs") - height(g,"hwr");
dep2(g)    = 0.5*dep(g);
test(v,g)  = (dep(g)  >= fgi(v)) and (dep(g)  < fgi(v+1));
test2(v,g) = (dep2(g) >= fgi(v)) and (dep2(g) < fgi(v+1));
fgc1(g,v)  = (fgo(v) + (dep2(g)-fgi(v))*(fgo(v+1)-fgo(v))/(fgi(v+1)-fgi(v)))$test2(v,g);
fga1(g,v)  = (fgo(v) + (dep(g)-fgi(v))*(fgo(v+1)-fgo(v))/(fgi(v+1)-fgi(v)))$test(v,g);
fgc(g)     = sum(v, fgc1(g,v));
fga(g)     = sum(v,fga1(g,v));

efr(g,m)    = (1.0 - drc - wl("fd",g,m))*wpara1("rain",g,m)*1000/areac(g,"at");
subirr(g,m) = wpara1("pevap",g,m)*(1-the1)*fga(g);
wn(g,c,t,s,w,m) = watreq(c,t,s,w,m) - efr(g,m) - subirr(g,m);
wn(g,c,t,s,w,m)$(wn(g,c,t,s,w,m) < 0) = 0.0;

* calculations below are to compute annual underflow from neighbouring polygons qggw(g)
Set
   gn 'neighbouring polygons' / poly-17+19, poly-18, pl15, pl16, pl22                      /
   gg(g,gn)                   / poly-17+19.(poly-18,pl15,pl22),  poly-18.(poly-17+19,pl16) /;

Table wpara2(*,* ) 'height of water table above sea level and transmissibility coeff.'
*        (in  feet)              (.001 sq.ft per sec)
         poly-17+19   poly-18    pl15    pl16   pl22
   hwr       280.58       279   197.1   181.9    347
   trn       458.8        448   493     370      412;

Table dist(*,g,gn) 'node to-node length (lnn in feet) perpendicular bisector (lpb in feet)'
                      pl15      pl16      pl22    poly-18    poly-17+19
   lnn.poly-17+19   604058              309840     104148
   lnn.poly-18                601455                             104148
   lpb.poly-17+19    67696              174448     343688
   lpb.poly-18                 70300                             343688;

Parameter
    qgf(gn,g) 'underflows from polygon gn to polygon g in kaf'
    qggw(g)   'total inflow from neighbouring polygons kaf';

qgf(gn,g) = ((dist("lpb",g,gn)*((wpara2("trn",gn) + wpara2("trn",g))/2.0/1000)
            *(wpara2("hwr",gn)-wpara2("hwr",g))*nrc)/dist("lnn",g,gn))$gg(g,gn)/1000;
qggw(g)   = sum(gn$gg(g,gn), qgf(gn,g));

* annual recharge from canals, rain, govt. tubewells, evaporation from ground water pumpage by govt. tubwells
* and net rercharge

Parameter
    efs(g,m)    'effective seepage of rainfall from slack culturable land in feet'
    seepcgw(g)  'canal water seepage to the groundwater in kaf'
    seeprain(g) 'seepage of rain to groundwater in kaf'
    seepgtw(g)  'seepage from the publick tubelwells kaf'
    trseep(g)   'seepage in transit by polygon kaf'
    etgw(g)     'evapotraspiration from the groundwater along the canal and irrigated area'
    delgw(g)    'annual change in groudwater thousand acre feet'
    ntw(g)      'number of existing tubewells'
    ntr(g)      'number of existing tractors'
    gtw(g,m)    'government tubewell pumping                     (af per month)'
    gtw1(g,m)   'government tubewell pumping routed to root zone (af per month)';

efs(g,m)    = (1.0-drc)*0.2*wpara1("rain",g,m)/areac(g,"at");
etgw(g)     = sum(m, wpara1("pevap",g,m)*(wa(g)*fgc(g) + areac(g,"alg")*fga(g)))/1000;
seeprain(g) = sum(m, (1.0-drc)*((ws("fd",g,m)+0.2)*(areac(g,"at")-areac(g,"alg"))
                   + ws("fd",g,m)*areac(g,"alg"))*wpara1("rain",g,m)/areac(g,"at"));
trseep(usp) = sum((m,i,g), wdiv(i,g,m)*1000*wsu("mb",i,g,m)$dsp(g,i));
seepcgw(g)  = sum(m, sum(i, wdiv(i,g,m)*1000*(1-wlu("mb",i,g,m))*( ws("mb",g,m) + (1-wl("mb",g,m))*(ws("dm",g,m)
                   + (1-wl("dm",g,m))*(ws("wc",g,m) + (1-wl("wc",g,m))*ws("fd",g,m)))))) + trseep(g);

ntw(g)     = 0.0;
ntr(g)     = 0.0;
gtw(g,m)   = 0.0;
delgw(g)   = 0.0;
gtw1(g,m)  = gtw(g,m)*wcdeleff(g,m) ;
seepgtw(g) = sum(m, gtw(g,m)*(ws("wc",g,m) + (1-wl("wc",g,m))*ws("fd",g,m)))/1000;

display wsu, weu, wlu, trseep;

* report in effficiencies, and seepage to groundwater
Parameter
   rep4(g,*,u,*) 'report on water losses'
   rep5(g,*,*)   'water requirements and avalability'
   rep6(g,*)     'groundwater inflows and outflows';

rep4(g,"seepage",u,m) = ws(u,g,m);
rep4(g,"evap",u,m)    = we(u,g,m);
rep4(g,"total",u,m)   = wl(u,g,m);
rep5(g,"rain",m)      = efr(g,m);

rep5(g,"rain","total")     = sum(m, efr(g,m));
rep5(g,"subirr",m)         = subirr(g,m);
rep5(g,"subirr","total")   = sum(m, subirr(g,m));
rep5(g,"cnldeleff",m)      = cnldeleff(g,m);
rep5(g,"canal",m)          = sum(i, wdiv(i,g,m)*1000);
rep5(g,"canal","total")    = sum(m, rep5(g,"canal",m));
rep5(g,"canal-rt",m)       = wr(g,m);
rep5(g,"canal-rt","total") = sum(m, wr(g,m));

rep6(g,"seepcgw")   = seepcgw(g);
rep6(g,"seeprain")  = seeprain(g);
rep6(g,"rivseep")   = sum(m, wpara1("rivseep",g,m))/1000;
rep6(g,"seepgtw")   = seepgtw(g);
rep6(g,"underflow") = qggw(g);
rep6(g,"gtw")       = sum(m, gtw(g,m) )/1000;
rep6(g,"etgw  ")    = etgw(g);
rep6(g,"net-seep")  =   seepcgw(g) + seeprain(g) + rep6(g,"rivseep") + rep6(g,"seepgtw")
                      + qggw(g) - rep6(g,"gtw") - etgw(g);
display rep4;

display "  rain        rain use by crops in feet "
        "  subirr      water used by crop by cappilary action from the gw acre feet/acre "
        "  cnldeleff   canal delivery efficiency from canal head to the root zone "
        "              excluding the losses in transit through other polygons     "
        "  canal       canal water availability at canal head kaf "
        "  canal-rt    canal water availability at root zone kaf ";
display rep5;
display "  seepcgw     seepage from canal water to the groundwater (gw) in kaf "
        "  seeprain    seepage from rain to gw kaf              "
        "  rivseep     seepage from rivers to gw in kaf"
        "  seepgtw     seepage from public tubewells to gw in kaf "
        "  underflow   flow  from neighbouring polygons to gw kaf "
        "  etgw        evaporation from groundwater kaf          "
        "  gtw         public tubewell pumpage kaf  "
        "  net-seep    net seeapge to gw from sources mentioned above kaf ";
display rep6;

$sTitle Equations and Variables
Variable
   yfa(g)        'farm income thousand rupees'
   yva(g)        'normal farm income thousand rupees'
   mad(g)        'mean absolute revenue deviation of cropping pattern thousand rupees'
   dr(g)         'public drainage from groundwater in kaf'
   inj(g)        'public injection to groundwater in kaf'
   tw(g,m)       'private tubewell water use in month m kaf'
   ts(g,m)       'private tractor services use m month thousand hrs'
   scc(g,c)      'sales of crop commodities c in thousand lbs'
   ccc(g,c)      'on farm consumption of crops commodities c in thousand lbs'
   pcc(g,c)      'farm purchases of crop commodities c in thousand lbs'
   slc(g,q)      'sales of livestock commodities q  in thousand lbs or liters'
   clc(g,q)      'farm consumption of livestock commodities q in thousand lbs or liter'
   plc(g,q)      'farm purchases of livestock commodity q in thousand lbs or liters'
   acost(g)      'farm cost  in thousand rupees'
*                 operating cost + annualized cost of capital services
   x(g,c,t,s,w)  'area of crop c in thousand acres using a technology'
   xca(g,c)      'area of crop c by polygon'
   animal(g,l)   'production of livestock type l thousands'
*                 thouands of bullocks pairs, cattle or buffaloes
   ppc           'purchases of protein concentrates thousand lbs'
   esl(g,m)      'employment of seasonal labor in month m  in thousand man-hrs'
   itw(g)        'investment in increased private tubewell capacity kaf per month'
   itr(g)        'investment in increased tractor capacity'
*                 thousand tractor-hours per month
   efl(g,m)      'employment of family labor in month m in k man-hrs'
   pdev(g,y)     'positive revenue deviation in thousand rupees'
   ndev(g,y)     'negative revenue deviation in thousand rupees'
   utl           'utility of income million rupees'
   slkland(g,m)  'slack land in k acres'
   slkwater(g,m) 'slack water in kaf';

Positive Variable pdev , ndev, mad, dr    , inj, tw , ts , scc, ccc, pcc, slc, clc, plc
                  acost, x   , xca, animal, ppc, esl, itw, itr, efl, slkwater, slkland;

Equation
   objt  'objective function'
   inbl  'income accounting balance'
   nfin  'normal farm income constraint'
   ddev  'definition of revenue deviation year'
   cost  'annual farm cost'
   cmbc  'commodity balances for crops'
   cmbq  'commodity balances for livestock'
   cblc  'farm consumption of crop commodity'
   cblq  'farm consumption of livestock commodity'
   fdsp  'seasonal maintenance of fodder supplies'
   slsk  'protein requirements of livestock in season s'
   sgfd  'requirement for green fodder by livestock'
   bupw  'bullock draft power constraint'
   buca  'bullock reproduction constraint'
   trpw  'tractor draft power balance'
   trcp  'tractor capacity constraint'
   labr  'labor requirements constraint'
   landc 'land constraint'
   cacr  'crop acreage'
   watr  'water requirements for cropping activities'
   tbcp  'tubewell capacity constraint'
   gwbl  'annual groundwater balance';

* objective function
objt..      utl    =e=   sum(g, yfa(g) - r(g)*kl*sum(y, pdev(g,y) + ndev(g,y))/card(y)
                              - cdrwell*(dr(g) + inj(g))$sg(g)/drcap)/1000;

inbl(g)..   yfa(g) =e=   sum(c$cnf(c), psc(c)*scc(g,c)) + sum(c$cc(c), psc(c)*ccc(g,c) - pbc(c)*pcc(g,c))
                       + sum(q, psq(q)*slc(g,q) - pbq(q)*plc(g,q) + psq(q)*clc(g,q)) -acost(g);

nfin(g)..   yva(g) =e= yfa(g) - sum(c$cc(c), (pbc(c)-psc(c))*pcc(g,c)) + sum(q, (pbq(q)-psq(q))*plc(g,q));

ddev(g,y).. sum(c, dev(c,y)*xca(g,c)) =e= pdev(g,y) - ndev(g,y);

cost(g)..   acost(g) =e=   sum((c,t,s,w), cwcaptl(c,t,s,w)*x(g,c,t,s,w)$tech(c,t,s,w))
                         + sum(m, misc("twopc")*tw(g,m) + misc("tropc")*ts(g,m)+ wage(m)*esl(g,m))
                         + sum(l, lwcaptl(l)*animal(g,l)) + sum(sea, pp*ppc(g,sea))
                         + misc("twinvt")*(itw(g) + ntw(g))$fsg(g) + misc("trinvt")*(itr(g) + ntr(g));

cmbc(g,c)$cnf(c).. sum((t,s,w), yc(c,t,s,w)*fc*x(g,c,t,s,w)$tech(c,t,s,w))-scc(g,c)-(ccc(g,c)-pcc(g,c))$cc(c) =e= 0;

cmbq(g,q)..        sum(l, yq(l,q)*animal(g,l)) - slc(g,q) - clc(g,q) + plc(g,q) =e= 0.0;

cblc(g,c)$cc(c)..  ccc(g,c) =g= alphc(g,c) + betac(g,c)*yva(g);

cblq(g,q)..        clc(g,q) =g= alphq(g,q) + betaq(g,q)*yva(g);

fdsp(g,sea)..     sum(l, livio(l,"tn")*animal(g,l))
              =l= gfd(g,sea) + sum((c,t,s,w), (tdy(c,t,s,w,sea)+wtd(c,t,s,w,sea))*x(g,c,t,s,w)$tech(c,t,s,w));

slsk(g,sea)..     sum(l, livio(l,"pr")*animal(g,l))
              =l= ppc(g,sea) + gdp(g,sea) + sum((c,t,s,w), (dpy(c,t,s,w,sea)+wdp(c,t,s,w,sea))*x(g,c,t,s,w)$tech(c,t,s,w));

sgfd(g,sea)..     gr*sum(l, livio(l,"tn")*animal(g,l))
              =l= gfd(g,sea) + sum( (cf,t,s,w), tdy(cf,t,s,w,sea)*x(g,cf,t,s,w)$tech(cf,t,s,w))
               +  sum((c,t,s,w), wtd(c,t,s,w,sea)*x(g,c,t,s,w)$tech(c,t,s,w));

bupw(g,m)..   sum((c,t,s,w), bpr(c,t,s,w,m)*x(g,c,t,s,w)$tech(c,t,s,w)) =l= bp(m)*animal(g,"bullock");

buca(g)..     animal(g,"bullock") =l= repco*animal(g,"cattle");

trpw(g,m)..   sum((c,t,s,w), tr(c,t,s,w,m)*x(g,c,t,s,w)$tech(c,t,s,w)) - ts(g,m) =e= 0.0;

trcp(g,m)..   ts(g,m) =l= areac(g,"trg1")/1000 + trcap*itr(g);

labr(g,m)..   sum((c,t,s,w), labor(c,t,s,w,m)*x(g,c,t,s,w)$tech(c,t,s,w))
            + sum(l, lbq(l,m)*animal(g,l)) =e= efl(g,m) + esl(g,m);

landc(g,m)..  sum((c,t,s,w), land(c,t,s,w,m)*x(g,c,t,s,w)$tech(c,t,s,w)) + slkland(g,m) =e= areac(g,"alg")/1000;

cacr(g,c)..   xca(g,c) =e= sum((t,s,w) , x(g,c,t,s,w)$tech(c,t,s,w));

watr(g,m)..       sum((c,t,s,w), wn(g,c,t,s,w,m)*x(g,c,t,s,w)$tech(c,t,s,w)) + slkwater(g,m)
              =e= wr(g,m) + twdeleff(g,m)*tw(g,m)$fsg(g) + gtw1(g,m)/1000;

tbcp(g,m)$fsg(g).. tw(g,m) =l= areac(g, "twg")/1000 + twcap*itw(g);

* annual groundwater balance --- additional seepage of rain from slack land + seepage from tubewell water +
* seepage from canal water + underflow + seepage of rain from unirrigated area + seepage from rain in
* irrigated area + seepage from rice fields and river reaches  = e = private and government tubewell pumpage
* + evapotranspiration along canals and aquifer + drainage/injection + change storage

gwbl(g)..     sum(m, efs(g,m)*slkland(g,m) + (ws("pg",g,m) + (1-wl("pg",g,m))*ws("fd",g,m))*tw(g,m)$fsg(g) + wpara1("rivseep",g,m)/1000)
           +  sum(cri, xca(g,cri))*1.5 + seepcgw(g) + seepgtw(g) + qggw(g) + seeprain(g)
          =e= sum(m, tw(g,m)$fsg(g) + gtw(g,m)/1000) + etgw(g) + (dr(g) - inj(g))$sg(g) + delgw(g);

* solve the model using lp
xca.up(g,"sc-mill") = areac(g,"sra")/1000;
efl.up(g,m) = flab(g,"fh");
esl.up(g,m) = flab(g,"lh");

* itr.fx(g) = 0.0;
* itw.fx(g) = 0.0;

Model indus1 / objt , inbl, nfin, ddev, cost, cmbc, cmbq, cblc
               cblq , fdsp, slsk, sgfd, bupw, buca, trpw, labr
               landc, cacr, watr, tbcp, gwbl, trcp             /;

* indus2 is model without groundwater balance constraint
Model indus2 / objt, inbl,  nfin, ddev, cost, cmbc, cmbq, cblc
               cblq, fdsp,  slsk, sgfd, bupw, buca, trpw, trcp
               labr, landc, cacr, watr, tbcp                   /;

solve indus1 maximizing utl using lp;